plot_3NDvsFD_FHC.C
Go to the documentation of this file.
1 
2 #include "CAFAna/Analysis/Exposures.h"
3 #include "CAFAna/Analysis/Calcs.h"
4 #include "CAFAna/Core/Binning.h"
5 #include "CAFAna/Core/Cut.h"
6 #include "CAFAna/Core/HistAxis.h"
7 #include "CAFAna/Core/LoadFromFile.h"
8 #include "CAFAna/Core/Spectrum.h"
9 #include "CAFAna/Core/SpectrumLoader.h"
10 #include "CAFAna/Core/SpectrumLoaderBase.h"
11 #include "CAFAna/Cuts/Cuts.h"
12 #include "3FlavorAna/Cuts/NumuCuts.h"
13 #include "CAFAna/Cuts/NueCutsSecondAna.h"
14 #include "CAFAna/Cuts/SpillCuts.h"
15 #include "CAFAna/Decomp/ProportionalDecomp.h"
16 #include "3FlavorAna/Decomp/NumuDecomp.h"
17 #include "CAFAna/Extrap/ModularExtrap.h"
18 #include "CAFAna/Prediction/PredictionExtrap.h"
19 #include "CAFAna/Prediction/PredictionNoExtrap.h"
20 #include "CAFAna/Vars/GenieWeights.h"
21 #include "3FlavorAna/Vars/NueVars.h"
22 #include "CAFAna/Vars/Vars.h"
23 
24 #include "OscLib/OscCalcPMNSOpt.h"
25 #include "OscLib/OscCalcDumb.h"
26 #include "CAFAna/Analysis/Style.h"
27 
28 //#include "CAFAna/nue/SecondAna/draw_plots_util.h"
29 
30 using namespace ana;
31 
32 #include "TFile.h"
33 #include "TH1.h"
34 #include "TH2.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 
46 //make a pretty legend
47 void Legend(double x0, double y0, double x1, double y1)
48 {
49  // Doesn't handle log-y well
50  // TLegend* leg = AutoPlaceLegend(x1-x0, y1-y0);
51  TLegend* leg = new TLegend(x0, y0, x1, y1);
52  leg->SetFillStyle(0);
53  //leg->SetTextSize(2/30.);
54  leg->SetHeader("#splitline{#bf{Neutrino Beam}}{CORE sample: cos #theta reweighted}");
55  TH1* dummy = new TH1F("", "", 1, 0, 1);
56  TH1* dummyFill = new TH1F("", "", 1, 0, 1);
57  dummyFill->SetLineColor(kTotalMCColor);
58  dummyFill->SetFillColor(kTotalMCErrorBandColor);
59 
60  dummy->SetMarkerStyle(kFullCircle);
61 
62  dummy->SetLineColor(kGray+2);
63  dummy->SetMarkerColor(kGray+2);
64  dummy->SetLineColor(kTotalMCColor);
65  leg->AddEntry(dummy->Clone(), "ND #nu_{#mu} MC reweighted", "l");
66  dummy->SetLineColor(1);
67  leg->AddEntry(dummy->Clone(), "ND #nu_{#mu} MC nominal", "l");
68 
69  dummy->SetLineColor(kNueSignalColor);
70  leg->AddEntry(dummy->Clone(), "FD MC #nu_{e} + WS", "l");
71 
72  leg->Draw();
73 }
74 
75 void plot_3NDvsFD_FHC(const std::string weight_name, const std::string fname_NDnon, const std::string fname_ND, const std::string fname_FD)
76 {
77  gStyle->SetPadTickX(1);
78  gStyle->SetPadTickY(0);
79 
80  TH1::AddDirectory(kFALSE);
81 
82  struct HistDef
83  {
84  std::string name;
85  };
86 
87  //which variables do you actually want to look at?
88  const int kNumVars = 5;
89  const HistDef defs[kNumVars] = {
90  "trueQ2",
91  "recoQ2",
92  "trueW2",
93  "PtP",
94  "CosNumi"
95  };
96 
97  //which cut tiers do you want to look at?
98  const int kNumSels = 1;
99  const std::string selNamesND[kNumSels] = {
100  "numuND"
101  };
102 
103  const std::string selNamesFD[kNumSels] = {
104  "CORE"
105  };
106 
107  IPrediction* predsND[kNumSels][kNumVars];
108  IPrediction* predsNDnon[kNumSels][kNumVars];
109 
110  osc::NoOscillations noosc;
111  osc::OscCalcDumb mycalc;
112  //osc::OscCalcPMNSOpt mycalc;
113  //ResetOscCalcToDefault(&mycalc);
114 
115  //mycalc.SetDmsq32(+fabs(mycalc.GetDmsq32()));
116  //mycalc.SetdCP(3*TMath::Pi()/2);
117  //mycalc.SetTh23(asin(sqrt(0.404)));
118 
119  //pick up spectra/prediction objects from the input root file
120  for(int selIdx = 0; selIdx < kNumSels; ++selIdx){
121  const char* selName = selNamesND[selIdx].c_str();
122  const char* selNameFD = selNamesFD[selIdx].c_str();
123  const char* wName = weight_name.c_str();
124 
125  for(int varIdx = 0; varIdx < kNumVars; ++varIdx){
126  const char* varName = defs[varIdx].name.c_str();
127 
128  std::cout<<TString::Format("%s/spect_%s_%s", selName, wName, varName).Data()<<std::endl;
129 
130  predsND[selIdx][varIdx] = LoadFromFile<IPrediction>(fname_ND, TString::Format("%s/pred_%s_%s", selNameFD, wName, varName).Data()).release();
131  predsNDnon[selIdx][varIdx] = LoadFromFile<IPrediction>(fname_NDnon, TString::Format("%s/pred_noweight_%s", selNameFD, varName).Data()).release();
132 
133  if(selIdx == 0 && varIdx == 0)
134  std::cout << "ND MC POT: "
135  << predsNDnon[selIdx][varIdx]->Predict(&noosc).POT() << std::endl;
136 
137  } // end for varIdx
138  } // end for selIdx
139 
140 
141  TFile * file_no = new TFile (fname_FD.c_str(),"READ");
142  TDirectory * dpred_no[kNumSels][kNumVars];
143  for(int selIdx = 0; selIdx < kNumSels; ++selIdx){
144  const std::string selName = selNamesFD[selIdx];
145  for(int varIdx = 0; varIdx < kNumVars; ++varIdx){
146  const std::string varName = defs[varIdx].name;
147  if(varName == "recoQ2") {
148  dpred_no[selIdx][varIdx] = file_no->GetDirectory(("prediction/nue_pred_noextrap_"+selName+"_trueQ2").c_str());
149  } else if(varName == "CosNumi"){
150  dpred_no[selIdx][varIdx] = file_no->GetDirectory(("prediction/nue_pred_noextrap_"+selName+"_CosTheta").c_str());
151  } else {
152  dpred_no[selIdx][varIdx] = file_no->GetDirectory(("prediction/nue_pred_noextrap_"+selName+"_"+varName).c_str());
153  }
154  }
155  }
156 
157  TH1* hFullNDHists[kNumVars-1];
158  for(int selIdx = 0; selIdx < kNumSels; ++selIdx){
159  const std::string selName = selNamesFD[selIdx];
160  for(int varIdx = 0; varIdx < kNumVars; ++varIdx){
161 
162  const std::string varName = defs[varIdx].name;
163 
164  //play with the binning choice for given variables:
165  int rebinFactor=1;
166 
167  if(varName=="PtP") rebinFactor = 5;
168 
169  //create a canvas
170  TCanvas *c2 =new TCanvas("c2","c2",500,400);
171 
172  TPad *padAbs = new TPad("padAbs", "padAbs", 0, 0.3, 1, 1.0);
173  padAbs -> SetBottomMargin(0);
174  padAbs->SetFillStyle(0);
175  padAbs->SetFillColor(0);
176  padAbs -> SetLeftMargin(0.15);
177 
178  padAbs -> Draw();
179  padAbs -> cd();
180 
181  double legendxlow=0.45;
182  double legendxhigh=0.93;
183 
184  if(varName=="trueW2"||varName=="CosNumi"){legendxlow=0.17; legendxhigh=0.6;}
185  //if(varName=="CosNumi"){legendxlow=0.16; legendxhigh=0.56;}
186 
187  //Pull down the variable we want
188  IPrediction* predND = predsND[selIdx][varIdx];
189  IPrediction* predNDnon = predsNDnon[selIdx][varIdx];
190 
191  const double kPOT = predND->Predict(&noosc).POT();
192 
193  auto predFD = ana::LoadFrom<IPrediction> (dpred_no[selIdx][varIdx]);
194 
195  //Get components of the prediction and store them in histograms
196 
197  TH1* hFD = predFD->PredictComponent(&mycalc,
198  Flavors::kNuMuToNuE,
199  Current::kCC,
200  Sign::kBoth).ToTH1(kPOT);
201 
202  TH1* hAux = predND->Predict(&noosc).ToTH1(kPOT);
203  TH1* hAuxnon = predNDnon->Predict(&noosc).ToTH1(kPOT);
204 
205  TH1* hND = predND->Predict(&noosc).ToTH1(kPOT*hFD->Integral()/hAux->Integral());
206  TH1* hNDnon = predNDnon->Predict(&noosc).ToTH1(kPOT*hFD->Integral()/hAuxnon->Integral());
207 
208  hND->SetLineColor(kTotalMCColor);
209  hNDnon->SetLineColor(1);
210 
211  hFD->SetLineColor(kNueSignalColor);
212  hFD->Rebin(rebinFactor);
213  hND->Rebin(rebinFactor);
214  hNDnon->Rebin(rebinFactor);
215 
216  std::cout << selName << " " << varName << ": " << hND->Integral() << std::endl;
217 
218  hND->GetXaxis()->CenterTitle();
219  hND->GetYaxis()->CenterTitle();
220  hND->GetXaxis()->SetDecimals();
221  hND->GetYaxis()->SetDecimals();
222  hND->GetYaxis()->SetTitleOffset(.8);
223  hND->GetXaxis()->SetLabelSize(.0);
224  hND->GetYaxis()->SetLabelSize(.06);
225  hND->GetYaxis()->SetTitleSize(.06);
226  hND->SetMaximum(1.1*TMath::Max(hND->GetMaximum(), TMath::Max(hNDnon->GetMaximum(), hFD->GetMaximum())));
227  if(varName=="trueW2"){
228  hND->SetMaximum(1.2*hND->GetMaximum());
229  }
230 
231  hND->GetYaxis()->SetTitle(TString::Format("Events area normalized"));
232 
233  if(varName=="recoQ2" || varName=="trueQ2"){
234  hND->GetXaxis()->SetRangeUser(0,2);
235  }
236 
237  if(varName=="recoE"){
238  hND->GetXaxis()->SetTitle("Reconstructed Neutrino Energy (GeV)");
239  }
240 
241  if(varName=="CosNumi") hND->GetXaxis()->SetRangeUser(0.5,1.0);
242 
243  hND->GetXaxis()->SetNdivisions(405,kFALSE);
244  hND->SetMinimum(.001);
245 
246  hND ->DrawCopy("h");
247 
248  hNDnon->Draw("h same");
249 
250  hND->Draw("h same");
251  hFD->Draw("hist same");
252 
253  //Fix up the axis
254  padAbs->RedrawAxis();
255  padAbs->Update();
256 
257  //Add the legend and preliminary
258  Legend(legendxlow, .45, legendxhigh, .85);
259  std::string wTitle = "FD FHC CORE sample: Q^{2} reweighted";
260 
261  if(weight_name == "Cos") wTitle = "FD FHC CORE sample: Cos reweighted";
262  if(weight_name == "PtP") wTitle = "FD FHC CORE sample: p_{T}/p reweighted";
263  wTitle = "NOvA Simulation";
264 
265  TLatex* selTitle = new TLatex(.935, .95, wTitle.c_str());
266  selTitle->SetNDC();
267  selTitle->SetTextSize(2/30.);
268  selTitle->SetTextAlign(32);
269  selTitle->SetTextColor(kGray+1);
270  selTitle->Draw();
271 
272  c2->cd();
273  TPad *padRatio = new TPad("padRatio", "padRatio", 0, 0.0, 1, 0.33);
274  padRatio -> SetTopMargin(0.11);
275  padRatio -> SetBottomMargin(0.33);
276  padRatio -> SetFillStyle(0);
277  padRatio -> SetLeftMargin(0.15);
278 
279  padRatio -> Draw();
280  padRatio -> cd();
281 
282  TH1* hRatio = (TH1*)hND -> Clone("hRatio");
283  TH1* hRationon = (TH1*)hNDnon -> Clone("hRationon");
284 
285  double b=5;
286  double a=0;
287  if(varName=="trueW2") {b = 2.0;}
288  if(varName=="trueQ2") {b = 2.0;}
289  if(varName=="recoQ2") {b = 2.0;}
290  if(varName=="PtP" || varName=="CosNumi") {b = 1.0;}
291  if(varName=="CosNumi") a=0.5;
292 
293  TLine* lOne = new TLine(a,1.0,b,1.0);
294  lOne -> SetLineStyle(2);
295 
296  hRatio -> Divide(hFD);
297  hRationon -> Divide(hFD);
298 
299  hRatio->GetXaxis()->CenterTitle();
300  hRatio->GetYaxis()->CenterTitle();
301  hRatio->GetYaxis()->SetLabelSize(.12);
302  hRatio->GetXaxis()->SetLabelSize(.12);
303  hRatio->GetXaxis()->SetTitleOffset(1.);
304  hRatio->GetXaxis()->SetLabelOffset(0.03);
305  hRatio->GetXaxis()->SetTickSize(.1);
306  hRatio->GetYaxis()->SetTitleSize(.12);
307  hRatio->GetXaxis()->SetTitleSize(.12);
308  hRatio->GetYaxis()->SetRangeUser(0,2);
309  if(varName=="CosNumi"){
310  hRatio->GetYaxis()->SetRangeUser(0.5,2.5);
311  }
312  hRatio->SetMarkerSize(.1);
313  hRatio->GetXaxis()->SetDecimals();
314  hRatio->GetYaxis()->SetDecimals();
315  hRatio->GetYaxis()->SetTitleOffset(0.4);
316  hRatio->GetYaxis()->SetTitle("Ratio to FD #nu_{e}");
317 
318  if(varName=="recoE"){
319  hRatio->GetXaxis()->SetTitle("Reconstructed Neutrino Energy (GeV)");
320  }
321  if(varName=="recoQ2"){
322  hRatio->GetXaxis()->SetTitle("ND Reconstructed Q^{2} or FD True Q^{2} (GeV^{2})");
323  }
324  if(varName=="CosNumi"){
325  hRatio->GetXaxis()->SetTitle("cos #theta");
326  }
327  if(varName=="PtP"){
328  hRatio->GetXaxis()->SetTitle("p_{T}/p");
329  }
330 
331  if(varName=="recoQ2" || varName=="trueQ2"){
332  hRatio->GetXaxis()->SetRangeUser(0,2);
333  }
334  if(varName=="trueW2"){
335  hRatio->GetXaxis()->SetTitle("True W^{2} (GeV^{2})");
336  }
337 
338  hRatio->GetYaxis()->SetNdivisions(502,kFALSE);
339  hRatio->GetXaxis()->SetNdivisions(405,kFALSE);
340 
341  c2->Update();
342  hRatio-> DrawCopy("hist");
343  lOne -> Draw("same");
344  hRationon->Draw("hist same");
345  hRatio->Draw("hist same");
346  padRatio->RedrawAxis();
347  padRatio->Update();
348 
349  c2->Update();
350  c2->Print(("plots/FDvsND_REW_"+weight_name+"_"+selName+"_"+varName+".png").c_str());
351  c2->Print(("plots/FDvsND_REW_"+weight_name+"_"+selName+"_"+varName+".pdf").c_str());
352 
353  } // end for varIdx
354  } // end for selIdx
355 }
356 
Draw
tree Draw("slc.nhit")
osc::_NoOscillations
Pass neutrinos through unchanged.
Definition: IOscCalc.h:45
Divide
ratio_hxv Divide(hxv, goal_hxv)
ana::kNumVars
const int kNumVars
Definition: vars.h:14
ana
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
y1
Float_t y1[n_points_granero]
Definition: compare.C:5
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
ana::IPrediction::Predict
virtual Spectrum Predict(osc::IOscCalc *calc) const =0
x1
Float_t x1[n_points_granero]
Definition: compare.C:5
SetLineStyle
hmean SetLineStyle(2)
plot_validation_datamc.Clone
def Clone(obj)
Definition: plot_validation_datamc.py:15
ana::HistDef::name
std::string name
Definition: NuePlotLists.h:12
ana::kTotalMCErrorBandColor
const Color_t kTotalMCErrorBandColor
Definition: Style.h:17
SetFillStyle
ntuple SetFillStyle(1001)
ana::kNueSignalColor
const Color_t kNueSignalColor
Definition: Style.h:19
Legend
void Legend(double x0, double y0, double x1, double y1)
Definition: plot_3NDvsFD_FHC.C:47
demo5.c2
c2
Definition: demo5.py:33
ana::Current::kCC
Charged-current interactions.
Definition: IPrediction.h:39
a
const double a
Definition: gUpMuFluxGen.cxx:164
osc::OscCalcDumb
Simple oscillation probability calculator that has no solar term or mass hierarchy or delta so it&#39;s s...
Definition: OscCalcDumb.h:16
ana::defs
const HistDef defs[kNumVars]
Definition: vars.h:15
ana::LoadFrom< IPrediction >
std::unique_ptr< IPrediction > LoadFrom< IPrediction >(TDirectory *dir, const std::string &label)
Definition: IPrediction.cxx:18
ana::kNumSels
const int kNumSels
Definition: vars.h:43
noosc
osc::NoOscillations noosc
Definition: NueSystFuncs.h:44
ana::Spectrum::POT
double POT() const
Definition: Spectrum.h:227
om::cout
OStream cout
Definition: OStream.cxx:6
plot_3NDvsFD_FHC
void plot_3NDvsFD_FHC(const std::string weight_name, const std::string fname_NDnon, const std::string fname_ND, const std::string fname_FD)
Definition: plot_3NDvsFD_FHC.C:75
b
const hit & b
Definition: hits.cxx:21
ana::Sign::kBoth
Both neutrinos and antineutrinos.
Definition: IPrediction.h:52
ana::IPrediction
Standard interface to all prediction techniques.
Definition: IPrediction.h:57
genie::utils::style::Format
void Format(TGraph *gr, int lcol, int lsty, int lwid, int mcol, int msty, double msiz)
Definition: Style.cxx:154
ana::kTotalMCColor
const Color_t kTotalMCColor
Definition: Style.h:16
cd
c cd(1)
string
enum BeamMode string
Generated on Wed Feb 24 2021 14:11:37 for NOvA by   doxygen 1.8.11