#ifndef NUFROMGAMMAHH
#define NUFROMGAMMAHH


#include <iostream>
#include <cstdio>
#include <fstream>
#include <sstream>
#include <vector>
#include <string>
#include <cmath>
#pragma GCC diagnostic push
#pragma GCC diagnostic ignored "-Wall"
#include "TROOT.h"
#include "TChain.h"
#include "TFile.h"
#include "TTree.h"
#include "TBranch.h"
#include "TH1.h"
#include "TF1.h"
#include "TCanvas.h"
#include "TGraph.h"
#include "TLegend.h"
#include "TSpline.h"
#pragma GCC diagnostic pop


/***************************************************************************************************
	
 From emission gamma spectrum (flux) calculata the neutrino emission spectrum following formulas (30-33)
 F.L. Villante, F. Vissani
 https://arxiv.org/abs/0807.4151v1

 R. Coniglione 3-9-2011 
 B. Caiffi, M. Sanguineti, V, Kulikovskiy 28-12-2020 adaptation for a aanet framework
***************************************************************************************************
Analytic gamma spectrum (input)
* spettro proporzionale VelaJunior HESS 2016 observation  
* 32.2e-11*pow(x[0]/1e3),-1.81)*exp(-x[0]/6.7e3);  // in GeV-1 m^-2 s^-1
* eHWC J1825-134
* "2.12e-12*pow(x[0]/1e4,-2.12)*exp(-x[0]/6.1e04) 
* eHWC J1907-063
* "0.95e-12*pow(x[0]/1e4,-2.46-0.11*ln(x[0]/1.e4))
* eHWC J2019-368
* "0.45e-12*pow(x[0]/1e4,-2.08-0.26*ln(x[0]/1.e4))  
*
*
* Usage: 
* nufromgamma nuflux(gammafluxeq); // where gammafluxeq are the strings shown above 
* nuflux.F_nu(energy);             //get flux[GeV^-1 m^-2 s^-1] at energy [GeV]
* nuflux.F_antnu(energy);          //antinu flux
* nuflux.F_nu_tot(energy);         //nu+antinu flux
*
*/ 

struct nufromgamma : public TObject {
       TF1 f_gamma;
       TSpline3 f_nu_tot;
       TSpline3 f_nu;
       TSpline3 f_antinu;
       //std::string typeOFconversion;
       int typeOFconversion;

       //failed to work in python... no idea why
       //static constexpr double logemin = 2;
       //static constexpr double logemax = 8;
       //static constexpr unsigned int nbinloge = 1e4;

       nufromgamma(): typeOFconversion(1) {
	       //cout << "DEBUG: calling empty constructor" << endl;
	       	}
  nufromgamma(const std::string& gammafluxeq,int IndexModConv) {
				 f_gamma = TF1("fgammaflux",gammafluxeq.c_str(),100.,10000000.);
                                 typeOFconversion = IndexModConv;
				   f_nu= to_spline("numu", IndexModConv);
				   f_antinu= to_spline("antinumu", IndexModConv);
				   f_nu_tot= to_spline("tot", IndexModConv);
       }

	     nufromgamma(const nufromgamma& master){
	       //cout << "calling copy constructor of nufromgamma" << endl;
         f_gamma = TF1("fgammaflux",master.f_gamma.GetExpFormula().Data(),100.,10000000.);
	 typeOFconversion = master.typeOFconversion;
				   f_nu= to_spline("numu", master.typeOFconversion );
				   f_antinu= to_spline("antinumu",  master.typeOFconversion );
				   f_nu_tot= to_spline("tot",  master.typeOFconversion );
         }	

	     nufromgamma& operator= (const nufromgamma& master) {
         //cout << "calling copy assignment operator of nufromgamma" << endl;
	 this->typeOFconversion = master.typeOFconversion;
         this->f_gamma = TF1("fgammaflux",master.f_gamma.GetExpFormula().Data(),100.,10000000.);
				   this->f_nu= to_spline("numu", master.typeOFconversion) ;
				   this->f_antinu= to_spline("antinumu", master.typeOFconversion) ;
				   this->f_nu_tot= to_spline("tot", master.typeOFconversion) ;
	       return *this;
	     }

       nufromgamma (nufromgamma&& master) {
	  //cout << "calling copy by reference" << endl;
         f_gamma = TF1("fgammaflux",master.f_gamma.GetExpFormula().Data(),100.,10000000.);
	 typeOFconversion = master.typeOFconversion;
				 f_nu= to_spline("numu", master.typeOFconversion );
				 f_antinu= to_spline("antinumu", master.typeOFconversion );
				 f_nu_tot= to_spline("tot", master.typeOFconversion );
       }

  
  TSpline3 to_spline(const std::string& type, int typeOFconversion )
	{
	  using namespace std;
	  
           const double logemin = 2;
           const double logemax = 8;
           const double nbinloge = 1e4;

           vector<double> xs_up;
	   vector<double> ys_up;

	   double E;

	   for(double loge=logemin;loge<logemax;loge+=(logemax-logemin)/nbinloge){
	    E=pow(10,loge);
	    xs_up.push_back(loge);
            if (typeOFconversion == 1){
	      if (type=="tot")   ys_up.push_back( F_nu_tot_up(E));
	      if (type=="numu")   ys_up.push_back( F_nu_up(E));
	      if (type=="antinumu")   ys_up.push_back( F_antinu_up(E));
            }
            if (typeOFconversion == 2){
	      if (type=="tot")   ys_up.push_back( F_nu_tot(E));
	      if (type=="numu")   ys_up.push_back( F_nu(E));
	      if (type=="antinumu")   ys_up.push_back( F_antinu(E));
            }
           if (typeOFconversion == 3) {
		   if (type=="tot")   ys_up.push_back( F_nu_KAB(E));
		   else ys_up.push_back(F_nu_KAB(E)/2.);
	   }
    	   }

	   xs_up.push_back( logemax );
	   ys_up.push_back( 0 );

	   TSpline3 spline_up(type.c_str(), &(xs_up[0]), &(ys_up[0]), xs_up.size() , "e1b1"  );
	   spline_up.SetLineColor(4);
	   spline_up.SetLineWidth(2);
           return spline_up;
	}




 
  
  


        // ***************************************

        // neutrno kernel (eq 31)
        Double_t K_nu(Double_t x){
            double k_nu=0;
            float r_k =0.0458;  
            float r_pi=0.573;
            if(x <= r_k)               k_nu = x*x*(15.34-28.93*x);
            if((r_k < x) && (x< r_pi)) k_nu = 0.0165 + 0.1193*x + 3.747*x*x - 3.981*x*x*x;
            if(x >= r_pi)              k_nu = (1.-x)*(1.-x)*(-0.6698+ 6.588*x);
            return k_nu;
        }

        // antineutrino kernel (eq 33)
        Double_t K_antinu(Double_t x){
            double k_antinu = 0;
            float r_k =0.0458;  
            float r_pi=0.573;
            if(x <= r_k)               k_antinu = x*x*(18.48-25.33*x);
            if((r_k < x) && (x< r_pi)) k_antinu = 0.0251 + 0.0826*x + 3.697*x*x - 3.548*x*x*x;
            if(x >= r_pi)              k_antinu = (1.-x)*(1.-x)*(0.0351+ 5.864*x);
            return k_antinu;
        }
  
        // neutrno kernel update 
        Double_t K_nu_up(Double_t x){
            double k_nu_up=0;
            float r_k =0.0458;  
            float r_pi=0.573;
            if(x <= r_k)               k_nu_up = x*x*(22.405-59.43*x);
            if((r_k < x) && (x< r_pi)) k_nu_up = 0.0146 + 0.2838*x + 6.89*x*x - 8.117*x*x*x;
            if(x >= r_pi)              k_nu_up = (1.-x)*(1.-x)*(-2.693+ 13.436*x);
            return k_nu_up;
        }

        // antineutrino kernel update
        Double_t K_antinu_up(Double_t x){
            double k_antinu_up = 0;
            float r_k =0.0458;  
            float r_pi=0.573;
            if(x <= r_k)               k_antinu_up = x*x*(23.791-46.799*x);
            if((r_k < x) && (x< r_pi)) k_antinu_up = 0.0244 + 0.1965*x + 6.042*x*x - 6.602*x*x*x;
            if(x >= r_pi)              k_antinu_up = (1.-x)*(1.-x)*(-1.4267+ 10.904*x);
            return k_antinu_up;
        }

        //  KAB electric neutrino kernel
        Double_t K_e_KAB(Double_t x){
            double k_e_KAB=0;  
            float r= 0.573;
            if(x< r)  k_e_KAB =  2/(3*(1.-r)*(1.-r))*((1.-r)*(6.-7.*r+11.*r*r-4.*r*r*r)+6.*r*log(r))+2.*(r-x)/(3.*r*r)*(7.*r*r - 4.*r*r*r+7.*x*r-4.*x*r*r-2.*x*x-4.*x*x*r);
            if(x >= r)  k_e_KAB = 2*(1-x)/(3*(1-r)*(1-r))*((1.-x)*(6*(1-x)*(1-x)+r*(5+5*x-4*x*x))+6*r*log(x));
            return k_e_KAB;
        }

        //  KAB muon neutrino kernel
        Double_t K_mu_KAB(Double_t x){
            double k_mu_KAB = 0; 
            float r=0.573;
            if(x< r)   k_mu_KAB = (3-2*r)/(9*(1-r)*(1-r))*(9*r*r-6*log(r)-4*r*r*r-5)+(1+2*r)*(r-x)/(9*r*r)*(9*(r+x)-4*(r*r+r*x+x*x));
            if(x >= r) k_mu_KAB = (3-2*r)/(9*(1-r)*(1-r))*(9*x*x-6*log(x)-4*x*x*x-5);
            return k_mu_KAB;
        }
  
        // nu integral (last part in eq 30)

        Double_t Int_nu(Double_t Energy){
            Double_t x_start = 0.;
            Double_t x_end   = 1.;
            Double_t dx = 0.001;
            Double_t Integral_nu = 0.;
            Double_t x= x_start + dx/2.;
            while (x < x_end){
                Integral_nu += K_nu(x) * f_gamma.Eval(Energy/x) * (1./x) * dx;
                x += dx;
            }
            return Integral_nu;
        }

        // nbar integral (last part in eq 32)

        Double_t Int_antinu(Double_t Energy){
            Double_t x_start = 0.;
            Double_t x_end   = 1.;
            Double_t dx = 0.001;
            Double_t Integral_antinu = 0.;
            Double_t x= x_start + dx/2.;
            while (x < x_end){
                Integral_antinu += K_antinu(x) * f_gamma.Eval(Energy/x) * (1./x) * dx;
                x += dx;
            }
            return Integral_antinu;
        }
  //integrals for udated Villante Visani method
        Double_t Int_nu_up(Double_t Energy){
            Double_t x_start = 0.;
            Double_t x_end   = 1.;
            Double_t dx = 0.001;
            Double_t Integral_nu_up = 0.;
            Double_t x= x_start + dx/2.;
            while (x < x_end){
                Integral_nu_up += K_nu_up(x) * f_gamma.Eval(Energy/x) * (1./x) * dx;
                x += dx;
            }
            return Integral_nu_up;
        }


        Double_t Int_antinu_up(Double_t Energy){
            Double_t x_start = 0.;
            Double_t x_end   = 1.;
            Double_t dx = 0.001;
            Double_t Integral_antinu_up = 0.;
            Double_t x= x_start + dx/2.;
            while (x < x_end){
                Integral_antinu_up += K_antinu_up(x) * f_gamma.Eval(Energy/x) * (1./x) * dx;
                x += dx;
            }
            return Integral_antinu_up;
        }

  //derivative and integral for Kelner Aharonian Bugayov method 

        Double_t derr_f(Double_t Energy){
	    Double_t derr = (f_gamma.Eval(Energy+Energy/1000)-f_gamma.Eval(Energy-Energy/1000))/(2*Energy/1000);
            return derr;
        }


        Double_t Int_nu_mu_KAB(Double_t Energy){
            Double_t x_start = 0.;
            Double_t x_end   = 1.;
            Double_t dx = 0.001;
            Double_t Integral_nu_mu_KAB = 0.;
            Double_t x= x_start + dx/2.;
            while (x < x_end){
	      Integral_nu_mu_KAB +=  (-1)*Energy*K_mu_KAB(x) * derr_f(Energy/x) * (1./(x*x)) * dx;
                x += dx;
            }
            return Integral_nu_mu_KAB;
        }


        Double_t Int_nu_e_KAB(Double_t Energy){
            Double_t x_start = 0.;
            Double_t x_end   = 1.;
            Double_t dx = 0.001;
            Double_t Integral_nu_e_KAB = 0.;
            Double_t x= x_start + dx/2.;
            while (x < x_end){
	      Integral_nu_e_KAB +=  (-1)*Energy*K_e_KAB(x) * derr_f(Energy/x) * (1./(x*x)) * dx;
                x += dx;
            }
            return Integral_nu_e_KAB;
        }


  
        //nu completed (eq 30) 
        Double_t F_nu(Double_t Energy){      
            float r_pi=0.573;
            float r_k =0.0458;
            Double_t Phi_nu = 0.380*f_gamma.Eval(Energy/(1.-r_pi))+0.0130*f_gamma.Eval(Energy/(1.-r_k)) + Int_nu(Energy);              
            return Phi_nu; 
        }


        Double_t F_nuS(double loge) {
          const double logemin = 2;
          const double logemax = 8;
          if (loge < logemin || loge > logemax) return 0; //
          double val = f_nu.Eval(loge);
          if(val<0 ) return 0.0; //spline fault
          return f_nu.Eval(loge);
        }

        //nbar completed (eq 33)  
        Double_t F_antinu(Double_t Energy){      
            float r_pi=0.573;
            float r_k =0.0458;
            Double_t Phi_antinu = 0.278*f_gamma.Eval(Energy/(1.-r_pi))+0.0090*f_gamma.Eval(Energy/(1.-r_k)) + Int_antinu(Energy);
            return Phi_antinu; 
        }


        Double_t F_antinuS(double loge) {
          const double logemin = 2;
          const double logemax = 8;
          if (loge < logemin || loge > logemax) return 0; //
          double val = f_antinu.Eval(loge);
          if(val<0 ) return 0.0; //spline fault
          return f_antinu.Eval(loge);
        }

        // nu+nbar

        Double_t F_nu_tot(Double_t Energy){
            Double_t Phi_nu_tot = F_nu(Energy) + F_antinu(Energy);
	    Phi_nu_tot = Phi_nu_tot;
            return Phi_nu_tot;
        }
        Double_t F_nu_totS(double loge) {
          const double logemin = 2;
          const double logemax = 8;
          if (loge < logemin || loge > logemax) return 0; //
          double val = f_nu_tot.Eval(loge);
          if(val<0 ) return 0.0; //spline fault
          return f_nu_tot.Eval(loge);
        }


        Double_t F_nu_up(Double_t Energy){      
            float r_pi=0.573;
            float r_k =0.0458;
            Double_t Phi_nu_up = 0.3935*f_gamma.Eval(Energy/(1.-r_pi))+0.013455*f_gamma.Eval(Energy/(1.-r_k)) + Int_nu_up(Energy);              
            return Phi_nu_up; 
        }


        //nbar completed (eq 33)  
        Double_t F_antinu_up(Double_t Energy){      
            float r_pi=0.573;
            float r_k =0.0458;
            Double_t Phi_antinu_up = 0.2879*f_gamma.Eval(Energy/(1.-r_pi))+0.009345*f_gamma.Eval(Energy/(1.-r_k)) + Int_antinu_up(Energy);
            return Phi_antinu_up; 
        }


        // nu+nbar

        Double_t F_nu_tot_up(Double_t Energy){
            Double_t Phi_nu_tot_up = F_nu_up(Energy) + F_antinu_up(Energy);
            return Phi_nu_tot_up;
        }

        Double_t F_nu_KAB(Double_t Energy){      
            float r=0.573;
            float P_mumu = 0.36;
	    float P_emu = 0.26;
            Double_t Phi_nu_KAB = P_mumu*((1./(1-r))*f_gamma.Eval(Energy/(1.-r))+ Int_nu_mu_KAB(Energy))+P_emu*Int_nu_e_KAB(Energy);              
            return Phi_nu_KAB; 
        }
  
  
  TSpline3 get_f_nu(){
    return f_nu;

  }
  
  TSpline3 get_f_antinu(){
    return f_antinu;

  }


  TSpline3 get_f_nu_tot(){
    return f_nu_tot;

  }


  
     
//****************** Helpers **************************
    //gamma *E**2
        Double_t F_gammaE(Double_t Energy){
          Double_t Phi_gammaE =  pow(Energy,2) * f_gamma.Eval(Energy);
          return Phi_gammaE;
        }
    // nu *E**2
        Double_t F_nuE(Double_t Energy){
          Double_t Phi_nuE =  pow(Energy,2) * F_nu(Energy);
          return Phi_nuE;
        }
    // nbar *E**2
        Double_t F_antinuE(Double_t Energy){
          Double_t Phi_antinuE =  pow(Energy,2) * F_antinu(Energy);
          return Phi_antinuE;
        }
    // nu+nbar *E**2
        Double_t F_nu_totE(Double_t Energy){
          Double_t Phi_nu_totE =  pow(Energy,2) * F_nu_tot(Energy);
          return Phi_nu_totE;
        }
    // gamma/nu_tot ratio 
        Double_t F_ratioE(Double_t Energy){
          Double_t Phi_ratioE =  F_gammaE(Energy) / F_nu_totE(Energy);
          return Phi_ratioE;
        }
	ClassDef ( nufromgamma, 1 )
};

#endif