//____________________________________________________________________________
/*!

\class   genie::SeaNtpWriter

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

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

\created November 25, 2015

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

*/
//____________________________________________________________________________



#include "OutputWriters/SeaNtpWriter.h"

using namespace std;
using namespace genie;

//___________________________________________________________________________
SeaNtpWriter::SeaNtpWriter(GenParam * GenPar, string Format)
{

  fGenPar = GenPar;
  fSeaEvt = new GSeaEvent();

  fTreeFormat=Format;

  if((fTreeFormat.compare("EventTree")==0) || fTreeFormat.compare("SimpleTree")==0)
    LOG("SeaNtpWriter", pINFO) << "Using Tree Format "<<fTreeFormat;
  else {
    LOG("SeaNtpWriter", pFATAL) << "Unknowing Tree Format "<<fTreeFormat;
    exit(1);
  }

  this->Initialize(fGenPar->EvFilePrefix);
}
//___________________________________________________________________________
SeaNtpWriter::~SeaNtpWriter(){

  fTHeader->Write();
  fTEvent->Write();

  fOutFile->Close();

  if(fSeaEvt)delete fSeaEvt;
  if(fOutFile)delete fOutFile;

}
//___________________________________________________________________________
void SeaNtpWriter::Initialize(string EvFilePrefix){

  string ftstr;
  if(fTreeFormat.compare("EventTree")==0)ftstr=".et";
  if(fTreeFormat.compare("SimpleTree")==0)ftstr=".st";

  ostringstream fnstr;
  fnstr << EvFilePrefix  << "."
        << fGenPar->RunNu
        << ftstr
        << ".root";

  fOutFileName=fnstr.str();

  this->OpenFile();
  this->CreateTreeHeader();
  this->CreateTreeEvent();
  this->CreateEventBranch();

  return;
}
//___________________________________________________________________________
void SeaNtpWriter::OpenFile(){
  fOutFile = new TFile(fOutFileName.c_str(),"recreate");
}
//___________________________________________________________________________
void SeaNtpWriter::CreateTreeEvent(void){
  fTEvent = new TTree("Events","gSeaGen events");
  fTEvent->SetAutoSave(10000000);
}
//___________________________________________________________________________
void SeaNtpWriter::CreateEventBranch(){

  int BufSize= 32000;

  if(fTreeFormat.compare("EventTree")==0){
    fTEvent->Branch("Events","GSeaEvents",fSeaEvt);
  }
  else if(fTreeFormat.compare("SimpleTree")==0){
    fTEvent->Branch("iEvt",             &fSeaEvt->iEvt,                         "iEvt/I"                ,BufSize);
    fTEvent->Branch("PScale",           &fSeaEvt->PScale,                       "PScale/D"              ,BufSize);
    fTEvent->Branch("TargetA",          &fSeaEvt->TargetA,                      "TargetA/I"             ,BufSize);
    fTEvent->Branch("TargetZ",          &fSeaEvt->TargetZ,                      "TargetZ/I"             ,BufSize);
    fTEvent->Branch("InterId",          &fSeaEvt->InterId,                      "InterId/I"             ,BufSize);
    fTEvent->Branch("ScattId",          &fSeaEvt->ScattId,                      "ScattId/I"             ,BufSize);
    fTEvent->Branch("Bx",               &fSeaEvt->Bx,                           "Bx/D"                  ,BufSize);
    fTEvent->Branch("By",               &fSeaEvt->By,                           "By/D"                  ,BufSize);
    if(fGenPar->GenMode.compare("POINT")==0){
      fTEvent->Branch("LST",            &fSeaEvt->LST,                          "LST/D"                 ,BufSize);
      fTEvent->Branch("MJD",            &fSeaEvt->MJD,                          "LST/D"                 ,BufSize);
    }
    fTEvent->Branch("VerInCan",         &fSeaEvt->VerInCan,                     "VerInCan/I"            ,BufSize);
    fTEvent->Branch("LepInCan",         &fSeaEvt->LepInCan,                     "LepInCan/O"            ,BufSize);
    fTEvent->Branch("WaterXSec",        &fSeaEvt->WaterXSec,                    "WaterXSec/D"           ,BufSize);
    fTEvent->Branch("WaterIntLen",      &fSeaEvt->WaterIntLen,                  "WaterIntLen/D"         ,BufSize);
    fTEvent->Branch("PEarth",           &fSeaEvt->PEarth,                       "PEarth/D"              ,BufSize);
    fTEvent->Branch("ColumnDepth",      &fSeaEvt->ColumnDepth,                  "ColumnDepth/D"         ,BufSize);
    fTEvent->Branch("XSecMean",         &fSeaEvt->XSecMean,                     "XSecMean/D"            ,BufSize);
    fTEvent->Branch("Agen",             &fSeaEvt->Agen,                         "Agen/D"                ,BufSize);
    fTEvent->Branch("WGen",             &fSeaEvt->GenWeight,                    "WGen/D"                ,BufSize);
    fTEvent->Branch("WEvt",             &fSeaEvt->EvtWeight,                    "WEvt/D"                ,BufSize);

    fTEvent->Branch("E_nu",             &fSeaEvt->Neutrino.E,                   "E_nu/D"                ,BufSize);
    fTEvent->Branch("Pdg_nu",           &fSeaEvt->Neutrino.PdgCode,             "Pdg_nu/I"              ,BufSize);
    fTEvent->Branch("Vx_nu",            &fSeaEvt->Neutrino.V1,                  "Vx_nu/D"               ,BufSize);
    fTEvent->Branch("Vy_nu",            &fSeaEvt->Neutrino.V2,                  "Vy_nu/D"               ,BufSize);
    fTEvent->Branch("Vz_nu",            &fSeaEvt->Neutrino.V3,                  "Vz_nu/D"               ,BufSize);
    fTEvent->Branch("Dx_nu",            &fSeaEvt->Neutrino.D1,                  "Dx_nu/D"               ,BufSize);
    fTEvent->Branch("Dy_nu",            &fSeaEvt->Neutrino.D2,                  "Dy_nu/D"               ,BufSize);
    fTEvent->Branch("Dz_nu",            &fSeaEvt->Neutrino.D3,                  "Dz_nu/D"               ,BufSize);
    fTEvent->Branch("T_nu",             &fSeaEvt->Neutrino.T,                   "T_nu/D"                ,BufSize);


    if(fGenPar->PropMode){
      fTEvent->Branch("NInter",         &fSeaEvt->NInter,                       "NInter/I"              ,BufSize);
      fTEvent->Branch("E_nui",          &fSeaEvt->EarthNeutrino.E,              "E_nui/D"               ,BufSize);
      fTEvent->Branch("Pdg_nui",        &fSeaEvt->EarthNeutrino.PdgCode,        "Pdg_nui/I"             ,BufSize);
      fTEvent->Branch("Vx_nui",         &fSeaEvt->EarthNeutrino.V1,             "Vx_nui/D"              ,BufSize);
      fTEvent->Branch("Vy_nui",         &fSeaEvt->EarthNeutrino.V2,             "Vy_nui/D"              ,BufSize);
      fTEvent->Branch("Vz_nui",         &fSeaEvt->EarthNeutrino.V3,             "Vz_nui/D"              ,BufSize);
      fTEvent->Branch("Dx_nui",         &fSeaEvt->EarthNeutrino.D1,             "Dx_nui/D"              ,BufSize);
      fTEvent->Branch("Dy_nui",         &fSeaEvt->EarthNeutrino.D2,             "Dy_nui/D"              ,BufSize);
      fTEvent->Branch("Dz_nui",         &fSeaEvt->EarthNeutrino.D3,             "Dz_nui/D"              ,BufSize);
      fTEvent->Branch("T_nui",          &fSeaEvt->EarthNeutrino.T,              "T_nui/D"               ,BufSize);
    }


    fTEvent->Branch("E_pl",             &fSeaEvt->PrimaryLepton.E,              "E_pl/D"                ,BufSize);
    fTEvent->Branch("Pdg_pl",           &fSeaEvt->PrimaryLepton.PdgCode,        "Pdg_pl/I"              ,BufSize);
    fTEvent->Branch("Vx_pl",            &fSeaEvt->PrimaryLepton.V1,             "Vx_pl/D"               ,BufSize);
    fTEvent->Branch("Vy_pl",            &fSeaEvt->PrimaryLepton.V2,             "Vy_pl/D"               ,BufSize);
    fTEvent->Branch("Vz_pl",            &fSeaEvt->PrimaryLepton.V3,             "Vz_pl/D"               ,BufSize);
    fTEvent->Branch("Dx_pl",            &fSeaEvt->PrimaryLepton.D1,             "Dx_pl/D"               ,BufSize);
    fTEvent->Branch("Dy_pl",            &fSeaEvt->PrimaryLepton.D2,             "Dy_pl/D"               ,BufSize);
    fTEvent->Branch("Dz_pl",            &fSeaEvt->PrimaryLepton.D3,             "Dz_pl/D"               ,BufSize);
    fTEvent->Branch("T_pl",             &fSeaEvt->PrimaryLepton.T,              "T_pl/D"                ,BufSize);

    fTEvent->Branch("NTracks",          &fNTracks,                              "NTracks/I"             ,BufSize);

    fTEvent->Branch("Id_tr",            fTrackId,                               "Id_tr[NTracks]/I"      ,BufSize);
    fTEvent->Branch("Pdg_tr",           fTrackPdg,                              "Pdg_tr[NTracks]/I"     ,BufSize);
    if(fGenPar->WriteParticlesMode!=0){
      fTEvent->Branch("St_tr",          fTrackSt,                               "St_tr[NTracks]/I"      ,BufSize);
      fTEvent->Branch("MotherId_tr",    fTrackMotherId,                         "MotherId_tr[NTracks]/I",BufSize);
    }
    fTEvent->Branch("E_tr",             fTrackE,                                "E_tr[NTracks]/D"       ,BufSize);
    fTEvent->Branch("Vx_tr",            fTrackVx,                               "Vx_tr[NTracks]/D"      ,BufSize);
    fTEvent->Branch("Vy_tr",            fTrackVy,                               "Vy_tr[NTracks]/D"      ,BufSize);
    fTEvent->Branch("Vz_tr",            fTrackVz,                               "Vz_tr[NTracks]/D"      ,BufSize);
    fTEvent->Branch("Dx_tr",            fTrackDx,                               "Dx_tr[NTracks]/D"      ,BufSize);
    fTEvent->Branch("Dy_tr",            fTrackDy,                               "Dy_tr[NTracks]/D"      ,BufSize);
    fTEvent->Branch("Dz_tr",            fTrackDz,                               "Dz_tr[NTracks]/D"      ,BufSize);
    fTEvent->Branch("T_tr",             fTrackT,                                "T_tr[NTracks]/D"       ,BufSize);

    if(gSystem->Getenv("SYSTLIST")){
      fTEvent->Branch("NSysWgt",        &fNSysWgt,                              "NSysWgt/I"             ,BufSize);
      fTEvent->Branch("WSys",           fSysWgt,                                "WSys[NSysWgt]/D"       ,BufSize);

      for(unsigned i=0; i<fGenPar->SystParams.SingleParams.size(); i++){
        if(fGenPar->SystParams.SingleParams[i].name.compare("NormCCQE")==0){
          fTEvent->Branch("NSysWgt_NormCCQE",           &fNSysWgt_NormCCQE,      "NSysWgt_NormCCQE/I"                   ,BufSize);
          fTEvent->Branch("WSys_NormCCQE",              fSysWgt_NormCCQE,        "WSys_NormCCQE[NSysWgt_NormCCQE]/D"    ,BufSize);
        }
        else if(fGenPar->SystParams.SingleParams[i].name.compare("MaCCQEshape")==0){
          fTEvent->Branch("NSysWgt_MaCCQEshape",        &fNSysWgt_MaCCQEshape,    "NSysWgt_MaCCQEshape/I"                       ,BufSize);
          fTEvent->Branch("WSys_MaCCQEshape",           fSysWgt_MaCCQEshape,      "WSys_MaCCQEshape[NSysWgt_MaCCQEshape]/D"     ,BufSize);
        }
        else if(fGenPar->SystParams.SingleParams[i].name.compare("NormCCRES")==0){
          fTEvent->Branch("NSysWgt_NormCCRES",          &fNSysWgt_NormCCRES,    "NSysWgt_NormCCRES/I"                   ,BufSize);
          fTEvent->Branch("WSys_NormCCRES",             fSysWgt_NormCCRES,      "WSys_NormCCRES[NSysWgt_NormCCRES]/D"   ,BufSize);
        }
        else if(fGenPar->SystParams.SingleParams[i].name.compare("VecFFCCQEshape")==0){
          fTEvent->Branch("NSysWgt_VecFFCCQEshape",     &fNSysWgt_VecFFCCQEshape,    "NSysWgt_VecFFCCQEshape/I"                 ,BufSize);
          fTEvent->Branch("WSys_VecFFCCQEshape",        fSysWgt_VecFFCCQEshape,      "WSys_VecFFCCQEshape[NSysWgt_VecFFCCQEshape]/D"    ,BufSize);
        }
        else if(fGenPar->SystParams.SingleParams[i].name.compare("MaCCRESshape")==0){
          fTEvent->Branch("NSysWgt_MaCCRESshape",       &fNSysWgt_MaCCRESshape,    "NSysWgt_MaCCRESshape/I"                     ,BufSize);
          fTEvent->Branch("WSys_MaCCRESshape",          fSysWgt_MaCCRESshape,      "WSys_MaCCRESshape[NSysWgt_MaCCRESshape]/D"  ,BufSize);
        }
        else if(fGenPar->SystParams.SingleParams[i].name.compare("MvCCRESshape")==0){
          fTEvent->Branch("NSysWgt_MvCCRESshape",       &fNSysWgt_MvCCRESshape,    "NSysWgt_MvCCRESshape/I"                     ,BufSize);
          fTEvent->Branch("WSys_MvCCRESshape",          fSysWgt_MvCCRESshape,      "WSys_MvCCRESshape[NSysWgt_MvCCRESshape]/D"  ,BufSize);
        }
        else if(fGenPar->SystParams.SingleParams[i].name.compare("NormNCRES")==0){
          fTEvent->Branch("NSysWgt_NormNCRES",          &fNSysWgt_NormNCRES,    "NSysWgt_NormNCRES/I"                   ,BufSize);
          fTEvent->Branch("WSys_NormNCRES",             fSysWgt_NormNCRES,      "WSys_NormNCRES[NSysWgt_NormNCRES]/D"   ,BufSize);
        }
        else if(fGenPar->SystParams.SingleParams[i].name.compare("MaNCRESshape")==0){
          fTEvent->Branch("NSysWgt_MaNCRESshape",       &fNSysWgt_MaNCRESshape,    "NSysWgt_MaNCRESshape/I"                     ,BufSize);
          fTEvent->Branch("WSys_MaNCRESshape",          fSysWgt_MaNCRESshape,      "WSys_MaNCRESshape[NSysWgt_MaNCRESshape]/D"  ,BufSize);
        }
        else if(fGenPar->SystParams.SingleParams[i].name.compare("MvNCRESshape")==0){
          fTEvent->Branch("NSysWgt_MvNCRESshape",       &fNSysWgt_MvNCRESshape,    "NSysWgt_MvNCRESshape/I"                     ,BufSize);
          fTEvent->Branch("WSys_MvNCRESshape",          fSysWgt_MvNCRESshape,      "WSys_MvNCRESshape[NSysWgt_MvNCRESshape]/D"  ,BufSize);
        }
        else if(fGenPar->SystParams.SingleParams[i].name.compare("NonRESBGvpCC1pi")==0){
          fTEvent->Branch("NSysWgt_NonRESBGvpCC1pi",    &fNSysWgt_NonRESBGvpCC1pi,    "NSysWgt_NonRESBGvpCC1pi/I"                       ,BufSize);
          fTEvent->Branch("WSys_NonRESBGvpCC1pi",       fSysWgt_NonRESBGvpCC1pi,      "WSys_NormNCRES[NSysWgt_NonRESBGvpCC1pi]/D"       ,BufSize);
        }
        else if(fGenPar->SystParams.SingleParams[i].name.compare("NonRESBGvnCC1pi")==0){
          fTEvent->Branch("NSysWgt_NonRESBGvnCC1pi",     &fNSysWgt_NonRESBGvnCC1pi,    "NSysWgt_NonRESBGvnCC1pi/I"                      ,BufSize);
          fTEvent->Branch("WSys_NonRESBGvnCC1pi",       fSysWgt_NonRESBGvnCC1pi,      "WSys_NonRESBGvnCC1pi[NSysWgt_NonRESBGvnCC1pi]/D" ,BufSize);
        }
        else if(fGenPar->SystParams.SingleParams[i].name.compare("MaCOHpi")==0){
          fTEvent->Branch("NSysWgt_MaCOHpi",            &fNSysWgt_MaCOHpi,    "NSysWgt_MaCOHpi/I"                       ,BufSize);
          fTEvent->Branch("WSys_MaCOHpi",               fSysWgt_MaCOHpi,      "WSys_MaCOHpi[NSysWgt_MaCOHpi]/D"         ,BufSize);
        }
        else if(fGenPar->SystParams.SingleParams[i].name.compare("MFP_pi")==0){
          fTEvent->Branch("NSysWgt_MFP_pi",             &fNSysWgt_MFP_pi,    "NSysWgt_MFP_pi/I"                 ,BufSize);
          fTEvent->Branch("WSys_MFP_pi",                fSysWgt_MFP_pi,      "WSys_MFP_pi[NSysWgt_MFP_pi]/D"    ,BufSize);
        }
        else if(fGenPar->SystParams.SingleParams[i].name.compare("MFP_N")==0){
          fTEvent->Branch("NSysWgt_MFP_N",              &fNSysWgt_MFP_N,    "NSysWgt_MFP_N/I"                   ,BufSize);
          fTEvent->Branch("WSys_MFP_N",                 fSysWgt_MFP_N,      "WSys_MFP_N[NSysWgt_MFP_N]/D"       ,BufSize);
        }
        else if(fGenPar->SystParams.SingleParams[i].name.compare("FrPiProd_pi")==0){
          fTEvent->Branch("NSysWgt_FrPiProd_pi",        &fNSysWgt_FrPiProd_pi,    "NSysWgt_FrPiProd_pi/I"                       ,BufSize);
          fTEvent->Branch("WSys_FrPiProd_pi",           fSysWgt_FrPiProd_pi,      "WSys_FrPiProd_pi[NSysWgt_FrPiProd_pi]/D"     ,BufSize);
        }
        else if(fGenPar->SystParams.SingleParams[i].name.compare("AGKYxF1pi")==0){
          fTEvent->Branch("NSysWgt_AGKYxF1pi",          &fNSysWgt_AGKYxF1pi,    "NSysWgt_AGKYxF1pi/I"                   ,BufSize);
          fTEvent->Branch("WSys_AGKYxF1pi",             fSysWgt_AGKYxF1pi,      "WSys_AGKYxF1pi[NSysWgt_AGKYxF1pi]/D"   ,BufSize);
        }
        else if(fGenPar->SystParams.SingleParams[i].name.compare("AGKYpT1pi")==0){
          fTEvent->Branch("NSysWgt_AGKYpT1pi",          &fNSysWgt_AGKYpT1pi,    "NSysWgt_AGKYpT1pi/I"                   ,BufSize);
          fTEvent->Branch("WSys_AGKYpT1pi",             fSysWgt_AGKYpT1pi,      "WSys_AGKYpT1pi[NSysWgt_AGKYpT1pi]/D"   ,BufSize);
        }

        else if(fGenPar->SystParams.SingleParams[i].name.compare("FormZone")==0){
          fTEvent->Branch("NSysWgt_FormZone",           &fNSysWgt_FormZone,    "NSysWgt_FormZone/I"                     ,BufSize);
          fTEvent->Branch("WSys_FormZone",              fSysWgt_FormZone,      "WSys_FormZone[NSysWgt_FormZone]/D"      ,BufSize);
        }
        else if(fGenPar->SystParams.SingleParams[i].name.compare("Theta_Delta2Npi")==0){
          fTEvent->Branch("NSysWgt_Theta_Delta2Npi",    &fNSysWgt_Theta_Delta2Npi,    "NSysWgt_Theta_Delta2Npi/I"                       ,BufSize);
          fTEvent->Branch("WSys_Theta_Delta2Npi",       fSysWgt_Theta_Delta2Npi,      "WSys_Theta_Delta2Npi[NSysWgt_Theta_Delta2Npi]/D" ,BufSize);
        }
      }

    }
  }

  return;
}
//___________________________________________________________________________
void SeaNtpWriter::CreateTreeHeader(void){

  fTHeader = new TTree("Header"," ");

  if(fTreeFormat.compare("EventTree")==0){
    fTHeader->Branch("Header","GenParam",&fGenPar);
  }
  else if(fTreeFormat.compare("SimpleTree")==0){

    int BufSize= 32000;
    int NCHAR=100;

    char gSeaGenVer[NCHAR], GenieVer[NCHAR], InpXSecFile[NCHAR],EventGeneratorList[NCHAR];
    char PropCode[NCHAR], GenMode[NCHAR];
    sprintf(gSeaGenVer,"%s",fGenPar->gSeaGenVer.c_str());
    sprintf(GenieVer,"%s",fGenPar->GenieVer.c_str());

    sprintf(EventGeneratorList,"%s",fGenPar->EventGeneratorList.c_str());
    sprintf(InpXSecFile,"%s",fGenPar->InpXSecFile.c_str());
    sprintf(PropCode,"%s",fGenPar->PropCode.c_str());
    sprintf(GenMode,"%s",fGenPar->GenMode.c_str());

    char* rundate=asctime(localtime (&fGenPar->RunTime));

    fTHeader->Branch("gSeaGenVer",      gSeaGenVer,                     "gSeaGenVar/C"          ,BufSize);
    fTHeader->Branch("GenieVer",        GenieVer,                       "GenieVer/C"            ,BufSize);
    fTHeader->Branch("RunTime",         rundate,                        "RunTime/C"             ,BufSize);
    fTHeader->Branch("RunNu",           &fGenPar->RunNu,                "RunNu/I"               ,BufSize);
    fTHeader->Branch("RanSeed",         &fGenPar->RanSeed,              "RanSeed/I"             ,BufSize);

    fTHeader->Branch("EventGeneratorList", EventGeneratorList,          "EventGeneratorList/C"  ,BufSize);

    fTHeader->Branch("InpXSecFile",     InpXSecFile,                    "InpXSecFile/C"         ,BufSize);
    fTHeader->Branch("NTot",            &fGenPar->NTot,                 "NTot/D"                ,BufSize);
    fTHeader->Branch("EvMin",           &fGenPar->EvMin,                "EvMin/D"               ,BufSize);
    fTHeader->Branch("EvMax",           &fGenPar->EvMax,                "EvMax/D"               ,BufSize);
    fTHeader->Branch("CtMin",           &fGenPar->CtMin,                "CtMin/D"               ,BufSize);
    fTHeader->Branch("CtMax",           &fGenPar->CtMax,                "CtMax/D"               ,BufSize);
    fTHeader->Branch("Alpha",           &fGenPar->Alpha,                "Alpha/D"               ,BufSize);
    if ( fGenPar->GenMode.compare("BIN")==0 ){
      fTHeader->Branch("Primary",           &fGenPar->Primary,                "Primary/I"               ,BufSize);
      fTHeader->Branch("EPrimMin",          &fGenPar->EPrimMin,               "EPrimMin/D"              ,BufSize);
      fTHeader->Branch("EPrimMax",          &fGenPar->EPrimMax,               "EPrimMax/D"              ,BufSize);
      fTHeader->Branch("PrimCtMin",         &fGenPar->PrimCtMin,              "PrimCtMin/D"             ,BufSize);
      fTHeader->Branch("PrimCtMax",         &fGenPar->PrimCtMax,              "PrimCtMax/D"             ,BufSize);
      fTHeader->Branch("CRModel",           &fGenPar->CRModel,                "CRModel/O"               ,BufSize);
      fTHeader->Branch("EmuMin",            &fGenPar->EmuMin,                 "EmuMin/D"                ,BufSize);
    }
    fTHeader->Branch("NBin",            &fGenPar->NBin,                 "NBin/I"                ,BufSize);
    if(fGenPar->detector.size()>0){
      char detector[NCHAR];
      sprintf(detector,"%s",fGenPar->detector.c_str());
      fTHeader->Branch("detector",        detector,                     "detector/C"            ,BufSize);
    }
    fTHeader->Branch("Can1",            &fGenPar->Can1,                 "Can1/D"                ,BufSize);
    fTHeader->Branch("Can2",            &fGenPar->Can2,                 "Can2/D"                ,BufSize);
    fTHeader->Branch("Can3",            &fGenPar->Can3,                 "Can3/D"                ,BufSize);
    fTHeader->Branch("OriginX",         &fGenPar->OriginX,              "OriginX/D"             ,BufSize);
    fTHeader->Branch("OriginY",         &fGenPar->OriginY,              "OriginY/D"             ,BufSize);
    fTHeader->Branch("OriginZ",         &fGenPar->OriginZ,              "OriginZ/D"             ,BufSize);
    fTHeader->Branch("HRock",           &fGenPar->HRock,                "HRock/D"               ,BufSize);
    fTHeader->Branch("HSeaWater",       &fGenPar->HSeaWater,            "HSeaWater/D"           ,BufSize);
    fTHeader->Branch("RVol",            &fGenPar->RVol,                 "RVol/D"                ,BufSize);
#ifdef _ANTARES_ENABLED__
    fTHeader->Branch("HBedRock",        &fGenPar->HBedRock,             "HBedRock/D"            ,BufSize);
    fTHeader->Branch("NAbsLength",      &fGenPar->NAbsLength,           "NAbsLength/D"          ,BufSize);
    fTHeader->Branch("AbsLength",       &fGenPar->AbsLength,            "AbsLength/D"           ,BufSize);
#endif
    fTHeader->Branch("SiteDepth",       &fGenPar->SiteDepth,            "SiteDepth/D"           ,BufSize);
    fTHeader->Branch("SiteLatitude",    &fGenPar->SiteLatitude,         "SiteLatitude/D"        ,BufSize);
    fTHeader->Branch("SiteLongitude",   &fGenPar->SiteLongitude,        "SiteLongitude/D"       ,BufSize);
    fTHeader->Branch("SeaBottomRadius", &fGenPar->SeaBottomRadius,       "SeaBottomRadius/D"    ,BufSize);
    fTHeader->Branch("GlobalGenWeight", &fGenPar->GlobalGenWeight,       "GlobalGenWeight/D"    ,BufSize);
    fTHeader->Branch("RhoSW",           &fGenPar->RhoSW,                "RhoSW/D"               ,BufSize);
    fTHeader->Branch("RhoSR",           &fGenPar->RhoSR,                "RhoSR/D"               ,BufSize);
    fTHeader->Branch("TGen",            &fGenPar->TGen,                 "TGen/D"                ,BufSize);

    if(fGenPar->TrackEvts>0){
      fTHeader->Branch("PropCode",        PropCode,                     "PropCode/C"            ,BufSize);
    }

    fTHeader->Branch("GenMode",         GenMode,                        "GenMode/C"             ,BufSize);
    if(fGenPar->GenMode.compare("POINT")==0){
      char SourceFile[NCHAR], SourceName[NCHAR];
      sprintf(SourceFile,"%s",fGenPar->SourceFile.c_str());
      sprintf(SourceName,"%s",fGenPar->SourceName.c_str());
      fTHeader->Branch("SourceFile",      SourceFile,                   "SourceFile/C"          ,BufSize);
      fTHeader->Branch("SourceName",      SourceName,                   "SourceName/C"          ,BufSize);
      fTHeader->Branch("Declination",     &fGenPar->Declination,        "Declination/D"         ,BufSize);
      fTHeader->Branch("RightAscension",  &fGenPar->RightAscension,     "RightAscension/D"      ,BufSize);
      fTHeader->Branch("SourceRadius",    &fGenPar->SourceRadius,       "SourceRadius/D"        ,BufSize);
      fTHeader->Branch("GenMJD",          &fGenPar->GenMJD,             "GenMJD/O"              ,BufSize);
      if(fGenPar->GenMJD){
        fTHeader->Branch("MJDStart",      &fGenPar->MJDStart,           "MJDStart/D"            ,BufSize);
        fTHeader->Branch("MJDStop",       &fGenPar->MJDStop,            "MJDStop/D"             ,BufSize);
      }
    }
    fTHeader->Branch("PropMode",          &fGenPar->PropMode,           "PropMode/O"            ,BufSize);


    int NNu=fGenPar->NuList.size();
    fTHeader->Branch("NNu",               &NNu,                         "NNu/I"                 ,BufSize);

    int NuList[6] = {0,0,0,0,0,0};
    for(unsigned i=0; i<fGenPar->NuList.size(); i++)NuList[i]=fGenPar->NuList[i];
    fTHeader->Branch("NuList",            NuList,                       "NuList[NNu]/I"         ,BufSize);

    char Flux[6][1000];
    for(unsigned iFlux=0; iFlux<fGenPar->FluxFiles.size(); iFlux++){
      ostringstream fline;
      fline<<fGenPar->FluxFiles[iFlux].NuType<<" "<<fGenPar->FluxFiles[iFlux].FluxSimul<<" ";
      char neu[10];
      sprintf(neu,"Flux%d",iFlux+1);
      string neu1(neu);
      string neu2=neu1;
      neu2.append("/C");
      char flux[1000];
      for(unsigned iFile=0; iFile<fGenPar->FluxFiles[iFlux].FileNames.size(); iFile++){
        fline<<fGenPar->FluxFiles[iFlux].FileNames[iFile]<<" ";
        sprintf(flux,"%s",fline.str().c_str());
      }
      strcpy(Flux[iFlux],flux);
      fTHeader->Branch(neu1.c_str(),      Flux[iFlux],                  neu2.c_str()            ,BufSize);
    }

    fTHeader->Branch("RangeSW",         &fGenPar->RangeSW);
    fTHeader->Branch("RangeSR",         &fGenPar->RangeSR);

  }

  fTHeader->Fill();
}
//________________________________________________________________________________________
void SeaNtpWriter::WriteEvent(GSeaEvent * SeaEvt){

  fSeaEvt->Copy(SeaEvt);

  if(fTreeFormat.compare("SimpleTree")==0){

    fNTracks=0;

    for(int j=0; j<NTOTMAX; j++) {
      fTrackId[j]        = 0;
      fTrackSt[j]        = 0;
      fTrackPdg[j]       = 0;
      fTrackMotherId[j]  = 0;
      fTrackE[j]         = 0.;
      fTrackVx[j]        = 0.;
      fTrackVy[j]        = 0.;
      fTrackVz[j]        = 0.;
      fTrackDx[j]        = 0.;
      fTrackDy[j]        = 0.;
      fTrackDz[j]        = 0.;
      fTrackT[j]         = 0.;
    }

    for(unsigned iTr=0; iTr<fSeaEvt->Tracks.size();iTr++){
      // we write only particles at the can (or in the bedrock if defined)
      if(fSeaEvt->Tracks[iTr].PdgCode<1000000000 && SeaEvt->Tracks[iTr].Status>0 ){
        fTrackId[fNTracks]        = fSeaEvt->Tracks[iTr].Id;
        fTrackSt[fNTracks]        = fSeaEvt->Tracks[iTr].Status;
        fTrackPdg[fNTracks]       = fSeaEvt->Tracks[iTr].PdgCode;
        fTrackMotherId[fNTracks]  = fSeaEvt->Tracks[iTr].MotherId;
        fTrackE[fNTracks]         = fSeaEvt->Tracks[iTr].E;
        fTrackVx[fNTracks]        = fSeaEvt->Tracks[iTr].V1;
        fTrackVy[fNTracks]        = fSeaEvt->Tracks[iTr].V2;
        fTrackVz[fNTracks]        = fSeaEvt->Tracks[iTr].V3;
        fTrackDx[fNTracks]        = fSeaEvt->Tracks[iTr].D1;
        fTrackDy[fNTracks]        = fSeaEvt->Tracks[iTr].D2;
        fTrackDz[fNTracks]        = fSeaEvt->Tracks[iTr].D3;
        fTrackT[fNTracks]         = fSeaEvt->Tracks[iTr].T;
        fNTracks++;
      }
    }



    if(gSystem->Getenv("SYSTLIST")){
      fNSysWgt=fGenPar->SystParams.ListParams.size();
      for(int j=0; j<NSYSMAX; j++) {
        fSysWgt[j] = 0.;
        fSysWgt_NormCCQE[j] = 0.;
        fSysWgt_MaCCQEshape[j] = 0.;
        fSysWgt_NormCCRES[j] = 0.;
        fSysWgt_VecFFCCQEshape[j] = 0.;
        fSysWgt_MaCCRESshape[j] = 0.;
        fSysWgt_MvCCRESshape[j] = 0.;
        fSysWgt_NormNCRES[j] = 0.;
        fSysWgt_MaNCRESshape[j] = 0.;
        fSysWgt_MvNCRESshape[j] = 0.;
        fSysWgt_NonRESBGvpCC1pi[j] = 0.;
        fSysWgt_NonRESBGvnCC1pi[j] = 0.;
        fSysWgt_MaCOHpi[j] = 0.;
        fSysWgt_MFP_pi[j] = 0.;
        fSysWgt_MFP_N[j] = 0.;
        fSysWgt_FrPiProd_pi[j] = 0.;
        fSysWgt_AGKYxF1pi[j] = 0.;
        fSysWgt_AGKYpT1pi[j] = 0.;
        fSysWgt_FormZone[j] = 0.;
        fSysWgt_Theta_Delta2Npi[j] = 0.;
      }
      for(int j=0; j<fNSysWgt; j++) {
        fSysWgt[j] = fSeaEvt->SysWgt.WList[j];
      }

      for(unsigned i=0; i<SeaEvt->SysWgt.WParam.size(); i++){
        if(fSeaEvt->SysWgt.WParam[i].name.compare("NormCCQE")==0){
          fNSysWgt_NormCCQE=fSeaEvt->SysWgt.WParam[i].wght.size();
          for(unsigned j=0; j<SeaEvt->SysWgt.WParam[i].wght.size(); j++)fSysWgt_NormCCQE[j]=SeaEvt->SysWgt.WParam[i].wght[j];
        }
        else if(fSeaEvt->SysWgt.WParam[i].name.compare("MaCCQEshape")==0){
          fNSysWgt_MaCCQEshape=fSeaEvt->SysWgt.WParam[i].wght.size();
          for(unsigned j=0; j<SeaEvt->SysWgt.WParam[i].wght.size(); j++)fSysWgt_MaCCQEshape[j]=SeaEvt->SysWgt.WParam[i].wght[j];
        }
        else if(fSeaEvt->SysWgt.WParam[i].name.compare("NormCCRES")==0){
          fNSysWgt_NormCCRES=fSeaEvt->SysWgt.WParam[i].wght.size();
          for(unsigned j=0; j<SeaEvt->SysWgt.WParam[i].wght.size(); j++)fSysWgt_NormCCRES[j]=SeaEvt->SysWgt.WParam[i].wght[j];
        }
        else if(fSeaEvt->SysWgt.WParam[i].name.compare("VecFFCCQEshape")==0){
          fNSysWgt_VecFFCCQEshape=fSeaEvt->SysWgt.WParam[i].wght.size();
          for(unsigned j=0; j<SeaEvt->SysWgt.WParam[i].wght.size(); j++)fSysWgt_VecFFCCQEshape[j]=SeaEvt->SysWgt.WParam[i].wght[j];
        }
        else if(fSeaEvt->SysWgt.WParam[i].name.compare("MaCCRESshape")==0){
          fNSysWgt_MaCCRESshape=fSeaEvt->SysWgt.WParam[i].wght.size();
          for(unsigned j=0; j<SeaEvt->SysWgt.WParam[i].wght.size(); j++)fSysWgt_MaCCRESshape[j]=SeaEvt->SysWgt.WParam[i].wght[j];
        }
        else if(fSeaEvt->SysWgt.WParam[i].name.compare("MvCCRESshape")==0){
          fNSysWgt_MvCCRESshape=fSeaEvt->SysWgt.WParam[i].wght.size();
          for(unsigned j=0; j<SeaEvt->SysWgt.WParam[i].wght.size(); j++)fSysWgt_MvCCRESshape[j]=SeaEvt->SysWgt.WParam[i].wght[j];
        }
        else if(fSeaEvt->SysWgt.WParam[i].name.compare("NormNCRES")==0){
          fNSysWgt_NormNCRES=fSeaEvt->SysWgt.WParam[i].wght.size();
          for(unsigned j=0; j<SeaEvt->SysWgt.WParam[i].wght.size(); j++)fSysWgt_NormNCRES[j]=SeaEvt->SysWgt.WParam[i].wght[j];
        }
        else if(fSeaEvt->SysWgt.WParam[i].name.compare("MaNCRESshape")==0){
          fNSysWgt_MaNCRESshape=fSeaEvt->SysWgt.WParam[i].wght.size();
          for(unsigned j=0; j<SeaEvt->SysWgt.WParam[i].wght.size(); j++)fSysWgt_MaNCRESshape[j]=SeaEvt->SysWgt.WParam[i].wght[j];
        }
        else if(fSeaEvt->SysWgt.WParam[i].name.compare("MvNCRESshape")==0){
          fNSysWgt_MvNCRESshape=fSeaEvt->SysWgt.WParam[i].wght.size();
          for(unsigned j=0; j<SeaEvt->SysWgt.WParam[i].wght.size(); j++)fSysWgt_MvNCRESshape[j]=SeaEvt->SysWgt.WParam[i].wght[j];
        }
        else if(fSeaEvt->SysWgt.WParam[i].name.compare("NonRESBGvpCC1pi")==0){
          fNSysWgt_NonRESBGvpCC1pi=fSeaEvt->SysWgt.WParam[i].wght.size();
          for(unsigned j=0; j<SeaEvt->SysWgt.WParam[i].wght.size(); j++)fSysWgt_NonRESBGvpCC1pi[j]=SeaEvt->SysWgt.WParam[i].wght[j];
        }
        else if(fSeaEvt->SysWgt.WParam[i].name.compare("NonRESBGvnCC1pi")==0){
          fNSysWgt_NonRESBGvnCC1pi=fSeaEvt->SysWgt.WParam[i].wght.size();
          for(unsigned j=0; j<SeaEvt->SysWgt.WParam[i].wght.size(); j++)fSysWgt_NonRESBGvnCC1pi[j]=SeaEvt->SysWgt.WParam[i].wght[j];
        }
        else if(fSeaEvt->SysWgt.WParam[i].name.compare("MaCOHpi")==0){
          fNSysWgt_MaCOHpi=fSeaEvt->SysWgt.WParam[i].wght.size();
          for(unsigned j=0; j<SeaEvt->SysWgt.WParam[i].wght.size(); j++)fSysWgt_MaCOHpi[j]=SeaEvt->SysWgt.WParam[i].wght[j];
        }
        else if(fSeaEvt->SysWgt.WParam[i].name.compare("MFP_pi")==0){
          fNSysWgt_MFP_pi=fSeaEvt->SysWgt.WParam[i].wght.size();
          for(unsigned j=0; j<SeaEvt->SysWgt.WParam[i].wght.size(); j++)fSysWgt_MFP_pi[j]=SeaEvt->SysWgt.WParam[i].wght[j];
        }
        else if(fSeaEvt->SysWgt.WParam[i].name.compare("MFP_N")==0){
          fNSysWgt_MFP_N=fSeaEvt->SysWgt.WParam[i].wght.size();
          for(unsigned j=0; j<SeaEvt->SysWgt.WParam[i].wght.size(); j++)fSysWgt_MFP_N[j]=SeaEvt->SysWgt.WParam[i].wght[j];
        }
        else if(fSeaEvt->SysWgt.WParam[i].name.compare("FrPiProd_pi")==0){
          fNSysWgt_FrPiProd_pi=fSeaEvt->SysWgt.WParam[i].wght.size();
          for(unsigned j=0; j<SeaEvt->SysWgt.WParam[i].wght.size(); j++)fSysWgt_FrPiProd_pi[j]=SeaEvt->SysWgt.WParam[i].wght[j];
        }
        else if(fSeaEvt->SysWgt.WParam[i].name.compare("AGKYxF1pi")==0){
          fNSysWgt_AGKYxF1pi=fSeaEvt->SysWgt.WParam[i].wght.size();
          for(unsigned j=0; j<SeaEvt->SysWgt.WParam[i].wght.size(); j++)fSysWgt_AGKYxF1pi[j]=SeaEvt->SysWgt.WParam[i].wght[j];
        }
        else if(fSeaEvt->SysWgt.WParam[i].name.compare("AGKYpT1pi")==0){
          fNSysWgt_AGKYpT1pi=fSeaEvt->SysWgt.WParam[i].wght.size();
          for(unsigned j=0; j<SeaEvt->SysWgt.WParam[i].wght.size(); j++)fSysWgt_AGKYpT1pi[j]=SeaEvt->SysWgt.WParam[i].wght[j];
        }
        else if(fSeaEvt->SysWgt.WParam[i].name.compare("FormZone")==0){
          fNSysWgt_FormZone=fSeaEvt->SysWgt.WParam[i].wght.size();
          for(unsigned j=0; j<SeaEvt->SysWgt.WParam[i].wght.size(); j++)fSysWgt_FormZone[j]=SeaEvt->SysWgt.WParam[i].wght[j];
        }
        else if(fSeaEvt->SysWgt.WParam[i].name.compare("Theta_Delta2Npi")==0){
          fNSysWgt_Theta_Delta2Npi=fSeaEvt->SysWgt.WParam[i].wght.size();
          for(unsigned j=0; j<SeaEvt->SysWgt.WParam[i].wght.size(); j++)fSysWgt_Theta_Delta2Npi[j]=SeaEvt->SysWgt.WParam[i].wght[j];
        }
      }

    }
  }
  fTEvent->Fill();

  return;
}