//________________________________________________________________________________________
/*!

\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 "definitions.h"
// #include <gperftools/profiler.h>

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

#include <boost/algorithm/string.hpp>
#include <boost/filesystem.hpp>

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


// GENIE

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


// gSeaGen

#include "SeaNuDrivers/GGenerateEventGenie.h"
#ifdef _KM3NET_ENABLED__
#include "SeaNuDrivers/GGenerateEventExternal.h"
#endif

#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"



#ifdef _PROPOSAL_ENABLED__
#include "PROPOSAL/PROPOSAL.h"
#endif

using namespace std;

using namespace genie;
using namespace genie::geometry;

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

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

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

void                    LoadAstroInfo      (void);

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

void			InitializeGenie	   (int argc, char ** argv);
void                    PrintGenieConfig   (void);



// User-specified options:

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

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

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

  // Calculate the UTM convergence angle
  gGenPar.UTMConvergeAngle = GAstro::CalcUTMConvergeAngle(gGenPar.SiteLatitude, gGenPar.SiteLongitude);

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

  // Define the generator
  GGenerateEventGenie * gGenEvtGenie=0;
#ifdef _KM3NET_ENABLED__
  GGenerateEventExternal * gGenEvtExt=0;
#endif


  GGenerateEvent * gGenEvt=0;

  if(gGenPar.GenCode.compare("GENIE")==0 || gGenPar.GenCode.compare("CORSIKA")==0){
    InitializeGenie(argc,argv);
    gGenEvtGenie = new GGenerateEventGenie(&gGenPar);
    gGenEvt = dynamic_cast <GGenerateEvent*> (gGenEvtGenie);
    PrintGenieConfig();
  }
#ifdef _KM3NET_ENABLED__
  else if(gGenPar.GenCode.compare("EXTERNAL")==0){
    gGenEvtExt = new GGenerateEventExternal(&gGenPar);
    gGenEvt = dynamic_cast <GGenerateEvent*> (gGenEvtExt);
  }
#endif

#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);
  if (gGenPar.WriteEvt) {
    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 = 0;
  if (gGenPar.WriteKM3NeT) {
    km3net_writer = new KM3NeTWrite( &gGenPar );
    km3net_writer->WriteHeader();
    SLOG("gSeaNuEvGen", pINFO) << "km3net writer initialized";
  }
#endif
  SeaNtpWriter * ntprt = 0;
  if (gGenPar.WriteNative) 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();
    // ProfilerStart("test_gseagen");
    // generate next event at the can
    while (1){
      
      gGenEvt->GenerateEvent();
      if(gGenEvt->GetNFluxNeutrinos()>VecBinPar[iBin].NEvBinScal || gGenEvt->EoF) break;

      if (gGenPar.WriteNative)
        ntprt->WriteEvent(gGenEvt->fSeaEvent);
      
      if (gGenPar.NativeOutput)
         gGenEvt->WriteNative(gGenEvt->fSeaEvent);
      
#ifdef _ANTARES_ENABLED__
      if (gGenPar.WriteEvt) {
        evtwrt->WriteEvent(gGenEvt->fSeaEvent);
      }
#endif
#ifdef _KM3NET_ENABLED__
      if (gGenPar.WriteKM3NeT) {
        km3net_writer->WriteEvent(gGenEvt->fSeaEvent);
      }
#endif
      VecBinPar[iBin].NEvBinGen++;
      if ( (gGenPar.SilentCORSIKA) && !(gGenEvt->fNEvt%100000) ) LOG("gSeaNuEvGen", pNOTICE) << "Generated events "<<gGenEvt->fNEvt;
      else if( !(gGenEvt->fNEvt%500) ) LOG("gSeaNuEvGen", pNOTICE) << "Generated events "<<gGenEvt->fNEvt;
      gGenEvt->CleanEvent();  // clean-up
    }
    // ProfilerStop();
    if( !(gGenPar.GenMode=="BIN" || gGenPar.GenMode=="EVT") ) {
      PrintStat();
      LOG("gSeaNuEvGen", pNOTICE) << "Done bin: " << iBin+1;
    }
  }
  // clean-up

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

#ifdef _KM3NET_ENABLED__
  if (gGenPar.WriteKM3NeT) {
    if (km3net_writer) delete km3net_writer;
  }
#endif
  if (gGenPar.WriteNative) {
    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 InitializeGenie(int argc, char ** argv){

  RunOpt::Instance()->EnableBareXSecPreCalc(true);
  RunOpt::Instance()->ReadFromCommandLine(argc,argv);


#ifdef _GENIEVERSIONGTEQ3__

//  RunOpt::Instance()->SetEventGeneratorList(gGenPar.EventGeneratorList);
//  RunOpt::Instance()->SetTuneName(gGenPar.GenTune);

  if ( ! RunOpt::Instance()->Tune() ) {
    LOG("gevgen", pFATAL) << " No TuneId in RunOption";
    exit(-1);
  }
  RunOpt::Instance()->BuildTune();

#endif

  // Initialization of random number generators, cross-section table,
  // messenger thresholds, cache file
  utils::app_init::MesgThresholds(RunOpt::Instance()->MesgThresholdFiles());
  utils::app_init::CacheFile(RunOpt::Instance()->CacheFile());
  utils::app_init::RandGen((long int)gGenPar.RanSeed);

  // in case of tune GVLE18_01a_00_000 we calculate xsec run time if file does not exist
  if(gGenPar.GenTune.compare("GVLE18_01a_00_000")==0){
    if(access(gGenPar.InpXSecFile.c_str(), F_OK) != -1) utils::app_init::XSecTable(gGenPar.InpXSecFile, false);
  }
  else
    utils::app_init::XSecTable(gGenPar.InpXSecFile, false);

  // Set GHEP print level
  GHepRecord::SetPrintLevel(RunOpt::Instance()->EventRecordPrintLevel());
}
//________________________________________________________________________________________
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 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();
  tRandom->SetSeed((long int)gGenPar.RanSeed);

  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.ZminVol = SeaGeom->GetZminVol();
    param.ZmaxVol = SeaGeom->GetZmaxVol();
    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 {
        double temp = (VecBinPar[iBin].RT/VecBinPar[VecBinPar.size()-1].RT);
        VecBinPar[iBin].NEvBinScal = tRandom->Poisson(VecBinPar[iBin].NEvBin*temp*temp);
      }
    }
    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.ZminVol   = VecBinPar[VecBinPar.size()-1].ZminVol;
  gGenPar.ZmaxVol   = VecBinPar[VecBinPar.size()-1].ZmaxVol;
  gGenPar.SeaCenter = VecBinPar[VecBinPar.size()-1].SeaCenter;
  gGenPar.RockCenter = VecBinPar[VecBinPar.size()-1].RockCenter;

  if( !(gGenPar.GenMode=="BIN" || gGenPar.GenMode=="EVT") ) 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(" ");

  std::printf(GSEAGEN_INFO.c_str(), ver.c_str());

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



  cout << ARGS_INFO;

}
#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_min-ZOM-OffSet);
    gGenPar.Can2 = min(gGenPar.SiteDepth, z_max+ZOM+OffSet);
    gGenPar.Can3 = radius+ROM+OffSet;


    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_min-OffSet);
    gGenPar.Can2 = min(gGenPar.SiteDepth, z_max+OffSet);
    gGenPar.Can3 = radius+OffSet;


    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.*0.5*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 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]);
      }
/*
      if(opt.compare("--write-evt")==0){
        i++;
        LOG("gSeaNuEvGen",pDEBUG) << "Write .evt output file";
        gGenPar.WriteEvt = (bool) atoi(argv[i]);
      }
      if(opt.compare("--write-km3net")==0){
        i++;
        LOG("gSeaNuEvGen",pDEBUG) << "Write km3net .root output file";
        gGenPar.WriteKM3NeT = (bool) atoi(argv[i]);
      }
      if(opt.compare("--write-native")==0){
        i++;
        LOG("gSeaNuEvGen",pDEBUG) << "Write native .root output file";
        gGenPar.WriteNative = (bool) 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 iRot=0; iRot<3; iRot++) gGenPar.CustomRotation[iRot] = atof(coords[iRot].c_str());
        if(rot_unit.compare("DEG")==0) for (int iRot=0; iRot<3; iRot++) gGenPar.CustomRotation[iRot] *= 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("--corsika-max-nrecords")==0){
        i++;
        gGenPar.MaxNRec = atof(argv[i]);
        if (gGenPar.MaxNRec==0) LOG("gSeaNuEvGen",pDEBUG) << "All records will be read from CORSIKA files";
        else LOG("gSeaNuEvGen",pDEBUG) << "At most " << gGenPar.MaxNRec << " records will be read from CORSIKA files";
      }
      if(opt.compare("--corsika-less-verbose")==0){
        i++;
        gGenPar.SilentCORSIKA = atof(argv[i]);
        assert( (gGenPar.SilentCORSIKA==0) || (gGenPar.SilentCORSIKA==1) );
        if (gGenPar.SilentCORSIKA==1) LOG("gSeaNuEvGen",pDEBUG) << "CORSIKA LOG messages will be less verbose";
      }
      if(opt.compare("--corsika-only-convert")==0){
        i++;
        gGenPar.NoPropagationCORSIKA = atof(argv[i]);
        assert( (gGenPar.NoPropagationCORSIKA==0) || (gGenPar.NoPropagationCORSIKA==1) );
        if (gGenPar.NoPropagationCORSIKA==1) LOG("gSeaNuEvGen",pINFO) << "CORSIKA input will NOT be propagated, only converted";
      }
      if(opt.compare("--muon-range-tolerance")==0){
        i++;
        gGenPar.MuonRangeTolerance = atof(argv[i]);
        assert( gGenPar.MuonRangeTolerance>=0 );
        if (gGenPar.MuonRangeTolerance>0) LOG("gSeaNuEvGen",pDEBUG) << "The muon range tolerance is set to "<< gGenPar.MuonRangeTolerance <<" meters";
      }
      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];
  
        boost::to_upper(gGenPar.GenCode);

        if(gGenPar.GenCode.compare("GENIE")==0)LOG("gSeaNuEvGen", pNOTICE) << "Neutrino interactions simulated with GENIE";
#ifdef _KM3NET_ENABLED__
        else if(gGenPar.GenCode.compare("EXTERNAL")==0)LOG("gSeaNuEvGen", pNOTICE) << "Neutrino interactions simulated with an external generator";
#endif
        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);
	boost::to_lower(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){
	  gGenPar.PropCodeVer=PROPOSAL_VERSION;
	  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;
      }
      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";
          }
        }
      }
#ifdef _ANTARES_ENABLED__
      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";
            gGenPar.GenCode ="CORSIKA";
          }

          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("-format")==0){              // tree format
        i++;
        gGenPar.WriteFormat.clear();
        gGenPar.WriteFormat = utils::str::Split(argv[i], ",");
      }
      
#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("--propagation-step")==0){
        i++;
        gGenPar.PropagationStep = atof(argv[i]);
        LOG("gSeaNuEvGen", pINFO) << "Reading the propagation step size (in meters): " << gGenPar.PropagationStep;
        assert(gGenPar.PropagationStep>0);
      }
      if(opt.compare("--approach-dist")==0){
        i++;
        gGenPar.DistAfterFirstStep = atof(argv[i]);
        LOG("gSeaNuEvGen", pINFO) << "Reading the distance from the can to which the first propagation step will be done (in meters): " << gGenPar.PropagationStep;
        assert(gGenPar.DistAfterFirstStep>0);
      }
      if(opt.compare("-h")==0){         // help
        help = true;
      }
    }
  }

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


// activate the writing of output files
  for(unsigned iFormat=0; iFormat<gGenPar.WriteFormat.size(); iFormat++){

    if(gGenPar.WriteFormat[iFormat].compare("et")==0){
      gGenPar.WriteNative = true;
      gOptTreeFormat="EventTree";
    }
#if defined _ANTARES_ENABLED__
    if(gGenPar.WriteFormat[iFormat].compare("evt")==0){
      gGenPar.WriteEvt = true;
    }
#endif
#if defined _KM3NET_ENABLED__
    if(gGenPar.WriteFormat[iFormat].compare("km3")==0){
      gGenPar.WriteKM3NeT = true;
    }
#endif
    if(gGenPar.WriteFormat[iFormat].compare("gen")==0){
      gGenPar.NativeOutput = true;
    }
  }
  
  // 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 || gGenPar.EventGeneratorList.compare("GLRES-Tau")==0)) {
      if ( TauDecCh!="muon" && TauDecCh!="shower" ) {
        LOG("gSeaNuEvGen",pFATAL) << "You are running "<<gGenPar.EventGeneratorList<<". 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 || gGenPar.EventGeneratorList.compare("GLRES-Tau")==0)) {
      LOG("gSeaNuEvGen",pFATAL) << "You are running "<<gGenPar.EventGeneratorList<<". 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.EventGeneratorList.compare("GLRES-Mu")==0 ) gGenPar.TrackEvts = 1;
#ifdef _HETAU_ENABLED__
      else if ( gGenPar.EventGeneratorList.compare("GLRESAtomic-Tau")==0 || gGenPar.EventGeneratorList.compare("GLRES-Tau")==0) {
        if      (TauDecCh=="muon")   gGenPar.TrackEvts = 2;
        else if (TauDecCh=="shower") gGenPar.TrackEvts = 3;
      }
#endif
      else if ( gGenPar.EventGeneratorList.compare("GLRESAtomic")==0 || gGenPar.EventGeneratorList.compare("GLRES")==0) gGenPar.TrackEvts = 1;
    }
    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);
  }


  if(gGenPar.NativeOutput && !(gGenPar.GenMode=="POINT" || gGenPar.GenMode=="DIFFUSE")){
    LOG("gSeaNuEvGen", pFATAL) << "GENIE native output cannot be written if not neutrino generation! Remove -w option.";
    gAbortingInErr = true;
    exit(1);
  }


  // check if the can is ok: can1>=0; can2<=Depth; can2>=can1

  if(gGenPar.Can1<0){
    LOG("gSeaNuEvGen", pFATAL) << "CAN Zmin value must be greater than zero. See -c option or file input with DEFPARAM environmental parameter.";
    gAbortingInErr = true;
    exit(1);
  }

  if(gGenPar.Can2>gGenPar.SiteDepth){
    LOG("gSeaNuEvGen", pFATAL) << "CAN Zmax value must be lower than the site depth. See -c option or file input with DEFPARAM environmental parameter.";
    gAbortingInErr = true;
    exit(1);
  }

  if(gGenPar.Can2<gGenPar.Can1){
    LOG("gSeaNuEvGen", pFATAL) << "CAN Zmax cannot be lower than Zmin. See -c option or file input with DEFPARAM environmental parameter.";
    gAbortingInErr = true;
    exit(1);
  }



  // check if hbedrock is ok: can1>=0; hbedrock>=0 and hbedrock=0 if can1>0

  if(gGenPar.HBedRock<0){
    LOG("gSeaNuEvGen", pFATAL) << "HBedRock cannot be lower than zero. See -bedrock option or file input with DEFPARAM environmental parameter.";
    gAbortingInErr = true;
    exit(1);
  }

  if(gGenPar.Can1>0 && gGenPar.HBedRock>0){
    LOG("gSeaNuEvGen", pFATAL) << "HBedRock must be zero if CAN bottom in not on the seabed (Zmin>0). See -bedrock option or file input with DEFPARAM environmental parameter.";
    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<<" ";
      fluxinfo<<gGenPar.FluxFiles[iFlux].NuType<<" "<<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";
  }


}


///////////////////////////////////////////////////////////////7
void PrintGenieConfig(){


  string fullname;
  string file;
  string command;
  
  file = "UnstableParticleDecayer.xml";
  fullname = utils::xml::GetXMLFilePath(file);
  
  LOG("gSeaNuEvGen", pINFO)<< "Using GENIE config file: "<< fullname;
        
#ifdef _KM3NET_ENABLED__
  
  command = "cat " + fullname;
  system(command.c_str());

#endif
  
  file = "CommonDecay.xml";
  fullname = utils::xml::GetXMLFilePath(file);
  
  LOG("gSeaNuEvGen", pINFO)<< "Using GENIE config file: "<< fullname;
  
#ifdef _KM3NET_ENABLED__
  
  command = "cat " + fullname;
  system(command.c_str());

#endif
  
  
    
  return;
}