plot_nueFD_Signal_REWvsNOM.C
Go to the documentation of this file.
1 // Macro to make table of nominal predictions at oscillation values
2 
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"
14 #include "CAFAna/Cuts/NueCutsSecondAna.h"
15 #include "CAFAna/Cuts/SpillCuts.h"
23 #include "CAFAna/Vars/Vars.h"
24 
25 #include "OscLib/OscCalcPMNSOpt.h"
26 #include "OscLib/OscCalcDumb.h"
27 #include "CAFAna/Analysis/Style.h"
28 
29 //#include "CAFAna/nue/SecondAna/draw_plots_util.h"
30 
31 #include "TFile.h"
32 #include "TH1.h"
33 #include "TCanvas.h"
34 #include "TGaxis.h"
35 #include "TLatex.h"
36 #include "TLegend.h"
37 #include "TLine.h"
38 #include "TSystem.h"
39 #include <iostream>
40 #include <iomanip>
41 #include <fstream>
42 #include <tgmath.h>
43 
44 using namespace ana;
45 
46 void Legend(double x0, double y0, double x1, double y1)
47 {
48  // Doesn't handle log-y well
49  // TLegend* leg = AutoPlaceLegend(x1-x0, y1-y0);
50  TLegend* leg = new TLegend(x0, y0, x1, y1);
51  leg->SetFillStyle(0);
52 
53  TH1* dummy = new TH1F("", "", 1, 0, 1);
54  TH1* dummyFill = new TH1F("", "", 1, 0, 1);
55  dummyFill->SetLineColor(kTotalMCColor);
56  dummyFill->SetFillColor(kTotalMCErrorBandColor);
57 
58  dummy->SetLineColor(kNueSignalColor);
59  leg->AddEntry(dummy->Clone(), "#nu_{e} signal NOM.", "l");
60  dummy->SetLineStyle(3);
61  leg->AddEntry(dummy->Clone(), "trueQ2 REW.", "l");
62  dummy->SetLineStyle(8);
63  leg->AddEntry(dummy->Clone(), "PtP REW.", "l");
64  dummy->SetLineStyle(9);
65  leg->AddEntry(dummy->Clone(), "CosTheta REW.", "l");
66 
67  leg->Draw();
68 }
69 
70 void plot_nueFD_Signal_REWvsNOM(const std::string filename_no, const std::string filename_re
71  )
72 
73 {
74  const int kNumSels = 1;
75  const int kNumVars = 1;
76 
77  struct HistDef
78  {
80  };
81 
82  const HistDef defs[kNumVars] = {
83  "recoE",
84  };
85 
86  const std::string selNames[kNumSels] = {
87  "AllSamples"
88  };
89 
90  TFile * file_no = new TFile (filename_no.c_str(),"READ");
91  if(file_no->IsZombie()){std::cerr << "Can't find " << filename_no.c_str() << std::endl; return;}
92  TFile * file_re = new TFile (filename_re.c_str(),"READ");
93  if(file_re->IsZombie()){std::cerr << "Can't find " << filename_re.c_str() << std::endl; return;}
94 
95  osc::OscCalcPMNSOpt mycalc;
96  ResetOscCalcToDefault(&mycalc);
97 
98  mycalc.SetDmsq32(+fabs(mycalc.GetDmsq32()));
99  mycalc.SetdCP(3*TMath::Pi()/2);
100  mycalc.SetTh23(asin(sqrt(0.404)));
101 
102  std::ofstream tfile;
103  const int kNumTabs = 1;
104  const std::string tabNames[kNumTabs] = {
105  "AllSamples"
106  };
107 
108 
109  tfile.open("plots/FDextrap-TableOut.tex", ios::out);
110 
111  tfile << "\\begin{tabular}{lllll} \n";
112  tfile << "\\multicolumn{5}{c}{prop FD $\nu_e$ signal ALL SAMPLES} \\\\ \n";
113  tfile << "Sample & Sig.+Bkg. & /NOMINAL & Signal & /NOMINAL \\\\ \n";
114  tfile << "\\hline \n";
115 
116  double intTotalFull, intAllFull, intSigFull, intNCFull, intCCFull, intBeamFull, intTauFull;
117 
118  for(int selIdx = 0; selIdx < kNumSels; ++selIdx){
119  const std::string selName = selNames[selIdx];
120  for(int varIdx = 0; varIdx < kNumVars; ++varIdx){
121  bool pidaxes = varIdx;
122  const std::string varName = defs[varIdx].name;
123  TDirectory * dpred_no = file_no->GetDirectory("prediction/nue_pred_AllSamples_recoE");
124  TDirectory * dpred_re = file_re->GetDirectory("prediction/nue_pred_trueQ2_AllSamples_recoE");
125  TDirectory * dpred_re_PtP = file_re->GetDirectory("prediction/nue_pred_PtP_AllSamples_recoE");
126  TDirectory * dpred_re_Cos = file_re->GetDirectory("prediction/nue_pred_Cos_AllSamples_recoE");
127 
128  std::cout << "HELLO" << endl;
129  //FD Plots
130  auto pred_no = ana::LoadFrom<IPrediction> (dpred_no);
131  auto pred_re = ana::LoadFrom<IPrediction> (dpred_re);
132  auto pred_re_PtP = ana::LoadFrom<IPrediction> (dpred_re_PtP);
133  auto pred_re_Cos = ana::LoadFrom<IPrediction> (dpred_re_Cos);
134 
135  std::cout << "HEllO" << endl;
136  if(selIdx == 0 && varIdx == 0)
137  std::cout << "MC POT: "
138  << pred_no->Predict(&mycalc).POT() << std::endl;
139 
140  //play with the binning choice for given variables:
141  int rebinFactor=1;
142 
143  //create a canvas
144  TCanvas *c2 =new TCanvas("c2","c2",500,400);
145 
146  TPad *padAbs = new TPad("padAbs", "padAbs", 0, 0.3, 1, 1.0);
147  padAbs -> SetBottomMargin(0);
148  padAbs->SetFillStyle(0);
149  padAbs -> SetLeftMargin(0.15);
150 
151  padAbs -> Draw();
152  padAbs -> cd();
153 
154  // POT set by hand from total POT of ND sample
155  const double kPOTFD = 8.09e20;
156 
157  //Backgrounds
158  TH1* hNC = pred_no.get()->PredictComponent(&mycalc,
160  Current::kNC,
161  Sign::kBoth).ToTH1(kPOTFD);
162 
163  TH1* hCC = pred_no.get()->PredictComponent(&mycalc,
165  Current::kCC,
166  Sign::kBoth).ToTH1(kPOTFD);
167 
168  TH1* hBeam = pred_no.get()->PredictComponent(&mycalc,
170  Current::kCC,
171  Sign::kBoth).ToTH1(kPOTFD);
172 
173  TH1* hAll = (TH1*)hNC->Clone(UniqueName().c_str());
174  hAll->Add(hCC);
175  hAll->Add(hBeam);
176 
177  //Signals
178  TH1* hSig_no = pred_no.get()->PredictComponent(&mycalc,
180  Current::kCC,
181  Sign::kBoth).ToTH1(kPOTFD);
182  hSig_no->SetLineColor(kNueSignalColor);
183  hSig_no->Rebin(rebinFactor);
184 
185  TH1* hSig_re = pred_re.get()->PredictComponent(&mycalc,
187  Current::kCC,
188  Sign::kBoth).ToTH1(kPOTFD);
189  hSig_re->SetLineColor(kNueSignalColor);
190  hSig_re->SetLineStyle(3);
191  hSig_re->Rebin(rebinFactor);
192 
193  TH1* hSig_re_PtP = pred_re_PtP.get()->PredictComponent(&mycalc,
195  Current::kCC,
196  Sign::kBoth).ToTH1(kPOTFD);
197  hSig_re_PtP->SetLineColor(kNueSignalColor);
198  hSig_re_PtP->SetLineStyle(8);
199  hSig_re_PtP->Rebin(rebinFactor);
200 
201  TH1* hSig_re_Cos = pred_re_Cos.get()->PredictComponent(&mycalc,
203  Current::kCC,
204  Sign::kBoth).ToTH1(kPOTFD);
205  hSig_re_Cos->SetLineColor(kNueSignalColor);
206  hSig_re_Cos->SetLineStyle(9);
207  hSig_re_Cos->Rebin(rebinFactor);
208 
209  double intAll = hAll->Integral();
210  double intSig_no = hSig_no->Integral();
211  double intSig_re = hSig_re->Integral();
212  double intSig_re_PtP = hSig_re_PtP->Integral();
213  double intSig_re_Cos = hSig_re_Cos->Integral();
214 
215  tfile << std::fixed;
216  tfile.precision(2);
217  tfile << "NOMINAL" << " & " << round(100.0*(intSig_no+intAll))/100.0 << " & 0.00 & ";
218  tfile << round(100.0*intSig_no)/100.0 << " & 0.00 \\\\ \n";
219  tfile << "trueQ2 REW." << " & " << round(100.0*(intSig_re+intAll))/100.0 << " & " << round(1000.0*((intSig_re+intAll)/(intSig_no+intAll)-1.0))/10.0 << " \\% & ";
220  tfile << round(100.0*intSig_re)/100.0 << " & "<< round(1000.0*((intSig_re/intSig_no)-1.0))/10.0 << " \\% \\\\ \n";
221  tfile << "PtP REW." << " & " << round(100.0*(intSig_re_PtP+intAll))/100.0 << " & " << round(1000.0*((intSig_re_PtP+intAll)/(intSig_no+intAll)-1.0))/10.0 << " \\% & ";
222  tfile << round(100.0*intSig_re_PtP)/100.0 << " & " << round(1000.0*((intSig_re_PtP/intSig_no)-1.0))/10.0 << " \\% \\\\ \n";
223  tfile << "Cos REW." << " & " << round(100.0*(intSig_re_Cos+intAll))/100.0 << " & " << round(1000.0*((intSig_re_Cos+intAll)/(intSig_no+intAll)-1.0))/10.0 << " \\% & ";
224  tfile << round(100.0*intSig_re_Cos)/100.0 << " & " << round(1000.0*((intSig_re_Cos/intSig_no))-1.0)/10.0 << " \\% \\\\ \n";
225  tfile << "\\hline \n";
226  tfile << "\\end{tabular}";
227  tfile.close();
228 
229  hSig_re->GetXaxis()->CenterTitle();
230  hSig_re->GetYaxis()->CenterTitle();
231  hSig_re->GetXaxis()->SetDecimals();
232  hSig_re->GetYaxis()->SetDecimals();
233  hSig_re->GetYaxis()->SetTitleOffset(1.15);
234  hSig_re->GetYaxis()->SetRangeUser(.001,1.005*TMath::Max(hSig_re->GetMaximum(),hSig_no->GetMaximum()));
235  hSig_re->GetYaxis()->SetLabelSize(3.7/70.);
236  hSig_re->GetXaxis()->SetLabelSize(.0);
237  hSig_re->SetMaximum(1.1*TMath::Max(hSig_re->GetMaximum(),hSig_no->GetMaximum()));
238  hSig_re->GetYaxis()->SetTitle(TString::Format("Events / %.02lf #times 10^{20} POT", kPOTFD/1E20));
239 
240  if(varName=="recoE"){
241  hSig_re->GetXaxis()->SetTitle("Reconstructed neutrino energy (GeV)");
242  }
243 
244  hSig_re->GetXaxis()->SetNdivisions(405,kFALSE);
245 
246  hSig_re->DrawCopy("hist");
247  hSig_no->Draw("same hist");
248  hSig_re_PtP->Draw("same hist");
249  hSig_re_Cos->Draw("same hist");
250 
251  //Fix up the axis
252  padAbs->RedrawAxis();
253  padAbs->Update();
254 
255  double legendxlow=0.5;
256  double legendxhigh=0.95;
257 
258  Legend(legendxlow, .55, legendxhigh, .85);
259 
260  TLatex* selTitle = new TLatex(.935, .95, "FD #nu_{e} signal reweighted");
261  selTitle->SetNDC();
262  selTitle->SetTextSize(2/30.);
263  selTitle->SetTextAlign(32);
264  selTitle->Draw();
265 
266  c2->cd();
267  TPad *padRatio = new TPad("padRatio", "padRatio", 0, 0.05, 1, 0.3);
268  padRatio -> SetTopMargin(0);
269  padRatio -> SetBottomMargin(0.25);
270  padRatio -> SetFillStyle(0);
271  padRatio -> SetLeftMargin(0.15);
272 
273  padRatio -> Draw();
274  padRatio -> cd();
275 
276  TH1* hRatio = (TH1*)hSig_re -> Clone("hRatio");
277  TH1* hRatio_PtP = (TH1*)hSig_re_PtP -> Clone("hRatio_PtP");
278  TH1* hRatio_Cos = (TH1*)hSig_re_Cos -> Clone("hRatio_Cos");
279  TLine* lOne = new TLine(0,1.0,5.0,1.0);
280  lOne -> SetLineStyle(2);
281 
282  hRatio -> Divide(hSig_no);
283  hRatio_PtP -> Divide(hSig_no);
284  hRatio_Cos -> Divide(hSig_no);
285 
286  hRatio->GetXaxis()->CenterTitle();
287  hRatio->GetYaxis()->CenterTitle();
288  hRatio->GetYaxis()->SetLabelSize(3.7/30.);
289  hRatio->GetXaxis()->SetLabelSize(3.7/30.);
290  hRatio->GetXaxis()->SetLabelOffset(.04);
291  hRatio->GetXaxis()->SetTickSize(.06);
292  hRatio->GetYaxis()->SetTitleSize(4.5/30.);
293  hRatio->GetXaxis()->SetTitleSize(4.5/30.);
294  hRatio->GetYaxis()->SetRangeUser(0.7,1.3);
295  hRatio->SetMarkerSize(.1);
296  hRatio->GetXaxis()->SetDecimals();
297  hRatio->GetYaxis()->SetDecimals();
298  hRatio->GetYaxis()->SetTitleOffset(0.4);
299  hRatio->GetYaxis()->SetTitle("REW./NOM.");
300 
301  if(varName=="recoE"){
302  hRatio->GetXaxis()->SetTitle("Reconstructed neutrino energy (GeV)");
303  }
304 
305  hRatio->GetYaxis()->SetNdivisions(502,kFALSE);
306  hRatio->GetXaxis()->SetNdivisions(405,kFALSE);
307 
308  hRatio -> Draw("hist");
309  hRatio_PtP -> Draw("hist same");
310  hRatio_Cos -> Draw("hist same");
311  lOne -> Draw("same");
312 
313  c2->Print(("plots/FDnue_pred_recoE_abs.png"));
314  } // end for varIdx
315  } // end for selIdx
316 } // end
tree Draw("slc.nhit")
ratio_hxv Divide(hxv, goal_hxv)
const int kNumVars
Definition: vars.h:14
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
void ResetOscCalcToDefault(osc::IOscCalcAdjustable *calc)
Reset calculator to default assumptions for all parameters.
Definition: Calcs.cxx:23
fvar< T > fabs(const fvar< T > &x)
Definition: fabs.hpp:15
Float_t y1[n_points_granero]
Definition: compare.C:5
(&#39; appearance&#39;)
Definition: IPrediction.h:18
Float_t x1[n_points_granero]
Definition: compare.C:5
(&#39;beam &#39;)
Definition: IPrediction.h:15
hmean SetLineStyle(2)
T sqrt(T number)
Definition: d0nt_math.hpp:156
std::string name
Definition: NuePlotLists.h:12
OStream cerr
Definition: OStream.cxx:7
const Color_t kTotalMCErrorBandColor
Definition: Style.h:17
fvar< T > round(const fvar< T > &x)
Definition: round.hpp:23
ntuple SetFillStyle(1001)
const Color_t kNueSignalColor
Definition: Style.h:19
c2
Definition: demo5.py:33
Charged-current interactions.
Definition: IPrediction.h:39
Optimized version of OscCalcPMNS.
Definition: StanTypedefs.h:31
const HistDef defs[kNumVars]
Definition: vars.h:15
void SetTh23(const T &th23) override
virtual T GetDmsq32() const
Definition: IOscCalc.h:91
std::unique_ptr< IPrediction > LoadFrom< IPrediction >(TDirectory *dir, const std::string &label)
Definition: IPrediction.cxx:18
const int kNumSels
Definition: vars.h:43
OStream cout
Definition: OStream.cxx:6
void SetDmsq32(const T &dmsq32) override
void SetdCP(const T &dCP) override
Neutral-current interactions.
Definition: IPrediction.h:40
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
All neutrinos, any flavor.
Definition: IPrediction.h:26
const std::string selNames[kNumSels]
Definition: vars.h:46
void Legend(double x0, double y0, double x1, double y1)
std::string UniqueName()
Return a different string each time, for creating histograms.
Definition: Utilities.cxx:29
c cd(1)
T asin(T number)
Definition: d0nt_math.hpp:60
void plot_nueFD_Signal_REWvsNOM(const std::string filename_no, const std::string filename_re)
enum BeamMode string