make_nue_filesyst_pred.C
Go to the documentation of this file.
4 #include "CAFAna/Core/HistAxis.h"
5 #include "CAFAna/Core/Spectrum.h"
8 #include "CAFAna/Cuts/Cuts.h"
9 #include "3FlavorAna/Cuts/NueCutsSecondAna.h"
12 #include "CAFAna/Cuts/SpillCuts.h"
14 #include "CAFAna/Decomp/IDecomp.h"
21 #include "CAFAna/Vars/HistAxes.h"
25 #include "CAFAna/Vars/Vars.h"
26 
28 
29 #include "TCanvas.h"
30 #include "TFile.h"
31 #include "TH1.h"
32 #include "TLegend.h"
33 #include "TPad.h"
34 
35 #include <iostream>
36 #include <memory>
37 
38 using namespace ana;
39 
40 
41 void make_nue_filesyst_pred(int systtype)
42 {
43 
44  std::vector <std::shared_ptr< Prod3DataLoaders> > systloaders;
45  ECAFType ltype = kNueConcat;
46 
47  /// A function to switch between different systematics loader class
48  switch (systtype)
49  {
50  case 0 :
51  std::cout << "Using NueNomLoaders2017" << std::endl;
52  systloaders.push_back(std::make_shared <Prod3NomLoaders> (ltype));
53  break;
54 
55  case 1 :
56  std::cout << "Using NueAbsCalibLoaders2017" << std::endl;
57  systloaders.push_back(std::make_shared <Prod3AbsCalibLoaders> (ltype, 1.0));
58  systloaders.push_back(std::make_shared <Prod3AbsCalibLoaders> (ltype, -1.0));
59  break;
60  case 2 :
61  std::cout<< "Using NueLightLevelLoaders2017" <<std::endl;
62  systloaders.push_back(std::make_shared <Prod3LightLevelLoaders> (ltype, 1.0));
63  systloaders.push_back(std::make_shared <Prod3LightLevelLoaders> (ltype, -1.0));
64  break;
65  case 3 :
66  std::cout<< "Using CherenkovLoaders2017" <<std::endl;
67  systloaders.push_back(std::make_shared <Prod3CherenkovLoaders> (ltype));
68  break;
69  case 4 :
70  std::cout<< "Using CalibShapeLoaders2017" <<std::endl;
71  systloaders.push_back(std::make_shared <Prod3CalibShapeLoaders> (ltype));
72  break;
73  default: std::cout << "Run me like <cafe 'make_nue_filesyst_pred.C+(0/1/2/3/4)'> for systematic type you want to study \n"
74  << "0 is for nominal, 1 is for calib-syst, 2 is for lightlevel-syst, 3 is for Cherenkov, 4 is for calib-shape \n"
75  << "Exiting(0)" << std::endl;
76  return;
77  };
78 
79  // Needs to be hardcoded for now, since gdb doesn't handle an array of variable size lenght elements very well.
80  const int kNumPreds = 2;
81 
82  const int kNumShifts = 2; // Needs to be hardcoded, same reason
83  const std::string shiftNames[kNumShifts] = {"plusOne", "minusOne"};
84 
85  const std::string predNames[kNumPreds] = {"pred_propDecomp", "pred_comboDecomp"};
86 
87 
88  //..........THIS SNIPPET WILL BE ELEGANT WAY TO ADD NOMINAL...................
89  //..........but gdb freaks out seeing "size()" although it is a const........
90 
91  // const int kNumShifts = systloaders.size();
92  //std::vector<std::string> shiftNames;
93  //if(kNumShifts==1) shiftNames.push_back("noShift")
94  // else
95  // {
96  // shiftNames.push_back("plusOne");
97  // shiftNames.push_back("minusOne");
98  // }
99  //...............leaving the piece here to work on later.....................
100 
101 
102  struct HistDef
103  {
104  const std::string name;
105  HistAxis axis;
106  };
107 
108  std::vector <HistDef> histdef;
109 
110  // Which variables and binning for nue and ND numu do we want to look at?
111  histdef.push_back({"shwCalE", {"Shower calE (GeV)", Binning::Simple(30.0, 0.0, 3.0), kShwCalE}});
112  histdef.push_back({"hadCalE", {"calE_{had} (GeV)", Binning::Simple(30.0, 0.0, 3.0), kHadCalE}});
113  histdef.push_back({"recoNuE", {"reco E_{#nu} (GeV)", Binning::Simple(50.0, 0.0, 5.0), kNueEnergy2017}});
114  histdef.push_back({"id_cvne", {"CVN e", Binning::Simple(40.0, 0.70, 1.1), kCVNe}});
115  histdef.push_back({"EnergyCVN2D", {"Nue Energy / Selection Bin", kNue2017Binning, kNue2017AnaBin}});
116 
117 
118  const int kNumVars = 5; // Needs to be hardcoded, same reason
119 
120  const HistAxis& axisNumu = kNumuNonQEAxisFirstAna;
121 
122 
123  // Setup a nue and a numu decomposition, extrapolation and a prediction based on the extrapolation
124  NumuDecomp* decompNumu [kNumShifts];
125  ProportionalDecomp* decompProp [kNumPreds][kNumShifts][kNumVars];
126  MichelDecomp* decompMichel[kNumPreds][kNumShifts][kNumVars];
127  BENDecomp* decompBen [kNumPreds][kNumShifts][kNumVars];
128  ModularExtrap extrapProp [kNumPreds][kNumShifts][kNumVars];
129  ModularExtrap extrapMichel[kNumPreds][kNumShifts][kNumVars];
130  PredictionExtrap* predProp [kNumPreds][kNumShifts][kNumVars];
131  PredictionExtrap* predMichel [kNumPreds][kNumShifts][kNumVars];
132 
133 
134 
135  // Loop over the 2 prediction types, 2 shifts and variables;
136  for(int predIdx = 0; predIdx < kNumPreds; ++predIdx)
137  {
138  for(int shiftIdx = 0; shiftIdx < kNumShifts; ++shiftIdx)
139  {
140  if ((systtype == 0 || systtype == 3 || systtype == 4) && shiftIdx > 0) continue;
141 
142  decompNumu[shiftIdx] = new NumuDecomp (*systloaders[shiftIdx], axisNumu, kNumuND,
144 
145  for(int varIdx = 0; varIdx < kNumVars; ++varIdx)
146  {
147  const HistAxis& histaxis = histdef[varIdx].axis;
148 
149  if (predIdx == 0)
150  {
151  decompProp[predIdx][shiftIdx][varIdx] = new ProportionalDecomp (*systloaders[shiftIdx],
152  histaxis, kNue2017NDCVNSsb, kNoShift,
153  kNoShift, kPPFXFluxCVWgt * kXSecCVWgt2017);
154 
155 
156  extrapProp[predIdx][shiftIdx][varIdx] = NueExtrap(*systloaders[shiftIdx],
157  *decompProp[predIdx][shiftIdx][varIdx],
158  *decompNumu[shiftIdx], histaxis, axisNumu,
160  kNoShift, kPPFXFluxCVWgt * kXSecCVWgt2017);
161 
162 
163  predProp[predIdx][shiftIdx][varIdx] = new PredictionExtrap (&extrapProp[predIdx][shiftIdx][varIdx]);
164 
165  }
166 
167  else
168  {
169  decompBen[predIdx][shiftIdx][varIdx] = new BENDecomp (*systloaders[shiftIdx],
170  histaxis, kNue2017NDCVNSsb, kNoShift,
171  kNoShift, kPPFXFluxCVWgt * kXSecCVWgt2017);
172 
173 
174  decompMichel[predIdx][shiftIdx][varIdx] = new MichelDecomp (*systloaders[shiftIdx],
175  histaxis, kNue2017NDCVNSsb,
176  decompBen[predIdx][shiftIdx][varIdx], kNoShift,
177  kNoShift, kPPFXFluxCVWgt * kXSecCVWgt2017);
178 
179 
180  extrapMichel[predIdx][shiftIdx][varIdx] = NueExtrap (*systloaders[shiftIdx],
181  *decompMichel[predIdx][shiftIdx][varIdx],
182  *decompNumu[shiftIdx], histaxis, axisNumu,
184  kNoShift, kPPFXFluxCVWgt * kXSecCVWgt2017);
185 
186 
187  predMichel[predIdx][shiftIdx][varIdx] = new PredictionExtrap (&extrapMichel[predIdx][shiftIdx][varIdx]);
188 
189  }
190 
191 
192  }// end of varIdx loop
193 
194  }// end of shiftIdx loop
195 
196  }// end of predIdx loop
197 
198 
199 
200  // Fill the prediction!
201  for(auto &it : systloaders)
202  {
203  it->SetSpillCut(kStandardSpillCuts);
204  it->Go();
205  }
206 
207 
208  // Saving predictions in a structural hierarchy under different pred-analysis
209  std::cout<<"making output file \n";
210  TFile* fout = new TFile(("fout_make_nue_"+systloaders[0]->GetSystString()+"_pred.root").c_str(), "RECREATE");
211 
212  for(int predIdx = 0; predIdx < 1; ++predIdx)
213  {
214  TDirectory* pred_dir = fout->mkdir(predNames[predIdx].c_str());
215  pred_dir->cd();
216 
217  TDirectory* syst_dir = pred_dir->mkdir(systloaders[predIdx]->GetSystString().c_str());
218  syst_dir->cd();
219 
220  for(int shiftIdx = 0; shiftIdx < kNumShifts; ++shiftIdx)
221  {
222  if ((systtype == 0 || systtype == 3 || systtype == 4) && shiftIdx > 0) continue;
223 
224  TDirectory* shift_dir;
225  if (systtype == 0)
226  {
227  shift_dir = syst_dir->mkdir("noShift");
228  shift_dir->cd();
229  }
230  else
231  {
232  shift_dir = syst_dir->mkdir(shiftNames[shiftIdx].c_str());
233  shift_dir->cd();
234  }
235 
236  for (int varIdx = 0; varIdx < kNumVars; ++varIdx)
237  {
238  if (predIdx == 0)
239  predProp[predIdx][shiftIdx][varIdx]->SaveTo(shift_dir, ("nue_pred_"+histdef[varIdx].name).c_str());
240 
241  else
242  predMichel[predIdx][shiftIdx][varIdx]->SaveTo(syst_dir, ("nue_pred_"+histdef[varIdx].name).c_str());
243 
244  } // end of varIdx loop
245 
246  } // end of shiftIdx loop
247 
248  }// end of predIdx loop
249 
250 }
const XML_Char * name
Definition: expat.h:151
const int kNumVars
Definition: vars.h:14
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
set< int >::iterator it
const Var kCVNe
PID
Definition: Vars.cxx:35
const Var kPPFXFluxCVWgt
Definition: PPFXWeights.h:16
const Var kShwCalE
std::string name
Definition: NuePlotLists.h:12
const Cut kNue2017NDCVNSsb
Definition: NueCuts2017.h:302
const Binning kNue2017Binning
const Var kHadCalE([](const caf::SRProxy *sr){if(!sr->vtx.elastic.IsValid) return float(sr->slc.calE);if(sr->vtx.elastic.fuzzyk.nshwlid< 1) return float(sr->slc.calE);return((sr->slc.calE- sr->vtx.elastic.fuzzyk.png[0].shwlid.calE));})
Definition: NueVarsExtra.h:40
void make_nue_filesyst_pred(int systtype)
Uses MC for NC and CC components, assigns remainder of data to CC.
Definition: NumuDecomp.h:10
const Var kNueEnergy2017([](const caf::SRProxy *sr){if(!sr->vtx.elastic.IsValid) return-1.0;if(sr->vtx.elastic.fuzzyk.npng< 1) return-1.0;return NueRecoE_2017FDFit(kCVNemE(sr), kCVNhadE(sr));})
Definition: NueEnergy2017.h:11
const Cut kNumuND
Definition: NumuCuts.h:55
HistAxis axis
Definition: NuePlotLists.h:13
const SystShifts kNoShift
Definition: SystShifts.cxx:21
OStream cout
Definition: OStream.cxx:6
void NueExtrap(string beam="fhc", string cvntype="oldpresel")
Definition: NueExtrap.C:28
const Var kNue2017AnaBin([](const caf::SRProxy *sr){int selBin=kNue2017SelectionBin(sr);float nuE=kNueEnergy2017(sr);int nuEBin=nuE/0.5;assert(nuEBin<=8 &&"An event with nuE > 4.5 should never happen");int anaBin=9 *selBin+nuEBin;return anaBin;})
Definition: NueCuts2017.h:311
virtual void SaveTo(TDirectory *dir, const std::string &name) const override
Splits Data proportionally according to MC.
const SpillCut kStandardSpillCuts
Apply this unless you&#39;re doing something special.
Definition: SpillCuts.h:49
const HistAxis axisNumu
Take the output of an extrapolation and oscillate it as required.
const Var kXSecCVWgt2017
Definition: XsecTunes.h:37
Extrapolate each component using a separate ModularExtrapComponent.
Definition: ModularExtrap.h:23
static Binning Simple(int n, double lo, double hi, const std::vector< std::string > &labels={})
Definition: Binning.cxx:107
const Cut kNue2017FDAllSamples
Our FD selection including all samples, for making predictions, etc.
Definition: NueCuts2017.h:155
ECAFType
Definition: Loaders.h:19
enum BeamMode string