plot_nueFD_Signal_REWvsNOM_pTExtrap_FHC.C
Go to the documentation of this file.
1 // Macro to make table of nominal predictions at oscillation values
5 #include "CAFAna/Core/Binning.h"
6 #include "CAFAna/Core/Cut.h"
7 #include "CAFAna/Core/HistAxis.h"
9 #include "CAFAna/Core/Spectrum.h"
12 #include "CAFAna/Cuts/Cuts.h"
16 #include "CAFAna/Cuts/NueCutsSecondAna.h"
17 #include "CAFAna/Cuts/SpillCuts.h"
26 #include "CAFAna/Vars/Vars.h"
27 
28 #include "OscLib/OscCalcPMNSOpt.h"
29 #include "OscLib/OscCalcDumb.h"
30 #include "CAFAna/Analysis/Style.h"
32 
33 #include "TFile.h"
34 #include "TH1.h"
35 #include "TCanvas.h"
36 #include "TGaxis.h"
37 #include "TLatex.h"
38 #include "TLegend.h"
39 #include "TLine.h"
40 #include "TSystem.h"
41 #include "TStyle.h"
42 #include <iostream>
43 #include <iomanip>
44 #include <fstream>
45 #include <tgmath.h>
46 
47 using namespace ana;
48 
49 void Legend(double x0, double y0, double x1, double y1)
50 {
51  // Doesn't handle log-y well
52  // TLegend* leg = AutoPlaceLegend(x1-x0, y1-y0);
53  TLegend* leg = new TLegend(x0, y0, x1, y1);
54  leg->SetFillStyle(0);
55  leg->SetTextSize(2/30.);
56 
57  // leg->SetHeader("#splitline{#bf{Neutrino Beam}}{#splitline{FD prediction}{#nu_{e} signal+WS}}");
58  leg->SetHeader("#splitline{#bf{Neutrino Beam}}{#nu_{e} signal+WS}");
59  TH1* dummy = new TH1F("", "", 1, 0, 1);
60  TH1* dummyFill = new TH1F("", "", 1, 0, 1);
61  leg->AddEntry(dummy->Clone(), "","");
62  dummyFill->SetLineColor(kTotalMCColor);
63  dummyFill->SetFillColor(kTotalMCErrorBandColor);
64 
65  dummy->SetLineColor(kNueSignalColor);
66  leg->AddEntry(dummy->Clone(), "Nominal", "l");
67  dummy->SetLineStyle(9);
68  dummy->SetLineColor(kOrange+1);
69  leg->AddEntry(dummy->Clone(), "Q^{2} reweighted", "l");
70  dummy->SetLineStyle(8);
71  dummy->SetLineColor(kGray+2);
72  leg->AddEntry(dummy->Clone(), "p_{T}/p reweighted", "l");
73  dummy->SetLineStyle(3);
74  // dummy->SetLineColor(kGray+2);
75  leg->AddEntry(dummy->Clone(), "cos #theta reweighted", "l");
76 
77  leg->Draw();
78 }
79 
80 TH1* create_coreperi_TH1(const TH1* coreTH1, const TH1* periTH1)
81 {
82  TH1* coreperiTH1 = (TH1*)coreTH1->Clone("coreperi");
83  for(int i = 18; i < 28; ++i){
84  coreperiTH1->SetBinContent(i, periTH1->GetBinContent(i));
85  }
86 
87  return coreperiTH1;
88 }
89 
90 void plot_nueFD_Signal_REWvsNOM_pTExtrap_FHC(const std::string filename_no, const std::string filename_re_core, const std::string filename_re_peri)
91 
92 {
93  gStyle->SetPadTickX(1);
94  gStyle->SetPadTickY(0);
95 
96  TFile * file_no = new TFile (filename_no.c_str(),"READ");
97  if(file_no->IsZombie()){std::cerr << "Can't find " << filename_no.c_str() << std::endl; return;}
98  TFile * file_re_core = new TFile (filename_re_core.c_str(),"READ");
99  if(file_re_core->IsZombie()){std::cerr << "Can't find " << filename_re_core.c_str() << std::endl; return;}
100  TFile * file_re_peri = new TFile (filename_re_peri.c_str(), "READ");
101  if(file_re_peri->IsZombie()){std::cerr << "Can't find " << filename_re_peri.c_str() << std::endl; return;}
102 
103  osc::OscCalcDumb mycalc;
104 
105  TFile* file = new TFile("FD_KinematicsCorrection_pTExtrap_FHC.root","RECREATE");
106 
107  TDirectory* dpred_noextrap = file_no->GetDirectory("NOReweight/nueAxis_NoExtrap");
108  PredictionNoExtrap* pred_noextrap = ana::LoadFrom<PredictionNoExtrap>(dpred_noextrap).release();
109 
110  //Save the predictions for each quant for nom, core and peri
111 
112  std::vector<const PredictionExtendToPeripheral* > pred_no_quant;
113  std::vector<const PredictionExtendToPeripheral* > pred_re_core_q2_quant;
114  std::vector<const PredictionExtendToPeripheral* > pred_re_core_PtP_quant;
115  std::vector<const PredictionExtendToPeripheral* > pred_re_core_Cos_quant;
116  std::vector<const PredictionExtendToPeripheral* > pred_re_peri_q2_quant;
117  std::vector<const PredictionExtendToPeripheral* > pred_re_peri_PtP_quant;
118  std::vector<const PredictionExtendToPeripheral* > pred_re_peri_Cos_quant;
119 
120  const unsigned int nQuants = 3;
121 
122  for ( unsigned int quant = 0; quant < nQuants; ++quant ){
123  //nominal
124  TDirectory* dpred_no_quant = file_no->GetDirectory(Form("NOReweight/nueAxis_NueSignalExtrap_Quant%i",quant+1));
125  auto extend_nom = new PredictionExtendToPeripheral((ana::LoadFrom<PredictionExtrap>(dpred_no_quant).release()), pred_noextrap, kNue2020Binning, false);
126  pred_no_quant.push_back(extend_nom);
127 
128  //core
129  TDirectory* dpred_re_core_q2_quant = file_re_core->GetDirectory(Form("CORE/nueAxis_NueSignalExtrap_Q2Weight_Quant%i", quant+1));
130  auto extend_core_q2 = new PredictionExtendToPeripheral((ana::LoadFrom<PredictionExtrap>(dpred_re_core_q2_quant).release()), pred_noextrap, kNue2020Binning, false);
131  pred_re_core_q2_quant.push_back(extend_core_q2);
132 
133  TDirectory* dpred_re_core_PtP_quant = file_re_core->GetDirectory(Form("CORE/nueAxis_NueSignalExtrap_PtPWeight_Quant%i", quant+1));
134  auto extend_core_PtP = new PredictionExtendToPeripheral((ana::LoadFrom<PredictionExtrap>(dpred_re_core_PtP_quant).release()), pred_noextrap, kNue2020Binning, false);
135  pred_re_core_PtP_quant.push_back(extend_core_PtP);
136 
137  TDirectory* dpred_re_core_Cos_quant = file_re_core->GetDirectory(Form("CORE/nueAxis_NueSignalExtrap_CosWeight_Quant%i", quant+1));
138  auto extend_core_Cos = new PredictionExtendToPeripheral((ana::LoadFrom<PredictionExtrap>(dpred_re_core_Cos_quant).release()), pred_noextrap, kNue2020Binning, false);
139  pred_re_core_Cos_quant.push_back(extend_core_Cos);
140 
141  //peri
142  TDirectory* dpred_re_peri_q2_quant = file_re_peri->GetDirectory(Form("PERI/nueAxis_NueSignalExtrap_Q2Weight_Quant%i", quant+1));
143  auto extend_peri_q2 = new PredictionExtendToPeripheral((ana::LoadFrom<PredictionExtrap>(dpred_re_peri_q2_quant).release()), pred_noextrap, kNue2020Binning, false);
144  pred_re_peri_q2_quant.push_back(extend_peri_q2);
145 
146  TDirectory* dpred_re_peri_PtP_quant = file_re_peri->GetDirectory(Form("PERI/nueAxis_NueSignalExtrap_PtPWeight_Quant%i", quant+1));
147  auto extend_peri_PtP = new PredictionExtendToPeripheral((ana::LoadFrom<PredictionExtrap>(dpred_re_peri_PtP_quant).release()), pred_noextrap, kNue2020Binning, false);
148  pred_re_peri_PtP_quant.push_back(extend_peri_PtP);
149 
150  TDirectory* dpred_re_peri_Cos_quant = file_re_peri->GetDirectory(Form("PERI/nueAxis_NueSignalExtrap_CosWeight_Quant%i", quant+1));
151  auto extend_peri_Cos = new PredictionExtendToPeripheral((ana::LoadFrom<PredictionExtrap>(dpred_re_peri_Cos_quant).release()), pred_noextrap, kNue2020Binning, false);
152  pred_re_peri_Cos_quant.push_back(extend_peri_Cos);
153  }
154 
155  //Add quants to form a single prediction for nom, core and peri
156  //add to form Spectrums (so mycalc already implemented)
157 
158  Spectrum spec_nom = pred_no_quant.front()->PredictComponent(&mycalc,
160  Current::kCC,
161  Sign::kBoth);
162  spec_nom.Clear();
163 
164  Spectrum spec_core_q2 = pred_re_core_q2_quant.front()->PredictComponent(&mycalc,
166  Current::kCC,
167  Sign::kBoth);
168  spec_core_q2.Clear();
169 
170  Spectrum spec_core_PtP = pred_re_core_PtP_quant.front()->PredictComponent(&mycalc,
172  Current::kCC,
173  Sign::kBoth);
174  spec_core_PtP.Clear();
175 
176  Spectrum spec_core_Cos = pred_re_core_Cos_quant.front()->PredictComponent(&mycalc,
178  Current::kCC,
179  Sign::kBoth);
180  spec_core_Cos.Clear();
181 
182  Spectrum spec_peri_q2 = pred_re_peri_q2_quant.front()->PredictComponent(&mycalc,
184  Current::kCC,
185  Sign::kBoth);
186  spec_peri_q2.Clear();
187 
188  Spectrum spec_peri_PtP = pred_re_peri_PtP_quant.front()->PredictComponent(&mycalc,
190  Current::kCC,
191  Sign::kBoth);
192  spec_peri_PtP.Clear();
193 
194  Spectrum spec_peri_Cos = pred_re_peri_Cos_quant.front()->PredictComponent(&mycalc,
196  Current::kCC,
197  Sign::kBoth);
198  spec_peri_Cos.Clear();
199 
200 
201 
202  for (unsigned int quant = 0; quant < nQuants; ++quant ){
203  spec_nom += pred_no_quant[quant]->PredictComponent(&mycalc,
205  Current::kCC,
206  Sign::kBoth);
207 
208  spec_core_q2 += pred_re_core_q2_quant[quant]->PredictComponent(&mycalc,
210  Current::kCC,
211  Sign::kBoth);
212 
213  spec_core_PtP += pred_re_core_PtP_quant[quant]->PredictComponent(&mycalc,
215  Current::kCC,
216  Sign::kBoth);
217 
218  spec_core_Cos += pred_re_core_Cos_quant[quant]->PredictComponent(&mycalc,
220  Current::kCC,
221  Sign::kBoth);
222 
223  spec_peri_q2 += pred_re_peri_q2_quant[quant]->PredictComponent(&mycalc,
225  Current::kCC,
226  Sign::kBoth);
227 
228  spec_peri_PtP += pred_re_peri_PtP_quant[quant]->PredictComponent(&mycalc,
230  Current::kCC,
231  Sign::kBoth);
232 
233  spec_peri_Cos += pred_re_peri_Cos_quant[quant]->PredictComponent(&mycalc,
235  Current::kCC,
236  Sign::kBoth);
237 
238  }
239 
240  //form histograms and merge core+peri to form a single prediction for FD
241  // POT set by hand from total POT of ND sample
242  const double kPOTFD = 13.5e20;
243 
244  TH1* hAll_no = spec_nom.ToTH1(kPOTFD);
245 
246  TH1* hAll_core_q2 = spec_core_q2.ToTH1(kPOTFD);
247  TH1* hAll_peri_q2 = spec_peri_q2.ToTH1(kPOTFD);
248  TH1* hAll_re_q2 = (TH1*)(create_coreperi_TH1(hAll_core_q2, hAll_peri_q2))->Clone("hAll_re_q2");
249 
250  TH1* hAll_core_PtP = spec_core_PtP.ToTH1(kPOTFD);
251  TH1* hAll_peri_PtP = spec_peri_PtP.ToTH1(kPOTFD);
252  TH1* hAll_re_PtP = (TH1*)(create_coreperi_TH1(hAll_core_PtP, hAll_peri_PtP))->Clone("hAll_re_PtP");
253 
254  TH1* hAll_core_Cos = spec_core_Cos.ToTH1(kPOTFD);
255  TH1* hAll_peri_Cos = spec_peri_Cos.ToTH1(kPOTFD);
256  TH1* hAll_re_Cos = (TH1*)(create_coreperi_TH1(hAll_core_Cos, hAll_peri_Cos))->Clone("hAll_re_Cos");
257 
258  //Set hist styles
259  hAll_re_q2->SetLineColor(kOrange+1);
260  hAll_re_q2->SetLineStyle(9);
261  // hAll_re_q2->SetLineWidth(1);
262 
263  hAll_re_PtP->SetLineColor(kGray+2);
264  hAll_re_PtP->SetLineStyle(8);
265  hAll_re_PtP->SetLineWidth(1);
266 
267  hAll_re_Cos->SetLineColor(kGray+2);
268  hAll_re_Cos->SetLineStyle(3);
269  hAll_re_Cos->SetLineWidth(1);
270 
271  //create a canvas
272  TCanvas *c2 =new TCanvas("c2","c2",500,400);
273 
274  TPad *padAbs = new TPad("padAbs", "padAbs", 0, 0.3, 1, 1.0);
275  padAbs-> SetBottomMargin(0);
276  padAbs->SetFillStyle(0);
277  padAbs-> SetLeftMargin(0.15);
278 
279  padAbs-> Draw();
280  padAbs-> cd();
281 
282  hAll_re_q2->GetXaxis()->CenterTitle();
283  hAll_re_q2->GetYaxis()->CenterTitle();
284  hAll_re_q2->GetXaxis()->SetDecimals();
285  hAll_re_q2->GetYaxis()->SetDecimals();
286  hAll_re_q2->GetYaxis()->SetTitleOffset(.8);
287  hAll_re_q2->GetYaxis()->SetTitleSize(.06);
288  hAll_re_q2->GetYaxis()->SetLabelSize(.06);
289  hAll_re_q2->GetXaxis()->SetLabelSize(.0);
290  hAll_re_q2->SetMaximum(1.1*hAll_re_q2->GetMaximum());
291  hAll_re_q2->GetYaxis()->SetTitle(TString::Format("Events / %.02lf #times 10^{21} POT", kPOTFD/1E21));
292 
293  hAll_no->Draw("hist");
294  hAll_re_PtP->Draw("same hist");
295  hAll_re_Cos->Draw("same hist");
296  hAll_re_q2->Draw("same hist");
297 
298  Nue2018ThreeBinAxis(hAll_re_q2, true, false, false);
299  hAll_re_q2->GetYaxis()->SetRangeUser(.001,20);
300 
301  //Fix up the axis
302  padAbs->RedrawAxis();
303  padAbs->Update();
304 
305  double legendxlow=0.55;
306  double legendxhigh=0.95;
307 
308  Legend(legendxlow, .5, legendxhigh, .85);
309 
310  TLatex* selTitle = new TLatex(.9, .95, "NOvA Simulation");
311  selTitle->SetNDC();
312  selTitle->SetTextColor(kGray+1);
313  selTitle->SetTextSize(2/30.);
314  selTitle->SetTextAlign(32);
315  selTitle->Draw();
316 
317  c2->cd();
318  TPad *padRatio = new TPad("padRatio", "padRatio", 0, 0.0, 1, 0.33);
319  padRatio -> SetTopMargin(0.11);
320  padRatio -> SetBottomMargin(0.33);
321  padRatio -> SetFillStyle(0);
322  padRatio -> SetLeftMargin(0.15);
323 
324  padRatio -> Draw();
325  padRatio -> cd();
326 
327  TH1* hRatio_q2 = (TH1*)hAll_re_q2 -> Clone("hRatio_q2");
328  TH1* hRatio_PtP = (TH1*)hAll_re_PtP -> Clone("hRatio_PtP");
329  TH1* hRatio_Cos = (TH1*)hAll_re_Cos -> Clone("hRatio_Cos");
330  TLine* lOne = new TLine(0,1.0,27.0,1.0);
331  lOne -> SetLineStyle(2);
332 
333  hRatio_q2 -> Divide(hAll_no);
334  hRatio_PtP -> Divide(hAll_no);
335  hRatio_Cos -> Divide(hAll_no);
336 
337  hRatio_q2->GetXaxis()->CenterTitle();
338  hRatio_q2->GetYaxis()->CenterTitle();
339  hRatio_q2->GetYaxis()->SetLabelSize(.12);
340  hRatio_q2->GetXaxis()->SetLabelSize(.12);
341  hRatio_q2->GetXaxis()->SetTitleOffset(1.0);
342  hRatio_q2->GetXaxis()->SetTickSize(.1);
343  hRatio_q2->GetYaxis()->SetTitleSize(.12);
344  hRatio_q2->GetXaxis()->SetTitleSize(.12);
345  hRatio_q2->SetMarkerSize(.1);
346  hRatio_q2->GetXaxis()->SetDecimals();
347  hRatio_q2->GetYaxis()->SetDecimals();
348  hRatio_q2->GetYaxis()->SetTitleOffset(0.4);
349  hRatio_q2->GetYaxis()->SetTitle("#splitline{Ratio to}{nominal}");
350 
351  hRatio_q2->GetYaxis()->SetNdivisions(502,kFALSE);
352  hRatio_q2 -> Draw("hist");
353  hRatio_PtP -> Draw("same hist");
354  hRatio_Cos -> Draw("same hist");
355  lOne -> Draw("same");
356  hRatio_q2 -> SetMinimum(0.8);
357  Nue2018ThreeBinAxis(hRatio_q2, true, false, false);
358  hRatio_q2->GetYaxis()->SetRangeUser(0.95,1.05);
359 
360  c2->Print(("plots/FDnue_pred_recoE_pTExtrap_FHC.png"));
361  c2->SaveAs("plots/FDnue_pred_recoE_pTExtrap_FHC.pdf","pdf");
362 
363  hRatio_q2-> Write("trueQ2_FD_KinematicsCorrection_FHC");
364  hRatio_PtP-> Write("PtP_FD_KinematicsCorrection_FHC");
365  hRatio_Cos-> Write("Cos_FD_KinematicsCorrection_FHC");
366  // file->Close();
367 
368  //Save the integral stats in a tex file
369  double intAll_no = hAll_no->Integral();
370  double intAll_re_q2 = hAll_re_q2->Integral();
371  double intAll_re_PtP = hAll_re_PtP->Integral();
372  double intAll_re_Cos = hAll_re_Cos->Integral();
373 
374  std::cout << intAll_no << std::endl;
375  std::cout << intAll_re_q2 << std::endl;
376  std::cout << intAll_re_PtP << std::endl;
377  std::cout << intAll_re_Cos << std::endl;
378 
379  std::ofstream tfile;
380 
381  tfile.open("plots/FDextrap-TableOut-pTExtrap-FHC.tex", ios::out);
382 
383  tfile << "\\begin{tabular}{lll} \n";
384  tfile << "\\multicolumn{3}{c}{prop FD $\nu_e$ signal + WS} \\\\ \n";
385  tfile << "Sample & Signal+WS & /NOMINAL \\\\ \n";
386  tfile << "\\hline \n";
387 
388  tfile << std::fixed;
389  tfile.precision(2);
390  tfile << "NOMINAL" << " & " << round(100.0*intAll_no)/100.0 << " & \\\\ \n";
391  tfile << "trueQ2 REW." << " & " << round(100.0*intAll_re_q2)/100.0 << " & " << round(1000.0*((intAll_re_q2/intAll_no)-1.0))/10.0 << " \\% \\\\ \n";
392  tfile << "PtP REW." << " & " << round(100.0*intAll_re_PtP)/100.0 << " & " << round(1000.0*((intAll_re_PtP/intAll_no)-1.0))/10.0 << " \\% \\\\ \n";
393  tfile << "Cos REW." << " & " << round(100.0*intAll_re_Cos)/100.0 << " & " << round(1000.0*((intAll_re_Cos/intAll_no)-1.0))/10.0 << " \\% \\\\ \n";
394  tfile << "\\hline \n";
395  tfile << "\\end{tabular}";
396  tfile.close();
397 
398 } // end
399 
400 
401 
void Nue2018ThreeBinAxis(THStack *axes, bool drawLabels, bool merged, bool coreOnly)
tree Draw("slc.nhit")
enum BeamMode kOrange
ratio_hxv Divide(hxv, goal_hxv)
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
Float_t y1[n_points_granero]
Definition: compare.C:5
(&#39; appearance&#39;)
Definition: IPrediction.h:18
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
Float_t x1[n_points_granero]
Definition: compare.C:5
hmean SetLineStyle(2)
TH1 * create_coreperi_TH1(const TH1 *coreTH1, const TH1 *periTH1)
OStream cerr
Definition: OStream.cxx:7
void Clear()
Definition: Spectrum.cxx:361
const Color_t kTotalMCErrorBandColor
Definition: Style.h:17
fvar< T > round(const fvar< T > &x)
Definition: round.hpp:23
ntuple SetFillStyle(1001)
void plot_nueFD_Signal_REWvsNOM_pTExtrap_FHC(const std::string filename_no, const std::string filename_re_core, const std::string filename_re_peri)
Representation of a spectrum in any variable, with associated POT.
Definition: Spectrum.h:40
const Color_t kNueSignalColor
Definition: Style.h:19
c2
Definition: demo5.py:33
Charged-current interactions.
Definition: IPrediction.h:39
void Legend(double x0, double y0, double x1, double y1)
const Binning kNue2020Binning
Simple oscillation probability calculator that has no solar term or mass hierarchy or delta so it&#39;s s...
Definition: OscCalcDumb.h:16
OStream cout
Definition: OStream.cxx:6
TFile * file
Definition: cellShifts.C:17
Both neutrinos and antineutrinos.
Definition: IPrediction.h:52
void Format(TGraph *gr, int lcol, int lsty, int lwid, int mcol, int msty, double msiz)
Definition: Style.cxx:154
const Color_t kTotalMCColor
Definition: Style.h:16
Prediction that just uses FD MC, with no extrapolation.
histo SetMinimum(1.0e-9)
c cd(1)
gm Write()
enum BeamMode string