multiverse_efficiency_macro.C
Go to the documentation of this file.
1 /*
2 #include "CAFAna/Core/Binning.h"
3 #include "CAFAna/Cuts/Cuts.h"
4 
5 #include "CAFAna/Core/Spectrum.h"
6 #include "CAFAna/Core/SpectrumLoader.h"
7 
8 #include "CAFAna/Core/Var.h"
9 #include "CAFAna/Core/MultiVar.h"
10 #include "CAFAna/Core/HistAxis.h"
11 #include "CAFAna/Core/SystShifts.h"
12 
13 #include "CAFAna/Vars/Vars.h"
14 #include "3FlavorAna/Vars/NueVarsExtra.h"
15 
16 #include "NDAna/numucc_1Pi/NumuCC1PiVars.h"
17 #include "NDAna/numucc_1Pi/NumuCC1PiCuts.h"
18 
19 #include "G4RwgtPlots/Geant4WeightVars.h"
20 
21 #include "CAFAna/XSec/Multiverse.h"
22 #include "CAFAna/XSec/SystematicDef.h"
23 //#include "CAFAna/XSec/CutOptimization.h"
24 
25 // for POT
26 #include "CAFAna/Analysis/Exposures.h"
27 
28 // for plotting
29 #include "CAFAna/Analysis/Plots.h"
30 
31 // to get datetime in plot names
32 #include "/nova/app/users/csweeney/xsec/code/Utilities/misc.h"
33 
34 #include "StandardRecord/Proxy/SRProxy.h"
35 
36 #include <vector>
37 #include <string>
38 
39 #include "TFile.h"
40 #include "TH1.h"
41 #include "TCanvas.h"
42 #include "TGraphAsymmErrors.h"
43 #include "TString.h"
44 */
45 
46 #include "CAFAna/Core/Binning.h"
47 #include "CAFAna/Cuts/Cuts.h"
48 
49 #include "CAFAna/Core/Spectrum.h"
50 #include "CAFAna/Core/SpectrumLoader.h"
51 
52 #include "CAFAna/Core/Var.h"
53 
54 #include "CAFAna/Vars/Vars.h"
55 
56 #include "CAFAna/XSec/Multiverse.h"
57 #include "CAFAna/XSec/SystematicDef.h"
58 
59 #include "StandardRecord/Proxy/SRProxy.h"
60 
61 #include <vector>
62 #include <string>
63 
64 #include "TFile.h"
65 #include "TH1.h"
66 #include "TCanvas.h"
67 #include "TGraphAsymmErrors.h"
68 #include "TString.h"
69 
70 #include "NDAna/numucc_1Pi/NumuCC1PiVars.h"
71 #include "NDAna/numucc_1Pi/NumuCC1PiCuts.h"
72 
73 #include "G4RwgtPlots/Geant4WeightVars.h"
74 
75 // for POT
76 #include "CAFAna/Analysis/Exposures.h"
77 
78 // for plotting
79 #include "CAFAna/Analysis/Plots.h"
80 
81 // to get datetime in plot names
82 #include "/nova/app/users/csweeney/xsec-TR/Utilities/misc.h"
83 
84 // for kVtxX
85 #include "3FlavorAna/Vars/NueVarsExtra.h"
86 
87 
88 
89 namespace ana
90 {
91 
92  Var kUnity([](const caf::SRProxy *sr){
93  // To be used as "common_weight" for Multiverse constructor
94  return 1.f;
95  });
96 
97 }
98 
99 
100 using namespace ana;
101 
102 const Var kTrueE_nu([](const caf::SRProxy* sr){
103  if(sr->mc.nnu < 1) return -5.f;
104  return float(sr->mc.nu[0].E);
105  });
106 
107 
108 void multiverse_efficiency_macro()
109 {
110 
111 
112  //std::string inName = "/pnfs/nova/scratch/users/csweeney/cached_grid_prod5_g4rwgt_1000/*.caf.root";
113  //std::string inName = "/pnfs/nova/scratch/users/csweeney/grid_prod5_g4rwgt_miniprod/*.caf.root";
114  std::string inName = "/nova/ana/users/csweeney/grid_prod5_g4rwgt_miniprod/*.caf.root";
115  std::string outName = "10_noMuonID_outFile_vtxX.root";
116 
117  SpectrumLoader* loader = new SpectrumLoader(inName);
118 
119  //std::string outName = "test.root";
120  //TFile *outFile = new TFile(outName.c_str(), "RECREATE");
121 
122  const SystShifts * shift = new SystShifts(kNoShift);
123  const Var * common_weight = new Var(kUnity);
124 
125 
126 
127  //~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
128  std::vector<std::vector<Var>> weights_vec{
129  GetkGeantPiplusWeights()//, GetkGeantPiplusWeights(), GetkGeantPiplusWeights(),
130  //GetkGeantPiminusWeights(), GetkGeantPiminusWeights(), GetkGeantPiminusWeights(),
131  //GetkGeantAllPionsWeights(), GetkGeantAllPionsWeights(), GetkGeantAllPionsWeights()
132  };
133 
134  std::vector<std::string> myTitle_vec{
135  "Piplus 100 universes"//, "Piplus 100 universes", "Piplus 100 universes",
136  //"Piminus 100 universes", "Piminus 100 universes", "Piminus 100 universes",
137  //"Total 100 universes", "Total 100 universes", "Total 100 universes"
138  };
139 
140 
141 
142 
143  std::vector<std::string> myLabel_vec{
144  "Vertex in X Direction (cm)"//, "Vertex in Y Direction (cm)", "Vertex in Z Direction (cm)",
145  //"Vertex in X Direction (cm)", "Vertex in Y Direction (cm)", "Vertex in Z Direction (cm)",
146  //"Vertex in X Direction (cm)", "Vertex in Y Direction (cm)", "Vertex in Z Direction (cm)"
147  };
148  std::vector<std::string> myName_vec = {
149  "piplus_vtxX"//, "eff_piplus_vtxY", "eff_piplus_vtxZ",
150  // "eff_piminus_vtxX", "eff_piminus_vtxY", "eff_piminus_vtxZ",
151  //"eff_total_vtxX", "eff_total_vtxY", "eff_total_vtxZ"
152  };
153  std::vector<Binning> myBins_vec{
154  Binning::Simple(10, -200, 200)//, Binning::Simple(10, -200, 200), Binning::Simple(10, 0, 1270),
155  //Binning::Simple(10, -200, 200), Binning::Simple(10, -200, 200), Binning::Simple(10, 0, 1270),
156  //Binning::Simple(10, -200, 200), Binning::Simple(10, -200, 200), Binning::Simple(10, 0, 1270)
157  };
158  std::vector<const Var *> kMyVar_vec{
159  new Var(kVtxX)//, new Var(kVtxY), new Var(kVtxZ),
160  //new Var(kVtxX), new Var(kVtxY), new Var(kVtxZ),
161  //new Var(kVtxX), new Var(kVtxY), new Var(kVtxZ)
162  };
163 
164 
165  /*
166  std::vector<std::string> myLabel_vec{ "True neutrino energy (GeV)" };
167  std::vector<std::string> myName_vec = { "piplus_nuE" };
168  std::vector<Binning> myBins_vec{Binning::Simple(10, 0, 5)};
169  std::vector<const Var *> kMyVar_vec{new Var(kTrueE_nu)};
170  */
171 
172  const Cut * kSignal = new Cut(CutFromNuTruthCut(kIsNumuCC1Pi_NT && kTrueDetector_NT));
173 
174  const Cut * kSelection = new Cut(
175  kFiducialRecoVtx &&
176  kProngContainment &&
177  kMyProngQuality //&&
178  //(kHighestMuonCVN > 0.3)
179  );
180 
181 
182  //-----
183 
184  std::vector<Multiverse *> mySig_MV_vec;
185  std::vector<Spectrum *> mySig_Spec_vec;
186 
187  std::vector<Multiverse *> mySelsig_MV_vec;
188  std::vector<Spectrum *> mySelsig_Spec_vec;
189 
190 
191 
192  for(uint i=0; i<myName_vec.size(); ++i){
193 
194  std::string myLabel = myLabel_vec[i];
195  Binning myBins = myBins_vec[i];
196  const Var * kMyVar = kMyVar_vec[i];
197  std::vector<Var> weights = weights_vec[i];
198 
199  // std::string myName = myName_vec[i];
200 
201 
202 
203  const HistAxis * axis = new HistAxis(myLabel,
204  myBins,
205  *kMyVar);
206 
207 
208  Multiverse * mySig_MV = new Multiverse(*loader,
209  *axis,
210  *kSignal,
211  *shift,
212  weights,
213  *common_weight);
214 
215 
216  Spectrum * sSig = new Spectrum(myLabel, myBins,
217  *loader, *kMyVar,
218  *kSignal);
219  mySig_MV_vec.push_back(mySig_MV);
220  mySig_Spec_vec.push_back(sSig);
221 
222  //~~~~~~~
223 
224  Multiverse * mySelsig_MV = new Multiverse(*loader,
225  *axis,
226  *kSignal && *kSelection,
227  *shift,
228  weights,
229  *common_weight);
230 
231 
232  Spectrum * sSelsig = new Spectrum(myLabel, myBins,
233  *loader, *kMyVar,
234  *kSignal && *kSelection);
235  mySelsig_MV_vec.push_back(mySelsig_MV);
236  mySelsig_Spec_vec.push_back(sSelsig);
237  }
238 
239  //~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
240 
241 
242 
243 
244  loader->Go();
245 
246 
247  TFile * outFile = new TFile( outName.c_str(), "recreate");
248 
249  for(uint i=0; i< myName_vec.size(); ++i){
250 
251  TDirectory* base_dir = outFile->mkdir(myName_vec[i].c_str());
252 
253  mySig_Spec_vec[i]->SaveTo( base_dir, "spec_Sig" );
254  mySig_MV_vec[i]->SaveTo(base_dir, "mv_Sig");
255 
256  mySelsig_Spec_vec[i]->SaveTo( base_dir, "spec_Selsig" );
257  mySelsig_MV_vec[i]->SaveTo(base_dir, "mv_Selsig");
258  }
259 
260 
261 }
262 
ana::kMyProngQuality
const Cut kMyProngQuality([](const caf::SRProxy *sr){if(sr->vtx.elastic.IsValid==false) return false;int nPngs=sr->vtx.elastic.fuzzyk.npng;bool isQuality=(nPngs > 0 && sr->slc.nhit > 20 && sr->slc.ncontplanes > 4);return isQuality;})
ana::Binning
Represent the binning of a Spectrum&#39;s x-axis.
Definition: Binning.h:16
ana
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
ana::SystShifts
Simple record of shifts applied to systematic parameters.
Definition: SystShifts.h:20
ana::kProngContainment
const Cut kProngContainment([](const caf::SRProxy *sr) {if(sr->vtx.elastic.IsValid==false) return false;int iMuonPng=kBestMuonPng(sr); for(unsigned int i=0;i< sr->vtx.elastic.fuzzyk.nshwlid;++i){const TVector3 start=sr->vtx.elastic.fuzzyk.png[i].shwlid.start;const TVector3 stop=sr->vtx.elastic.fuzzyk.png[i].shwlid.stop;if(std::min(start.X(), stop.X())< -180.0) return false;if(std::max(start.X(), stop.X()) > 180.0) return false;if(std::min(start.Y(), stop.Y())< -180.0) return false;if(std::max(start.Y(), stop.Y()) > 180.0) return false;if(std::min(start.Z(), stop.Z())< 20.0) return false;if(std::max(start.Z(), stop.Z()) > 1565.0) return false;}int nPngs=sr->vtx.elastic.fuzzyk.npng;if(nPngs< 1) return false;for(int i=0;i< nPngs;++i){if(i==iMuonPng) continue; else if(sr->vtx.elastic.fuzzyk.png[i].shwlid.start.Z() > 1275|| sr->vtx.elastic.fuzzyk.png[i].shwlid.stop.Z() > 1275) return false;} if(iMuonPng< 0||iMuonPng > nPngs) return false;auto &muon_prong=sr->vtx.elastic.fuzzyk.png[iMuonPng]; if(muon_prong.bpf.muon.IsValid==false) return false;auto &muon_bpf=muon_prong.bpf.muon;bool ret=((muon_bpf.stop.z< 1275||muon_bpf.trkyposattrans< 55) &&muon_bpf.trkfwdcellnd > 5 &&muon_bpf.trkbakcellnd > 10);return ret;})
caf::Proxy< caf::StandardRecord >
Proxy for caf::StandardRecord.
Definition: SRProxy.h:2126
caf::Proxy< caf::SRTruthBranch >::nu
caf::Proxy< std::vector< caf::SRNeutrino > > nu
Definition: SRProxy.h:618
cmf::Spectrum
std::vector< double > Spectrum
Definition: Constants.h:746
caf::Proxy< caf::SRTruthBranch >::nnu
caf::Proxy< short int > nnu
Definition: SRProxy.h:617
ana::kVtxX
const Var kVtxX
Definition: NueVarsExtra.cxx:31
ana::Spectrum
Representation of a spectrum in any variable, with associated POT.
Definition: Spectrum.h:39
multiverse_efficiency_macro
void multiverse_efficiency_macro()
Definition: multiverse_efficiency_macro.C:108
ana::kTrueDetector_NT
const NuTruthCut kTrueDetector_NT([](const caf::SRNeutrinoProxy *nu){bool fid=(nu->vtx.X()< 192 && nu->vtx.X() >-191 && nu->vtx.Y()< 194 && nu->vtx.Y() >-187 && nu->vtx.Z()< 1270 && nu->vtx.Z() > 0);return fid; })
outFile
TFile * outFile
Definition: PlotXSec.C:135
ana::Var
_Var< caf::SRProxy > Var
Definition: Var.h:7
weights
Var weights
Definition: getFlatWeightSpectra.C:30
sr
caf::StandardRecord * sr
Definition: FillTrainingTrees.h:18
ana::SpectrumLoader::Go
virtual void Go() override
Load all the registered spectra.
Definition: SpectrumLoader.cxx:66
demo0.loader
loader
Definition: demo0.py:10
ana::kIsNumuCC1Pi_NT
const NuTruthCut kIsNumuCC1Pi_NT([](const caf::SRNeutrinoProxy *nu){int nPions=0;int nMuons=0;int nPrims=nu->prim.size();for(int prim_idx=0;prim_idx< nPrims;prim_idx++){auto &prim=nu->prim[prim_idx];int pdg=prim.pdg;if(abs(pdg)==211)++nPions;else if(pdg==13)++nMuons;}bool is_numucc=nu->iscc &&nu->pdg==14;bool ret=is_numucc &&(nPions==1)&&(nMuons==1);return ret;})
ana::kNoShift
const SystShifts kNoShift
Definition: SystShifts.cxx:22
caf::Proxy< caf::StandardRecord >::mc
caf::Proxy< caf::SRTruthBranch > mc
Definition: SRProxy.h:2138
ana::kSignal
const Cut kSignal
Definition: SINCpi0_Cuts.h:383
ana::_HistAxis< Var >
ana::kFiducialRecoVtx
const Cut kFiducialRecoVtx([](const caf::SRProxy *sr){if(sr->vtx.elastic.IsValid==false) return false;bool isFid=(sr->vtx.elastic.vtx.x< my_vtxmax.X()&& sr->vtx.elastic.vtx.x > my_vtxmin.X()&& sr->vtx.elastic.vtx.y< my_vtxmax.Y()&& sr->vtx.elastic.vtx.y > my_vtxmin.Y()&& sr->vtx.elastic.vtx.z< my_vtxmax.Z()&& sr->vtx.elastic.vtx.z > my_vtxmin.Z());return isFid;})
ana::SpectrumLoader
Collaborates with Spectrum and OscillatableSpectrum to fill spectra from CAF files.
Definition: SpectrumLoader.h:21
ana::Cut
_Cut< caf::SRProxy, caf::SRSpillProxy > Cut
Definition: Cut.h:9
ana::kTrueE_nu
const Var kTrueE_nu([](const caf::SRProxy *sr){if(sr->mc.nnu< 1) return-5.f;return float(sr->mc.nu[0].E);})
ana::_Cut< caf::SRProxy, caf::SRSpillProxy >
ana::_Var< caf::SRProxy >
ana::kUnity
Var kUnity([](const caf::SRProxy *sr){return 1.f;})
ana::Binning::Simple
static Binning Simple(int n, double lo, double hi, const std::vector< std::string > &labels={})
Definition: Binning.cxx:107
ana::HistAxis
_HistAxis< Var > HistAxis
Definition: HistAxis.h:10
ana::CutFromNuTruthCut
Cut CutFromNuTruthCut(const NuTruthCut &stc)
Definition: Cut.cxx:7
string
enum BeamMode string
uint
unsigned int uint
Definition: DAQChannelMapBaseOLD.h:22
Generated on Thu May 6 2021 14:10:20 for NOvA by   doxygen 1.8.11