preselection_cutflow.C
Go to the documentation of this file.
1 #include "CAFAna/Core/Spectrum.h"
3 
6 
11 
12 #include "TCanvas.h"
13 #include "TLegend.h"
14 
16 
18 
19 using namespace nuebarccinc;
20 
21 
22 std::vector<std::string> cut_labels = {"Data Quality",
23  "NHits",
24  "Front Planes",
25  "Fiducial",
26  "Containment",
27  "MuonID"};
28 
29 struct order_cuts {
30  bool operator() (const std::string & a, const std::string & b)
31  {
32  auto aidx = std::find(cut_labels.begin(), cut_labels.end(), a);
33  auto bidx = std::find(cut_labels.begin(), cut_labels.end(), b);
34  return aidx < bidx;
35  }
36 };
37 
38 struct order_chns {
39  bool operator() (const std::string & a, const std::string & b)
40  {
41  std::vector<std::string> styles;
42  for(auto i = 0u; i < kNumChns; i++) styles.push_back(chns_style[i].latex);
43 
44  auto aidx = std::find(styles.begin(), styles.end(), a);
45  auto bidx = std::find(styles.begin(), styles.end(), b);
46  return aidx < bidx;
47  }
48 };
49 
50 template<class Compare = less<std::string> >
51 void Plot(std::map<std::string, TH1D*, Compare> hists, std::string name, bool show_integral);
52 
53 void Plot(TH1D* hist, std::string label, std::string name, bool show_integral);
54 void preselection_cutflow(std::string input_file_name = "", std::string plot_dump = ".")
55 {
56 
57  std::vector<ana::Cut> cut_flow = {kDQ,
58  kDQ && kNHitsCut,
59  kDQ && kNHitsCut && kFrontPlanes,
60  kDQ && kNHitsCut && kFrontPlanes && kNuebarCCIncFiducial,
61  kDQ && kNHitsCut && kFrontPlanes && kNuebarCCIncFiducial && kNuebarCCIncContainment,
62  kDQ && kNHitsCut && kFrontPlanes && kNuebarCCIncFiducial && kNuebarCCIncContainment && kMubarIDCut};
63 
64 
66  const ana::Var & kCVWeight = ana::kPPFXFluxCVWgt * ana::kXSecCVWgt2020;
67 
68  std::map<std::string, const ana::HistAxis*> axes;
69  axes["nu_e"] = new ana::HistAxis("True Neutrino Energy (GeV)",
70  ana::Binning::Simple(40, 0, 10),
72 // axes["electron_e"] = new ana::HistAxis("True Electron Energy (GeV)",
73 // ana::Binning::Simple(40, 0, 10),
74 // kTrueElectronE);
75 // axes["electron_costheta"] = new ana::HistAxis("True Electron cos#theta",
76 // ana::Binning::Simple(40, -1, 1),
77 // kTrueElectronCosTheta);
78 
79 
80  if(input_file_name == "") {
83 
84 
85  std::map<std::string, std::map<std::string, ana::Spectrum* , order_cuts> > spectra_selected ;
86  std::map<std::string, std::map<std::string, ana::Spectrum* , order_cuts> > spectra_selected_signal;
87  std::map<std::string, std::map<std::string, ana::Spectrum* , order_cuts> > spectra_selected_bkgd ;
88  std::map<std::string, ana::Spectrum* > spectra_signal;
89  std::map<std::string, std::map<std::string, std::map<std::string, ana::Spectrum* > , order_cuts> > spectra_selected_components;
90 
91 
92 
93  for(auto ax_it = axes.begin(); ax_it != axes.end(); ax_it++) {
94  for(auto icut = 0u; icut < cut_flow.size(); icut++) {
95  for(auto ichn = 0; ichn < kNumChns; ichn++) {
96 
97  spectra_selected_components[ax_it->first][cut_labels[icut]][chns[ichn].name] =
98  new ana::Spectrum(loader,
99  *ax_it->second,
100  cut_flow[icut] && chns[ichn].cut,
102  kCVWeight);
103  }
104  spectra_selected_bkgd[ax_it->first][cut_labels[icut]] = new ana::Spectrum(loader,
105  *ax_it->second,
106  !kSignal && cut_flow[icut],
108  kCVWeight);
109  spectra_selected_signal[ax_it->first][cut_labels[icut]] = new ana::Spectrum(loader,
110  *ax_it->second,
111  kSignal && cut_flow[icut],
113  kCVWeight);
114  spectra_selected[ax_it->first][cut_labels[icut]] = new ana::Spectrum(loader,
115  *ax_it->second,
116  cut_flow[icut],
118  kCVWeight);
119  }
120  spectra_signal[ax_it->first] = new ana::Spectrum(loader,
121  *ax_it->second,
122  kSignal,
124  kCVWeight);
125  }
126 
127  loader.Go();
128 
129  TFile * output = new TFile("preselection_cutflow.root", "recreate");
130 
131 // for(auto ichn = 0; ichn < kNumChns; ichn++) {
132 // muonid_components [ichn]->SaveTo(output, ("muonid_" + chns[ichn].name).c_str());
133 // muonid_components_contained[ichn]->SaveTo(output, ("muonid_contained_" + chns[ichn].name).c_str());
134 // }
135  for(auto icut = 0u; icut < cut_flow.size(); icut++) {
136  for(auto ichn = 0; ichn < kNumChns; ichn++) {
137  spectra_selected_components[icut][ichn]->SaveTo(output, (cut_labels[icut] + "_" + chns[ichn].name).c_str());
138  }
139  spectra_selected_bkgd [icut]->SaveTo(output, ("background_and_" + cut_labels[icut]).c_str());
140  spectra_selected_signal[icut]->SaveTo(output, ("signal_and_" + cut_labels[icut]).c_str());
141  spectra_selected [icut]->SaveTo(output, cut_labels[icut] .c_str());
142  }
143  spectrum_signal->SaveTo(output, "signal");
144  } // end fill spectra
145  else {
146  int colors[] = {kRed, kMagenta+2, kBlue-4, kCyan, kGreen, kOrange, kGray+2, kYellow+2};
147 
148  TFile * input_file = TFile::Open(input_file_name.c_str(), "read");
149  std::map<std::string, std::map<std::string, ana::Spectrum* , order_cuts> > spectra_selected ;
150  std::map<std::string, std::map<std::string, ana::Spectrum* , order_cuts> > spectra_selected_signal;
151  std::map<std::string, std::map<std::string, ana::Spectrum* , order_cuts> > spectra_selected_bkgd ;
152  std::map<std::string, ana::Spectrum* > spectra_signal;
153  std::map<std::string, std::map<std::string, std::map<std::string, ana::Spectrum* >, order_cuts > > spectra_selected_components;
154 
155  for(auto ax_it = axes.begin(); ax_it != axes.end(); ax_it++) {
156  for(auto icut = 0u; icut < cut_labels.size(); icut++) {
157  for(auto ichn = 0; ichn < kNumChns; ichn++) {
158  auto name = TString::Format("%s_%s_%s",
159  ax_it->first.c_str(),
160  cut_labels[icut].c_str(),
161  chns[ichn].name.c_str());
162  spectra_selected_components[ax_it->first][cut_labels[icut]][chns_style[ichn].latex] = ana::Spectrum::LoadFrom(input_file->GetDirectory(name)).release();
163 
164  ana::Spectrum * spectrum_signal = ana::Spectrum::LoadFrom(input_file, "signal")).release(;
165 // for(auto ichn = 0; ichn < kNumChns; ichn++) {
166 // muonid_components .push_back(ana::Spectrum::LoadFrom(input_file, ("muonid_" + chns[ichn].name).c_str())).release();
167 // muonid_components_contained.push_back(ana::Spectrum::LoadFrom(input_file, ("muonid_contained_" + chns[ichn].name).c_str())).release();
168 // }
169  for(auto icut = 0u; icut < cut_labels.size(); icut++) {
170  for(auto ichn = 0; ichn < kNumChns; ichn++) {
171  spectra_selected_components[icut][ichn] = ana::Spectrum::LoadFrom(input_file, (cut_labels[icut] + "_" + chns[ichn].name).c_str())).release(;
172  }
173  spectra_selected [icut] = ana::Spectrum::LoadFrom(input_file, cut_labels[icut].c_str())).release(;
174  spectra_selected_signal[icut] = ana::Spectrum::LoadFrom(input_file, ("signal_and_" + cut_labels[icut]).c_str())).release(;
175  spectra_selected_bkgd [icut] = ana::Spectrum::LoadFrom(input_file, ("background_and_" + cut_labels[icut]).c_str())).release(;
176  }
177 
178 
179  std::map<std::string, std::map<std::string, TH1D* , order_cuts> > hists_selected ;
180  std::map<std::string, std::map<std::string, TH1D* , order_cuts> > hists_selected_signal;
181  std::map<std::string, std::map<std::string, TH1D* , order_cuts> > hists_selected_bkgd ;
182  std::map<std::string, TH1D* > hists_signal;
183  std::map<std::string, std::map<std::string, std::map<std::string, TH1D* , order_chns> > > hists_selected_components;
184 
185  for(auto ax_it = axes.begin(); ax_it != axes.end(); ax_it++) {
186  hists_signal[ax_it->first] = spectra_signal.at(ax_it->first)->ToTH1(ana::kAna2019RHCPOT);
187  hists_signal.at(ax_it->first)->SetTitle("Signal");
188  for(auto icut = 0u; icut < cut_labels.size(); icut++) {
189  for(auto ichn = 0u; ichn < spectra_selected_components.at(ax_it->first).at(cut_labels[icut]).size(); ichn++) {
190  hists_selected_components[ax_it->first][cut_labels[icut]][chns_style[ichn].latex] = spectra_selected_components.at(ax_it->first).at(cut_labels[icut]).at(chns_style[ichn].latex)->ToTH1(ana::kAna2019RHCPOT);
191  hists_selected_components.at(ax_it->first).at(cut_labels[icut]).at(chns_style[ichn].latex)->SetLineColor(chns_style[ichn].color);
192  hists_selected_components.at(ax_it->first).at(cut_labels[icut]).at(chns_style[ichn].latex)->SetTitle(("Selected: " + cut_labels[icut]).c_str());
193  }
194 
195  hists_selected [ax_it->first][cut_labels[icut]] = spectra_selected .at(ax_it->first).at(cut_labels[icut])->ToTH1(ana::kAna2019RHCPOT);
196  hists_selected_signal[ax_it->first][cut_labels[icut]] = spectra_selected_signal.at(ax_it->first).at(cut_labels[icut])->ToTH1(ana::kAna2019RHCPOT);
197  hists_selected_bkgd [ax_it->first][cut_labels[icut]] = spectra_selected_bkgd .at(ax_it->first).at(cut_labels[icut])->ToTH1(ana::kAna2019RHCPOT);
198 
199  hists_selected .at(ax_it->first).at(cut_labels[icut])->SetLineColor(colors[icut]);
200  hists_selected_signal.at(ax_it->first).at(cut_labels[icut])->SetLineColor(colors[icut]);
201  hists_selected_bkgd .at(ax_it->first).at(cut_labels[icut])->SetLineColor(colors[icut]);
202  hists_selected .at(ax_it->first).at(cut_labels[icut])->SetTitle("Selected" );
203  hists_selected_signal.at(ax_it->first).at(cut_labels[icut])->SetTitle("Selected Signal");
204  hists_selected_bkgd .at(ax_it->first).at(cut_labels[icut])->SetTitle("Selected Background");
205  }
206  }
207  std::cout << "Hists are created" << std::endl;
208  std::map<std::string, std::map<std::string, TH1D*> > efficiencies;
209  std::map<std::string, std::map<std::string, TH1D*> > purities ;
210  for(auto ax_it = axes.begin(); ax_it != axes.end(); ax_it++) {
211  for(auto icut = 0u; icut < cut_labels.size(); icut++) {
212  efficiencies[ax_it->first][cut_labels[icut]] = (TH1D*) hists_selected_signal.at(ax_it->first).at(cut_labels[icut])->Clone(ana::UniqueName().c_str());
213 
214  efficiencies.at(ax_it->first).at(cut_labels[icut])->Divide(hists_signal.at(ax_it->first));
215  efficiencies.at(ax_it->first).at(cut_labels[icut])->SetTitle("Efficiency");
216  efficiencies.at(ax_it->first).at(cut_labels[icut])->GetYaxis()->SetTitle("");
217 
218  purities [ax_it->first][cut_labels[icut]] = (TH1D*) hists_selected_signal.at(ax_it->first).at(cut_labels[icut])->Clone(ana::UniqueName().c_str());
219  purities .at(ax_it->first).at(cut_labels[icut])->Divide(hists_selected.at(ax_it->first).at(cut_labels[icut]));
220  purities .at(ax_it->first).at(cut_labels[icut])->SetTitle("Purity");
221  purities .at(ax_it->first).at(cut_labels[icut])->GetYaxis()->SetTitle("");
222  }
223  }
224 
225  for(auto ax_it = axes.begin(); ax_it != axes.end(); ax_it++) {
226  for(auto icut = 0u; icut < cut_labels.size(); icut++) {
227  Plot(hists_selected_components.at(ax_it->first).at(cut_labels[icut]), plot_dump + "/" + ax_it->first + "_preselection_cut_flow_bkgd_components_" + cut_labels[icut] + ".pdf", false);
228  }
229 
230  Plot(efficiencies .at(ax_it->first), plot_dump + "/" + ax_it->first + "_preselection_cut_flow_efficiencies.pdf" , false);
231  Plot(purities .at(ax_it->first), plot_dump + "/" + ax_it->first + "_preselection_cut_flow_purities.pdf" , false);
232  Plot(hists_selected .at(ax_it->first), plot_dump + "/" + ax_it->first + "_preselection_cut_flow_selected.pdf" , true);
233  Plot(hists_selected_signal.at(ax_it->first), plot_dump + "/" + ax_it->first + "_preselection_cut_flow_selected_signal.pdf" , true);
234  Plot(hists_selected_bkgd .at(ax_it->first), plot_dump + "/" + ax_it->first + "_preselection_cut_flow_selected_bkgd.pdf" , true);
235  Plot(hists_signal .at(ax_it->first), "Signal" , plot_dump + "/" + ax_it->first + "_preselection_cut_flow_signal.pdf", true);
236 
237  }
238  for(auto ax_it = axes.begin(); ax_it != axes.end(); ax_it++) {
239  std::cout << "Cut Flow (N fiducial nuebar)" << std::endl;
240  for(auto icut = 0u; icut < cut_labels.size(); icut++) {
241  auto nnubar = hists_selected_components.at(ax_it->first).at(cut_labels[icut]).at("#bar{#nu}_{e}")->Integral();
242  auto nnu = hists_selected_components.at(ax_it->first).at(cut_labels[icut]).at("#nu_{e}")->Integral();
243  std::cout << cut_labels[icut] << "\t" << nnubar + nnu << std::endl;
244  }
245  }
246  }
247 }
248 
249 void Plot(TH1D* hist, std::string label, std::string name, bool show_integral)
250 {
251  TCanvas * c = new TCanvas();
252  TLegend * leg = new TLegend();
253  double integral = hist->Integral();
254  if(show_integral)
255  leg->AddEntry(hist, TString::Format("%.2f: %s", integral, label.c_str()), "l");
256  else
257  leg->AddEntry(hist, TString::Format("%s", label.c_str()), "l");
258 
259  hist->Draw("hist");
260  leg->Draw();
261  c->Print(name.c_str());
262 }
263 
264 template<class Compare = less<std::string> >
265 void Plot(std::map<std::string, TH1D*, Compare> hists, std::string name, bool show_integral)
266 {
267  TCanvas * c = new TCanvas();
268  TLegend * leg = new TLegend(0.3, 0.3);
269  double maxval = 0;
270 
271  for(auto h_it = hists.begin(); h_it != hists.end(); h_it++) {
272  maxval = std::max(maxval, h_it->second->GetMaximum());
273  }
274  for(auto h_it = hists.begin(); h_it != hists.end(); h_it++) {
275  double integral = h_it->second->Integral();
276  if(show_integral)
277  leg->AddEntry(h_it->second, TString::Format("%-12d: %s", (int)integral, h_it->first.c_str()), "l");
278  else
279  leg->AddEntry(h_it->second, TString::Format("%s", h_it->first.c_str()), "l");
280  if(h_it == hists.begin()) {
281  h_it->second->SetMaximum(maxval * 1.2);
282  h_it->second->Draw("hist");
283  }
284  else {
285  h_it->second->Draw("hist same");
286  }
287  }
288  // leg->SetTextSize(0.03);
289  leg->Draw();
290  c->Print(name.c_str());
291 }
T max(const caf::Proxy< T > &a, T b)
const XML_Char * name
Definition: expat.h:151
enum BeamMode kOrange
ofstream output
void latex(double x, double y, std::string txt, double ang=0, int align=22, double size=0.042)
const ana::NuTruthCut kNuebarCCST([](const caf::SRNeutrinoProxy *truth){return(truth->iscc && truth->pdg==-12 && truth->pdgorig==-12);})
enum BeamMode kRed
_HistAxis< Var > HistAxis
Definition: HistAxis.h:103
const SelDef cut_flow[knumcutflow]
const Var kPPFXFluxCVWgt
Definition: PPFXWeights.h:16
const ana::Cut kFrontPlanes
const ana::NuTruthCut kTrueDetectorST([](const caf::SRNeutrinoProxy *sr){return(sr->vtx.X()< tvtxmax.X()&&sr->vtx.X() > tvtxmin.X()&&sr->vtx.Y() > tvtxmin.Y()&&sr->vtx.Y()< tvtxmax.Y()&&sr->vtx.Z() > tvtxmin.Z()&&sr->vtx.Z()< tvtxmax.Z());})
const ana::NuTruthCut kNueCCST([](const caf::SRNeutrinoProxy *truth){return(truth->iscc && truth->pdg==12 && truth->pdgorig==12);})
void SetSpillCut(const SpillCut &cut)
TString hists[nhists]
Definition: bdt_com.C:3
std::vector< std::string > cut_labels
Representation of a spectrum in any variable, with associated POT.
Definition: Spectrum.h:40
const SelDef chns[kNumChns]
const char * label
int colors[6]
Definition: tools.h:1
const int kNumChns
static std::unique_ptr< Spectrum > LoadFrom(TDirectory *dir, const std::string &name)
Definition: Spectrum.cxx:535
const std::string PROD5_MC_RHC_NOMINAL
const double a
virtual void Go() override
Load all the registered spectra.
loader
Definition: demo0.py:10
std::vector< float > Spectrum
Definition: Constants.h:759
const ana::Cut kNuebarCCIncFiducial
const SystShifts kNoShift
Definition: SystShifts.cxx:22
OStream cout
Definition: OStream.cxx:6
const double kAna2019RHCPOT
Definition: Exposures.h:224
const ana::Cut kNHitsCut
const ana::Cut kNuebarCCIncContainment
Collaborates with Spectrum and OscillatableSpectrum to fill spectra from CAF files.
const Style chns_style[kNumChns]
const ana::Var kTrueNeutrinoE
const hit & b
Definition: hits.cxx:21
const ana::Cut kDQ([](const caf::SRProxy *sr){if(sr->sel.nuecosrej.hitsperplane >=8) return false;if(!sr->vtx.elastic.IsValid) return false;if(sr->vtx.elastic.fuzzyk.npng==0) return false;if(sr->vtx.elastic.fuzzyk.nshwlid==0) return false;if(sr->vtx.elastic.fuzzyk.png[0].shwlid.nhitx<=5||sr->vtx.elastic.fuzzyk.png[0].shwlid.nhity<=5) return false;if(sr->vtx.elastic.fuzzyk.png[0].shwlid.gap > 100) return false; if(sr->vtx.elastic.fuzzyk.nshwlid >=2){if((sr->vtx.elastic.fuzzyk.png[0].shwlid.dir). Dot(sr->vtx.elastic.fuzzyk.png[1].shwlid.dir)< -0.95) return false;}float xyhitdiff=abs(sr->vtx.elastic.fuzzyk.png[0].shwlid.nhitx-sr->vtx.elastic.fuzzyk.png[0].shwlid.nhity);float xyhitsum=sr->vtx.elastic.fuzzyk.png[0].shwlid.nhitx+sr->vtx.elastic.fuzzyk.png[0].shwlid.nhity;if(xyhitdiff/xyhitsum >=0.4) return false;float nhitall=0.0;for(unsigned int j=0;j< sr->vtx.elastic.fuzzyk.nshwlid;j++){nhitall+=sr->vtx.elastic.fuzzyk.png[j].shwlid.nhit;}if(nhitall/sr->slc.nhit<=0.7) return false;return true;})
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 ana::Cut kMubarIDCut
const Var kMuonID(muonid_classifier::MuonID)
void preselection_cutflow(std::string input_file_name="", std::string plot_dump=".")
void Plot(std::map< std::string, TH1D *, Compare > hists, std::string name, bool show_integral)
enum BeamMode kGreen
string input_file
parseSTDOUT
enum BeamMode kBlue
static Binning Simple(int n, double lo, double hi, const std::vector< std::string > &labels={})
Definition: Binning.cxx:107
Definition: styles.py:1
std::string UniqueName()
Return a different string each time, for creating histograms.
Definition: Utilities.cxx:29
Cut CutFromNuTruthCut(const NuTruthCut &stc)
Definition: Cut.cxx:7
const Var kXSecCVWgt2020
Definition: XsecTunes.h:105
enum BeamMode string