fd_eff_vs_pos.C
Go to the documentation of this file.
2 #include "CAFAna/Core/Ratio.h"
8 using namespace ana;
9 
11 
12 #include "TCanvas.h"
13 #include "TFile.h"
14 #include "TH2.h"
15 
16 const Cut kHasVtx([](const caf::SRProxy* sr){return sr->vtx.elastic.IsValid;});
17 
18 // OK to have no checks because presel will guarantee vertex exists
19 const Var kVtxX([](const caf::SRProxy* sr){return sr->vtx.elastic.vtx.X();});
20 const Var kVtxY([](const caf::SRProxy* sr){return sr->vtx.elastic.vtx.Y();});
21 
22 const Var kNHits = SIMPLEVAR(slc.nhit);
23 
24 void fill()
25 {
26  SpectrumLoader loaderNonswap("prod_caf_S15-05-22_fd_genie_fhc_nonswap_matchedruns");
27  SpectrumLoader loaderSwap("prod_caf_S15-05-22_fd_genie_fhc_fluxswap_matchedruns");
28  // SpectrumLoader loaderNonswap("prod_decaf_S15-05-22a_fd_genie_fhc_nonswap_nue_contain_matchedruns");
29  // SpectrumLoader loaderSwap("prod_decaf_S15-05-22a_fd_genie_fhc_fluxswap_nue_contain_matchedruns");
30  loaderNonswap.SetSpillCut(kStandardSpillCuts);
31  loaderSwap.SetSpillCut(kStandardSpillCuts);
32 
33  const Binning bins1D = Binning::Simple(160, -800, +800);
34  const Binning bins2D = Binning::Simple(80*80, 0, 80*80);
35 
36  const Var kVtxXY = Var2D(kVtxX, 80, -800, +800, kVtxY, 80, -800, +800);
37 
38  std::map<std::string, IPrediction*> preds;
39 
40  for(int pid = 0; pid < 2; ++pid){
41  for(int view = 0; view < 3; ++view){
42  std::string title, varStr;
44  if(view == 0){title = "X (cm)"; var = kVtxX; varStr = "x";}
45  if(view == 1){title = "Y (cm)"; var = kVtxY; varStr = "y";}
46  if(view == 2){title = "XY"; var = kVtxXY; varStr = "xy";}
47  Binning bins = bins1D;
48  if(view == 2) bins = bins2D;
49  for(int selIdx = 0; selIdx < 3; ++selIdx){
50  Cut cut = kHasVtx;
51  if(pid == 0 && selIdx == 1) cut = kNueFirstAnaFullLEMPresel;
52  if(pid == 1 && selIdx == 1) cut = kNueFirstAnaFullLIDPresel;
53  if(pid == 0 && selIdx == 2) cut = kNueFDLEM;
54  if(pid == 1 && selIdx == 2) cut = kNueFDLID;
55  std::string selStr;
56  if(selIdx == 0) selStr = "nosel";
57  if(selIdx == 1) selStr = "presel";
58  if(selIdx == 2) selStr = "pidcut";
59 
60  for(int varIdx = 0; varIdx < 3; ++varIdx){
61  Var wei = kUnweighted;
62  std::string weiStr = "count";
63  if(varIdx == 1){wei = kNHits; weiStr = "nhits";}
64  if(varIdx == 2){wei = kCaloE; weiStr = "calE";}
65  const std::string key =
66  TString::Format("%s_%s_%s_%s",
67  (pid == 0) ? "lem" : "lid",
68  selStr.c_str(),
69  varStr.c_str(),
70  weiStr.c_str()).Data();
71  preds[key] = new PredictionNoExtrap(loaderNonswap, loaderSwap,
72  title, bins,
73  var, cut, kNoShift, wei);
74  }
75  }
76  }
77  }
78 
79  loaderNonswap.Go();
80  loaderSwap.Go();
81 
82  TFile fout("fd_eff_vs_pos.root", "RECREATE");
83 
84  for(auto it: preds) it.second->SaveTo(fout.mkdir(it.first.c_str()));
85 }
86 
87 void ratio(std::string prefix, std::string suffix, bool preselDenom)
88 {
89  const std::string fname = "fd_eff_vs_pos.root";
90 
91  IPrediction* predPresel = LoadFromFile<IPrediction>(fname, prefix+(preselDenom ? "_presel_" : "_nosel_")+suffix+"_count").release();
92  IPrediction* predPID = LoadFromFile<IPrediction>(fname, prefix+"_pidcut_"+suffix+"_count").release();
93 
95 
96  TH1* numSig = predPID->PredictComponent(&calc, Flavors::kNuMuToNuE, Current::kCC, Sign::kBoth).ToTH1(1e20);
97  TH1* denomSig = predPresel->PredictComponent(&calc, Flavors::kNuMuToNuE, Current::kCC, Sign::kBoth).ToTH1(1e20);
98 
99  TH1* numBkg = predPID->Predict(&calc).ToTH1(1e20);
100  TH1* denomBkg = predPresel->Predict(&calc).ToTH1(1e20);
101 
102  numBkg->Add(numSig, -1);
103  denomBkg->Add(denomSig, -1);
104 
105  std::cout << prefix << " " << suffix << std::endl;
106  // std::cout << numSig->GetXaxis()->GetBinCenter(21) << " "
107  // << numSig->GetXaxis()->GetBinCenter(80) << " "
108  // << numSig->GetXaxis()->GetBinCenter(81) << " "
109  // << numSig->GetXaxis()->GetBinCenter(140) << std::endl;
110  const double effNom = numSig->Integral(81, 140)/denomSig->Integral(81, 140);
111  const double effLow = numSig->Integral(21, 80)/denomSig->Integral(21, 80);
112  std::cout << effNom << " " << effLow << std::endl;
113 
114  numSig->Rebin(5);
115  denomSig->Rebin(5);
116 
117  numBkg->Rebin(5);
118  denomBkg->Rebin(5);
119 
120  numSig->Divide(numSig, denomSig, 1, 1, "B");
121  numSig->Scale(100);
122  numSig->GetYaxis()->SetTitle("#nu_{e} signal selection efficiency (%)");
123  numSig->GetYaxis()->SetRangeUser(0, preselDenom ? 55 : 30);
124  numSig->GetXaxis()->CenterTitle();
125  numSig->GetYaxis()->CenterTitle();
126 
127  numBkg->Divide(numBkg, denomBkg, 1, 1, "B");
128  numBkg->Scale(100);
129  numBkg->GetYaxis()->SetTitle("Background selection efficiency (%)");
130  numBkg->GetYaxis()->SetRangeUser(0, 6);
131  numBkg->GetXaxis()->CenterTitle();
132  numBkg->GetYaxis()->CenterTitle();
133 
134  if(!preselDenom) suffix += "_over_nocut";
135 
136  new TCanvas;
137  numSig->Draw("ep");
138  gPad->Print(("sig_eff_"+prefix+"_"+suffix+".png").c_str());
139  gPad->Print(("sig_eff_"+prefix+"_"+suffix+".eps").c_str());
140 
141  new TCanvas;
142  numBkg->Draw("ep");
143  gPad->Print(("bkg_eff_"+prefix+"_"+suffix+".png").c_str());
144  gPad->Print(("bkg_eff_"+prefix+"_"+suffix+".eps").c_str());
145 }
146 
147 void ratio2D(std::string prefix, bool preselDenom)
148 {
149  const std::string fname = "fd_eff_vs_pos.root";
150 
151  IPrediction* predPresel = LoadFromFile<IPrediction>(fname, prefix+(preselDenom ? "_presel_xy_count" : "_nosel_xy_count")).release();
152  IPrediction* predPID = LoadFromFile<IPrediction>(fname, prefix+"_pidcut_xy_count").release();
153 
155 
156  const Binning binsX = Binning::Simple(80, -800, +800);
157  const Binning binsY = Binning::Simple(80, -800, +800);
158 
159  TH2* numSig = ToTH2(predPID->PredictComponent(&calc, Flavors::kNuMuToNuE, Current::kCC, Sign::kBoth), 1e20, kPOT, binsX, binsY);
160  TH2* denomSig = ToTH2(predPresel->PredictComponent(&calc, Flavors::kNuMuToNuE, Current::kCC, Sign::kBoth), 1e20, kPOT, binsX, binsY);
161 
162  std::cout << prefix << " xy" << std::endl;
163  // std::cout << numSig->GetXaxis()->GetBinCenter(11) << " "
164  // << numSig->GetXaxis()->GetBinCenter(40) << " "
165  // << numSig->GetXaxis()->GetBinCenter(41) << " "
166  // << numSig->GetXaxis()->GetBinCenter(70) << std::endl;
167  const double effNom = numSig->Integral(41, 70, 41, 70)/denomSig->Integral(41, 70, 41, 70);
168  const double effLow = numSig->Integral(11, 40, 11, 40)/denomSig->Integral(11, 40, 11, 40);
169  const double effMed1 = numSig->Integral(11, 40, 41, 70)/denomSig->Integral(11, 40, 41, 70);
170  const double effMed2 = numSig->Integral(41, 70, 11, 40)/denomSig->Integral(41, 70, 11, 40);
171  std::cout << effNom << " " << effLow << " " << effMed1 << " " << effMed2 << std::endl;
172 
173  numSig->Rebin2D(2, 2);
174  denomSig->Rebin2D(2, 2);
175 
176  numSig->Divide(denomSig);
177  numSig->Scale(100);
178 
179  numSig->GetXaxis()->SetTitle("X (cm)");
180  numSig->GetYaxis()->SetTitle("Y (cm)");
181 
182  numSig->GetXaxis()->CenterTitle();
183  numSig->GetYaxis()->CenterTitle();
184 
185  numSig->GetZaxis()->SetRangeUser(preselDenom ? 30 : 15,
186  preselDenom ? 55 : 30);
187 
188  std::string suffix = (preselDenom ? "": "_over_nocut");
189  new TCanvas;
190  numSig->Draw("colz");
191  gPad->Print(("sig_eff_"+prefix+"_xy"+suffix+".png").c_str());
192  gPad->Print(("sig_eff_"+prefix+"_xy"+suffix+".eps").c_str());
193 }
194 
196 {
197  const std::string fname = "fd_eff_vs_pos.root";
198 
199  IPrediction* predNHit = LoadFromFile<IPrediction>(fname, prefix+(presel ? "_presel_" : "_pidcut_")+suffix+"_nhits").release();
200  IPrediction* predCount = LoadFromFile<IPrediction>(fname, prefix+(presel ? "_presel_" : "_pidcut_")+suffix+"_count").release();
201 
203 
204  TH1* num = predNHit->Predict(&calc).ToTH1(1e20);
205  TH1* denom = predCount->Predict(&calc).ToTH1(1e20);
206 
207  num->Rebin(5);
208  denom->Rebin(5);
209 
210  num->Divide(denom);
211  num->GetYaxis()->SetTitle("Average number of hits");
212  num->GetYaxis()->SetRangeUser(50, 90);
213 
214  new TCanvas;
215  num->Draw("ep");
216  gPad->Print(("nhits_"+prefix+"_"+suffix+(presel ? "_presel" : "")+".png").c_str());
217  gPad->Print(("nhits_"+prefix+"_"+suffix+(presel ? "_presel" : "")+".eps").c_str());
218 }
219 
221 {
222  const std::string fname = "fd_eff_vs_pos.root";
223 
224  IPrediction* predNHit = LoadFromFile<IPrediction>(fname, prefix+(presel ? "_presel_" : "_pidcut_")+"xy_nhits").release();
225  IPrediction* predCount = LoadFromFile<IPrediction>(fname, prefix+(presel ? "_presel_" : "_pidcut_")+"xy_count").release();
226 
228 
229  const Binning binsX = Binning::Simple(80, -800, +800);
230  const Binning binsY = Binning::Simple(80, -800, +800);
231 
232  TH2* num = ToTH2(predNHit->Predict(&calc), 1e20, binsX, binsY);
233  TH2* denom = ToTH2(predCount->Predict(&calc), 1e20, binsX, binsY);
234 
235  num->Rebin2D(2, 2);
236  denom->Rebin2D(2, 2);
237 
238  num->Divide(denom);
239  num->GetXaxis()->SetTitle("X (cm)");
240  num->GetYaxis()->SetTitle("Y (cm)");
241  num->GetZaxis()->SetRangeUser(50, 90);
242 
243  new TCanvas;
244  num->Draw("colz");
245  gPad->Print(("nhits_"+prefix+"_xy"+(presel ? "_presel" : "")+".png").c_str());
246  gPad->Print(("nhits_"+prefix+"_xy"+(presel ? "_presel" : "")+".eps").c_str());
247 }
248 
249 void plot()
250 {
251  ratio("lem", "x", true);
252  ratio("lem", "y", true);
253  ratio("lem", "x", false);
254  ratio("lem", "y", false);
255  ratio2D("lem", true);
256  ratio2D("lem", false);
257 
258  ratio("lid", "x", true);
259  ratio("lid", "y", true);
260  ratio("lid", "x", false);
261  ratio("lid", "y", false);
262  ratio2D("lid", true);
263  ratio2D("lid", false);
264 
265  nhits("lem", "x", false);
266  nhits("lem", "y", false);
267  nhits("lid", "x", false);
268  nhits("lid", "y", false);
269  nhits2D("lem", false);
270  nhits2D("lid", false);
271 
272  nhits("lem", "x", true);
273  nhits("lem", "y", true);
274  nhits("lid", "x", true);
275  nhits("lid", "y", true);
276  nhits2D("lem", true);
277  nhits2D("lid", true);
278 }
279 
280 void fd_eff_vs_pos(bool doFill)
281 {
282  if(doFill) fill();
283  plot();
284 }
Preselection Object.
Definition: FillPIDs.h:19
void fill()
Definition: fd_eff_vs_pos.C:24
SRElasticProxy elastic
Definition: SRProxy.h:2229
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.
set< int >::iterator it
void fd_eff_vs_pos(bool doFill)
(&#39; appearance&#39;)
Definition: IPrediction.h:18
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
virtual Spectrum PredictComponent(osc::IOscCalculator *calc, Flavors::Flavors_t flav, Current::Current_t curr, Sign::Sign_t sign) const =0
const Var kNHits([](const caf::SRProxy *sr){float nhit=sr->slc.nhit;return nhit;})
float X() const
Definition: SRProxy.h:951
void nhits(std::string prefix, std::string suffix, bool presel)
var
void plot()
Proxy for StandardRecord.
Definition: SRProxy.h:2237
void SetSpillCut(const SpillCut &cut)
float Y() const
Definition: SRProxy.h:952
const Var kVtxX
virtual Spectrum Predict(osc::IOscCalculator *calc) const =0
const Cut kNueFirstAnaFullLIDPresel
From Tian&#39;s docdb 12673.
Charged-current interactions.
Definition: IPrediction.h:39
const Var kCaloE
Summed calorimetric energy of all hits in slice, uncorrected.
Definition: Vars.cxx:50
Simple oscillation probability calculator that has no solar term or mass hierarchy or delta so it&#39;s s...
Proxy< bool > IsValid
Definition: SRProxy.h:2216
const Cut kNueFDLEM
virtual void Go() override
Load all the registered spectra.
const XML_Char * prefix
Definition: expat.h:380
TH2 * ToTH2(const Spectrum &s, double exposure, ana::EExposureType expotype, const Binning &binsx, const Binning &binsy, ana::EBinType bintype)
Definition: Utilities.cxx:374
void ratio2D(std::string prefix, bool preselDenom)
const SystShifts kNoShift
Definition: SystShifts.h:112
OStream cout
Definition: OStream.cxx:6
#define SIMPLEVAR(CAFNAME)
For Vars where literally all you need is a single CAF variable.
Definition: Var.h:94
const Binning bins
Definition: NumuCC_CPiBin.h:8
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
const Cut cut
Definition: exporter_fd.C:30
const Var kVtxY
Collaborates with Spectrum and OscillatableSpectrum to fill spectra from CAF files.
int num
Definition: f2_nu.C:119
Both neutrinos and antineutrinos.
Definition: IPrediction.h:52
Standard interface to all prediction techniques.
Definition: IPrediction.h:57
const SpillCut kStandardSpillCuts
Apply this unless you&#39;re doing something special.
Definition: SpillCuts.h:28
void Format(TGraph *gr, int lcol, int lsty, int lwid, int mcol, int msty, double msiz)
Definition: Style.cxx:154
Prediction that just uses FD MC, with no extrapolation.
void nhits2D(std::string prefix, bool presel)
GenericVar< T > Var2D(const GenericVar< T > &a, const Binning &binsa, const GenericVar< T > &b, const Binning &binsb)
Variable formed from two input variables.
Definition: Var.cxx:245
const Cut kNueFDLID
loaderSwap
Definition: demo4.py:9
static Binning Simple(int n, double lo, double hi, const std::vector< std::string > &labels={})
Definition: Binning.cxx:38
SRVector3DProxy vtx
Definition: SRProxy.h:2215
SRVertexBranchProxy vtx
Definition: SRProxy.h:2249
const Var kUnweighted
The simplest possible Var, always 1. Used as a default weight.
Definition: Var.h:100
const Cut kHasVtx([](const caf::SRProxy *sr){return sr->vtx.elastic.IsValid;})
static constexpr Double_t sr
Definition: Munits.h:164
void ratio(std::string prefix, std::string suffix, bool preselDenom)
Definition: fd_eff_vs_pos.C:87
const Cut kNueFirstAnaFullLEMPresel
From Tian&#39;s docdb 12673.