Nus18SystsTauLoad.C
Go to the documentation of this file.
1 /////////////////////////////////////////////////////////////////////
2 // Tau 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 /// Run systematics for Tau scale
35 void Nus18SystsTauLoad(TString opt) {
36 
37  TH1::AddDirectory(0);
38  const Var kReweight = kXSecCVWgt2018*kPPFXFluxCVWgt;
39 
40  // Get options
41  ECAFType cafType;
42  Loaders::FluxType fluxType;
43  bool extrap = false;
44  if (opt.Contains("caf", TString::kIgnoreCase))
45  cafType = ECAFType::kFullCAF;
46  else if (opt.Contains("concat", TString::kIgnoreCase))
47  cafType = ECAFType::kNusConcat;
48  else {
49  std::cout << "ECAFType required: caf or concat" << std::endl;
50  abort();
51  }
52  if (opt.Contains("fhc", TString::kIgnoreCase))
53  fluxType = Loaders::FluxType::kFHC;
54  else if (opt.Contains("rhc", TString::kIgnoreCase))
55  fluxType = Loaders::FluxType::kRHC;
56  else {
57  std::cout << "FluxType required: fhc or rhc" << std::endl;
58  abort();
59  }
60  if (opt.Contains("extrap", TString::kIgnoreCase))
61  extrap = true;
62 
63  // Set up the loaders
64  Loaders* loaders = new Prod4NomLoaders(cafType, fluxType);
66 
67  // Vector of sigmas, levels to set each systematic
68  const std::vector<int> sigmas{-1, +1};
69 
70  // Get the labels for Tau systematics
71  std::map<std::string, std::string> shiftLabels;
72  shiftLabels["Tau"] = "Tau";
73  const Var kTauP = kReweight*kFDTauWgtP;
74  const Var kTauM = kReweight*kFDTauWgtM;
75 
76  // Define map of systematic objects, indexed by same labels as previous map
77  std::map<std::string, std::map<int, const Var*> > systematics;
78  systematics["Tau"][-1] = &kTauM;
79  systematics["Tau"][+1] = &kTauP;
80 
81  // Define a map of samples and cuts/axes to set them up
82  // Each sample, essentially, is a way to set up the analysis
83  HistAxis AxisRecoE = kNus18AxisE;
84  std::map<std::string, std::pair<HistAxis*, std::pair<Cut*, Cut*> > > cut_samples;
85  cut_samples["Nus18"] = std::make_pair(&AxisRecoE, std::make_pair(new Cut(kNus18ND), new Cut(kNus18FD)));
86 
87  // Set up maps to the decompositions/predictions (each flavor component)
88  // Nominal maps are indexed purely by sample label
89  // Shifted maps are indexed by sample label, systematic label, then sigma of the shift
90  std::map<std::string, NDPredictionSterile*> predsND_nominal;
91  std::map<std::string, FDPredictionSterile*> predsFD_nominal;
92  std::map<std::string, PredictionSterile*> predsExtrap_nominal;
93  std::map<std::string, std::map<std::string, std::map<int, NDPredictionSterile*> > > predsND_shifted;
94  std::map<std::string, std::map<std::string, std::map<int, FDPredictionSterile*> > > predsFD_shifted;
95  std::map<std::string, std::map<std::string, std::map<int, PredictionSterile*> > > predsExtrap_shifted;
96 
97  // Set up the actual decompositions/predictions
98  // sample.first = the sample label
99  // shiftLabel.first = the systematic shift label
100  // sample.second.first = the sample HistAxis
101  // sample.second.second = the sample Cut
102  // syst.first = the integer number of sigmas for the shift
103  // syst.second = the systematic being shifted
104  for (const auto& sample : cut_samples) {
105 
106  // const NDPredictionGenerator gennd(*sample.second.first, *sample.second.second.first, kNoShift, kReweight);
107  // const FDPredictionGenerator genfd(*sample.second.first, *sample.second.second.second, kNoShift, kReweight);
108  SterileGenerator* genExtrap;
109  if (extrap)
110  genExtrap = new SterileGenerator(*sample.second.first, kNus18BinsNumuCCAxis,
111  *sample.second.second.second, *sample.second.second.first,
112  kNumuCutND2018, kNoShift, kReweight);
113 
114  // Set up the nominal decompositions
115  // Only one is needed per sample
116  // predsND_nominal[sample.first] = (NDPredictionSterile*)gennd.Generate(*loaders).release();
117  // predsFD_nominal[sample.first] = (FDPredictionSterile*)genfd.Generate(*loaders).release();
118  if (extrap)
119  predsExtrap_nominal[sample.first] = (PredictionSterile*)genExtrap->Generate(*loaders).release();
120 
121  if (extrap)
122  delete genExtrap;
123 
124  // Set up the shifted decompositions
125  // One is needed for every sample, systematic, and sigma level
126  for(const auto& shiftLabel : shiftLabels) {
127  for(const auto& syst : systematics[shiftLabel.first]) {
128 
129  // Create new prediction generators with correct reweights
130  // const NDPredictionGenerator gennd_tau
131  // (*sample.second.first, *sample.second.second.first, kNoShift, *syst.second);
132  // const FDPredictionGenerator genfd_tau
133  // (*sample.second.first, *sample.second.second.second, kNoShift, *syst.second);
134  SterileGenerator* genExtrap_tau;
135  if (extrap)
136  genExtrap_tau = new SterileGenerator(*sample.second.first, kNus18BinsNumuCCAxis,
137  *sample.second.second.second, *sample.second.second.first,
138  kNumuCutND2018, kNoShift, *syst.second);
139 
140  // Shift the spectra
141  // predsND_shifted[sample.first][shiftLabel.first][syst.first]
142  // = (NDPredictionSterile*)gennd_tau.Generate(*loaders).release();
143  // predsFD_shifted[sample.first][shiftLabel.first][syst.first]
144  // = (FDPredictionSterile*)genfd_tau.Generate(*loaders).release();
145  if (extrap)
146  predsExtrap_shifted[sample.first][shiftLabel.first][syst.first]
147  = (PredictionSterile*)genExtrap_tau->Generate(*loaders).release();
148 
149  if (extrap)
150  delete genExtrap_tau;
151 
152  }
153  }
154 
155  }
156 
157  loaders->Go();
158 
159  std::string folder = "./";
160  std::string filenm = "SystsTau18";
161 
162  std::string fullLocation = folder + filenm + ".root";
163  TFile* rootF = new TFile(fullLocation.c_str(), "RECREATE");
164  // SaveMaps(rootF, predsND_nominal, predsND_shifted);
165  // SaveMaps(rootF, predsFD_nominal, predsFD_shifted);
166  if (extrap)
167  SaveMaps(rootF, predsExtrap_nominal, predsExtrap_shifted);
168  rootF->Close();
169 
170  // delete loaders;
171 
172  return;
173 
174 }
175 
176 // FDExtrap* extrapfds = new FDExtrap
177 // (FDExtrap::FDExtrap_c(
178 // loaders.GetLoader(caf::kFARDET, Loaders::kMC, kBeam, Loaders::kNonSwap),
179 // loaders.GetLoader(caf::kFARDET, Loaders::kMC, kBeam, Loaders::kFluxSwap),
180 // loaders.GetLoader(caf::kFARDET, Loaders::kMC, kBeam, Loaders::kTauSwap),
181 // *sample.second.first , kNus18FD, kNoShift, kWgtP3));
182 // predsFD_nominal[sample.first] = new FDPredictionSterile( extrapfds );
183 
184 // for(const auto& shiftlabel : shiftlabels) {
185 // for(const auto& syst : systematics[shiftlabel.first]) {
186 
187 // FDExtrap* extrapfdst = new FDExtrap
188 // (FDExtrap::FDExtrap_c
189 // (
190 // loaders.GetLoader(caf::kFARDET, Loaders::kMC, kBeam, Loaders::kNonSwap),
191 // loaders.GetLoader(caf::kFARDET, Loaders::kMC, kBeam, Loaders::kFluxSwap),
192 // loaders.GetLoader(caf::kFARDET, Loaders::kMC, kBeam, Loaders::kTauSwap),
193 // *sample.second.first , kNus18FD, kNoShift, *syst.second));
194 // predsFD_shifted[sample.first][shiftlabel.first][syst.first] = new FDPredictionSterile( extrapfdst );
195 // }
196 // }
197 // }
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 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 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
const Cut kNus18ND
Full Nus18 ND analysis selection.
Definition: NusCuts18.h:137
const SystShifts kNoShift
Definition: SystShifts.cxx:21
OStream cout
Definition: OStream.cxx:6
void Nus18SystsTauLoad(TString opt)
Run systematics for Tau scale.
const Var kXSecCVWgt2018
Definition: XsecTunes.h:49
std::vector< Loaders * > loaders
Definition: syst_header.h:386
const Var kFDTauWgtP
Definition: NusVars.h:99
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
const Var kFDTauWgtM
Definition: NusVars.h:96
For nominal spectra and reweighting systs (xsec/flux)
Definition: Prod4Loaders.h:96
ECAFType
Definition: Loaders.h:19
enum BeamMode string