//____________________________________________________________________________
/*!

\class   genie::BinParam, genie::GenParam

\brief   Data structure for gSeaGen

\author  Carla Distefano <distefano_c@lns.infn.it>
         LNS-INFN, Catania
         Included modifications by Alfonso Garcia <alfonsog@nikhef.nl>
         NIKHEF, Amsterdam (January 2019)

\created May 13, 2012

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

*/
//____________________________________________________________________________



#include "SeaEvent/GBinParam.h"

using namespace std;

GenParam::GenParam(){

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

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


  this->SetDefParams();
  this->LoadDefParams();
  this->LoadMediaComp();
};
//____________________________________________________________________________

GenParam::~GenParam(){};
//____________________________________________________________________________

void GenParam::SetDefParams(void){

// Set default parameters, this function must be called before the GetCommandLineArgs or at beginning of it

  //running time
  time(&RunTime);
  
  RunNu = 100000000;

  NTot = 0;
  Rotation.SetToIdentity();

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


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

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

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

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

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


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

  RanSeed=12345;

  PropMode=false;
  Drawing="surface";

  GenMode="DIFFUSE";
  TGen = 1.;             	// 1 second

  EvFilePrefix = "gseagen";     // event file prefix (override with -o)
  NativeOutput = false;
  WriteEvt = false;
  WriteKM3NeT = false;
  WriteNative = false;

#if defined _KM3NET_ENABLED__	  
  WriteFormat.push_back("km3");
#else
  WriteFormat.push_back("et");
#endif




  PropCode  = "PropaMuon";
  PropCodeVer="";
  PropagationStep=5.;
  DistAfterFirstStep=100.;
  MuonRangeTolerance=0.;
  PropagationTolerance=5.;  // the default tolerance is 5m

  GenCode ="GENIE";

  RTOpt = "can";

  TrackEvts = 0;

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

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

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


  EventGeneratorList="Default";

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

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

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

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

  return;

};

//____________________________________________________________________________
void GenParam::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...";
      cout<< filename <<" does not exist --> check env. parameter DEFPARAM; exiting..."<<endl;
      exit(1);
    }
    else
////      LOG("gSeaNuEvGen", pINFO) <<"Reading new default parameters from file "<< filename;
      cout<< "Reading new default parameters from file "<< filename<<endl;

    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... ";
          cout<<" Invalid root element in filename " << filename <<"; exiting... "<<endl;
	  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"){
          RunNu=atoi((const char *)value);
        }
        else if(selem=="RanSeed"){
          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);
            Can1 = atof(cansize[0].c_str());
            Can2 = atof(cansize[1].c_str());
            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);
            OriginX = atof(origsize[0].c_str());
            OriginY = atof(origsize[1].c_str());
            OriginZ = atof(origsize[2].c_str());
          }
        }
        else if(selem=="SiteLatitude"){
          SiteLatitude=atof((const char *)value);
        }
        else if(selem=="SiteLongitude"){
          SiteLongitude=atof((const char *)value);
        }
        else if(selem=="SiteDepth"){
          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);
            EvMin = atof(rangesize[0].c_str());
            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);
            CtMin = atof(rangesize[0].c_str());
            CtMax = atof(rangesize[1].c_str());
          }
        }
        else if(selem=="NBin"){
          NBin=atoi((const char *)value);
        }
        else if(selem=="Alpha"){
        Alpha=atof((const char *)value);
        }
        else if(selem=="RTOpt"){
          RTOpt=(const char *)value;
          RTOpt.erase(remove(RTOpt.begin(),RTOpt.end(),' '),RTOpt.end());
          RTOpt.erase(remove(RTOpt.begin(),RTOpt.end(),'\t'),RTOpt.end());
        }
        else if(selem=="NAbsLength"){
          NAbsLength=atof((const char *)value);
        }
        else if(selem=="AbsLength"){
          AbsLength=atof((const char *)value);
        }
        else if(selem=="HBedRock"){
          HBedRock=atof((const char *)value);
        }
        else if(selem=="PropCode"){
          PropCode=(const char *)value;
          PropCode.erase(remove(PropCode.begin(),PropCode.end(),' '),PropCode.end());
          PropCode.erase(remove(PropCode.begin(),PropCode.end(),'\t'),PropCode.end());
          PropCode = genie::utils::str::ToLower(PropCode);
          if(PropCode.compare("propamuon")==0)PropCode = "PropaMuon";
          if(PropCode.compare("music")==0){
#ifdef _MUSIC_ENABLED__
            PropCode = "MUSIC";
#else
////            LOG("gSeaNuEvGen", pFATAL) <<"Propagation code: "<<PropCode<<", is not enabled, change file "<<filename <<"; exiting...";
////            gAbortingInErr = true;
            cout<<"Propagation code: "<<PropCode<<", is not enabled, change file "<<filename <<"; exiting..."<<endl;
            exit(1);
#endif
          }
          if(PropCode.compare("proposal")==0){
#ifdef _PROPOSAL_ENABLED__
            PropCode = "PROPOSAL";
#else
////            LOG("gSeaNuEvGen", pFATAL) <<"Propagation code: "<<PropCode<<", is not enabled, change file "<<filename <<"; exiting...";
////            gAbortingInErr = true;
            cout<<"Propagation code: "<<PropCode<<", is not enabled, change file "<<filename <<"; exiting..."<<endl;
            exit(1);
#endif
          }
          if(PropCode.compare("jpp")==0){
#ifdef _JPP_ENABLED__
            PropCode = "JPP";
#else
////            LOG("gSeaNuEvGen", pFATAL) <<"Propagation code: "<<PropCode<<", is not enabled, change file "<<filename <<"; exiting...";
////            gAbortingInErr = true;
            cout<<"Propagation code: "<<PropCode<<", is not enabled, change file "<<filename <<"; exiting..."<<endl;
            exit(1);
#endif
          }
        }
        else if(selem=="WriteParticlesMode"){
          WriteParticlesMode=atoi((const char *)value);
        }
        else if(selem=="WeightOpt"){
        WeightOpt=atoi((const char *)value);
        }
        else if(selem=="EvFilePrefix"){
          EvFilePrefix=(const char *)value;
          EvFilePrefix.erase(remove(EvFilePrefix.begin(),EvFilePrefix.end(),' '),EvFilePrefix.end());
          EvFilePrefix.erase(remove(EvFilePrefix.begin(),EvFilePrefix.end(),'\t'),EvFilePrefix.end());
        }
        else if(selem=="InpXSecFile"){
          InpXSecFile=(const char *)value;
          InpXSecFile.erase(remove(InpXSecFile.begin(),InpXSecFile.end(),' '),InpXSecFile.end());
          InpXSecFile.erase(remove(InpXSecFile.begin(),InpXSecFile.end(),'\t'),InpXSecFile.end());
        }
        else if(selem=="EventGeneratorList"){
          EventGeneratorList=(const char *)value;
          EventGeneratorList.erase(remove(EventGeneratorList.begin(),EventGeneratorList.end(),' '),EventGeneratorList.end());
          EventGeneratorList.erase(remove(EventGeneratorList.begin(),EventGeneratorList.end(),'\t'),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)UseUTMSystem = true;
          else if(CheckIfUTM.compare("false")==0)UseUTMSystem = false;
          else {
////            LOG("gSeaNuEvGen", pFATAL) <<"Error reading new default parameters: Param UseUTMSystem must be true or false; exiting...";
////            gAbortingInErr = true;
            cout<<"Error reading new default parameters: Param UseUTMSystem must be true or false; exiting..."<<endl;
            exit(1);
          }
        }
      }

      xmlFree(name);
      xmlFree(value);
    }

    }

  return;

};

//____________________________________________________________________________
void GenParam::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...";
      cout << filename <<" does not exist --> check env. parameter MEDIACOMP; exiting..."<<endl;
      exit(1);
    }
    else
//      LOG("gSeaNuEvGen", pFATAL) << filename <<" does not exist --> check env. parameter MEDIACOMP; exiting...";
      cout << filename <<" does not exist --> check env. parameter MEDIACOMP; exiting..."<<endl;

    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... ";
        cout << " invalid root element in filename " << filename <<"; exiting... "<<endl; 
        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;
          SeaWaterComp.clear();
//          LOG("gSeaNuEvGen", pINFO) <<"Defining new SeaWater composition";
          cout <<"Defining new SeaWater composition"<<endl;
          if(xdensity!=NULL){
            sdensity = utils::str::TrimSpaces((const char*)xdensity);
//            LOG("gSeaNuEvGen",pINFO) <<"Defining new SeaWater density";
            cout <<"Defining new SeaWater density"<<endl;
            RhoSW=(double)atof(sdensity.c_str());
          }
        }
        if(smedia=="Rock" && CheckRock){
          CheckRock=false;
          RockComp.clear();
//          LOG("gSeaNuEvGen", pINFO) <<"Defining new Rock composition";
          cout <<"Defining new Rock composition"<<endl;
          if(xdensity!=NULL){
            sdensity = utils::str::TrimSpaces((const char*)xdensity);
//            LOG("gSeaNuEvGen",pINFO) <<"Defining new Rock density";
            cout <<"Defining new Rock density"<<endl;
            RhoSR=(double)atof(sdensity.c_str());
          }
        }
        if(smedia=="Mantle" && CheckMantle){
          CheckMantle=false;
          MantleComp.clear();
//          LOG("gSeaNuEvGen", pINFO) <<"Defining new Mantle composition";
          cout <<"Defining new Mantle composition"<<endl;
//          if(xdensity!=NULL) LOG("gSeaNuEvGen",pINFO) <<"Defining new Mantle density, ignoring...";
          if(xdensity!=NULL) cout <<"Defining new Mantle density, ignoring..."<<endl;        
	}
        if(smedia=="Core" && CheckCore){
          CheckCore=false;
          CoreComp.clear();
//          LOG("gSeaNuEvGen", pINFO) <<"Defining new Core composition";
          cout <<"Defining new Core composition"<<endl;
//          if(xdensity!=NULL) LOG("gSeaNuEvGen",pINFO) <<"Defining new Core density, ignoring...";
          if(xdensity!=NULL) cout <<"Defining new Core density, ignoring..."<<endl;
        }

        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... ";
          cout << smedia <<": element "<< selem << " not present in the list..., exiting... "<<endl;
          exit(1);
        }
        if(smedia=="SeaWater")  SeaWaterComp.insert(map<int, double>::value_type(pdg,atof((const char *)value)));
        if(smedia=="Rock")      RockComp.insert(map<int, double>::value_type(pdg,atof((const char *)value)));
        if(smedia=="Mantle")    MantleComp.insert(map<int, double>::value_type(pdg,atof((const char *)value)));
        if(smedia=="Core")      CoreComp.insert(map<int, double>::value_type(pdg,atof((const char *)value)));
      }
      xmlFree(name);
      xmlFree(value);
    }
  }

  return;
}