PredictionSystJoint2018.cxx
Go to the documentation of this file.
2 
6 #include "3FlavorAna/Systs/DummySystStorage.h" // norm systs
7 
9 #include "3FlavorAna/Cuts/NueCuts2018.h" // kNue2018Binning
10 
11 #include "OscLib/IOscCalc.h"
12 
13 #include "TObjString.h"
14 #include "TROOT.h"
15 
16 #include <iostream>
17 #include <unistd.h>
18 
19 namespace ana
20 {
21  REGISTER_LOADFROM("PredictionSystJoint2018", IPrediction, PredictionSystJoint2018);
22 
23  //----------------------------------------------------------------------
24  // Nue constructor
27  const std::vector<const ISyst*>& systs,
29  bool isFakeData,
30  bool mergePeripheral,
31  const std::string fname,
32  bool minimizeMemory)
33  : fBeam(beam)
34  {
35  gROOT->SetMustClean(false);
36 
37  fOscOrigin = osc->Copy();
38 
39  fSplitBySign = IsFHC();
40 
41  std::cout << "mergePeripheral is "
42  << (mergePeripheral ? "true" : "false") << std::endl;
43 
44  std::cout << "use FakeNDData predictions is "
45  << (isFakeData ? "true" : "false") << std::endl;
46  std::string extrapStr = "pred_";
47  if (extrap == kProportional) extrapStr += "xp_prop";
48  if (extrap == kCombo) extrapStr += "xp_combo";
49  if (extrap == kNoExtrap) extrapStr += "nxp";
50 
51  if (extrap == kFake) extrapStr += "xp_prop";
52 
53  extrapStr += "_";
54 
55  std::cout << "Loading nominals..." << std::endl;
56 
57  IPrediction *nom;
58 
59  std::string fnameSyst;
60 
61  if (fname != ""){
62  std::cout << "Loading predictions from unofficial file..." << std::endl;
63  std::cout << "\t--> " << fname << std::endl;
64  fnameSyst = fname;
65  }
66 
67  else if (IsFHC() && !isFakeData)
68  fnameSyst = "/nova/ana/nu_e_ana/Ana2018/Predictions/provisional_singles/NueProd4SystematicsOnRealNDData_2018-05-22/pred_allxp_Nue2018Axis_full_FHCAllSyst_nueconcat_realNDData_2018-05-22.root";
69  else if (!IsFHC() && !isFakeData)
70  fnameSyst = "/nova/ana/nu_e_ana/Ana2018/Predictions/provisional_singles/NueProd4SystematicsOnRealNDData_2018-05-22/pred_allxp_Nue2018Axis_full_RHCAllSyst_nueconcat_realNDData_2018-05-22.root";
71  else if (IsFHC() && isFakeData)
72  fnameSyst = "/nova/ana/nu_e_ana/Ana2018/Predictions/provisional_singles/NueProd4SystematicsOnFakeNDData_2018-05-23/pred_allxp_Nue2018Axis_full_FHCAllSysts_nueconcat_fakeNDData_05_22_2018_pass1_merged.root";
73  else if (!IsFHC() && isFakeData)
74  fnameSyst = "/nova/ana/nu_e_ana/Ana2018/Predictions/provisional_singles/NueProd4SystematicsOnFakeNDData_2018-05-23/pred_allxp_Nue2018Axis_full_RHCAllSysts_nueconcat_fakeNDData_05_22_2018_pass1_merged.root";
75 
76  TFile* fSyst = TFile::Open(fnameSyst.c_str(), "read");
77  if(fSyst->IsZombie()){
78  std::cout << "Couldn't open " << fnameSyst << std::endl;
79  abort();
80  }
81 
82  std::string varname ="/nue_pred_Nue2018Axis";
83 
84  // make PredictionExtendOwning
85 
86  IPrediction* nomCore = ana::LoadFrom<IPrediction>(fSyst, extrapStr + "Nominal/noShift" + varname).release();
87  PredictionNoExtrap* nomNoExtrap = ana::LoadFrom<PredictionNoExtrap>(fSyst, "pred_nxp_Nominal/noShift" + varname).release();
88 
89  fPredNom = std::make_unique<PredictExtendOwning>(nomCore, nomNoExtrap, kNue2018Binning, true);
90  nom = new PredictExtendOwning(nomCore, nomNoExtrap, kNue2018Binning, true);
91 
92 
93  // We have exactly 1 syst for the 2018 that is non-standard is some way
94  // The light level syst uses an alternative nominal. Grab that now.
95 
96  IPrediction* nomLightLevel;
97 
99  fSyst, extrapStr + "Lightlevel/noShift"+varname).release();
100  PredictionNoExtrap* nomLLNoExtrap = ana::LoadFrom<PredictionNoExtrap>(
101  fSyst, "pred_nxp_Lightlevel/noShift"+varname).release();
102 
103  nomLightLevel = new PredictExtendOwning(
104  nomLLCore, nomLLNoExtrap, true);
105 
106  for(const ISyst* syst: systs){
107 
108  if(syst == &kAna2018NormFHC ||
109  syst == &kAna2018NormRHC) continue;
110 
111  std::string systName = syst->ShortName();
112 
113  std::cout << "Loading " << systName << "..." << std::endl;
114 
115  ShiftedPreds sp;
116  sp.systName = systName;
117 
118  for(int sigma: GetDummyShifts(syst))
119  {
120  sp.shifts.push_back(sigma);
121 
122  if (sigma == 0)
123  {
124  if(syst == &kAnaLightlevelSyst ||
125  syst == &kAnaCherenkovSyst)
126  {
127  std::cout << "\tLooking at systName " << systName << ", so pushing back nomLightLevel " << std::endl;
128  sp.preds.push_back(nomLightLevel);
129  }
130  else
131  sp.preds.push_back(nom);
132  continue;
133  }
134 
135  const std::string sigmaStr = "/"+SigmaToString(sigma);
136 
137  IPrediction* shCore = ana::LoadFrom<IPrediction>(fSyst, extrapStr + systName + sigmaStr +varname).release();
138  PredictionNoExtrap* shNoExtrap = ana::LoadFrom<PredictionNoExtrap>(fSyst, "pred_nxp_" + systName + sigmaStr +varname).release();
139 
140  sp.preds.push_back(new PredictExtendOwning(shCore, shNoExtrap, kNue2018Binning, true));
141 
142  } // end for sigma
143  // std::cout<<" test find "<<systName<<" "<<syst->find(systName)<<std::endl;
144  fPreds.emplace(syst, sp);
145  } // end for syst
146 
147  // Calculate all the ratios etc
148  InitFits();
149 
150  AddNormSysts(systs);
151 
152  // Clean up. Don't use if you need to make debug plots
153  if(minimizeMemory){
154  MinimizeMemory();
155  }
156 
157  delete fSyst;
158  }
159 
160  //----------------------------------------------------------------------
161  // Numu constructor
164  const std::string beam,
165  int WhichQuant,
166  const std::vector<const ISyst*>& systs,
167  const std::string fname,
168  bool minimizeMemory)
169  : fBeam(beam)
170  {
171  gROOT->SetMustClean(false);
172 
173  fOscOrigin = osc->Copy();
174 
175  fSplitBySign = IsFHC();
176 
177  std::string extrapStr = "pred_";
178  if (extrap == kNuMu) extrapStr += "xp_numu";
179  if (extrap == kNuMuNoExtrap) extrapStr += "nxp";
180  extrapStr += "_";
181 
182  std::string VarName = "numu_pred_Energy_Quant"+std::to_string(WhichQuant);
183  if (WhichQuant == -1) VarName = "numu_pred_NumuEnergy";
184 
185  std::cout << "\nLooking at VarName " << VarName << " --> Loading nominals now..." << std::endl;
186 
187  TFile* fin = TFile::Open(fname.c_str(), "read");
188 
189  IPrediction *nom;
190  std::string NomVar = extrapStr+"Nominal/noShift/"+VarName;
192  nom = ana::LoadFrom<IPrediction>(fin, NomVar).release();
193  // --- The light level syst uses an alternative nominal. Grab that now.
194  std::string LNomVar = extrapStr+"Lightlevel/noShift/"+VarName;
195  IPrediction* nomLightLevel = ana::LoadFrom<IPrediction>(fin, LNomVar).release();
196 
197  for(const ISyst* syst: systs){
198 
199  if(syst == &kAna2018NormFHC ||
200  syst == &kAna2018NormRHC ||
201  syst->ShortName() == "RPAShapeRES2018Test" // || // As we aren't using it in Ana18, and so not in the Pred. files.
202  // -- When using Ana18 Predictions, need to comment out the following XSec Systs.
203  //syst->ShortName() == "RPAShapeenh2019" ||
204  //syst->ShortName() == "RPAShapesupp2019" ||
205  //syst->ShortName() == "RPAShapeRES2019" ||
206  //syst->ShortName() == "MECShape2018RPAFixNu" ||
207  //syst->ShortName() == "MECShape2018RPAFixAntiNu"
208  ) continue;
209 
210  std::string systName = syst->ShortName();
211 
212  std::cout << "Loading " << systName << "..." << std::endl;
213 
214  ShiftedPreds sp;
215  sp.systName = systName;
216 
217  for(int sigma: GetDummyShifts(syst)){
218  sp.shifts.push_back(sigma);
219 
220  if (sigma == 0){
221  if (systName == "Lightlevel" || systName == "Cherenkov") {
222  std::cout << "\tLooking at systName " << systName << ", so pushing back nomLightLevel " << std::endl;
223  sp.preds.push_back(nomLightLevel);
224  } else {
225  sp.preds.push_back(nom);
226  }
227  continue;
228  }
229 
230  const std::string sigmaStr = "/"+SigmaToString(sigma)+"/";
231 
232  auto shCore = ana::LoadFrom<PredictionExtrap> (fin, extrapStr + systName + sigmaStr +VarName).release();
233  sp.preds.push_back(shCore);
234 
235  } // end for sigma
236  fPreds.emplace(syst, sp);
237  } // end for syst
238 
239  // Calculate all the ratios etc
240  InitFits();
241 
242  AddNormSysts(systs);
243 
244  if(minimizeMemory){
245  MinimizeMemory();
246  }
247 
248  delete fin;
249  }
250 
251  //----------------------------------------------------------------------
253  {
254  }
255 
256  //----------------------------------------------------------------------
258  AddNormSysts(const std::vector<const ISyst*>& systs)
259  {
260  // DocDb-27230-v7 and DocDb-27821 >v1
261  // leave horn currents uncorrelated
262  // e.g. FHC norm applies shift to fhc,
263  // no shift to rhc, and vice versa
264  const double kNormFHC = 0.0144;
265  const double kNormRHC = 0.0064;
266 
267  // Include syst for wrong beam type to keep the fit happy, but set it to
268  // zero magnitude
269  if(IsFHC()){
270  AddNormSyst(&kAna2018NormFHC, systs, kNormFHC);
271  AddNormSyst(&kAna2018NormRHC, systs, 0);
272  }
273  else{
274  AddNormSyst(&kAna2018NormFHC, systs, 0);
275  AddNormSyst(&kAna2018NormRHC, systs, kNormRHC);
276  }
277  }
278 
279  //----------------------------------------------------------------------
281  AddNormSyst(const ISyst* syst,
282  const std::vector<const ISyst*>& systs,
283  double frac)
284  {
285  // Bail out if we're not actually requested to include this systematic
286  if(std::find(systs.begin(), systs.end(), syst) == systs.end()) return;
287 
288  TH1* h = Predict(fOscOrigin).ToTH1(1);
289  const int nBins = h->GetNbinsX();
290  delete h;
291 
292  ShiftedPreds normsp;
293  normsp.systName = syst->ShortName();
294  normsp.shifts = {0};
295  normsp.preds = {fPredNom.get()};
296  normsp.fits.resize(kNCoeffTypes);
297  if(fSplitBySign) normsp.fitsNubar.resize(kNCoeffTypes);
298 
299  for(int coeff = 0; coeff < kNCoeffTypes; ++coeff){
300  normsp.fits[coeff].resize(nBins+2);
301  if(fSplitBySign) normsp.fitsNubar[coeff].resize(nBins+2);
302  for(int binIdx = 0; binIdx < nBins+2; ++binIdx){
303  normsp.fits[coeff][binIdx].emplace_back(0, 0, frac, 1);
304  if(fSplitBySign) normsp.fitsNubar[coeff][binIdx].emplace_back(0, 0, frac, 1);
305  }
306  }
307 
308  normsp.FillRemaps();
309 
310  normsp.nCoeffs = 1;
311 
312  fPreds.emplace(syst, normsp);
313  }
314 
315  //----------------------------------------------------------------------
316  std::unique_ptr<PredictionSystJoint2018> PredictionSystJoint2018::LoadFrom(TDirectory* dir, const std::string& name)
317  {
318  dir = dir->GetDirectory(name.c_str()); // switch to subdir
319  assert(dir);
320 
321  TObjString* tag = (TObjString*)dir->Get("type");
322  assert(tag);
323 
324  assert(tag->GetString() == "PredictionSystJoint2018");
325  delete tag;
326 
327  std::unique_ptr<PredictionSystJoint2018> ret(new PredictionSystJoint2018);
328  LoadFromBody(dir, ret.get(), {&kAna2018NormFHC, &kAna2018NormRHC});
329 
330  const TObjString* beam = (TObjString*)dir->Get("beam");
331  assert(beam);
332  ret->fBeam = beam->GetString();
333 
334  ret->fSplitBySign = ret->IsFHC();
335 
336  // Just enable all norm systs unconditionally. Way too hard to do otherwise
337  // NB fBeam was set above. That's necessary for correctness.
338  ret->AddNormSysts({&kAna2018NormFHC, &kAna2018NormRHC});
339 
340  delete dir;
341 
342  return ret;
343  }
344 
345  //----------------------------------------------------------------------
346  void PredictionSystJoint2018::SaveTo(TDirectory* dir, const std::string& name) const
347  {
348  PredictionInterp::SaveTo(dir, name);
349 
350  TDirectory* tmp = gDirectory;
351 
352  dir = (TDirectory*)dir->GetDirectory(name.c_str());
353  dir->cd();
354  // Overwrite "IPrediction"
355  TObjString("PredictionSystJoint2018").Write("type");
356 
357  TObjString(fBeam.c_str()).Write("beam");
358 
359  dir->Write();
360  delete dir;
361 
362  tmp->cd();
363  }
364 
365  //----------------------------------------------------------------------
367  {
368  if(sigma == -2) return "minusTwo";
369  if(sigma == -1) return "minusOne";
370  if(sigma == 0) return "";
371  if(sigma == +1) return "plusOne";
372  if(sigma == +2) return "plusTwo";
373 
374  std::cout << "PredictionSystJoint2018: unknown sigma " << sigma << std::endl;
375  abort();
376  }
377 
378  //----------------------------------------------------------------------
380  {
381  if(fBeam == "fhc" ) return true;
382  if(fBeam == "rhc") return false;
383 
384  std::cout << "PredictionSystJoint2018: Bad beam string: "
385  << fBeam << std::endl;
386  abort();
387  }
388 
389 
390 } // namespace
TString fin
Definition: Style.C:24
const XML_Char * name
Definition: expat.h:151
void AddNormSysts(const std::vector< const ISyst * > &systs)
Must set fBeam first.
int nBins
Definition: plotROC.py:16
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
std::vector< SystGroupDef > systs
Definition: syst_header.h:385
bool mergePeripheral
std::string systName
What systematic we&#39;re interpolating.
void SaveTo(TDirectory *dir, const std::string &name) const override
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
virtual const std::string & ShortName() const final
The name printed out to the screen.
Definition: ISyst.h:27
Loads shifted spectra from files.
General interface to oscillation calculators.
Definition: StanTypedefs.h:22
std::vector< std::vector< std::vector< Coeffs > > > fitsNubar
Will be filled if signs are separated, otherwise not.
std::vector< int > GetDummyShifts(const ISyst *s)
Float_t tmp
Definition: plot.C:36
Encapsulate code to systematically shift a caf::SRProxy.
Definition: ISyst.h:14
std::unique_ptr< IPrediction > fPredNom
The nominal prediction.
const DummyAnaSyst kAnaLightlevelSyst("Lightlevel","Lightlevel")
static std::unique_ptr< PredictionSystJoint2018 > LoadFrom(TDirectory *dir, const std::string &name)
virtual void SaveTo(TDirectory *dir, const std::string &name) const override
const DummyAnaSyst kAnaCherenkovSyst("Cherenkov","Cherenkov")
std::vector< double > shifts
Shift values sampled.
std::vector< IPrediction * > preds
double frac(double x)
Fractional part.
const Binning kNue2018Binning
Definition: NueCuts2018.cxx:92
std::unique_ptr< IPrediction > LoadFrom< IPrediction >(TDirectory *dir, const std::string &label)
Definition: IPrediction.cxx:18
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
virtual _IOscCalc * Copy() const =0
Spectrum Predict(osc::IOscCalc *calc) const override
bool IsFHC() const
Interprets fBeam.
const DummyAnaSyst kAna2018NormRHC("NormRHC2018","RHC. Norm.")
TDirectory * dir
Definition: macro.C:5
REGISTER_LOADFROM("BENDecomp", IDecomp, BENDecomp)
assert(nhit_max >=nhit_nbins)
Standard interface to all prediction techniques.
Definition: IPrediction.h:57
std::string SigmaToString(int sigma) const
Prediction that just uses FD MC, with no extrapolation.
std::string to_string(ModuleType mt)
Definition: ModuleType.h:32
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={})
void AddNormSyst(const ISyst *syst, const std::vector< const ISyst * > &systs, double frac)
const DummyAnaSyst kAna2018NormFHC("NormFHC2018","FHC. Norm.")
std::vector< std::vector< std::vector< Coeffs > > > fits
Indices: [type][histogram bin][shift bin].
gm Write()
enum BeamMode string