mre_blessed.C
Go to the documentation of this file.
1 #include <fstream>
2 #include "CAFAna/Core/LoadFromFile.h"
3 #include "CAFAna/Analysis/Plots.h"
4 #include "CAFAna/Core/Spectrum.h"
5 #include "CAFAna/Prediction/IPrediction.h"
6 #include "CAFAna/Analysis/Plots.h"
7 #include "CAFAna/Core/Ratio.h"
8 #include "CAFAna/Vars/Vars.h"
9 #include "3FlavorAna/Vars/MRVars.h"
10 #include "CAFAna/Analysis/Style.h"
11 #include "CAFAna/Analysis/Calcs.h"
12 #include "CAFAna/Cuts/NueCutsSecondAna.h"
13 #include "3FlavorAna/Cuts/NueCuts2020.h"
14 #include "CAFAna/Core/SpectrumLoader.h"
15 #include "CAFAna/Vars/GenieWeights.h"
16 #include "CAFAna/Vars/PPFXWeights.h"
17 #include "CAFAna/Cuts/SpillCuts.h"
18 #include "CAFAna/Vars/XsecTunes.h"
19 
20 #include "StandardRecord/Proxy/SRProxy.h"
21 
22 using namespace ana;
23 
24 #include "OscLib/IOscCalc.h"
25 #include "OscLib/OscCalcDumb.h"
26 #include "OscLib/OscCalcPMNSOpt.h"
27 
28 #include "TCanvas.h"
29 #include "TObject.h"
30 #include "TArrow.h"
31 #include "TColor.h"
32 #include "TFile.h"
33 #include "TLine.h"
34 #include "TH1.h"
35 #include "TText.h"
36 #include "TLegend.h"
37 #include "TLatex.h"
38 #include "TPad.h"
39 #include "TGraphAsymmErrors.h"
40 
41 void Prelim()
42 {
43  TLatex* prelim = new TLatex(.93, .9, "NOvA Preliminary");
44  prelim->SetTextColor(kBlue);
45  prelim->SetNDC();
46  prelim->SetTextSize(2/30.);
47  prelim->SetTextAngle(270);
48  prelim->SetTextAlign(12);
49  prelim->Draw();
50 }
51 
52 TCanvas *MakeMePlot(std::vector<TH1*> Histos, std::vector<TString> Names, bool isFHC, string yname = "none"){
53 
54  TCanvas *c = new TCanvas("c", "canvas", 800, 700);
55  TPad* p1 = new TPad("", "", 0, 0, 1, 1);
56  TPad* p2 = new TPad("", "", 0, 0, 1, 1);
57  p1->SetBottomMargin(.3);
58  p2->SetTopMargin(.7);
59 
60  for(auto p:{p1,p2}){
61  p->SetLeftMargin(0.1);
62  p->SetFillStyle(0);
63  p->Draw();
64  }
65  p1->cd();
66 
67  auto denom = (TH1*)Histos[0]->Clone();
68  denom->SetMaximum(denom->GetMaximum()*1.4);
69  if (yname == "Selection efficiency") denom->SetMaximum(1.35);
70  denom->SetMinimum(0);
71  denom->GetYaxis()->CenterTitle();
72  denom->GetXaxis()->CenterTitle();
73  denom->GetXaxis()->SetLabelOffset(999);
74  if(yname !="none") denom->GetYaxis()->SetTitle(yname.c_str());
75 
76  denom->Draw("hist E1");
77  for(int i=1; i<Histos.size(); i++){
78  Histos[i]->Draw("hist E1 same");
79  }
80  denom->Draw("hist E1 same");
81  p1->RedrawAxis();
82  gPad->Update();
83  TText *tt;
84  if(isFHC){
85  auto *tt = new TText(gPad->GetUxmax()*0.7, 1.02*gPad->GetUymax(), "FHC beam");
86  tt->SetTextColor(kBlue);
87  tt->SetTextSize(0.05);
88  tt->Draw();}
89  else{
90  auto *tt = new TText(gPad->GetUxmax()*0.7, 1.02*gPad->GetUymax(), "RHC beam");
91  tt->SetTextColor(kBlue);
92  tt->SetTextSize(0.05);
93  tt->Draw();}
94 
95  TLegend *leg = new TLegend(0.5,0.76,0.95,0.87);
96  leg->SetTextSize(0.025);
97  for(int i=0; i<Histos.size(); i++) leg->AddEntry(Histos[i], Names[i], "l");
98  leg->Draw();
99 
100  gPad->RedrawAxis();
101  p2->cd();
102  auto temp = (TH1*)Histos[1]->Clone();
103  temp->Divide(denom);
104  temp->Draw();
105  for(int i=2; i<Histos.size(); i++){
106  auto temp = (TH1*)Histos[i]->Clone();
107  temp->Divide(denom);
108  temp->Draw("same");
109  }
110  temp->GetYaxis()->SetRangeUser(0.8, 1.25);
111  temp->GetYaxis()->SetTitle("Ratio");
112  TLine *l = new TLine(-0.1,1,4.5,1);
113  l->SetLineColor(kBlack);
114  l->Draw();
115  p2->RedrawAxis();
116 
117  return c;
118 
119 }
120 
121 
122 void mre_blessed(bool plot = true, bool isFHC = true, string weight = "", TString drawing_options = "red"){
123 
124  const std::string mode_tag = isFHC?"fhc":"rhc";
125  std::string file_mc = "mre_pred_nd_cuts_"+mode_tag+"_mc.root";
126 
127  const Var kSlcE([](const caf::SRProxy *sr){
128  return sr->slc.calE;
129  });
130 
131  const Var kCalE1erProng([](const caf::SRProxy *sr){
132  if ( !sr->vtx.elastic.IsValid ) return -1000.f;
133  if (sr->vtx.elastic.fuzzyk.npng < 1) return -1000.f;
134  return float(sr->vtx.elastic.fuzzyk.png[0].calE);});
135 
136  const Var kHadCalENo1erProng = kSlcE - kCalE1erProng;
137 
138  const int kNumVars = 5;
139 
140  struct HistDef
141  {
142  std::string name;
143  HistAxis axis;
144  int lowedge;
145  int highedge;
146  int rebin;
147  };
148 
149  const HistDef defs[kNumVars] = {
150  {"energy", {"#nu_{e} Energy (GeV)", Binning::Simple(120,0,12), kNueEnergy2020}, 1, 45, 1},
151  {"caloE", {"Calorimetric energy (GeV)", Binning::Simple(120,0,12), kSlcE}, 1, 45, 1},
152  {"caloE_1pr", {"Calorimetric energy, 1st prong", Binning::Simple(120,0,12), kCalE1erProng}, 1, 45, 1},
153  {"calE_no1pr", {"Calorimetric energy, not 1st prong", Binning::Simple(120,0,12), kHadCalENo1erProng}, 1, 45, 1},
154  {"trueE", {"True Energy (GeV)", Binning::Simple(120,0,12), kMREParentNumuE}, 1, 45, 1},
155  //{"ResE", {"Energy Resolution", Binning::Simple(100, -1, 1), (kMREParentNumuE-kNueEnergy2020)/kMREParentNumuE}, 1, 100, 1},
156  };
157 
158  const int kNumSels = 2;
159 
160  const Cut sels[kNumSels] = {
161  kNue2020NDPresel && kNue2017MRParentSliceCut,
162  kNue2020ND && kNue2017MRParentSliceCut
163  //kNue2020NDFiducial && kNue2020NDContain && kNue2020NDNHits && kNue2020NDEnergy && kNue2020NDProngLength && kNue2017MRParentSliceCut,
164  //kNue2020NDFiducial && kNue2020NDContain && kNue2020NDNHits && kNue2020NDEnergy && kNue2020NDProngLength && kNue2017MRParentSliceCut && kNue2020PID
165  };
166 
167  const std::string selNames[kNumSels] = {
168  "nd-presel-mre",
169  "nd-full-mre"
170  //"nd-presel-wo-dq-mre",
171  //"nd-full-wo-dq-mre",
172  };
173 
174  if(!plot){
175 
176  std::string fname_data = "";
177  std::string fname = "";
178  std::string fname_calibxy_pos = "";
179  std::string fname_calibxy_neg = "";
180  std::string fname_calib = "";
181  std::string fname_light_up = "";
182  std::string fname_light_down = "";
183  std::string fname_ckv = "";
184 
185  if(isFHC){
186  fname_data = "defname: prod_mrecaf_R19-11-18-prod5reco.r_nd_numi_fhc_full_v1";
187  fname = "defname: prod_mrecaf_R19-11-18-prod5reco.r_nd_genie_N1810j0211a_nonswap_fhc_nova_v08_full_v1";
188  fname_calibxy_pos = "defname: prod_mrecaf_R19-11-18-prod5reco.r_nd_genie_N1810j0211a_nonswap_fhc_nova_v08_full_calibup_v1 with stride 1";
189  fname_calibxy_neg = "defname: prod_mrecaf_R19-11-18-prod5reco.r_nd_genie_N1810j0211a_nonswap_fhc_nova_v08_full_calibdown_v1 with stride 1";
190  fname_calib = "defname: prod_mrecaf_R19-11-18-prod5reco.r_nd_genie_N1810j0211a_nonswap_fhc_nova_v08_full_calibshape_v1 with stride 1";
191  fname_light_up = "defname: prod_mrecaf_R19-11-18-prod5reco.r_nd_genie_N1810j0211a_nonswap_fhc_nova_v08_full_llup_v1 with stride 1";
192  fname_light_down = "defname: prod_mrecaf_R19-11-18-prod5reco.r_nd_genie_N1810j0211a_nonswap_fhc_nova_v08_full_lldown_v1 with stride 1";
193  fname_ckv = "defname: prod_mrecaf_R19-11-18-prod5reco.r_nd_genie_N1810j0211a_nonswap_fhc_nova_v08_full_cherenkov_v1";
194 
195  }
196  if(!isFHC) {
197  fname_data = "defname: prod_mrecaf_R19-11-18-prod5reco.r_nd_numi_rhc_full_v1";
198  fname = "defname: prod_mrecaf_R19-11-18-prod5reco.r_nd_genie_N1810j0211a_nonswap_rhc_nova_v08_full_v1";
199  fname_calibxy_pos = "defname: prod_mrecaf_R19-11-18-prod5reco.r_nd_genie_N1810j0211a_nonswap_rhc_nova_v08_full_calibup_v1 with stride 1";
200  fname_calibxy_neg = "defname: prod_mrecaf_R19-11-18-prod5reco.r_nd_genie_N1810j0211a_nonswap_rhc_nova_v08_full_calibdown_v1 with stride 1";
201  fname_calib = "defname: prod_mrecaf_R19-11-18-prod5reco.r_nd_genie_N1810j0211a_nonswap_rhc_nova_v08_full_calibshape_v1 with stride 1";
202  fname_light_up = "defname: prod_mrecaf_R19-11-18-prod5reco.r_nd_genie_N1810j0211a_nonswap_rhc_nova_v08_full_llup_v1 with stride 1";
203  fname_light_down = "defname: prod_mrecaf_R19-11-18-prod5reco.r_nd_genie_N1810j0211a_nonswap_rhc_nova_v08_full_lldown_v1 with stride 1";
204  fname_ckv = "defname: prod_mrecaf_R19-11-18-prod5reco.r_nd_genie_N1810j0211a_nonswap_rhc_nova_v08_full_cherenkov_v1 with stride 1";
205  }
206 
207  SpectrumLoader loader_data(fname_data);
208  SpectrumLoader loader(fname);
209  SpectrumLoader loader_calibxy_pos(fname_calibxy_pos);
210  SpectrumLoader loader_calibxy_neg(fname_calibxy_neg);
211  SpectrumLoader loader_calib(fname_calib);
212  SpectrumLoader loader_light_up(fname_light_up);
213  SpectrumLoader loader_light_down(fname_light_down);
214  SpectrumLoader loader_ckv(fname_ckv);
215 
216  loader_data.SetSpillCut(kStandardSpillCuts);
217  loader.SetSpillCut(kStandardSpillCuts);
218  loader_calibxy_pos.SetSpillCut(kStandardSpillCuts);
219  loader_calibxy_neg.SetSpillCut(kStandardSpillCuts);
220  loader_calib.SetSpillCut(kStandardSpillCuts);
221  loader_light_up.SetSpillCut(kStandardSpillCuts);
222  loader_light_down.SetSpillCut(kStandardSpillCuts);
223  loader_ckv.SetSpillCut(kStandardSpillCuts);
224 
225  Spectrum* spec_data[kNumSels][kNumVars];
226  Spectrum* spec[kNumSels][kNumVars];
227  //Spectrum* spec_beam[kNumSels][kNumVars];
228  //Spectrum* spec_xsec[kNumSels][kNumVars];
229  //Spectrum* spec_now[kNumSels][kNumVars];
230  Spectrum* spec_calibxy_pos[kNumSels][kNumVars];
231  Spectrum* spec_calibxy_neg[kNumSels][kNumVars];
232  Spectrum* spec_calib[kNumSels][kNumVars];
233  Spectrum* spec_light_up[kNumSels][kNumVars];
234  Spectrum* spec_light_down[kNumSels][kNumVars];
235  Spectrum* spec_ckv[kNumSels][kNumVars];
236 
237  for(int selIdx = 0; selIdx < kNumSels; ++selIdx){
238  for (int varIdx=0; varIdx<kNumVars;varIdx++){
239  const HistAxis& axis = defs[varIdx].axis;
240  spec_data[selIdx][varIdx] = new Spectrum(loader_data, axis, sels[selIdx], kNoShift);
241  spec[selIdx][varIdx] = new Spectrum(loader, axis, sels[selIdx], kNoShift, kPPFXFluxCVWgt*kXSecCVWgt2020);
242  //spec_beam[selIdx][varIdx] = new Spectrum(loader, axis, sels[selIdx], kNoShift, kPPFXFluxCVWgt);
243  //spec_xsec[selIdx][varIdx] = new Spectrum(loader, axis, sels[selIdx], kNoShift, kXSecCVWgt2020);
244  //spec_now[selIdx][varIdx] = new Spectrum(loader, axis, sels[selIdx], kNoShift);
245  spec_calibxy_pos[selIdx][varIdx] = new Spectrum(loader_calibxy_pos, axis, sels[selIdx], kNoShift, kPPFXFluxCVWgt*kXSecCVWgt2020);
246  spec_calibxy_neg[selIdx][varIdx] = new Spectrum(loader_calibxy_neg, axis, sels[selIdx], kNoShift, kPPFXFluxCVWgt*kXSecCVWgt2020);
247  spec_calib[selIdx][varIdx] = new Spectrum(loader_calib, axis, sels[selIdx], kNoShift, kPPFXFluxCVWgt*kXSecCVWgt2020);
248  spec_light_up[selIdx][varIdx] = new Spectrum(loader_light_up, axis, sels[selIdx], kNoShift, kPPFXFluxCVWgt*kXSecCVWgt2020);
249  spec_light_down[selIdx][varIdx] = new Spectrum(loader_light_down, axis, sels[selIdx], kNoShift, kPPFXFluxCVWgt*kXSecCVWgt2020);
250  spec_ckv[selIdx][varIdx] = new Spectrum(loader_ckv, axis, sels[selIdx], kNoShift, kPPFXFluxCVWgt*kXSecCVWgt2020);
251  }
252  }
253 
254  loader_data.Go();
255  loader.Go();
256  loader_calibxy_pos.Go();
257  loader_calibxy_neg.Go();
258  loader_calib.Go();
259  loader_light_up.Go();
260  loader_light_down.Go();
261  loader_ckv.Go();
262 
263  TFile* file = new TFile(file_mc.c_str(),"recreate");
264 
265  for(int selIdx = 0; selIdx < kNumSels; ++selIdx){
266  TDirectory* d = file->mkdir(selNames[selIdx].c_str());
267  for (int varIdx=0; varIdx<kNumVars;varIdx++){
268  const char* name = defs[varIdx].name.c_str();
269  spec_data[selIdx][varIdx]->SaveTo(d, TString::Format("data_%s", name));
270  spec[selIdx][varIdx]->SaveTo(d, TString::Format("%s", name));
271  //spec_beam[selIdx][varIdx]->SaveTo(d, TString::Format("beam_%s", name));
272  //spec_xsec[selIdx][varIdx]->SaveTo(d, TString::Format("xsec_%s", name));
273  //spec_now[selIdx][varIdx]->SaveTo(d, TString::Format("now_%s", name));
274  spec_calibxy_pos[selIdx][varIdx]->SaveTo(d, TString::Format("calibxy_pos_%s", name));
275  spec_calibxy_neg[selIdx][varIdx]->SaveTo(d, TString::Format("calibxy_neg_%s", name));
276  spec_calib[selIdx][varIdx]->SaveTo(d, TString::Format("calib_%s", name));
277  spec_light_up[selIdx][varIdx]->SaveTo(d, TString::Format("light_up_%s", name));
278  spec_light_down[selIdx][varIdx]->SaveTo(d, TString::Format("light_down_%s", name));
279  spec_ckv[selIdx][varIdx]->SaveTo(d, TString::Format("ckv_%s", name));
280 
281  }
282  }
283  file->Close();
284  }
285 
286  if(plot){
287 
288  Spectrum *data[kNumSels][kNumVars];
289  Spectrum *mc[kNumSels][kNumVars];
290  Spectrum *mc_calibxy_pos[kNumSels][kNumVars];
291  Spectrum *mc_calibxy_neg[kNumSels][kNumVars];
292  Spectrum *mc_calib[kNumSels][kNumVars];
293  Spectrum *mc_light_up[kNumSels][kNumVars];
294  Spectrum *mc_light_down[kNumSels][kNumVars];
295  Spectrum *mc_ckv[kNumSels][kNumVars];
296 
297  for(int selIdx = 0; selIdx < kNumSels; ++selIdx){
298  const char* selName = selNames[selIdx].c_str();
299  for(int varIdx = 0; varIdx < kNumVars; ++varIdx){
300  const char* varName = defs[varIdx].name.c_str();
301  data[selIdx][varIdx] = LoadFromFile<Spectrum>(file_mc, TString::Format("%s/data_%s", selName, varName).Data()).release();
302  mc[selIdx][varIdx] = LoadFromFile<Spectrum>(file_mc, TString::Format("%s/%s%s", selName, weight.c_str(), varName).Data()).release();
303  mc_calibxy_pos[selIdx][varIdx] = LoadFromFile<Spectrum>(file_mc, TString::Format("%s/calibxy_pos_%s", selName, varName).Data()).release();
304  mc_calibxy_neg[selIdx][varIdx] = LoadFromFile<Spectrum>(file_mc, TString::Format("%s/calibxy_neg_%s", selName, varName).Data()).release();
305  mc_calib[selIdx][varIdx] = LoadFromFile<Spectrum>(file_mc, TString::Format("%s/calib_%s", selName, varName).Data()).release();
306  mc_light_up[selIdx][varIdx] = LoadFromFile<Spectrum>(file_mc, TString::Format("%s/light_up_%s", selName, varName).Data()).release();
307  mc_light_down[selIdx][varIdx] = LoadFromFile<Spectrum>(file_mc, TString::Format("%s/light_down_%s", selName, varName).Data()).release();
308  mc_ckv[selIdx][varIdx] = LoadFromFile<Spectrum>(file_mc, TString::Format("%s/ckv_%s", selName, varName).Data()).release();
309 
310  } // end for varIdx
311  } // end for selIdx
312 
313  TH1* hData[kNumSels][kNumVars];
314  TH1* hMC[kNumSels][kNumVars];
315  TH1* hMC_calibxy_pos[kNumSels][kNumVars];
316  TH1* hMC_calibxy_neg[kNumSels][kNumVars];
317  TH1* hMC_calib[kNumSels][kNumVars];
318  TH1* hMC_light_up[kNumSels][kNumVars];
319  TH1* hMC_light_down[kNumSels][kNumVars];
320  TH1* hMC_ckv[kNumSels][kNumVars];
321 
322  double POTmc = data[0][0]->POT();
323  cout<<POTmc<<endl;
324 
325  for(int selIdx = 0; selIdx < 2; ++selIdx){
326  const char* selName = selNames[selIdx].c_str();
327  for(int varIdx = 0; varIdx < kNumVars; ++varIdx){
328  const char* varName = defs[varIdx].name.c_str();
329  hData[selIdx][varIdx]=data[selIdx][varIdx]->ToTH1(POTmc, kBlack, kSolid, ana::EExposureType::kPOT);
330  hMC[selIdx][varIdx]=mc[selIdx][varIdx]->ToTH1(POTmc, kBlack, kSolid, ana::EExposureType::kPOT);
331  hMC_calibxy_pos[selIdx][varIdx]=mc_calibxy_pos[selIdx][varIdx]->ToTH1(POTmc, kBlue, kSolid, ana::EExposureType::kPOT);
332  hMC_calibxy_neg[selIdx][varIdx]=mc_calibxy_neg[selIdx][varIdx]->ToTH1(POTmc, kRed, kSolid, ana::EExposureType::kPOT);
333  hMC_calib[selIdx][varIdx]=mc_calib[selIdx][varIdx]->ToTH1(POTmc, kGreen+2, kSolid, ana::EExposureType::kPOT);
334  hMC_light_up[selIdx][varIdx]=mc_light_up[selIdx][varIdx]->ToTH1(POTmc, kBlue, kSolid, ana::EExposureType::kPOT);
335  hMC_light_down[selIdx][varIdx]=mc_light_down[selIdx][varIdx]->ToTH1(POTmc, kRed, kSolid, ana::EExposureType::kPOT);
336  hMC_ckv[selIdx][varIdx]=mc_ckv[selIdx][varIdx]->ToTH1(POTmc, kGreen+2, kSolid, ana::EExposureType::kPOT);
337 
338 
339  for(auto &h:{hData[selIdx][varIdx], hMC[selIdx][varIdx], hMC_calibxy_pos[selIdx][varIdx], hMC_calibxy_neg[selIdx][varIdx], hMC_calib[selIdx][varIdx], hMC_light_up[selIdx][varIdx], hMC_light_down[selIdx][varIdx], hMC_ckv[selIdx][varIdx]}){
340  //for(auto &h:{hData[selIdx][varIdx], hMC[selIdx][varIdx], hMC_calibxy_pos[selIdx][varIdx], hMC_calibxy_neg[selIdx][varIdx], hMC_calib[selIdx][varIdx], hMC_light_up[selIdx][varIdx], hMC_light_down[selIdx][varIdx]}){
341  if(defs[varIdx].rebin>1) h->Rebin(defs[varIdx].rebin);
342  h->GetXaxis()->SetRange(defs[varIdx].lowedge, defs[varIdx].highedge/defs[varIdx].rebin);
343  }
344 
345  double nom_mean = hMC[selIdx][varIdx]->GetMean();
346  double pos_mean = hMC_calibxy_pos[selIdx][varIdx]->GetMean();
347  double neg_mean = hMC_calibxy_neg[selIdx][varIdx]->GetMean();
348  double shape_mean = hMC_calib[selIdx][varIdx]->GetMean();
349  string nom_name = "Nominal (mean "+std::to_string(nom_mean)+")";
350  string pos_name = "Calib XY pos (mean "+std::to_string(pos_mean)+")";
351  string neg_name = "Calib XY neg (mean "+std::to_string(neg_mean)+")";
352  string shape_name = "Calib shape (mean "+std::to_string(shape_mean)+")";
353 
354  auto c = MakeMePlot({hMC[selIdx][varIdx], hMC_calibxy_pos[selIdx][varIdx], hMC_calibxy_neg[selIdx][varIdx], hMC_calib[selIdx][varIdx]}, {nom_name, pos_name, neg_name, shape_name}, isFHC);
355  c->SaveAs(("blpl/"+mode_tag+"_plot-calib-"+selName+"-"+varName+".pdf").c_str());
356  delete c;
357 
358  auto c1 = MakeMePlot({hMC[selIdx][varIdx], hMC_light_up[selIdx][varIdx], hMC_light_down[selIdx][varIdx], hMC_ckv[selIdx][varIdx]}, {"Nominal", "Light up", "Light down", "Cherenkov"}, isFHC);
359  c1->SaveAs(("blpl/"+mode_tag+"_plot-light-ckv-"+selName+"-"+varName+".pdf").c_str());
360  delete c1;
361 
362  if(selIdx>0 && (selIdx%2!=0)){
363 
364  TH1* hEff_mc = (TH1*)hMC[selIdx][varIdx]->Clone();
365 
366  TH1* hEff_mc_denom = (TH1*)hMC[selIdx-1][varIdx]->Clone();
367  TH1* hEff_data = (TH1*)hData[selIdx][varIdx]->Clone();
368  TH1* hEff_data_denom = (TH1*)hData[selIdx-1][varIdx]->Clone();
369  TH1* hEff_mc_calibxy_pos = (TH1*)hMC_calibxy_pos[selIdx][varIdx]->Clone();
370  TH1* hEff_mc_calibxy_pos_denom = (TH1*)hMC_calibxy_pos[selIdx-1][varIdx]->Clone();
371  TH1* hEff_mc_calibxy_neg = (TH1*)hMC_calibxy_neg[selIdx][varIdx]->Clone();
372  TH1* hEff_mc_calibxy_neg_denom = (TH1*)hMC_calibxy_neg[selIdx-1][varIdx]->Clone();
373  TH1* hEff_mc_calib = (TH1*)hMC_calib[selIdx][varIdx]->Clone();
374  TH1* hEff_mc_calib_denom = (TH1*)hMC_calib[selIdx-1][varIdx]->Clone();
375  TH1* hEff_mc_light_up = (TH1*)hMC_light_up[selIdx][varIdx]->Clone();
376  TH1* hEff_mc_light_up_denom = (TH1*)hMC_light_up[selIdx-1][varIdx]->Clone();
377  TH1* hEff_mc_light_down = (TH1*)hMC_light_down[selIdx][varIdx]->Clone();
378  TH1* hEff_mc_light_down_denom = (TH1*)hMC_light_down[selIdx-1][varIdx]->Clone();
379  TH1* hEff_mc_ckv = (TH1*)hMC_ckv[selIdx][varIdx]->Clone();
380  TH1* hEff_mc_ckv_denom = (TH1*)hMC_ckv[selIdx-1][varIdx]->Clone();
381 
382  hEff_mc->Divide(hEff_mc_denom);
383  hEff_data->Divide(hEff_data_denom);
384  hEff_mc_calibxy_pos->Divide(hEff_mc_calibxy_pos_denom);
385  hEff_mc_calibxy_neg->Divide(hEff_mc_calibxy_neg_denom);
386  hEff_mc_calib->Divide(hEff_mc_calib_denom);
387  hEff_mc_light_up->Divide(hEff_mc_light_up_denom);
388  hEff_mc_light_down->Divide(hEff_mc_light_down_denom);
389  hEff_mc_ckv->Divide(hEff_mc_ckv_denom);
390 
391  auto c2 = MakeMePlot({hEff_mc, hEff_mc_calibxy_pos, hEff_mc_calibxy_neg, hEff_mc_calib}, {"Nominal", "Calib XY pos", "Calib XY neg", "Calib shape"}, isFHC, "Selection efficiency");
392  c2->SaveAs(("blpl/"+mode_tag+"_plot-efficiency-calib-"+selName+"-"+varName+".pdf").c_str());
393  delete c2;
394 
395  auto c3 = MakeMePlot({hEff_mc, hEff_mc_light_up, hEff_mc_light_down, hEff_mc_ckv}, {"Nominal", "Light up", "Light down", "Cherenkov"}, isFHC, "Selection efficiency");
396  c3->SaveAs(("blpl/"+mode_tag+"_plot-efficiency-light-"+selName+"-"+varName+".pdf").c_str());
397  delete c3;
398 
399  hEff_mc_light_up->SetLineColor(kViolet);
400  hEff_mc_light_down->SetLineColor(kPink);
401  hEff_mc_calib->SetLineColor(kOrange);
402  auto c4 = MakeMePlot({hEff_mc, hEff_mc_light_up, hEff_mc_light_down, hEff_mc_ckv, hEff_mc_calibxy_pos, hEff_mc_calibxy_neg, hEff_mc_calib}, {"Nominal", "Light up", "Light down", "Cherenkov", "Calib XY pos", "Calib XY neg", "Calib shape"}, isFHC, "Selection efficiency");
403  //c4->SaveAs(("blpl/"+mode_tag+"_plot-efficiency-"+selName+"-"+varName+".pdf").c_str());
404  delete c4;
405 
406  auto tmp = MakeMePlot({hEff_data, hEff_mc}, {"Data", "MC"}, isFHC, "Selection efficiency");
407  tmp->SaveAs(("blpl/"+weight+mode_tag+"_plot-efficiency-"+selName+"-"+varName+".pdf").c_str());
408  delete tmp;
409 
410  // blessing plot - Data vs MC eff band + syst. band
411 
412  TCanvas *c5 = new TCanvas("c5", "canvas", 1000, 1000);
413  TPad* p1 = new TPad("", "", 0, 0, 1, 1);
414  TPad* p2 = new TPad("", "", 0, 0, 1, 1);
415  p1->SetBottomMargin(.3);
416  p2->SetTopMargin(.7);
417 
418  for(auto p:{p1,p2}){
419  p->SetLeftMargin(0.105);
420  p->SetFillStyle(0);
421  p->Draw();
422  }
423  p1->cd();
424 
425  /*Int_t ci = TColor::GetFreeColorIndex();
426  TColor *color = new TColor(ci, 155/255.0, 233/255.0, 233/255.0);
427  Int_t line = TColor::GetFreeColorIndex();
428  TColor *color_line = new TColor(line, 0/255.0, 39/255.0, 34/255.0);
429  Int_t line_data = TColor::GetFreeColorIndex();
430  TColor *color_data = new TColor(line_data, 255/255.0, 13/255.0, 109/255.0);*/ // too cheerful
431  Int_t ci = TColor::GetFreeColorIndex();
432  TColor *color = new TColor(ci, 216/255.0, 216/255.0, 216/255.0); // lighter gray
433  //ci = kGray;
434  Int_t line = kGray+2;
435  Int_t line_data = kBlack;
436  if(drawing_options.Contains("red")) {
437  ci= kRed-10;
438  line = kPink;//kRed;
439  line_data = kBlack;}
440 
441  hEff_data->SetLineColor(line_data);
442  hEff_mc->SetLineColor(line);
443  hEff_data->SetLineWidth(2);
444  hEff_mc->SetLineWidth(2);
445  hEff_mc->SetMaximum(1.2);
446  hEff_mc->SetMinimum(0);
447  hEff_mc->GetXaxis()->SetRange(0, 45);
448  hEff_mc->GetYaxis()->CenterTitle();
449  hEff_mc->GetXaxis()->CenterTitle();
450  hEff_mc->GetYaxis()->SetTitle("Selection efficiency");
451  hEff_mc->Draw("hist E1");
452  std::vector<TH1*> syst = {hEff_mc_light_up, hEff_mc_light_down, hEff_mc_ckv, hEff_mc_calibxy_pos, hEff_mc_calibxy_neg, hEff_mc_calib};
453  TGraphAsymmErrors* g = new TGraphAsymmErrors;
454  TH1F* hUp = (TH1F*)hEff_mc->Clone();
455  TH1F* hDown = (TH1F*)hEff_mc->Clone();
456  for(int binIdx = 0; binIdx < hEff_mc->GetNbinsX()+2; ++binIdx){
457  const double y = hEff_mc->GetBinContent(binIdx);
458  g->SetPoint(binIdx, hEff_mc->GetXaxis()->GetBinCenter(binIdx), y);
459  const double w = hEff_mc->GetXaxis()->GetBinWidth(binIdx);
460  double errUp = 0;
461  double errDn = 0;
462  for(unsigned int systIdx = 0; systIdx < syst.size(); ++systIdx){
463  double shift = syst[systIdx]->GetBinContent(binIdx)-y;
464  if(shift < 0) errDn += shift*shift;
465  if(shift >= 0) errUp+= shift*shift;
466  }
467  g->SetPointError(binIdx, w/2, w/2, sqrt(errDn), sqrt(errUp));
468  hUp->SetBinContent(binIdx, hEff_mc->GetBinContent(binIdx)+sqrt(errUp));
469  hDown->SetBinContent(binIdx, hEff_mc->GetBinContent(binIdx) - sqrt(errDn));
470  }
471  g->SetFillColor(ci);
472  g->Draw("e2 same");
473  hEff_mc->SetLineWidth(3);
474  hEff_mc->Draw("hist E1 same");
475  hEff_data->SetMarkerStyle(8);
476  hEff_data->SetMarkerSize(2);
477  hEff_data->Draw("E1 same");
478  c5->RedrawAxis();
479  gPad->Update();
480  TText *tt;
481  double systname_x = 0.05;
482  double leg_x = 0.63;
483  double leg_x_max = 0.85;
484  if (varIdx == 0){
485  systname_x = 0.55;
486  leg_x = 0.15;
487  leg_x_max = 0.38;
488  if (!isFHC){
489  leg_x = 0.15;
490  leg_x_max = 0.33;
491  }
492  }
493  if(isFHC){
494  auto *tt = new TText(gPad->GetUxmax()*0.55, 1.05*gPad->GetUymax(), "Neutrino beam, MRE");
495  tt->SetTextSize(0.04);
496  tt->Draw();
497  tt = new TText(gPad->GetUxmax()*systname_x, 0.9*gPad->GetUymax(), "Calibration + lightlevel syst.");
498  tt->SetTextSize(0.03);
499  tt->Draw();}
500  else{
501  auto *tt = new TText(gPad->GetUxmax()*0.5, 1.05*gPad->GetUymax(), "Antineutrino beam, MRE");
502  tt->SetTextSize(0.04);
503  tt->Draw();
504  tt = new TText(gPad->GetUxmax()*systname_x, 0.9*gPad->GetUymax(), "Calibration + lightlevel syst.");
505  tt->SetTextSize(0.03);
506  tt->Draw();}
507 
508  TLegend *leg = new TLegend(leg_x,0.77,leg_x_max,0.87);
509  leg->SetTextSize(0.033);
510  leg->AddEntry(hEff_data, "Data", "p");
511  leg->AddEntry(hEff_mc, "Monte Carlo", "l");
512  leg->AddEntry(g, "1#sigma syst. error", "f");
513  leg->Draw();
514 
515  Prelim();
516 
517  gPad->RedrawAxis();
518 
519  p2->cd();
520  auto temp = (TH1*)hEff_mc->Clone();
521  temp->Divide(hEff_data);
522  temp->Draw();
523  auto temp_up = (TH1*)hUp->Clone();
524  auto temp_down = (TH1*)hDown->Clone();
525  temp_up->Divide(hEff_data);
526  temp_down->Divide(hEff_data);
527 
528  TGraphAsymmErrors* g_ratio = new TGraphAsymmErrors;
529  for(int binIdx = 0; binIdx < temp->GetNbinsX()+2; ++binIdx){
530  const double w = temp->GetXaxis()->GetBinWidth(binIdx);
531  g_ratio->SetPoint(binIdx, temp->GetXaxis()->GetBinCenter(binIdx), temp->GetBinContent(binIdx));
532  g_ratio->SetPointError(binIdx, w/2, w/2, temp->GetBinContent(binIdx) - temp_down->GetBinContent(binIdx), temp_up->GetBinContent(binIdx) - temp->GetBinContent(binIdx));
533  }
534 
535  g_ratio->SetFillColor(ci);
536  g_ratio->Draw("e2 same");
537  temp->Draw("same");
538  TLine *l = new TLine(-0.1,1,4.5,1);
539  l->SetLineColor(kBlack);
540  l->Draw();
541  temp->GetYaxis()->SetRangeUser(0.5, 1.5);
542  temp->GetYaxis()->SetTitle("MC / Data");
543  p2->RedrawAxis();
544 
545  c5->SaveAs(("blpl/"+mode_tag+"_plot-efficiency-error-band-"+selName+"-"+varName+"_"+drawing_options+".pdf"));
546 
547  } // effic
548 
549  } // var
550  } // sel
551 
552  cout<<"Table for "<<mode_tag<<endl;
553  cout<<"Data "<<hData[0][0]->Integral()<<" "<<hData[1][0]->Integral()<<" "<<hData[1][0]->Integral()/hData[0][0]->Integral()<<endl;
554  cout<<"MC "<<hMC[0][0]->Integral()<<" "<<hMC[1][0]->Integral()<<" "<<hMC[1][0]->Integral()/hMC[0][0]->Integral()<<endl;
555  cout<<"difference "<<(hMC[1][0]->Integral()/hMC[0][0]->Integral() - hData[1][0]->Integral()/hData[0][0]->Integral())/(hMC[1][0]->Integral()/hMC[0][0]->Integral())<<endl;
556 
557  // bl plot with Data vs MC at presel and full sel.
558 
559  for(int varIdx = 0; varIdx < kNumVars; ++varIdx){
560  const char* varName = defs[varIdx].name.c_str();
561  TCanvas *c6 = new TCanvas("c6", "canvas", 1000, 800);
562  c6->SetLeftMargin(0.105);
563 
564  Int_t ci = TColor::GetFreeColorIndex();
565  TColor *color = new TColor(ci, 216/255.0, 216/255.0, 216/255.0); // lighter gray
566  //ci = kGray;
567  Int_t line = kGray+2;
568  Int_t line_data = kBlack;
569  if(drawing_options.Contains("red")) {
570  ci= kRed-10;
571  line = kPink;//kRed;
572  line_data = kBlack;
573  }
574 
575  auto hMC_presel = hMC[0][varIdx];
576  auto hMC_full = hMC[1][varIdx];
577  auto hData_presel = hData[0][varIdx];
578  auto hData_full = hData[1][varIdx];
579 
580  hMC_presel->SetLineStyle(7);
581  hData_presel->SetLineStyle(7);
582  hData_presel->SetLineColor(line);
583  hData_full->SetLineColor(line_data);
584  hMC_presel->SetLineColor(line);
585  hMC_full->SetLineColor(line_data);
586  for (auto &h:{hMC_presel, hMC_full, hData_presel, hData_full}){
587  h->SetLineWidth(2);
588  h->Scale(1/10000.0);
589  }
590  hData_presel->SetMarkerStyle(8);
591  hData_presel->SetMarkerColor(line);
592  hData_presel->SetMarkerSize(2);
593  hData_full->SetMarkerStyle(8);
594  hData_full->SetMarkerSize(2);
595  hData_presel->SetMaximum(hData_presel->GetMaximum()*1.2);
596  hData_presel->SetMinimum(0);
597  hMC_presel->GetXaxis()->SetRange(defs[varIdx].lowedge, defs[varIdx].highedge/defs[varIdx].rebin);;
598  hMC_presel->GetYaxis()->CenterTitle();
599  hMC_presel->GetXaxis()->CenterTitle();
600  hMC_presel->GetYaxis()->SetTitle(TString::Format("10^{4} events / %1.2g #times 10^{20} POT", POTmc/1E20).Data());
601  //hData_presel->Draw("E1");
602  hMC_presel->Draw("hist E1");
603  hData_presel->Draw("E1 same");
604  hData_full->Draw("E1 same");
605  hMC_full->Draw("hist E1 same");
606  c6->RedrawAxis();
607  gPad->Update();
608  TText *tt;
609  if(isFHC){
610  auto *tt = new TText(gPad->GetUxmax()*0.6, 1.05*gPad->GetUymax(), "Neutrino beam, MRE");
611  tt->SetTextSize(0.04);
612  tt->Draw();
613  }
614  else{
615  auto *tt = new TText(gPad->GetUxmax()*0.55, 1.05*gPad->GetUymax(), "Antineutrino beam, MRE");
616  tt->SetTextSize(0.04);
617  tt->Draw();
618  }
619 
620  TLegend *leg = new TLegend(0.55, 0.68, 0.85, 0.9);
621  leg->SetTextSize(0.033);
622  leg->AddEntry(hData_presel, "Data, preselected", "p");
623  leg->AddEntry(hMC_presel, "Monte Carlo, preselected", "l");
624  leg->AddEntry(hData_full, "Data, fully selected", "p");
625  leg->AddEntry(hMC_full, "Monte Carlo, fully selected", "l");
626  leg->Draw();
627 
628  Prelim();
629 
630  gPad->RedrawAxis();
631 
632  c6->SaveAs(("blpl/"+weight+mode_tag+"_plot-data-mc-presel-and-fullsel-"+varName+"_"+drawing_options+".pdf"));
633 
634  }
635 
636  }// is plot
637 }
638 
caf::Proxy< caf::SRFuzzyK >::npng
caf::Proxy< size_t > npng
Definition: SRProxy.h:1997
ana::kNue2017MRParentSliceCut
const Cut kNue2017MRParentSliceCut([](const caf::SRProxy *sr){const caf::SRMRCCParentProxy &parent=sr->parent.mrccpar;if(parent.eff< 0.9||parent.pur< 0.9) return false;if(parent.remid< 0.75) return false;if(parent.numuE< 0|| parent.nhit< 20|| parent.contplanes< 5) return false;if(parent.firstplane< 2 || parent.lastplane > 211|| parent.muonstop.Z() > 1275|| parent.muonfwdcell< 5 || parent.muonbkcell< 9 || parent.hadEinmucat > 0.03) return false;return true;})
Definition: NueCuts2017.h:336
plot_validation_datamc.p2
p2
Definition: plot_validation_datamc.py:199
name
const XML_Char * name
Definition: expat.h:151
plot_validation_datamc.fname
fname
Definition: plot_validation_datamc.py:180
ana::kNumVars
const int kNumVars
Definition: vars.h:14
caf::Proxy< caf::SRElastic >::fuzzyk
caf::Proxy< caf::SRFuzzyK > fuzzyk
Definition: SRProxy.h:2018
ana
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
ana::kPPFXFluxCVWgt
const Var kPPFXFluxCVWgt
Definition: PPFXWeights.h:16
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::weight
const Var weight
Definition: PIDCutTuningDefs.h:96
plot_validation_datamc_2018.color
color
Definition: plot_validation_datamc_2018.py:516
Integral
double Integral(const Spectrum &s, const double pot, int cvnregion)
Definition: syst_table_from_preds.C:46
ana::kMREParentNumuE
const Var kMREParentNumuE([](const caf::SRProxy *sr){if(std::isnan(1.*sr->parent.mrccpar.numuE)) return-1000.0f;return float(sr->parent.mrccpar.numuE);})
Definition: MRVars.h:13
make_syst_table_plots.line
line
Definition: make_syst_table_plots.py:268
p
const char * p
Definition: xmltok.h:285
caf::Proxy< caf::StandardRecord >
Proxy for caf::StandardRecord.
Definition: SRProxy.h:2085
std::sqrt
T sqrt(T number)
Definition: d0nt_math.hpp:156
Prelim
void Prelim()
Definition: mre_blessed.C:37
plot_validation_datamc.Clone
def Clone(obj)
Definition: plot_validation_datamc.py:15
ana::HistDef::name
std::string name
Definition: NuePlotLists.h:12
tmp
Float_t tmp
Definition: plot.C:36
ana::SpectrumLoaderBase::SetSpillCut
void SetSpillCut(const SpillCut &cut)
Definition: SpectrumLoaderBase.cxx:64
MakeMePlot
TCanvas * MakeMePlot(std::vector< TH1 * > Histos, std::vector< TString > Names, bool isFHC, string yname="none")
Definition: mre_blessed.C:48
line_data
Style_t line_data
Definition: saveFDMCHists.C:55
ana::kNue2020ND
const Cut kNue2020ND
Definition: NueCuts2020.h:178
ana::Spectrum
Representation of a spectrum in any variable, with associated POT.
Definition: Spectrum.h:40
submit_syst.y
y
Definition: submit_syst.py:249
data
const XML_Char const XML_Char * data
Definition: expat.h:268
caf::Proxy< caf::SRVertexBranch >::elastic
caf::Proxy< caf::SRElastic > elastic
Definition: SRProxy.h:2077
demo5.c2
c2
Definition: demo5.py:33
caf::Proxy< caf::SRFuzzyK >::png
caf::Proxy< std::vector< caf::SRFuzzyKProng > > png
Definition: SRProxy.h:2002
tt
Definition: type_traits.h:56
if
if(dump)
Definition: checkMeasured_Target.C:110
ana::sels
const Cut sels[kNumSels]
Definition: vars.h:44
ana::kNue2020NDPresel
const Cut kNue2020NDPresel
Definition: NueCuts2020.h:169
make_syst_table_plots.h
h
Definition: make_syst_table_plots.py:240
mre_blessed
void mre_blessed(bool plot=true, bool isFHC=true, TString drawing_options="red")
Definition: mre_blessed.C:114
plot_validation_datamc.c
c
Definition: plot_validation_datamc.py:198
d
Float_t d
Definition: plot.C:236
sr
caf::StandardRecord * sr
Definition: FillTrainingTrees.h:18
ana::SpectrumLoader::Go
virtual void Go() override
Load all the registered spectra.
Definition: SpectrumLoader.cxx:66
ana::Spectrum::SaveTo
void SaveTo(TDirectory *dir, const std::string &name) const
Definition: Spectrum.cxx:517
ana::defs
const HistDef defs[kNumVars]
Definition: vars.h:15
demo0.loader
loader
Definition: demo0.py:10
ana::kNumSels
const int kNumSels
Definition: vars.h:43
cmf::Spectrum
std::vector< float > Spectrum
Definition: Constants.h:570
submit_hadd.l
l
Definition: submit_hadd.py:89
prelim
TLatex * prelim
Definition: Xsec_final.C:133
ana::Spectrum::POT
double POT() const
Definition: Spectrum.h:219
isFHC
static bool isFHC
Definition: draw_plots_util.h:46
ana::HistDef::axis
HistAxis axis
Definition: NuePlotLists.h:13
ana::kNoShift
const SystShifts kNoShift
Definition: SystShifts.cxx:21
om::cout
OStream cout
Definition: OStream.cxx:6
ana::kNueEnergy2020
const Var kNueEnergy2020([](const caf::SRProxy *sr){if(!sr->vtx.elastic.IsValid) return-1.0;if(sr->vtx.elastic.fuzzyk.npng< 1) return-1.0;if(kIsRHC(sr)) return kNueEnergyRHC2020(sr);else return kNueEnergyFHC2020(sr);})
Nue energy with 3d prong info. only (old version of vars)
Definition: NueEnergy2020.h:38
xml_schema::string
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
caf::Proxy< caf::SRElastic >::IsValid
caf::Proxy< bool > IsValid
Definition: SRProxy.h:2017
caf::Proxy< caf::StandardRecord >::slc
caf::Proxy< caf::SRSlice > slc
Definition: SRProxy.h:2101
ana::_HistAxis< Var >
ana::SpectrumLoader
Collaborates with Spectrum and OscillatableSpectrum to fill spectra from CAF files.
Definition: SpectrumLoader.h:21
caf::Proxy< caf::SRSlice >::calE
caf::Proxy< float > calE
Definition: SRProxy.h:1248
mc
bool mc
Definition: attenuation_calibration_comparisons.C:33
file
TFile * file
Definition: cellShifts.C:17
demo5.c1
c1
Definition: demo5.py:24
ana::kStandardSpillCuts
const SpillCut kStandardSpillCuts
Apply this unless you&#39;re doing something special.
Definition: SpillCuts.h:49
genie::utils::style::Format
void Format(TGraph *gr, int lcol, int lsty, int lwid, int mcol, int msty, double msiz)
Definition: Style.cxx:154
art::to_string
std::string to_string(ModuleType mt)
Definition: ModuleType.h:32
ana::_Cut< caf::SRProxy >
caf::Proxy< caf::StandardRecord >::vtx
caf::Proxy< caf::SRVertexBranch > vtx
Definition: SRProxy.h:2105
ana::selNames
const std::string selNames[kNumSels]
Definition: vars.h:46
ana::_Var< caf::SRProxy >
w
Float_t w
Definition: plot.C:20
ana::Binning::Simple
static Binning Simple(int n, double lo, double hi, const std::vector< std::string > &labels={})
Definition: Binning.cxx:114
plot_validation_datamc.p1
p1
Definition: plot_validation_datamc.py:199
plot
void plot(std::string label, std::map< std::string, std::map< std::string, Spectrum * >> &plots, bool log)
Definition: plot_resolution.C:184
ana::kXSecCVWgt2020
const Var kXSecCVWgt2020
Definition: XsecTunes.h:106
Generated on Tue Oct 20 2020 05:11:25 for NOvA by   doxygen 1.8.11