PredictionSystNue2017.cxx
Go to the documentation of this file.
5 
8 
9 #include "OscLib/IOscCalc.h"
10 
11 #include "TObjString.h"
12 
13 #include <iostream>
14 
15 namespace ana
16 {
17 
18  REGISTER_LOADFROM("PredictionSystNue2017", IPrediction, PredictionSystNue2017);
19 
20  //----------------------------------------------------------------------
21 
22  const DummyNue2017Syst kNue2017MECq0ShapeSyst ("nue2017MECq0Shape",
23  "MECq0Shape",
24  {-2, -1, 0, +1, +2}, false);
25  const DummyNue2017Syst kNue2017MaCCQE_reducedSyst ("nue2017MaCCQE_reduced",
26  "MaCCQE_reduced",
27  {-2, -1, 0, +1, +2}, false);
28  const DummyNue2017Syst kNue2017RPAShapeRESSyst ("nue2017RPAShapeRES",
29  "RPAShapeRES",
30  {0, +1, +2}, false);
31  const DummyNue2017Syst kNue2017RPAShapeEnhSyst ("nue2017RPAShapeenh",
32  "RPAShapeEnh",
33  {-2, -1, 0, +1, +2}, false);
34  const DummyNue2017Syst kNue2017MaCCRESSyst ("nue2017MaCCRES",
35  "MaCCRES",
36  {-2, -1, 0, +1, +2}, false);
37  const DummyNue2017Syst kNue2017MaNCRESSyst ("nue2017MaNCRES",
38  "MaNCRES",
39  {-2, -1, 0, +1, +2}, false);
40  const DummyNue2017Syst kNue2017MvCCRESSyst ("nue2017MvCCRES",
41  "MvCCRES",
42  {-2, -1, 0, +1, +2}, false);
43  const DummyNue2017Syst kNue2017RadCorrNueSyst ("nue2017radcorrnue",
44  "RadCorrNue",
45  {-2, -1, 0, +1, +2}, false);
46  const DummyNue2017Syst kNue2017RadCorrNueBarSyst ("nue2017radcorrnuebar",
47  "RadCorrNueBar",
48  {-2, -1, 0, +1, +2}, false);
49  const DummyNue2017Syst kNue2017SecondClassCurrSyst ("nue20172ndclasscurr",
50  "SecondClassCurr",
51  {-2, -1, 0, +1, +2}, false);
52  const DummyNue2017Syst kNue2017DISvnCC1piSyst ("nue2017DISvnCC1pi",
53  "DISvnCC1pi",
54  {-2, -1, 0, +1, +2}, false);
55  const DummyNue2017Syst kNue2017AllSmallXSecSyst ("nue2017SumSmallXSecJoint2017",
56  "SumSmallXSecJoint2017",
57  {-2, -1, 0, +1, +2}, false);
58 
59 
60  const DummyNue2017Syst kNue2017AbsCalibSyst ("nue2017Calibration",
61  "AbsCalib",
62  {-1, 0, +1}, false);
63  const DummyNue2017Syst kNue2017RelCalibSyst ("nue2017RelativeCalib",
64  "RelCalib",
65  {-2, 0, +2}, false);
66  const DummyNue2017Syst kNue2017CalibShapeSyst ("nue2017CalibShape",
67  "CalibShape",
68  {0, +1}, false);
69  const DummyNue2017Syst kNue2017LightLevelSyst ("nue2017Lightlevel",
70  "LightLevel",
71  {-1, 0, +1}, false);
72  const DummyNue2017Syst kNue2017CherenkovSyst ("nue2017Cherenkov",
73  "Cherenkov",
74  {0, +1}, false);
75 
76  const DummyNue2017Syst kNue2017PPFXPC00Syst ("nue2017ppfx_pc00",
77  "PPFXPC00Sys",
78  {-2, -1, 0, +1, +2}, false);
79  const DummyNue2017Syst kNue2017PPFXPC01Syst ("nue2017ppfx_pc01",
80  "PPFXPC01",
81  {-2, -1, 0, +1, +2}, false);
82  const DummyNue2017Syst kNue2017PPFXPC02Syst ("nue2017ppfx_pc02",
83  "PPFXPC02",
84  {-2, -1, 0, +1, +2}, false);
85  const DummyNue2017Syst kNue2017PPFXPC03Syst ("nue2017ppfx_pc03",
86  "PPFXPC03",
87  {-2, -1, 0, +1, +2}, false);
88  const DummyNue2017Syst kNue2017PPFXPC04Syst ("nue2017ppfx_pc04",
89  "PPFXPC04",
90  {-2, -1, 0, +1, +2}, false);
91  const DummyNue2017Syst kNue2017BeamTransportSyst ("nue2017beamTransportComb",
92  "Combined Beam Transport",
93  {-2, -1, 0, +1, +2}, false);
94 
95  const DummyNue2017Syst kNue2017ExtrapSigSyst ("nue2017extrap_signalkin",
96  "Nue Signal Extrapolation",
97  {-2, -1, 0, +1, +2}, false);
98  const DummyNue2017Syst kNue2017ExtrapBkgSyst ("nue2017extrap_bkg",
99  "Nue Bkg Extrapolation",
100  {-2, -1, 0, +1, +2}, false);
101 
102  const DummyNue2017Syst kNue2017NormSystSig ("nue2017NueNormSig",
103  "Sig. Norm.",
104  {0, +1}, false);
105  const DummyNue2017Syst kNue2017NormSystBkg ("nue2017NueNormBkg",
106  "Bkg. Norm.",
107  {0, +1}, false);
108  const DummyNue2017Syst kNue2017NormSystBoth("nue2017NormBoth",
109  "Sig+Bkg Norm.",
110  {0, +1}, false);
111 
112 
113  //----------------------------------------------------------------------
116  const std::vector<const DummyNue2017Syst*>& systs,
117  const std::string fname)
118  : PredictionInterp()
119  {
120  fOscOrigin = osc->Copy();
121 
122  std::string extrapStr = "pred_";
123  if (extrap == ENue2017ExtrapType::kProportional) extrapStr += "xp_prop";
124  if (extrap == ENue2017ExtrapType::kCombo) extrapStr += "xp_combo";
125  if (extrap == ENue2017ExtrapType::kNoExtrap) extrapStr += "nxp";
126  extrapStr += "_";
127 
128  std::cout << "Loading nominals..." << std::endl;
129 
130  IPrediction *nom;
131 
132  //std::string varname ="";
133  std::string varname ="/nue_pred_EnergyCVN2D";
134  if (extrap == ENue2017ExtrapType::kNoExtrap){
135  fPredNom = LoadFromFile<IPrediction>
136  (fname, extrapStr + "Nominal/noShift" + varname);
137  nom = LoadFromFile<IPrediction>
138  (fname, extrapStr + "Nominal/noShift" + varname).release();
139  }
140  else{
141  // We need to pull out a PredictionExtrap and a PredictionNoExtrap
142  // to make the final PredictionExtendToPeripheral
143  PredictionExtrap* nomCore = LoadFromFile<PredictionExtrap>
144  (fname, extrapStr + "Nominal/noShift" + varname).release();
145  PredictionNoExtrap* nomNoExtrap = LoadFromFile<PredictionNoExtrap>
146  (fname, "pred_nxp_Nominal/noShift" + varname).release();
147  fPredNom = std::unique_ptr<PredictionExtendToPeripheral>
148  (new PredictionExtendToPeripheral(nomCore,nomNoExtrap,true));
149 
150  nom = new PredictionExtendToPeripheral(nomCore,nomNoExtrap,true);
151 
152  }
153 
154  // We have exactly 1 syst for the 2017 that is non-standard is some way
155  // The light level syst uses an alternative nominal. Grab that now.
156  IPrediction* nomLightLevel;
157 
158  if (extrap == ENue2017ExtrapType::kNoExtrap){
159  nomLightLevel = LoadFromFile<IPrediction>(fname, extrapStr + "Lightlevel/noShift"+
160  varname).release();
161  }
162  else{
163  PredictionExtrap* nomLLCore = LoadFromFile<PredictionExtrap>
164  (fname, extrapStr + "Lightlevel/noShift"+varname).release();
165  PredictionNoExtrap* nomLLNoExtrap = LoadFromFile<PredictionNoExtrap>
166  (fname,"pred_nxp_Nominal/noShift"+varname).release();
167  nomLightLevel = new PredictionExtendToPeripheral(nomLLCore,nomLLNoExtrap,true);
168  }
169 
170  for(const DummyNue2017Syst* syst: systs){
171  if(syst == &kNue2017NormSystBoth ||
172  syst == &kNue2017NormSystSig ||
173  syst == &kNue2017NormSystBkg) continue;
174 
175  // Hack. For sake of distinct names, every syst's fist 7 characters
176  // are nue2017
177  std::string systName = syst->ShortName();
178  systName.erase(0,7);
179 
180  std::cout << "Loading " << systName << "..." << std::endl;
181 
182  ShiftedPreds sp;
183  sp.systName = systName;
184 
185  for(int sigma: syst->Shifts()){
186  sp.shifts.push_back(sigma);
187 
188  if (sigma == 0){
189  if (systName == "Lightlevel" || systName == "Cherenkov") {
190  std::cout << "\tLooking at systName " << systName
191  << ", so pushing back nomLightLevel " << std::endl;
192  sp.preds.push_back(nomLightLevel);
193  } else {
194  sp.preds.push_back(nom);
195  }
196  continue;
197  }
198 
199  std::string sigmaStr = "/";
200  if(sigma == -2) sigmaStr += "minusTwo";
201  if(sigma == -1) sigmaStr += "minusOne";
202  if(sigma == +1) sigmaStr += "plusOne";
203  if(sigma == +2) sigmaStr += "plusTwo";
204 
205  //std::string varname ="";
206  std::string varname ="/nue_pred_EnergyCVN2D";
207  auto shCore = LoadFromFile<PredictionExtrap>
208  (fname, extrapStr + systName + sigmaStr +varname).release();
209  auto shNoExtrap = LoadFromFile<PredictionNoExtrap>
210  (fname, "pred_nxp_" + systName + sigmaStr +varname).release();
211  sp.preds.push_back(new PredictionExtendToPeripheral(shCore,shNoExtrap,true));
212  } // end for sigma
213 
214  fPreds.emplace(syst, sp);
215  } // end for syst
216 
217  // Calculate all the ratios etc
218  InitFits();
219 
220  AddNormSysts(systs);
221  }
222 
223  //----------------------------------------------------------------------
225  {
226  }
227 
228  //----------------------------------------------------------------------
229  void PredictionSystNue2017::AddNormSysts(const std::vector<const DummyNue2017Syst*>& systs)
230  {
231  TH1* h = Predict(fOscOrigin).ToTH1(1);
232  const int nBins = h->GetNbinsX();
233  delete h;
234 
235  // Got this 3.5% from Alex. TODO add reference
236  // leaving Sig and bkg hooks in place, but they are not used in the fit
237  const double normBothSyst = 3.5/100;
238  const double normBkgSyst = 3.5/100;
239  const double normSigSyst = 3.5/100;
240 
241  // TODO find a way to write these three norm systs as a loop
242  // if(std::count(systs.begin(), systs.end(), &kNue2017NormSystBoth)){
243  if(std::count_if(systs.begin(), systs.end(),
244  [](const ISyst* s){return s->ShortName() == "nue2017NormBoth";})){
245  ShiftedPreds normsp;
246  normsp.systName = "NormBoth";
247  normsp.shifts = {0};
248  normsp.preds = {fPredNom.get()};
249  for(int coeff = 0; coeff < kNCoeffTypes; ++coeff){
250  normsp.fits[coeff].resize(nBins+2);
251  for(int binIdx = 0; binIdx < nBins+2; ++binIdx){
252  normsp.fits[coeff][binIdx].push_back(Coeffs(0, 0, normBothSyst, 1));
253  }
254  }
255  normsp.FillRemaps();
256  fPreds.emplace(&kNue2017NormSystBoth, normsp);
257  }
258 
259 
260  // if(std::count(systs.begin(), systs.end(), &kNue2017NormSystBkg)){
261  if(std::count_if(systs.begin(), systs.end(),
262  [](const ISyst* s){return s->ShortName() == "nue2017NueNormBkg";})){
263  ShiftedPreds normbksp;
264  normbksp.systName = "NueNormBkg";
265  normbksp.shifts = {0};
266  normbksp.preds = {fPredNom.get()};
267  for(int coeff = 0; coeff < kNCoeffTypes; ++coeff){
268  normbksp.fits[coeff].resize(nBins+2);
269  for(int binIdx = 0; binIdx < nBins+2; ++binIdx){
270  normbksp.fits[coeff][binIdx].push_back(Coeffs(0, 0, (coeff != kNueApp) ? normBkgSyst : 0, 1));
271  }
272  }
273  normbksp.FillRemaps();
274  fPreds.emplace(&kNue2017NormSystBkg, normbksp);
275  }
276 
277 
278  // if(std::count(systs.begin(), systs.end(), &kNue2017NormSystSig)){
279  if(std::count_if(systs.begin(), systs.end(),
280  [](const ISyst* s){return s->ShortName() == "nue2017NueNormSig";})){
281  ShiftedPreds normsigsp;
282  normsigsp.systName = "NueNormSig";
283  normsigsp.shifts ={0};
284  normsigsp.preds = {fPredNom.get()};
285  for(int coeff = 0; coeff < kNCoeffTypes; ++coeff){
286  normsigsp.fits[coeff].resize(nBins+2);
287  for(int binIdx = 0; binIdx < nBins+2; ++binIdx){
288  normsigsp.fits[coeff][binIdx].push_back(Coeffs(0, 0, (coeff == kNueApp) ? normSigSyst : 0, 1));
289  }
290  }
291  normsigsp.FillRemaps();
292  fPreds.emplace(&kNue2017NormSystSig, normsigsp);
293  }
294  }
295 
296  //----------------------------------------------------------------------
297  std::unique_ptr<PredictionSystNue2017> PredictionSystNue2017::LoadFrom(TDirectory* dir, const std::string& name)
298  {
299  dir = dir->GetDirectory(name.c_str()); // switch to subdir
300  assert(dir);
301 
302  TObjString* tag = (TObjString*)dir->Get("type");
303  assert(tag);
304 
305  assert(tag->GetString() == "PredictionSystNue2017");
306  delete tag;
307 
308  std::unique_ptr<PredictionSystNue2017> ret(new PredictionSystNue2017);
310 
311  // Just enable all norm systs unconditionally. Way too hard to do otherwise
313 
314  delete dir;
315 
316  return ret;
317  }
318 
319  //----------------------------------------------------------------------
320  void PredictionSystNue2017::SaveTo(TDirectory* dir, const std::string& name) const
321  {
322  PredictionInterp::SaveTo(dir, name);
323 
324  TDirectory* tmp = gDirectory;
325 
326  dir = dir->GetDirectory(name.c_str());
327  dir->cd();
328  // Overwrite "IPrediction"
329  TObjString("PredictionSystNue2017").Write((name+"/type").c_str());
330 
331  dir->Write();
332  delete dir;
333 
334  tmp->cd();
335  }
336 
337 } // namespace
const DummyNue2017Syst kNue2017NormSystSig("nue2017NueNormSig","Sig. Norm.",{0,+1}, false)
const XML_Char * name
Definition: expat.h:151
const DummyNue2017Syst kNue2017CherenkovSyst("nue2017Cherenkov","Cherenkov",{0,+1}, false)
int nBins
Definition: plotROC.py:16
Implements systematic errors by interpolation between shifted templates.
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
const DummyNue2017Syst kNue2017PPFXPC04Syst("nue2017ppfx_pc04","PPFXPC04",{-2,-1, 0,+1,+2}, false)
const DummyNue2017Syst kNue2017NormSystBkg("nue2017NueNormBkg","Bkg. Norm.",{0,+1}, false)
std::vector< SystGroupDef > systs
Definition: syst_header.h:385
std::string systName
What systematic we&#39;re interpolating.
TH1D * ToTH1(double exposure, Color_t col=kBlack, Style_t style=kSolid, EExposureType expotype=kPOT, EBinType bintype=kBinContent) const
Histogram made from this Spectrum, scaled to some exposure.
Definition: Spectrum.cxx:148
void AddNormSysts(const std::vector< const DummyNue2017Syst * > &systs)
const DummyNue2017Syst kNue2017MaNCRESSyst("nue2017MaNCRES","MaNCRES",{-2,-1, 0,+1,+2}, false)
const DummyNue2017Syst kNue2017ExtrapSigSyst("nue2017extrap_signalkin","Nue Signal Extrapolation",{-2,-1, 0,+1,+2}, false)
const DummyNue2017Syst kNue2017MaCCRESSyst("nue2017MaCCRES","MaCCRES",{-2,-1, 0,+1,+2}, false)
const DummyNue2017Syst kNue2017ExtrapBkgSyst("nue2017extrap_bkg","Nue Bkg Extrapolation",{-2,-1, 0,+1,+2}, false)
General interface to oscillation calculators.
Definition: StanTypedefs.h:22
const DummyNue2017Syst kNue2017BeamTransportSyst("nue2017beamTransportComb","Combined Beam Transport",{-2,-1, 0,+1,+2}, false)
const DummyNue2017Syst kNue2017MvCCRESSyst("nue2017MvCCRES","MvCCRES",{-2,-1, 0,+1,+2}, false)
Float_t tmp
Definition: plot.C:36
const DummyNue2017Syst kNue2017PPFXPC02Syst("nue2017ppfx_pc02","PPFXPC02",{-2,-1, 0,+1,+2}, false)
const DummyNue2017Syst kNue2017CalibShapeSyst("nue2017CalibShape","CalibShape",{0,+1}, false)
Encapsulate code to systematically shift a caf::SRProxy.
Definition: ISyst.h:14
const DummyNue2017Syst kNue2017NormSystBoth("nue2017NormBoth","Sig+Bkg Norm.",{0,+1}, false)
const DummyNue2017Syst kNue2017RPAShapeRESSyst("nue2017RPAShapeRES","RPAShapeRES",{0,+1,+2}, false)
std::unique_ptr< IPrediction > fPredNom
The nominal prediction.
const DummyNue2017Syst kNue2017RadCorrNueBarSyst("nue2017radcorrnuebar","RadCorrNueBar",{-2,-1, 0,+1,+2}, false)
const XML_Char * s
Definition: expat.h:262
Loads shifted spectra from files.
virtual void SaveTo(TDirectory *dir, const std::string &name) const override
const DummyNue2017Syst kNue2017LightLevelSyst("nue2017Lightlevel","LightLevel",{-1, 0,+1}, false)
std::vector< double > shifts
Shift values sampled.
const DummyNue2017Syst kNue2017DISvnCC1piSyst("nue2017DISvnCC1pi","DISvnCC1pi",{-2,-1, 0,+1,+2}, false)
const DummyNue2017Syst kNue2017PPFXPC00Syst("nue2017ppfx_pc00","PPFXPC00Sys",{-2,-1, 0,+1,+2}, false)
const DummyNue2017Syst kNue2017AbsCalibSyst("nue2017Calibration","AbsCalib",{-1, 0,+1}, false)
const DummyNue2017Syst kNue2017MECq0ShapeSyst("nue2017MECq0Shape","MECq0Shape",{-2,-1, 0,+1,+2}, false)
const DummyNue2017Syst kNue2017PPFXPC01Syst("nue2017ppfx_pc01","PPFXPC01",{-2,-1, 0,+1,+2}, false)
std::vector< IPrediction * > preds
const DummyNue2017Syst kNue2017RelCalibSyst("nue2017RelativeCalib","RelCalib",{-2, 0,+2}, false)
double coeff(double W, int m, int l, bool real)
Oscillation probability calculators.
Definition: Calcs.h:5
double sigma(TH1F *hist, double percentile)
OStream cout
Definition: OStream.cxx:6
const DummyNue2017Syst kNue2017SecondClassCurrSyst("nue20172ndclasscurr","SecondClassCurr",{-2,-1, 0,+1,+2}, false)
virtual _IOscCalc * Copy() const =0
const DummyNue2017Syst kNue2017RadCorrNueSyst("nue2017radcorrnue","RadCorrNue",{-2,-1, 0,+1,+2}, false)
Spectrum Predict(osc::IOscCalc *calc) const override
static std::unique_ptr< PredictionSystNue2017 > LoadFrom(TDirectory *dir, const std::string &name)
TDirectory * dir
Definition: macro.C:5
REGISTER_LOADFROM("BENDecomp", IDecomp, BENDecomp)
void SaveTo(TDirectory *dir, const std::string &name) const override
assert(nhit_max >=nhit_nbins)
const DummyNue2017Syst kNue2017MaCCQE_reducedSyst("nue2017MaCCQE_reduced","MaCCQE_reduced",{-2,-1, 0,+1,+2}, false)
Standard interface to all prediction techniques.
Definition: IPrediction.h:57
const DummyNue2017Syst kNue2017PPFXPC03Syst("nue2017ppfx_pc03","PPFXPC03",{-2,-1, 0,+1,+2}, false)
Prediction that just uses FD MC, with no extrapolation.
Take the output of an extrapolation and oscillate it as required.
std::unordered_map< const ISyst *, ShiftedPreds > fPreds
osc::IOscCalc * fOscOrigin
The oscillation values we assume when evaluating the coefficients.
static void LoadFromBody(TDirectory *dir, PredictionInterp *ret, std::vector< const ISyst * > veto={})
const DummyNue2017Syst kNue2017AllSmallXSecSyst("nue2017SumSmallXSecJoint2017","SumSmallXSecJoint2017",{-2,-1, 0,+1,+2}, false)
const DummyNue2017Syst kNue2017RPAShapeEnhSyst("nue2017RPAShapeenh","RPAShapeEnh",{-2,-1, 0,+1,+2}, false)
std::vector< std::vector< std::vector< Coeffs > > > fits
Indices: [type][histogram bin][shift bin].
enum BeamMode string