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