plot_nominal_spectra.C
Go to the documentation of this file.
7 #include "CAFAna/Core/Spectrum.h"
15 #include "OscLib/IOscCalc.h"
16 
17 #include "TFile.h"
18 #include "TH1.h"
19 #include "TCanvas.h"
20 #include "TGaxis.h"
21 #include "TGraph.h"
22 #include "TLatex.h"
23 #include "TLegend.h"
24 #include "TLine.h"
25 #include "TSystem.h"
26 #include <fstream>
27 #include <iostream>
28 #include <iomanip>
29 #include <string>
30 #include <cstddef>
31 
32 using namespace ana;
33 
34 void PlotNuePredictionFourBins (const IPrediction* pred, osc::IOscCalc* calc, double fdpot);
35 void TableNuePredictionFourBins(const IPrediction* pred, osc::IOscCalc* calc, double fdpot);
36 void PlotNueDecompFourBins (MichelDecomp* combo = 0, ProportionalDecomp * prop = 0,
37  bool mccomponly = false, bool overlay = false);
38 void TableNueNDDataMC (MichelDecomp* combo);
40 
42 
43 
44  for(std::string decomp :{"noextrap","prop","combo"}){
46 
47  // FD plots
48  auto calc = DefaultOscCalc();
49  SetFakeCalc(calc, 0.402, 2.7e-3, 1.5);
50 
51  double fdpot = kAna2017POT;
52 
54  gPad->Print((TString)"plot_nue_fd_"+decomp+"_rock_spectrum.pdf");
55  std::cout << decomp ;
57  }
58  // ND Plots
59  //Add prop to the mix
60  auto pred_file = new TFile("/nova/ana/nu_e_ana/Ana2017/Predictions/provisional/pred_nue_reduced_v4.root", "read");
61  auto xp_dir = pred_file->GetDirectory((TString) "pred_" + "xp_combo" +
62  "_Nominal/noShift/");
63  auto xp_dir2 = pred_file->GetDirectory((TString) "pred_" + "xp_prop" +
64  "_Nominal/noShift/");
65  TString deco_dir ="extrap/EEextrap/Decomp";
66  auto combo = MichelDecomp::LoadFrom (xp_dir, deco_dir).release();
67  auto prop = ProportionalDecomp::LoadFrom(xp_dir2, deco_dir).release();
68 
69 
70  PlotNueDecompFourBins (combo,0,true,false);
71  gPad->Print("plot_nue_nd_nominal_mconly.pdf");
72 
73  PlotNueDecompFourBins (combo,0,false,true);
74  gPad->Print("plot_nue_nd_nominal_combo_vs_mc.pdf");
75 
76  PlotNueDecompFourBins (combo,prop,false,true);
77  gPad->Print("plot_nue_nd_nominal_prop_vs_mc.pdf");
78 
79 
80  TableNueNDDataMC(combo);
81 
82  TableNueNDComponents(combo,0);
83 
84  TableNueNDComponents(combo,prop);
85 
86 }
87 void PrintTableFourBins(std::vector <TH1D*> mcnom, std::vector <TString> labels){
88 
89  assert(mcnom.size()==labels.size());
90 
91  int N = 4;
92  int Nbins = mcnom[0]->GetNbinsX();
93  int binNbins = Nbins / N;
94  for(TString b:{"Bin1","Bin2","Bin3","Bin4","AllBins"}){
95  std::cout << "\t & \t " << b;
96  }
97  std::cout << "\t \\\\ \\\\hline \n";
98  for(int i = 0; i < (int) mcnom.size(); ++ i){
99  std::cout << labels[i] << "\t & \t";
100  for(int n = 0; n < 5; ++n){
101  int bin1 = 1 + (n * binNbins);
102  int bin2 = (n+1) * binNbins;
103  if(n==4) {bin1 = 1; bin2 = Nbins;}
104  std::cout << mcnom[i]->Integral(bin1,bin2)
105  << (n == 4 ? "\t \\\\\n" : "\t & \t ");
106  }//end bins
107  }//end rows
108  std::cout << "\n\n";
109 }
110 
111 
112 
113 std::vector <TH1D*> GetNueNDComponents(IDecomp* decomp, double pot)
114 {
115  std::vector <TH1D*> ret;
116  ret.push_back( (decomp->NueComponent() + decomp->AntiNueComponent() +
117  decomp->NumuComponent() + decomp->AntiNumuComponent() +
118  decomp->NCTotalComponent() ). ToTH1(pot,kTotalMCColor) );
119  ret.push_back( (decomp->NueComponent() + decomp->AntiNueComponent()).
120  ToTH1(pot,kBeamNueBackgroundColor) );
121  ret.push_back( (decomp->NCTotalComponent()).
122  ToTH1(pot,kNCBackgroundColor) );
123  ret.push_back( (decomp->NumuComponent() + decomp->AntiNumuComponent()).
124  ToTH1(pot,kNumuBackgroundColor) );
125  return ret;
126 }
127 
128 std::vector <TH1D*> GetNueNDComponentsMC(MichelDecomp* decomp, double pot)
129 {
130  std::vector <TH1D*> ret;
131  ret.push_back( (decomp->MC_NueComponent() + decomp->MC_AntiNueComponent() +
132  decomp->MC_NumuComponent() + decomp->MC_AntiNumuComponent() +
133  decomp->MC_NCTotalComponent() ). ToTH1(pot,kTotalMCColor) );
134  ret.push_back( (decomp->MC_NueComponent() + decomp->MC_AntiNueComponent()).
135  ToTH1(pot,kBeamNueBackgroundColor) );
136  ret.push_back( (decomp->MC_NCTotalComponent()).
137  ToTH1(pot,kNCBackgroundColor) );
138  ret.push_back( (decomp->MC_NumuComponent() + decomp->MC_AntiNumuComponent()).
139  ToTH1(pot,kNumuBackgroundColor) );
140  return ret;
141 }
142 
143 
145 {
146  auto data = combo->Data_Component();
147  auto ndpot = data.POT();
148  auto hdata = data.ToTH1(ndpot);
149  auto comp = GetNueNDComponentsMC (combo, ndpot);
150 
151  std::cout << "Data vs uncorrected MC\n";
152  comp.insert(comp.begin(),hdata);
153  std::vector <TString> labels = {"Data","Tot MC","Beam","NC","Numu"};
154  std::cout << std::setprecision(0) << std::fixed;
155  PrintTableFourBins(comp,labels);
156 }
157 
159 {
160  auto data = combo->Data_Component();
161  auto ndpot = data.POT();
162  std::vector <TH1D*> comp;
163  if(prop) comp = GetNueNDComponents (prop, ndpot);
164  else comp = GetNueNDComponents (combo, ndpot);
165  std::vector <TString> labels = {"Tot MC","Beam","NC","Numu"};
166  std::cout << std::setprecision(0) << std::fixed;
167  std::cout << (prop?"prop":"combo");
168  PrintTableFourBins(comp,labels);
169 }
170 
172  bool mccomponly, bool overlay)
173 {
174  auto data = combo->Data_Component();
175  auto ndpot = data.POT();
176  auto hdata = data.ToTH1(ndpot);
177  hdata->SetMarkerStyle(kFullCircle);
178  hdata->Draw("e1");
179 
180  std::vector <TH1D*> comp_decomp;
181  if(prop) comp_decomp = GetNueNDComponents(prop, ndpot);
182  else comp_decomp = GetNueNDComponents(combo, ndpot);
183 
184  if(!mccomponly) for (auto h:comp_decomp) h->Draw("hist same");
185 
186  std::vector <TH1D*> comp_mc;
187  if( mccomponly || overlay) comp_mc = GetNueNDComponentsMC(combo, ndpot);
188  if( overlay ) for (auto h:comp_mc) h->SetLineStyle(7);
189  if( mccomponly || overlay) for (auto h:comp_mc) h->Draw("hist same");
190  hdata->Draw("e1 same");
191 
192  hdata->GetYaxis()->SetTitle("Events / 0.5 GeV bin");
193  hdata->SetMaximum(1.2*hdata->GetMaximum());
194  Nue2017FourBinLabels(0.98, 0.8 * hdata->GetXaxis()->GetLabelSize());
196  Nue2017FourBinAxis(hdata);
197 
198  auto leg = new TLegend (0.70,0.2,0.89,0.75);
199  leg->SetTextSize(hdata->GetYaxis()->GetLabelSize());
200  leg->SetHeader(TString::Format("#splitline{NOvA ND}{%.2f#times10^{20} POT}",
201  ndpot/1.E20));
202  leg->AddEntry(hdata, "ND Data","epl");
203  leg->AddEntry(comp_decomp[0], "Total","l");
204  leg->AddEntry(comp_decomp[1], "Beam #nu_{e} CC","l");
205  leg->AddEntry(comp_decomp[2], "NC","l");
206  leg->AddEntry(comp_decomp[3], "#nu_{#mu} CC","l");
207  if(overlay){
208  auto dummy = new TGraph();
209  dummy->SetLineWidth(2);
210  dummy->SetLineColor(kGray+2);
211  dummy->SetLineStyle(7);
212  leg->AddEntry(dummy, "Uncorr. MC", "l");
213  }
214  leg->Draw();
215 }
216 
218 
219  auto htot = pred->Predict(calc).ToTH1(fdpot, kTotalMCColor);
221  ToTH1(fdpot, kNueSignalColor);
222  auto hbeam = pred->PredictComponent(calc, Flavors::kNuEToNuE, Current::kCC, Sign::kBoth).
223  ToTH1(fdpot, kBeamNueBackgroundColor);
224  auto hnc = pred->PredictComponent(calc, Flavors::kAll, Current::kNC, Sign::kBoth).
225  ToTH1(fdpot, kNCBackgroundColor);
226  auto hnumu = pred->PredictComponent(calc, Flavors::kAllNuMu, Current::kCC, Sign::kBoth).
227  ToTH1(fdpot, kNumuBackgroundColor);
228 
229  htot->Draw("hist");
230  for(auto h:{hsig,hbeam,hnc,hnumu})h->Draw("hist same");
231 
232  htot->GetYaxis()->SetTitle("Events / 0.5 GeV bin");
233  htot->SetMaximum(1.2*htot->GetMaximum());
234  Nue2017FourBinLabels(0.98, 0.8 * htot->GetXaxis()->GetLabelSize());
236  Nue2017FourBinAxis(htot);
237 
238  auto leg = new TLegend (0.12,0.4,0.35,0.75);
239  leg->SetTextSize(htot->GetYaxis()->GetLabelSize());
240  leg->SetHeader(TString::Format("NOvA FD, %.2f #times10^{20} POT eq.",
241  kAna2017FullDetEquivPOT/1.E20));
242  leg->AddEntry(htot, "Nominal Osc. Pred.","l");
243  leg->AddEntry(hsig, "Signal","l");
244  leg->AddEntry(hbeam, "Beam #nu_{e} CC","l");
245  leg->AddEntry(hnc, "NC","l");
246  leg->AddEntry(hnumu, "#nu_{#mu} CC","l");
247  leg->Draw();
248 }
249 
251 {
252  auto htot = pred->Predict(calc).ToTH1(fdpot, kTotalMCColor);
254  ToTH1(fdpot, kNueSignalColor);
255  auto hbkg = (TH1D*) htot->Clone(UniqueName().c_str());
256  hbkg->Add(hsig,-1);
257  auto hbeam = pred->PredictComponent(calc, Flavors::kNuEToNuE, Current::kCC, Sign::kBoth).
258  ToTH1(fdpot, kBeamNueBackgroundColor);
259  auto hnc = pred->PredictComponent(calc, Flavors::kAll, Current::kNC, Sign::kBoth).
260  ToTH1(fdpot, kNCBackgroundColor);
261  auto hnumu = pred->PredictComponent(calc, Flavors::kAllNuMu, Current::kCC, Sign::kBoth).
262  ToTH1(fdpot, kNumuBackgroundColor);
263  auto htau = pred->PredictComponent(calc, Flavors::kAllNuTau, Current::kCC, Sign::kBoth).
264  ToTH1(fdpot, kNumuBackgroundColor);
265 
266 
267  std::vector <TH1D*> mcnom = {htot,hsig,hbkg,hbeam,hnc,hnumu, htau};
268  std::vector <TString> labels = {"Total", "Signal","Bkg.","Beam","NC","Numu","Tau"};
269 
270  std::cout << std::setprecision(3) ;
271  PrintTableFourBins(mcnom,labels);
272 }
void Nue2017FourBinLabels(const double yNDC, const double textsize, const int color, const bool merged)
void Nue2017FourBinAxis(TH1 *axes, bool drawLabels, bool merged)
virtual Spectrum Data_Component() const override
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
virtual Spectrum PredictComponent(osc::IOscCalc *calc, Flavors::Flavors_t flav, Current::Current_t curr, Sign::Sign_t sign) const =0
Spectrum MC_AntiNumuComponent() const override
std::vector< TH1D * > GetNueNDComponentsMC(MichelDecomp *decomp, double pot)
virtual Spectrum AntiNueComponent() const =0
void TableNueNDComponents(MichelDecomp *combo, ProportionalDecomp *prop=0)
(&#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:209
void TableNuePredictionFourBins(const IPrediction *pred, osc::IOscCalc *calc, double fdpot)
virtual Spectrum NumuComponent() const =0
virtual Spectrum Predict(osc::IOscCalc *calc) const =0
void plot_nominal_spectra()
(&#39;beam &#39;)
Definition: IPrediction.h:15
const IPrediction * GetNuePrediction2017(std::string decomp, osc::IOscCalc *calc, bool corrSysts, bool GetFromUPS=false)
const double kAna2017POT
Definition: Exposures.h:174
TH1F * hsig
Definition: Xdiff_gwt.C:59
osc::OscCalcDumb calc
osc::IOscCalcAdjustable * DefaultOscCalc()
Create a new calculator with default assumptions for all parameters.
Definition: Calcs.cxx:49
const Color_t kNumuBackgroundColor
Definition: Style.h:30
static std::unique_ptr< ProportionalDecomp > LoadFrom(TDirectory *dir, const std::string &name)
const Color_t kNueSignalColor
Definition: Style.h:19
static std::unique_ptr< MichelDecomp > LoadFrom(TDirectory *dir, const std::string &name)
std::vector< TH1D * > GetNueNDComponents(IDecomp *decomp, double pot)
const XML_Char const XML_Char * data
Definition: expat.h:268
void PrintTableFourBins(std::vector< TH1D * > mcnom, std::vector< TString > labels)
Charged-current interactions.
Definition: IPrediction.h:39
Spectrum MC_NCTotalComponent() const override
const Color_t kBeamNueBackgroundColor
Definition: Style.h:24
#define pot
void PlotNueDecompFourBins(MichelDecomp *combo=0, ProportionalDecomp *prop=0, bool mccomponly=false, bool overlay=false)
virtual Spectrum AntiNumuComponent() const =0
void SetFakeCalc(osc::IOscCalcAdjustable *calc, double ssth23=-5, double dmsq32=-5, double dCP_Pi=-5)
double POT() const
Definition: Spectrum.h:231
void PlotNuePredictionFourBins(const IPrediction *pred, osc::IOscCalc *calc, double fdpot)
OStream cout
Definition: OStream.cxx:6
Spectrum MC_NumuComponent() const override
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
Splits Data proportionally according to MC.
const double kAna2017FullDetEquivPOT
Definition: Exposures.h:195
virtual Spectrum NCTotalComponent() const
Definition: IDecomp.h:18
Standard interface to all decomposition techniques.
Definition: IDecomp.h:13
const hit & b
Definition: hits.cxx:21
Neutral-current interactions.
Definition: IPrediction.h:40
assert(nhit_max >=nhit_nbins)
Both neutrinos and antineutrinos.
Definition: IPrediction.h:52
Standard interface to all prediction techniques.
Definition: IPrediction.h:57
void Format(TGraph *gr, int lcol, int lsty, int lwid, int mcol, int msty, double msiz)
Definition: Style.cxx:154
Spectrum MC_NueComponent() const override
const Color_t kTotalMCColor
Definition: Style.h:16
void Nue2017FourBinDivisions(const int color, const int style)
All neutrinos, any flavor.
Definition: IPrediction.h:26
Float_t e
Definition: plot.C:35
void TableNueNDDataMC(MichelDecomp *combo)
Spectrum MC_AntiNueComponent() const override
const Color_t kNCBackgroundColor
Definition: Style.h:22
std::string UniqueName()
Return a different string each time, for creating histograms.
Definition: Utilities.cxx:30
virtual Spectrum NueComponent() const =0
TH1F * hbkg
Definition: Xdiff_gwt.C:58