plotSystBandForTheseSystematics.C
Go to the documentation of this file.
1 
2 
3 #include <cstring>
4 
5 
7 {
8  Int_t niceColours[16]={kBlue-3,kRed-3,kGreen-3,
9  kMagenta-3,kCyan-3,kOrange-3,
10  kViolet-3,kTeal-3,
11  kBlue+1,kRed-3,kGreen+1,
12  kMagenta+1,kCyan+1,kOrange+1,
13  kViolet+1,kTeal+1};
14  if(index<16) return niceColours[index];
15 
16  Int_t offset=index/16;
17  Int_t colourInd=index%16;
18 
19  return (niceColours[colourInd]+offset);
20 
21 }
22 
23 
24 TGraphAsymmErrors *getErrorBand(TH1F *histNoShift, TH1F **histPlus, TH1F **histMinus, Int_t numShifts) ;
25 
26 void plotSystBandForTheseSystematics(char *rootFileList,char *legTagList=0) {
27  std::cout << rootFileList << "\n";
28  // return;
29  std::ifstream InFileList(rootFileList);
30  char fileNames[1024][1024];
31  int countFiles=0;
32  while(InFileList >> fileNames[countFiles]) {
33  std::cout << fileNames[countFiles] << "\n";
34  countFiles++;
35  }
36  std::cout << "Got " << countFiles << " files\n";
37 
38  char legTags[1024][100];
39  for(int i=0;i<1024;i++) {
40  sprintf(legTags[i],"Syst.");
41  }
42  if(legTagList) {
43  int countTags=0;
44  std::ifstream LegTagList(legTagList);
45  while(LegTagList >> legTags[countTags]) {
46  for(int si=0;si<strlen(legTags[countTags]);si++) {
47  // std::cout << si << "\t" << legTags[countTags][si] << endl;
48  if(legTags[countTags][si]=='~') {
49  legTags[countTags][si]=' ';
50  }
51  }
52  std::cout << legTags[countTags] << "\n";//t" << index(legTags,'~') << "\n";
53  countTags++;
54  }
55  std::cout << "Got " << countTags << " files\n";
56 
57 
58  }
59 
60 
61  // KEY: TH1F FDCosmicMuonCorrected;1
62 
63 
64 
65 
66  const char*ndFullHistName="plot/MultiExperiment/NOvA_NuMuSel/CorrectedStacks/NDMCCorrected";
67  const char*fdFullHistName="plot/MultiExperiment/NOvA_NuMuSel/CorrectedStacks/FDMCCorrected";
68  // const char*ndFullHistName="plot/MultiExperiment/NOvA_NuMuSel/UncorrectedStacks/NDMCUncorrected";
69  // const char*fdFullHistName="plot/MultiExperiment/NOvA_NuMuSel/UncorrectedStacks/FDMCUncorrected";
70  const char*ndNCHistName="plot/MultiExperiment/NOvA_NuMuSel/CorrectedStacks/NDmc_ncCorrected";
71  const char*fdNCHistName="plot/MultiExperiment/NOvA_NuMuSel/CorrectedStacks/FDmc_ncCorrected";
72 
73  // const char*ndNuTauCCHistName="plot/MultiExperiment/NOvA_NuMuSel/CorrectedStacks/NDNuTauCCCorrected";
74  const char*fdNuTauCCHistName="plot/MultiExperiment/NOvA_NuMuSel/CorrectedStacks/FDmc_mtCorrected";
75  // const char*ndNuTauBarCCHistName="plot/MultiExperiment/NOvA_NuMuSel/CorrectedStacks/NDNuTauBarCCCorrected";
76  const char*fdNuTauBarCCHistName="plot/MultiExperiment/NOvA_NuMuSel/CorrectedStacks/FDmc_mtbarCorrected";
77 
78  const char*ndNuECCHistName="plot/MultiExperiment/NOvA_NuMuSel/CorrectedStacks/NDmc_eeCorrected";
79  const char*fdNuECCHistName="plot/MultiExperiment/NOvA_NuMuSel/CorrectedStacks/FDmc_eeCorrected";
80  const char*ndNuEBarCCHistName="plot/MultiExperiment/NOvA_NuMuSel/CorrectedStacks/NDmc_eebarCorrected";
81  const char*fdNuEBarCCHistName="plot/MultiExperiment/NOvA_NuMuSel/CorrectedStacks/FDmc_eebarCorrected";
82 
83  const char*fdCosmicHistName="plot/MultiExperiment/NOvA_NuMuSel/CorrectedStacks/FDcosmicCorrected";
84 
85 
86  TFile *fpNoShift = TFile::Open(fileNames[0]); //0th file is the nominal
87  TH1F *histFullNDNoShift = (TH1F*)fpNoShift->Get(ndFullHistName);
88  histFullNDNoShift->SetName("histFullNDNoShift");
89  TH1F *histFullFDNoShift = (TH1F*)fpNoShift->Get(fdFullHistName);
90  histFullFDNoShift->SetName("histFullFDNoShift");
91 
92  // TH1F *histNCNDNoShift = (TH1F*)fpNoShift->Get(ndNCHistName);
93  // histNCNDNoShift->SetName("histNCNDNoShift");
94  // TH1F *histNCFDNoShift = (TH1F*)fpNoShift->Get(fdNCHistName);
95  // histNCFDNoShift->SetName("histNCFDNoShift");
96 
97 
98  // TH1F *histNuECCNDNoShift = (TH1F*)fpNoShift->Get(ndNuECCHistName);
99  // histNuECCNDNoShift->SetName("histNuECCNDNoShift");
100  // TH1F *histNuECCFDNoShift = (TH1F*)fpNoShift->Get(fdNuECCHistName);
101  // histNuECCFDNoShift->SetName("histNuECCFDNoShift");
102 
103  // TH1F *histNuEBarCCNDNoShift = (TH1F*)fpNoShift->Get(ndNuEBarCCHistName);
104  // histNuEBarCCNDNoShift->SetName("histNuEBarCCNDNoShift");
105  // TH1F *histNuEBarCCFDNoShift = (TH1F*)fpNoShift->Get(fdNuEBarCCHistName);
106  // histNuEBarCCFDNoShift->SetName("histNuEBarCCFDNoShift");
107  // histNuECCNDNoShift->Add(histNuEBarCCNDNoShift);
108  // histNuECCFDNoShift->Add(histNuEBarCCFDNoShift);
109 
110  // TH1F *histCosmicFDNoShift = (TH1F*)fpNoShift->Get("fdCosmicHistName");
111  // histCosmicFDNoShift->SetName("histCosmicFDNoShift");
112 
113 
114  // // TH1F *histNuTauCCNDNoShift = fpNoShift->Get(ndNuTauCCHistName);
115  // // histNuTauCCNDNoShift->SetName("histNuTauCCNDNoShift");
116  // TH1F *histNuTauCCFDNoShift = (TH1F*)fpNoShift->Get(fdNuTauCCHistName);
117  // histNuTauCCFDNoShift->SetName("histNuTauCCFDNoShift");
118 
119  // // TH1F *histNuTauBarCCNDNoShift = fpNoShift->Get(ndNuTauBarCCHistName);
120  // // histNuTauBarCCNDNoShift->SetName("histNuTauBarCCNDNoShift");
121  // TH1F *histNuTauBarCCFDNoShift = (TH1F*)fpNoShift->Get(fdNuTauBarCCHistName);
122  // histNuTauBarCCFDNoShift->SetName("histNuTauBarCCFDNoShift");
123  // // histNuTauCCNDNoShift->Add(histNuTauBarCCNDNoShift);
124  // histNuTauCCFDNoShift->Add(histNuTauBarCCFDNoShift);
125 
126 
127  int numSysts=(countFiles-1)/2;
128  std::cout << "Got " << numSysts << " systematics\n";
129 
130  char *legTag="Syst";
131 
132  TH1F *histFullNDPlus1[1024];
133  TH1F *histFullNDMinus1[1024];
134  TH1F *histFullFDPlus1[1024];
135  TH1F *histFullFDMinus1[1024];
136 
137  TH1F *histNCNDPlus1[1024];
138  TH1F *histNCNDMinus1[1024];
139  TH1F *histNCFDPlus1[1024];
140  TH1F *histNCFDMinus1[1024];
141 
142  TH1F *histNuTauCCNDPlus1[1024];
143  TH1F *histNuTauCCNDMinus1[1024];
144  TH1F *histNuTauCCFDPlus1[1024];
145  TH1F *histNuTauCCFDMinus1[1024];
146 
147  TH1F *histNuTauBarCCNDPlus1[1024];
148  TH1F *histNuTauBarCCNDMinus1[1024];
149  TH1F *histNuTauBarCCFDPlus1[1024];
150  TH1F *histNuTauBarCCFDMinus1[1024];
151 
152  TH1F *histNuECCNDPlus1[1024];
153  TH1F *histNuECCNDMinus1[1024];
154  TH1F *histNuECCFDPlus1[1024];
155  TH1F *histNuECCFDMinus1[1024];
156 
157  TH1F *histNuEBarCCNDPlus1[1024];
158  TH1F *histNuEBarCCNDMinus1[1024];
159  TH1F *histNuEBarCCFDPlus1[1024];
160  TH1F *histNuEBarCCFDMinus1[1024];
161 
162  TH1F *histCosmicFDPlus1[1024];
163  TH1F *histCosmicFDMinus1[1024];
164 
165  // numSysts=5;
166  char histName[180];
167  for(int systNum=0;systNum<numSysts;systNum++) {
168 
169  TFile *fpPlus1Shift = TFile::Open(fileNames[(systNum*2)+1]);
170  histFullNDPlus1[systNum] = (TH1F*) fpPlus1Shift->Get(ndFullHistName);
171  sprintf(histName,"histFullNDPlus1_%d",systNum);
172  histFullNDPlus1[systNum]->SetName(histName);
173  histFullFDPlus1[systNum] = (TH1F*) fpPlus1Shift->Get(fdFullHistName);
174  sprintf(histName,"histFullFDPlus1_%d",systNum);
175  histFullFDPlus1[systNum]->SetName(histName);
176 
177  histNCNDPlus1[systNum] = (TH1F*) fpPlus1Shift->Get(ndNCHistName);
178  sprintf(histName,"histNCNDPlus1_%d",systNum);
179  histNCNDPlus1[systNum]->SetName(histName);
180  histNCFDPlus1[systNum] = (TH1F*) fpPlus1Shift->Get(fdNCHistName);
181  sprintf(histName,"histNCFDPlus1_%d",systNum);
182  histNCFDPlus1[systNum]->SetName(histName);
183 
184  continue;
185  histCosmicFDPlus1[systNum] = (TH1F*) fpPlus1Shift->Get(fdCosmicHistName);
186  sprintf(histName,"histCosmicFDPlus1_%d",systNum);
187  histCosmicFDPlus1[systNum]->SetName(histName);
188 
189 
190  // histNuTauCCNDPlus1[systNum] = (TH1F*) fpPlus1Shift->Get(ndNuTauCCHistName);
191  // sprintf(histName,"histNuTauCCNDPlus1_%d",systNum);
192  // histNuTauCCNDPlus1[systNum]->SetName(histName);
193  histNuTauCCFDPlus1[systNum] = (TH1F*) fpPlus1Shift->Get(fdNuTauCCHistName);
194  sprintf(histName,"histNuTauCCFDPlus1_%d",systNum);
195  histNuTauCCFDPlus1[systNum]->SetName(histName);
196 
197  // histNuTauBarCCNDPlus1[systNum] = (TH1F*) fpPlus1Shift->Get(ndNuTauBarCCHistName);
198  // sprintf(histName,"histNuTauBarCCNDPlus1_%d",systNum);
199  // histNuTauBarCCNDPlus1[systNum]->SetName(histName);
200  histNuTauBarCCFDPlus1[systNum] = (TH1F*) fpPlus1Shift->Get(fdNuTauBarCCHistName);
201  sprintf(histName,"histNuTauBarCCFDPlus1_%d",systNum);
202  histNuTauBarCCFDPlus1[systNum]->SetName(histName);
203 
204  // histNuTauCCNDPlus1[systNum]->Add(histNuTauBarCCNDPlus1[systNum]);
205  histNuTauCCFDPlus1[systNum]->Add(histNuTauBarCCFDPlus1[systNum]);
206 
207 
208  histNuECCNDPlus1[systNum] = (TH1F*) fpPlus1Shift->Get(ndNuECCHistName);
209  sprintf(histName,"histNuECCNDPlus1_%d",systNum);
210  histNuECCNDPlus1[systNum]->SetName(histName);
211  histNuECCFDPlus1[systNum] = (TH1F*) fpPlus1Shift->Get(fdNuECCHistName);
212  sprintf(histName,"histNuECCFDPlus1_%d",systNum);
213  histNuECCFDPlus1[systNum]->SetName(histName);
214 
215  histNuEBarCCNDPlus1[systNum] = (TH1F*) fpPlus1Shift->Get(ndNuEBarCCHistName);
216  sprintf(histName,"histNuEBarCCNDPlus1_%d",systNum);
217  histNuEBarCCNDPlus1[systNum]->SetName(histName);
218  histNuEBarCCFDPlus1[systNum] = (TH1F*) fpPlus1Shift->Get(fdNuEBarCCHistName);
219  sprintf(histName,"histNuEBarCCFDPlus1_%d",systNum);
220  histNuEBarCCFDPlus1[systNum]->SetName(histName);
221 
222  histNuECCNDPlus1[systNum]->Add(histNuEBarCCNDPlus1[systNum]);
223  histNuECCFDPlus1[systNum]->Add(histNuEBarCCFDPlus1[systNum]);
224 
225 
226 
227 
228  TFile *fpMinus1Shift = TFile::Open(fileNames[(systNum*2)+2]);
229  histFullNDMinus1[systNum] = (TH1F*) fpMinus1Shift->Get(ndFullHistName);
230  sprintf(histName,"histFullNDMinus1_%d",systNum);
231  histFullNDMinus1[systNum]->SetName(histName);
232  histFullFDMinus1[systNum] = (TH1F*) fpMinus1Shift->Get(fdFullHistName);
233  sprintf(histName,"histFullFDMinus1_%d",systNum);
234  histFullFDMinus1[systNum]->SetName(histName);
235 
236  histNCNDMinus1[systNum] = (TH1F*) fpMinus1Shift->Get(ndNCHistName);
237  sprintf(histName,"histNCNDMinus1_%d",systNum);
238  histNCNDMinus1[systNum]->SetName(histName);
239  histNCFDMinus1[systNum] = (TH1F*) fpMinus1Shift->Get(fdNCHistName);
240  sprintf(histName,"histNCFDMinus1_%d",systNum);
241  histNCFDMinus1[systNum]->SetName(histName);
242 
243 
244  histCosmicFDMinus1[systNum] = (TH1F*) fpMinus1Shift->Get(fdCosmicHistName);
245  sprintf(histName,"histCosmicFDMinus1_%d",systNum);
246  histCosmicFDMinus1[systNum]->SetName(histName);
247 
248 
249 
250  // histNuTauCCNDMinus1[systNum] = (TH1F*) fpMinus1Shift->Get(ndNuTauCCHistName);
251  // sprintf(histName,"histNuTauCCNDMinus1_%d",systNum);
252  // histNuTauCCNDMinus1[systNum]->SetName(histName);
253  histNuTauCCFDMinus1[systNum] = (TH1F*) fpMinus1Shift->Get(fdNuTauCCHistName);
254  sprintf(histName,"histNuTauCCFDMinus1_%d",systNum);
255  histNuTauCCFDMinus1[systNum]->SetName(histName);
256 
257  // histNuTauBarCCNDMinus1[systNum] = (TH1F*) fpMinus1Shift->Get(ndNuTauBarCCHistName);
258  // sprintf(histName,"histNuTauBarCCNDMinus1_%d",systNum);
259  // histNuTauBarCCNDMinus1[systNum]->SetName(histName);
260  histNuTauBarCCFDMinus1[systNum] = (TH1F*) fpMinus1Shift->Get(fdNuTauBarCCHistName);
261  sprintf(histName,"histNuTauBarCCFDMinus1_%d",systNum);
262  histNuTauBarCCFDMinus1[systNum]->SetName(histName);
263 
264  // histNuTauCCNDMinus1[systNum]->Add(histNuTauBarCCNDMinus1[systNum]);
265  histNuTauCCFDMinus1[systNum]->Add(histNuTauBarCCFDMinus1[systNum]);
266 
267 
268  histNuECCNDMinus1[systNum] = (TH1F*) fpMinus1Shift->Get(ndNuECCHistName);
269  sprintf(histName,"histNuECCNDMinus1_%d",systNum);
270  histNuECCNDMinus1[systNum]->SetName(histName);
271  histNuECCFDMinus1[systNum] = (TH1F*) fpMinus1Shift->Get(fdNuECCHistName);
272  sprintf(histName,"histNuECCFDMinus1_%d",systNum);
273  histNuECCFDMinus1[systNum]->SetName(histName);
274 
275  histNuEBarCCNDMinus1[systNum] = (TH1F*) fpMinus1Shift->Get(ndNuEBarCCHistName);
276  sprintf(histName,"histNuEBarCCNDMinus1_%d",systNum);
277  histNuEBarCCNDMinus1[systNum]->SetName(histName);
278  histNuEBarCCFDMinus1[systNum] = (TH1F*) fpMinus1Shift->Get(fdNuEBarCCHistName);
279  sprintf(histName,"histNuEBarCCFDMinus1_%d",systNum);
280  histNuEBarCCFDMinus1[systNum]->SetName(histName);
281 
282  histNuECCNDMinus1[systNum]->Add(histNuEBarCCNDMinus1[systNum]);
283  histNuECCFDMinus1[systNum]->Add(histNuEBarCCFDMinus1[systNum]);
284 
285  }
286  // return;
287  char legText[180];
288  char canName[180];
289  char ratioName[180];
290  for(int systNum=0;systNum<numSysts;systNum++) {
291  sprintf(canName,"canND%d",systNum);
292  TCanvas *canND = new TCanvas(canName,canName,800,800);
293  canND->Divide(1,2);
294  canND->cd(1);
295  histFullNDMinus1[systNum]->SetLineColor(getFnexColour(1));
296  histFullNDMinus1[systNum]->SetLineStyle(2);
297  histFullNDMinus1[systNum]->Draw("hist");
298  histFullNDMinus1[systNum]->GetXaxis()->SetRangeUser(0,5);
299  histFullNDNoShift->Draw("hist same");
300  histFullNDPlus1[systNum]->SetLineColor(getFnexColour(1));
301  histFullNDPlus1[systNum]->Draw("hist same");
302  // histNCNDPlus1[systNum]->SetLineColor(getFnexColour(2));
303  // histNCNDPlus1[systNum]->Draw("hist same");
304  // histNCNDMinus1[systNum]->SetLineColor(getFnexColour(2));
305  // histNCNDMinus1[systNum]->SetLineStyle(2);
306  // histNCNDMinus1[systNum]->Draw("hist same");
307  // histNuECCNDPlus1[systNum]->SetLineColor(getFnexColour(4));
308  // histNuECCNDPlus1[systNum]->Draw("hist same");
309  // histNuECCNDMinus1[systNum]->SetLineColor(getFnexColour(4));
310  // histNuECCNDMinus1[systNum]->SetLineStyle(2);
311  // histNuECCNDMinus1[systNum]->Draw("hist same");
312  Double_t ndPlusMax=histFullNDMinus1[systNum]->GetMaximum();
313  if(ndPlusMax<histFullNDPlus1[systNum]->GetMaximum()) {
314  ndPlusMax=histFullNDPlus1[systNum]->GetMaximum();
315  }
316  if(ndPlusMax<histFullNDNoShift->GetMaximum()) {
317  ndPlusMax=histFullNDNoShift->GetMaximum();
318  }
319  histFullNDMinus1[systNum]->GetYaxis()->SetRangeUser(0,1.1*ndPlusMax);
320  canND->cd(2);
321  sprintf(ratioName,"ratioFullNDPlus1_%d",systNum);
322  TH1F *ratioFullNDPlus1 = (TH1F*)histFullNDPlus1[systNum]->Clone(ratioName);
323  ratioFullNDPlus1->Divide(histFullNDNoShift);
324  ratioFullNDPlus1->SetLineWidth(3);
325  ratioFullNDPlus1->Draw("hist");
326  ratioFullNDPlus1->GetXaxis()->SetRangeUser(0,5);
327  ratioFullNDPlus1->GetYaxis()->SetRangeUser(0.5,1.5);
328  ratioFullNDPlus1->SetYTitle("Ratio to nominal");
329  sprintf(ratioName,"ratioFullNDMinus1_%d",systNum);
330  TH1F *ratioFullNDMinus1 = (TH1F*)histFullNDMinus1[systNum]->Clone(ratioName);
331  ratioFullNDMinus1->Divide(histFullNDNoShift);
332  ratioFullNDMinus1->SetLineWidth(3);
333  ratioFullNDMinus1->Draw("same hist");
334  // TH1F *ratioNCNDPlus1 = (TH1F*)histNCNDPlus1[systNum]->Clone(ratioName);
335  // ratioNCNDPlus1->Divide(histNCNDNoShift);
336  // ratioNCNDPlus1->SetLineWidth(1);
337  // ratioNCNDPlus1->Draw("same hist");
338  // TH1F *ratioNCNDMinus1 = (TH1F*)histNCNDMinus1[systNum]->Clone(ratioName);
339  // ratioNCNDMinus1->Divide(histNCNDNoShift);
340  // ratioNCNDMinus1->SetLineWidth(1);
341  // ratioNCNDMinus1->Draw("same hist");
342  // TH1F *ratioNuECCNDPlus1 = (TH1F*)histNuECCNDPlus1[systNum]->Clone(ratioName);
343  // ratioNuECCNDPlus1->Divide(histNuECCNDNoShift);
344  // ratioNuECCNDPlus1->SetLineWidth(1);
345  // ratioNuECCNDPlus1->Draw("same hist");
346  // TH1F *ratioNuECCNDMinus1 = (TH1F*)histNuECCNDMinus1[systNum]->Clone(ratioName);
347  // ratioNuECCNDMinus1->Divide(histNuECCNDNoShift);
348  // ratioNuECCNDMinus1->SetLineWidth(1);
349  // ratioNuECCNDMinus1->Draw("same hist");
350 
351  canND->cd(1);
352  TLegend *leggyND = makeLegendBoxInTopRightCorner(0.3,0.3);
353  leggyND->AddEntry(histFullNDNoShift,"ND Nominal MC","l");
354  sprintf(legText,"ND -- %s +1#sigma",legTags[systNum+1]);
355  leggyND->AddEntry(histFullNDPlus1[systNum],legText,"l");
356  sprintf(legText,"ND -- %s -1#sigma",legTags[systNum+1]);
357  leggyND->AddEntry(histFullNDMinus1[systNum],legText,"l");
358  // sprintf(legText,"ND NC -- %s +1#sigma",legTags[systNum+1]);
359  // leggyND->AddEntry(histNCNDPlus1[systNum],legText,"l");
360  // // sprintf(legText,"ND NC -- %s -1#sigma",legTags[systNum+1]);
361  // // leggyND->AddEntry(histNCNDMinus1[systNum],legText,"l");
362  // sprintf(legText,"ND #nu_{e} CC -- %s +1#sigma",legTags[systNum+1]);
363  // leggyND->AddEntry(histNuECCNDPlus1[systNum],legText,"l");
364  leggyND->Draw();
365  PrintPDFandPNG(canND,canName);
366 
367 
368  sprintf(canName,"canFD%d",systNum);
369  TCanvas *canFD = new TCanvas(canName,canName,800,800);
370  canFD->Divide(1,2);
371  canFD->cd(1);
372  histFullFDMinus1[systNum]->SetLineColor(getFnexColour(1));
373  histFullFDMinus1[systNum]->SetLineStyle(2);
374  histFullFDMinus1[systNum]->Draw("hist");
375  histFullFDMinus1[systNum]->GetXaxis()->SetRangeUser(0,5);
376  histFullFDNoShift->Draw("hist same");
377  histFullFDPlus1[systNum]->SetLineColor(getFnexColour(1));
378  histFullFDPlus1[systNum]->Draw("hist same");
379  // histNCFDPlus1[systNum]->SetLineColor(getFnexColour(2));
380  // histNCFDPlus1[systNum]->Draw("hist same");
381  // histNCFDMinus1[systNum]->SetLineColor(getFnexColour(2));
382  // histNCFDMinus1[systNum]->SetLineStyle(2);
383  // histNCFDMinus1[systNum]->Draw("hist same");
384  // histNuTauCCFDPlus1[systNum]->SetLineColor(getFnexColour(3));
385  // histNuTauCCFDPlus1[systNum]->Draw("hist same");
386  // histNuTauCCFDMinus1[systNum]->SetLineColor(getFnexColour(3));
387  // histNuTauCCFDMinus1[systNum]->SetLineStyle(2);
388  // histNuTauCCFDMinus1[systNum]->Draw("hist same");
389  // histNuECCFDPlus1[systNum]->SetLineColor(getFnexColour(4));
390  // histNuECCFDPlus1[systNum]->Draw("hist same");
391  // histNuECCFDMinus1[systNum]->SetLineColor(getFnexColour(4));
392  // histNuECCFDMinus1[systNum]->SetLineStyle(2);
393  // histNuECCFDMinus1[systNum]->Draw("hist same");
394  // histCosmicFDPlus1[systNum]->SetLineColor(getFnexColour(5));
395  // histCosmicFDPlus1[systNum]->Draw("hist same");
396  // histCosmicFDMinus1[systNum]->SetLineColor(getFnexColour(5));
397  // histCosmicFDMinus1[systNum]->SetLineStyle(2);
398  // histCosmicFDMinus1[systNum]->Draw("hist same");
399  Double_t fdPlusMax=histFullFDMinus1[systNum]->GetMaximum();
400  if(fdPlusMax<histFullFDPlus1[systNum]->GetMaximum()) {
401  fdPlusMax=histFullFDPlus1[systNum]->GetMaximum();
402  }
403  if(fdPlusMax<histFullFDNoShift->GetMaximum()) {
404  fdPlusMax=histFullFDNoShift->GetMaximum();
405  }
406  histFullFDMinus1[systNum]->GetYaxis()->SetRangeUser(0,1.1*fdPlusMax);
407  canFD->cd(2);
408  sprintf(ratioName,"ratioFullFDPlus1_%d",systNum);
409  TH1F *ratioFullFDPlus1 = (TH1F*)histFullFDPlus1[systNum]->Clone(ratioName);
410  ratioFullFDPlus1->Divide(histFullFDNoShift);
411  ratioFullFDPlus1->SetLineWidth(3);
412  ratioFullFDPlus1->Draw("hist");
413  ratioFullFDPlus1->GetXaxis()->SetRangeUser(0,5);
414  ratioFullFDPlus1->GetYaxis()->SetRangeUser(0.5,1.5);
415  ratioFullFDPlus1->SetYTitle("Ratio to nominal");
416  sprintf(ratioName,"ratioFullFDMinus1_%d",systNum);
417  TH1F *ratioFullFDMinus1 = (TH1F*)histFullFDMinus1[systNum]->Clone(ratioName);
418  ratioFullFDMinus1->Divide(histFullFDNoShift);
419  ratioFullFDMinus1->SetLineWidth(3);
420  ratioFullFDMinus1->Draw("same hist");
421  // TH1F *ratioNCFDPlus1 = (TH1F*)histNCFDPlus1[systNum]->Clone(ratioName);
422  // ratioNCFDPlus1->Divide(histNCFDNoShift);
423  // ratioNCFDPlus1->SetLineWidth(1);
424  // ratioNCFDPlus1->Draw("same hist");
425  // TH1F *ratioNCFDMinus1 = (TH1F*)histNCFDMinus1[systNum]->Clone(ratioName);
426  // ratioNCFDMinus1->Divide(histNCFDNoShift);
427  // ratioNCFDMinus1->SetLineWidth(1);
428  // ratioNCFDMinus1->Draw("same hist");
429  // TH1F *ratioNuECCFDPlus1 = (TH1F*)histNuECCFDPlus1[systNum]->Clone(ratioName);
430  // ratioNuECCFDPlus1->Divide(histNuECCFDNoShift);
431  // ratioNuECCFDPlus1->SetLineWidth(1);
432  // ratioNuECCFDPlus1->Draw("same hist");
433  // TH1F *ratioNuECCFDMinus1 = (TH1F*)histNuECCFDMinus1[systNum]->Clone(ratioName);
434  // ratioNuECCFDMinus1->Divide(histNuECCFDNoShift);
435  // ratioNuECCFDMinus1->SetLineWidth(1);
436  // ratioNuECCFDMinus1->Draw("same hist");
437  // TH1F *ratioNuTauCCFDPlus1 = (TH1F*)histNuTauCCFDPlus1[systNum]->Clone(ratioName);
438  // ratioNuTauCCFDPlus1->Divide(histNuTauCCFDNoShift);
439  // ratioNuTauCCFDPlus1->SetLineWidth(1);
440  // ratioNuTauCCFDPlus1->Draw("same hist");
441  // TH1F *ratioNuTauCCFDMinus1 = (TH1F*)histNuTauCCFDMinus1[systNum]->Clone(ratioName);
442  // ratioNuTauCCFDMinus1->Divide(histNuTauCCFDNoShift);
443  // ratioNuTauCCFDMinus1->SetLineWidth(1);
444  // ratioNuTauCCFDMinus1->Draw("same hist");
445  // TH1F *ratioCosmicFDPlus1 = (TH1F*)histCosmicFDPlus1[systNum]->Clone(ratioName);
446  // ratioCosmicFDPlus1->Divide(histCosmicFDNoShift);
447  // ratioCosmicFDPlus1->Draw("same hist");
448  // TH1F *ratioCosmicFDMinus1 = (TH1F*)histCosmicFDMinus1[systNum]->Clone(ratioName);
449  // ratioCosmicFDMinus1->Divide(histCosmicFDNoShift);
450  // ratioCosmicFDMinus1->Draw("same hist");
451 
452  canFD->cd(1);
453  TLegend *leggyFD = makeLegendBoxInTopRightCorner(0.3,0.3);
454  leggyFD->AddEntry(histFullFDNoShift,"FD Nominal MC","l");
455  sprintf(legText,"FD -- %s +1#sigma",legTags[systNum+1]);
456  leggyFD->AddEntry(histFullFDPlus1[systNum],legText,"l");
457  sprintf(legText,"FD -- %s -1#sigma",legTags[systNum+1]);
458  leggyFD->AddEntry(histFullFDMinus1[systNum],legText,"l");
459  // sprintf(legText,"FD NC -- %s +1#sigma",legTags[systNum+1]);
460  // leggyFD->AddEntry(histNCFDPlus1[systNum],legText,"l");
461  // sprintf(legText,"FD #nu_{#tau} CC -- %s +1#sigma",legTags[systNum+1]);
462  // leggyFD->AddEntry(histNuTauCCFDPlus1[systNum],legText,"l");
463  // sprintf(legText,"FD #nu_{e} CC -- %s +1#sigma",legTags[systNum+1]);
464  // leggyFD->AddEntry(histNuECCFDPlus1[systNum],legText,"l");
465  // sprintf(legText,"FD Cosmic -- %s +1#sigma",legTags[systNum+1]);
466  // leggyFD->AddEntry(histCosmicFDPlus1[systNum],legText,"l");
467  // sprintf(legText,"FD NC -- %s -1#sigma",legTags[systNum+1]);
468  // leggyFD->AddEntry(histNCFDMinus1[systNum],legText,"l");
469  leggyFD->Draw();
470  PrintPDFandPNG(canFD,canName);
471 
472  }
473 
474  TCanvas *canErrorND = new TCanvas("canErrorND","canErrorND",800,800);
475  TGraphAsymmErrors *grND=getErrorBand(histFullNDNoShift,histFullNDPlus1,histFullNDMinus1,numSysts);
476  histFullNDNoShift->Draw();
477  histFullNDNoShift->SetMaximum(histFullNDNoShift->GetMaximum()*1.3);
478  histFullNDNoShift->GetXaxis()->SetTitle("Reco Neutrino Energy (Gev)");
479  grND->SetFillColor(getFnexColour(1));
480  grND->SetFillStyle(3001);
481  grND->Draw("e3");
482  histFullNDNoShift->DrawCopy("same");
483  histFullNDNoShift->GetXaxis()->SetRangeUser(0,5);
484  {
485  TLegend *leggyND = makeLegendBoxInTopRightCorner(0.3,0.3);
486  leggyND->AddEntry(histFullNDNoShift,"ND Nominal MC","l");
487  leggyND->AddEntry(grND,"ND Error Band","f");
488  leggyND->Draw();
489  }
490  TCanvas *canErrorFD = new TCanvas("canErrorFD","canErrorFD",800,800);
491  TGraphAsymmErrors *grFD=getErrorBand(histFullFDNoShift,histFullFDPlus1,histFullFDMinus1,numSysts);
492  histFullFDNoShift->Draw();
493  histFullFDNoShift->SetMaximum(histFullFDNoShift->GetMaximum()*1.3);
494  histFullFDNoShift->GetXaxis()->SetTitle("Reco Neutrino Energy (Gev)");
495  grFD->SetFillColor(getFnexColour(1));
496  grFD->SetFillStyle(3001);
497  grFD->Draw("e3");
498  histFullFDNoShift->DrawCopy("same");
499  histFullFDNoShift->GetXaxis()->SetRangeUser(0,5);
500  {
501  TLegend *leggyFD = makeLegendBoxInTopRightCorner(0.3,0.3);
502  leggyFD->AddEntry(histFullFDNoShift,"FD Nominal Prediction","l");
503  leggyFD->AddEntry(grFD,"FD Error Band","f");
504  leggyFD->Draw();
505  }
506 
507 
508 }
509 
510 TGraphAsymmErrors *getErrorBand(TH1F *histNoShift, TH1F **histPlus, TH1F **histMinus, Int_t numShifts)
511 {
512  char histName[180];
513  sprintf(histName,"%s_systband",histNoShift->GetName());
514  TGraphAsymmErrors *grErrorBand = new TGraphAsymmErrors();
515  grErrorBand->SetName(histName);
516  for(int bin=1;bin<=histNoShift->GetNbinsX();bin++) {
517  double plusError=0;
518  double minusError=0;
519  double nominal=histNoShift->GetBinContent(bin);
520  double binCent=histNoShift->GetBinCenter(bin);
521  double binLow=histNoShift->GetBinLowEdge(bin);
522  double binHigh=binLow+histNoShift->GetBinWidth(bin);
523  for(int shift=0;shift<numShifts;shift++) {
524  double plus=histPlus[shift]->GetBinContent(bin);
525  double minus=histMinus[shift]->GetBinContent(bin);
526  if(plus<minus) {
527  double temp=plus;
528  plus=minus;
529  minus=temp;
530  }
531  plusError+=(plus-nominal)*(plus-nominal);
532  minusError+=(minus-nominal)*(minus-nominal);
533  }
534  plusError=sqrt(plusError);
535  minusError=sqrt(minusError);
536  grErrorBand->SetPoint(bin-1,binCent,nominal);
537  grErrorBand->SetPointEXlow(bin-1,binCent-binLow);
538  grErrorBand->SetPointEXhigh(bin-1,binHigh-binCent);
539  grErrorBand->SetPointEYlow(bin-1,plusError);
540  grErrorBand->SetPointEYhigh(bin-1,minusError);
541  }
542 
543  return grErrorBand;
544 }
545 
546 
547 
548 
549 
550 
551 
552 
enum BeamMode kOrange
enum BeamMode kRed
void plotSystBandForTheseSystematics(char *rootFileList, char *legTagList=0)
T sqrt(T number)
Definition: d0nt_math.hpp:156
TGraphAsymmErrors * getErrorBand(TH1F *histNoShift, TH1F **histPlus, TH1F **histMinus, Int_t numShifts)
int getFnexColour(int index)
float bin[41]
Definition: plottest35.C:14
OStream cout
Definition: OStream.cxx:6
def countFiles(dir, testStrs=[])
enum BeamMode kViolet
enum BeamMode kGreen
enum BeamMode kBlue
T minus(const T &x)
Definition: minus.hpp:15