//____________________________________________________________________________
/*!

\class   genie::flux::GSeaAntaresFileFlux

\brief   A flux driver to read incident neutrinos from evt files

\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

*/
//____________________________________________________________________________

#ifdef _ANTARES_ENABLED__


#include <fstream>
#include <cassert>
#include <cstdio>
#include <iostream>
#include <fstream>

#include <TH2D.h>
#include <TH3D.h>
#include <TMath.h>
#include <TF1.h>
#include <TF3.h>
#include <TFile.h>

#ifdef _GENIEVERSIONGTEQ3__
#include "Framework/Conventions/Constants.h"
#include "Framework/Messenger/Messenger.h"
#include "Framework/ParticleData/PDGLibrary.h"
#else
#include "Conventions/Constants.h"
#include "Messenger/Messenger.h"
#include "PDG/PDGLibrary.h"
#endif

#include "SeaNuDrivers/GSeaAntaresFileFlux.h"


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




//___________________________________________________________________________
GSeaAntaresFileFlux::GSeaAntaresFileFlux(GenParam * GenPar) :
GSeaFileFlux(GenPar)
{

  fEvt = new event();

  if(fGenPar->FluxFiles.size()>1){
    LOG("GSeaAntaresFileFlux", pFATAL) << "Only one input file is allowed, exiting...";
    exit(1);
  }

  if(fGenPar->FluxFiles[0].FileNames.size()>1){
    LOG("GSeaAntaresFileFlux", pFATAL) << "Only one input file is allowed, exiting...";
    exit(1);
  }

  for(unsigned i=0; i<fGenPar->NuList.size(); i++)fPdgCList->push_back(fGenPar->NuList[i]);

  fInFile.open(fGenPar->FluxFiles[0].FileNames[0]);
  if(!fInFile.is_open()){
    LOG("GSeaAntaresFileFlux", pFATAL) << "Error in opening file "<<fGenPar->FluxFiles[0].FileNames[0]<<", exiting...";
    exit(1);
  }



  this->LoadHeaderFile();






  // Configure GFileFlux options
  // 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);


  LOG("GSeaAntaresFileFlux", pNOTICE)
       << "Instantiating the Antares file flux driver";


}
//___________________________________________________________________________
GSeaAntaresFileFlux::~GSeaAntaresFileFlux()
{

  fInFile.close();

}
//___________________________________________________________________________
void GSeaAntaresFileFlux::LoadHeaderFile()
{

  bool new_r, new_evt;
  unsigned int nRun, nEvt, Type;

  fEvt->read(fInFile);
  fEvt->info(new_r, new_evt, nRun, nEvt, Type);

  if(new_r){

    unsigned nLine;
    string line;
    int aux;

    nLine=fEvt->ndat("cut_seamuon");
    for(unsigned i=0; i<nLine; i++){
      line = fEvt->next("cut_seamuon");
//      sscanf(line.c_str(),"%lf %lf %lf %lf",&fGenPar->EvMin,&fGenPar->EvMax,&fGenPar->CtMin,&fGenPar->CtMax);

        // dobbiamo controllare che emin e emax sono validi....

    }

    nLine=fEvt->ndat("norma");
    line = fEvt->next("norma");
    sscanf(line.c_str(),"%d %lf ",&aux,&fGenPar->NTot);

  }
  else{
    LOG("GSeaAntaresFileFlux", pFATAL) << "Evt file does not have a header, exiting...";
    exit(1);
  }

  fMinEv = fGenPar->EvMin;
  fMaxEv = fGenPar->EvMax;
  fCosThetaMin = fGenPar->CtMin;
  fCosThetaMax = fGenPar->CtMax;
  fPhiMin = 0.;
  fPhiMax = 2.*kPi;

  this->SetMinEnergy(fGenPar->EvMin);
  this->SetMaxEnergy(fGenPar->EvMax);
  this->ResetNFluxNeutrinos();

  return;
}

//___________________________________________________________________________
vector<GSeaTrack> GSeaAntaresFileFlux::ReadFileEvent()
{

  vector<GSeaTrack> Tracks;
  int id = 0;

  bool new_r, new_evt;
  unsigned int nRun, nEvt, Type;

  unsigned nLine;
  string line;

  double RT=-1.;
  if     (fGenPar->RTOpt.compare("proj")==0) RT = 0.;
  else if(fGenPar->RTOpt.compare("can")==0)  RT = sqrt(pow((fGenPar->Can2-fGenPar->Can1),2)+pow(fGenPar->Can3*2,2))/2.;
  assert(RT>=0.);

  double Xc=0.;
  double Yc=0.;
  double Zc=-fGenPar->SeaBottomRadius+fGenPar->Can1;  // in m
  double R=fGenPar->SeaBottomRadius+fGenPar->SiteDepth;

  RandomGen * rnd = RandomGen::Instance();

  if(fEvt->read(fInFile)==0){

    fEvt->info(new_r, new_evt, nRun, nEvt, Type);

    if(!new_r){

      GSeaTrack AuxTrack,AuxMom,AuxGrandMom;
      fDistaMax=0.;

      nLine=fEvt->ndat("track_primary");
      line=fEvt->next("track_primary");
      sscanf(line.c_str(),"%*d %lf %lf %lf %lf %lf %lf %lf %lf %d ",&AuxTrack.V1,&AuxTrack.V2,&AuxTrack.V3,&AuxTrack.D1,&AuxTrack.D2,&AuxTrack.D3,&AuxTrack.E,&AuxTrack.T,&AuxTrack.PdgCode);
      AuxTrack.Id = id++;
      AuxTrack.MotherId=-1;
      AuxTrack.Status=0;
      AuxTrack.V1 = TMath::Abs(AuxTrack.V3/AuxTrack.D3)*-1*AuxTrack.D1;
      AuxTrack.V2 = TMath::Abs(AuxTrack.V3/AuxTrack.D3)*-1*AuxTrack.D2;
      AuxTrack.Length = 0.;
      Tracks.push_back(AuxTrack);

      nLine=fEvt->ndat("track_sealepton");
      for(unsigned i=0; i<nLine; i++){
        line = fEvt->next("track_sealepton");
        sscanf(line.c_str(),"%*d %lf %lf %lf %lf %lf %lf %lf %lf %d %lf %*f %d %d",&AuxTrack.V1,&AuxTrack.V2,&AuxTrack.V3,&AuxTrack.D1,&AuxTrack.D2,&AuxTrack.D3,&AuxTrack.E,&AuxTrack.T,&AuxTrack.PdgCode,&AuxTrack.Length,&AuxMom.PdgCode,&AuxGrandMom.PdgCode);
        if ( pdg::IsNeutrino(TMath::Abs(AuxTrack.PdgCode)) && (AuxTrack.E<fGenPar->EvMin || AuxTrack.E>fGenPar->EvMax) ) continue;
        AuxGrandMom.Id       = id++;
        AuxGrandMom.MotherId = Tracks[0].Id;
        AuxGrandMom.Status   = -21;
        AuxMom.Id            = id++;
        AuxMom.MotherId      = AuxGrandMom.Id;
        AuxMom.Status        = -21;
        AuxTrack.Id          = id++;
        AuxTrack.MotherId    = AuxMom.Id;
        AuxTrack.Status      = -1;
        Tracks.push_back(AuxGrandMom);
        Tracks.push_back(AuxMom);
        Tracks.push_back(AuxTrack);

        // only muons that have sufficient range are allowed to affect the generation area:
        double fDistaCan = fabs((fGenPar->SiteDepth-fGenPar->Can2-AuxTrack.V3)/AuxTrack.D3); // fGenPar->Can2 = Zmax
        double Range = pow(10.,fGenPar->RangeSW.Eval(log10(AuxTrack.E)))/(fGenPar->RhoSW);
        if (Range>=fDistaCan) fDistaMax = TMath::Max(fDistaMax,TMath::Sqrt(TMath::Power(AuxTrack.V1,2)+TMath::Power(AuxTrack.V2,2))*TMath::Abs(Tracks[0].D3));
      }

      fDistaMax += 100;

      // generate the primary vertex

      double b  = Tracks[0].D1*(fX0-Xc) + Tracks[0].D2*(fY0-Yc) + Tracks[0].D3*(fZ0-Zc);
      double c  = TMath::Power(fX0-Xc,2)+TMath::Power(fY0-Yc,2)+TMath::Power(fZ0-Zc,2)-R*R;
      double RL = TMath::Abs(-b-TMath::Sqrt(b*b-c));

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

      if (RT==0.) {
        double height = (fGenPar->Can2-fGenPar->Can1+fGenPar->HBedRock) + fDistaMax;
        double radius = fGenPar->Can3 + fDistaMax;
        double AreaTop  = kPi * TMath::Power(radius,2) * TMath::Abs(Tracks[0].D3) ;
        double AreaSide = 2 * height * radius * TMath::Sqrt( 1. - TMath::Power(Tracks[0].D3,2) );
        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 = Tracks[0].D3 > 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( Tracks[0].D2, Tracks[0].D1 );
          v.RotateZ( phidir );
          x = v.X();
          y = v.Y();
          z = v.Z();
        }
        x += -RL * Tracks[0].D1;
        y += -RL * Tracks[0].D2;
        z += -RL * Tracks[0].D3;
      }
      else {
        RT+=fDistaMax;
        fAgen = kPi*TMath::Power(RT,2);
        x += -RL * Tracks[0].D1;  // we have to multiply RL times the arrival direction cosines
        y += -RL * Tracks[0].D2;
        z += -RL * Tracks[0].D3;
        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)*RT*TMath::Cos(psi));
        dvec2.SetMag(TMath::Sqrt(random)*RT*TMath::Sin(psi));
        x += dvec1.X() + dvec2.X();
        y += dvec1.Y() + dvec2.Y();
        z += dvec1.Z() + dvec2.Z();
      }

      x += fX0;
      y += fY0;
      z += fZ0;

      nLine=fEvt->ndat("weights");
      for(unsigned i=0; i<nLine; i++){
        double waux;
        line = fEvt->next("weights");
        sscanf(line.c_str(),"%lf %lf %lf",&waux,&fGenWeight,&fEvtWeight);
        fGenWeight*=fAgen;
        fEvtWeight*=fAgen;

        fEvt->tagd("weights", i);
        char ch[1024];
        sprintf(ch, "%12.3E %12.3E %12.3E",fAgen,fGenWeight,fEvtWeight);
        fEvt->taga("weights",ch);

        sprintf(ch, "%f",fDistaMax);
        fEvt->taga("dista_max",ch);

        sprintf(ch, "%f",RT);
        fEvt->taga("RADD",ch);

      }

      for (auto& Track : Tracks) {
        Track.V1 += x;
        Track.V2 += y;
        Track.V3 += z;
        // for muons and neutrinos, we translated the starting point from the generation area up to the sealevel
        if ( pdg::IsNeutrino(TMath::Abs(Track.PdgCode)) || pdg::IsMuon(TMath::Abs(Track.PdgCode)) ) {
          double bmu  = Track.D1*(Track.V1-Xc)+ Track.D2*(Track.V2-Yc) +Track.D3*(Track.V3-Zc);
          double cmu  = TMath::Power(Track.V1-Xc,2)+TMath::Power(Track.V2-Yc,2)+TMath::Power(Track.V3-Zc,2)-R*R;
          double RLmu = -bmu-TMath::Sqrt(bmu*bmu-cmu);
          Track.V1+=RLmu*Track.D1;
          Track.V2+=RLmu*Track.D2;
          Track.V3+=RLmu*Track.D3;
        }

      }
    }
  }

  return Tracks;
}

#endif