SystsCalibAbsAna.C
Go to the documentation of this file.
4 
5 using namespace ana;
6 
7 void PlotSysts(TH1* nom, TH1* hup, TH1* hdn, TH1* hsx, TH1* hsy,
8  TDirectory* rootOUT, FILE* textOFS,
9  strings strs);
10 
12 {
13  TH1::AddDirectory(0);
14 
15  // A map of systematic labels
16  std::map<std::string, std::string> shiftlabels;
17  shiftlabels["CalFlat095"] = "Flat Miscalibration Down";
18  shiftlabels["CalFlat105"] = "Flat Miscalibration Up";
19  shiftlabels["CalSlopeX"] = "Sloped Miscalibration, X View";
20  shiftlabels["CalSlopeY"] = "Sloped Miscalibration, Y View";
21 
22  std::string labelRecoE = "Calorimetric Energy (GeV)";
23  std::vector<std::string> cut_samples;
24  cut_samples.push_back("Ana01");
25 
26  // Set up maps to the decompositions/predictions (each flavor component)
27  // Nominal maps are indexed purely by sample label
28  // Shifted maps are indexed by sample label, systematic label, then sigma of the shift
29  std::map<std::string, IDecomp*> decomps_nominal;
30  std::map<std::string, std::map<std::string, std::map<int, IDecomp*> > > decomps_shifted;
31  std::map<std::string, PredictionNoExtrap*> predsNE_nominal;
32  std::map<std::string, std::map<std::string, std::map<int, PredictionNoExtrap*> > > predsNE_shifted;
33  std::map<std::string, PredictionSterile*> predsSt_nominal;
34  std::map<std::string, std::map<std::string, std::map<int, PredictionSterile*> > > predsSt_shifted;
35 
36  std::string folder = "/nova/ana/steriles/Ana01/Systematics/";
37  std::string filenm = "SystsCalibAbs";
38 
39  std::string loadLocation = folder + filenm + ".root";
40  std::string saveLocation = folder + filenm + "Ana.root";
41  std::string textLocation = folder + filenm + "Ana.txt";
42 
43  TFile* rootL = new TFile(loadLocation.c_str(), "READ");
44  LoadMaps(rootL, &decomps_nominal, &decomps_shifted);
45  LoadMaps(rootL, &predsNE_nominal, &predsNE_shifted);
46  LoadMaps(rootL, &predsSt_nominal, &predsSt_shifted);
47  rootL->Close();
48 
49  TFile* rootF = new TFile(saveLocation.c_str(), "RECREATE");
50  FILE* textF;
51  textF = fopen(textLocation.c_str(), "w");
52 
53  for(const auto& sample : cut_samples) {
54  strings strs;
55  strs.fComponent = "NC";
56  strs.fDet = "ND";
58  strs.fSample = sample;
59  strs.fSystType = "CalibAbs";
60  strs.fXLabel = labelRecoE;
61 
62  InitializeSystText(textF, strs);
63 
64  TH1* hNomNC_decomp = GetNC(decomps_nominal[sample], NCSCALE);
65  TH1* hNomBG_decomp = GetBG(decomps_nominal[sample], NCSCALE);
66 
67  TH1* hNomNC_predNE = GetNC(predsNE_nominal[sample], NCSCALE);
68  TH1* hNomBG_predNE = GetBG(predsNE_nominal[sample], NCSCALE);
69 
70  TH1* hNomNC_predSt = GetNC(predsSt_nominal[sample], NCSCALE);
71  TH1* hNomBG_predSt = GetBG(predsSt_nominal[sample], NCSCALE);
72 
73  for(const auto& shifts : shiftlabels) {
74  std::cout << shifts.first << std::endl;
75  strs.fSystS = shifts.first;
76  strs.fSystL = shifts.second;
77 
78  // Flat systematics MUST be first,
79  // as they are handled with use_max boolean
80 
81  strs.fComponent = "NC";
82  strs.fDet = "ND";
83  if(shifts.first.find("Flat") != std::string::npos) {
84  PlotSyst(hNomNC_decomp,
85  decomps_shifted[sample][shifts.first][1],
86  rootF, textF, strs, NCSCALE, true, true);
87  }
88  else {
89  PlotSyst(hNomNC_decomp,
90  decomps_shifted[sample][shifts.first][1],
91  rootF, textF, strs, NCSCALE, true, false);
92  }
93 
94  strs.fComponent = "BG";
95  if(shifts.first.find("Flat") != std::string::npos) {
96  PlotSyst(hNomBG_decomp,
97  decomps_shifted[sample][shifts.first][1],
98  rootF, textF, strs, NCSCALE, false, true);
99  }
100  else {
101  PlotSyst(hNomBG_decomp,
102  decomps_shifted[sample][shifts.first][1],
103  rootF, textF, strs, NCSCALE, false, false);
104  }
105 
106  strs.fComponent = "NC";
107  strs.fDet = "FD";
108  if(shifts.first.find("Flat") != std::string::npos) {
109  PlotSyst(hNomNC_predNE,
110  predsNE_shifted[sample][shifts.first][1],
111  rootF, textF, strs, NCSCALE, true, true);
112  }
113  else {
114  PlotSyst(hNomNC_predNE,
115  predsNE_shifted[sample][shifts.first][1],
116  rootF, textF, strs, NCSCALE, true, false);
117  }
118 
119  strs.fComponent = "BG";
120  if(shifts.first.find("Flat") != std::string::npos) {
121  PlotSyst(hNomBG_predNE,
122  predsNE_shifted[sample][shifts.first][1],
123  rootF, textF, strs, NCSCALE, false, true);
124  }
125  else {
126  PlotSyst(hNomBG_predNE,
127  predsNE_shifted[sample][shifts.first][1],
128  rootF, textF, strs, NCSCALE, false, false);
129  }
130 
131  strs.fComponent = "NC";
132  strs.fDet = "EX";
133  if(shifts.first.find("Flat") != std::string::npos) {
134  PlotSyst(hNomNC_predSt,
135  predsSt_shifted[sample][shifts.first][1],
136  rootF, textF, strs, NCSCALE, true, true);
137  }
138  else {
139  PlotSyst(hNomNC_predSt,
140  predsSt_shifted[sample][shifts.first][1],
141  rootF, textF, strs, NCSCALE, true, false);
142  }
143 
144  strs.fComponent = "BG";
145  if(shifts.first.find("Flat") != std::string::npos) {
146  PlotSyst(hNomBG_predSt,
147  predsSt_shifted[sample][shifts.first][1],
148  rootF, textF, strs, NCSCALE, false, true);
149  }
150  else {
151  PlotSyst(hNomBG_predSt,
152  predsSt_shifted[sample][shifts.first][1],
153  rootF, textF, strs, NCSCALE, false, false);
154  }
155  }
156 
157  TH1* hFlUNC_decomp = GetNC(decomps_shifted[sample]["CalFlat105"][1], NCSCALE);
158  TH1* hFlUBG_decomp = GetBG(decomps_shifted[sample]["CalFlat105"][1], NCSCALE);
159  TH1* hFlDNC_decomp = GetNC(decomps_shifted[sample]["CalFlat095"][1], NCSCALE);
160  TH1* hFlDBG_decomp = GetBG(decomps_shifted[sample]["CalFlat095"][1], NCSCALE);
161  TH1* hSlXNC_decomp = GetNC(decomps_shifted[sample]["CalSlopeX"][1], NCSCALE);
162  TH1* hSlXBG_decomp = GetBG(decomps_shifted[sample]["CalSlopeX"][1], NCSCALE);
163  TH1* hSlYNC_decomp = GetNC(decomps_shifted[sample]["CalSlopeY"][1], NCSCALE);
164  TH1* hSlYBG_decomp = GetBG(decomps_shifted[sample]["CalSlopeY"][1], NCSCALE);
165 
166  TH1* hFlUNC_predNE = GetNC(predsNE_shifted[sample]["CalFlat105"][1], NCSCALE);
167  TH1* hFlUBG_predNE = GetBG(predsNE_shifted[sample]["CalFlat105"][1], NCSCALE);
168  TH1* hFlDNC_predNE = GetNC(predsNE_shifted[sample]["CalFlat095"][1], NCSCALE);
169  TH1* hFlDBG_predNE = GetBG(predsNE_shifted[sample]["CalFlat095"][1], NCSCALE);
170  TH1* hSlXNC_predNE = GetNC(predsNE_shifted[sample]["CalSlopeX"][1], NCSCALE);
171  TH1* hSlXBG_predNE = GetBG(predsNE_shifted[sample]["CalSlopeX"][1], NCSCALE);
172  TH1* hSlYNC_predNE = GetNC(predsNE_shifted[sample]["CalSlopeY"][1], NCSCALE);
173  TH1* hSlYBG_predNE = GetBG(predsNE_shifted[sample]["CalSlopeY"][1], NCSCALE);
174 
175  TH1* hFlUNC_predSt = GetNC(predsSt_shifted[sample]["CalFlat105"][1], NCSCALE);
176  TH1* hFlUBG_predSt = GetBG(predsSt_shifted[sample]["CalFlat105"][1], NCSCALE);
177  TH1* hFlDNC_predSt = GetNC(predsSt_shifted[sample]["CalFlat095"][1], NCSCALE);
178  TH1* hFlDBG_predSt = GetBG(predsSt_shifted[sample]["CalFlat095"][1], NCSCALE);
179  TH1* hSlXNC_predSt = GetNC(predsSt_shifted[sample]["CalSlopeX"][1], NCSCALE);
180  TH1* hSlXBG_predSt = GetBG(predsSt_shifted[sample]["CalSlopeX"][1], NCSCALE);
181  TH1* hSlYNC_predSt = GetNC(predsSt_shifted[sample]["CalSlopeY"][1], NCSCALE);
182  TH1* hSlYBG_predSt = GetBG(predsSt_shifted[sample]["CalSlopeY"][1], NCSCALE);
183 
184  strs.fComponent = "NC";
185  strs.fDet = "ND";
186  PlotSysts(hNomNC_decomp,
187  hFlUNC_decomp, hFlDNC_decomp,
188  hSlXNC_decomp, hSlYNC_decomp,
189  rootF, textF, strs);
190 
191  strs.fComponent = "BG";
192  PlotSysts(hNomBG_decomp,
193  hFlUBG_decomp, hFlDBG_decomp,
194  hSlXBG_decomp, hSlYBG_decomp,
195  rootF, textF, strs);
196 
197  strs.fComponent = "NC";
198  strs.fDet = "FD";
199  PlotSysts(hNomNC_predNE,
200  hFlUNC_predNE, hFlDNC_predNE,
201  hSlXNC_predNE, hSlYNC_predNE,
202  rootF, textF, strs);
203 
204  strs.fComponent = "BG";
205  PlotSysts(hNomBG_predNE,
206  hFlUBG_predNE, hFlDBG_predNE,
207  hSlXBG_predNE, hSlYBG_predNE,
208  rootF, textF, strs);
209 
210  strs.fComponent = "NC";
211  strs.fDet = "EX";
212  PlotSysts(hNomNC_predSt,
213  hFlUNC_predSt, hFlDNC_predSt,
214  hSlXNC_predSt, hSlYNC_predSt,
215  rootF, textF, strs);
216 
217  strs.fComponent = "BG";
218  PlotSysts(hNomBG_predSt,
219  hFlUBG_predSt, hFlDBG_predSt,
220  hSlXBG_predSt, hSlYBG_predSt,
221  rootF, textF, strs);
222  }
223 
224  fclose(textF);
225  rootF->Close();
226 
227  std::cout << "Finished. Ctrl+Z should be safe if the macro didn't quit." << std::endl;
228 }
229 
230 void PlotSysts(TH1* nom, TH1* hup, TH1* hdn, TH1* hsx, TH1* hsy,
231  TDirectory* rootOUT, FILE* textOFS,
232  strings strs)
233 {
234  double nNom = nom->Integral();
235  double n5Up = hup->Integral();
236  double n5Dn = hdn->Integral();
237  double nSlX = hsx->Integral();
238  double nSlY = hsy->Integral();
239 
240  fprintf(textOFS, "\nEvent counts for %.2s %.2s %s Systs:\n",
241  strs.fDet.c_str(), strs.fComponent.c_str(), strs.fSystType.c_str());
242  fprintf(textOFS, "NC: %6.2f, Flat Up: %6.2f, Flat Down: %6.2f, Slope X: %6.2f, Slope Y: %6.2f\n",
243  nNom, n5Up, n5Dn, nSlX, nSlY);
244  fprintf(textOFS, "Flat Up Diff: %6.2f, Flat Down Diff: %6.2f, Slope X Diff: %6.2f, Slope Y Diff: %6.2f",
245  100.*std::abs(n5Up/nNom - 1.), 100.*std::abs(n5Dn/nNom - 1.),
246  100.*std::abs(nSlX/nNom - 1.), 100.*std::abs(nSlY/nNom - 1.));
247 
248  std::string xlabel = strs.fXLabel;
249  std::string ylabel = "Events / " + strs.fPOT + " / 0.25 GeV";
250  std::string title = strs.fDet + " " + strs.fComponent +
251  " Errors for All " + strs.fSystType + " Systematics";
252 // std::string fullTitle = title + ";" + xlabel + ";" + ylabel;
253  std::string fullTitle = ";" + xlabel + ";" + ylabel;
254  std::string name = "c" + strs.fComponent + strs.fDet +
255  strs.fSystType + "SystsAll" + strs.fSample;
256  std::string ylabelRat = "Shifted / Nominal";
257  std::string fullTitleRat = ";" + xlabel + ";" + ylabelRat;
258 
259  TH1* hRat = (TH1*)nom->Clone();
260  TH1* hUpRat = (TH1*)hup->Clone();
261  TH1* hDnRat = (TH1*)hdn->Clone();
262  TH1* hSXRat = (TH1*)hsx->Clone();
263  TH1* hSYRat = (TH1*)hsy->Clone();
264 
265  int nBins = nom->GetNbinsX();
266  for(int i = 0; i <= nBins+1; ++i) {
267  hRat->SetBinContent(i, 1.);
268  }
269  hUpRat->Divide(nom);
270  hDnRat->Divide(nom);
271  hSXRat->Divide(nom);
272  hSYRat->Divide(nom);
273 
274  double maxval = nom->GetBinContent(nom->GetMaximumBin());
275  maxval = std::max(maxval, hup->GetBinContent(hup->GetMaximumBin()));
276  maxval = std::max(maxval, hdn->GetBinContent(hdn->GetMaximumBin()));
277  maxval = std::max(maxval, hsx->GetBinContent(hsx->GetMaximumBin()));
278  maxval = std::max(maxval, hsy->GetBinContent(hsy->GetMaximumBin()));
279 
280  double minval = nom->GetBinContent(nom->GetMinimumBin());
281  minval = std::min(minval, hup->GetBinContent(hup->GetMinimumBin()));
282  minval = std::min(minval, hdn->GetBinContent(hdn->GetMinimumBin()));
283  minval = std::min(minval, hsx->GetBinContent(hsx->GetMinimumBin()));
284  minval = std::min(minval, hsy->GetBinContent(hsy->GetMinimumBin()));
285 
286  if(maxval < minval) { std::swap(maxval, minval); }
287 
288  double maxvalRat = hRat->GetBinContent(hRat->GetMaximumBin());
289  maxvalRat = std::max(maxvalRat, hUpRat->GetBinContent(hUpRat->GetMaximumBin()));
290  maxvalRat = std::max(maxvalRat, hDnRat->GetBinContent(hDnRat->GetMaximumBin()));
291  maxvalRat = std::max(maxvalRat, hSXRat->GetBinContent(hSXRat->GetMaximumBin()));
292  maxvalRat = std::max(maxvalRat, hSYRat->GetBinContent(hSYRat->GetMaximumBin()));
293 
294  double minvalRat = hRat->GetBinContent(hRat->GetMinimumBin());
295  for(int i = 1; i <= nBins; ++i) {
296  if(hUpRat->GetBinContent(i) > 0.) {
297  minvalRat = std::min(minvalRat, hUpRat->GetBinContent(i));
298  }
299  if(hDnRat->GetBinContent(i) > 0.) {
300  minvalRat = std::min(minvalRat, hDnRat->GetBinContent(i));
301  }
302  if(hSXRat->GetBinContent(i) > 0.) {
303  minvalRat = std::min(minvalRat, hSXRat->GetBinContent(i));
304  }
305  if(hSYRat->GetBinContent(i) > 0.) {
306  minvalRat = std::min(minvalRat, hSYRat->GetBinContent(i));
307  }
308  }
309 
310  if(std::abs(maxvalRat - 1.) > std::abs(1. - minvalRat)) {
311  minvalRat = 1 - std::abs(maxvalRat - 1.);
312  maxvalRat = 1 + std::abs(maxvalRat - 1.);
313  }
314  else {
315  minvalRat = 1 - std::abs(1. - minvalRat);
316  maxvalRat = 1 + std::abs(1. - minvalRat);
317  }
318 
319  CenterTitles(nom);
320  nom->SetLineColor(kBlack);
321  nom->SetMaximum(1.1*maxval);
322  nom->SetMinimum(0);
323  nom->SetTitle(fullTitle.c_str());
324 
325  CenterTitles(hup);
326  hup->SetLineColor(kRed);
327  hup->SetLineStyle(2);
328  hup->SetMaximum(1.1*maxval);
329  hup->SetMinimum(0);
330  hup->SetTitle(fullTitle.c_str());
331 
332  CenterTitles(hdn);
333  hdn->SetLineColor(kRed);
334  hdn->SetLineStyle(5);
335  hdn->SetMaximum(1.1*maxval);
336  hdn->SetMinimum(0);
337  hdn->SetTitle(fullTitle.c_str());
338 
339  CenterTitles(hsx);
340  hsx->SetLineColor(kBlue);
341  hsx->SetLineStyle(2);
342  hsx->SetMaximum(1.1*maxval);
343  hsx->SetMinimum(0);
344  hsx->SetTitle(fullTitle.c_str());
345 
346  CenterTitles(hsy);
347  hsy->SetLineColor(kBlue);
348  hsy->SetLineStyle(5);
349  hsy->SetMaximum(1.1*maxval);
350  hsy->SetMinimum(0);
351  hsy->SetTitle(fullTitle.c_str());
352 
353  CenterTitles(hRat);
354  hRat->SetLineColor(kBlack);
355  hRat->SetMaximum(1.1*maxvalRat);
356  hRat->SetMinimum(0.9*minvalRat);
357  hRat->SetTitle(fullTitleRat.c_str());
358 
359  CenterTitles(hUpRat);
360  hUpRat->SetLineColor(kRed);
361  hUpRat->SetLineStyle(2);
362  hUpRat->SetMaximum(1.1*maxvalRat);
363  hUpRat->SetMinimum(0.9*minvalRat);
364  hUpRat->SetTitle(fullTitleRat.c_str());
365 
366  CenterTitles(hDnRat);
367  hDnRat->SetLineColor(kRed);
368  hDnRat->SetLineStyle(5);
369  hDnRat->SetMaximum(1.1*maxvalRat);
370  hDnRat->SetMinimum(0.9*minvalRat);
371  hDnRat->SetTitle(fullTitleRat.c_str());
372 
373  CenterTitles(hSXRat);
374  hSXRat->SetLineColor(kBlue);
375  hSXRat->SetLineStyle(2);
376  hSXRat->SetMaximum(1.1*maxvalRat);
377  hSXRat->SetMinimum(0.9*minvalRat);
378  hSXRat->SetTitle(fullTitleRat.c_str());
379 
380  CenterTitles(hSYRat);
381  hSYRat->SetLineColor(kBlue);
382  hSYRat->SetLineStyle(5);
383  hSYRat->SetMaximum(1.1*maxvalRat);
384  hSYRat->SetMinimum(0.9*minvalRat);
385  hSYRat->SetTitle(fullTitleRat.c_str());
386 
387  double xL = 0.6, xR = 0.85;
388  double yB = 0.6, yT = 0.85;
389  TLegend* leg = new TLegend(xL, yB, xR, yT);
390  leg->AddEntry(nom, "Nominal", "l");
391  leg->AddEntry(hup, "Flat Scale Up", "l");
392  leg->AddEntry(hdn, "Flat Scale Down", "l");
393  leg->AddEntry(hsx, "Slope X View", "l");
394  leg->AddEntry(hsy, "Slope Y View", "l");
395 
396  TCanvas* c = new TCanvas(name.c_str(), title.c_str(), 800, 800);
397  TPad* pSpecs = new TPad("pSpecs", "", 0., 0.375, 1., 1.);
398  TPad* pRatio = new TPad("pRatio", "", 0., 0., 1., 0.375);
399  pSpecs->Draw();
400  pRatio->Draw();
401 
402  pSpecs->cd();
403  nom->Draw("hist");
404  hup->Draw("hist same");
405  hdn->Draw("hist same");
406  hsx->Draw("hist same");
407  hsy->Draw("hist same");
408  leg->Draw();
409  Simulation();
410 
411  pRatio->cd();
412  hRat->Draw("hist");
413  hUpRat->Draw("hist same");
414  hDnRat->Draw("hist same");
415  hSXRat->Draw("hist same");
416  hSYRat->Draw("hist same");
417 
418  c->Update();
419  rootOUT->WriteTObject(c);
420  return;
421 }
void SystsCalibAbsAna()
void Simulation()
Definition: tools.h:16
T max(const caf::Proxy< T > &a, T b)
const XML_Char * name
Definition: expat.h:151
int nBins
Definition: plotROC.py:16
std::string fDet
Definition: PPFXHelper.h:104
enum BeamMode kRed
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
std::string fComponent
Definition: PPFXHelper.h:103
void PlotSyst(TH1 *h, TH1 *hp1, TH1 *hm1, TH1 *hp2, TH1 *hm2, TDirectory *out, FILE *text, strings strs)
Definition: PPFXHelper.h:662
void PlotSysts(TH1 *nom, TH1 *hup, TH1 *hdn, TH1 *hsx, TH1 *hsy, TDirectory *rootOUT, FILE *textOFS, strings strs)
void CenterTitles(TH1 *histo)
Definition: Plots.cxx:1481
float abs(float number)
Definition: d0nt_math.hpp:39
void swap(art::HLTGlobalStatus &lhs, art::HLTGlobalStatus &rhs)
TH1 * GetNC(IDecomp *specs, double POT)
Definition: PPFXHelper.h:211
fclose(fg1)
void InitializeSystText(FILE *text, strings strs)
Print some initial text about a systematic – the systematic type and sample.
Definition: PPFXHelper.h:263
std::string fSystL
Definition: PPFXHelper.h:109
std::string fPOT
Definition: PPFXHelper.h:105
OStream cout
Definition: OStream.cxx:6
std::string fXLabel
Definition: PPFXHelper.h:110
std::string fSystType
Definition: PPFXHelper.h:107
std::string fSystS
Definition: PPFXHelper.h:108
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
T min(const caf::Proxy< T > &a, T b)
TH1 * GetBG(IDecomp *specs, double POT)
Definition: PPFXHelper.h:236
enum BeamMode kBlue
std::string StringFromDouble(double pot)
Definition: PPFXHelper.h:128
A helper structure to contain a group of string for plotting.
Definition: PPFXHelper.h:101
std::string fSample
Definition: PPFXHelper.h:106
enum BeamMode string