draw_michel_plots.C
Go to the documentation of this file.
1 // Macro to draw 2018 nue FHC analysis MichelDecomp plot
2 
3 #include "OscLib/OscCalcDumb.h"
4 
5 #include "CAFAna/Vars/Vars.h"
8 
9 #include "CAFAna/Cuts/Cuts.h"
10 #include "CAFAna/Cuts/SpillCuts.h"
12 #include "CAFAna/Cuts/TruthCuts.h"
13 #include "CAFAna/Cuts/NueCutsSecondAna.h"
14 #include "CAFAna/Core/Spectrum.h"
17 
18 #include "CAFAna/Analysis/Style.h"
19 #include "CAFAna/Analysis/SALoaders.h"
29 
30 #include "CAFAna/Systs/Systs.h"
31 #include "CAFAna/Core/SystShifts.h"
32 
33 #include "Utilities/func/MathUtil.h"
34 
35 #include "TROOT.h"
36 #include "TCanvas.h"
37 #include "TH2.h"
38 #include "TFile.h"
39 #include "TPad.h"
40 #include "TLegend.h"
41 #include "TProfile.h"
42 #include "TRandom3.h"
43 #include "THStack.h"
44 #include "TGaxis.h"
45 #include "TLine.h"
46 #include "TStyle.h"
47 #include "TLatex.h"
48 #include "TGraph.h"
49 
50 #include "TColor.h"
51 
52 
53 #include <cmath>
54 #include <vector>
55 #include <string>
56 #include <iostream>
57 #include <sstream>
58 #include <iomanip>
59 
60 
61 using namespace ana;
62 
63 void FormatRatio(TH1 *hrat, std::string xtitle, std::string ytitle)
64 {
65  gPad->SetTopMargin(0);
66  gPad->SetBottomMargin(0.2);
67  hrat->GetXaxis()->SetTitleSize(0.10);
68  hrat->GetXaxis()->SetLabelSize(0.08);
69  hrat->GetXaxis()->SetTitleOffset(0.8);
70  hrat->GetYaxis()->SetTitleSize(0.10);
71  hrat->GetYaxis()->SetLabelSize(0.08);
72  hrat->GetYaxis()->SetTitleOffset(0.5);
73  hrat->SetTitle("");
74  hrat->GetYaxis()->SetRangeUser(0.6501,1.3499);
75  hrat->GetXaxis()->SetTitle(xtitle.c_str());
76  hrat->GetYaxis()->SetTitle(ytitle.c_str());
77  hrat->GetYaxis()->CenterTitle();
78  hrat->GetXaxis()->CenterTitle();
79 }
80 
81 
82 void MakeStackNME(std::vector<const TH2F*> hs, double potrat, int idx1, int idx2,
83  std::vector<const TH1F*> scales,
85  std::string yrattitle, double &tot, bool applyScales = false)
86 {
87  TCanvas *can = new TCanvas(UniqueName().c_str(),"",1500,700);
88  can->Divide(3,1);
89 
90 
91  TH1D *hnm = (TH1D*)(hs[0]->ProjectionY(UniqueName().c_str(),idx1,idx2))->Clone(UniqueName().c_str());
92  TH1D *hnc = (TH1D*)(hs[1]->ProjectionY(UniqueName().c_str(),idx1,idx2))->Clone(UniqueName().c_str());
93  TH1D *hne = (TH1D*)(hs[2]->ProjectionY(UniqueName().c_str(),idx1,idx2))->Clone(UniqueName().c_str());
94  TH1D *hda = (TH1D*)(hs[3]->ProjectionY(UniqueName().c_str(),idx1,idx2))->Clone(UniqueName().c_str());
95 
96  TH1D *g0 = (TH1D*)(hs[0]->ProjectionY(UniqueName().c_str(),idx1,idx2))->Clone(UniqueName().c_str());
97  TH1D *g1 = (TH1D*)(hs[1]->ProjectionY(UniqueName().c_str(),idx1,idx2))->Clone(UniqueName().c_str());
98  TH1D *g2 = (TH1D*)(hs[2]->ProjectionY(UniqueName().c_str(),idx1,idx2))->Clone(UniqueName().c_str());
99 
100  TH1D *f0 = (TH1D*)(hs[0]->ProjectionY(UniqueName().c_str(),idx1,idx2))->Clone(UniqueName().c_str());
101  TH1D *f1 = (TH1D*)(hs[1]->ProjectionY(UniqueName().c_str(),idx1,idx2))->Clone(UniqueName().c_str());
102  TH1D *f2 = (TH1D*)(hs[2]->ProjectionY(UniqueName().c_str(),idx1,idx2))->Clone(UniqueName().c_str());
103 
104 
105  hnm->Scale(potrat*1e-3);
106  hnc->Scale(potrat*1e-3);
107  hne->Scale(potrat*1e-3);
108 
109  g0->Scale(potrat*1e-3);
110  g1->Scale(potrat*1e-3);
111  g2->Scale(potrat*1e-3);
112 
113  f0->Scale(potrat*1e-3);
114  f1->Scale(potrat*1e-3);
115  f2->Scale(potrat*1e-3);
116 
117 
118  for (int i = 1; i <= hda->GetNbinsX(); i++){
119  hda->SetBinContent(i,0.001*hda->GetBinContent(i));
120  hda->SetBinError (i,0.001*hda->GetBinError(i));
121  }
122 
123 
124  TH1D *hmc = (TH1D*)(hs[0]->ProjectionY(UniqueName().c_str(),idx1,idx2))->Clone(UniqueName().c_str());
125  hmc->Scale(potrat*1e-3);
126  hmc->Add(hnc);
127  hmc->Add(hne);
128 
129  TH1D *fmc = (TH1D*)(hs[0]->ProjectionY(UniqueName().c_str(),idx1,idx2))->Clone(UniqueName().c_str());
130  fmc->Scale(potrat*1e-3);
131 
132 
133  hnm->SetFillColor(kGreen+2);
134  hnc->SetFillColor(kAzure);
135  hne->SetFillColor(kPink+9);
136  hda->SetMinimum(1e-3);
137 
138  int EIdx = idx1%10;
139  std::string ERange = "1.0 < E_{#nu} < 1.5 GeV";
140  if (EIdx == 1) ERange = "0.0 < E_{#nu} < 0.5 GeV";
141  if (EIdx == 2) ERange = "0.5 < E_{#nu} < 1.0 GeV";
142  if (EIdx == 3) ERange = "1.0 < E_{#nu} < 1.5 GeV";
143  if (EIdx == 4) ERange = "1.5 < E_{#nu} < 2.0 GeV";
144  if (EIdx == 5) ERange = "2.0 < E_{#nu} < 2.5 GeV";
145  if (EIdx == 6) ERange = "2.5 < E_{#nu} < 3.0 GeV";
146  if (EIdx == 7) ERange = "3.0 < E_{#nu} < 3.5 GeV";
147  if (EIdx == 8) ERange = "3.5 < E_{#nu} < 4.0 GeV";
148  if (EIdx == 9) ERange = "4.0 < E_{#nu} < 4.5 GeV";
149 
150  ERange += " and ";
151 
152  int CIdx = idx1/10;
153  std::string CRange = "0.84 <= CVN < 0.96";
154  if (CIdx == 1) CRange = "CVN >= 0.96";
155 
156  hda->SetTitle((ERange + CRange).c_str());
157  hmc->GetYaxis()->SetTitleSize(0.04);
158  hmc->GetYaxis()->SetTitleOffset(1.4);
159  hda->GetYaxis()->SetTitleSize(0.04);
160  hda->GetYaxis()->SetTitleOffset(1.4);
161 
162 
163  THStack *st = new THStack();
164  st->Add(hne);
165  st->Add(hnc);
166  st->Add(hnm);
167 
168  can->cd(1);
169 
170 
171  double max = std::max(hmc->GetMaximum(),hda->GetMaximum());
172  hmc->GetYaxis()->SetRangeUser(1e-3,1.2*max);
173  hda->GetYaxis()->SetRangeUser(1e-3,1.2*max);
174 
175  hmc->GetYaxis()->SetTitle("10^{3}Events / 8#times10^{20}");
176  hmc->GetYaxis()->CenterTitle();
177  hda->GetYaxis()->SetTitle("10^{3}Events / 8#times10^{20}");
178  hda->GetYaxis()->CenterTitle();
179 
180  g0->Divide(hmc);
181  g1->Divide(hmc);
182  g2->Divide(hmc);
183 
184  g0->SetFillColor(kGreen+2);
185  g1->SetFillColor(kAzure);
186  g2->SetFillColor(kPink+9);
187 
188  THStack *stg = new THStack();
189  stg->Add(g2);
190  stg->Add(g1);
191  stg->Add(g0);
192 
193  g0->GetXaxis()->SetTitle("N_{ME}");
194  g0->GetXaxis()->CenterTitle();
195  g0->GetYaxis()->SetTitle("Fraction of Events");
196  g0->GetYaxis()->CenterTitle();
197  g0->GetYaxis()->SetTitleOffset(1.2);
198  gPad->SetLeftMargin(0.13);
199 
200  g0->GetYaxis()->SetRangeUser(0,1);
201  g0->Draw("hist");
202  stg->Draw("hist,same");
203 
204  can->cd(2);
205 
206  TPad *pRat2 = new TPad(UniqueName().c_str(),"",0,0,1,0.32);
207  TPad *pSpec2= new TPad(UniqueName().c_str(),"",0,0.32,1,1);
208 
209  pRat2->Draw();
210  pSpec2->Draw();
211 
212  pSpec2->cd();
213  gPad->SetBottomMargin(0);
214 
215  hda->GetXaxis()->SetTitle("");
216 
217  hda->Draw("E");
218  st->Draw("hist,same");
219  hda->Draw("E,same");
220 
221  TLegend *leg = new TLegend(0.6,0.55,0.8,0.8);
222  leg->AddEntry(hda,"Data","el");
223  leg->AddEntry(hnm,"#nu_{#mu} CC","f");
224  leg->AddEntry(hnc,"NC","f");
225  leg->AddEntry(hne,"#nu_{e} CC","f");
226  leg->Draw();
227 
228  pRat2->cd();
229 
230  TH1F *hrat = (TH1F*)hda->Clone(UniqueName().c_str());
231  hrat->Divide(hmc);
232 
233  FormatRatio(hrat,"N_{ME}","Ratio");
234 
235  for (int i = 1; i <=3; i++)
236  hrat->SetBinError(i,1./sqrt(1000*hda->GetBinContent(i)));
237  //hrat->SetBinError(i,1./sqrt(hda->GetBinContent(i)));
238  hrat->Draw("E");
239 
240  TLine *l = new TLine(0,1,3,1);
241  l->SetLineStyle(7);
242  l->Draw("same");
243 
244  can->cd(3);
245 
246  TPad *pRat3 = new TPad(UniqueName().c_str(),"",0,0,1,0.32);
247  TPad *pSpec3= new TPad(UniqueName().c_str(),"",0,0.32,1,1);
248 
249  pRat3->Draw();
250  pSpec3->Draw();
251 
252  pSpec3->cd();
253  gPad->SetBottomMargin(0);
254 
255  f0->SetFillColor(kGreen+2);
256  f1->SetFillColor(kAzure);
257  f2->SetFillColor(kPink+9);
258 
259  f0->Scale(scales[0]->GetBinContent(idx1));
260  f1->Scale(scales[1]->GetBinContent(idx1));
261  f2->Scale(scales[2]->GetBinContent(idx1));
262 
263  fmc->Scale(scales[0]->GetBinContent(idx1));
264  fmc->Add(f1);
265  fmc->Add(f2);
266 
267  THStack *st3 = new THStack();
268  st3->Add(f2);
269  st3->Add(f1);
270  st3->Add(f0);
271 
272  TH1F *hda2 = (TH1F*)hda->Clone(UniqueName().c_str());
273  hda2->SetTitle("");
274 
275  hda2->Draw("E");
276  st3->Draw("hist,same");
277  hda2->Draw("E,same");
278 
279  TH1F *frat = (TH1F*)hda2->Clone(UniqueName().c_str());
280  frat->Divide(fmc);
281 
282  double chisq = 0;
283  for (int i = 1; i <=3; i++){
284  frat->SetBinError(i,1./sqrt(1000*hda->GetBinContent(i)));
285  chisq += util::sqr(fmc->GetBinContent(i)-hda->GetBinContent(i))/
286  hda->GetBinContent(i);
287  }
288  chisq *= 1000;
289  tot += chisq;
290  std::cout << chisq << std::endl;
291 
292  TLatex latex;
293  latex.SetTextSize(0.09);
294 
295  std::ostringstream out;
296  out << std::setprecision(3) << chisq;
297 
298  latex.DrawLatexNDC(0.46,0.65,("#chi^{2} = "+out.str()).c_str());
299 
300  pRat3->cd();
301 
302  FormatRatio(frat,"N_{ME}","Ratio");
303 
304  frat->Draw("E");
305  l->Draw("same");
306 
307  can->SaveAs(("Fig_AnaBin_"+std::to_string(idx1)+".png").c_str());
308 }
309 
310 void MakeFracPlot(std::vector<const TH1F*> hs)
311 {
312  TH1F *f0 = (TH1F*)hs[0]->Clone(UniqueName().c_str());
313  TH1F *f1 = (TH1F*)hs[1]->Clone(UniqueName().c_str());
314  TH1F *f2 = (TH1F*)hs[2]->Clone(UniqueName().c_str());
315  TH1F *sum= (TH1F*)hs[2]->Clone(UniqueName().c_str());
316  sum->Add(f1);
317  sum->Add(f0);
318 
319  f0->Divide(sum);
320  f1->Divide(sum);
321  f2->Divide(sum);
322 
323  f0->SetFillColor(kGreen+2);
324  f1->SetFillColor(kAzure);
325  f2->SetFillColor(kPink+9);
326 
327 
328  TCanvas *can = new TCanvas();
329 
330  TH2F *h = new TH2F(UniqueName().c_str(),";NueEnergy_2018/CVN;Fraction of Events",1,0,9.5,1,0,1);
331 
332  h->Draw();
333 
334  THStack *st = new THStack();
335  st->Add(f2);
336  st->Add(f1);
337  st->Add(f0);
338  st->Draw("hist same");
339 
340  Nue2018ThreeBinAxis(h,true,true,true);
341 
342  TLine *l = new TLine(0,0.8,19,0.8);
343  l->Draw("same");
344 
345  can->SaveAs("Fracs.png");
346 }
347 
348 void MakeRawPlot(std::vector<const TH1F*> hs, double potrat)
349 {
350  TH1F *hrat = (TH1F*)hs[3]->Clone(UniqueName().c_str());
351  TH1F *hmc = (TH1F*)hs[0]->Clone(UniqueName().c_str());
352  TH1F *hda = (TH1F*)hs[3]->Clone(UniqueName().c_str());
353 
354  TH1F *hs0 = (TH1F*)hs[0]->Clone(UniqueName().c_str());
355  TH1F *hs1 = (TH1F*)hs[1]->Clone(UniqueName().c_str());
356  TH1F *hs2 = (TH1F*)hs[2]->Clone(UniqueName().c_str());
357 
358  hmc->Add(hs1);
359  hmc->Add(hs2);
360  hmc->Scale(potrat*1e-3);
361  hrat->Scale(1e-3);
362  hrat->Divide(hmc);
363 
364  hs0->Scale(potrat*1e-3);
365  hs1->Scale(potrat*1e-3);
366  hs2->Scale(potrat*1e-3);
367 
368  hs0->SetFillColor(kGreen+2);
369  hs1->SetFillColor(kAzure);
370  hs2->SetFillColor(kPink+9);
371 
372  TCanvas *can = new TCanvas();
373 
374  TPad *pRat1 = new TPad("rat1","",0,0,1,0.35);
375  TPad *pSpec1= new TPad("spec1","",0,0.35,1,1);
376 
377  pRat1->Draw();
378  pSpec1->Draw();
379 
380  pSpec1->cd();
381  gPad->SetBottomMargin(0);
382 
383  TH2F *h = new TH2F(UniqueName().c_str(),";NueEnergy2018/CVN;10^{3} Events",1,0,9.5,20000,1e-3,4.8);
384  h->GetYaxis()->CenterTitle();
385  h->GetXaxis()->CenterTitle();
386  h->SetMinimum(1e-3);
387  h->Draw();
388 
389  Nue2018ThreeBinAxis(h,true,true,true);
390  Nue2018ThreeBinLabels(0.95,0.08,kGray+3,true);
391 
392  THStack *st = new THStack();
393  st->Add(hs2);
394  st->Add(hs1);
395  st->Add(hs0);
396 
397  st->Draw("hist,same");
398  hda->SetMarkerStyle(8);
399  hda->Scale(1e-3);
400  hda->Draw("E,same");
401 
402  TLegend *leg = new TLegend(0.36,0.47,0.6,0.77);
403  leg->AddEntry(hs2,"Beam #nu_{e} CC","f");
404  leg->AddEntry(hs0,"#nu_{#mu} CC","f");
405  leg->AddEntry(hs1,"NC","f");
406  leg->AddEntry(hda,"Data");
407  leg->Draw();
408 
409  pRat1->cd();
410  gPad->SetTopMargin(0);
411  gPad->SetBottomMargin(0.2);
412  hrat->GetXaxis()->SetTitleSize(0.10);
413  hrat->GetXaxis()->SetLabelSize(0.08);
414  hrat->GetXaxis()->SetTitleOffset(1);
415  hrat->GetXaxis()->SetRangeUser(0,20);
416  hrat->GetYaxis()->SetTitleSize(0.10);
417  hrat->GetYaxis()->SetLabelSize(0.08);
418  hrat->GetYaxis()->SetTitleOffset(0.5);
419  hrat->SetLineColor(kBlack);
420  hrat->SetMaximum(1.25001);
421  hrat->SetMinimum(0.74999);
422  hrat->GetYaxis()->SetTitle("Data/MC");
423  hrat->GetXaxis()->SetTitle("NueEnergy2018/CVN");
424  hrat->GetYaxis()->CenterTitle();
425  hrat->GetXaxis()->CenterTitle();
426  hrat->Draw("E");
427 
428  Nue2018ThreeBinAxis(hrat,true,true,true);
429  std::cout << "Number of events Uncorr." << endl;
430  std::cout << "Data: " << hda->Integral() << endl;
431  std::cout << "Numu: " << hs0->Integral() << endl;
432  std::cout << "NC: " << hs1->Integral() << endl;
433  std::cout << "Beam: " << hs2->Integral() << endl;
434 
435  TLine *l = new TLine(0,1,19,1);
436  l->SetLineStyle(7);
437  l->Draw("same");
438 
439  can->SaveAs("RawDaMC.png");
440  can->SaveAs("RawDaMC.pdf");
441 }
442 
443 void MakeScalesPlot(std::vector<const TH1F*> hs)
444 {
445  TH1F *hs0 = (TH1F*)hs[0]->Clone(UniqueName().c_str());
446  TH1F *hs1 = (TH1F*)hs[1]->Clone(UniqueName().c_str());
447  TH1F *hs2 = (TH1F*)hs[2]->Clone(UniqueName().c_str());
448 
449  hs0->SetLineColor(kGreen+2);
450  hs1->SetLineColor(kAzure);
451  hs2->SetLineColor(kPink+9);
452 
453  hs0->GetYaxis()->SetRangeUser(0,2);
454  hs0->GetXaxis()->SetRangeUser(0,20);
455  hs0->GetXaxis()->CenterTitle();
456  hs0->GetXaxis()->SetTitle("NueEnergy2018/CVN");
457  hs0->GetYaxis()->SetTitle("Decomp Scale Factor");
458  hs0->GetYaxis()->CenterTitle();
459 
460  TCanvas *can = new TCanvas();
461  hs0->Draw();
462  hs1->Draw("same");
463  hs2->Draw("same");
464 
465  TLegend *leg = new TLegend(0.55,0.64,0.66,0.87);
466  leg->AddEntry(hs0,"#nu_{#mu} CC","le");
467  leg->AddEntry(hs1,"NC","le");
468  leg->AddEntry(hs2,"#nu_{e} CC","le");
469  leg->Draw();
470 
471  TLine *l1 = new TLine(1,0,1,2);
472  l1->Draw("same");
473  TLine *l2 = new TLine(6,0,6,2);
474  l2->Draw("same");
475 
476  can->SaveAs("DecompScales.png");
477 }
478 
479 void MakeResultPlot(MichelDecomp *m, TH1 *hda, std::vector<const TH1F*> ht, std::vector<const TH1F*> hs, double pot, double potrat)
480 {
481  TCanvas *can = new TCanvas();
482 
483  TH1F *ht0 = (TH1F*)ht[0]->Clone(UniqueName().c_str());
484  TH1F *ht1 = (TH1F*)ht[1]->Clone(UniqueName().c_str());
485  TH1F *ht2 = (TH1F*)ht[2]->Clone(UniqueName().c_str());
486 
487  m->NumuComponent();
488  TH1* estNM=(m->NumuComponent() + m->AntiNumuComponent()).ToTH1(pot);
489  TH1* estNE=(m->NueComponent() + m->AntiNueComponent()).ToTH1(pot);
490  TH1* estNC= m->NCTotalComponent().ToTH1(pot);
491 
492  TH1F *ratNM = (TH1F*)hs[0]->Clone(UniqueName().c_str());
493  TH1F *ratNC = (TH1F*)hs[1]->Clone(UniqueName().c_str());
494  TH1F *ratNE = (TH1F*)hs[2]->Clone(UniqueName().c_str());
495 
496  TH1* mcNM=(m->MC_NumuComponent()+m->MC_AntiNumuComponent()).ToTH1(pot);
497  TH1* mcNE=(m->MC_NueComponent() +m->MC_AntiNueComponent()).ToTH1(pot);
498  TH1* mcNC= m->MC_NCTotalComponent().ToTH1(pot);
499 
500  for (int i = 1; i <= hda->GetNbinsX(); i++){
501  hda->SetBinContent(i,1e-3*hda->GetBinContent(i));
502  hda->SetBinError(i,1e-3*hda->GetBinError(i));
503  }
504 
505  mcNM->Scale(1e-3);
506  mcNE->Scale(1e-3);
507  mcNC->Scale(1e-3);
508  estNM->Scale(1e-3);
509  estNE->Scale(1e-3);
510  estNC->Scale(1e-3);
511 
512  TH1* nn2 = (TH1*)estNE->Clone();
513 
514  TH1* nn1 = (TH1*)estNE->Clone();
515  nn1->Add(estNC);
516 
517  TH1* nn0 = (TH1*)estNE->Clone();
518  nn0->Add(estNC);
519  nn0->Add(estNM);
520 
521  ht0->Scale(potrat*1e-3);
522  ht1->Scale(potrat*1e-3);
523  ht2->Scale(potrat*1e-3);
524 
525  estNM->SetFillColor(kGreen-9);
526  estNC->SetFillColor(kAzure-9);
527  estNE->SetFillColor(kPink+1);
528  ratNM->SetLineColor(kGreen-9);
529  ratNC->SetLineColor(kAzure-9);
530  ratNE->SetLineColor(kPink+1);
531  ratNM->SetMarkerColor(kGreen-9);
532  ratNC->SetMarkerColor(kAzure-9);
533  ratNE->SetMarkerColor(kPink+1);
534 
535  ht0->SetLineStyle(2);
536  ht1->SetLineStyle(2);
537  ht2->SetLineStyle(2);
538 
539  ht1->Add(ht2);
540 
541  ht0->Add(ht1);
542 
543  TH1* htdum=(TH1*)ht0->Clone();
544 
545  ht0->SetLineColor(kGreen+4);
546  ht1->SetLineColor(kAzure+3);
547  ht2->SetLineColor(kPink+4);
548 
549  THStack *st = new THStack();
550  st->Add(estNE);
551  st->Add(estNC);
552  st->Add(estNM);
553 
554 
555  hda->SetMarkerStyle(8);
556  hda->GetXaxis()->SetTitle("");
557  hda->GetXaxis()->CenterTitle();
558  hda->GetYaxis()->SetTitle("10^{3} Events");
559  hda->GetYaxis()->CenterTitle();
560  hda->SetTitle("");
561 
562  hda->SetMaximum(4);
563  hda->SetMinimum(1e-3);
564 
565  hda->GetXaxis()->SetRangeUser(0,20);
566  hda->Draw("E");
567  st ->Draw("hist,same");
568  ht0->Draw("hist,same");
569  ht1->Draw("hist,same");
570  ht2->Draw("hist,same");
571 
572  nn0->Draw("hist,same");
573  nn1->Draw("hist,same");
574  nn2->Draw("hist,same");
575  hda->Draw("E,same");
576 
577  Nue2018ThreeBinAxis(hda,true,true,true);
578  Nue2018ThreeBinLabels(0.95,0.08,kGray+3,true);
579 
580  TLegend *leg = new TLegend(0.3,0.47,0.6,0.77);
581  leg->AddEntry(estNE,"Beam #nu_{e} CC","f");
582  leg->AddEntry(estNC,"NC","f");
583  leg->AddEntry(estNM,"#nu_{#mu} CC","f");
584  leg->AddEntry(hda, "Data","le");
585  leg->AddEntry(htdum, "Uncorr.","l");
586  leg->Draw();
587 
588  std::cout << "Number of events after BEN+MICHEL" << endl;
589  std::cout << "Data: " << hda->Integral() << endl;
590  std::cout << "Numu: " << estNM->Integral() << endl;
591  std::cout << "NC: " << estNC->Integral() << endl;
592  std::cout << "Beam: " << estNE->Integral() << endl;
593 
594  can->SaveAs("Fig_DCMP.png");
595  can->SaveAs("Fig_DCMP.pdf");
596 }
597 
598 
600 {
601  TGaxis::SetMaxDigits(3);
602 
603 
604  std::string fname = "nue_dcmp.root";
605 
606  std::string dname="nue_decomps/michelcheatDCMP";
607 
608  MichelDecomp *mdcmp = LoadFromFile<MichelDecomp>(fname,dname).release();
609  // Uncomment this for the first run to generate 'MDCMP.root'
610  //TFile *out = new TFile("MDCMP.root","recreate");
611  //mdcmp->SavePlots(out->mkdir("mdcmp"));
612  //out->Close();
613  //return;
614 
615  TFile *in = new TFile(fname.c_str(),"read");
616 
617  TH2F *tnm = (TH2F*)(in->Get((dname+"/TempNumu/hist").c_str()))->Clone(UniqueName().c_str());
618  TH2F *tnc = (TH2F*)(in->Get((dname+"/TempNC/hist").c_str()))->Clone(UniqueName().c_str());
619  TH2F *tne = (TH2F*)(in->Get((dname+"/TempNue/hist").c_str()))->Clone(UniqueName().c_str());
620  TH2F *tda = (TH2F*)(in->Get((dname+"/TempData/hist").c_str()))->Clone(UniqueName().c_str());
621  const std::vector<const TH2F*> hs = {tnm,tnc,tne,tda};
622 
623  TH1F *nm = (TH1F*)(in->Get((dname+"/Numu/hist").c_str()))->Clone(UniqueName().c_str());
624  TH1F *anm= (TH1F*)(in->Get((dname+"/AntiNumu/hist").c_str()))->Clone(UniqueName().c_str());
625  nm->Add(anm);
626  TH1F *ne = (TH1F*)(in->Get((dname+"/Nue/hist").c_str()))->Clone(UniqueName().c_str());
627  TH1F *ane= (TH1F*)(in->Get((dname+"/AntiNue/hist").c_str()))->Clone(UniqueName().c_str());
628  ne->Add(ane);
629  TH1F *nc = (TH1F*)(in->Get((dname+"/NC/hist").c_str()))->Clone(UniqueName().c_str());
630  TH1F *anc = (TH1F*)(in->Get((dname+"/AntiNC/hist").c_str()))->Clone(UniqueName().c_str());
631  nc->Add(anc);
632  TH1F *da = (TH1F*)(in->Get((dname+"/Data/hist").c_str()))->Clone(UniqueName().c_str());
633 
634  TH1D *potda = (TH1D*)(in->Get((dname+"/TempData/pot").c_str()))->Clone(UniqueName().c_str());
635  TH1D *potmc = (TH1D*)(in->Get((dname+"/TempNumu/pot").c_str()))->Clone(UniqueName().c_str());
636  std::cout << "POT " << potda->Integral() << " " << potmc->Integral() << std::endl;
637  const std::vector<const TH1F*> bs = {nm,nc,ne,da};
638 
639  TFile *in2 = new TFile("MDCMP.root", "read");
640  TH1F *hscc = (TH1F*)(in2->Get("mdcmp/numu/numu_ratio"))->Clone(UniqueName().c_str());
641  TH1F *hsnc = (TH1F*)(in2->Get("mdcmp/nc/nc_ratio"))->Clone(UniqueName().c_str());
642  TH1F *hsne = (TH1F*)(in2->Get("mdcmp/nue/nue_ratio"))->Clone(UniqueName().c_str());
643  const std::vector<const TH1F*> scales = {hscc,hsnc,hsne};
644 
645  std::vector<int> idxs = {1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18};
646  double tot = 0;
647 
648  for (int idx : idxs)
649  //MakeStackNME(hs,potda->Integral()/potmc->Integral(),idx,idx,scales,
650  MakeStackNME(hs,potda->Integral()/potmc->Integral(), idx, idx, scales, "NueEnergy/CVN","Events / 8#times10^{20}","Ratio",tot);
651 
652  std::cout << tot << std::endl;
653 
654  MakeScalesPlot(scales);
655 
656  MakeRawPlot(bs,potda->Integral()/potmc->Integral());
657 
658  MakeFracPlot(bs);
659 
660  //MakeResultPlot(mdcmp, da, bs, scales, potda->Integral(), potda->Integral()/potmc->Integral());
661 }
void Nue2018ThreeBinAxis(THStack *axes, bool drawLabels, bool merged, bool coreOnly)
T max(const caf::Proxy< T > &a, T b)
void latex(double x, double y, std::string txt, double ang=0, int align=22, double size=0.042)
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
Spectrum MC_AntiNumuComponent() const override
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
void MakeStackNME(std::vector< const TH2F * > hs, double potrat, int idx1, int idx2, std::vector< const TH1F * > scales, std::string xtitle, std::string ytitle, std::string yrattitle, double &tot, bool applyScales=false)
void MakeFracPlot(std::vector< const TH1F * > hs)
T sqrt(T number)
Definition: d0nt_math.hpp:156
virtual Spectrum AntiNumuComponent() const override
static constexpr Double_t nm
Definition: Munits.h:133
void FormatRatio(TH1 *hrat, std::string xtitle, std::string ytitle)
TFile * fmc
Definition: bdt_com.C:14
Float_t f2
virtual Spectrum NumuComponent() const override
T sqr(T x)
More efficient square function than pow(x,2)
Definition: MathUtil.h:23
virtual Spectrum AntiNueComponent() const override
Spectrum MC_NCTotalComponent() const override
virtual Spectrum NueComponent() const override
#define pot
Float_t f1
void draw_michel_plots()
OStream cout
Definition: OStream.cxx:6
Spectrum MC_NumuComponent() const override
void MakeResultPlot(MichelDecomp *m, TH1 *hda, std::vector< const TH1F * > ht, std::vector< const TH1F * > hs, double pot, double potrat)
ifstream in
Definition: comparison.C:7
Eigen::ArrayXd ProjectionY(const Eigen::MatrixXd &mat)
Helper for WeightingVariable.
void Nue2018ThreeBinLabels(const double yNDC, const double textsize, const int color, const bool nd)
enum BeamMode nc
Spectrum MC_NueComponent() const override
std::string to_string(ModuleType mt)
Definition: ModuleType.h:32
Double_t sum
Definition: plot.C:31
Float_t e
Definition: plot.C:35
void MakeScalesPlot(std::vector< const TH1F * > hs)
Spectrum MC_AntiNueComponent() const override
T max(sqlite3 *const db, std::string const &table_name, std::string const &column_name)
Definition: statistics.h:68
enum BeamMode kGreen
std::string UniqueName()
Return a different string each time, for creating histograms.
Definition: Utilities.cxx:29
const double potmc
virtual Spectrum NCTotalComponent() const override
void MakeRawPlot(std::vector< const TH1F * > hs, double potrat)
enum BeamMode string