plot_recoE_numu.C
Go to the documentation of this file.
1 //------------------------------------------------------------------------------
2 // Plotting macro to make the numu box opening/blessed plots based on older
3 // versions from previous analyses. This was meant to be general purpose but
4 // as requests for more plots came in it quickly got more entangled
5 //
6 // Options are pretty self explanatory but you can refer to docDB 45809
7 // to see the plots that are produced.
8 //
9 // The error bars on the data are generated using fns in CAFAna/Analysis/Plots.h
10 // which use TFeldmanCousins::Calculate{Upper,Lower}Limit set to 1sigma
11 //
12 // That function takes the number of observed (data) points and an estimate of
13 // the background. In the past we've done mixed things sometimes using a bkgd
14 // of 0 and sometimes using the estimate from MC. I've decided to use the MC
15 // estimate in all cases here. This has the effect of slightly increasing the
16 // error bars in the bins with significant bkgd
17 //
18 // Andrew Sutton - ats2mk@virginia.edu
19 //------------------------------------------------------------------------------
20 
21 #include "plotting_header.h"
22 #include "TRatioPlot.h"
23 
24 using namespace ana;
25 
27  bool drawData = true,
28  bool drawPred = true,
29  bool drawSystBand = true,
30  bool useBF = true,
31  bool makeUnoscPlots = true,
32  bool makeUnoscRatio = true,
33  bool isPreliminary = true)
34 {
35  //----------------------------------------------------------------------------
36  // Options setup
37  //----------------------------------------------------------------------------
38  uint nQuantiles = 4;
39  double maxPredScale;
40  double quantPredScale;
41  double pot, lt;
42  bool isFHC;
44  std::string HORN;
45  Sign::Sign_t WS;
46  std::string sWSLabel;
47  if ( horn == "fhc" ) {
48  maxPredScale = 1.5;
49  quantPredScale = maxPredScale * 2.;
50  pot = kAna2020FHCPOT;
52  isFHC = true;
53  Beam = BeamType2020::kFHC;
54  HORN = "FHC";
55  WS = Sign::kAntiNu;
56  sWSLabel = "Wrong Sign: #bar{#nu}_{#mu}CC";
57  } else {
58  maxPredScale = 1.75;
59  quantPredScale = maxPredScale * 2.2;
60  pot = kAna2020RHCPOT;
62  isFHC = false;
63  Beam = BeamType2020::kRHC;
64  HORN = "RHC";
65  WS = Sign::kNu;
66  sWSLabel = "Wrong Sign: #nu_{#mu}CC";
67  }
68 
69  if (!drawPred && drawSystBand) {
70  std::cout << "Not drawing oscillated predicion but asking to draw the syst bands...no" << std::endl;
71  std::abort();
72  }
73 
74  // -- Set the fle append structure for outputs
75  std::string sDataPred = "_DataPred";
76  std::string sSystBand = "";
77  if (drawSystBand) sSystBand = "_SystBand";
78  if (!drawData) sDataPred = "_PredOnly";
79  if (!drawPred) sDataPred = "_DataOnly";
80 
81  //----------------------------------------------------------------------------
82  // Dir and file names
83  //----------------------------------------------------------------------------
84  std::string sOutDir = "";
85  std::string sHieOct = "NH_UO";
86  std::string sPredDir = "/nova/ana/3flavor/Ana2020/Predictions/";
87  std::string sCosDir = "/nova/ana/3flavor/Ana2020/Predictions/";
88  std::string sDataDir = "/nova/ana/3flavor/Ana2020/Data/";
89  std::string sFitDir = "/nova/ana/3flavor/Ana2020/Fits/BestFits/";
90  std::string sPredFile = sPredDir + "pred_NumuEnergy_AllQuantiles_full_" + HORN + "AllSysts_numuconcat_20200428.root";
91  std::string sCosFile = sCosDir + "fd_cosmic_pred_numu_" + horn + "_v2.root";
92  std::string sDataFile = sDataDir + "all_2020_data.root";
93  std::string sFitFile = sFitDir + "bestfits_joint_realData_both_systs.root";
94 
95 
96  //----------------------------------------------------------------------------
97  // Initial setup, calcs and systs
98  //----------------------------------------------------------------------------
99  std::vector< const ISyst * > systs = {};
100  systs = get3FlavorAna2020AllSysts(EAnaType2020::kNumuAna, true, Beam, false, true);
101 
102  osc::NoOscillations base_NoOsc;
103  auto calcNO = base_NoOsc.Copy();
104  auto calc = DefaultOscCalc();
105 
106  // Best fit stuff
107  // Pull BF calc and shifts from file
108  // or you can set the calc above
109  SystShifts bfshifts = SystShifts();
110  if (useBF) {
111  TFile *fBestFit = new TFile(sFitFile.c_str(), "READ");
112  calc = (osc::IOscCalcAdjustable*)(LoadFrom<osc::IOscCalc>(fBestFit, (sHieOct + "_calc").c_str())).release();
113 
114  // Only used the BF systs that we loaded in using get3Flavor...
115  auto systsBF = SystShifts::LoadFrom(fBestFit, (sHieOct + "_systs").c_str()).release();
116  for (auto bfsyst : systsBF->ActiveSysts())
117  if (std::find(systs.begin(), systs.end(), bfsyst) != systs.end())
118  bfshifts.SetShift(bfsyst, systsBF->GetShift(bfsyst));
119 
120  fBestFit->Close();
121  }// endif useBF
122 
123  //----------------------------------------------------------------------------
124  // Load things in
125  //----------------------------------------------------------------------------
126  // Cosmics and Data
127  TFile *fCosmics = new TFile(sCosFile.c_str() , "READ");
128  TFile *fData = 0;
129  if (drawData) fData = new TFile(sDataFile.c_str(), "READ");
130 
131  std::vector<TH1*> hCosm;
132  std::vector<TH1*> hData;
133  std::vector<TH1*> hDataUnscaled;
134  for (uint q = 0; q < nQuantiles+1; ++q) {
135  if (drawData) {
136  Spectrum* data = Spectrum::LoadFrom( fData, ("numu_" + horn + "_q" + std::to_string(q)).c_str()) .release();
137  hData.push_back( data->ToTH1(data->POT()) );
138  hDataUnscaled.push_back( data->ToTH1(data->POT()) );
139  }
140 
141  Spectrum* cosm = Spectrum::LoadFrom( fCosmics, Form("numu_cos_q%i",q)) .release();
142  hCosm.push_back( cosm->ToTH1( lt, kLivetime ) );
143  }
144  fCosmics->Close();
145  if (drawData) fData->Close();
146 
147  // Predictions
148  std::vector<PredictionSyst3Flavor2020*> vPred;
149  for (uint q = 0; q < nQuantiles; ++q)
150  vPred.push_back( new PredictionSyst3Flavor2020(ENu2020ExtrapType::kNuMuPtExtrap, calc, horn, q+1, systs, false, "") );
151 
152  //----------------------------------------------------------------------------
153  // Turn into hists
154  // have to do some fancy footwork to make element 0 the sum of all quants
155  //----------------------------------------------------------------------------
156  std::vector<TH1*> hNoOsc;
157  std::vector<TH1*> hNOBkg;
158  std::vector<TH1*> hNONumu;
159  std::vector<TH1*> hPred;
160  std::vector<TH1*> hBkg;
161  std::vector<TH1*> hNumu;
162  std::vector<TH1*> hWS;
163 
164  // Initialize the 0th element with Q1
165  std::string WSH = "WS" + HORN;
166  hNoOsc.push_back ( vPred[0]->PredictComponentSyst(calcNO, {} , predSels.at("All").flav , predSels.at("All").curr , predSels.at("All").sign ).ToTH1(pot) );
167  hNOBkg.push_back ( vPred[0]->PredictComponentSyst(calcNO, {} , predSels.at("All").flav , predSels.at("All").curr , predSels.at("All").sign ).ToTH1(pot) );
168  hNONumu.push_back( vPred[0]->PredictComponentSyst(calcNO, {} , predSels.at("Numu").flav, predSels.at("Numu").curr, predSels.at("Numu").sign).ToTH1(pot) );
169  hPred.push_back ( vPred[0]->PredictComponentSyst(calc , bfshifts, predSels.at("All").flav , predSels.at("All").curr , predSels.at("All").sign ).ToTH1(pot) );
170  hBkg.push_back ( vPred[0]->PredictComponentSyst(calc , bfshifts, predSels.at("All").flav , predSels.at("All").curr , predSels.at("All").sign ).ToTH1(pot) );
171  hNumu.push_back ( vPred[0]->PredictComponentSyst(calc , bfshifts, predSels.at("Numu").flav, predSels.at("Numu").curr, predSels.at("Numu").sign).ToTH1(pot) );
172  hWS.push_back ( vPred[0]->PredictComponentSyst(calc , bfshifts, predSels.at(WSH).flav , predSels.at(WSH).curr , predSels.at(WSH).sign ).ToTH1(pot) );
173 
174  // Now do the rest of the quants
175  for (uint q = 0; q < nQuantiles; ++q) {
176  hNoOsc.push_back ( vPred[q]->PredictComponentSyst(calcNO, {} , predSels.at("All").flav , predSels.at("All").curr , predSels.at("All").sign ).ToTH1(pot) );
177  hNOBkg.push_back ( vPred[q]->PredictComponentSyst(calcNO, {} , predSels.at("All").flav , predSels.at("All").curr , predSels.at("All").sign ).ToTH1(pot) );
178  hNONumu.push_back( vPred[q]->PredictComponentSyst(calcNO, {} , predSels.at("Numu").flav, predSels.at("Numu").curr, predSels.at("Numu").sign).ToTH1(pot) );
179  hPred.push_back ( vPred[q]->PredictComponentSyst(calc , bfshifts, predSels.at("All").flav , predSels.at("All").curr , predSels.at("All").sign ).ToTH1(pot) );
180  hBkg.push_back ( vPred[q]->PredictComponentSyst(calc , bfshifts, predSels.at("All").flav , predSels.at("All").curr , predSels.at("All").sign ).ToTH1(pot) );
181  hNumu.push_back ( vPred[q]->PredictComponentSyst(calc , bfshifts, predSels.at("Numu").flav, predSels.at("Numu").curr, predSels.at("Numu").sign).ToTH1(pot) );
182  hWS.push_back ( vPred[q]->PredictComponentSyst(calc , bfshifts, predSels.at(WSH).flav , predSels.at(WSH).curr , predSels.at(WSH).sign ).ToTH1(pot) );
183 
184  // Tack onto element 0 for other quants
185  if (q > 0){
186  hNoOsc [0]->Add(hNoOsc [q+1]);
187  hNOBkg [0]->Add(hNOBkg [q+1]);
188  hNONumu[0]->Add(hNONumu[q+1]);
189  hPred [0]->Add(hPred [q+1]);
190  hBkg [0]->Add(hBkg [q+1]);
191  hNumu [0]->Add(hNumu [q+1]);
192  hWS [0]->Add(hWS [q+1]);
193  }
194  }// end quant loop
195 
196  // Make the necessary combos for the custom HistStack then scale
197  // HistStack from bottom to top is cosm, bkg, ws (draw them in reverse order)
198  std::vector<TGraph*> gData;
199  for (uint q = 0; q < nQuantiles+1; ++q) {
200  hBkg [q]->Add(hNumu [q], -1); // remove the signal for bkg
201  hNOBkg[q]->Add(hNONumu[q], -1); // remove the signal for bkg
202 
203  std::cout << "\n*************************";
204  if (q == 0) std::cout << "\nAll Quantiles";
205  else std::cout << "\nQuantile " << q;
206  std::cout << "\n NoOsc : " << hNoOsc[q]->Integral();
207  if (drawData) {
208  std::cout << "\n Data : " << hData[q]->Integral();
209  }
210  std::cout << "\n Pred : " << hPred[q]->Integral()
211  << "\n Signal : " << hNumu[q]->Integral()
212  << "\n Cosmics : " << hCosm[q]->Integral()
213  << "\n Other Bkgd: " << hBkg [q]->Integral()
214  << "\n Wrong Sign: " << hWS [q]->Integral()
215  << "\n*************************\n"
216  << std::endl;
217 
218  // Add backgrounds
219  hBkg [q]->Add(hCosm[q], +1); // add on the cosmics for stack
220  hWS [q]->Add(hBkg [q], +1); // add bkg to WS for stack
221  hPred [q]->Add(hCosm[q], +1); // add cosmics to full prediction
222  hNoOsc[q]->Add(hCosm[q], +1); // add cosmics to NoOsc prediction
223  hNOBkg[q]->Add(hCosm[q], +1); // add cosmics to NoOsc prediction
224 
225 
226  // Scale everything
227  if (drawData) {
228  hData [q]->Scale(0.1, "width");
229  }
230 
231  hCosm [q]->Scale(0.1, "width");
232  hNoOsc [q]->Scale(0.1, "width");
233  hNOBkg [q]->Scale(0.1, "width");
234  hNONumu[q]->Scale(0.1, "width");
235  hPred [q]->Scale(0.1, "width");
236  hBkg [q]->Scale(0.1, "width");
237  hNumu [q]->Scale(0.1, "width");
238  hWS [q]->Scale(0.1, "width");
239 
240  // Pimp it out
241  if (q == 0) {
242  if (drawData)
243  PimpHist(hData [q], kSolid , kBlack, 2, kFullCircle , kBlack , 1);
244 
245  PimpHist(hCosm [q], kSolid, kCosmicBackgroundColor, 1, kFullCircle, kCosmicBackgroundColor, 1);
246  PimpHist(hNoOsc[q], kSolid, kNoOscColor , 3, kFullCircle, kNoOscColor , 1);
247  PimpHist(hPred [q], kSolid, kFullPredColor , 3, kFullCircle, kFullPredColor , 1);
248  PimpHist(hBkg [q], kSolid, kTotalBkgColor , 1, kFullCircle, kTotalBkgColor , 1);
249  PimpHist(hWS [q], kSolid, kNumuWSColor , 1, kFullCircle, kNumuWSColor , 1);
250  } else {
251  if (drawData)
252  PimpHist(hData [q], kSolid , kBlack, 2, kFullCircle , kBlack , 2);
253 
254  PimpHist(hCosm [q], kSolid, kCosmicBackgroundColor, 2, kFullCircle, kCosmicBackgroundColor, 2);
255  PimpHist(hNoOsc[q], kSolid, kNoOscColor , 2, kFullCircle, kNoOscColor , 2);
256  PimpHist(hPred [q], kSolid, kFullPredColor , 2, kFullCircle, kFullPredColor , 2);
257  PimpHist(hBkg [q], kSolid, kTotalBkgColor , 2, kFullCircle, kTotalBkgColor , 2);
258  PimpHist(hWS [q], kSolid, kNumuWSColor , 2, kFullCircle, kNumuWSColor , 2);
259  }
260 
261  FillWithDimColor(hWS [q], 0);
262  FillWithDimColor(hCosm[q], 0);
263  FillWithDimColor(hBkg [q], 0);
264 
265  // Set the maximum values for plotting
266  if (q == 0) {
267  hNoOsc[0]->SetMaximum(hNoOsc[0]->GetMaximum() * maxPredScale);
268  hPred [0]->SetMaximum(hPred [0]->GetMaximum() * maxPredScale);
269  if (drawData)
270  hData [0]->SetMaximum(hPred [0]->GetMaximum());
271  } else if (q == 1) {
272  hNoOsc[1]->SetMaximum(hNoOsc[1]->GetMaximum() * maxPredScale);
273  hPred [1]->SetMaximum(hPred [1]->GetMaximum() * quantPredScale);
274  if (drawData)
275  hData [1]->SetMaximum(hPred [1]->GetMaximum());
276  } else {
277  hNoOsc[q]->SetMaximum(hNoOsc[1]->GetMaximum());
278  hPred [q]->SetMaximum(hPred [1]->GetMaximum());
279  if (drawData)
280  hData [q]->SetMaximum(hPred [1]->GetMaximum());
281  }
282 
283  // Make graph with error bars for data
284  if (drawData) {
285  gData.push_back( graphAsymmErrorWithBkgScaled(hData[q], hBkg[q]) );
286  gData.back()->GetXaxis()->SetLimits(0., 5.);
287  }
288  }
289 
290 
291  //----------------------------------------------------------------------------
292  // Get Syst Bands
293  //----------------------------------------------------------------------------
294  std::vector<TH1*> ups[nQuantiles+1], dns[nQuantiles+1];
295  std::vector<TH1*> upsNO[nQuantiles+1], dnsNO[nQuantiles+1];
296  if (drawSystBand) {
297  // Initialize the 0th element wuth Q1
298  GetSystBands(calcNO, vPred[0], systs, upsNO[0], dnsNO[0], pot);
299  GetBFSystBands(calc, vPred[0], systs, bfshifts, ups[0], dns[0], pot);
300 
301  // Now do the rest of the quants
302  for (uint q = 0; q < nQuantiles; ++q) {
303  GetSystBands(calcNO, vPred[q], systs, upsNO[q+1], dnsNO[q+1], pot);
304  GetBFSystBands(calc, vPred[q], systs, bfshifts, ups[q+1], dns[q+1], pot);
305 
306  // Tack onto element 0 for all other quants
307  if (q > 0) {
308  for (uint s = 0; s < systs.size(); ++s) {
309  ups [0][s]->Add(ups [q+1][s]);
310  dns [0][s]->Add(dns [q+1][s]);
311  upsNO[0][s]->Add(upsNO[q+1][s]);
312  dnsNO[0][s]->Add(dnsNO[q+1][s]);
313  }
314  }//end adding to AllQuants
315 
316  }// end quant loop
317 
318  // One more set of loops to add in cosmics, could
319  // do this above but it's neater this way
320  for (uint q = 0; q < nQuantiles+1; ++q) {
321  for (uint s = 0; s < systs.size(); ++s) {
322  ups [q][s]->Scale(0.1, "width");
323  dns [q][s]->Scale(0.1, "width");
324  upsNO[q][s]->Scale(0.1, "width");
325  dnsNO[q][s]->Scale(0.1, "width");
326 
327  ups [q][s]->Add(hCosm[q]);
328  dns [q][s]->Add(hCosm[q]);
329  upsNO[q][s]->Add(hCosm[q]);
330  dnsNO[q][s]->Add(hCosm[q]);
331  }
332  }
333  }// endif drawSystBand
334 
335 
336  //----------------------------------------------------------------------------
337  // Plotting Setup
338  //----------------------------------------------------------------------------
339  TFile* fOut = new TFile((sOutDir + "File_" + HORN + sDataPred + sSystBand + ".root").c_str(), "RECREATE");
340 
341  TLatex* xTitle = new TLatex(0.5, 0.03, "Reconstructed neutrino energy (GeV)");
342  xTitle->SetNDC();
343  xTitle->SetTextAlign(22);
344  TLatex* yTitle = new TLatex(0.02, 0.5, "Events / 0.1 GeV");
345  yTitle->SetNDC();
346  yTitle->SetTextAlign(22);
347  yTitle->SetTextAngle(90.);
348 
349  float xSize = xTitle->GetTextSize();
350  float ySize = yTitle->GetTextSize();
351 
352  //----------------------------------------------------------------------------
353  // Plotting Oscillated only
354  //----------------------------------------------------------------------------
355  std::string sTitle = HORN + sDataPred + sSystBand;
356  std::ofstream f((sTitle + ".txt").c_str());
357  std::ofstream f_ByQ((sTitle + "_ByQ.txt").c_str());
358  TCanvas* c = new TCanvas((sTitle).c_str(), (sTitle).c_str(), 960, 800);
359  TCanvas* c_ByQ;
360  TPad *p1, *p2, *p3, *p4;
361  SplitCanvasQuant(c_ByQ, p1, p2, p3, p4);
362  c_ByQ->SetTitle((sTitle + "_ByQ").c_str());
363 
364  // setup the legends with options for drawing data and bands
365  TLegend *leg = new TLegend(0.55, 0.49, 0.78, 0.85);
366  TLegend *leg_l = new TLegend(0.12,0.75,0.27,0.88);
367  TLegend *leg_r = new TLegend(0.55,0.75,0.70,0.88);
368 
369  TLegendEntry *eBand;
370  if (drawData) {
371  leg ->AddEntry(gData[0], "FD data", "LEP");
372  leg_l->AddEntry(gData[0], "FD data", "LEP");
373  }
374 
375  if (drawPred) {
376  leg->AddEntry(hPred[0], "2020 Best-fit", "L");
377  if (drawSystBand)
378  leg_r->AddEntry(hPred[0], "2020 Best-fit", "L");
379  else
380  leg_l->AddEntry(hPred[0], "2020 Best-fit", "L");
381 
382  if (drawSystBand){
383  eBand = leg->AddEntry("error_l","1-#sigma syst. range", "BF");
384  eBand->SetFillStyle(1001);
385  eBand->SetFillColor(kPredErrColor);
386  eBand->SetLineColor(kPredErrColor);
387 
388  eBand = leg_r->AddEntry("error_l","1-#sigma syst. range", "BF");
389  eBand->SetFillStyle(1001);
390  eBand->SetFillColor(kPredErrColor);
391  eBand->SetLineColor(kPredErrColor);
392  }
393  leg->AddEntry(hBkg[0] , "Background", "BF"); leg_l->AddEntry(hBkg[0] , "Background", "BF");
394  // leg->AddEntry(hWS[0] , sWSLabel.c_str(), "BF"); leg_r->AddEntry(hWS[0] , sWSLabel.c_str(), "BF");
395  // leg->AddEntry(hCosm[0], "Cosmic bkgd." , "BF"); leg_r->AddEntry(hCosm[0], "Cosmic bkgd." , "BF");
396  }
397 
398  leg ->SetTextSize(0.05); leg ->SetBorderSize(0); leg ->SetFillStyle (0);
399  leg_l->SetTextSize(0.03); leg_l->SetBorderSize(0); leg_l->SetFillStyle (0);
400  leg_r->SetTextSize(0.03); leg_r->SetBorderSize(0); leg_r->SetFillStyle (0);
401 
402  // Quantile loop to actually plot things
403  for (uint q = 0; q < nQuantiles+1; ++q) {
404  MakeHistCanvasReady_Quant(q, hPred[q], hPred[q]->GetMaximum());
405  if (drawData)
406  MakeHistCanvasReady_Quant(q, gData[q]->GetHistogram(), hData[q]->GetMaximum());
407 
408  if (q == 0) c ->cd();
409  if (q == 1) p1->cd();
410  if (q == 2) p2->cd();
411  if (q == 3) p3->cd();
412  if (q == 4) p4->cd();
413 
414  // Can use this for AllQuants and if ups/dns vectors
415  // are empty then they just aren't drawn
416  if (drawPred) {
417  PlotWithSystErrorBand_Quant(q, hPred[q], ups[q], dns[q],
420  hPred[q]->GetMaximum(), true);
421 
422 
423  // hWS [q]->Draw("hist same");
424  hBkg [q]->Draw("hist same");
425  // hCosm[q]->Draw("hist same");
426 
427  if (drawData) gData[q]->Draw("P");
428 
429  } else if (drawData) {
430  gData[q]->Draw("AP");
431  }
432 
433  DrawQuantLabel(q);
434  }
435 
436  // Finish up and save
437  c->cd();
438  leg->Draw();
439  xTitle->SetTextSize(xSize);
440  yTitle->SetTextSize(ySize);
441  xTitle->Draw();
442  yTitle->Draw();
443  if (isPreliminary) Preliminary();
444  DrawBeamLabel(isFHC);
445  c->Update();
446  if (!drawData) leg->SetY1NDC(0.58);
447  if (!drawSystBand) leg->SetY1NDC(0.58);
448  if (!drawData && !drawSystBand) leg->SetY1NDC(0.67);
449  c->Modified();
450  c->SaveAs((TString)c->GetTitle() + ".png");
451  c->SaveAs((TString)c->GetTitle() + ".pdf");
452  c->SaveAs((TString)c->GetTitle() + ".eps");
453  c->Write();
454  f << HORN
455  << " reconstructed energy spectra.";
456  if (drawPred) f << " The prediction is at the 2020 joint analysis best fit"
457  << " point with the associated systematic pulls applied.";
458  if (drawData) f << " FD data points shown with Poisson errors.";
459  if (drawSystBand) f << " +/- 1 sigma systematic bands shown on the prediction.";
460 
461  f << std::endl;
462  f.close();
463 
464  c_ByQ->cd();
465  if (isPreliminary) Preliminary();
466  DrawBeamLabel(isFHC);
467  xTitle->SetTextSize(0.04);
468  yTitle->SetTextSize(0.04);
469  xTitle->Draw();
470  yTitle->Draw();
471  p1->cd();
472  leg_l->Draw();
473  p2->cd();
474  leg_r->Draw();
475  c_ByQ->Update();
476  c_ByQ->SaveAs((TString)c_ByQ->GetTitle() + ".png");
477  c_ByQ->SaveAs((TString)c_ByQ->GetTitle() + ".pdf");
478  c_ByQ->SaveAs((TString)c_ByQ->GetTitle() + ".eps");
479  c_ByQ->Write();
480  f_ByQ << HORN
481  << " reconstructed energy spectra split by hadronic energy fraction quartile.";
482  if (drawPred) f_ByQ << " The prediction is at the 2020 joint analysis best fit"
483  << " point with the associated systematic pulls applied.";
484  if (drawData) f_ByQ << " FD data points shown with Poisson errors.";
485  if (drawSystBand) f_ByQ << " +/- 1 sigma systematic bands shown on the prediction.";
486 
487 
488  f_ByQ << std::endl;
489  f_ByQ.close();
490 
491  //----------------------------------------------------------------------------
492  // Plotting Unoscillated overlays
493  //----------------------------------------------------------------------------
494  if (makeUnoscPlots) {
495  std::string sTitleNO = HORN + sDataPred + sSystBand + "_NoOsc";
496  std::ofstream fNO((sTitleNO + ".txt").c_str());
497  std::ofstream fNO_ByQ((sTitleNO + "_ByQ.txt").c_str());
498  TCanvas* cNO = new TCanvas((sTitleNO).c_str(), (sTitleNO).c_str(), 960, 800);
499  TCanvas* cNO_ByQ;
500  TPad *p1NO, *p2NO, *p3NO, *p4NO;
501  SplitCanvasQuant(cNO_ByQ, p1NO, p2NO, p3NO, p4NO);
502  cNO_ByQ->SetTitle((sTitleNO + "_ByQ").c_str());
503 
504  // setup the legends
505  TLegend *legNO = new TLegend(0.55, 0.40, 0.78, 0.85);
506  TLegend *legNO_l = new TLegend(0.30, 0.62, 0.43, 0.8);
507  TLegend *legNO_r = new TLegend(0.73, 0.67, 0.86, 0.8);
508  TLegendEntry *eBandNO;
509  if (drawData) {
510  legNO->AddEntry(gData[0], "FD data", "LEP");
511  legNO_l->AddEntry(gData[0], "FD data", "LEP");
512  }
513 
514  legNO->AddEntry(hNoOsc[0], "No oscillation", "L");
515  legNO_l->AddEntry(hNoOsc[0], "No oscillation", "L");
516 
517  if (drawPred) {
518  legNO->AddEntry(hPred [0], "2020 Best-fit" , "L");
519  if (drawSystBand || drawData)
520  legNO_r->AddEntry(hPred [0], "2020 Best-fit" , "L");
521  else
522  legNO_l->AddEntry(hPred [0], "2020 Best-fit" , "L");
523 
524  if (drawSystBand){
525  eBandNO = legNO->AddEntry("errorNO","1-#sigma syst. range", "BF");
526  eBandNO->SetFillStyle(1001);
527  eBandNO->SetFillColor(kNoOscBandColor);
528  eBandNO->SetLineColor(kNoOscBandColor);
529 
530  eBandNO = legNO_l->AddEntry("errorNO_l","1-#sigma syst. range", "BF");
531  eBandNO->SetFillStyle(1001);
532  eBandNO->SetFillColor(kNoOscBandColor);
533  eBandNO->SetLineColor(kNoOscBandColor);
534  }
535  legNO->AddEntry(hBkg[0], "Background","BF" );
536  if (!drawSystBand && !drawData)
537  legNO_l->AddEntry(hBkg[0], "Background","BF" );
538  else
539  legNO_r->AddEntry(hBkg[0], "Background","BF" );
540  }
541 
542  legNO ->SetTextSize(0.05); legNO ->SetBorderSize(0); legNO ->SetFillStyle (0);
543  legNO_l->SetTextSize(0.03); legNO_l->SetBorderSize(0); legNO_l->SetFillStyle (0);
544  legNO_r->SetTextSize(0.03); legNO_r->SetBorderSize(0); legNO_r->SetFillStyle (0);
545 
546  // Quantile loop to actually plot things
547  for (uint q = 0; q < nQuantiles+1; ++q) {
548  MakeHistCanvasReady_Quant(q, hNoOsc[q], hNoOsc[q]->GetMaximum());
549 
550  if (q == 0) cNO ->cd();
551  if (q == 1) p1NO->cd();
552  if (q == 2) p2NO->cd();
553  if (q == 3) p3NO->cd();
554  if (q == 4) p4NO->cd();
555 
556  // Can use this for AllQuants and if ups/dns vectors
557  // are empty then the just aren't drawn
558  PlotWithSystErrorBand_Quant(q, hNoOsc[q], upsNO[q], dnsNO[q],
559  kNoOscColor,
561  hNoOsc[q]->GetMaximum(), true);
562 
563  if (drawPred) {
564  PlotWithSystErrorBand_Quant(q, hPred[q], ups[q], dns[q],
567  hNoOsc[q]->GetMaximum(), false);
568  hBkg[q]->Draw("hist same");
569  }
570 
571  if (drawData) gData[q]->Draw("P");
572 
573  DrawQuantLabel(q);
574  }
575 
576  // Finish up and save
577  cNO->cd();
578  legNO->Draw();
579  xTitle->SetTextSize(xSize);
580  yTitle->SetTextSize(ySize);
581  xTitle->Draw();
582  yTitle->Draw();
583  if (isPreliminary) Preliminary();
584  DrawBeamLabel(isFHC);
585  cNO->Update();
586  cNO->SaveAs((TString)cNO->GetTitle() + ".png");
587  cNO->SaveAs((TString)cNO->GetTitle() + ".pdf");
588  cNO->SaveAs((TString)cNO->GetTitle() + ".eps");
589  cNO->Write();
590  fNO << HORN
591  << " reconstructed energy spectra with the unoscillated prediction.";
592  if (drawPred) fNO << " The oscillated prediction is at the 2020 joint analysis best"
593  << " fit point with the associated systematic pulls applied.";
594  if (drawData) fNO << " FD data points shown with Poisson errors.";
595  if (drawSystBand) fNO << " +/- 1 sigma systematic bands shown on the prediction.";
596 
597 
598  fNO << std::endl;
599  fNO.close();
600 
601  cNO_ByQ->cd();
602  if (isPreliminary) Preliminary();
603  DrawBeamLabel(isFHC);
604  xTitle->SetTextSize(0.04);
605  yTitle->SetTextSize(0.04);
606  xTitle->Draw();
607  yTitle->Draw();
608  p1NO ->cd();
609  legNO_l->Draw();
610  p2NO ->cd();
611  legNO_r->Draw();
612  cNO_ByQ->Update();
613  cNO_ByQ->SaveAs((TString)cNO_ByQ->GetTitle() + ".png");
614  cNO_ByQ->SaveAs((TString)cNO_ByQ->GetTitle() + ".pdf");
615  cNO_ByQ->SaveAs((TString)cNO_ByQ->GetTitle() + ".eps");
616  cNO_ByQ->Write();
617  fNO_ByQ << HORN
618  << " reconstructed energy spectra split by hadronic energy fraction"
619  << " quartile with the unoscillated prediction.";
620  if (drawPred) fNO_ByQ << " The oscillated prediction is at the 2020 joint analysis best"
621  << " fit point with the associated systematic pulls applied.";
622  if (drawData) fNO_ByQ << " FD data points shown with Poisson errors.";
623  if (drawSystBand) fNO_ByQ << " +/- 1 sigma systematic bands shown on the prediction.";
624 
625 
626  fNO_ByQ << std::endl;
627  fNO_ByQ.close();
628 
629  }
630 
631  //----------------------------------------------------------------------------
632  // Plotting Ratio to NoOsc
633  //----------------------------------------------------------------------------
634  if (makeUnoscRatio) {
635  // Remove backgrounds
636  for (uint q = 0; q < nQuantiles+1; ++q) {
637  hNoOsc[q]->Add(hNOBkg[q], -1);
638  if (drawData)
639  hData[q]->Add(hBkg[q], -1);
640  if(drawPred)
641  hPred[q]->Add(hBkg[q], -1);
642 
643  if (drawSystBand)
644  for (uint s = 0; s < systs.size(); ++s) {
645  ups[q][s]->Add(hBkg[q], -1);
646  dns[q][s]->Add(hBkg[q], -1);
647  }
648 
649  }
650 
651  // Make the ratio graph
652  std::vector<TGraph*> gDataRat;
653  if (drawData){
654  for (uint q = 0; q < nQuantiles+1; ++q) {
655  gDataRat.push_back(RatioAsymmErrorScaled(hData[q], hNoOsc[q]));
656  gDataRat.back()->GetHistogram()->SetMaximum(1.4);
657  gDataRat.back()->GetXaxis()->SetLimits(0., 5.);
658  }
659  }
660 
661  std::string sTitleRat = HORN + sDataPred + sSystBand + "_Ratio";
662 
663  std::ofstream fRat((sTitleRat + ".txt").c_str());
664  std::ofstream fRat_ByQ((sTitleRat + "_ByQ.txt").c_str());
665  TCanvas* cRat = new TCanvas((sTitleRat).c_str(), (sTitleRat).c_str(), 960, 800);
666  TCanvas* cRat_ByQ;
667  TPad *p1Rat, *p2Rat, *p3Rat, *p4Rat;
668  SplitCanvasQuant(cRat_ByQ, p1Rat, p2Rat, p3Rat, p4Rat);
669  cRat_ByQ->SetTitle((sTitleRat + "_ByQ").c_str());
670 
671  TLatex* yTitleRat = new TLatex(0.02, 0.5, "Ratio to no oscillation");
672  yTitleRat->SetNDC();
673  yTitleRat->SetTextAlign(22);
674  yTitleRat->SetTextAngle(90.);
675 
676  TLegend *legRat = new TLegend(0.43, 0.73, 0.66, 0.90);
677  TLegend *legRat_l = new TLegend(0.12,0.75,0.27,0.88);
678 
679  TLine *unitline = new TLine(0, 1, 5, 1);
680  unitline->SetLineColor(kGray+1);
681  unitline->SetLineStyle(9);
682  unitline->SetLineWidth(4);
683 
684  if (drawData) {
685  legRat ->AddEntry(gDataRat[0], "FD data", "LEP");
686  legRat_l->AddEntry(gDataRat[0], "FD data", "LEP");
687  }
688  if (drawPred) {
689  legRat ->AddEntry(hPred[0], "2020 Best-fit", "L");
690  legRat_l->AddEntry(hPred[0], "2020 Best-fit", "L");
691 
692  if (drawSystBand) {
693  eBand = legRat->AddEntry("error_l","1-#sigma syst. range", "BF");
694  eBand->SetFillStyle(1001);
695  eBand->SetFillColor(kPredErrColor);
696  eBand->SetLineColor(kPredErrColor);
697 
698  eBand = legRat_l->AddEntry("error_l","1-#sigma syst. range", "BF");
699  eBand->SetFillStyle(1001);
700  eBand->SetFillColor(kPredErrColor);
701  eBand->SetLineColor(kPredErrColor);
702  }
703  }
704  legRat ->SetTextSize(0.05); legRat ->SetBorderSize(0); legRat ->SetFillStyle (0);
705  legRat_l->SetTextSize(0.03); legRat_l->SetBorderSize(0); legRat_l->SetFillStyle (0);
706 
707  for (uint q = 0; q < nQuantiles+1; ++q) {
708  if (q == 0) cRat ->cd();
709  if (q == 1) p1Rat->cd();
710  if (q == 2) p2Rat->cd();
711  if (q == 3) p3Rat->cd();
712  if (q == 4) p4Rat->cd();
713 
714  // Get a graph to start with
715 
716  if (drawPred) {
717  if (drawSystBand) {
718  TGraphAsymmErrors* gPred =
719  PlotWithSystErrorBand_Quant(q, hPred[q], ups[q], dns[q],
722  hPred[q]->GetMaximum(), true);
723 
724  for (int bin = 1; bin <= 19; ++bin) {
725 
726  double x, y;
727  gPred->GetPoint(bin, x, y);
728  double den = hNoOsc[q]->GetBinContent(bin);
729 
730  double errup = 0, errdn = 0;
731  errup = gPred->GetErrorYhigh(bin) / den;
732  errdn = gPred->GetErrorYlow(bin) / den;
733 
734  gPred->SetPoint(bin, x, y / den);
735  gPred->SetPointEYhigh(bin, errup);
736  gPred->SetPointEYlow(bin, errdn);
737  }
738 
739  gPred->GetYaxis()->SetTitle("");
740  gPred->GetXaxis()->SetLimits(0., 5.);
741  gPred->GetHistogram()->SetMaximum(1.4);
742  gPred->SetFillColor(kPredErrColor);
743  gPred->Draw("A2");
744  }
745 
746  hPred[q]->Divide((TH1*)hNoOsc[q]);
747  MakeHistCanvasReady_Quant(q, hPred[q], 1.4);
748  hPred[q]->GetYaxis()->SetTitle("Ratio to no oscillation");
749  hPred[q]->Draw("axis same ][");
750  hPred[q]->Draw("hist same ][");
751 
752  if (drawData) gDataRat[q]->Draw("P");
753  } else if (drawData) gDataRat[q]->Draw("AP");
754 
755  unitline->Draw();
756 
757  DrawQuantLabel(q);
758  }
759 
760  cRat->cd();
761  legRat->Draw();
762  xTitle->SetTextSize(xSize);
763  xTitle->Draw();
764  yTitleRat->Draw();
765  if (isPreliminary) Preliminary();
766  DrawBeamLabel(isFHC);
767  cRat->Update();
768  cRat->SaveAs((TString)cRat->GetTitle() + ".png");
769  cRat->SaveAs((TString)cRat->GetTitle() + ".pdf");
770  cRat->SaveAs((TString)cRat->GetTitle() + ".eps");
771  cRat->Write();
772  fRat << HORN
773  << " reconstructed energy ratio spectra with the unoscillated"
774  << " prediction as the denominator. Backgrounds are removed from"
775  << " the numerator and denomenator before the ratio is taken.";
776  if (drawPred) fRat << " The prediction is at the 2020 joint analysis best fit point with"
777  << " the associated systematic pulls applied.";
778  if (drawData) fRat << " Ratio with FD data shown with Poisson errors evaluated after backgrounds are removed.";
779  if (drawSystBand) fRat << " +/- 1 sigma systematic bands shown on the prediction. The bands"
780  << " only consider the systematic shifts at the BF oscillation point."
781  << " The oscillation point is held constant and the ratio is the hi/lo"
782  << " syst band of the oscillated spectrum divided by the unoscillated"
783  << " spectrum. No systematic errors are considered for the unsosc pred.";
784 
785 
786  fRat << std::endl;
787  fRat.close();
788 
789  cRat_ByQ->cd();
790  if (isPreliminary) Preliminary();
791  DrawBeamLabel(isFHC);
792  xTitle->SetTextSize(0.04);
793  xTitle->Draw();
794  yTitleRat->SetTextSize(0.04);
795  yTitleRat->Draw();
796  p1Rat->cd();
797  legRat_l->Draw();
798  cRat_ByQ->Update();
799  cRat_ByQ->SaveAs((TString)cRat_ByQ->GetTitle() + ".png");
800  cRat_ByQ->SaveAs((TString)cRat_ByQ->GetTitle() + ".pdf");
801  cRat_ByQ->SaveAs((TString)cRat_ByQ->GetTitle() + ".eps");
802  cRat_ByQ->Write();
803  fRat_ByQ << HORN
804  << " reconstructed energy ratio spectra split by hadronic energy"
805  << " fraction quartiles with the unoscillated prediction as the denominator."
806  << " Backgrounds are removed from the numerator and denomenator before the ratio is taken.";
807  if (drawPred) fRat_ByQ << " The prediction is at the 2020 joint analysis best fit point with"
808  << " the associated systematic pulls applied.";
809  if (drawData) fRat_ByQ << " Ratio with FD data shown with Poisson errors.";
810  if (drawSystBand) fRat_ByQ << " +/- 1 sigma systematic bands shown on the prediction. The bands"
811  << " only consider the systematic shifts at the BF oscillation point."
812  << " The oscillation point is held constant and the ratio is the hi/lo"
813  << " syst band of the oscillated spectrum divided by the unoscillated"
814  << " spectrum. No systematic errors are considered for the unsosc pred.";
815 
816  fRat_ByQ << std::endl;
817  fRat_ByQ.close();
818 
819  }
820 
821  fOut->Close();
822 }
Pass neutrinos through unchanged.
Definition: IOscCalc.h:40
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
double lt
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
Loads shifted spectra from files.
const Color_t kCosmicBackgroundColor
Definition: Style.h:41
TH1D * ToTH1(double exposure, Color_t col=kBlack, Style_t style=kSolid, EExposureType expotype=kPOT, EBinType bintype=kBinContent) const
Histogram made from this Spectrum, scaled to some exposure.
Definition: Spectrum.cxx:148
const Color_t kNoOscBandColor
Definition: Style.h:35
Simple record of shifts applied to systematic parameters.
Definition: SystShifts.h:20
const std::map< std::string, PredSel > predSels
Float_t den
Definition: plot.C:36
void plot_recoE_numu(std::string horn="fhc", bool drawData=true, bool drawPred=true, bool drawSystBand=true, bool useBF=true, bool makeUnoscPlots=true, bool makeUnoscRatio=true, bool isPreliminary=true)
TPaveText * DrawQuantLabel(int quant)
Definition: Plots.cxx:1574
const Color_t kNoOscColor
Definition: Style.h:34
osc::OscCalcDumb calc
osc::IOscCalcAdjustable * DefaultOscCalc()
Create a new calculator with default assumptions for all parameters.
Definition: Calcs.cxx:49
void GetSystBands(osc::IOscCalc *calc, const IPrediction *pred, std::vector< const ISyst * > systs, std::vector< TH1 * > &hUps, std::vector< TH1 * > &hDns, double pot, Flavors::Flavors_t flav, Current::Current_t curr, Sign::Sign_t sign)
Definition: Plots.cxx:244
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
const XML_Char const XML_Char * data
Definition: expat.h:268
const Color_t kFullPredColor
Definition: Style.h:36
const XML_Char * s
Definition: expat.h:262
static std::unique_ptr< SystShifts > LoadFrom(TDirectory *dir, const std::string &name)
Definition: SystShifts.cxx:256
static std::unique_ptr< Spectrum > LoadFrom(TDirectory *dir, const std::string &name)
Definition: Spectrum.cxx:546
TPaveText * DrawBeamLabel(bool isFHC)
Put the standardized beam label in the left corner of the active canvas.
Definition: Plots.cxx:1555
const Color_t kPredErrColor
Definition: Style.h:37
virtual _IOscCalc< T > * Copy() const override
Definition: IOscCalc.h:44
#define pot
const double kAna2020FHCLivetime
Definition: Exposures.h:237
const Color_t kTotalBkgColor
Definition: Style.h:39
const double kAna2020FHCPOT
Definition: Exposures.h:233
float bin[41]
Definition: plottest35.C:14
double POT() const
Definition: Spectrum.h:219
const double kAna2020RHCPOT
Definition: Exposures.h:235
static bool isFHC
void Preliminary()
Neutrinos-only.
Definition: IPrediction.h:49
OStream cout
Definition: OStream.cxx:6
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
TGraph * graphAsymmErrorWithBkgScaled(TH1 *data, TH1 *bkgd, double overallScale)
Definition: Plots.cxx:1020
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
const double kAna2020RHCLivetime
Definition: Exposures.h:238
void GetBFSystBands(osc::IOscCalc *calc, const IPrediction *pred, std::vector< const ISyst * > systs, const SystShifts &bfshifts, std::vector< TH1 * > &hUps, std::vector< TH1 * > &hDns, double pot, Flavors::Flavors_t flav, Current::Current_t curr, Sign::Sign_t sign)
Definition: Plots.cxx:271
const Color_t kNumuWSColor
Definition: Style.h:38
std::vector< const ISyst * > get3FlavorAna2020AllSysts(const EAnaType2020 ana, const bool smallgenie, const BeamType2020 beam, const bool isFit, const bool ptExtrap)
void SplitCanvasQuant(TCanvas *&canvas, TPad *&pad1, TPad *&pad2, TPad *&pad3, TPad *&pad4)
Definition: Plots.cxx:1437
TGraph * RatioAsymmErrorScaled(TH1 *data, TH1 *pred, double overallScale)
Definition: Plots.cxx:1105
void Beam(bool isRHC)
std::string to_string(ModuleType mt)
Definition: ModuleType.h:32
void SetShift(const ISyst *syst, double shift, bool force=false)
Definition: SystShifts.cxx:80
void FillWithDimColor(TH1 *h, bool usealpha, float dim)
unsigned int uint