specprod_systematics.C

Go to the documentation of this file.

//////////////////////////////////////////////////////////////////
/// \brief Macro to create a given xsec systematic for the 
///        numucc inc analysis
/// \author Connor Johnson
/// \date June 2019
//////////////////////////////////////////////////////////////////

#include "CAFAna/Vars/Vars.h"
#include "CAFAna/Core/Binning.h"
#include "CAFAna/Cuts/Cuts.h"
#include "CAFAna/Cuts/SpillCuts.h"
#include "CAFAna/Core/Spectrum.h"
#include "CAFAna/Core/SpectrumLoader.h"
#include "CAFAna/Core/SpectrumLoaderBase.h"
#include "3FlavorAna/Vars/NumuVars.h"
#include "CAFAna/Systs/Systs.h"
#include "CAFAna/Systs/XSecSystLists.h"
#include "CAFAna/Systs/BeamSysts.h"
#include "CAFAna/Systs/EnergySysts2018.h"
#include "CAFAna/XSec/GenieMultiverseSyst.h"

#include "CAFAna/Vars/GenieWeights.h"
#include "CAFAna/Vars/PPFXWeights.h"

#include "NDAna/numucc_inc/NumuCCIncAnalysis.h"
#include "NDAna/numucc_inc/NumuCCIncBins.h"
#include "NDAna/numucc_inc/NumuCCIncVars.h"
#include "NDAna/numucc_inc/NumuCCIncCuts.h"
#include "NDAna/numucc_inc/NDXSecMuonPID.h"
#include "NDAna/numucc_inc/UnfoldingVariable.h"

#include "TFile.h"
#include "TH1.h"
#include "TH2.h"

#include <iostream>
#include <cmath>
#include <sstream>
#include <string>
#include <fstream>
#include <iomanip>
#include <map>
#include <algorithm>

//////////////////////////////////////////////////////////////////
/// List of systematics to allow studies over. Each section treated differently
//////////////////////////////////////////////////////////////////

/// Systematics defined by a separate hist axis
std::map<string, const vector<const ana::HistAxis* > > histaxis_systematics = {
  {"muoneup",                  {&ana::kRecoMuKEVsCosVsEavailStandardAxisUp,      &ana::kRecoEStandardAxis_MuonEUp, &ana::kRecoQ2StandardAxis_MuonEUp}},
  {"muonedw",                  {&ana::kRecoMuKEVsCosVsEavailStandardAxisDw,      &ana::kRecoEStandardAxis_MuonEDw, &ana::kRecoQ2StandardAxis_MuonEDw}},
  {"angleup",                  {&ana::kRecoMuKEVsCosVsEavailStandardAxisAngleUp, &ana::kRecoEStandardAxis        , &ana::kRecoQ2StandardAxis_AngleUp}},
  {"angledw",                  {&ana::kRecoMuKEVsCosVsEavailStandardAxisAngleDw, &ana::kRecoEStandardAxis        , &ana::kRecoQ2StandardAxis_AngleDw}}
};

/// Systematics defined by a BeamSyst
std::map<string, const ana::BeamSyst*> beamSystematics = {
  {"2kA",                      &ana::kBeamHornCurrent},
  {"02mmBeamSpotSize",         &ana::kBeamSpotSize},
  {"1mmBeamShiftX",            &ana::kBeamPosX},
  {"1mmBeamShiftY",            &ana::kBeamPosY},
  {"3mmHorn1X",                &ana::kBeamH1PosX},
  {"3mmHorn1Y",                &ana::kBeamH1PosY},
  {"3mmHorn2X",                &ana::kBeamH2PosX},
  {"3mmHorn2Y",                &ana::kBeamH2PosY},
  {"7mmTargetZ",               &ana::kBeamTarget},
  {"MagneticFieldinDecayPipe", &ana::kBeamMagField},
  {"1mmHornWater",             &ana::kBeamGeomWater}
};

/// Other pre-defined systematics
std::map<string, const ana::ISyst*> otherSystematics = {
  {"neutron",              &ana::kNeutronVisEScalePrimariesSyst2018}
};

/// Systematics that run on the nominal dataset
std::vector<std::string> standardVars = 
{
  "cv",
  "nominal",
  "lightdown",
  "lightup",
  "ckv",
  "calibneg",
  "calibpos",
  "calibshift",
  "neutron_nom",
  "nom_good_seventh"
};

/// Systematics that use a different dataset
std::map<string, const std::string> dataset_to_samdef = {
  {"nominal",          "nd_fhc_remid-hotfix_caf_minus_muonid_training_minus_fakedata"},
  {"nom_good_seventh", "nd_fhc_remid-hotfix_goodseventh_common_subruns"},
  {"lightdown",        "prod_caf_R17-11-14-prod4reco.remid-hotfix.b_nd_genie_nonswap_fhc_nova_v08_full_lightmodel-lightdown-calibup_v1_good_seventh"},
  {"lightup",          "prod_caf_R17-11-14-prod4reco.remid-hotfix.b_nd_genie_nonswap_fhc_nova_v08_full_lightmodel-lightup-calibdown_v1_good_seventh"},
  {"ckv",              "prod_caf_R17-11-14-prod4reco.remid-hotfix.b_nd_genie_nonswap_fhc_nova_v08_full_ckv-proton-shift-down_v1"},
  {"calibneg",         "prod_caf_R17-11-14-prod4reco.remid-hotfix.b_nd_genie_nonswap_fhc_nova_v08_full_calib-shift-nd-xyview-neg-offset_v1_good_seventh"},
  {"calibpos",         "prod_caf_R17-11-14-prod4reco.remid-hotfix.b_nd_genie_nonswap_fhc_nova_v08_full_calib-shift-nd-xyview-pos-offset_v1_good_seventh"},
  {"calibshift",       "prod_caf_R17-11-14-prod4reco.remid-hotfix.b_nd_genie_nonswap_fhc_nova_v08_full_calib-shift-nd-func_v1_good_seventh"},
  {"neutron",          "prod_caf_R17-11-14-prod4reco.neutron-respin.c_nd_genie_nonswap_fhc_nova_v08_full_v1"},
  {"neutron_nom",      "prod_caf_R17-11-14-prod4reco.neutron-respin.c_nd_genie_nonswap_fhc_nova_v08_full_v1"}
};

//////////////////////////////////////////////////////////////////
/// Override folder naming in certain cases
//////////////////////////////////////////////////////////////////
std::map<string, const std::string> syst_to_folder = {
  {"muoneup",                 "MuESUp"},
  {"muonedw",                 "MuESDw"}
};

//////////////////////////////////////////////////////////////////
/// Define the weights. Build off of each other, and default to "1".
//////////////////////////////////////////////////////////////////

const ana::NuTruthVar wgt_xsST = ana::kXSecCVWgt2018_smallerDISScale_NT;
const ana::Var wgt_xs          = ana::VarFromNuTruthVar(wgt_xsST, 1);
const ana::NuTruthVar wgtST    = ana::kPPFXFluxCVWgtST*wgt_xsST;
const ana::Var wgt             = ana::VarFromNuTruthVar(wgtST, 1);

//////////////////////////////////////////////////////////////////
/// Process all \a systs input, and store output in \a outfile
//////////////////////////////////////////////////////////////////
void specprod_systematics(const std::vector<std::string> systs, const std::string outfile = "fSystSpec.root"){
  std::map<std::string, ana::SpectrumLoaderBase*> vloaders; ///< Map of loaders
  std::map<std::string, std::map<std::string, ana::xsec::UnfoldingVariable*>>   vUnfoldingVars;  ///< Nested map of UnfoldingVariables ([systName][variableName])

  // Loop through systematics
  for (unsigned int isyst = 0; isyst < systs.size(); isyst++){
    std::string syst = systs[isyst];
    std::string dataset;
    if (dataset_to_samdef.find(syst) != dataset_to_samdef.end())
      dataset = syst;
    else
      dataset = "nominal";

    cout << "Creating spectrums for spec: " << syst << " on dataset: " << dataset << endl;

    /// Create loader, only one per dataset
    ana::SpectrumLoaderBase * loader;
    if (vloaders.find(dataset) != vloaders.end())
      loader = vloaders[dataset];
    else{
      loader = new ana::SpectrumLoader(dataset_to_samdef[dataset]);
      loader->SetSpillCut(ana::kStandardSpillCuts);
      vloaders[dataset] = loader;
    }

    // // CV, and the systematic datasets
    if (std::find(standardVars.begin(), standardVars.end(), syst) != standardVars.end()){
      vUnfoldingVars[syst]["EAvail"] = new ana::xsec::UnfoldingVariable(loader, 
        &ana::kTrueMuKEVsCosVsEavailStandardAxisST, &ana::kRecoMuKEVsCosVsEavailStandardAxis,
        &ana::kAllNumuCCCuts, &ana::kIsTrueSigST, &ana::kNoShift, &wgtST);

      vUnfoldingVars[syst]["ENu"]    = new ana::xsec::UnfoldingVariable(loader, 
        &ana::kTrueEStandardAxisST, &ana::kRecoEStandardAxis,
        &ana::kAllNumuCC1DCuts, &ana::kIsTrueSig1DST, &ana::kNoShift, &wgtST);

      vUnfoldingVars[syst]["Q2"]     = new ana::xsec::UnfoldingVariable(loader,
        &ana::kTrueQ2StandardAxisST, &ana::kRecoQ2StandardAxis,
        &ana::kAllNumuCC1DCuts, &ana::kIsTrueSig1DST, &ana::kNoShift, &wgtST);
    }

    // HistAxis systematics (muon energy, angle systematics)
    else if (histaxis_systematics.find(syst) != histaxis_systematics.end()){
      // const ana::HistAxis * eavail_histaxis = histaxis_systematics[syst][0];
      const ana::HistAxis * enu_histaxis    = histaxis_systematics[syst][1];
      const ana::HistAxis * q2_histaxis     = histaxis_systematics[syst][2];
      vUnfoldingVars[syst]["EAvail"] = new ana::xsec::UnfoldingVariable(loader,
        &ana::kTrueMuKEVsCosVsEavailStandardAxisST, eavail_histaxis,
        &ana::kAllNumuCCCuts, &ana::kIsTrueSigST, &ana::kNoShift, &wgtST);

      vUnfoldingVars[syst]["ENu"]    = new ana::xsec::UnfoldingVariable(loader,
        &ana::kTrueEStandardAxisST, enu_histaxis,
        &ana::kAllNumuCC1DCuts, &ana::kIsTrueSig1DST, &ana::kNoShift, &wgtST);

      vUnfoldingVars[syst]["Q2"]     = new ana::xsec::UnfoldingVariable(loader,
        &ana::kTrueQ2StandardAxisST, q2_histaxis,
        &ana::kAllNumuCC1DCuts, &ana::kIsTrueSig1DST, &ana::kNoShift, &wgtST);
    }

    // Beam Systematics
    else if (beamSystematics.find(syst) != beamSystematics.end()){
      std::string bs_up = syst + "_Up";
      std::string bs_dw = syst + "_Dw";

      // First the up-shifts
      const ana::BeamSyst* bshift = beamSystematics[syst];
      ana::SystShifts * bsystshiftup = new ana::SystShifts(bshift, +1);
      vUnfoldingVars[bs_up]["EAvail"] = new ana::xsec::UnfoldingVariable(loader, 
        &ana::kTrueMuKEVsCosVsEavailStandardAxisST, &ana::kRecoMuKEVsCosVsEavailStandardAxis,
        &ana::kAllNumuCCCuts, &ana::kIsTrueSigST, bsystshiftup, &wgtST);
      vUnfoldingVars[bs_up]["ENu"]    = new ana::xsec::UnfoldingVariable(loader, 
        &ana::kTrueEStandardAxisST, &ana::kRecoEStandardAxis,
        &ana::kAllNumuCC1DCuts, &ana::kIsTrueSig1DST, bsystshiftup, &wgtST);
      vUnfoldingVars[bs_up]["Q2"]     = new ana::xsec::UnfoldingVariable(loader,
        &ana::kTrueQ2StandardAxisST, &ana::kRecoQ2StandardAxis,
        &ana::kAllNumuCC1DCuts, &ana::kIsTrueSig1DST, bsystshiftup, &wgtST);

      // Next the down-shifts
      ana::SystShifts * bsystshiftdw = new ana::SystShifts(bshift, -1);
      vUnfoldingVars[bs_dw]["EAvail"] = new ana::xsec::UnfoldingVariable(loader, 
        &ana::kTrueMuKEVsCosVsEavailStandardAxisST, &ana::kRecoMuKEVsCosVsEavailStandardAxis,
        &ana::kAllNumuCCCuts, &ana::kIsTrueSigST, bsystshiftdw, &wgtST);
      vUnfoldingVars[bs_dw]["ENu"]    = new ana::xsec::UnfoldingVariable(loader, 
        &ana::kTrueEStandardAxisST, &ana::kRecoEStandardAxis,
        &ana::kAllNumuCC1DCuts, &ana::kIsTrueSig1DST, bsystshiftdw, &wgtST);
      vUnfoldingVars[bs_dw]["Q2"]     = new ana::xsec::UnfoldingVariable(loader,
        &ana::kTrueQ2StandardAxisST, &ana::kRecoQ2StandardAxis,
        &ana::kAllNumuCC1DCuts, &ana::kIsTrueSig1DST, bsystshiftdw, &wgtST);
    }
    else if (otherSystematics.find(syst) != otherSystematics.end()){
      std::string syst_up = syst + "_Up";
      std::string syst_dw = syst + "_Dw";

      // First the up-shifts
      const ana::ISyst* isyst = otherSystematics[syst];
      ana::SystShifts * systshiftup = new ana::SystShifts(isyst, +1);
      vUnfoldingVars[syst_up]["EAvail"] = new ana::xsec::UnfoldingVariable(loader, 
        &ana::kTrueMuKEVsCosVsEavailStandardAxisST, &ana::kRecoMuKEVsCosVsEavailStandardAxis,
        &ana::kAllNumuCCCuts, &ana::kIsTrueSigST, systshiftup, &wgtST);
      vUnfoldingVars[syst_up]["ENu"]    = new ana::xsec::UnfoldingVariable(loader, 
        &ana::kTrueEStandardAxisST, &ana::kRecoEStandardAxis,
        &ana::kAllNumuCC1DCuts, &ana::kIsTrueSig1DST, systshiftup, &wgtST);
      vUnfoldingVars[syst_up]["Q2"]     = new ana::xsec::UnfoldingVariable(loader,
        &ana::kTrueQ2StandardAxisST, &ana::kRecoQ2StandardAxis,
        &ana::kAllNumuCC1DCuts, &ana::kIsTrueSig1DST, systshiftup, &wgtST);

      // Next the down-shifts
      ana::SystShifts * systshiftdw = new ana::SystShifts(isyst, -1);
      vUnfoldingVars[syst_dw]["EAvail"] = new ana::xsec::UnfoldingVariable(loader, 
        &ana::kTrueMuKEVsCosVsEavailStandardAxisST, &ana::kRecoMuKEVsCosVsEavailStandardAxis,
        &ana::kAllNumuCCCuts, &ana::kIsTrueSigST, systshiftdw, &wgtST);
      vUnfoldingVars[syst_dw]["ENu"]    = new ana::xsec::UnfoldingVariable(loader, 
        &ana::kTrueEStandardAxisST, &ana::kRecoEStandardAxis,
        &ana::kAllNumuCC1DCuts, &ana::kIsTrueSig1DST, systshiftdw, &wgtST);
      vUnfoldingVars[syst_dw]["Q2"]     = new ana::xsec::UnfoldingVariable(loader,
        &ana::kTrueQ2StandardAxisST, &ana::kRecoQ2StandardAxis,
        &ana::kAllNumuCC1DCuts, &ana::kIsTrueSig1DST, systshiftdw, &wgtST);
    }
    else {
      cerr << "Attempt to run on unknown systematic. Syst: " << syst << "Not found. Skipping." << endl;
      exit(1);
    }
  } // end systematic loop

  // Call go on all loaders (no particular order, could be dumb. Hmm...)
  for(pair<std::string, ana::SpectrumLoaderBase*> loader : vloaders){
    loader.second->Go();
  }

  // Store it all.
  TFile* fOut = new TFile(outfile.c_str(), "RECREATE");
  fOut->cd();
  TDirectory* dir;

  // Loop through all systematics
  for (std::pair<std::string, std::map<std::string, ana::xsec::UnfoldingVariable*> > systVars : vUnfoldingVars){
    std::string syst = systVars.first;
    std::map<std::string, ana::xsec::UnfoldingVariable*> unfoldVars = systVars.second;

    std::string foldername = "";
    if (syst_to_folder.find(syst) != syst_to_folder.end())
      foldername = syst_to_folder[syst];
    else
      foldername = syst;

    dir = fOut->mkdir(syst.c_str());
    dir->cd();
    for (std::pair<std::string, ana::xsec::UnfoldingVariable*> unfoldVar : unfoldVars){
      std::string varName     = unfoldVar.first;
      ana::xsec::UnfoldingVariable * var = unfoldVar.second;

      TDirectory * varDir = dir->mkdir(varName.c_str());
      var->SaveSpectrums(varDir);
    }
  }
  fOut->Close();
  return;
}

//////////////////////////////////////////////////////////////////
/// Process a single \a syst and store it in \a outfile
//////////////////////////////////////////////////////////////////
void specprod_systematics(const std::string syst, std::string outfile = ""){
  if (outfile.empty())
    outfile = "fSystSpec_" + syst + ".root";
  if (syst == "all")
    specprod_systematics({
      "cv", "muoneup", "muonedw", "angleup", "angledw", "2kA",                     
      "02mmBeamSpotSize",
      "1mmBeamShiftX",
      "1mmBeamShiftY",
      "3mmHorn1X",
      "3mmHorn1Y",
      "3mmHorn2X",
      "3mmHorn2Y",
      "7mmTargetZ",
      "MagneticFieldinDecayPipe",
      "1mmHornWater",
      "neutron",
      "lightdown", "lightup", "ckv",
      "calibneg", "calibpos", "calibshift",
      "neutron_nom", "nom_good_seventh"
      }, outfile);
  else
    specprod_systematics(std::vector<string> {syst}, outfile);
  return;
}