fiducial_accounting.C
Go to the documentation of this file.
1 #include <bitset>
5 #include "CAFAna/Core/HistAxis.h"
6 #include "CAFAna/Core/Loaders.h"
7 #include "CAFAna/Vars/Vars.h"
8 #include "CAFAna/Core/Binning.h"
9 #include "CAFAna/Core/Cut.h"
10 #include "CAFAna/Cuts/Cuts.h"
12 #include "CAFAna/Cuts/NueCutsSecondAna.h"
16 #include "CAFAna/Decomp/NCDecomp.h"
21 #include "TFile.h"
22 #include "OscLib/OscCalcDumb.h"
23 #include "CAFAna/Core/Spectrum.h"
24 #include "TF1.h"
25 #include "TH1.h"
26 #include "TH2.h"
27 #include "TEfficiency.h"
29 #include "TCanvas.h"
31 #include "CAFAna/Vars/FitVars.h"
32 #include "CAFAna/Core/Utilities.h"
33 #include "Utilities/func/MathUtil.h"
35 
36 using namespace ana;
37 
38 const Cut kBoxCutIdeal(
39  [](const caf::SRProxy* sr)
40  {
41  if (sr->mc.nnu == 0) return false;
42  std::bitset<14> binary(sr->hdr.dibmask);
43  std::pair<int,int> planesA = calcFirstLastLivePlane(sr->slc.firstplane,
44  binary);
45  std::pair<int,int> planesB = calcFirstLastLivePlane(sr->slc.lastplane,
46  binary);
47  if ((planesA.first != planesB.first) || (planesA.second != planesB.second))
48  return false;
49  if ((planesA.second - planesA.first + 1)/64 < 14) return false;
50  if (sr->mc.nu[0].vtx.x < - 100) return false;
51  if (sr->mc.nu[0].vtx.x > 100) return false;
52  if (sr->mc.nu[0].vtx.y < - 100) return false;
53  if (sr->mc.nu[0].vtx.y > 100) return false;
54  if (sr->mc.nu[0].vtx.z < 1000) return false;
55  if (sr->mc.nu[0].vtx.z > 3500) return false;
56  return true;
57  });
58 
59 
60 Cut kBoxCut( [](const caf::SRProxy* sr)
61  {
62  if (sr->mc.nnu == 0) return false;
63  if (sr->mc.nu[0].vtx.x < - 762) return false;
64  if (sr->mc.nu[0].vtx.x > 763) return false;
65  if (sr->mc.nu[0].vtx.y < - 762) return false;
66  if (sr->mc.nu[0].vtx.y > 763) return false;
67  std::bitset<14> binary(sr->hdr.dibmask);
68  std::pair<int,int> planesA = calcFirstLastLivePlane(sr->slc.firstplane,
69  binary);
70  std::pair<int,int> planesB = calcFirstLastLivePlane(sr->slc.lastplane,
71  binary);
72  if ((planesA.first != planesB.first) || (planesA.second != planesB.second))
73  return false;
74  if ((planesA.second - planesA.first + 1)/64 < 4) return false;
75  int first = planesA.first/64;
76  int last = (planesA.second +1)/64;
77  if (sr->mc.nu[0].vtx.z < ((first)*426.0)) return false;
78  if (sr->mc.nu[0].vtx.z > ((last)*426.0)) return false;
79  return true;
80  });
81 
82 const Var kBlock( [](const caf::SRProxy* sr)
83  {
84  int size = -1;
85  std::bitset<14> binary(sr->hdr.dibmask);
86  std::pair<int,int> planesA = calcFirstLastLivePlane(sr->slc.firstplane,
87  binary);
88  std::pair<int,int> planesB = calcFirstLastLivePlane(sr->slc.lastplane,
89  binary);
90  if ((planesA.first != planesB.first) || (planesA.second != planesB.second))
91  return size;
92  size = (planesA.second - planesA.first + 1)/64;
93  return size;
94  });
95 
97 {
98  TH1::SetDefaultSumw2();
99  TH2::SetDefaultSumw2();
100  double kPOTFD = all->Predict(calc).POT();
101  TH1* hSig = all->PredictComponent(calc,
103  Current::kCC,
104  Sign::kBoth).ToTH1(kPOTFD );
105  hSig->SetLineColor(kRed-10);
106 
107  TH1* hSigSel = sel->PredictComponent(calc,
109  Current::kCC,
110  Sign::kBoth).ToTH1(kPOTFD );
111 
112  TH1D* hN = new TH1D("hN_"+title,"Nue CC selection efficiency;diblocks;efficiency",11,3.5,14.5);
113  TH1D* hD = new TH1D("hD_"+title,"Nue CC selection efficiency;diblocks;efficiency",11,3.5,14.5);
114  for(unsigned int i=5;i<16;++i){
115  for(unsigned int j=0;j<hSigSel->GetBinContent(i);++j){
116  hN->Fill(i-1);
117  }
118  for(unsigned int j=0;j<hSig->GetBinContent(i);++j){
119  hD->Fill(i-1);
120  }
121  }
122  hN->Write();
123  hD->Write();
124 
125 }
126 
128 {
129 
130  TH1::SetDefaultSumw2();
131  TH2::SetDefaultSumw2();
132 
133  const Cut kSelectionCut = kNue2017CorePresel;
134 
135  // SpectrumLoader loaderFDSwap(
136  // "prod_caf_R17-03-01-prod3reco.l_fd_genie_fluxswap_fhc_nova_v08_full_v1" );
137  //
138  // NullLoader loaderFDNonSwap;
139  Loaders loader_fdmc;
140  loader_fdmc.SetLoaderPath("prod_caf_R17-03-01-prod3reco.l_fd_genie_fluxswap_fhc_nova_v08_full_v1",
142 
143  const HistAxis axis(
144  "diblocks", Binning::Simple(15,-0.5,14.5), kBlock );
145 
146  // make prediction for all true neutrinos in fiducial volume
147  PredictionNoExtrap* predAll = new PredictionNoExtrap(
148  loader_fdmc,
149  axis, kBoxCut
150  );
151  //prediction for all true neutrinos in fiducial volume that pass preselection
152  PredictionNoExtrap* predSel = new PredictionNoExtrap(
153  loader_fdmc,
154  axis, kSelectionCut&&kBoxCut
155  );
156  // make prediction for all true neutrinos in small box at center of detector
157  PredictionNoExtrap* predIdealAll = new PredictionNoExtrap(
158  loader_fdmc,
159  axis, kBoxCutIdeal
160  );
161  //prediction for all true neutrinos in fiducial volume that pass preselection in small box at center of detector
162  PredictionNoExtrap* predIdealSel = new PredictionNoExtrap(
163  loader_fdmc,
164  axis, kSelectionCut&&kBoxCutIdeal
165  );
166 
167  loader_fdmc.SetSpillCut(kStandardSpillCuts);
168 
169 
170 
171  // Load from files
172  // loaderFDNonSwap.Go();
173  // loaderFDSwap.Go();
174  loader_fdmc.Go();
175 
176  // Open output file
177  TFile file("fd_fiducial_mass_efficiency.root","RECREATE");
178  // Save things in output file
179 
180 
181 
183 
184  //caclulate preselection efficieny as a function of diblock
185  GetExtrap(predAll,predSel, &calc,"FidMass");
186  //caclulate preselection efficieny for small window in ideal detector
187  GetExtrap(predIdealAll,predIdealSel, &calc,"SmallWindow");
188 
189  file.Close();
190 
191 }
Spectrum PredictComponent(osc::IOscCalc *calc, Flavors::Flavors_t flav, Current::Current_t curr, Sign::Sign_t sign) const override
enum BeamMode kRed
Far Detector at Ash River.
Definition: SREnums.h:11
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
(&#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:148
caf::Proxy< caf::SRHeader > hdr
Definition: SRProxy.h:2137
Proxy for caf::StandardRecord.
Definition: SRProxy.h:2126
Collection of SpectrumLoaders for many configurations.
Definition: Loaders.h:23
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
osc::OscCalcDumb calc
void Go()
Call Go() on all the loaders.
Definition: Loaders.cxx:162
void SetSpillCut(const SpillCut &cut)
Definition: Loaders.cxx:121
Charged-current interactions.
Definition: IPrediction.h:39
if(dump)
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
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
caf::Proxy< unsigned int > lastplane
Definition: SRProxy.h:1309
const double j
Definition: BetheBloch.cxx:29
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
void fiducial_accounting()
caf::Proxy< caf::SRTruthBranch > mc
Definition: SRProxy.h:2138
caf::Proxy< caf::SRSlice > slc
Definition: SRProxy.h:2142
void GetExtrap(PredictionNoExtrap *all, PredictionNoExtrap *sel, osc::IOscCalc *calc, TString title)
std::pair< int, int > calcFirstLastLivePlane(int plane, std::bitset< 14 > binary, caf::Det_t det)
const Cut kNue2017CorePresel
Full preselection for 2017 "core" sample with Cut-based CosRej.
Definition: NueCuts2017.h:105
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.
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;})
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:107