//____________________________________________________________________________
/*!

\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


*/
//____________________________________________________________________________


#include <cassert>

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


#include "SeaNuDrivers/GAstro.h"
#include "SeaNuDrivers/GSeaFileFlux.h"

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

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

  fInitialized = 0;
  fNewNeutrino = true;

  fGenPar = GenPar;

  this->Initialize();

}
//___________________________________________________________________________
GSeaFileFlux::~GSeaFileFlux()
{
  this->CleanUp();
}
//___________________________________________________________________________
void GSeaFileFlux::SetNewNeutrino(bool NewNeutrino){
  fNewNeutrino=NewNeutrino;
}
//___________________________________________________________________________
double GSeaFileFlux::MaxEnergy(void)
{
  return TMath::Min(fMaxEv, fMaxEvCut);
}
//___________________________________________________________________________

bool GSeaFileFlux::GenerateNext(void)
{

  this->GenerateNext_1try();

  return true;
}
//___________________________________________________________________________
bool GSeaFileFlux::GenerateNext_1try(void)
{

  // Must have run intitialization
  assert(fInitialized);

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

  return true;
}
//___________________________________________________________________________
void GSeaFileFlux::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;

  // Report and exit
  LOG("GSeaFileFlux", pDEBUG)
       << "Set neutrino: "
       << "\n pdg-code: " << fgPdgC
       << "\n p4: " << utils::print::P4AsShortString(&fgP4)
       << "\n x4: " << utils::print::X4AsString(&fgX4);
}
//___________________________________________________________________________
void GSeaFileFlux::ResetNFluxNeutrinos(void)
{
  fNNeutrinos=0;
  return;
}
//___________________________________________________________________________
long int GSeaFileFlux::NFluxNeutrinos(void) const
{
  return fNNeutrinos;
}
//___________________________________________________________________________
void GSeaFileFlux::SetMinEnergy(double emin)
{
  emin = TMath::Max(0., emin);
  fMinEvCut = emin;
}
//___________________________________________________________________________
void GSeaFileFlux::SetMaxEnergy(double emax)
{
  emax = TMath::Max(0., emax);
  fMaxEvCut = emax;
}
//___________________________________________________________________________
void GSeaFileFlux::SetMinCosTheta(double ctmin)
{
  ctmin = TMath::Max(-1.,ctmin);
  fMinCtCut = ctmin;
}
//___________________________________________________________________________
void GSeaFileFlux::SetMaxCosTheta(double ctmax)
{
  ctmax = TMath::Min(1.,ctmax);
  fMaxCtCut = ctmax;
}

//___________________________________________________________________________
void GSeaFileFlux::GenerateWeighted(bool gen_weighted)
{
  fGenWeighted = gen_weighted;
}

//___________________________________________________________________________
void GSeaFileFlux::InitGenBin(double Emin, double Emax, double RL, double RT, double X0, double Y0, double Z0)
{
  this->SetMinEnergy(Emin);
  this->SetMaxEnergy(Emax);

  fRl=RL;
  fRt=RT;
  fX0=X0;
  fY0=Y0;
  fZ0=Z0;
}
//___________________________________________________________________________
void GSeaFileFlux::Clear(Option_t * opt)
{
// Dummy clear method needed to conform to GFluxI interface
//
  LOG("GSeaFileFlux", pERROR) << "No clear method implemented for option:"<< opt;
}

//___________________________________________________________________________
void GSeaFileFlux::SetUserCoordSystem(TRotation & rotation)
{
  fRotTHz2User = rotation;
}
//___________________________________________________________________________
void GSeaFileFlux::Initialize(void)
{
  LOG("GSeaFileFlux", pNOTICE) << "Initializing file flux driver";

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


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

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

}
//___________________________________________________________________________
void GSeaFileFlux::ComputeWeights(GSeaEvent * SeaEvt){

  SeaEvt->DistaMax=fDistaMax;
  SeaEvt->Agen=fAgen;
  SeaEvt->GenWeight=fGenWeight;
  SeaEvt->EvtWeight=fEvtWeight;


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


  return;
}
//___________________________________________________________________________
void GSeaFileFlux::GetEvtTime(unsigned int * TimeStampSec, unsigned int * TimeStampTick) {

  if(fGenPar->GenMJD){
    RandomGen * rnd = RandomGen::Instance();
    double MJD = fGenPar->MJDStart+rnd->RndFlux().Rndm()*(fGenPar->MJDStop-fGenPar->MJDStart);
    GAstro::CalcTimeStamp(MJD, TimeStampSec, TimeStampTick);
  }
  return;
}
//___________________________________________________________________________
void GSeaFileFlux::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.);

  fDistaMax  = 0.;
  fAgen      = 0.;
  fGenWeight = 0.;
  fEvtWeight = 0.;

}
//___________________________________________________________________________
void GSeaFileFlux::CleanUp(void)
{
  LOG("GSeaFileFlux", pNOTICE) << "Cleaning up...";
  if (fPdgCList ) delete fPdgCList;
}