Nus18SystsPPFXLoad.C
Go to the documentation of this file.
1 /////////////////////////////////////////////////////////////////////
2 // PPFX systematics for Nus18
3 //
4 // Author: Sijith Edayath
5 // Modified by Mike Wallbank (wallbank@fnal.gov)
6 /////////////////////////////////////////////////////////////////////
7 
9 #include "CAFAna/Core/Loaders.h"
10 #include "CAFAna/Core/SystShifts.h"
11 #include "CAFAna/Core/Utilities.h"
12 #include "CAFAna/Cuts/Cuts.h"
13 #include "NuXAna/Cuts/NusCuts18.h"
15 #include "CAFAna/Cuts/SpillCuts.h"
23 #include "CAFAna/Systs/BeamSysts.h"
24 #include "CAFAna/Systs/Systs.h"
26 #include "CAFAna/Vars/HistAxes.h"
27 #include "NuXAna/Vars/HistAxes.h"
29 
31 
32 using namespace ana;
33 
34 void Nus18SystsPPFXLoad(TString opt) {
35 
36  TH1::AddDirectory(0);
37  const Var kReweight = kXSecCVWgt2018*kPPFXFluxCVWgt;
38 
39  // Get options
40  ECAFType cafType;
41  Loaders::FluxType fluxType;
42  bool extrap = false;
43  if (opt.Contains("caf", TString::kIgnoreCase))
44  cafType = ECAFType::kFullCAF;
45  else if (opt.Contains("concat", TString::kIgnoreCase))
46  cafType = ECAFType::kNusConcat;
47  else {
48  std::cout << "ECAFType required: caf or concat" << std::endl;
49  abort();
50  }
51  if (opt.Contains("fhc", TString::kIgnoreCase))
52  fluxType = Loaders::FluxType::kFHC;
53  else if (opt.Contains("rhc", TString::kIgnoreCase))
54  fluxType = Loaders::FluxType::kRHC;
55  else {
56  std::cout << "FluxType required: fhc or rhc" << std::endl;
57  abort();
58  }
59  if (opt.Contains("extrap", TString::kIgnoreCase))
60  extrap = true;
61 
62  // Set up the loaders
63  Loaders* loaders = new Prod4NomLoaders(cafType, fluxType);
65 
66  // Define a map of samples and cuts/axes to set them up
67  // Each sample, essentially, is a way to set up the analysis
68  HistAxis AxisRecoE = kNus18AxisE;
69  std::map<std::string, std::pair<HistAxis*, std::pair<Cut*, Cut*> > > cut_samples;
70  cut_samples["Nus18"] = std::make_pair(&AxisRecoE, std::make_pair(new Cut(kNus18ND), new Cut(kNus18FD)));
71 
72  // Set up maps to the decompositions/predictions (each flavor component)
73  // Nominal maps are indexed purely by sample label
74  // Shifted maps are indexed by sample label, systematic label, then sigma of the shift
75  std::map<std::string, NDPredictionSterile*> predsND_nominal;
76  std::map<std::string, FDPredictionSterile*> predsFD_nominal;
77  std::map<std::string, PredictionSterile*> predsExtrap_nominal;
78  std::map<std::string, std::map<std::string, std::map<int, NDPredictionSterile*> > > predsND_shifted;
79  std::map<std::string, std::map<std::string, std::map<int, FDPredictionSterile*> > > predsFD_shifted;
80  std::map<std::string, std::map<std::string, std::map<int, PredictionSterile*> > > predsExtrap_shifted;
81 
82  // Get vector of reweighted PPFX universes
83  std::vector<Var> kPPFXFluxUnivWgts;
84  for (unsigned int UnivIdx = 0; UnivIdx < 100; UnivIdx++) {
85  const NuTruthVar tempPPFXWgt_NT([UnivIdx](const caf::SRNeutrinoProxy* nu) {
86  if(nu->rwgt.ppfx.vuniv[UnivIdx] <= 0)
87  return 1.f;
88  return (float)nu->rwgt.ppfx.vuniv[UnivIdx];
89  });
90  const Var kPPFXFluxUnivW = VarFromNuTruthVar(tempPPFXWgt_NT, 1);
91  kPPFXFluxUnivWgts.push_back(kPPFXFluxUnivW);
92  }
93 
94 
95  // Set up the actual decompositions/predictions
96  // sample.first = the sample label
97  // shiftlabel.first = the systematic shift label
98  // sample.second.first = the sample HistAxis
99  // sample.second.second = the sample Cut
100  for (const auto& sample : cut_samples) {
101 
102  // Get prediction generators
103  const NDPredictionGenerator gennd(*sample.second.first, *sample.second.second.first, kNoShift, kReweight);
104  const FDPredictionGenerator genfd(*sample.second.first, *sample.second.second.second, kNoShift, kReweight);
105  SterileGenerator* genExtrap;
106  if (extrap)
107  genExtrap = new SterileGenerator(*sample.second.first, kNus18BinsNumuCCAxis,
108  *sample.second.second.second, *sample.second.second.first,
109  kNumuCutND2018, kNoShift, kReweight);
110 
111  // Set up the nominal decompositions
112  // Only one is needed per sample
113  predsND_nominal[sample.first] = (NDPredictionSterile*)gennd.Generate(*loaders).release();
114  predsFD_nominal[sample.first] = (FDPredictionSterile*)genfd.Generate(*loaders).release();
115  if (extrap)
116  predsExtrap_nominal[sample.first] = (PredictionSterile*)genExtrap->Generate(*loaders).release();
117 
118  // if (extrap)
119  // delete genExtrap;
120 
121  // Set up the shifted decompositions
122  // One is needed for every sample and PPFX universe
123  for (unsigned int Idx = 0; Idx < 100; ++Idx) {
124 
125  // Get new weights based on random PPFX universe
126  const Var kPPFXFluxWgt = kPPFXFluxUnivWgts[Idx];
127  const Var kMCUnivWeight = kPPFXFluxWgt*kXSecCVWgt2018;
128 
129  // Get new prediction generators
130  const NDPredictionGenerator gennd_ppfx
131  (*sample.second.first, *sample.second.second.first, kNoShift, kMCUnivWeight);
132  const FDPredictionGenerator genfd_ppfx
133  (*sample.second.first, *sample.second.second.second, kNoShift, kMCUnivWeight);
134  SterileGenerator* genExtrap_ppfx;
135  if (extrap)
136  genExtrap_ppfx = new SterileGenerator(*sample.second.first, kNus18BinsNumuCCAxis,
137  *sample.second.second.second, *sample.second.second.first,
138  kNumuCutND2018, kNoShift, kMCUnivWeight);
139 
140  // Get new predictions
141  predsND_shifted[sample.first][Form("%d",Idx)][1]
142  = (NDPredictionSterile*)gennd_ppfx.Generate(*loaders).release();
143  predsFD_shifted[sample.first][Form("%d",Idx)][1]
144  = (FDPredictionSterile*)genfd_ppfx.Generate(*loaders).release();
145  if (extrap)
146  predsExtrap_shifted[sample.first][Form("%d",Idx)][1]
147  = (PredictionSterile*)genExtrap_ppfx->Generate(*loaders).release();
148 
149  // if (extrap)
150  // delete genExtrap_ppfx;
151 
152  }
153 
154  }
155 
156  loaders->Go();
157 
158  std::string folder = "./";
159  std::string filenm = "SystsPPFX18";
160 
161  std::string fullLocation = folder + filenm + ".root";
162  TFile* rootF = new TFile(fullLocation.c_str(), "RECREATE");
163  SaveMaps(rootF, predsND_nominal, predsND_shifted);
164  SaveMaps(rootF, predsFD_nominal, predsFD_shifted);
165  if (extrap)
166  SaveMaps(rootF, predsExtrap_nominal, predsExtrap_shifted);
167  rootF->Close();
168 
169  // delete loaders;
170 
171  return;
172 
173 }
caf::Proxy< std::vector< float > > vuniv
Definition: SRProxy.h:490
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
const Var kPPFXFluxCVWgt
Definition: PPFXWeights.h:16
std::unique_ptr< IPrediction > Generate(Loaders &loaders, const SystShifts &shiftMC=kNoShift) const override
std::unique_ptr< IPrediction > Generate(Loaders &loaders, const SystShifts &shiftMC=kNoShift) const override
Proxy for caf::SRNeutrino.
Definition: SRProxy.h:510
const Cut kNumuCutND2018
Definition: NumuCuts2018.h:41
Collection of SpectrumLoaders for many configurations.
Definition: Loaders.h:23
void SaveMaps(TDirectory *out, 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)
Save all of the objects in the input maps to the out directory/file.
Definition: PPFXHelper.h:1077
caf::Proxy< caf::SRMCReweight > rwgt
Definition: SRProxy.h:560
std::pair< Spectrum *, CheatDecomp * > make_pair(SpectrumLoaderBase &loader_data, SpectrumLoaderBase &loader_mc, HistAxis *axis, Cut *cut, const SystShifts &shift, const Var &wei)
Definition: DataMCLoad.C:336
Generates extrapolated NC predictions using ProportionalDecomp.
const Cut kNus18FD
Definition: NusCuts18.h:100
caf::Proxy< caf::SRFluxWeights > ppfx
Definition: SRProxy.h:506
const HistAxis kNus18AxisE("Energy Deposited in Scintillator (GeV)", kNus18EnergyBinning, kNus18Energy)
Axes used in Nus18 analysis by nus group.
Definition: HistAxes.h:16
void Go()
Call Go() on all the loaders.
Definition: Loaders.cxx:162
void Nus18SystsPPFXLoad(TString opt)
void SetSpillCut(const SpillCut &cut)
Definition: Loaders.cxx:121
_Cut< caf::SRProxy > Cut
Representation of a cut (selection) to be applied to a caf::StandardRecord object.
Definition: Cut.h:96
const Cut kNus18ND
Full Nus18 ND analysis selection.
Definition: NusCuts18.h:137
Take the output of an extrapolation and oscillate it as required.
const SystShifts kNoShift
Definition: SystShifts.cxx:21
OStream cout
Definition: OStream.cxx:6
Var VarFromNuTruthVar(const NuTruthVar &stv, double _default)
Definition: Var.cxx:7
const Var kXSecCVWgt2018
Definition: XsecTunes.h:49
std::vector< Loaders * > loaders
Definition: syst_header.h:386
A prediction object compatible with sterile oscillations.
const HistAxis kNus18BinsNumuCCAxis("Reconstructed Neutrino Energy (GeV)", kNus18EnergyBinning, kCCE)
Definition: HistAxes.h:17
std::unique_ptr< IPrediction > Generate(Loaders &loaders, const SystShifts &shiftMC=kNoShift) const override
Take the output of an extrapolation and oscillate it as required.
const SpillCut kStandardSpillCuts
Apply this unless you&#39;re doing something special.
Definition: SpillCuts.h:49
For nominal spectra and reweighting systs (xsec/flux)
Definition: Prod4Loaders.h:96
Generates Near Detector predictions.
ECAFType
Definition: Loaders.h:19
enum BeamMode string