plot_NDvsFD_weights_RHC.C
Go to the documentation of this file.
1 
4 #include "CAFAna/Core/Binning.h"
5 #include "CAFAna/Core/Cut.h"
6 #include "CAFAna/Core/HistAxis.h"
8 #include "CAFAna/Core/Spectrum.h"
11 #include "CAFAna/Cuts/Cuts.h"
13 #include "CAFAna/Cuts/NueCutsSecondAna.h"
14 #include "CAFAna/Cuts/SpillCuts.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 <iostream>
42 #include <iomanip>
43 #include <fstream>
44 
45 //make a pretty legend
46 void Legend(double x0, double y0, double x1, double y1)
47 {
48  // Doesn't handle log-y well
49  // TLegend* leg = AutoPlaceLegend(x1-x0, y1-y0);
50  TLegend* leg = new TLegend(x0, y0, x1, y1);
51  leg->SetFillStyle(0);
52 
53  TH1* dummy = new TH1F("", "", 1, 0, 1);
54  TH1* dummyFill = new TH1F("", "", 1, 0, 1);
55  dummyFill->SetLineColor(kTotalMCColor);
56  dummyFill->SetFillColor(kTotalMCErrorBandColor);
57  dummy->SetLineColor(kGray+2);
58  dummy->SetMarkerColor(kGray+2);
59  dummy->SetLineColor(1);
60  leg->AddEntry(dummy->Clone(), "ND #bar{#nu}_{#mu} Total MC", "l");
61  dummy->SetLineColor(kNueSignalColor);
62  leg->AddEntry(dummy->Clone(), "FD signal #bar{#nu}_{e} + WS", "l");
63 
64  leg->Draw();
65 }
66 
67 void plot_NDvsFD_weights_RHC(const std::string fname_ND = "datamc_ND_numu_kinematics.root", const std::string fname_FD = "FDprediction_kinematics.root")
68 {
69  struct HistDef
70  {
72  };
73 
74  //which variables do you actually want to look at?
75  const int kNumVars = 5;
76  const HistDef defs[kNumVars] = {
77  "trueQ2",
78  "recoQ2",
79  "trueW2",
80  "PtP",
81  "CosNumi"
82  };
83 
84  //which cut tiers do you want to look at?
85  const int kNumSels = 3;
86  const std::string selNamesND[kNumSels] = {
87  "numuND",
88  "numuND",
89  "numuND"
90  };
91  const std::string selNamesFD[kNumSels] = {
92  "CORE",
93  "PERI",
94  "ALL"
95  };
96 
97  IPrediction* predsND[kNumSels][kNumVars];
98 
100  osc::OscCalcDumb mycalc;
101 
102  //pick up spectra/prediction objects from the input ND root file
103  for(int selIdx = 0; selIdx < kNumSels; ++selIdx){
104  const char* selName = selNamesND[selIdx].c_str();
105  for(int varIdx = 0; varIdx < kNumVars; ++varIdx){
106  const char* varName = defs[varIdx].name.c_str();
107 
108  std::cout<<TString::Format("%s/spect_%s", selName, varName).Data()<<std::endl;
109  predsND[selIdx][varIdx] = LoadFromFile<IPrediction>(fname_ND, TString::Format("%s/pred_%s", selName, varName).Data()).release();
110  if(selIdx == 0 && varIdx == 0)
111  std::cout << "ND MC POT: "
112  << predsND[selIdx][varIdx]->Predict(&noosc).POT() << std::endl;
113 
114  } // end for varIdx
115  } // end for selIdx
116 
117 
118  TFile * file_no = new TFile (fname_FD.c_str(),"READ");
119  TDirectory * dpred_no[kNumSels][kNumVars];
120  for(int selIdx = 0; selIdx < kNumSels; ++selIdx){
121  const std::string selName = selNamesFD[selIdx];
122  for(int varIdx = 0; varIdx < kNumVars; ++varIdx){
123  const std::string varName = defs[varIdx].name;
124 
125  // need to pick equivalent variables
126  if(varName == "recoQ2") {
127  dpred_no[selIdx][varIdx] = file_no->GetDirectory(("prediction/nue_pred_noextrap_"+selName+"_trueQ2").c_str());
128  } else if(varName == "CosNumi"){
129  dpred_no[selIdx][varIdx] = file_no->GetDirectory(("prediction/nue_pred_noextrap_"+selName+"_CosTheta").c_str());
130  } else {
131  dpred_no[selIdx][varIdx] = file_no->GetDirectory(("prediction/nue_pred_noextrap_"+selName+"_"+varName).c_str());
132  }
133  }
134  }
135 
136  //loop over the different selector/cut levels to plot each one
137  for(int selIdx = 0; selIdx < kNumSels; ++selIdx){
138  const std::string selName = selNamesFD[selIdx];
139  //TFile * file = new TFile(("Weights_"+selName+"_RHC.root").c_str(),"RECREATE");
140  for(int varIdx = 0; varIdx < kNumVars; ++varIdx){
141 
142  const std::string varName = defs[varIdx].name;
143 
144  //create a canvas
145  TCanvas *c2 =new TCanvas("c2","c2",500,400);
146 
147  //absolute pad
148  TPad *padAbs = new TPad("padAbs", "padAbs", 0, 0.3, 1, 1.0);
149  padAbs -> SetBottomMargin(0);
150  padAbs->SetFillStyle(0);
151  padAbs -> SetLeftMargin(0.15);
152 
153  padAbs -> Draw();
154  padAbs -> cd();
155 
156  double legendxlow=0.5;
157  double legendxhigh=0.9;
158 
159  if(varName=="trueW2"||varName=="CosNumi"){legendxlow=0.2; legendxhigh=0.5;}
160 
161  //Pull down the variable we want
162 
163  IPrediction* predND = predsND[selIdx][varIdx];
164  auto predFD = ana::LoadFrom<IPrediction> (dpred_no[selIdx][varIdx]);
165 
166  //Get components of the prediction and store them in histograms
167 
168  const double kPOT = predND->Predict(&noosc).POT();
169 
170  TH1* hFD = predFD->PredictComponent(&mycalc,
172  Current::kCC,
173  Sign::kBoth).ToTH1(kPOT);
174 
175  TH1* hAux = predND->Predict(&noosc).ToTH1(kPOT);
176 
177  TH1* hND = predND->Predict(&noosc).ToTH1(kPOT*hFD->Integral()/hAux->Integral());
178 
179  hND->SetLineColor(1);
180  hFD->SetLineColor(kNueSignalColor);
181 
182  std::cout << selName << " " << varName << ": " << hND->Integral() << std::endl;
183 
184  hND->GetXaxis()->CenterTitle();
185  hND->GetYaxis()->CenterTitle();
186  hND->GetXaxis()->SetDecimals();
187  hND->GetYaxis()->SetDecimals();
188  hND->GetYaxis()->SetTitleOffset(1.0);
189  hND->GetXaxis()->SetLabelSize(.0);
190  hND->GetYaxis()->SetLabelSize(.05);
191  hND->SetMaximum(1.1*hND->GetMaximum());
192  hND->GetYaxis()->SetTitle(TString::Format("Events"));
193 
194  if(varName=="recoQ2" || varName=="trueQ2"){
195  hND->GetXaxis()->SetRangeUser(0,2);
196  }
197 
198  hND->GetXaxis()->SetNdivisions(405,kFALSE);
199  hND->SetMinimum(.001);
200 
201  hND ->DrawCopy("hist");
202  hND->Draw("hist same");
203  hFD->Draw("hist same");
204 
205  //Fix up the axis
206  padAbs->RedrawAxis();
207  padAbs->Update();
208 
209  //Add the legend and preliminary
210  Legend(legendxlow, .70, legendxhigh, .85);
211  //Preliminary();
212 
213  TLatex* selTitle = new TLatex(.935, .95, ("nominal "+selName).c_str());
214  selTitle->SetNDC();
215  selTitle->SetTextSize(2/30.);
216  selTitle->SetTextAlign(32);
217  selTitle->Draw();
218 
219  c2->cd();
220 
221  //ratio pad
222  TPad *padRatio = new TPad("padRatio", "padRatio", 0, 0.02, 1, 0.33);
223  padRatio -> SetTopMargin(0.1);
224  padRatio -> SetBottomMargin(0.25);
225  padRatio -> SetFillStyle(0);
226  padRatio -> SetLeftMargin(0.15);
227 
228  padRatio -> Draw();
229  padRatio -> cd();
230 
231  TH1* hRatio = (TH1*)hND -> Clone("hRatio");
232 
233  double b=5;
234  if(varName=="trueW2") {b = 2.0;}
235  if(varName=="trueQ2") {b = 2.0;}
236  if(varName=="recoQ2") {b = 2.0;}
237  if(varName=="PtP" || varName=="CosNumi") {b = 1.0;}
238 
239  TLine* lOne = new TLine(0,1.0,b,1.0);
240  lOne -> SetLineStyle(2);
241 
242  hRatio -> Divide(hFD);
243 
244  TH1* hWeight = (TH1*)hFD -> Clone("hWeight");
245  hWeight -> Divide(hND);
246 
247  hRatio->GetXaxis()->CenterTitle();
248  hRatio->GetYaxis()->CenterTitle();
249  hRatio->GetYaxis()->SetLabelSize(4.0/30.);
250  hRatio->GetXaxis()->SetLabelSize(4.0/30.);
251  hRatio->GetXaxis()->SetLabelOffset(.04);
252  hRatio->GetXaxis()->SetTickSize(.06);
253  hRatio->GetYaxis()->SetTitleSize(3.5/30.);
254  hRatio->GetXaxis()->SetTitleSize(3.5/30.);
255  hRatio->GetYaxis()->SetRangeUser(0,2.0);
256  hRatio->SetMarkerSize(.1);
257  hRatio->GetXaxis()->SetDecimals();
258  hRatio->GetYaxis()->SetDecimals();
259  hRatio->GetYaxis()->SetTitleOffset(0.4);
260  hRatio->GetXaxis()->SetTitleOffset(1.05);
261  hRatio->GetYaxis()->SetTitle("#bar{#nu}_{#mu} ND/Sig.+WS FD");
262 
263  if(varName=="recoQ2"){
264  hRatio->GetXaxis()->SetTitle("ND reco Q^{2} vs. FD true Q^{2} (GeV^{2})");
265  }
266  if(varName=="recoQ2" || varName=="trueQ2"){
267  hRatio->GetXaxis()->SetRangeUser(0,2);
268  }
269 
270  hRatio->GetYaxis()->SetNdivisions(502,kFALSE);
271  hRatio->GetXaxis()->SetNdivisions(405,kFALSE);
272 
273  hRatio -> Draw("hist");
274  lOne -> Draw("same");
275  padRatio -> RedrawAxis();
276  padRatio -> Update();
277 
278  c2->Print(("plots/FDvsND_"+selName+"_"+varName+".png").c_str());
279 
280  // save further weights in weights file
281  /*if(varName=="trueQ2"){
282  hWeight->Write(("trueQ2_weight_"+selName).c_str());
283  }
284  if(varName=="PtP"){
285  hWeight->Write(("PtP_weight_"+selName).c_str());
286  }
287  if(varName=="CosNumi"){
288  hWeight->Write(("CosNumi_weight_"+selName).c_str());
289  }*/
290 
291  } // end for varIdx
292  //file->Close();
293  } // end for selIdx
294 
295  // close weights file
296 
297 }
298 
void plot_NDvsFD_weights_RHC(const std::string fname_ND="datamc_ND_numu_kinematics.root", const std::string fname_FD="FDprediction_kinematics.root")
tree Draw("slc.nhit")
Pass neutrinos through unchanged.
Definition: IOscCalc.h:45
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
c1 Update()
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:148
virtual Spectrum Predict(osc::IOscCalc *calc) const =0
Float_t x1[n_points_granero]
Definition: compare.C:5
hmean SetLineStyle(2)
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
Charged-current interactions.
Definition: IPrediction.h:39
Simple oscillation probability calculator that has no solar term or mass hierarchy or delta so it&#39;s s...
Definition: OscCalcDumb.h:16
const HistDef defs[kNumVars]
Definition: vars.h:15
std::unique_ptr< IPrediction > LoadFrom< IPrediction >(TDirectory *dir, const std::string &label)
Definition: IPrediction.cxx:18
const int kNumSels
Definition: vars.h:43
osc::NoOscillations noosc
Definition: NueSystFuncs.h:44
double POT() const
Definition: Spectrum.h:227
OStream cout
Definition: OStream.cxx:6
const hit & b
Definition: hits.cxx:21
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
c cd(1)
void Legend(double x0, double y0, double x1, double y1)
enum BeamMode string