plot_rationoosc.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 = false;
14 
15 double maxy = 0;
16 double maxy_one = 4.999;
17 double maxy_two = 5.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 = 35;
22 double maxy_fhc = 90;
23 
25  std::string anaName,
26  bool nh)
27 {
28 
29  double sinHere = sin_rhcfhc_joint;
30  double deltaHere = delta_rhcfhc_joint;
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_pad0/4.;
45  }
46  else{
47  pot = pot_two;
48  livetime = livetime_two;
49  wrongs = ws_two;
50  maxy = maxy_pad0/4.;
51  }
52 
53 }// end HornSettings
54 
55 
56 
57 void plot_rationoosc(std::string anaName = "fhc",
58  std::string calculator = "2018calc",
59  bool nh = true,
60  bool fake = true)
61 {
62 
63  ///////////////////////////////////////////////////////////////////////
64  // REMEMBER: //
65  // remove bkg from unoscillated ratio //
66  // remove bkg and cosmics from data for ratio without oscillations //
67  ///////////////////////////////////////////////////////////////////////
68 
69  bool preliminary = false;
70  std::string anaName2 ="";
71  if(anaName=="rhc"){
72  anaName2="Antineutrino beam";
73  totquant = 4;
74  twobeams = false;
76  }
77  // if fhc, change settings so the first 4 quantiles are fhc instead of rhc
78  // dont do anything if rhc as the first 4 quantiles are rhc already
79  if(anaName=="fhc"){
80  anaName2="Neutrino beam";
81  horn_one = "fhc";
82  period_one = "full";
83  pot_one = pot_fhc;
86  maxy_one = 5.999;
87  twobeams = false;
88  totquant = 4;
90  }
91 
92  // calculator
93  std::string calcName = "2018calc";
94  std::string hierarchy = "";
96  RestartCalculator(calc, anaName, nh);
98  auto calcNoOsc = osc.Copy();
99 
100 
101  // print on screen /////////////////////////////////////////////////////////
102  std::cout << "\n\nmake table with predicted number of events\n" << std::endl;
103  std::cout << "ana name: " << anaName << std::endl;
104  std::cout << "calculator: " << calcName << ", hierarchy: " << hierarchy << std::endl;
105 
106 
107  //inputs
108  TString dout_name = "../plots_blessing/";
109  std::string ddata_name = "/nova/ana/nu_mu_ana/Ana2018/Data/";
110  std::string dcosm_name = "/nova/ana/nu_mu_ana/Ana2018/Cosmics/";
111  std::string dpred_name = "/nova/ana/nu_mu_ana/Ana2018/Predictions/";
112  std::string fpred_name = "";
113  std::string fpred_suffix = "";
114 
115  TFile* fdata_one = new TFile((ddata_name + "fd_dataspectrum_" + horn_one + "_full__numu2018.root").c_str());
116  TFile* fdata_two = new TFile((ddata_name + "fd_dataspectrum_" + horn_two + "_full__numu2018.root").c_str());
117  std::vector<Spectrum*> s_realdata;
118  std::cout << "\nloading data files and spectra" << std::endl;
119  for(int quant=1; quant<totquant+1; quant++){
120  if(quant<5) s_realdata.push_back(LoadFrom<Spectrum>(fdata_one, Form("fd_data_q%d", quant) ).release());
121  else s_realdata.push_back(LoadFrom<Spectrum>(fdata_two, Form("fd_data_q%d", quant-4) ).release());
122  }
123 
124 
125  double pot1here = s_realdata[0]->POT();
126  std::cout << "\n\n pot for this sample = " << pot1here << "\n\n" << std::endl;
127 
128  // cosmics
129  std::string cosmics = "nocosmics";
130  std::vector<Spectrum> s_cosmics;
131  std::vector<TH1*> hCosmics;
132  std::cout << "\nloading cosmics" << std::endl;
133  std::cout << "" << horn_one << std::endl;
134  std::string fcosmics_name_one = dcosm_name + "cosmics_" + horn_one + "__numu2018.root";
135  TFile fcosmics_one(fcosmics_name_one.c_str());
136  hCosmics.push_back((TH1*)fcosmics_one.Get("cosmics_q0")); // remember, quant0 = all
137  hCosmics.push_back((TH1*)fcosmics_one.Get("cosmics_q1"));
138  hCosmics.push_back((TH1*)fcosmics_one.Get("cosmics_q2"));
139  hCosmics.push_back((TH1*)fcosmics_one.Get("cosmics_q3"));
140  hCosmics.push_back((TH1*)fcosmics_one.Get("cosmics_q4"));
141  auto *cosmics_quant1_one = (TH1*)fcosmics_one.Get("cosmics_q1");
142  auto *cosmics_quant2_one = (TH1*)fcosmics_one.Get("cosmics_q2");
143  auto *cosmics_quant3_one = (TH1*)fcosmics_one.Get("cosmics_q3");
144  auto *cosmics_quant4_one = (TH1*)fcosmics_one.Get("cosmics_q4");
145  s_cosmics.push_back(Spectrum (cosmics_quant1_one, pot_one, livetime_one));
146  s_cosmics.push_back(Spectrum (cosmics_quant2_one, pot_one, livetime_one));
147  s_cosmics.push_back(Spectrum (cosmics_quant3_one, pot_one, livetime_one));
148  s_cosmics.push_back(Spectrum (cosmics_quant4_one, pot_one, livetime_one));
149  if(twobeams){
150  std::string fcosmics_name_two = dcosm_name + "cosmics_" + horn_two + "__numu2018.root";
151  TFile fcosmics_two(fcosmics_name_two.c_str());
152  hCosmics.push_back((TH1*)fcosmics_two.Get("cosmics_q1"));
153  hCosmics.push_back((TH1*)fcosmics_two.Get("cosmics_q2"));
154  hCosmics.push_back((TH1*)fcosmics_two.Get("cosmics_q3"));
155  hCosmics.push_back((TH1*)fcosmics_two.Get("cosmics_q4"));
156  auto *cosmics_quant1_two = (TH1*)fcosmics_two.Get("cosmics_q1");
157  auto *cosmics_quant2_two = (TH1*)fcosmics_two.Get("cosmics_q2");
158  auto *cosmics_quant3_two = (TH1*)fcosmics_two.Get("cosmics_q3");
159  auto *cosmics_quant4_two = (TH1*)fcosmics_two.Get("cosmics_q4");
160  s_cosmics.push_back(Spectrum (cosmics_quant1_two, pot_two, livetime_two));
161  s_cosmics.push_back(Spectrum (cosmics_quant2_two, pot_two, livetime_two));
162  s_cosmics.push_back(Spectrum (cosmics_quant3_two, pot_two, livetime_two));
163  s_cosmics.push_back(Spectrum (cosmics_quant4_two, pot_two, livetime_two));
164  }
165  if(twobeams){
166  for(int quantId = 5; quantId < totquant+1; quantId++){
167  hCosmics[0]->Add(hCosmics[quantId]);
168  }
169  }
170 
171 
172  std::cout << "\nloading predictions" << std::endl;
173  const std::string dpred_name_fhc= dpred_name+"pred_numuconcat_fhc__numu2018.root";
174  const std::string dpred_name_rhc= dpred_name+"pred_numuconcat_rhc__numu2018.root";
175  std::string dpred_name_one = dpred_name_rhc;
176  std::string dpred_name_two = dpred_name_fhc;
177  if(anaName=="fhc") dpred_name_one = dpred_name_fhc;
178 
179  std::vector< const ISyst* > systs;
180  systs = getAllAna2018Systs(kNumuAna2018,true);
181  systs={};
182  int systSize = systs.size();
183 
184  std::cout << "\nloading predictions" << std::endl;
185  std::vector<PredictionSystJoint2018*> predictions;
186  ENu2018ExtrapType extrap = kNuMu;
187  for(int quant = 1; quant < totquant+1; quant++){
188  std::cout << "quantile " << quant << std::endl;
189  if(quant<5) predictions.push_back(new PredictionSystJoint2018(extrap, calc, horn_one, quant, systs, dpred_name_one));
190  else predictions.push_back(new PredictionSystJoint2018(extrap, calc, horn_two, quant-4, systs, dpred_name_two));
191  }
192 
193 
194  std::vector<TH1*> hDataUnscaled;
195  std::vector<TH1*> hData;
196  std::vector<TH1*> hNC;
197  std::vector<TH1*> hCC;
198  std::vector<TH1*> hNuMu;
199  std::vector<TH1*> hNoOsc;
200  std::vector<TH1*> hPred;
201  std::vector<TH1*> hBkg;
202  std::vector<TH1*> hAllBkg;
203  std::vector<TH1*> hWS;
204 
205 
207  hData.push_back(s_realdata[0]->ToTH1(s_realdata[0]->POT()));
208  hDataUnscaled.push_back(s_realdata[0]->ToTH1(s_realdata[0]->POT()));
209  hNC.push_back(predictions[0]->PredictComponent(calc, Flavors::kAll, Current::kNC, Sign::kBoth).ToTH1(pot));
210  hCC.push_back(predictions[0]->PredictComponent(calc, Flavors::kAll, Current::kCC, Sign::kBoth).ToTH1(pot));
211  hNuMu.push_back(predictions[0]->PredictComponent(calc, Flavors::kAllNuMu, Current::kCC, Sign::kBoth).ToTH1(pot));
212  hBkg.push_back(predictions[0]->Predict(calc).ToTH1(pot));
213  hAllBkg.push_back(predictions[0]->Predict(calc).ToTH1(pot));
214  hPred.push_back(predictions[0]->Predict(calc).ToTH1(pot));
215  hNoOsc.push_back(predictions[0]->Predict(calcNoOsc).ToTH1(pot));
216  hWS.push_back(predictions[0]->PredictComponent(calc, Flavors::kAllNuMu, Current::kCC, wrongs).ToTH1(pot));
217 
218  for(int quantId = 0; quantId < totquant; quantId++){
219  HornSettings(quantId, pot, livetime, wrongs, maxy);
220  hData.push_back(s_realdata[quantId]->ToTH1(s_realdata[0]->POT()));
221  hDataUnscaled.push_back(s_realdata[quantId]->ToTH1(s_realdata[0]->POT()));
222  hNC.push_back(predictions[quantId]->PredictComponent(calc, Flavors::kAll, Current::kNC, Sign::kBoth).ToTH1(pot));
223  hCC.push_back(predictions[quantId]->PredictComponent(calc, Flavors::kAll, Current::kCC, Sign::kBoth).ToTH1(pot));
224  hNuMu.push_back(predictions[quantId]->PredictComponent(calc, Flavors::kAllNuMu, Current::kCC, Sign::kBoth).ToTH1(pot));
225  hBkg.push_back(predictions[quantId]->Predict(calc).ToTH1(pot));
226  hAllBkg.push_back(predictions[quantId]->Predict(calc).ToTH1(pot));
227  hPred.push_back(predictions[quantId]->Predict(calc).ToTH1(pot));
228  hNoOsc.push_back(predictions[quantId]->Predict(calcNoOsc).ToTH1(pot));
229  hWS.push_back(predictions[quantId]->PredictComponent(calc, Flavors::kAllNuMu, Current::kCC, wrongs).ToTH1(pot));
230  if(quantId > 0){ // already added the first quantile
231  hDataUnscaled[0] ->Add(s_realdata[quantId]->ToTH1(s_realdata[0]->POT()));
232  hData[0] ->Add(s_realdata[quantId]->ToTH1(s_realdata[0]->POT()));
233  hNC[0] ->Add(predictions[quantId]->PredictComponent(calc, Flavors::kAll, Current::kNC, Sign::kBoth).ToTH1(pot));
234  hCC[0] ->Add(predictions[quantId]->PredictComponent(calc, Flavors::kAll, Current::kCC, Sign::kBoth).ToTH1(pot));
235  hNuMu[0] ->Add(predictions[quantId]->PredictComponent(calc, Flavors::kAllNuMu, Current::kCC, Sign::kBoth).ToTH1(pot));
236  hBkg[0] ->Add(predictions[quantId]->Predict(calc).ToTH1(pot));
237  hAllBkg[0]->Add(predictions[quantId]->Predict(calc).ToTH1(pot));
238  hPred[0] ->Add(predictions[quantId]->Predict(calc).ToTH1(pot));
239  hNoOsc[0] ->Add(predictions[quantId]->Predict(calcNoOsc).ToTH1(pot));
240  hWS[0] ->Add(predictions[quantId]->PredictComponent(calc, Flavors::kAllNuMu, Current::kCC, wrongs).ToTH1(pot));
241  }
242  }
243  // add bkg and cosmics to ws
244  // add cosmics to bkg
245  for(int quantId = 0; quantId < totquant+1; quantId++){
246  // we want to compare to pure CC beam
247  // in this case, ws doesnt count as bkg
248  hBkg[quantId]->Add(hNuMu[quantId], -1);
249  hBkg[quantId]->Add(hWS[quantId], -1);
250  hAllBkg[quantId]->Add(hNuMu[quantId],-1);
251  hAllBkg[quantId]->Add(hWS[quantId],-1);
252  hAllBkg[quantId]->Add(hCosmics[quantId]);
253  hPred[quantId]->Add(hCosmics[quantId]);
254  hNoOsc[quantId]->Add(hCosmics[quantId]);
255  }
256 
257 
258  //////////////////////////////////////
259  // actually do sth for no oscilations
260  /////////////////////////////////////
261  std::vector<TH1*> hNoOscDenom;
262  std::vector<TH1*> hPredRatio;
263  std::vector<TGraphAsymmErrors*> hDataRatio;
264 
265  for(int quantId = 0; quantId < totquant+1; quantId++){
266  // predictions ratio
267  TH1* tempdenom = (TH1*)hNoOsc[quantId]->Clone();
268  TH1* tempratio = (TH1*)hPred[quantId]->Clone();
269  hNoOscDenom.push_back(tempdenom);
270  hPredRatio.push_back(tempratio);
271  hNoOscDenom[quantId]->Add(hAllBkg[quantId], -1);
272  hPredRatio[quantId]->Add(hAllBkg[quantId], -1);
273  hPredRatio[quantId]->Divide(hNoOscDenom[quantId]);
274  hPredRatio[quantId]->GetYaxis()->SetTitle("Reconstructed Neutrino Energy (GeV)");
275  hPredRatio[quantId]->GetYaxis()->SetTitle("Ratio to no oscillation");
276  hPredRatio[quantId]->GetYaxis()->CenterTitle();
277  hPredRatio[quantId]->GetXaxis()->CenterTitle();
278 
279  // data ratio
280  TFeldmanCousins fc(0.6827);
281  hDataRatio.push_back(GraphWithPoissonErrors(hData[quantId]));
282  hDataRatio[quantId]->SetMarkerStyle(kFullCircle);
283  hDataRatio[quantId]->SetLineWidth(2);
284 
285  for(int bin = 0; bin < hData[quantId]->GetNbinsX(); bin++) {
286  double oldx, oldy, newy, bkg;
287  hDataRatio[quantId]->GetPoint(bin,oldx,oldy);
288  bkg = hAllBkg[quantId]->GetBinContent(bin+1);
289  newy = oldy - bkg;
290  hDataRatio[quantId]->SetPoint(bin,oldx,newy/hNoOscDenom[quantId]->GetBinContent(bin+1));
291  hDataRatio[quantId]->SetPointEYhigh(bin,(fc.CalculateUpperLimit(oldy,bkg)-newy)/hNoOscDenom[quantId]->GetBinContent(bin+1));
292  hDataRatio[quantId]->SetPointEYlow(bin,(newy-fc.CalculateLowerLimit(oldy,bkg))/hNoOscDenom[quantId]->GetBinContent(bin+1));
293  if(bin==0) {
294  hDataRatio[quantId]->SetPoint(bin,oldx,-50);
295  hDataRatio[quantId]->SetPointEYhigh(bin,0);
296  hDataRatio[quantId]->SetPointEYlow(bin,0);
297  }
298  } // for bins
299 
300  }// for quantiles
301 
302 
303 
304  // scale everything, everything!
305  std::cout << "\nScaling histograms to 0.1 GeV" << std::endl;
306  for(int quantId = 0; quantId < totquant+1; quantId++){
307  hData[quantId] ->Scale(0.1, "width");
308  hCosmics[quantId]->Scale(0.1, "width");
309  hNC[quantId] ->Scale(0.1, "width");
310  hCC[quantId] ->Scale(0.1, "width");
311  hBkg[quantId] ->Scale(0.1, "width");
312  hAllBkg[quantId] ->Scale(0.1, "width");
313  hNuMu[quantId] ->Scale(0.1, "width");
314  hPred[quantId] ->Scale(0.1, "width");
315  hNoOsc[quantId] ->Scale(0.1, "width");
316  hWS[quantId] ->Scale(0.1, "width");
317  }
318 
319  // data graphs
320  std::cout << "\nDefine data graphs\n" << std::endl;
321  std::vector<TGraph*> grData;
322  for(int quantId = 0; quantId < totquant+1; quantId++){
323  grData.push_back(graphAsymmErrorScaled(hData[quantId], hDataUnscaled[quantId]));
324  }
325 
326 
327  // pimping
328  for(int quantId = 0; quantId< 1; quantId++){
329  PimpHist(hData[quantId], line_data, color_data, 3, marker_data, color_data, 1);
330  PimpHist(hWS[quantId], line_ws, color_ws, 3, marker_ws, color_ws, 1);
331  PimpHist(hNC[quantId], line_nc, color_nc, 3, marker_nc, color_nc, 1);
332  PimpHist(hPred[quantId], line_pred, color_pred, 3, marker_pred, color_pred, 1);
333  PimpHist(hCosmics[quantId], line_cos, color_cos, 3, marker_cos, color_cos, 1);
334  PimpHist(hAllBkg[quantId], line_bkg, color_bkg, 3, marker_bkg, color_bkg, 1);
335  PimpHist(hBkg[quantId], line_bkg, color_bkg, 3, marker_bkg, color_bkg, 1);
336  PimpHist(hNoOsc[quantId], line_noosc, color_noosc, 3, marker_noosc, color_noosc, 1);
337  PimpHist(hPredRatio[quantId], line_pred, color_pred, 3, marker_pred, color_pred, 1);
338  FillWithDimColor(hWS[quantId], 0);
339  FillWithDimColor(hCosmics[quantId], 0);
340  FillWithDimColor(hAllBkg[quantId], 0);
341  grData[quantId]->SetMarkerSize(1.);
342  }// for quantId
343  //root is not behaving for the quantiles
344  //set different line width and marker size
345  for(int quantId = 1; quantId< totquant+1; quantId++){
346  PimpHist(hData[quantId], line_data, color_data, 3, marker_data, color_data, 1.);
347  PimpHist(hWS[quantId], line_ws, color_ws, 3, marker_ws, color_ws, 3);
348  PimpHist(hNC[quantId], line_nc, color_nc, 3, marker_nc, color_nc, 3);
349  PimpHist(hPred[quantId], line_pred, color_pred, 3, marker_pred, color_pred, 3);
350  PimpHist(hCosmics[quantId], line_cos, color_cos, 3, marker_cos, color_cos, 3);
351  PimpHist(hAllBkg[quantId], line_bkg, color_bkg, 3, marker_bkg, color_bkg, 3);
352  PimpHist(hBkg[quantId], line_bkg, color_bkg, 3, marker_bkg, color_bkg, 3);
353  PimpHist(hNoOsc[quantId], line_noosc, color_noosc, 3, marker_noosc, color_noosc, 3);
354  PimpHist(hPredRatio[quantId], line_pred, color_pred, 3, marker_pred, color_pred, 3);
355  FillWithDimColor(hWS[quantId], 0);
356  FillWithDimColor(hCosmics[quantId], 0);
357  FillWithDimColor(hAllBkg[quantId], 0);
358  grData[quantId]->SetMarkerSize(1.);
359  }// for quantId
360 
361 
362  // for later use in the 4pad plots
363  TLatex* xTitle = new TLatex(0.5, 0.03, "Reconstructed Energy (GeV)");
364  xTitle->SetNDC();
365  xTitle->SetTextAlign(22);
366  xTitle->SetTextSize(0.055 /paperscale);
367  TLatex* yTitle = new TLatex(0.02, 0.5, "Events / 0.1 GeV");
368  yTitle->SetNDC();
369  yTitle->SetTextAlign(22);
370  yTitle->SetTextAngle(90.);
371  yTitle->SetTextSize(0.055 / paperscale);
372 
373 
374  // ratios (spectra later)
375  if(DianaIsLazy){
376  std::ofstream file;
377  TString outName = dout_name + "fd_rationoosc_spectra__" + anaName;
378  std::string idName = anaName2;
379  int quantId = 0;
380  TCanvas *canvas = new TCanvas("canvas","canvas", 800, 500);
381  canvas->cd();
382  ana::CenterTitles(hPred[0]);
383  hPredRatio[0]->GetYaxis()->SetTitle("Ratio to no oscillation");
384  hPredRatio[0]->GetYaxis()->SetRangeUser(-0.1, 1.4);
385  hPredRatio[0]->GetXaxis()->SetLabelSize(hPredRatio[0]->GetXaxis()->GetLabelSize());
386  hPredRatio[0]->GetYaxis()->SetLabelSize(hPredRatio[0]->GetYaxis()->GetLabelSize());
387  hPredRatio[0]->Draw("hist");
388  hDataRatio[quantId]->Draw("p0 same");
389 
390 
391  TPaveText *pText0 = new TPaveText(0.35, 0.70, 0.45, 0.85, "NDC");
392  pText0->SetBorderSize(0);
393  pText0->SetFillStyle(0);
394  TText *text = pText0->AddText("All Quartiles");
395  text->SetTextAlign(22);
396  text->SetTextSize(0.05);
397  pText0->Draw();
398  //TLegend *legend = new TLegend(0.75,0.73,0.89,0.88);
399  TLegend *legend = new TLegend(0.32,0.58,0.45,0.73);
400  legend->AddEntry(hData[quantId], "Data","lep");
401  legend->AddEntry(hPred[quantId], "Prediction","l");
402  legend->SetTextSize(0.04); //no border for legend
403  legend->SetBorderSize(0); //no border for legend
404  legend->SetFillStyle(0); //fill colour is white
405  legend->Draw();
406  CornerLabel(idName);
408  //drawQuantLabel(0);
409  gPad->RedrawAxis();
410  gPad->Update();
411  canvas->cd();
412  canvas->SaveAs(outName + "__"+calcName+".eps");
413  canvas->SaveAs(outName + "__"+calcName+".png");
414  canvas->SaveAs(outName + "__"+calcName+".pdf");
415  file.open (outName + "__"+calcName+".txt");
416  file<< "Spectrum ratio of data over unoscillated prediction and best fit prediction over unoscillated prediction for the NuMI beam in " << anaName << " mode.";
417  file.close();
418  }
419 
420 
421 
422  const int cx = 800, cy = 960;
423  // all quantiles together
424  // equals quant 0
425  if(DianaIsLazy){
426  std::ofstream file;
427  TString outName = dout_name + "fd_dataprediction_noosc_spectra__" + anaName;
428  std::string idName = anaName2;
429  int quantId = 0;
430  TCanvas *canvas = new TCanvas("canvas","canvas", cx, cy);
431  canvas->cd();
432  ana::CenterTitles(hPred[0]);
433  hPred[0]->GetYaxis()->SetTitle("Events / 0.1 GeV");
434  hPred[0]->GetYaxis()->SetRangeUser(0.0, maxy_pad0);
435  hPred[0]->GetXaxis()->SetLabelSize(hPred[0]->GetXaxis()->GetLabelSize());
436  hPred[0]->GetYaxis()->SetLabelSize(hPred[0]->GetYaxis()->GetLabelSize());
437  hPred[0]->Draw("hist");
439  hPred[quantId]->Draw("hist same");
440  hAllBkg[quantId]->Draw("hist same");
441  grData[quantId]->Draw("p0 same");
442  grData[quantId]->Draw("p0 same");
443  hNoOsc[quantId]->Draw("hist same");
444  grData[quantId]->Draw("p0 same");
445  TLegend *legend = new TLegend(0.58,0.69,0.89,0.89);
446  legend->AddEntry(hData[quantId], "FD Data","lep");
447  legend->AddEntry(hPred[quantId], "Prediction","l");
448  legend->AddEntry(hNoOsc[quantId], "No oscillation","l");
449  legend->AddEntry(hAllBkg[quantId], "Total bkg.","l");
450  legend->SetTextSize(0.04); //no border for legend
451  legend->SetBorderSize(0); //no border for legend
452  legend->SetFillStyle(0); //fill colour is white
453  legend->Draw();
454  CornerLabel(idName);
456  drawQuantLabel(0);
457  gPad->RedrawAxis();
458  gPad->Update();
459  canvas->cd();
460 
461  canvas->SaveAs(outName + "__"+calcName+".eps");
462  canvas->SaveAs(outName + "__"+calcName+".png");
463  canvas->SaveAs(outName + "__"+calcName+".pdf");
464  file.open (outName + "__"+calcName+".txt");
465  file<< "FD spectra showing the prediction without oscillations, and beam background for the NuMI beam in " << anaName << " mode.";
466  file.close();
467 
468  } // plot all quantiles
469 
470 
471 
472 
473  if(DianaIsLazy){
474 
475  std::ofstream file;
476  TCanvas *canvas;
477  TPad *pad1, *pad2, *pad3, *pad4;
478  SplitCanvasQuant(canvas, pad1, pad2, pad3, pad4);
479  gStyle->SetLineScalePS(3./paperscale);
480  canvas->SetCanvasSize(paperscale * canvas->GetWw(), paperscale * canvas->GetWh());
481  TString outName = dout_name + "fd_dataprediction_noosc_spectra__" + anaName;
482  std::string idName = anaName2;
483  for(int quantId:{3,4,1,2}){
484  if(quantId==1) pad1->cd();
485  if(quantId==2) pad2->cd();
486  if(quantId==3) pad3->cd();
487  if(quantId==4) pad4->cd();
488  HornSettings(quantId-1, pot, livetime, wrongs, maxy);
489  MakeHistCanvasReady_Quant(quantId, hPred[quantId], maxy);
490  // hack to avoid 0-5 overlapping in the x-axis
491  hPred[quantId]->GetXaxis()->SetLabelSize(hPred[quantId]->GetXaxis()->GetLabelSize()/paperscale);
492  hPred[quantId]->GetYaxis()->SetLabelSize(hPred[quantId]->GetYaxis()->GetLabelSize()/paperscale);
493  hPred[quantId]->Draw("hist");
494  if(quantId==3){
495  TGaxis* gax1 = new TGaxis(0., 0., 4., 0., 0.,4.,4, "");
496  gax1->SetLabelSize(hPred[quantId]->GetXaxis()->GetLabelSize());
497  gax1->SetLabelOffset(hPred[quantId]->GetXaxis()->GetLabelOffset());
498  gax1->SetLabelFont(hPred[quantId]->GetXaxis()->GetLabelFont());
499  hPred[quantId]->GetXaxis()->SetLabelSize(0.);
500  hPred[quantId]->Draw("hist");
501  gax1->Draw("same");
502  }
503  hPred[quantId]->Draw("hist same");
504  hAllBkg[quantId]->Draw("hist same");
505  grData[quantId]->Draw("p0 same");
506  hNoOsc[quantId]->Draw("hist same");
507  drawQuantLabel(quantId);
508  grData[quantId]->Draw("p0 same");
509  if(quantId==2){
510  pad2->cd();
511  //TLegend *legend = new TLegend(0.09,0.75,0.27,0.89,"","NDC");
512  TLegend *legend = new TLegend(0.75,0.62,0.89,0.75,"","NDC");
513  legend->AddEntry(hData[quantId], "FD Data","lep");
514  legend->AddEntry(hPred[quantId], "Prediction","l");
515  legend->AddEntry(hNoOsc[quantId], "No oscillation","l");
516  legend->AddEntry(hAllBkg[quantId], "Total bkg.","l");
517  legend->SetTextSize(0.025); //no border for legend
518  legend->SetBorderSize(0); //no border for legend
519  legend->SetFillStyle(0); //fill colour is white
520  legend->Draw();
521  }
522  CornerLabel(idName);
524  gPad->Update();
525  }// for quantiles
526 
527  for(auto & p:{pad3,pad4,pad1,pad2}){
528  p->RedrawAxis();
529  }
530  canvas->Update();
531  canvas->cd();
532  xTitle->Draw();
533  yTitle->Draw();
534  canvas->SaveAs(outName + "_quant14__"+calcName+".png");
535  canvas->SaveAs(outName + "_quant14__"+calcName+".eps");
536  canvas->SaveAs(outName + "_quant14__"+calcName+".pdf");
537  file.open (outName + "_quant14__"+calcName+".txt");
538  file<< "FD spectra showing data, prediction at best fit, predictions without oscillations, and beam background for the NuMI beam in " << anaName << " mode, separated per quartiles.";
539  file.close();
540 
541  } // plot quantiles 1 to 4
542 
543 
544 
545 
546  ////////////////////////////////
547  ///// no oscillation only //
548  ///////////////////////////////
549  // all quantiles together
550  // equals quant 0
551  if(DianaIsLazy){
552  std::ofstream file;
553  TString outName = dout_name + "fd_prediction_noosc_spectra__" + anaName;
554  std::string idName = anaName2;
555  int quantId = 0;
556  TCanvas *canvas = new TCanvas("canvas","canvas", cx, cy);
557  canvas->cd();
558  ana::CenterTitles(hNoOsc[0]);
559  hNoOsc[0]->GetYaxis()->SetTitle("Events / 0.1 GeV");
560  hNoOsc[0]->GetYaxis()->SetRangeUser(0.0, maxy_pad0);
561  hNoOsc[0]->GetXaxis()->SetLabelSize(hNoOsc[0]->GetXaxis()->GetLabelSize());
562  hNoOsc[0]->GetYaxis()->SetLabelSize(hNoOsc[0]->GetYaxis()->GetLabelSize());
563  hNoOsc[0]->Draw("hist same");
564  hPred[0]->Draw("hist same");
566  hNoOsc[quantId]->Draw("hist same");
567  TLegend *legend = new TLegend(0.58,0.79,0.89,0.89);//(.89-.69)->4
568  legend->AddEntry(hPred[quantId], "Prediction","l");
569  legend->AddEntry(hNoOsc[quantId], "No oscillation","l");
570  legend->SetTextSize(0.04);
571  legend->SetBorderSize(0);
572  legend->SetFillStyle(0);
573  legend->Draw();
574  CornerLabel(idName);
576  //drawQuantLabel(0);
577  gPad->RedrawAxis();
578  gPad->Update();
579  canvas->cd();
580  canvas->SaveAs(outName + "__"+calcName+".eps");
581  canvas->SaveAs(outName + "__"+calcName+".png");
582  canvas->SaveAs(outName + "__"+calcName+".pdf");
583  file.open (outName + "__"+calcName+".txt");
584  file<< "FD spectra showing the prediction with neutrino oscillations (violet) and without oscillations (green) for the NuMI beam in " << anaName << " mode for the 2018 exposure.";
585  file.close();
586 
587  } // plot all quantiles
588 
589 
590 
591 
592 
593  if(DianaIsLazy){
594 
595  std::ofstream file;
596  TCanvas *canvas;
597  TPad *pad1, *pad2, *pad3, *pad4;
598  SplitCanvasQuant(canvas, pad1, pad2, pad3, pad4);
599  gStyle->SetLineScalePS(3./paperscale);
600  canvas->SetCanvasSize(paperscale * canvas->GetWw(), paperscale * canvas->GetWh());
601  TString outName = dout_name + "fd_prediction_noosc_spectra__" + anaName;
602  std::string idName = anaName2;
603  for(int quantId:{3,4,1,2}){
604  if(quantId==1) pad1->cd();
605  if(quantId==2) pad2->cd();
606  if(quantId==3) pad3->cd();
607  if(quantId==4) pad4->cd();
608  HornSettings(quantId-1, pot, livetime, wrongs, maxy);
609  MakeHistCanvasReady_Quant(quantId, hNoOsc[quantId], maxy);
610  // hack to avoid 0-5 overlapping in the x-axis
611  hNoOsc[quantId]->GetXaxis()->SetLabelSize(hNoOsc[quantId]->GetXaxis()->GetLabelSize()/paperscale);
612  hNoOsc[quantId]->GetYaxis()->SetLabelSize(hNoOsc[quantId]->GetYaxis()->GetLabelSize()/paperscale);
613  hNoOsc[quantId]->Draw("hist same");
614  hPred[quantId]->Draw("hist same");
615  if(quantId==3){
616  TGaxis* gax1 = new TGaxis(0., 0., 4., 0., 0.,4.,4, "");
617  gax1->SetLabelSize(hNoOsc[quantId]->GetXaxis()->GetLabelSize());
618  gax1->SetLabelOffset(hNoOsc[quantId]->GetXaxis()->GetLabelOffset());
619  gax1->SetLabelFont(hNoOsc[quantId]->GetXaxis()->GetLabelFont());
620  hNoOsc[quantId]->GetXaxis()->SetLabelSize(0.);
621  hNoOsc[quantId]->Draw("hist same");
622  hPred[quantId]->Draw("hist same");
623  gax1->Draw("same");
624  }
625 
626  hNoOsc[quantId]->Draw("hist same");
627  hPred[quantId]->Draw("hist same");
628  drawQuantLabel(quantId);
629  if(quantId==2){
630  pad2->cd();
631  //TLegend *legend = new TLegend(0.09,0.75,0.27,0.89,"","NDC");
632  TLegend *legend = new TLegend(0.75,0.685,0.89,0.75,"","NDC"); // (75-62)->4
633  legend->AddEntry(hPred[quantId], "Prediction","l");
634  legend->AddEntry(hNoOsc[quantId], "No oscillation","l");
635  legend->SetTextSize(0.025); //no border for legend
636  legend->SetBorderSize(0); //no border for legend
637  legend->SetFillStyle(0); //fill colour is white
638  legend->Draw();
639  }
640  CornerLabel(idName);
642  gPad->Update();
643  }// for quantiles
644 
645 
646  for(auto & p:{pad3,pad4,pad1,pad2}){
647  p->RedrawAxis();
648  }
649 
650  canvas->Update();
651  canvas->cd();
652  xTitle->Draw();
653  yTitle->Draw();
654  canvas->SaveAs(outName + "_quant14__"+calcName+".png");
655  canvas->SaveAs(outName + "_quant14__"+calcName+".eps");
656  canvas->SaveAs(outName + "_quant14__"+calcName+".pdf");
657  file.open (outName + "_quant14__"+calcName+".txt");
658  file<< "FD spectra showing the prediction with neutrino oscillations (violet) and without oscillations (green) for the NuMI beam in " << anaName << " mode, separated per quartiles, for the 2018 exposure.";
659  file.close();
660 
661  } // plot quantiles 1 to 4
662 
663 
664 
665 
666 } // End of function
667 
void PreliminaryNoOsc()
Definition: numu_tools.h:206
Pass neutrinos through unchanged.
Definition: IOscCalc.h:45
Style_t marker_cos
Definition: saveFDMCHists.C:63
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
correl_xv GetYaxis() -> SetDecimals()
Color_t color_noosc
Definition: settings.h:64
Loads shifted spectra from files.
double pot_two
Definition: saveFDMCHists.C:20
int totquant
Definition: saveFDMCHists.C:40
void PreliminaryNoOscQuant()
Definition: numu_tools.h:223
const char * p
Definition: xmltok.h:285
T sqrt(T number)
Definition: d0nt_math.hpp:156
Color_t color_pred
Definition: saveFDMCHists.C:47
void CornerLabel(std::string &s)
Definition: numu_tools.h:145
Style_t marker_data
Definition: saveFDMCHists.C:61
const double livetime_fhc
Definition: saveFDMCHists.C:12
Style_t line_bkg
Definition: saveFDMCHists.C:58
void CenterTitles(TH1 *histo)
Definition: Plots.cxx:1481
TGraphAsymmErrors * GraphWithPoissonErrors(const TH1 *h, bool noErrorsXaxis, bool drawEmptyBins)
Calculate statistical errors appropriate for small Poisson numbers.
Definition: Plots.cxx:910
Style_t line_data
Definition: saveFDMCHists.C:55
Color_t color_data
Definition: saveFDMCHists.C:50
osc::OscCalcDumb calc
osc::IOscCalcAdjustable * DefaultOscCalc()
Create a new calculator with default assumptions for all parameters.
Definition: Calcs.cxx:49
double pot_one
Definition: saveFDMCHists.C:19
virtual void SetDmsq32(const T &dmsq32)=0
bool twobeams
Definition: saveFDMCHists.C:39
Representation of a spectrum in any variable, with associated POT.
Definition: Spectrum.h:40
Style_t line_pred
Definition: saveFDMCHists.C:56
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
Style_t line_cos
Definition: saveFDMCHists.C:57
double maxy_rhc
double livetime_two
Definition: saveFDMCHists.C:23
std::vector< const ISyst * > getAllAna2018Systs(const EAnaType2018 ana, const bool smallgenie, const BeamType2018 beam, bool isFit)
Charged-current interactions.
Definition: IPrediction.h:39
TGraph * graphAsymmErrorScaled(TH1 *histScaled, TH1 *hist, double overallScale)
Definition: Plots.cxx:978
Style_t marker_noosc
Definition: settings.h:78
correl_xv GetXaxis() -> SetDecimals()
bool DianaIsLazy
virtual _IOscCalc< T > * Copy() const override
Definition: IOscCalc.h:49
Color_t color_ws
Definition: saveFDMCHists.C:48
#define pot
void HornSettings(int quantile, double &pot, double &livetime, Sign::Sign_t &wrongs, double &maxy)
Sign::Sign_t ws_one
Definition: saveFDMCHists.C:28
double paperscale
Int_t preliminary
Definition: SimpleIterate.C:63
bool errorBand
float bin[41]
Definition: plottest35.C:14
Color_t color_nc
Definition: saveFDMCHists.C:49
Sign::Sign_t ws_two
Definition: saveFDMCHists.C:29
Oscillation probability calculators.
Definition: Calcs.h:5
Style_t marker_pred
Definition: saveFDMCHists.C:62
void plot_rationoosc(std::string anaName="fhc", std::string calculator="2018calc", bool nh=true, bool fake=true)
double maxy
unsigned int quantile(std::vector< std::vector< double >> v, float E, float hadEfrac)
Definition: Toy_analyses.C:480
double maxy_fhc
OStream cout
Definition: OStream.cxx:6
std::vector< double > POT
double maxy_pad0
const double pot_fhc
Definition: saveFDMCHists.C:10
std::string extrapName
double maxy_one
Style_t line_noosc
Definition: settings.h:71
void PreliminaryNoOscRatio()
Definition: numu_tools.h:240
double livetime
Definition: saveFDMCHists.C:21
double livetime_one
Definition: saveFDMCHists.C:22
Sign::Sign_t wrongs
Definition: saveFDMCHists.C:27
TPaveText * drawQuantLabel(int quantId=0)
Definition: numu_tools.h:4
virtual void SetTh23(const T &th23)=0
const double delta_rhcfhc_joint
Definition: settings.h:26
Color_t color_cos
Definition: saveFDMCHists.C:51
Neutral-current interactions.
Definition: IPrediction.h:40
Style_t marker_bkg
Definition: saveFDMCHists.C:64
TFile * file
Definition: cellShifts.C:17
Both neutrinos and antineutrinos.
Definition: IPrediction.h:52
void RestartCalculator(osc::IOscCalcAdjustable *calc, std::string anaName, bool nh)
TPad * pad3
Definition: analysis.C:13
const double sin_rhcfhc_joint
Definition: settings.h:25
Style_t line_nc
Definition: saveFDMCHists.C:54
void SplitCanvasQuant(TCanvas *&canvas, TPad *&pad1, TPad *&pad2, TPad *&pad3, TPad *&pad4)
Definition: Plots.cxx:1437
TPad * pad2
Definition: analysis.C:13
Style_t marker_ws
Definition: saveFDMCHists.C:59
All neutrinos, any flavor.
Definition: IPrediction.h:26
double maxy_two
std::string period_one
Definition: saveFDMCHists.C:35
Style_t marker_nc
Definition: saveFDMCHists.C:60
void FillWithDimColor(TH1 *h, bool usealpha, float dim)
std::string horn_one
Definition: saveFDMCHists.C:33
Style_t line_ws
Definition: saveFDMCHists.C:53
Color_t color_bkg
Definition: saveFDMCHists.C:52
std::string horn_two
Definition: saveFDMCHists.C:34
TPad * pad1
Definition: analysis.C:13
T asin(T number)
Definition: d0nt_math.hpp:60
enum BeamMode string