Functions | Variables
plot_datapredictions.C File Reference
#include "3FlavorAna/Ana2018/numu/includes.h"
#include "../settings.h"
#include "../numu_tools.h"
#include "3FlavorAna/Plotting/NuePlotStyle.h"

Go to the source code of this file.

Functions

void RestartCalculator (osc::IOscCalcAdjustable *calc, std::string anaName, bool nh)
 
void HornSettings (int quantile, double &pot, double &livetime, Sign::Sign_t &wrongs, double &maxy)
 
void plot_datapredictions (std::string anaName="fhc", std::string calculator="2018calc", bool nh=true)
 

Variables

std::string extrapName = "extrap"
 
bool DianaIsLazy = true
 
bool errorBand = true
 
double maxy = 0
 
double maxy_one = 4.999
 
double maxy_two = 5.999
 
double paperscale = 0.8/0.5
 
double maxy_pad0 = 0
 
double maxy_rhc = 9
 
double maxy_fhc = 13
 

Function Documentation

void HornSettings ( int  quantile,
double &  pot,
double &  livetime,
Sign::Sign_t wrongs,
double &  maxy 
)

Definition at line 57 of file plot_datapredictions.C.

References livetime_one, livetime_two, maxy_one, maxy_two, pot_one, pot_two, ws_one, and ws_two.

Referenced by plot_datapredictions().

57  {
58  if(quantile<4){
59  pot = pot_one;
61  wrongs = ws_one;
62  maxy = maxy_one;
63  }
64  else{
65  pot = pot_two;
67  wrongs = ws_two;
68  maxy = maxy_two;
69  }
70 
71 }// end HornSettings
double maxy_two
double pot_two
Definition: saveFDMCHists.C:20
double maxy
double pot_one
Definition: saveFDMCHists.C:19
double maxy_one
double livetime_two
Definition: saveFDMCHists.C:23
#define pot
Sign::Sign_t ws_one
Definition: saveFDMCHists.C:28
Sign::Sign_t ws_two
Definition: saveFDMCHists.C:29
unsigned int quantile(std::vector< std::vector< double >> v, float E, float hadEfrac)
Definition: Toy_analyses.C:480
double livetime
Definition: saveFDMCHists.C:21
double livetime_one
Definition: saveFDMCHists.C:22
Sign::Sign_t wrongs
Definition: saveFDMCHists.C:27
void plot_datapredictions ( std::string  anaName = "fhc",
std::string  calculator = "2018calc",
bool  nh = true 
)

Definition at line 75 of file plot_datapredictions.C.

References calc, make_mec_shifts_plots::canvas, ana::CenterTitles(), color_band, color_bkg, color_cos, color_data, color_nc, color_pred, color_ws, osc::_NoOscillations< T >::Copy(), CornerLabel(), make_associated_cosmic_defs::cosmics, om::cout, ana::DefaultOscCalc(), DianaIsLazy, drawQuantLabel(), allTimeWatchdog::endl, HTMLTools::entry(), errorBand, file, ana::FillWithDimColor(), ana::getAllAna2018Systs(), GetXaxis(), GetYaxis(), ana::graphAsymmErrorScaled(), horn_one, horn_two, HornSettings(), MECModelEnuComparisons::i, ana::Flavors::kAll, ana::Flavors::kAllNuMu, ana::Sign::kAntiNu, ana::Sign::kBoth, ana::Current::kCC, ana::Current::kNC, ana::kNuMu, ana::kNumuAna2018, make_mec_shifts_plots::legend, line_bkg, line_cos, line_data, line_nc, line_pred, line_ws, livetime, livetime_fhc, livetime_one, livetime_two, ana::MakeHistCanvasReady_Quant(), marker_bkg, marker_cos, marker_data, marker_nc, marker_pred, marker_ws, maxy, maxy_fhc, maxy_one, maxy_pad0, maxy_rhc, runNovaSAM::outName, pad1, pad2, pad3, paperscale, period_one, ana::PimpHist(), ana::PlotWithSystErrorBand_Quant(), POT, pot, pot_fhc, pot_one, pot_two, preliminary, PreliminaryBoxOpening(), runNovaSAM::release, RestartCalculator(), ana::SystShifts::SetShift(), ana::SplitCanvasQuant(), systs, totquant, twobeams, wrongs, and ws_one.

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

Definition at line 24 of file plot_datapredictions.C.

References std::asin(), delta_rhcfhc_joint, ana::ResetOscCalcToDefault(), osc::_IOscCalcAdjustable< T >::SetDmsq32(), osc::_IOscCalcAdjustable< T >::SetTh23(), sin_rhcfhc_joint, and std::sqrt().

Referenced by plot_datapredictions().

27 {
28 
29  // take different values from ../settings.h
30  // either numu only or joint BF
31  // we use the joint BF for neutrino
32  double sinHere = sin_rhcfhc_joint;
33  double deltaHere = delta_rhcfhc_joint;
34  // naively changing hierarchy
35  // this isnt totally right
36  /*
37  double deltaSign = 1.;
38  if(nh==false) deltaSign = -1.;
39 
40  if(anaName=="rhc"){
41  sinHere = sin_rhcfhc;
42  deltaHere = delta_rhcfhc*deltaSign;
43  }
44  if(anaName=="fhc"){
45  sinHere = sin_rhcfhc;
46  deltaHere = delta_rhcfhc*deltaSign;
47  }
48  */
50  calc->SetTh23(asin(sqrt(sinHere)));
51  calc->SetDmsq32(deltaHere);
52 
53 }
void ResetOscCalcToDefault(osc::IOscCalcAdjustable *calc)
Reset calculator to default assumptions for all parameters.
Definition: Calcs.cxx:23
T sqrt(T number)
Definition: d0nt_math.hpp:156
virtual void SetDmsq32(const T &dmsq32)=0
virtual void SetTh23(const T &th23)=0
const double delta_rhcfhc_joint
Definition: settings.h:26
const double sin_rhcfhc_joint
Definition: settings.h:25
T asin(T number)
Definition: d0nt_math.hpp:60

Variable Documentation

bool DianaIsLazy = true

Definition at line 12 of file plot_datapredictions.C.

Referenced by plot_datapredictions().

bool errorBand = true

Definition at line 13 of file plot_datapredictions.C.

Referenced by plot_datapredictions().

std::string extrapName = "extrap"
double maxy = 0
double maxy_fhc = 13

Definition at line 22 of file plot_datapredictions.C.

Referenced by plot_datapredictions().

double maxy_one = 4.999

Definition at line 16 of file plot_datapredictions.C.

Referenced by HornSettings(), and plot_datapredictions().

double maxy_pad0 = 0

Definition at line 20 of file plot_datapredictions.C.

Referenced by plot_datapredictions().

double maxy_rhc = 9

Definition at line 21 of file plot_datapredictions.C.

Referenced by plot_datapredictions().

double maxy_two = 5.999

Definition at line 17 of file plot_datapredictions.C.

Referenced by HornSettings().

double paperscale = 0.8/0.5

Definition at line 18 of file plot_datapredictions.C.

Referenced by plot_datapredictions().