Nus18SystsKaonLoad.C
Go to the documentation of this file.
1 /////////////////////////////////////////////////////////////////////
2 // Kaon 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"
18 #include "CAFAna/Extrap/ExtrapSterile.h"
24 #include "CAFAna/Systs/BeamSysts.h"
25 #include "CAFAna/Systs/Systs.h"
27 #include "CAFAna/Vars/HistAxes.h"
28 #include "NuXAna/Vars/HistAxes.h"
30 
32 
33 using namespace ana;
34 
35 /// Run systematics for Kaons
36 void Nus18SystsKaonLoad(TString opt) {
37 
38  TH1::AddDirectory(0);
39  const Var kReweight = kXSecCVWgt2018*kPPFXFluxCVWgt;
40 
41  // Get options
42  ECAFType cafType;
43  Loaders::FluxType fluxType;
44  bool extrap = false;
45  if (opt.Contains("caf", TString::kIgnoreCase))
46  cafType = ECAFType::kFullCAF;
47  else if (opt.Contains("concat", TString::kIgnoreCase))
48  cafType = ECAFType::kNusConcat;
49  else {
50  std::cout << "ECAFType required: caf or concat" << std::endl;
51  abort();
52  }
53  if (opt.Contains("fhc", TString::kIgnoreCase))
54  fluxType = Loaders::FluxType::kFHC;
55  else if (opt.Contains("rhc", TString::kIgnoreCase))
56  fluxType = Loaders::FluxType::kRHC;
57  else {
58  std::cout << "FluxType required: fhc or rhc" << std::endl;
59  abort();
60  }
61  if (opt.Contains("extrap", TString::kIgnoreCase))
62  extrap = true;
63 
64  // Set up the loaders
65  Loaders* loaders = new Prod4NomLoaders(cafType, fluxType);
67 
68  // Get the labels for Kaon systematics
69  std::map<std::string, std::string> shiftLabels;
70  shiftLabels["Kaon"] = "Kaon";
71  const Var kKaonP = kReweight*kNus18BeamWeightP;
72  const Var kKaonM = kReweight*kNus18BeamWeightM;
73 
74  // Define map of systematic objects, indexed by same labels as previous map
75  std::map<std::string, std::map<int, const Var*> > systematics;
76  systematics["Kaon"][-1] = &kKaonM;
77  systematics["Kaon"][+1] = &kKaonP;
78 
79  // Define a map of samples and cuts/axes to set them up
80  // Each sample, essentially, is a way to set up the analysis
81  HistAxis AxisRecoE = kNus18AxisE;
82  std::map<std::string, std::pair<HistAxis*, std::pair<Cut*, Cut*> > > cut_samples;
83  cut_samples["Nus18"] = std::make_pair(&AxisRecoE, std::make_pair(new Cut(kNus18ND), new Cut(kNus18FD)));
84 
85  // Set up maps to the decompositions/predictions (each flavor component)
86  // Nominal maps are indexed purely by sample label
87  // Shifted maps are indexed by sample label, systematic label, then sigma of the shift
88  std::map<std::string, NDPredictionSterile*> predsND_nominal;
89  std::map<std::string, FDPredictionSterile*> predsFD_nominal;
90  std::map<std::string, PredictionSterile*> predsExtrap_nominal;
91  std::map<std::string, std::map<std::string, std::map<int, NDPredictionSterile*> > > predsND_shifted;
92  std::map<std::string, std::map<std::string, std::map<int, FDPredictionSterile*> > > predsFD_shifted;
93  std::map<std::string, std::map<std::string, std::map<int, PredictionSterile*> > > predsExtrap_shifted;
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  // syst.first = the integer number of sigmas for the shift
101  // syst.second = the systematic being shifted
102  for (const auto& sample : cut_samples) {
103 
104  // const NDPredictionGenerator gennd(*sample.second.first, *sample.second.second.first, kNoShift, kReweight);
105  // const FDPredictionGenerator genfd(*sample.second.first, *sample.second.second.second, kNoShift, kReweight);
106  SterileGenerator* genExtrap;
107  // if (extrap) {
108  genExtrap = new SterileGenerator(*sample.second.first, kNus18BinsNumuCCAxis,
109  *sample.second.second.second, *sample.second.second.first,
110  kNumuCutND2018, kNoShift, kReweight);
111  // ProportionalDecomp* decompNC_nominal = new ProportionalDecomp
112  // (loaders->GetLoader(caf::Det_t::kNEARDET, Loaders::DataMC::kMC),
113  // loaders->GetLoader(caf::Det_t::kNEARDET, Loaders::DataMC::kData),
114  // *sample.second.first, *sample.second.second.first,
115  // kNoShift, kNoShift, kReweight);
116  // ProportionalDecomp* decompNuMu_nominal = new ProportionalDecomp
117  // (loaders->GetLoader(caf::Det_t::kNEARDET, Loaders::DataMC::kMC),
118  // loaders->GetLoader(caf::Det_t::kNEARDET, Loaders::DataMC::kData),
119  // kNus18BinsNumuCCAxis, kNumuCutND2018,
120  // kNoShift, kNoShift, kReweight);
121  // ModularExtrapSterile* extrap_nominal = new ModularExtrapSterile
122  // (ModularExtrapSterile::NCDisappearance
123  // (loaders->GetLoader(caf::Det_t::kNEARDET, Loaders::DataMC::kMC),
124  // loaders->GetLoader(caf::Det_t::kFARDET, Loaders::DataMC::kMC, DataSource::kBeam, Loaders::SwappingConfig::kFluxSwap),
125  // loaders->GetLoader(caf::Det_t::kFARDET, Loaders::DataMC::kMC, DataSource::kBeam, Loaders::SwappingConfig::kNonSwap),
126  // loaders->GetLoader(caf::Det_t::kFARDET, Loaders::DataMC::kMC, DataSource::kBeam, Loaders::SwappingConfig::kTauSwap),
127  // *decompNC_nominal, *decompNuMu_nominal,
128  // *sample.second.first, kNus18BinsNumuCCAxis,
129  // *sample.second.second.second, *sample.second.second.first, kNumuCutND2018,
130  // kNoShift, kReweight));
131  // predsExtrap_nominal[sample.first] = new PredictionSterile(extrap_nominal);
132  // }
133 
134  // Set up the nominal decompositions
135  // Only one is needed per sample
136  // predsND_nominal[sample.first] = (NDPredictionSterile*)gennd.Generate(*loaders).release();
137  // predsFD_nominal[sample.first] = (FDPredictionSterile*)genfd.Generate(*loaders).release();
138  // if (extrap) {
139  predsExtrap_nominal[sample.first] = (PredictionSterile*)genExtrap->Generate(*loaders).release();
140  // }
141 
142  // if (extrap)
143  // delete genExtrap;
144 
145  // Set up the shifted decompositions
146  // One is needed for every sample, systematic, and sigma level
147  for(const auto& shiftLabel : shiftLabels) {
148  for(const auto& syst : systematics[shiftLabel.first]) {
149 
150  // // Create new prediction generators with correct reweights
151  // const NDPredictionGenerator gennd_kaon
152  // (*sample.second.first, *sample.second.second.first, kNoShift, *syst.second);
153  // const FDPredictionGenerator genfd_kaon
154  // (*sample.second.first, *sample.second.second.second, kNoShift, *syst.second);
155  SterileGenerator* genExtrap_kaon;
156  // if (extrap) {
157  genExtrap_kaon = new SterileGenerator(*sample.second.first, kNus18BinsNumuCCAxis,
158  *sample.second.second.second, *sample.second.second.first,
159  kNumuCutND2018, kNoShift, *syst.second);
160  // ProportionalDecomp* decompNC_shifted = new ProportionalDecomp
161  // (loaders->GetLoader(caf::Det_t::kNEARDET, Loaders::DataMC::kMC),
162  // loaders->GetLoader(caf::Det_t::kNEARDET, Loaders::DataMC::kData),
163  // *sample.second.first, *sample.second.second.first,
164  // kNoShift, kNoShift, *syst.second);
165  // ProportionalDecomp* decompNuMu_shifted = new ProportionalDecomp
166  // (loaders->GetLoader(caf::Det_t::kNEARDET, Loaders::DataMC::kMC),
167  // loaders->GetLoader(caf::Det_t::kNEARDET, Loaders::DataMC::kData),
168  // kNus18BinsNumuCCAxis, kNumuCutND2018,
169  // kNoShift, kNoShift, *syst.second);
170  // ModularExtrapSterile* extrap_shifted = new ModularExtrapSterile
171  // (ModularExtrapSterile::NCDisappearance
172  // (loaders->GetLoader(caf::Det_t::kNEARDET, Loaders::DataMC::kMC),
173  // loaders->GetLoader(caf::Det_t::kFARDET, Loaders::DataMC::kMC, DataSource::kBeam, Loaders::SwappingConfig::kFluxSwap),
174  // loaders->GetLoader(caf::Det_t::kFARDET, Loaders::DataMC::kMC, DataSource::kBeam, Loaders::SwappingConfig::kNonSwap),
175  // loaders->GetLoader(caf::Det_t::kFARDET, Loaders::DataMC::kMC, DataSource::kBeam, Loaders::SwappingConfig::kTauSwap),
176  // *decompNC_shifted, *decompNuMu_shifted,
177  // *sample.second.first, kNus18BinsNumuCCAxis,
178  // *sample.second.second.second, *sample.second.second.first, kNumuCutND2018,
179  // kNoShift, *syst.second));
180  // predsExtrap_shifted[sample.first][shiftLabel.first][syst.first] = new PredictionSterile(extrap_shifted);
181  // }
182 
183  // Shift the spectra
184  // predsND_shifted[sample.first][shiftLabel.first][syst.first]
185  // = (NDPredictionSterile*)gennd_kaon.Generate(*loaders).release();
186  // predsFD_shifted[sample.first][shiftLabel.first][syst.first]
187  // = (FDPredictionSterile*)genfd_kaon.Generate(*loaders).release();
188  // if (extrap) {
189  predsExtrap_shifted[sample.first][shiftLabel.first][syst.first]
190  = (PredictionSterile*)genExtrap_kaon->Generate(*loaders).release();
191  // }
192 
193  // if (extrap)
194  // delete genExtrap_kaon;
195 
196  }
197  }
198 
199  }
200 
201  loaders->Go();
202 
203  std::string folder = "./";
204  std::string filenm = "SystsKaon18";
205 
206  std::string fullLocation = folder + filenm + ".root";
207  TFile* rootF = new TFile(fullLocation.c_str(), "RECREATE");
208  // SaveMaps(rootF, predsND_nominal, predsND_shifted);
209  // SaveMaps(rootF, predsFD_nominal, predsFD_shifted);
210  // if (extrap)
211  SaveMaps(rootF, predsExtrap_nominal, predsExtrap_shifted);
212  rootF->Close();
213 
214  // delete loaders;
215 
216  return;
217 
218 }
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
const Var kNus18BeamWeightM
Definition: NusVars.h:90
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
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
const Var kNus18BeamWeightP
Definition: NusVars.h:87
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 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
void Nus18SystsKaonLoad(TString opt)
Run systematics for Kaons.
const Cut kNus18ND
Full Nus18 ND analysis selection.
Definition: NusCuts18.h:137
const SystShifts kNoShift
Definition: SystShifts.cxx:22
OStream cout
Definition: OStream.cxx:6
const Var kXSecCVWgt2018
Definition: XsecTunes.h:49
std::vector< Loaders * > loaders
Definition: syst_header.h:386
bool systematics
Definition: fnexvscaf.C:31
A prediction object compatible with sterile oscillations.
const HistAxis kNus18BinsNumuCCAxis("Reconstructed Neutrino Energy (GeV)", kNus18EnergyBinning, kCCE)
Definition: HistAxes.h:17
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
ECAFType
Definition: Loaders.h:19
enum BeamMode string