//____________________________________________________________________________
/*!

\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

*/
//____________________________________________________________________________


#ifndef _GATMO_FLUX_H_
#define _GATMO_FLUX_H_

#include <string>
#include <map>
#include <vector>
#include <cstdio>
#include <iostream>
#include <fstream>
#include <cmath>

#include <TLorentzVector.h>
#include <TRotation.h>

#include "SeaNuDrivers/GAstro.h"
#include "SeaNuDrivers/GPEarth.h"
#include "SeaEvent/GSeaEvent.h"
#include "SeaEvent/GBinParam.h"

#ifdef _GENIEVERSIONGTEQ3__
#include "Framework/Conventions/Constants.h"
#include "Framework/EventGen/GFluxI.h"
#ifndef __CINT__
#include "Framework/Messenger/Messenger.h"
#endif
#include "Framework/Numerical/RandomGen.h"
#include "Framework/ParticleData/PDGCodeList.h"
#include "Framework/Utils/PrintUtils.h"
#else
#include "Conventions/Constants.h"
#include "EVGDrivers/GFluxI.h"
#ifndef __CINT__
#include "Messenger/Messenger.h"
#endif
#include "Numerical/RandomGen.h"
#include "PDG/PDGCodeList.h"
#include "Utils/PrintUtils.h"
#endif



class TF1;
class TF2;
class TF3;
class TH1D;
class TH2D;
class TH3D;



using namespace std;
using namespace genie;


namespace genie {
namespace flux  {

struct FluxFunc1D {TF1* Flux1D; double Range[2];};
struct FluxFunc2D {TF2* Flux2D; double Range[4];};
struct FluxFunc3D {TF3* Flux3D; double Range[6];};

struct FluxHisto1D {TH1D* Flux1D; double Range[2];};
struct FluxHisto2D {TH2D* Flux2D; double Range[4];};
struct FluxHisto3D {TH3D* Flux3D; double Range[6];};


class GSeaAtmoFlux: public GFluxI {

public :
  virtual ~GSeaAtmoFlux();

  // methods implementing the GENIE GFluxI interface
//  virtual const PDGCodeList &    FluxParticles (void) { return *fPdgCList; }
  virtual PDGCodeList &          FluxParticles (void);  // Alfonso Garcia
  virtual double                 MaxEnergy     (void);
  virtual bool                   GenerateNext  (void);
  virtual int                    PdgCode       (void) { return fgPdgC;     }
  virtual double                 Weight        (void) { return fWeight;    }
  virtual const TLorentzVector & Momentum      (void) { return fgP4;       }
  virtual const TLorentzVector & Position      (void) { return fgX4;       }
  virtual bool                   End           (void) { return false;      }
  virtual long int               Index         (void) { return -1;         }
  virtual void                   Clear            (Option_t * opt);
  virtual void                   GenerateWeighted (bool gen_weighted);


  // get neutrino energy/direction of generated events
  double Enu        (void) { return fgP4.Energy(); }
  double Energy     (void) { return fgP4.Energy(); }
  double CosTheta   (void) { return -fgP4.Pz()/fgP4.Energy(); }
  double CosZenith  (void) { return -fgP4.Pz()/fgP4.Energy(); }


  // set neutrino to simulate a new interaction (PropMode)
  void   SetNeutrino(TLorentzVector * VecEne, TLorentzVector * VecPos, int pdg, bool NewNeutrino);
  void   SetNewNeutrino     (bool NewNeutrino);


  // methods specific to the atmospheric flux drivers
  long int NFluxNeutrinos     (void) const; ///< Number of flux nu's generated. Not the same as the number of nu's thrown towards the geometry (if there are cuts).
  void     ResetNFluxNeutrinos(void);
  void     SetMinEnergy       (double emin);
  void     SetMaxEnergy       (double emax);
  void     SetMinCosTheta     (double ctmin);
  void     SetMaxCosTheta     (double ctmax);
  void     SetSpectralIndex   (double index);
  void     SetRadii           (double Rlongitudinal, double Rtransverse);
  void     InitGenBin         (double Emin, double Emax, double RL, double RT, double X0, double Y0, double Z0);
  void     SetUserCoordSystem (TRotation & rotation); ///< Rotation: Topocentric Horizontal -> User-defined Topocentric Coord System.
  double   InitializeWeight   ();
  double   GenWeight          (void); //return the generation weight
  double   GetGlobalGenWeight (void) {return fGlobalGenWeight;}
  double   GetAgen            (void) {return fAgen;}
  double   GetPEarth          (void);
  double   GetColumnDepth     (void);
  double   GetSigmaMean       (void);
  double   GetLST             (void);
  double   GetMJD             (void);
  void     GetEvtTime         (unsigned int * TimeStampSec, unsigned int * TimeStampTick);

  void     ComputeWeights     (void);
  void     ComputeWeights     (double X0, double Y0, double Z0);
  void     ComputeWeights     (GSeaEvent * SeaEvt);

  double   GetFlux            (int flavour, double LogEne, double Theta, double Phi);


protected:

  // abstract class, ctor hidden
  GSeaAtmoFlux(GenParam * GenPar);

  // protected methods
  bool    GenerateNext_1try (void);
  void    Initialize        (void);
  void    CleanUp           (void);
  void    ResetSelection    (void);
  double  MinEnergy         (void) { return fMinEvCut; }
  double  MaxCosTheta       (void) { return fMaxCtCut; }
  double  MinCosTheta       (void) { return fMinCtCut; }


  TH2D * CreateFluxHisto2D (string name, string title);
  TH3D * CreateFluxHisto3D (string name, string title);

  void FillFuncFlux(unsigned int iFlux, string FuncName);


  // protected data members

  GAstro  * fAstro;
  GPEarth * PEarth;

  double           fMinEv;            ///< minimum energy in input flux files
  double           fMaxEv;            ///< maximum energy in input flux files
  double           fCosThetaMin;      ///< minimum cosTheta in input flux files
  double           fCosThetaMax;      ///< maximum cosTheta in input flux files
  double           fPhiMin;           ///< minimum phi in input flux files
  double           fPhiMax;           ///< maximum phi in input flux files
  unsigned int     fNumLogEvBins;     ///< number of log-e bins in input flux data files
  unsigned int     fNumLogEvBinsPerDecade;  ///< number of log-e bins per decade in input flux data files
  unsigned int     fNumCosThetaBins;  ///< number of cos(theta) bins in input flux data files
  unsigned int     fNumPhiBins;       ///< number of phi bins in input flux data files
  double           fDeltaLogEv;       ///< energy bin width in log-e

  PDGCodeList *    fPdgCList;         ///< input list of neutrino pdg-codes
  int              fgPdgC;            ///< current nu pdg-code
  TLorentzVector   fgP4;              ///< current nu 4-momentum
  TLorentzVector   fgX4;              ///< current nu 4-position
  long int         fNNeutrinos;       ///< number of flux neutrinos thrown so far
  double           fMaxEvCut;         ///< user-defined cut: maximum energy
  double           fMinEvCut;         ///< user-defined cut: minimum energy
  double           fMaxCtCut;         ///< user-difined cut: maximum neutrino cosTheta
  double           fMinCtCut;         ///< user-difined cut: minimum neutrino cosTheta
  double           fRl;               ///< defining flux neutrino generation surface: longitudinal radius
  double           fRt;               ///< defining flux neutrino generation surface: transverse radius

  double           fX0;               ///< defining flux nuetrino generation surface: X coordinate of fRl radius sphere
  double           fY0;               ///< defining flux nuetrino generation surface: Y coordinate of fRl radius sphere
  double           fZ0;               ///< defining flux nuetrino generation surface: Z coordinate of fRl radius sphere


  double           fEv;               ///< generated energy
  double           fTheta;            ///< generated coming direction
  double           fPhi;              ///< generated coming direction

  TRotation        fRotTHz2User;      ///< coord. system rotation: THZ -> Topocentric user-defined
  bool             fGenWeighted;      ///< generate a weighted or unweighted flux?
  double           fSpectralIndex;    ///< power law function used for weighted flux
  bool             fInitialized;      ///< flag to check that initialization is run
  double           fGlobalGenWeight;  ///< generation weight
  double           fGenWeight;        ///< generation weight (fGlobalGenWeight*pow(Ev,fSpectralIndex))
  double           fWeight;           ///< current generated nu weight
  double           fAgen;             ///< area used to shoot the neutrino

  double           fMJD;
  double           fLST;

  GenParam *       fGenPar;

  map<int, vector<FluxFunc1D> >     fFlux1FRange;  ///< flux = f(Ev)                    for each neutrino species
  map<int, vector<FluxFunc2D> >     fFlux2FRange;  ///< flux = f(Ev,cos_theta)          for each neutrino species
  map<int, vector<FluxFunc3D> >     fFlux3FRange;  ///< flux = f(Ev,cos_theta,phi)      for each neutrino species

  map<int, vector<FluxHisto1D> >    fFlux1DRange;  ///< flux = f(Ev)                    for each neutrino species
  map<int, vector<FluxHisto2D> >    fFlux2DRange;  ///< flux = f(Ev,cos_theta)          for each neutrino species
  map<int, vector<FluxHisto3D> >    fFlux3DRange;  ///< flux = f(Ev,cos_theta,phi)      for each neutrino species


////  double       fGCut;              ///< Neutrino direction precut dist

  bool             fNewNeutrino;
  int              fNInt;

  TLorentzVector   fgP4I;              ///< generated nu 4-momentum
  TLorentzVector   fgX4I;              ///< generated nu 4-position
  int              fgPdgCI;            ///< generated nu pdg-code



};

} // flux namespace
} // genie namespace

#endif // _GATMO_FLUX_H_