syst_header.h

Go to the documentation of this file.

#pragma once

#include "CAFAna/Core/Loaders.h"
#include "CAFAna/Analysis/Prod4Loaders.h"
#include "CAFAna/Core/SystShifts.h"
#include "CAFAna/Cuts/Cuts.h"
#include "3FlavorAna/Cuts/NueCuts2018.h"
#include "3FlavorAna/Cuts/NumuCuts2018.h"
#include "3FlavorAna/Cuts/QuantileCuts.h"
#include "3FlavorAna/Prediction/PredictionGenerator3Flavor.h"
#include "3FlavorAna/Systs/NueExtrapSysts.h"
#include "CAFAna/Systs/XSecSysts.h"
#include "CAFAna/Systs/XSecSystLists.h"
#include "CAFAna/Systs/BeamSysts.h"
#include "CAFAna/Vars/GenieWeights.h"
#include "3FlavorAna/Vars/HistAxes.h"
#include "CAFAna/Vars/PPFXWeights.h"
#include "3FlavorAna/Systs/EnergySysts.h"
#include "3FlavorAna/Systs/EnergySysts2017.h"
#include "3FlavorAna/Systs/JointAna2017Systs.h"
#include "3FlavorAna/Systs/GeniePCASyst.h"
#include "3FlavorAna/Systs/DummySystStorage.h"
#include "3FlavorAna/Systs/NueAcceptSysts.h"
#include "3FlavorAna/Systs/GeniePCASyst.h"
#include "3FlavorAna/Systs/MichelTaggingSyst.h"

#include <TFile.h>
#include <TString.h>
#include <vector>
#include <iostream>


using namespace ana;

// Current structure of header (could change depending on taste and usage)
// SystGroupDef defines a group of systs with their valid shifts. 
// Right now I'm using it to define groups for shifts which have their own loaders like CalibUp etc 
// Also a Flux group which doesn't require separate loaders is also implemented

typedef std::map<TString, std::map<TString, SystShifts>> Group;

struct SystGroupDef
{
  Loaders* loader;
  TString syst_group_name;
  Group group_shifts;
};

Loaders* GetLoaders2018(const TString option,
                        const TString period="full")
{
  Loaders * loaders = new Loaders();

  // Define CAF type
  auto caftype = ana::ECAFType::kFullCAF;
  if (option.Contains("nue2018_concat")) caftype = ana::ECAFType::kNueConcat;
  if (option.Contains("numu2018_concat")) caftype = ana::ECAFType::kNumuConcat;

  // Get flux (FHC or RHC) from option
  if (!(option.Contains("FHC") || option.Contains("RHC")))
  {
    std:: cerr << "Please mention FHC or RHC in option" << std::endl;
    abort(); // unable to determine if FHC or RHC
  }
  auto flux = option.Contains("RHC") ? Loaders::kRHC : Loaders::kFHC;

  if(option.Contains("CalibUp"))
    loaders = new Prod4AbsCalibLoaders(caftype, flux, +1,
                                       period.Data(), period.Data());
  else if(option.Contains("CalibDown"))
    loaders = new Prod4AbsCalibLoaders(caftype, flux, -1,
                                       period.Data(), period.Data());
  else if(option.Contains("CalibShape"))
    loaders = new Prod4CalibShapeLoaders(caftype, flux,
                                         period.Data(), period.Data());
  else if(option.Contains("Cherenkov"))
    loaders = new Prod4CherenkovLoaders(caftype, flux,
                                        period.Data(), period.Data());
  else if(option.Contains("LightUp"))
    loaders = new Prod4LightLevelLoaders(caftype, flux, +1,
                                         period.Data(), period.Data());
  else if(option.Contains("LightDown"))
    loaders = new Prod4LightLevelLoaders(caftype, flux, -1,
                                         period.Data(), period.Data());
  else if(option.Contains("LightNom"))  // FD ONLY
  {
    loaders = new Prod4LightLevelLoaders(caftype, flux, 0,
                                         period.Data(), period.Data());

    auto temploaders = new Prod4NomLoaders(caftype, flux,
                                           period.Data(), "full");
    auto nominalNDMC = temploaders->GetLoaderPath(
      caf::kNEARDET, Loaders::kMC, ana::kBeam, Loaders::kNonSwap);
    delete temploaders;
  
    std::cout << "LightNom ONLY for FD" << std::endl
      << " --> adding loader path for standard ND nominal\n\n";
  
    loaders->SetLoaderPath(
      nominalNDMC, caf::kNEARDET, Loaders::kMC, ana::kBeam, Loaders::kNonSwap);
  
    std::cout << "Now: " << loaders->GetLoaderPath(
      caf::kNEARDET, Loaders::kMC, ana::kBeam, Loaders::kNonSwap)
      << std::endl << std::endl;
  }
  else if(option.Contains("RelCalib"))
  {
    int sign = (option.Contains("Up") ? +1:-1);
    loaders = new Prod4AbsCalibLoaders(caftype, flux, sign,
                                       period.Data(), period.Data());
    auto temploaders = new Prod4AbsCalibLoaders(
      caftype, flux, - sign, period.Data(), period.Data());
    auto oppositeNDMC = temploaders->GetLoaderPath(
      caf::kNEARDET, Loaders::kMC, ana::kBeam, Loaders::kNonSwap);
    delete temploaders;
    
    loaders->SetLoaderPath(oppositeNDMC, caf::kNEARDET,
                           Loaders::kData, ana::kBeam, Loaders::kNonSwap);
 
    std::cout << "Swapped calibration loader" << oppositeNDMC << std::endl;
  }
  else if(option.Contains("Normal")) 
    loaders = new Prod4NomLoaders(caftype, flux,
                                    period.Data(), period.Data());

  // TODO: Do we still need this?
  if (option.Contains("NoTau"))
  {
    loaders->DisableLoader(caf::kFARDET, Loaders::kMC,
                          ana::kBeam, Loaders::kTauSwap);
  }
  if(option.Contains("FakeData"))
  {
    auto temploaders = new Prod4NomLoaders (caftype, flux, period.Data(), "full");
    auto nominalNDMC = temploaders->GetLoaderPath(
      caf::kNEARDET, Loaders::kMC, ana::kBeam, Loaders::kNonSwap);
    delete temploaders;
 
    std::cout << "\n\nSwapping ND Data for fake data \n\n";
    std::cout << "Before " << loaders->GetLoaderPath(
      caf::kNEARDET, Loaders::kData, ana::kBeam, Loaders::kNonSwap) << "\n";
 
    loaders->SetLoaderPath(
      nominalNDMC, caf::kNEARDET, Loaders::kData, ana::kBeam, Loaders::kNonSwap);
 
    std::cout << "After " << loaders->GetLoaderPath(
      caf::kNEARDET, Loaders::kData, ana::kBeam, Loaders::kNonSwap)
      << std::endl << std::endl;
  }
  return loaders;
}
                    
void AddSystGroup(Loaders* loader,
                  TString syst_group_name, 
                  Group group_shifts, 
                  std::vector<SystGroupDef> &syst_group)
{
  syst_group.push_back({loader, syst_group_name, group_shifts});
}

// helper class for syst groups
class SystGroupHelper
{
public: 
  SystGroupHelper(std::vector<const ISyst*> systs, std::vector<int> shifts, bool dummy)
    : fSysts(systs), fShifts(shifts), fDummy(dummy) {}
  
  SystGroupHelper(std::vector<const ISyst*> systs)
    : fSysts(systs), fShifts({+2, +1, -1, -2}), fDummy(false) {}
  
  Group GetShifts()
  {
    Group ret;
    for(auto syst: fSysts){
      std::map<TString, SystShifts> systshifts;
      for(auto sigma: fShifts)
        systshifts.insert({fShiftNames[sigma], 
                           fDummy ? kNoShift : SystShifts(syst, sigma)
                          });
      ret.insert({TString(syst->ShortName()), systshifts});
    }
    return ret;
  }
private:
  std::map<int, TString> fShiftNames = 
    {{+1, "PlusOne"}, {-1, "MinusOne"}, 
    {0, "NoShift"}, 
    {+2, "PlusTwo"}, {-2, "MinusTwo"}};
  std::vector<const ISyst*> fSysts;
  std::vector<int> fShifts;
  bool fDummy;
};
    
// group all analysis systematics
// The names of these functions are IMPORTANT for now. The python script takes in a syst group argument like 'X' to call GetXGroupShifts(). I want to do something more airtight later. 
Group GetCalibUpGroupShifts(){return SystGroupHelper({&kAnaCalibrationSyst}, {+1}, true).GetShifts();}
Group GetCalibDownGroupShifts(){return SystGroupHelper({&kAnaCalibrationSyst}, {-1}, true).GetShifts();}
Group GetRelCalibUpGroupShifts(){return SystGroupHelper({&kAnaRelativeCalibSyst}, {+2}, true).GetShifts();}
Group GetRelCalibDownGroupShifts(){return SystGroupHelper({&kAnaRelativeCalibSyst}, {-2}, true).GetShifts();}
Group GetCalibShapeGroupShifts(){return SystGroupHelper({&kAnaCalibShapeSyst}, {+1}, true).GetShifts();}

Group GetLightUpGroupShifts(){return SystGroupHelper({&kAnaLightlevelSyst}, {+1}, true).GetShifts();}
Group GetLightDownGroupShifts(){return SystGroupHelper({&kAnaLightlevelSyst}, {-1}, true).GetShifts();}
Group GetLightNomGroupShifts(){return SystGroupHelper({&kAnaLightlevelSyst}, {0}, true).GetShifts();}
Group GetCherenkovGroupShifts(){return SystGroupHelper({&kAnaCherenkovSyst}, {+1}, true).GetShifts();}


Group GetFluxGroupShifts(){
  std::vector<const ISyst*> flux_systs;
  for(int i = 0; i < 5; i++)
    flux_systs.push_back(GetFluxPrincipals2018(i));
  return SystGroupHelper(flux_systs).GetShifts();
}

Group GetGeniePCAGroupShifts(){
  std::vector<const ISyst*> geniepca_systs;
  for(int i = 0; i < 5; i++)
    geniepca_systs.push_back(GetGeniePrincipals2018Small(i));
  return SystGroupHelper(geniepca_systs).GetShifts();
}

std::vector<const ISyst*> XSecSystListHelper(int minsyst = 0, int maxsyst = -1)
{
  std::vector<const ISyst*> geniesysts = getAllXsecSysts_2018_RPAFix();
  geniesysts.push_back(&kTauScaleSyst);
  int lastsyst = maxsyst;
  if(maxsyst < 0) lastsyst = geniesysts.size();
  return std::vector<const ISyst*>(geniesysts.begin() + minsyst, geniesysts.begin() + lastsyst);
}

Group GetXSecAllGroupShifts(){return SystGroupHelper(XSecSystListHelper()).GetShifts();}
// break XSec systs into 6 groups
Group GetXSec1GroupShifts(){return SystGroupHelper(XSecSystListHelper(0, 14)).GetShifts();}
Group GetXSec2GroupShifts(){return SystGroupHelper(XSecSystListHelper(14, 28)).GetShifts();}
Group GetXSec3GroupShifts(){return SystGroupHelper(XSecSystListHelper(28, 42)).GetShifts();}
Group GetXSec4GroupShifts(){return SystGroupHelper(XSecSystListHelper(42, 56)).GetShifts();}
Group GetXSec5GroupShifts(){return SystGroupHelper(XSecSystListHelper(56, 70)).GetShifts();}
Group GetXSec6GroupShifts(){return SystGroupHelper(XSecSystListHelper(70, -1)).GetShifts();}

// Logic for which ones are selected handled in python script
Group GetMuEnergyJointGroupShifts(){return SystGroupHelper({&kMuEScaleSyst2017}).GetShifts();}
Group GetMuEnergyNumuGroupShifts(){
  return SystGroupHelper({&kRelMuEScaleSyst2017, &kDirectRelHadEScaleSyst2017, &kDirectHadEScaleSyst2017
                        }).GetShifts();
}

Group GetNeutronGroupShifts(){return SystGroupHelper({&kNeutronVisEScalePrimariesSyst2018}).GetShifts();}

Group GetNueAcceptanceGroupShifts(){
  return SystGroupHelper({&kNueAcceptSystBkg2018FHC, &kNueAcceptSystSignalKin2018FHC,
                          &kNueAcceptSystBkg2018RHC, &kNueAcceptSystSignalKin2018RHC
                         }).GetShifts();
}

Group GetJointSmallGroupShifts(){return SystGroupHelper({getAna2018SummedSmallXsecSysts(kJointAna2018)}).GetShifts();}
Group GetNumuSmallGroupShifts(){return SystGroupHelper({getAna2018SummedSmallXsecSysts(kNumuAna2018)}).GetShifts();}

Group GetNueMichelGroupShifts(){return SystGroupHelper({&kMichelTaggingSyst2018}).GetShifts();}

// other analysis components
std::vector<Cut> GetQuantileCuts(std::string f_quantile,
                                 const HistAxis axisNumu,
                                 const HistAxis axisHadEFrac)
{
  // TFile* inFile = TFile::Open( pnfs2xrootd(f_quantile).c_str() );
  TFile* inFile = TFile::Open( f_quantile.c_str() );
  // Check that the file is valid.
  if (inFile->IsZombie()) {
    std::cout << "Problem with file " << f_quantile << std::endl;
    abort();
  }
  TH2 *FDSpec2D = (TH2*)inFile->FindObjectAny("FDSpec2D");
  const unsigned int nquantiles = 4;
  std::vector<Cut> HadEFracQuantCuts = QuantileCutsFromTH2( FDSpec2D, axisNumu, axisHadEFrac, nquantiles );
  return HadEFracQuantCuts;
}

// Helper functions to add analysis generators
// Like GetXGroupShifts, AddXGenerator is called whenever 'X' is passed to the python script
void AddNumuExtrapGenerator(std::vector<Cut> HadEFracQuantCuts,
                            const HistAxis axisNumu,
                            const Cut cutNumuFD,
                            const Cut cutNumuND,
                            Var kMCWeight,
                            std::map<TString, IPredictionGenerator*> &gens) 
{
  gens.insert({"Numu_Extrap_Quant0", 
              new NumuExtrapGenerator(axisNumu, cutNumuFD, cutNumuND, kNoShift, kMCWeight) 
              });

  for(int i = 1; i <= (int)HadEFracQuantCuts.size(); i++)
    gens.insert({TString("Numu_Extrap_Quant"+std::to_string(i)),
                new NumuExtrapGenerator(axisNumu, 
                                        cutNumuFD && HadEFracQuantCuts[i-1], 
                                        cutNumuND && HadEFracQuantCuts[i-1], 
                                        kNoShift, kMCWeight) 
                });

}

void AddNumuNoExtrapGenerator(std::vector<Cut> HadEFracQuantCuts,
                              const HistAxis axisNumu,
                              const Cut cutNumuFD,
                              Var kMCWeight, 
                              std::map<TString, IPredictionGenerator*> &gens) 
{
  gens.insert({"Numu_NoExtrap_Quant0", 
              new NoExtrapGenerator(axisNumu, cutNumuFD, kMCWeight) });
  for(int i = 1; i <= (int)HadEFracQuantCuts.size(); i++)
    gens.insert({TString("Numu_NoExtrap_Quant"+std::to_string(i)),
                new NoExtrapGenerator(axisNumu, cutNumuFD && HadEFracQuantCuts[i-1], kMCWeight) });

}

void AddNueComboExtrapGenerator(const HistAxis axisNue,
                                const HistAxis axisNumuForNue,
                                const Cut cutNueFD,
                                const Cut cutNueND,
                                const Cut cutNumuForNueND,
                                Var kMCWeight, 
                                std::map<TString, IPredictionGenerator*> &gens) 
{
  gens.insert({"Nue_ComboExtrap", 
              new NueComboExtrapGenerator(axisNue, 
                                          axisNumuForNue, 
                                          cutNueFD, 
                                          cutNueND, 
                                          cutNumuForNueND, 
                                          kNoShift, kMCWeight) 
              });

}

void AddNuePropExtrapGenerator(const HistAxis axisNue,
                               const HistAxis axisNumuForNue,
                               const Cut cutNueFD,
                               const Cut cutNueND,
                               const Cut cutNumuForNueND,
                               Var kMCWeight, 
                               std::map<TString, IPredictionGenerator*> &gens) 
{
  gens.insert({"Nue_PropExtrap", 
              new NuePropExtrapGenerator(axisNue, 
                                         axisNumuForNue, 
                                         cutNueFD, 
                                         cutNueND, 
                                         cutNumuForNueND, 
                                         kNoShift, kMCWeight) 
              });

}

void AddNueNoExtrapGenerator(const HistAxis axisNue,
                             const Cut cutNueFD,
                             Var kMCWeight,
                             std::map<TString, IPredictionGenerator*> &gens) 
{
  gens.insert({"Nue_NoExtrap", 
              new NoExtrapGenerator(axisNue, cutNueFD, kMCWeight)});

}
// Cuts
// NuMu
const Cut cutNumuND = kNumuCutND2018;
const Cut cutNumuFD = kNumuCutFD2018;

// Nue
const Cut cutNumuForNueND = kNumuCutND2018;
const Cut cutNueND = kNue2018NDCVNSsb;
const Cut cutNueFD = kNue2018FDAllSamples;  

// Axes
// NuMu
const HistAxis axisNumu = kNumuCCOptimisedAxis;
const HistAxis axisHadEFrac = kHadEFracAxis;

// Nue
const HistAxis axisNue = kNue2018Axis;
const HistAxis axisNumuForNue = kNumuCCOptimisedAxis;

// Weight
Var kMCWeight = kPPFXFluxCVWgt * kXSecCVWgt2018RPAFix;

// Create objects used in make_predictions_systs.C
std::vector<SystGroupDef> systs;
std::vector<Loaders*> loaders;
std::map<TString, IPredictionGenerator*> gens;