nus17_fiducial_accounting.C
Go to the documentation of this file.
1 #include <bitset>
2 #include <string>
3 #include <vector>
4 
8 #include "CAFAna/Core/HistAxis.h"
9 #include "CAFAna/Core/Spectrum.h"
10 #include "CAFAna/Vars/Vars.h"
11 #include "CAFAna/Core/Binning.h"
12 #include "CAFAna/Core/Cut.h"
13 #include "CAFAna/Core/Utilities.h"
14 #include "CAFAna/Cuts/Cuts.h"
16 #include "NuXAna/Cuts/NusCuts17.h"
18 #include "OscLib/OscCalcDumb.h"
19 
20 #include "TFile.h"
21 #include "TF1.h"
22 #include "TH1.h"
23 #include "TH2.h"
24 #include "TEfficiency.h"
25 #include "TCanvas.h"
26 
27 
28 using namespace ana;
29 
30 const Cut kBoxCutIdeal(
31  [](const caf::SRProxy* sr)
32  {
33  if(sr->mc.nnu == 0) return false;
34  std::bitset<14> binary(sr->hdr.dibmask);
35  std::pair<int,int> planesA = calcFirstLastLivePlane(sr->slc.firstplane,
36  binary);
37  std::pair<int,int> planesB = calcFirstLastLivePlane(sr->slc.lastplane,
38  binary);
39  if((planesA.first != planesB.first) || (planesA.second != planesB.second))
40  return false;
41  if((planesA.second - planesA.first + 1)/64 < 14) return false;
42  if(sr->mc.nu[0].vtx.X() < -100) return false;
43  if(sr->mc.nu[0].vtx.X() > 100) return false;
44  if(sr->mc.nu[0].vtx.Y() < -100) return false;
45  if(sr->mc.nu[0].vtx.Y() > 100) return false;
46  if(sr->mc.nu[0].vtx.Z() < 1000) return false;
47  if(sr->mc.nu[0].vtx.Z() > 3500) return false;
48  return true;
49  });
50 
51 
52 const Cut kBoxCut(
53  [](const caf::SRProxy* sr)
54  {
55  if(sr->mc.nnu == 0) return false;
56  if(sr->mc.nu[0].vtx.X() < -762) return false;
57  if(sr->mc.nu[0].vtx.X() > 763) return false;
58  if(sr->mc.nu[0].vtx.Y() < -762) return false;
59  if(sr->mc.nu[0].vtx.Y() > 763) return false;
60  std::bitset<14> binary(sr->hdr.dibmask);
61  std::pair<int,int> planesA = calcFirstLastLivePlane(sr->slc.firstplane,
62  binary);
63  std::pair<int,int> planesB = calcFirstLastLivePlane(sr->slc.lastplane,
64  binary);
65  if((planesA.first != planesB.first) || (planesA.second != planesB.second))
66  return false;
67  if((planesA.second - planesA.first + 1)/64 < 4) return false;
68  int first = planesA.first/64;
69  int last = (planesA.second +1)/64;
70  if(sr->mc.nu[0].vtx.Z() < ((first)*426.0)) return false;
71  if(sr->mc.nu[0].vtx.Z() > ((last)*426.0)) return false;
72  return true;
73  });
74 
75 const Var kBlock(
76  [](const caf::SRProxy* sr)
77  {
78  int size = -1;
79  std::bitset<14> binary(sr->hdr.dibmask);
80  std::pair<int,int> planesA = calcFirstLastLivePlane(sr->slc.firstplane,
81  binary);
82  std::pair<int,int> planesB = calcFirstLastLivePlane(sr->slc.lastplane,
83  binary);
84  if((planesA.first != planesB.first) || (planesA.second != planesB.second))
85  return size;
86  size = (planesA.second - planesA.first + 1)/64;
87  return size;
88  });
89 
91  PredictionNoExtrap* sel,
93  TString title)
94 {
95  TH1::SetDefaultSumw2();
96  TH2::SetDefaultSumw2();
97  double kPOTFD = all->Predict(calc).POT();
98 
99  TH1* hSig = all->PredictComponent(calc,
101  Current::kNC,
102  Sign::kBoth).ToTH1(kPOTFD);
103  hSig->SetLineColor(kRed-10);
104 
105  TH1* hSigSel = sel->PredictComponent(calc,
107  Current::kNC,
108  Sign::kBoth).ToTH1(kPOTFD);
109 
110 
111  TH1D* hSigEff = new TH1D("hSigEff_"+title,
112  "Nus NC Selection Efficiency;Diblocks;Efficiency", 11, 3.5, 14.5);
113  TH1D* hN = new TH1D("hN_"+title,
114  "Nus NC Selection Efficiency;Diblocks;Efficiency", 11, 3.5, 14.5);
115  TH1D* hD = new TH1D("hD_"+title,
116  "Nus NC Selection Efficiency;Diblocks;Efficiency", 11, 3.5, 14.5);
117  for(unsigned int i = 5; i < 16; ++i){
118  for(unsigned int j = 0; j < hSigSel->GetBinContent(i); ++j){
119  hN->Fill(i-1);
120  }
121  for(unsigned int j = 0; j < hSig->GetBinContent(i); ++j){
122  hD->Fill(i-1);
123  }
124  }
125  hN->Write();
126  hD->Write();
127  hSigEff->Divide(hN, hD, 1.0, 1.0, "b");
128 
129  std::cout<<"Calculating efficiency for: "<<title<<std::endl;
130  for(unsigned int i = 1; i < 12; ++i){
131  float sigN = hN->GetBinContent(i);
132  float sigD = hD->GetBinContent(i);
133  float sigE = hSigEff->GetBinContent(i);
134  std::cout << "Signal Eff, DB: " << i+3
135  << " Selected: " << sigN
136  << " All: " << sigD
137  << " Efficiency: " << sigE << std::endl;
138  }
139  hSigEff->Write();
140 }
141 
143 {
144 
145  TH1::SetDefaultSumw2();
146  TH2::SetDefaultSumw2();
147 
148  SpectrumLoader loaderFDSwap("prod_caf_R17-03-01-prod3reco.l_fd_genie_fluxswap_fhc_nova_v08_full_v1");
149 
150  SpectrumLoader loaderFDNonSwap("prod_caf_R17-03-01-prod3reco.l_fd_genie_nonswap_fhc_nova_v08_full_v1");
151 
152  const Cut kNusEff = kNus17FDPresel && kNHit >= 25;
153 
154  const HistAxis axis("Diblocks", Binning::Simple(15, -0.5, 14.5), kBlock);
155 
156  // Prediction for all true neutrinos in fiducial volume
157  PredictionNoExtrap* predAll = new PredictionNoExtrap(
158  loaderFDNonSwap, loaderFDSwap,
159  axis.GetLabels()[0], axis.GetBinnings()[0], axis.GetVars()[0], kBoxCut
160  );
161 
162  // Prediction for all true neutrinos in fiducial volume that pass preselection
163  PredictionNoExtrap* predSel = new PredictionNoExtrap(
164  loaderFDNonSwap, loaderFDSwap,
165  axis.GetLabels()[0], axis.GetBinnings()[0], axis.GetVars()[0], kNusEff && kBoxCut
166  );
167 
168  // Pprediction for all true neutrinos in small box at center of detector
169  PredictionNoExtrap* predIdealAll = new PredictionNoExtrap(
170  loaderFDNonSwap, loaderFDSwap,
171  axis.GetLabels()[0], axis.GetBinnings()[0], axis.GetVars()[0], kBoxCutIdeal
172  );
173  // Prediction for all true neutrinos in fiducial volume that pass preselection
174  // in small box at center of detector
175  PredictionNoExtrap* predIdealSel = new PredictionNoExtrap(
176  loaderFDNonSwap, loaderFDSwap,
177  axis.GetLabels()[0], axis.GetBinnings()[0], axis.GetVars()[0], kNusEff && kBoxCutIdeal
178  );
179 
180  loaderFDSwap.SetSpillCut(kStandardSpillCuts);
181  loaderFDNonSwap.SetSpillCut(kStandardSpillCuts);
182 
183  // Load from files
184  loaderFDNonSwap.Go();
185  loaderFDSwap.Go();
186 
187  // Open output file
188  TFile file(outfile.c_str(),"RECREATE");
189 
191 
192  // Caclulate preselection efficieny as a function of diblock
193  GetExtrap(predAll, predSel, &calc, "FidMass");
194  // Caclulate preselection efficieny for small window in ideal detector
195  GetExtrap(predIdealAll, predIdealSel, &calc, "SmallWindow");
196 
197  file.Close();
198 }
Spectrum PredictComponent(osc::IOscCalc *calc, Flavors::Flavors_t flav, Current::Current_t curr, Sign::Sign_t sign) const override
enum BeamMode kRed
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
const std::vector< T > & GetVars() const
Definition: HistAxis.h:92
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
caf::Proxy< caf::SRHeader > hdr
Definition: SRProxy.h:2137
Proxy for caf::StandardRecord.
Definition: SRProxy.h:2126
General interface to oscillation calculators.
Definition: StanTypedefs.h:22
caf::Proxy< std::vector< caf::SRNeutrino > > nu
Definition: SRProxy.h:618
caf::Proxy< short int > nnu
Definition: SRProxy.h:617
void SetSpillCut(const SpillCut &cut)
osc::OscCalcDumb calc
void nus17_fiducial_accounting(std::string outfile)
Spectrum Predict(osc::IOscCalc *calc) const override
const Cut kBoxCut([](const caf::SRProxy *sr){if(sr->mc.nnu==0) return false;if(sr->mc.nu[0].vtx.x< -762) return false;if(sr->mc.nu[0].vtx.x > 763) return false;if(sr->mc.nu[0].vtx.y< -762) return false;if(sr->mc.nu[0].vtx.y > 763) return false;std::bitset< 14 > binary(sr->hdr.dibmask);std::pair< int, int > planesA=calcFirstLastLivePlane(sr->slc.firstplane, binary);std::pair< int, int > planesB=calcFirstLastLivePlane(sr->slc.lastplane, binary);if((planesA.first!=planesB.first)||(planesA.second!=planesB.second)) return false;if((planesA.second-planesA.first+1)/64< 4) return false;int first=planesA.first/64;int last=(planesA.second+1)/64;if(sr->mc.nu[0].vtx.z< ((first)*426.0)) return false;if(sr->mc.nu[0].vtx.z >((last)*426.0)) return false;return true;})
Definition: AnalysisMasks.h:25
const Var kNHit
Definition: Vars.cxx:71
void GetExtrap(PredictionNoExtrap *all, PredictionNoExtrap *sel, osc::IOscCalc *calc, TString title)
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.
caf::Proxy< unsigned int > lastplane
Definition: SRProxy.h:1309
const double j
Definition: BetheBloch.cxx:29
const Var kBlock([](const caf::SRProxy *sr){int size=-1;std::bitset< 14 > binary(sr->hdr.dibmask);std::pair< int, int > planesA=calcFirstLastLivePlane(sr->slc.firstplane, binary);std::pair< int, int > planesB=calcFirstLastLivePlane(sr->slc.lastplane, binary);if((planesA.first!=planesB.first)||(planesA.second!=planesB.second)) return size;size=(planesA.second-planesA.first+1)/64;return size;})
double POT() const
Definition: Spectrum.h:227
const Cut kBoxCutIdeal([](const caf::SRProxy *sr){if(sr->mc.nnu==0) return false;std::bitset< 14 > binary(sr->hdr.dibmask);std::pair< int, int > planesA=calcFirstLastLivePlane(sr->slc.firstplane, binary);std::pair< int, int > planesB=calcFirstLastLivePlane(sr->slc.lastplane, binary);if((planesA.first!=planesB.first)||(planesA.second!=planesB.second)) return false;if((planesA.second-planesA.first+1)/64< 14) return false;if(sr->mc.nu[0].vtx.x< -100) return false;if(sr->mc.nu[0].vtx.x > 100) return false;if(sr->mc.nu[0].vtx.y< -100) return false;if(sr->mc.nu[0].vtx.y > 100) return false;if(sr->mc.nu[0].vtx.z< 1000) return false;if(sr->mc.nu[0].vtx.z > 3500) return false;return true;})
A cut to select only events in a very small central region of the detector.
Definition: AnalysisMasks.h:21
OStream cout
Definition: OStream.cxx:6
caf::Proxy< caf::SRTruthBranch > mc
Definition: SRProxy.h:2138
const std::vector< Binning > & GetBinnings() const
Definition: LabelsAndBins.h:69
caf::Proxy< caf::SRSlice > slc
Definition: SRProxy.h:2142
Collaborates with Spectrum and OscillatableSpectrum to fill spectra from CAF files.
std::pair< int, int > calcFirstLastLivePlane(int plane, std::bitset< 14 > binary, caf::Det_t det)
Neutral-current interactions.
Definition: IPrediction.h:40
TFile * file
Definition: cellShifts.C:17
Both neutrinos and antineutrinos.
Definition: IPrediction.h:52
caf::Proxy< unsigned int > firstplane
Definition: SRProxy.h:1305
const SpillCut kStandardSpillCuts
Apply this unless you&#39;re doing something special.
Definition: SpillCuts.h:49
caf::Proxy< short unsigned int > dibmask
Definition: SRProxy.h:235
Prediction that just uses FD MC, with no extrapolation.
All neutrinos, any flavor.
Definition: IPrediction.h:26
const Cut kNus17FDPresel
The Nus17 preselection cuts for the Far Detector from docdb 21113.
Definition: NusCuts17.h:38
static Binning Simple(int n, double lo, double hi, const std::vector< std::string > &labels={})
Definition: Binning.cxx:107
FILE * outfile
Definition: dump_event.C:13
const std::vector< std::string > & GetLabels() const
Definition: LabelsAndBins.h:68
enum BeamMode string