SystsCalibRelAna17.C
Go to the documentation of this file.
4 
5 using namespace ana;
6 
7 void PlotSysts(TH1* nom, TH1* hfu, TH1* hfd, TH1* hnu, TH1* hnd,
8  TDirectory* rootOUT, FILE* textOFS,
9  strings strs);
10 
12 {
13  TH1::AddDirectory(0);
14  double NCSCALE =8.85e20;
15  // A map of systematic labels
16  std::map<std::string, std::string> shiftlabels;
17  shiftlabels["CalFlat095ND"] = "Flat Miscalibration Down at ND";
18  shiftlabels["CalFlat105ND"] = "Flat Miscalibration Up at ND";
19  shiftlabels["CalFlat095FD"] = "Flat Miscalibration Down at FD";
20  shiftlabels["CalFlat105FD"] = "Flat Miscalibration Up at FD";
21 
22  std::string labelRecoE = "Energy Deposited in Scintillator (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  std::string folder1="./Output_Ana/";
36  std::string folder = "/nova/ana/users/sedayath/CAFAna_work_place/job/NUS17/Syst/";
37  // std::string filenm = "Cal_Rel";
38  std::string filenm = "Cal_Rel1";
39  std::string loadLocation = folder + filenm + ".root";
40  std::string saveLocation = folder1 + filenm + "Ana.root";
41  std::string textLocation = folder1 + 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";
57  strs.fPOT = StringFromDouble(NCSCALE);
58  strs.fSample = sample;
59  strs.fSystType = "CalibRel";
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  strs.fSystS = shifts.first;
75  strs.fSystL = shifts.second;
76 
77  strs.fComponent = "NC";
78  strs.fDet = "ND";
79  PlotSyst(hNomNC_decomp,
80  decomps_shifted[sample][shifts.first][1],
81  rootF, textF, strs, NCSCALE, true, true);
82 
83  strs.fComponent = "BG";
84  PlotSyst(hNomBG_decomp,
85  decomps_shifted[sample][shifts.first][1],
86  rootF, textF, strs, NCSCALE, false, true);
87 
88  strs.fComponent = "NC";
89  strs.fDet = "FD";
90  PlotSyst(hNomNC_predNE,
91  predsNE_shifted[sample][shifts.first][1],
92  rootF, textF, strs, NCSCALE, true, true);
93 
94  strs.fComponent = "BG";
95  PlotSyst(hNomBG_predNE,
96  predsNE_shifted[sample][shifts.first][1],
97  rootF, textF, strs, NCSCALE, false, true);
98 
99  strs.fComponent = "NC";
100  strs.fDet = "EX";
101  PlotSyst(hNomNC_predSt,
102  predsSt_shifted[sample][shifts.first][1],
103  rootF, textF, strs, NCSCALE, true, true);
104 
105  strs.fComponent = "BG";
106  PlotSyst(hNomBG_predSt,
107  predsSt_shifted[sample][shifts.first][1],
108  rootF, textF, strs, NCSCALE, false, true);
109  }
110 
111  TH1* hFDUNC_decomp = GetNC(decomps_shifted[sample]["CalFlat105FD"][1], NCSCALE);
112  TH1* hFDUBG_decomp = GetBG(decomps_shifted[sample]["CalFlat105FD"][1], NCSCALE);
113  TH1* hFDDNC_decomp = GetNC(decomps_shifted[sample]["CalFlat095FD"][1], NCSCALE);
114  TH1* hFDDBG_decomp = GetBG(decomps_shifted[sample]["CalFlat095FD"][1], NCSCALE);
115  TH1* hNDUNC_decomp = GetNC(decomps_shifted[sample]["CalFlat105ND"][1], NCSCALE);
116  TH1* hNDUBG_decomp = GetBG(decomps_shifted[sample]["CalFlat105ND"][1], NCSCALE);
117  TH1* hNDDNC_decomp = GetNC(decomps_shifted[sample]["CalFlat095ND"][1], NCSCALE);
118  TH1* hNDDBG_decomp = GetBG(decomps_shifted[sample]["CalFlat095ND"][1], NCSCALE);
119 
120  TH1* hFDUNC_predNE = GetNC(predsNE_shifted[sample]["CalFlat105FD"][1], NCSCALE);
121  TH1* hFDUBG_predNE = GetBG(predsNE_shifted[sample]["CalFlat105FD"][1], NCSCALE);
122  TH1* hFDDNC_predNE = GetNC(predsNE_shifted[sample]["CalFlat095FD"][1], NCSCALE);
123  TH1* hFDDBG_predNE = GetBG(predsNE_shifted[sample]["CalFlat095FD"][1], NCSCALE);
124  TH1* hNDUNC_predNE = GetNC(predsNE_shifted[sample]["CalFlat105ND"][1], NCSCALE);
125  TH1* hNDUBG_predNE = GetBG(predsNE_shifted[sample]["CalFlat105ND"][1], NCSCALE);
126  TH1* hNDDNC_predNE = GetNC(predsNE_shifted[sample]["CalFlat095ND"][1], NCSCALE);
127  TH1* hNDDBG_predNE = GetBG(predsNE_shifted[sample]["CalFlat095ND"][1], NCSCALE);
128 
129  TH1* hFDUNC_predSt = GetNC(predsSt_shifted[sample]["CalFlat105FD"][1], NCSCALE);
130  TH1* hFDUBG_predSt = GetBG(predsSt_shifted[sample]["CalFlat105FD"][1], NCSCALE);
131  TH1* hFDDNC_predSt = GetNC(predsSt_shifted[sample]["CalFlat095FD"][1], NCSCALE);
132  TH1* hFDDBG_predSt = GetBG(predsSt_shifted[sample]["CalFlat095FD"][1], NCSCALE);
133  TH1* hNDUNC_predSt = GetNC(predsSt_shifted[sample]["CalFlat105ND"][1], NCSCALE);
134  TH1* hNDUBG_predSt = GetBG(predsSt_shifted[sample]["CalFlat105ND"][1], NCSCALE);
135  TH1* hNDDNC_predSt = GetNC(predsSt_shifted[sample]["CalFlat095ND"][1], NCSCALE);
136  TH1* hNDDBG_predSt = GetBG(predsSt_shifted[sample]["CalFlat095ND"][1], NCSCALE);
137 
138  strs.fComponent = "NC";
139  strs.fDet = "ND";
140  PlotSysts(hNomNC_decomp,
141  hFDUNC_decomp, hFDDNC_decomp, hNDUNC_decomp, hNDDNC_decomp,
142  rootF, textF, strs);
143 
144  strs.fComponent = "BG";
145  PlotSysts(hNomBG_decomp,
146  hFDUBG_decomp, hFDDBG_decomp, hNDUBG_decomp, hNDDBG_decomp,
147  rootF, textF, strs);
148 
149  strs.fComponent = "NC";
150  strs.fDet = "FD";
151  PlotSysts(hNomNC_predNE,
152  hFDUNC_predNE, hFDDNC_predNE, hNDUNC_predNE, hNDDNC_predNE,
153  rootF, textF, strs);
154 
155  strs.fComponent = "BG";
156  PlotSysts(hNomBG_predNE,
157  hFDUBG_predNE, hFDDBG_predNE, hNDUBG_predNE, hNDDBG_predNE,
158  rootF, textF, strs);
159 
160  strs.fComponent = "NC";
161  strs.fDet = "EX";
162  PlotSysts(hNomNC_predSt,
163  hFDUNC_predSt, hFDDNC_predSt, hNDUNC_predSt, hNDDNC_predSt,
164  rootF, textF, strs);
165 
166  strs.fComponent = "BG";
167  PlotSysts(hNomBG_predSt,
168  hFDUBG_predSt, hFDDBG_predSt, hNDUBG_predSt, hNDDBG_predSt,
169  rootF, textF, strs);
170  }
171 
172  fclose(textF);
173  rootF->Close();
174 
175  std::cout << "Finished. Ctrl+Z should be safe if the macro didn't quit." << std::endl;
176 }
177 
178 void PlotSysts(TH1* nom, TH1* hfu, TH1* hfd, TH1* hnu, TH1* hnd,
179  TDirectory* rootOUT, FILE* textOFS,
180  strings strs)
181 {
182  double nNom = nom->Integral();
183  double nFDU = hfu->Integral();
184  double nFDD = hfd->Integral();
185  double nNDU = hnu->Integral();
186  double nNDD = hnd->Integral();
187 
188  ZeroError(nom);
189  ZeroError(hfu);
190  ZeroError(hfd);
191  ZeroError(hnu);
192  ZeroError(hnd);
193 
194  fprintf(textOFS, "\nEvent counts for %.2s %.2s %s Systs:\n",
195  strs.fDet.c_str(), strs.fComponent.c_str(), strs.fSystType.c_str());
196  fprintf(textOFS, "NC: %6.2f, FD Up: %6.2f, FD Down: %6.2f, ND Up: %6.2f, ND Down: %6.2f\n",
197  nNom, nFDU, nFDD, nNDU, nNDD);
198  fprintf(textOFS, "FD Up Diff: %6.2f, FD Down Diff: %6.2f, ND Up Diff: %6.2f, ND Down Diff: %6.2f",
199  100.*std::abs(nFDU/nNom - 1.), 100.*std::abs(nFDD/nNom - 1.),
200  100.*std::abs(nNDU/nNom - 1.), 100.*std::abs(nNDD/nNom - 1.));
201 
202  double diffs[] = {std::abs(nFDU - nNom),
203  std::abs(nFDD - nNom),
204  std::abs(nNDU - nNom),
205  std::abs(nNDD - nNom)};
206  if( *std::max_element(diffs, diffs+4) == diffs[0]) {
207  AddErrorInQuadrature(nom, hfu, true);
208  }
209  else if(*std::max_element(diffs, diffs+4) == diffs[1]){
210  AddErrorInQuadrature(nom, hfd, true);
211  }
212  else if(*std::max_element(diffs, diffs+4) == diffs[2]){
213  AddErrorInQuadrature(nom, hnu, true);
214  }
215  else {
216  AddErrorInQuadrature(nom, hnd, true);
217  }
218 
219  std::string xlabel = strs.fXLabel;
220  std::string ylabel = "Events / " + strs.fPOT + " / 0.25 GeV";
221  std::string title = strs.fDet + " " + strs.fComponent +
222  " Errors for All " + strs.fSystType + " Systematics";
223 // std::string fullTitle = title + ";" + xlabel + ";" + ylabel;
224  std::string fullTitle = ";" + xlabel + ";" + ylabel;
225  std::string name = "c" + strs.fComponent + strs.fDet +
226  strs.fSystType + "SystsAll" + strs.fSample;
227  std::string ylabelRat = "Shifted / Nominal";
228  std::string fullTitleRat = ";" + xlabel + ";" + ylabelRat;
229 
230  TH1* hRat = (TH1*)nom->Clone();
231  TH1* hFDURat = (TH1*)hfu->Clone();
232  TH1* hFDDRat = (TH1*)hfd->Clone();
233  TH1* hNDURat = (TH1*)hnu->Clone();
234  TH1* hNDDRat = (TH1*)hnd->Clone();
235 
236  int nBins = nom->GetNbinsX();
237  for(int i = 0; i <= nBins+1; ++i) {
238  hRat->SetBinContent(i, 1.);
239  }
240  hFDURat->Divide(nom);
241  hFDDRat->Divide(nom);
242  hNDURat->Divide(nom);
243  hNDDRat->Divide(nom);
244 
245  double maxval = nom->GetBinContent(nom->GetMaximumBin());
246  maxval = std::max(maxval, hfu->GetBinContent(hfu->GetMaximumBin()));
247  maxval = std::max(maxval, hfd->GetBinContent(hfd->GetMaximumBin()));
248  maxval = std::max(maxval, hnu->GetBinContent(hnu->GetMaximumBin()));
249  maxval = std::max(maxval, hnd->GetBinContent(hnd->GetMaximumBin()));
250 
251  double minval = nom->GetBinContent(nom->GetMinimumBin());
252  minval = std::min(minval, hfu->GetBinContent(hfu->GetMinimumBin()));
253  minval = std::min(minval, hfd->GetBinContent(hfd->GetMinimumBin()));
254  minval = std::min(minval, hnu->GetBinContent(hnu->GetMinimumBin()));
255  minval = std::min(minval, hnd->GetBinContent(hnd->GetMinimumBin()));
256 
257  if(maxval < minval) { std::swap(maxval, minval); }
258 
259  double maxvalRat = hRat->GetBinContent(hRat->GetMaximumBin());
260  maxvalRat = std::max(maxvalRat, hFDURat->GetBinContent(hFDURat->GetMaximumBin()));
261  maxvalRat = std::max(maxvalRat, hFDDRat->GetBinContent(hFDDRat->GetMaximumBin()));
262  maxvalRat = std::max(maxvalRat, hNDURat->GetBinContent(hNDURat->GetMaximumBin()));
263  maxvalRat = std::max(maxvalRat, hNDDRat->GetBinContent(hNDDRat->GetMaximumBin()));
264 
265  double minvalRat = hRat->GetBinContent(hRat->GetMinimumBin());
266  for(int i = 1; i <= nBins; ++i) {
267  if(hFDURat->GetBinContent(i) > 0.) {
268  minvalRat = std::min(minvalRat, hFDURat->GetBinContent(i));
269  }
270  if(hFDDRat->GetBinContent(i) > 0.) {
271  minvalRat = std::min(minvalRat, hFDDRat->GetBinContent(i));
272  }
273  if(hNDURat->GetBinContent(i) > 0.) {
274  minvalRat = std::min(minvalRat, hNDURat->GetBinContent(i));
275  }
276  if(hNDDRat->GetBinContent(i) > 0.) {
277  minvalRat = std::min(minvalRat, hNDDRat->GetBinContent(i));
278  }
279  }
280 
281  if(std::abs(maxvalRat - 1.) > std::abs(1. - minvalRat)) {
282  minvalRat = 1 - std::abs(maxvalRat - 1.);
283  maxvalRat = 1 + std::abs(maxvalRat - 1.);
284  }
285  else {
286  minvalRat = 1 - std::abs(1. - minvalRat);
287  maxvalRat = 1 + std::abs(1. - minvalRat);
288  }
289 
290  CenterTitles(nom);
291  nom->SetLineColor(kBlack);
292  nom->SetMaximum(1.1*maxval);
293  nom->SetMinimum(0);
294  nom->SetTitle(fullTitle.c_str());
295 
296  CenterTitles(hfu);
297  hfu->SetLineColor(kRed);
298  hfu->SetLineStyle(2);
299  hfu->SetMaximum(1.1*maxval);
300  hfu->SetMinimum(0);
301  hfu->SetTitle(fullTitle.c_str());
302 
303  CenterTitles(hfd);
304  hfd->SetLineColor(kBlue);
305  hfd->SetLineStyle(2);
306  hfd->SetMaximum(1.1*maxval);
307  hfd->SetMinimum(0);
308  hfd->SetTitle(fullTitle.c_str());
309 
310  CenterTitles(hnu);
311  hnu->SetLineColor(kGreen+2);
312  hnu->SetLineStyle(2);
313  hnu->SetMaximum(1.1*maxval);
314  hnu->SetMinimum(0);
315  hnu->SetTitle(fullTitle.c_str());
316 
317  CenterTitles(hnd);
318  hnd->SetLineColor(kViolet-5);
319  hnd->SetLineStyle(2);
320  hnd->SetMaximum(1.1*maxval);
321  hnd->SetMinimum(0);
322  hnd->SetTitle(fullTitle.c_str());
323 
324  CenterTitles(hRat);
325  hRat->SetLineColor(kBlack);
326  hRat->SetMaximum(1.1*maxvalRat);
327  hRat->SetMinimum(0.9*minvalRat);
328  hRat->SetTitle(fullTitleRat.c_str());
329 
330  CenterTitles(hFDURat);
331  hFDURat->SetLineColor(kRed);
332  hFDURat->SetLineStyle(2);
333  hFDURat->SetMaximum(1.1*maxvalRat);
334  hFDURat->SetMinimum(0.9*minvalRat);
335  hFDURat->SetTitle(fullTitleRat.c_str());
336 
337  CenterTitles(hFDDRat);
338  hFDDRat->SetLineColor(kBlue);
339  hFDDRat->SetLineStyle(2);
340  hFDDRat->SetMaximum(1.1*maxvalRat);
341  hFDDRat->SetMinimum(0.9*minvalRat);
342  hFDDRat->SetTitle(fullTitleRat.c_str());
343 
344  CenterTitles(hNDURat);
345  hNDURat->SetLineColor(kGreen+2);
346  hNDURat->SetLineStyle(2);
347  hNDURat->SetMaximum(1.1*maxvalRat);
348  hNDURat->SetMinimum(0.9*minvalRat);
349  hNDURat->SetTitle(fullTitleRat.c_str());
350 
351  CenterTitles(hNDDRat);
352  hNDDRat->SetLineColor(kViolet-5);
353  hNDDRat->SetLineStyle(2);
354  hNDDRat->SetMaximum(1.1*maxvalRat);
355  hNDDRat->SetMinimum(0.9*minvalRat);
356  hNDDRat->SetTitle(fullTitleRat.c_str());
357 
358  double xL = 0.6, xR = 0.85;
359  double yB = 0.6, yT = 0.85;
360  TLegend* leg = new TLegend(xL, yB, xR, yT);
361  leg->AddEntry(nom, "Nominal", "l");
362  leg->AddEntry(hfu, "FD Scale Up", "l");
363  leg->AddEntry(hfd, "FD Scale Down", "l");
364  leg->AddEntry(hnu, "ND Scale Up", "l");
365  leg->AddEntry(hnd, "ND Scale Down", "l");
366 
367  TCanvas* c = new TCanvas(name.c_str(), title.c_str(), 800, 800);
368  TPad* pSpecs = new TPad("pSpecs", "", 0., 0.375, 1., 1.);
369  TPad* pRatio = new TPad("pRatio", "", 0., 0., 1., 0.375);
370  pSpecs->Draw();
371  pRatio->Draw();
372 
373  pSpecs->cd();
374  nom->Draw("hist");
375  hfu->Draw("hist same");
376  hfd->Draw("hist same");
377  hnu->Draw("hist same");
378  hnd->Draw("hist same");
379  leg->Draw();
380  Simulation();
381 
382  pRatio->cd();
383  hRat->Draw("hist");
384  hFDURat->Draw("hist same");
385  hFDDRat->Draw("hist same");
386  hNDURat->Draw("hist same");
387  hNDDRat->Draw("hist same");
388 
389  c->Update();
390  rootOUT->WriteTObject(c);
391  return;
392 }
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 ZeroError(TH1 *nom)
Set the errors of the input histogram to 0.
Definition: PPFXHelper.h:149
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 AddErrorInQuadrature(TH1 *nom, TH1 *up, bool use_max)
Definition: PPFXHelper.h:163
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
void SystsCalibRelAna17()
std::string fPOT
Definition: PPFXHelper.h:105
OStream cout
Definition: OStream.cxx:6
std::string fXLabel
Definition: PPFXHelper.h:110
enum BeamMode kViolet
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 kGreen
enum BeamMode kBlue
void PlotSysts(TH1 *nom, TH1 *hfu, TH1 *hfd, TH1 *hnu, TH1 *hnd, TDirectory *rootOUT, FILE *textOFS, strings strs)
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