Functions
plot_3NDvsFD.C File Reference
#include "CAFAna/Analysis/Exposures.h"
#include "CAFAna/Analysis/Calcs.h"
#include "CAFAna/Core/Binning.h"
#include "CAFAna/Core/Cut.h"
#include "CAFAna/Core/HistAxis.h"
#include "CAFAna/Core/LoadFromFile.h"
#include "CAFAna/Core/Spectrum.h"
#include "CAFAna/Core/SpectrumLoader.h"
#include "CAFAna/Core/SpectrumLoaderBase.h"
#include "CAFAna/Cuts/Cuts.h"
#include "3FlavorAna/Cuts/NumuCuts.h"
#include "CAFAna/Cuts/NueCutsSecondAna.h"
#include "CAFAna/Cuts/SpillCuts.h"
#include "CAFAna/Decomp/ProportionalDecomp.h"
#include "3FlavorAna/Decomp/NumuDecomp.h"
#include "CAFAna/Extrap/ModularExtrap.h"
#include "CAFAna/Prediction/PredictionExtrap.h"
#include "CAFAna/Prediction/PredictionNoExtrap.h"
#include "CAFAna/Vars/GenieWeights.h"
#include "3FlavorAna/Vars/NueVars.h"
#include "CAFAna/Vars/Vars.h"
#include "OscLib/OscCalcPMNSOpt.h"
#include "OscLib/OscCalcDumb.h"
#include "CAFAna/Analysis/Style.h"
#include "TFile.h"
#include "TH1.h"
#include "TH2.h"
#include "TCanvas.h"
#include "TGaxis.h"
#include "TLatex.h"
#include "TLegend.h"
#include "TLine.h"
#include "TSystem.h"
#include <iostream>
#include <iomanip>
#include <fstream>

Go to the source code of this file.

Functions

void Legend (double x0, double y0, double x1, double y1)
 
void plot_3NDvsFD (const std::string weight_name, const std::string fname_NDnon, const std::string fname_ND, const std::string fname_FD)
 

Function Documentation

void Legend ( double  x0,
double  y0,
double  x1,
double  y1 
)

Definition at line 45 of file plot_3NDvsFD.C.

References APDHVSetting::dummy, ana::kNueSignalColor, ana::kTotalMCColor, ana::kTotalMCErrorBandColor, and MECModelEnuComparisons::leg.

Referenced by plot_3NDvsFD().

46 {
47  // Doesn't handle log-y well
48  // TLegend* leg = AutoPlaceLegend(x1-x0, y1-y0);
49  TLegend* leg = new TLegend(x0, y0, x1, y1);
50  leg->SetFillStyle(0);
51 
52  TH1* dummy = new TH1F("", "", 1, 0, 1);
53  TH1* dummyFill = new TH1F("", "", 1, 0, 1);
54  dummyFill->SetLineColor(kTotalMCColor);
55  dummyFill->SetFillColor(kTotalMCErrorBandColor);
56 
57  dummy->SetMarkerStyle(kFullCircle);
58 
59  dummy->SetLineColor(kGray+2);
60  dummy->SetMarkerColor(kGray+2);
61  dummy->SetLineColor(kTotalMCColor);
62  leg->AddEntry(dummy->Clone(), "ND #nu_{#mu} Total MC REW.", "l");
63  dummy->SetLineColor(1);
64  leg->AddEntry(dummy->Clone(), "ND #nu_{#mu} Total MC NOM.", "l");
65 
66  dummy->SetLineColor(kNueSignalColor);
67  leg->AddEntry(dummy->Clone(), "FD signal #nu_{e}", "l");
68 
69  leg->Draw();
70 }
y1
Float_t y1[n_points_granero]
Definition: compare.C:5
x1
Float_t x1[n_points_granero]
Definition: compare.C:5
ana::kTotalMCErrorBandColor
const Color_t kTotalMCErrorBandColor
Definition: Style.h:17
ana::kNueSignalColor
const Color_t kNueSignalColor
Definition: Style.h:19
ana::kTotalMCColor
const Color_t kTotalMCColor
Definition: Style.h:16
void plot_3NDvsFD ( const std::string  weight_name,
const std::string  fname_NDnon,
const std::string  fname_ND,
const std::string  fname_FD 
)

Definition at line 72 of file plot_3NDvsFD.C.

References std::asin(), b, demo5::c2, cd(), plot_validation_datamc::Clone(), om::cout, ana::defs, Divide(), Draw(), allTimeWatchdog::endl, stan::math::fabs(), genie::utils::style::Format(), osc::_IOscCalcAdjustable< T >::GetDmsq32(), ana::Sign::kBoth, ana::Current::kCC, ana::kNueSignalColor, ana::kNumSels, ana::Flavors::kNuMuToNuE, ana::kNumVars, ana::kPOT, ana::kTotalMCColor, Legend(), ana::LoadFrom< IPrediction >(), ana::HistDef::name, noosc, ana::Spectrum::POT(), ana::IPrediction::Predict(), runNovaSAM::release, ana::ResetOscCalcToDefault(), osc::_OscCalcPMNSOpt< T >::SetdCP(), osc::_OscCalcPMNSOpt< T >::SetDmsq32(), SetFillStyle(), SetLineStyle(), osc::_OscCalcPMNSOpt< T >::SetTh23(), std::sqrt(), string, and ana::Spectrum::ToTH1().

73 {
74  const double kPOT = 8.09e20;
75 
76  struct HistDef
77  {
78  std::string name;
79  };
80 
81  //which variables do you actually want to look at?
82  const int kNumVars = 5;
83  const HistDef defs[kNumVars] = {
84  "trueQ2",
85  "recoQ2",
86  "trueW2",
87  "PtP",
88  "CosNumi"
89  };
90 
91  //which cut tiers do you want to look at?
92  const int kNumSels = 1;
93  const std::string selNamesND[kNumSels] = {
94  "numuND"
95  };
96 
97  const std::string selNamesFD[kNumSels] = {
98  "AllSamples"
99  };
100 
101  IPrediction* predsND[kNumSels][kNumVars];
102  IPrediction* predsNDnon[kNumSels][kNumVars];
103 
104  osc::NoOscillations noosc;
105  osc::OscCalcPMNSOpt mycalc;
106  ResetOscCalcToDefault(&mycalc);
107 
108  mycalc.SetDmsq32(+fabs(mycalc.GetDmsq32()));
109  mycalc.SetdCP(3*TMath::Pi()/2);
110  mycalc.SetTh23(asin(sqrt(0.404)));
111 
112  //pick up spectra/prediction objects from the input root file
113  for(int selIdx = 0; selIdx < kNumSels; ++selIdx){
114  const char* selName = selNamesND[selIdx].c_str();
115  const char* wName = weight_name.c_str();
116 
117  for(int varIdx = 0; varIdx < kNumVars; ++varIdx){
118  const char* varName = defs[varIdx].name.c_str();
119 
120  std::cout<<TString::Format("%s/spect_%s_%s", selName, wName, varName).Data()<<std::endl;
121 
122  predsND[selIdx][varIdx] = LoadFromFile<IPrediction>(fname_ND, TString::Format("%s/pred_%s_%s", selName, wName, varName).Data()).release();
123  predsNDnon[selIdx][varIdx] = LoadFromFile<IPrediction>(fname_NDnon, TString::Format("%s/pred_%s", selName, varName).Data()).release();
124 
125  if(selIdx == 0 && varIdx == 0)
126  std::cout << "ND MC POT: "
127  << predsNDnon[selIdx][varIdx]->Predict(&noosc).POT() << std::endl;
128 
129  } // end for varIdx
130  } // end for selIdx
131 
132 
133  TFile * file_no = new TFile (fname_FD.c_str(),"READ");
134  TDirectory * dpred_no[kNumSels][kNumVars];
135  for(int selIdx = 0; selIdx < kNumSels; ++selIdx){
136  const std::string selName = selNamesFD[selIdx];
137  for(int varIdx = 0; varIdx < kNumVars; ++varIdx){
138  const std::string varName = defs[varIdx].name;
139  if(varName == "recoQ2") {
140  dpred_no[selIdx][varIdx] = file_no->GetDirectory(("prediction/nue_pred_noextrap_"+selName+"_trueQ2").c_str());
141  } else if(varName == "CosNumi"){
142  dpred_no[selIdx][varIdx] = file_no->GetDirectory(("prediction/nue_pred_noextrap_"+selName+"_CosTheta").c_str());
143  } else {
144  dpred_no[selIdx][varIdx] = file_no->GetDirectory(("prediction/nue_pred_noextrap_"+selName+"_"+varName).c_str());
145  }
146  }
147  }
148 
149  TH1* hFullNDHists[kNumVars-1];
150  for(int selIdx = 0; selIdx < kNumSels; ++selIdx){
151  const std::string selName = selNamesFD[selIdx];
152  for(int varIdx = 0; varIdx < kNumVars; ++varIdx){
153 
154  const std::string varName = defs[varIdx].name;
155 
156  //play with the binning choice for given variables:
157  int rebinFactor=1;
158 
159  //create a canvas
160  TCanvas *c2 =new TCanvas("c2","c2",500,400);
161 
162 
163  TPad *padAbs = new TPad("padAbs", "padAbs", 0, 0.3, 1, 1.0);
164  padAbs -> SetBottomMargin(0);
165  padAbs->SetFillStyle(0);
166  padAbs -> SetLeftMargin(0.15);
167 
168  padAbs -> Draw();
169  padAbs -> cd();
170 
171  double legendxlow=0.45;
172  double legendxhigh=0.85;
173 
174  if(varName=="trueW2"||varName=="CosNumi"){legendxlow=0.2; legendxhigh=0.6;}
175 
176  //Pull down the variable we want
177 
178  IPrediction* predND = predsND[selIdx][varIdx];
179  IPrediction* predNDnon = predsNDnon[selIdx][varIdx];
180 
181  auto predFD = ana::LoadFrom<IPrediction> (dpred_no[selIdx][varIdx]);
182 
183  //Get components of the prediction and store them in histograms
184 
185  TH1* hFD = predFD->PredictComponent(&mycalc,
186  Flavors::kNuMuToNuE,
187  Current::kCC,
188  Sign::kBoth).ToTH1(kPOT);
189 
190  TH1* hAux = predND->Predict(&noosc).ToTH1(kPOT);
191  TH1* hAuxnon = predNDnon->Predict(&noosc).ToTH1(kPOT);
192 
193  TH1* hND = predND->Predict(&noosc).ToTH1(kPOT*hFD->Integral()/hAux->Integral());
194  TH1* hNDnon = predNDnon->Predict(&noosc).ToTH1(kPOT*hFD->Integral()/hAuxnon->Integral());
195 
196  hND->SetLineColor(kTotalMCColor);
197  hNDnon->SetLineColor(1);
198 
199  hFD->SetLineColor(kNueSignalColor);
200  hND->Rebin(rebinFactor);
201  hNDnon->Rebin(rebinFactor);
202 
203  std::cout << selName << " " << varName << ": " << hND->Integral() << std::endl;
204 
205  hND->GetXaxis()->CenterTitle();
206  hND->GetYaxis()->CenterTitle();
207  hND->GetXaxis()->SetDecimals();
208  hND->GetYaxis()->SetDecimals();
209  hND->GetYaxis()->SetTitleOffset(1.15);
210  hND->GetXaxis()->SetLabelSize(.0);
211  hND->GetYaxis()->SetLabelSize(.05);
212  hND->SetMaximum(1.1*TMath::Max(hND->GetMaximum(), TMath::Max(hNDnon->GetMaximum(), hFD->GetMaximum())));
213  if(varName=="trueW2"){
214  hND->SetMaximum(1.2*hND->GetMaximum());
215  }
216 
217  hND->GetYaxis()->SetTitle(TString::Format("Events"));
218 
219  if(varName=="recoQ2" || varName=="trueQ2"){
220  hND->GetXaxis()->SetRangeUser(0,2);
221  }
222 
223  if(varName=="recoE"){
224  hND->GetXaxis()->SetTitle("Reconstructed neutrino energy (GeV)");
225  }
226  hND->GetXaxis()->SetNdivisions(405,kFALSE);
227  hND->SetMinimum(.001);
228 
229  hND ->DrawCopy("hist");
230 
231  hNDnon->Draw("hist same");
232 
233  hND->Draw("hist same");
234  hFD->Draw("hist same");
235 
236  //Fix up the axis
237  padAbs->RedrawAxis();
238  padAbs->Update();
239 
240  //Add the legend and preliminary
241  Legend(legendxlow, .63, legendxhigh, .85);
242  std::string wTitle = weight_name+" reweighted";
243 
244  TLatex* selTitle = new TLatex(.935, .95, wTitle.c_str());
245  selTitle->SetNDC();
246  selTitle->SetTextSize(2/30.);
247  selTitle->SetTextAlign(32);
248  selTitle->Draw();
249 
250  c2->cd();
251  TPad *padRatio = new TPad("padRatio", "padRatio", 0, 0.05, 1, 0.3);
252  padRatio -> SetTopMargin(0);
253  padRatio -> SetBottomMargin(0.25);
254  padRatio -> SetFillStyle(0);
255  padRatio -> SetLeftMargin(0.15);
256 
257  padRatio -> Draw();
258  padRatio -> cd();
259 
260  TH1* hRatio = (TH1*)hND -> Clone("hRatio");
261  TH1* hRationon = (TH1*)hNDnon -> Clone("hRationon");
262 
263  double b=5;
264  if(varName=="trueW2") {b = 2.0;}
265  if(varName=="trueQ2") {b = 2.0;}
266  if(varName=="recoQ2") {b = 2.0;}
267  if(varName=="PtP" || varName=="CosNumi") {b = 1.0;}
268 
269  TLine* lOne = new TLine(0,1.0,b,1.0);
270  lOne -> SetLineStyle(2);
271 
272  hRatio -> Divide(hFD);
273  hRationon -> Divide(hFD);
274 
275  hRatio->GetXaxis()->CenterTitle();
276  hRatio->GetYaxis()->CenterTitle();
277  hRatio->GetYaxis()->SetLabelSize(4.0/30.);
278  hRatio->GetXaxis()->SetLabelSize(4.0/30.);
279  hRatio->GetXaxis()->SetLabelOffset(.04);
280  hRatio->GetXaxis()->SetTickSize(.06);
281  hRatio->GetYaxis()->SetTitleSize(4.5/30.);
282  hRatio->GetXaxis()->SetTitleSize(4.5/30.);
283  hRatio->GetYaxis()->SetRangeUser(0,2);
284  if(varName=="CosNumi"){
285  hRatio->GetYaxis()->SetRangeUser(0.5,2.5);
286  }
287  hRatio->SetMarkerSize(.1);
288  hRatio->GetXaxis()->SetDecimals();
289  hRatio->GetYaxis()->SetDecimals();
290  hRatio->GetYaxis()->SetTitleOffset(0.4);
291  hRatio->GetYaxis()->SetTitle("#nu_{#mu} ND/Sig. FD");
292 
293  if(varName=="recoE"){
294  hRatio->GetXaxis()->SetTitle("Reconstructed neutrino energy (GeV)");
295  }
296  if(varName=="recoQ2"){
297  hRatio->GetXaxis()->SetTitle("ND reco Q^{2} or FD true Q^{2} (GeV^{2})");
298  }
299 
300  if(varName=="recoQ2" || varName=="trueQ2"){
301  hRatio->GetXaxis()->SetRangeUser(0,2);
302  }
303 
304  hRatio->GetYaxis()->SetNdivisions(502,kFALSE);
305  hRatio->GetXaxis()->SetNdivisions(405,kFALSE);
306 
307  hRatio-> DrawCopy("hist");
308  lOne -> Draw("same");
309  hRationon->Draw("hist same");
310  hRatio->Draw("hist same");
311 
312  c2->Print(("plots/FDvsND_REW_"+weight_name+"_"+selName+"_"+varName+".png").c_str());
313 
314  } // end for varIdx
315  } // end for selIdx
316 }
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::ResetOscCalcToDefault
void ResetOscCalcToDefault(osc::IOscCalcAdjustable *calc)
Reset calculator to default assumptions for all parameters.
Definition: Calcs.cxx:23
stan::math::fabs
fvar< T > fabs(const fvar< T > &x)
Definition: fabs.hpp:15
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
SetLineStyle
hmean SetLineStyle(2)
std::sqrt
T sqrt(T number)
Definition: d0nt_math.hpp:156
plot_validation_datamc.Clone
def Clone(obj)
Definition: plot_validation_datamc.py:15
ana::HistDef::name
std::string name
Definition: NuePlotLists.h:12
SetFillStyle
ntuple SetFillStyle(1001)
ana::kNueSignalColor
const Color_t kNueSignalColor
Definition: Style.h:19
demo5.c2
c2
Definition: demo5.py:33
ana::Current::kCC
Charged-current interactions.
Definition: IPrediction.h:39
osc::_OscCalcPMNSOpt
Optimized version of OscCalcPMNS.
Definition: StanTypedefs.h:31
ana::defs
const HistDef defs[kNumVars]
Definition: vars.h:15
osc::_OscCalcPMNSOpt::SetTh23
void SetTh23(const T &th23) override
Definition: OscCalcPMNSOpt.h:42
osc::_IOscCalcAdjustable::GetDmsq32
virtual T GetDmsq32() const
Definition: IOscCalc.h:91
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
osc::_OscCalcPMNSOpt::SetDmsq32
void SetDmsq32(const T &dmsq32) override
Definition: OscCalcPMNSOpt.h:39
osc::_OscCalcPMNSOpt::SetdCP
void SetdCP(const T &dCP) override
Definition: OscCalcPMNSOpt.h:43
b
const hit & b
Definition: hits.cxx:21
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
Legend
void Legend(double x0, double y0, double x1, double y1)
Definition: plot_3NDvsFD.C:45
cd
c cd(1)
std::asin
T asin(T number)
Definition: d0nt_math.hpp:60
string
enum BeamMode string
Generated on Fri Apr 9 2021 14:13:53 for NOvA by   doxygen 1.8.11