SystsGENIEAna.C

Go to the documentation of this file.

#include "NuXAna/Prediction/PredictionSterile.h"
#include "CAFAna/Systs/Systs.h"
#include "NuXAna/macros/NuSPlotFunctions.h"
#include "NuXAna/macros/NuSSystFunctions.h"

using namespace ana;

void SystsGENIEAna()
{
  TH1::AddDirectory(0);

  // A map of systematic labels
  std::map<std::string, rwgt::ReweightLabel_t> shiftlabels; // Easier to get label information like this...
  // NCEL parameters:
  shiftlabels["ReweightMaNCEL"]              = rwgt::fReweightMaNCEL;
  shiftlabels["ReweightEtaNCEL"]             = rwgt::fReweightEtaNCEL;
  //shiftlabels["ReweightMaCCQE"]              = rwgt::fReweightMaCCQE;
  shiftlabels["ReweightMaCCQEshape"]         = rwgt::fReweightMaCCQEshape;
  shiftlabels["ReweightNormCCQE"]            = rwgt::fReweightNormCCQE;
  //shiftlabels["ReweightNormCCQEenu"]         = rwgt::fReweightNormCCQEenu;
  shiftlabels["ReweightCCQEPauliSupViaKF"]   = rwgt::fReweightCCQEPauliSupViaKF;
  shiftlabels["ReweightVecCCQEshape"]        = rwgt::fReweightVecCCQEshape;
  shiftlabels["ReweightCCQEMomDistroFGtoSF"] = rwgt::fReweightCCQEMomDistroFGtoSF;
  //shiftlabels["ReweightNormCCRES"]           = rwgt::fReweightNormCCRES;
  //shiftlabels["ReweightNormNCRES"]           = rwgt::fReweightNormNCRES;
  //shiftlabels["ReweightMaCCRESshape"]        = rwgt::fReweightMaCCRESshape;
  //shiftlabels["ReweightMvCCRESshape"]        = rwgt::fReweightMvCCRESshape;
  //shiftlabels["ReweightMaNCRESshape"]        = rwgt::fReweightMaNCRESshape;
  //shiftlabels["ReweightMvNCRESshape"]        = rwgt::fReweightMvNCRESshape;
  shiftlabels["ReweightMaCCRES"]             = rwgt::fReweightMaCCRES;
  shiftlabels["ReweightMvCCRES"]             = rwgt::fReweightMvCCRES;
  shiftlabels["ReweightMaNCRES"]             = rwgt::fReweightMaNCRES;
  shiftlabels["ReweightMvNCRES"]             = rwgt::fReweightMvNCRES;
  shiftlabels["ReweightMaCOHpi"]             = rwgt::fReweightMaCOHpi;
  shiftlabels["ReweightR0COHpi"]             = rwgt::fReweightR0COHpi;
  shiftlabels["ReweightRvpCC1pi"]            = rwgt::fReweightRvpCC1pi;
  shiftlabels["ReweightRvpCC2pi"]            = rwgt::fReweightRvpCC2pi;
  shiftlabels["ReweightRvnCC1pi"]            = rwgt::fReweightRvnCC1pi;
  shiftlabels["ReweightRvnCC2pi"]            = rwgt::fReweightRvnCC2pi;
  shiftlabels["ReweightRvpNC1pi"]            = rwgt::fReweightRvpNC1pi;
  shiftlabels["ReweightRvpNC2pi"]            = rwgt::fReweightRvpNC2pi;
  shiftlabels["ReweightRvnNC1pi"]            = rwgt::fReweightRvnNC1pi;
  shiftlabels["ReweightRvnNC2pi"]            = rwgt::fReweightRvnNC2pi;
  shiftlabels["ReweightRvbarpCC1pi"]         = rwgt::fReweightRvbarpCC1pi;
  shiftlabels["ReweightRvbarpCC2pi"]         = rwgt::fReweightRvbarpCC2pi;
  shiftlabels["ReweightRvbarnCC1pi"]         = rwgt::fReweightRvbarnCC1pi;
  shiftlabels["ReweightRvbarnCC2pi"]         = rwgt::fReweightRvbarnCC2pi;
  shiftlabels["ReweightRvbarpNC1pi"]         = rwgt::fReweightRvbarpNC1pi;
  shiftlabels["ReweightRvbarpNC2pi"]         = rwgt::fReweightRvbarpNC2pi;
  shiftlabels["ReweightRvbarnNC1pi"]         = rwgt::fReweightRvbarnNC1pi;
  shiftlabels["ReweightRvbarnNC2pi"]         = rwgt::fReweightRvbarnNC2pi;
  shiftlabels["ReweightAhtBY"]               = rwgt::fReweightAhtBY;
  shiftlabels["ReweightBhtBY"]               = rwgt::fReweightBhtBY;
  shiftlabels["ReweightCV1uBY"]              = rwgt::fReweightCV1uBY;
  shiftlabels["ReweightCV2uBY"]              = rwgt::fReweightCV2uBY;
  //shiftlabels["ReweightAhtBYshape"]          = rwgt::fReweightAhtBYshape;
  //shiftlabels["ReweightBhtBYshape"]          = rwgt::fReweightBhtBYshape;
  //shiftlabels["ReweightCV1uBYshape"]         = rwgt::fReweightCV1uBYshape;
  //shiftlabels["ReweightCV2uBYshape"]         = rwgt::fReweightCV2uBYshape;
  shiftlabels["ReweightNormDISCC"]           = rwgt::fReweightNormDISCC;
  shiftlabels["ReweightRnubarnuCC"]          = rwgt::fReweightRnubarnuCC;
  shiftlabels["ReweightDISNuclMod"]          = rwgt::fReweightDISNuclMod;
  shiftlabels["ReweightAGKY_xF1pi"]          = rwgt::fReweightAGKY_xF1pi;
  shiftlabels["ReweightAGKY_pT1pi"]          = rwgt::fReweightAGKY_pT1pi;
  shiftlabels["ReweightFormZone"]            = rwgt::fReweightFormZone;
  shiftlabels["ReweightTheta_Delta2Npi"]     = rwgt::fReweightTheta_Delta2Npi;
  shiftlabels["ReweightBR1gamma"]            = rwgt::fReweightBR1gamma;
  shiftlabels["ReweightBR1eta"]              = rwgt::fReweightBR1eta;
  shiftlabels["ReweightMFP_N"]               = rwgt::fReweightMFP_N;
  shiftlabels["ReweightFrCEx_N"]             = rwgt::fReweightFrCEx_N;
  shiftlabels["ReweightFrElas_N"]            = rwgt::fReweightFrElas_N;
  shiftlabels["ReweightFrInel_N"]            = rwgt::fReweightFrInel_N;
  shiftlabels["ReweightFrAbs_N"]             = rwgt::fReweightFrAbs_N;
  shiftlabels["ReweightFrPiProd_N"]          = rwgt::fReweightFrPiProd_N;
  shiftlabels["ReweightMFP_pi"]              = rwgt::fReweightMFP_pi;
  shiftlabels["ReweightFrCEx_pi"]            = rwgt::fReweightFrCEx_pi;
  shiftlabels["ReweightFrElas_pi"]           = rwgt::fReweightFrElas_pi;
  shiftlabels["ReweightFrInel_pi"]           = rwgt::fReweightFrInel_pi;
  shiftlabels["ReweightFrAbs_pi"]            = rwgt::fReweightFrAbs_pi;
  shiftlabels["ReweightFrPiProd_pi"]         = rwgt::fReweightFrPiProd_pi;
  //shiftlabels["ReweightNC"]                  = rwgt::fReweightNC;
  // These two shifts are GENIE systs, but not in the form of reweights as above
  // Unfortunately they need a rwgt label, so give them the null one.
  shiftlabels["mecScale"] = rwgt::kReweightNull;
  shiftlabels["RPA"]      = rwgt::kReweightNull;

  std::string labelRecoE = "Calorimetric Energy (GeV)";
  std::vector<std::string> cut_samples;
  cut_samples.push_back("Ana01");

  // Set up maps to the decompositions/predictions (each flavor component)
  // Nominal maps are indexed purely by sample label
  // Shifted maps are indexed by sample label, systematic label, then sigma of the shift
  std::map<std::string,                                     IDecomp*>     decomps_nominal;
  std::map<std::string, std::map<std::string, std::map<int, IDecomp*> > > decomps_shifted;
  std::map<std::string,                                     PredictionNoExtrap*>     predsNE_nominal;
  std::map<std::string, std::map<std::string, std::map<int, PredictionNoExtrap*> > > predsNE_shifted;
  std::map<std::string,                                     PredictionSterile*>      predsSt_nominal;
  std::map<std::string, std::map<std::string, std::map<int, PredictionSterile*> > >  predsSt_shifted;

  std::string folder = "/nova/ana/steriles/Ana01/Systematics/";
  std::string filenm = "SystsGENIE";

  std::string loadLocation = folder + filenm + ".root";
  std::string saveLocation = folder + filenm + "Ana.root";
  std::string textLocation = folder + filenm + "Ana.txt";

  TFile* rootL = new TFile(loadLocation.c_str(), "READ");
  LoadMaps(rootL, &decomps_nominal, &decomps_shifted);
  LoadMaps(rootL, &predsNE_nominal, &predsNE_shifted);
  LoadMaps(rootL, &predsSt_nominal, &predsSt_shifted);
  rootL->Close();

  TFile* rootF = new TFile(saveLocation.c_str(), "RECREATE");
  FILE* textF;
  textF = fopen(textLocation.c_str(), "w");

  for(const auto& sample : cut_samples) {
    strings strs;
    strs.fComponent = "NC";
    strs.fDet = "ND";
    strs.fPOT = StringFromDouble(NCSCALE);
    strs.fSample = sample;
    strs.fSystType = "GENIE";
    strs.fXLabel = labelRecoE;

    InitializeSystText(textF, strs);

    TH1* hNomNC_decomp = GetNC(decomps_nominal[sample], NCSCALE);
    TH1* hNomBG_decomp = GetBG(decomps_nominal[sample], NCSCALE);

    TH1* hNomNC_predNE = GetNC(predsNE_nominal[sample], NCSCALE);
    TH1* hNomBG_predNE = GetBG(predsNE_nominal[sample], NCSCALE);

    TH1* hNomNC_predSt = GetNC(predsSt_nominal[sample], NCSCALE);
    TH1* hNomBG_predSt = GetBG(predsSt_nominal[sample], NCSCALE);
    
    for(const auto& shifts : shiftlabels) {
      if(shifts.first.compare("RPA") != 0 &&
         shifts.first.compare("mecScale") != 0) {
        ISyst* temp = new GenieRwgtTableSyst(shifts.second);
        strs.fSystS = temp->ShortName();
        strs.fSystL = temp->LatexName();
      }
      else {
        if(shifts.first.compare("RPA") == 0) {
          strs.fSystS = "RPA";
          strs.fSystL = "RPA Weight";
        }
        else {
          strs.fSystS = "mecScale";
          strs.fSystL = "MEC Scale";
        }
      }

      strs.fComponent = "NC";
      strs.fDet = "ND";
      PlotSyst(hNomNC_decomp,
               decomps_shifted[sample][shifts.first],
               rootF, textF, strs, NCSCALE, true, false);

      strs.fComponent = "BG";
      PlotSyst(hNomBG_decomp,
               decomps_shifted[sample][shifts.first],
               rootF, textF, strs, NCSCALE, false, false);

      strs.fComponent = "NC";
      strs.fDet = "FD";
      PlotSyst(hNomNC_predNE,
               predsNE_shifted[sample][shifts.first],
               rootF, textF, strs, NCSCALE, true, false);

      strs.fComponent = "BG";
      PlotSyst(hNomBG_predNE,
               predsNE_shifted[sample][shifts.first],
               rootF, textF, strs, NCSCALE, false, false);

      strs.fComponent = "NC";
      strs.fDet = "EX";
      PlotSyst(hNomNC_predSt,
               predsSt_shifted[sample][shifts.first],
               rootF, textF, strs, NCSCALE, true, false);

      strs.fComponent = "BG";
      PlotSyst(hNomBG_predSt,
               predsSt_shifted[sample][shifts.first],
               rootF, textF, strs, NCSCALE, false, false);
    }

    strs.fComponent = "NC";
    strs.fDet = "ND";
    PlotSystBand(hNomNC_decomp, rootF, textF, strs);

    strs.fComponent = "BG";
    PlotSystBand(hNomBG_decomp, rootF, textF, strs);

    strs.fComponent = "NC";
    strs.fDet = "FD";
    PlotSystBand(hNomNC_predNE, rootF, textF, strs);

    strs.fComponent = "BG";
    PlotSystBand(hNomBG_predNE, rootF, textF, strs);

    strs.fComponent = "NC";
    strs.fDet = "EX";
    PlotSystBand(hNomNC_predSt, rootF, textF, strs);

    strs.fComponent = "BG";
    PlotSystBand(hNomBG_predSt, rootF, textF, strs);
  }

  fclose(textF);
  rootF->Close();

  std::cout << "Finished. Ctrl+Z should be safe if the macro didn't quit." << std::endl;
}