//________________________________________________________________________________________
/*!

\program gSeaNuEvGen

\brief   A GENIE application to simulate neutrino events in neutrino telescopes

\created May 5, 2012

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

\cpright  Copyright (c) 2012-2021, The KM3NeT Collaboration


*/
//_________________________________________________________________________________________

#include <cassert>
#include <cstdlib>
#include <cctype>
#include <string>
#include <vector>
#include <sstream>
#include <map>
#include <cstdio>
#include <iostream>
#include <fstream>
#include <cmath>
#include <algorithm>

#include "libxml/parser.h"
#include "libxml/xmlmemory.h"
#include "libxml/xmlreader.h"

#include <boost/preprocessor/stringize.hpp>

#define STR_(X) #X
#define STR(X) STR_(X)


// GENIE

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


// gSeaGen

#include "SeaNuDrivers/GGenerateEventGenie.h"
#include "SeaNuDrivers/GSeaGeometry.h"



#ifdef _ANTARES_ENABLED__
#include "OutputWriters/EvtWrite.h"
#endif
#ifdef _KM3NET_ENABLED__
#include "OutputWriters/KM3NeTWrite.h"
#endif
#include "OutputWriters/SeaNtpWriter.h"

#include "SeaEvent/GBinParam.h"


using namespace std;

using namespace genie;
using namespace genie::geometry;

void                    PrintName          (void);
void                    PrintSyntax        (void);

void                    CreaBin            (void);
void                    PrintStat          (void);

void                    SetDefParams       (void);
void                    LoadDefParams      (void);

#ifdef _ANTARES_ENABLED__
bool                    LoadGeometryDet    (string filename);
bool                    LoadGeometryDetx   (string filemane);
#endif

void                    CalcUTMConvergeAngle(void);

void                    LoadMediaComp      (void);

void                    LoadAstroInfo      (void);

void                    DecodeCommandLine  (int argc, char * argv[]);



// User-specified options:

bool             gOptPDG       = false;
string           gOptTreeFormat= "SimpleTree";

vector<BinParam> VecBinPar;
GenParam         gGenPar;

tm *gDate;


//________________________________________________________________________________________
int main(int argc, char** argv)
{

  // get gSeaGen version
  string VerFile = gSystem->Getenv("GSEAGEN");
  VerFile.append("/VERSION");
  ifstream file_stream(VerFile.c_str(), ios::in);
  file_stream>>gGenPar.gSeaGenVer;
  file_stream.close();

  // get GENIE version
  gGenPar.GenieVer=BOOST_PP_STRINGIZE( __GENIE_RELEASE__);
  while(1){
    size_t pos = gGenPar.GenieVer.find("\"");
    if(pos==string::npos)break;
    gGenPar.GenieVer.erase(pos,1);
  }

  PrintName();
  if(argc==1){
    PrintSyntax();
    exit(0);
  }

  if(!gSystem->Getenv("GSEAGEN")){
    LOG("gSeaNuEvGen", pFATAL) << "env. parameter GSEAGEN does not exist; exiting...";
    gAbortingInErr = true;
    exit(1);
  }

  // Before parsing command line arguments, we set default parameters
  SetDefParams();
  LoadDefParams();

  // Load and re-define media composition if request by the user
  LoadMediaComp();

  // Parse command line arguments
  DecodeCommandLine(argc,argv);

  // Calculate the UTM convergence angle
  CalcUTMConvergeAngle();

  // Create the bins and their geometries
  CreaBin();

  // Define the generator
  GGenerateEventGenie * gGenEvtGenie=0;

  GGenerateEvent * gGenEvt=0;

  if(gGenPar.GenCode.compare("GENIE")==0){
    gGenEvtGenie = new GGenerateEventGenie(&gGenPar);
    gGenEvt = dynamic_cast <GGenerateEvent*> (gGenEvtGenie);
  }

#ifdef _ANTARES_ENABLED__
  // initialize an evt writer (to be done after the instance of the flux driver...)
  EvtWrite * evtwrt = 0;
  if(gGenPar.GenMode=="EVT") evtwrt = new EvtWrite(&gGenPar,gOptPDG,gGenEvt->fAntaresEvt);
  else                       evtwrt = new EvtWrite(&gGenPar,gOptPDG);
  evtwrt->WriteHeader();
  SLOG("gSeaNuEvGen", pINFO) << "evt writer initialized";
#endif
#ifdef _KM3NET_ENABLED__
  // initialize an aaevt writer (to be done after the instance of the flux driver...)
  KM3NeTWrite * km3net_writer = new KM3NeTWrite( &gGenPar );
  km3net_writer->WriteHeader();
  SLOG("gSeaNuEvGen", pINFO) << "km3net writer initialized";
#endif
  SeaNtpWriter * ntprt = 0;
  if(gGenPar.GenMode=="DIFFUSE" || gGenPar.GenMode=="POINT" ) ntprt = new SeaNtpWriter(&gGenPar, gOptTreeFormat);
  for(int iBin=0; iBin<gGenPar.NBin; iBin++){
    VecBinPar[iBin].NEvBinGen=0;
    gGenEvt->InitGenBin(VecBinPar[iBin].EneMin*units::GeV,VecBinPar[iBin].EneMax*units::GeV,VecBinPar[iBin].RL,VecBinPar[iBin].RT,VecBinPar[iBin].X0,VecBinPar[iBin].Y0,VecBinPar[iBin].Z0);
    if(iBin==0 || gGenPar.TrackEvts>0 ) gGenEvt->SetGeometry(VecBinPar[iBin].geometry);
    gGenEvt->Configure(VecBinPar[iBin].NEvBinScal);
    VecBinPar[iBin].ProbScale = gGenEvt->GlobProbScale();

    // generate next event at the can
    while (1){
      gGenEvt->GenerateEvent();
      if(gGenEvt->GetNFluxNeutrinos()>VecBinPar[iBin].NEvBinScal || gGenEvt->EoF) break;
#ifdef _ANTARES_ENABLED__
      evtwrt->WriteEvent(gGenEvt->fSeaEvent);
#endif
#ifdef _KM3NET_ENABLED__
      km3net_writer->WriteEvent(gGenEvt->fSeaEvent);
#endif
      if(gGenPar.GenMode=="POINT" || gGenPar.GenMode=="DIFFUSE") ntprt->WriteEvent(gGenEvt->fSeaEvent);
      gGenEvt->WriteNative(gGenEvt->fSeaEvent);
      VecBinPar[iBin].NEvBinGen++;
      if( !(gGenEvt->fNEvt%500) ) LOG("gSeaNuEvGen", pNOTICE) << "Generated events "<<gGenEvt->fNEvt;
      // clean-up
      gGenEvt->CleanEvent();

    }

    PrintStat();
    LOG("gSeaNuEvGen", pNOTICE) << "Done bin: " << iBin+1;
  }

  // clean-up

  if(gGenPar.GenCode.compare("GENIE")==0){
    if(gGenEvtGenie)delete gGenEvtGenie;
  }
#ifdef _ANTARES_ENABLED__
  if(evtwrt)delete evtwrt;
#endif

#ifdef _KM3NET_ENABLED__
  if(km3net_writer)delete km3net_writer;
#endif

  if(gGenPar.GenMode=="POINT" || gGenPar.GenMode=="DIFFUSE")delete ntprt;

  system("rm geometry-*.root");
#ifdef _MUSIC_ENABLED__
  if(gGenPar.TrackEvts>0 && gGenPar.PropCode.compare("MUSIC")==0)system("rm music*.dat");

#endif

#ifdef _HETAU_ENABLED__
  system("rm tausic*.dat");
#endif

  return 0;
}
//________________________________________________________________________________________
void SetDefParams(void){
// Set default parameters, this function must be called before the GetCommandLineArgs or at beginning of it

  //running time
  time(&gGenPar.RunTime);
  gDate = localtime (&gGenPar.RunTime);

  gGenPar.RunNu = 100000000;

  gGenPar.NTot = 0;
  gGenPar.Rotation.SetToIdentity();

  gGenPar.Can1 = 0.000;
  gGenPar.Can2 = 659.862;
  gGenPar.Can3 = 309.450;


  gGenPar.OriginX = 0.;
  gGenPar.OriginY = 0.;
  gGenPar.OriginZ = 0.;

  gGenPar.WriteParticlesMode = 0;
  gGenPar.SaveMu = true;
  gGenPar.SaveNu = false;
  gGenPar.ChancesPerShower = 0;
  for (int i=0; i<3; i++) gGenPar.CustomRotation[i] = 0.;
  gGenPar.SaveOther = false;
  gGenPar.WeightOpt = 0;

  gGenPar.Alpha=0.;
  gGenPar.Primary=-1;
  gGenPar.CRModel="GST3";
  gGenPar.EmuMin=0.;
  gGenPar.EPrimMin=0.;
  gGenPar.EPrimMax=0.;
  gGenPar.PrimCtMin=0.;
  gGenPar.PrimCtMax=0.;

  gGenPar.RhoSW = 1.03975;              // g/cm3
  gGenPar.RhoSR = 2.65;                 // g/cm3

  gGenPar.NAbsLength = 3.;
  gGenPar.AbsLength = 70.;              // m
  gGenPar.SiteDepth = 2475.;            // m
  gGenPar.HBedRock = 0.;                // m
  gGenPar.SiteLatitude = 0.74700;       // in rad
  gGenPar.SiteLongitude = 0.10763;      // in rad
  gGenPar.SiteDeclination = 2.28*kPi/180.; // in rad ; default value (declination of ORCA site on 12.01.2021 according to WMM2020)
  gGenPar.UTMConvergeAngle = 0.;        // in rad, calculate at run-time
  gGenPar.UseUTMSystem = false;
  gGenPar.SeaBottomRadius=6368.E3;      // m (Earth radius at the sea bottom level)


  gGenPar.NBin = 1;
  gGenPar.Alpha = 1.4;
  gGenPar.EvMin = 1.;
  gGenPar.EvMax = 1000.;
  gGenPar.CtMin = 0.;
  gGenPar.CtMax = 1.;

  gGenPar.RanSeed=12345;

  gGenPar.PropMode=false;
  gGenPar.Drawing="surface";

  gGenPar.GenMode="DIFFUSE";
//  gGenPar.TGen = 31556926.;             // 1 yr in seconds
  gGenPar.TGen = 1.;             	// 1 second



  gGenPar.EvFilePrefix = "gseagen";     // event file prefix (override with -o)
  gGenPar.NativeOutput = false;
  gGenPar.PropCode  = "PropaMuon";
  gGenPar.PropCodeVer="";
  
  gGenPar.GenCode ="GENIE";

  gGenPar.RTOpt = "can";

  gGenPar.TrackEvts = 0;

  gGenPar.SourceName="";
  gGenPar.Declination=0.;
  gGenPar.RightAscension=0.;
  gGenPar.SourceRadius=0.;

  gGenPar.TimeStampStart = 0;
  gGenPar.TimeStampStop = 0;
  gGenPar.MJDStart=0.;
  gGenPar.MJDStop=0.;
  gGenPar.GenMJD=false;

  gGenPar.InpXSecFile=("gxspl-FNALsmall.xml");  // standard genie official xsec file name


  gGenPar.EventGeneratorList="Default";

  // media composition: SeaWater, Rock (Crust), Mantle, Core
  // SeaWater composition
  gGenPar.SeaWaterComp.insert(map<int, double>::value_type(1000080160,0.8584));         // O16
  gGenPar.SeaWaterComp.insert(map<int, double>::value_type(1000010010,0.1082));         // H1
  gGenPar.SeaWaterComp.insert(map<int, double>::value_type(1000170350,1.94E-2));        // Cl35
  gGenPar.SeaWaterComp.insert(map<int, double>::value_type(1000110230,1.08E-2));        // Na23
  gGenPar.SeaWaterComp.insert(map<int, double>::value_type(1000120240,0.1292E-2));      // Mg24
  gGenPar.SeaWaterComp.insert(map<int, double>::value_type(1000160320,0.091E-2));       // S32
  gGenPar.SeaWaterComp.insert(map<int, double>::value_type(1000200400,0.04E-2));        // Ca40
  gGenPar.SeaWaterComp.insert(map<int, double>::value_type(1000190390,0.04E-2));        // K39
  gGenPar.SeaWaterComp.insert(map<int, double>::value_type(1000350800,0.0067E-2));      // Br80
  gGenPar.SeaWaterComp.insert(map<int, double>::value_type(1000060120,0.0028E-2));      // C12

  // Rock composition: Crust
  gGenPar.RockComp.insert(map<int, double>::value_type(1000080160,0.463));              // O16
  gGenPar.RockComp.insert(map<int, double>::value_type(1000140280,0.282));              // Si28
  gGenPar.RockComp.insert(map<int, double>::value_type(1000130270,0.0823));             // Al27
  gGenPar.RockComp.insert(map<int, double>::value_type(1000260560,0.0563));             // Fe56
  gGenPar.RockComp.insert(map<int, double>::value_type(1000200400,0.0415));             // Ca40
  gGenPar.RockComp.insert(map<int, double>::value_type(1000110230,0.0236));             // Na23
  gGenPar.RockComp.insert(map<int, double>::value_type(1000120240,0.0233));             // Mg24
  gGenPar.RockComp.insert(map<int, double>::value_type(1000190390,0.0209));             // K39
  gGenPar.RockComp.insert(map<int, double>::value_type(1000220480,0.0057));             // Ti48
  gGenPar.RockComp.insert(map<int, double>::value_type(1000010010,0.0014));             // H1

  // Mantle composition
  gGenPar.MantleComp.insert(map<int, double>::value_type(1000080160,0.4522));           // O16
  gGenPar.MantleComp.insert(map<int, double>::value_type(1000120240,0.2283));           // Mg24
  gGenPar.MantleComp.insert(map<int, double>::value_type(1000140280,0.2149));           // Si28
  gGenPar.MantleComp.insert(map<int, double>::value_type(1000260560,0.0597));           // Fe56
  gGenPar.MantleComp.insert(map<int, double>::value_type(1000130270,0.0225));           // Al27
  gGenPar.MantleComp.insert(map<int, double>::value_type(1000200400,0.0224));           // Ca40

  // Core composition
  gGenPar.CoreComp.insert(map<int, double>::value_type(1000260560,0.9));                // Fe56
  gGenPar.CoreComp.insert(map<int, double>::value_type(1000280580,0.1));                // Ni58


  return;
}
//________________________________________________________________________________________
void LoadAstroInfo(void){

  xmlTextReaderPtr reader;

  string filename=gGenPar.SourceFile;
  string Equinox="J2000"; // default

  // astro file input, read info from file

  bool CheckDecl=false;
  bool CheckRA=false;
  bool CheckLat=false;
  bool CheckLong=false;

  LOG("gSeaNuEvGen", pINFO) <<"Reading source info from file "<< filename;
  reader = xmlNewTextReaderFilename(filename.c_str());

  string smedia,selem;

  //  string DateStart,DateStop;
  double GalacticLatitude=0.,GalacticLongitude=0.;

  while(xmlTextReaderRead(reader)==1){

    xmlChar * name  = xmlTextReaderName     (reader);
    xmlChar * value = xmlTextReaderValue    (reader);
    int       type  = xmlTextReaderNodeType (reader);
    int       depth = xmlTextReaderDepth    (reader);

    if(depth==0 && type==1){
      if(xmlStrcmp(name, (const xmlChar *) "astro_source")){
         LOG("gSeaNuEvGen", pFATAL) << " Invalid root element in filename " << filename <<"; exiting... ";
         exit(1);
       }
    }
    if(!xmlStrcmp(name, (const xmlChar*) "param_set" ) && type==1){
      xmlChar * xmedia   = xmlTextReaderGetAttribute(reader,(const xmlChar*)"source_name");
      smedia = utils::str::TrimSpaces((const char *)xmedia);
      gGenPar.SourceName=smedia;
      xmlFree(xmedia);
    }
    if(!xmlStrcmp(name, (const xmlChar *) "param" ) && type==1) {
      xmlChar * xelem   = xmlTextReaderGetAttribute(reader,(const xmlChar*)"name");
      selem = utils::str::TrimSpaces((const char *)xelem);
      xmlFree(xelem);
    }
    if(depth==3){
      if(selem=="Declination"){
        gGenPar.Declination=atof((const char *)value);
        CheckDecl=true;
      }
      else if(selem=="RightAscension"){
        gGenPar.RightAscension=atof((const char *)value);
        CheckRA=true;
      }
      else if(selem=="Radius"){
        gGenPar.SourceRadius=atof((const char *)value);
        gGenPar.SourceRadius*=kPi/180.;
      }
      else if(selem=="GalacticLatitude"){
        GalacticLatitude=atof((const char *)value);
        CheckLat=true;
      }
      else if(selem=="GalacticLongitude"){
        GalacticLongitude=atof((const char *)value);
        CheckLong=true;
      }
      else if(selem=="Equinox"){
        Equinox=(const char *)value;
        Equinox.erase(remove_if(Equinox.begin(), Equinox.end(), ::isspace), Equinox.end());  // erase blanc spaces
      }
    }
    xmlFree(name);
    xmlFree(value);
  }

  // check info validity and some calculations

  bool Equatorial=false;

  if(CheckDecl){
    if(fabs(gGenPar.Declination)>90.){
      LOG("gSeaNuEvGen", pFATAL) << "Source declination must be in the range [-90.,+90.] degrees, exiting...";
      gAbortingInErr = true;
      exit(1);
    }
  }
  if(CheckRA){
    if(gGenPar.RightAscension<0. and gGenPar.RightAscension>360.){
      LOG("gSeaNuEvGen", pFATAL) << "Source Right Ascension must be in the range [0,360] degrees, exiting...";
      gAbortingInErr = true;
      exit(1);
    }
  }
  if(CheckLat){
    if(fabs(GalacticLatitude)>90.){
      LOG("gSeaNuEvGen", pFATAL) << "Source Galactic Latitude must be in the range [-90.,+90.] degrees, exiting...";
      gAbortingInErr = true;
      exit(1);
    }
  }
  if(CheckLong){
    if(GalacticLongitude<0. and GalacticLongitude>360.){
      LOG("gSeaNuEvGen", pFATAL) << "Source Galactic Longitude must be in the range [0,360] degrees, exiting...";
      gAbortingInErr = true;
      exit(1);
    }
  }

  if((CheckDecl && !CheckRA) || (!CheckDecl && CheckRA)){
    LOG("gSeaNuEvGen", pFATAL) << "Both Horizontal Coordinates must be defined, exiting...";
    gAbortingInErr = true;
    exit(1);
  }
  else if(CheckDecl && CheckRA) {
    Equatorial=true;
    LOG("gSeaNuEvGen", pINFO) << "Point-source declination: "<<gGenPar.Declination;
    LOG("gSeaNuEvGen", pINFO) << "Point-source right ascension: "<<gGenPar.RightAscension;
    gGenPar.Declination*=kPi/180.;
    gGenPar.RightAscension*=kPi/180.;
  }
  else if((CheckLat && !CheckLong) || (!CheckLat && CheckLong)){
    LOG("gSeaNuEvGen", pFATAL) << "Both Galactic Coordinates must be defined, exiting...";
    gAbortingInErr = true;
    exit(1);
  }
  else {
     LOG("gSeaNuEvGen", pINFO) << "Point-source galactic latitude: "<<GalacticLatitude;
     LOG("gSeaNuEvGen", pINFO) << "Point-source galactic longitude: "<<GalacticLongitude;
     GalacticLatitude*=kPi/180.;
     GalacticLongitude*=kPi/180.;
  }


  if(Equinox.compare("J2000")!=0 && Equinox.compare("B1950")!=0){
    LOG("gSeaNuEvGen", pFATAL)<< "Unknown Equinox, it must be J2000 or B1950, exiting... "<<Equinox;
    gAbortingInErr = true;
    exit(1);
  }


  if(!Equatorial  || Equinox.compare("B1950")==0){
    // we have to convert coordinates to euqtorial at epoch J2000.0

    GAstro gAstro(gGenPar.SiteLatitude, gGenPar.SiteLongitude);

    if(Equatorial){
      gAstro.CalcPrecession(33281.924, 51544.5, gGenPar.RightAscension, gGenPar.Declination, &gGenPar.RightAscension, &gGenPar.Declination);
    }
    else{
      gAstro.Gal2Eq(Equinox, GalacticLatitude, GalacticLongitude, &gGenPar.RightAscension, &gGenPar.Declination);

    }
  }

  if(gGenPar.SourceRadius<0.){
    LOG("gSeaNuEvGen", pFATAL) << "Source radius must be equal or greater than zero, exiting...";
    gAbortingInErr = true;
    exit(1);
  }


  return;
}
//________________________________________________________________________________________
void LoadDefParams(void){
  // Load and re-define default parameters if request by the user

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

    xmlTextReaderPtr reader;

    string filename=gSystem->Getenv("DEFPARAM");

    if(!(access(filename.c_str(), F_OK) != -1)){
      LOG("gSeaNuEvGen", pFATAL) << filename <<" does not exist --> check env. parameter DEFPARAM; exiting...";
      exit(1);
    }
    else
      LOG("gSeaNuEvGen", pINFO) <<"Reading new default parameters from file "<< filename;


    reader = xmlNewTextReaderFilename(filename.c_str());

    string smedia,selem,sdensity;




    while(xmlTextReaderRead(reader)==1){

      xmlChar * name  = xmlTextReaderName     (reader);
      xmlChar * value = xmlTextReaderValue    (reader);
      int       type  = xmlTextReaderNodeType (reader);
      int       depth = xmlTextReaderDepth    (reader);

      if(depth==0 && type==1){
        if(xmlStrcmp(name, (const xmlChar *) "default_params")){
          LOG("gSeaNuEvGen", pFATAL) << " Invalid root element in filename " << filename <<"; exiting... ";
          exit(1);
        }
      }


      if(!xmlStrcmp(name, (const xmlChar *) "param" ) && type==1) {
        xmlChar * xelem   = xmlTextReaderGetAttribute(reader,(const xmlChar*)"name");
        selem = utils::str::TrimSpaces((const char *)xelem);
        xmlFree(xelem);
      }
      if(depth==2){
        if(selem=="RunNu"){
          gGenPar.RunNu=atoi((const char *)value);
        }
        else if(selem=="RanSeed"){
          gGenPar.RanSeed=atoi((const char *)value);
        }
        else if(selem=="Can"){
          string incan=(const char *)value;
          if(incan.find(",") != string::npos){
            // split the comma separated list
            vector<string> cansize = utils::str::Split(incan, ",");
            assert(cansize.size() == 3);
            gGenPar.Can1 = atof(cansize[0].c_str());
            gGenPar.Can2 = atof(cansize[1].c_str());
            gGenPar.Can3 = atof(cansize[2].c_str());
          }
        }
        else if(selem=="Origin"){
          string inorig=(const char *)value;
          if(inorig.find(",") != string::npos){
            // split the comma separated list
            vector<string> origsize = utils::str::Split(inorig, ",");
            assert(origsize.size() == 3);
            gGenPar.OriginX = atof(origsize[0].c_str());
            gGenPar.OriginY = atof(origsize[1].c_str());
            gGenPar.OriginZ = atof(origsize[2].c_str());
          }
        }
        else if(selem=="SiteLatitude"){
          gGenPar.SiteLatitude=atof((const char *)value);
        }
        else if(selem=="SiteLongitude"){
          gGenPar.SiteLongitude=atof((const char *)value);
        }
        else if(selem=="SiteDepth"){
          gGenPar.SiteDepth=atof((const char *)value);
        }
        else if(selem=="EvRange"){
          string inrange=(const char *)value;
          if(inrange.find(",") != string::npos){
            // split the comma separated list
            vector<string> rangesize = utils::str::Split(inrange, ",");
            assert(rangesize.size() == 2);
            gGenPar.EvMin = atof(rangesize[0].c_str());
            gGenPar.EvMax = atof(rangesize[1].c_str());
          }
        }
        else if(selem=="CtRange"){
          string inrange=(const char *)value;
          if(inrange.find(",") != string::npos){
            // split the comma separated list
            vector<string> rangesize = utils::str::Split(inrange, ",");
            assert(rangesize.size() == 2);
            gGenPar.CtMin = atof(rangesize[0].c_str());
            gGenPar.CtMax = atof(rangesize[1].c_str());
          }
        }
        else if(selem=="NBin"){
          gGenPar.NBin=atoi((const char *)value);
        }
        else if(selem=="Alpha"){
        gGenPar.Alpha=atof((const char *)value);
        }
        else if(selem=="RTOpt"){
          gGenPar.RTOpt=(const char *)value;
          gGenPar.RTOpt.erase(remove(gGenPar.RTOpt.begin(),gGenPar.RTOpt.end(),' '),gGenPar.RTOpt.end());
          gGenPar.RTOpt.erase(remove(gGenPar.RTOpt.begin(),gGenPar.RTOpt.end(),'\t'),gGenPar.RTOpt.end());
        }
        else if(selem=="NAbsLength"){
          gGenPar.NAbsLength=atof((const char *)value);
        }
        else if(selem=="AbsLength"){
          gGenPar.AbsLength=atof((const char *)value);
        }
        else if(selem=="HBedRock"){
          gGenPar.HBedRock=atof((const char *)value);
        }
        else if(selem=="PropCode"){
          gGenPar.PropCode=(const char *)value;
          gGenPar.PropCode.erase(remove(gGenPar.PropCode.begin(),gGenPar.PropCode.end(),' '),gGenPar.PropCode.end());
          gGenPar.PropCode.erase(remove(gGenPar.PropCode.begin(),gGenPar.PropCode.end(),'\t'),gGenPar.PropCode.end());
          gGenPar.PropCode = genie::utils::str::ToLower(gGenPar.PropCode);
          if(gGenPar.PropCode.compare("propamuon")==0)gGenPar.PropCode = "PropaMuon";
          if(gGenPar.PropCode.compare("music")==0){
#ifdef _MUSIC_ENABLED__
            gGenPar.PropCode = "MUSIC";
#else
            LOG("gSeaNuEvGen", pFATAL) <<"Propagation code: "<<gGenPar.PropCode<<", is not enabled, change file "<<filename <<"; exiting...";
            gAbortingInErr = true;
            exit(1);
#endif
          }
          if(gGenPar.PropCode.compare("proposal")==0){
#ifdef _PROPOSAL_ENABLED__
            gGenPar.PropCode = "PROPOSAL";
#else
            LOG("gSeaNuEvGen", pFATAL) <<"Propagation code: "<<gGenPar.PropCode<<", is not enabled, change file "<<filename <<"; exiting...";
            gAbortingInErr = true;
            exit(1);
#endif
          }
          if(gGenPar.PropCode.compare("jpp")==0){
#ifdef _JPP_ENABLED__
            gGenPar.PropCode = "JPP";
#else
            LOG("gSeaNuEvGen", pFATAL) <<"Propagation code: "<<gGenPar.PropCode<<", is not enabled, change file "<<filename <<"; exiting...";
            gAbortingInErr = true;
            exit(1);
#endif
          }
        }
//        else if(selem=="TGen"){
//          gGenPar.TGen=atof((const char *)value);
//        }
        else if(selem=="WriteParticlesMode"){
          gGenPar.WriteParticlesMode=atoi((const char *)value);
        }
        else if(selem=="WeightOpt"){
        gGenPar.WeightOpt=atoi((const char *)value);
        }
        else if(selem=="EvFilePrefix"){
          gGenPar.EvFilePrefix=(const char *)value;
          gGenPar.EvFilePrefix.erase(remove(gGenPar.EvFilePrefix.begin(),gGenPar.EvFilePrefix.end(),' '),gGenPar.EvFilePrefix.end());
          gGenPar.EvFilePrefix.erase(remove(gGenPar.EvFilePrefix.begin(),gGenPar.EvFilePrefix.end(),'\t'),gGenPar.EvFilePrefix.end());
        }
        else if(selem=="InpXSecFile"){
          gGenPar.InpXSecFile=(const char *)value;
          gGenPar.InpXSecFile.erase(remove(gGenPar.InpXSecFile.begin(),gGenPar.InpXSecFile.end(),' '),gGenPar.InpXSecFile.end());
          gGenPar.InpXSecFile.erase(remove(gGenPar.InpXSecFile.begin(),gGenPar.InpXSecFile.end(),'\t'),gGenPar.InpXSecFile.end());
        }
        else if(selem=="EventGeneratorList"){
          gGenPar.EventGeneratorList=(const char *)value;
          gGenPar.EventGeneratorList.erase(remove(gGenPar.EventGeneratorList.begin(),gGenPar.EventGeneratorList.end(),' '),gGenPar.EventGeneratorList.end());
          gGenPar.EventGeneratorList.erase(remove(gGenPar.EventGeneratorList.begin(),gGenPar.EventGeneratorList.end(),'\t'),gGenPar.EventGeneratorList.end());
        }

        else if(selem=="UseUTMSystem"){
          string CheckIfUTM = (const char *)value;
          CheckIfUTM.erase(remove(CheckIfUTM.begin(),CheckIfUTM.end(),' '),CheckIfUTM.end());
          CheckIfUTM.erase(remove(CheckIfUTM.begin(),CheckIfUTM.end(),'\t'),CheckIfUTM.end());
          CheckIfUTM = genie::utils::str::ToLower(CheckIfUTM);
          if(CheckIfUTM.compare("true")==0)gGenPar.UseUTMSystem = true;
          else if(CheckIfUTM.compare("false")==0)gGenPar.UseUTMSystem = false;
          else {
            LOG("gSeaNuEvGen", pFATAL) <<"Error reading new default parameters: Param UseUTMSystem must be true or false; exiting...";
            gAbortingInErr = true;
            exit(1);
          }
        }
      }

      xmlFree(name);
      xmlFree(value);
    }

    }

  return;
  }
//________________________________________________________________________________________
void LoadMediaComp(void){
  // Load and re-define media composition if request by the user

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


    bool CheckSeaWater=true;
    bool CheckRock=true;
    bool CheckMantle=true;
    bool CheckCore=true;

    xmlTextReaderPtr reader;

    string filename=gSystem->Getenv("MEDIACOMP");

    if(!(access(filename.c_str(), F_OK) != -1)){
      LOG("gSeaNuEvGen", pFATAL) << filename <<" does not exist --> check env. parameter MEDIACOMP; exiting...";
      exit(1);
    }
    else
      LOG("gSeaNuEvGen", pINFO) <<"Reading new media composition from file "<< filename;

    reader = xmlNewTextReaderFilename(filename.c_str());

    string smedia,selem,sdensity;
    int pdg=0;

    while(xmlTextReaderRead(reader)==1){

      xmlChar * name  = xmlTextReaderName     (reader);
      xmlChar * value = xmlTextReaderValue    (reader);
      int       type  = xmlTextReaderNodeType (reader);
      int       depth = xmlTextReaderDepth    (reader);

      if(depth==0 && type==1){
        if(xmlStrcmp(name, (const xmlChar *) "media_comp")){
        LOG("gSeaNuEvGen", pFATAL) << " invalid root element in filename " << filename <<"; exiting... ";
        exit(1);
       }
      }
      if(!xmlStrcmp(name, (const xmlChar*) "param_set" ) && type==1){
        xmlChar * xmedia   = xmlTextReaderGetAttribute(reader,(const xmlChar*)"media");
        smedia = utils::str::TrimSpaces((const char *)xmedia);

        xmlChar * xdensity = xmlTextReaderGetAttribute(reader,(const xmlChar*)"density");

        if(smedia=="SeaWater" && CheckSeaWater){
          CheckSeaWater=false;
          gGenPar.SeaWaterComp.clear();
          LOG("gSeaNuEvGen", pINFO) <<"Defining new SeaWater composition";
          if(xdensity!=NULL){
            sdensity = utils::str::TrimSpaces((const char*)xdensity);
            LOG("gSeaNuEvGen",pINFO) <<"Defining new SeaWater density";
            gGenPar.RhoSW=(double)atof(sdensity.c_str());
          }
        }
        if(smedia=="Rock" && CheckRock){
          CheckRock=false;
          gGenPar.RockComp.clear();
          LOG("gSeaNuEvGen", pINFO) <<"Defining new Rock composition";
          if(xdensity!=NULL){
            sdensity = utils::str::TrimSpaces((const char*)xdensity);
            LOG("gSeaNuEvGen",pINFO) <<"Defining new Rock density";
            gGenPar.RhoSR=(double)atof(sdensity.c_str());
          }
        }
        if(smedia=="Mantle" && CheckMantle){
          CheckMantle=false;
          gGenPar.MantleComp.clear();
          LOG("gSeaNuEvGen", pINFO) <<"Defining new Mantle composition";
          if(xdensity!=NULL) LOG("gSeaNuEvGen",pINFO) <<"Defining new Mantle density, ignoring...";
        }
        if(smedia=="Core" && CheckCore){
          CheckCore=false;
          gGenPar.CoreComp.clear();
          LOG("gSeaNuEvGen", pINFO) <<"Defining new Core composition";
          if(xdensity!=NULL) LOG("gSeaNuEvGen",pINFO) <<"Defining new Core density, ignoring...";
        }

        xmlFree(xmedia);
      }

      if(!xmlStrcmp(name, (const xmlChar *) "param" ) && type==1) {
        xmlChar * xelem   = xmlTextReaderGetAttribute(reader,(const xmlChar*)"name");
        selem = utils::str::TrimSpaces((const char *)xelem);
        xmlFree(xelem);
      }

      if(depth==3){
        if(selem=="H")       pdg=1000010010;            // H1
        else if(selem=="C")  pdg=1000060120;            // C12
        else if(selem=="O")  pdg=1000080160;            // O16
        else if(selem=="Na") pdg=1000110230;            // Na23
        else if(selem=="Mg") pdg=1000120240;            // Mg24
        else if(selem=="Al") pdg=1000130270;            // Al27
        else if(selem=="Si") pdg=1000140280;            // Si28
        else if(selem=="S")  pdg=1000160320;            // S32
        else if(selem=="Cl") pdg=1000170350;            // Cl35
        else if(selem=="K")  pdg=1000190390;            // K39
        else if(selem=="Ca") pdg=1000200400;            // Ca40
        else if(selem=="Ti") pdg=1000220480;            // Ti48
        else if(selem=="Fe") pdg=1000260560;            // Fe56
        else if(selem=="Ni") pdg=1000280580;            // Ni58
        else if(selem=="Br") pdg=1000350800;            // Br80

        else {
          LOG("gSeaNuEvGen", pFATAL) << smedia <<": element "<< selem << " not present in the list..., exiting... ";
          exit(1);
        }
        if(smedia=="SeaWater")  gGenPar.SeaWaterComp.insert(map<int, double>::value_type(pdg,atof((const char *)value)));
        if(smedia=="Rock")      gGenPar.RockComp.insert(map<int, double>::value_type(pdg,atof((const char *)value)));
        if(smedia=="Mantle")    gGenPar.MantleComp.insert(map<int, double>::value_type(pdg,atof((const char *)value)));
        if(smedia=="Core")      gGenPar.CoreComp.insert(map<int, double>::value_type(pdg,atof((const char *)value)));
      }
      xmlFree(name);
      xmlFree(value);
    }
  }

  return;
}

//________________________________________________________________________________________
void PrintStat(void)
{

  char run[1024];
  sprintf(run,"%d",gGenPar.RunNu);

  string filename=gGenPar.EvFilePrefix;
  filename.append(".");
  filename.append(run);
  filename.append(".binstat");

  FILE *fstat;
  fstat=fopen(filename.c_str(),"w");

  fprintf(fstat,"Bin      EneMin        EneMax       NEvBin NEvBinScal  NEvBinGen         RT         RL         X          Y          Z        RVol  HSeaWater      HRock    ProbScale \n");

  string format = "%2d %12.2f  %12.2f %12.2E %10d %10d %10.3f %10.3f %10.3f %10.3f %10.3f %10.3f %10.3f %10.3f %12.2E\n";
  for(int iBin=0; iBin<gGenPar.NBin; iBin++){
    fprintf(fstat,format.c_str(),iBin+1,VecBinPar[iBin].EneMin,VecBinPar[iBin].EneMax, VecBinPar[iBin].NEvBin,VecBinPar[iBin].NEvBinScal, VecBinPar[iBin].NEvBinGen,
         VecBinPar[iBin].RT,VecBinPar[iBin].RL,VecBinPar[iBin].X0,VecBinPar[iBin].Y0,VecBinPar[iBin].Z0,VecBinPar[iBin].RVol,VecBinPar[iBin].HSeaWater,VecBinPar[iBin].HRock,VecBinPar[iBin].ProbScale);
  }



  fclose(fstat);

  return;
}
//________________________________________________________________________________________
void CreaBin(void)
{
  BinParam param;
  TRandom * tRandom = new TRandom();

  double DeltaE=(log10(gGenPar.EvMax)-log10(gGenPar.EvMin))/gGenPar.NBin;
  double IE,IETOT;
  char bin[1024];
  if(gGenPar.Alpha==1)
    IETOT=log(gGenPar.EvMax/gGenPar.EvMin);
  else
    IETOT=(pow(gGenPar.EvMax,(1.-gGenPar.Alpha))-pow(gGenPar.EvMin,(1.-gGenPar.Alpha)))/(1.-gGenPar.Alpha);
  GSeaGeometry * SeaGeom = new GSeaGeometry(&gGenPar);

  for (int iBin=0; iBin<gGenPar.NBin; iBin++){

    param.EneMin = log10(gGenPar.EvMin)+iBin*DeltaE;
    param.EneMax = param.EneMin+DeltaE;

    param.EneMin = pow(10,param.EneMin);
    param.EneMax = pow(10,param.EneMax);

    if(gGenPar.Alpha==1)
      IE=log(param.EneMax/param.EneMin);
    else
      IE=(pow(param.EneMax,(1.-gGenPar.Alpha))-pow(param.EneMin,(1.-gGenPar.Alpha)))/(1.-gGenPar.Alpha);

    if(gGenPar.TrackEvts==0 && iBin==0){
      param.geometry = "geometry-el.root";
      SeaGeom->BuildGeometryElectron(param.geometry);
    }
    else if(gGenPar.TrackEvts==1){
      param.geometry = "geometry-mu-";
      sprintf(bin,"%d",iBin+1);
      param.geometry.append(bin);
      param.geometry.append(".root");
      SeaGeom->BuildGeometryMuon(param.EneMax, param.geometry);
    }
#ifdef _HETAU_ENABLED__
    else if(gGenPar.TrackEvts==2){
      param.geometry = "geometry-tau-muon-";
      sprintf(bin,"%d",iBin+1);
      param.geometry.append(bin);
      param.geometry.append(".root");
      SeaGeom->BuildGeometryTauTrack(param.EneMax, param.geometry);
    }
    else if(gGenPar.TrackEvts==3){
      param.geometry = "geometry-tau-shower-";
      sprintf(bin,"%d",iBin+1);
      param.geometry.append(bin);
      param.geometry.append(".root");
      SeaGeom->BuildGeometryTauShower(param.EneMax, param.geometry);
    }
#endif

    if(gGenPar.PropMode)param.geometry = "geometry-Earth.root";

    SeaGeom->CalcRadii();

    param.RT = SeaGeom->GetRT();
    param.RL = SeaGeom->GetRL();


    param.X0 = SeaGeom->GetX0();
    param.Y0 = SeaGeom->GetY0();
    param.Z0 = SeaGeom->GetZ0();

    param.RVol = SeaGeom->GetRVol();
    param.HSeaWater = SeaGeom->GetHSeaWater();
    param.HRock = SeaGeom->GetHRock();
    param.RockCenter = SeaGeom->GetRockCenter();
    param.SeaCenter = SeaGeom->GetSeaCenter();

    gGenPar.RangeSW = SeaGeom->GetRangeSW();
    gGenPar.RangeSR = SeaGeom->GetRangeSR();

    if(gGenPar.GenMode=="POINT" || gGenPar.GenMode=="DIFFUSE")
      param.NEvBin = gGenPar.NTot*(IE/IETOT);
    else
      param.NEvBin = 0;


    param.ProbScale = 0.;
    param.ProbScaleSW = 0.;
    param.ProbScaleSR = 0.;
    param.NEvBinGenSW = 0;
    param.NEvBinGenSR = 0;
    param.NEvBinGen = 0;

    param.NEvBinScal = 0;
    param.NEvBinScalSW = 0;
    param.NEvBinScalSR = 0;

    VecBinPar.push_back(param);

  }

  if(gGenPar.PropMode)gGenPar.SeaBottomRadius=SeaGeom->GetSeaBottomRadius()*1E3;

  for (unsigned iBin=0; iBin<VecBinPar.size(); iBin++)
    if(gGenPar.GenMode=="POINT" || gGenPar.GenMode=="DIFFUSE")
      if (VecBinPar[iBin].RT<=0.)
        VecBinPar[iBin].NEvBinScal = tRandom->Poisson(VecBinPar[iBin].NEvBin);
      else
        VecBinPar[iBin].NEvBinScal = tRandom->Poisson(VecBinPar[iBin].NEvBin*pow((VecBinPar[iBin].RT/VecBinPar[VecBinPar.size()-1].RT),2));
    else
      VecBinPar[iBin].NEvBinScal = 0;


  gGenPar.HRock     = VecBinPar[VecBinPar.size()-1].HRock;
  gGenPar.HSeaWater = VecBinPar[VecBinPar.size()-1].HSeaWater;
  gGenPar.RVol      = VecBinPar[VecBinPar.size()-1].RVol;
  gGenPar.SeaCenter = VecBinPar[VecBinPar.size()-1].SeaCenter;
  gGenPar.RockCenter = VecBinPar[VecBinPar.size()-1].RockCenter;

  PrintStat();

  delete SeaGeom;
  delete tRandom;

  return;
}

//________________________________________________________________________________________
void PrintName(void)
{

  string ver=gGenPar.gSeaGenVer;
  int iBlanc=9-ver.size();
  for(int i=0; i<iBlanc; i++)ver.append(" ");

  cout
   << "\n"
   << "\n"
   << "**********************************************************************************************************\n"
   << "**********************************************************************************************************\n"
   << "**                                                                                                      **\n"
   << "**                                                                                                      **\n"
   << "**                                                                                                      **\n"
   << "**                        .oooooo..o                       .oooooo.                                     **\n"
   << "**                       d8P'    `Y8                      d8P'  `Y8b                                    **\n"
   << "**             .oooooooo Y88bo.       .ooooo.   .oooo.   888            .ooooo.  ooo. .oo.              **\n"
   << "**            888' `88b   `\"Y8888o.  d88' `88b `P  )88b  888           d88' `88b `888P\"Y88b             **\n"
   << "**            888   888       `\"Y88b 888ooo888  .oP\"888  888     ooooo 888ooo888  888   888             **\n"
   << "**            `88bod8P'  oo     .d8P 888    .o d8(  888  `88.    .88'  888    .o  888   888             **\n"
   << "**            `8oooooo.  8\"\"88888P'  `Y8bod8P' `Y888\"\"8o  `Y8bood8P'   `Y8bod8P' o888o o888o            **\n"
   << "**            d\"     YD                                                                                 **\n"
   << "**            \"Y88888P'                                                                                 **\n"
   << "**                                                                                                      **\n"
   << "**                               a GENIE-based code for neutrino telescopes                             **\n"
   << "**                                                                                                      **\n"
   << "**                                            version  "<<ver<<"                                        **\n"
   << "**                                                                                                      **\n"
   << "**   The gSeaGen code is a neutrino event generator for neutrino telescopes. It simulates particles     **\n"
   << "**   created in all-flavour neutrino interactions, which may produce detectable Cherenkov light.        **\n"
   << "**   gSeaGen uses GENIE (http://www.genie-mc.org) to simulate the neutrino interactions.                **\n"
   << "**                                                                                                      **\n"
   << "**   Developers: Carla Distefano, INFN-LNS, Via S.Sofia 62, 95123, Catania, Italy                       **\n"
   << "**               e-mail: distefano_c@lns.infn.it                                                        **\n"
   << "**               Alfonso Andres Garcia Soto, Nikhef, Science Park 904, Amsterdam 1098 XH, Netherlands   **\n"
   << "**               e-mail: alfonsog@nikhef.nl                                                             **\n"
   << "**               Piotr Kalaczynski, 47National Centre for Nuclear Research, 02-093 Warsaw, Poland       **\n"
   << "**               e-mail: pkalaczynski@km3net.de                                                         **\n"
   << "**                                                                                                      **\n"
   << "**   Repository: https://git.km3net.de/simulation/gseagen                                               **\n"
   << "**                                                                                                      **\n"
   << "**   Reference:  S. Aiello et al. (the KM3NeT Collaboration),                                           **\n"
   << "**               Computer Physics Communications 256, 107477 (2020)                                     **\n"
   << "**                                                                                                      **\n"
   << "**   Web page:   https://git.km3net.de/simulation/gseagen/wikis                                         **\n"
   << "**                                                                                                      **\n"
   << "**   Copyright (C) 2012-2021  The KM3NeT Collaboration                                                  **\n"
   << "**                                                                                                      **\n"
   << "**   This program is free software: you can redistribute it and/or modify it under the terms of the     **\n"
   << "**   GNU General Public License as published by the Free Software Foundation, either version 3 of       **\n"
   << "**   the License, or (at your option) any later version.                                                **\n"
   << "**   This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY;          **\n"
   << "**   without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.          **\n"
   << "**   See the GNU General Public License for more details.                                               **\n"
   << "**   You should have received a copy of the GNU General Public License along with this program.         **\n"
   << "**   If not, see <http://www.gnu.org/licenses/>.                                                        **\n"
   << "**                                                                                                      **\n"
   << "**********************************************************************************************************\n"
   << "**********************************************************************************************************\n"
   << "\n";

  sleep(1);
}
//________________________________________________________________________________________
void PrintSyntax(void)
{

#ifdef _MUSIC_ENABLED__
  string MUSICPATH = BOOST_PP_STRINGIZE(_MUSICDATA__);
  while(1){
    size_t pos = MUSICPATH.find("\"");
    if(pos==string::npos)break;
    MUSICPATH.erase(pos,1);
  }
#endif

  cout
   << "\n"
   << "\n **Syntax**"
   << "\n gSeaNuEvGen "
   << "\n           [-a generation_spectral_index]"
   << "\n           [-b n_of_bins]"
   << "\n           [-bedrock bedrock_height]"
   << "\n           [-c can_definition]"
   << "\n           [-C]"
   << "\n           [--coord site_coordinate_definition]"
#if defined _ANTARES_ENABLED__  || defined _KM3NET_ENABLED__
   << "\n           [--site-dec site_declination]"
#endif
   << "\n           [--cross-sections XSecName]"
   << "\n           [-d depth]"
   << "\n           [-e min_energy,max_energy]"
   << "\n           [--event-generator-list list_name]"
   << "\n           [ -f simulation:flux[neutrino_code],...]"
#ifdef _ANTARES_ENABLED__
   << "\n           [-l light_absorption_lenght,[n_light_absorption_lenght]]"
#endif
   << "\n           [ -n n_of_events]"
   << "\n           [-o output_event_file_prefix]"
   << "\n           [-origin X0,Y0,Z0]"
   << "\n           [-p]"
   << "\n           [-point source_info]"
   << "\n           [-prop prop_code]"
   << "\n           [-r run_number]"
#ifdef _ANTARES_ENABLED__
   << "\n           [-R coordinate_rotation_matrix]"
   << "\n           [-rt RT_option]"
#endif
   << "\n           [--seed random_number_seed]"
   << "\n           [-utm]"
   << "\n           [-t min_costheta,max_costheta]"
   << "\n           [-time time_format:StartTime,StopTime]"
   << "\n           [-tree tree_format]"
   << "\n           [--tune generator_tune]"
   << "\n           [-w]"
#if defined _ANTARES_ENABLED__  || defined _KM3NET_ENABLED__
   << "\n           [-write write_particles_mode]"
#endif
   << "\n           [-save saved_particles]"
   << "\n           [-chances chances_per_shower]"
   << "\n           [--rot-showers UNIT:yaw,pitch,roll]"
   << "\n           [-wgt weighting scheme]"
   << "\n           [--cr-flux CR flux model used to weight CORSIKA events]"
#ifdef _HETAU_ENABLED__
   << "\n           [-tau decay_mode]"
#endif
   << "\n"
   << "\n **Code Options**"
   << "\n"
   << "\n [] Denotes an optional argument."
   << "\n"
   << "\n -a        Specifies the generation spectral index [default: 1.4]."
   << "\n"
   << "\n -b        Specifies the number of energy bins [default: 1]. "
   << "\n"
   << "\n -bedrock  Specifies the height of the bedrock just below the seabottom (in meters) [default: 0.]. "
   << "\n"
   << "\n -c        Specifies input for the detector can definition. "
   << "\n           Possible choices are: "
   << "\n             1)  '-c CAN:Zmin,Zmax,Rad' "
   << "\n               where Zmin and Zmax are bottom and top Z coordinates of the can "
   << "\n               and Rad is its radius in meters. eg. 'CAN:0.000,659.862,309.450' [default]"
   << "\n               The site depth will be set to the default value if not specified using -d option."
   << "\n               Take care that the depth site enters in the definition of the interaction volume for down-going neutrinos!"
   << "\n               The detector coordinate sites are set to the default values if not specified using -coord option."
   << "\n"
#ifdef _ANTARES_ENABLED__
   << "\n             2)  '-c DET:/path/geometry.det' "
   << "\n               where geometry.det is the detector geometry file in the ANTARES format. "
   << "\n               In this case, the can will be defined from the detector size adding n [default n=3]"
   << "\n               times the light absorption lenght [default La=70 m] (see -l option). "
   << "\n               The site depth and coordinates are read from geometry.det file. Site declination is"
   << "\n               set to the default value of 2.28*kPi/180. RAD"
   << "\n"
   << "\n             3)  '-c DETX:/path/geometry.detx' "
   << "\n               where geometry.detx is the detector geometry file in the KM3NeT format (vesion V2)."
   << "\n               In this case, the can will be defined from the detector size adding n [default n=3]"
   << "\n               times the light absorption lenght [default La=70 m] (see -l option). "
   << "\n               The site depth is computed from the UTM reference point, the site coordinates"
   << "\n               are computed from the detector barycenter UTM coordinates. Site declination is"
   << "\n               set to the default value of 2.28*kPi/180. RAD"
   << "\n"
#endif
   << "\n -C        Force the interaction volume to be the detector CAN. This option has  "
   << "\n           an effect only in case of muon and tau CC interaction. In the other "
   << "\n           cases the interaction volume is already the CAN."
   << "\n"
   << "\n --coord   Specifies the detector site coordinates. "
   << "\n           Possible choices are: "
   << "\n             1)  --coord UNIT:latitude,longitude [,declination]  "
   << "\n               where detector site latitude and longitude are in radians (UNIT=RAD) or in degrees (UNIT=DEG) [default RAD:0.74700,0.10763]."
   << "\n               The longitude is positive eastwards from the meridian of Greenwich, and negative to the West."
   << "\n               Additionally, magnetic declination of the detector site may be specified in radians (UNIT=RAD) or in degrees (UNIT=DEG) [default RAD:2.28*kPi/180.]"
   << "\n"
   << "\n             2)  --coord SITE:site_name  "
   << "\n               detector site coordinates and magnetic declination are set by site name. Currently the predefined sites are: ARCA, ORCA, IceCube, GVD."
   << "\n"
#if defined _ANTARES_ENABLED__  || defined _KM3NET_ENABLED__
   << "\n --site-dec site_declination "
   << "\n           Set detector site declination by hand. May be used to overwrite the default value used when DET or DETX files are given in -c option"
   << "\n           Input ignored if option -c CAN:Zmin,Zmax,Rad is used "
#endif
   << "\n"
   << "\n --cross-sections XSecName "
   << "\n           Name of an XML file with pre-computed"
   << "\n           cross-section values used for constructing splines [default: gxspl-FNALsmall.xml]"
   << "\n"
   << "\n           eg. '--cross-sections gxspl-FNALsmall.xml' "
   << "\n"
   << "\n           Despite their names, files in $GSEAGEN/dat directory contain cross sections also for Rock, Mantle and Crust media! "
   << "\n"
   << "\n -d        Specifies the site depth in meters [default: 2475]. "
   << "\n           This option has an effect only if the detector can is defined by hand (see -c option)."
   << "\n"
   << "\n -e        Specifies the neutrino energy in GeV [default: 1,1000]. "
   << "\n           Must be a comma-separated pair of numbers, eg. `-e 0.3,70' "
   << "\n"
   << "\n --event-generator-list list_name"
   << "\n           List of event generators to load. Valid settings are the XML block names appearing in"
   << "\n           config/EventGeneratorListAssembler.xml [default: 'Default']."
   << "\n"
   << "\n           eg. '--event-generator-list  CC' "
   << "\n           eg. '--event-generator-list  NC' "
   << "\n"
   << "\n -f        Specifies the neutrino types to be generated and the neutrino flux models used to weight the events."
   << "\n           The general syntax is: `-f simulation:flux,...[neutrino_code],flux...[neutrino_code];simulation:...`"
   << "\n           where the `simulation' string specifies the model, 'flux' is the neutrino flux or the file containing it"
   << "\n           and the 'neutrino_code' is the neutrino type in the PDG code (12: nu_e, -12: anti_nu_e, 14: nu_mu,"
   << "\n           -14: anti_nu_mu, 16: nu_tau, -16: anti_nu_tau)."
   << "\n"
   << "\n           Alternatively, if gSeaGen is not used to generate particles, but only propagate the ones generated"
   << "\n           with an external code (CORSIKA), -f keyword may be used in 2 ways:"
   << "\n           `-f EVT:corant_file.evt` for files preprocessed with corant (https://git.km3net.de/simulation/corant)"
   << "\n           and"
   << "\n           `-f BIN:corsika_binary_file` for binary CORSIKA output files"
   << "\n"
   << "\n           Notes: "
   << "\n            - Possible values for `simulation' are: FLUKA, BARTOL, HONDA or FUNCX (where X=1,2,3)."
   << "\n              In case of FLUKA, BARTOL or HONDA, flux is the file name (incl. full path)."
   << "\n              See:"
   << "\n              Bartol fluxes: http://www-pnp.physics.ox.ac.uk/~barr/fluxfiles"
   << "\n              Fluka fluxes:  http://pcbat1.mi.infn.it/~battist/neutrino.html"
   << "\n              Honda fluxes:  http://www.icrr.u-tokyo.ac.jp/~mhonda/nflx2014/index.html"
   << "\n"
   << "\n              eg. '-f BARTOL:dat/f210_3_z.kam_num[14],dat/f210_3_z.kam_nbm[-14]' "
   << "\n"
   << "\n            - More than one file can be defined for each neutrino code. "
   << "\n"
   << "\n              eg. '-f BARTOL:dat/fmin10_0401z.kam_nue,dat/f210_3_z.kam_nue[12]' "
   << "\n"
   << "\n              In case of `FUNC1', flux is a mathematical function (c++ syntax) depending on x, which represents "
   << "\n              the neutrino energy in units of GeV. The flux is in unit of 1/(GeV m2 s sr). "
   << "\n"
   << "\n              eg. '-f FUNC1:1E-5*pow(x,-2)[14]' "
   << "\n"
   << "\n              Two-dimensional  functions  depending  on  x  and  y  can  be  provided to define a dependence on the  "
   << "\n              neutrino arrival direction.  In this case the tag simul has value FUNC2 and the variable y represents cos(theta).  "
   << "\n              Similarly three-dimensional functions depending also on the variable z (representing the angle phi in radians)  "
   << "\n              can be provided by setting the tag simul to FUNC3. "
   << "\n"
   << "\n            - More than one 'simulation' input can be defined. "
   << "\n              Fluxes corresponding to the same neutrino code will added. "
   << "\n"
   << "\n              eg. '-f BARTOL:dat/f210_3_z.kam_num[14];FUNC1:x[14]' "
   << "\n"
   << "\n            - More than one -f option is allowed. "
   << "\n              Fluxes corresponding to the same neutrino code will added. "
   << "\n"
   << "\n              eg. '-f BARTOL:dat/f210_3_z.kam_nue[12] -f BGLRS:dat/f210_3_z.kam_nbm[-12]' "
   << "\n"
   << "\n            - Values for `simulation' and 'flux' are not mandatory. In this case the code will set by default: "
   << "\n              simulation= FUNC1 and flux= 1E-5*pow(x,-2) "
   << "\n"
   << "\n              eg. '-f [14]' "
   << "\n              eg. '-f [12,14]' "
   << "\n              eg. '-f [16,-16]' "
   << "\n"
#ifdef _ANTARES_ENABLED__
   << "\n -l light_absorption_lenght,[n_light_absorption_lenght]       "
   << "\n"
   << "\n           Specifies the light absorption length in meters and the number of light absorption lengths to calculate "
   << "\n           the can size from a detector geometry file (see option -c)  [default: 70,3]. "
   << "\n           It is possible to input only the light absorption length "
   << "\n"
   << "\n              eg. '-l 60.' "
   << "\n"
   << "\n           In this case default value of n_light_absorption_lenght=3 is used. "
#endif
   << "\n"
   << "\n -n        Specifies how many neutrinos to generate. "
   << "\n"
   << "\n -o        Sets the prefix of the output event file [default: gseagen]. "
   << "\n           The output filename is built as: "
   << "\n           [prefix].[run_number].evt "
   << "\n"
   << "\n -origin X0,Y0,Z0 "
   << "\n           Sets the origin position coordinates (in meters) in the detector rest frame [default: 0.,0.,0.]. "
   << "\n"
   << "\n -p        Activates the neutrino propagation through the Earth [default: no propagation]."
   << "\n"
   << "\n -prop     Set the muon propagation code."
   << "\n           Possible choices are: "
   << "\n             1)  -prop 'PropaMuon' [internal code, default]"
#ifdef _MUSIC_ENABLED__
   << "\n             2)  -prop 'MUSIC' "
#endif
#ifdef _PROPOSAL_ENABLED__
   << "\n             3)  -prop 'PROPOSAL' "
#endif
#ifdef _JPP_ENABLED__
   << "\n             3)  -prop 'JPP' "
#endif
   << "\n"
   << "\n -point    Activates the generation from point-like (and extended) sources. "
   << "\n           Possible choices are: "
   << "\n             1)  -point 'DEC:src_decl'"
   << "\n               where src_decl is the source declination in degrees. The source right ascension will be RA=0; "
//   << "\n               the generated interval time will be TGen = 31556926 sec (1 yr as in the diffuse mode); "
   << "\n               the Local Sidereal Time will be generated between 0h and 24h."
   << "\n"
   << "\n             2)  -point 'FILE:AstroSource.xml'"
   << "\n                info about the source are read from file AstroSource.xml"
   << "\n                The user has the possibility to set:"
   <<" \n                   - the source coordinates using equatorial or galactic coordinates;"
//   << "\n                   - the generated interval time in MJDs or dates: TGen = (MJDStop-MJDStart)*86400 sec"
   << "\n                   - the source angular radius (extended source)"
   << "\n                See $GSEAGEN/dat/AstroSource.xml for an example."
   << "\n"
   << "\n           If point/extended mode is activated, only tag FUNC for -f option is available."
   << "\n"
   << "\n -r        Specifies the MC run number [default: 100000000]. "
   << "\n"
#ifdef _ANTARES_ENABLED__
   << "\n -R        Input rotation matrix for transforming the flux neutrino coordinates "
   << "\n           from the default Topocentric Horizontal (see GENIE manual) coordinate "
   << "\n           system to the user-defined topocentric coordinate system.  "
   << "\n           The rotation is specified by the 3 Euler angles (phi, theta, psi). "
   << "\n           The Euler angles are input as a comma separated list as: "
   << "\n           `-R <convention>:phi,theta,psi', "
   << "\n           where <convention> is either X (for X-convention), Y (for Y-convention), "
   << "\n           X^-1 or Y^-1 (as previously, but using the inverse matrix). "
   << "\n           By default, the X-convention (rotation about Z-axis, then about the  "
   << "\n           new X-axis, then about the Z-axis) is used.  "
   << "\n           Notes:"
   << "\n             - (Extract from TRotation documentation) "
   << "\n             'Euler angles usually define the rotation of the new coordinate  "
   << "\n              system with respect to the original system, however, the TRotation  "
   << "\n              class specifies the rotation of the object in the original system  "
   << "\n              (an active rotation). To recover the usual Euler rotations (ie.  "
   << "\n              rotate the system not the object), you must take the inverse of  "
   << "\n              the rotation.' "
   << "\n           Examples: "
   << "\n             1. To set the Euler angles phi=3.14, theta=1.28, psi=1.0 using the "
   << "\n             X-convention, type: `-R 3.14,1.28,1.0', or `-R X:3.14,1.28,1.0' "
   << "\n             2. To set the Euler angles phi=3.14, theta=1.28, psi=1.0 using the "
   << "\n             Y-convention, type: `-R Y:3.14,1.28,1.0' "
   << "\n             3. To set the Euler angles phi=3.14, theta=1.28, psi=1.0 using the "
   << "\n             Y-convention, and then use the inverse rotation matrix, type: "
   << "\n             `-R Y^-1:3.14,1.28,1.0'  "
   << "\n"
   << "\n -rt       Set option to define the generation area radius RT"
   << "\n           Possible choices are: "
   << "\n             1)  -rt 'can' "
   << "\n               the generation surface will be a circle covering the detector can projected onto a plane perpendicular to each neutrino direction."
   << "\n             2)  -rt 'vol' "
   << "\n               the generation surface will be a circle covering the interaction volume projected onto a plane perpendicular to each neutrino direction."
   << "\n             3)  -rt proj "
   << "\n               the generation surface will be the detector can projected onto a plane perpendicular to each neutrino direction."
   << "\n             4)  -rt value "
   << "\n               the generation surface radius will RT = value (in meters)."
   << "\n"
#endif
   << "\n --seed random_number_seed [default: 12345]"
   << "\n"
   << "\n -utm      Use UTM detector system  [default: false]"
   << "\n"
   << "\n -t        Specifies the neutrino angular range. [default: 0.,1.] (up-going neutrinos)"
   << "\n           Must be a comma-separated pair of numbers, eg. `-t 0.,1.' "
   << "\n"
   << "\n -time     Specifies the generation time interval. "
   << "\n           The general syntax is: `-time time_format:StartTime,StopTime"
   << "\n           where the time_format string specifies the format for StartTime and StopTime."
   << "\n"
   << "\n           Notes: "
   << "\n            - Possible values for time_format are: "
   << "\n              TS: StartTime and StopTime are in Unix Time Stamp."
   << "\n              MJD: StartTime and StopTime are in Modified Julian Date."
   << "\n              DATE: StartTime and StopTime are in date in format dd-mm-yyyy;hh:mm:ss"
   << "\n"
   << "\n            - Value for time_format is not mandatory. In this case the code will set by default:"
   << "\n              time_format= TS "
   << "\n"
   << "\n -tree     Specifies the output ROOT file format. "
   << " [default: st]"
   << "\n           Possible choices are: "
   << "\n           1) st    --> SimpleTree format: tree branches are instances of simple variables"
   << "\n           2) et    --> EventTree format: tree branches are instances of gSeaGen classes"
   << "\n"
   << "\n --tune    Specifies the generator tune, required for GENIE version >=3.00.00 [default: G18_02a_00_000]"
   << "\n           Example: "
   << "\n           --tune G18_02a_00_000 "
   << "\n"
   << "\n -w        Writes the native GENIE output in [prefix].[run_number].root "
   << "\n"
#if defined _ANTARES_ENABLED__  || defined _KM3NET_ENABLED__
   << "\n -write    Choose the particles that are stored in the output file [default: 0] "
   << "\n           Possible choices are: "
   << "\n             1)  -write 0 "
   << "\n               only write events with lepton in can and store final tracks."
   << "\n             2)  -write 1 "
   << "\n               only write events with lepton in can and store intermidiate and final tracks."
   << "\n             3)  -write 2 "
   << "\n               write all generated events and store intermidiate and final tracks."
   << "\n"
#endif
   << "\n -save     Choose the particles produced in CORSIKA that are saved for propagation  [default: 'mu']"
   << "\n           Must be used with -write 2 except for -save 'mu', which will work with any -write option"
   << "\n           Possible choices are: "
   << "\n             1)  -save 'mu' "
   << "\n               only saves muons."
   << "\n             2)  -save 'nu' "
   << "\n               only saves neutrinos."
   << "\n             2)  -save 'lep' "
   << "\n               saves both muons and neutrinos."
   << "\n             3)  -save 'all' "
   << "\n               saves muons, neutrinos and all other particles that reached the sea level."
   << "\n"
   << "\n -chances  Number of extra chances each CORSIKA shower will get to hit the can [default: 0]"
   << "\n"
   << "\n --rot-showers     Specifies the custom rotation of directional cosines for all the CORSIKA showers at the sea level"
   << "\n                   in radians (UNIT=RAD) or in degrees (UNIT=DEG) [default: RAD:0.0,0.0,0.0]. The rotation is specified"
   << "\n                   in terms of yaw (around z-axis), pitch (around y) and roll (around x) angle: --rot-showers UNIT:yaw,pitch,roll"
   << "\n"
   << "\n -wgt      Weight the events accoirding to interaction probabilty to increase statistic [default: 0]"
   << "\n           Possible choices are: "
   << "\n             1)  -wgt 0 "
   << "\n               events are not weighted."
   << "\n             2)  -wgt 1 "
   << "\n               events are weighted according to a maximum probability."
#ifdef _KM3NET_ENABLED__
   << "\n             3)  -wgt 2 "
   << "\n               events are weighted because they are force to interact in the volume."
   << "\n"
   << "\n --cr-flux "
   << "\n           CR flux model needed to compute event weights for CORSIKA showers."
   << "\n           The default (and so far the only supported value) is 'GST3'."
#endif
#ifdef _HETAU_ENABLED__
   << "\n"
   << "\n -tau decay_mode"
   << "\n           Choose the decay mode of tau. Needed for HEDIS NUTAU-CC channel"
   << "\n           Possible choices are: "
   << "\n             1)  -tau 'muon' "
   << "\n               tau will decay into muon."
   << "\n             2)  -tau 'shower' "
   << "\n               tau will decay into shower."
   << "\n"
#endif
   << "\n"
   << "\v -h        Display this help and exit."
   << "\n"
   << "\n"
   << "\n **Environmental Parameters**"
   << "\n"
   << "\n GSeaGen environmental parameters: "
   << "\n"
   << "\n GSEAGEN:     pointing to the top directory of the code [mandatory]. "
   << "\n"
#ifdef _MUSIC_ENABLED__
   << "\n MUSICPATH:   pointing to the directory with music input files [default: MUSICPATH = "<<MUSICPATH<<"]. "
   << "\n"
#endif
   << "\n GMURNGFL:    pointing to file with max muon ranges [default: GMURNGFL = $GSEAGEN/dat/muon_rmax_music.dat]. "
   << "\n"
   << "\n MEDIACOMP:   overriding the default composition for target media: SeaWater, Rock, Mantle, Crust."
   << "\n                  Eg. MEDIACOMP = $GSEAGEN/dat/MediaComposition.xml "
   << "\n"
   << "\n PROPOSALCONF: pointing to the PROPOSAL configuration file. Default inputs are used if the env. parameter does not exist."
   << "\n                  Eg. PROPOSALCONF = proposaldir/resources/config.json"
   << "\n"
   << "\n SYSTLIST:    pointing to the file defining the systematic parameters. The calculation is activated if the env. parameter exists."
   << "\n                  Eg. SYSTLIST = $GSEAGEN/dat/SystParams.xml "
   << "\n"
   << "\n DEFPARAM:    pointing to the file defining the new default valus of input systematic parameters. Values will be changed with running options if used."
   << "\n                  Eg. SYSTLIST = $GSEAGEN/dat/SystParams.xml "
   << "\n"
   << "\n Based on GENIE, gSeaGen accepts all standard GENIE environmental parameters. For examples: "
   << "\n"
   << "\n GALGCONF:    overriding the default configuration folder [default: GALGCONF = $GENIE/config]"
   << "\n                  Eg. GALGCONF = $GENIE/config "
   << "\n              to use the default GENIE configuration folder. "
   << "\n"
   << "\n GXMLPATH:    to get user XML input files using GXMLPATH + default locations."
   << "\n"
   << "\n"
   << "\n"
   << " Please also read the GENIE documentation at http://www.genie-mc.org"
   << "\n";

}
#ifdef _ANTARES_ENABLED__
//___________________________________________________________________________
bool LoadGeometryDet(string filename)
{

  bool is_accessible = ! (gSystem->AccessPathName(filename.c_str()));

  if(is_accessible) {
    SLOG("gSeaNuEvGen", pINFO) << "Loading geometry file: "<< filename;
    ifstream geofile(filename.c_str(), ios::in);

    event geo;
    int ntags, ierr, id;
    double x,y,z,x_cg=0.,y_cg=0.,z_cg=0.;
    double radius=0.,r;
    double ROM=0.,ZOM=0.;
    double GroundZ;
    double OffSet = gGenPar.NAbsLength*gGenPar.AbsLength;
    double x0,y0,z0,z_max=0.,z_min=0.;
    vector<double> PosX,PosY,PosZ;
    string line;

    double MagneticDeclination;

    ierr = geo.read(geofile);

    if(ierr != 0){
      LOG("gSeaNuEvGen", pFATAL) << "Error "<<ierr<<" reading detector file";
      gAbortingInErr = true;
      exit(1);
      return false;
    }

    if (!geo.run_header()){
      LOG("gSeaNuEvGen", pFATAL) << "No run header found in the detector file";
      gAbortingInErr = true;
      exit(1);
      return false;
    }

    // reads the environment
    ntags=geo.ndat("environment");
    if(ntags==1){
      line=geo.next("environment");
      istringstream iline(line.c_str());
      iline >> gGenPar.SiteLatitude >> gGenPar.SiteLongitude >> MagneticDeclination>> gGenPar.SiteDepth>> GroundZ;
      gGenPar.SiteLongitude*=-1.;
    }
    else {
      if(ntags==0)LOG("gSeaNuEvGen", pFATAL) << " No environment defined in the detector file";
      if(ntags>1)LOG("gSeaNuEvGen", pFATAL) << "More than one environment defined in the detector file";
      gAbortingInErr = true;
      exit(1);
      return false;
    }

    // reads OM positions in the cluster
    ntags=geo.ndat("OM_position");
    if(ntags==0){
      LOG("gSeaNuEvGen", pFATAL) << " No OM_position defined in the detector file";
      gAbortingInErr = true;
      exit(1);
      return false;
    }
    else {
      for(int i=0; i<ntags; i++){
        line=geo.next("OM_position");
        istringstream iline(line.c_str());
        iline >> id >> x >> y >> z;
        r = sqrt(x*x+y*y);
        if(r>ROM)ROM=r;
        if(fabs(z)>ZOM)ZOM=fabs(z);
      }
    }

    ierr = geo.findev(geofile, 1);

    if(ierr != 0){
      LOG("gSeaNuEvGen", pFATAL) << "Error "<<ierr<<" reading detector file event";
      gAbortingInErr = true;
      exit(1);
      return false;
    }

    // reads OM_cluster positions
    ntags=geo.ndat("OM_cluster_data");
    if(ntags==0){
      LOG("gSeaNuEvGen", pFATAL) << " No OM_cluster_data defined in the detector file";
      gAbortingInErr = true;
      exit(1);
      return false;
    }
    else {
      for(int i=0; i<ntags; i++){
        line=geo.next("OM_cluster_data");
        istringstream iline(line.c_str());
        iline >> id >> x >> y >> z;
        PosX.push_back(x);
        PosY.push_back(y);
        PosZ.push_back(z);
      }
    }

    for(unsigned i=0; i<PosX.size(); i++){
      x_cg += PosX[i];
      y_cg += PosY[i];
      z_cg += PosZ[i];
    }

    x_cg /= PosX.size();
    y_cg /= PosX.size();
    z_cg /= PosX.size();

    for(unsigned i=0; i<PosX.size(); i++){
      x0=PosX[i]-x_cg;
      y0=PosY[i]-y_cg;
      z0=PosZ[i]-z_cg;
      r = sqrt(x0*x0+y0*y0);
      if(r>radius)radius=r;
      if(z0>z_max)z_max=z0;
      if(z0<z_min)z_min=z0;
    }

    gGenPar.Can1 = max(GroundZ-z_cg, z_min-ZOM-OffSet);
    gGenPar.Can2 = min(gGenPar.SiteDepth-z_cg, z_max+ZOM+OffSet);
    gGenPar.Can3 = radius+ROM+OffSet;

    // Can1 and Can2 are with respect the center of mass. We want them with respect the seabottom
    gGenPar.Can1 += z_cg;
    gGenPar.Can2 += z_cg;


    geofile.close();
  } else {
    LOG("gSeaNuEvGen", pFATAL)
       << "Specified detector geometry file [" << filename << "] is not accessible!";
    gAbortingInErr = true;
    exit(1);
    return false;
  }

  LOG("gSeaNuEvGen", pINFO)<< "Detector can defined from file";
  LOG("gSeaNuEvGen", pINFO)<< "Can: "<< gGenPar.Can1 <<" "<< gGenPar.Can2<< " " << gGenPar.Can3 <<" SiteDepth: "<<gGenPar.SiteDepth<<" Latitude: "<<gGenPar.SiteLatitude<<" Longitude: "<<gGenPar.SiteLongitude;

  return true;
}
//___________________________________________________________________________
bool LoadGeometryDetx(string filename)
{


  int global_det_id;
  string format_version;

  double UTC_validity_from,UTC_validity_to;

  string UTM_ref_grid,ZoneStr;
  double UTM_ref_easting,UTM_ref_northing,UTM_ref_z;

  int ndoms;
  int dom_id,line_id,floor_id,npmts;

  int pmt_id_global;
  double x, y, z, dx, dy, dz, t0;
  double x_cg=0.,y_cg=0.,z_cg=0.;
  double radius=0.,r=0.;
  double x0,y0,z0;
  double z_min=0.,z_max=0.;

  vector<double> PosX,PosY,PosZ;
  double GroundZ=0.;
  double OffSet = gGenPar.NAbsLength*gGenPar.AbsLength;

  bool is_accessible = ! (gSystem->AccessPathName(filename.c_str()));

  if(is_accessible) {
    SLOG("gSeaNuEvGen", pINFO) << "Loading geometry file: "<< filename;
    ifstream geofile(filename.c_str(), ios::in);

    //skip lines with comments at beginning of file
    string auxline;
    while ( std::getline(geofile, auxline) ) {
      if ( auxline.compare(0,1,"#")!=0 ) break;
    }
    istringstream iss(auxline);
    iss >> global_det_id;
    iss >> format_version;

    geofile>>UTC_validity_from>>UTC_validity_to;
    geofile>>UTM_ref_grid;

    string utm;
    geofile>>utm;
    UTM_ref_grid.append(" ");
    UTM_ref_grid.append(utm);
    geofile>>utm;
    UTM_ref_grid.append(" ");
    UTM_ref_grid.append(utm);


    ZoneStr=utm.substr(0,2);

    geofile>>UTM_ref_easting>>UTM_ref_northing>>UTM_ref_z;
    gGenPar.SiteDepth=fabs(UTM_ref_z);



    geofile>>ndoms;


    for(int idom=0; idom<ndoms; idom++){

      geofile>>dom_id>>line_id>>floor_id>>npmts;

      for(int ipmt=0; ipmt<npmts; ipmt++){
        geofile>>pmt_id_global>>x>>y>>z>>dx>>dy>>dz>>t0;
        PosX.push_back(x);
        PosY.push_back(y);
        PosZ.push_back(z);
      }
    }


    for(unsigned i=0; i<PosX.size(); i++){
      x_cg += PosX[i];
      y_cg += PosY[i];
      z_cg += PosZ[i];
    }

    x_cg /= PosX.size();
    y_cg /= PosX.size();
    z_cg /= PosX.size();

    for(unsigned i=0; i<PosX.size(); i++){
      x0=PosX[i]-x_cg;
      y0=PosY[i]-y_cg;
      z0=PosZ[i]-z_cg;
      r = sqrt(x0*x0+y0*y0);
      if(r>radius)radius=r;
      if(z0>z_max)z_max=z0;
      if(z0<z_min)z_min=z0;
    }

    gGenPar.Can1 = max(GroundZ-z_cg, z_min-OffSet);
    gGenPar.Can2 = min(gGenPar.SiteDepth-z_cg, z_max+OffSet);
    gGenPar.Can3 = radius+OffSet;


// Can1 and Can2 are with respect the center of mass. We want them with respect the seabottom

    gGenPar.Can1 += z_cg;
    gGenPar.Can2 += z_cg;


    geofile.close();
  } else {
    LOG("gSeaNuEvGen", pFATAL)
       << "Specified detector geometry file [" << filename << "] is not accessible!";
    gAbortingInErr = true;
    exit(1);
    return false;
  }

  double N=UTM_ref_northing+y_cg;
  double E=UTM_ref_easting+x_cg;

  int Zone;
  sscanf (ZoneStr.c_str(), "%d",&Zone);

  double N0=0.; // Northen Hemisphere
  double k0=0.9996;
  double E_0=500000.;//  m

  double a=6378.137*1E3;
  double f=1./298.257223563;
  double n=f/(2.-f);

  double A=a/(1+n)*(1+n*n/4.+n*n*n*n/64);

  double xi=(N-N0)/(k0*A);
  double eta=(E-E_0)/(k0*A);

  double beta_1=1./2.*n-2./3.*n*n+37./96.*n*n*n;
  double beta_2=1./48.*n*n-1./15.*n*n*n;
  double beta_3=17./480.*n*n*n;

  double delta_1=2*n-2./3.*n*n-2.*n*n*n;
  double delta_2=7./3.*n*n-8./5.*n*n*n;
  double delta_3=56./15.*n*n*n;

  double xip=xi-(beta_1*sin(2*1*xi)*cosh(2*1*eta)+beta_2*sin(2*2*xi)*cosh(2*2*eta)+beta_3*sin(2*3*xi)*cosh(2*3*eta));
  double etap=eta-(beta_1*cos(2*1*xi)*sinh(2*1*eta)+beta_2*cos(2*2*xi)*sinh(2*2*eta)+beta_3*cos(2*3*xi)*sinh(2*3*eta));
  double chi=asin(sin(xip)/cosh(etap));

  double phi=(chi+(delta_1*sin(2*1*chi)+delta_2*sin(2*2*chi)+delta_3*sin(2*3*chi)));
  double lambda0=((double)Zone*6.-183.)*kPi/180.;  // 6 e 183 in gradi
  double lambda=(lambda0+atan2(sinh(etap),cos(xi)));

  gGenPar.SiteLatitude=phi;
  gGenPar.SiteLongitude=lambda;

  LOG("gSeaNuEvGen", pINFO)<< "Detector can defined from file";
  LOG("gSeaNuEvGen", pINFO)<< "Can: "<< gGenPar.Can1 <<" "<< gGenPar.Can2<< " " << gGenPar.Can3 <<" SiteDepth: "<<gGenPar.SiteDepth<<" Latitude: "<<gGenPar.SiteLatitude<<" Longitude: "<<gGenPar.SiteLongitude;

  return true;
}
#endif
//___________________________________________________________________________
void CalcUTMConvergeAngle(void){

  // Calculate the UTM convergence angle

  double zone=1+int( 60.* (  gGenPar.SiteLongitude*180./kPi +180.)/360.);
  double lambda0=(-180.+(zone-1.)*6.+3)*kPi/180.; // in rad
//  gGenPar.UTMConvergeAngle=atan(tan(gGenPar.SiteLongitude-lambda0)*sin(gGenPar.SiteLatitude));

  double lambda  = gGenPar.SiteLongitude;
  double phi     = gGenPar.SiteLatitude;
  double omega   = lambda-lambda0;

  // parameters of the Earth ellipsoid
  double a  = 6378137.;      // semi-major axis in meters (WGS84)
  double e2 = 0.0066943800;  // eccentricity (WGS84)

  double rho = a*(1.-e2)/pow(1.-e2*pow(sin(phi),2),3./2.);
  double nu  = a/sqrt(1-e2*pow(sin(phi),2.));
  double psi = nu/rho;
  double t = tan(phi);

  // power series for convergence angle
  gGenPar.UTMConvergeAngle = sin(phi) * omega
    - sin(phi) * pow(omega,3.)/3.  * pow(cos(phi),2.) * (2.*pow(psi,2.)-psi)
    - sin(phi) * pow(omega,5.)/15. * pow(cos(phi),4.) * (pow(psi,4.)*(11.-24.*pow(t,2.))-
							 pow(psi,3.)*(11.-36.*pow(t,2.))+
							 2.*pow(psi,2.)*(1.-7.*pow(t,2.))+
							 psi*pow(t,2.)
							 )
    - sin(phi) * pow(omega,7.)/315. * pow(cos(phi),6.) * (17.-26.*pow(t,2.)+2.*pow(t,4.));

  LOG("gSeaNuEvGen", pINFO)<< "UTM Convergence Angle: "<< gGenPar.UTMConvergeAngle*180./kPi << " deg";

  return;
}
//___________________________________________________________________________
void DecodeCommandLine(int argc, char * argv[]){

  // Get the command line arguments
  LOG("gSeaNuEvGen", pNOTICE) << "Parsing command line arguments";

  bool help=false;

  // Decode command line
  vector<string> VecArg;

  size_t apos;
  string opt,value;

  bool defcan=false; // true if can set by hand
  bool setcoord=false; // true if detector coordinates set by hand
  double depth=0., latitude=0., longitude=0., declination=gGenPar.SiteDeclination;

  bool DefDepth=false;

  bool CheckFluxOpt=false;

  bool PointFile=false;

  bool GenEl  = false;
  bool GenMu  = false;
  bool GenTau = false;
  bool GenCan = false;

#ifdef _HETAU_ENABLED__
  string TauDecCh = "";
#endif

  string cantype;

  for (int i = 1; i < argc; i++)
  {
    string arg(argv[i]);
    if (arg.compare(0,1,"-")==0){
      apos=arg.find(" ");
      opt=arg.substr(0,apos);

      if(opt.compare("-r")==0){         // run number
        i++;
        LOG("gSeaNuEvGen", pDEBUG) << "Reading MC run number";
        gGenPar.RunNu = atoi(argv[i]);
      }
      if(opt.compare("-w")==0){         // write the genie native output
        LOG("gSeaNuEvGen", pDEBUG) << "Write the generator native output";
        gGenPar.NativeOutput = true;
      }
#if defined _ANTARES_ENABLED__  || defined _KM3NET_ENABLED__
      if(opt.compare("-write")==0){
        i++;
        LOG("gSeaNuEvGen",pDEBUG) << "Write particle mode";
        gGenPar.WriteParticlesMode = atoi(argv[i]);
      }
#endif
      if(opt.compare("-save")==0){
        i++;
        LOG("gSeaNuEvGen",pDEBUG) << "Choosing the sea-level particles to be saved";

        if ( (gGenPar.WriteParticlesMode!=2) && (((string) argv[i]).compare("mu")!=0) ) {
          LOG("gSeaNuEvGen",pFATAL) << "For -save options other than 'mu', -write 2 must be used, exiting";
          gAbortingInErr = true;
          exit(1);
        }
        if(((string) argv[i]).compare("mu")==0) gGenPar.SaveMu = true;
        else if(((string) argv[i]).compare("nu")==0) {gGenPar.SaveMu = false; gGenPar.SaveNu = true;}
        else if(((string) argv[i]).compare("lep")==0) {gGenPar.SaveMu = true; gGenPar.SaveNu = true;}
        else if(((string) argv[i]).compare("all")==0) gGenPar.SaveOther = true;
        else {
          LOG("gSeaNuEvGen",pFATAL) << "Invalid particles specified for the -save option, exiting";
          gAbortingInErr = true;
          exit(1);
        }
      }
      if(opt.compare("-chances")==0){
        i++;
        LOG("gSeaNuEvGen",pDEBUG) << "Number of extra chances for CORSIKA showers";
        gGenPar.ChancesPerShower = atoi(argv[i]);
      }
      if(opt.compare("--rot-showers")==0){
        i++;
        LOG("gSeaNuEvGen", pINFO) << "Reading the angles by which the CORSIKA showers at sea level should be rotated";
        string rot_str = argv[i];
        string::size_type column = rot_str.find_first_of(":",0);
        if(column==string::npos) { // check if units are specified
          LOG("gSeaNuEvGen",pFATAL) << "You need to specify the rotation units (DEG or RAD)!";
          gAbortingInErr = true;
          exit(1);
        }
        string rot_unit = rot_str.substr(0,column);
        rot_str.erase(0,column+1);

        string::size_type comma = rot_str.find_first_of(",",0);
        if(comma==string::npos){ // check if the inputs are comma-separated
          LOG("gSeaNuEvGen",pFATAL) << "Error in detector coordinates input ";
          gAbortingInErr = true;
          exit(1);
        }
        vector<string> coords = utils::str::Split(rot_str, ",");
        assert( coords.size() == 3 ); // 3 arguments must be given
        for (int i=0; i<3; i++) gGenPar.CustomRotation[i] = atof(coords[i].c_str());
        if(rot_unit.compare("DEG")==0) for (int i=0; i<3; i++) gGenPar.CustomRotation[i] *= kPi/180.;
        else if(rot_unit.compare("RAD")!=0){ // check if a valid unit was used
          LOG("gSeaNuEvGen",pFATAL) << "Error in rotation angle units, they must be DEG or RAD!";
          gAbortingInErr = true;
          exit(1);
        }
      }
      if(opt.compare("-wgt")==0){
        i++;
        LOG("gSeaNuEvGen",pDEBUG) << "Weight events";
        gGenPar.WeightOpt = atoi(argv[i]);
      }
      if(opt.compare("--cr-flux")==0){
        i++;
        gGenPar.CRModel = (string) argv[i];
        LOG("gSeaNuEvGen",pDEBUG) << "Use " << gGenPar.CRModel << " CR flux model";
      }
      if(opt.compare("-o")==0){         // event file prefix
        i++;
        LOG("gSeaNuEvGen", pDEBUG) << "Reading the event filename prefix";
        gGenPar.EvFilePrefix = (string) argv[i];
      }
      if(opt.compare("-a")==0){         // generation spectral index
        i++;
        LOG("gSeaNuEvGen", pDEBUG) << "Reading the generation spectral index";
        gGenPar.Alpha = atof(argv[i]);
      }
      if(opt.compare("-b")==0){         // number of energy bins
        i++;
        LOG("gSeaNuEvGen", pDEBUG) << "Reading the number of energy bins";
        gGenPar.NBin = atoi(argv[i]);
      }
#ifdef _ANTARES_ENABLED__
      if(opt.compare("-l")==0){         // light abs lenght for the can definition
        i++;
        LOG("gSeaNuEvGen", pINFO) << "Reading the light absorption lenght";
        string sla = argv[i];
        // may be a comma separated set of values
        if(sla.find(",") != string::npos) {
          // split the comma separated list
          vector<string> vla = utils::str::Split(sla, ",");
          assert(vla.size() == 2);
          gGenPar.AbsLength  = atof(vla[0].c_str());
          gGenPar.NAbsLength = atof(vla[1].c_str());
        }
        else{
          gGenPar.AbsLength = atof(argv[i]);
        }
        cout<<gGenPar.AbsLength<<" "<<gGenPar.NAbsLength<<endl;
      }
#endif
      if(opt.compare("--seed")==0){             // random generation seed
        i++;
        LOG("gSeaNuEvGen", pDEBUG) << "Reading random number seed";
        gGenPar.RanSeed = atoi(argv[i]);
      }
      if(opt.compare("-point")==0){             // point-mode
        i++;
        LOG("gSeaNuEvGen", pINFO) << "Point-mode activated";
        gGenPar.GenMode="POINT";
        gGenPar.SourceName="POINT";
        string source_str = argv[i];

        string::size_type column = source_str.find_first_of(":",0);
        if(column!=string::npos) {

          if(source_str.find("DEC")!=string::npos){
            source_str.erase(0,column+1);
            gGenPar.Declination=atof(source_str.c_str());
            if(fabs(gGenPar.Declination)>90.){
              LOG("gSeaNuEvGen",pFATAL) << "Source declination must be in the range [-90.,+90.] degrees, exiting...";
              gAbortingInErr = true;
              exit(1);
            }
          }
          else if(source_str.find("FILE")!=string::npos){
            source_str.erase(0,column+1);
            gGenPar.SourceFile=source_str;
            struct stat buffer;
            if(stat (gGenPar.SourceFile.c_str(), &buffer) != 0){
              LOG("gSeaNuEvGen",pFATAL) << "Source info file does not exits, exiting...";
              gAbortingInErr = true;
              exit(1);
            }
            PointFile=true;
          }
        }
        else {
          LOG("gSeaNuEvGen",pFATAL) << "Error in point-mode source info, exiting...";
          gAbortingInErr = true;
          exit(1);
        }
      }
      if(opt.compare("--cross-sections")==0){   // input cross-section file
        i++;
        LOG("gSeaNuEvGen", pINFO) << "Reading cross-section file";
        gGenPar.InpXSecFile = argv[i];
      }
      if(opt.compare("--event-generator-list")==0){     // input event generator list N.B. MUST BE BEFORE PARSING OPT -F!!!!
        i++;
        LOG("gSeaNuEvGen", pINFO) << "Reading event-generator-list";
        gGenPar.EventGeneratorList = argv[i];
      }
      if(opt.compare("--tune")==0){     // input event generator tune!!!!
        i++;
        LOG("gSeaNuEvGen", pINFO) << "Reading event-generator-tune";
        gGenPar.GenTune = argv[i];
      }
      if(opt.compare("-PDG")==0){               // use PDG code
        gOptPDG = true;
      }
      if(opt.compare("-p")==0){
        gGenPar.PropMode=true;
        LOG("gSeaNyGen", pINFO) << "Propagation through the Earth activated";
      }
      if(opt.compare("-n")==0){         //  number of neutrinos to generate
        i++;
        LOG("gSeaNuEvGen", pDEBUG) << "Reading number of events to generate";
        gGenPar.NTot = atof(argv[i]);
      }
      if(opt.compare("-e")==0){         // neutrino energy range
        i++;
        LOG("gSeaNuEvGen", pINFO) << "Reading neutrino energy range";
        string nue = argv[i];
        // must be a comma separated set of values
        if(nue.find(",") != string::npos) {
          // split the comma separated list
          vector<string> nurange = utils::str::Split(nue, ",");
          assert(nurange.size() == 2);
          double emin = atof(nurange[0].c_str());
          double emax = atof(nurange[1].c_str());
          if(emin==emax)emax+=0.001;
          assert(emax>emin && emin>=0.);
          gGenPar.EvMin = emin;
          gGenPar.EvMax = emax;
        }
      }
      if(opt.compare("-t")==0){         // neutrino cosTheta range
        i++;
        LOG("gSeaNuEvGen", pINFO) << "Reading neutrino cosTheta range";
        string nuct = argv[i];
        // must be a comma separated set of values
        if(nuct.find(",") != string::npos) {
          // split the comma separated list
          vector<string> nuctrange = utils::str::Split(nuct, ",");
          assert(nuctrange.size() == 2);
          double ctmin = atof(nuctrange[0].c_str());
          double ctmax = atof(nuctrange[1].c_str());
          assert(ctmax>=ctmin && ctmin>=-1. && ctmax<=1.);
          gGenPar.CtMin = ctmin;
          gGenPar.CtMax = ctmax;
        }
      }
      // can
      if(opt.compare("-c")==0){         // detector can
        i++;
        string incan = argv[i];
        // get input source (DET, DETX or CAN)
        string::size_type canend = incan.find_first_of(":",0);
        if(canend==string::npos) {
#ifdef _ANTARES_ENABLED__
          LOG("gSeaNuEvGen",pFATAL) << "You need to specify the can input type: CAN, DET or DETX";
#endif
#ifndef _ANTARES_ENABLED__
          LOG("gSeaNuEvGen",pFATAL) << "You need to specify the can input type: CAN";
#endif
          gAbortingInErr = true;
          exit(1);
        }
        cantype = incan.substr(0,canend);
        incan.erase(0,canend+1);
        if(cantype.compare("CAN")==0) {
          defcan=true;
          LOG("gSeaNuEvGen",pINFO)<< "Reading detector can size";
          if(incan.find(",") != string::npos){
            // split the comma separated list
            vector<string> cansize = utils::str::Split(incan, ",");
            assert(cansize.size() == 3);
            gGenPar.Can1 = atof(cansize[0].c_str());
            gGenPar.Can2 = atof(cansize[1].c_str());
            gGenPar.Can3 = atof(cansize[2].c_str());
            LOG("gSeaNuEvGen", pINFO)<< "Can: "<< gGenPar.Can1 <<" "<< gGenPar.Can2<< " " << gGenPar.Can3;
          }
          else {
            LOG("gSeaNuEvGen", pFATAL)<< "Invalid can size. Use `-c can1,can2,can3', eg. `-c -278.15,341.47,266.11'";
            gAbortingInErr = true;
            exit(1);
          }
        }
#ifdef _ANTARES_ENABLED__
        else if(cantype.compare("DET")==0) gGenPar.detector = incan; // the can size will be calculated after option reading
#endif
        else if(cantype.compare("DETX")==0) gGenPar.detector = incan; // the can size will be calculated after option reading
        else if (true){
#ifdef _ANTARES_ENABLED__
          LOG("gSeaNuEvGen",pFATAL) << "You made an error in specifying the can input type: choose CAN or DET";
#endif
#ifndef _ANTARES_ENABLED__
          LOG("gSeaNuEvGen",pFATAL) << "You made an error in specifying the can input type: choose CAN";
#endif
          gAbortingInErr = true;
          exit(1);
        }
        // setting site declination by hand if Needed
/*
        if(opt.compare("--site-dec")==0){
          i++;
          declination = atof(argv[i]);
          setcoord=true;
          if(depth<=0.){
            LOG("gSeaNuEvGen", pFATAL)<< "The site depth must be a positive number', eg. `-d 3500'";
            gAbortingInErr = true;
            exit(1);
          }
        }
*/
        if((cantype.compare("DET")==0) || (cantype.compare("DETX")==0)){
          // setting site declination by hand if Needed
          if(opt.compare("--site-dec")==0){
            i++;
            declination = atof(argv[i]);
            assert(abs(declination)<=0.5*kPi);
            setcoord=true;
          }
        }

      } // can
      else if(opt.compare("--coord")==0){      // detector coordinates (if the site name was not specified)
        i++;
        setcoord=true;
        string coord_str = argv[i];
        string::size_type column = coord_str.find_first_of(":",0);
        if(column==string::npos) {
          LOG("gSeaNuEvGen",pFATAL) << "You need to specify the detector coordinate units (DEG or RAD) or site name ";
          gAbortingInErr = true;
          exit(1);
        }
        string coord_unit = coord_str.substr(0,column);
        coord_str.erase(0,column+1);

        if(coord_unit.compare("SITE")==0){ // Experimental site (predefined location)
          LOG("gSeaNuEvGen", pINFO) << "Reading site location (by name)";

          // all declinations are according to WMM-2020 on 12.01.2021 (https://www.ngdc.noaa.gov/geomag/calculators/magcalc.shtml)
          if (coord_str.compare("ARCA")==0){ // KM3NeT-ARCA site (Mediterranean Sea, close to Sicily)
            latitude  = 36.266667*kPi/180.;
            longitude = 16.100000*kPi/180.;
            declination = 3.69*kPi/180.; // ± 0.33°  changing by  0.11° E per year
          }
          else if (coord_str.compare("ORCA")==0){ // KM3NeT-ORCA site (Mediterranean Sea, close to Toulon)
            latitude  = 42.800000*kPi/180.;
            longitude = 6.033333*kPi/180.;
            declination = 2.28*kPi/180.; // ± 0.35°  changing by  0.15° E per year
          }
          else if (coord_str.compare("IceCube")==0){ // IceCube site (South Pole)
            latitude  = 89.99*kPi/180.;
            longitude = -63.453056*kPi/180.;
            declination = 29.01*kPi/180.; //  ± 0.41°  changing by  0.11° W per year
          }
          else if (coord_str.compare("GVD")==0){ // GVD site (Lake Baikal)
            latitude  = 51.7763311*kPi/180.;
            longitude = 104.3736644*kPi/180.;
            declination = -4.35*kPi/180.; // ± 0.40°  changing by  0.15° W per year
          }
          else {
            LOG("gSeaNuEvGen",pFATAL) << "You must specify a valid detector site";
            gAbortingInErr = true;
            exit(1);
          }
        }
        else{
          string::size_type comma = coord_str.find_first_of(",",0);
          if(comma==string::npos){
            LOG("gSeaNuEvGen",pFATAL) << "Error in detector coordinates input ";
            gAbortingInErr = true;
            exit(1);
          }
          vector<string> coords = utils::str::Split(coord_str, ",");
          assert((coords.size() == 2) || (coords.size() == 3)); // 2 or 3 arguments
          latitude  = atof(coords[0].c_str());
          longitude = atof(coords[1].c_str());
          if (coords.size()==3) declination = atof(coords[2].c_str());
          if(coord_unit.compare("DEG")==0){
            latitude  *= kPi/180.;
            longitude *= kPi/180.;
            if (coords.size()==3) declination*= kPi/180.;
          }
          else if(coord_unit.compare("RAD")!=0){
            LOG("gSeaNuEvGen",pFATAL) << "Error in detector coordinate units; they must be DEG or RAD (or you must give site name with SITE) ";
            gAbortingInErr = true;
            exit(1);
          }
        }
      }
      if ((opt.compare("EVT")>0 || opt.compare("BIN")>0) && !setcoord) {
        LOG("gSeaNuEvGen",pFATAL) << "Error: detector declination MUST be specified for CORSIKA input (use --coord)";
        gAbortingInErr = true;
        exit(1);
      }
      if(opt.compare("-origin")==0){            // origin position in the detector rest frame
        i++;
        LOG("gSeaNuEvGen",pINFO)<< "Reading origin position in the detector rest frame";
        string origin = argv[i];

        if(origin.find(",") != string::npos){
            // split the comma separated list
            vector<string> coord = utils::str::Split(origin, ",");
            if(coord.size()==3){
              gGenPar.OriginX = atof(coord[0].c_str());
              gGenPar.OriginY = atof(coord[1].c_str());
              gGenPar.OriginZ = atof(coord[2].c_str());
              LOG("gSeaNuEvGen", pINFO)<< "Origin Position: "<< gGenPar.OriginX <<" "<< gGenPar.OriginY<< " " << gGenPar.OriginZ;
            }
            else {
              LOG("gSeaNuEvGen", pFATAL)<< "Invalid Origin Position. Use `-origin X0,Y0,Z0'";
              gAbortingInErr = true;
              exit(1);
            }
          }
          else {
            LOG("gSeaNuEvGen", pFATAL)<< "Invalid Origin Position. Use `-origin X0,Y0,Z0'";
            gAbortingInErr = true;
            exit(1);
          }
      }
      if(opt.compare("-bedrock")==0){           // set the bedrock height
        i++;
        LOG("gSeaNuEvGen", pDEBUG) << "Reading the bedrock height";
        gGenPar.HBedRock = atof(argv[i]);
      }


      if(opt.compare("-g")==0){         // set the generator code
        i++;

        gGenPar.GenCode = argv[i];
        gGenPar.GenCode = genie::utils::str::ToLower(gGenPar.GenCode);

        if(gGenPar.GenCode.compare("genie")==0)gGenPar.GenCode = "GENIE";

        if(gGenPar.GenCode.compare("GENIE")==0)LOG("gSeaNuEvGen", pNOTICE) << "Neutrino interaction simulated with GENIE";

        else {
          LOG("gSeaNuEvGen", pFATAL) << "Unknown neutrino interaction code: "<<gGenPar.GenCode<<", exiting...";
          gAbortingInErr = true;
          exit(1);
        }
      }

      if(opt.compare("-prop")==0){              // set the propagation code
        i++;

        string PropCode = argv[i];
        PropCode = genie::utils::str::ToLower(PropCode);

        if(PropCode.compare("propamuon")==0)gGenPar.PropCode = "PropaMuon";
        if(PropCode.compare("music")==0)gGenPar.PropCode = "MUSIC";
        if(PropCode.compare("proposal")==0)gGenPar.PropCode = "PROPOSAL";
        if(PropCode.compare("jpp")==0)gGenPar.PropCode= "JPP";
        if(gGenPar.PropCode.compare("PropaMuon")==0)LOG("gSeaNuEvGen", pNOTICE) << "Muon propagation with PropaMuon activated";  
#ifdef _MUSIC_ENABLED__
        else if(gGenPar.PropCode.compare("MUSIC")==0){
          gGenPar.PropCodeVer="2004";
	  LOG("gSeaNuEvGen", pNOTICE)     << "Muon propagation with MUSIC activated";
	}
#endif
#ifdef _PROPOSAL_ENABLED__
        else if(gGenPar.PropCode.compare("PROPOSAL")==0){
	  char ch[100];
	  sprintf(ch, "%s.%s.%s", STR(_PROPOSAL_VERSION_MAJOR__),STR(_PROPOSAL_VERSION_MINOR__),STR(_PROPOSAL_VERSION_PATCH__));
          gGenPar.PropCodeVer=ch;
	  LOG("gSeaNuEvGen", pNOTICE)     << "Muon propagation with PROPOSAL "<<gGenPar.PropCodeVer<<" activated";
	}
#endif
#ifdef _JPP_ENABLED__
        else if(gGenPar.PropCode.compare("JPP")==0)LOG("gSeaNuEvGen", pNOTICE)     << "Muon propagation with JPP activated";
#endif
        else {
          LOG("gSeaNuEvGen", pFATAL) << "Unknown Muon propagation code: "<<gGenPar.PropCode<<", exiting...";
          gAbortingInErr = true;
          exit(1);
        }
      }

      if(opt.compare("-d")==0){
        i++;
        depth = atof(argv[i]);
        DefDepth=true;
        if(depth<=0.){
          LOG("gSeaNuEvGen", pFATAL)<< "The site depth must be a positive number', eg. `-d 3500'";
          gAbortingInErr = true;
          exit(1);
        }
      }
      if(opt.compare("-C")==0){         // use the can as interaction volume
        LOG("gSeaNuEvGen",pINFO) << "Setting the detector can as interation volume";
        GenCan = true;
      }
#ifdef _ANTARES_ENABLED__
      if(opt.compare("-rt")==0){
        i++;
        string rtopt = argv[i];
        gGenPar.RTOpt = genie::utils::str::ToLower(rtopt);
        if(gGenPar.RTOpt.compare("can")==0){
          LOG("gSeaNuEvGen", pNOTICE) << "RTOpt = 'can': Generation surface will cover the detector can";
        }
        else if(gGenPar.RTOpt.compare("vol")==0){
          LOG("gSeaNuEvGen", pNOTICE) << "RTOpt = 'vol': Generation surface will cover the interaction volume";
        }
        else if(gGenPar.RTOpt.compare("proj")==0){
          LOG("gSeaNuEvGen", pNOTICE) << "RTOpt = 'proj': Generation surface will cover the projected detector can";
        }
        else {
          double RTValue=atof(gGenPar.RTOpt.c_str());
          ostringstream ss;
          ss << RTValue;
          string sch(ss.str());
          if(gGenPar.RTOpt.size()!=sch.size() || RTValue<=0){
            LOG("gSeaNuEvGen", pNOTICE) << "Error in setting RT: it must be can, vol or a positive number";
            gAbortingInErr = true;
            exit(1);
          }
          else {
            LOG("gSeaNuEvGen", pNOTICE) << "Generation surface radius set by hand: RT = "<<RTValue<<" m";
          }
        }
      }


      if(opt.compare("-R")==0){         // Coordinate rotation matrix
        i++;
        string rotarg = argv[i];
        //get convention
        string::size_type j = rotarg.find_first_of(":",0);
        string convention = "";
        if(j==string::npos) { convention = "X"; }
        else                { convention = rotarg.substr(0,j); }
        //get angles phi,theta,psi
        rotarg.erase(0,j+1);
        vector<string> euler_angles = utils::str::Split(rotarg,",");
        if(euler_angles.size() != 3) {
          LOG("gSeaNuEvGen", pFATAL) << "You didn't specify all 3 Euler angles using the -R option";
          gAbortingInErr = true;
          exit(1);
        }
        double phi   = atof(euler_angles[0].c_str());
        double theta = atof(euler_angles[1].c_str());
        double psi   = atof(euler_angles[2].c_str());
        //set Euler angles using appropriate convention
        if(convention.find("X")!=string::npos || convention.find("x")!=string::npos) {
          LOG("gSeaNuEvGen", pNOTICE) << "Using X-convention for input Euler angles";
          gGenPar.Rotation.SetXEulerAngles(phi,theta,psi);
        } else
        if(convention.find("Y")!=string::npos || convention.find("y")!=string::npos) {
          LOG("gSeaNuEvGen", pNOTICE) << "Using Y-convention for input Euler angles";
          gGenPar.Rotation.SetYEulerAngles(phi,theta,psi);
        } else {
          LOG("gSeaNuEvGen", pFATAL) << "Unknown Euler angle convention. Please use the X- or Y-convention";
          gAbortingInErr = true;
          exit(1);
        }
        //invert?
        if(convention.find("^-1")!=string::npos) {
          LOG("gSeaNuEvGen", pNOTICE) << "Inverting rotation matrix";
          gGenPar.Rotation.Invert();
        }
      }
#endif


      if(opt.compare("-time")==0){         // Generation interval time
        i++;

        string timeformat="TS";
        string intime = argv[i];

        // get input time format
        string::size_type timeend = intime.find_first_of(":",0);
        if(timeend!=string::npos) {
          timeformat = intime.substr(0,timeend);
          timeformat = genie::utils::str::ToUpper(timeformat);
          intime.erase(0,timeend+1);
        }

        vector<string> timeinterval = utils::str::Split(intime, ",");

        if(timeinterval.size() != 2){
            LOG("gSeaNuEvGen", pFATAL)<< "Error in -time option. Use `-time TimeFormat:TimeStart,TimeStop'";
            gAbortingInErr = true;
            exit(1);
        }

        if(timeformat.compare("TS")==0){
          gGenPar.TimeStampStart = atoi(timeinterval[0].c_str());
          gGenPar.TimeStampStop = atoi(timeinterval[1].c_str());
          gGenPar.MJDStart = GAstro::CalcMJDFromTimeStamp(gGenPar.TimeStampStart,0);
          gGenPar.MJDStop = GAstro::CalcMJDFromTimeStamp(gGenPar.TimeStampStop,0);
        }
        else if(timeformat.compare("MJD")==0 || timeformat.compare("DATE")==0){
          if(timeformat.compare("MJD")==0){
            gGenPar.MJDStart = atof(timeinterval[0].c_str());
            gGenPar.MJDStop = atof(timeinterval[1].c_str());
          }
          else {
            gGenPar.MJDStart = GAstro::CalcMJD(timeinterval[0].c_str());
            gGenPar.MJDStop = GAstro::CalcMJD(timeinterval[1].c_str());
          }
          unsigned int TimeStampSec, TimeStampTick;
          GAstro::CalcTimeStamp(gGenPar.MJDStart, &TimeStampSec, &TimeStampTick);
          gGenPar.TimeStampStart = TimeStampSec;
          GAstro::CalcTimeStamp(gGenPar.MJDStop, &TimeStampSec, &TimeStampTick);
          gGenPar.TimeStampStop = TimeStampSec;
        }
        else {
          LOG("gSeaNuEvGen", pFATAL) << "Error in -time option. Unknown time format: "<<timeformat<<", exiting...";
          gAbortingInErr = true;
          exit(1);
        }

        gGenPar.GenMJD=true;
//        gGenPar.TGen=(gGenPar.MJDStop-gGenPar.MJDStart)*86400.; // otherwise default value for TGen is used

        string start,stop;
        GAstro::CalcDate(gGenPar.MJDStart, &start);
        GAstro::CalcDate(gGenPar.MJDStop, &stop);

        LOG("gSeaNuEvGen", pINFO) << "Generation interval time: "<<start<<" - "<<stop;

      }


      FluxFileStr fluxstr;

      if(opt.compare("-f")==0){         // neutrinos fluxes and types
        i++;
        CheckFluxOpt=true;
        LOG("gSeaNuEvGen", pDEBUG) << "Getting input flux files";
        string fluxtot = argv[i];


        if(fluxtot.find("[")==0){       // case -f [nu,nub] --> simulate with the default FUNC
          string neutrinos=fluxtot;
          fluxtot="FUNC1:1E-5*pow(x,-2)";
          fluxtot.append(neutrinos);
        }

        // find FluxSim
        vector<string::size_type> nColon;
        vector<string::size_type> nSemiColon;
        string::size_type pos=-1;
        while(1){
          pos=fluxtot.find(":",pos+1);
          if(pos==string::npos)break;
          nColon.push_back(pos);
        }
        pos=-1;
        while(1){
          pos=fluxtot.find(";",pos+1);
          if(pos==string::npos)break;
          nSemiColon.push_back(pos);
        }
        nSemiColon.push_back(string::npos);
        if(nSemiColon.size()!=nColon.size()){
          LOG("gSeaNuEvGen", pFATAL) << "You made an error in specifying the flux info";
          gAbortingInErr = true;
          exit(1);
        }

        pos=0;
        for(unsigned iFluxSim=0; iFluxSim<nColon.size(); iFluxSim++){
        // get flux simulation info (FLUKA,BGLRS,...) so as to instantiate the
        // appropriate flux driver
          fluxstr.FluxSimul = fluxtot.substr(pos,nColon[iFluxSim]-pos);
          for(string::size_type iFlux=0; iFlux<fluxstr.FluxSimul.size(); iFlux++) {
            fluxstr.FluxSimul[iFlux] = toupper(fluxstr.FluxSimul[iFlux]);
          }

          if(fluxstr.FluxSimul.compare("BGLRS")==0)fluxstr.FluxSimul = "BARTOL";

          if(fluxstr.FluxSimul.compare("EVT")==0 || fluxstr.FluxSimul.compare("BIN")==0){
            if(gGenPar.GenMode.compare("POINT")==0){
              LOG("gSeaNuEvGen", pINFO)<< "POINT mode not compatible with reading incident neutrinos from file, exiting...";
              LOG("gSeaNuEvGen", pINFO)<< "Remove -point option or modify -f option";
              exit(1);
            }
            if     (fluxstr.FluxSimul.compare("EVT")==0) gGenPar.GenMode = "EVT";
            else if(fluxstr.FluxSimul.compare("BIN")==0) gGenPar.GenMode = "BIN";
          }

          string flux = fluxtot.substr(nColon[iFluxSim]+1,nSemiColon[iFluxSim]-nColon[iFluxSim]-1);

          // now get the list of input files and the corresponding neutrino codes.

          if(gGenPar.GenMode=="BIN" || gGenPar.GenMode=="EVT"){

              string::size_type open_bracket   = flux.find("[",0);
              if(open_bracket!=string::npos)flux = flux.substr(0,open_bracket);
              fluxstr.NuType= 0;
              fluxstr.FileNames.push_back(flux);

              //we cannot propagate neutrinos yet from CORSIKA or CORANT files
              /*gGenPar.NuList.push_back(12);
              gGenPar.NuList.push_back(-12);
              gGenPar.NuList.push_back(14);
              gGenPar.NuList.push_back(-14);
              gGenPar.NuList.push_back(16);
              gGenPar.NuList.push_back(-16);*/

              gGenPar.FluxFiles.push_back(fluxstr);

          }
          else {

            while(1){
              string::size_type open_bracket   = flux.find("[",0);
              string::size_type close_bracket  = flux.find("]",0);
              if (open_bracket ==string::npos || close_bracket==string::npos){
                LOG("gSeaNuEvGen", pFATAL) << "You made an error in specifying the flux info";
                gAbortingInErr = true;
                exit(1);
              }

              string neutrino_pdg = flux.substr(open_bracket+1,close_bracket-open_bracket-1);
              vector<int> vec_neutrino_pdg;
              int pdg_id;
              string::size_type comnu=neutrino_pdg.find(",");
              if(comnu==string::npos){
                pdg_id=atoi(neutrino_pdg.c_str());
                vec_neutrino_pdg.push_back(pdg_id);
              }
              else {
                string nu_pdg1,nu_pdg2;
                nu_pdg1=neutrino_pdg.substr(0,comnu);
                nu_pdg2=neutrino_pdg;
                nu_pdg2.erase(0,comnu+1);
                pdg_id=atoi(nu_pdg1.c_str());
                vec_neutrino_pdg.push_back(pdg_id);
                pdg_id=atoi(nu_pdg2.c_str());
                vec_neutrino_pdg.push_back(pdg_id);
              }

              for(unsigned iSize=0; iSize<vec_neutrino_pdg.size(); iSize++){
                fluxstr.NuType = vec_neutrino_pdg[iSize];

                string filelist = flux.substr(0,open_bracket);
                string::size_type comma;
                fluxstr.FileNames.clear();
                if(fluxstr.FluxSimul.find("FUNC")!=std::string::npos){ // commas are part of the math formula
                  if(fluxstr.FluxSimul=="FUNC")fluxstr.FluxSimul="FUNC1"; // copatibility with old versions
                  fluxstr.FileNames.push_back(filelist);

                }
                else {
                   // commas separate files for the same neutrino types
                  while(1){
                    comma=filelist.find(",",0);
                    string filename = filelist.substr(0,comma);
                    filelist.erase(0,comma+1);
                    fluxstr.FileNames.push_back(filename);
                    if(filename.size()==0){
                      LOG("gSeaNuEvGen", pFATAL) << "No neutrino flux file or function is defined";
                      gAbortingInErr = true;
                      exit(1);
                    }
                    if(comma==string::npos)break;
                  }
                }

                if(abs(fluxstr.NuType)==12)GenEl=true;
                if(abs(fluxstr.NuType)==14)GenMu=true;
                if(abs(fluxstr.NuType)==16)GenTau=true;

                gGenPar.NuList.push_back(fluxstr.NuType);
                gGenPar.FluxFiles.push_back(fluxstr);
              }

              flux.erase(0,close_bracket+1);
              if(flux.compare(0,1,",")==0)flux.erase(0,1);
              else break;

            }
            pos=nSemiColon[iFluxSim]+1;
          }
          sort(gGenPar.NuList.begin(),gGenPar.NuList.end());
          for(unsigned int iList=1; iList<gGenPar.NuList.size(); iList++){
            if(gGenPar.NuList[iList]==gGenPar.NuList[iList-1])gGenPar.NuList.erase(gGenPar.NuList.begin()+iList);
          }
        }

      } // end -f option

      if(opt.compare("-utm")==0){         // use UTM detector system
        LOG("gSeaNuEvGen", pINFO) << "Use UTM detector coordinate system";
        gGenPar.UseUTMSystem = true;
      }
      if(opt.compare("-tree")==0){              // tree format
        i++;
        string format=argv[i];
        if(format.compare("st")==0)gOptTreeFormat="SimpleTree";
        else if(format.compare("et")==0)gOptTreeFormat="EventTree";
        else{
          LOG("gSeaNuEvGen", pFATAL) << "Unknown output ROOT file format "<<format;
          gAbortingInErr = true;
          exit(1);
        }

        LOG("gSeaNuEvGen",pINFO)<< "Reading output ROOT file format "<<gOptTreeFormat;
      }
#ifdef _HETAU_ENABLED__
      if (opt.compare("-tau")==0) {
        i++;
        TauDecCh = argv[i];
        LOG("gSeaNuEvGen", pINFO) << "Simulating tau decay only with: " << TauDecCh;
      }
#endif
      if(opt.compare("-h")==0){         // help
        help = true;
      }
    }
  }

  if(help || argc==1){
    PrintSyntax();
    exit(0);
  }


  // load detector geometry file to set the can size
#ifdef _ANTARES_ENABLED__
  if(cantype.compare("DETX")==0){
    LoadGeometryDetx(gGenPar.detector);
  }
  else if(cantype.compare("DET")==0) {
    LoadGeometryDet(gGenPar.detector);
  }
#endif


  // set site depth to the values given with -d option and/or detector coordinates given with option -coord
  // but only if the can is set by hand, otherwise the -d option is ignored
  if(defcan){
    if(DefDepth)gGenPar.SiteDepth=depth;
    if(setcoord){
      gGenPar.SiteLatitude=latitude;
      gGenPar.SiteLongitude=longitude;
      assert(abs(declination)<=0.5*kPi);
      gGenPar.SiteDeclination=declination;
    }
    LOG("gSeaNuEvGen", pINFO)<< "Detector depth: "<< gGenPar.SiteDepth;
    LOG("gSeaNuEvGen", pINFO)<< "Detector coordinates Latitude and Longitude: "<< gGenPar.SiteLatitude <<" "<< gGenPar.SiteLongitude;
  }


  // set point source info... but after definition of site coordinates (in gal to eq coord transf. is required)
  if(gGenPar.GenMode.compare("POINT")==0){
    if(PointFile)LoadAstroInfo();
    else {
     LOG("gSeaNuEvGen", pINFO) << "Point-source declination: "<<gGenPar.Declination;
     gGenPar.Declination*=kPi/180.;
    }
  }


#ifdef _HETAU_ENABLED__
    if ( GenTau && gGenPar.EventGeneratorList.compare("CCHEDIS")==0 ){
      if ( TauDecCh!="muon" && TauDecCh!="shower" ) {
        LOG("gSeaNuEvGen",pFATAL) << "You are running NUTAU-CCHEDIS. Thus, you have to set -tau option to 'muon' or 'shower'";
        gAbortingInErr = true;
        exit(1);
      }
    }
    else if ( GenEl && gGenPar.EventGeneratorList.compare("GLRESAtomic-Tau")==0 ) {
      if ( TauDecCh!="muon" && TauDecCh!="shower" ) {
        LOG("gSeaNuEvGen",pFATAL) << "You are running NUEBAR-GLRESAtomic-Tau. Thus, you have to set -tau option to 'muon' or 'shower'";
        gAbortingInErr = true;
        exit(1);
      }
    }
    else {
      if ( TauDecCh=="muon" || TauDecCh=="shower" ) {
        LOG("gSeaNuEvGen",pFATAL) << "You are not running NUTAU-CCHEDIS. Thus, do not set -tau option";
        gAbortingInErr = true;
        exit(1);
      }
    }
#else
    if ( GenTau && gGenPar.EventGeneratorList.compare("CCHEDIS")==0 ){
      LOG("gSeaNuEvGen",pFATAL) << "You are running NUTAU-CCHEDIS. You need to link GSEAGEN to TAUOLA & TAUSIC in order to run this channel";
      gAbortingInErr = true;
      exit(1);
    }
    else if ( GenEl && gGenPar.EventGeneratorList.compare("GLRESAtomic-Tau")==0 ) {
      LOG("gSeaNuEvGen",pFATAL) << "You are running NUEBAR-GLRESAtomic-Tau. You need to link GSEAGEN to TAUOLA & TAUSIC in order to run this channel";
      gAbortingInErr = true;
      exit(1);
    }
#endif

  // check wich geometry should be used
  if(gGenPar.GenMode=="EVT" || gGenPar.GenMode=="BIN") {
    if(GenCan) {
      LOG("gSeaNuEvGen", pFATAL) << "Generate events in can (-C) is not compatible with -f FILE";
      gAbortingInErr = true;
      exit(1);
    }
    gGenPar.TrackEvts = 1;
  }
  else {
    if      ( GenCan )                                           gGenPar.TrackEvts = 0;
    else if ( gGenPar.EventGeneratorList.compare("NC")==0 )      gGenPar.TrackEvts = 0;
    else if ( gGenPar.EventGeneratorList.compare("NCHEDIS")==0 ) gGenPar.TrackEvts = 0;
    else if ( GenMu )                                            gGenPar.TrackEvts = 1;
    else if ( GenEl ) {
      gGenPar.TrackEvts = 0;
      if      ( gGenPar.EventGeneratorList.compare("GLRESAtomic-Mu")==0 ) gGenPar.TrackEvts = 1;
#ifdef _HETAU_ENABLED__
      else if ( gGenPar.EventGeneratorList.compare("GLRESAtomic-Tau")==0 ) {
        if      (TauDecCh=="muon")   gGenPar.TrackEvts = 2;
        else if (TauDecCh=="shower") gGenPar.TrackEvts = 3;
      }
#endif
    }
    else if ( GenTau ) {
      gGenPar.TrackEvts = 1;
#ifdef _HETAU_ENABLED__
      if      (TauDecCh=="muon")   gGenPar.TrackEvts = 2;
      else if (TauDecCh=="shower") gGenPar.TrackEvts = 3;
#endif
    }
  }

  if( (gGenPar.GenMode=="POINT" || gGenPar.GenMode=="DIFFUSE") && gGenPar.NTot == 0){
    LOG("gSeaNuEvGen", pFATAL) << "You must specify the number of events to generate";
    gAbortingInErr = true;
    exit(1);
  }
  if(gGenPar.EvMin<0 || gGenPar.EvMax<0 || gGenPar.EvMin>gGenPar.EvMax) {
    LOG("gSeaNuEvGen", pFATAL) << "Invalid energy range. Use `-e emin,emax', eg `-e 0.5,100.";
    gAbortingInErr = true;
    exit(1);
  }
  if(gGenPar.CtMin<-1 || gGenPar.CtMax<-1 || gGenPar.CtMin>gGenPar.CtMax) {
    LOG("gSeaNuEvGen", pFATAL) << "Invalid angular range. Use `-t ctmin,ctmax', eg `-t 0.,1.";
    gAbortingInErr = true;
    exit(1);
  }
  if(!CheckFluxOpt){
    LOG("gSeaNuEvGen", pFATAL) << "No flux info was specified! Use the -f option.";
    gAbortingInErr = true;
    exit(1);
  }




  //
  // print-out summary
  //


  ostringstream rotation;
  rotation << "\t| " <<  gGenPar.Rotation.XX() << "  " << gGenPar.Rotation.XY() << "  " << gGenPar.Rotation.XZ() << " |\n";
  rotation << "\t| " <<  gGenPar.Rotation.YX() << "  " << gGenPar.Rotation.YY() << "  " << gGenPar.Rotation.YZ() << " |\n";
  rotation << "\t| " <<  gGenPar.Rotation.ZX() << "  " << gGenPar.Rotation.ZY() << "  " << gGenPar.Rotation.ZZ() << " |\n";

  if(gGenPar.GenMode=="BIN" || gGenPar.GenMode=="EVT"){

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

    if(gGenPar.NBin>1){
      LOG("gSeaNuEvGen", pNOTICE) << "Number of bins must be 1 when reading neutrinos from file. User input ignored...";
      gGenPar.NBin=1;
    }


    ostringstream fluxinfo;
    fluxinfo << "Neutrinos and muons from file: "<<gGenPar.FluxFiles[0].FileNames[0]<<"\n";

    LOG("gSeaNuEvGen", pNOTICE)
       << "\n"
       << "\n ****** MC Job (" << gGenPar.RunNu << ") Settings ****** "
       << "\n @@ Flux"
       << "\n\t" << fluxinfo.str()

       << "\n\t Using number of energy bins = " << gGenPar.NBin

       << "\n @@ Coordinate transformation (Rotation THZ -> User-defined coordinate system)"
       << "\n" << rotation.str()
       << "\n\n";
  }
  else {

    PDGLibrary * pdglib = PDGLibrary::Instance();
    ostringstream fluxinfo;
    fluxinfo << "Using real fluxes: \n";

    for(unsigned iFlux=0; iFlux<gGenPar.FluxFiles.size(); iFlux++){
      fluxinfo<<pdglib->Find(gGenPar.FluxFiles[iFlux].NuType)->GetName()<<" "<<gGenPar.FluxFiles[iFlux].FluxSimul<<" ";
      for(unsigned iFile=0; iFile<gGenPar.FluxFiles[iFlux].FileNames.size(); iFile++)
        fluxinfo<<gGenPar.FluxFiles[iFlux].FileNames[iFile]<<" ";
      fluxinfo<<"\n";
    }

    ostringstream expinfo;
    if(gGenPar.NTot > 0)            { expinfo << gGenPar.NTot            << " events";   }

    LOG("gSeaNuEvGen", pNOTICE)
       << "\n"
       << "\n ****** MC Job (" << gGenPar.RunNu << ") Settings ****** "
       << "\n @@ Flux"
       << "\n\t" << fluxinfo.str()
       << "\n @@ Exposure"
       << "\n\t" << expinfo.str()
       << "\n @@ Cuts"
       << "\n\t Using energy range = (" << gGenPar.EvMin << " GeV, " << gGenPar.EvMax << " GeV)"
       << "\n\t Using cosTheta range = (" << gGenPar.CtMin << ", " << gGenPar.CtMax << ")"
       << "\n\t Using generation spectral index = " << gGenPar.Alpha
       << "\n\t Using number of energy bins = " << gGenPar.NBin
       << "\n @@ Coordinate transformation (Rotation THZ -> User-defined coordinate system)"
       << "\n" << rotation.str()
       << "\n\n";
  }


}