//____________________________________________________________________________
/*!

\class   genie::flux::GSeaRealAtmoFlux

\brief   A flux driver for the Atmospheric Neutrino Flux

\author  Carla Distefano <distefano_c@lns.infn.it>
         LNS-INFN, Catania

\created September 13, 2012

\cpright  Copyright (c) 2012-2019, The KM3NeT Collaboration
          For the full text of the license see $GSEAGEN/LICENSE

*/
//____________________________________________________________________________


#include <cassert>

#include <TF1.h>
#include <TF2.h>
#include <TF3.h>


#include <TH2D.h>
#include <TH3D.h>
#include <TMath.h>
#include <TGraph.h>
#include <TSystem.h>



#include "SeaNuDrivers/GSeaAtmoFlux.h"

using namespace std;
using namespace genie;
using namespace genie::flux;
using namespace genie::constants;

//____________________________________________________________________________
GSeaAtmoFlux::GSeaAtmoFlux(GenParam * GenPar)
{

  fInitialized = 0;
  fNewNeutrino = true;

  fGenPar = GenPar;

  fAstro = new GAstro(fGenPar->SiteLatitude,fGenPar->SiteLongitude,fGenPar->UTMConvergeAngle,fGenPar->UseUTMSystem);

  if(!fGenPar->PropMode)PEarth = new GPEarth(fGenPar);


  this->Initialize();

  // set generation spectral index
  this->SetSpectralIndex(fGenPar->Alpha);

  // Configure GAtmoFlux options (common to all concrete atmospheric flux drivers)
  // set min/max energy:
  this->SetMinEnergy (fGenPar->EvMin * units::GeV);
  this->SetMaxEnergy (fGenPar->EvMax * units::GeV);

  // set min/max costheta:
  this->SetMinCosTheta (fGenPar->CtMin);
  this->SetMaxCosTheta (fGenPar->CtMax);

  // after setting og energy and angular ranges
  this->InitializeWeight();

}
//___________________________________________________________________________
GSeaAtmoFlux::~GSeaAtmoFlux()
{
  this->CleanUp();
  if(PEarth){PEarth=0; delete PEarth;}
  if(fAstro){fAstro=0; delete fAstro;}
}
//___________________________________________________________________________
void GSeaAtmoFlux::SetNewNeutrino(bool NewNeutrino){
  fNewNeutrino=NewNeutrino;
}
//___________________________________________________________________________
double GSeaAtmoFlux::MaxEnergy(void)
{
  return TMath::Min(fMaxEv, fMaxEvCut);
}
//___________________________________________________________________________
PDGCodeList & GSeaAtmoFlux::FluxParticles(void)
{
  if (fGenPar->PropMode) {
    PDGCodeList * fList = new PDGCodeList(false);
    fList->push_back(12);
    fList->push_back(-12);
    fList->push_back(14);
    fList->push_back(-14);
    fList->push_back(16);
    fList->push_back(-16);
    return *fList;
  }
  return *fPdgCList;
}//___________________________________________________________________________

bool GSeaAtmoFlux::GenerateNext(void)
{

  if(fNewNeutrino){
    LOG("GSeaAtmoFlux", pDEBUG) << "flux: new neutrino";
    this->GenerateNext_1try();
    return true;
  }
  else{
    LOG("GSeaAtmoFlux", pDEBUG) << "flux: new interaction?";
    fNewNeutrino=true;

    bool AcceptUp = fgP4.Pz()>0 && fgX4.Z()<fGenPar->Can2;
    bool AcceptDown = fgP4.Pz()<0 && fgX4.Z()>fGenPar->Can1;

    if(AcceptUp || AcceptDown){
     return true;

    }
    else {
      this->GenerateNext_1try();
      return true;
    }
   }

   return false;

}
//___________________________________________________________________________
bool GSeaAtmoFlux::GenerateNext_1try(void)
{

  // Must have run intitialization
  assert(fInitialized);


  // Reset previously generated neutrino code / 4-p / 4-x / wights
  this->ResetSelection();

  // Get a RandomGen instance
  RandomGen * rnd = RandomGen::Instance();

  // Generate the neutrino pdg code
  unsigned int nnu = fPdgCList->size();
  unsigned int inu = rnd->RndFlux().Integer(nnu);
  fgPdgC = (*fPdgCList)[inu];

  // Generate neutrino energy
  if(fSpectralIndex==1){
    fEv=exp(log(fMinEvCut)+rnd->RndFlux().Rndm()*log(fMaxEvCut/fMinEvCut));
  }
  else {
    double emin = TMath::Power(fMinEvCut,1.-fSpectralIndex);
    double emax = TMath::Power(fMaxEvCut,1.-fSpectralIndex);
    fEv = TMath::Power(emin+(emax-emin)*rnd->RndFlux().Rndm(),1./(1.-fSpectralIndex));
  }
  // Generate source neutrino direction

  double Dx=0.,Dy=0.,Dz=0.; // direction of coming point (corresponding to fTheta and fPhi)

  if(fGenPar->GenMJD)fMJD = fGenPar->MJDStart+rnd->RndFlux().Rndm()*(fGenPar->MJDStop-fGenPar->MJDStart);

  if(fGenPar->GenMode.compare("DIFFUSE")==0){
    Dz = fMinCtCut+(fMaxCtCut-fMinCtCut)*rnd->RndFlux().Rndm();
    fTheta = acos(Dz);
    fPhi = 2.*kPi* rnd->RndFlux().Rndm();
    Dx = TMath::Sin(fTheta) * TMath::Sin(fPhi);
    Dy = TMath::Sin(fTheta) * TMath::Cos(fPhi);
  }
  else if(fGenPar->GenMode.compare("POINT")==0){

    if(fGenPar->GenMJD){
      fAstro->CalcSrcDirMJD(fMJD,fGenPar->RightAscension,fGenPar->Declination,&Dx,&Dy,&Dz);
    }
    else{
      fLST = 2.*kPi* rnd->RndFlux().Rndm();
      fAstro->CalcSrcDirLST(fLST,fGenPar->RightAscension,fGenPar->Declination,&Dx,&Dy,&Dz);
    }
    fTheta = acos(Dz);
    fPhi = atan2(Dy,Dx);

    if(fGenPar->SourceRadius>0.){  // extended source
      double cosTh=Dz;
      double sinTh=sin(fTheta);
      double cosPh=cos(fPhi);
      double sinPh=sin(fPhi);
      double cosThExt=cos(fGenPar->SourceRadius)+(1.-cos(fGenPar->SourceRadius))*rnd->RndFlux().Rndm();
      double sinThExt=sqrt(1-cosThExt*cosThExt);
      double PhExt=2.*kPi* rnd->RndFlux().Rndm();
      double cosPhExt=cos(PhExt);
      double sinPhExt=sin(PhExt);

      Dx = cosTh*cosPh*cosPhExt*sinThExt-sinPh*sinPhExt*sinThExt+cosPh*sinTh*cosThExt;
      Dy = cosTh*sinPh*cosPhExt*sinThExt+cosPh*sinPhExt*sinThExt+sinTh*sinPh*cosThExt;
      Dz = -sinTh*cosPhExt*sinThExt+cosTh*cosThExt;
      fTheta = acos(Dz);
      fPhi = atan2(Dy,Dx);
    }
  }
  else{
    LOG("GSeaAtmoFlux", pFATAL) << " invalid generation mode; exiting... ";
    exit(1);
  }

  // Compute the neutrino momentum and set the correspondig 4-vector
  double pz = -1.* fEv * Dz;
  double py = -1.* fEv * Dy;
  double px = -1.* fEv * Dx;

  fgP4.SetPxPyPzE(px,py,pz,fEv);

  double x = 0.;
  double y = 0.;
  double z = 0.;

  if (fRt<=0.) {

//------------------------------------------------------------------------------------------
//new method -> shooting neutrinos using projected surface (Agen is variable)
//------------------------------------------------------------------------------------------

    //RT was saved as -fRTPlus, so we have to convert it back to positive
    double height = (fGenPar->Can2-fGenPar->Can1+fGenPar->HBedRock) - fRt;
    double radius = fGenPar->Can3 - fRt;

    double AreaTop  = kPi * TMath::Power(radius,2) * TMath::Abs(Dz) ;
    double AreaSide = 2 * height * radius * TMath::Sqrt( 1 - Dz*Dz );
    fAgen = AreaTop + AreaSide;

    if ( rnd->RndFlux().Rndm() < AreaTop / fAgen ) {
      double r   = radius * TMath::Sqrt( rnd->RndFlux().Rndm() );
      double phi = 2 * kPi * rnd->RndFlux().Rndm();
      x = r*TMath::Sin(phi);
      y = r*TMath::Cos(phi);
      z = Dz < 0 ? -height/2. : height/2.;
    }
    else {
      double zaux = height * (rnd->RndFlux().Rndm()-0.5);
      double yaux = 2 * radius * (rnd->RndFlux().Rndm()-0.5);
      double xaux = -TMath::Sqrt( radius*radius - yaux*yaux );
      TVector3 v(xaux,yaux,zaux);
      double phidir = TMath::ATan2( -Dy, -Dx );
      v.RotateZ( phidir );
      x = v.X();
      y = v.Y();
      z = v.Z();
    }

    if( fRl>0.0 ){
      z += fRl * Dz;
      y += fRl * Dy;
      x += fRl * Dx;
    }

  }
  else {

//------------------------------------------------------------------------------------------
//old method -> shooting neutrinos from a circular sphere (Agen is constant)
//------------------------------------------------------------------------------------------

    // Shift the neutrino position onto the flux generation surface.
    // The position is computed at the surface of a sphere with R=fRl
    // at the topocentric horizontal (THZ) coordinate system.
    if( fRl>0.0 ){
      z += fRl * Dz;
      y += fRl * Dy;
      x += fRl * Dx;
    }

    // If the position is left as is, then all generated neutrinos
    // would point towards the origin.
    // Displace the position randomly on the surface that is
    // perpendicular to the selected point P(xo,yo,zo) on the sphere
    if( fRt>0.0 ){
      fAgen = kPi*TMath::Power(fRt,2);
      TVector3 vec(x,y,z);         // vector towards selected point
      TVector3 dvec1 = vec.Orthogonal(); // orthogonal vector
      TVector3 dvec2 = dvec1;        // second orthogonal vector
      dvec2.Rotate(-kPi/2.0,vec);        // rotate second vector by 90deg, now forming a new orthogonal cartesian coordinate system
      double psi = 2.*kPi* rnd->RndFlux().Rndm(); // rndm angle [0,2pi]
      double random = rnd->RndFlux().Rndm();  // rndm number  [0,1]
      dvec1.SetMag(TMath::Sqrt(random)*fRt*TMath::Cos(psi));
      dvec2.SetMag(TMath::Sqrt(random)*fRt*TMath::Sin(psi));
      x += dvec1.X() + dvec2.X();
      y += dvec1.Y() + dvec2.Y();
      z += dvec1.Z() + dvec2.Z();
    }

  }

  // Set the neutrino position 4-vector with values
  // calculated at previous steps.
  fgX4.SetXYZT(x+fX0,y+fY0,z+fZ0,0.);

  // Apply user-defined rotation from THZ -> user-defined topocentric coordinate system
  if(!fRotTHz2User.IsIdentity()){
    TVector3 tx3 = fgX4.Vect();
    TVector3 tp3 = fgP4.Vect();
    tx3 = fRotTHz2User * tx3;
    tp3 = fRotTHz2User * tp3;
    fgX4.SetVect(tx3);
    fgP4.SetVect(tp3);
  }

  // Increment flux neutrino counter used for sample normalization purposes.
  fNNeutrinos++;

  fNInt=1;  // it is the first interaction

  fgP4I=fgP4;
  fgX4I=fgX4;
  fgPdgCI=fgPdgC;

  // Report and exit
  LOG("GSeaAtmoFlux", pDEBUG)
       << "Generated neutrino: "
       << "\n pdg-code: " << fgPdgC
       << "\n p4: " << utils::print::P4AsShortString(&fgP4)
       << "\n x4: " << utils::print::X4AsString(&fgX4);

  return true;
}
//___________________________________________________________________________
void GSeaAtmoFlux::SetNeutrino(TLorentzVector * VecEne, TLorentzVector * VecPos, int pdg, bool NewNeutrino){

// fEv, fTheta and fPhi are not set. They are used to comupute the weights and must be the generated ones.

  fgPdgC = pdg;

  fgX4.SetXYZT(VecPos->X(),VecPos->Y(),VecPos->Z(),VecPos->T());
  fgP4.SetPxPyPzE(VecEne->Px(),VecEne->Py(),VecEne->Pz(),VecEne->Energy());


  fNewNeutrino = NewNeutrino;
  fNInt++;                      // increment the number of interactions

  // Report and exit
  LOG("GSeaAtmoFlux", pDEBUG)
       << "Set neutrino: "
       << "\n pdg-code: " << fgPdgC
       << "\n p4: " << utils::print::P4AsShortString(&fgP4)
       << "\n x4: " << utils::print::X4AsString(&fgX4);


}
//___________________________________________________________________________
void GSeaAtmoFlux::ResetNFluxNeutrinos(void)
{
  fNNeutrinos=0;
  return;
}
//___________________________________________________________________________
long int GSeaAtmoFlux::NFluxNeutrinos(void) const
{
  return fNNeutrinos;
}
//___________________________________________________________________________
void GSeaAtmoFlux::SetMinEnergy(double emin)
{
  emin = TMath::Max(0., emin);
  fMinEvCut = emin;
}
//___________________________________________________________________________
void GSeaAtmoFlux::SetMaxEnergy(double emax)
{
  emax = TMath::Max(0., emax);
  fMaxEvCut = emax;
}
//___________________________________________________________________________
void GSeaAtmoFlux::SetMinCosTheta(double ctmin)
{
  ctmin = TMath::Max(-1.,ctmin);
  fMinCtCut = ctmin;
}
//___________________________________________________________________________
void GSeaAtmoFlux::SetMaxCosTheta(double ctmax)
{
  ctmax = TMath::Min(1.,ctmax);
  fMaxCtCut = ctmax;
}
//___________________________________________________________________________
void GSeaAtmoFlux::InitGenBin(double Emin, double Emax, double RL, double RT, double X0, double Y0, double Z0)
{
  this->SetMinEnergy(Emin);
  this->SetMaxEnergy(Emax);
  this->ResetNFluxNeutrinos();
  this->SetRadii(RL,RT);

  fX0=X0;
  fY0=Y0;
  fZ0=Z0;

}

//___________________________________________________________________________
void GSeaAtmoFlux::Clear(Option_t * opt)
{
// Dummy clear method needed to conform to GFluxI interface
//
  LOG("GSeaAtmoFlux", pERROR) << "No clear method implemented for option:"<< opt;
}
//___________________________________________________________________________
void GSeaAtmoFlux::GenerateWeighted(bool gen_weighted)
{
  fGenWeighted = gen_weighted;
}
//___________________________________________________________________________
void GSeaAtmoFlux:: SetSpectralIndex(double index)
{
  fSpectralIndex = index;
  LOG("GSeaAtmoFlux", pNOTICE) << "Using Spectral Index = " << -index;
}
//___________________________________________________________________________
void GSeaAtmoFlux::SetUserCoordSystem(TRotation & rotation)
{
  fRotTHz2User = rotation;
}
//___________________________________________________________________________
void GSeaAtmoFlux::Initialize(void)
{
  LOG("GSeaAtmoFlux", pNOTICE) << "Initializing atmospheric flux driver";

  bool allow_dup = false;
  fPdgCList = new PDGCodeList(allow_dup);

  for(unsigned i=0; i<fGenPar->NuList.size(); i++){

    vector<FluxFunc1D> Vec1F;
    fFlux1FRange.insert(map<int, vector<FluxFunc1D> >::value_type(fGenPar->NuList[i],Vec1F));

    vector<FluxFunc2D> Vec2F;
    fFlux2FRange.insert(map<int, vector<FluxFunc2D> >::value_type(fGenPar->NuList[i],Vec2F));

    vector<FluxFunc3D> Vec3F;
    fFlux3FRange.insert(map<int, vector<FluxFunc3D> >::value_type(fGenPar->NuList[i],Vec3F));

    vector<FluxHisto1D> Vec1D;
    fFlux1DRange.insert(map<int, vector<FluxHisto1D> >::value_type(fGenPar->NuList[i],Vec1D));

    vector<FluxHisto2D> Vec2D;
    fFlux2DRange.insert(map<int, vector<FluxHisto2D> >::value_type(fGenPar->NuList[i],Vec2D));

    vector<FluxHisto3D> Vec3D;
    fFlux3DRange.insert(map<int, vector<FluxHisto3D> >::value_type(fGenPar->NuList[i],Vec3D));

    fPdgCList->push_back(fGenPar->NuList[i]);
  }

  // Default radii
  fRl = 0.;
  fRt = 0.;

  // Default detector coord system: Topocentric Horizontal Coordinate system
  fRotTHz2User.SetToIdentity();

  // Reset `current' selected flux neutrino
  this->ResetSelection();

  // Reset number of neutrinos thrown so far
  fNNeutrinos = 0;

  // Reset the generation weight
  fGlobalGenWeight = -1.;

  // Reset the LST and MJD
  fMJD=0.;
  fLST=0.;

  // Done!
  fInitialized = 1;
  LOG("GSeaAtmoFlux", pNOTICE) << "Initializing atmospheric flux driver... done";

}
//___________________________________________________________________________
double GSeaAtmoFlux::InitializeWeight()
{
  LOG("GSeaAtmoFlux", pNOTICE) << "Initializing generation weight";

  double IE=1.;

  if(fSpectralIndex==1){
    IE=log(fMaxEvCut/fMinEvCut);
  }
  else {
    IE=(pow(fMaxEvCut,(1.-fSpectralIndex))-pow(fMinEvCut,(1.-fSpectralIndex)))/(1.-fSpectralIndex);
  }

  double ITheta=1.;
  if(fGenPar->GenMode.compare("DIFFUSE")==0)ITheta=2.*kPi*(fMaxCtCut-fMinCtCut);
//  if(fGenPar->GenMode.compare("POINT")==0 && fGenPar->SourceRadius>0)ITheta=kPi*TMath::Power(fGenPar->SourceRadius,2); --> 16 sept 2019 no solid angle in extended sources

  fGlobalGenWeight=IE*ITheta*fGenPar->TGen*fGenPar->NuList.size();

#ifdef _HETAU_ENABLED__
  //branching ratios of the tau decays
  if      (fGenPar->TrackEvts==2) fGlobalGenWeight *= 0.174;    //muon
  else if (fGenPar->TrackEvts==3) fGlobalGenWeight *= 1.-0.174; //others
#endif

  return fGlobalGenWeight;
}
//___________________________________________________________________________
void GSeaAtmoFlux::ComputeWeights(void){

  double X0=fgX4.X();
  double Y0=fgX4.Y();
  double Z0=fgX4.Z();

  this->ComputeWeights(X0, Y0, Z0);

  return;
}
//___________________________________________________________________________
void GSeaAtmoFlux::ComputeWeights(double X0, double Y0, double Z0){

  fWeight = this->GetFlux(fgPdgC,fEv,fTheta,fPhi);

  fGenWeight = fGlobalGenWeight*fAgen*TMath::Power(fEv,fSpectralIndex);

  if(!fGenPar->PropMode)fGenWeight *= PEarth->EarthAbsorption(X0,Y0,Z0,fgP4.Px()/fEv,fgP4.Py()/fEv,fgP4.Pz()/fEv,fEv, fgPdgC);

  return;
}
//___________________________________________________________________________
void GSeaAtmoFlux::ComputeWeights(GSeaEvent * SeaEvt){


  this->ComputeWeights(SeaEvt->Vertex.X(), SeaEvt->Vertex.Y(), SeaEvt->Vertex.Z());

  SeaEvt->NInter=fNInt;
  SeaEvt->ColumnDepth=this->GetColumnDepth();
  SeaEvt->PEarth=this->GetPEarth();
  SeaEvt->XSecMean=this->GetSigmaMean();
  SeaEvt->Agen=this->GetAgen();

  SeaEvt->GenWeight=this->GenWeight()*SeaEvt->PScale*fNInt;
  SeaEvt->EvtWeight=SeaEvt->GenWeight*fWeight;

  SeaEvt->LST=this->GetLST();
  SeaEvt->MJD=this->GetMJD();

  this->GetEvtTime(SeaEvt->EvtTime,SeaEvt->EvtTime+1);

  return;
}
//___________________________________________________________________________
double GSeaAtmoFlux::GenWeight(void)
{
  if(fGlobalGenWeight<0){
    LOG ("GSeaAtmoFlux", pWARN)
     << "Computation of generation weight not initialized";
    fGenWeight=0.;
  }

  return fGenWeight;
}
//___________________________________________________________________________
void GSeaAtmoFlux::ResetSelection(void)
{
  // initializing running neutrino pdg-code, 4-position, 4-momentum, weights

  fgPdgC = 0;

  fgP4.SetPxPyPzE (0.,0.,0.,0.);
  fgX4.SetXYZT    (0.,0.,0.,0.);

  fGenWeight = 0.;
  fWeight = 0.;
  fAgen = 0.;
}
//___________________________________________________________________________
void GSeaAtmoFlux::CleanUp(void)
{
  LOG("GSeaAtmoFlux", pNOTICE) << "Cleaning up...";


  map<int,vector<FluxFunc1D> >::iterator MapEntry1F;
  vector<FluxFunc1D> Vec1F;
  for ( MapEntry1F=fFlux1FRange.begin() ; MapEntry1F != fFlux1FRange.end(); MapEntry1F++ ){
     Vec1F=MapEntry1F->second;
     for(unsigned i=0; i< Vec1F.size(); i++){
       TF1 * flux_histogram = MapEntry1F->second[i].Flux1D;
       if(flux_histogram){
       delete flux_histogram;
       flux_histogram = 0;
       }
     }
   }
  fFlux1FRange.clear();

  map<int,vector<FluxFunc2D> >::iterator MapEntry2F;
  vector<FluxFunc2D> Vec2F;
  for ( MapEntry2F=fFlux2FRange.begin() ; MapEntry2F != fFlux2FRange.end(); MapEntry2F++ ){
     Vec2F=MapEntry2F->second;
     for(unsigned i=0; i< Vec2F.size(); i++){
       TF1 * flux_histogram = MapEntry2F->second[i].Flux2D;
       if(flux_histogram){
       delete flux_histogram;
       flux_histogram = 0;
       }
     }
   }
  fFlux2FRange.clear();

  map<int,vector<FluxFunc3D> >::iterator MapEntry3F;
  vector<FluxFunc3D> Vec3F;
  for ( MapEntry3F=fFlux3FRange.begin() ; MapEntry3F != fFlux3FRange.end(); MapEntry3F++ ){
     Vec3F=MapEntry3F->second;
     for(unsigned i=0; i< Vec3F.size(); i++){
       TF1 * flux_histogram = MapEntry3F->second[i].Flux3D;
       if(flux_histogram){
       delete flux_histogram;
       flux_histogram = 0;
       }
     }
   }
  fFlux3FRange.clear();

  map<int,vector<FluxHisto1D> >::iterator MapEntry1D;
  vector<FluxHisto1D> Vec1D;
  for ( MapEntry1D=fFlux1DRange.begin() ; MapEntry1D != fFlux1DRange.end(); MapEntry1D++ ){
     Vec1D=MapEntry1D->second;
     for(unsigned i=0; i< Vec1D.size(); i++){
       TH1D * flux_histogram = MapEntry1D->second[i].Flux1D;
       if(flux_histogram){
       delete flux_histogram;
       flux_histogram = 0;
       }
     }
   }
  fFlux1DRange.clear();

  map<int,vector<FluxHisto2D> >::iterator MapEntry2D;
  vector<FluxHisto2D> Vec2D;
  for ( MapEntry2D=fFlux2DRange.begin() ; MapEntry2D != fFlux2DRange.end(); MapEntry2D++ ){
     Vec2D=MapEntry2D->second;
     for(unsigned i=0; i< Vec2D.size(); i++){
       TH2D * flux_histogram = MapEntry2D->second[i].Flux2D;
       if(flux_histogram){
       delete flux_histogram;
       flux_histogram = 0;
       }
     }
   }
  fFlux2DRange.clear();

  map<int,vector<FluxHisto3D> >::iterator MapEntry3;
  vector<FluxHisto3D> Vec3;
  for ( MapEntry3=fFlux3DRange.begin() ; MapEntry3 != fFlux3DRange.end(); MapEntry3++ ){
     Vec3=MapEntry3->second;
     for(unsigned i=0; i< Vec3.size(); i++){
       TH3D * flux_histogram = MapEntry3->second[i].Flux3D;
       if(flux_histogram){
       delete flux_histogram;
       flux_histogram = 0;
       }
     }
   }
  fFlux3DRange.clear();


  if (fPdgCList ) delete fPdgCList;
}
//___________________________________________________________________________
void GSeaAtmoFlux::SetRadii(double Rlongitudinal, double Rtransverse)
{
  LOG ("GSeaAtmoFlux", pNOTICE) << "Setting R[longitudinal] = " << Rlongitudinal;
  LOG ("GSeaAtmoFlux", pNOTICE) << "Setting R[transverse]   = " << Rtransverse;

  fRl = Rlongitudinal;
  fRt = Rtransverse;
}
//___________________________________________________________________________
double GSeaAtmoFlux::GetFlux(int flavour, double energy, double Theta, double phi)
{
  double flux=0.;
  double Emin,Emax,cosThMin=-1,cosThMax=1,PhiMin=0., PhiMax=2.*kPi;
  double cosTheta = TMath::Cos(Theta);

  if(fGenPar->UseUTMSystem) phi -= fGenPar->UTMConvergeAngle + kPi/2.;


  map<int,vector<FluxFunc1D> >::iterator MapEntry1F = fFlux1FRange.find(flavour);
  for(unsigned i=0; i<MapEntry1F->second.size(); i++){
    Emin=MapEntry1F->second[i].Range[0];
    Emax=MapEntry1F->second[i].Range[1];
    if(energy>=Emin && energy<=Emax)
      flux += MapEntry1F->second[i].Flux1D->Eval(energy);
  }

  map<int,vector<FluxFunc2D> >::iterator MapEntry2F = fFlux2FRange.find(flavour);
  for(unsigned i=0; i<MapEntry2F->second.size(); i++){
    Emin=MapEntry2F->second[i].Range[0];
    Emax=MapEntry2F->second[i].Range[1];
    cosThMin=MapEntry2F->second[i].Range[2];
    cosThMax=MapEntry2F->second[i].Range[3];
    if(energy>=Emin && energy<=Emax && cosTheta>=cosThMin && cosTheta<=cosThMax)
      flux += MapEntry2F->second[i].Flux2D->Eval(energy,cosTheta);
  }

  map<int,vector<FluxFunc3D> >::iterator MapEntry3F = fFlux3FRange.find(flavour);
  for(unsigned i=0; i<MapEntry3F->second.size(); i++){
    Emin=MapEntry3F->second[i].Range[0];
    Emax=MapEntry3F->second[i].Range[1];
    cosThMin=MapEntry3F->second[i].Range[2];
    cosThMax=MapEntry3F->second[i].Range[3];
    PhiMin=MapEntry3F->second[i].Range[4];
    PhiMax=MapEntry3F->second[i].Range[5];
    if(energy>=Emin && energy<=Emax && cosTheta>=cosThMin && cosTheta<=cosThMax && phi>=PhiMin && phi<=PhiMax)
      flux += MapEntry3F->second[i].Flux3D->Eval(energy,cosTheta,phi);
  }

  map<int,vector<FluxHisto1D> >::iterator MapEntry1D = fFlux1DRange.find(flavour);
  for(unsigned i=0; i<MapEntry1D->second.size(); i++){
    Emin=MapEntry1D->second[i].Range[0];
    Emax=MapEntry1D->second[i].Range[1];
    if(energy>=Emin && energy<=Emax)
      flux += MapEntry1D->second[i].Flux1D->Interpolate(energy);
  }


  map<int,vector<FluxHisto2D> >::iterator MapEntry2D = fFlux2DRange.find(flavour);
  for(unsigned i=0; i<MapEntry2D->second.size(); i++){
    Emin=MapEntry2D->second[i].Range[0];
    Emax=MapEntry2D->second[i].Range[1];
    cosThMin=MapEntry2D->second[i].Range[2];
    cosThMax=MapEntry2D->second[i].Range[3];
    if(energy>=Emin && energy<=Emax && cosTheta>=cosThMin && cosTheta<=cosThMax)
      flux += MapEntry2D->second[i].Flux2D->Interpolate(energy,cosTheta);
  }

  map<int,vector<FluxHisto3D> >::iterator MapEntry3D = fFlux3DRange.find(flavour);
  for(unsigned i=0; i<MapEntry3D->second.size(); i++){
    Emin=MapEntry3D->second[i].Range[0];
    Emax=MapEntry3D->second[i].Range[1];
    cosThMin=MapEntry3D->second[i].Range[2];
    cosThMax=MapEntry3D->second[i].Range[3];
    PhiMin=MapEntry3D->second[i].Range[4];
    PhiMax=MapEntry3D->second[i].Range[5];

    if(energy>=Emin && energy<=Emax && cosTheta>=cosThMin && cosTheta<=cosThMax && phi>=PhiMin && phi<=PhiMax){

     // we take into account that int TH3 interpolation,
     // the given values (x,y,z) must be between first bin center and last bin center for each coordinate

      if(energy<MapEntry3D->second[i].Flux3D->GetXaxis()->GetBinCenter(1))
        energy=MapEntry3D->second[i].Flux3D->GetXaxis()->GetBinCenter(1);
      if(energy>MapEntry3D->second[i].Flux3D->GetXaxis()->GetBinCenter(MapEntry3D->second[i].Flux3D->GetXaxis()->GetNbins()))
        energy=(MapEntry3D->second[i].Flux3D->GetXaxis()->GetBinCenter(MapEntry3D->second[i].Flux3D->GetXaxis()->GetNbins())+MapEntry3D->second[i].Flux3D->GetXaxis()->GetBinLowEdge(MapEntry3D->second[i].Flux3D->GetXaxis()->GetNbins()))/2.;

      if(cosTheta<MapEntry3D->second[i].Flux3D->GetYaxis()->GetBinCenter(1))
        cosTheta=MapEntry3D->second[i].Flux3D->GetYaxis()->GetBinCenter(1);
      if(cosTheta>MapEntry3D->second[i].Flux3D->GetYaxis()->GetBinCenter(MapEntry3D->second[i].Flux3D->GetYaxis()->GetNbins()))
        cosTheta=(MapEntry3D->second[i].Flux3D->GetYaxis()->GetBinCenter(MapEntry3D->second[i].Flux3D->GetYaxis()->GetNbins())+MapEntry3D->second[i].Flux3D->GetYaxis()->GetBinLowEdge(MapEntry3D->second[i].Flux3D->GetYaxis()->GetNbins()))/2.;

      if(phi<MapEntry3D->second[i].Flux3D->GetZaxis()->GetBinCenter(1))
        phi=MapEntry3D->second[i].Flux3D->GetZaxis()->GetBinCenter(1);
      if(phi>MapEntry3D->second[i].Flux3D->GetZaxis()->GetBinCenter(MapEntry3D->second[i].Flux3D->GetZaxis()->GetNbins()))
        phi=(MapEntry3D->second[i].Flux3D->GetZaxis()->GetBinCenter(MapEntry3D->second[i].Flux3D->GetZaxis()->GetNbins())+MapEntry3D->second[i].Flux3D->GetZaxis()->GetBinLowEdge(MapEntry3D->second[i].Flux3D->GetZaxis()->GetNbins()))/2.;

      flux += MapEntry3D->second[i].Flux3D->Interpolate(energy,cosTheta,phi);

    }
  }

  return flux;
}
//___________________________________________________________________________
double GSeaAtmoFlux::GetPEarth(void) {


  if(!fGenPar->PropMode)return PEarth->GetPEarth();

  return 0.;
}
//___________________________________________________________________________
double GSeaAtmoFlux::GetColumnDepth(void) {

  if(!fGenPar->PropMode)return PEarth->GetColumnDepth();
  else return 0.;

}
//___________________________________________________________________________
double GSeaAtmoFlux::GetSigmaMean(void) {
  if(!fGenPar->PropMode)return PEarth->GetSigmaMean();
  else return 0.;
}
//___________________________________________________________________________
double GSeaAtmoFlux::GetMJD(void) {
  if(fGenPar->GenMJD)return fMJD;
  else return 0.;
}
//___________________________________________________________________________
double GSeaAtmoFlux::GetLST(void) {
  if(fGenPar->GenMode.compare("POINT")==0)return fLST;
  else return 0.;
}
//___________________________________________________________________________
void GSeaAtmoFlux::GetEvtTime(unsigned int * TimeStampSec, unsigned int * TimeStampTick) {

  if(fGenPar->GenMJD)fAstro->CalcTimeStamp(fMJD, TimeStampSec, TimeStampTick);

  return;
}
//___________________________________________________________________________
TH2D * GSeaAtmoFlux::CreateFluxHisto2D(string name, string title)
{
  LOG("GSeaAtmoFlux", pNOTICE) << "Instantiating histogram: [" << name << "]";

  double EnergyBins[fNumLogEvBins+1];
  double LogEmin=TMath::Log10(fMinEv);

  for(unsigned int i=0; i<=fNumLogEvBins; i++)
    EnergyBins[i]=TMath::Power(10.,LogEmin+i*fDeltaLogEv);

  TH2D * h2 = new TH2D(name.c_str(), title.c_str(),fNumLogEvBins, EnergyBins, fNumCosThetaBins, fCosThetaMin,fCosThetaMax);

  return h2;
}

//___________________________________________________________________________
TH3D * GSeaAtmoFlux::CreateFluxHisto3D(string name, string title)
{
  LOG("GSeaAtmoFlux", pNOTICE) << "Instantiating histogram: [" << name << "]";

  double EnergyBins[fNumLogEvBins+1];
  double LogEmin=TMath::Log10(fMinEv);

  for(unsigned int i=0; i<=fNumLogEvBins; i++)
    EnergyBins[i]=TMath::Power(10.,LogEmin+((double)i)*fDeltaLogEv);

  double CosThetaBins[fNumCosThetaBins+1];
  for(unsigned int i=0; i<=fNumCosThetaBins; i++)
    CosThetaBins[i]=fCosThetaMin+((double)i)*(fCosThetaMax-fCosThetaMin)/(double)fNumCosThetaBins;

  double PhiBins[fNumPhiBins+1];
  for(unsigned int i=0; i<=fNumPhiBins; i++)
    PhiBins[i]=fPhiMin+((double)i)*(fPhiMax-fPhiMin)/(double)fNumPhiBins;

  TH3D * h3 = new TH3D(name.c_str(), title.c_str(),fNumLogEvBins, EnergyBins, fNumCosThetaBins, CosThetaBins, fNumPhiBins, PhiBins);

  return h3;
}

//___________________________________________________________________________
void GSeaAtmoFlux::FillFuncFlux(unsigned int iFlux, string FuncName){

  if(fGenPar->FluxFiles[iFlux].FluxSimul == "FUNC1"){

    TF1 *  input_func = new TF1(FuncName.c_str(), fGenPar->FluxFiles[iFlux].FileNames[0].c_str(), fMinEv, fMaxEv);
    FluxFunc1D flux;
    flux.Flux1D=input_func;
    flux.Range[0]=fMinEv;
    flux.Range[1]=fMaxEv;

    map<int,vector<FluxFunc1D> >::iterator MapEntry = fFlux1FRange.find(fGenPar->FluxFiles[iFlux].NuType);
    MapEntry->second.push_back(flux);

  }
  else if(fGenPar->FluxFiles[iFlux].FluxSimul == "FUNC2"){

    TF2 *  input_func = new TF2(FuncName.c_str(), fGenPar->FluxFiles[iFlux].FileNames[0].c_str(), fMinEv, fMaxEv,fCosThetaMin,fCosThetaMax);
    FluxFunc2D flux;
    flux.Flux2D=input_func;
    flux.Range[0]=fMinEv;
    flux.Range[1]=fMaxEv;
    flux.Range[2]=fCosThetaMin;
    flux.Range[3]=fCosThetaMax;

    map<int,vector<FluxFunc2D> >::iterator MapEntry = fFlux2FRange.find(fGenPar->FluxFiles[iFlux].NuType);
    MapEntry->second.push_back(flux);

  }
  else if(fGenPar->FluxFiles[iFlux].FluxSimul == "FUNC3"){

    TF3 *  input_func = new TF3(FuncName.c_str(), fGenPar->FluxFiles[iFlux].FileNames[0].c_str(), fMinEv, fMaxEv,fCosThetaMin,fCosThetaMax,fPhiMin,fPhiMax);
    FluxFunc3D flux;
    flux.Flux3D=input_func;
    flux.Range[0]=fMinEv;
    flux.Range[1]=fMaxEv;
    flux.Range[2]=fCosThetaMin;
    flux.Range[3]=fCosThetaMax;
    flux.Range[4]=fPhiMin;
    flux.Range[5]=fPhiMax;

    map<int,vector<FluxFunc3D> >::iterator MapEntry = fFlux3FRange.find(fGenPar->FluxFiles[iFlux].NuType);
    MapEntry->second.push_back(flux);

  }

  return;
}