numu_cut_flow.C
Go to the documentation of this file.
1 // Fills spectra for systematic error band
2 #ifdef __CINT__
3 void numu_cut_flow()
4 {
5  std::cout << "Sorry, you must run in compiled mode" << std::endl;
6 }
7 #else
8 
9 
10 
11 #include "CAFAna/Cuts/Cuts.h"
12 #include "CAFAna/Cuts/NumuCuts.h"
13 #include "CAFAna/Core/Binning.h"
14 #include "CAFAna/Core/Spectrum.h"
16 #include "CAFAna/Systs/Systs.h"
17 #include "CAFAna/Systs/NumuSysts.h"
18 #include "CAFAna/Cuts/SpillCuts.h"
19 #include "CAFAna/Analysis/Plots.h"
20 #include "CAFAna/Vars/HistAxes.h"
22 #include "CAFAna/Core/Spectrum.h"
23 #include "OscLib/IOscCalc.h"
24 #include "OscLib/OscCalcDumb.h"
25 
26 
27 #include "TCanvas.h"
28 #include "TFile.h"
29 #include "TGraph.h"
30 #include "TH1.h"
31 #include "TH2.h"
32 #include "TMath.h"
33 #include "TGaxis.h"
34 #include "TMultiGraph.h"
35 #include "TLegend.h"
36 #include "TLatex.h"
37 #include "TStyle.h"
38 #include "THStack.h"
39 #include "TPaveText.h"
40 #include <cmath>
41 #include <iostream>
42 #include <vector>
43 #include <list>
44 #include <sstream>
45 #include <string>
46 #include <fstream>
47 #include <iomanip>
48 #include <stdexcept>
49 
50 
51 using namespace ana;
52 
53 class CanMan {
54 
55 public:
56  CanMan():
57  fCans()
58  {}
59 
61  {
62  fCans.push_back(new TCanvas(name.c_str(), title.c_str()));
63  return fCans.back();
64  }
65 
66  TCanvas* NewCan(std::string nameTitle)
67  {
68  return NewCan(nameTitle, nameTitle);
69  }
70 
71  TCanvas* SaveCans(std::string prefix="", std::string suffix=".pdf")
72  {
73  for(const auto& can:fCans)
74  {
75  std::string name = prefix + std::string(can->GetName()) + suffix;
76  can->SaveAs(name.c_str());
77  }
78  }
79 
80  std::vector<TCanvas*> fCans;
81 
82 };
83 
84 
86 {
87  if (opt.find("same") == std::string::npos &&
88  opt.find("SAME") == std::string::npos)
89  opt += " SAME ";
90 
91 }
92 
93 
94 class Sample
95 {
96 public:
97 
99  :
100  fShortName(shortName),
101  fLongName(longName),
102  fCut(cut)
103  {};
104 
108 
109 };
110 
111 class CutLevel
112 {
113 public:
114 
116  SpectrumLoader& loader, const HistAxis& axis, int color)
117  :
118  fSpec(loader, axis, cut),
119  fCut(cut),
120  fShortName(shortName),
121  fLongName(longName),
122  fColor(color)
123  {}
124 
125  void Draw(float pot, std::string opt,
126  osc::IOscCalc* calc, int from, int to, TLegend* leg,
127  int style=1)
128  {
129  TH1* h = fSpec.Oscillated(calc, from, to).ToTH1(pot);
130  h->SetLineColor(fColor);
131  h->SetLineStyle(style);
132  h->Draw(opt.c_str());
133  leg->AddEntry(h, fLongName.c_str(), "l");
134 
135  };
136 
138  osc::IOscCalc* calc, int from, int to, TH1* denominator,
139  TLegend* leg)
140  {
141  TH1* h = fSpec.Oscillated(calc, from, to).ToTH1(pot);
142  h->GetYaxis()->SetTitle("Efficiency");
143  h->Divide(denominator);
144  h->SetLineColor(fColor);
145  h->SetLineStyle(fStyle);
146  h->Draw(opt.c_str());
147  leg->AddEntry(h, fLongName.c_str(), "l");
148  };
149 
150  TH1* ToTH1(float pot, osc::IOscCalc* calc, int from, int to)
151  {
152  return fSpec.Oscillated(calc, from, to).ToTH1(pot);
153  }
154 
155 
156 
161  int fColor;
162  int fStyle;
163 };
164 
165 
166 class CutFlow
167 {
168 public:
170  std::vector<Sample> samples):
171  fName(name),
172  fLoader(loader),
173  fAxis(axis)
174  {
175 
176  fSamples.reserve(samples.size());
177  for(const auto& sample:samples)
178  {
179  fSamples.emplace_back(sample.fShortName, sample.fLongName,
180  sample.fCut,
181  fLoader, fAxis, kBlack);
182  }
183  fLevels.resize(fSamples.size());
184  };
185 
187  {
188  for(size_t iSample = 0; iSample < fSamples.size(); ++iSample)
189  {
190  fLevels[iSample].emplace_back(shortName, longName,
191  fSamples[iSample].fCut && cut,
192  fLoader, fAxis, color);
193  }
194 
195  };
196 
197 
198  void Reserve(size_t n)
199  {
200 
201  for(auto& samp:fLevels)
202  samp.reserve(n);
203  };
204 
205 
206  int findSample(std::string sampleName)
207  {
208  int i = -1;
209  for(size_t iSample = 0; iSample < fSamples.size(); ++iSample){
210  if(!sampleName.compare(fSamples[iSample].fShortName))
211  {
212  i = iSample;
213  }
214  }
215  if(i == -1) throw std::invalid_argument("Sample not found");
216  return i;
217  }
218 
219  void Draw(float pot, std::string sample, std::string opt,
220  osc::IOscCalc* calc, int from, int to, int style=1)
221  {
222  int i = findSample(sample);
223 
224  std::string name = fName + "_" + fSamples[i].fShortName;
225  std::string currOpt(opt);
226  TLegend* leg = new TLegend(0.55, 0.55, 0.85, 0.85);
227  fSamples[i].Draw(pot, currOpt, calc, from, to, leg, style);
228  addSame(currOpt);
229  for(size_t iLevel = 0; iLevel < fLevels[i].size(); ++iLevel){
230  fLevels[i][iLevel].Draw(pot, currOpt, calc, from, to, leg, style);
231  }
232  leg->Draw();
233 
234  };
235 
236 
238  osc::IOscCalc* calc, int from, int to)
239  {
240 
241  int i = findSample(sample);
242  std::string name = fName + "_" + fSamples[i].fShortName;
243  std::string currOpt(opt);
244  TLegend* leg = AutoPlaceLegend(0.3, 0.23, 0.33);
245 
246  TH1* denominator = fSamples[i].ToTH1(pot, calc, from, to);
247 
248  for(size_t iLevel = 0; iLevel < fLevels[i].size(); ++iLevel){
249  fLevels[i][iLevel].DrawEfficiency(pot, currOpt,
250  calc, from, to, denominator, leg);
251  addSame(currOpt);
252  }
253  leg->Draw();
254 
255  };
256 
257 
258 
262  std::vector<CutLevel> fSamples;
263 
264  std::vector<std::vector<CutLevel> > fLevels;
265 
266 
267 };
268 
269 
270 
272 {
273 
274 
275  Cut truthContain({"sel.contain.missE", "mc.nu.E", "mc.nnu"},
276  [](const caf::StandardRecord* sr)
277  {
278  if(sr->mc.nnu == 0) return false;
279  return sr->sel.contain.missE / sr->mc.nu[0].E < 0.01;
280  });
281 
282  Cut cvn({"sel.cvnnumu.output", "sel.cvnnumu.noutput"},
283  [](const caf::StandardRecord* sr)
284  {
285  if(!sr->sel.cvnnumu.noutput) return false;
286  return sr->sel.cvnnumu.output[0] > 0.35;
287  });
288 
289  Cut cvn2({"sel.cvnnumu.output", "sel.cvnnumu.noutput"},
290  [](const caf::StandardRecord* sr)
291  {
292  if(!sr->sel.cvnnumu.noutput) return false;
293  return sr->sel.cvnnumu.output[0] > 0.5;
294  });
295 
296 
297  Var cvnOut({"sel.cvnnumu.output", "sel.cvnnumu.noutput"},
298  [](const caf::StandardRecord* sr)
299  {
300  if(!sr->sel.cvnnumu.noutput) return -5.f;
301  return sr->sel.cvnnumu.output[0];
302  });
303 
304  HistAxis cvnAxis("; CVN #nu_{#mu};Events", Binning::Simple(40, 0,1), cvnOut);
305 
306  SpectrumLoader nonswap(
307  "/cloud/login/rocco/development/results/2015-10-20_caf_nonswap_numu/fardet_genie_fhc_nonswap_ideal_3000_r01000001_s*");
309  CutFlow remid("remid", nonswap,
310  kNumuCCAxis,
311  {Sample("sig", "Truth Containment",
312  kIsNumuCC && truthContain),
313  Sample("bkg", "Truth Containment",
314  !kIsNumuCC && truthContain)});
315 
316 
317  remid.Reserve(10);
318 
319 
320  remid.AddLevel("quality", "Quality", truthContain && kNumuQuality,
321  kGreen+2);
322  remid.AddLevel("containment", "Containment",
323  truthContain && kNumuQuality && kNumuContainFD, kBlue);
324  remid.AddLevel("NCRej", "CVN Loose",
325  truthContain && kNumuQuality && kNumuContainFD && cvn,
326  kMagenta-3);
327 
328  remid.AddLevel("NCRej", "CVN Tight",
329  truthContain && kNumuQuality && kNumuContainFD && cvn2,
330  kCyan-3);
331  remid.AddLevel("NCRej", "ReMId",
332  truthContain && kNumuQuality && kNumuContainFD && kNumuNCRej,
333  kRed);
334 
335 
336  OscillatableSpectrum cvnSig(nonswap, cvnAxis, kIsNumuCC && truthContain);
337  OscillatableSpectrum cvnBkg(nonswap, cvnAxis, !kIsNumuCC && truthContain);
338 
339 
340 
341  nonswap.Go();
342 
343  CanMan canMan;
344 
346  canMan.NewCan("sig");
347  remid.Draw(3e20, "sig","H", calc, 14, 14);
348 
349  canMan.NewCan("sig_eff");
350  remid.DrawEfficiency(3e20, "sig","H", calc, 14, 14);
351 
352  canMan.NewCan("bkg");
353  remid.Draw(3e20, "bkg","H", calc, 14, 14);
354 
355  canMan.NewCan("bkg_eff");
356  remid.DrawEfficiency(3e20, "bkg","H", calc, 14, 14);
357 
358  canMan.NewCan("cvn_dist");
359  TH1* hSig = cvnSig.Oscillated(calc, 14, 14).ToTH1(3e20);
360  TH1* hBkg = cvnBkg.Oscillated(calc, 14, 14).ToTH1(3e20);
361  hSig->SetLineColor(kRed);
362  hBkg->SetLineColor(kBlue);
363  hSig->Draw("HIST ");
364  hBkg->Draw("HIST SAME");
365 
366  TLegend* leg = AutoPlaceLegend(0.35, 0.15);
367  leg->AddEntry(hSig, "Contained #nu_{#mu} Signal", "l");
368  leg->AddEntry(hBkg, "Contained Background", "l");
369  leg->Draw();
370 
371 
372 
373  delete calc;
374 
375  osc::IOscCalc* calcOsc = new osc::OscCalcDumb();
376  canMan.NewCan("sig_osc");
377  remid.Draw(3e20, "sig","H", calcOsc, 14, 14);
378 
379  canMan.NewCan("sig_eff_osc");
380  remid.DrawEfficiency(3e20, "sig","H", calcOsc, 14, 14);
381 
382  canMan.NewCan("bkg_osc");
383  remid.Draw(3e20, "bkg","H", calcOsc, 14, 14);
384 
385  canMan.NewCan("bkg_eff_osc");
386  remid.DrawEfficiency(3e20, "bkg","H", calcOsc, 14, 14);
387 
388  canMan.NewCan("cvn_dist_osc");
389  TH1* hSigOsc = cvnSig.Oscillated(calcOsc, 14, 14).ToTH1(3e20);
390  TH1* hBkgOsc = cvnBkg.Oscillated(calcOsc, 14, 14).ToTH1(3e20);
391  hSigOsc->SetLineColor(kRed);
392  hBkgOsc->SetLineColor(kBlue);
393  hBkgOsc->Draw("HIST");
394  hSigOsc->Draw("HIST SAME ");
395 
396  TLegend* legOsc = AutoPlaceLegend(0.35, 0.15);
397  legOsc->AddEntry(hSigOsc, "Contained #nu_{#mu} Signal", "l");
398  legOsc->AddEntry(hBkgOsc, "Contained Background", "l");
399  legOsc->Draw();
400 
401  canMan.SaveCans("cvn_plots/", ".pdf");
402 
403 }
404 
405 
406 #endif
void DrawEfficiency(float pot, std::string opt, osc::IOscCalc *calc, int from, int to, TH1 *denominator, TLegend *leg)
std::vector< TCanvas * > fCans
Definition: numu_cut_flow.C:80
std::string fShortName
const XML_Char * name
Definition: expat.h:151
const Cut kIsNumuCC(CCFlavSel(14, 14))
Select CC .
Definition: canMan.py:1
void Reserve(size_t n)
TH1 * ToTH1(float pot, osc::IOscCalc *calc, int from, int to)
enum BeamMode kRed
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
float missE
sum of energy of particles that leave detector. example cut: missE/trueE > 0.01 is truly uncontained ...
Definition: SRContain.h:21
Sample(std::string shortName, std::string longName, Cut cut)
Definition: numu_cut_flow.C:98
std::vector< std::vector< CutLevel > > fLevels
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
OscillatableSpectrum fSpec
TCanvas * SaveCans(std::string prefix="", std::string suffix=".pdf")
Definition: numu_cut_flow.C:71
std::string fLongName
General interface to oscillation calculators.
Definition: StanTypedefs.h:22
std::string fName
std::string fShortName
const Cut kNumuContainFD([](const caf::SRProxy *sr){ std::pair< int, int > planes=calcFirstLastLivePlane(sr->slc.firstplane, std::bitset< 14 >(sr->hdr.dibmask));int planestofront=sr->slc.firstplane-planes.first;int planestoback=planes.second-sr->slc.lastplane;return( sr->slc.ncellsfromedge > 1 &&planestofront > 1 &&planestoback > 1 &&sr->sel.contain.kalfwdcell > 10 &&sr->sel.contain.kalbakcell > 10 &&sr->sel.contain.cosfwdcell > 0 &&sr->sel.contain.cosbakcell > 0);})
Definition: NumuCuts.h:20
void SetSpillCut(const SpillCut &cut)
HistAxis fAxis
Defines an enumeration for prong classification.
osc::OscCalcDumb calc
int findSample(std::string sampleName)
_NoOscillations< double > NoOscillations
Definition: IOscCalc.h:67
void AddLevel(std::string shortName, std::string longName, Cut cut, int color)
SpectrumLoader & fLoader
const HistAxis kNumuCCAxis("Reconstructed Neutrino Energy (GeV)", kNumuEnergyBinning, kCCE)
Definition: HistAxes.h:8
TCanvas * NewCan(std::string name, std::string title)
Definition: numu_cut_flow.C:60
#define pot
Simple oscillation probability calculator that has no solar term or mass hierarchy or delta so it&#39;s s...
Definition: OscCalcDumb.h:16
caf::StandardRecord * sr
virtual void Go() override
Load all the registered spectra.
std::vector< CutLevel > fSamples
void invalid_argument(const char *function, const char *name, const T &y, const char *msg1, const char *msg2)
const Cut kNumuNCRej([](const caf::SRProxy *sr){return(sr->sel.remid.pid >0.75);})
Definition: NumuCuts.h:24
loader
Definition: demo0.py:10
short nnu
Number of neutrinos in nu vector (0 or 1)
Definition: SRTruthBranch.h:37
const XML_Char * prefix
Definition: expat.h:380
OStream cout
Definition: OStream.cxx:6
The StandardRecord is the primary top-level object in the Common Analysis File trees.
const Cut cut
Definition: exporter_fd.C:30
Spectrum Oscillated(osc::IOscCalc *calc, int from, int to) const
Collaborates with Spectrum and OscillatableSpectrum to fill spectra from CAF files.
SRIDBranch sel
Selector (PID) branch.
void numu_cut_flow()
string shortName
THUMBNAIL BLOCK: We need to make a thumbnail for each.
SRTruthBranch mc
Truth branch for MC: energy, flavor, etc.
void DrawEfficiency(float pot, std::string sample, std::string opt, osc::IOscCalc *calc, int from, int to)
TLegend * AutoPlaceLegend(double dx, double dy, double yPin)
Create a legend, maximizing distance from all histograms.
Definition: Plots.cxx:717
const SpillCut kStandardSpillCuts
Apply this unless you&#39;re doing something special.
Definition: SpillCuts.h:49
A PID for muons.
Definition: FillPIDs.h:10
void addSame(std::string &opt)
Definition: numu_cut_flow.C:85
const Cut kNumuQuality
Definition: NumuCuts.h:18
CutLevel(std::string shortName, std::string longName, Cut cut, SpectrumLoader &loader, const HistAxis &axis, int color)
void Draw(float pot, std::string sample, std::string opt, osc::IOscCalc *calc, int from, int to, int style=1)
void Draw(float pot, std::string opt, osc::IOscCalc *calc, int from, int to, TLegend *leg, int style=1)
Spectrum with true energy information, allowing it to be oscillated
enum BeamMode kGreen
enum BeamMode kBlue
static Binning Simple(int n, double lo, double hi, const std::vector< std::string > &labels={})
Definition: Binning.cxx:107
SRContain contain
Output from SRContain (containment related variables)
Definition: SRIDBranch.h:51
std::vector< SRNeutrino > nu
implemented as a vector to maintain mc.nu structure, i.e. not a pointer, but with 0 or 1 entries...
Definition: SRTruthBranch.h:25
TCanvas * NewCan(std::string nameTitle)
Definition: numu_cut_flow.C:66
CutFlow(std::string name, SpectrumLoader &loader, const HistAxis &axis, std::vector< Sample > samples)
std::string fLongName
enum BeamMode string