fiducial_opt.C
Go to the documentation of this file.
2 #include "CAFAna/Core/HistAxis.h"
3 #include "CAFAna/Core/Binning.h"
4 #include "CAFAna/Core/Var.h"
5 
9 
10 #include "CAFAna/Vars/XsecTunes.h"
13 
15 
16 #include "TFile.h"
17 
22 
23 // using namespace ana::xsec::numubar;
24 
25 const ana::Var vtx_x([](const caf::SRProxy* sr) -> float
26 {
27  if (sr->vtx.elastic.IsValid)
28  return sr->vtx.elastic.vtx.x;
29  return -999.0;
30 });
31 
32 const ana::Var vtx_y([](const caf::SRProxy* sr) -> float
33 {
34  if (sr->vtx.elastic.IsValid)
35  return sr->vtx.elastic.vtx.y;
36  return -999.0;
37 });
38 
39 const ana::Var vtx_z([](const caf::SRProxy* sr) -> float
40 {
41  if (sr->vtx.elastic.IsValid)
42  return sr->vtx.elastic.vtx.z;
43  return -999.0;
44 });
45 
46 void fiducial_opt(string outFilename = "fiducial_opt_hists.root")
47 {
48  // Shared variables
50  ana::HistAxis * vtx_x_axis = new ana::HistAxis("Vtx_X", ana::Binning::Simple(4000, -200, 200), vtx_x);
51  ana::HistAxis * vtx_y_axis = new ana::HistAxis("Vtx_Y", ana::Binning::Simple(4000, -200, 200), vtx_y);
52  ana::HistAxis * vtx_z_axis = new ana::HistAxis("Vtx_Z", ana::Binning::Simple(16000, 0, 1600), vtx_y);
53 
54  // Cuts
56  const ana::Cut signalCut = CutFromNuTruthCut(signalCut_NT);
57  const ana::Cut selectionCut = ana::xsec::numubarcc::kQualityCut && ana::xsec::numubarcc::kContainmentCut; // No vertex applied before fiducial opt
58 
59  // MuonID Optimization object
60  ana::CutOptimization vtx_x_optimization(vtx_x_axis, &signalCut, &selectionCut);
61  ana::CutOptimization vtx_y_optimization(vtx_y_axis, &signalCut, &selectionCut);
62  ana::CutOptimization vtx_z_optimization(vtx_z_axis, &signalCut, &selectionCut);
63 
64  std::map<ana::CutOptimization*, ana::HistAxis*> optimize_axes
65  {
66  {&vtx_x_optimization, vtx_x_axis},
67  {&vtx_y_optimization, vtx_y_axis},
68  {&vtx_z_optimization, vtx_z_axis},
69  };
70 
71  std::map<std::string, ana::CutOptimization*> optimize_labels
72  {
73  {"x", &vtx_x_optimization},
74  {"y", &vtx_y_optimization},
75  {"z", &vtx_z_optimization},
76  };
77 
78  // Systematic Loaders
79  ana::SpectrumLoader * cherenkov_loader = new ana::SpectrumLoader(ana::xsec::numubarcc::syst_datasets["cherenkov"].c_str());
80  ana::SpectrumLoader * lightup_loader = new ana::SpectrumLoader(ana::xsec::numubarcc::syst_datasets["lightup"].c_str());
81  ana::SpectrumLoader * lightdw_loader = new ana::SpectrumLoader(ana::xsec::numubarcc::syst_datasets["lightdw"].c_str());
82  ana::SpectrumLoader * calibup_loader = new ana::SpectrumLoader(ana::xsec::numubarcc::syst_datasets["calibup"].c_str());
83  ana::SpectrumLoader * calibdw_loader = new ana::SpectrumLoader(ana::xsec::numubarcc::syst_datasets["calibdw"].c_str());
84  ana::SpectrumLoader * calibshape_loader = new ana::SpectrumLoader(ana::xsec::numubarcc::syst_datasets["calibshape"].c_str());
85 
86  for(std::pair<std::string, ana::CutOptimization*> optimize_label : optimize_labels)
87  {
88  std::string label = optimize_label.first;
89  ana::CutOptimization * optimization = optimize_label.second;
90  ana::HistAxis* hist_axis = optimize_axes.at(optimization); // Should fail on not-found
91 
92  // Nominal
93  ana::SystematicDef * nominal_syst = new ana::SystematicDef("nominal_" + label, nominal_loader, hist_axis, &ana::xsec::numubarcc::std_wgt);
94  optimization->DefineNominal(nominal_syst);
95 
96  // Shifted samples
97  optimization->DefineSystematic(
98  new ana::SystematicDef("cherenkov_" + label,
99  cherenkov_loader, hist_axis, &ana::xsec::numubarcc::std_wgt));
100  optimization->DefineSystematic(
101  new ana::SystematicDef("lightlevel_" + label,
102  lightup_loader, lightdw_loader, hist_axis, &ana::xsec::numubarcc::std_wgt));
103  optimization->DefineSystematic(
104  new ana::SystematicDef("calibration_" + label,
105  calibup_loader, calibdw_loader, hist_axis, &ana::xsec::numubarcc::std_wgt));
106  optimization->DefineSystematic(
107  new ana::SystematicDef("calibshape_" + label,
108  calibshape_loader, hist_axis, &ana::xsec::numubarcc::std_wgt));
109 
110  // GENIE multiverse
111  ana::SystematicDef * genie_syst = new ana::SystematicDef("genie_" + label, nominal_loader, hist_axis, &ana::xsec::numubarcc::std_wgt);
113  optimization->DefineMultiverseSystematic(genie_syst, genie_univs.GetSystShifts());
114 
115  // PPFX Multiverse
116  ana::SystematicDef * ppfx_syst = new ana::SystematicDef("ppfx_" + label, nominal_loader, hist_axis, &ana::xsec::numubarcc::xsec_wgt);
117  optimization->DefineMultiverseSystematic(ppfx_syst, ana::GetkPPFXFluxUnivWgt());
118  }
119 
120  vtx_x_optimization.Go();
121  vtx_y_optimization.Go();
122  vtx_z_optimization.Go();
123 
124  TFile * fout = TFile::Open(outFilename.c_str(), "RECREATE");
125  vtx_x_optimization.SaveTo(fout, "vtx_x_optimization");
126  vtx_y_optimization.SaveTo(fout, "vtx_y_optimization");
127  vtx_z_optimization.SaveTo(fout, "vtx_z_optimization");
128  fout->Close();
129 }
std::map< const std::string, const std::string > syst_datasets
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 genie_syst(int ndatafiles=INT_MAX, int nmcfiles=INT_MAX, string dataflistpn="/nova/ana/users/slin/sam_definition_file_lists/prod_caf_R16-12-20-prod3recopreview.b_nd_numi_fhc_period3_v1_goodruns.txt", string mcflistpn="/nova/ana/users/slin/sam_definition_file_lists/prod_caf_R16-12-20-prod3recopreview.d_nd_genie_nonswap_fhc_nova_v08_full_v3.txt")
Definition: genie_syst.C:38
void DefineNominal(SystematicDef *nominal)
std::vector< Var > GetkPPFXFluxUnivWgt()
Definition: PPFXWeights.cxx:53
const char * label
caf::Proxy< caf::SRElastic > elastic
Definition: SRProxy.h:2118
std::vector< const ISyst * > getAllXsecSysts_2020()
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
void DefineMultiverseSystematic(SystematicDef *syst, std::vector< Var > mv_weights)
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 std::string nominal_dataset
const ana::Var std_wgt
const ana::Var vtx_y([](const caf::SRProxy *sr) -> float{if(sr->vtx.elastic.IsValid) return sr->vtx.elastic.vtx.y;return-999.0;})
const ana::Var xsec_wgt
caf::Proxy< bool > IsValid
Definition: SRProxy.h:2058
const NuTruthCut kIsNumubarCC_NT([](const caf::SRNeutrinoProxy *nu){return(nu->iscc &&nu->pdg==-14);})
Collaborates with Spectrum and OscillatableSpectrum to fill spectra from CAF files.
caf::Proxy< float > y
Definition: SRProxy.h:107
void DefineSystematic(SystematicDef *syst)
std::vector< SystShifts > GetSystShifts()
caf::Proxy< caf::SRVector3D > vtx
Definition: SRProxy.h:2073
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
void fiducial_opt(string outFilename="fiducial_opt_hists.root")
Definition: fiducial_opt.C:46
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
const ana::Var vtx_z([](const caf::SRProxy *sr) -> float{if(sr->vtx.elastic.IsValid) return sr->vtx.elastic.vtx.z;return-999.0;})
enum BeamMode string