Plotting_DataAndPrediction.C
Go to the documentation of this file.
1 #pragma once
2 
6 
7 using namespace ana;
8 
9 //STYLING.
10 Color_t color_band = kViolet-9;
11 Color_t color_pred = kViolet+1;
12 Color_t color_ws = kGreen+2;
13 Color_t color_nc = kBlue-7;
14 Color_t color_data = kBlack;
15 Color_t color_cos = kAzure+1;
16 Color_t color_bkg = kGray+1;
17 Color_t color_noosc = kRed+1;
18 Color_t color_noosc_b = kRed-9;
19 Style_t line_ws = kSolid;
20 Style_t line_nc = kSolid;
21 Style_t line_data = kSolid;
22 Style_t line_pred = kSolid;
23 Style_t line_cos = kSolid;
24 Style_t line_bkg = kSolid;
25 Style_t marker_ws = kFullCircle;
26 Style_t marker_nc = kFullSquare;
27 Style_t marker_data = kFullCircle;
28 Style_t marker_pred = kFullCircle;
29 Style_t marker_cos = kFullTriangleUp;
30 Style_t marker_bkg = kFullCircle;
31 
32 void Plotting_DataAndPrediction(std::string s_HornCurrent = "fhc",
33  bool drawSystBand = true,
34  bool makeUnoscPlots = true,
35  std::string s_FileAppend = "",
36  double dmsq = 0.,
37  double th23 = 0.,
38  bool isPreliminary = true)
39 {
40  std::string s_OutDir = "";
41  std::string s_Calc = "/nova/ana/nu_e_ana/Ana2019/Results/BestFits/bestfits_joint_realData_both_systs.root";
42  std::string s_HiOct = "NH_UO_calc";
43  std::string s_CosDir = "/nova/ana/nu_mu_ana/Ana2019/Cosmics/";
44  std::string s_PredDir = "/nova/ana/nu_mu_ana/Ana2018/Predictions/";
45  std::string s_DataDir = "/nova/ana/nu_mu_ana/Ana2019/Data/";
46  std::string s_DataHist = s_DataDir;
47 
48  s_FileAppend += "_" + s_HornCurrent;
49  if(drawSystBand)
50  {
51  s_FileAppend += "_SystBand";
52  }
53 
54  unsigned int nQuantiles = 4;
55  double pot;
56  double livetime;
57  Sign::Sign_t WS;
58 
59  if(s_HornCurrent=="fhc")
60  {
61  pot = kAna2019FHCPOT;
62  livetime = kAna2019FHCLivetime;
63  WS = Sign::kAntiNu;
64  s_DataDir += "fd_spectra_TGraphs_FHC2019.root";
65  s_DataHist += "fd_dataspectrum_fhc_full__numu2019.root";
66  }
67  else
68  {
69  pot = kAna2019RHCPOT;
70  livetime = kAna2019RHCLivetime;
71  WS = Sign::kNu;
72  s_DataDir += "fd_spectra_TGraphs_RHC2019.root";
73  s_DataHist += "fd_dataspectrum_rhc_full__numu2019.root";
74  }
75 
77 
78  auto calc_Def = DefaultOscCalc();
79 
80  osc::NoOscillations base_NoOsc;
81  auto calc_NoOsc = base_NoOsc.Copy();
82 
84  calc_User->SetTh23(th23);
85  calc_User->SetDmsq32(dmsq);
86 
87  TFile *f_Calc = new TFile(s_Calc.c_str(), "OPEN");
88  auto calc_File = (osc::IOscCalcAdjustable*)(LoadFrom<osc::IOscCalc>(f_Calc, s_HiOct.c_str())).release();
89  auto bfshiftsFile = LoadFrom<SystShifts>(f_Calc, "NH_UO_systs").release();
90 
91  auto calc = ((dmsq==0. && th23==0.) ? calc_File->Copy() : calc_User->Copy());
92 
93  std::string s_Year = "2019";
94  std::string fcosmics_name_one = s_HornCurrent + "cosmics_" + s_HornCurrent + "__numu" + s_Year + ".root";
95  TFile *f_Cosmics = new TFile((s_CosDir+"cosmics_" + s_HornCurrent + "__numu" + s_Year + ".root").c_str(), "READ");
96  std::vector<TH1*> hCosmics;
97  for(unsigned int i = 0; i < nQuantiles+1; i++)
98  {
99  hCosmics.push_back((TH1*)f_Cosmics->Get(Form("cosmics_q%i", i)));
100  }
101 
102  std::string s_Pred = s_PredDir + "pred_numuconcat_" + s_HornCurrent + "__numu2018.root";
103 
104  SystShifts bfshifts;
105  for(auto syst: bfshiftsFile->ActiveSysts()){
106  if(std::find(systs.begin(), systs.end(), syst) != systs.end())
107  bfshifts.SetShift(syst, bfshiftsFile->GetShift(syst));
108  }
109 
110  std::vector<PredictionSystJoint2018*> vec_Pred;
111  ENu2018ExtrapType extrapType = kNuMu;
112  for(unsigned int i = 0; i < nQuantiles; i++)
113  {
114  vec_Pred.push_back(new PredictionSystJoint2018(extrapType, calc_Def, s_HornCurrent, i+1, systs, s_Pred));
115  }
116 
117  std::vector<TH1*> hCC;
118  std::vector<TH1*> hNuMu;
119  std::vector<TH1*> hNC;
120  std::vector<TH1*> hWS;
121  std::vector<TH1*> hNoOsc;
122  std::vector<TH1*> hPred;
123  std::vector<TH1*> hBkg;
124  std::vector<TH1*> hAllBkg;
125  std::vector<TH1*> hNoOscNC;
126  std::vector<TH1*> hNoOscNumu;
127 
128  for(unsigned int i = 0; i < nQuantiles; i++)
129  {
130  hCC.push_back (vec_Pred[i]->PredictComponentSyst(calc, bfshifts, Flavors::kAll, Current::kCC, Sign::kBoth).ToTH1(pot));
131  hNuMu.push_back (vec_Pred[i]->PredictComponentSyst(calc, bfshifts, Flavors::kAllNuMu, Current::kCC, Sign::kBoth).ToTH1(pot));
132  hNC.push_back (vec_Pred[i]->PredictComponentSyst(calc, bfshifts, Flavors::kAll, Current::kNC, Sign::kBoth).ToTH1(pot));
133  hWS.push_back (vec_Pred[i]->PredictComponentSyst(calc, bfshifts, Flavors::kAllNuMu, Current::kCC, WS).ToTH1(pot));
134  hNoOsc.push_back (vec_Pred[i]->Predict(calc_NoOsc).ToTH1(pot));
135  hPred.push_back (vec_Pred[i]->PredictSyst(calc, bfshifts).ToTH1(pot));
136  hBkg.push_back (vec_Pred[i]->PredictSyst(calc, bfshifts).ToTH1(pot));
137  hAllBkg.push_back(vec_Pred[i]->PredictSyst(calc, bfshifts).ToTH1(pot));
138 
139  hNoOscNC.push_back (vec_Pred[i]->PredictComponent(calc_NoOsc, Flavors::kAll, Current::kNC, Sign::kBoth).ToTH1(pot));
140  hNoOscNumu.push_back(vec_Pred[i]->PredictComponent(calc_NoOsc, Flavors::kAllNuMu, Current::kCC, Sign::kBoth).ToTH1(pot));
141  if(i==0)
142  {
143  hCC.push_back (vec_Pred[i]->PredictComponentSyst(calc, bfshifts, Flavors::kAll, Current::kCC, Sign::kBoth).ToTH1(pot));
144  hNuMu.push_back (vec_Pred[i]->PredictComponentSyst(calc, bfshifts, Flavors::kAllNuMu, Current::kCC, Sign::kBoth).ToTH1(pot));
145  hNC.push_back (vec_Pred[i]->PredictComponentSyst(calc, bfshifts, Flavors::kAll, Current::kNC, Sign::kBoth).ToTH1(pot));
146  hWS.push_back (vec_Pred[i]->PredictComponentSyst(calc, bfshifts, Flavors::kAllNuMu, Current::kCC, WS).ToTH1(pot));
147  hNoOsc.push_back (vec_Pred[i]->Predict(calc_NoOsc).ToTH1(pot));
148  hPred.push_back (vec_Pred[i]->PredictSyst(calc, bfshifts).ToTH1(pot));
149  hBkg.push_back (vec_Pred[i]->PredictSyst(calc, bfshifts).ToTH1(pot));
150  hAllBkg.push_back(vec_Pred[i]->PredictSyst(calc, bfshifts).ToTH1(pot));
151 
152  hNoOscNC.push_back (vec_Pred[i]->PredictComponent(calc_NoOsc, Flavors::kAll, Current::kNC, Sign::kBoth).ToTH1(pot));
153  hNoOscNumu.push_back(vec_Pred[i]->PredictComponent(calc_NoOsc, Flavors::kAllNuMu, Current::kCC, Sign::kBoth).ToTH1(pot));
154  }
155  else
156  {
157  hCC [0]->Add(vec_Pred[i]->PredictComponentSyst(calc, bfshifts, Flavors::kAll, Current::kCC, Sign::kBoth).ToTH1(pot));
158  hNuMu [0]->Add(vec_Pred[i]->PredictComponentSyst(calc, bfshifts, Flavors::kAllNuMu, Current::kCC, Sign::kBoth).ToTH1(pot));
159  hNC [0]->Add(vec_Pred[i]->PredictComponentSyst(calc, bfshifts, Flavors::kAll, Current::kNC, Sign::kBoth).ToTH1(pot));
160  hWS [0]->Add(vec_Pred[i]->PredictComponentSyst(calc, bfshifts, Flavors::kAllNuMu, Current::kCC, WS).ToTH1(pot));
161  hNoOsc [0]->Add(vec_Pred[i]->Predict(calc_NoOsc).ToTH1(pot));
162  hPred [0]->Add(vec_Pred[i]->PredictSyst(calc, bfshifts).ToTH1(pot));
163  hBkg [0]->Add(vec_Pred[i]->PredictSyst(calc, bfshifts).ToTH1(pot));
164  hAllBkg[0]->Add(vec_Pred[i]->PredictSyst(calc, bfshifts).ToTH1(pot));
165 
166  hNoOscNC [0]->Add(vec_Pred[i]->PredictComponent(calc_NoOsc, Flavors::kAll, Current::kNC, Sign::kBoth).ToTH1(pot));
167  hNoOscNumu[0]->Add(vec_Pred[i]->PredictComponent(calc_NoOsc, Flavors::kAllNuMu, Current::kCC, Sign::kBoth).ToTH1(pot));
168  }
169  }
170 
171  //FAKE STACK.
172  for(unsigned int i = 0; i < hBkg.size(); i++){
173  hBkg[i] ->Add(hNuMu [i], -1);
174  hAllBkg[i]->Add(hNuMu [i], -1);
175  hAllBkg[i]->Add(hCosmics[i]);
176  hWS[i] ->Add(hAllBkg [i]);
177  hPred[i] ->Add(hCosmics[i]);
178  hNoOsc[i] ->Add(hCosmics[i]);
179  }
180 
181  for(unsigned int i = 0; i < hBkg.size(); i++)
182  {
183  hCosmics[i]->Scale(0.1, "width");
184  hCC[i] ->Scale(0.1, "width");
185  hNuMu[i] ->Scale(0.1, "width");
186  hNC[i] ->Scale(0.1, "width");
187  hWS[i] ->Scale(0.1, "width");
188  hNoOsc[i] ->Scale(0.1, "width");
189  hPred[i] ->Scale(0.1, "width");
190  hBkg[i] ->Scale(0.1, "width");
191  hAllBkg[i] ->Scale(0.1, "width");
192  }
193 
194  TFile *f_Data = new TFile(s_DataDir.c_str(), "OPEN");
195  TFile *f_DataHist = new TFile(s_DataHist.c_str(), "OPEN");
196  std::vector<TGraphAsymmErrors*> vec_Data;
197  std::vector<Spectrum> spec_Pred;
198  std::vector<Spectrum> spec_NoOsc;
199  for(unsigned int i = 0; i < hPred.size(); i++)
200  {
201  Spectrum spec(hPred[i], pot, livetime);
202  spec_Pred.push_back(spec);
203  Spectrum spec_NO(hNoOsc[i], pot, livetime);
204  spec_NoOsc.push_back(spec_NO);
205  vec_Data.push_back((TGraphAsymmErrors*)f_Data->Get(Form("g1_fd_spectra_Q%i", i)));
206  vec_Data.at(i)->SetMarkerStyle(kFullCircle);
207  }
208 
209  std::vector<Spectrum> ups[nQuantiles+1], dns[nQuantiles+1];
210  std::vector<Spectrum> upsNoOsc[nQuantiles+1], dnsNoOsc[nQuantiles+1];
211 
212  if(drawSystBand)
213  {
214  for(unsigned int i = 0; i < hPred.size(); i++)
215  {
216  for(const ISyst* syst: systs){
217  SystShifts shifts;
218  if(i==0)
219  {
220  shifts.SetShift(syst, +1);
221  TH1 * htempup = vec_Pred[0]->PredictSyst(calc, shifts).ToTH1(pot);
222  htempup ->Add(vec_Pred[0]->PredictSyst(calc, SystShifts::Nominal()).ToTH1(pot), -1);
223  htempup ->Add(vec_Pred[0]->PredictSyst(calc, bfshifts).ToTH1(pot));
224  TH1 * htempupNoOsc = vec_Pred[0]->PredictSyst(calc_NoOsc, shifts).ToTH1(pot);
225  for(unsigned int j = 1; j < nQuantiles; j++)
226  {
227  htempup ->Add(vec_Pred[j]->PredictSyst(calc, shifts).ToTH1(pot));
228  htempup ->Add(vec_Pred[j]->PredictSyst(calc, SystShifts::Nominal()).ToTH1(pot), -1);
229  htempup ->Add(vec_Pred[j]->PredictSyst(calc, bfshifts).ToTH1(pot));
230  htempupNoOsc->Add(vec_Pred[j]->PredictSyst(calc_NoOsc, shifts).ToTH1(pot));
231  }
232  htempup->Scale(0.1, "width");
233  htempup->Add(hCosmics[0]);
234  htempupNoOsc->Scale(0.1, "width");
235 
236  Spectrum specup(htempup, pot, livetime);
237  ups[i].push_back(specup);
238  Spectrum specupNoOsc(htempupNoOsc, pot, livetime);
239  upsNoOsc[i].push_back(specupNoOsc);
240 
241  shifts.SetShift(syst, -1);
242  TH1 * htempdn = vec_Pred[0]->PredictSyst(calc, shifts).ToTH1(pot);
243  htempdn ->Add(vec_Pred[0]->PredictSyst(calc, SystShifts::Nominal()).ToTH1(pot), -1);
244  htempdn ->Add(vec_Pred[0]->PredictSyst(calc, bfshifts).ToTH1(pot));
245  TH1 * htempdnNoOsc = vec_Pred[0]->PredictSyst(calc_NoOsc, shifts).ToTH1(pot);
246  for(unsigned int j = 1; j < nQuantiles; j++)
247  {
248  htempdn ->Add(vec_Pred[j]->PredictSyst(calc, shifts).ToTH1(pot));
249  htempdn ->Add(vec_Pred[j]->PredictSyst(calc, SystShifts::Nominal()).ToTH1(pot), -1);
250  htempdn ->Add(vec_Pred[j]->PredictSyst(calc, bfshifts).ToTH1(pot));
251  htempdnNoOsc->Add(vec_Pred[j]->PredictSyst(calc_NoOsc, shifts).ToTH1(pot));
252  }
253  htempdn->Scale(0.1, "width");
254  htempdn->Add(hCosmics[0]);
255  htempdnNoOsc->Scale(0.1, "width");
256 
257  Spectrum specdn(htempdn, pot, livetime);
258  dns[i].push_back(specdn);
259  Spectrum specdnNoOsc(htempdnNoOsc, pot, livetime);
260  dnsNoOsc[i].push_back(specdnNoOsc);
261  }
262  else
263  {
264  shifts.SetShift(syst, +1);
265  TH1* htempup = vec_Pred[i-1]->PredictSyst(calc, shifts).ToTH1(pot);
266  htempup ->Add(vec_Pred[i-1]->PredictSyst(calc, SystShifts::Nominal()).ToTH1(pot), -1);
267  htempup ->Add(vec_Pred[i-1]->PredictSyst(calc, bfshifts).ToTH1(pot));
268  TH1* htempupNoOsc = vec_Pred[i-1]->PredictSyst(calc_NoOsc, shifts).ToTH1(pot);
269  htempup->Scale(0.1, "width");
270  htempup->Add(hCosmics[i]);
271  htempupNoOsc->Scale(0.1, "width");
272  Spectrum specup(htempup, pot, livetime);
273  ups[i].push_back(specup);
274  Spectrum specupNoOsc(htempupNoOsc, pot, livetime);
275  upsNoOsc[i].push_back(specupNoOsc);
276  shifts.SetShift(syst, -1);
277  TH1 * htempdn = vec_Pred[i-1]->PredictSyst(calc, shifts).ToTH1(pot);
278  htempdn ->Add(vec_Pred[i-1]->PredictSyst(calc, SystShifts::Nominal()).ToTH1(pot), -1);
279  htempdn ->Add(vec_Pred[i-1]->PredictSyst(calc, bfshifts).ToTH1(pot));
280  TH1 * htempdnNoOsc = vec_Pred[i-1]->PredictSyst(calc_NoOsc, shifts).ToTH1(pot);
281  htempdn->Scale(0.1, "width");
282  htempdn->Add(hCosmics[i]);
283  htempdnNoOsc->Scale(0.1, "width");
284  Spectrum specdn(htempdn, pot, livetime);
285  dns[i].push_back(specdn);
286  Spectrum specdnNoOsc(htempdnNoOsc, pot, livetime);
287  dnsNoOsc[i].push_back(specdnNoOsc);
288  }
289  }
290  }
291  }
292 
293  DecorateHist(hWS [0], line_ws, color_ws, 1, marker_ws, color_ws, 1);
294  DecorateHist(hNC [0], line_nc, color_nc, 1, marker_nc, color_nc, 1);
296  DecorateHist(hCosmics[0], line_cos, color_cos, 1, marker_cos, color_cos, 1);
297  DecorateHist(hAllBkg [0], line_bkg, color_bkg, 1, marker_bkg, color_bkg, 1);
298  DecorateHist(hBkg [0], line_bkg, color_bkg, 1, marker_bkg, color_bkg, 1);
300  FillWithDimColor(hWS [0], 0);
301  FillWithDimColor(hCosmics[0], 0);
302  FillWithDimColor(hAllBkg [0], 0);
303 
304  double maxPredByQ = 0.;
305  double maxDataByQ = 0.;
306  double maxNoOscByQ = 0.;
307  double maxQ = 0.;
308  for(unsigned int i = 1; i < hBkg.size(); i++)
309  {
310  DecorateHist(hWS [i], line_ws, color_ws, 2, marker_ws, color_ws, 2);
311  DecorateHist(hNC [i], line_nc, color_nc, 2, marker_nc, color_nc, 2);
313  DecorateHist(hCosmics[i], line_cos, color_cos, 2, marker_cos, color_cos, 2);
314  DecorateHist(hAllBkg [i], line_bkg, color_bkg, 2, marker_bkg, color_bkg, 2);
315  DecorateHist(hBkg [i], line_bkg, color_bkg, 2, marker_bkg, color_bkg, 2);
317  FillWithDimColor(hWS [i], 0);
318  FillWithDimColor(hCosmics[i], 0);
319  FillWithDimColor(hAllBkg [i], 0);
320  if(hPred[i]->GetMaximum()>maxPredByQ)
321  {
322  maxPredByQ = hPred[i]->GetMaximum();
323  }
324  if(hNoOsc[i]->GetMaximum()>maxNoOscByQ)
325  {
326  maxNoOscByQ = hNoOsc[i]->GetMaximum();
327  }
328  for(int j = 0; j < vec_Data[i]->GetN(); j++)
329  {
330  double x, y;
331  vec_Data[i]->GetPoint(j, x, y);
332  double ey = vec_Data[i]->GetErrorYhigh(j);
333  if((y+ey)>maxDataByQ)
334  {
335  maxDataByQ = y + ey;
336  }
337  }
338  }
339 
340  maxQ = std::max(maxPredByQ, maxDataByQ);
341 
342  TFile *f_Out = new TFile((s_OutDir+"dataAndPrediction_" + s_FileAppend + ".root").c_str(), "RECREATE");
343  TCanvas *c = new TCanvas("c_dataAndPrediction_", "c_dataAndPrediction_", 600, 500);
344  TCanvas *c_ByQ = new TCanvas("c_dataAndPrediction_ByQ","c_dataAndPrediction_ByQ", 600, 500);
345  TPad *p1, *p2, *p3, *p4;
346  SplitCanvasQuant(c_ByQ, p1, p2, p3, p4);
347 
348  TLatex* xTitle = new TLatex(0.5, 0.03, "Reconstructed Neutrino Energy (GeV)");
349  xTitle->SetNDC();
350  xTitle->SetTextAlign(22);
351  TLatex* yTitle = new TLatex(0.02, 0.5, "Events / 0.1 GeV");
352  yTitle->SetNDC();
353  yTitle->SetTextAlign(22);
354  yTitle->SetTextAngle(90.);
355 
356  std::string s_HornLabel = s_HornCurrent=="rhc" ? "Antineutrino beam" : "Neutrino beam";
357 
358  TLegend *leg = new TLegend(0.58,0.50,0.78,0.9);
359  leg->AddEntry(vec_Data[0], "FD Data", "LEP");
360  leg->AddEntry(hPred [0], "Prediction", "L" );
361  TLegendEntry *ent_Band;
362  if(drawSystBand)
363  {
364  ent_Band =leg->AddEntry("error","1-#sigma syst. range", "BF");
365  ent_Band->SetFillStyle(1001);
366  ent_Band->SetFillColor(color_band);
367  ent_Band->SetLineColor(color_band);
368  }
369  leg->AddEntry(hWS [0], (s_HornCurrent=="rhc" ? "Wrong Sign: #nu_{#mu}CC" : "Wrong Sign: #bar{#nu}_{#mu}CC"),"BF");
370  leg->AddEntry(hAllBkg [0], "Total bkgd.", "BF" );
371  leg->AddEntry(hCosmics[0], "Cosmic bkgd.", "BF" );
372  leg->SetTextSize(0.05);
373  leg->SetBorderSize(0);
374  leg->SetFillStyle(0);
375 
376  c->cd();
377  ana::CenterTitles(hPred[0]);
378  if(drawSystBand)
379  {
380  TH1 *nom = spec_Pred[0].ToTH1(pot);
381  nom->GetYaxis()->SetTitle("Events / 0.1 GeV");
382  std::vector<TH1*> upsh, dnsh;
383  for(const auto& upShift: ups[0]) upsh.push_back(upShift.ToTH1(pot));
384  for(const auto& downShift: dns[0]) dnsh.push_back(downShift.ToTH1(pot));
385  PlotWithSystErrorBand(nom, upsh, dnsh, color_pred, color_band, 1.6);
386  }
387  else
388  {
389  hPred[0]->Draw("HIST");
390  hPred[0]->GetYaxis()->SetRangeUser(0, hPred[0]->GetMaximum()*1.6);
391  hPred[0]->GetYaxis()->SetTitle("Events / 0.1 GeV");
392  }
393  hWS [0]->Draw("HIST SAME");
394  hAllBkg [0]->Draw("HIST SAME");
395  hCosmics[0]->Draw("HIST SAME");
396  vec_Data[0]->Draw("P");
397 
399  MyCornerLabel(s_HornLabel);
400  drawQuantLabel(0);
401  leg->Draw();
402 
403  c->SetName(c->GetName()+(TString)"Spectra");
404  c->SaveAs((TString)s_OutDir + c->GetName() + (TString)s_FileAppend + ".pdf");
405  c->Write();
406 
407  c_ByQ->cd();
408 
409  for(unsigned int i = 1; i < hBkg.size(); i++)
410  {
411  if(i==1) p1->cd();
412  if(i==2) p2->cd();
413  if(i==3) p3->cd();
414  if(i==4) p4->cd();
415  if(drawSystBand)
416  {
417  TH1 *nom = spec_Pred[i].ToTH1(pot);
418  if(i==3)
419  {
420  nom->GetXaxis()->ChangeLabel(6,-1,0);
421  }
422  std::vector<TH1*> upsh, dnsh;
423  for(const auto& upShift: ups[i]) upsh.push_back(upShift.ToTH1(pot));
424  for(const auto& downShift: dns[i]) dnsh.push_back(downShift.ToTH1(pot));
425  //PlotWithSystErrorBand_Quant(i, nom, upsh, dnsh, color_pred, color_band, maxQ*1.25);
426  PlotWithSystErrorBand_Quant(i, nom, upsh, dnsh, color_pred, color_band, 5.99);
427  }
428  else
429  {
430  MakeHistCanvasReady_Quant(i, hPred[i], 5.99);
431  hPred[i]->Draw("HIST");
432  if(i==3)
433  {
434  hPred[i]->GetXaxis()->ChangeLabel(6,-1,0);
435  }
436  }
437  hWS [i]->Draw("HIST SAME");
438  hAllBkg [i]->Draw("HIST SAME");
439  hCosmics[i]->Draw("HIST SAME");
440  vec_Data[i]->Draw("P");
441  drawQuantLabel(i);
442  leg->Draw();
443  }
444 
445  c_ByQ->cd();
447  MyCornerLabel(s_HornLabel);
448  xTitle->Draw();
449  yTitle->Draw();
450 
451  c_ByQ->SetName("c_dataAndPrediction_ByQ_Spectra");
452  c_ByQ->SaveAs((TString)s_OutDir + c_ByQ->GetName() + (TString)s_FileAppend + ".pdf");
453  c_ByQ->Write();
454 
455  if(makeUnoscPlots)
456  {
457  TCanvas *c_NoOsc = new TCanvas("c_dataAndPrediction_NoOsc", "c_dataAndPrediction_NoOsc", 600, 500);
458  TCanvas *c_ByQ_NoOsc = new TCanvas("c_dataAndPrediction_ByQ_NoOsc","c_dataAndPrediction_ByQ_NoOsc", 600, 500);
459  TPad *p1_NoOsc, *p2_NoOsc, *p3_NoOsc, *p4_NoOsc;
460  SplitCanvasQuant(c_ByQ_NoOsc, p1_NoOsc, p2_NoOsc, p3_NoOsc, p4_NoOsc);
461 
462  TLegend *leg_NoOsc = new TLegend(0.55, 0.40, 0.78, 0.85);
463  leg_NoOsc->AddEntry(vec_Data[0], "FD Data", "LEP");
464  leg_NoOsc->AddEntry(hPred [0], "Prediction", "L" );
465  TLegendEntry *ent_Band;
466  TLegendEntry *ent_Band_NoOsc;
467  if(drawSystBand)
468  {
469  /*
470  ent_Band =leg_NoOsc->AddEntry("error","1-#sigma syst. range", "BF");
471  ent_Band->SetFillStyle(1001);
472  ent_Band->SetFillColor(color_band);
473  ent_Band->SetLineColor(color_band);
474  */
475  leg_NoOsc->AddEntry(hNoOsc[0], "No oscillation", "L" );
476  ent_Band_NoOsc =leg_NoOsc->AddEntry("errorNoOsc","1-#sigma syst. range", "BF");
477  ent_Band_NoOsc->SetFillStyle(1001);
478  ent_Band_NoOsc->SetFillColor(color_noosc_b);
479  ent_Band_NoOsc->SetLineColor(color_noosc_b);
480  }
481  else
482  {
483  leg_NoOsc->AddEntry(hNoOsc[0], "No oscillation", "L" );
484  }
485  leg_NoOsc->AddEntry(hAllBkg[0], "Total bkgd.", "BF");
486 
487  leg_NoOsc->SetTextSize(0.05);
488  leg_NoOsc->SetBorderSize(0);
489  leg_NoOsc->SetFillStyle(0);
490 
491  c_NoOsc->cd();
492  ana::CenterTitles(hNoOsc[0]);
493  if(drawSystBand)
494  {
495  TH1 *nomNoOsc = spec_NoOsc[0].ToTH1(pot);
496  nomNoOsc->GetYaxis()->SetTitle("Events / 0.1 GeV");
497  std::vector<TH1*> upshNoOsc, dnshNoOsc;
498  for(const auto& upShift: upsNoOsc[0]) upshNoOsc.push_back(upShift .ToTH1(pot));
499  for(const auto& downShift: dnsNoOsc[0]) dnshNoOsc.push_back(downShift.ToTH1(pot));
500  PlotWithSystErrorBand(nomNoOsc, upshNoOsc, dnshNoOsc, color_noosc, color_noosc_b, 1.3);
501  TH1 *nom = spec_Pred[0].ToTH1(pot);
502  std::vector<TH1*> upsh, dnsh;
503  for(const auto& upShift: ups[0]) upsh.push_back(upShift.ToTH1(pot));
504  for(const auto& downShift: dns[0]) dnsh.push_back(downShift.ToTH1(pot));
505  PlotWithSystErrorBand(nom, upsh, dnsh, color_pred, color_band, 1., false);
506  }
507  else
508  {
509  hNoOsc[0]->Draw("HIST");
510  hNoOsc[0]->GetYaxis()->SetRangeUser(0, hNoOsc[0]->GetMaximum()*1.15);
511  hNoOsc[0]->GetYaxis()->SetTitle("Events / 0.1 GeV");
512  }
513  hAllBkg [0]->Draw("HIST SAME");
514  hPred [0]->Draw("HIST SAME");
515  vec_Data[0]->Draw("P");
516 
518  MyCornerLabel(s_HornLabel);
519  drawQuantLabel(0);
520  leg_NoOsc->Draw();
521 
522  c_NoOsc->SetName(c_NoOsc->GetName()+(TString)"_Spectra");
523  c_NoOsc->SaveAs((TString)s_OutDir + c_NoOsc->GetName() + (TString)s_FileAppend + ".pdf");
524  c_NoOsc->Write();
525 
526  c_ByQ_NoOsc->cd();
527 
528  for(unsigned int i = 1; i < hBkg.size(); i++)
529  {
530  if(i==1) p1_NoOsc->cd();
531  if(i==2) p2_NoOsc->cd();
532  if(i==3) p3_NoOsc->cd();
533  if(i==4) p4_NoOsc->cd();
534  if(drawSystBand)
535  {
536  TH1 *nom = spec_Pred [i].ToTH1(pot);
537  TH1 *nomNoOsc = spec_NoOsc[i].ToTH1(pot);
538  if(i==3)
539  {
540  nomNoOsc->GetXaxis()->ChangeLabel(6,-1,0);
541  }
542  std::vector<TH1*> upsh, dnsh;
543  for(const auto& upShift: ups[i]) upsh.push_back(upShift.ToTH1(pot));
544  for(const auto& downShift: dns[i]) dnsh.push_back(downShift.ToTH1(pot));
545  std::vector<TH1*> upshNoOsc, dnshNoOsc;
546  for(const auto& upShift: upsNoOsc[i]) upshNoOsc.push_back(upShift .ToTH1(pot));
547  for(const auto& downShift: dnsNoOsc[i]) dnshNoOsc.push_back(downShift.ToTH1(pot));
548  PlotWithSystErrorBand_Quant(i, nomNoOsc, upshNoOsc, dnshNoOsc, color_noosc, color_noosc_b, 15.99);
549  PlotWithSystErrorBand_Quant(i, nom, upsh, dnsh, color_pred, color_band, 15.99, false);
550  //PlotWithSystErrorBand_Quant(i, nomNoOsc, upshNoOsc, dnshNoOsc, color_noosc, color_noosc_b, maxNoOscByQ*1.3);
551  //PlotWithSystErrorBand_Quant(i, nom, upsh, dnsh, color_pred, color_band, maxNoOscByQ*1.3, false);
552  }
553  else
554  {
555  //MakeHistCanvasReady_Quant(i, hNoOsc[i], maxNoOscByQ*1.14);
556  MakeHistCanvasReady_Quant(i, hNoOsc[i], 27.1);
557  hNoOsc[i]->Draw("HIST");
558  if(i==3)
559  {
560  hNoOsc[i]->GetXaxis()->ChangeLabel(6,-1,0);
561  }
562  }
563  hPred [i]->Draw("HIST SAME");
564  hAllBkg [i]->Draw("HIST SAME");
565  vec_Data[i]->Draw("P");
566  drawQuantLabel(i);
567  leg_NoOsc->Draw();
568  }
569 
570  c_ByQ_NoOsc->cd();
572  MyCornerLabel(s_HornLabel);
573  xTitle->Draw();
574  yTitle->Draw();
575 
576  c_ByQ_NoOsc->SetName("c_dataAndPrediction_ByQ_NoOsc_Spectra");
577  c_ByQ_NoOsc->SaveAs((TString)s_OutDir + c_ByQ_NoOsc->GetName() + (TString)s_FileAppend + ".pdf");
578  c_ByQ_NoOsc->Write();
579 
580  TCanvas *c_Rat = new TCanvas("c_dataAndPrediction_RatioNoOsc", "c_dataAndPrediction_RatioNoOsc", 600, 500);
581  TCanvas *c_ByQ_Rat = new TCanvas("c_dataAndPrediction_ByQ_RatioNoOsc","c_dataAndPrediction_ByQ_RatioNoOsc", 600, 500);
582  TPad *p1_Rat, *p2_Rat, *p3_Rat, *p4_Rat;
583  SplitCanvasQuant(c_ByQ_Rat, p1_Rat, p2_Rat, p3_Rat, p4_Rat);
584 
585  std::vector<TH1*> hDataUnscaled;
586  for(unsigned int i = 0; i <= nQuantiles; i++)
587  {
588  std::unique_ptr<Spectrum> spec = Spectrum::LoadFrom(f_DataHist, Form("fd_data_q%i", i));
589  hDataUnscaled.push_back(std::move(spec)->ToTH1(std::move(spec)->POT()));
590  hDataUnscaled[i]->Add(hBkg[i], -1);
591  hDataUnscaled[i]->Add(hCosmics[i], -1);
592  }
593 
594  std::vector<TH1*> hNoOscDenom;
595  std::vector<TH1*> hPredRatio;
596  std::vector<TGraphAsymmErrors*> vec_DataRatio = vec_Data;
597 
598  std::vector<TH1*> upsNoOscRatio[nQuantiles+1], dnsNoOscRatio[nQuantiles+1];
599  for(unsigned int i = 0; i < nQuantiles+1; i++)
600  {
601  for(unsigned int j = 0; j < upsNoOsc[i].size(); j++)
602  {
603  Ratio ratioups(ups[i][j], upsNoOsc[i][j]);
604  upsNoOscRatio[i].push_back(ratioups.ToTH1());
605  Ratio ratiodns(dns[i][j], dnsNoOsc[i][j]);
606  dnsNoOscRatio[i].push_back(ratiodns.ToTH1());
607  }
608  }
609 
610  for(unsigned int i = 0; i < hNoOsc.size(); i++)
611  {
612  TH1* tempdenom = (TH1*)hNoOsc[i]->Clone();
613  TH1* tempratio = (TH1*)hPred[i]->Clone();
614  hNoOscDenom.push_back(tempdenom);
615  hPredRatio .push_back(tempratio);
616  hNoOscDenom[i]->Add(hBkg[i], -1);
617  hPredRatio [i]->Add(hBkg[i], -1);
618  hPredRatio [i]->Divide(hNoOscDenom[i]);
619 
620  hPredRatio[i]->GetXaxis()->CenterTitle();
621  hPredRatio[i]->GetYaxis()->CenterTitle();
622  hPredRatio[i]->GetXaxis()->SetTitle("Reconstructed Neutrino Energy (GeV)");
623  hPredRatio[i]->GetYaxis()->SetTitle("Ratio To No Oscillation");
624 
625  TFeldmanCousins fc(0.6827);
626  for(int j = 0; j < vec_DataRatio[i]->GetN(); j++)
627  {
628  double x, y, yUnscaled;
629  vec_DataRatio[i]->GetPoint(j, x, y);
630 
631  int bin = hCosmics[i]->FindBin(x);
632  yUnscaled = hDataUnscaled[i]->GetBinContent(bin);
633  double scaleFactor = 0.1*hCosmics[i]->GetBinWidth(bin);
634  y -= hCosmics[i]->GetBinContent(bin) - hBkg[i]->GetBinContent(bin);
635  vec_DataRatio[i]->SetPoint (j, x, y/hNoOscDenom[i]->GetBinContent(bin));
636  if(yUnscaled<0)
637  {
638  vec_DataRatio[i]->SetPointEYlow (j, 0.);
639  }
640  else
641  {
642  vec_DataRatio[i]->SetPointEYlow (j,(yUnscaled-fc.CalculateLowerLimit(yUnscaled, 0))/hNoOscDenom[i]->GetBinContent(bin));
643  }
644  vec_DataRatio[i]->SetPointEYhigh(j,(fc.CalculateUpperLimit(yUnscaled, 0)-yUnscaled)/hNoOscDenom[i]->GetBinContent(bin));
645  /*
646  std::cout << i << " " << bin << " " << y << " "
647  << (fc.CalculateLowerLimit(y, 0)) << " " << fc.CalculateUpperLimit(y, 0) << " " << hNoOscDenom[i]->GetBinContent(bin) << " "
648  << (y-fc.CalculateLowerLimit(y, 0))/hNoOscDenom[i]->GetBinContent(bin) << " " << vec_DataRatio[i]->GetErrorYlow(j) << " "
649  << (fc.CalculateUpperLimit(y, 0) -y)/hNoOscDenom[i]->GetBinContent(bin) << std::endl;
650  */
651  }
652  }
653 
654  TLegend *leg_Rat = new TLegend(0.22, 0.65, 0.57, 0.8);
655  //TLegend *leg_Rat = new TLegend(0.26, 0.7, 0.6, 0.85);
656  leg_Rat->AddEntry(vec_Data [0], "FD Data", "LEP");
657  leg_Rat->AddEntry(hPredRatio[0], "Prediction", "L" );
658  TLegendEntry *ent_Band_Rat;
659  if(drawSystBand)
660  {
661  ent_Band_Rat = leg_Rat->AddEntry("errorRat","1-#sigma syst. range", "BF");
662  ent_Band_Rat->SetFillStyle(1001);
663  ent_Band_Rat->SetFillColor(color_band);
664  ent_Band_Rat->SetLineColor(color_band);
665  }
666  leg_Rat->SetTextSize(0.05);
667  leg_Rat->SetBorderSize(0);
668  leg_Rat->SetFillStyle(0);
669 
670  c_Rat->cd();
671  PlotWithSystErrorBand(hPredRatio[0], upsNoOscRatio[0], dnsNoOscRatio[0], color_pred, color_band, 1.6);
672  TLine *l_At1 = new TLine(hPredRatio[0]->GetXaxis()->GetXmin(), 1., hPredRatio[0]->GetXaxis()->GetXmax(), 1.);
673  TLine *l_At0 = new TLine(hPredRatio[0]->GetXaxis()->GetXmin(), 0., hPredRatio[0]->GetXaxis()->GetXmax(), 0.);
674  l_At1->SetLineColor(kGray);
675  l_At1->SetLineStyle(kDashed);
676  l_At1->Draw();
677  l_At0->SetLineColor(kGray);
678  l_At0->SetLineStyle(kDashed);
679  l_At0->Draw();
680  vec_DataRatio[0]->Draw("P0");
681  hPredRatio [0]->GetYaxis()->SetRangeUser(-0.1, 1.6);
682 
684  MyCornerLabel(s_HornLabel);
685  //drawQuantLabel(0);
686  gPad->RedrawAxis();
687 
688  leg_Rat->Draw();
689  c_Rat->SaveAs((TString)s_OutDir + c_Rat->GetName() + (TString)s_FileAppend + ".pdf");
690  c_Rat->Write();
691 
692  c_ByQ_Rat->cd();
693 
694  for(unsigned int i = 1; i < hBkg.size(); i++)
695  {
696  if(i==1) p1_Rat->cd();
697  if(i==2) p2_Rat->cd();
698  if(i==3) p3_Rat->cd();
699  if(i==4) p4_Rat->cd();
700  //MakeHistCanvasReady_Quant(i, hPredRatio[i], 1.4);
701  PlotWithSystErrorBand_Quant(i, hPredRatio[i], upsNoOscRatio[i], dnsNoOscRatio[i], color_pred, color_band, 1.6);
702  //hPredRatio [i]->Draw("HIST");
703  vec_DataRatio[i]->Draw("P0");
704  if(i==3)
705  {
706  hPredRatio[i]->GetXaxis()->ChangeLabel(6,-1,0);
707  }
708  hPredRatio[i]->GetYaxis()->SetRangeUser(-0.15, 1.6);
709  leg_Rat->Draw();
710  drawQuantLabel(i);
711  }
712 
713  c_ByQ_Rat->cd();
715  MyCornerLabel(s_HornLabel);
716  xTitle->Draw();
717  yTitle->Draw();
718 
719  c_ByQ_Rat->SetName("c_dataAndPrediction_ByQ_RatioNoOsc");
720  c_ByQ_Rat->SaveAs((TString)s_OutDir + c_ByQ_Rat->GetName() + (TString)s_FileAppend + ".pdf");
721  c_ByQ_Rat->Write();
722 
723 
724  /*
725  TCanvas *c_RatSpec = new TCanvas("c_RatSpec_", "c_RatSpec_", 1000, 1000);
726  TPad *p1_RatSpec = new TPad("p1_RatSpec", "p1_RatSpec", 0, 0.3, 1, 1);
727  TPad *p2_RatSpec = new TPad("p2_RatSpec", "p2_RatSpec", 0, 0, 1, 0.3);
728  p1_RatSpec->SetBottomMargin(0.);
729  p2_RatSpec->SetTopMargin(0.);
730 
731  p1_RatSpec->cd();
732  if(drawSystBand)
733  {
734  TH1 *nomNoOsc = spec_NoOsc[0].ToTH1(pot);
735  nomNoOsc->GetYaxis()->SetTitle("Events / 0.1 GeV");
736  std::vector<TH1*> upshNoOsc, dnshNoOsc;
737  for(const auto& upShift: upsNoOsc[0]) upshNoOsc.push_back(upShift .ToTH1(pot));
738  for(const auto& downShift: dnsNoOsc[0]) dnshNoOsc.push_back(downShift.ToTH1(pot));
739  PlotWithSystErrorBand(nomNoOsc, upshNoOsc, dnshNoOsc, color_noosc, color_noosc_b, 1.3);
740  TH1 *nom = spec_Pred[0].ToTH1(pot);
741  std::vector<TH1*> upsh, dnsh;
742  for(const auto& upShift: ups[0]) upsh.push_back(upShift.ToTH1(pot));
743  for(const auto& downShift: dns[0]) dnsh.push_back(downShift.ToTH1(pot));
744  PlotWithSystErrorBand(nom, upsh, dnsh, color_pred, color_band, 1., false);
745  }
746  else
747  {
748  hNoOsc[0]->Draw("HIST");
749  hNoOsc[0]->GetYaxis()->SetRangeUser(0, hNoOsc[0]->GetMaximum()*1.15);
750  hNoOsc[0]->GetYaxis()->SetTitle("Events / 0.1 GeV");
751  }
752  hAllBkg [0]->Draw("HIST SAME");
753  hPred [0]->Draw("HIST SAME");
754  vec_Data[0]->Draw("P");
755 
756  //drawQuantLabel(0);
757  //leg_NoOsc->Draw();
758 
759  p2_RatSpec->cd();
760  // PlotWithSystErrorBand(hPredRatio[0], upsNoOscRatio[0], dnsNoOscRatio[0], color_pred, color_band, 1.6);
761  vec_DataRatio[0]->Draw("P0");
762  l_At1->Draw();
763  hPredRatio[0]->GetYaxis()->SetRangeUser(-0.1, 1.6);
764 
765  leg_Rat->Draw();
766  gPad->RedrawAxis();
767 
768 
769  //c_RatSpec->cd();
770  //PreliminaryBoxOpening();
771  //MyCornerLabel(s_HornLabel);
772 
773  c_RatSpec->cd();
774  p1_RatSpec->Draw();
775  p2_RatSpec->Draw();
776  c_RatSpec->SetName("c_dataAndPrediction_RatioNoOsc_RatSpec");
777  c_RatSpec->SaveAs((TString)s_OutDir + c_RatSpec->GetName() + (TString)s_FileAppend + ".pdf");
778  c_RatSpec->Write();
779  */
780 
781  }
782 
783  f_Out->Close();
784 
785 }
786 
T max(const caf::Proxy< T > &a, T b)
Color_t color_band
Pass neutrinos through unchanged.
Definition: IOscCalc.h:45
enum BeamMode kRed
virtual _IOscCalcAdjustable< T > * Copy() const =0
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
Antineutrinos-only.
Definition: IPrediction.h:50
TH1D * ToTH1(Color_t col=kBlack, Style_t style=kSolid) const
Definition: Ratio.cxx:67
std::vector< SystGroupDef > systs
Definition: syst_header.h:385
void MakeHistCanvasReady_Quant(const int quant, TH1 *hist, double ymax)
Definition: Plots.cxx:1365
double th23
Definition: runWimpSim.h:98
const double kAna2019RHCLivetime
Definition: Exposures.h:227
Loads shifted spectra from files.
Style_t line_pred
Simple record of shifts applied to systematic parameters.
Definition: SystShifts.h:20
TGraphAsymmErrors * PlotWithSystErrorBand(IPrediction *pred, const std::vector< const ISyst * > &systs, osc::IOscCalc *calc, double pot, int col, int errCol, float headroom, bool newaxis, EBinType bintype, double alpha)
Plot prediction with +/-1sigma error band.
Definition: Plots.cxx:304
Style_t marker_data
Style_t marker_bkg
Style_t marker_cos
void CenterTitles(TH1 *histo)
Definition: Plots.cxx:1483
static SystShifts Nominal()
Definition: SystShifts.h:34
osc::OscCalcDumb calc
osc::IOscCalcAdjustable * DefaultOscCalc()
Create a new calculator with default assumptions for all parameters.
Definition: Calcs.cxx:49
Encapsulate code to systematically shift a caf::SRProxy.
Definition: ISyst.h:14
virtual void SetDmsq32(const T &dmsq32)=0
Representation of a spectrum in any variable, with associated POT.
Definition: Spectrum.h:39
Style_t line_ws
Style_t marker_ws
void Plotting_DataAndPrediction(std::string s_HornCurrent="fhc", bool drawSystBand=true, bool makeUnoscPlots=true, std::string s_FileAppend="", double dmsq=0., double th23=0., bool isPreliminary=true)
std::vector< const ISyst * > getAllAna2018Systs(const EAnaType2018 ana, const bool smallgenie, const BeamType2018 beam, bool isFit)
Charged-current interactions.
Definition: IPrediction.h:39
Color_t color_ws
Color_t color_noosc
static std::unique_ptr< Spectrum > LoadFrom(TDirectory *dir, const std::string &name)
Definition: Spectrum.cxx:563
Style_t marker_nc
correl_xv GetXaxis() -> SetDecimals()
virtual _IOscCalc< T > * Copy() const override
Definition: IOscCalc.h:49
#define pot
const double j
Definition: BetheBloch.cxx:29
const double kAna2019FHCLivetime
Definition: Exposures.h:226
void PreliminaryBoxOpening()
Definition: numu_tools.h:187
float bin[41]
Definition: plottest35.C:14
Color_t color_cos
void MyCornerLabel(std::string Str)
Represent the ratio between two spectra.
Definition: Ratio.h:8
Neutrinos-only.
Definition: IPrediction.h:49
std::vector< double > POT
const double kAna2019RHCPOT
Definition: Exposures.h:224
Color_t color_bkg
double livetime
Definition: saveFDMCHists.C:21
TGraphAsymmErrors * PlotWithSystErrorBand_Quant(const int quant, IPrediction *pred, const std::vector< const ISyst * > &systs, osc::IOscCalc *calc, double pot, int col, int errCol, float maxy, bool newaxis)
Definition: Plots.cxx:543
enum BeamMode kViolet
Color_t color_noosc_b
Style_t line_bkg
TPaveText * drawQuantLabel(int quantId=0)
Definition: numu_tools.h:4
virtual void SetTh23(const T &th23)=0
Neutral-current interactions.
Definition: IPrediction.h:40
Style_t line_nc
Color_t color_data
const double kAna2019FHCPOT
Definition: Exposures.h:223
Both neutrinos and antineutrinos.
Definition: IPrediction.h:52
Color_t color_pred
void SplitCanvasQuant(TCanvas *&canvas, TPad *&pad1, TPad *&pad2, TPad *&pad3, TPad *&pad4)
Definition: Plots.cxx:1439
All neutrinos, any flavor.
Definition: IPrediction.h:26
Color_t color_nc
enum BeamMode kGreen
void DecorateHist(TH1 *h, Style_t lineStyle, Color_t lineColor, int lineWidth, Style_t markerStyle, Color_t markerColor, double markerSize)
Definition: PlottingTools.h:94
enum BeamMode kBlue
void SetShift(const ISyst *syst, double shift, bool force=false)
Definition: SystShifts.cxx:81
void FillWithDimColor(TH1 *h, bool usealpha, float dim)
Style_t line_data
Style_t line_cos
Style_t marker_pred
enum BeamMode string