vertex_fiducial_optimization.C
Go to the documentation of this file.
1 // Author: Prabhjot Singh (prabhjot@fnal.gov)
2 // Date: 21 Feb 2021
3 // Purpose: Macro for fiducial optimization
4 // - This macro is inspired from NDAna/numubarcc_inc/Analysis/FiducialOpt/fiducial_opt.C
5 
6 // CAFAna includes
8 #include "CAFAna/Core/HistAxis.h"
9 #include "CAFAna/Core/Binning.h"
10 #include "CAFAna/Core/Var.h"
11 
15 
16 #include "CAFAna/Vars/XsecTunes.h"
19 
21 
22 // NDAna includes
27 
28 // ROOT includes
29 #include "TFile.h"
30 
31 const ana::Var vtx_x([](const caf::SRProxy* sr) -> float
32 {
33  if (sr->vtx.elastic.IsValid)
34  return sr->vtx.elastic.vtx.x;
35  return -999.0;
36 });
37 
38 const ana::Var vtx_y([](const caf::SRProxy* sr) -> float
39 {
40  if (sr->vtx.elastic.IsValid)
41  return sr->vtx.elastic.vtx.y;
42  return -999.0;
43 });
44 
45 const ana::Var vtx_z([](const caf::SRProxy* sr) -> float
46 {
47  if (sr->vtx.elastic.IsValid)
48  return sr->vtx.elastic.vtx.z;
49  return -999.0;
50 });
51 
52 void vertex_fiducial_optimization(string outFilename = "vertex_fiducial_optimization_hists.root")
53 {
54  // Shared variables
55  // Prod 5.1. Flatsumdecaf RHC ND MC dataset. This dataset has 400 files
56  ana::SpectrumLoader * nominal_loader = new ana::SpectrumLoader("prod_flatsumdecaf_development_nd_genie_N1810j0211a_nonswap_rhc_nova_v08_full_ndphysics_contain_v1");
57  ana::HistAxis * vtx_x_axis = new ana::HistAxis("Vtx_X", ana::Binning::Simple(4000, -200, 200), vtx_x);
58  ana::HistAxis * vtx_y_axis = new ana::HistAxis("Vtx_Y", ana::Binning::Simple(4000, -200, 200), vtx_y);
59  ana::HistAxis * vtx_z_axis = new ana::HistAxis("Vtx_Z", ana::Binning::Simple(16000, 0, 1600), vtx_y);
60 
61  // Cuts
63  const ana::Cut signalCut = CutFromNuTruthCut(signalCut_NT);
65 
66  // Vertex Optimization object
67  ana::CutOptimization vtx_x_optimization(vtx_x_axis, &signalCut, &selectionCut);
68  ana::CutOptimization vtx_y_optimization(vtx_y_axis, &signalCut, &selectionCut);
69  ana::CutOptimization vtx_z_optimization(vtx_z_axis, &signalCut, &selectionCut);
70 
71  std::map<ana::CutOptimization*, ana::HistAxis*> optimize_axes
72  {
73  {&vtx_x_optimization, vtx_x_axis},
74  {&vtx_y_optimization, vtx_y_axis},
75  {&vtx_z_optimization, vtx_z_axis},
76  };
77 
78  std::map<std::string, ana::CutOptimization*> optimize_labels
79  {
80  {"x", &vtx_x_optimization},
81  {"y", &vtx_y_optimization},
82  {"z", &vtx_z_optimization},
83  };
84 
85  // Systematic Loaders
86  // None at the moment
87 
88  for(std::pair<std::string, ana::CutOptimization*> optimize_label : optimize_labels)
89  {
90  std::string label = optimize_label.first;
91  ana::CutOptimization * optimization = optimize_label.second;
92  ana::HistAxis* hist_axis = optimize_axes.at(optimization);
93 
94  // Nominal
95  ana::SystematicDef * nominal_syst = new ana::SystematicDef("nominal_" + label, nominal_loader, hist_axis, &ana::xsec::numubarcc::std_wgt);
96  optimization->DefineNominal(nominal_syst);
97 
98  // Shifted samples
99  // None at the moment
100 
101  // GENIE multiverse
102  // None at the moment
103 
104  // PPFX Multiverse
105  // None at the moment
106  }
107 
108  vtx_x_optimization.Go();
109  vtx_y_optimization.Go();
110  vtx_z_optimization.Go();
111 
112  TFile * fout = TFile::Open(outFilename.c_str(), "RECREATE");
113  vtx_x_optimization.SaveTo(fout, "vtx_x_optimization");
114  vtx_y_optimization.SaveTo(fout, "vtx_y_optimization");
115  vtx_z_optimization.SaveTo(fout, "vtx_z_optimization");
116  fout->Close();
117 }
void vertex_fiducial_optimization(string outFilename="vertex_fiducial_optimization_hists.root")
Proxy for caf::StandardRecord.
Definition: SRProxy.h:2126
const ana::Var vtx_x([](const caf::SRProxy *sr) -> float{if(sr->vtx.elastic.IsValid) return sr->vtx.elastic.vtx.x;return-999.0;})
void DefineNominal(SystematicDef *nominal)
const char * label
caf::Proxy< caf::SRElastic > elastic
Definition: SRProxy.h:2118
const Cut kQualityCut([](const caf::SRProxy *sr){return(sr->trk.kalman.ntracks > 0 &&sr->slc.nhit > 20 &&sr->slc.ncontplanes > 4);})
Preselection.
caf::Proxy< float > z
Definition: SRProxy.h:108
caf::Proxy< float > x
Definition: SRProxy.h:106
void SaveTo(TDirectory *dir, const std::string &name) const
const Cut kContainmentCut([](const caf::SRProxy *sr){if(!sr->vtx.elastic.IsValid) return false;int ibesttrk=kBestMuonIDIndex(sr);if(ibesttrk< 0||ibesttrk >=int(sr->trk.kalman.ntracks)) return false;for(const caf::SRFuzzyKProngProxy &prong:sr->vtx.elastic.fuzzyk.png) if(!VtxInBounds(&prong.shwlid.start, containLow, containHigh)||!VtxInBounds(&prong.shwlid.stop, containLow, containHigh)) return false;if(sr->trk.kalman.ntracks< 1) return false;const unsigned short muon_catcher_edge=1275;for(unsigned int i=0;i< sr->trk.kalman.ntracks;++i){if(int(i)==ibesttrk) continue;if(sr->trk.kalman.tracks[i].start.Z() > muon_catcher_edge|| sr->trk.kalman.tracks[i].stop.Z() > muon_catcher_edge) return false;}const caf::SRKalmanTrackProxy &besttrack=sr->trk.kalman.tracks[ibesttrk];return((besttrack.stop.Z()< muon_catcher_edge||besttrack.trkyposattrans< 55) &&besttrack.trkfwdcellnd > 5 &&besttrack.trkbakcellnd > 10);})
GenericSystematicDef< Cut, Var > SystematicDef
caf::StandardRecord * sr
const ana::Var std_wgt
caf::Proxy< bool > IsValid
Definition: SRProxy.h:2058
const NuTruthCut kIsNumubarCC_NT([](const caf::SRNeutrinoProxy *nu){return(nu->iscc &&nu->pdg==-14);})
const ana::Var vtx_y([](const caf::SRProxy *sr) -> float{if(sr->vtx.elastic.IsValid) return sr->vtx.elastic.vtx.y;return-999.0;})
Collaborates with Spectrum and OscillatableSpectrum to fill spectra from CAF files.
caf::Proxy< float > y
Definition: SRProxy.h:107
caf::Proxy< caf::SRVector3D > vtx
Definition: SRProxy.h:2073
const ana::Var vtx_z([](const caf::SRProxy *sr) -> float{if(sr->vtx.elastic.IsValid) return sr->vtx.elastic.vtx.z;return-999.0;})
const NuTruthCut kTrueVtxLooseCut_NT([](const caf::SRNeutrinoProxy *nu){return VtxInBounds(&nu->vtx, loose_vtx_min, loose_vtx_max);})
Template for Cuts applied to any type of object.
Definition: Cut.h:15
caf::Proxy< caf::SRVertexBranch > vtx
Definition: SRProxy.h:2146
static Binning Simple(int n, double lo, double hi, const std::vector< std::string > &labels={})
Definition: Binning.cxx:107
_HistAxis< Var > HistAxis
Definition: HistAxis.h:10
Cut CutFromNuTruthCut(const NuTruthCut &stc)
Definition: Cut.cxx:7
enum BeamMode string