//____________________________________________________________________________
/*!

\class   genie::EvtWrite

\brief   A class to write gSeaGen output in a EVT file

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

\created May 22, 2012

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

*/
//____________________________________________________________________________


#ifdef _ANTARES_ENABLED__

#include "OutputWriters/EvtWrite.h"

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


//___________________________________________________________________________
EvtWrite::EvtWrite(GenParam * GenPar, bool OptPDG)
{

  fGenPar = GenPar;

  fOptPDG = OptPDG;

  fNewEvt = true;

  ostringstream fnstr;
  fnstr << fGenPar->EvFilePrefix  << "."
        << fGenPar->RunNu
        << ".evt";

  string OutFileName=fnstr.str();
  fOutFile.open(OutFileName.c_str());

  fPdg = new TDatabasePDG();

  fEvt = new event();

}

//___________________________________________________________________________
EvtWrite::EvtWrite(GenParam * GenPar, bool OptPDG, event *evt)
{
  fGenPar = GenPar;

  fOptPDG = OptPDG;

  fNewEvt = true;

  ostringstream fnstr;
  fnstr << fGenPar->EvFilePrefix  << "."
        << fGenPar->RunNu
        << ".evt";

  string OutFileName=fnstr.str();
  fOutFile.open(OutFileName.c_str());

  fPdg = new TDatabasePDG();

  fEvt = evt;
  fNewEvt = false;

}

//___________________________________________________________________________
EvtWrite::~EvtWrite(){
  fOutFile.close();
  if(fPdg) delete fPdg;
  if(fEvt) delete fEvt;
}

//___________________________________________________________________________
void EvtWrite::WriteHeader(void)
{

  tm *date;
  date = localtime (&fGenPar->RunTime);

  char ch[1024];

  if(fNewEvt)fEvt->defrun(fGenPar->RunNu);
      

  if ( (fGenPar->GenMode.compare("BIN")==0) || (fGenPar->GenMode.compare("EVT")==0) ){
    sprintf(ch,"CORSIKA %5.3f %2.2d  %2.2d%2.2d ", fGenPar->CORSIKAVer,fGenPar->CORSIKADate,0,0);
    fEvt->taga("simul",ch);
  }

  sprintf(ch, "gSeaGen %s %2.2d%2.2d%2.2d  %2.2d%2.2d",fGenPar->gSeaGenVer.c_str(),
      date->tm_year-100,date->tm_mon+1,date->tm_mday,date->tm_hour,date->tm_min);

  fEvt->taga("simul",ch);

  if ( fGenPar->GenMode.compare("POINT")==0 || fGenPar->GenMode.compare("DIFFUSE")==0 ) {

    sprintf(ch, "GENIE %s",fGenPar->GenieVer.c_str());
    fEvt->taga("physics",ch);

#ifdef _GENIEVERSIONGTEQ3__
    string ReWgtVer;
    string VerFile = gSystem->Getenv("GENIE_REWEIGHT");
    VerFile.append("/VERSION");
    ifstream file_stream(VerFile.c_str(), ios::in);
    file_stream>>ReWgtVer;
    file_stream.close();

//    sprintf(ch, "GENIE_REWEIGHT %s %2.2d%2.2d%2.2d  %2.2d%2.2d%2.2d",ReWgtVer.c_str(),
//      date->tm_year-100,date->tm_mon+1,date->tm_mday,date->tm_hour,date->tm_min,date->tm_sec);
//    fEvt->taga("simul",ch);
    
    sprintf(ch, "GENIE_REWEIGHT %s",ReWgtVer.c_str());
    fEvt->taga("physics",ch);
    
    sprintf(ch, "GENIE_TUNE %s %s", fGenPar->GenTune.c_str(), fGenPar->EventGeneratorList.c_str());
    fEvt->taga("physics",ch);
    
        
#endif

    sprintf(ch,"%10.2E %10.2E %10.2E %10.2E",fGenPar->EvMin,fGenPar->EvMax,fGenPar->CtMin,fGenPar->CtMax);
    fEvt->taga("cut_nu",ch);

    if(fGenPar->InpXSecFile.size()>0){
      sprintf(ch, "%s", fGenPar->InpXSecFile.c_str());
      fEvt->taga("XSecFile",ch);
    }

  }
  
  else if ( (fGenPar->GenMode.compare("BIN")==0) || (fGenPar->GenMode.compare("EVT")==0) ) {
    sprintf(ch, "CORSIKA %s %s",fGenPar->CORSIKAHadronicModelLE.c_str(),fGenPar->CORSIKAHadronicModelHE.c_str());
    fEvt->taga("physics",ch);  
  }  
  
  


  if(gSystem->Getenv("GXMLPATH")){
    string gxmlpath=gSystem->Getenv("GXMLPATH");
    sprintf(ch, "%s",gxmlpath.c_str());
    fEvt->taga("GXMLPATH",ch);
  }

  fEvt->taga("drawing",fGenPar->Drawing);

  sprintf(ch, "%+4.2f", -fGenPar->Alpha);
  if(fNewEvt)fEvt->taga("spectrum",ch);

  sprintf(ch, "gSeaGen %d", fGenPar->RanSeed);
  fEvt->taga("seed",ch);


  if(fGenPar->detector.size()>0){
    sprintf(ch, "gSeaGen %s", fGenPar->detector.c_str());
    fEvt->taga("detector",ch);
  }


  if(fGenPar->TrackEvts>0){
    sprintf(ch, "gSeaGen %s %s", fGenPar->PropCode.c_str(),fGenPar->PropCodeVer.c_str());
    fEvt->taga("physics",ch);
  }

//  sprintf(ch, "%f", fGenPar->TGen);
//  fEvt->taga("tgen",ch);

  if ( fGenPar->GenMode.compare("BIN")==0 ){
    sprintf(ch,"%10d", fGenPar->Primary );
    fEvt->taga("primary",ch);
  }
  else{
    ostringstream nulist;
    for(unsigned iNuType=0; iNuType<fGenPar->NuList.size(); iNuType++)nulist << fGenPar->NuList[iNuType]  << " ";
    string snulist=nulist.str();
    fEvt->taga("primary",snulist);      
  }
  
  if ( (fGenPar->GenMode.compare("BIN")==0) || (fGenPar->GenMode.compare("EVT")==0) ) sprintf(ch, "FILE"); 
  else sprintf(ch, "%s", fGenPar->GenMode.c_str());
  fEvt->taga("source_mode",ch);
  
  if ( fGenPar->GenMode == "POINT") {
    if ( fGenPar->SourceName != "" ){
      ostringstream name;
      name<<"Dec:";
      name.precision(3);
      name<<fixed;    
      name << fGenPar->Declination*180./kPi;
      cout << name.str();
    }
    sprintf(ch,"%f %f %f %s",fGenPar->RightAscension,fGenPar->Declination,fGenPar->SourceRadius*180./kPi,fGenPar->SourceName.c_str());
    fEvt->taga("source_info",ch);
  }    
  

  if(fGenPar->GenMJD){
    sprintf(ch, "%d %d", fGenPar->TimeStampStart,fGenPar->TimeStampStop);
    fEvt->taga("time_interval",ch);
  }




  sprintf(ch, "%10.3f %10.3f %10.3f %10.3f %10.3f",
             fGenPar->OriginX,fGenPar->OriginY,fGenPar->Can1-fGenPar->OriginZ,fGenPar->Can2-fGenPar->OriginZ,fGenPar->Can3);
  fEvt->taga("fixedcan",ch);



  sprintf(ch, "%10.6f %10.6f",
             fGenPar->SiteLatitude,fGenPar->SiteLongitude);
  fEvt->taga("coord_site",ch);


  sprintf(ch, "%10.3f", fGenPar->SiteDepth);
  fEvt->taga("depth",ch);

  if(fGenPar->HBedRock>0){
    sprintf(ch,"%10.3f", fGenPar->HBedRock );
    fEvt->taga("bedrock",ch);
  }   


  double IntVol=constants::kPi*pow(fGenPar->RVol,2);
  IntVol *= fGenPar->HSeaWater*fGenPar->RhoSW+fGenPar->HRock*fGenPar->RhoSR;

  double Zmin=fGenPar->Can1-fGenPar->HRock+fGenPar->OriginZ;
  double Zmax=fGenPar->Can1+fGenPar->HSeaWater+fGenPar->OriginZ;
  sprintf(ch,"%10.3f %10.3f %10.3f %12.2E %12.2E",Zmin,Zmax,fGenPar->RVol,IntVol,fGenPar->NTot);
  fEvt->taga("genvol",ch);

//  sprintf(ch,"%d %d", 0,(int)fGenPar->NTot );
//  fEvt->taga("norma",ch);

  for(unsigned iFlux=0; iFlux<fGenPar->FluxFiles.size(); iFlux++){
    ostringstream fline;
    fline<<fGenPar->FluxFiles[iFlux].NuType<<" "<<fGenPar->FluxFiles[iFlux].FluxSimul<<" ";
    for(unsigned iFile=0; iFile<fGenPar->FluxFiles[iFlux].FileNames.size(); iFile++)
      fline<<fGenPar->FluxFiles[iFlux].FileNames[iFile]<<" ";

    string flux=fline.str();
    fEvt->taga("flux",flux);
  }

  fEvt->write(fOutFile);

  return;
}
//________________________________________________________________________________________

void EvtWrite::WriteEvent(GSeaEvent * SeaEvt){


  char ch[1024];

  int Geant3Code;
  int PdgCode;
  int WriteCode;

  if(fNewEvt)fEvt->defeve(SeaEvt->iEvt,1);

  if(SeaEvt->Neutrino.Id!=0){

  sprintf(ch, "%d %12.3f %12.3f %12.3f %12.6f %12.6f %12.6f %12.4E %12.6f %12.6f %12.6f %d %d %d",
               SeaEvt->Neutrino.Id,
               SeaEvt->Neutrino.V1,
               SeaEvt->Neutrino.V2,
               SeaEvt->Neutrino.V3,
               SeaEvt->Neutrino.D1,
               SeaEvt->Neutrino.D2,
               SeaEvt->Neutrino.D3,
               SeaEvt->Neutrino.E,
               SeaEvt->Neutrino.T,
               SeaEvt->Bx,SeaEvt->By,SeaEvt->ScattId,
               SeaEvt->Neutrino.PdgCode,
               SeaEvt->InterId);
  fEvt->taga("neutrino",ch);

  // if neutral current, neutrinos have pdg code in the primarylepton code (and not geant3 code)
  PdgCode=SeaEvt->PrimaryLepton.PdgCode;
  Geant3Code=fPdg->ConvertPdgToGeant3(PdgCode);

  // if neutral current, neutrinos have pdg code in the primarylepton code (and not geant3 code)
  if(fOptPDG || SeaEvt->InterId==3)WriteCode = PdgCode;
  else WriteCode = Geant3Code;

  sprintf(ch, "%d %12.3f %12.3f %12.3f %12.6f %12.6f %12.6f %12.4E %12.6f %d",
               SeaEvt->PrimaryLepton.Id,
               SeaEvt->PrimaryLepton.V1,
               SeaEvt->PrimaryLepton.V2,
               SeaEvt->PrimaryLepton.V3,
               SeaEvt->PrimaryLepton.D1,
               SeaEvt->PrimaryLepton.D2,
               SeaEvt->PrimaryLepton.D3,
               SeaEvt->PrimaryLepton.E,
               SeaEvt->PrimaryLepton.T,
               WriteCode);

  fEvt->taga("primarylepton",ch);

  sprintf(ch, "%d %d",SeaEvt->TargetZ,SeaEvt->TargetA);
  fEvt->taga("target",ch);

  sprintf(ch, "%d %d", SeaEvt->VerInCan, (int)SeaEvt->LepInCan);
  fEvt->taga("vertex",ch);

  sprintf(ch, "%12.3E %12.3E %12.3E", SeaEvt->Agen, SeaEvt->GenWeight, SeaEvt->EvtWeight);
  fEvt->taga("weights",ch);

  //only save the first 11 entries in w2list_gseagen.csv
  sprintf(ch, "%12.3E %12.3E %12.3E %12.3E %12.4E %12.3E %12.3E",
              fGenPar->GlobalGenWeight,pow(SeaEvt->Neutrino.E,fGenPar->Alpha),SeaEvt->XSecMean,SeaEvt->ColumnDepth,SeaEvt->PEarth,SeaEvt->WaterIntLen,SeaEvt->PScale);
  fEvt->taga("w2list",ch);

  }





  //we start from 1 because the first entry is the neutrino (or the cosmic ray) which has been alrady stored
  for(unsigned iTrack=1; iTrack<SeaEvt->Tracks.size(); iTrack++){

    PdgCode=SeaEvt->Tracks[iTrack].PdgCode;

    if(fOptPDG)WriteCode = PdgCode;
    else {
      // all neutrinos and anti-neutrinos have Geant3Code=4
      if(TMath::Abs(PdgCode)==12 || TMath::Abs(PdgCode)==14 || TMath::Abs(PdgCode)==16)Geant3Code=4;
      else Geant3Code=fPdg->ConvertPdgToGeant3(PdgCode);
      if(Geant3Code==0 && SeaEvt->Tracks[iTrack].Status==1 && PdgCode<1000000000)cout<<"No Geant3 Code for PDG "<<PdgCode<<endl;
      WriteCode = Geant3Code;
    }

    sprintf(ch, "%3d %12.3f %12.3f %12.3f %12.6f %12.6f %12.6f %12.4E %15.6f %d %12.6f",
                 SeaEvt->Tracks[iTrack].Id,
                 SeaEvt->Tracks[iTrack].V1,
                 SeaEvt->Tracks[iTrack].V2,
                 SeaEvt->Tracks[iTrack].V3,
                 SeaEvt->Tracks[iTrack].D1,
                 SeaEvt->Tracks[iTrack].D2,
                 SeaEvt->Tracks[iTrack].D3,
                 SeaEvt->Tracks[iTrack].E,
                 SeaEvt->Tracks[iTrack].T,
                 WriteCode,
                 SeaEvt->Tracks[iTrack].Length);

    if      (SeaEvt->Tracks[iTrack].Status<-1000)               fEvt->taga("track_earthlepton",ch); //propagated particles
    else if (SeaEvt->Tracks[iTrack].Status==1 && WriteCode!=0 ) fEvt->taga("track_in",ch);          //partciles that produce light with geant3 code
    else                                                        fEvt->taga("track_aux",ch);         //auxiluary particles

  }




  if(gSystem->Getenv("SYSTLIST")){
    ostringstream fline;
    fline<<SeaEvt->SysWgt.WList.size()<<" ";
    for(unsigned i=0; i<SeaEvt->SysWgt.WList.size(); i++)fline<<SeaEvt->SysWgt.WList[i]<<" ";
    fEvt->taga("wSysGen",fline.str());
    for(unsigned i=0; i<SeaEvt->SysWgt.WParam.size(); i++){
      fline.str("");
      fline<<SeaEvt->SysWgt.WParam[i].wght.size()<<" ";
      string tag="wSysGen_";
      for(unsigned j=0; j<SeaEvt->SysWgt.WParam[i].wght.size(); j++)fline<<SeaEvt->SysWgt.WParam[i].wght[j]<<" ";
        tag.append(SeaEvt->SysWgt.WParam[i].name);
        fEvt->taga(tag.c_str(),fline.str());
      }
  }



  if(fGenPar->GenMode.compare("POINT")==0){
    sprintf(ch, "%f", SeaEvt->LST);
    fEvt->taga("LST",ch);
    if(fGenPar->GenMJD){
      sprintf(ch, "%f ", SeaEvt->MJD);
      fEvt->taga("MJD",ch);
    }
  }

  sprintf(ch,"%d %d",SeaEvt->EvtTime[0],SeaEvt->EvtTime[1]);
  fEvt->taga("eventtime",ch);


  fEvt->write(fOutFile);

  return;
}

#endif