#ifndef NUFROMGAMMAHH
#define NUFROMGAMMAHH


#include <iostream>
#include <cstdio>
#include <fstream>
#include <sstream>
#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 <fstream>
#include <sstream>
#include <vector>
#include <string>
#include "TSpline.h"
using namespace std;

/***************************************************************************************************
	
 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;

       //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() {	}
       nufromgamma(string gammafluxeq) {
					//cout << "Debug: nufromgamma constructor called with: " << gammafluxeq << endl;
				 f_gamma = TF1("fgammaflux",gammafluxeq.c_str(),100.,10000000.);
				 f_nu= to_spline("numu");
				 f_antinu= to_spline("antinumu");
				 f_nu_tot= to_spline("tot");
       }

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

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

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

       //    not the right position!!!!!!!!!!!

	TSpline3 to_spline(string type )
	{
           const double logemin = 2;
           const double logemax = 8;
           const double nbinloge = 1e4;

           vector<double> xs;
	   vector<double> ys;

	   double E;

	   for(double loge=logemin;loge<logemax;loge+=(logemax-logemin)/nbinloge){
	    E=pow(10,loge);
	    xs.push_back(loge);
	    if (type=="tot")   ys.push_back( F_nu_tot(E));
	    if (type=="numu")   ys.push_back( F_nu(E));
	    if (type=="antinumu")   ys.push_back( F_antinu(E));

    	   }

	   xs.push_back( logemax );
	   ys.push_back( 0 );

	   TSpline3 spline(type.c_str(), &(xs[0]), &(ys[0]), xs.size() , "e1b1"  );
	   spline.SetLineColor(4);
	   spline.SetLineWidth(2);
           return spline;
	}


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

        // 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;
        }

        // 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;
        }

        //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);
        }

  
  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