plot_NDvsFD_weights.C
Go to the documentation of this file.
3 #include "CAFAna/Core/Binning.h"
4 #include "CAFAna/Core/Cut.h"
5 #include "CAFAna/Core/HistAxis.h"
7 #include "CAFAna/Core/Spectrum.h"
10 #include "CAFAna/Cuts/Cuts.h"
12 #include "CAFAna/Cuts/NueCutsSecondAna.h"
13 #include "CAFAna/Cuts/SpillCuts.h"
21 #include "CAFAna/Vars/Vars.h"
22 
23 #include "OscLib/OscCalcPMNSOpt.h"
24 #include "OscLib/OscCalcDumb.h"
25 #include "CAFAna/Analysis/Style.h"
26 
27 //#include "CAFAna/nue/SecondAna/draw_plots_util.h"
28 
29 using namespace ana;
30 
31 #include "TFile.h"
32 #include "TH1.h"
33 #include "TH2.h"
34 #include "TCanvas.h"
35 #include "TGaxis.h"
36 #include "TLatex.h"
37 #include "TLegend.h"
38 #include "TLine.h"
39 #include "TSystem.h"
40 #include <iostream>
41 #include <iomanip>
42 #include <fstream>
43 
44 //make a pretty legend
45 void Legend(double x0, double y0, double x1, double y1)
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  //leg->AddEntry(dummy->Clone(), "ND data", "lep");
60  dummy->SetLineColor(kGray+2);
61  dummy->SetMarkerColor(kGray+2);
62  // leg->AddEntry(dummy->Clone(), "ND data stagger", "lep");
63  dummy->SetLineColor(1);
64  leg->AddEntry(dummy->Clone(), "ND #nu_{#mu} Total MC", "l");
65  dummy->SetLineColor(kNueSignalColor);
66  leg->AddEntry(dummy->Clone(), "FD signal #nu_{e}", "l");
67 
68  leg->Draw();
69 }
70 
71 void plot_NDvsFD_weights(const std::string fname_ND = "datamc_ND_numu_kinematics.root", const std::string fname_FD = "FDprediction_kinematics.root")
72 {
73  const double kPOT = 8.09e20;
74 
75  struct HistDef
76  {
78  };
79 
80  //which variables do you actually want to look at?
81  const int kNumVars = 5;
82  const HistDef defs[kNumVars] = {
83  "trueQ2",
84  "recoQ2",
85  "trueW2",
86  "PtP",
87  "CosNumi"
88  };
89 
90  //which cut tiers do you want to look at?
91  const int kNumSels = 1;
92  const std::string selNamesND[kNumSels] = {
93  "numuND"
94  };
95  const std::string selNamesFD[kNumSels] = {
96  "AllSamples"
97  };
98 
99  IPrediction* predsND[kNumSels][kNumVars];
100 
102  osc::OscCalcPMNSOpt mycalc;
103  ResetOscCalcToDefault(&mycalc);
104 
105  mycalc.SetDmsq32(+fabs(mycalc.GetDmsq32()));
106  mycalc.SetdCP(3*TMath::Pi()/2);
107  mycalc.SetTh23(asin(sqrt(0.404)));
108 
109  //pick up spectra/prediction objects from the input ND root file
110  for(int selIdx = 0; selIdx < kNumSels; ++selIdx){
111  const char* selName = selNamesND[selIdx].c_str();
112  for(int varIdx = 0; varIdx < kNumVars; ++varIdx){
113  const char* varName = defs[varIdx].name.c_str();
114 
115  std::cout<<TString::Format("%s/spect_%s", selName, varName).Data()<<std::endl;
116 
117  predsND[selIdx][varIdx] = LoadFromFile<IPrediction>(fname_ND, TString::Format("%s/pred_%s", selName, varName).Data()).release();
118 
119  if(selIdx == 0 && varIdx == 0)
120  std::cout << "ND MC POT: "
121  << predsND[selIdx][varIdx]->Predict(&noosc).POT() << std::endl;
122 
123  } // end for varIdx
124  } // end for selIdx
125 
126 
127  TFile * file_no = new TFile (fname_FD.c_str(),"READ");
128  TDirectory * dpred_no[kNumSels][kNumVars];
129  for(int selIdx = 0; selIdx < kNumSels; ++selIdx){
130  const std::string selName = selNamesFD[selIdx];
131  for(int varIdx = 0; varIdx < kNumVars; ++varIdx){
132  const std::string varName = defs[varIdx].name;
133 
134  // need to pick equivalent variables
135  if(varName == "recoQ2") {
136  dpred_no[selIdx][varIdx] = file_no->GetDirectory(("prediction/nue_pred_noextrap_"+selName+"_trueQ2").c_str());
137  } else if(varName == "CosNumi"){
138  dpred_no[selIdx][varIdx] = file_no->GetDirectory(("prediction/nue_pred_noextrap_"+selName+"_CosTheta").c_str());
139  } else {
140  dpred_no[selIdx][varIdx] = file_no->GetDirectory(("prediction/nue_pred_noextrap_"+selName+"_"+varName).c_str());
141  }
142  }
143  }
144 
145  TH1* hFullNDHists[kNumVars-1];
146  TFile * file = new TFile("Weights.root","RECREATE");
147 
148  //loop over the different selector/cut levels to plot each one
149  for(int selIdx = 0; selIdx < kNumSels; ++selIdx){
150  const std::string selName = selNamesFD[selIdx];
151  for(int varIdx = 0; varIdx < kNumVars; ++varIdx){
152 
153  const std::string varName = defs[varIdx].name;
154 
155  //play with the binning choice for given variables:
156  int rebinFactor=1;
157 
158  //create a canvas
159  TCanvas *c2 =new TCanvas("c2","c2",500,400);
160 
161  //absolute pad
162  TPad *padAbs = new TPad("padAbs", "padAbs", 0, 0.3, 1, 1.0);
163  padAbs -> SetBottomMargin(0);
164  padAbs->SetFillStyle(0);
165  padAbs -> SetLeftMargin(0.15);
166 
167  padAbs -> Draw();
168  padAbs -> cd();
169 
170  double legendxlow=0.55;
171  double legendxhigh=0.95;
172 
173  if(varName=="trueW2"||varName=="CosNumi"){legendxlow=0.2; legendxhigh=0.5;}
174 
175  //Pull down the variable we want
176 
177  IPrediction* predND = predsND[selIdx][varIdx];
178  auto predFD = ana::LoadFrom<IPrediction> (dpred_no[selIdx][varIdx]);
179 
180  //Get components of the prediction and store them in histograms
181 
182  TH1* hFD = predFD->PredictComponent(&mycalc,
184  Current::kCC,
185  Sign::kBoth).ToTH1(kPOT);
186 
187  TH1* hAux = predND->Predict(&noosc).ToTH1(kPOT);
188 
189  TH1* hND = predND->Predict(&noosc).ToTH1(kPOT*hFD->Integral()/hAux->Integral());
190 
191  hND->SetLineColor(1);
192  hFD->SetLineColor(kNueSignalColor);
193 
194  hND->Rebin(rebinFactor);
195 
196  std::cout << selName << " " << varName << ": " << hND->Integral() << std::endl;
197 
198  hND->GetXaxis()->CenterTitle();
199  hND->GetYaxis()->CenterTitle();
200  hND->GetXaxis()->SetDecimals();
201  hND->GetYaxis()->SetDecimals();
202  hND->GetYaxis()->SetTitleOffset(1.15);
203  hND->GetXaxis()->SetLabelSize(.0);
204  hND->GetYaxis()->SetLabelSize(.05);
205  hND->SetMaximum(1.1*hND->GetMaximum());
206  hND->GetYaxis()->SetTitle(TString::Format("Events"));
207 
208  if(varName=="recoQ2" || varName=="trueQ2"){
209  hND->GetXaxis()->SetRangeUser(0,2);
210  }
211 
212  if(varName=="recoE"){
213  hND->GetXaxis()->SetTitle("Reconstructed neutrino energy (GeV)");
214  }
215  hND->GetXaxis()->SetNdivisions(405,kFALSE);
216  hND->SetMinimum(.001);
217 
218  hND ->DrawCopy("hist");
219  hND->Draw("hist same");
220  hFD->Draw("hist same");
221 
222  //Fix up the axis
223  padAbs->RedrawAxis();
224  padAbs->Update();
225 
226  //Add the legend and preliminary
227  Legend(legendxlow, .70, legendxhigh, .85);
228  //Preliminary();
229 
230  TLatex* selTitle = new TLatex(.935, .95, "nominal");
231  selTitle->SetNDC();
232  selTitle->SetTextSize(2/30.);
233  selTitle->SetTextAlign(32);
234  selTitle->Draw();
235 
236  c2->cd();
237 
238  //ratio pad
239  TPad *padRatio = new TPad("padRatio", "padRatio", 0, 0.05, 1, 0.3);
240  padRatio -> SetTopMargin(0);
241  padRatio -> SetBottomMargin(0.25);
242  padRatio -> SetFillStyle(0);
243  padRatio -> SetLeftMargin(0.15);
244 
245  padRatio -> Draw();
246  padRatio -> cd();
247 
248  TH1* hRatio = (TH1*)hND -> Clone("hRatio");
249 
250  double b=5;
251  if(varName=="trueW2") {b = 2.0;}
252  if(varName=="trueQ2") {b = 2.0;}
253  if(varName=="recoQ2") {b = 2.0;}
254  if(varName=="PtP" || varName=="CosNumi") {b = 1.0;}
255 
256  TLine* lOne = new TLine(0,1.0,b,1.0);
257  lOne -> SetLineStyle(2);
258 
259  hRatio -> Divide(hFD);
260 
261  TH1* hWeight = (TH1*)hFD -> Clone("hWeight");
262  hWeight -> Divide(hND);
263 
264  hRatio->GetXaxis()->CenterTitle();
265  hRatio->GetYaxis()->CenterTitle();
266  hRatio->GetYaxis()->SetLabelSize(4.0/30.);
267  hRatio->GetXaxis()->SetLabelSize(4.0/30.);
268  hRatio->GetXaxis()->SetLabelOffset(.04);
269  hRatio->GetXaxis()->SetTickSize(.06);
270  hRatio->GetYaxis()->SetTitleSize(4.5/30.);
271  hRatio->GetXaxis()->SetTitleSize(4.5/30.);
272  hRatio->GetYaxis()->SetRangeUser(0,2.0);
273  hRatio->SetMarkerSize(.1);
274  hRatio->GetXaxis()->SetDecimals();
275  hRatio->GetYaxis()->SetDecimals();
276  hRatio->GetYaxis()->SetTitleOffset(0.4);
277  hRatio->GetYaxis()->SetTitle("#nu_{#mu} ND/Sig. FD");
278 
279  if(varName=="recoE"){
280  hRatio->GetXaxis()->SetTitle("Reconstructed neutrino energy (GeV)");
281  }
282  if(varName=="recoQ2"){
283  hRatio->GetXaxis()->SetTitle("ND reco Q^{2} vs. FD true Q^{2} (GeV^{2})");
284  }
285  if(varName=="recoQ2" || varName=="trueQ2"){
286  hRatio->GetXaxis()->SetRangeUser(0,2);
287  }
288 
289  hRatio->GetYaxis()->SetNdivisions(502,kFALSE);
290  hRatio->GetXaxis()->SetNdivisions(405,kFALSE);
291 
292  hRatio -> Draw("hist");
293  lOne -> Draw("same");
294 
295  c2->Print(("plots/FDvsND_"+selName+"_"+varName+".png").c_str());
296 
297  // save further weights in weights file
298  if(varName=="trueQ2"){
299  hWeight->Write(("trueQ2_weight_"+selName).c_str());
300  }
301  if(varName=="PtP"){
302  hWeight->Write(("PtP_weight_"+selName).c_str());
303  }
304  if(varName=="CosNumi"){
305  hWeight->Write(("CosNumi_weight_"+selName).c_str());
306  }
307 
308  } // end for varIdx
309  } // end for selIdx
310 
311  // close weights file
312  file->Close();
313 }
tree Draw("slc.nhit")
Pass neutrinos through unchanged.
Definition: IOscCalc.h:40
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
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
virtual Spectrum Predict(osc::IOscCalc *calc) const =0
Float_t x1[n_points_granero]
Definition: compare.C:5
hmean SetLineStyle(2)
T sqrt(T number)
Definition: d0nt_math.hpp:156
std::string name
Definition: NuePlotLists.h:12
const Color_t kTotalMCErrorBandColor
Definition: Style.h:17
ntuple SetFillStyle(1001)
const Color_t kNueSignalColor
Definition: Style.h:19
c2
Definition: demo5.py:33
void Legend(double x0, double y0, double x1, double y1)
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:86
std::unique_ptr< IPrediction > LoadFrom< IPrediction >(TDirectory *dir, const std::string &label)
Definition: IPrediction.cxx:28
const int kNumSels
Definition: vars.h:43
osc::NoOscillations noosc
Definition: NueSystFuncs.h:44
double POT() const
Definition: Spectrum.h:231
OStream cout
Definition: OStream.cxx:6
void SetDmsq32(const T &dmsq32) override
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
void SetdCP(const T &dCP) override
const hit & b
Definition: hits.cxx:21
TFile * file
Definition: cellShifts.C:17
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
const Color_t kTotalMCColor
Definition: Style.h:16
void plot_NDvsFD_weights(const std::string fname_ND="datamc_ND_numu_kinematics.root", const std::string fname_FD="FDprediction_kinematics.root")
c cd(1)
T asin(T number)
Definition: d0nt_math.hpp:60