//_____________________________________________________________________________
/*!

\class   genie::GSyst

\brief   A class to evaluate the systematic weights using the GENIE class GReWeight

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

\created May 14, 2015

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

*/
//____________________________________________________________________________


#include <cassert>

#include "TSystem.h"

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



#include "SeaNuDrivers/GSeaSyst.h"

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

using namespace std;
using namespace genie;

//____________________________________________________________________________
GSeaSyst::GSeaSyst()
{
  this->SetSystParamList();
  this->Configure();
  this->BuildSystInput();
}
//___________________________________________________________________________
GSeaSyst::~GSeaSyst(){

  for(unsigned iParam=0; iParam<fSystParam.size(); iParam++){
    for(unsigned iTwk=0; iTwk<fSystParam[iParam].ntwk; iTwk++){
       fSystParam[iParam].ReWeight[iTwk]=0;
       delete fSystParam[iParam].ReWeight[iTwk];
    }
  }

  for(unsigned iList=0; iList<fSystList.size(); iList++){
    fSystList[iList].ReWeight = 0;
    delete fSystList[iList].ReWeight;
  }

}
//___________________________________________________________________________
SystParam GSeaSyst::SetParam(string name, int ntwk){
  SystParam P;
  P.name=name;
  P.param=GSyst::FromString(name);
    if(P.param==0){
    LOG("GSeaSyst", pFATAL) <<"systematic parameter "<< name <<" does not exist..., exiting";
    exit(1);
  }

  double twk_dial_min  = -1.0;
  double twk_dial_max  =  1.0;
  double twk_dial_step;

  if(ntwk % 2 == 0)ntwk+=1;
  if(ntwk < 3)ntwk=3;
  twk_dial_step=(twk_dial_max-twk_dial_min)/(ntwk-1);

  double tweaking_value;

  for(int itwk=0; itwk<ntwk; itwk++){
    tweaking_value=twk_dial_min+(double)itwk*twk_dial_step;
    P.tweaking_value.push_back(tweaking_value);
  }

  P.ntwk=P.tweaking_value.size();

  return P;
}
//___________________________________________________________________________
ListParam GSeaSyst::SetParam(string name, double tweaking_value){
  ListParam P;
  P.name=name;
  P.param=GSyst::FromString(name);
    if(P.param==0){
    LOG("GSeaSyst", pFATAL) <<"systematic parameter "<< name <<" does not exist..., exiting";
    exit(1);
  }
  P.tweaking_value=tweaking_value;
  return P;
}



//___________________________________________________________________________
void GSeaSyst::SetSystParamList(){

  ListParam P;
  SystList List;

  if(gSystem->Getenv("SYSTLIST")){  // systematic parameter list from file

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

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

      this->ReadParamListFile(filename);

  }
  else { // default systematic parameter lists

    P=this->SetParam("NormCCQE",        +1.0);
    List.params.push_back(P);
    P=this->SetParam("MaCCQEshape",     +1.0);
    List.params.push_back(P);
    P=this->SetParam("NormCCRES",       -1.0);
    List.params.push_back(P);
    P=this->SetParam("VecFFCCQEshape",  -1.0);
    List.params.push_back(P);
    P=this->SetParam("MaCCRESshape",    -1.0);
    List.params.push_back(P);
    P=this->SetParam("MvCCRESshape",    +0.5);
    List.params.push_back(P);
    P=this->SetParam("NormNCRES",       +1.0);
    List.params.push_back(P);
    P=this->SetParam("MaNCRESshape",    -0.7);
    List.params.push_back(P);
    P=this->SetParam("MvNCRESshape",    +0.3);
    List.params.push_back(P);
    P=this->SetParam("NonRESBGvpCC1pi", +0.5);
    List.params.push_back(P);
    P=this->SetParam("NonRESBGvnCC1pi", +0.5);
    List.params.push_back(P);
    P=this->SetParam("MaCOHpi",         -0.5);
    List.params.push_back(P);
    P=this->SetParam("MFP_pi",          +1.0);
    List.params.push_back(P);
    P=this->SetParam("MFP_N",           -1.0);
    List.params.push_back(P);
    P=this->SetParam("FrPiProd_pi",     -0.7);
    List.params.push_back(P);
    P=this->SetParam("AGKYxF1pi",       -1.0);
    List.params.push_back(P);
    P=this->SetParam("AGKYpT1pi",       +1.0);
    List.params.push_back(P);
    P=this->SetParam("FormZone",        +1.0);
//    List.params.push_back(P);
    P=this->SetParam("Theta_Delta2Npi", +1.0);
    List.params.push_back(P);

    List.npar=List.params.size();

    fSystList.push_back(List);
    List.params.clear();

}

  return;
}

//___________________________________________________________________________
void GSeaSyst::ReadParamListFile(string filename){

  ListParam P;
  SystList List;

  SystParam Param;

  fSystList.clear();
  fSystParam.clear();


  xmlTextReaderPtr reader;
  reader = xmlNewTextReaderFilename(filename.c_str());


  string stype,sparam,stweak;
  int ntwk;
  double tweaking_value;

  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 *) "syst_params")){
      LOG("GSeaSyst", pFATAL) << "root element in filename " << filename <<" must be syst_param; exiting... ";
      exit(1);
     }
    }

    if(!xmlStrcmp(name, (const xmlChar*) "param_set" ) && type==1){
      xmlChar * xtype   = xmlTextReaderGetAttribute(reader,(const xmlChar*)"set_type");
      stype = utils::str::TrimSpaces((const char *)xtype);
      size_t pos=stype.find(";");
      if(pos==string::npos)stweak="simple";
      else {
        stweak=stype.substr(pos+1,stype.size()-pos);
        stype.erase(pos,stype.size()-pos);
      }
      xmlFree(xtype);
    }


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

    if(depth==3){
      if(stype=="parameter"){
        ntwk=atoi((const char *)value);
        Param=this->SetParam(sparam,ntwk);
        fSystParam.push_back(Param);
      }
      else if(stype=="list"){
        tweaking_value=(double)atof((const char *)value);
        P=this->SetParam(sparam,tweaking_value);
        List.params.push_back(P);
      }
      else {
        LOG("GSeaSyst", pFATAL) << "error in filename " << filename <<"; check attribute set_type; exiting... ";
        exit(1);
      }

    }

    if(!xmlStrcmp(name, (const xmlChar*) "param_set" ) && stype=="list" && type!=1){
      List.npar=List.params.size();
      fSystList.push_back(List);
      if(stweak=="mirror"){
        for(unsigned iParam=0; iParam<List.params.size(); iParam++)List.params[iParam].tweaking_value *= -1;
        fSystList.push_back(List);
      }
      List.params.clear();
    }

    xmlFree(name);
    xmlFree(value);

  }

  return;
}
//___________________________________________________________________________
void GSeaSyst::Configure(void){

  for(unsigned iParam=0; iParam<fSystParam.size(); iParam++){
    for(unsigned iTwk=0; iTwk<fSystParam[iParam].ntwk; iTwk++){

      GReWeight * ReWeight;
      fSystParam[iParam].ReWeight.push_back(ReWeight);
      fSystParam[iParam].ReWeight[iTwk] = new GReWeight;

      // add calculators
      fSystParam[iParam].ReWeight[iTwk]->AdoptWghtCalc( "xsec_ncel",       new GReWeightNuXSecNCEL );
      fSystParam[iParam].ReWeight[iTwk]->AdoptWghtCalc( "xsec_ccqe",       new GReWeightNuXSecCCQE );
      fSystParam[iParam].ReWeight[iTwk]->AdoptWghtCalc( "xsec_ccqe_vec",   new GReWeightNuXSecCCQEvec );
      fSystParam[iParam].ReWeight[iTwk]->AdoptWghtCalc( "xsec_ccres",      new GReWeightNuXSecCCRES );
      fSystParam[iParam].ReWeight[iTwk]->AdoptWghtCalc( "xsec_ncres",      new GReWeightNuXSecNCRES );
      fSystParam[iParam].ReWeight[iTwk]->AdoptWghtCalc( "xsec_nonresbkg",  new GReWeightNonResonanceBkg );
      fSystParam[iParam].ReWeight[iTwk]->AdoptWghtCalc( "xsec_dis",        new GReWeightNuXSecDIS );
      fSystParam[iParam].ReWeight[iTwk]->AdoptWghtCalc( "xsec_coh",        new GReWeightNuXSecCOH );
      fSystParam[iParam].ReWeight[iTwk]->AdoptWghtCalc( "nuclear_qe",      new GReWeightFGM );
      fSystParam[iParam].ReWeight[iTwk]->AdoptWghtCalc( "nuclear_dis",     new GReWeightDISNuclMod );
      fSystParam[iParam].ReWeight[iTwk]->AdoptWghtCalc( "hadro_res_decay", new GReWeightResonanceDecay );
      fSystParam[iParam].ReWeight[iTwk]->AdoptWghtCalc( "hadro_fzone",     new GReWeightFZone );
      fSystParam[iParam].ReWeight[iTwk]->AdoptWghtCalc( "hadro_intranuke", new GReWeightINuke );
      fSystParam[iParam].ReWeight[iTwk]->AdoptWghtCalc( "hadro_agky",      new GReWeightAGKY );

      GSystSet & syst = fSystParam[iParam].ReWeight[iTwk]->Systematics();

      // Fine-tune weight calculators

      if(fSystParam[iParam].param == kXSecTwkDial_MaCCQE) {
        // By default GReWeightNuXSecCCQE is in `NormAndMaShape' mode
        // where Ma affects the shape of dsigma/dQ2 and a different param affects the normalization
        // If the input is MaCCQE, switch the weight calculator to `Ma' mode
        GReWeightNuXSecCCQE * rwccqe = dynamic_cast<GReWeightNuXSecCCQE *> (fSystParam[iParam].ReWeight[iTwk]->WghtCalc("xsec_ccqe"));
        rwccqe->SetMode(GReWeightNuXSecCCQE::kModeMa);
      }
      if(fSystParam[iParam].param == kXSecTwkDial_MaCCRES || fSystParam[iParam].param == kXSecTwkDial_MvCCRES) {
        // As above, but for the GReWeightNuXSecCCRES weight calculator
        GReWeightNuXSecCCRES * rwccres = dynamic_cast<GReWeightNuXSecCCRES *> (fSystParam[iParam].ReWeight[iTwk]->WghtCalc("xsec_ccres"));
        rwccres->SetMode(GReWeightNuXSecCCRES::kModeMaMv);
      }
      if(fSystParam[iParam].param == kXSecTwkDial_MaNCRES || fSystParam[iParam].param == kXSecTwkDial_MvNCRES) {
        // As above, but for the GReWeightNuXSecNCRES weight calculator
        GReWeightNuXSecNCRES * rwncres = dynamic_cast<GReWeightNuXSecNCRES *> (fSystParam[iParam].ReWeight[iTwk]->WghtCalc("xsec_ncres"));
        rwncres->SetMode(GReWeightNuXSecNCRES::kModeMaMv);
      }
      if(fSystParam[iParam].param == kXSecTwkDial_AhtBYshape  || fSystParam[iParam].param == kXSecTwkDial_BhtBYshape ||
        fSystParam[iParam].param == kXSecTwkDial_CV1uBYshape || fSystParam[iParam].param == kXSecTwkDial_CV2uBYshape ) {
        // Similarly for the GReWeightNuXSecDIS weight calculator.
        // There the default behaviour is for the Aht, Bht, CV1u and CV2u Bodek-Yang
        // params to affects both normalization and dsigma/dxdy shape.
        // Switch mode if a shape-only param is specified.
        GReWeightNuXSecDIS * rwdis = dynamic_cast<GReWeightNuXSecDIS *> (fSystParam[iParam].ReWeight[iTwk]->WghtCalc("xsec_dis"));
        rwdis->SetMode(GReWeightNuXSecDIS::kModeABCV12uShape);
      }

      syst.Set(fSystParam[iParam].param,fSystParam[iParam].tweaking_value[iTwk]);
      fSystParam[iParam].ReWeight[iTwk]->Reconfigure();

    }
  }


  for(unsigned iList=0; iList<fSystList.size(); iList++){

    fSystList[iList].ReWeight = new GReWeight;

      // add calculators
    fSystList[iList].ReWeight->AdoptWghtCalc( "xsec_ncel",       new GReWeightNuXSecNCEL );
    fSystList[iList].ReWeight->AdoptWghtCalc( "xsec_ccqe",       new GReWeightNuXSecCCQE );
    fSystList[iList].ReWeight->AdoptWghtCalc( "xsec_ccqe_vec",   new GReWeightNuXSecCCQEvec );
    fSystList[iList].ReWeight->AdoptWghtCalc( "xsec_ccres",      new GReWeightNuXSecCCRES );
    fSystList[iList].ReWeight->AdoptWghtCalc( "xsec_ncres",      new GReWeightNuXSecNCRES );
    fSystList[iList].ReWeight->AdoptWghtCalc( "xsec_nonresbkg",  new GReWeightNonResonanceBkg );
    fSystList[iList].ReWeight->AdoptWghtCalc( "xsec_dis",        new GReWeightNuXSecDIS );
    fSystList[iList].ReWeight->AdoptWghtCalc( "xsec_coh",        new GReWeightNuXSecCOH );
    fSystList[iList].ReWeight->AdoptWghtCalc( "nuclear_qe",      new GReWeightFGM );
    fSystList[iList].ReWeight->AdoptWghtCalc( "nuclear_dis",     new GReWeightDISNuclMod );
    fSystList[iList].ReWeight->AdoptWghtCalc( "hadro_res_decay", new GReWeightResonanceDecay );
    fSystList[iList].ReWeight->AdoptWghtCalc( "hadro_fzone",     new GReWeightFZone );
    fSystList[iList].ReWeight->AdoptWghtCalc( "hadro_intranuke", new GReWeightINuke );
    fSystList[iList].ReWeight->AdoptWghtCalc( "hadro_agky",      new GReWeightAGKY );

    GSystSet & syst = fSystList[iList].ReWeight->Systematics();

    for(unsigned iParam=0; iParam<fSystList[iList].params.size(); iParam++){

      // Fine-tune weight calculators
      if(fSystList[iList].params[iParam].param == kXSecTwkDial_MaCCQE) {
        // By default GReWeightNuXSecCCQE is in `NormAndMaShape' mode
        // where Ma affects the shape of dsigma/dQ2 and a different param affects the normalization
        // If the input is MaCCQE, switch the weight calculator to `Ma' mode
        GReWeightNuXSecCCQE * rwccqe = dynamic_cast<GReWeightNuXSecCCQE *> (fSystList[iList].ReWeight->WghtCalc("xsec_ccqe"));
        rwccqe->SetMode(GReWeightNuXSecCCQE::kModeMa);
      }
      if(fSystList[iList].params[iParam].param == kXSecTwkDial_MaCCRES || fSystList[iList].params[iParam].param == kXSecTwkDial_MvCCRES) {
        // As above, but for the GReWeightNuXSecCCRES weight calculator
        GReWeightNuXSecCCRES * rwccres = dynamic_cast<GReWeightNuXSecCCRES *> (fSystList[iList].ReWeight->WghtCalc("xsec_ccres"));
        rwccres->SetMode(GReWeightNuXSecCCRES::kModeMaMv);
      }
      if(fSystList[iList].params[iParam].param == kXSecTwkDial_MaNCRES || fSystList[iList].params[iParam].param == kXSecTwkDial_MvNCRES) {
        // As above, but for the GReWeightNuXSecNCRES weight calculator
        GReWeightNuXSecNCRES * rwncres = dynamic_cast<GReWeightNuXSecNCRES *> (fSystList[iList].ReWeight->WghtCalc("xsec_ncres"));
        rwncres->SetMode(GReWeightNuXSecNCRES::kModeMaMv);
      }
      if(fSystList[iList].params[iParam].param == kXSecTwkDial_AhtBYshape  || fSystList[iList].params[iParam].param == kXSecTwkDial_BhtBYshape ||
        fSystList[iList].params[iParam].param == kXSecTwkDial_CV1uBYshape || fSystList[iList].params[iParam].param == kXSecTwkDial_CV2uBYshape ) {
        // Similarly for the GReWeightNuXSecDIS weight calculator.
        // There the default behaviour is for the Aht, Bht, CV1u and CV2u Bodek-Yang
        // params to affects both normalization and dsigma/dxdy shape.
        // Switch mode if a shape-only param is specified.
        GReWeightNuXSecDIS * rwdis = dynamic_cast<GReWeightNuXSecDIS *> (fSystList[iList].ReWeight->WghtCalc("xsec_dis"));
        rwdis->SetMode(GReWeightNuXSecDIS::kModeABCV12uShape);
      }

      syst.Set(fSystList[iList].params[iParam].param,fSystList[iList].params[iParam].tweaking_value);
    }

    fSystList[iList].ReWeight->Reconfigure();
  }

  return;
}
//_________________________________________________________________________________
SystWeights GSeaSyst::ComputeWeights(EventRecord * event){

  double wg=1.;
  EventRecord & evt = * event;
  SystWeights Weights;

  double weight=evt.Weight();

  evt.SetWeight(1);

  for(unsigned iParam=0; iParam<fSystParam.size(); iParam++){
    StrParamWeight pWght;
    pWght.name=fSystParam[iParam].name;
    for(unsigned iTwk=0; iTwk<fSystParam[iParam].ntwk; iTwk++){
      wg=fSystParam[iParam].ReWeight[iTwk]->CalcWeight(evt);
      pWght.wght.push_back(wg);
    }
    Weights.WParam.push_back(pWght);
  }

  for(unsigned iList=0; iList<fSystList.size(); iList++){
    wg=fSystList[iList].ReWeight->CalcWeight(evt);
    Weights.WList.push_back(wg);
  }

  evt.SetWeight(weight);

  return Weights;
}

//_________________________________________________________________________________
void GSeaSyst::BuildSystInput(){

  InputParam param;
  vector<InputParam> vecparam;

  for(unsigned i=0; i<fSystParam.size(); i++){
    param.name=fSystParam[i].name;
    param.ntwk=fSystParam[i].ntwk;
    param.tweaking_value=fSystParam[i].tweaking_value;
    fSystInput.SingleParams.push_back(param);
  }

  for(unsigned i=0; i<fSystList.size(); i++){

    vecparam.clear();
    for(unsigned j=0; j<fSystList[i].params.size(); j++){
      param.name=fSystList[i].params[j].name;
      param.ntwk=1;
      param.tweaking_value.clear();
      param.tweaking_value.push_back(fSystList[i].params[j].tweaking_value);
      vecparam.push_back(param);
    }

    fSystInput.ListParams.push_back(vecparam);

  }

  return;
}