efficiencySA.C
Go to the documentation of this file.
1 // Fills spectra for systematic error band
2 
3 #include "CAFAna/Cuts/Cuts.h"
4 #include "CAFAna/Cuts/NumuCuts.h"
5 #include "CAFAna/Core/Binning.h"
6 #include "CAFAna/Core/Loaders.h"
7 #include "CAFAna/Core/Spectrum.h"
8 #include "CAFAna/Core/Var.h"
11 #include "CAFAna/Systs/Systs.h"
12 #include "CAFAna/Systs/NumuSysts.h"
13 #include "CAFAna/Vars/NumuVars.h"
15 #include "CAFAna/Cuts/SpillCuts.h"
16 #include "CAFAna/Analysis/Calcs.h"
17 #include "CAFAna/Analysis/Plots.h"
18 #include "CAFAna/Analysis/Style.h"
19 #include "CAFAna/Vars/HistAxes.h"
20 
22 
23 #include "TCanvas.h"
24 #include "TFile.h"
25 #include "TGraph.h"
26 #include "TGraphAsymmErrors.h"
27 #include "TH1.h"
28 #include "TH2.h"
29 #include "TMath.h"
30 #include "TGaxis.h"
31 #include "TMultiGraph.h"
32 #include "TLegend.h"
33 #include "TLatex.h"
34 #include "TStyle.h"
35 #include "THStack.h"
36 #include "TPaveText.h"
37 #include "TList.h"
38 #include "TGaxis.h"
39 #include "TAttLine.h"
40 #include "TAttMarker.h"
41 
42 #include <cmath>
43 #include <iostream>
44 #include <vector>
45 #include <list>
46 #include <sstream>
47 #include <string>
48 #include <sstream>
49 #include <fstream>
50 #include <iomanip>
51 
52 const double pot = 6e20;
53 
54 using namespace ana;
55 
57 {
58 
59  std::string fd_nonswap = "bzamoran_prod_caf_R16-03-03-prod2reco.f_fd_genie_nonswap_fhc_nova_v08_epoch1-3c_v1_prod2-snapshot_2pc";
60 
63 
65 
68 
69  const Binning kEnergyBinning = Binning::Simple(20,0,5);
70 
71  //// ------ Now the cuts ------ /////
72 
73  const Cut kAny({"mc.nnu"}, // Any mode
74  [](const caf::StandardRecord* sr)
75  {
76  return (sr->mc.nnu > 0);
77  });
78  const Cut kQE({"mc.nnu", "mc.nu.*"}, // QE
79  [](const caf::StandardRecord* sr)
80  {
81  return sr->mc.nnu > 0 && sr->mc.nu[0].mode == 0;
82  });
83  const Cut kRes({"mc.nnu", "mc.nu.*"}, // Res
84  [](const caf::StandardRecord* sr)
85  {
86  return sr->mc.nnu > 0 && sr->mc.nu[0].mode == 1;
87  });
88  const Cut kDIS({"mc.nnu", "mc.nu.*"}, // DIS
89  [](const caf::StandardRecord* sr)
90  {
91  return sr->mc.nnu > 0 && sr->mc.nu[0].mode == 2;
92  });
93  const Cut kCoh({"mc.nnu", "mc.nu.*"}, // Coh
94  [](const caf::StandardRecord* sr)
95  {
96  return sr->mc.nnu > 0 && sr->mc.nu[0].mode == 3;
97  });
98  const Cut kMEC({"mc.nnu", "mc.nu.*"}, // MEC
99  [](const caf::StandardRecord* sr)
100  {
101  return sr->mc.nnu > 0 && sr->mc.nu[0].mode == 10;
102  });
103 
104  const Cut kTruthContainFD({"mc.nnu","mc.nu.*","sel.contain.missE"}, // only approx
105  [](const caf::StandardRecord* sr)
106  {
107  return sr->mc.nnu > 0 && sr->mc.nu[0].E > 0 && sr->sel.contain.missE / sr->mc.nu[0].E < 0.01;
108  });
109 
110  struct Modes
111  {
113  Cut cut;
114  };
115 
116  const int kNumModes = 6;
117 
118  Modes modes[kNumModes] = {
119  {"All", kAny},
120  {"QE", kQE},
121  {"Res", kRes},
122  {"DIS", kDIS},
123  {"Coh", kCoh},
124  {"MEC", kMEC}
125  };
126 
127  struct Plot
128  {
131  Binning bins;
132  Var var;
133  };
134 
135  const int kNumPlots = 2;
136 
137  Plot plots[kNumPlots] = {
138  {";Reconstructed Neutrino Energy (GeV);Efficiency", "kCCE", kEnergyBinning, kCCE},
139  {";True Neutrino Energy (GeV);Efficiency", "kTrueE", kEnergyBinning, kTrueE}
140  };
141 
142  struct Selection
143  {
146  Cut cut;
147  };
148 
149  const int kNumCuts = 5;
150 
151  const Cut kBasicCut = kIsNumuCC && kTruthContainFD;
152 
153  Selection cuts[kNumCuts] = {
154  {"No cuts","NoCut", kBasicCut},
155  {"Data quality","DQ", kBasicCut && kNumuQuality},
156  {"+ containment","DQ_contain", kBasicCut && kNumuQuality && kNumuContainFD},
157  {"+ PID","DQ_contain_NCRej", kBasicCut && kNumuQuality && kNumuContainFD && kNumuNCRej},
158  {"+ cosmic rej.","NumuFD", kBasicCut &&kNumuQuality && kNumuContainFD && kNumuNCRej && kNumuCosmicRej}
159  };
160 
161  int colours[kNumCuts] = {kBlack, kBlue, kGreen+2, kOrange-2, kRed};
162 
163  PredictionNoExtrap* sSpect[kNumPlots][kNumCuts][kNumModes];
164 
165  for(int k = 0; k < kNumModes; k++)
166  {
167  Cut mode = modes[k].cut;
168 
169  for(int i = 0; i < kNumPlots; i++)
170  {
171  Plot p = plots[i];
172 
173  for(int j = 0; j < kNumCuts; j++)
174  {
175  Cut cut = cuts[j].cut;
176 
177  sSpect[i][j][k] = new PredictionNoExtrap(loaders, p.label, p.bins, p.var, cut && mode, kNoShift, kTuftsWeightCC );
178  } // loop over cuts
179  } // loop over plots
180  } // loop over modes
181 
182  loaders.Go();
183 
184  TH1* hHist[kNumPlots][kNumCuts][kNumModes];
185 
186  for(int k = 0; k < kNumModes; k++)
187  {
188  for(int i = 0; i < kNumPlots; i++)
189  {
190  for(int j = 0; j < kNumCuts; j++)
191  {
192  std::string histName = plots[i].fname+"_"+cuts[j].fname;
193 
194  // hHist[i][j][k] = sSpect[i][j][k]->Predict(&calc).ToTH1(pot);
195  hHist[i][j][k] = sSpect[i][j][k]->PredictUnoscillated().ToTH1(pot);
196  hHist[i][j][k]->SetName(histName.c_str());
197 
198  hHist[i][j][k]->SetMarkerColor(colours[j]);
199  hHist[i][j][k]->SetLineColor(colours[j]);
200  } // end loop over cuts
201  } // end loop over plots
202  } // end loop over modes
203 
204  for(int k = 0; k < kNumModes; k++)
205  {
206  for(int i = 0; i < kNumPlots; i++)
207  {
208  std::string canvasName = plots[i].fname+"_"+modes[k].fname;
209  TCanvas* c1 = new TCanvas("",canvasName.c_str());
210 
211  TLegend* leg = new TLegend(0.65,0.18,0.85,0.40);
212  leg->SetFillColor(10);
213  leg->SetLineWidth(0);
214  leg->SetLineColor(10);
215 
216  for (int j = 1; j < kNumCuts; j++)
217  {
218  hHist[i][j][k]->Divide(hHist[i][0][k]);
219 
220  std::cout << modes[k].fname.c_str() << " "
221  << plots[i].fname.c_str() << " "
222  << cuts[j].fname.c_str() << " : "
223  << hHist[i][j][k]->Integral() << std::endl;
224 
225  if(j == 1) hHist[i][j][k]->Draw("hist");
226  else hHist[i][j][k]->Draw("hist same");
227  leg->AddEntry(hHist[i][j][k], cuts[j].label.c_str(), "l");
228  } // loop over cuts
229 
230  leg->Draw();
231  c1->SaveAs(("effSA/"+canvasName+".pdf").c_str());
232  c1->SaveAs(("effSA/"+canvasName+".root").c_str());
233  } // loop over plots
234  } // loop over modes
235 
236 }
const Cut kIsNumuCC(CCFlavSel(14, 14))
Select CC .
Far Detector at Ash River.
Definition: SREnums.h:11
osc::OscCalculatorDumb calc
Represent the binning of a Spectrum&#39;s x-axis.
Definition: Binning.h:16
Oscillation analysis framework, runs over CAF files outside of ART.
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:553
std::string label
Definition: CutFlow_Data.C:28
std::string fname
Definition: getTimePeak.C:83
Optimized version of OscCalculatorPMNS.
Definition: StanTypedefs.h:28
const char * p
Definition: xmltok.h:285
Collection of SpectrumLoaders for many configurations.
Definition: Loaders.h:17
virtual Spectrum PredictUnoscillated() const
Definition: IPrediction.cxx:82
const Cut kNumuCosmicRej([](const caf::SRProxy *sr){return(sr->sel.cosrej.anglekal > 0.5 && sr->sel.cosrej.numucontpid2019 > 0.535 && sr->slc.nhit< 400);})
Definition: NumuCuts.h:31
var
void ResetOscCalcToDefault(osc::IOscCalculatorAdjustable *calc)
Reset calculator to default assumptions for all parameters.
Definition: Calcs.cxx:23
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:19
void Go()
Call Go() on all the loaders.
Definition: Loaders.cxx:170
const char * label
void SetSpillCut(const SpillCut &cut)
Definition: Loaders.cxx:121
====================================================================== ///
Definition: CutFlow_Data.C:27
void efficiencySA()
Definition: efficiencySA.C:56
const Var kTrueE([](const caf::SRProxy *sr){return(sr->mc.nnu==0)?0.:float(sr->mc.nu[0].E);})
Definition: Vars.cxx:87
const double pot
Definition: efficiencySA.C:52
const int kNumPlots
Definition: GetSpectra.h:1
const Var kCCE
Definition: Vars.h:90
const Binning kEnergyBinning
const double j
Definition: BetheBloch.cxx:29
const Cut kNumuNCRej([](const caf::SRProxy *sr){return(sr->sel.remid.pid >0.75);})
Definition: NumuCuts.h:23
const SystShifts kNoShift
Definition: SystShifts.h:112
OStream cout
Definition: OStream.cxx:6
const std::vector< PlotDef > plots
const Binning bins
Definition: NumuCC_CPiBin.h:8
std::string fname
Definition: CutFlow_Data.C:29
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
The StandardRecord is the primary top-level object in the Common Analysis File trees.
const Cut cut
Definition: exporter_fd.C:30
const Var kTuftsWeightCC
Definition: XsecTunes.h:32
std::vector< Loaders * > loaders
Definition: syst_header.h:385
c1
Definition: demo5.py:24
const SpillCut kStandardSpillCuts
Apply this unless you&#39;re doing something special.
Definition: SpillCuts.h:28
Binning bins
Definition: CutFlow_Data.C:30
Prediction that just uses FD MC, with no extrapolation.
const Cut kNumuQuality
Definition: NumuCuts.h:17
Var var
Definition: CutFlow_Data.C:31
void SetLoaderPath(const std::string &path, caf::Det_t det, DataMC datamc, DataSource src=kBeam, SwappingConfig swap=kNonSwap)
Configure loader via wildcard path.
Definition: Loaders.cxx:25
static Binning Simple(int n, double lo, double hi, const std::vector< std::string > &labels={})
Definition: Binning.cxx:38
static constexpr Double_t sr
Definition: Munits.h:164