PlotXSec.C
Go to the documentation of this file.
1 // PlotXSec.C
2 
5 
6 // Non-standard Multiverse include
7 #include "NDAna/numucc_inc/Multiverse.h"
8 
9 #include "TH1.h"
10 #include "TH2.h"
11 #include "TH3.h"
12 #include "TFile.h"
13 #include "TAxis.h"
14 #include "TCanvas.h"
15 #include "TLegend.h"
16 
17 #include <string>
18 #include <vector>
19 #include <map>
20 #include <iostream>
21 
22 map<string, string> systTitles
23 {
24  {"cv", "Central Value"},
25  {"light", "Light Level"},
26  {"calib", "Calibration"},
27  {"calibshift", "Calib Shape"},
28  {"ckv", "Cherenkov"},
29  {"muone", "Muon Energy"},
30  {"neutron", "Neutron Syst"},
31  {"angle", "Muon Angle"},
32  {"neutron_nom","Nominal (Neutron syst)"},
33  {"nom_good_seventh", "Nominal (calib / light levels)"},
34  {"genie", "XSec and FSI"},
35  {"ppfx", "Beam - HP"},
36  {"stat", "Statistical"},
37  {"beam", "Beam - Focusing"}
38 };
39 
40 vector<string> standardSysts
41 {
42  // "cv",
43  "calibshift",
44  "ckv"
45 };
46 
47 map<string, pair<string, string>> standardSystsUpDown
48 {
49  {"light", {"lightup", "lightdown"}},
50  {"calib", {"calibpos", "calibneg"}},
51  {"muone", {"muoneup", "muonedw"}},
52  {"angle", {"angleup", "angledw"}},
53  {"neutron", {"neutron_Up", "neutron_Dw"}}
54 };
55 
56 const map<string, string> otherNomSysts
57 {
58  {"neutron", "neutron_nom"},
59  {"light", "nom_good_seventh"},
60  {"calib", "nom_good_seventh"},
61  {"calibshift", "nom_good_seventh"}
62 };
63 
64 vector<string> beamSysts
65 {
66  "2kA",
67  "02mmBeamSpotSize",
68  "1mmBeamShiftX",
69  "1mmBeamShiftY",
70  "3mmHorn1X",
71  "3mmHorn1Y",
72  "3mmHorn2X",
73  "3mmHorn2Y",
74  "7mmTargetZ",
75  "MagneticFieldinDecayPipe",
76  "1mmHornWater"
77 };
78 
79 map<string, Color_t> systColors
80 {
81  {"total", kBlack},
82  {"beam", kBlue},
83  {"ckv", kOrange},
84  {"neutron", kGreen},
85  {"angle", kRed},
86  {"muone", kViolet},
87  {"stat", kGray},
88  {"genie", kCyan},
89  {"ppfx", kOrange}
90 };
91 
92 // Slice map of what which bins to present
93 map<int, double> maxTMuBins
94 {
95  {2, 1.1}, // 0.50-0.56
96  {3, 1.1}, // 0.56-0.62
97  {4, 1.1}, // 0.62-0.68
98  {5, 1.4}, // 0.68-0.74
99  {6, 1.4}, // 0.74-0.80
100  {7, 1.4}, // 0.80-0.85
101  {8, 1.8}, // 0.85-0.88
102  {9, 1.9}, // 0.88-0.91
103  {10, 2.2}, // 0.91-0.94
104 };
105 
106 const string var3DName = "EAvail";
107 string uncertLabel = "_uncert";
108 
109 const vector<string> vars3D {
110  "EAvail"
111 };
112 
113 // Nominal systs
114 TH3D * nom3DXSec;
115 TH2D * nom2DXSec;
116 
117 // Non-standard nominal (neutron, calibration, light-levels, etc)
118 string defaultNominal = "cv";
119 map<string, TH3D*> nominalXSecs3D;
120 map<string, TH2D*> nominalXSecs2D;
121 map<string, TH1D*> nominalXSecsENu;
122 map<string, TH1D*> nominalXSecsQ2;
123 
124 // Hold every TH2D we make
125 map<string, TH3D*> all3DXSecs;
126 map<string, TH2D*> xSecs2D;
127 map<string, TH1D*> xSecsENu;
128 map<string, TH1D*> xSecsQ2;
129 map<string, TH2D*> uncerts2D;
130 map<string, TH1D*> uncertsENu;
131 map<string, TH1D*> uncertsQ2;
132 
133 // File handlers
134 TFile * inputFile;
135 TFile * outFile;
136 
137 // Canvas
138 TCanvas * c1;
139 
140 // Define the visible axis to use.
141 const TAxis angleAxis(1, 0.5, 1.0);
142 const TAxis muoneAxis(1, 0.5, 2.5);
143 const TAxis uncertAxis(1, -0.2, 0.2);
144 const TAxis enuAxis(1, 0.0, 5.0);
145 const TAxis q2Axis(1, 0.0, 2.8);
146 
147 // Genie and PPFX
148 const unsigned int N_GENIE = 1000;
149 const unsigned int N_PPFX = 100;
150 
152 {
153  hist->GetXaxis()->SetRangeUser(q2Axis.GetXmin(), q2Axis.GetXmax());
154  // string histName = hist->GetName();
155  // if (histName.find(uncertLabel) != string::npos)
156  // hist->GetYaxis()->SetRangeUser(uncertAxis.GetXmin(), uncertAxis.GetXmax());
157 }
158 
160 {
161  hist->GetXaxis()->SetRangeUser(enuAxis.GetXmin(), enuAxis.GetXmax());
162  // string histName = hist->GetName();
163  // if (histName.find(uncertLabel) != string::npos)
164  // hist->GetYaxis()->SetRangeUser(uncertAxis.GetXmin(), uncertAxis.GetXmax());
165 }
166 
167 
169 {
170  hist->GetXaxis()->SetRangeUser(angleAxis.GetXmin(), angleAxis.GetXmax());
171  hist->GetYaxis()->SetRangeUser(muoneAxis.GetXmin(), muoneAxis.GetXmax());
172  // string histName = hist->GetName();
173  // if (histName.find(uncertLabel) != string::npos)
174  // hist->GetZaxis()->SetRangeUser(uncertAxis.GetXmin(), uncertAxis.GetXmax());
175 }
176 
177 void AddUpDwUncerts(TH1* destHist, TH1* upHist, TH1* dwHist)
178 {
179  assert(destHist->GetNbinsX() == upHist->GetNbinsX() && destHist->GetNbinsX() == dwHist->GetNbinsX());
180  for(int ibinx = 0; ibinx <= destHist->GetNbinsX() + 1; ++ibinx)
181  {
182  double upErr = abs(upHist->GetBinContent(ibinx));
183  double dwErr = abs(dwHist->GetBinContent(ibinx));
184  double addErr = max(upErr, dwErr);
185  double currentErr = destHist->GetBinContent(ibinx);
186  destHist->SetBinContent(ibinx, sqrt( pow(currentErr, 2) + pow(addErr, 2)));
187  }
188 }
189 
190 void AddUpDwUncerts2D(TH2* destHist, TH2* upHist, TH2* dwHist)
191 {
192  assert(destHist->GetNbinsX() == upHist->GetNbinsX() && destHist->GetNbinsX() == dwHist->GetNbinsX());
193  assert(destHist->GetNbinsY() == upHist->GetNbinsY() && destHist->GetNbinsY() == dwHist->GetNbinsY());
194  for(int ibinx = 0; ibinx <= destHist->GetNbinsX() + 1; ++ibinx)
195  for (int ibiny = 0; ibiny <= destHist->GetNbinsY() + 1; ++ibiny)
196  {
197  double upErr = abs(upHist->GetBinContent(ibinx, ibiny));
198  double dwErr = abs(dwHist->GetBinContent(ibinx, ibiny));
199  double addErr = max(upErr, dwErr);
200  double currentErr = destHist->GetBinContent(ibinx, ibiny);
201  destHist->SetBinContent(ibinx, ibiny, sqrt( pow(currentErr, 2) + pow(addErr, 2)));
202  }
203 }
204 
205 // Get a 3D XSec, loading from file only if necessary
206 TH3D* Load3DXSec(string syst)
207 {
208  if (all3DXSecs.find(syst) != all3DXSecs.end())
209  return all3DXSecs.at(syst);
210 
211  TH3D* xsec3D = (TH3D*) inputFile->Get((syst + "/" + var3DName + "/xsec_" + syst + "_" + var3DName).c_str());
212  all3DXSecs[syst] = xsec3D;
213  if (!xsec3D)
214  return NULL;
215  xsec3D->SetTitle((systTitles[syst] + " XSec").c_str());
216  xsec3D->Write();
217  return xsec3D;
218 }
219 
220 // Get 2D XSec
221 TH2D* XSec2D(string syst)
222 {
223  if (xSecs2D.find(syst) != xSecs2D.end())
224  return xSecs2D.at(syst);
225 
226  TH3D * syst3DXSec = Load3DXSec(syst);
227  TH2D * syst2DXSec;
228  if (!syst3DXSec)
229  {
230  xSecs2D[syst] = NULL;
231  return NULL;
232  }
233  syst2DXSec = (TH2D*) syst3DXSec->Project3D("yx");
234  syst2DXSec->Scale(1, "width"); // Area-normalize the xsec
235  FixAxisRange2D(syst2DXSec);
236  xSecs2D[syst] = syst2DXSec;
237  syst2DXSec->SetName(("xsec_" + syst + "_2DProj").c_str());
238  syst2DXSec->SetTitle((systTitles[syst] + " XSec").c_str());
239  syst2DXSec->Write();
240  return syst2DXSec;
241 }
242 
243 // XSec ENu
244 TH1D* XSecENu(string syst)
245 {
246  if (xSecsENu.find(syst) != xSecsENu.end())
247  return xSecsENu.at(syst);
248 
249  TH1D * xSecENu = (TH1D*) inputFile->Get((syst + "/ENu/xsec_" + syst + "_ENu").c_str());
250  FixAxisRangeENu(xSecENu);
251  xSecsENu[syst] = xSecENu;
252  xSecENu->SetTitle((systTitles[syst] + " XSec").c_str());
253  xSecENu->Write();
254  return xSecENu;
255 }
256 
257 // XSec Q2
258 TH1D* XSecQ2(string syst)
259 {
260  if (xSecsQ2.find(syst) != xSecsQ2.end())
261  return xSecsQ2.at(syst);
262 
263  TH1D * xSecQ2 = (TH1D*) inputFile->Get((syst + "/Q2/xsec_" + syst + "_Q2").c_str());
264  xSecQ2->Scale(1, "width");
265  FixAxisRangeQ2(xSecQ2);
266  xSecsQ2[syst] = xSecQ2;
267  xSecQ2->SetTitle((systTitles[syst] + " XSec").c_str());
268  xSecQ2->Write();
269  return xSecQ2;
270 }
271 
272 TH2D* XSec2DUncert(string syst)
273 {
274  if (uncerts2D.find(syst) != uncerts2D.end())
275  return uncerts2D.at(syst);
276 
277  assert(xSecs2D.find(syst) != xSecs2D.end());
278  TH2D * systXSec = xSecs2D.at(syst);
279  if (!systXSec){
280  uncerts2D[syst] = NULL;
281  return NULL;
282  }
283 
284  TH2D * thisNom = nominalXSecs2D["cv"];
285  if (otherNomSysts.find(syst) != otherNomSysts.end()){
286  cout << "Using non-standard nominal : " << otherNomSysts.at(syst) << " for syst " << syst << endl;
287  thisNom = nominalXSecs2D[otherNomSysts.at(syst)];
288  }
289  TH2D * syst2DUncert = new TH2D(*systXSec - *thisNom);
290  syst2DUncert->SetName((syst + "_2d_uncert").c_str());
291  syst2DUncert->SetTitle((syst + " Uncert").c_str());
292  syst2DUncert->Divide(thisNom);
293  FixAxisRange2D(syst2DUncert);
294  // syst2DUncert->GetZaxis()->SetRangeUser(uncertAxis.GetXmin(), uncertAxis.GetXmax());
295  uncerts2D[syst] = syst2DUncert;
296  syst2DUncert->Write();
297  return syst2DUncert;
298 }
299 
300 TH1D* XSecENuUncert(string syst)
301 {
302  if (uncertsENu.find(syst) != uncertsENu.end())
303  return uncertsENu.at(syst);
304 
305  assert(xSecsENu.find(syst) != xSecsENu.end());
306  TH1D * systXSec = xSecsENu.at(syst);
307 
308  TH1D * thisNom = nominalXSecsENu["cv"];
309  if (otherNomSysts.find(syst) != otherNomSysts.end()){
310  cout << "Using non-standard nominal : " << otherNomSysts.at(syst) << " for syst " << syst << endl;
311  thisNom = nominalXSecsENu[otherNomSysts.at(syst)];
312  }
313  TH1D * systUncert = new TH1D(*systXSec - *thisNom);
314  systUncert->SetName((syst + "_enu_uncert").c_str());
315  systUncert->SetTitle((syst + " Uncert").c_str());
316  systUncert->Divide(thisNom);
317  FixAxisRangeENu(systUncert);
318  // systUncert->GetZaxis()->SetRangeUser(uncertAxis.GetXmin(), uncertAxis.GetXmax());
319  uncertsENu[syst] = systUncert;
320  systUncert->Write();
321  return systUncert;
322 }
323 
324 TH1D* XSecQ2Uncert(string syst)
325 {
326  if (uncertsQ2.find(syst) != uncertsQ2.end())
327  return uncertsQ2.at(syst);
328 
329  assert(xSecsQ2.find(syst) != xSecsQ2.end());
330  TH1D * systXSec = xSecsQ2.at(syst);
331 
332  TH1D * thisNom = nominalXSecsQ2["cv"];
333  if (otherNomSysts.find(syst) != otherNomSysts.end()){
334  cout << "Using non-standard nominal : " << otherNomSysts.at(syst) << " for syst " << syst << endl;
335  thisNom = nominalXSecsQ2[otherNomSysts.at(syst)];
336  }
337  TH1D * systUncert = new TH1D(*systXSec - *thisNom);
338  systUncert->SetName((syst + "_q2_uncert").c_str());
339  systUncert->SetTitle((syst + " Uncert").c_str());
340  systUncert->Divide(thisNom);
341  FixAxisRangeQ2(systUncert);
342  // systUncert->GetZaxis()->SetRangeUser(uncertAxis.GetXmin(), uncertAxis.GetXmax());
343  uncertsQ2[syst] = systUncert;
344  systUncert->Write();
345  return systUncert;
346 }
347 
348 void LoadNominal(string syst)
349 {
350  if (nominalXSecs3D.find(syst) != nominalXSecs3D.end())
351  return;
352 
353  TDirectory * dir = outFile->GetDirectory(syst.c_str());
354  if (!dir)
355  dir = outFile->mkdir(syst.c_str());
356  dir->cd();
357 
358  TH3D * nomXSec = Load3DXSec(syst);
359  TH2D * nomXSec2D = XSec2D(syst);
360  TH1D * nomXSecENu = XSecENu(syst);
361  TH1D * nomXSecQ2 = XSecQ2(syst);
362 
363  nominalXSecs3D[syst] = nomXSec;
364  nominalXSecs2D[syst] = nomXSec2D;
365  nominalXSecsENu[syst] = nomXSecENu;
366  nominalXSecsQ2 [syst] = nomXSecQ2;
367 
368  outFile->cd();
369 }
370 
371 void LoadSyst(string syst)
372 {
373  TDirectory * dir = outFile->GetDirectory(syst.c_str());
374  if (!dir)
375  dir = outFile->mkdir(syst.c_str());
376  dir->cd();
377  TH3D * xSec3D = Load3DXSec(syst);
378  TH2D * xSec2D = XSec2D(syst);
379  TH1D * xSecENu = XSecENu(syst);
380  TH1D * xSecQ2 = XSecQ2(syst);
381 
382  outFile->cd();
383 }
384 
385 void LoadUncerts(string syst)
386 {
387  TDirectory * dir = outFile->GetDirectory(syst.c_str());
388  if (!dir)
389  dir = outFile->mkdir(syst.c_str());
390  dir->cd();
391  TH2 * uncert2D = XSec2DUncert(syst);
392  TH1 * uncertENu = XSecENuUncert(syst);
393  TH1 * uncertQ2 = XSecQ2Uncert(syst);
394 
395  outFile->cd();
396 }
397 
398 void LoadStatUncert(string nom)
399 {
400  if (uncerts2D. find("stat") != uncerts2D.end())
401  return; // Already loaded
402 
403  // Create Stat Uncerts
404  TH2D * statUncert2D = NULL;
405  if (nominalXSecs2D.find(nom) != nominalXSecs2D.end() && nominalXSecs2D[nom])
406  statUncert2D = (TH2D*) nominalXSecs2D[nom]->Clone("stat_uncert_2D");
407  for (int ibinx = 0; statUncert2D && ibinx <= statUncert2D->GetNbinsX() + 1; ++ibinx)
408  for (int ibiny = 0; statUncert2D && ibiny <= statUncert2D->GetNbinsY() + 1; ++ibiny){
409  double content = statUncert2D->GetBinContent(ibinx, ibiny);
410  double uncert = statUncert2D->GetBinError (ibinx, ibiny);
411  statUncert2D->SetTitle("Stat Uncert;#cos{#theta_{#mu}};T_{#mu};;");
412  statUncert2D->SetBinContent(ibinx, ibiny, content == 0.0 ? 0.0 : uncert / content);
413  statUncert2D->SetBinError(ibinx, ibiny, 0);
414  }
415  TH1D * statUncertENu = (TH1D*) nominalXSecsENu[nom]->Clone("stat_uncert_ENu");
416  for (int ibinx = 0; ibinx <= statUncertENu->GetNbinsX() + 1; ++ibinx){
417  double content = statUncertENu->GetBinContent(ibinx);
418  double uncert = statUncertENu->GetBinError(ibinx);
419  statUncertENu->SetTitle("Statistical Uncertainty;E_{nu} (GeV);;;");
420  statUncertENu->SetBinContent(ibinx, content == 0.0 ? 0.0 : uncert / content);
421  statUncertENu->SetBinError(ibinx, 0);
422  }
423  TH1D * statUncertQ2 = (TH1D*) nominalXSecsQ2[nom]->Clone("stat_uncert_Q2");
424  for (int ibinx = 0; ibinx <= statUncertQ2->GetNbinsX() + 1; ++ibinx){
425  double content = statUncertQ2->GetBinContent(ibinx);
426  double uncert = statUncertQ2->GetBinError(ibinx);
427  statUncertQ2->SetTitle("Statistical Uncertainty;Q^{2} (GeV^{2});;;");
428  statUncertQ2->SetBinContent(ibinx, content == 0.0 ? 0.0 : uncert / content);
429  statUncertQ2->SetBinError(ibinx, 0);
430  }
431  if (statUncert2D)
432  uncerts2D ["stat"] = statUncert2D;
433  uncertsENu["stat"] = statUncertENu;
434  uncertsQ2 ["stat"] = statUncertQ2;
435 }
436 
437 void AddQuadrature(TH1* dest, TH1* hist)
438 {
439  assert(dest && hist && dest->GetNbinsX() == hist->GetNbinsX());
440  for(int ibin = 0; ibin <= hist->GetNbinsX() + 1; ++ibin)
441  dest->SetBinContent(ibin, sqrt(pow(dest->GetBinContent(ibin), 2) + pow(hist->GetBinContent(ibin), 2)));
442 }
443 
444 void AddQuadrature(TH1* dest, vector<TH1*> hists)
445 {
446  assert(!hists.empty());
447  for (TH1 * hist : hists)
448  AddQuadrature(dest, hist);
449 }
450 
451 void AddQuadrature2D(TH2* dest, TH2* hist)
452 {
453  assert(dest && hist && dest->GetNbinsX() == hist->GetNbinsX() && dest->GetNbinsY() == hist->GetNbinsY());
454  for(int ibinx = 0; ibinx <= dest->GetNbinsX() + 1; ++ibinx)
455  for(int ibiny = 0; ibiny <= dest->GetNbinsY() + 1; ++ibiny)
456  dest->SetBinContent(ibinx, ibiny, sqrt(pow(dest->GetBinContent(ibinx, ibiny), 2) + pow(hist->GetBinContent(ibinx, ibiny), 2)));
457  // for(int ibinx = 0; ibinx <= dest->GetNbinsX() + 1; ++ibinx)
458  // for (int ibiny = 0; ibiny <= dest->GetNbinsY() + 1; ++ibiny)
459  // {
460  // double addErr = hist->GetBinContent(ibinx, ibiny);
461  // double currentErr = dest->GetBinContent(ibinx, ibiny);
462  // dest->SetBinContent(ibinx, ibiny, sqrt( pow(currentErr, 2) + pow(addErr, 2)));
463  // }
464 }
465 
466 void AddQuadrature2D(TH2* dest, vector<TH2*> hists)
467 {
468  assert(!hists.empty());
469  for (TH2 * hist : hists)
470  AddQuadrature2D(dest, hist);
471 }
472 
473 TH1 * GetMaxUpDown(TH1 * upHist, TH1 * dwHist)
474 {
475  TH1 * upDwError = (TH1*) upHist->Clone(ana::UniqueName().c_str());
476  upDwError->Reset();
477  for (int ibin = 0; ibin <= upHist->GetNbinsX() + 1; ++ibin)
478  upDwError->SetBinContent(ibin, max(abs(upHist->GetBinContent(ibin)), abs(dwHist->GetBinContent(ibin))));
479  return upDwError;
480 }
481 
482 TH2 * GetMaxUpDown2D(TH2 * upHist, TH2 * dwHist)
483 {
484  TH2 * upDwError = (TH2*) upHist->Clone(ana::UniqueName().c_str());
485  upDwError->Reset();
486  for (int ibinx = 0; ibinx <= upHist->GetNbinsX() + 1; ++ibinx)
487  for (int ibiny = 0; ibiny <= upHist->GetNbinsY() + 1; ++ibiny)
488  upDwError->SetBinContent(ibinx, ibiny, max(abs(upHist->GetBinContent(ibinx, ibiny)), abs(dwHist->GetBinContent(ibinx, ibiny))));
489  return upDwError;
490 }
491 
492 // Absolute for 1D
493 TH1* GetAbsolute(TH1* hist)
494 {
495  TH1 * histAbs = (TH1*) hist->Clone(ana::UniqueName().c_str());
496  for(int ibin = 0; ibin <= hist->GetNbinsX(); ++ibin)
497  histAbs->SetBinContent(ibin, abs(hist->GetBinContent(ibin)));
498  return histAbs;
499 }
500 
501 // Absolute for 2D
502 TH2* GetAbsolute2D(TH2* hist)
503 {
504  TH2 * histAbs = (TH2*) hist->Clone(ana::UniqueName().c_str());
505  for(int ibinx = 0; ibinx <= hist->GetNbinsX(); ++ibinx)
506  for(int ibiny = 0; ibiny <= hist->GetNbinsY(); ++ibiny)
507  histAbs->SetBinContent(ibinx, ibiny, abs(hist->GetBinContent(ibinx, ibiny)));
508  return histAbs;
509 }
510 
511 // Get slice
512 TH1 * GetSliceHist(TH2* hist, int ibinx)
513 {
514  // TH1 * histSlice = hist->ProjectionY((hist->GetName() + string("_") + to_string(ibinx)).c_str(), ibinx, ibinx);
515  TH1 * histSlice = hist->ProjectionY((ana::UniqueName()).c_str(), ibinx, ibinx);
516  TH1 * absSlice = GetAbsolute(histSlice);
517  absSlice->SetName((histSlice->GetName() + string("_abs")).c_str());
518  absSlice->SetTitle((to_string(histSlice->GetXaxis()->GetXmin()) + " < cos(#theta_{#mu}) < " + to_string(histSlice->GetXaxis()->GetXmax())).c_str());
519  return absSlice;
520 }
521 
522 // Combine the beam logic
524 {
525  TH2D * totalBeamUncert = (TH2D*) nominalXSecs2D["cv"]->Clone("totalBeamUncert2D");
526  totalBeamUncert->Reset();
527  for (string beamSyst : beamSysts)
528  {
529  TH2 * upSyst = uncerts2D[beamSyst + "_Up"];
530  TH2 * dwSyst = uncerts2D[beamSyst + "_Dw"];
531  AddUpDwUncerts2D(totalBeamUncert, upSyst, dwSyst);
532  }
533  uncerts2D["totalBeam"] = totalBeamUncert;
534  return totalBeamUncert;
535 }
536 
538 {
539  TH1D * totalBeamUncert = (TH1D*) nominalXSecsENu["cv"]->Clone("totalBeamUncertENu");
540  totalBeamUncert->Reset();
541  for (string beamSyst : beamSysts)
542  {
543  TH1 * upSyst = uncertsENu[beamSyst + "_Up"];
544  TH1 * dwSyst = uncertsENu[beamSyst + "_Dw"];
545  AddUpDwUncerts(totalBeamUncert, upSyst, dwSyst);
546  }
547  uncertsENu["totalBeam"] = totalBeamUncert;
548  return totalBeamUncert;
549 }
550 
552 {
553  TH1D * totalBeamUncert = (TH1D*) nominalXSecsQ2["cv"]->Clone("totalBeamUncertQ2");
554  totalBeamUncert->Reset();
555  for (string beamSyst : beamSysts)
556  {
557  TH1 * upSyst = uncertsQ2[beamSyst + "_Up"];
558  TH1 * dwSyst = uncertsQ2[beamSyst + "_Dw"];
559  AddUpDwUncerts(totalBeamUncert, upSyst, dwSyst);
560  }
561  uncertsQ2["totalBeam"] = totalBeamUncert;
562  return totalBeamUncert;
563 }
564 
565 // Plot only the phase space we want
566 void LimitTMuBins(TH1* hist, int ibiny)
567 {
568  if (maxTMuBins.find(ibiny) != maxTMuBins.end())
569  for (int ibin = 1; ibin <= hist->GetNbinsX(); ++ibin)
570  if (hist->GetXaxis()->GetBinCenter(ibin) >= maxTMuBins.at(ibiny)){
571  hist->SetBinContent(ibin, 0);
572  hist->SetBinError(ibin, 0);
573  }
574 }
575 
576 ////////////////////////////////
577 // Plotting functions
578 ////////////////////////////////
579 void MakeSystematicPlots(vector<string> systs, string systFilename)
580 {
581  TH2 * totalUncert2D = NULL;
582  if (nominalXSecs2D.find("cv") != nominalXSecs2D.end() && nominalXSecs2D["cv"])
583  totalUncert2D = (TH2*) nominalXSecs2D["cv"] ->Clone("totalUncert2D");
584  TH1 * totalUncertENu = (TH1*) nominalXSecsENu["cv"]->Clone("totalUncertENu");
585  TH1 * totalUncertQ2 = (TH1*) nominalXSecsQ2["cv"] ->Clone("totalUncertQ2");
586 
587  // Switch the
588  if (totalUncert2D)
589  totalUncert2D->Reset();
590  totalUncertENu ->Reset();
591  totalUncertQ2 ->Reset();
592 
593  map<string, TH2*> uncert2DPlots;
594  map<string, TH1*> uncertENuPlots;
595  map<string, TH1*> uncertQ2Plots;
596 
597  for (string syst : systs)
598  {
599  TH2 * uncert2D = NULL;
600  TH1 * uncertENu = NULL;
601  TH1 * uncertQ2 = NULL;
602  if(std::find(standardSysts.begin(), standardSysts.end(), syst) != standardSysts.end() || syst == "stat")
603  {
604  if (uncerts2D.find(syst) != uncerts2D.end() && uncerts2D[syst]){
605  cout << "Found 2D uncert for: " << syst << endl;
606  uncert2D = GetAbsolute2D(uncerts2D[syst]);
607  }
608  uncertENu = GetAbsolute(uncertsENu[syst]);
609  uncertQ2 = GetAbsolute(uncertsQ2[syst]);
610  }
611  else if (standardSystsUpDown.find(syst) != standardSystsUpDown.end())
612  {
613  string upSyst = standardSystsUpDown[syst].first;
614  string dwSyst = standardSystsUpDown[syst].second;
615  if (uncerts2D.find(upSyst) != uncerts2D.end() && uncerts2D.find(dwSyst) != uncerts2D.end() && uncerts2D[upSyst] && uncerts2D[dwSyst])
616  uncert2D = GetMaxUpDown2D(uncerts2D[upSyst], uncerts2D[dwSyst]);
617  uncertENu = GetMaxUpDown(uncertsENu[upSyst], uncertsENu[dwSyst]);
618  uncertQ2 = GetMaxUpDown(uncertsQ2[upSyst], uncertsQ2[dwSyst]);
619  }
620  else if (syst == "beam")
621  {
622  if (uncerts2D.find("2kA_Up") != uncerts2D.end() && uncerts2D["2kA_Up"])
623  uncert2D = GetTotalBeamUncert2D();
624  uncertENu = GetTotalBeamUncertENu();
625  uncertQ2 = GetTotalBeamUncertQ2();
626  }
627  else if (syst == "genie")
628  {
629  vector<TH1*> genieUnivs2D;
630  vector<TH1*> genieUnivsENu;
631  vector<TH1*> genieUnivsQ2;
632  for (unsigned int iuniv = 0; iuniv < N_GENIE; ++iuniv)
633  {
634  string syst = "genie" + to_string(iuniv);
635  if(xSecs2D.find(syst) != xSecs2D.end() && uncerts2D[syst])
636  genieUnivs2D. push_back(xSecs2D[syst]);
637  genieUnivsENu.push_back(xSecsENu[syst]);
638  genieUnivsQ2. push_back(xSecsQ2[syst]);
639  }
640  ana::Multiverse * genieVerse2D = NULL;
641  if (!genieUnivs2D.empty())
642  genieVerse2D = new ana::Multiverse(genieUnivs2D, 2);
643  ana::Multiverse * genieVerseENu = new ana::Multiverse(genieUnivsENu, 1);
644  ana::Multiverse * genieVerseQ2 = new ana::Multiverse(genieUnivsQ2, 1);
645  if (genieVerse2D)
646  uncert2D = GetMaxUpDown2D(
647  new TH2D(*(TH2D*)genieVerse2D->GetPlusOneSigmaShift(nominalXSecs2D["cv"]) - *nominalXSecs2D["cv"]),
648  new TH2D(*(TH2D*)genieVerse2D->GetMinusOneSigmaShift(nominalXSecs2D["cv"]) - *nominalXSecs2D["cv"]));
649  uncertENu = GetMaxUpDown(
650  new TH1D(*(TH1D*)genieVerseENu->GetPlusOneSigmaShift(nominalXSecsENu["cv"]) - *nominalXSecsENu["cv"]),
651  new TH1D(*(TH1D*)genieVerseENu->GetMinusOneSigmaShift(nominalXSecsENu["cv"]) - *nominalXSecsENu["cv"]));
652  uncertQ2 = GetMaxUpDown(
653  new TH1D(*(TH1D*)genieVerseQ2->GetPlusOneSigmaShift(nominalXSecsQ2["cv"]) - *nominalXSecsQ2["cv"]),
654  new TH1D(*(TH1D*)genieVerseQ2->GetMinusOneSigmaShift(nominalXSecsQ2["cv"]) - *nominalXSecsQ2["cv"]));
655 
656  if (uncert2D)
657  uncert2D->Divide(nominalXSecs2D["cv"]);
658  uncertENu->Divide(nominalXSecsENu["cv"]);
659  uncertQ2->Divide(nominalXSecsQ2["cv"]);
660  }
661  else if (syst == "ppfx")
662  {
663  vector<TH1*> ppfxUnivs2D;
664  vector<TH1*> ppfxUnivsENu;
665  vector<TH1*> ppfxUnivsQ2;
666  for (unsigned int iuniv = 0; iuniv < N_PPFX; ++iuniv)
667  {
668  string syst = "ppfx" + to_string(iuniv);
669  if(xSecs2D.find(syst) != xSecs2D.end() && uncerts2D[syst])
670  ppfxUnivs2D. push_back(xSecs2D[syst]);
671  ppfxUnivsENu.push_back(xSecsENu[syst]);
672  ppfxUnivsQ2. push_back(xSecsQ2[syst]);
673  }
674  ana::Multiverse * ppfxVerse2D = NULL;
675  if (!ppfxUnivs2D.empty())
676  ppfxVerse2D = new ana::Multiverse(ppfxUnivs2D, 2);
677  ana::Multiverse * ppfxVerseENu = new ana::Multiverse(ppfxUnivsENu, 1);
678  ana::Multiverse * ppfxVerseQ2 = new ana::Multiverse(ppfxUnivsQ2, 1);
679  if (ppfxVerse2D)
680  uncert2D = GetMaxUpDown2D(
681  new TH2D(*(TH2D*)ppfxVerse2D->GetPlusOneSigmaShift(nominalXSecs2D["cv"]) - *nominalXSecs2D["cv"]),
682  new TH2D(*(TH2D*)ppfxVerse2D->GetMinusOneSigmaShift(nominalXSecs2D["cv"]) - *nominalXSecs2D["cv"]));
683  uncertENu = GetMaxUpDown(
684  new TH1D(*(TH1D*)ppfxVerseENu->GetPlusOneSigmaShift(nominalXSecsENu["cv"]) - *nominalXSecsENu["cv"]),
685  new TH1D(*(TH1D*)ppfxVerseENu->GetMinusOneSigmaShift(nominalXSecsENu["cv"]) - *nominalXSecsENu["cv"]));
686  uncertQ2 = GetMaxUpDown(
687  new TH1D(*(TH1D*)ppfxVerseQ2->GetPlusOneSigmaShift(nominalXSecsQ2["cv"]) - *nominalXSecsQ2["cv"]),
688  new TH1D(*(TH1D*)ppfxVerseQ2->GetMinusOneSigmaShift(nominalXSecsQ2["cv"]) - *nominalXSecsQ2["cv"]));
689 
690  if (uncert2D)
691  uncert2D->Divide(nominalXSecs2D["cv"]);
692  uncertENu->Divide(nominalXSecsENu["cv"]);
693  uncertQ2->Divide(nominalXSecsQ2["cv"]);
694  }
695  if (totalUncert2D && uncert2D)
696  AddQuadrature2D(totalUncert2D, uncert2D);
697  AddQuadrature(totalUncertENu, uncertENu);
698  AddQuadrature(totalUncertQ2, uncertQ2);
699  if (uncert2D) uncert2D ->SetLineColor(systColors[syst]);
700  if (uncertENu) uncertENu->SetLineColor(systColors[syst]);
701  if (uncertQ2) uncertQ2 ->SetLineColor(systColors[syst]);
702  uncert2DPlots[syst] = uncert2D;
703  uncertENuPlots[syst] = uncertENu;
704  uncertQ2Plots[syst] = uncertQ2;
705  }
706 
707  // Plotting
708  TCanvas * systCanv = new TCanvas();
709  systCanv->SetLeftMargin(0.12);
710  systCanv->SetBottomMargin(0.12);
711  systCanv->SaveAs((systFilename + "[").c_str());
712  TLegend * leg = new TLegend();
713 
714  // 2D Slices
715  if (totalUncert2D){
716  for (int ibinx = 2; ibinx <= totalUncert2D->GetNbinsX(); ++ibinx)
717  {
718  TH1 * totalUncertSlice = GetSliceHist(totalUncert2D, ibinx);
719  LimitTMuBins(totalUncertSlice, ibinx);
720  totalUncertSlice->SetTitle(("Uncertainties: " +
721  to_string(totalUncert2D->GetXaxis()->GetBinLowEdge(ibinx)) + " < cos(#theta_{#mu}) < " +
722  to_string(totalUncert2D->GetXaxis()->GetBinUpEdge(ibinx))).c_str());
723  totalUncertSlice->GetYaxis()->SetRangeUser(0.0, 0.38);
724  totalUncertSlice->Draw("hist");
725  vector<pair<TH1*, string> > histsForLeg;
726  for (string syst : systs){
727  TH1 * histSlice = GetSliceHist(uncert2DPlots[syst], ibinx);
728  LimitTMuBins(histSlice, ibinx);
729  histSlice->SetLineColor(systColors[syst]);
730  histSlice->Draw("hist same");
731  histsForLeg.push_back({histSlice, systTitles[syst]});
732  }
733  leg = ana::AutoPlaceLegend(0.2, 0.35);
734  leg->AddEntry(totalUncertENu, "Total Uncert");
735  for (pair<TH1*, string> hists : histsForLeg)
736  leg->AddEntry(hists.first, hists.second.c_str());
737  leg->Draw();
738  systCanv->SaveAs(systFilename.c_str());
739  systCanv->SaveAs(("plots/SliceSyst_" + to_string(ibinx) + ".png").c_str());
740  leg->Clear();
741  delete leg;
742  }
743  }
744  else{
745  // ENu Plots
746  totalUncertENu->SetTitle("Uncertainties - ENu");
747  totalUncertENu->GetYaxis()->SetTitle("#frac{#delta #sigma}{#sigma}");
748  totalUncertENu->GetYaxis()->CenterTitle();
749  totalUncertENu->Draw("hist");
750  for (string syst : systs){
751  TH1 * hist = uncertENuPlots[syst];
752  hist->SetLineColor(systColors[syst]);
753  hist->Draw("hist same");
754  }
755  leg = ana::AutoPlaceLegend(0.2, 0.35);
756  leg->AddEntry(totalUncertENu, "Total Uncert");
757  for (string syst : systs)
758  leg->AddEntry(uncertENuPlots[syst], systTitles[syst].c_str());
759  leg->Draw();
760  systCanv->SaveAs(systFilename.c_str());
761  leg->Clear();
762  delete leg;
763 
764  // Q2 Plots
765  totalUncertQ2->SetTitle("Uncertainties - Q2");
766  totalUncertQ2->GetYaxis()->SetTitle("#frac{#delta #sigma}{#sigma}");
767  totalUncertQ2->GetYaxis()->CenterTitle();
768  totalUncertQ2->Draw("hist");
769  for (string syst : systs){
770  TH1 * hist = uncertQ2Plots[syst];
771  hist->SetLineColor(systColors[syst]);
772  hist->Draw("hist same");
773  }
774  leg = ana::AutoPlaceLegend(0.2, 0.35);
775  leg->AddEntry(totalUncertQ2, "Total Uncert");
776  for (string syst : systs)
777  leg->AddEntry(uncertQ2Plots[syst], systTitles[syst].c_str());
778  leg->Draw();
779  systCanv->SaveAs(systFilename.c_str());
780  leg->Clear();
781  delete leg;
782  }
783 
784  systCanv->SaveAs((systFilename + "]").c_str());
785 
786  systCanv->Clear();
787  delete systCanv;
788 }
789 
790 
791 ////////////////////////////////
792 // Main
793 ////////////////////////////////
794 
795 
796 void PlotXSec(string filename, string outputFolder)
797 {
798  inputFile = new TFile(filename.c_str(), "READ");
799  outFile = new TFile((outputFolder + "/xsecUncertResults.root").c_str(), "RECREATE");
800 
801  c1 = new TCanvas();
802  string outputPDF = outputFolder + "/allUncertainties.pdf";
803  c1->SetTopMargin(0.15);
804  c1->SetBottomMargin(0.15);
805  c1->SetRightMargin(0.15);
806  c1->SaveAs((outputPDF + "[").c_str());
807 
808  string output1D = outputFolder + "/1DXSecs.pdf";
809 
810  string doubleDiffTitle = "#frac{d^{2}#sigma}{dcos(#theta)dTmu} (cm^{2} / GeV / nucleon)";
811  string uncertTitle = "#frac{#delta#sigma}{#sigma}";
812 
813  LoadNominal("cv");
814  for(pair<string, string> nonNominal : otherNomSysts)
815  {
816  string otherNom = nonNominal.second;
817  if (otherNomSysts.find(otherNom) == otherNomSysts.end())
818  LoadNominal(otherNom); // Don't repeat for the same systematic
819  }
820  if (nominalXSecs2D["cv"]){
821  nominalXSecs2D["cv"]->Draw("COLZ");
822  c1->SaveAs(outputPDF.c_str());
823  }
824  LoadStatUncert("cv");
825 
826  // Standard XSecs
827  for (string syst : standardSysts)
828  {
829  LoadSyst(syst);
830  LoadUncerts(syst);
831  }
832 
833  // Standard XSecs
834  for (pair<string, pair<string, string> > systDef : standardSystsUpDown)
835  {
836  string systBase = systDef.first;
837  string systUp = systDef.second.first;
838  string systDw = systDef.second.second;
839 
840  LoadSyst(systUp);
841  LoadUncerts(systUp);
842 
843  LoadSyst(systDw);
844  LoadUncerts(systDw);
845  }
846 
847  // Beam XSecs
848  for (string beamSyst : beamSysts)
849  {
850  string systUp = beamSyst + "_Up";
851  string systDw = beamSyst + "_Dw";
852 
853  LoadSyst(systUp);
854  LoadUncerts(systUp);
855 
856  LoadSyst(systDw);
857  LoadUncerts(systDw);
858  }
859 
860  // Handle Genie and PPFX
861  for (unsigned int igenie = 0; igenie < N_GENIE; ++igenie)
862  {
863  string syst = "genie" + to_string(igenie);
864  LoadSyst(syst);
865  LoadUncerts(syst);
866  }
867  for (unsigned int ippfx = 0; ippfx < N_PPFX; ++ippfx)
868  {
869  string syst = "ppfx" + to_string(ippfx);
870  LoadSyst(syst);
871  LoadUncerts(syst);
872  }
873 
874  c1->SaveAs((outputPDF + "]").c_str());
875 
877  {"muone", "angle", "neutron", "genie", "ppfx", "beam", "stat"},
878  outputFolder + "/systematicPlots.pdf");
879  c1->Clear();
880 
881  nominalXSecsENu["cv"]->GetYaxis()->SetRangeUser(0.0, nominalXSecsENu["cv"]->GetMaximum() * 1.5);
882  nominalXSecsENu["cv"]->SetLineColor(kRed);
883  nominalXSecsENu["cv"]->Draw("hist");
884  for (unsigned int igenie = 0; igenie < N_GENIE; ++igenie){
885  TH1 * hist = xSecsENu["genie" + to_string(igenie)];
886  hist->SetLineColor(kBlue);
887  hist->Draw("hist same");
888  }
889  nominalXSecsENu["cv"]->Draw("hist same");
890  c1->SaveAs((outputFolder + "/enuGenie.png").c_str());
891 
892  nominalXSecsQ2["cv"]->GetYaxis()->SetRangeUser(0.0, nominalXSecsQ2["cv"]->GetMaximum() * 1.5);
893  nominalXSecsQ2["cv"]->SetLineColor(kRed);
894  nominalXSecsQ2["cv"]->Draw("hist");
895  for (unsigned int ippfx = 0; ippfx < N_PPFX; ++ippfx){
896  TH1 * hist = xSecsQ2["ppfx" + to_string(ippfx)];
897  hist->SetLineColor(kBlue);
898  hist->Draw("hist same");
899  }
900  nominalXSecsQ2["cv"]->Draw("hist same");
901  c1->SaveAs((outputFolder + "/q2Genie.png").c_str());
902 
903  outFile->Close();
904 }
TH3D * Load3DXSec(string syst)
Definition: PlotXSec.C:206
TH2 * GetMaxUpDown2D(TH2 *upHist, TH2 *dwHist)
Definition: PlotXSec.C:482
TH2D * XSec2D(string syst)
Definition: PlotXSec.C:221
enum BeamMode kOrange
TFile * inputFile
Definition: PlotXSec.C:134
map< string, TH1D * > uncertsENu
Definition: PlotXSec.C:130
enum BeamMode kRed
std::vector< SystGroupDef > systs
Definition: syst_header.h:385
Spectrum * GetPlusOneSigmaShift(const Spectrum *)
Definition: Multiverse.cxx:443
map< string, TH1D * > nominalXSecsQ2
Definition: PlotXSec.C:122
void MakeSystematicPlots(vector< string > systs, string systFilename)
Definition: PlotXSec.C:579
TH2D * XSec2DUncert(string syst)
Definition: PlotXSec.C:272
TH1 * GetMaxUpDown(TH1 *upHist, TH1 *dwHist)
Definition: PlotXSec.C:473
vector< string > standardSysts
Definition: PlotXSec.C:41
T sqrt(T number)
Definition: d0nt_math.hpp:156
constexpr T pow(T x)
Definition: pow.h:75
TH1D * XSecQ2Uncert(string syst)
Definition: PlotXSec.C:324
TH2 * GetAbsolute2D(TH2 *hist)
Definition: PlotXSec.C:502
map< string, TH3D * > all3DXSecs
Definition: PlotXSec.C:125
map< string, TH1D * > xSecsENu
Definition: PlotXSec.C:127
TH2D * nom2DXSec
Definition: PlotXSec.C:115
map< int, double > maxTMuBins
Definition: PlotXSec.C:94
void abs(TH1 *hist)
string filename
Definition: shutoffs.py:106
const vector< string > vars3D
Definition: PlotXSec.C:109
TH2 * GetTotalBeamUncert2D()
Definition: PlotXSec.C:523
TH1 * GetAbsolute(TH1 *hist)
Definition: PlotXSec.C:493
void LimitTMuBins(TH1 *hist, int ibiny)
Definition: PlotXSec.C:566
TString hists[nhists]
Definition: bdt_com.C:3
void LoadSyst(string syst)
Definition: PlotXSec.C:371
TH1 * GetTotalBeamUncertQ2()
Definition: PlotXSec.C:551
Spectrum * GetMinusOneSigmaShift(const Spectrum *)
Definition: Multiverse.cxx:450
const TAxis uncertAxis(1,-0.2, 0.2)
void LoadNominal(string syst)
Definition: PlotXSec.C:348
TH1 * GetTotalBeamUncertENu()
Definition: PlotXSec.C:537
map< string, string > systTitles
Definition: PlotXSec.C:23
TCanvas * c1
Definition: PlotXSec.C:138
void AddUpDwUncerts2D(TH2 *destHist, TH2 *upHist, TH2 *dwHist)
Definition: PlotXSec.C:190
TFile * outFile
Definition: PlotXSec.C:135
TH1D * XSecENu(string syst)
Definition: PlotXSec.C:244
map< string, TH2D * > uncerts2D
Definition: PlotXSec.C:129
const unsigned int N_PPFX
Definition: PlotXSec.C:149
const TAxis angleAxis(1, 0.5, 1.0)
base_types push_back(int_type())
vector< string > beamSysts
Definition: PlotXSec.C:65
const unsigned int N_GENIE
Definition: PlotXSec.C:148
map< string, Color_t > systColors
Definition: PlotXSec.C:80
void AddUpDwUncerts(TH1 *destHist, TH1 *upHist, TH1 *dwHist)
Definition: PlotXSec.C:177
const TAxis enuAxis(1, 0.0, 5.0)
map< string, TH1D * > uncertsQ2
Definition: PlotXSec.C:131
OStream cout
Definition: OStream.cxx:6
const string var3DName
Definition: PlotXSec.C:106
map< string, TH2D * > nominalXSecs2D
Definition: PlotXSec.C:120
void PlotXSec(string filename, string outputFolder)
Definition: PlotXSec.C:796
const map< string, string > otherNomSysts
Definition: PlotXSec.C:57
void AddQuadrature(TH1 *dest, TH1 *hist)
Definition: PlotXSec.C:437
TDirectory * dir
Definition: macro.C:5
TH1 * GetSliceHist(TH2 *hist, int ibinx)
Definition: PlotXSec.C:512
void FixAxisRange2D(TH2 *hist)
Definition: PlotXSec.C:168
string defaultNominal
Definition: PlotXSec.C:118
map< string, TH2D * > xSecs2D
Definition: PlotXSec.C:126
enum BeamMode kViolet
void LoadUncerts(string syst)
Definition: PlotXSec.C:385
map< string, TH3D * > nominalXSecs3D
Definition: PlotXSec.C:119
assert(nhit_max >=nhit_nbins)
map< string, TH1D * > xSecsQ2
Definition: PlotXSec.C:128
const TAxis muoneAxis(1, 0.5, 2.5)
const TAxis q2Axis(1, 0.0, 2.8)
TLegend * AutoPlaceLegend(double dx, double dy, double yPin)
Create a legend, maximizing distance from all histograms.
Definition: Plots.cxx:715
std::string to_string(ModuleType mt)
Definition: ModuleType.h:32
void FixAxisRangeENu(TH1 *hist)
Definition: PlotXSec.C:159
T max(sqlite3 *const db, std::string const &table_name, std::string const &column_name)
Definition: statistics.h:68
void LoadStatUncert(string nom)
Definition: PlotXSec.C:398
enum BeamMode kGreen
enum BeamMode kBlue
TH1D * XSecENuUncert(string syst)
Definition: PlotXSec.C:300
TH1D * XSecQ2(string syst)
Definition: PlotXSec.C:258
void FixAxisRangeQ2(TH1 *hist)
Definition: PlotXSec.C:151
map< string, pair< string, string > > standardSystsUpDown
Definition: PlotXSec.C:48
std::string UniqueName()
Return a different string each time, for creating histograms.
Definition: Utilities.cxx:29
TH3D * nom3DXSec
Definition: PlotXSec.C:114
void AddQuadrature2D(TH2 *dest, TH2 *hist)
Definition: PlotXSec.C:451
string uncertLabel
Definition: PlotXSec.C:107
map< string, TH1D * > nominalXSecsENu
Definition: PlotXSec.C:121
enum BeamMode string