print_nus17_fd_cut_tables.C
Go to the documentation of this file.
1 // This macro is used to make the nue energy blessed plots. (see doc 11789)
2 // Oscillation parameters are same as those in nue_ana_basic macro
3 
4 #include "CAFAna/Analysis/Calcs.h"
5 #include "CAFAna/Analysis/Exposures.h"
6 #include "CAFAna/Cuts/Cuts.h"
7 #include "CAFAna/Cuts/SpillCuts.h"
8 #include "CAFAna/Cuts/NueCutsFirstAna.h"
9 #include "CAFAna/Core/Utilities.h"
10 #include "CAFAna/Core/SpectrumLoader.h"
11 #include "CAFAna/Core/Spectrum.h"
12 #include "CAFAna/Vars/Vars.h"
13 #include "CAFAna/Cuts/TimingCuts.h"
14 #include "CAFAna/Analysis/Plots.h"
15 #include "CAFAna/Analysis/Style.h"
16 #include "3FlavorAna/Plotting/NuePlotStyle.h"
17 #include "CAFAna/Core/LoadFromFile.h"
18 
19 #include "OscLib/OscCalcDumb.h"
20 #include "OscLib/OscCalcSterile.h"
21 #include "CAFAna/Prediction/PredictionNoExtrap.h"
22 #include "NuXAna/Prediction/PredictionSterile.h"
23 
24 #include "NuXAna/macros/Nus17/nus17_fd_arrays.h"
25 
26 #include "TH1D.h"
27 #include "TH2.h"
28 #include "TCanvas.h"
29 #include "TColor.h"
30 #include "TGraph.h"
31 #include "THStack.h"
32 #include "TMultiGraph.h"
33 #include "TLine.h"
34 #include "TLegend.h"
35 #include "TString.h"
36 //#include "Utilities/rootlogon.C"
37 
38 using namespace ana;
39 
40 
41 const double pot_val = kAna2017POT; // POT
42 const double livetime = kAna2017Livetime; // seconds
43 
44 
45 void GetSpectra(IPrediction* pred, osc::IOscCalc* calc,
46  TH1*& hSig, TH1*& hNue, TH1*& hCC, TH1*& hTau, TH1*& hTotBkg)
47 {
48  hSig = pred->PredictComponent(calc,
49  Flavors::kAll,
50  Current::kNC,
51  Sign::kBoth).ToTH1(pot_val);
52  hSig->SetLineColor(kNCBackgroundColor);
53 
54 
55  hNue = pred->PredictComponent(calc,
56  Flavors::kAllNuE,
57  Current::kCC,
58  Sign::kBoth).ToTH1(pot_val);
59  hNue->SetLineColor(kNueSignalColor);
60 
61  hCC = pred->PredictComponent(calc,
62  Flavors::kAllNuMu,
63  Current::kCC,
64  Sign::kBoth).ToTH1(pot_val);
65  hCC->SetLineColor(kNumuBackgroundColor);
66 
67  hTau = pred->PredictComponent(calc,
68  Flavors::kAllNuTau,
69  Current::kCC,
70  Sign::kBoth).ToTH1(pot_val);
71  hTau->SetLineColor(kMagenta);
72 
73  hTotBkg = (TH1*)hNue->Clone(UniqueName().c_str());
74  hTotBkg->Add(hCC);
75  hTotBkg->Add(hTau);
76 
77  hTotBkg->SetLineColor(kTotalMCColor);
78 }
79 void GetSpectra2D(IPrediction* pred, osc::IOscCalc* calc,
80  TH2*& hSig, TH2*&hTotBkg)
81 {
82  hSig = pred->PredictComponent(calc,
83  Flavors::kAll,
84  Current::kNC,
85  Sign::kBoth).ToTH2(pot_val);
86 
87  TH1* hNue = pred->PredictComponent(calc,
88  Flavors::kNuMuToNuE,
89  Current::kCC,
90  Sign::kBoth).ToTH2(pot_val);
91 
92  TH1* hCC = pred->PredictComponent(calc,
93  Flavors::kAllNuMu,
94  Current::kCC,
95  Sign::kBoth).ToTH2(pot_val);
96 
97  TH1* hBeam = pred->PredictComponent(calc,
98  Flavors::kNuEToNuE,
99  Current::kCC,
100  Sign::kBoth).ToTH2(pot_val);
101 
102  hTotBkg = (TH2*)hNue->Clone(UniqueName().c_str());
103  hTotBkg->Add(hCC);
104  hTotBkg->Add(hBeam);
105 
106 }
107 
108 
109 //------------------------------------------------------------
110 
111 void print_tables(const std::string* cutnames,
112  const std::string pred_basename,
113  const std::string spec_basename,
114  const std::string fname,
115  const double cos_scaling,
116  const int ncuts,
117  bool savePlots=false)
118 {
119  osc::OscCalcSterile* osc = DefaultSterileCalc(3);
120 
121 
122  for(int iv = 0; iv < nvars; iv++){
123  for(int ic = 0; ic < ncuts; ic++){
124 
125 
126  char name[50];
127 
128  sprintf(name, "%s_%s_%s",
129  pred_basename.c_str(),
130  cutnames[ic].c_str(),
131  axisarray_names_opt[iv].c_str());
132 
133  IPrediction* pred =
134  LoadFromFile<IPrediction>(fname, name).release();
135 
136  sprintf(name, "%s_%s_%s",
137  spec_basename.c_str(),
138  cutnames[ic].c_str(),
139  axisarray_names_opt[iv].c_str());
140 
141  Spectrum* spec =
142  LoadFromFile<Spectrum>(fname, name).release();
143 
144  TH1 *sighist, *nuehist, *numuhist, *nutauhist, *bghist;
145  GetSpectra(pred, osc, sighist, nuehist, numuhist, nutauhist, bghist);
146  TH1* coshist = spec->ToTH1(cos_scaling, kBlack, kSolid,
147  ana::EExposureType::kLivetime);
148 
149 
150 
151  if(axisarray_names[iv] == "nusE"){ // nusE/cvn
152  std::cout<<cutnames[ic]<<" & "
153  <<sighist->Integral()<<" & "
154  <<numuhist->Integral()<<" & "
155  <<nuehist->Integral()<<" & "
156  <<nutauhist->Integral()<<" & "
157  <<bghist->Integral()<<" & "
158  <<coshist->Integral()<<" & "
159  <<sighist->Integral()/sqrt(sighist->Integral() +
160  bghist->Integral() +
161  coshist->Integral())
162 
163  << " "<<std::endl;
164  }
165 
166  const Color_t kCosmicBkgdColor = kOrange+7;
167 
168  new TCanvas();
169  gPad->SetFillStyle(0);
170 
171  double maxy = std::max(coshist->GetMaximum(),
172  sighist->GetMaximum());
173 
174  sighist->SetLineColor(kNCBackgroundColor);
175  bghist->SetLineColor(kNumuBackgroundColor);
176  coshist->SetLineColor(kCosmicBkgdColor);
177 
178  if(cutnames[ic].find("nocut") != std::string::npos) gPad->SetLogy();
179  if(cutnames[ic].find("cvnnc") != std::string::npos) gPad->SetLogy();
180  if(cutnames[ic].find("eq17") != std::string::npos) gPad->SetLogy();
181  //gPad->SetLogy();
182 
183  if(axisarray_names[iv] == "nusE"){
184  coshist->GetXaxis()->SetRangeUser(0,12);
185  bghist->GetXaxis()->SetRangeUser(0,12);
186  sighist->GetXaxis()->SetRangeUser(0,12);
187  }
188 
189  if(coshist->Integral() > sighist->Integral()){
190  coshist->Draw("hist");
191  sighist->Draw("hist same");
192  bghist->Draw("hist same");
193  coshist->GetYaxis()->SetTitle("Events / 8.85 #times 10^{20} POT");
194  coshist->GetXaxis()->SetTitle(axisarray_names_xaxis[iv].c_str());
195  CenterTitles(coshist);
196  }
197  else if((bghist->Integral() > sighist->Integral()) &&
198  (bghist->Integral() > coshist->Integral())){
199  bghist->Draw("hist");
200  sighist->Draw("hist same");
201  coshist->Draw("hist same");
202  bghist->GetYaxis()->SetTitle("Events / 8.85 #times 10^{20} POT");
203  bghist->GetXaxis()->SetTitle(axisarray_names_xaxis[iv].c_str());
204  CenterTitles(sighist);
205  }
206  else{
207  sighist->Draw("hist");
208  bghist->Draw("hist same");
209  coshist->Draw("hist same");
210  sighist->GetYaxis()->SetTitle("Events / 8.85 #times 10^{20} POT");
211  sighist->GetXaxis()->SetTitle(axisarray_names_xaxis[iv].c_str());
212  CenterTitles(sighist);
213  }
214 
215  TLegend *leg;
216  if(axisarray_names[iv] == "nusE") leg = new TLegend(0.5, 0.65, 0.8, 0.85);
217  if(axisarray_names[iv] == "cvn") leg = new TLegend(0.3, 0.65, 0.7, 0.85);
218  if(axisarray_names[iv] == "hits") leg = new TLegend(0.5, 0.65, 0.85, 0.85);
219  else leg = new TLegend(0.5, 0.65, 0.8, 0.85);
220  leg->SetBorderSize(0);
221  leg->SetFillColor(kWhite);
222  leg->SetFillStyle(0);
223  leg->SetFillStyle(0);
224  leg->SetTextFont(42);
225  leg->SetTextSize(0.037);
226  leg->AddEntry(sighist, "NC Signal", "l");
227  leg->AddEntry(bghist, "CC Background", "l");
228  leg->AddEntry(coshist,"Cosmic Background", "l");
229  leg->Draw();
230  //Preliminary();
231 
232  if(savePlots){
233  gPad->SaveAs(("plots/final_"+spec_basename+cutnames[ic]+"_"+axisarray_names[iv]+".png").c_str());
234  gPad->SaveAs(("plots/final_"+spec_basename+cutnames[ic]+"_"+axisarray_names[iv]+".pdf").c_str());
235  gPad->SaveAs(("plots/final_"+spec_basename+cutnames[ic]+"_"+axisarray_names[iv]+".eps").c_str());
236  }
237 
238 
239  new TCanvas();
240  gPad->SetFillStyle(0);
241 
242  if(cutnames[ic].find("nocut") != std::string::npos) gPad->SetLogy();
243  if(cutnames[ic].find("cvnnc") != std::string::npos) gPad->SetLogy();
244 
245  TH1 *cosclone, *bgclone, *sigclone;
246  if(sighist->Integral() > coshist->Integral()){
247  cosclone = (TH1*)coshist->Clone();
248  bgclone = (TH1*)cosclone->Clone();
249  bgclone->Add(bghist);
250  sigclone = (TH1*)bgclone->Clone();
251  sigclone->Add(sighist);
252  }
253  else{
254  bgclone = (TH1*)bghist->Clone();
255  sigclone = (TH1*)bgclone->Clone();
256  sigclone->Add(sighist);
257  cosclone = (TH1*)sigclone->Clone();
258  cosclone->Add(coshist);
259  }
260 
261  sigclone->SetLineColor(kNCBackgroundColor);
262  sigclone->SetFillColor(kNCBackgroundColor);
263  FillWithDimColor(sigclone);
264 
265  bgclone->SetLineColor(kNumuBackgroundColor);
266  bgclone->SetFillColor(kNumuBackgroundColor);
267  FillWithDimColor(bgclone);
268  cosclone->SetLineColor(kCosmicBkgdColor);
269  cosclone->SetFillColor(kCosmicBkgdColor);
270  FillWithDimColor(cosclone);
271 
272  if(axisarray_names[iv] == "nusE"){
273  cosclone->GetXaxis()->SetRangeUser(0,12);
274  bgclone->GetXaxis()->SetRangeUser(0,12);
275  sigclone->GetXaxis()->SetRangeUser(0,12);
276  }
277 
278  if(sighist->Integral() > coshist->Integral()){
279  sigclone->Draw("hist");
280  bgclone->Draw("hist same");
281  cosclone->Draw("hist same");
282  sigclone->Draw("histsameaxis");
283  cosclone->GetYaxis()->SetTitle("Events / 8.85 #times 10^{20} POT");
284  cosclone->GetXaxis()->SetTitle(axisarray_names_xaxis[iv].c_str());
285  sigclone->GetYaxis()->SetTitle("Events / 8.85 #times 10^{20} POT");
286  sigclone->GetXaxis()->SetTitle(axisarray_names_xaxis[iv].c_str());
287  bgclone->GetYaxis()->SetTitle("Events / 8.85 #times 10^{20} POT");
288  bgclone->GetXaxis()->SetTitle(axisarray_names_xaxis[iv].c_str());
289 
290  CenterTitles(cosclone);
291  CenterTitles(bgclone);
292  CenterTitles(sigclone);
293  }
294  else{
295  cosclone->Draw("hist");
296  sigclone->Draw("hist same");
297  bgclone->Draw("hist same");
298  sigclone->GetYaxis()->SetTitle("Events / 8.85 #times 10^{20} POT");
299  sigclone->GetXaxis()->SetTitle(axisarray_names_xaxis[iv].c_str());
300  bgclone->GetYaxis()->SetTitle("Events / 8.85 #times 10^{20} POT");
301  bgclone->GetXaxis()->SetTitle(axisarray_names_xaxis[iv].c_str());
302  cosclone->GetYaxis()->SetTitle("Events / 8.85 #times 10^{20} POT");
303  bgclone->GetXaxis()->SetTitle(axisarray_names_xaxis[iv].c_str());
304 
305  CenterTitles(sigclone);
306  CenterTitles(cosclone);
307  CenterTitles(bgclone);
308  }
309  TLegend *legStack;
310  if(axisarray_names[iv] == "nusE") legStack = new TLegend(0.5, 0.65, 0.8, 0.85);
311  if(axisarray_names[iv] == "cvn") legStack = new TLegend(0.3, 0.65, 0.7, 0.85);
312  if(axisarray_names[iv] == "hits") legStack = new TLegend(0.5, 0.65, 0.85, 0.85);
313  else legStack = new TLegend(0.5, 0.65, 0.8, 0.85);
314  legStack->SetBorderSize(0);
315  legStack->SetFillColor(kWhite);
316  legStack->SetFillStyle(0);
317  legStack->SetFillStyle(0);
318  legStack->SetTextFont(42);
319  legStack->SetTextSize(0.037);
320 
321  legStack->AddEntry(sigclone, "NC Signal", "f");
322  legStack->AddEntry(bgclone, "CC Background", "f");
323  legStack->AddEntry(cosclone,"Cosmic Background", "f");
324  legStack->Draw();
325 
326  //Preliminary();
327 
328  if(savePlots){
329  gPad->SaveAs(("plots/final_stack_"+spec_basename+cutnames[ic]+"_"+axisarray_names[iv]+".png").c_str());
330  gPad->SaveAs(("plots/final_stack_"+spec_basename+cutnames[ic]+"_"+axisarray_names[iv]+".root").c_str());
331  gPad->SaveAs(("plots/final_stack_"+spec_basename+cutnames[ic]+"_"+axisarray_names[iv]+".pdf").c_str());
332  }
333  }// end loop over vars
334  }// end loop over cuts
335 
336  // for(int ic = 0; ic < ncuts; ic++){
337 
338  // char name[50];
339 
340  // sprintf(name, "%s_%s_%s",
341  // pred_basename.c_str(),
342  // cutnames[ic].c_str(),
343  // "2Dptp");
344 
345  // std::cout << "name: "<< name << std::endl;
346  // IPrediction* pred =
347  // LoadFromFile<IPrediction>(fname, name).release();
348 
349  // sprintf(name, "%s_%s_%s",
350  // spec_basename.c_str(),
351  // cutnames[ic].c_str(),
352  // "2Dptp");
353  // std::cout << "name: "<< name << std::endl;
354  // Spectrum* spec =
355  // LoadFromFile<Spectrum>(fname, name).release();
356 
357  // TH2 *sighist, *bghist;
358  // GetSpectra2D(pred, &osc, sighist, bghist);
359  // std::cout << "here: "<< std::endl;
360  // TH2* coshist = spec->ToTH2(cos_scaling,
361  // ana::EExposureType::kLivetime);
362 
363  // new TCanvas;
364  // coshist->Draw("colz");
365  // gPad->SaveAs(("plots/"+spec_basename+cutnames[ic]+"2Dcos.png").c_str());
366  // //gPad->SaveAs(("plots/"+spec_basename+cutnames[ic]+"2Dcos.root").c_str());
367 
368  // new TCanvas;
369  // sighist->Draw("colz");
370  // gPad->SaveAs(("plots/"+spec_basename+cutnames[ic]+"2Dsig.png").c_str());
371  // //gPad->SaveAs(("plots/"+spec_basename+cutnames[ic]+"2Dsig.root").c_str());
372 
373  // new TCanvas;
374  // bghist->Draw("colz");
375  // gPad->SaveAs(("plots/"+spec_basename+cutnames[ic]+"2Dbg.png").c_str());
376  // //gPad->SaveAs(("plots/"+spec_basename+cutnames[ic]+"2Dbg.root").c_str());
377  // }
378 
379 }
380 
381 //------------------------------------------------------------
382 
383 void print_nus17_fd_cut_tables(bool savePlots=false)
384 {
385 
386  //std::string fname = "/nova/ana/users/gsdavies/nus/nus17/FD_cutTable_extendedE.root";
387  std::string fname = "/nova/ana/users/gsdavies/nus/nus17/nus17FDCutTablePlots.root";
388  //std::string fname = "mcv1_cutTable.root";
389 
390  // cosmic livetime scaling
391  TFile fin(fname.c_str());
392  if(fin.IsZombie()){
393  std::cerr << "File: " << fname << " does not exist. Check the name." << std::endl;
394  exit(0);
395  }
396 
397 
398  TH1D* time = (TH1D*)gDirectory->Get("time_cos");
399  double scaling = 440.;
400 
401  char name[15];
402 
403 
404  std::cout<<"\n--------------------------------------"
405  <<"\n---------------N-1 Table -------------\n";
406 
407  print_tables(cutnames_nminus1,
408  "pred_nminus1",
409  "spec_nminus1",
410  fname, scaling, ncuts_nminus1, savePlots);
411 
412  std::cout<<"\n------------------------------------------"
413  <<"\n--------------Cut Flow Table -------------\n";
414 
415  print_tables(cutnames,
416  "pred",
417  "spec",
418  fname, scaling, ncuts, savePlots);
419 
420 }
std::max
T max(const caf::Proxy< T > &a, T b)
Definition: BasicTypesProxy.h:414
fin
TString fin
Definition: Style.C:24
name
const XML_Char * name
Definition: expat.h:151
plot_validation_datamc.fname
fname
Definition: plot_validation_datamc.py:180
kOrange
enum BeamMode kOrange
plot_validation_datamc.pred
pred
Definition: plot_validation_datamc.py:175
ana::Spectrum::ToTH2
TH2 * ToTH2(double exposure, EExposureType expotype=kPOT, EBinType bintype=kBinContent) const
Spectrum must be 2D to obtain TH2.
Definition: Spectrum.cxx:165
ana
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
ana::IPrediction::PredictComponent
virtual Spectrum PredictComponent(osc::IOscCalc *calc, Flavors::Flavors_t flav, Current::Current_t curr, Sign::Sign_t sign) const =0
axisarray_names_opt
const std::string axisarray_names_opt[nvars]
Definition: nus17_fd_arrays.h:187
ana::Flavors::kNuMuToNuE
(&#39; appearance&#39;)
Definition: IPrediction.h:18
ana::Spectrum::ToTH1
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
maxy
double maxy
Definition: plot_datapredictions.C:15
ana::Flavors::kNuEToNuE
(&#39;beam &#39;)
Definition: IPrediction.h:15
std::sqrt
T sqrt(T number)
Definition: d0nt_math.hpp:156
osc::_IOscCalc
General interface to oscillation calculators.
Definition: StanTypedefs.h:22
osc::OscCalcSterile
Adapt the PMNS_Sterile calculator to standard interface.
Definition: OscCalcSterile.h:23
om::cerr
OStream cerr
Definition: OStream.cxx:7
ana::kAna2017POT
const double kAna2017POT
Definition: Exposures.h:174
pot_val
const double pot_val
Definition: print_nus17_fd_cut_tables.C:41
livetime
const double livetime
Definition: print_nus17_fd_cut_tables.C:42
ana::CenterTitles
void CenterTitles(TH1 *histo)
Definition: Plots.cxx:1481
print_nus17_fd_cut_tables
void print_nus17_fd_cut_tables(bool savePlots=false)
Definition: print_nus17_fd_cut_tables.C:383
calc
osc::OscCalcDumb calc
Definition: make_nueFDprediction_kinematics.C:70
ana::kNumuBackgroundColor
const Color_t kNumuBackgroundColor
Definition: Style.h:30
ana::Spectrum
Representation of a spectrum in any variable, with associated POT.
Definition: Spectrum.h:40
ana::kNueSignalColor
const Color_t kNueSignalColor
Definition: Style.h:19
scaling
void scaling(TH1D *hIn, const double shape_scale)
Definition: Calculate_covariances.C:431
nus17_fd_arrays.h
ana::Current::kCC
Charged-current interactions.
Definition: IPrediction.h:39
ana::DefaultSterileCalc
osc::OscCalcSterile * DefaultSterileCalc(int nflavors)
Create a sterile calculator with default assumptions for all parameters.
Definition: Calcs.cxx:97
ana::kAna2017Livetime
const double kAna2017Livetime
Definition: Exposures.h:200
cutnames
const std::string cutnames[ncuts]
Definition: nus17_fd_arrays.h:41
GetSpectra2D
void GetSpectra2D(IPrediction *pred, osc::IOscCalc *calc, TH2 *&hSig, TH2 *&hTotBkg)
Definition: print_nus17_fd_cut_tables.C:79
cutnames_nminus1
const std::string cutnames_nminus1[ncuts_nminus1]
Definition: nus17_fd_arrays.h:56
print_tables
void print_tables(const std::string *cutnames, const std::string pred_basename, const std::string spec_basename, const std::string fname, const double cos_scaling, const int ncuts, bool savePlots=false)
Definition: print_nus17_fd_cut_tables.C:111
osc
Oscillation probability calculators.
Definition: Calcs.h:5
om::cout
OStream cout
Definition: OStream.cxx:6
ncuts_nminus1
const int ncuts_nminus1
Definition: nus17_fd_arrays.h:54
nvars
const int nvars
Definition: nus17_fd_arrays.h:168
exit
exit(0)
GetSpectra
void GetSpectra(IPrediction *pred, osc::IOscCalc *calc, TH1 *&hSig, TH1 *&hNue, TH1 *&hCC, TH1 *&hTau, TH1 *&hTotBkg)
Definition: print_nus17_fd_cut_tables.C:45
ana::Current::kNC
Neutral-current interactions.
Definition: IPrediction.h:40
ana::Sign::kBoth
Both neutrinos and antineutrinos.
Definition: IPrediction.h:52
ana::IPrediction
Standard interface to all prediction techniques.
Definition: IPrediction.h:57
ana::kTotalMCColor
const Color_t kTotalMCColor
Definition: Style.h:16
ana::Flavors::kAll
All neutrinos, any flavor.
Definition: IPrediction.h:26
ncuts
const int ncuts
Definition: nus17_fd_arrays.h:9
ana::kNCBackgroundColor
const Color_t kNCBackgroundColor
Definition: Style.h:22
ana::UniqueName
std::string UniqueName()
Return a different string each time, for creating histograms.
Definition: Utilities.cxx:29
ana::FillWithDimColor
void FillWithDimColor(TH1 *h, bool usealpha, float dim)
Definition: NuePlotStyle.cxx:55
axisarray_names
const std::string axisarray_names[nvars]
Definition: nus17_fd_arrays.h:178
string
enum BeamMode string
axisarray_names_xaxis
const std::string axisarray_names_xaxis[nvars]
Definition: nus17_fd_arrays.h:211
Generated on Mon Nov 30 2020 05:14:58 for NOvA by   doxygen 1.8.11