syst_header.h
Go to the documentation of this file.
1 #pragma once
2 
3 #include "CAFAna/Core/Loaders.h"
6 #include "CAFAna/Cuts/Cuts.h"
12 #include "CAFAna/Systs/XSecSysts.h"
14 #include "CAFAna/Systs/BeamSysts.h"
26 
27 #include <TFile.h>
28 #include <TString.h>
29 #include <vector>
30 #include <iostream>
31 
32 
33 using namespace ana;
34 
35 // Current structure of header (could change depending on taste and usage)
36 // SystGroupDef defines a group of systs with their valid shifts.
37 // Right now I'm using it to define groups for shifts which have their own loaders like CalibUp etc
38 // Also a Flux group which doesn't require separate loaders is also implemented
39 
40 typedef std::map<TString, std::map<TString, SystShifts>> Group;
41 
43 {
45  TString syst_group_name;
47 };
48 
49 Loaders* GetLoaders2018(const TString option,
50  const TString period="full")
51 {
52  Loaders * loaders = new Loaders();
53 
54  // Define CAF type
55  auto caftype = ana::ECAFType::kFullCAF;
56  if (option.Contains("nue2018_concat")) caftype = ana::ECAFType::kNueConcat;
57  if (option.Contains("numu2018_concat")) caftype = ana::ECAFType::kNumuConcat;
58 
59  // Get flux (FHC or RHC) from option
60  if (!(option.Contains("FHC") || option.Contains("RHC")))
61  {
62  std:: cerr << "Please mention FHC or RHC in option" << std::endl;
63  abort(); // unable to determine if FHC or RHC
64  }
65  auto flux = option.Contains("RHC") ? Loaders::kRHC : Loaders::kFHC;
66 
67  if(option.Contains("CalibUp"))
68  loaders = new Prod4AbsCalibLoaders(caftype, flux, +1,
69  period.Data(), period.Data());
70  else if(option.Contains("CalibDown"))
71  loaders = new Prod4AbsCalibLoaders(caftype, flux, -1,
72  period.Data(), period.Data());
73  else if(option.Contains("CalibShape"))
74  loaders = new Prod4CalibShapeLoaders(caftype, flux,
75  period.Data(), period.Data());
76  else if(option.Contains("Cherenkov"))
77  loaders = new Prod4CherenkovLoaders(caftype, flux,
78  period.Data(), period.Data());
79  else if(option.Contains("LightUp"))
80  loaders = new Prod4LightLevelLoaders(caftype, flux, +1,
81  period.Data(), period.Data());
82  else if(option.Contains("LightDown"))
83  loaders = new Prod4LightLevelLoaders(caftype, flux, -1,
84  period.Data(), period.Data());
85  else if(option.Contains("LightNom")) // FD ONLY
86  {
87  loaders = new Prod4LightLevelLoaders(caftype, flux, 0,
88  period.Data(), period.Data());
89 
90  auto temploaders = new Prod4NomLoaders(caftype, flux,
91  period.Data(), "full");
92  auto nominalNDMC = temploaders->GetLoaderPath(
94  delete temploaders;
95 
96  std::cout << "LightNom ONLY for FD" << std::endl
97  << " --> adding loader path for standard ND nominal\n\n";
98 
99  loaders->SetLoaderPath(
101 
102  std::cout << "Now: " << loaders->GetLoaderPath(
104  << std::endl << std::endl;
105  }
106  else if(option.Contains("RelCalib"))
107  {
108  int sign = (option.Contains("Up") ? +1:-1);
109  loaders = new Prod4AbsCalibLoaders(caftype, flux, sign,
110  period.Data(), period.Data());
111  auto temploaders = new Prod4AbsCalibLoaders(
112  caftype, flux, - sign, period.Data(), period.Data());
113  auto oppositeNDMC = temploaders->GetLoaderPath(
115  delete temploaders;
116 
117  loaders->SetLoaderPath(oppositeNDMC, caf::kNEARDET,
119 
120  std::cout << "Swapped calibration loader" << oppositeNDMC << std::endl;
121  }
122  else if(option.Contains("Normal"))
123  loaders = new Prod4NomLoaders(caftype, flux,
124  period.Data(), period.Data());
125 
126  // TODO: Do we still need this?
127  if (option.Contains("NoTau"))
128  {
131  }
132  if(option.Contains("FakeData"))
133  {
134  auto temploaders = new Prod4NomLoaders (caftype, flux, period.Data(), "full");
135  auto nominalNDMC = temploaders->GetLoaderPath(
137  delete temploaders;
138 
139  std::cout << "\n\nSwapping ND Data for fake data \n\n";
140  std::cout << "Before " << loaders->GetLoaderPath(
142 
143  loaders->SetLoaderPath(
145 
146  std::cout << "After " << loaders->GetLoaderPath(
148  << std::endl << std::endl;
149  }
150  return loaders;
151 }
152 
154  TString syst_group_name,
155  Group group_shifts,
156  std::vector<SystGroupDef> &syst_group)
157 {
158  syst_group.push_back({loader, syst_group_name, group_shifts});
159 }
160 
161 // helper class for syst groups
163 {
164 public:
165  SystGroupHelper(std::vector<const ISyst*> systs, std::vector<int> shifts, bool dummy)
166  : fSysts(systs), fShifts(shifts), fDummy(dummy) {}
167 
168  SystGroupHelper(std::vector<const ISyst*> systs)
169  : fSysts(systs), fShifts({+2, +1, -1, -2}), fDummy(false) {}
170 
172  {
173  Group ret;
174  for(auto syst: fSysts){
175  std::map<TString, SystShifts> systshifts;
176  for(auto sigma: fShifts)
177  systshifts.insert({fShiftNames[sigma],
178  fDummy ? kNoShift : SystShifts(syst, sigma)
179  });
180  ret.insert({TString(syst->ShortName()), systshifts});
181  }
182  return ret;
183  }
184 private:
185  std::map<int, TString> fShiftNames =
186  {{+1, "PlusOne"}, {-1, "MinusOne"},
187  {0, "NoShift"},
188  {+2, "PlusTwo"}, {-2, "MinusTwo"}};
189  std::vector<const ISyst*> fSysts;
190  std::vector<int> fShifts;
191  bool fDummy;
192 };
193 
194 // group all analysis systematics
195 // 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.
196 Group GetCalibUpGroupShifts(){return SystGroupHelper({&kAnaCalibrationSyst}, {+1}, true).GetShifts();}
197 Group GetCalibDownGroupShifts(){return SystGroupHelper({&kAnaCalibrationSyst}, {-1}, true).GetShifts();}
200 Group GetCalibShapeGroupShifts(){return SystGroupHelper({&kAnaCalibShapeSyst}, {+1}, true).GetShifts();}
201 
202 Group GetLightUpGroupShifts(){return SystGroupHelper({&kAnaLightlevelSyst}, {+1}, true).GetShifts();}
203 Group GetLightDownGroupShifts(){return SystGroupHelper({&kAnaLightlevelSyst}, {-1}, true).GetShifts();}
204 Group GetLightNomGroupShifts(){return SystGroupHelper({&kAnaLightlevelSyst}, {0}, true).GetShifts();}
205 Group GetCherenkovGroupShifts(){return SystGroupHelper({&kAnaCherenkovSyst}, {+1}, true).GetShifts();}
206 
207 
209  std::vector<const ISyst*> flux_systs;
210  for(int i = 0; i < 5; i++)
211  flux_systs.push_back(GetFluxPrincipals2018(i));
212  return SystGroupHelper(flux_systs).GetShifts();
213 }
214 
216  std::vector<const ISyst*> geniepca_systs;
217  for(int i = 0; i < 5; i++)
218  geniepca_systs.push_back(GetGeniePrincipals2018Small(i));
219  return SystGroupHelper(geniepca_systs).GetShifts();
220 }
221 
222 std::vector<const ISyst*> XSecSystListHelper(int minsyst = 0, int maxsyst = -1)
223 {
224  std::vector<const ISyst*> geniesysts = getAllXsecSysts_2018_RPAFix();
225  geniesysts.push_back(&kTauScaleSyst);
226  int lastsyst = maxsyst;
227  if(maxsyst < 0) lastsyst = geniesysts.size();
228  return std::vector<const ISyst*>(geniesysts.begin() + minsyst, geniesysts.begin() + lastsyst);
229 }
230 
232 // break XSec systs into 6 groups
239 
240 // Logic for which ones are selected handled in python script
244  }).GetShifts();
245 }
246 
248 
252  }).GetShifts();
253 }
254 
257 
259 
260 // other analysis components
261 std::vector<Cut> GetQuantileCuts(std::string f_quantile,
262  const HistAxis axisNumu,
263  const HistAxis axisHadEFrac)
264 {
265  // TFile* inFile = TFile::Open( pnfs2xrootd(f_quantile).c_str() );
266  TFile* inFile = TFile::Open( f_quantile.c_str() );
267  // Check that the file is valid.
268  if (inFile->IsZombie()) {
269  std::cout << "Problem with file " << f_quantile << std::endl;
270  abort();
271  }
272  TH2 *FDSpec2D = (TH2*)inFile->FindObjectAny("FDSpec2D");
273  const unsigned int nquantiles = 4;
274  std::vector<Cut> HadEFracQuantCuts = QuantileCutsFromTH2( FDSpec2D, axisNumu, axisHadEFrac, nquantiles );
275  return HadEFracQuantCuts;
276 }
277 
278 // Helper functions to add analysis generators
279 // Like GetXGroupShifts, AddXGenerator is called whenever 'X' is passed to the python script
280 void AddNumuExtrapGenerator(std::vector<Cut> HadEFracQuantCuts,
281  const HistAxis axisNumu,
282  const Cut cutNumuFD,
283  const Cut cutNumuND,
284  Var kMCWeight,
285  std::map<TString, IPredictionGenerator*> &gens)
286 {
287  gens.insert({"Numu_Extrap_Quant0",
288  new NumuExtrapGenerator(axisNumu, cutNumuFD, cutNumuND, kNoShift, kMCWeight)
289  });
290 
291  for(int i = 1; i <= (int)HadEFracQuantCuts.size(); i++)
292  gens.insert({TString("Numu_Extrap_Quant"+std::to_string(i)),
293  new NumuExtrapGenerator(axisNumu,
294  cutNumuFD && HadEFracQuantCuts[i-1],
295  cutNumuND && HadEFracQuantCuts[i-1],
296  kNoShift, kMCWeight)
297  });
298 
299 }
300 
301 void AddNumuNoExtrapGenerator(std::vector<Cut> HadEFracQuantCuts,
302  const HistAxis axisNumu,
303  const Cut cutNumuFD,
304  Var kMCWeight,
305  std::map<TString, IPredictionGenerator*> &gens)
306 {
307  gens.insert({"Numu_NoExtrap_Quant0",
308  new NoExtrapGenerator(axisNumu, cutNumuFD, kMCWeight) });
309  for(int i = 1; i <= (int)HadEFracQuantCuts.size(); i++)
310  gens.insert({TString("Numu_NoExtrap_Quant"+std::to_string(i)),
311  new NoExtrapGenerator(axisNumu, cutNumuFD && HadEFracQuantCuts[i-1], kMCWeight) });
312 
313 }
314 
316  const HistAxis axisNumuForNue,
317  const Cut cutNueFD,
318  const Cut cutNueND,
319  const Cut cutNumuForNueND,
320  Var kMCWeight,
321  std::map<TString, IPredictionGenerator*> &gens)
322 {
323  gens.insert({"Nue_ComboExtrap",
324  new NueComboExtrapGenerator(axisNue,
325  axisNumuForNue,
326  cutNueFD,
327  cutNueND,
328  cutNumuForNueND,
329  kNoShift, kMCWeight)
330  });
331 
332 }
333 
335  const HistAxis axisNumuForNue,
336  const Cut cutNueFD,
337  const Cut cutNueND,
338  const Cut cutNumuForNueND,
339  Var kMCWeight,
340  std::map<TString, IPredictionGenerator*> &gens)
341 {
342  gens.insert({"Nue_PropExtrap",
343  new NuePropExtrapGenerator(axisNue,
344  axisNumuForNue,
345  cutNueFD,
346  cutNueND,
347  cutNumuForNueND,
348  kNoShift, kMCWeight)
349  });
350 
351 }
352 
354  const Cut cutNueFD,
355  Var kMCWeight,
356  std::map<TString, IPredictionGenerator*> &gens)
357 {
358  gens.insert({"Nue_NoExtrap",
359  new NoExtrapGenerator(axisNue, cutNueFD, kMCWeight)});
360 
361 }
362 // Cuts
363 // NuMu
366 
367 // Nue
371 
372 // Axes
373 // NuMu
376 
377 // Nue
380 
381 // Weight
383 
384 // Create objects used in make_predictions_systs.C
385 std::vector<SystGroupDef> systs;
386 std::vector<Loaders*> loaders;
387 std::map<TString, IPredictionGenerator*> gens;
Near Detector underground.
Definition: SREnums.h:10
Group GetRelCalibDownGroupShifts()
Definition: syst_header.h:199
Far Detector at Ash River.
Definition: SREnums.h:11
std::map< TString, IPredictionGenerator * > gens
Definition: syst_header.h:387
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
const NueAcceptSystSignalKin2018RHC kNueAcceptSystSignalKin2018RHC
Group group_shifts
Definition: syst_header.h:46
std::vector< SystGroupDef > systs
Definition: syst_header.h:385
Group GetXSec4GroupShifts()
Definition: syst_header.h:236
Group GetXSec2GroupShifts()
Definition: syst_header.h:234
const HistAxis axisHadEFrac
Definition: syst_header.h:375
const SummedSyst * getAna2018SummedSmallXsecSysts(const EAnaType2018 ana)
const Var kPPFXFluxCVWgt
Definition: PPFXWeights.h:16
TString syst_group_name
Definition: syst_header.h:45
Group GetLightUpGroupShifts()
Definition: syst_header.h:202
void AddNueComboExtrapGenerator(const HistAxis axisNue, const HistAxis axisNumuForNue, const Cut cutNueFD, const Cut cutNueND, const Cut cutNumuForNueND, Var kMCWeight, std::map< TString, IPredictionGenerator * > &gens)
Definition: syst_header.h:315
void AddNumuExtrapGenerator(std::vector< Cut > HadEFracQuantCuts, const HistAxis axisNumu, const Cut cutNumuFD, const Cut cutNumuND, Var kMCWeight, std::map< TString, IPredictionGenerator * > &gens)
Definition: syst_header.h:280
std::vector< Cut > GetQuantileCuts(std::string f_quantile, const HistAxis axisNumu, const HistAxis axisHadEFrac)
Definition: syst_header.h:261
Simple record of shifts applied to systematic parameters.
Definition: SystShifts.h:20
const Cut kNumuCutND2018
Definition: NumuCuts2018.h:41
Collection of SpectrumLoaders for many configurations.
Definition: Loaders.h:23
void AddNumuNoExtrapGenerator(std::vector< Cut > HadEFracQuantCuts, const HistAxis axisNumu, const Cut cutNumuFD, Var kMCWeight, std::map< TString, IPredictionGenerator * > &gens)
Definition: syst_header.h:301
const Cut kNue2018NDCVNSsb
Definition: NueCuts2018.h:153
const NuTruthSystComponentScale kTauScaleSyst("NuTauScale","#nu_{#tau} Scale", kIsTau_NT &&!kIsNC_NT, 0.6, NuTruthSystComponentScale::kLinear)
100% uncertainty scale on taus
Definition: Systs.h:176
const HistAxis kHadEFracAxis("E_{had.} / E_{#nu}", Binning::Simple(200, 0, 1), kHadEFrac)
HistAxis that implements Hadronic Energy fraction binning used by L.Vinton to derive Hadronic Energy ...
Definition: HistAxes.h:30
const Cut cutNumuND
Definition: syst_header.h:364
Group GetGeniePCAGroupShifts()
Definition: syst_header.h:215
OStream cerr
Definition: OStream.cxx:7
void DisableLoader(caf::Det_t det, DataMC datamc, DataSource src=kBeam, SwappingConfig swap=kNonSwap)
Definition: Loaders.cxx:65
const Cut cutNumuForNueND
Definition: syst_header.h:368
Loaders for Cherenkov paths/definitions.
Definition: Prod4Loaders.h:172
Group GetCherenkovGroupShifts()
Definition: syst_header.h:205
const NeutronVisEScaleSyst2018 kNeutronVisEScalePrimariesSyst2018(false)
void AddNuePropExtrapGenerator(const HistAxis axisNue, const HistAxis axisNumuForNue, const Cut cutNueFD, const Cut cutNueND, const Cut cutNumuForNueND, Var kMCWeight, std::map< TString, IPredictionGenerator * > &gens)
Definition: syst_header.h:334
const DummyAnaSyst kAnaRelativeCalibSyst("RelativeCalib","RelCalib")
Generates FD-only predictions (no extrapolation)
Loaders::FluxType flux
BeamSyst * GetFluxPrincipals2018(int PCIdx)
Definition: BeamSysts.cxx:80
Group GetNumuSmallGroupShifts()
Definition: syst_header.h:256
Generates extrapolated Numu predictions.
Group GetMuEnergyNumuGroupShifts()
Definition: syst_header.h:242
const MuEScaleSyst2017 kMuEScaleSyst2017(0.0074, 0.0012)
Loaders for calibration shape paths/definitions.
Definition: Prod4Loaders.h:196
void AddSystGroup(Loaders *loader, TString syst_group_name, Group group_shifts, std::vector< SystGroupDef > &syst_group)
Definition: syst_header.h:153
Group GetLightDownGroupShifts()
Definition: syst_header.h:203
const DummyAnaSyst kAnaLightlevelSyst("Lightlevel","Lightlevel")
Group GetJointSmallGroupShifts()
Definition: syst_header.h:255
std::string GetLoaderPath(caf::Det_t det, DataMC datamc, DataSource src=kBeam, SwappingConfig swap=kNonSwap) const
Definition: Loaders.cxx:89
Group GetRelCalibUpGroupShifts()
Definition: syst_header.h:198
ifstream inFile
Definition: AnaPlotMaker.h:34
Loaders for absolute calibration paths/definitions.
Definition: Prod4Loaders.h:119
Group GetFluxGroupShifts()
Definition: syst_header.h:208
const DummyAnaSyst kAnaCalibrationSyst("Calibration","AbsCalib")
Loaders for light level paths/definitions.
Definition: Prod4Loaders.h:145
Generates extrapolated Nue predictions using ProportionalDecomp.
const DummyAnaSyst kAnaCherenkovSyst("Cherenkov","Cherenkov")
std::vector< const ISyst * > fSysts
Definition: syst_header.h:189
Loaders * GetLoaders2018(const TString option, const TString period="full")
Definition: syst_header.h:49
Group GetCalibShapeGroupShifts()
Definition: syst_header.h:200
const RelMuEScaleSyst2017 kRelMuEScaleSyst2017(0.0045, 10.5)
Group GetNueAcceptanceGroupShifts()
Definition: syst_header.h:249
const Cut cutNumuFD
Definition: syst_header.h:365
Group GetXSec1GroupShifts()
Definition: syst_header.h:233
const NueAcceptSystBkg2018RHC kNueAcceptSystBkg2018RHC
RHC.
Group GetNeutronGroupShifts()
Definition: syst_header.h:247
Group GetXSec6GroupShifts()
Definition: syst_header.h:238
const NueAcceptSystSignalKin2018FHC kNueAcceptSystSignalKin2018FHC
const DirectHadEScaleSyst2017 kDirectHadEScaleSyst2017(0.05)
const Cut cutNueND
Definition: syst_header.h:369
void AddNueNoExtrapGenerator(const HistAxis axisNue, const Cut cutNueFD, Var kMCWeight, std::map< TString, IPredictionGenerator * > &gens)
Definition: syst_header.h:353
const Cut cutNueFD
Definition: syst_header.h:370
const Var kXSecCVWgt2018RPAFix
Definition: XsecTunes.h:53
loader
Definition: demo0.py:10
const NueAcceptSystBkg2018FHC kNueAcceptSystBkg2018FHC
FHC.
SystGroupHelper(std::vector< const ISyst * > systs)
Definition: syst_header.h:168
const HistAxis kNumuCCOptimisedAxis("Reconstructed Neutrino Energy (GeV)", kNumuCCEOptimisedBinning, kCCE)
HistAxis that implements optimised numuCCE from L. Vinton. See docdb 16332. This was close to &#39;custC&#39;...
Definition: HistAxes.h:24
double sigma(TH1F *hist, double percentile)
GeniePCASyst * GetGeniePrincipals2018Small(int PCIdx)
const HistAxis axisNumuForNue
Definition: syst_header.h:379
const SystShifts kNoShift
Definition: SystShifts.h:115
OStream cout
Definition: OStream.cxx:6
const HistAxis axisNue
Definition: syst_header.h:378
std::map< TString, std::map< TString, SystShifts > > Group
Definition: syst_header.h:40
SystGroupHelper(std::vector< const ISyst * > systs, std::vector< int > shifts, bool dummy)
Definition: syst_header.h:165
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
Group GetXSecAllGroupShifts()
Definition: syst_header.h:231
const DirectRelHadEScaleSyst2017 kDirectRelHadEScaleSyst2017(0.05)
std::vector< const ISyst * > getAllXsecSysts_2018_RPAFix()
Loaders * loader
Definition: syst_header.h:44
Group GetShifts()
Definition: syst_header.h:171
Group GetXSec5GroupShifts()
Definition: syst_header.h:237
std::vector< Loaders * > loaders
Definition: syst_header.h:386
Group GetCalibUpGroupShifts()
Definition: syst_header.h:196
Generates extrapolated Nue predictions using Michel+BEN decomposition.
const MichelTaggingSyst2018 kMichelTaggingSyst2018
const Cut kNue2018FDAllSamples
Definition: NueCuts2018.h:77
std::vector< const ISyst * > XSecSystListHelper(int minsyst=0, int maxsyst=-1)
Definition: syst_header.h:222
std::vector< Cut > QuantileCutsFromTH2(TH2 *quantileHist, const HistAxis &independentAxis, const HistAxis &quantileAxis, const unsigned int &numQuantiles, const bool verbose)
: Do the same as the QuantileCuts function but taking in the TH2 instead of making it...
Group GetCalibDownGroupShifts()
Definition: syst_header.h:197
const Cut kNumuCutFD2018
Definition: NumuCuts2018.h:39
Group GetNueMichelGroupShifts()
Definition: syst_header.h:258
const HistAxis axisNumu
Definition: syst_header.h:374
For nominal spectra and reweighting systs (xsec/flux)
Definition: Prod4Loaders.h:96
void SetLoaderPath(const std::string &path, caf::Det_t det, DataMC datamc, DataSource src=kBeam, SwappingConfig swap=kNonSwap)
Configure loader via wildcard path.
Definition: Loaders.cxx:25
Var kMCWeight
Definition: syst_header.h:382
const HistAxis kNue2018Axis("NuE Energy / Analysis Bin", kNue2018Binning, kNue2018AnaBin)
Use this Axis for Ana2018, official Axis.
std::vector< int > fShifts
Definition: syst_header.h:190
Group GetLightNomGroupShifts()
Definition: syst_header.h:204
Group GetMuEnergyJointGroupShifts()
Definition: syst_header.h:241
const DummyAnaSyst kAnaCalibShapeSyst("CalibShape","CalibShape")
def sign(x)
Definition: canMan.py:197
Group GetXSec3GroupShifts()
Definition: syst_header.h:235