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