NueSystHelper.cxx

Go to the documentation of this file.

#include "3FlavorAna/Systs/NueSystHelper.h"

#include "3FlavorAna/Systs/JointAna2018Systs.h"
#include "CAFAna/Systs/Systs.h"
#include "CAFAna/Systs/XSecSystLists.h"
#include "CAFAna/Systs/DISSysts.h"
#include "CAFAna/Systs/MECSysts.h"
#include "CAFAna/Systs/RPASysts.h"
#include "CAFAna/Systs/BeamSysts.h"
#include "3FlavorAna/Systs/EnergySysts.h"
#include "3FlavorAna/Systs/EnergySysts2017.h"
#include "3FlavorAna/Systs/EnergySysts2018.h"
#include "3FlavorAna/Systs/NueExtrapSysts.h"
#include "3FlavorAna/Systs/GeniePCASyst.h" // GetGeniePrincipals2018
#include "3FlavorAna/Systs/NueAcceptSysts.h"
#include "3FlavorAna/Systs/MichelTaggingSyst.h"//MicheclSystematic
#include "3FlavorAna/Systs/WrongSignSysts.h"

#include "3FlavorAna/Systs/DummySystStorage.h"

#include <iostream>

namespace ana
{

// Functions compied from syst_variations

  std::vector <const ISyst*> GetSystematics2018(const TString option)
  {
    std::vector <const ISyst*> ret;
    std::vector<const ISyst*> tmp;

    if (option.Contains("RelativeCalib") || option.Contains("All"))
      ret.push_back(&kAnaRelativeCalibSyst);
    if (option.Contains("Calibration") || option.Contains("All"))
      ret.push_back(&kAnaCalibrationSyst);
    if (option.Contains("CalibShape") || option.Contains("All"))
      ret.push_back(&kAnaCalibShapeSyst);
    if (option.Contains("LightLevel") || option.Contains("All"))
      ret.push_back(&kAnaLightlevelSyst);
    if (option.Contains("Cherenkov") || option.Contains("All"))
      ret.push_back(&kAnaCherenkovSyst);

    if(option.Contains("BeamTransport"))
    {
      std::cerr << "Option 'BeamTransport' is deprecated." << std::endl;
      std::cerr << "'FluxSyst' includes both 'BeamTransport' and 'PPFXSyst'"
        << std::endl;
      abort(); // outdated option
    }
    if(option.Contains("BeamAllTransport"))
    {
      std::cerr << "Option 'BeamAllTransport' is deprecated." << std::endl;
      std::cerr << "'FluxSyst' includes both 'BeamTransport' and 'PPFXSyst'"
        << std::endl;
      abort(); // outdated option
    }
    if(option.Contains("PPFXSyst"))
    {
      std::cerr << "Option 'PPFXSyst' is deprecated." << std::endl;
      std::cerr << "'FluxSyst' includes both 'BeamTransport' and 'PPFXSyst'"
        << std::endl;
      abort(); // outdated option
    }

    if(option.Contains("FluxSyst") || option.Contains("All"))
    {
      // Selecting 8 principal components to include in the big PDF
      // In the final fit we expect this number to be smaller ~4-5
      // TODO: Add functionality to specify desired number of components at
      // runtime.
      for (int i = 0; i < 8; i++)
        ret.push_back(GetFluxPrincipals2018(i));
    }
    if(option.Contains("GENIEPCASyst") || option.Contains("All"))
    {
      // Selecting 5 principal components to include in the big PDF
      // This matches the number we think we will include in the final fit
      // TODO: Add functionality to specify desired number of components at
      // runtime.
      for (int i = 0; i < 5; i++)
        ret.push_back(GetGeniePrincipals2018Small(i));
    }
    if(option.Contains("XSecSyst") || option.Contains("All"))
    {
      std::vector<const ISyst*> geniesysts = getAllXsecSysts_2018();
      // After 2018 analysis some new systematics were added. This returns the
      // 2018 systematics list plus 5 new systematics:
      // - kRPACCQEEnhSyst2019
      // - kRPACCQESuppSyst2019
      // - kRPARESSyst2019
      // - kMECShapeSyst2018RPAFixNu
      // - kMECShapeSyst2018RPAFixAntiNu
      // ONLY use if you are also applying the RPAFix weights
      if (option.Contains("RPAFix")){
        std::cout<<"Replacing 2018RPA systs by RPAFix..."<<std::endl;
        geniesysts.clear();
        geniesysts = getAllXsecSysts_2018_RPAFix();
      }
      // The tau systematic description is in the Nue executive summary,
      // p. 25 of doc-26699-v9
      geniesysts.push_back(&kTauScaleSyst);
      int minsyst=0, maxsyst=geniesysts.size();
      if(option.Contains("XSecSyst_Frac"))
      {
        minsyst = 0;
        maxsyst = 3;
        std::cout << "\n\nOnly looking at a small range: " << std::endl;
      }
      if(option.Contains("XSecSyst_RPATest"))
      {
        std::cout << "\n\nWe are not using kRPARESSyst2018Test" << std::endl;
      }
      if(option.Contains("XSecSyst_COH"))
      {
        minsyst = 53;
        maxsyst = 55;
        std::cout << "\n\nOnly looking at COHCC and COHNC (which replaced MaCOHpi and R0COHpi): " << std::endl;
      }
      if(option.Contains("XSecSyst_1of6")){ minsyst= 0; maxsyst=14;}
      if(option.Contains("XSecSyst_2of6")){ minsyst=14; maxsyst=28;}
      if(option.Contains("XSecSyst_3of6")){ minsyst=28; maxsyst=42;}
      if(option.Contains("XSecSyst_4of6")){ minsyst=42; maxsyst=56;}
      if(option.Contains("XSecSyst_5of6")){ minsyst=56; maxsyst=70;}
      if(option.Contains("XSecSyst_6of6")){ minsyst=70; maxsyst=87;}
      if(option.Contains("All")){ minsyst=0; maxsyst=87;}

      std::cout << minsyst << " -- " << maxsyst << "\n\n";

      for(int i = minsyst; i < maxsyst; ++i)
      {
        if (i >= (int)geniesysts.size()) break;
        ret.push_back(geniesysts[i]);
      }
    }
    if (option.Contains("MuEnergySyst") || (option.Contains("All")))
    {
      // 1.0% abs (FD + ND) muon E scale applied to muon track length
      ret.push_back(&kMuEScaleSyst2017);  // Nue and Numu
      if(option.Contains("numu", TString::kIgnoreCase))  // Numu only!
      {
        // .76% rel (FD only) muon E scale applied to muon track length
        ret.push_back(&kRelMuEScaleSyst2017);
        // 5% hadronic energy scale shift.
        ret.push_back(&kDirectHadEScaleSyst2017);
        // The relative hadronic energy scale shift
        ret.push_back(&kDirectRelHadEScaleSyst2017);
      }
    }

    //Neutron systematic calculation
    if(option.Contains("NeutronSyst")||(option.Contains("All")))
      {
        ret.push_back(&kNeutronVisEScalePrimariesSyst2018);
      }


    if(option.Contains ("NueExtrapSyst") ||
       (option.Contains("All") && option.Contains("nue")))
    {
      tmp = {
        //&kNueExtrapSystSignalKin2017,
        //&kNueExtrapSystBkg2017
      };
      for (auto entry:tmp)
        ret.push_back(entry);
    }
    if (option.Contains("SummedSmall"))
    {
      const ISyst* small  = getAna2018SummedSmallXsecSysts(kJointAna2018);
      ret.push_back(small);
    }
    if (option.Contains("NuMuSmall"))
    {
      const ISyst* small = getAna2018SummedSmallXsecSysts(kNumuAna2018);
      ret.push_back(small);
    }
    // Add acceptance sytematics
    if (option.Contains("AcceptanceSyst")|| (option.Contains("All") && option.Contains("nue")))
      {
        tmp = {
          &kNueAcceptSystBkg2018FHC,  // background, FHC
          &kNueAcceptSystSignalKin2018FHC,  // signal, FHC
          &kNueAcceptSystBkg2018RHC,  // background, RHC
          &kNueAcceptSystSignalKin2018RHC,  // signal, RHC
        };
        for (auto entry:tmp)
          ret.push_back(entry);
      }

    // Add MichelTagging sytematics
    if (option.Contains("MichelTaggingSyst")|| (option.Contains("All") && option.Contains("nue")))
      {
        ret.push_back(&kMichelTaggingSyst2018);
      }

    // Wrong sign scale systematic
    if(option.Contains("WrongSignScale"))
      {
        ret.push_back(&kWrongSignScale);
        ret.push_back(&kWrongSignScale100);
      }

    // Wrong sign scale by energy bin
    if(option.Contains("WrongSignEnergyBin"))
      {
        ret.push_back(&kWrongSignEnergyBin);
        ret.push_back(&kWrongSignEnergyBinRHC);
      }


    return ret;
  }// std::vector <const ISyst*> GetSystematics2018

  std::vector <ShiDef> GetShifts2018(const TString option)
  {
    std::vector <ShiDef> ret;
    TString topstr, lowstr;

    // Add nominal
    if (option.Contains("All") || option.Contains("Nominal"))
    {
      topstr = "Nominal"; lowstr= "noShift";
      ret.push_back({topstr, lowstr, kNoShift});
    }
    if(option.Contains("FluxSyst")           ||
       option.Contains("GENIEPCASyst")       ||
       option.Contains("XSecSyst")           ||
       option.Contains("MuEnergySyst")       ||
       option.Contains("NueExtrapSyst")      ||
       option.Contains("AllSmall")           ||
       option.Contains("NuMuSmall")          ||
       option.Contains("AcceptanceSyst")     ||
       option.Contains("NeutronSyst")        ||
       option.Contains("MichelTaggingSyst")  ||
       option.Contains("WrongSignScale")     ||
       option.Contains("WrongSignEnergyBin") ||
       option.Contains("All"))
    {
      auto systs = GetSystematics2018(option);

      for (auto syst:systs){
        for (auto sigma:{+2,+1,-1,-2})
        {
          TString sig_str = (sigma>0?"plus":"minus");
          if(abs(sigma)==1) sig_str+="One";
          if(abs(sigma)==2) sig_str+="Two";
          ret.push_back({syst->ShortName(),sig_str,SystShifts(syst,sigma)});
          std::cout << (TString) "Added " + syst->ShortName() + "_" + sig_str
            << std::endl;
        }//end sigmas
        std::cout << std::endl;
      }//end systs
    }//if FluxSyst || GENIEPCASyst || XSecSyst || NumuEnergySyst ||
     //NueExtrapSyst || AllSmall || NuMuSmall || WrongSign || DoubleWS

    if (option.Contains("CalibrationUp") || option.Contains("All"))
    {
      topstr = "Calibration"; lowstr ="plusOne";
      ret.push_back({topstr, lowstr, kNoShift});
    }
    if(option.Contains("CalibrationDown") || option.Contains("All"))
    {
      topstr = "Calibration"; lowstr ="minusOne";
      ret.push_back({topstr, lowstr, kNoShift});
    }
    if(option.Contains("RelativeCalibUp") || option.Contains("All"))
    {
      topstr = "RelativeCalib"; lowstr ="plusTwo";
      ret.push_back({topstr, lowstr, kNoShift});
    }
    if(option.Contains("RelativeCalibDown") || option.Contains("All"))
    {
      topstr = "RelativeCalib"; lowstr ="minusTwo";
      ret.push_back({topstr, lowstr, kNoShift});
    }
    if(option.Contains("CalibShape") || option.Contains("All"))
    {
      topstr = "CalibShape"; lowstr = "plusOne";
      ret.push_back({topstr, lowstr, kNoShift});
    }
    if(option.Contains("Cherenkov") || option.Contains("All"))
    {
      topstr = "Cherenkov"; lowstr = "plusOne";
      ret.push_back({topstr, lowstr, kNoShift});
    }
    if(option.Contains("LightLevelUp") || option.Contains("All"))
    {
      topstr = "Lightlevel"; lowstr = "plusOne";
      ret.push_back({topstr, lowstr, kNoShift});
    }
    if(option.Contains("LightLevelDown") || option.Contains("All"))
    {
      topstr = "Lightlevel"; lowstr = "minusOne";
      ret.push_back({topstr, lowstr, kNoShift});
    }
    if(option.Contains("LightLevelNom") || option.Contains("All"))
    {
      topstr = "Lightlevel"; lowstr = "noShift";
      ret.push_back({topstr, lowstr, kNoShift});
    }
    if(option.Contains("CalibrationUp_2sigmaNoTau") ||
      option.Contains("All"))
    {
      topstr = "Calibration"; lowstr = "plusTwo";
      ret.push_back({topstr, lowstr, kNoShift});
    }
    if(option.Contains("CalibrationDown_2sigmaNoTau") ||
      option.Contains("All"))
    {
      topstr = "Calibration"; lowstr = "minusTwo";
      ret.push_back({topstr, lowstr, kNoShift});
    }
    return ret;
  } //std::vector <ShiDef> GetShifts2018

  std::vector<SNameDef> GetShiftNames2018(const TString option)
  {
    std::vector <SNameDef> names;

    if (option.Contains("All"))
    {
      names.push_back({"SumSmallXSecJoint2018",
                      {"plusOne","minusOne","plusTwo","minusTwo"}});
      names.push_back({"extrap_signalkin",
                      {"plusOne","minusOne","plusTwo","minusTwo"}});
      names.push_back({"extrap_bkg",
                      {"plusOne","minusOne","plusTwo","minusTwo"}});
    }
    if((option == "Calibration") || option.Contains("All"))
      names.push_back({"Calibration", {"plusOne","minusOne"}});
    if((option == "RelativeCalib") || option.Contains("All"))
      names.push_back({"RelativeCalib",{"plusTwo","minusTwo"}});
    if((option == "CalibShape") || option.Contains("All"))
      names.push_back({"CalibShape", {"plusOne"}});
    if((option == "Cherenkov") || option.Contains("All"))
      names.push_back({"Cherenkov", {"plusOne"}});
    if((option == "Lightlevel") || option.Contains("All"))
      names.push_back({"Lightlevel", {"plusOne","minusOne"}});

    if(option.Contains("All") && option.Contains("Reduced") )
    {
      std::vector <const ISyst*> systs;
      if(!option.Contains("Super") )
      {
        systs = {&kMECq0ShapeSyst2017,
                 &kRPARESSyst2017,
                 GetGenieKnobSyst(rwgt::fReweightMaCCRES),
                 GetGenieKnobSyst(rwgt::fReweightMaNCRES),
                 GetGenieKnobSyst(rwgt::fReweightMvCCRES),
                 &kRPACCQEEnhSyst2017,
                 &kDISvnCC1pi_2017,
                 &kMAQEGenieReducedSyst2017,
                 &kRadCorrNue,
                 &kRadCorrNuebar,
                 &k2ndClassCurrs,
                 getAna2018SummedSmallXsecSysts(kJointAna2018)};
        for (int i =0; i<5; i++)
        {
          systs.push_back(GetPPFXPrincipals(i));
        }
      }
      // Some nue specific ones...
      if(option.Contains("Nue") )
      {
        systs.push_back(&kNueExtrapSystSignalKin2017);
        systs.push_back(&kNueExtrapSystBkg2017);
      }
      if(option.Contains("Numu") )
      {
        systs.push_back(&kMuEScaleSyst2017);
        systs.push_back(&kRelMuEScaleSyst2017);
        systs.push_back(&kDirectHadEScaleSyst2017);
        systs.push_back(&kDirectRelHadEScaleSyst2017);
        systs.push_back(getAna2018SummedSmallXsecSysts(kNumuAna2018));
      }
      // Want kBeamAll
      systs.push_back(&kBeamAllTransport);
      // Fill names with the systs we've just added..
      for (auto const & syst:systs)
        names.push_back(
          {syst->ShortName(), {"plusOne","minusOne","plusTwo","minusTwo"}});
    }// if "All" && "Reduced"
    if (option.Contains("All"))
    {
      for(auto c:{"FluxSyst", "NumuEnergySyst",
                  "NueExtrapSyst","AllSmall","XSecSyst"})
      { //nue wants XSec to be last
        auto systs = GetSystematics2018(c);
        for (auto const & syst:systs)
        names.push_back(
          {syst->ShortName(),{"plusOne","minusOne","plusTwo","minusTwo"}});
      }//end beam or xsec
    }
    else
    {
      auto systs = GetSystematics2018(option);
      for (auto const & syst:systs)
        names.push_back(
          {syst->ShortName(), {"plusOne","minusOne","plusTwo","minusTwo"}});
    }
    return names;
  }


}