#include "Flux.hh"
#include "PremModel.h"
#include "PMNS_Fast.h"
#include "km3net-dataformat/definitions/w2list_gseagen.hh"

struct nu_summary
{
	double bx;             // Bjorken x
	double by;             // 'Bjorken' y
	double ichan;          // ichan from the mc production 
	double cc;             // cc flag from the mc production
	double atm_flux_numu;  // flux of (anti) numu at zenith and energy of the neutrino [GeV-1 s-1 sr-1 m-2]
	double atm_flux_nue;   // flux of (anti) nue at zenith and energy of the neutrino [GeV-1 s-1 sr-1 m-2]
	double posc_from_numu; // oscillation probabiliy from numu into the flavour of the neutrino 
	double posc_from_nue;  // oscillation probabiliy from nue  into the flavour of the neutrino 
	double atm_flux;       // total flux of flavour the neutrino after oscillations [GeV-1 s-1 sr-1 m-2]	
        double path_length;    // Neutrino path length according to PREM model (km)
};




OscProb::PMNS_Fast& oscprob_pmns()
{
	static OscProb::PMNS_Fast o;
	return o;
}



inline double gusr( AAObject& o, string key , double defa = -1 ) 
{
	if ( ! o.haveusr(key) ) {
		cout << "ERROR: did not find usr key " << key << endl;
		return defa;
	}
	return o.getusr(key);
}



inline bool summarize_nu ( Evt& evt, nu_summary& s )
{
	Trk* nup = evt.neutrino();
	if (!nup)
	{
		print("ERROR: failed to find neutrino");
		return false;
	}

	Trk& nu = *nup;

	// check if they are stored in w2list or usr
	if (evt.w2list.size()>9) {
        s.bx    = evt.w2list[W2LIST_GSEAGEN_BX];
        s.by    = evt.w2list[W2LIST_GSEAGEN_BY];
        s.ichan = evt.w2list[W2LIST_GSEAGEN_ICHAN];
        s.cc    = evt.w2list[W2LIST_GSEAGEN_CC];
	}
	else {
		// failure of these is allowed (default = -1 )
        s.bx    = gusr( nu , "bx" );
        s.by    = gusr( nu , "by" );
        s.ichan = gusr( nu , "ichan" );
        s.cc    = gusr( nu , "cc" );
	}
	// It's either CC, or not CC, right ?
	if      (s.cc == 2) {s.cc = 1;}
	else if (s.cc == 3) {s.cc = 0;}

    static Flux_Atmospheric F;

    int sign = nu.type > 0 ? 1 : -1;

    s.atm_flux_numu = F.dNdEdOmega( sign * 14, log10( nu.E ), nu.dir.z );
    s.atm_flux_nue  = F.dNdEdOmega( sign * 12, log10( nu.E ), nu.dir.z );

    OscProb::PMNS_Fast& pmns = oscprob_pmns();
    static OscProb::PremModel prem;
    
    // static auto init = [] {
    //  	cout << " init oscprob" << endl;
    // pmns.SetDeltaMsqrs(7.400e-5, 2.494e-3);
    // pmns.SetMix(5.868e-1,8.238e-1,1.491e-1,234.*TMath::Pi()/180.);
    //     return true;
    // 	}(); // javascript ftw!


    

    int op_flav = ( abs(nu.type)-12 ) / 2; // 0 = nue, 1 = numu, 2 = nutau
    prem.FillPath( -1.*nu.dir.z );
    pmns.SetPath( prem.GetNuPath() );
    s.path_length = prem.GetTotalL( -1.*nu.dir.z);
    pmns.SetIsNuBar(sign<0);
    s.posc_from_nue = pmns.Prob( 0, op_flav , nu.E);
    s.posc_from_numu= pmns.Prob( 1, op_flav , nu.E);

    if (s.cc) {
      s.atm_flux = s.posc_from_nue * s.atm_flux_nue + s.posc_from_numu * s.atm_flux_numu;
    } else {
      s.atm_flux = s.atm_flux_nue + s.atm_flux_numu;
    }
 
    return true;
}