demo_cut_optimization.C
Go to the documentation of this file.
1 #pragma once
2 
4 
5 // analysis cuts/vars
9 
10 // for data pot
12 
13 // Multiverse includes
15 #include "CAFAna/Vars/XsecTunes.h"
19 
20 using namespace ana;
21 
22 
23 void demo_cut_optimization(std::string input_file_name = "")
24 {
25  if(input_file_name == "" ) {
26  std::map<TString, SpectrumLoader *> loaders;
27  loaders["kNominal"] = new SpectrumLoader("prod_caf_R17-11-14-prod4reco.remid-hotfix.b_nd_genie_nonswap_fhc_nova_v08_full_v1");
28  loaders["kCalibNeg"] = new SpectrumLoader("prod_caf_R17-11-14-prod4reco.remid-hotfix.b_nd_genie_nonswap_fhc_nova_v08_full_calib-shift-nd-xyview-neg-offset_v1_good_seventh");
29  loaders["kCalibPos"] = new SpectrumLoader("prod_caf_R17-11-14-prod4reco.remid-hotfix.b_nd_genie_nonswap_fhc_nova_v08_full_calib-shift-nd-xyview-pos-offset_v1_good_seventh");
30  loaders["kCalibShape"] = new SpectrumLoader("prod_caf_R17-11-14-prod4reco.remid-hotfix.b_nd_genie_nonswap_fhc_nova_v08_full_calib-shift-nd-func_v1_good_seventh");
31  loaders["kLightUpCalibDown"] = new SpectrumLoader("prod_caf_R17-11-14-prod4reco.remid-hotfix.b_nd_genie_nonswap_fhc_nova_v08_full_lightmodel-lightup-calibdown_v1_good_seventh");
32  loaders["kLightDownCalibUp"] = new SpectrumLoader("prod_caf_R17-11-14-prod4reco.remid-hotfix.b_nd_genie_nonswap_fhc_nova_v08_full_lightmodel-lightdown-calibup_v1_good_seventh");
33  loaders["kCherenkov"] = new SpectrumLoader("prod_caf_R17-11-14-prod4reco.remid-hotfix.b_nd_genie_nonswap_fhc_nova_v08_full_ckv-proton-shift-down_v1_good_seventh");
34 
35  // define a hist axis to optimize
36  const HistAxis * opti_axis = new HistAxis("Vertex in X Direction (cm)",
37  Binning::Simple(50, -200, 200),
38  kVtxX);
39  // Define selection and signal cuts
41  // This demo is for fiducial cut optimization so the selection
42  // should be all preselection minus fiducial constraints
43  const Cut * kSelection = new Cut(kcDQ &&
45  kcFrontPlanes &&
46  kcNHits);
47 
48  // central value weight to be applied to all samples
49  const Var * kCVWeight = new Var(kPPFXFluxCVWgt * kXSecCVWgt2018);
50 
51  // nominal definition
52  SystematicDef * sNominal =
53  new SystematicDef("sNominal",
54  loaders["kNominal"],
55  opti_axis,
56  kCVWeight);
57 
58  // systematic definitions
59  SystematicDef * sLightlevelUpDown =
60  new SystematicDef("sLightlevelUpDown",
61  loaders["kLightUpCalibDown"],
62  loaders["kLightDownCalibUp"],
63  opti_axis,
64  kCVWeight);
65  SystematicDef * sCalibPosNeg =
66  new SystematicDef("sCalibPosNeg",
67  loaders["kCalibPos"],
68  loaders["kCalibNeg"],
69  opti_axis,
70  kCVWeight);
71  SystematicDef * sCalibShape =
72  new SystematicDef("sCalibShape",
73  loaders["kCalibShape"],
74  opti_axis,
75  kCVWeight);
76  SystematicDef * sCherenkov =
77  new SystematicDef("sCherenkov",
78  loaders["kCherenkov"],
79  opti_axis,
80  kCVWeight);
81 
82  // Setup the genie multiverses
83  // We take flux uncertainty to be constant, so no need for PPFX universe
84  // for the sake of example, just do 10 universes
85  GenieMultiverseParameters genie_multiverse(10, getAllXsecSysts_2018());
86  std::vector<SystShifts> genie_shifts = genie_multiverse.GetSystShifts();
87 
88  SystematicDef * sGenieMultiverse =
89  new SystematicDef("sGenieMultiverse",
90  loaders["kNominal"],
91  opti_axis,
92  kCVWeight);
93 
94  // create the optimizer.
95  CutOptimization * optimize = new CutOptimization(opti_axis,
96  kSignal,
97  kSelection);
98  optimize->DefineNominal(sNominal);
99  optimize->DefineSystematic(sLightlevelUpDown);
100  optimize->DefineSystematic(sCalibPosNeg);
101  optimize->DefineSystematic(sCalibShape);
102  optimize->DefineSystematic(sCherenkov);
103  optimize->DefineMultiverseSystematic(sGenieMultiverse, genie_shifts);
104 
105  // optionally set spill cuts for all loaders
106  optimize->SetSpillCuts(kStandardSpillCuts);
107 
108  // Let the loaders go
109  optimize->Go();
110 
111  // save to file
112  TFile * output = new TFile("demo_cut_optimization.root", "recreate");
113  optimize->SaveTo(output->mkdir("vtx_cut_optimize"));
114  output->Close();
115  }
116  else {
117  // load input file provided in argument
118  TFile * input = TFile::Open(input_file_name.c_str(), "read");
119  CutOptimization * optimize = CutOptimization::LoadFrom(input->GetDirectory("vtx_cut_optimize")).release();
120  input->Close();
121 
122 
123  // save resulting histograms back to file.
124  TFile * results = new TFile("demo_cut_optimization_results.root", "recreate");
125  results->mkdir("vtx_cut_optimize/results");
126  TDirectory * save_dir = results->GetDirectory("vtx_cut_optimize/results");
127 
128  // directory to save plots
129  std::string plot_dump = "demo_cut_optimization_plots";
130 
131  optimize->Optimize(CutOptimization::kBinByBin, // OptimizationMethod_t
133  kAna2018FHCPOT, // scale to data pot
134  save_dir,
135  plot_dump,
136  true);
137  delete optimize;
138  }
139 
140 }
const NuTruthCut kIsNueCCST([](const caf::SRNeutrinoProxy *truth){return(truth->iscc &&truth->pdg==12);})
Definition: NumuCCIncCuts.h:40
const Cut kcNueCCIncContainment([](const caf::SRProxy *sr){if(sr->vtx.elastic[0].fuzzyk.nshwlid< 1) return false;TVector3 start=sr->vtx.elastic[0].fuzzyk.png[0].shwlid.start;TVector3 stop=sr->vtx.elastic[0].fuzzyk.png[0].shwlid.stop;for(uint ix=0;ix< sr->vtx.elastic[0].fuzzyk.nshwlid;ix++){TVector3 start_muoncatcher=sr->vtx.elastic[0].fuzzyk.png[ix].shwlid.start;TVector3 stop_muoncatcher=sr->vtx.elastic[0].fuzzyk.png[ix].shwlid.stop;if(std::max(start_muoncatcher.Z(), stop_muoncatcher.Z()) > 1250.0) return false;}if(sr->sel.nuecosrej.distallpngtop< 50) return false;if(sr->sel.nuecosrej.distallpngbottom< 30) return false;if(sr->sel.nuecosrej.distallpngeast< 50) return false;if(sr->sel.nuecosrej.distallpngwest< 30) return false;if(sr->sel.nuecosrej.distallpngfront< 150) return false;return true;})
ofstream output
const Cut kcFrontPlanes
Definition: NueCCIncCuts.h:117
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
const Cut kcNHits
Definition: NueCCIncCuts.h:118
std::vector< const ISyst * > getAllXsecSysts_2018()
Get master XSec syst list for 2018 analyses.
const Var kPPFXFluxCVWgt
Definition: PPFXWeights.h:16
GenericVar< caf::SRProxy > Var
Representation of a variable to be retrieved from a caf::StandardRecord object.
Definition: Var.h:76
const NuTruthCut kTrueDetectorST([](const caf::SRNeutrinoProxy *sr){return(sr->vtx.X()< tvtxmax.X()&& sr->vtx.X() > tvtxmin.X()&& sr->vtx.Y() > tvtxmin.Y()&& sr->vtx.Y()< tvtxmax.Y()&& sr->vtx.Z() > tvtxmin.Z()&& sr->vtx.Z()< tvtxmax.Z());})
Definition: NumuCCIncCuts.h:32
static std::unique_ptr< CutOptimization > LoadFrom(TDirectory *dir, const std::string &name)
GenericCut< caf::SRProxy > Cut
Representation of a cut (selection) to be applied to a caf::StandardRecord object.
Definition: Cut.h:96
GenericSystematicDef< caf::SRProxy > SystematicDef
const Var kVtxX
void DefineNominal(SystematicDef *nominal)
const Cut kcDQ
Definition: NueCCIncCuts.h:114
void SaveTo(TDirectory *dir, const std::string &name) const
void DefineMultiverseSystematic(SystematicDef *syst, std::vector< Var > mv_weights)
void demo_cut_optimization(std::string input_file_name="")
void SetSpillCuts(const SpillCut &)
const Cut kSelection
Definition: SINCpi0_Cuts.h:328
const Cut kSignal
Definition: SINCpi0_Cuts.h:325
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
Collaborates with Spectrum and OscillatableSpectrum to fill spectra from CAF files.
void DefineSystematic(SystematicDef *syst)
const Var kXSecCVWgt2018
Definition: XsecTunes.h:49
std::vector< Loaders * > loaders
Definition: syst_header.h:386
std::vector< SystShifts > GetSystShifts()
GenericHistAxis< Var > HistAxis
Definition: HistAxis.h:111
const SpillCut kStandardSpillCuts
Apply this unless you&#39;re doing something special.
Definition: SpillCuts.h:46
const double kAna2018FHCPOT
Definition: Exposures.h:207
static Binning Simple(int n, double lo, double hi, const std::vector< std::string > &labels={})
Definition: Binning.cxx:114
Cut CutFromNuTruthCut(const NuTruthCut &stc)
Definition: Cut.cxx:134