muonid_optimization.C
Go to the documentation of this file.
1 #pragma once
2 
4 
5 // analysis cuts/vars
11 
12 // for data pot
14 
15 // Multiverse includes
17 #include "CAFAna/Vars/XsecTunes.h"
21 
22 
23 
24 using namespace ana;
25 
27 {
28 
29  std::map<TString, SpectrumLoader *> loaders;
30  loaders["kNominal"] = new SpectrumLoader(nuebarccinc::PROD5_MC_RHC_NOMINAL );
31  loaders["kCalibNeg"] = new SpectrumLoader(nuebarccinc::PROD5_MC_RHC_CALIB_DOWN );
32  loaders["kCalibPos"] = new SpectrumLoader(nuebarccinc::PROD5_MC_RHC_CALIB_UP );
33  loaders["kCalibShape"] = new SpectrumLoader(nuebarccinc::PROD5_MC_RHC_CALIB_SHAPE);
34  loaders["kLightUpCalibDown"] = new SpectrumLoader(nuebarccinc::PROD5_MC_RHC_LIGHT_UP );
35  loaders["kLightDownCalibUp"] = new SpectrumLoader(nuebarccinc::PROD5_MC_RHC_LIGHT_DOWN );
36  loaders["kCherenkov"] = new SpectrumLoader(nuebarccinc::PROD5_MC_RHC_CHERENKOV );
37 
38  // define a hist axis to optimize
39  const HistAxis * opti_axis = new HistAxis("MuonID",
40  Binning::Simple(50, -1, 1),
42  // Define selection and signal cuts
46  // This demo is for fiducial cut optimization so the selection
47  // should be all preselection minus fiducial constraints
48  const Cut * kSelection = new Cut(nuebarccinc::kDQ &&
54 
55  // central value weight to be applied to all samples
56  const Var * kCVWeight = new Var(kPPFXFluxCVWgt * kXSecCVWgt2020);
57 
58  // nominal definition
59  SystematicDef * sNominal =
60  new SystematicDef("sNominal",
61  loaders["kNominal"],
62  opti_axis,
63  kCVWeight);
64 
65  // systematic definitions
66  SystematicDef * sLightlevelUpDown =
67  new SystematicDef("sLightlevelUpDown",
68  loaders["kLightUpCalibDown"],
69  loaders["kLightDownCalibUp"],
70  opti_axis,
71  kCVWeight);
72  SystematicDef * sCalibPosNeg =
73  new SystematicDef("sCalibPosNeg",
74  loaders["kCalibPos"],
75  loaders["kCalibNeg"],
76  opti_axis,
77  kCVWeight);
78  SystematicDef * sCalibShape =
79  new SystematicDef("sCalibShape",
80  loaders["kCalibShape"],
81  opti_axis,
82  kCVWeight);
83  SystematicDef * sCherenkov =
84  new SystematicDef("sCherenkov",
85  loaders["kCherenkov"],
86  opti_axis,
87  kCVWeight);
88 
89  // Setup the genie multiverses
90  // We take flux uncertainty to be constant, so no need for PPFX universe
91  GenieMultiverseParameters genie_multiverse(300, getAllXsecSysts_2020());
92  std::vector<SystShifts> genie_shifts = genie_multiverse.GetSystShifts();
93 
94  SystematicDef * sGenieMultiverse =
95  new SystematicDef("sGenieMultiverse",
96  loaders["kNominal"],
97  opti_axis,
98  kCVWeight);
99 
100  // create the optimizer.
101  CutOptimization * optimize = new CutOptimization(opti_axis,
102  kSignal,
103  kSelection);
104 
105  optimize->DefineNominal(sNominal);
106  optimize->DefineSystematic(sLightlevelUpDown);
107  optimize->DefineSystematic(sCalibPosNeg);
108  optimize->DefineSystematic(sCalibShape);
109  optimize->DefineSystematic(sCherenkov);
110  optimize->DefineMultiverseSystematic(sGenieMultiverse, genie_shifts);
111 
112  // optionally set spill cuts for all loaders
113  optimize->SetSpillCuts(kStandardSpillCuts);
114 
115  // Let the loaders go
116  optimize->Go();
117 
118  // save to file
119  TFile * output = new TFile("muonid_optimization.root", "recreate");
120  optimize->SaveTo(output, "muonid_cut_optimize");
121  output->Close();
122 
123 }
ofstream output
const ana::Cut kCVNNueIDCut
const ana::NuTruthCut kNuebarCCST([](const caf::SRNeutrinoProxy *truth){return(truth->iscc && truth->pdg==-12 && truth->pdgorig==-12);})
_HistAxis< Var > HistAxis
Definition: HistAxis.h:103
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
const ana::Var kNHits([](const caf::SRProxy *sr){return sr->slc.nhit;})
const Var kPPFXFluxCVWgt
Definition: PPFXWeights.h:16
void muonid_optimization()
const ana::Cut kFrontPlanes
const std::string PROD5_MC_RHC_LIGHT_UP
const ana::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());})
const ana::NuTruthCut kNueCCST([](const caf::SRNeutrinoProxy *truth){return(truth->iscc && truth->pdg==12 && truth->pdgorig==12);})
const std::string PROD5_MC_RHC_CALIB_DOWN
GenericSystematicDef< caf::SRProxy > SystematicDef
const std::string PROD5_MC_RHC_LIGHT_DOWN
void DefineNominal(SystematicDef *nominal)
_Var< caf::SRProxy > Var
Representation of a variable to be retrieved from a caf::StandardRecord object.
Definition: Var.h:74
std::vector< const ISyst * > getAllXsecSysts_2020()
const ana::Cut kNuebarCCIncContainmentLoose([](const caf::SRProxy *sr){if(!sr->vtx.elastic.IsValid) return false;if(sr->vtx.elastic.fuzzyk.nshwlid< 1) return false;TVector3 start=sr->vtx.elastic.fuzzyk.png[0].shwlid.start;TVector3 stop=sr->vtx.elastic.fuzzyk.png[0].shwlid.stop;for(uint ix=0;ix< sr->vtx.elastic.fuzzyk.nshwlid;ix++){TVector3 start_muoncatcher=sr->vtx.elastic.fuzzyk.png[ix].shwlid.start;TVector3 stop_muoncatcher=sr->vtx.elastic.fuzzyk.png[ix].shwlid.stop;if(std::max(start_muoncatcher.Z(), stop_muoncatcher.Z()) > 1250.0) return false;}if(sr->sel.nuecosrej.distallpngtop< 30) return false;if(sr->sel.nuecosrej.distallpngbottom< 10) return false;if(sr->sel.nuecosrej.distallpngeast< 30) return false;if(sr->sel.nuecosrej.distallpngwest< 10) return false;if(sr->sel.nuecosrej.distallpngfront< 100) return false;return true;})
const std::string PROD5_MC_RHC_NOMINAL
void SaveTo(TDirectory *dir, const std::string &name) const
void DefineMultiverseSystematic(SystematicDef *syst, std::vector< Var > mv_weights)
const std::string PROD5_MC_RHC_CALIB_UP
const std::string PROD5_MC_RHC_CALIB_SHAPE
_Cut< caf::SRProxy > Cut
Representation of a cut (selection) to be applied to a caf::StandardRecord object.
Definition: Cut.h:96
void SetSpillCuts(const SpillCut &)
const std::string PROD5_MC_RHC_CHERENKOV
const Cut kSelection
Definition: SINCpi0_Cuts.h:404
const Cut kSignal
Definition: SINCpi0_Cuts.h:400
Collaborates with Spectrum and OscillatableSpectrum to fill spectra from CAF files.
void DefineSystematic(SystematicDef *syst)
std::vector< Loaders * > loaders
Definition: syst_header.h:386
std::vector< SystShifts > GetSystShifts()
const ana::Cut kDQ([](const caf::SRProxy *sr){if(sr->sel.nuecosrej.hitsperplane >=8) return false;if(!sr->vtx.elastic.IsValid) return false;if(sr->vtx.elastic.fuzzyk.npng==0) return false;if(sr->vtx.elastic.fuzzyk.nshwlid==0) return false;if(sr->vtx.elastic.fuzzyk.png[0].shwlid.nhitx<=5||sr->vtx.elastic.fuzzyk.png[0].shwlid.nhity<=5) return false;if(sr->vtx.elastic.fuzzyk.png[0].shwlid.gap > 100) return false; if(sr->vtx.elastic.fuzzyk.nshwlid >=2){if((sr->vtx.elastic.fuzzyk.png[0].shwlid.dir). Dot(sr->vtx.elastic.fuzzyk.png[1].shwlid.dir)< -0.95) return false;}float xyhitdiff=abs(sr->vtx.elastic.fuzzyk.png[0].shwlid.nhitx-sr->vtx.elastic.fuzzyk.png[0].shwlid.nhity);float xyhitsum=sr->vtx.elastic.fuzzyk.png[0].shwlid.nhitx+sr->vtx.elastic.fuzzyk.png[0].shwlid.nhity;if(xyhitdiff/xyhitsum >=0.4) return false;float nhitall=0.0;for(unsigned int j=0;j< sr->vtx.elastic.fuzzyk.nshwlid;j++){nhitall+=sr->vtx.elastic.fuzzyk.png[j].shwlid.nhit;}if(nhitall/sr->slc.nhit<=0.7) return false;return true;})
const SpillCut kStandardSpillCuts
Apply this unless you&#39;re doing something special.
Definition: SpillCuts.h:49
const Var kMuonID(muonid_classifier::MuonID)
static Binning Simple(int n, double lo, double hi, const std::vector< std::string > &labels={})
Definition: Binning.cxx:107
const ana::Cut kNuebarCCIncFiducialLoose([](const caf::SRProxy *sr){if(!sr->vtx.elastic.IsValid) return false;if(sr->vtx.elastic.vtx.X()< -165.0) return false;if(sr->vtx.elastic.vtx.X() > 165.0) return false;if(sr->vtx.elastic.vtx.Y()< -165.0) return false;if(sr->vtx.elastic.vtx.Y() > 165.0) return false;if(sr->vtx.elastic.vtx.Z()< 100.0) return false;if(sr->vtx.elastic.vtx.Z() > 1000.0) return false;return true;})
Cut CutFromNuTruthCut(const NuTruthCut &stc)
Definition: Cut.cxx:7
const Var kXSecCVWgt2020
Definition: XsecTunes.h:105