mre_comp_split.C
Go to the documentation of this file.
2 #include "CAFAna/Core/Spectrum.h"
4 #include "CAFAna/Core/Binning.h"
6 #include "CAFAna/Cuts/Cuts.h"
7 #include "CAFAna/Vars/Vars.h"
9 
10 #include "CAFAna/Cuts/NueCutsSecondAna.h"
13 
19 
20 #include "TCanvas.h"
21 #include "TFile.h"
22 #include "TH1.h"
23 #include "TROOT.h"
24 #include "TRint.h"
25 #include "TColor.h"
26 #include "TFile.h"
27 #include "TLine.h"
28 #include "TText.h"
29 #include "TLegend.h"
30 #include "TLatex.h"
31 #include "TPad.h"
32 
33 #include <fstream>
34 #include <iomanip>
35 #include <iostream>
36 
37 #include "Utilities/func/MathUtil.h"
38 
39 void Prelim()
40 {
41  TLatex* prelim = new TLatex(.93, .9, "NOvA Preliminary");
42  prelim->SetTextColor(kBlue);
43  prelim->SetNDC();
44  prelim->SetTextSize(2/30.);
45  prelim->SetTextAngle(270);
46  prelim->SetTextAlign(12);
47  prelim->Draw();
48 }
49 
50 using namespace ana;
51 
52 void mre_comp_split(bool plot = true, bool isFHC = true){
53 
54  const std::string mode_tag = isFHC?"fhc":"rhc";
55  std::string file_name = "mre_pred_no_extrap_nd_cuts_"+mode_tag+"_mc_vs_data.root";
56 
57  const int kNumVars = 3;
58 
59  struct HistDef
60  {
62  HistAxis axis;
63  };
64 
65  const Var kSlcE([](const caf::SRProxy *sr){
66  return sr->slc.calE;
67  });
68 
69  const HistDef defs[kNumVars] = {
70  {"energy", {"#nu_{e} Energy(GeV)", Binning::Simple(120,0,12), kNueEnergy2018}},
71  {"caloE", {"Calorimetric energy", Binning::Simple(120,0,12), kSlcE}},
72  {"cvne", {"CVNe", Binning::Simple(120,-0.1,1.1), kCVNSSe}},
73  };
74 
75  const int kNumSels=2;
76 
77  const Cut sels[kNumSels] = {
79  kNue2018NDCVNSsb && kNue2017MRParentSliceCut
80  };
81 
82  const std::string selNames[kNumSels] = {
83  "nd-presel-mre",
84  "nd-full-mre"
85  };
86 
87  if (!plot){
88 
89  std::string fname = "";
90  std::string fname_data = "";
91 
92  // FHC MC
93  if(isFHC){
94  fname = "defname: prod_mrecaf_R17-11-14-prod4mre.b_nd_genie_fhc_full_v1 with stride 1";
95  fname_data = "defname: prod_mrecaf_R17-11-14-prod4mre.b_nd_numi_fhc_full_v1_goodruns";
96  }
97  // RHC MC
98  if(!isFHC) {
99  fname = "defname: prod_mrecaf_R17-11-14-prod4mre.b_nd_genie_rhc_full_v1 with stride 1";
100  fname_data = "defname: prod_mrecaf_R17-11-14-prod4mre.b_nd_numi_rhc_full_v1_goodruns";
101  }
102 
103 
104  SpectrumLoader loader(fname);
105  SpectrumLoader loader_data(fname_data);
106 
108  loader_data.SetSpillCut(kStandardSpillCuts);
109 
110 
111  Spectrum* spec_data[kNumSels][kNumVars];
112  Spectrum* spec[kNumSels][kNumVars];
113  Spectrum* spec_dis[kNumSels][kNumVars];
114  Spectrum* spec_res[kNumSels][kNumVars];
115  Spectrum* spec_qe[kNumSels][kNumVars];
116  Spectrum* spec_mec[kNumSels][kNumVars];
117  Spectrum* spec_coh[kNumSels][kNumVars];
118 
119  for(int selIdx = 0; selIdx < kNumSels; ++selIdx){
120  for (int varIdx=0; varIdx<kNumVars;varIdx++){
121  const HistAxis& axis = defs[varIdx].axis;
122  spec[selIdx][varIdx] = new Spectrum(loader, axis, sels[selIdx], kNoShift, kPPFXFluxCVWgt*kXSecCVWgt2018); // All int
123  spec_dis[selIdx][varIdx] = new Spectrum(loader, axis, sels[selIdx] && kIsDIS, kNoShift, kPPFXFluxCVWgt*kXSecCVWgt2018); // DIS
124  spec_res[selIdx][varIdx] = new Spectrum(loader, axis, sels[selIdx] && kIsRes, kNoShift, kPPFXFluxCVWgt*kXSecCVWgt2018); // RES
125  spec_qe[selIdx][varIdx] = new Spectrum(loader, axis, sels[selIdx] && kIsQE, kNoShift, kPPFXFluxCVWgt*kXSecCVWgt2018); // QE
126  spec_mec[selIdx][varIdx] = new Spectrum(loader, axis, sels[selIdx] && kIsDytmanMEC, kNoShift, kPPFXFluxCVWgt*kXSecCVWgt2018); // MEC
127  spec_coh[selIdx][varIdx] = new Spectrum(loader, axis, sels[selIdx] && kIsCoh, kNoShift, kPPFXFluxCVWgt*kXSecCVWgt2018); // COH
128  spec_data[selIdx][varIdx] = new Spectrum(loader_data, axis, sels[selIdx]);
129  }
130  }
131 
132  loader.Go();
133  loader_data.Go();
134 
135  TFile* file = new TFile(file_name.c_str(),"recreate");
136 
137  for(int selIdx = 0; selIdx < kNumSels; ++selIdx){
138  TDirectory* d = file->mkdir(selNames[selIdx].c_str());
139  for (int varIdx=0; varIdx<kNumVars;varIdx++){
140  const char* name = defs[varIdx].name.c_str();
141  spec[selIdx][varIdx]->SaveTo(d, TString::Format("%s", name));
142  spec_dis[selIdx][varIdx]->SaveTo(d, TString::Format("dis_%s", name));
143  spec_res[selIdx][varIdx]->SaveTo(d, TString::Format("res_%s", name));
144  spec_qe[selIdx][varIdx]->SaveTo(d, TString::Format("qe_%s", name));
145  spec_mec[selIdx][varIdx]->SaveTo(d, TString::Format("mec_%s", name));
146  spec_coh[selIdx][varIdx]->SaveTo(d, TString::Format("coh_%s", name));
147  spec_data[selIdx][varIdx] ->SaveTo(d, TString::Format("data_%s", name));
148  }
149  }
150  file->Close();
151  } // end of !plot
152 
153  if (plot){
154  Spectrum* spec_data[kNumSels][kNumVars];
155  Spectrum* spec[kNumSels][kNumVars];
156  Spectrum* spec_dis[kNumSels][kNumVars];
157  Spectrum* spec_res[kNumSels][kNumVars];
158  Spectrum* spec_qe[kNumSels][kNumVars];
159  Spectrum* spec_mec[kNumSels][kNumVars];
160  Spectrum* spec_coh[kNumSels][kNumVars];
161 
162  TH1* hMC[kNumSels][kNumVars];
163  TH1* hMC_dis[kNumSels][kNumVars];
164  TH1* hMC_res[kNumSels][kNumVars];
165  TH1* hMC_qe[kNumSels][kNumVars];
166  TH1* hMC_mec[kNumSels][kNumVars];
167  TH1* hMC_coh[kNumSels][kNumVars];
168  TH1* hData[kNumSels][kNumVars];
169 
170  double POT = -1;
171 
172  for(int selIdx = 0; selIdx < kNumSels; ++selIdx){
173  const char* selName = selNames[selIdx].c_str();
174  for(int varIdx = 0; varIdx < kNumVars; ++varIdx){
175  const char* varName = defs[varIdx].name.c_str();
176 
177  spec[selIdx][varIdx] = LoadFromFile<Spectrum>(file_name, TString::Format("%s/%s", selName, varName).Data()).release();
178  spec_dis[selIdx][varIdx] = LoadFromFile<Spectrum>(file_name, TString::Format("%s/dis_%s", selName, varName).Data()).release();
179  spec_res[selIdx][varIdx] = LoadFromFile<Spectrum>(file_name, TString::Format("%s/res_%s", selName, varName).Data()).release();
180  spec_qe[selIdx][varIdx] = LoadFromFile<Spectrum>(file_name, TString::Format("%s/qe_%s", selName, varName).Data()).release();
181  spec_mec[selIdx][varIdx] = LoadFromFile<Spectrum>(file_name, TString::Format("%s/mec_%s", selName, varName).Data()).release();
182  spec_coh[selIdx][varIdx] = LoadFromFile<Spectrum>(file_name, TString::Format("%s/coh_%s", selName, varName).Data()).release();
183  spec_data[selIdx][varIdx] = LoadFromFile<Spectrum>(file_name, TString::Format("%s/data_%s", selName, varName).Data()).release();
184 
185  POT = spec_data[selIdx][varIdx]->POT();
186  cout<<"data pot is "<<POT<<endl;
187 
188  hMC[selIdx][varIdx]=spec[selIdx][varIdx]->ToTH1(POT, kPink, kSolid, ana::EExposureType::kPOT);
189  hMC_dis[selIdx][varIdx]=spec_dis[selIdx][varIdx]->ToTH1(POT, kGray+1, kSolid, ana::EExposureType::kPOT);
190  hMC_res[selIdx][varIdx]=spec_res[selIdx][varIdx]->ToTH1(POT, kGreen+2, kSolid, ana::EExposureType::kPOT);
191  hMC_qe[selIdx][varIdx]=spec_qe[selIdx][varIdx]->ToTH1(POT, kAzure+7, kSolid, ana::EExposureType::kPOT);
192  hMC_mec[selIdx][varIdx]=spec_mec[selIdx][varIdx]->ToTH1(POT, kOrange-3, kSolid, ana::EExposureType::kPOT);
193  hMC_coh[selIdx][varIdx]=spec_coh[selIdx][varIdx]->ToTH1(POT, kViolet-1, kSolid, ana::EExposureType::kPOT);
194  hData[selIdx][varIdx]=spec_data[selIdx][varIdx]->ToTH1(POT, kBlack, kSolid, ana::EExposureType::kPOT);
195 
196  } // end for varIdx
197  } // end for selIdx
198 
199 
200  // now let's make one plot for blessing - CVNe for presel with different int. types
201  int selIdx = 0;
202  int varIdx = 2;
203 
204  cout<<"integrals "<<hMC[selIdx][varIdx]->Integral()<<" "<<hMC_dis[selIdx][varIdx]->Integral()<<" "<<hMC_res[selIdx][varIdx]->Integral()<<" "<<hMC_qe[selIdx][varIdx]->Integral()<<" "<<hMC_mec[selIdx][varIdx]->Integral()<<" "<<hMC_coh[selIdx][varIdx]->Integral()<<" "<<hData[selIdx][varIdx]->Integral()<<endl;
205 
206  TCanvas *c = new TCanvas("c", "canvas", 1000, 1000);
207  TPad* p1 = new TPad("", "", 0, 0, 1, 1);
208  TPad* p2 = new TPad("", "", 0, 0, 1, 1);
209  p1->SetBottomMargin(.3);
210  p2->SetTopMargin(.7);
211 
212  for(auto p:{p1,p2}){
213  p->SetLeftMargin(0.105);
214  p->SetFillStyle(0);
215  p->Draw();
216  }
217  p1->cd();
218 
219  for (auto &h:{hMC[selIdx][varIdx], hMC_dis[selIdx][varIdx], hMC_res[selIdx][varIdx], hMC_qe[selIdx][varIdx], hMC_mec[selIdx][varIdx], hMC_coh[selIdx][varIdx], hData[selIdx][varIdx]}){
220  //h->SetLineWidth(2);
221  h->Rebin(2);
222  h->Scale(1/10000.0);
223  }
224  hMC[selIdx][varIdx]->SetMaximum(hMC[selIdx][varIdx]->GetMaximum()*1.2);
225  hMC[selIdx][varIdx]->SetMinimum(0);
226  hMC[selIdx][varIdx]->GetYaxis()->CenterTitle();
227  hMC[selIdx][varIdx]->GetXaxis()->CenterTitle();
228  hMC[selIdx][varIdx]->GetYaxis()->SetTitleSize(0.045);
229  hMC[selIdx][varIdx]->GetYaxis()->SetTitleOffset(1.17);
230  hMC[selIdx][varIdx]->GetYaxis()->SetTitle(TString::Format("10^{4} events / %1.2g #times 10^{20} POT", POT/1E20).Data());
231 
232  hMC[selIdx][varIdx]->Draw("hist");
233  hMC_dis[selIdx][varIdx]->Draw("hist same");
234  hMC_res[selIdx][varIdx]->Draw("hist same");
235  hMC_qe[selIdx][varIdx]->Draw("hist same");
236  hMC_mec[selIdx][varIdx]->Draw("hist same");
237  hMC_coh[selIdx][varIdx]->Draw("hist same");
238  hData[selIdx][varIdx]->SetMarkerStyle(20);
239  hData[selIdx][varIdx]->Draw("E1 same");
240 
241  c->RedrawAxis();
242  gPad->Update();
243 
244  TText *tt;
245  TText *tt2;
246  TLine *l2;
247  if(isFHC){
248  auto *tt = new TText(gPad->GetUxmax()*0.53, 1.05*gPad->GetUymax(), "Neutrino beam, MRE");
249  tt->SetTextSize(0.04);
250  tt->Draw();
251  auto *l2 = new TLine(0.84,0,0.84,gPad->GetUymax());
252  l2->SetLineColor(kBlack);
253  l2->Draw();
254  auto *tt2 = new TText(0.82, 0.6*gPad->GetUymax(), "Cut position");
255  tt2->SetTextAngle(90.0);
256  tt2->SetTextSize(0.035);
257  tt2->Draw();
258  }
259  else{
260  auto *tt = new TText(gPad->GetUxmax()*0.46, 1.05*gPad->GetUymax(), "Antineutrino beam, MRE");
261  tt->SetTextSize(0.04);
262  tt->Draw();
263  auto *l2 = new TLine(0.89,0,0.89,gPad->GetUymax());
264  l2->SetLineColor(kBlack);
265  l2->Draw();
266  auto *tt2 = new TText(0.87, 0.6*gPad->GetUymax(), "Cut position");
267  tt2->SetTextAngle(90.0);
268  tt2->SetTextSize(0.035);
269  tt2->Draw();
270  }
271 
272  TLegend *leg = new TLegend(0.2,0.6,0.5,0.85);
273  leg->SetTextSize(0.03);
274  leg->AddEntry(hData[selIdx][varIdx], "Data", "p");
275  leg->AddEntry(hMC[selIdx][varIdx], "Monte Carlo, Total", "l");
276  leg->AddEntry(hMC_qe[selIdx][varIdx], "Monte Carlo, QE", "l");
277  leg->AddEntry(hMC_res[selIdx][varIdx], "Monte Carlo, RES", "l");
278  leg->AddEntry(hMC_mec[selIdx][varIdx], "Monte Carlo, MEC", "l");
279  leg->AddEntry(hMC_dis[selIdx][varIdx], "Monte Carlo, DIS", "l");
280  leg->AddEntry(hMC_coh[selIdx][varIdx], "Monte Carlo, COH", "l");
281  leg->Draw();
282 
283  Prelim();
284 
285  p2->cd();
286  auto temp = (TH1*)hData[selIdx][varIdx]->Clone();
287  temp->Divide(hMC[selIdx][varIdx]);
288  temp->SetLineColor(kPink);
289  temp->SetMarkerStyle(1);
290  temp->GetXaxis()->SetTitle("#nu_{e} CVN classifier");
291  temp->Draw();
292 
293  TLine *l = new TLine(-0.1,1,1.1,1);
294  l->SetLineColor(kBlack);
295  l->Draw();
296  temp->GetYaxis()->SetRangeUser(0.75, 1.25);
297  temp->GetYaxis()->CenterTitle();
298  temp->GetXaxis()->CenterTitle();
299  temp->GetYaxis()->SetTitle("Data / MC");
300  temp->GetYaxis()->SetTitleSize(0.04);
301  temp->GetYaxis()->SetTitleOffset(1.2);
302  p2->RedrawAxis();
303 
304  c->SaveAs(("blpl/"+mode_tag+"_plot-data-mc_cvne_presel.pdf").c_str());
305  }
306 }
const Cut kIsQE
Definition: TruthCuts.cxx:104
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
void Prelim()
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
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
const Cut kIsRes
Definition: TruthCuts.cxx:111
const char * p
Definition: xmltok.h:285
Proxy for caf::StandardRecord.
Definition: SRProxy.h:2085
const Cut kNue2018NDCVNSsb
Definition: NueCuts2018.h:153
std::string name
Definition: NuePlotLists.h:12
const Cut kIsDIS
Definition: TruthCuts.cxx:118
void SetSpillCut(const SpillCut &cut)
const Cut kNue2018NDPresel
Definition: NueCuts2018.h:145
Representation of a spectrum in any variable, with associated POT.
Definition: Spectrum.h:40
const Cut kIsDytmanMEC
Definition: TruthCuts.cxx:187
void mre_comp_split(bool plot=true, bool isFHC=true)
Definition: type_traits.h:56
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
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
std::vector< double > POT
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
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
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
const std::string selNames[kNumSels]
Definition: vars.h:46
static Binning Simple(int n, double lo, double hi, const std::vector< std::string > &labels={})
Definition: Binning.cxx:114
const Cut kIsCoh
Definition: TruthCuts.cxx:133
void plot(std::string label, std::map< std::string, std::map< std::string, Spectrum * >> &plots, bool log)