//____________________________________________________________________________
/*!

\class   genie::GGenerateEvent.cxx

\brief   A base generation driver for an under-water neutrino telescope

\author  Carla Distefano <distefano_c@lns.infn.it>
         LNS-INFN, Catania
         Included modifications by Alfonso Garcia <alfonsog@nikhef.nl>
         NIKHEF, Amsterdam (January 2019)
         Included JPP propagation by Victor Carretero <victor.carretero@ific.uv.es>
         IFIC, Valencia (October 2020)
\created July 19, 2017

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

*/
//____________________________________________________________________________


#include <sys/types.h>
#include <sys/stat.h>

#ifdef _GENIEVERSIONGTEQ3__
#include "Framework/Conventions/Constants.h"
#else
#include "Conventions/Constants.h"
#endif



#include "SeaNuDrivers/GGenerateEvent.h"


#if (_PROPOSAL_VERSION_MAJOR__  == 5 && !(_PROPOSAL_VERSION_MINOR__ == 0 && _PROPOSAL_VERSION_PATCH__ == 0)) || (_PROPOSAL_VERSION_MAJOR__ == 6 && _PROPOSAL_VERSION_MINOR__ == 0 && _PROPOSAL_VERSION_PATCH__ == 0 )
#include "PROPOSAL/json.hpp"
#endif

#include <boost/property_tree/json_parser.hpp>



#ifdef _MUSIC_ENABLED__
extern "C"              { void initialize_music_sr_(int *muin);}
extern "C"              { void initialize_music_sw_(int *muin);}
extern "C"              { void muon_transport_sr_(double *VX,double *VY,double *VZ,double *CX,double *CY,double *CZ,double *E,double *Depth,double *T,double *Rho,int *flag1,int *flag2,double *VXF,double *VYF,double *VZF,double *CXF,double *CYF,double *CZF,double *EF,double *DepthF,double *TF);}
extern "C"              { void muon_transport_sw_(double *VX,double *VY,double *VZ,double *CX,double *CY,double *CZ,double *E,double *Depth,double *T,double *Rho,int *flag1,int *flag2,double *VXF,double *VYF,double *VZF,double *CXF,double *CYF,double *CZF,double *EF,double *DepthF,double *TF);}
#endif

#ifdef _HETAU_ENABLED__
#include "Tauola/Tauola.h"
#include "Tauola/TauolaHEPEVTParticle.h"
using namespace Tauolapp;
int TAUMODEL=1;
int ITFLAG=1;
double lifetime = 2.91e-13; //s
double cSpeed = kLightSpeed/(units::meter/units::second);
extern "C"              { void initialize_tausic_sr_(int *tauin);}
extern "C"              { void initialize_tausic_sw_(int *tauin);}
extern "C"              { void tau_transport_sr_(double *VX,double *VY,double *VZ,double *CX,double *CY,double *CZ,double *E,double *Depth,double *T,double *Rho,int *flag1,int *flag2,double *VXF,double *VYF,double *VZF,double *CXF,double *CYF,double *CZF,double *EF,double *DepthF,double *TF, int *IDEC, int *ITFLAG, double *TAUTI, double *TAUTF);}
extern "C"              { void tau_transport_sw_(double *VX,double *VY,double *VZ,double *CX,double *CY,double *CZ,double *E,double *Depth,double *T,double *Rho,int *flag1,int *flag2,double *VXF,double *VYF,double *VZF,double *CXF,double *CYF,double *CZF,double *EF,double *DepthF,double *TF, int *IDEC, int *ITFLAG, double *TAUTI, double *TAUTF);}
#endif


//___________________________________________________________________________
void copy( const GSeaTrack& from , PropaTrack& to )
{
  to.Id             = from.Id;
  to.Status         = from.Status;
  to.Pdg            = from.PdgCode;
  to.MotherId       = from.MotherId;
  to.D1             = from.D1;
  to.D2             = from.D2;
  to.D3             = from.D3;
  to.V1             = from.V1;
  to.V2             = from.V2;
  to.V3             = from.V3;
  to.E              = from.E;
  to.Length         = from.Length;
  to.T              = from.T;
}

//___________________________________________________________________________
void copy( const  PropaTrack& from, GSeaTrack& to)
{
  to.Id             = from.Id;
  to.Status         = from.Status;
  to.PdgCode        = from.Pdg;
  to.MotherId       = from.MotherId;
  to.D1             = from.D1;
  to.D2             = from.D2;
  to.D3             = from.D3;
  to.V1             = from.V1;
  to.V2             = from.V2;
  to.V3             = from.V3;
  to.E              = from.E;
  to.Length         = from.Length;
  to.T              = from.T;
}



//___________________________________________________________________________
GGenerateEvent::GGenerateEvent(GenParam * GenPar)
{
  fGenPar = GenPar;
  fSeaEvent = new GSeaEvent();
  fNEvt=0;
  fNTot=0;

  EoF = false;

  // get flux driver
  flux_driver_diffuse=0;
  flux_driver_point=0;
  fFluxDriver=0;
  fFileDriver=0;

#ifdef _ANTARES_ENABLED__
  fAntaresEvt=0;
#endif

  this->CalcMediaParams();

  if(fGenPar->GenMode.compare("DIFFUSE")==0){
    flux_driver_diffuse = this->GetFluxDiffuse();
    fFluxDriver=dynamic_cast<GSeaAtmoFlux *>(flux_driver_diffuse);
  }
  else if(fGenPar->GenMode.compare("POINT")==0) {
    flux_driver_point = this->GetFluxPoint();
    fFluxDriver=dynamic_cast<GSeaAtmoFlux *>(flux_driver_point);
  }
  else if (fGenPar->GenMode.compare("BIN")==0) {
    LOG("GGenerateEvent", pINFO) <<"Generating an event from CORSIKA input";
    GSeaCORSIKAFileFlux * corsika_file_driver = new GSeaCORSIKAFileFlux(fGenPar);
    fFileDriver = dynamic_cast<GSeaFileFlux *>(corsika_file_driver);
  }
#ifdef _ANTARES_ENABLED__
  else if(fGenPar->GenMode.compare("EVT")==0){
    GSeaAntaresFileFlux * antares_file_driver = new GSeaAntaresFileFlux(fGenPar);
    fFileDriver = dynamic_cast<GSeaFileFlux *>(antares_file_driver);
    fAntaresEvt=antares_file_driver->GetEvent();
  }
#endif
  else {
    LOG("GGenerateEvent", pFATAL) << "Unknown generation mode "<<fGenPar->GenMode<<", exiting...";
    exit(1);
  }


  if(fGenPar->TrackEvts>0){     // we initialize the muon propagator only if we use them
    if(fGenPar->PropCode.compare("PropaMuon")==0){
      PropaMuonSW = new PropaMuon(fGenPar->SeaWaterComp,fGenPar->RhoSW);
      PropaMuonSR = new PropaMuon(fGenPar->RockComp,fGenPar->RhoSR);
    }
#ifdef _MUSIC_ENABLED__
   else if(fGenPar->PropCode.compare("MUSIC")==0){
     string music_path;
     if(gSystem->Getenv("MUSICPATH"))music_path = gSystem->Getenv("MUSICPATH");
     else {
       music_path = gSystem->Getenv("GSEAGEN");
       music_path.append("/dat");

       music_path = BOOST_PP_STRINGIZE(_MUSICDATA__);
       while(1){
         size_t pos = music_path.find("\"");
         if(pos==string::npos)break;
         music_path.erase(pos,1);
       }

     }
     string cp_command("cp ");
     cp_command.append(music_path);
     cp_command.append("/music*.dat .");
     system(cp_command.c_str());

     int muin =2;
     initialize_music_sr_(&muin);
     initialize_music_sw_(&muin);
  }
#endif
#ifdef _PROPOSAL_ENABLED__
   else if(fGenPar->PropCode.compare("PROPOSAL")==0){
     this->ConfigProposal();
   }
#endif
  }

#ifdef _HETAU_ENABLED__
  if (fGenPar->TrackEvts==2) {
    Tauola::setSameParticleDecayMode(2);
    Tauola::setOppositeParticleDecayMode(2);
  }
  Tauola::initialize();
  if (fGenPar->TrackEvts==3) Tauola::setTauBr(2,0.);
  string tausic_path;
  if(gSystem->Getenv("TAUSICPATH"))tausic_path = gSystem->Getenv("TAUSICPATH");
  else {
   tausic_path = gSystem->Getenv("GSEAGEN");
   tausic_path.append("/dat");
   tausic_path = BOOST_PP_STRINGIZE(_TAUSICDATA__);
   while(1){
     size_t pos = tausic_path.find("\"");
     if(pos==string::npos)break;
     tausic_path.erase(pos,1);
   }
  }
  string cp_command("cp ");
  cp_command.append(tausic_path);
  cp_command.append("/tausic*.dat .");
  system(cp_command.c_str());
  int tauin =2;
  initialize_tausic_sr_(&tauin);
  initialize_tausic_sw_(&tauin);
#endif


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


  if(fFluxDriver)delete fFluxDriver;
#ifdef _ANTARES_ENABLED__
  if(fFileDriver)delete fFileDriver;
#endif

  if(fGenPar->TrackEvts){
    if(fGenPar->PropCode.compare("PropaMuon")==0){
      if(PropaMuonSW)delete PropaMuonSW;
      if(PropaMuonSR)delete PropaMuonSR;
    }
#ifdef _PROPOSAL_ENABLED__
    if(fGenPar->PropCode.compare("PROPOSAL")==0){
      if(ProposalMu)delete ProposalMu;
    }
#endif
  }

  if(fSeaEvent)delete fSeaEvent;

}
//___________________________________________________________________________
void GGenerateEvent::CalcMediaParams(void){
// reads the medium compositions and calculate the molar masses, the number of nucleons and electrons of the media
// fill the ComponentSW and ComponentSR if Proposal if enabled at compiling and choosen at run-times

  map<int, double> AtomicMass;
  map<int, string> ElementName;

// load isotope table
  string filename = string(gSystem->Getenv("GENIE")) +
                  string("/data/evgen/catalogues/iso/natural-isotopes.data");


  bool is_accessible = ! (gSystem->AccessPathName( filename.c_str() ));
  if (!is_accessible) {
     LOG("gSeaGen", pFATAL) << "Can not read file: " << filename <<", exiting... ";
     exit(1);
  }

  // load the natural isotopes .txt file
  string input_buf;
  std::ifstream input(filename.c_str());
  if (input.is_open()){

    //skip first 8 lines (comments)
    for(int i=0; i<8; i++){
      string buffer;
      getline(input, buffer);
    }

    int Z = -1, Z_previous = -1, nelements = 0, pdgcode = 0;
    double atomicmass = 0, abundance = 0;
    string elementname, subelementname;

    while( !input.eof() ) {

      //read in naturally occuring element info
      input >> Z;
      input >> elementname;
      input >> nelements;

      // check not re-reading same element
      if(Z!=Z_previous){
        LOG("NatIsotop", pDEBUG) << "Reading entry for Z = " << Z;
        for(int n=0 ; n < nelements; n++){
          input >> subelementname;
          input >> pdgcode;
          input >> atomicmass;
          input >> abundance;

          AtomicMass.insert(map<int, double>::value_type(pdgcode,atomicmass));
          ElementName.insert(map<int, string>::value_type(pdgcode,elementname));
        }
      }
      Z_previous = Z;
    } //!eof

  }
  else {
    LOG("gSeaGen", pFATAL) << "Can not open file: " << filename <<", exiting... ";
    exit(1);
  }

  // Br80 is not in the table but used in the default compositions, we add by hand
  if(AtomicMass.find(1000350800) == AtomicMass.end()){
    AtomicMass.insert(map<int, double>::value_type(1000350800,79.9185293));
    ElementName.insert(map<int, string>::value_type(1000350800,"Br"));
  }


  // SeaWater

  int TgtPdgCode;
  double TgtMassFrac;
  double AtomInMolecule;
  double MoleMin,mass;
  int nAtoms;


  double TgtMassFracMax=0.;
  double MoliFracMass=1.;

  vector<double> AtomMass,Mole,MassFrac;
  vector<int> Z,A;
  vector<string> AtomName;

#ifdef _JPP_ENABLED__
  JPHYSICS::JRadiationSource::source_type EErad(&JPHYSICS::JRadiation::TotalCrossSectionEErad, &JPHYSICS::JRadiation::EfromEErad);
  JPHYSICS::JRadiationSource::source_type Brems(&JPHYSICS::JRadiation::TotalCrossSectionBrems, &JPHYSICS::JRadiation::EfromBrems);
  JPHYSICS::JRadiationSource::source_type GNrad(&JPHYSICS::JRadiation::TotalCrossSectionGNrad, &JPHYSICS::JRadiation::EfromGNrad);
  JPHYSICS::JRadiationSource::source_type Ion  (&JPHYSICS::JRadiation::CalculateACoeff       , NULL);
#endif

  map<int,double>::iterator iter = fGenPar->SeaWaterComp.begin();

  for ( ; iter != fGenPar->SeaWaterComp.end(); ++iter) {
    TgtPdgCode=iter->first;
    TgtMassFrac=iter->second;
    if(AtomicMass.find(TgtPdgCode) != AtomicMass.end()){

      if(TgtMassFrac>TgtMassFracMax){
        TgtMassFracMax=TgtMassFrac;
        MoliFracMass=TgtMassFrac/AtomicMass.find(TgtPdgCode)->second;
      }

      MassFrac.push_back(TgtMassFrac);
      mass=AtomicMass.find(TgtPdgCode)->second;
      AtomMass.push_back(AtomicMass.find(TgtPdgCode)->second);
      AtomName.push_back(ElementName.find(TgtPdgCode)->second);
      A.push_back(pdg::IonPdgCodeToA(TgtPdgCode));
      Z.push_back(pdg::IonPdgCodeToZ(TgtPdgCode));
      Mole.push_back(TgtMassFrac/mass);
    }
    else {
      LOG("gSeaGen", pFATAL) << "Element "<< TgtPdgCode<<" not found in file " << filename <<", exiting... ";
      exit(1);
    }
  }

  MoleMin=*std::min_element(Mole.begin(),Mole.end());

  fGenPar->MolarMassSW=0.;
  fGenPar->NNucleonsSW=0;
  fGenPar->NElectronsSW=0;


  for(unsigned iComp=0; iComp<Mole.size(); iComp++){
    nAtoms=round(Mole[iComp]/MoleMin);
    fGenPar->MolarMassSW+=(double)nAtoms*AtomMass[iComp];
    fGenPar->NNucleonsSW+=nAtoms*A[iComp];
    fGenPar->NElectronsSW+=nAtoms*Z[iComp];

    AtomInMolecule=MassFrac[iComp]/AtomMass[iComp]/MoliFracMass;
#ifdef _PROPOSAL_ENABLED__
      if(fGenPar->PropCode.compare("PROPOSAL")==0)
#if _PROPOSAL_VERSION_MAJOR__  == 5 || (_PROPOSAL_VERSION_MAJOR__ == 6 && _PROPOSAL_VERSION_MINOR__ == 0 && _PROPOSAL_VERSION_PATCH__ == 0 )
        ComponentSW.push_back(shared_ptr<PROPOSAL::Components::Component>(new PROPOSAL::Components::Component(AtomName[iComp].c_str(),Z[iComp],AtomMass[iComp],AtomInMolecule)));
#else
        ComponentSW.push_back(PROPOSAL::Components::Component(AtomName[iComp].c_str(),Z[iComp],AtomMass[iComp],AtomInMolecule));
#endif
#endif
#ifdef _JPP_ENABLED__
      if(fGenPar->PropCode.compare("JPP")==0){
        JPHYSICS::JRadiation jppelement(Z[iComp],A[iComp],40,0.01,0.1,0.1); // (Z,A,integration steps, minimal energy for brems, minimal energy for epair, minimal energy for nuclear interaction)
        JLANG::JSharedPointer<JPHYSICS::JRadiation> jppEl(new JPHYSICS::JRadiationFunction(jppelement, 300, 0.2, 1.0e11));                      //(element, number of energy bins, min energy, max energy)
        radiationSW.push_back( new JPHYSICS::JRadiationSource(jppEl, fGenPar->RhoSW * MassFrac[iComp],  EErad));
        radiationSW.push_back( new JPHYSICS::JRadiationSource(jppEl, fGenPar->RhoSW * MassFrac[iComp],  Brems));
        radiationSW.push_back( new JPHYSICS::JRadiationSource(jppEl, fGenPar->RhoSW * MassFrac[iComp],  GNrad));
        ionizationSW.push_back(new JPHYSICS::JRadiationSource(jppEl, fGenPar->RhoSW * MassFrac[iComp],  Ion));
      }
#endif
  }

  // Rock


  TgtMassFracMax=0.;
  MoliFracMass=1.;


  AtomMass.clear(); Mole.clear(); MassFrac.clear();
  Z.clear(); A.clear();
  AtomName.clear();

  iter = fGenPar->RockComp.begin();


  for ( ; iter != fGenPar->RockComp.end(); ++iter) {
    TgtPdgCode=iter->first;
    TgtMassFrac=iter->second;
    if(AtomicMass.find(TgtPdgCode) != AtomicMass.end()){

      if(TgtMassFrac>TgtMassFracMax){
        TgtMassFracMax=TgtMassFrac;
        MoliFracMass=TgtMassFrac/AtomicMass.find(TgtPdgCode)->second;
      }

      MassFrac.push_back(TgtMassFrac);
      mass=AtomicMass.find(TgtPdgCode)->second;
      AtomMass.push_back(AtomicMass.find(TgtPdgCode)->second);
      AtomName.push_back(ElementName.find(TgtPdgCode)->second);
      A.push_back(pdg::IonPdgCodeToA(TgtPdgCode));
      Z.push_back(pdg::IonPdgCodeToZ(TgtPdgCode));
      Mole.push_back(TgtMassFrac/mass);
    }
    else {
      LOG("gSeaGen", pFATAL) << "Element "<< TgtPdgCode<<" not found in file " << filename <<", exiting... ";
      exit(1);
    }
  }

  MoleMin=*std::min_element(Mole.begin(),Mole.end());

  fGenPar->MolarMassSR=0.;
  fGenPar->NNucleonsSR=0;
  fGenPar->NElectronsSR=0;


  for(unsigned iComp=0; iComp<Mole.size(); iComp++){
    nAtoms=round(Mole[iComp]/MoleMin);
    fGenPar->MolarMassSR+=(double)nAtoms*AtomMass[iComp];
    fGenPar->NNucleonsSR+=nAtoms*A[iComp];
    fGenPar->NElectronsSR+=nAtoms*Z[iComp];

    AtomInMolecule=MassFrac[iComp]/AtomMass[iComp]/MoliFracMass;
#ifdef _PROPOSAL_ENABLED__
      if(fGenPar->PropCode.compare("PROPOSAL")==0)


#if _PROPOSAL_VERSION_MAJOR__  == 5 || (_PROPOSAL_VERSION_MAJOR__ == 6 && _PROPOSAL_VERSION_MINOR__ == 0 && _PROPOSAL_VERSION_PATCH__ == 0 )
        ComponentSR.push_back(shared_ptr<PROPOSAL::Components::Component>(new PROPOSAL::Components::Component(AtomName[iComp].c_str(),Z[iComp],AtomMass[iComp],AtomInMolecule)));
#else
        ComponentSR.push_back(PROPOSAL::Components::Component(AtomName[iComp].c_str(),Z[iComp],AtomMass[iComp],AtomInMolecule));
#endif
#endif

#ifdef _JPP_ENABLED__
      if(fGenPar->PropCode.compare("JPP")==0){
        JPHYSICS::JRadiation jppelement(Z[iComp],A[iComp],40,0.01,0.1,0.1); // (Z,A,integration steps, minimal energyfor brems, minimal energy for epair, minimal energy for nuclear interaction)
        JLANG::JSharedPointer<JPHYSICS::JRadiation> jppEl(new JPHYSICS::JRadiationFunction(jppelement, 300, 0.2, 1.0e11));                      //(element, number of energy bins, min energy, max energy)
        radiationSR.push_back( new JPHYSICS::JRadiationSource(jppEl, fGenPar->RhoSR * MassFrac[iComp],  EErad));
        radiationSR.push_back( new JPHYSICS::JRadiationSource(jppEl, fGenPar->RhoSR * MassFrac[iComp],  Brems));
        radiationSR.push_back( new JPHYSICS::JRadiationSource(jppEl, fGenPar->RhoSR * MassFrac[iComp],  GNrad));
        ionizationSR.push_back(new JPHYSICS::JRadiationSource(jppEl, fGenPar->RhoSR * MassFrac[iComp],  Ion));
      }
#endif
  }


  // Mantle

  AtomMass.clear(); Mole.clear(); MassFrac.clear();
  Z.clear(); A.clear();
  AtomName.clear();

  iter = fGenPar->MantleComp.begin();
  for ( ; iter != fGenPar->MantleComp.end(); ++iter) {
    TgtPdgCode=iter->first;
    TgtMassFrac=iter->second;
    if(AtomicMass.find(TgtPdgCode) != AtomicMass.end()){

      if(TgtMassFrac>TgtMassFracMax){
        TgtMassFracMax=TgtMassFrac;
        MoliFracMass=TgtMassFrac/AtomicMass.find(TgtPdgCode)->second;
      }

      MassFrac.push_back(TgtMassFrac);
      mass=AtomicMass.find(TgtPdgCode)->second;
      AtomMass.push_back(AtomicMass.find(TgtPdgCode)->second);
      AtomName.push_back(ElementName.find(TgtPdgCode)->second);
      A.push_back(pdg::IonPdgCodeToA(TgtPdgCode));
      Z.push_back(pdg::IonPdgCodeToZ(TgtPdgCode));
      Mole.push_back(TgtMassFrac/mass);
    }
    else {
      LOG("gSeaGen", pFATAL) << "Element "<< TgtPdgCode<<" not found in file " << filename <<", exiting... ";
      exit(1);
    }
  }

  MoleMin=*std::min_element(Mole.begin(),Mole.end());

  fGenPar->MolarMassMantle=0.;
  fGenPar->NNucleonsMantle=0;
  fGenPar->NElectronsMantle=0;


  for(unsigned iComp=0; iComp<Mole.size(); iComp++){
    nAtoms=round(Mole[iComp]/MoleMin);
    fGenPar->MolarMassMantle+=(double)nAtoms*AtomMass[iComp];
    fGenPar->NNucleonsMantle+=nAtoms*A[iComp];
    fGenPar->NElectronsMantle+=nAtoms*Z[iComp];
  }


  // Core


  AtomMass.clear(); Mole.clear(); MassFrac.clear();
  Z.clear(); A.clear();
  AtomName.clear();

  iter = fGenPar->CoreComp.begin();
  for ( ; iter != fGenPar->CoreComp.end(); ++iter) {
    TgtPdgCode=iter->first;
    TgtMassFrac=iter->second;
    if(AtomicMass.find(TgtPdgCode) != AtomicMass.end()){

      if(TgtMassFrac>TgtMassFracMax){
        TgtMassFracMax=TgtMassFrac;
        MoliFracMass=TgtMassFrac/AtomicMass.find(TgtPdgCode)->second;
      }

      MassFrac.push_back(TgtMassFrac);
      mass=AtomicMass.find(TgtPdgCode)->second;
      AtomMass.push_back(AtomicMass.find(TgtPdgCode)->second);
      AtomName.push_back(ElementName.find(TgtPdgCode)->second);
      A.push_back(pdg::IonPdgCodeToA(TgtPdgCode));
      Z.push_back(pdg::IonPdgCodeToZ(TgtPdgCode));
      Mole.push_back(TgtMassFrac/mass);
    }
    else {
      LOG("gSeaGen", pFATAL) << "Element "<< TgtPdgCode<<" not found in file " << filename <<", exiting... ";
      exit(1);
    }
  }

  MoleMin=*std::min_element(Mole.begin(),Mole.end());

  fGenPar->MolarMassCore=0.;
  fGenPar->NNucleonsCore=0;
  fGenPar->NElectronsCore=0;


  for(unsigned iComp=0; iComp<Mole.size(); iComp++){
    nAtoms=round(Mole[iComp]/MoleMin);
    fGenPar->MolarMassCore+=(double)nAtoms*AtomMass[iComp];
    fGenPar->NNucleonsCore+=nAtoms*A[iComp];
    fGenPar->NElectronsCore+=nAtoms*Z[iComp];
  }

  AtomMass.clear(); Mole.clear(); MassFrac.clear();
  Z.clear(); A.clear();
  AtomName.clear();

  return;
}
//___________________________________________________________________________
void GGenerateEvent::ConfigProposal(void){
// configure Proposal if enabled at compiling and choosen at run-times
#ifdef _PROPOSAL_ENABLED__

  double global_ecut_inside  = 500;
  double global_ecut_infront = -1;
  double global_ecut_behind  = -1;

  double global_vcut_inside  = -1;
  double global_vcut_infront = 0.001;
  double global_vcut_behind  = -1;

  bool global_cont_inside  = false;
  bool global_cont_infront = true;
  bool global_cont_behind = false;

  bool do_interpolation = true;

  string path_to_tables_readonly=PROPOSALRES;
  string path_to_tables=PROPOSALRES;

  if(!gSystem->Getenv("PROPOSALCONF")){
    struct stat info;
    if( stat( path_to_tables.c_str(), &info ) != 0 ){
      LOG("GGenerateEvent", pFATAL) <<"Error: invalid PROPOSAL table path "<<path_to_tables<<"; exiting...";
      exit(1);
    }
    else if( !(info.st_mode & S_IFDIR) ) {
      LOG("GGenerateEvent", pFATAL) <<"Error: invalid PROPOSAL table path "<<path_to_tables<<"; exiting...";
      exit(1);
    }
  }

///////////////////////////////////////////////////////////


  PROPOSAL::RandomGenerator::Get().SetSeed(fGenPar->RanSeed);

  PROPOSAL::InterpolationDef interpolation_def;

  interpolation_def.path_to_tables=path_to_tables;

#if !(_PROPOSAL_VERSION_MAJOR__ == 5 && _PROPOSAL_VERSION_MINOR__ == 0 && _PROPOSAL_VERSION_PATCH__ == 0 ) && !(_PROPOSAL_VERSION_MAJOR__ == 5 && _PROPOSAL_VERSION_MINOR__ == 1 && _PROPOSAL_VERSION_PATCH__ == 0 )
  interpolation_def.path_to_tables_readonly=path_to_tables_readonly;
#endif

  unique_ptr<PROPOSAL::Sector::Definition> sec_def_global(new PROPOSAL::Sector::Definition());

  if(gSystem->Getenv("PROPOSALCONF")){

    string config_file = gSystem->Getenv("PROPOSALCONF");

    if(!(access(config_file.c_str(), F_OK) != -1)){
      LOG("GGenerateEvent", pFATAL) << config_file <<" does not exist --> check env. parameter PROPOSALCONF; exiting...";
      exit(1);
    }
    else
      LOG("GGenerateEvent", pINFO) <<"Reading PROPOSAL inputs from file "<< config_file;

// check if table paths exist before proposal parses the file to avoid abort

    struct stat info;
    string pathname;
    boost::property_tree::ptree json_config_check;
    boost::property_tree::json_parser::read_json(config_file, json_config_check);
    vector<string> option;

#if (_PROPOSAL_VERSION_MAJOR__ == 5 && _PROPOSAL_VERSION_MINOR__ == 0 && _PROPOSAL_VERSION_PATCH__ == 0 )
    option.push_back("global.path_to_tables");
#else
    option.push_back("global.interpolation.path_to_tables");
    option.push_back("global.interpolation.path_to_tables_readonly");
#endif
    boost::optional<string> optional_param;
    for(unsigned i=0; i<option.size(); i++){
      optional_param = json_config_check.get_optional<string>(option[i]);
      if(optional_param){
        pathname = optional_param.get();
        if( stat( pathname.c_str(), &info ) != 0 ){
          LOG("GGenerateEvent", pFATAL) <<"Error in PROPOSAL input file "<< config_file<<": invalid path "<<pathname<<"; exiting...";
          exit(1);
        }
        else if( !(info.st_mode & S_IFDIR) ) {
          LOG("GGenerateEvent", pFATAL) <<"Error in PROPOSAL input file "<< config_file<<": invalid path "<<pathname<<"; exiting...";
          exit(1);
        }
      }
    }


#if (_PROPOSAL_VERSION_MAJOR__ == 5 && _PROPOSAL_VERSION_MINOR__ == 0 && _PROPOSAL_VERSION_PATCH__ == 0 )

    boost::property_tree::ptree json_config;
    boost::property_tree::json_parser::read_json(config_file, json_config);

    if(json_config.get_optional<bool>("global.interpolate")) do_interpolation = json_config.get_optional<bool>("global.interpolate").get();
    if(json_config.get_optional<bool>("global.raw")) interpolation_def.raw = json_config.get_optional<bool>("global.raw").get();
    if(json_config.get_optional<string>("global.path_to_tables")) interpolation_def.path_to_tables = json_config.get_optional<string>("global.path_to_tables").get();

    if(json_config.get_optional<double>("global.brems_multiplier")) sec_def_global->utility_def.brems_def.multiplier = json_config.get_optional<double>("global.brems_multiplier").get();
    if(json_config.get_optional<double>("global.photo_multiplier")) sec_def_global->utility_def.photo_def.multiplier = json_config.get_optional<double>("global.photo_multiplier").get();
    if(json_config.get_optional<double>("global.ioniz_multiplier")) sec_def_global->utility_def.ioniz_def.multiplier = json_config.get_optional<double>("global.ioniz_multiplier").get();
    if(json_config.get_optional<double>("global.epair_multiplier")) sec_def_global->utility_def.epair_def.multiplier = json_config.get_optional<double>("global.epair_multiplier").get();

    if(json_config.get_optional<bool>("global.lpm")) sec_def_global->utility_def.epair_def.lpm_effect = json_config.get_optional<bool>("global.lpm").get();
    if(json_config.get_optional<bool>("global.lpm")) sec_def_global->utility_def.brems_def.lpm_effect = json_config.get_optional<bool>("global.lpm").get();

    if(json_config.get_optional<bool>("global.exact_time")) sec_def_global->do_exact_time_calculation = json_config.get_optional<bool>("global.exact_time").get();
    if(json_config.get_optional<bool>("global.continous_loss_output")) sec_def_global->do_continuous_energy_loss_output = json_config.get_optional<bool>("global.continous_loss_output").get();

    if(json_config.get_optional<bool>("global.stopping_decay")) sec_def_global->stopping_decay = json_config.get_optional<bool>("global.stopping_decay").get();
    if(json_config.get_optional<bool>("global.do_stochastic_loss_weighting")) sec_def_global->do_stochastic_loss_weighting = json_config.get_optional<bool>("global.do_stochastic_loss_weighting").get();
    if(json_config.get_optional<double>("global.stochastic_loss_weighting")) sec_def_global->stochastic_loss_weighting = json_config.get_optional<double>("global.stochastic_loss_weighting").get();

    if(json_config.get_optional<string>("global.scattering")) {
      string scattering = json_config.get_optional<string>("global.scattering").get();
      sec_def_global->scattering_model =  PROPOSAL::ScatteringFactory::Get().GetEnumFromString(scattering);
    }

    if(json_config.get_optional<string>("global.brems")) {
      string brems = json_config.get_optional<string>("global.brems").get();
      sec_def_global->utility_def.brems_def.parametrization =  PROPOSAL::BremsstrahlungFactory::Get().GetEnumFromString(brems);
    }

    if(json_config.get_optional<string>("global.epair")) {
      string epair = json_config.get_optional<string>("global.epair").get();
      sec_def_global->utility_def.epair_def.parametrization =  PROPOSAL::EpairProductionFactory::Get().GetEnumFromString(epair);
    }

    if(json_config.get_optional<string>("global.photo")) {
      string photo = json_config.get_optional<string>("global.photo").get();
      sec_def_global->utility_def.photo_def.parametrization =  PROPOSAL::PhotonuclearFactory::Get().GetEnumFromString(photo);
    }

    if(json_config.get_optional<string>("global.photo_shadow")) {
      string photo_shadow = json_config.get_optional<string>("global.photo_shadow").get();
      sec_def_global->utility_def.photo_def.shadow =  PROPOSAL::PhotonuclearFactory::Get().GetShadowEnumFromString(photo_shadow);
    }

    if(json_config.get_optional<bool>("global.photo_hard_component")) sec_def_global->utility_def.photo_def.hard_component = json_config.get_optional<bool>("global.photo_hard_component").get();


    if(json_config.get_optional<double>("global.global_ecut_inside"))  global_ecut_inside  = json_config.get_optional<double>("global.global_ecut_inside").get();
    if(json_config.get_optional<double>("global.global_ecut_infront")) global_ecut_infront = json_config.get_optional<double>("global.global_ecut_infront").get();
    if(json_config.get_optional<double>("global.global_ecut_behind"))  global_ecut_behind  = json_config.get_optional<double>("global.global_ecut_behind").get();
    if(json_config.get_optional<double>("global.global_vcut_inside"))  global_vcut_inside  = json_config.get_optional<double>("global.global_vcut_inside").get();
    if(json_config.get_optional<double>("global.global_vcut_infront")) global_vcut_infront = json_config.get_optional<double>("global.global_vcut_infront").get();
    if(json_config.get_optional<double>("global.global_vcut_behind"))  global_vcut_behind  = json_config.get_optional<double>("global.global_vcut_behind").get();
    if(json_config.get_optional<bool>("global.global_cont_inside"))    global_cont_inside  = json_config.get_optional<bool>("global.global_cont_inside").get();
    if(json_config.get_optional<bool>("global.global_cont_infront"))   global_cont_infront = json_config.get_optional<bool>("global.global_cont_infront").get();
    if(json_config.get_optional<bool>("global.global_cont_behind"))    global_cont_behind  = json_config.get_optional<bool>("global.global_cont_behind").get();


#elif _PROPOSAL_VERSION_MAJOR__  == 5 || (_PROPOSAL_VERSION_MAJOR__ == 6 && _PROPOSAL_VERSION_MINOR__ == 0 && _PROPOSAL_VERSION_PATCH__ == 0 )

    // Create the json parser
    nlohmann::json json_config;
    std::ifstream infilestream(config_file);
    infilestream >> json_config;

    nlohmann::json json_global = json_config["global"];

    if(json_global.find("continous_loss_output") != json_global.end()) sec_def_global->do_continuous_energy_loss_output = json_global["continous_loss_output"].get<bool>();
#if (_PROPOSAL_VERSION_MAJOR__ == 6 && _PROPOSAL_VERSION_MINOR__ == 0 && _PROPOSAL_VERSION_PATCH__ == 0 )
    if(json_global.find("only_loss_inside_detector") != json_global.end()) sec_def_global->only_loss_inside_detector = json_global["only_loss_inside_detector"].get<bool>();
#endif

    if (json_global.find("interpolation") != json_global.end()){
      if(json_global["interpolation"].find("do_interpolation") != json_global["interpolation"].end()) do_interpolation = json_global["interpolation"]["do_interpolation"];
      if(json_global["interpolation"].find("path_to_tables") != json_global["interpolation"].end()) interpolation_def.path_to_tables = json_global["interpolation"]["path_to_tables"];
      if(json_global["interpolation"].find("do_binary_tables") != json_global["interpolation"].end()) interpolation_def.do_binary_tables = json_global["interpolation"]["do_binary_tables"];

#if !(_PROPOSAL_VERSION_MAJOR__ == 5 && _PROPOSAL_VERSION_MINOR__ == 1 && _PROPOSAL_VERSION_PATCH__ == 0 )
      if(json_global["interpolation"].find("nodes_propagate") != json_global["interpolation"].end()) interpolation_def.nodes_propagate = json_global["interpolation"]["nodes_propagate"];
      if(json_global["interpolation"].find("nodes_continous_randomization") != json_global["interpolation"].end()) interpolation_def.nodes_continous_randomization = json_global["interpolation"]["nodes_continous_randomization"];
      if(json_global["interpolation"].find("nodes_cross_section") != json_global["interpolation"].end()) interpolation_def.nodes_cross_section = json_global["interpolation"]["nodes_cross_section"];
      if(json_global["interpolation"].find("max_node_energy") != json_global["interpolation"].end()) interpolation_def.max_node_energy = json_global["interpolation"]["max_node_energy"];
      if(json_global["interpolation"].find("path_to_tables_readonly") != json_global["interpolation"].end()) interpolation_def.path_to_tables_readonly = json_global["interpolation"]["path_to_tables_readonly"];
#endif
#if (_PROPOSAL_VERSION_MAJOR__ == 6 && _PROPOSAL_VERSION_MINOR__ == 0 && _PROPOSAL_VERSION_PATCH__ == 0 )
      if(json_global["interpolation"].find("just_use_readonly_path") != json_global["interpolation"].end()) interpolation_def.just_use_readonly_path = json_global["interpolation"]["just_use_readonly_path"];
#endif
    }

    if(json_global.find("exact_time") != json_global.end()) sec_def_global->do_exact_time_calculation = json_global["exact_time"].get<bool>();
    if(json_global.find("stopping_decay") != json_global.end()) sec_def_global->stopping_decay = json_global["stopping_decay"].get<bool>();
    if(json_global.find("scattering") != json_global.end()) {
      string scattering = json_global["scattering"].get<std::string>();
      sec_def_global->scattering_model =  PROPOSAL::ScatteringFactory::Get().GetEnumFromString(scattering);
    }

    if(json_global.find("brems_multiplier") != json_global.end()) sec_def_global->utility_def.brems_def.multiplier = json_global["brems_multiplier"].get<double>();
    if(json_global.find("epair_multiplier") != json_global.end()) sec_def_global->utility_def.epair_def.multiplier = json_global["epair_multiplier"].get<double>();
    if(json_global.find("ioniz_multiplier") != json_global.end()) sec_def_global->utility_def.ioniz_def.multiplier = json_global["ioniz_multiplier"].get<double>();
    if(json_global.find("photo_multiplier") != json_global.end()) sec_def_global->utility_def.photo_def.multiplier = json_global["photo_multiplier"].get<double>();
    if(json_global.find("brems") != json_global.end()) {
      string brems = json_global["brems"].get<std::string>();
      sec_def_global->utility_def.brems_def.parametrization =  PROPOSAL::BremsstrahlungFactory::Get().GetEnumFromString(brems);
    }
    if(json_global.find("epair") != json_global.end()) {
      string epair = json_global["epair"].get<std::string>();
      sec_def_global->utility_def.epair_def.parametrization =  PROPOSAL::EpairProductionFactory::Get().GetEnumFromString(epair);
    }
    if(json_global.find("photo") != json_global.end()) {
      string photo = json_global["photo"].get<std::string>();
      sec_def_global->utility_def.photo_def.parametrization =  PROPOSAL::PhotonuclearFactory::Get().GetEnumFromString(photo);
    }
    if(json_global.find("photo_shadow") != json_global.end()) {
      string photo_shadow = json_global["photo_shadow"].get<std::string>();
      sec_def_global->utility_def.photo_def.shadow =  PROPOSAL::PhotonuclearFactory::Get().GetShadowEnumFromString(photo_shadow);
    }
    if(json_global.find("photo_hard_component") != json_global.end()) sec_def_global->utility_def.photo_def.hard_component = json_global["photo_hard_component"].get<bool>();
    if(json_global.find("lpm") != json_global.end()) sec_def_global->utility_def.epair_def.lpm_effect = json_global["lpm"].get<bool>();
    if(json_global.find("lpm") != json_global.end()) sec_def_global->utility_def.brems_def.lpm_effect = json_global["lpm"].get<bool>();


#if !(_PROPOSAL_VERSION_MAJOR__ == 5 && _PROPOSAL_VERSION_MINOR__ == 1 && _PROPOSAL_VERSION_PATCH__ == 0 )
    if(json_global.find("mupair_multiplier") != json_global.end()) sec_def_global->utility_def.mupair_def.multiplier = json_global["mupair_multiplier"].get<double>();
    if(json_global.find("mupair_particle_output") != json_global.end()) sec_def_global->utility_def.mupair_def.particle_output = json_global["mupair_particle_output"].get<bool>();
#endif
#if (_PROPOSAL_VERSION_MAJOR__ == 5 && _PROPOSAL_VERSION_MINOR__ == 2 && _PROPOSAL_VERSION_PATCH__ == 0 )
    if(json_global.find("mupair_enable") != json_global.end()) sec_def_global->utility_def.mupair_def.mupair_enable = json_global["mupair_enable"].get<bool>();
    if(json_global.find("weak_enable") != json_global.end()) sec_def_global->utility_def.weak_def.weak_enable = json_global["weak_enable"].get<bool>();
#endif


#if (_PROPOSAL_VERSION_MAJOR__ == 6 && _PROPOSAL_VERSION_MINOR__ == 0 && _PROPOSAL_VERSION_PATCH__ == 0 )

    if(json_global.find("weak_multiplier") != json_global.end()) sec_def_global->utility_def.weak_def.multiplier = json_global["weak_multiplier"].get<double>();
    if(json_global.find("annihilation_multiplier") != json_global.end()) sec_def_global->utility_def.annihilation_def.multiplier = json_global["annihilation_multiplier"].get<double>();

    if(json_global.find("mupair") != json_global.end()) {
      string mupair = json_global["mupair"].get<std::string>();
      sec_def_global->utility_def.mupair_def.parametrization =  PROPOSAL::MupairProductionFactory::Get().GetEnumFromString(mupair);
    }

    if(json_global.find("weak") != json_global.end()) {
      string weak = json_global["weak"].get<std::string>();
      sec_def_global->utility_def.weak_def.parametrization =  PROPOSAL::WeakInteractionFactory::Get().GetEnumFromString(weak);
    }


    if(json_global.find("ioniz") != json_global.end()) {
      string ioniz = json_global["ioniz"].get<std::string>();
      sec_def_global->utility_def.ioniz_def.parametrization =  PROPOSAL::IonizationFactory::Get().GetEnumFromString(ioniz);
    }

    if(json_global.find("annihilation") != json_global.end()) {
      string annihilation = json_global["annihilation"].get<std::string>();
      sec_def_global->utility_def.annihilation_def.parametrization =  PROPOSAL::AnnihilationFactory::Get().GetEnumFromString(annihilation);
    }
#endif

    if (json_global.find("cuts_infront") != json_global.end()){
      if(json_global["cuts_infront"].find("e_cut") != json_global["cuts_infront"].end()) global_ecut_infront = json_global["cuts_infront"]["e_cut"];
      if(json_global["cuts_infront"].find("v_cut") != json_global["cuts_infront"].end()) global_vcut_infront = json_global["cuts_infront"]["v_cut"];
      if(json_global["cuts_infront"].find("cont_rand") != json_global["cuts_infront"].end()) global_cont_infront = json_global["cuts_infront"]["cont_rand"];
    }
    if (json_global.find("cuts_inside") != json_global.end()){
      if(json_global["cuts_inside"].find("e_cut") != json_global["cuts_inside"].end()) global_ecut_inside = json_global["cuts_inside"]["e_cut"];
      if(json_global["cuts_inside"].find("v_cut") != json_global["cuts_inside"].end()) global_vcut_inside = json_global["cuts_inside"]["v_cut"];
      if(json_global["cuts_inside"].find("cont_rand") != json_global["cuts_inside"].end()) global_cont_inside = json_global["cuts_inside"]["cont_rand"];
    }
    if (json_global.find("cuts_behind") != json_global.end()){
      if(json_global["cuts_behind"].find("e_cut") != json_global["cuts_behind"].end()) global_ecut_behind = json_global["cuts_behind"]["e_cut"];
      if(json_global["cuts_behind"].find("v_cut") != json_global["cuts_behind"].end()) global_vcut_behind = json_global["cuts_behind"]["v_cut"];
      if(json_global["cuts_behind"].find("cont_rand") != json_global["cuts_behind"].end()) global_cont_behind = json_global["cuts_behind"]["cont_rand"];
    }

    if (json_global.find("stochastic_loss_weighting") != json_global.end()){
      sec_def_global->stochastic_loss_weighting = json_global["stochastic_loss_weighting"].get<double>();
      sec_def_global->do_stochastic_loss_weighting = 1e-10 < std::abs(sec_def_global->stochastic_loss_weighting);
    }


#else
    // Create the json parser
    nlohmann::json json_config;
    std::ifstream infilestream(config_file);
    infilestream >> json_config;

    nlohmann::json json_global;

    if(json_config.contains("global")){
      json_global = json_config["global"];

      if (json_global.contains("interpolation")){
        nlohmann::json json_interpol = json_global["interpolation"];
        interpolation_def = PROPOSAL::InterpolationDef(json_interpol);
        do_interpolation = json_interpol.value("do_interpolation", do_interpolation);
      }
      sec_def_global.reset(new PROPOSAL::Sector::Definition(json_global));

    }
#endif
  }



#if _PROPOSAL_VERSION_MAJOR__  == 5 || (_PROPOSAL_VERSION_MAJOR__ == 6 && _PROPOSAL_VERSION_MINOR__ == 0 && _PROPOSAL_VERSION_PATCH__ == 0 )

  PROPOSAL::Medium *med_sw = new PROPOSAL::Medium("gSeaWater",
           1,
           75.0,    // I
           -3.5017, // C
           0.09116, // a
           3.4773,  // m
           0.2400,  // X0
           2.8004,  // X1
           0,       // d0
           fGenPar->RhoSW, // massDensity
           ComponentSW);


  PROPOSAL::Medium *med_sr = new PROPOSAL::Medium("gRock",
           1,
           149.0,   // I
           -5.053,  // C
           0.078,   // a
           3.645,   // m
           0.288,   // X0
           3.196,   // X1
           0,       // d0
           fGenPar->RhoSR, // massDensity
           ComponentSR);

  PROPOSAL::Cylinder geoCan;
  geoCan.SetPosition(PROPOSAL::Vector3D(0, 0, (fGenPar->Can1+fGenPar->Can2)/2.*100.));
  geoCan.SetRadius(fGenPar->Can3*100.);
  geoCan.SetInnerRadius(0.*100.);
  geoCan.SetZ((fGenPar->Can2-fGenPar->Can1)*100.);

  PROPOSAL::Geometry* geometrySW = PROPOSAL::GeometryFactory::Get().CreateGeometry("cylinder");
  PROPOSAL::Cylinder* cylinderSW = dynamic_cast<PROPOSAL::Cylinder*>(geometrySW);

  cylinderSW->SetPosition(PROPOSAL::Vector3D(0, 0, fGenPar->SeaCenter*100.));
  cylinderSW->SetRadius(fGenPar->RVol*100.);
  cylinderSW->SetInnerRadius(0.*100.);
  cylinderSW->SetZ(fGenPar->HSeaWater*100.);
  PROPOSAL::Geometry* geoSeaWater = cylinderSW;
  geoSeaWater->SetHierarchy(1);

  PROPOSAL::Geometry* geometrySR = PROPOSAL::GeometryFactory::Get().CreateGeometry("cylinder");
  PROPOSAL::Cylinder* cylinderSR = dynamic_cast<PROPOSAL::Cylinder*>(geometrySR);

  cylinderSR->SetPosition(PROPOSAL::Vector3D(0, 0, fGenPar->RockCenter*100.));
  cylinderSR->SetRadius(fGenPar->RVol*100.);
  cylinderSR->SetInnerRadius(0.*100.);
  cylinderSR->SetZ(fGenPar->HRock*100.);
  PROPOSAL::Geometry* geoRock = cylinderSR;
  geoRock->SetHierarchy(1);


#else


  shared_ptr<const PROPOSAL::Medium> med_sw = make_shared<const PROPOSAL::Medium>("gSeaWater",
           1,
           75.0,    // I
           -3.5017, // C
           0.09116, // a
           3.4773,  // m
           0.2400,  // X0
           2.8004,  // X1
           0,       // d0
           fGenPar->RhoSW, // massDensity
           ComponentSW);

  shared_ptr<const PROPOSAL::Medium> med_sr = make_shared<const PROPOSAL::Medium>("gRock",
           1,
           149.0,   // I
           -5.053,  // C
           0.078,   // a
           3.645,   // m
           0.288,   // X0
           3.196,   // X1
           0,       // d0
           fGenPar->RhoSR, // massDensity
           ComponentSR);


  shared_ptr<const PROPOSAL::Geometry> geoCan;
  geoCan = make_shared<PROPOSAL::Cylinder>(PROPOSAL::Vector3D(0, 0, (fGenPar->Can1+fGenPar->Can2)/2.*100.), fGenPar->Can3*100., 0., (fGenPar->Can2-fGenPar->Can1)*100.);


  shared_ptr<PROPOSAL::Geometry> geometrySW;
  geometrySW = make_shared<PROPOSAL::Cylinder>(PROPOSAL::Vector3D(0, 0, fGenPar->SeaCenter*100.), fGenPar->RVol*100., 0., fGenPar->HSeaWater*100.);
  geometrySW->SetHierarchy(1);
  shared_ptr<const PROPOSAL::Geometry> geoSeaWater=geometrySW;


  shared_ptr<PROPOSAL::Geometry> geometrySR;
  geometrySR = make_shared<PROPOSAL::Cylinder>(PROPOSAL::Vector3D(0, 0, fGenPar->RockCenter*100.), fGenPar->RVol*100., 0., fGenPar->HRock*100.);
  geometrySR->SetHierarchy(1);
  shared_ptr<const PROPOSAL::Geometry> geoRock=geometrySR;

#endif


  vector<PROPOSAL::Sector::Definition> sectors;

  PROPOSAL::Sector::Definition sec_def_sector=*sec_def_global;

  sec_def_sector.location = PROPOSAL::Sector::ParticleLocation::InfrontDetector;

#if _PROPOSAL_VERSION_MAJOR__ > 5 && _PROPOSAL_VERSION_MINOR__ > 0
  sec_def_sector.SetMedium(med_sw);
  sec_def_sector.SetGeometry(geoSeaWater);
#else
  sec_def_sector.SetMedium(*med_sw);
  sec_def_sector.SetGeometry(*geoSeaWater);
#endif

  sec_def_sector.cut_settings.SetEcut(global_ecut_infront);
  sec_def_sector.cut_settings.SetVcut(global_vcut_infront);
  sec_def_sector.do_continuous_randomization = global_cont_infront;

  sectors.push_back(sec_def_sector);

#if _PROPOSAL_VERSION_MAJOR__ > 5 && _PROPOSAL_VERSION_MINOR__ > 0
  sec_def_sector.SetMedium(med_sr);
  sec_def_sector.SetGeometry(geoRock);
#else
  sec_def_sector.SetMedium(*med_sr);
  sec_def_sector.SetGeometry(*geoRock);
#endif

  sectors.push_back(sec_def_sector);

#if _PROPOSAL_VERSION_MAJOR__ > 5 && _PROPOSAL_VERSION_MINOR__ > 0
  sec_def_sector.SetMedium(med_sw);
  sec_def_sector.SetGeometry(geoSeaWater);
#else
  sec_def_sector.SetMedium(*med_sw);
  sec_def_sector.SetGeometry(*geoSeaWater);
#endif
  sec_def_sector.location = PROPOSAL::Sector::ParticleLocation::InsideDetector;
  sec_def_sector.cut_settings.SetEcut(global_ecut_inside);
  sec_def_sector.cut_settings.SetVcut(global_vcut_inside);
  sec_def_sector.do_continuous_randomization = global_cont_inside;

  sectors.push_back(sec_def_sector);

#if _PROPOSAL_VERSION_MAJOR__ > 5 && _PROPOSAL_VERSION_MINOR__ > 0
  sec_def_sector.SetMedium(med_sr);
  sec_def_sector.SetGeometry(geoRock);
#else
  sec_def_sector.SetMedium(*med_sr);
  sec_def_sector.SetGeometry(*geoRock);
#endif

  sectors.push_back(sec_def_sector);

#if _PROPOSAL_VERSION_MAJOR__ > 5 && _PROPOSAL_VERSION_MINOR__ > 0
  sec_def_sector.SetMedium(med_sw);
  sec_def_sector.SetGeometry(geoSeaWater);
#else
  sec_def_sector.SetMedium(*med_sw);
  sec_def_sector.SetGeometry(*geoSeaWater);
#endif
  sec_def_sector.location = PROPOSAL::Sector::ParticleLocation::BehindDetector;
  sec_def_sector.cut_settings.SetEcut(global_ecut_behind);
  sec_def_sector.cut_settings.SetVcut(global_vcut_behind);
  sec_def_sector.do_continuous_randomization = global_cont_behind;

  sectors.push_back(sec_def_sector);

#if _PROPOSAL_VERSION_MAJOR__ > 5 && _PROPOSAL_VERSION_MINOR__ > 0
  sec_def_sector.SetMedium(med_sr);
  sec_def_sector.SetGeometry(geoRock);
#else
  sec_def_sector.SetMedium(*med_sr);
  sec_def_sector.SetGeometry(*geoRock);
#endif

  sectors.push_back(sec_def_sector);

  if(do_interpolation)
    ProposalMu = new PROPOSAL::Propagator(PROPOSAL::MuMinusDef::Get(), sectors, geoCan,interpolation_def);
  else
    ProposalMu = new PROPOSAL::Propagator(PROPOSAL::MuMinusDef::Get(), sectors, geoCan);

#endif





  return;
}
//________________________________________________________________________________________
GSeaRealAtmoFlux * GGenerateEvent::GetFluxDiffuse(void)
{

  // Instantiate the diffuse flux driver

  GSeaRealAtmoFlux * atmo_flux_driver = new GSeaRealAtmoFlux(fGenPar);

  // set rotation for coordinate tranformation from the topocentric horizontal
  // system to a user-defined coordinate system:

  if(!fGenPar->Rotation.IsIdentity()) {
     atmo_flux_driver->SetUserCoordSystem(fGenPar->Rotation);
  }

  fGenPar->GlobalGenWeight=atmo_flux_driver->GetGlobalGenWeight();

  return atmo_flux_driver;
}
//________________________________________________________________________________________
GSeaRealPointFlux* GGenerateEvent::GetFluxPoint(void)
{
  // Instantiate the point flux driver

  GSeaRealPointFlux * point_flux_driver = new GSeaRealPointFlux(fGenPar);

  // set rotation for coordinate tranformation from the topocentric horizontal
  // system to a user-defined coordinate system:
  if(!fGenPar->Rotation.IsIdentity()) {
     point_flux_driver->SetUserCoordSystem(fGenPar->Rotation);
  }

  fGenPar->GlobalGenWeight=point_flux_driver->GetGlobalGenWeight();

  return point_flux_driver;
}
//________________________________________________________________________________________
int GGenerateEvent::PosInCan(double VertX, double VertY, double VertZ){
  // check if the interaction vertex is inside the detector can

  int VerInCan = 0;

  if (TMath::Sqrt(VertX*VertX+VertY*VertY)<=fGenPar->Can3) {
    if(VertZ>= fGenPar->Can1 && VertZ<=fGenPar->Can2)VerInCan = 1;
    else if(VertZ>= fGenPar->Can1-fGenPar->HBedRock && VertZ<fGenPar->Can1)VerInCan = 2;

  }

  return VerInCan;
}
//________________________________________________________________________________________
void GGenerateEvent::FillEvent()
{

  // just copy the thing.
  fSeaEvent->Tracks = gWritenTracks;
  fSeaEvent->NTracks=fSeaEvent->Tracks.size();

}
//________________________________________________________________________________________
bool GGenerateEvent::PropagationMu(GSeaTrack *t)
{
  // wrapper function to call propagation with regular GSeaTrack

  PropaTrack pt;
  copy( *t, pt );
  bool r = PropagationMu ( &pt );
  copy( pt, *t );
  return r;
}
//________________________________________________________________________________________
bool GGenerateEvent::PropagationMu(PropaTrack * trackf){

  if(trackf->E<0.16)return false; // muon at rest (no propagation)

  PropaTrack track=*trackf;

  int iFlagCan;
  bool LepCan=false;
  double DepthW,DepthR;

  double range=pow(10.,fGenPar->RangeSW.Eval(log10(track.E)))/(fGenPar->RhoSW);

#ifdef _MUSIC_ENABLED__
  int MUMODEL2=1;
  double DepthsF;
  double RhoSW=fGenPar->RhoSW;
  double RhoSR=fGenPar->RhoSR;
#endif

// We want to know if the muon track intercepts the can

  iFlagCan=DistToCan(track);

  if(iFlagCan==1)LepCan = true;                         // the track is inside the can (no propagation)
  else if(iFlagCan==2 || iFlagCan==-1)LepCan = false;   // the track does not intercept the can (no propagation)
  else if(range<fDistaCan)LepCan = false;               // muon does not have chance to reach the can (no propagation)
  else {
#ifdef _JPP_ENABLED__
    double JPP_Emin=1; //GeV Minimal energy, if energy is below 1 GeV, only ionization is computed until stop.
#endif
    while(1){
#ifdef _JPP_ENABLED__
      JSIRENE::JVertex vertex(0, track.T, track.E); //The "Z" to propagate in jpp is the first component.
#endif
      bool InRock=false;
      DepthR=fabs((fGenPar->Can1-track.V3)/(track.D3));
      DepthW=fDistaCan;

      if(track.V3<=fGenPar->Can1){
        if(DepthR>0.01)InRock=true;
        else InRock=false;
      }

      if(InRock){                       // I'm in rock

        // Transport in rock the muon
        if(fGenPar->PropCode.compare("PropaMuon")==0){
          if(!PropaMuonSR->Propagate(track,DepthR)){LepCan=false; break;}
           *trackf = PropaMuonSR->GetPropaTrack();
        }
#ifdef _MUSIC_ENABLED__
        else if(fGenPar->PropCode.compare("MUSIC")==0){
          DepthR*=100.;
          track.V1*=100.;
          track.V2*=100.;
          track.V3*=100.;
          muon_transport_sr_(&track.V1,&track.V2,&track.V3,&track.D1,&track.D2,&track.D3,&track.E,&DepthR,&track.T,&RhoSR,&MUMODEL2,&MUMODEL2,&trackf->V1,&trackf->V2,&trackf->V3,&trackf->D1,&trackf->D2,&trackf->D3,&trackf->E,&DepthsF,&trackf->T);
          trackf->V1/=100.;
          trackf->V2/=100.;
          trackf->V3/=100.;
        }
#endif
#ifdef _JPP_ENABLED__
        else if(fGenPar->PropCode.compare("JPP")==0){
          while (vertex.getE() >= JPP_Emin && vertex.getZ() < DepthR) {
            const int Nr = radiationSR.size();
            double li[Nr];                                                 // inverse interaction lengths
            double ls = 1.0e-5;                                           // minimal total inverse interaction length (100 km)^-1

            for (int i = 0; i != Nr; ++i) {
              ls += li[i] = radiationSR[i]->getInverseInteractionLength(vertex.getE());
            }

            const int Ni = ionizationSR.size();

            double Ai[Ni];
            double As=0;

            for (int i = 0; i != Ni; ++i) {
              As +=Ai[i]= ionizationSR[i]->getA(vertex.getE());
            }
            const double ds = min(DepthR-vertex.getZ(),min(gRandom->Exp(1.0) / ls, vertex.getRange(As)));
            vertex.step(ds,As); //Ionization continuous losses
            if (vertex.getE() >= JPP_Emin) {
              double Es = 0.0;                                            // shower energy [GeV]
              double y  = gRandom->Uniform(ls);
              for (int i = 0; i != Nr; ++i) {
                y -= li[i];
                if (y < 0.0) {
                  Es = radiationSR[i]->getEnergyOfShower(vertex.getE());    // shower energy [GeV]
                  break;
                }
              }
              vertex.applyEloss(Es);
            }
          }
          // add muon end point
          const int Ni = ionizationSR.size();

          double Ai[Ni];
          double As=0;

          for (int i = 0; i != Ni; ++i) {
            As +=Ai[i]= ionizationSR[i]->getA(vertex.getE());
          }
          if (vertex.getE() < JPP_Emin && vertex.getRange(As) > 0.0 && vertex.getZ()  < DepthR) {
            vertex.step(min(vertex.getRange(),DepthR-vertex.getZ()),As);
          }
          trackf->E=vertex.getE();
          //Projecting the distance traveled into the track direction.
          trackf->V1+=vertex.getZ()*track.D1;
          trackf->V2+=vertex.getZ()*track.D2;
          trackf->V3+=vertex.getZ()*track.D3;
          trackf->T=vertex.getT();

        }
#endif
#ifdef _PROPOSAL_ENABLED__
        else if(fGenPar->PropCode.compare("PROPOSAL")==0){
	
	  PROPOSAL::Vector3D direction(track.D1,track.D2,track.D3);
	  direction.CalculateSphericalCoordinates();

#if _PROPOSAL_VERSION_MAJOR__ > 5
          PROPOSAL::DynamicData particle_mu(PROPOSAL::MuMinusDef::Get().particle_type);

          particle_mu.SetEnergy(track.E*1E3); // [MeV]
          particle_mu.SetPropagatedDistance(0);
          particle_mu.SetPosition(PROPOSAL::Vector3D(track.V1*1E2,track.V2*1E2,track.V3*1E2));
          particle_mu.SetDirection(direction);

          particle_mu.SetTime(track.T); // in ?

          PROPOSAL::Secondaries sec = ProposalMu->Propagate(particle_mu,DepthW*100.);

          trackf->E=sec.GetEnergy().back()/1.E3;
          trackf->V1=sec.GetPosition().back().GetX()/1.E2;
          trackf->V2=sec.GetPosition().back().GetY()/1.E2;
          trackf->V3=sec.GetPosition().back().GetZ()/1.E2;

          trackf->D1=sec.GetDirection().back().GetX();
          trackf->D2=sec.GetDirection().back().GetY();
          trackf->D3=sec.GetDirection().back().GetZ();

          trackf->T=sec.GetTime().back();
#else

          ProposalMu->GetParticle().SetEnergy(track.E*1E3);     // in Mev

          ProposalMu->GetParticle().SetDirection(direction);
          ProposalMu->GetParticle().SetPosition(PROPOSAL::Vector3D(track.V1*1E2,track.V2*1E2,track.V3*1E2));
          ProposalMu->GetParticle().SetTime(track.T);   // in ?

          ProposalMu->GetParticle().SetPropagatedDistance(0);
          ProposalMu->Propagate(DepthW*100.);

          trackf->E=ProposalMu->GetParticle().GetEnergy()/1.E3;
          trackf->V1=ProposalMu->GetParticle().GetPosition().GetX()/1.E2;
          trackf->V2=ProposalMu->GetParticle().GetPosition().GetY()/1.E2;
          trackf->V3=ProposalMu->GetParticle().GetPosition().GetZ()/1.E2;

          trackf->D1=ProposalMu->GetParticle().GetDirection().GetX();
          trackf->D2=ProposalMu->GetParticle().GetDirection().GetY();
          trackf->D3=ProposalMu->GetParticle().GetDirection().GetZ();

          trackf->T=ProposalMu->GetParticle().GetTime();
#endif

        }
#endif

	trackf->DirectionCosines();

      }
      else{

        // Transport in sea water the muon
        if(fGenPar->PropCode.compare("PropaMuon")==0){
          if(!PropaMuonSW->Propagate(track,DepthW)){LepCan=false; break;}
           *trackf = PropaMuonSW->GetPropaTrack();
        }
#ifdef _MUSIC_ENABLED__
        else if(fGenPar->PropCode.compare("MUSIC")==0){
          DepthW*=100.;
          track.V1*=100.;
          track.V2*=100.;
          track.V3*=100.;
          muon_transport_sw_(&track.V1,&track.V2,&track.V3,&track.D1,&track.D2,&track.D3,&track.E,&DepthW,&track.T,&RhoSW,&MUMODEL2,&MUMODEL2,&trackf->V1,&trackf->V2,&trackf->V3,&trackf->D1,&trackf->D2,&trackf->D3,&trackf->E,&DepthsF,&trackf->T);
          trackf->V1/=100.;
          trackf->V2/=100.;
          trackf->V3/=100.;
        }
#endif
#ifdef _JPP_ENABLED__
        else if(fGenPar->PropCode.compare("JPP")==0){
          while (vertex.getE() >= JPP_Emin && vertex.getZ()  < DepthW) {
            const int Nr = radiationSW.size();
            double li[Nr];                                                 // inverse interaction lengths
            double ls = 1.0e-5;                                           // minimal total inverse interaction length (100 km)^-1

            for (int i = 0; i != Nr; ++i) {
              ls += li[i] = radiationSW[i]->getInverseInteractionLength(vertex.getE());
            }

            const int Ni = ionizationSW.size();

            double Ai[Ni];
            double As=0;

            for (int i = 0; i != Ni; ++i) {
              As +=Ai[i]= ionizationSW[i]->getA(vertex.getE());
            }
            const double ds = min(DepthW-vertex.getZ(),min(gRandom->Exp(1.0) / ls, vertex.getRange(As)));
            vertex.step(ds,As); //Ionization continuous losses
            if (vertex.getE() >= JPP_Emin) {

              double Es = 0.0;                                            // shower energy [GeV]
              double y  = gRandom->Uniform(ls);

              for (int i = 0; i != Nr; ++i) {

                y -= li[i];
                if (y < 0.0) {
                  Es = radiationSW[i]->getEnergyOfShower(vertex.getE());    // shower energy [GeV]
                  break;
                }
              }

              vertex.applyEloss(Es);
            }
          }
          // add muon end point
          const int Ni = ionizationSW.size();

          double Ai[Ni];
          double As=0;

          for (int i = 0; i != Ni; ++i) {
            As +=Ai[i]= ionizationSW[i]->getA(vertex.getE());
          }
          if (vertex.getE() < JPP_Emin && vertex.getRange(As) > 0.0 && vertex.getZ()< DepthW) {
            vertex.step(min(DepthW-vertex.getZ(),vertex.getRange()),As);
          }
          trackf->E=vertex.getE();
          //Projecting the distance traveled into the track direction.
          trackf->V1+=vertex.getZ()*track.D1;
          trackf->V2+=vertex.getZ()*track.D2;
          trackf->V3+=vertex.getZ()*track.D3;
          trackf->T=vertex.getT();
        }
#endif
#ifdef _PROPOSAL_ENABLED__
        else if(fGenPar->PropCode.compare("PROPOSAL")==0){
	  
	  PROPOSAL::Vector3D direction(track.D1,track.D2,track.D3);
	  direction.CalculateSphericalCoordinates();

#if _PROPOSAL_VERSION_MAJOR__ > 5
          PROPOSAL::DynamicData particle_mu(PROPOSAL::MuMinusDef::Get().particle_type);

          particle_mu.SetEnergy(track.E*1E3); // [MeV]
          particle_mu.SetPropagatedDistance(0);
          particle_mu.SetPosition(PROPOSAL::Vector3D(track.V1*1E2,track.V2*1E2,track.V3*1E2));
          particle_mu.SetDirection(direction);

          particle_mu.SetTime(track.T); // in ?

          PROPOSAL::Secondaries sec = ProposalMu->Propagate(particle_mu,DepthW*100.);

          trackf->E=sec.GetEnergy().back()/1.E3;
          trackf->V1=sec.GetPosition().back().GetX()/1.E2;
          trackf->V2=sec.GetPosition().back().GetY()/1.E2;
          trackf->V3=sec.GetPosition().back().GetZ()/1.E2;

          trackf->D1=sec.GetDirection().back().GetX();
          trackf->D2=sec.GetDirection().back().GetY();
          trackf->D3=sec.GetDirection().back().GetZ();

          trackf->T=sec.GetTime().back();
#else

          ProposalMu->GetParticle().SetEnergy(track.E*1E3);     // in Mev

          ProposalMu->GetParticle().SetDirection(direction);
          ProposalMu->GetParticle().SetPosition(PROPOSAL::Vector3D(track.V1*1E2,track.V2*1E2,track.V3*1E2));
          ProposalMu->GetParticle().SetTime(track.T);   // in ?

          ProposalMu->GetParticle().SetPropagatedDistance(0);
          ProposalMu->Propagate(DepthW*100.);

          trackf->E=ProposalMu->GetParticle().GetEnergy()/1.E3;
          trackf->V1=ProposalMu->GetParticle().GetPosition().GetX()/1.E2;
          trackf->V2=ProposalMu->GetParticle().GetPosition().GetY()/1.E2;
          trackf->V3=ProposalMu->GetParticle().GetPosition().GetZ()/1.E2;

          trackf->D1=ProposalMu->GetParticle().GetDirection().GetX();
          trackf->D2=ProposalMu->GetParticle().GetDirection().GetY();
          trackf->D3=ProposalMu->GetParticle().GetDirection().GetZ();

          trackf->T=ProposalMu->GetParticle().GetTime();
#endif

        }
#endif

	trackf->DirectionCosines();
      }

      if(trackf->E<0.16){LepCan=false; break;}
      track=*trackf;

      iFlagCan=DistToCan(track);
      if(fDistaCan>=0 && fDistaCan<0.01){LepCan=true; break;}
      if(iFlagCan==1){LepCan=true; break;} // the track is inside the can (no propagation)
      else if(iFlagCan==2 || iFlagCan==-1){LepCan=false; break;}

    }  // close while(1)

  }
  if(LepCan) trackf->Status=1;

  return LepCan;

  }
//________________________________________________________________________________________
#ifdef _HETAU_ENABLED__
double GGenerateEvent::PropagationTau(GSeaTrack *t)
{
  // wrapper function to call propagation with regular GSeaTrack

  PropaTrack pt;
  copy( *t, pt );
  double r = PropagationTau ( &pt );
  copy( pt, *t );
  return r;
}
double GGenerateEvent::PropagationTau(PropaTrack * trackf){

  if(trackf->E<1.8) return -1; // tau at rest (no propagation)

  int iFlagCan = DistToCan(*trackf); // We want to know if the tau track intercepts the can
  assert(iFlagCan!=1); //tau should never be inside the can because we check that before

  RandomGen * rnd = RandomGen::Instance();

  if (iFlagCan==2 || iFlagCan==-1) return -1; // the track does not intercept the can (no propagation)
  else { //the track is outside the can and pointing

    double tauti = -TMath::Log( rnd->RndGen().Rndm() ) * lifetime*1e9; //ns for tausic

    PropaTrack track=*trackf;
    double RhoSW=fGenPar->RhoSW;
    double RhoSR=fGenPar->RhoSR;
    double DepthW,DepthR;
    double DepthsF;
    int idec;
    double tautf = -1.;
    bool InRock=false;

    while(1){
      DepthR=fabs((fGenPar->Can1-track.V3)/(track.D3));
      if(track.V3<=fGenPar->Can1) InRock = (DepthR>0.01) ? true : false;
      if (tautf>0.) tauti = tautf;
      if(InRock){     // I'm in rock
        DepthR*=100.;
        track.V1*=100.; track.V2*=100.; track.V3*=100.;
        tau_transport_sr_(&track.V1,&track.V2,&track.V3,&track.D1,&track.D2,&track.D3,&track.E,&DepthR,&track.T,&RhoSR,&TAUMODEL,&TAUMODEL,&trackf->V1,&trackf->V2,&trackf->V3,&trackf->D1,&trackf->D2,&trackf->D3,&trackf->E,&DepthsF,&trackf->T,&idec,&ITFLAG,&tauti,&tautf);
        trackf->V1/=100.; trackf->V2/=100.; trackf->V3/=100.;
      }
      else{
        DepthW=fDistaCan*100.;
        track.V1*=100.; track.V2*=100.; track.V3*=100.;
        tau_transport_sw_(&track.V1,&track.V2,&track.V3,&track.D1,&track.D2,&track.D3,&track.E,&DepthW,&track.T,&RhoSW,&TAUMODEL,&TAUMODEL,&trackf->V1,&trackf->V2,&trackf->V3,&trackf->D1,&trackf->D2,&trackf->D3,&trackf->E,&DepthsF,&trackf->T,&idec,&ITFLAG,&tauti,&tautf);
        trackf->V1/=100.; trackf->V2/=100.; trackf->V3/=100.;
      }
      track=*trackf;
      iFlagCan = DistToCan(track);
      if (track.E<1.8)                   return 0.;    //tau at rest before reaching the can
      if (idec==1)                       return 0.;    //tau has decayed before reaching the can
      if (fDistaCan>0 && fDistaCan<0.01) { trackf->Status=1; return tautf; } //tau has reached the can before decaying
      if (iFlagCan==1)                   { trackf->Status=1; return tautf; } //tau has reached the can before decaying
      if (iFlagCan==2 || iFlagCan==-1)   return -1;    //tau not pointing at detector
    }  // close while(1)

  }

  return -1;

}
//________________________________________________________________________________________
vector<GSeaTrack> GGenerateEvent::DecayTau(GSeaTrack tau, double tauti) {

    vector<GSeaTrack> decayed;

    GSeaTrack propagated_tau = tau;
    if (tauti>0) {
      double RhoSW=fGenPar->RhoSW;
      double Depth=1e10;
      double DepthsF;
      int idec;
      double tautf = -1.;
      tau.V1*=100.; tau.V2*=100.; tau.V3*=100.;
      tau_transport_sw_(&tau.V1,&tau.V2,&tau.V3,&tau.D1,&tau.D2,&tau.D3,&tau.E,&Depth,&tau.T,&RhoSW,&TAUMODEL,&TAUMODEL,&propagated_tau.V1,&propagated_tau.V2,&propagated_tau.V3,&propagated_tau.D1,&propagated_tau.D2,&propagated_tau.D3,&propagated_tau.E,&DepthsF,&propagated_tau.T,&idec,&ITFLAG,&tauti,&tautf);
      propagated_tau.V1/=100.; propagated_tau.V2/=100.; propagated_tau.V3/=100.;
    }

    //now we decay the propagated tau
    double momf = TMath::Sqrt( propagated_tau.E*propagated_tau.E - kTauMass2 );
    TauolaHEPEVTEvent * Tauola_evt = new TauolaHEPEVTEvent();
    TauolaHEPEVTParticle *Tauola_tau = new TauolaHEPEVTParticle( propagated_tau.PdgCode, 1, propagated_tau.D1*momf, propagated_tau.D2*momf, propagated_tau.D3*momf, propagated_tau.E, kTauMass, -1, -1, -1, -1 );
    Tauola_evt->addParticle(Tauola_tau);
    Tauola::decayOne(Tauola_tau,true,0.,0.,1);
    for ( int sec=1; sec<Tauola_evt->getParticleCount(); sec++ ) {
      GSeaTrack DecTrack;
      DecTrack.Status   = Tauola_evt->getParticle(sec)->getStatus();
      DecTrack.PdgCode  = Tauola_evt->getParticle(sec)->getPdgID();
      DecTrack.E        = Tauola_evt->getParticle(sec)->getE();
      double mom = TMath::Sqrt(TMath::Power(Tauola_evt->getParticle(sec)->getPx(),2)+TMath::Power(Tauola_evt->getParticle(sec)->getPy(),2)+TMath::Power(Tauola_evt->getParticle(sec)->getPz(),2));
      DecTrack.D1       = Tauola_evt->getParticle(sec)->getPx()/mom;
      DecTrack.D2       = Tauola_evt->getParticle(sec)->getPy()/mom;
      DecTrack.D3       = Tauola_evt->getParticle(sec)->getPz()/mom;
      DecTrack.V1       = propagated_tau.V1;
      DecTrack.V2       = propagated_tau.V2;
      DecTrack.V3       = propagated_tau.V3;
      DecTrack.T        = propagated_tau.T;
      DecTrack.MotherId = propagated_tau.Id;
      DecTrack.Length   = 0.;
      if (this->PosInCan(DecTrack.V1,DecTrack.V2,DecTrack.V3)==0) DecTrack.Status *= -1;
      decayed.push_back(DecTrack);
    }

    return decayed;

}
#endif
//___________________________________________________________________________
int GGenerateEvent::DistToCan(PropaTrack track)
{

// iFlag = -1 something has gone wrong!
// iFlag =  0 if the track intercepts the can AND the range is positive
// iFlag =  1 if the point is inside the can
// iFlag =  2 if the track never intercepts the can

  int iFlag=-1;
  fDistaCan = -1;

  double Zmin= fGenPar->Can1;
  double Zmax= fGenPar->Can2;
  double Radius = fGenPar->Can3;

  double Vx = track.V1;
  double Vy = track.V2;
  double Vz = track.V3;

  double Dx = track.D1;
  double Dy = track.D2;
  double Dz = track.D3;

  double utit=sqrt(Dx*Dx+Dy*Dy+Dz*Dz);

  Dx/=utit;
  Dy/=utit;
  Dz/=utit;

  fDistaX=0.; fDistaY=0.; fDistaZ=0.;

  if(Vz<=Zmax && Vz>=Zmin && Vx*Vx+Vy*Vy<=Radius*Radius){
    iFlag = 1;
    fDistaCan = 0.;
  }
  else if(fabs(Dz)==1){
    if (sqrt(Vx*Vx+Vy*Vy) > Radius) iFlag = 2;
    else if(Dz==1 && Vz>Zmax) iFlag = 2;
    else if(Dz==-1 && Vz<Zmin) iFlag = 2;
    else{
      iFlag =  0;
      if(Dz==1){
        fDistaX=Vx;
        fDistaY=Vy;
        fDistaZ=Zmin;
        fDistaCan = fabs(Zmin-Vz);
      }
      else {
        fDistaX=Vx;
        fDistaY=Vy;
        fDistaZ=Zmax;
        fDistaCan = fabs(Vz-Zmax);
      }
    }
  }
  else {
    double a = Dx*Dx+Dy*Dy;
    double b = 2.*(Vx*Dx+Vy*Dy);
    double c = Vx*Vx+Vy*Vy-Radius*Radius;
    double delta=b*b-4.*a*c;
    if(delta<0)iFlag = 2;
    else {
      delta = sqrt(delta);
      double t1=(-b-delta)/(2.*a);
      double t2=(-b+delta)/(2.*a);
      double z1=Vz+Dz*t1;
      double z2=Vz+Dz*t2;
      if((z1<Zmin && z2<Zmin) || (z1>Zmax && z2>Zmax)) iFlag = 2;
      else if(z1>=Zmin && z1<=Zmax){
        if(t1>=0){
          iFlag = 0;
          fDistaX=Vx+Dx*t1;
          fDistaY=Vy+Dy*t1;
          fDistaZ=Vz+Dz*t1;
          fDistaCan = fabs(t1);
        } else iFlag = 2;       // intersection before the vertex
      }
      else if(z1<Zmin && Dz>0){  // up-going
        iFlag = 0;
        t1=(Zmin-Vz)/Dz;
        if(t1>=0){
          fDistaX=Vx+Dx*t1;
          fDistaY=Vy+Dy*t1;
          fDistaZ=Vz+Dz*t1;
          fDistaCan=fabs((Zmin-Vz)/Dz);
        } else iFlag = 2;       // intersection before the vertex
      }
      else if(z1>Zmax && Dz<0){  // down-going
        iFlag = 0;
        t1=(Zmax-Vz)/Dz;
        if(t1>=0){
          fDistaX=Vx+Dx*t1;
          fDistaY=Vy+Dy*t1;
          fDistaZ=Vz+Dz*t1;
          fDistaCan=fabs((Zmax-Vz)/Dz);
        } else iFlag = 2;       // intersection before the vertex
      }
    }
  }

  return iFlag;
}
//________________________________________________________________________________________________________
bool GGenerateEvent::PropagateTracks(vector<GSeaTrack> GeneratorTracks) {

//Move all the muons to the can and add them to gWriteTracks

  gPropagatedTracks.clear();

  RandomGen * rnd = RandomGen::Instance();

  int nGenTracks = GeneratorTracks.size();

  //first propagate the taus (we need to propagate all of them)
#ifdef _HETAU_ENABLED__
  for ( int i=0; i<nGenTracks; i++ ) {
    if( abs(GeneratorTracks[i].PdgCode)!=15 ) continue;
    vector<GSeaTrack> decayed;
    if( GeneratorTracks[i].Status==1 ) { //propagate all final state taus inside can
      decayed = DecayTau(GeneratorTracks[i],-TMath::Log( rnd->RndGen().Rndm() )*lifetime*1e9);
      GeneratorTracks[i].Length = TMath::Sqrt(TMath::Power(decayed[0].V1-GeneratorTracks[i].V1,2)+TMath::Power(decayed[0].V2-GeneratorTracks[i].V2,2)+TMath::Power(decayed[0].V3-GeneratorTracks[i].V3,2));
    }
    else if( GeneratorTracks[i].Status==-1) { //outside the can
      auto propagated_tau = GeneratorTracks[i];
      double tauti = PropagationTau(&propagated_tau);
      if ( tauti==0 ) { // decayed before got to the can
        GeneratorTracks[i].Status -= 2000;
        decayed = DecayTau(propagated_tau,tauti);
        GeneratorTracks[i].Length = TMath::Sqrt(TMath::Power(decayed[0].V1-GeneratorTracks[i].V1,2)+TMath::Power(decayed[0].V2-GeneratorTracks[i].V2,2)+TMath::Power(decayed[0].V3-GeneratorTracks[i].V3,2));
      }
      else if ( tauti>0  ) { // didnt decay before got to the can
        GeneratorTracks[i].Status -= 1000;
        GeneratorTracks[i].Length = TMath::Sqrt(TMath::Power(propagated_tau.V1-GeneratorTracks[i].V1,2)+TMath::Power(propagated_tau.V2-GeneratorTracks[i].V2,2)+TMath::Power(propagated_tau.V3-GeneratorTracks[i].V3,2));
        propagated_tau.MotherId = GeneratorTracks[i].Id;
        propagated_tau.Id       = GeneratorTracks.size(); // post-increment
        decayed = DecayTau(propagated_tau,tauti);
        propagated_tau.Length = TMath::Sqrt(TMath::Power(decayed[0].V1-propagated_tau.V1,2)+TMath::Power(decayed[0].V2-propagated_tau.V2,2)+TMath::Power(decayed[0].V3-propagated_tau.V3,2));
        GeneratorTracks.push_back( propagated_tau );
      }
    }
    for ( unsigned int k=0; k<decayed.size(); k++ ) {
      decayed[k].Id = GeneratorTracks.size();
      GeneratorTracks.push_back( decayed[k] );
    }
  }
  nGenTracks = GeneratorTracks.size();
#endif

  //now go to other particles
  bool LepInCan = false;
  for ( int i=0; i<nGenTracks; i++ ) {
    if( abs(GeneratorTracks[i].PdgCode)==13 && GeneratorTracks[i].Status==-1 ) { //propagate final state muons that are outside can
      auto propagated_mu = GeneratorTracks[i];
      if ( PropagationMu( &propagated_mu ) ){ //mu has reached the boundary of the can
        GeneratorTracks[i].Status -= 1000; //Status < 1000 means: propagated track that got to the can
        GeneratorTracks[i].Length = TMath::Sqrt(TMath::Power(propagated_mu.V1-GeneratorTracks[i].V1,2)+TMath::Power(propagated_mu.V2-GeneratorTracks[i].V2,2)+TMath::Power(propagated_mu.V3-GeneratorTracks[i].V3,2));
        propagated_mu.MotherId = GeneratorTracks[i].Id;
        propagated_mu.Length   = 0.;
        propagated_mu.Id       = GeneratorTracks.size(); // post-increment
        GeneratorTracks.push_back( propagated_mu );
        LepInCan = true;
      }
    }
    else if ( GeneratorTracks[i].Status==1 ) {
      LepInCan = true;
      if ( TMath::Abs(GeneratorTracks[i].PdgCode)==311 ) { // K0 (K0Bar) is detected as K0S or KOL
        GeneratorTracks[i].PdgCode = (rnd->RndFlux().Rndm()>0.5) ? 130 : 310;
      }
    }
  }

  gPropagatedTracks = GeneratorTracks;

  return LepInCan;

}
//________________________________________________________________________________________________________
bool GGenerateEvent::WriteTracks(bool LepInCan) {

//Add particles to the gWriteTracks depending on WriteParticlesMode.

  gWritenTracks.clear();

  if      ( fGenPar -> WriteParticlesMode == 2 ) { gWritenTracks = gPropagatedTracks; }
  else if ( fGenPar -> WriteParticlesMode == 1 ) { if (LepInCan) gWritenTracks = gPropagatedTracks; }
  else if ( fGenPar -> WriteParticlesMode == 0 ) {
    if (LepInCan) {
      for ( auto& track : gPropagatedTracks ) {
        if      ( track.Id==0        ) gWritenTracks.push_back(track); //we always store the primary particle (neutrino or cosmic ray)
        else if ( track.Status<-1000 ) gWritenTracks.push_back(track); //we always store the propagated particles
        else if ( track.Status==1    ) gWritenTracks.push_back(track); //we store the final state particles inside the can
      }
    }
  }

  return gWritenTracks.size() > 0;

}
//___________________________________________________________________________
double GGenerateEvent::GetXSecWater(int neutrino_pdg, double Energy){

  double XSec=0.;

  map<int,TGraph>::iterator iter = fGenPar->XSecWater.begin();
  for ( ; iter != fGenPar->XSecWater.end(); ++iter) {
    if(neutrino_pdg==iter->first){
       TGraph GrSecTot = iter->second;
       XSec=GrSecTot.Eval(Energy);
       break;
    }
  }

  return XSec;
}
//___________________________________________________________________________
double GGenerateEvent::GetXSecSW(int neutrino_pdg, double Energy){

  double XSec=0.;

  map<int,TGraph>::iterator iter = fGenPar->XSecTotSW.begin();
  for ( ; iter != fGenPar->XSecTotSW.end(); ++iter) {
    if(neutrino_pdg==iter->first){
       TGraph GrSecTot = iter->second;
       XSec=GrSecTot.Eval(Energy);
       break;
    }
  }

  return XSec;
}