plotNueSystBandForTheseSystematics.C
Go to the documentation of this file.
1 
2 #include <iostream>
3 #include <cstring>
4 
5 #include "TStyle.h"
6 
7  const int NUM_PID_BINS=3;
8 
9 
11 {
12  Int_t niceColours[16]={kBlue-3,kRed-3,kGreen-3,
13  kMagenta-3,kCyan-3,kOrange-3,
14  kViolet-3,kTeal-3,
15  kBlue+1,kRed-3,kGreen+1,
16  kMagenta+1,kCyan+1,kOrange+1,
17  kViolet+1,kTeal+1};
18  if(index<16) return niceColours[index];
19 
20  Int_t offset=index/16;
21  Int_t colourInd=index%16;
22 
23  return (niceColours[colourInd]+offset);
24 
25 }
26 
27 TGraphAsymmErrors *getFractionalErrorBand(TH1F *histNoShift[3], TH1F *histPlus[1024][3], TH1F *histMinus[1024][3], Int_t numShifts, Int_t pidInd) ;
28 TGraphAsymmErrors *getErrorBand(TH1F *histNoShift[3], TH1F *histPlus[1024][3], TH1F *histMinus[1024][3], Int_t numShifts, Int_t pidInd) ;
29 TCanvas *getCanvasWithHists(int systNum,const char *ndOrFD,const char *legTag,TH1F *histFullNoShift[3],TH1F *histFullPlus1[3],TH1F *histFullMinus1[3],TH1F *histSignalNoShift[3],TH1F *histSignalPlus1[3],TH1F *histSignalMinus1[3],TH1F *histBackgroundNoShift[3],TH1F *histBackgroundPlus1[3],TH1F *histBackgroundMinus1[3] );
30 
31 void plotNueSystBandForTheseSystematics(char *rootFileList,char *legTagList=0) {
32  std::cout << rootFileList << "\n";
33  gStyle->SetOptStat(00000000);
34  // return;
35  std::ifstream InFileList(rootFileList);
36  char fileNames[1024][1024];
37  int countFiles=0;
38  while(InFileList >> fileNames[countFiles]) {
39  std::cout << fileNames[countFiles] << "\n";
40  countFiles++;
41  }
42  std::cout << "Got " << countFiles << " files\n";
43 
44  char legTags[1024][100];
45  for(int i=0;i<1024;i++) {
46  sprintf(legTags[i],"Syst.");
47  }
48  if(legTagList) {
49  int countTags=0;
50  std::ifstream LegTagList(legTagList);
51  while(LegTagList >> legTags[countTags]) {
52  for(int si=0;si<strlen(legTags[countTags]);si++) {
53  // std::cout << si << "\t" << legTags[countTags][si] << endl;
54  if(legTags[countTags][si]=='~') {
55  legTags[countTags][si]=' ';
56  }
57  }
58  std::cout << legTags[countTags] << "\n";//t" << index(legTags,'~') << "\n";
59  countTags++;
60  }
61  std::cout << "Got " << countTags << " files\n";
62 
63 
64  }
65 
66 
67  // KEY: TH1F FDCosmicMuonCorrected;1
68 
69 
70 
71  const char *histNamePreFactor="plot/MultiExperiment";
72  const char *nueSelectionLabels[NUM_PID_BINS]={"NOvA_NuESel_LowPID_NuMuSel","NOvA_NuESel_MidPID_NuMuSel","NOvA_NuESel_HighPID_NuMuSel"};
73 
74 
75 
76 
77  const char *ndFullHistName="CorrectedStacks/NDMCCorrected";
78  const char *fdFullHistName="CorrectedStacks/FDMCCorrected";
79  const char *fdSignalHistName="CorrectedStacks/FDmc_meCorrected";
80  const char *fdSignalHistName2="CorrectedStacks/FDmc_mebarCorrected";
81  // const char*ndFullHistName="plot/MultiExperiment/NOvA_NuESel_HighPID_NuMuSel/UncorrectedStacks/NDMCUncorrected";
82  // const char*fdFullHistName="plot/MultiExperiment/NOvA_NuESel_HighPID_NuMuSel/UncorrectedStacks/FDMCUncorrected";
83  // const char*ndNCHistName="plot/MultiExperiment/NOvA_NuESel_HighPID_NuMuSel/CorrectedStacks/NDmc_ncCorrected";
84  // const char*fdNCHistName="plot/MultiExperiment/NOvA_NuESel_HighPID_NuMuSel/CorrectedStacks/FDmc_ncCorrected";
85 
86  // const char*ndNuTauCCHistName="plot/MultiExperiment/NOvA_NuESel_HighPID_NuMuSel/CorrectedStacks/NDNuTauCCCorrected";
87  // const char*fdNuTauCCHistName="plot/MultiExperiment/NOvA_NuESel_HighPID_NuMuSel/CorrectedStacks/FDmc_mtCorrected";
88  // const char*ndNuTauBarCCHistName="plot/MultiExperiment/NOvA_NuESel_HighPID_NuMuSel/CorrectedStacks/NDNuTauBarCCCorrected";
89  // const char*fdNuTauBarCCHistName="plot/MultiExperiment/NOvA_NuESel_HighPID_NuMuSel/CorrectedStacks/FDmc_mtbarCorrected";
90 
91  // const char*ndNuECCHistName="plot/MultiExperiment/NOvA_NuESel_HighPID_NuMuSel/CorrectedStacks/NDmc_eeCorrected";
92  // const char*fdNuECCHistName="plot/MultiExperiment/NOvA_NuESel_HighPID_NuMuSel/CorrectedStacks/FDmc_eeCorrected";
93  // const char*ndNuEBarCCHistName="plot/MultiExperiment/NOvA_NuESel_HighPID_NuMuSel/CorrectedStacks/NDmc_eebarCorrected";
94  // const char*fdNuEBarCCHistName="plot/MultiExperiment/NOvA_NuESel_HighPID_NuMuSel/CorrectedStacks/FDmc_eebarCorrected";
95  // const char*fdCosmicHistName="plot/MultiExperiment/NOvA_NuESel_HighPID_NuMuSel/CorrectedStacks/FDcosmicCorrected";
96 
97 
98  char histName[FILENAME_MAX];
99  TFile *fpNoShift = TFile::Open(fileNames[0]); //0th file is the nominal
100  TH1F *histFullNDNoShift[NUM_PID_BINS];
101  TH1F *histFullFDNoShift[NUM_PID_BINS];
102  TH1F *histFDSignalNoShift[NUM_PID_BINS];
103  TH1F *histFDBackgroundNoShift[NUM_PID_BINS];
104 
105  Double_t fdSignalTotal=0;
106  Double_t fdBackgroundTotal=0;
107 
108 
109  for(int pidInd=0;pidInd<NUM_PID_BINS;pidInd++) {
110  sprintf(histName,"%s/%s/%s",histNamePreFactor,nueSelectionLabels[pidInd],ndFullHistName);
111  histFullNDNoShift[pidInd]= (TH1F*)fpNoShift->Get(histName);
112  sprintf(histName,"%s/%s/%s",histNamePreFactor,nueSelectionLabels[pidInd],fdFullHistName);
113  histFullFDNoShift[pidInd]= (TH1F*)fpNoShift->Get(histName);
114  sprintf(histName,"%s/%s/%s",histNamePreFactor,nueSelectionLabels[pidInd],fdSignalHistName);
115  histFDSignalNoShift[pidInd]= (TH1F*)fpNoShift->Get(histName);
116  sprintf(histName,"%s/%s/%s",histNamePreFactor,nueSelectionLabels[pidInd],fdSignalHistName2);
117  TH1F *htemp= (TH1F*)fpNoShift->Get(histName);
118  histFDSignalNoShift[pidInd]->Add(htemp);
119  sprintf(histName,"histBackground%d",pidInd);
120  histFDBackgroundNoShift[pidInd]=(TH1F*)histFullFDNoShift[pidInd]->Clone(histName);
121  histFDBackgroundNoShift[pidInd]->Add(histFDSignalNoShift[pidInd],-1);
122 
123  fdSignalTotal+=histFDSignalNoShift[pidInd]->Integral();
124  fdBackgroundTotal+=histFDBackgroundNoShift[pidInd]->Integral();
125 
126  }
127 
128 
129  int numSysts=(countFiles-1)/2;
130  std::cout << "Got " << numSysts << " systematics\n";
131  std::cout << "Signal: " << fdSignalTotal << "\tBackground: " << fdBackgroundTotal << "\n";
132 
133  const char *legTag="Syst";
134 
135  TH1F *histFullNDPlus1[1024][NUM_PID_BINS];
136  TH1F *histFullNDMinus1[1024][NUM_PID_BINS];
137  TH1F *histFullFDPlus1[1024][NUM_PID_BINS];
138  TH1F *histFullFDMinus1[1024][NUM_PID_BINS];
139  TH1F *histFDSignalPlus1[1024][NUM_PID_BINS];
140  TH1F *histFDSignalMinus1[1024][NUM_PID_BINS];
141  TH1F *histFDBackgroundPlus1[1024][NUM_PID_BINS];
142  TH1F *histFDBackgroundMinus1[1024][NUM_PID_BINS];
143 
144  Double_t fdSignalTotalSyst[1024][2]={{0}};
145  Double_t fdBackgroundTotalSyst[1024][2]={{0}};
146 
147 
148  // numSysts=1;
149  for(int systNum=0;systNum<numSysts;systNum++) {
150  std::cout << "Here " << systNum << "\n";
151  TFile *fpPlus1Shift = TFile::Open(fileNames[(systNum*2)+1]);
152  for(int pidInd=0;pidInd<NUM_PID_BINS;pidInd++) {
153  sprintf(histName,"%s/%s/%s",histNamePreFactor,nueSelectionLabels[pidInd],ndFullHistName);
154  histFullNDPlus1[systNum][pidInd]= (TH1F*)fpPlus1Shift->Get(histName);
155  sprintf(histName,"%s/%s/%s",histNamePreFactor,nueSelectionLabels[pidInd],fdFullHistName);
156  histFullFDPlus1[systNum][pidInd]= (TH1F*)fpPlus1Shift->Get(histName);
157  sprintf(histName,"%s/%s/%s",histNamePreFactor,nueSelectionLabels[pidInd],fdSignalHistName);
158  histFDSignalPlus1[systNum][pidInd]= (TH1F*)fpPlus1Shift->Get(histName);
159  sprintf(histName,"%s/%s/%s",histNamePreFactor,nueSelectionLabels[pidInd],fdSignalHistName2);
160  TH1F *htemp= (TH1F*)fpPlus1Shift->Get(histName);
161  histFDSignalPlus1[systNum][pidInd]->Add(htemp);
162  sprintf(histName,"histBackground%d",pidInd);
163  histFDBackgroundPlus1[systNum][pidInd]=(TH1F*)histFullFDPlus1[systNum][pidInd]->Clone(histName);
164  histFDBackgroundPlus1[systNum][pidInd]->Add(histFDSignalPlus1[systNum][pidInd],-1);
165 
166  fdSignalTotalSyst[systNum][0]+=histFDSignalPlus1[systNum][pidInd]->Integral();
167  fdBackgroundTotalSyst[systNum][0]+=histFDBackgroundPlus1[systNum][pidInd]->Integral();
168 
169  }
170 
171 
172  TFile *fpMinus1Shift = TFile::Open(fileNames[(systNum*2)+2]);
173  for(int pidInd=0;pidInd<NUM_PID_BINS;pidInd++) {
174  sprintf(histName,"%s/%s/%s",histNamePreFactor,nueSelectionLabels[pidInd],ndFullHistName);
175  histFullNDMinus1[systNum][pidInd]= (TH1F*)fpMinus1Shift->Get(histName);
176  sprintf(histName,"%s/%s/%s",histNamePreFactor,nueSelectionLabels[pidInd],fdFullHistName);
177  histFullFDMinus1[systNum][pidInd]= (TH1F*)fpMinus1Shift->Get(histName);
178  sprintf(histName,"%s/%s/%s",histNamePreFactor,nueSelectionLabels[pidInd],fdSignalHistName);
179  histFDSignalMinus1[systNum][pidInd]= (TH1F*)fpMinus1Shift->Get(histName);
180  sprintf(histName,"%s/%s/%s",histNamePreFactor,nueSelectionLabels[pidInd],fdSignalHistName2);
181  TH1F *htemp= (TH1F*)fpMinus1Shift->Get(histName);
182  histFDSignalMinus1[systNum][pidInd]->Add(htemp);
183  sprintf(histName,"histBackground%d",pidInd);
184  histFDBackgroundMinus1[systNum][pidInd]=(TH1F*)histFullFDMinus1[systNum][pidInd]->Clone(histName);
185  histFDBackgroundMinus1[systNum][pidInd]->Add(histFDSignalMinus1[systNum][pidInd],-1);
186  fdSignalTotalSyst[systNum][1]+=histFDSignalMinus1[systNum][pidInd]->Integral();
187  fdBackgroundTotalSyst[systNum][1]+=histFDBackgroundMinus1[systNum][pidInd]->Integral();
188  }
189 
190  }
191  // return;
192 
193 
194  for(int systNum=0;systNum<numSysts;systNum++) {
195  TCanvas *canND = getCanvasWithHists(systNum,"ND",legTags[systNum+1],histFullNDNoShift,histFullNDPlus1[systNum],histFullNDMinus1[systNum],NULL,NULL,NULL,NULL,NULL,NULL);
196  TCanvas *canFDAll = getCanvasWithHists(systNum,"FD",legTags[systNum+1],histFullFDNoShift,histFullFDPlus1[systNum],histFullFDMinus1[systNum],histFDSignalNoShift,histFDSignalPlus1[systNum],histFDSignalMinus1[systNum],histFDBackgroundNoShift,histFDBackgroundPlus1[systNum],histFDBackgroundMinus1[systNum]);
197  // TCanvas *canFDSignal = getCanvasWithHists(systNum,"FD_Signal",legTags[systNum+1],histFDSignalNoShift,histFDSignalPlus1[systNum],histFDSignalMinus1[systNum]);
198 
199  }
200 
201  for(int systNum=0;systNum<numSysts;systNum++) {
202  std::cout << legTags[systNum+1] << "\t" << fdSignalTotalSyst[systNum][0] << "\t" << fdSignalTotalSyst[systNum][1] << "\t"
203  << fdBackgroundTotalSyst[systNum][0] << "\t" << fdBackgroundTotalSyst[systNum][1] << "\n";
204 
205  }
206 
207  Int_t pidInd=0;
208 
209  TCanvas *canErrorND = new TCanvas("canErrorND","canErrorND",800,800);
210  TGraphAsymmErrors *grND=getErrorBand(histFullNDNoShift,histFullNDPlus1,histFullNDMinus1,numSysts,pidInd);
211  histFullNDNoShift[pidInd]->Draw();
212  histFullNDNoShift[pidInd]->SetMaximum(histFullNDNoShift[pidInd]->GetMaximum()*1.3);
213  histFullNDNoShift[pidInd]->GetXaxis()->SetTitle("Reco Neutrino Energy (Gev)");
214  grND->SetFillColor(getFnexColour(1));
215  grND->SetFillStyle(3001);
216  grND->Draw("e3");
217  histFullNDNoShift[pidInd]->DrawCopy("same");
218  histFullNDNoShift[pidInd]->GetXaxis()->SetRangeUser(0,5);
219  {
220  TLegend *leggyND = makeLegendBoxInTopRightCorner(0.3,0.3);
221  leggyND->AddEntry(histFullNDNoShift[pidInd],"ND Nominal MC","l");
222  leggyND->AddEntry(grND,"ND Error Band","f");
223  leggyND->Draw();
224  }
225  TCanvas *canErrorFD = new TCanvas("canErrorFD","canErrorFD",800,800);
226  TGraphAsymmErrors *grFD=getErrorBand(histFullFDNoShift,histFullFDPlus1,histFullFDMinus1,numSysts,pidInd);
227  histFullFDNoShift[pidInd]->Draw();
228  histFullFDNoShift[pidInd]->SetMaximum(histFullFDNoShift[pidInd]->GetMaximum()*1.3);
229  histFullFDNoShift[pidInd]->GetXaxis()->SetTitle("Reco Neutrino Energy (Gev)");
230  grFD->SetFillColor(getFnexColour(1));
231  grFD->SetFillStyle(3001);
232  grFD->Draw("e4");
233  histFullFDNoShift[pidInd]->DrawCopy("same");
234  histFullFDNoShift[pidInd]->GetXaxis()->SetRangeUser(0,5);
235  {
236  TLegend *leggyFD = makeLegendBoxInTopRightCorner(0.3,0.3);
237  leggyFD->AddEntry(histFullFDNoShift[pidInd],"FD Nominal Prediction","l");
238  leggyFD->AddEntry(grFD,"FD Error Band","f");
239  leggyFD->Draw();
240  }
241 
242 
243  TCanvas *canFrac = new TCanvas("canFrac","canFrac");
244  canFrac->Divide(1,2);
245  canFrac->cd(1);
246  TGraphAsymmErrors *grFracND=getFractionalErrorBand(histFullNDNoShift,histFullNDPlus1,histFullNDMinus1,numSysts,pidInd);
247  grFracND->SetFillColor(getFnexColour(1));
248  grFracND->SetFillStyle(3001);
249  grFracND->Draw("ap");
250  grFracND->GetXaxis()->SetRangeUser(0,5);
251  canFrac->cd(2);
252  TGraphAsymmErrors *grFracFD=getFractionalErrorBand(histFullFDNoShift,histFullFDPlus1,histFullFDMinus1,numSysts,pidInd);
253  grFracFD->SetFillColor(getFnexColour(1));
254  grFracFD->SetFillStyle(3001);
255  grFracFD->Draw("ap");
256  grFracFD->GetXaxis()->SetRangeUser(0,5);
257 
258 
259 
260 
261 }
262 
263 TGraphAsymmErrors *getErrorBand(TH1F *histNoShift[3], TH1F *histPlus[1024][3], TH1F *histMinus[1024][3], Int_t numShifts, Int_t pidInd)
264 
265 {
266  char histName[180];
267  sprintf(histName,"%s_systband",histNoShift[pidInd]->GetName());
268  TGraphAsymmErrors *grErrorBand = new TGraphAsymmErrors();
269  grErrorBand->SetName(histName);
270  for(int bin=1;bin<=histNoShift[pidInd]->GetNbinsX();bin++) {
271  double plusError=0;
272  double minusError=0;
273  double nominal=histNoShift[pidInd]->GetBinContent(bin);
274  double binCent=histNoShift[pidInd]->GetBinCenter(bin);
275  double binLow=histNoShift[pidInd]->GetBinLowEdge(bin);
276  double binHigh=binLow+histNoShift[pidInd]->GetBinWidth(bin);
277  for(int shift=0;shift<numShifts;shift++) {
278  double plus=histPlus[shift][pidInd]->GetBinContent(bin);
279  double minus=histMinus[shift][pidInd]->GetBinContent(bin);
280  if(plus<minus) {
281  double temp=plus;
282  plus=minus;
283  minus=temp;
284  }
285  plusError+=(plus-nominal)*(plus-nominal);
286  minusError+=(minus-nominal)*(minus-nominal);
287  }
288  plusError=sqrt(plusError);
289  minusError=sqrt(minusError);
290  grErrorBand->SetPoint(bin-1,binCent,nominal);
291  grErrorBand->SetPointEXlow(bin-1,binCent-binLow);
292  grErrorBand->SetPointEXhigh(bin-1,binHigh-binCent);
293  grErrorBand->SetPointEYlow(bin-1,plusError);
294  grErrorBand->SetPointEYhigh(bin-1,minusError);
295  }
296 
297 
298  return grErrorBand;
299 }
300 
301 
302 TGraphAsymmErrors *getFractionalErrorBand(TH1F *histNoShift[3], TH1F *histPlus[1024][3], TH1F *histMinus[1024][3], Int_t numShifts, Int_t pidInd)
303 
304 {
305  char histName[180];
306  sprintf(histName,"%s_systband",histNoShift[pidInd]->GetName());
307  TGraphAsymmErrors *grErrorBand = new TGraphAsymmErrors();
308  grErrorBand->SetName(histName);
309  for(int bin=1;bin<=histNoShift[pidInd]->GetNbinsX();bin++) {
310  double plusError=0;
311  double minusError=0;
312  double nominal=histNoShift[pidInd]->GetBinContent(bin);
313  double binCent=histNoShift[pidInd]->GetBinCenter(bin);
314  double binLow=histNoShift[pidInd]->GetBinLowEdge(bin);
315  double binHigh=binLow+histNoShift[pidInd]->GetBinWidth(bin);
316  for(int shift=0;shift<numShifts;shift++) {
317  double plus=histPlus[shift][pidInd]->GetBinContent(bin);
318  double minus=histMinus[shift][pidInd]->GetBinContent(bin);
319  if(plus<minus) {
320  double temp=plus;
321  plus=minus;
322  minus=temp;
323  }
324  plusError+=(plus-nominal)*(plus-nominal);
325  minusError+=(minus-nominal)*(minus-nominal);
326  }
327  plusError=sqrt(plusError);
328  minusError=sqrt(minusError);
329  plusError/=nominal;
330  minusError/=nominal;
331  grErrorBand->SetPoint(bin-1,binCent,0);
332  grErrorBand->SetPointEXlow(bin-1,binCent-binLow);
333  grErrorBand->SetPointEXhigh(bin-1,binHigh-binCent);
334  grErrorBand->SetPointEYlow(bin-1,plusError);
335  grErrorBand->SetPointEYhigh(bin-1,minusError);
336  }
337 
338  return grErrorBand;
339 }
340 
341 
342 TCanvas *getCanvasWithHists(int systNum,const char *ndOrFD,const char *legTag,TH1F *histFullNoShift[3],TH1F *histFullPlus1[3],TH1F *histFullMinus1[3],TH1F *histSignalNoShift[3],TH1F *histSignalPlus1[3],TH1F *histSignalMinus1[3],TH1F *histBackgroundNoShift[3],TH1F *histBackgroundPlus1[3],TH1F *histBackgroundMinus1[3] ) {
343 
344  char legText[180];
345  char canName[180];
346  char padName[180];
347  char ratioName[180];
348  sprintf(canName,"can%s%d",ndOrFD,systNum);
349  TCanvas *can = new TCanvas(canName,canName,1200,800);
350  can->SetTopMargin(0);
351  can->SetRightMargin(0);
352  can->SetLeftMargin(0);
353  can->SetBottomMargin(0);
354 
355  Double_t eps=0.005;
356  TPad *subPad[6];
357 
358 
359  Double_t xVals[4]={0,0.33,0.66,1.0};
360  Double_t yBottom[3]={0.5,0};
361  Double_t yTop[3]={1.0,0.49};
362 
363  for(int row=0;row<2;row++) {
364  for(int col=0;col<3;col++) {
365  sprintf(padName,"can_%d_%d_%d",systNum,row,col);
366  subPad[col+3*row]= new TPad(padName,padName,xVals[col],yBottom[row],xVals[col+1],yTop[row],0);
367  subPad[col+3*row]->Draw();
368  if(row==0)
369  subPad[col+3*row]->SetBottomMargin(0);
370  else
371  subPad[col+3*row]->SetTopMargin(0);
372 
373  if(col!=0)
374  subPad[col+3*row]->SetLeftMargin(0);
375  if(col!=2)
376  subPad[col+3*row]->SetRightMargin(0);
377 
378  }
379  }
380 
381 
382  Double_t histMax=0;
383 
384  for(int pidInd=0;pidInd<NUM_PID_BINS;pidInd++) {
385  if(histFullMinus1[pidInd]->GetMaximum()>histMax)
386  histMax=histFullMinus1[pidInd]->GetMaximum();
387  if(histMax<histFullPlus1[pidInd]->GetMaximum()) {
388  histMax=histFullPlus1[pidInd]->GetMaximum();
389  }
390  if(histMax<histFullNoShift[pidInd]->GetMaximum()) {
391  histMax=histFullNoShift[pidInd]->GetMaximum();
392  }
393  }
394 
395 
396 
397  for(int pidInd=0;pidInd<NUM_PID_BINS;pidInd++) {
398  subPad[pidInd]->cd();
399  histFullMinus1[pidInd]->SetLineColor(getFnexColour(1));
400  histFullMinus1[pidInd]->SetLineStyle(2);
401  histFullMinus1[pidInd]->Draw("hist");
402  histFullMinus1[pidInd]->GetXaxis()->SetRangeUser(0,5);
403  histFullMinus1[pidInd]->GetXaxis()->SetTitle("Reco E. (GeV)");
404  histFullMinus1[pidInd]->GetYaxis()->SetTitle("#Events");
405 
406  histFullNoShift[pidInd]->Draw("hist same");
407  histFullPlus1[pidInd]->SetLineColor(getFnexColour(1));
408  histFullPlus1[pidInd]->Draw("hist same");
409  histFullMinus1[pidInd]->GetYaxis()->SetRangeUser(0,1.1*histMax);
410 
411  if(histSignalNoShift) {
412  histSignalNoShift[pidInd]->SetLineColor(getFnexColour(0));
413  histSignalNoShift[pidInd]->Draw("hist same");
414  histSignalMinus1[pidInd]->SetLineColor(getFnexColour(0));
415  histSignalMinus1[pidInd]->SetLineStyle(2);
416  histSignalMinus1[pidInd]->Draw("same");
417  histSignalPlus1[pidInd]->SetLineColor(getFnexColour(0));
418  histSignalPlus1[pidInd]->SetLineStyle(3);
419  histSignalPlus1[pidInd]->Draw("same");
420  }
421 
422 
423  subPad[3+pidInd]->cd();
424  sprintf(ratioName,"ratio%sFullPlus1_%d_%d",ndOrFD,systNum,pidInd);
425  TH1F *ratioFullPlus1 = (TH1F*)histFullPlus1[pidInd]->Clone(ratioName);
426  ratioFullPlus1->Divide(histFullNoShift[pidInd]);
427  ratioFullPlus1->SetLineWidth(3);
428  ratioFullPlus1->Draw("hist");
429  ratioFullPlus1->GetXaxis()->SetRangeUser(0,5);
430  ratioFullPlus1->GetYaxis()->SetRangeUser(0.5,1.5);
431  ratioFullPlus1->GetXaxis()->SetTitle("Reco. E. (GeV)");
432  ratioFullPlus1->SetYTitle("Ratio to nominal");
433  sprintf(ratioName,"ratio%sFullMinus1_%d",ndOrFD,systNum);
434  TH1F *ratioFullMinus1 = (TH1F*)histFullMinus1[pidInd]->Clone(ratioName);
435  ratioFullMinus1->Divide(histFullNoShift[pidInd]);
436  ratioFullMinus1->SetLineWidth(3);
437  ratioFullMinus1->Draw("same hist");
438 
439  if(histSignalNoShift) {
440  sprintf(ratioName,"ratio%sSignalPlus1_%d_%d",ndOrFD,systNum,pidInd);
441  TH1F *ratioSignalPlus1 = (TH1F*)histSignalPlus1[pidInd]->Clone(ratioName);
442  ratioSignalPlus1->Divide(histSignalNoShift[pidInd]);
443  ratioSignalPlus1->SetLineWidth(3);
444  ratioSignalPlus1->Draw("hist same");
445  sprintf(ratioName,"ratio%sSignalMinus1_%d",ndOrFD,systNum);
446  TH1F *ratioSignalMinus1 = (TH1F*)histSignalMinus1[pidInd]->Clone(ratioName);
447  ratioSignalMinus1->Divide(histSignalNoShift[pidInd]);
448  ratioSignalMinus1->SetLineWidth(3);
449  ratioSignalMinus1->Draw("same hist");
450 
451  }
452 
453 
454 
455  // can->cd(1+pidInd);
456  if(pidInd==0) {
457  subPad[pidInd]->cd();
458  TLegend *leggy = makeLegendBoxInTopRightCorner(0.8,0.5);
459  sprintf(legText,"%s Nominal MC",ndOrFD);
460 
461  leggy->AddEntry(histFullNoShift[pidInd],legText,"l");
462  sprintf(legText,"%s -- %s +1#sigma",ndOrFD,legTag);
463  leggy->AddEntry(histFullPlus1[pidInd],legText,"l");
464  sprintf(legText,"%s -- %s -1#sigma",ndOrFD,legTag);
465  leggy->AddEntry(histFullMinus1[pidInd],legText,"l");
466 
467  if(histSignalNoShift) {
468  sprintf(legText,"%s Nominal #nu_{#mu}#rightarrow#nu_{e}",ndOrFD);
469  leggy->AddEntry(histSignalNoShift[pidInd],legText,"l");
470  sprintf(legText,"%s -- %s +1#sigma",ndOrFD,legTag);
471  leggy->AddEntry(histSignalPlus1[pidInd],legText,"l");
472  sprintf(legText,"%s -- %s -1#sigma",ndOrFD,legTag);
473  leggy->AddEntry(histSignalMinus1[pidInd],legText,"l");
474  }
475 
476 
477 
478  leggy->Draw();
479  }
480  }
481  PrintPDFandPNG(can,canName);
482  return can;
483 
484 }
485 
486 
487 
488 
489 
490 
enum BeamMode kOrange
enum BeamMode kRed
TGraphAsymmErrors * getErrorBand(TH1F *histNoShift[3], TH1F *histPlus[1024][3], TH1F *histMinus[1024][3], Int_t numShifts, Int_t pidInd)
T sqrt(T number)
Definition: d0nt_math.hpp:156
void plotNueSystBandForTheseSystematics(char *rootFileList, char *legTagList=0)
const int NUM_PID_BINS
std::string GetName(int i)
Int_t col[ntarg]
Definition: Style.C:29
TGraphAsymmErrors * getFractionalErrorBand(TH1F *histNoShift[3], TH1F *histPlus[1024][3], TH1F *histMinus[1024][3], Int_t numShifts, Int_t pidInd)
int getFnexColour(int index)
float bin[41]
Definition: plottest35.C:14
OStream cout
Definition: OStream.cxx:6
def countFiles(dir, testStrs=[])
enum BeamMode kViolet
TCanvas * getCanvasWithHists(int systNum, const char *ndOrFD, const char *legTag, TH1F *histFullNoShift[3], TH1F *histFullPlus1[3], TH1F *histFullMinus1[3], TH1F *histSignalNoShift[3], TH1F *histSignalPlus1[3], TH1F *histSignalMinus1[3], TH1F *histBackgroundNoShift[3], TH1F *histBackgroundPlus1[3], TH1F *histBackgroundMinus1[3])
enum BeamMode kGreen
enum BeamMode kBlue
T minus(const T &x)
Definition: minus.hpp:15