//____________________________________________________________________________
/*!

\class   genie::flux::GSeaFileFlux

\brief   A flux driver to read probing neutrinos from file

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

\created March 18, 2019

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

*/
//____________________________________________________________________________


#ifndef _GFILE_FLUX_H_
#define _GFILE_FLUX_H_

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

#include <TLorentzVector.h>
#include <TRotation.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

using namespace std;
using namespace genie;


namespace genie {
namespace flux  {

class GSeaFileFlux: public GFluxI {

public :

  virtual ~GSeaFileFlux();

  // methods implementing the GENIE GFluxI interface
  PDGCodeList &    FluxParticles (void) { return *fPdgCList; }

  virtual double                 MaxEnergy     (void);
  virtual bool                   GenerateNext  (void);
  virtual int                    PdgCode       (void) { return fgPdgC;     }
  virtual double                 Weight        (void) { return 1.;    }

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


  virtual vector<GSeaTrack>      ReadFileEvent(void) = 0;

  // 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
  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     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.

  void     ComputeWeights     (GSeaEvent * SeaEvt);
  void     GetEvtTime         (unsigned int * TimeStampSec, unsigned int * TimeStampTick);

  // needed for CORSIKA file driver
  string   GetGenCode         (void) { return fGenPar->GenCode; }

protected:

  // abstract class, ctor hidden
  GSeaFileFlux(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; }


  // protected data members


  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

  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?
  bool             fInitialized;      ///< flag to check that initialization is run
  double           fAgen;             ///< area used to shoot [w1]
  double           fGenWeight;        ///< current generated weight [w2]
  double           fEvtWeight;        ///< current event weight [w3]
  double           fDistaMax;        ///< current event weight [w3]

  GenParam *       fGenPar;

  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_