NusSystMaker.C
Go to the documentation of this file.
1 // Based on results from individual systematic studies,
2 // create a file containing histograms of shift/nominal ratios.
3 // These histograms are used by NusSystFromHist
4 // to create Nus group ISyst objects.
5 // There are three types of syst groups:
6 // Errors added in quadrature -- these have a std::vector of labels
7 // On/off errors -- these have a single label
8 // An up/down systematic -- these have a std::pair of labels
9 
13 
14 using namespace ana;
15 
16 // Fill a second map with clones of histograms in the first map
17 void CloneMap(std::map<std::string, std::map<std::string, TH1*> >& noms,
18  std::map<std::string, std::map<std::string, TH1*> >& clone);
19 
20 // Open the file systfile
21 // Fill the decomp, predNE, predSt maps from the file
22 // Put the nominal spectra from the maps into noms
23 void FillMaps(std::string folder, std::string systfile, std::string sample,
24  std::map<std::string, std::map<std::string, TH1*> >& noms,
25  std::map<std::string, IDecomp*>& decomps_nominal,
26  std::map<std::string, std::map<std::string, std::map<int, IDecomp*> > >& decomps_shifted,
27  std::map<std::string, PredictionNoExtrap*>& predsNE_nominal,
28  std::map<std::string, std::map<std::string, std::map<int, PredictionNoExtrap*> > >& predsNE_shifted,
29  std::map<std::string, PredictionSterile*>& predsSt_nominal,
30  std::map<std::string, std::map<std::string, std::map<int, PredictionSterile*> > >& predsSt_shifted);
31 
32 // Make a standardized name for the syst ratio histogram
35 // Make a standardized title for the syst ratio histogram
38 
39 // Make nominal ratio histograms, a line at 1 corresponding to 0 sigma
40 void NominalRat(std::map<std::string, std::map<std::string, TH1*> >& noms,
41  std::map<std::string, std::map<std::string, TH1*> >& line);
42 
43 // Add errors in quadrature for all systs in the shiftlabels vector
44 // Set the errors for the hists in noms to these errors
45 void QuadErrors(const std::vector<std::string>& shiftlabels, std::string sample,
46  std::map<std::string, std::map<std::string, TH1*> >& noms,
47  std::map<std::string, std::map<std::string, std::map<int, IDecomp*> > >& decomps_shifted,
48  std::map<std::string, std::map<std::string, std::map<int, PredictionNoExtrap*> > >& predsNE_shifted,
49  std::map<std::string, std::map<std::string, std::map<int, PredictionSterile*> > >& predsSt_shifted);
50 
51 // Save the histograms in hists to out
52 void SaveMaps(std::map<std::string, std::map<std::string, TH1*> >& hists,
53  TDirectory* out, std::string syst, std::string sigma);
54 
55 // Set bins in shft to shift/nominal using the shiftlabel syst
56 void SetBinsOnOffShift(std::string shiftlabel, std::string sample,
57  std::map<std::string, std::map<std::string, TH1*> >& noms,
58  std::map<std::string, std::map<std::string, TH1*> >& shft,
59  std::map<std::string, std::map<std::string, std::map<int, IDecomp*> > >& decomps_shifted,
60  std::map<std::string, std::map<std::string, std::map<int, PredictionNoExtrap*> > >& predsNE_shifted,
61  std::map<std::string, std::map<std::string, std::map<int, PredictionSterile*> > >& predsSt_shifted);
62 
63 // Perform the on/off syst type method for a syst that has an up and down
64 // These are shifts that are not added in quadrature
65 // The +1 sigma shift should be the first parameter in shiftlabels
66 void SetBinsPlusMinusOne(std::pair<std::string, std::string> shiftlabels, std::string sample,
67  std::map<std::string, std::map<std::string, TH1*> >& noms,
68  std::map<std::string, std::map<std::string, TH1*> >& sfp1,
69  std::map<std::string, std::map<std::string, TH1*> >& sfm1,
70  std::map<std::string, std::map<std::string, std::map<int, IDecomp*> > >& decomps_shifted,
71  std::map<std::string, std::map<std::string, std::map<int, PredictionNoExtrap*> > >& predsNE_shifted,
72  std::map<std::string, std::map<std::string, std::map<int, PredictionSterile*> > >& predsSt_shifted);
73 
74 // Set bins in sfp1 and sfm1 to (nominal +/- error)/nominal
75 // Errors should have been calculated by QuadErrors
76 void SetBinsPlusMinusOne(std::map<std::string, std::map<std::string, TH1*> >& noms,
77  std::map<std::string, std::map<std::string, TH1*> >& sfp1,
78  std::map<std::string, std::map<std::string, TH1*> >& sfm1);
79 
80 //------------------------------------------------------------------------------
82 {
83  TH1::AddDirectory(0);
84 
85  // Choose token for systematic samples
86  std::string sample = "Ana01";
87 
88  // Lists of labels to use from each syst group
89  std::vector<std::string> labelsBeam;
90  labelsBeam.push_back("BeamHornCurrent");
91  labelsBeam.push_back("BeamSpotSize");
92  labelsBeam.push_back("BeamPosX");
93  labelsBeam.push_back("BeamPosY");
94  labelsBeam.push_back("BeamH1PosXY");
95  labelsBeam.push_back("BeamH2PosXY");
96  labelsBeam.push_back("BeamTarget");
97  labelsBeam.push_back("BeamExp");
98  labelsBeam.push_back("BeamGeomWater");
99  labelsBeam.push_back("BeamFlukaVersion");
100  labelsBeam.push_back("BeamG4");
101  labelsBeam.push_back("BeamNA49");
102 
103  std::string labelBirks = "BirksC";
104 
105  std::pair<std::string, std::string> labelsFlat =
106  std::make_pair("CalFlat105", "CalFlat095");
107  std::string labelSlopeX = "CalSlopeX";
108  std::string labelSlopeY = "CalSlopeY";
109 
110  std::pair<std::string, std::string> labelsCalibRel =
111  std::make_pair("CalFlat105ND", "CalFlat095ND");
112 
113  std::string labelNumu = "Numu";
114  std::pair<std::string, std::string> labelsNue =
115  std::make_pair("NueUp", "NueDn");
116 
117  // Commented systs are >1% effect on NC signal (None!)
118  std::vector<std::string> labelsGENIE;
119  labelsGENIE.push_back("ReweightMaNCEL");
120  labelsGENIE.push_back("ReweightEtaNCEL");
121  labelsGENIE.push_back("ReweightMaCCQEshape");
122  labelsGENIE.push_back("ReweightNormCCQE");
123  labelsGENIE.push_back("ReweightCCQEPauliSupViaKF");
124  labelsGENIE.push_back("ReweightVecCCQEshape");
125  labelsGENIE.push_back("ReweightCCQEMomDistroFGtoSF");
126  labelsGENIE.push_back("ReweightMaCCRES");
127  labelsGENIE.push_back("ReweightMvCCRES");
128  labelsGENIE.push_back("ReweightMaNCRES");
129  labelsGENIE.push_back("ReweightMvNCRES");
130  labelsGENIE.push_back("ReweightMaCOHpi");
131  labelsGENIE.push_back("ReweightR0COHpi");
132  labelsGENIE.push_back("ReweightRvpCC1pi");
133  labelsGENIE.push_back("ReweightRvpCC2pi");
134  labelsGENIE.push_back("ReweightRvnCC1pi");
135  labelsGENIE.push_back("ReweightRvnCC2pi");
136  labelsGENIE.push_back("ReweightRvpNC1pi");
137  labelsGENIE.push_back("ReweightRvpNC2pi");
138  labelsGENIE.push_back("ReweightRvnNC1pi");
139  labelsGENIE.push_back("ReweightRvnNC2pi");
140  labelsGENIE.push_back("ReweightRvbarpCC1pi");
141  labelsGENIE.push_back("ReweightRvbarpCC2pi");
142  labelsGENIE.push_back("ReweightRvbarnCC1pi");
143  labelsGENIE.push_back("ReweightRvbarnCC2pi");
144  labelsGENIE.push_back("ReweightRvbarpNC1pi");
145  labelsGENIE.push_back("ReweightRvbarpNC2pi");
146  labelsGENIE.push_back("ReweightRvbarnNC1pi");
147  labelsGENIE.push_back("ReweightRvbarnNC2pi");
148  labelsGENIE.push_back("ReweightAhtBY");
149  labelsGENIE.push_back("ReweightBhtBY");
150  labelsGENIE.push_back("ReweightCV1uBY");
151  labelsGENIE.push_back("ReweightCV2uBY");
152  labelsGENIE.push_back("ReweightNormDISCC");
153  labelsGENIE.push_back("ReweightRnubarnuCC");
154  labelsGENIE.push_back("ReweightDISNuclMod");
155  labelsGENIE.push_back("ReweightAGKY_xF1pi");
156  labelsGENIE.push_back("ReweightAGKY_pT1pi");
157  labelsGENIE.push_back("ReweightFormZone");
158  labelsGENIE.push_back("ReweightTheta_Delta2Npi");
159  labelsGENIE.push_back("ReweightBR1gamma");
160  labelsGENIE.push_back("ReweightBR1eta");
161  labelsGENIE.push_back("ReweightMFP_N");
162  labelsGENIE.push_back("ReweightFrCEx_N");
163  labelsGENIE.push_back("ReweightFrElas_N");
164  labelsGENIE.push_back("ReweightFrInel_N");
165  labelsGENIE.push_back("ReweightFrAbs_N");
166  labelsGENIE.push_back("ReweightFrPiProd_N");
167  labelsGENIE.push_back("ReweightMFP_pi");
168  labelsGENIE.push_back("ReweightFrCEx_pi");
169  labelsGENIE.push_back("ReweightFrElas_pi");
170  labelsGENIE.push_back("ReweightFrInel_pi");
171  labelsGENIE.push_back("ReweightFrAbs_pi");
172  labelsGENIE.push_back("ReweightFrPiProd_pi");
173  labelsGENIE.push_back("mecScale");
174  labelsGENIE.push_back("RPA");
175 
176  std::string labelNDRock = "NDRock";
177 
178  // Set up maps to the decompositions/predictions (each flavor component)
179  // Nominal maps are indexed purely by sample label
180  // Shifted maps are indexed by sample label, systematic label, then sigma of the shift
181  std::map<std::string, IDecomp*> decomps_nominal;
182  std::map<std::string, std::map<std::string, std::map<int, IDecomp*> > > decomps_shifted;
183  std::map<std::string, PredictionNoExtrap*> predsNE_nominal;
184  std::map<std::string, std::map<std::string, std::map<int, PredictionNoExtrap*> > > predsNE_shifted;
185  std::map<std::string, PredictionSterile*> predsSt_nominal;
186  std::map<std::string, std::map<std::string, std::map<int, PredictionSterile*> > > predsSt_shifted;
187 
188  // Folder where all systematic files live
189  std::string folder = "/nova/ana/steriles/Ana01/Systematics/";
190 
191  // Output file name that will contain histograms for ISyst objects
192  std::string saveLocation = folder + "NusSystsAna01.root";
193 
194  TFile* rootF = new TFile(saveLocation.c_str(), "RECREATE");
195 
196  std::cout << "Beam" << std::endl; // Inform the user of progress
197  // Create maps to contain nominal spectra,
198  // -1, 0, +1 sigma ratio histograms
199  // First index is channel (NC/BG), second is det (ND/FD/EX)
200  std::map<std::string, std::map<std::string, TH1*> > nomsBeam;
201  std::map<std::string, std::map<std::string, TH1*> > sfp1Beam;
202  std::map<std::string, std::map<std::string, TH1*> > sfm1Beam;
203  std::map<std::string, std::map<std::string, TH1*> > lineBeam;
204 
205  // Open Beam systematic file, fill maps that are ready
206  // Second input is base file name without extension
207  FillMaps(folder, "SystsBeam", sample, nomsBeam,
208  decomps_nominal, decomps_shifted,
209  predsNE_nominal, predsNE_shifted,
210  predsSt_nominal, predsSt_shifted);
211 
212  // Add beam errors in quadrature
213  QuadErrors(labelsBeam, sample, nomsBeam,
214  decomps_shifted, predsNE_shifted, predsSt_shifted);
215 
216  // Populate +/-1 sigma maps, make sure they have the correct axis
217  CloneMap(nomsBeam, sfp1Beam);
218  CloneMap(nomsBeam, sfm1Beam);
219  // Populate 0 sigma map, set these histograms to be lines at 1
220  NominalRat(nomsBeam, lineBeam);
221 
222  // Set the +/-1 sigma histogram bins to appropiate ratio values
223  SetBinsPlusMinusOne(nomsBeam, sfp1Beam, sfm1Beam);
224 
225  // Save the histograms to file
226  SaveMaps(sfp1Beam, rootF, "Beam", "+1");
227  SaveMaps(sfm1Beam, rootF, "Beam", "-1");
228  SaveMaps(lineBeam, rootF, "Beam", "0");
229 
230  // Reset these maps as they are reused
231  decomps_nominal.clear();
232  decomps_shifted.clear();
233  predsNE_nominal.clear();
234  predsNE_shifted.clear();
235  predsSt_nominal.clear();
236  predsSt_shifted.clear();
237 
238  std::cout << "Birks" << std::endl;
239  std::map<std::string, std::map<std::string, TH1*> > nomsBirks;
240  std::map<std::string, std::map<std::string, TH1*> > shftBirks;
241  std::map<std::string, std::map<std::string, TH1*> > lineBirks;
242 
243  FillMaps(folder, "SystsBirks", sample, nomsBirks,
244  decomps_nominal, decomps_shifted,
245  predsNE_nominal, predsNE_shifted,
246  predsSt_nominal, predsSt_shifted);
247 
248  CloneMap(nomsBirks, shftBirks);
249  NominalRat(nomsBirks, lineBirks);
250 
251  SetBinsOnOffShift(labelBirks, sample, nomsBirks, shftBirks,
252  decomps_shifted, predsNE_shifted, predsSt_shifted);
253 
254  SaveMaps(shftBirks, rootF, "Birks", "+1");
255  SaveMaps(shftBirks, rootF, "Birks", "-1");
256  SaveMaps(lineBirks, rootF, "Birks", "0");
257 
258  decomps_nominal.clear();
259  decomps_shifted.clear();
260  predsNE_nominal.clear();
261  predsNE_shifted.clear();
262  predsSt_nominal.clear();
263  predsSt_shifted.clear();
264 
265  std::cout << "Calib, Both Dets" << std::endl;
266  std::map<std::string, std::map<std::string, TH1*> > nomsCalAbs;
267  std::map<std::string, std::map<std::string, TH1*> > sfp1CalFlat;
268  std::map<std::string, std::map<std::string, TH1*> > sfm1CalFlat;
269  std::map<std::string, std::map<std::string, TH1*> > lineCalFlat;
270  std::map<std::string, std::map<std::string, TH1*> > shftCalSlopeX;
271  std::map<std::string, std::map<std::string, TH1*> > lineCalSlopeX;
272  std::map<std::string, std::map<std::string, TH1*> > shftCalSlopeY;
273  std::map<std::string, std::map<std::string, TH1*> > lineCalSlopeY;
274 
275  FillMaps(folder, "SystsCalibAbs", sample, nomsCalAbs,
276  decomps_nominal, decomps_shifted,
277  predsNE_nominal, predsNE_shifted,
278  predsSt_nominal, predsSt_shifted);
279 
280  CloneMap(nomsCalAbs, sfp1CalFlat);
281  CloneMap(nomsCalAbs, sfm1CalFlat);
282  NominalRat(nomsCalAbs, lineCalFlat);
283  CloneMap(nomsCalAbs, shftCalSlopeX);
284  NominalRat(nomsCalAbs, lineCalSlopeX);
285  CloneMap(nomsCalAbs, shftCalSlopeY);
286  NominalRat(nomsCalAbs, lineCalSlopeY);
287 
288  SetBinsPlusMinusOne(labelsFlat, sample, nomsCalAbs, sfp1CalFlat, sfm1CalFlat,
289  decomps_shifted, predsNE_shifted, predsSt_shifted);
290  SetBinsOnOffShift(labelSlopeX, sample, nomsCalAbs, shftCalSlopeX,
291  decomps_shifted, predsNE_shifted, predsSt_shifted);
292  SetBinsOnOffShift(labelSlopeY, sample, nomsCalAbs, shftCalSlopeY,
293  decomps_shifted, predsNE_shifted, predsSt_shifted);
294 
295  SaveMaps(sfp1CalFlat, rootF, "CalFlat", "+1");
296  SaveMaps(sfm1CalFlat, rootF, "CalFlat", "-1");
297  SaveMaps(lineCalFlat, rootF, "CalFlat", "0");
298  SaveMaps(shftCalSlopeX, rootF, "CalSlopeX", "+1");
299  SaveMaps(shftCalSlopeX, rootF, "CalSlopeX", "-1");
300  SaveMaps(lineCalSlopeX, rootF, "CalSlopeX", "0");
301  SaveMaps(shftCalSlopeY, rootF, "CalSlopeY", "+1");
302  SaveMaps(shftCalSlopeY, rootF, "CalSlopeY", "-1");
303  SaveMaps(lineCalSlopeY, rootF, "CalSlopeY", "0");
304 
305  decomps_nominal.clear();
306  decomps_shifted.clear();
307  predsNE_nominal.clear();
308  predsNE_shifted.clear();
309  predsSt_nominal.clear();
310  predsSt_shifted.clear();
311 
312  std::cout << "Calib, Single Det" << std::endl;
313  std::map<std::string, std::map<std::string, TH1*> > nomsCalRel;
314  std::map<std::string, std::map<std::string, TH1*> > sfp1CalRel;
315  std::map<std::string, std::map<std::string, TH1*> > sfm1CalRel;
316  std::map<std::string, std::map<std::string, TH1*> > lineCalRel;
317 
318  FillMaps(folder, "SystsCalibRel", sample, nomsCalRel,
319  decomps_nominal, decomps_shifted,
320  predsNE_nominal, predsNE_shifted,
321  predsSt_nominal, predsSt_shifted);
322 
323  CloneMap(nomsCalRel, sfp1CalRel);
324  CloneMap(nomsCalRel, sfm1CalRel);
325  NominalRat(nomsCalRel, lineCalRel);
326 
327  SetBinsPlusMinusOne(labelsCalibRel, sample, nomsCalRel, sfp1CalRel, sfm1CalRel,
328  decomps_shifted, predsNE_shifted, predsSt_shifted);
329 
330  SaveMaps(sfp1CalRel, rootF, "CalRel", "+1");
331  SaveMaps(sfm1CalRel, rootF, "CalRel", "-1");
332  SaveMaps(lineCalRel, rootF, "CalRel", "0");
333 
334  decomps_nominal.clear();
335  decomps_shifted.clear();
336  predsNE_nominal.clear();
337  predsNE_shifted.clear();
338  predsSt_nominal.clear();
339  predsSt_shifted.clear();
340 
341  std::cout << "Decomp / CC Background / Data/MC Diff" << std::endl;
342  std::map<std::string, std::map<std::string, TH1*> > nomsDecomp;
343  std::map<std::string, std::map<std::string, TH1*> > shftNumu;
344  std::map<std::string, std::map<std::string, TH1*> > lineNumu;
345  std::map<std::string, std::map<std::string, TH1*> > sfp1Nue;
346  std::map<std::string, std::map<std::string, TH1*> > sfm1Nue;
347  std::map<std::string, std::map<std::string, TH1*> > lineNue;
348 
349  FillMaps(folder, "SystsDecomp", sample, nomsDecomp,
350  decomps_nominal, decomps_shifted,
351  predsNE_nominal, predsNE_shifted,
352  predsSt_nominal, predsSt_shifted);
353 
354  CloneMap(nomsDecomp, shftNumu);
355  NominalRat(nomsDecomp, lineNumu);
356  CloneMap(nomsDecomp, sfp1Nue);
357  CloneMap(nomsDecomp, sfm1Nue);
358  NominalRat(nomsDecomp, lineNue);
359 
360  SetBinsOnOffShift(labelNumu, sample, nomsDecomp, shftNumu,
361  decomps_shifted, predsNE_shifted, predsSt_shifted);
362  SetBinsPlusMinusOne(labelsNue, sample, nomsDecomp, sfp1Nue, sfm1Nue,
363  decomps_shifted, predsNE_shifted, predsSt_shifted);
364 
365  SaveMaps(shftNumu, rootF, "NumuCC", "+1");
366  SaveMaps(shftNumu, rootF, "NumuCC", "-1");
367  SaveMaps(lineNumu, rootF, "NumuCC", "0");
368  SaveMaps(sfp1Nue, rootF, "NueCC", "+1");
369  SaveMaps(sfm1Nue, rootF, "NueCC", "-1");
370  SaveMaps(lineNue, rootF, "NueCC", "0");
371 
372  decomps_nominal.clear();
373  decomps_shifted.clear();
374  predsNE_nominal.clear();
375  predsNE_shifted.clear();
376  predsSt_nominal.clear();
377  predsSt_shifted.clear();
378 
379  std::cout << "Small GENIE" << std::endl;
380  std::map<std::string, std::map<std::string, TH1*> > nomsGENIE;
381  std::map<std::string, std::map<std::string, TH1*> > sfp1GENIE;
382  std::map<std::string, std::map<std::string, TH1*> > sfm1GENIE;
383  std::map<std::string, std::map<std::string, TH1*> > lineGENIE;
384 
385  FillMaps(folder, "SystsGENIE", sample, nomsGENIE,
386  decomps_nominal, decomps_shifted,
387  predsNE_nominal, predsNE_shifted,
388  predsSt_nominal, predsSt_shifted);
389 
390  QuadErrors(labelsGENIE, sample, nomsGENIE,
391  decomps_shifted, predsNE_shifted, predsSt_shifted);
392 
393  CloneMap(nomsGENIE, sfp1GENIE);
394  CloneMap(nomsGENIE, sfm1GENIE);
395  NominalRat(nomsGENIE, lineGENIE);
396 
397  SetBinsPlusMinusOne(nomsGENIE, sfp1GENIE, sfm1GENIE);
398 
399  SaveMaps(sfp1GENIE, rootF, "GENIESm", "+1");
400  SaveMaps(sfm1GENIE, rootF, "GENIESm", "-1");
401  SaveMaps(lineGENIE, rootF, "GENIESm", "0");
402 
403  decomps_nominal.clear();
404  decomps_shifted.clear();
405  predsNE_nominal.clear();
406  predsNE_shifted.clear();
407  predsSt_nominal.clear();
408  predsSt_shifted.clear();
409 
410  std::cout << "ND Rock" << std::endl;
411  std::map<std::string, std::map<std::string, TH1*> > nomsRock;
412  std::map<std::string, std::map<std::string, TH1*> > shftRock;
413  std::map<std::string, std::map<std::string, TH1*> > lineRock;
414 
415  FillMaps(folder, "SystsNDRock", sample, nomsRock,
416  decomps_nominal, decomps_shifted,
417  predsNE_nominal, predsNE_shifted,
418  predsSt_nominal, predsSt_shifted);
419 
420  CloneMap(nomsRock, shftRock);
421  NominalRat(nomsRock, lineRock);
422 
423  SetBinsOnOffShift(labelNDRock, sample, nomsRock, shftRock,
424  decomps_shifted, predsNE_shifted, predsSt_shifted);
425 
426  SaveMaps(shftRock, rootF, "NDRock", "+1");
427  SaveMaps(shftRock, rootF, "NDRock", "-1");
428  SaveMaps(lineRock, rootF, "NDRock", "0");
429 
430  decomps_nominal.clear();
431  decomps_shifted.clear();
432  predsNE_nominal.clear();
433  predsNE_shifted.clear();
434  predsSt_nominal.clear();
435  predsSt_shifted.clear();
436 
437  rootF->Close();
438 
439  std::cout << "Finished. Ctrl+Z should be safe if the macro didn't quit." << std::endl;
440 }
441 
442 //------------------------------------------------------------------------------
443 void CloneMap(std::map<std::string, std::map<std::string, TH1*> >& noms,
444  std::map<std::string, std::map<std::string, TH1*> >& clone)
445 {
446  for(const auto& stringmap : noms) {
447  for(const auto& stringhist: stringmap.second) {
448  clone[stringmap.first][stringhist.first] = (TH1*)stringhist.second->Clone();
449  }
450  }
451 
452  return;
453 }
454 
455 //------------------------------------------------------------------------------
457  std::map<std::string, std::map<std::string, TH1*> >& noms,
458  std::map<std::string, IDecomp*>& decomps_nominal,
459  std::map<std::string, std::map<std::string, std::map<int, IDecomp*> > >& decomps_shifted,
460  std::map<std::string, PredictionNoExtrap*>& predsNE_nominal,
461  std::map<std::string, std::map<std::string, std::map<int, PredictionNoExtrap*> > >& predsNE_shifted,
462  std::map<std::string, PredictionSterile*>& predsSt_nominal,
463  std::map<std::string, std::map<std::string, std::map<int, PredictionSterile*> > >& predsSt_shifted)
464 {
465  std::string loadLocation = folder + systfile + ".root";
466  TFile* rootL = new TFile(loadLocation.c_str(), "READ");
467  LoadMaps(rootL, &decomps_nominal, &decomps_shifted);
468  LoadMaps(rootL, &predsNE_nominal, &predsNE_shifted);
469  LoadMaps(rootL, &predsSt_nominal, &predsSt_shifted);
470  rootL->Close();
471 
472  noms["NC"]["ND"] = GetNC(decomps_nominal[sample], NCSCALE);
473  noms["BG"]["ND"] = GetBG(decomps_nominal[sample], NCSCALE);
474  noms["NC"]["FD"] = GetNC(predsNE_nominal[sample], NCSCALE);
475  noms["BG"]["FD"] = GetBG(predsNE_nominal[sample], NCSCALE);
476  noms["NC"]["EX"] = GetNC(predsSt_nominal[sample], NCSCALE);
477  noms["BG"]["EX"] = GetBG(predsSt_nominal[sample], NCSCALE);
478 
479  return;
480 }
481 
482 //------------------------------------------------------------------------------
485 {
486  // This pattern must match that in NusSystFromHist::LoadHists
487  std::string ret = "h" + chan + "_" + det + "_" + syst + "_" + sigma;
488  return ret;
489 }
490 
491 //------------------------------------------------------------------------------
494 {
495  std::string ret = chan + " " + det + " " + syst + " Systematic, "
496  + sigma + "#sigma";
497  return ret;
498 }
499 
500 //------------------------------------------------------------------------------
501 void NominalRat(std::map<std::string, std::map<std::string, TH1*> >& noms,
502  std::map<std::string, std::map<std::string, TH1*> >& line)
503 {
504  CloneMap(noms, line); // Make sure axis matches
505 
507  for(int i = 0; i <= nbins; ++i) {
508  for(const auto& stringmap : line) {
509  std::string chan = stringmap.first;
510 
511  for(const auto& stringhist : stringmap.second) {
512  std::string det = stringhist.first;
513 
514  line[chan][det]->SetBinContent(i, 1.);
515  }
516  }
517  }
518 }
519 
520 //------------------------------------------------------------------------------
521 void QuadErrors(const std::vector<std::string>& shiftlabels, std::string sample,
522  std::map<std::string, std::map<std::string, TH1*> >& noms,
523  std::map<std::string, std::map<std::string, std::map<int, IDecomp*> > >& decomps_shifted,
524  std::map<std::string, std::map<std::string, std::map<int, PredictionNoExtrap*> > >& predsNE_shifted,
525  std::map<std::string, std::map<std::string, std::map<int, PredictionSterile*> > >& predsSt_shifted)
526 {
527  for(const auto& shift : shiftlabels) {
529  noms["NC"]["ND"],
530  GetNC(decomps_shifted[sample][shift][+1], NCSCALE),
531  GetNC(decomps_shifted[sample][shift][-1], NCSCALE)
532  );
534  noms["BG"]["ND"],
535  GetBG(decomps_shifted[sample][shift][+1], NCSCALE),
536  GetBG(decomps_shifted[sample][shift][-1], NCSCALE)
537  );
539  noms["NC"]["FD"],
540  GetNC(predsNE_shifted[sample][shift][+1], NCSCALE),
541  GetNC(predsNE_shifted[sample][shift][-1], NCSCALE)
542  );
544  noms["BG"]["FD"],
545  GetBG(predsNE_shifted[sample][shift][+1], NCSCALE),
546  GetBG(predsNE_shifted[sample][shift][-1], NCSCALE)
547  );
549  noms["NC"]["EX"],
550  GetNC(predsSt_shifted[sample][shift][+1], NCSCALE),
551  GetNC(predsSt_shifted[sample][shift][-1], NCSCALE)
552  );
554  noms["BG"]["EX"],
555  GetBG(predsSt_shifted[sample][shift][+1], NCSCALE),
556  GetBG(predsSt_shifted[sample][shift][-1], NCSCALE)
557  );
558  }
559 
560  // Don't forget to take the square root!
561  SqrtError(noms["NC"]["ND"]);
562  SqrtError(noms["BG"]["ND"]);
563  SqrtError(noms["NC"]["FD"]);
564  SqrtError(noms["BG"]["FD"]);
565  SqrtError(noms["NC"]["EX"]);
566  SqrtError(noms["BG"]["EX"]);
567 
568  return;
569 }
570 
571 //------------------------------------------------------------------------------
572 void SaveMaps(std::map<std::string, std::map<std::string, TH1*> >& hists,
573  TDirectory* out, std::string syst, std::string sigma)
574 {
575  for(const auto& stringmap : hists) {
576  for(const auto& stringhist: stringmap.second) {
577  std::string chan = stringmap.first;
578  std::string det = stringhist.first;
579 
580  hists[chan][det]->SetName( MakeName( chan, det, syst, sigma).c_str());
581  hists[chan][det]->SetTitle(MakeTitle(chan, det, syst, sigma).c_str());
582  ZeroError(hists[chan][det]);
583  out->WriteTObject(hists[chan][det]);
584  }
585  }
586 
587  return;
588 }
589 
590 //------------------------------------------------------------------------------
591 void SetBinsOnOffShift(std::string shiftlabel, std::string sample,
592  std::map<std::string, std::map<std::string, TH1*> >& noms,
593  std::map<std::string, std::map<std::string, TH1*> >& shft,
594  std::map<std::string, std::map<std::string, std::map<int, IDecomp*> > >& decomps_shifted,
595  std::map<std::string, std::map<std::string, std::map<int, PredictionNoExtrap*> > >& predsNE_shifted,
596  std::map<std::string, std::map<std::string, std::map<int, PredictionSterile*> > >& predsSt_shifted)
597 {
599  double nom = 0., sft = 0.;
600  std::string chan = "", det = "";
601 
602  std::map<std::string, std::map<std::string, TH1*> > sfts;
603  sfts["NC"]["ND"] = GetNC(decomps_shifted[sample][shiftlabel][1], NCSCALE);
604  sfts["BG"]["ND"] = GetBG(decomps_shifted[sample][shiftlabel][1], NCSCALE);
605  sfts["NC"]["FD"] = GetNC(predsNE_shifted[sample][shiftlabel][1], NCSCALE);
606  sfts["BG"]["FD"] = GetBG(predsNE_shifted[sample][shiftlabel][1], NCSCALE);
607  sfts["NC"]["EX"] = GetNC(predsSt_shifted[sample][shiftlabel][1], NCSCALE);
608  sfts["BG"]["EX"] = GetBG(predsSt_shifted[sample][shiftlabel][1], NCSCALE);
609 
610  for(int i = 0; i <= nbins; ++i) {
611  for(const auto& stringmap : noms) {
612  chan = stringmap.first;
613 
614  for(const auto& stringhist : stringmap.second) {
615  det = stringhist.first;
616 
617  nom = noms[chan][det]->GetBinContent(i);
618  sft = sfts[chan][det]->GetBinContent(i);
619 
620  // On/off systematics are just a straight ratio of shift/nominal
621  // Protect against negative bin values
622  shft[chan][det]->SetBinContent(i, std::max(0., sft/nom));
623  }
624  }
625  }
626 
627  return;
628 }
629 
630 //------------------------------------------------------------------------------
631 void SetBinsPlusMinusOne(std::pair<std::string, std::string> shiftlabels, std::string sample,
632  std::map<std::string, std::map<std::string, TH1*> >& noms,
633  std::map<std::string, std::map<std::string, TH1*> >& sfp1,
634  std::map<std::string, std::map<std::string, TH1*> >& sfm1,
635  std::map<std::string, std::map<std::string, std::map<int, IDecomp*> > >& decomps_shifted,
636  std::map<std::string, std::map<std::string, std::map<int, PredictionNoExtrap*> > >& predsNE_shifted,
637  std::map<std::string, std::map<std::string, std::map<int, PredictionSterile*> > >& predsSt_shifted)
638 {
639  // This syst type essentially consists of two on/off systematics,
640  // where one is +1 sigma and the other -1 sigma
641  SetBinsOnOffShift(shiftlabels.first, sample, noms, sfp1,
642  decomps_shifted, predsNE_shifted, predsSt_shifted);
643  SetBinsOnOffShift(shiftlabels.second, sample, noms, sfm1,
644  decomps_shifted, predsNE_shifted, predsSt_shifted);
645 }
646 
647 //------------------------------------------------------------------------------
648 void SetBinsPlusMinusOne(std::map<std::string, std::map<std::string, TH1*> >& noms,
649  std::map<std::string, std::map<std::string, TH1*> >& sfp1,
650  std::map<std::string, std::map<std::string, TH1*> >& sfm1)
651 {
653  double nom = 0., err = 0.;
654  std::string chan = "", det = "";
655 
656  for(int i = 0; i <= nbins; ++i) {
657  for(const auto& stringmap : noms) {
658  chan = stringmap.first;
659 
660  for(const auto& stringhist : stringmap.second) {
661  det = stringhist.first;
662 
663  nom = noms[chan][det]->GetBinContent(i);
664  err = noms[chan][det]->GetBinError(i);
665 
666  // This type of systematic has been added in quadrature
667  // The expected format is that the relative shift is the bin error
668  sfp1[chan][det]->SetBinContent(i, std::max(0., (nom+err)/nom));
669  sfm1[chan][det]->SetBinContent(i, std::max(0., (nom-err)/nom));
670  }
671  }
672  }
673 
674  return;
675 }
T max(const caf::Proxy< T > &a, T b)
std::string MakeName(std::string chan, std::string det, std::string syst, std::string sigma)
Definition: NusSystMaker.C:483
void SetBinsPlusMinusOne(std::pair< std::string, std::string > shiftlabels, std::string sample, std::map< std::string, std::map< std::string, TH1 * > > &noms, std::map< std::string, std::map< std::string, TH1 * > > &sfp1, std::map< std::string, std::map< std::string, TH1 * > > &sfm1, std::map< std::string, std::map< std::string, std::map< int, IDecomp * > > > &decomps_shifted, std::map< std::string, std::map< std::string, std::map< int, PredictionNoExtrap * > > > &predsNE_shifted, std::map< std::string, std::map< std::string, std::map< int, PredictionSterile * > > > &predsSt_shifted)
Definition: NusSystMaker.C:631
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
void SqrtError(TH1 *nom)
Set the errors of the input histogram to be their square roots.
Definition: PPFXHelper.h:139
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
void NominalRat(std::map< std::string, std::map< std::string, TH1 * > > &noms, std::map< std::string, std::map< std::string, TH1 * > > &line)
Definition: NusSystMaker.C:501
std::string MakeTitle(std::string chan, std::string det, std::string syst, std::string sigma)
Definition: NusSystMaker.C:492
void ZeroError(TH1 *nom)
Set the errors of the input histogram to 0.
Definition: PPFXHelper.h:149
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
TString hists[nhists]
Definition: bdt_com.C:3
TH1 * GetNC(IDecomp *specs, double POT)
Definition: PPFXHelper.h:211
const int nbins
Definition: cellShifts.C:15
void AddErrorInQuadrature(TH1 *nom, TH1 *up, bool use_max)
Definition: PPFXHelper.h:163
double sigma(TH1F *hist, double percentile)
OStream cout
Definition: OStream.cxx:6
void CloneMap(std::map< std::string, std::map< std::string, TH1 * > > &noms, std::map< std::string, std::map< std::string, TH1 * > > &clone)
Definition: NusSystMaker.C:443
int NBins() const
Definition: Binning.h:29
void QuadErrors(const std::vector< std::string > &shiftlabels, std::string sample, std::map< std::string, std::map< std::string, TH1 * > > &noms, std::map< std::string, std::map< std::string, std::map< int, IDecomp * > > > &decomps_shifted, std::map< std::string, std::map< std::string, std::map< int, PredictionNoExtrap * > > > &predsNE_shifted, std::map< std::string, std::map< std::string, std::map< int, PredictionSterile * > > > &predsSt_shifted)
Definition: NusSystMaker.C:521
void LoadMaps(TDirectory *dir, std::map< std::string, IDecomp * > *nominal, std::map< std::string, std::map< std::string, std::map< int, IDecomp * > > > *shifted)
Definition: PPFXHelper.h:273
static double NCSCALE
TH1 * GetBG(IDecomp *specs, double POT)
Definition: PPFXHelper.h:236
const Binning kNCDisappearanceEnergyBinning
Energy binnings used in Ana01 for nus extrapolation.
Definition: Binning.cxx:68
void FillMaps(std::string folder, std::string systfile, std::string sample, std::map< std::string, std::map< std::string, TH1 * > > &noms, 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)
Definition: NusSystMaker.C:456
void SetBinsOnOffShift(std::string shiftlabel, std::string sample, std::map< std::string, std::map< std::string, TH1 * > > &noms, std::map< std::string, std::map< std::string, TH1 * > > &shft, std::map< std::string, std::map< std::string, std::map< int, IDecomp * > > > &decomps_shifted, std::map< std::string, std::map< std::string, std::map< int, PredictionNoExtrap * > > > &predsNE_shifted, std::map< std::string, std::map< std::string, std::map< int, PredictionSterile * > > > &predsSt_shifted)
Definition: NusSystMaker.C:591
void NusSystMaker()
Definition: NusSystMaker.C:81
enum BeamMode string