plot_predictions.C
Go to the documentation of this file.
1 #pragma once
2 
7 
8 using namespace ana;
9 
11 
12 bool DianaIsLazy = true;
13 bool errorBand = true;
14 
15 double maxy = 0;
16 double maxy_one = 1.999;
17 double maxy_two = 3.999;
18 double paperscale = 0.8/0.5; //proportions in the latex document, bc reasons
19 
20 double maxy_pad0 = 0;
21 double maxy_rhc = 5;
22 double maxy_fhc = 12;
23 
25  std::string anaName,
26  bool nh)
27 {
28 
29  double sinHere = 0.558;
30  double deltaHere = 0.00244;
32  calc->SetTh23(asin(sqrt(sinHere)));
33  calc->SetDmsq32(deltaHere);
34 
35 }
36 
37 
38 
39 void HornSettings(int quantile, double &pot, double &livetime, Sign::Sign_t &wrongs, double &maxy){
40  if(quantile<4){
41  pot = pot_one;
42  livetime = livetime_one;
43  wrongs = ws_one;
44  maxy = maxy_one;
45  }
46  else{
47  pot = pot_two;
48  livetime = livetime_two;
49  wrongs = ws_two;
50  maxy = maxy_two;
51  }
52 
53 }// end HornSettings
54 
55 
56 
57 void plot_predictions(std::string anaName = "fhc",
58  std::string calculator = "2018calc",
59  bool nh = true,
60  bool fake = true)
61 {
62 
63 
64  bool preliminary = false;
65  std::string anaName2 ="";
66  if(anaName=="rhc"){
67  anaName2="Antineutrino beam";
68  totquant = 4;
69  twobeams = false;
71  }
72  // if fhc, change settings so the first 4 quantiles are fhc instead of rhc
73  // dont do anything if rhc as the first 4 quantiles are rhc already
74  if(anaName=="fhc"){
75  anaName2="Neutrino beam";
76  horn_one = "fhc";
77  period_one = "full";
78  pot_one = pot_fhc;
81  maxy_one = 5.999;
82  twobeams = false;
83  totquant = 4;
85  }
86 
87  // calculator
88  std::string calcName = "2018calc";
89  std::string hierarchy = "";
91  RestartCalculator(calc, anaName, nh);
93  auto calcNoOsc = osc.Copy();
94 
95 
96  // print on screen /////////////////////////////////////////////////////////
97  std::cout << "\n\nmake table with predicted number of events\n" << std::endl;
98  std::cout << "ana name: " << anaName << std::endl;
99  std::cout << "calculator: " << calcName << ", hierarchy: " << hierarchy << std::endl;
100 
101 
102  //inputs
103  TString dout_name = "";
104  std::string ddata_name = "/nova/ana/nu_mu_ana/Ana2018/Data/";
105  std::string dcosm_name = "/nova/ana/nu_mu_ana/Ana2018/Cosmics/";
106  std::string dpred_name = "/nova/ana/nu_mu_ana/Ana2018/Predictions/";
107  std::string fpred_name = "";
108  std::string fpred_suffix = "";
109 
110  TFile* fdata_one = new TFile((ddata_name + "fd_dataspectrum_" + horn_one + "_full__numu2018.root").c_str());
111  TFile* fdata_two = new TFile((ddata_name + "fd_dataspectrum_" + horn_two + "_full__numu2018.root").c_str());
112  std::vector<Spectrum*> s_realdata;
113  std::cout << "\nloading data files and spectra" << std::endl;
114  for(int quant=1; quant<totquant+1; quant++){
115  if(quant<5) s_realdata.push_back(LoadFrom<Spectrum>(fdata_one, Form("fd_data_q%d", quant) ).release());
116  else s_realdata.push_back(LoadFrom<Spectrum>(fdata_two, Form("fd_data_q%d", quant-4) ).release());
117  }
118 
119 
120  double pot1here = s_realdata[0]->POT();
121  std::cout << "\n\n pot for this sample = " << pot1here << "\n\n" << std::endl;
122 
123  // cosmics
124  std::string cosmics = "nocosmics";
125  std::vector<Spectrum> s_cosmics;
126  std::vector<TH1*> hCosmics;
127  std::cout << "\nloading cosmics" << std::endl;
128  std::cout << "" << horn_one << std::endl;
129  std::string fcosmics_name_one = dcosm_name + "cosmics_" + horn_one + "__numu2018.root";
130  TFile fcosmics_one(fcosmics_name_one.c_str());
131  hCosmics.push_back((TH1*)fcosmics_one.Get("cosmics_q0")); // remember, quant0 = all
132  hCosmics.push_back((TH1*)fcosmics_one.Get("cosmics_q1"));
133  hCosmics.push_back((TH1*)fcosmics_one.Get("cosmics_q2"));
134  hCosmics.push_back((TH1*)fcosmics_one.Get("cosmics_q3"));
135  hCosmics.push_back((TH1*)fcosmics_one.Get("cosmics_q4"));
136  auto *cosmics_quant1_one = (TH1*)fcosmics_one.Get("cosmics_q1");
137  auto *cosmics_quant2_one = (TH1*)fcosmics_one.Get("cosmics_q2");
138  auto *cosmics_quant3_one = (TH1*)fcosmics_one.Get("cosmics_q3");
139  auto *cosmics_quant4_one = (TH1*)fcosmics_one.Get("cosmics_q4");
140  s_cosmics.push_back(Spectrum (cosmics_quant1_one, pot_one, livetime_one));
141  s_cosmics.push_back(Spectrum (cosmics_quant2_one, pot_one, livetime_one));
142  s_cosmics.push_back(Spectrum (cosmics_quant3_one, pot_one, livetime_one));
143  s_cosmics.push_back(Spectrum (cosmics_quant4_one, pot_one, livetime_one));
144  if(twobeams){
145  std::string fcosmics_name_two = dcosm_name + "cosmics_" + horn_two + "__numu2018.root";
146  TFile fcosmics_two(fcosmics_name_two.c_str());
147  hCosmics.push_back((TH1*)fcosmics_two.Get("cosmics_q1"));
148  hCosmics.push_back((TH1*)fcosmics_two.Get("cosmics_q2"));
149  hCosmics.push_back((TH1*)fcosmics_two.Get("cosmics_q3"));
150  hCosmics.push_back((TH1*)fcosmics_two.Get("cosmics_q4"));
151  auto *cosmics_quant1_two = (TH1*)fcosmics_two.Get("cosmics_q1");
152  auto *cosmics_quant2_two = (TH1*)fcosmics_two.Get("cosmics_q2");
153  auto *cosmics_quant3_two = (TH1*)fcosmics_two.Get("cosmics_q3");
154  auto *cosmics_quant4_two = (TH1*)fcosmics_two.Get("cosmics_q4");
155  s_cosmics.push_back(Spectrum (cosmics_quant1_two, pot_two, livetime_two));
156  s_cosmics.push_back(Spectrum (cosmics_quant2_two, pot_two, livetime_two));
157  s_cosmics.push_back(Spectrum (cosmics_quant3_two, pot_two, livetime_two));
158  s_cosmics.push_back(Spectrum (cosmics_quant4_two, pot_two, livetime_two));
159  }
160  if(twobeams){
161  for(int quantId = 5; quantId < totquant+1; quantId++){
162  hCosmics[0]->Add(hCosmics[quantId]);
163  }
164  }
165 
166 
167  std::cout << "\nloading predictions" << std::endl;
168  const std::string dpred_name_fhc= dpred_name+"pred_numuconcat_fhc__numu2018.root";
169  const std::string dpred_name_rhc= dpred_name+"pred_numuconcat_rhc__numu2018.root";
170  std::string dpred_name_one = dpred_name_rhc;
171  std::string dpred_name_two = dpred_name_fhc;
172  if(anaName=="fhc") dpred_name_one = dpred_name_fhc;
173 
174  std::vector< const ISyst* > systs;
175  systs = getAllAna2018Systs(kNumuAna2018,true);
176  int systSize = systs.size();
177 
178  std::cout << "\nloading predictions" << std::endl;
179  std::vector<PredictionSystJoint2018*> predictions;
180  ENu2018ExtrapType extrap = kNuMu;
181  for(int quant = 1; quant < totquant+1; quant++){
182  std::cout << "quantile " << quant << std::endl;
183  if(quant<5) predictions.push_back(new PredictionSystJoint2018(extrap, calc, horn_one, quant, systs, dpred_name_one));
184  else predictions.push_back(new PredictionSystJoint2018(extrap, calc, horn_two, quant-4, systs, dpred_name_two));
185  }
186 
187 
188  std::vector<TH1*> hDataUnscaled;
189  std::vector<TH1*> hData;
190  std::vector<TH1*> hNC;
191  std::vector<TH1*> hCC;
192  std::vector<TH1*> hNuMu;
193  std::vector<TH1*> hNoOsc;
194  std::vector<TH1*> hPred;
195  std::vector<TH1*> hBkg;
196  std::vector<TH1*> hAllBkg;
197  std::vector<TH1*> hWS;
198 
199 
201  hData.push_back(s_realdata[0]->ToTH1(s_realdata[0]->POT()));
202  hDataUnscaled.push_back(s_realdata[0]->ToTH1(s_realdata[0]->POT()));
203  hNC.push_back(predictions[0]->PredictComponent(calc, Flavors::kAll, Current::kNC, Sign::kBoth).ToTH1(pot));
204  hCC.push_back(predictions[0]->PredictComponent(calc, Flavors::kAll, Current::kCC, Sign::kBoth).ToTH1(pot));
205  hNuMu.push_back(predictions[0]->PredictComponent(calc, Flavors::kAllNuMu, Current::kCC, Sign::kBoth).ToTH1(pot));
206  hBkg.push_back(predictions[0]->Predict(calc).ToTH1(pot));
207  hAllBkg.push_back(predictions[0]->Predict(calc).ToTH1(pot));
208  hPred.push_back(predictions[0]->Predict(calc).ToTH1(pot));
209  hNoOsc.push_back(predictions[0]->Predict(calcNoOsc).ToTH1(pot));
210  hWS.push_back(predictions[0]->PredictComponent(calc, Flavors::kAllNuMu, Current::kCC, wrongs).ToTH1(pot));
211 
212  for(int quantId = 0; quantId < totquant; quantId++){
213  HornSettings(quantId, pot, livetime, wrongs, maxy);
214  hData.push_back(s_realdata[quantId]->ToTH1(s_realdata[0]->POT()));
215  hDataUnscaled.push_back(s_realdata[quantId]->ToTH1(s_realdata[0]->POT()));
216  hNC.push_back(predictions[quantId]->PredictComponent(calc, Flavors::kAll, Current::kNC, Sign::kBoth).ToTH1(pot));
217  hCC.push_back(predictions[quantId]->PredictComponent(calc, Flavors::kAll, Current::kCC, Sign::kBoth).ToTH1(pot));
218  hNuMu.push_back(predictions[quantId]->PredictComponent(calc, Flavors::kAllNuMu, Current::kCC, Sign::kBoth).ToTH1(pot));
219  hBkg.push_back(predictions[quantId]->Predict(calc).ToTH1(pot));
220  hAllBkg.push_back(predictions[quantId]->Predict(calc).ToTH1(pot));
221  hPred.push_back(predictions[quantId]->Predict(calc).ToTH1(pot));
222  hNoOsc.push_back(predictions[quantId]->Predict(calcNoOsc).ToTH1(pot));
223  hWS.push_back(predictions[quantId]->PredictComponent(calc, Flavors::kAllNuMu, Current::kCC, wrongs).ToTH1(pot));
224  if(quantId > 0){ // already added the first quantile
225  hDataUnscaled[0] ->Add(s_realdata[quantId]->ToTH1(s_realdata[0]->POT()));
226  hData[0] ->Add(s_realdata[quantId]->ToTH1(s_realdata[0]->POT()));
227  hNC[0] ->Add(predictions[quantId]->PredictComponent(calc, Flavors::kAll, Current::kNC, Sign::kBoth).ToTH1(pot));
228  hCC[0] ->Add(predictions[quantId]->PredictComponent(calc, Flavors::kAll, Current::kCC, Sign::kBoth).ToTH1(pot));
229  hNuMu[0] ->Add(predictions[quantId]->PredictComponent(calc, Flavors::kAllNuMu, Current::kCC, Sign::kBoth).ToTH1(pot));
230  hBkg[0] ->Add(predictions[quantId]->Predict(calc).ToTH1(pot));
231  hAllBkg[0]->Add(predictions[quantId]->Predict(calc).ToTH1(pot));
232  hPred[0] ->Add(predictions[quantId]->Predict(calc).ToTH1(pot));
233  hNoOsc[0] ->Add(predictions[quantId]->Predict(calcNoOsc).ToTH1(pot));
234  hWS[0] ->Add(predictions[quantId]->PredictComponent(calc, Flavors::kAllNuMu, Current::kCC, wrongs).ToTH1(pot));
235  }
236  }
237  // fake stack: draw bkg, cosmics, ws in this order
238  // add bkg and cosmics to ws
239  // add cosmics to bkg
240  for(int quantId = 0; quantId < totquant+1; quantId++){
241  hBkg[quantId]->Add(hNuMu[quantId], -1);
242  hAllBkg[quantId]->Add(hNuMu[quantId], -1);
243  hAllBkg[quantId]->Add(hCosmics[quantId]);
244  hWS[quantId]->Add(hAllBkg[quantId]);
245  hPred[quantId]->Add(hCosmics[quantId]);
246  }
247 
248 
249  std::cout << "\n\nIntegrated events = " << hData[0]->Integral() << "\n\n" << std::endl;
250 
251  // scale everything, everything!
252  std::cout << "\nScaling histograms to 0.1 GeV" << std::endl;
253  for(int quantId = 0; quantId < totquant+1; quantId++){
254  hData[quantId] ->Scale(0.1, "width");
255  hCosmics[quantId]->Scale(0.1, "width");
256  hNC[quantId] ->Scale(0.1, "width");
257  hCC[quantId] ->Scale(0.1, "width");
258  hBkg[quantId] ->Scale(0.1, "width");
259  hAllBkg[quantId] ->Scale(0.1, "width");
260  hNuMu[quantId] ->Scale(0.1, "width");
261  hPred[quantId] ->Scale(0.1, "width");
262  hNoOsc[quantId] ->Scale(0.1, "width");
263  hWS[quantId] ->Scale(0.1, "width");
264  }
265 
266  // data graphs
267  std::cout << "\nDefine data graphs\n" << std::endl;
268  std::vector<TGraph*> grData;
269  for(int quantId = 0; quantId < totquant+1; quantId++){
270  grData.push_back(graphAsymmErrorScaled(hData[quantId], hDataUnscaled[quantId]));
271  }
272 
273 
274 
275  // predictions again?
276  std::vector<Spectrum> sPrediction;
277  for(int quantId = 0; quantId < totquant+1; quantId++){
278  if(quantId==0){
280  Spectrum prediction(hPred[quantId], pot, livetime);
281  sPrediction.push_back(prediction);
282  }
283  else{
284  HornSettings(quantId-1, pot, livetime, wrongs, maxy);
285  Spectrum prediction(hPred[quantId], pot, livetime);
286  sPrediction.push_back(prediction);
287  }
288  }
289 
290 
291  // save up and dn shifts for each quantile
292  // looks messy but really makes it less difficult
293  std::vector<Spectrum> ups[totquant+1], dns[totquant+1];
294  std::cout << "\nFilling ups and dns\n" << std::endl;
295 
296  if(errorBand){
297  for(int quantId = 0; quantId < totquant+1; quantId++){
298 
299  std::cout << "quantile " << quantId << std::endl;
300  // spectrums for syst band best fit
301  for(const ISyst* syst: systs){
302  SystShifts shifts;
303 
304  if(quantId==0){
306  shifts.SetShift(syst, +1);
307  TH1 * htempup = predictions[0]->PredictSyst(calc, shifts).ToTH1(pot);
308  for(int i = 1; i < totquant; i++){htempup->Add(predictions[i]->PredictSyst(calc, shifts).ToTH1(pot));}
309  htempup->Scale(0.1, "width");
310  htempup->Add(hCosmics[0]);
311  Spectrum specup(htempup, pot, livetime);
312  ups[quantId].push_back(specup);
313 
314  shifts.SetShift(syst, -1);
315  TH1 * htempdn = predictions[0]->PredictSyst(calc, shifts).ToTH1(pot);
316  for(int i = 1; i < totquant; i++){htempdn->Add(predictions[i]->PredictSyst(calc, shifts).ToTH1(pot));}
317  htempdn->Scale(0.1, "width");
318  htempdn->Add(hCosmics[0]);
319  Spectrum specdn(htempdn, pot, livetime);
320  dns[quantId].push_back(specdn);
321  }// all quantiles
322 
323  else{
324  HornSettings(quantId-1, pot, livetime, wrongs, maxy);
325  shifts.SetShift(syst, +1);
326  TH1* htempup = predictions[quantId-1]->PredictSyst(calc, shifts).ToTH1(pot);
327  htempup->Scale(0.1, "width");
328  htempup->Add(hCosmics[quantId]);
329  Spectrum specup(htempup, pot, livetime);
330  ups[quantId].push_back(specup);
331  shifts.SetShift(syst, -1);
332  TH1 * htempdn = predictions[quantId-1]->PredictSyst(calc, shifts).ToTH1(pot);
333  htempdn->Scale(0.1, "width");
334  htempdn->Add(hCosmics[quantId]);
335  Spectrum specdn(htempdn, pot, livetime);
336  dns[quantId].push_back(specdn);
337  }
338  } // shifts
339  }// for quantiles
340  }
341 
342 
343 
344  // pimping
345  for(int quantId = 0; quantId< 1; quantId++){
346  PimpHist(hData[quantId], line_data, color_data, 3, marker_data, color_data, 1);
347  PimpHist(hWS[quantId], line_ws, color_ws, 1, marker_ws, color_ws, 1);
348  PimpHist(hNC[quantId], line_nc, color_nc, 1, marker_nc, color_nc, 1);
349  PimpHist(hPred[quantId], line_pred, color_pred, 3, marker_pred, color_pred, 1);
350  PimpHist(hCosmics[quantId], line_cos, color_cos, 1, marker_cos, color_cos, 1);
351  PimpHist(hAllBkg[quantId], line_bkg, color_bkg, 1, marker_bkg, color_bkg, 1);
352  PimpHist(hBkg[quantId], line_bkg, color_bkg, 1, marker_bkg, color_bkg, 1);
353  FillWithDimColor(hWS[quantId], 0);
354  FillWithDimColor(hCosmics[quantId], 0);
355  FillWithDimColor(hAllBkg[quantId], 0);
356  grData[quantId]->SetMarkerSize(1.);
357  }// for quantId
358  //root is not behaving for the quantiles
359  //set different line width and marker size
360  for(int quantId = 1; quantId< totquant+1; quantId++){
361  PimpHist(hData[quantId], line_data, color_data, 2, marker_data, color_data, 3);
362  PimpHist(hWS[quantId], line_ws, color_ws, 2, marker_ws, color_ws, 3);
363  PimpHist(hNC[quantId], line_nc, color_nc, 2, marker_nc, color_nc, 3);
364  PimpHist(hPred[quantId], line_pred, color_pred, 2, marker_pred, color_pred, 3);
365  PimpHist(hCosmics[quantId], line_cos, color_cos, 2, marker_cos, color_cos, 3);
366  PimpHist(hAllBkg[quantId], line_bkg, color_bkg, 2, marker_bkg, color_bkg, 3);
367  PimpHist(hBkg[quantId], line_bkg, color_bkg, 2, marker_bkg, color_bkg, 3);
368  FillWithDimColor(hWS[quantId], 0);
369  FillWithDimColor(hCosmics[quantId], 0);
370  FillWithDimColor(hAllBkg[quantId], 0);
371  grData[quantId]->SetMarkerSize(3.);
372  }// for quantId
373 
374 
375  // for later use in the 4pad plots
376  TLatex* xTitle = new TLatex(0.5, 0.03, "Reconstructed Energy (GeV)");
377  xTitle->SetNDC();
378  xTitle->SetTextAlign(22);
379  xTitle->SetTextSize(0.055 /paperscale);
380  TLatex* yTitle = new TLatex(0.02, 0.5, "Events / 0.1 GeV");
381  yTitle->SetNDC();
382  yTitle->SetTextAlign(22);
383  yTitle->SetTextAngle(90.);
384  yTitle->SetTextSize(0.055 / paperscale);
385 
386  const int cx = 960, cy = 800;
387  // all quantiles together
388  // equals quant 0
389  if(DianaIsLazy){
390  std::ofstream file;
391  TString outName = dout_name + "fd_prediction_spectra__" + anaName + "__2017calc";
392  std::string idName = anaName2;
393  int quantId = 0;
394  TCanvas *canvas = new TCanvas("canvas","canvas", cx, cy);
395  canvas->cd();
396  ana::CenterTitles(hPred[0]);
397  hPred[0]->GetYaxis()->SetTitle("Events / 0.1 GeV");
398  hPred[0]->GetYaxis()->SetRangeUser(0.0, maxy_pad0);
399  hPred[0]->GetXaxis()->SetLabelSize(hPred[0]->GetXaxis()->GetLabelSize());
400  hPred[0]->GetYaxis()->SetLabelSize(hPred[0]->GetYaxis()->GetLabelSize());
401  hPred[0]->Draw("hist");
403  if(errorBand) PlotWithSystErrorBand_Quant(quantId, sPrediction[quantId], ups[quantId], dns[quantId], pot, color_pred, color_band, maxy_pad0, false);
404  else hPred[quantId]->Draw("hist same");
405  hWS[quantId]->Draw("hist same");
406  hAllBkg[quantId]->Draw("hist same");
407  hCosmics[quantId]->Draw("hist same");
408  hPred[quantId]->Draw("hist same");
409  TLegend *legend = new TLegend(0.55,0.55,0.88,0.88);
410  legend->AddEntry(hPred[quantId], "2017 Best Fit Pred.","l");
411  TLegendEntry *entry=legend->AddEntry("error","1-#sigma syst. range","bf");
412  entry->SetFillStyle(1001);
413  entry->SetFillColor(color_band);
414  entry->SetLineColor(color_band);
415  if(anaName=="fhc")legend->AddEntry(hWS[quantId], "Wrong Sign:#bar{#nu}_{#mu}CC","bf");
416  else legend->AddEntry(hWS[quantId], "Wrong Sign:#nu_{#mu}CC","bf");
417  legend->AddEntry(hAllBkg[quantId], "Total bkg.","bf");
418  legend->AddEntry(hCosmics[quantId], "Cosmic bkg.","bf");
419  legend->SetTextSize(0.04); //no border for legend
420  legend->SetBorderSize(0); //no border for legend
421  legend->SetFillStyle(0); //fill colour is white
422  legend->Draw();
423  CornerLabel(idName);
425  drawQuantLabel(0);
426  gPad->RedrawAxis();
427  gPad->Update();
428  canvas->cd();
429 
430  canvas->SaveAs(outName+".eps");
431  canvas->SaveAs(outName+".png");
432  canvas->SaveAs(outName+".pdf");
433  file.open (outName+".txt");
434  file<< "FD spectra showing the prediction at the 2017 best fit with systematics error band, and beam and cosmic backgrounds for the NuMI beam in " << anaName << " mode.";
435  file.close();
436 
437  } // plot all quantiles
438 
439 
440  if(DianaIsLazy){
441 
442  std::ofstream file;
443  TCanvas *canvas;
444  TPad *pad1, *pad2, *pad3, *pad4;
445  SplitCanvasQuant(canvas, pad1, pad2, pad3, pad4);
446  gStyle->SetLineScalePS(3./paperscale);
447  canvas->SetCanvasSize(paperscale * canvas->GetWw(), paperscale * canvas->GetWh());
448  TString outName = dout_name + "fd_prediction_spectra__" + anaName+ "__2017calc";
449  std::string idName = anaName2;
450  for(int quantId:{3,4,1,2}){
451  if(quantId==1) pad1->cd();
452  if(quantId==2) pad2->cd();
453  if(quantId==3) pad3->cd();
454  if(quantId==4) pad4->cd();
455  HornSettings(quantId-1, pot, livetime, wrongs, maxy);
456  MakeHistCanvasReady_Quant(quantId, hPred[quantId], maxy);
457  // hack to avoid 0-5 overlapping in the x-axis
458  hPred[quantId]->GetXaxis()->SetLabelSize(hPred[quantId]->GetXaxis()->GetLabelSize()/paperscale);
459  hPred[quantId]->GetYaxis()->SetLabelSize(hPred[quantId]->GetYaxis()->GetLabelSize()/paperscale);
460  hPred[quantId]->Draw("hist");
461  if(quantId==3){
462  TGaxis* gax1 = new TGaxis(0., 0., 4., 0., 0.,4.,4, "");
463  gax1->SetLabelSize(hPred[quantId]->GetXaxis()->GetLabelSize());
464  gax1->SetLabelOffset(hPred[quantId]->GetXaxis()->GetLabelOffset());
465  gax1->SetLabelFont(hPred[quantId]->GetXaxis()->GetLabelFont());
466  hPred[quantId]->GetXaxis()->SetLabelSize(0.);
467  hPred[quantId]->Draw("hist");
468  gax1->Draw("same");
469  }
470 
471  if(errorBand) PlotWithSystErrorBand_Quant(quantId, sPrediction[quantId], ups[quantId], dns[quantId], pot, color_pred, color_band, maxy, false);
472  else hPred[quantId]->Draw("hist same");
473  hWS[quantId]->Draw("hist same");
474  hAllBkg[quantId]->Draw("hist same");
475  hCosmics[quantId]->Draw("hist same");
476  hPred[quantId]->Draw("hist same");
477  drawQuantLabel(quantId);
478  if(quantId==1){
479  pad1->cd();
480  TLegend *legend = new TLegend(0.09,0.73,0.27,0.87,"","NDC");
481  //legend->AddEntry(hData[quantId], "FD Data","lep");
482  legend->AddEntry(hPred[quantId], "2017 Best Fit Pred.","l");
483  TLegendEntry *entry=legend->AddEntry("error","1-#sigma syst. range","bf");
484  entry->SetFillStyle(1001);
485  entry->SetFillColor(color_band);
486  entry->SetLineColor(color_band);
487  if(anaName=="fhc")legend->AddEntry(hWS[quantId], "Wrong Sign:#bar{#nu}_{#mu}CC","bf");
488  else legend->AddEntry(hWS[quantId], "Wrong Sign:#nu_{#mu}CC","bf");
489  legend->AddEntry(hAllBkg[quantId], "Total bkg.","bf");
490  legend->AddEntry(hCosmics[quantId], "Cosmic bkg.","bf");
491  legend->SetTextSize(0.025); //no border for legend
492  legend->SetBorderSize(0); //no border for legend
493  legend->SetFillStyle(0); //fill colour is white
494  legend->Draw();
495  CornerLabel(idName);
497  gPad->Update();
498  }
499  }// for quantiles
500 
501 
502  for(auto & p:{pad3,pad4,pad1,pad2}){
503  p->RedrawAxis();
504  }
505 
506  canvas->Update();
507  canvas->cd();
508  xTitle->Draw();
509  yTitle->Draw();
510  canvas->SaveAs(outName + "_quant14.png");
511  canvas->SaveAs(outName + "_quant14.eps");
512  canvas->SaveAs(outName + "_quant14.pdf");
513  file.open (outName + "_quant14.txt");
514  file<< "FD spectra showing the prediction at the 2017 best fit with systematics error band, and beam and cosmic backgrounds for the NuMI beam in " << anaName << " mode, separated per quartiles.";
515  file.close();
516 
517  } // plot quantiles 1 to 4
518 
519 
520 
521  if(DianaIsLazy && twobeams){
522 
523  std::ofstream file;
524  TCanvas *canvas;
525  TPad *pad1, *pad2, *pad3, *pad4;
526  SplitCanvasQuant(canvas, pad1, pad2, pad3, pad4);
527  gStyle->SetLineScalePS(3./paperscale);
528  canvas->SetCanvasSize(paperscale * canvas->GetWw(), paperscale * canvas->GetWh());
529  TString outName = dout_name + "fd_prediction_spectra__" + anaName+ "__2017calc";
530  std::string idName = anaName2;
531  for(int quantId:{7,8,5,6}){
532  if(quantId==5) pad1->cd();
533  if(quantId==6) pad2->cd();
534  if(quantId==7) pad3->cd();
535  if(quantId==8) pad4->cd();
536  HornSettings(quantId-1, pot, livetime, wrongs, maxy);
537  MakeHistCanvasReady_Quant(quantId, hPred[quantId], maxy);
538  if(errorBand) PlotWithSystErrorBand_Quant(quantId-4, sPrediction[quantId], ups[quantId], dns[quantId], pot, color_pred, color_band, maxy, true);
539  else hPred[quantId]->Draw("hist same");
540  hAllBkg[quantId]->Draw("hist same");
541  hCosmics[quantId]->Draw("hist same");
542  hPred[quantId]->Draw("hist same");
543  //grData[quantId]->Draw("p0 same");
544  drawQuantLabel(quantId-4);
545  if(quantId==1){
546  pad1->cd();
547  TLegend *legend = new TLegend(0.09,0.75,0.27,0.87,"","NDC");
548  legend->AddEntry(hPred[quantId], "2017 Best Fit Pred.","l");
549  TLegendEntry *entry=legend->AddEntry("error","1-#sigma syst. range","bf");
550  entry->SetFillStyle(1001);
551  entry->SetFillColor(color_band);
552  entry->SetLineColor(color_band);
553  legend->AddEntry(hAllBkg[quantId], "Total bkg.","bf");
554  legend->AddEntry(hCosmics[quantId], "Cosmic bkg.","bf");
555  legend->SetTextSize(0.025); //no border for legend
556  legend->SetBorderSize(0); //no border for legend
557  legend->SetFillStyle(0); //fill colour is white
558  legend->Draw();
559  CornerLabel(idName);
561  gPad->Update();
562  }
563  }// for quantiles
564 
565  for(auto & p:{pad3,pad4,pad1,pad2}){
566  p->RedrawAxis();
567  }
568 
569  canvas->Update();
570  canvas->cd();
571  xTitle->Draw();
572  yTitle->Draw();
573  canvas->SaveAs(outName + "_quant58.png");
574  canvas->SaveAs(outName + "_quant58.eps");
575  canvas->SaveAs(outName + "_quant58.pdf");
576  file.open (outName + "_quant58.txt");
577  file<< "FD spectra showing the prediction at the 2017 best fit with systematics error band, and beam and cosmic backgrounds for the NuMI beam in " << anaName << " mode, separated per quartiles.";
578  file.close();
579 
580  } // plot quantiles 1 to 4
581 
582 
583 
584 } // End of function
585 
Style_t marker_pred
Color_t color_data
Color_t color_pred
Pass neutrinos through unchanged.
Definition: IOscCalc.h:45
const double pot_fhc
void RestartCalculator(osc::IOscCalcAdjustable *calc, std::string anaName, bool nh)
double maxy_fhc
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
void ResetOscCalcToDefault(osc::IOscCalcAdjustable *calc)
Reset calculator to default assumptions for all parameters.
Definition: Calcs.cxx:23
Antineutrinos-only.
Definition: IPrediction.h:50
std::vector< SystGroupDef > systs
Definition: syst_header.h:385
void MakeHistCanvasReady_Quant(const int quant, TH1 *hist, double ymax)
Definition: Plots.cxx:1363
double maxy
correl_xv GetYaxis() -> SetDecimals()
int totquant
std::string extrapName
Loads shifted spectra from files.
Sign::Sign_t ws_one
Style_t line_data
Simple record of shifts applied to systematic parameters.
Definition: SystShifts.h:20
double pot_one
const double livetime_fhc
const char * p
Definition: xmltok.h:285
bool twobeams
T sqrt(T number)
Definition: d0nt_math.hpp:156
void CornerLabel(std::string &s)
Definition: numu_tools.h:145
Color_t color_bkg
double paperscale
Color_t color_band
void HornSettings(int quantile, double &pot, double &livetime, Sign::Sign_t &wrongs, double &maxy)
Style_t line_cos
void CenterTitles(TH1 *histo)
Definition: Plots.cxx:1481
double maxy_one
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
double maxy_rhc
Style_t line_bkg
Representation of a spectrum in any variable, with associated POT.
Definition: Spectrum.h:40
void PimpHist(TH1 *hist, Style_t linestyle, Color_t linecolor, int linewidth, Style_t markerstyle, Color_t markercolor, double markersize)
Pimp histogram once and for all.
Definition: Plots.cxx:1406
double pot_two
double livetime_one
void preliminary()
std::vector< const ISyst * > getAllAna2018Systs(const EAnaType2018 ana, const bool smallgenie, const BeamType2018 beam, bool isFit)
Charged-current interactions.
Definition: IPrediction.h:39
Style_t line_nc
double maxy_two
TGraph * graphAsymmErrorScaled(TH1 *histScaled, TH1 *hist, double overallScale)
Definition: Plots.cxx:978
double livetime
Color_t color_nc
correl_xv GetXaxis() -> SetDecimals()
virtual _IOscCalc< T > * Copy() const override
Definition: IOscCalc.h:49
Style_t marker_bkg
Style_t line_pred
std::string horn_two
void PreliminaryBoxOpening()
Definition: numu_tools.h:187
Style_t marker_ws
Oscillation probability calculators.
Definition: Calcs.h:5
unsigned int quantile(std::vector< std::vector< double >> v, float E, float hadEfrac)
Definition: Toy_analyses.C:480
OStream cout
Definition: OStream.cxx:6
std::vector< double > POT
Style_t line_ws
Style_t marker_nc
Color_t color_ws
Style_t marker_data
double livetime_two
Sign::Sign_t ws_two
std::string period_one
bool errorBand
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:541
TPaveText * drawQuantLabel(int quantId=0)
Definition: numu_tools.h:4
virtual void SetTh23(const T &th23)=0
Neutral-current interactions.
Definition: IPrediction.h:40
TFile * file
Definition: cellShifts.C:17
std::string horn_one
Both neutrinos and antineutrinos.
Definition: IPrediction.h:52
void plot_predictions(std::string anaName="fhc", std::string calculator="2018calc", bool nh=true, bool fake=true)
TPad * pad3
Definition: analysis.C:13
void SplitCanvasQuant(TCanvas *&canvas, TPad *&pad1, TPad *&pad2, TPad *&pad3, TPad *&pad4)
Definition: Plots.cxx:1437
double pot
bool DianaIsLazy
TPad * pad2
Definition: analysis.C:13
Color_t color_cos
All neutrinos, any flavor.
Definition: IPrediction.h:26
Style_t marker_cos
def entry(str)
Definition: HTMLTools.py:26
Sign::Sign_t wrongs
void SetShift(const ISyst *syst, double shift, bool force=false)
Definition: SystShifts.cxx:80
void FillWithDimColor(TH1 *h, bool usealpha, float dim)
double maxy_pad0
TPad * pad1
Definition: analysis.C:13
T asin(T number)
Definition: d0nt_math.hpp:60
enum BeamMode string