containment_optimization.C
Go to the documentation of this file.
1 #pragma once
2 
4 
5 // analysis cuts/vars
10 
11 // for data pot
13 
14 // Multiverse includes
16 #include "CAFAna/Vars/XsecTunes.h"
20 
21 
22 // temporary measure to use prod4 trained muonid
23 #include "NDAna/muonid/NDXSecMuonPID.h"
24 namespace nuebarccinc {
25  const ana::Cut kMuonIDProd4Cut = ana::kMuonIDProd4 < -0.2;
26 }
27 
29 {
30  const int ret = gSystem->Load("libNDAnamuonid");
31  if(ret != 0) exit(1);
32 }
33 
34 using namespace ana;
36 {
38  // axis = 0, 1, 2, 3, 4, for top, bottom, east, west, front, respectively
39  std::map<TString, SpectrumLoader *> loaders;
40  loaders["kNominal"] = new SpectrumLoader(nuebarccinc::PROD5_MC_RHC_NOMINAL );
41  loaders["kCalibNeg"] = new SpectrumLoader(nuebarccinc::PROD5_MC_RHC_CALIB_DOWN );
42  loaders["kCalibPos"] = new SpectrumLoader(nuebarccinc::PROD5_MC_RHC_CALIB_UP );
43  loaders["kCalibShape"] = new SpectrumLoader(nuebarccinc::PROD5_MC_RHC_CALIB_SHAPE);
44  loaders["kLightUpCalibDown"] = new SpectrumLoader(nuebarccinc::PROD5_MC_RHC_LIGHT_UP );
45  loaders["kLightDownCalibUp"] = new SpectrumLoader(nuebarccinc::PROD5_MC_RHC_LIGHT_DOWN );
46  loaders["kCherenkov"] = new SpectrumLoader(nuebarccinc::PROD5_MC_RHC_CHERENKOV );
47 
48  std::vector<std::string> opti_labels = {"contain_top", "contain_bottom", "contain_east",
49  "contain_west", "contain_front"};
50  // define a hist axis to optimize
51  std::vector<const HistAxis *> opti_axes = {new HistAxis("Min. Prong Distance from Top (cm)",
52  Binning::Simple(40, 0, 400),
53  kDistAllTop),
54  new HistAxis("Min. Prong Distance from Bottom (cm)",
55  Binning::Simple(40, 0, 400),
57  new HistAxis("Min. Prong Distance from East (cm)",
58  Binning::Simple(40, 0, 400),
59  kDistAllEast),
60  new HistAxis("Min. Prong Distance from West (cm)",
61  Binning::Simple(40, 0, 400),
62  kDistAllWest),
63  new HistAxis("Min. Prong Distance from Front (cm)",
64  Binning::Simple(80, 0, 1200),
65  kDistAllFront)};
66 
67  // Define selection and signal cuts
70  // This demo is for fiducial cut optimization so the selection
71  // should be all preselection minus fiducial constraints
72  const Cut * kSelection = new Cut(nuebarccinc::kDQ &&
79 
80 
81  // central value weight to be applied to all samples
82  const Var * kCVWeight = new Var(kPPFXFluxCVWgt * kXSecCVWgt2020);
83 
84  // nominal definition
85  SystematicDef * sNominal =
86  new SystematicDef("sNominal",
87  loaders["kNominal"],
88  opti_axes[axis],
89  kCVWeight);
90 
91  // systematic definitions
92  SystematicDef * sLightlevelUpDown =
93  new SystematicDef("sLightlevelUpDown",
94  loaders["kLightUpCalibDown"],
95  loaders["kLightDownCalibUp"],
96  opti_axes[axis],
97  kCVWeight);
98  SystematicDef * sCalibPosNeg =
99  new SystematicDef("sCalibPosNeg",
100  loaders["kCalibPos"],
101  loaders["kCalibNeg"],
102  opti_axes[axis],
103  kCVWeight);
104  SystematicDef * sCalibShape =
105  new SystematicDef("sCalibShape",
106  loaders["kCalibShape"],
107  opti_axes[axis],
108  kCVWeight);
109  SystematicDef * sCherenkov =
110  new SystematicDef("sCherenkov",
111  loaders["kCherenkov"],
112  opti_axes[axis],
113  kCVWeight);
114 
115  // Setup the genie multiverses
116  // We take flux uncertainty to be constant, so no need for PPFX universe
117  GenieMultiverseParameters genie_multiverse(300, getAllXsecSysts_2020());
118  std::vector<SystShifts> genie_shifts = genie_multiverse.GetSystShifts();
119 
120  SystematicDef * sGenieMultiverse =
121  new SystematicDef("sGenieMultiverse",
122  loaders["kNominal"],
123  opti_axes[axis],
124  kCVWeight);
125 
126  // create the optimizer.
128  optimizer = new CutOptimization(opti_axes[axis],
129  kSignal,
130  kSelection);
131  optimizer->DefineNominal(sNominal);
132  optimizer->DefineSystematic(sLightlevelUpDown);
133  optimizer->DefineSystematic(sCalibPosNeg );
134  optimizer->DefineSystematic(sCalibShape );
135  optimizer->DefineSystematic(sCherenkov );
136  optimizer->DefineMultiverseSystematic(sGenieMultiverse, genie_shifts);
137  // optionally set spill cuts for all loaders
138  optimizer->SetSpillCuts(kStandardSpillCuts);
139 
140 
141  // Let the loaders go
142  optimizer->Go();
143 
144  // save to file
145  TFile * output = new TFile(("containment_optimization_" + opti_labels[axis] + ".root").c_str(), "recreate");
146  optimizer->SaveTo(output, opti_labels[axis].c_str());
147  output->Close();
148 
149 }
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 kDistAllBottom([](const caf::SRProxy *sr){if(std::isnan(1.*sr->sel.nuecosrej.distallpngbottom)) return-1000.0f;return float(sr->sel.nuecosrej.distallpngbottom);})
Distance of all showers in slice from the bottom edge of detector.
Definition: NueVars.h:33
const Var kPPFXFluxCVWgt
Definition: PPFXWeights.h:16
const Var kDistAllWest([](const caf::SRProxy *sr){if(std::isnan(1.*sr->sel.nuecosrej.distallpngwest)) return-1000.0f;return float(sr->sel.nuecosrej.distallpngwest);})
Distance of all showers in slice from the west edge of detector.
Definition: NueVars.h:36
const Var kDistAllTop([](const caf::SRProxy *sr){if(std::isnan(1.*sr->sel.nuecosrej.distallpngtop)) return-1000.0f;return float(sr->sel.nuecosrej.distallpngtop);})
Distance of all showers in slice from the top edge of detector.
Definition: NueVars.h:30
const ana::Cut kFrontPlanes
void load_libs_muonid()
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
const Var kDistAllEast([](const caf::SRProxy *sr){if(std::isnan(1.*sr->sel.nuecosrej.distallpngeast)) return-1000.0f;return float(sr->sel.nuecosrej.distallpngeast);})
Distance of all showers in slice from the east edge of detector.
Definition: NueVars.h:39
std::vector< const ISyst * > getAllXsecSysts_2020()
void Cut(double x)
Definition: plot_outliers.C:1
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
void SetSpillCuts(const SpillCut &)
const std::string PROD5_MC_RHC_CHERENKOV
const Cut kSelection
Definition: SINCpi0_Cuts.h:404
const Var kDistAllFront([](const caf::SRProxy *sr){if(std::isnan(1.*sr->sel.nuecosrej.distallpngfront)) return-1000.0f;return float(sr->sel.nuecosrej.distallpngfront);})
Distance of all showers in slice from the front edge of detector.
Definition: NueVars.h:45
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()
exit(0)
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 ana::Cut kMuonIDProd4Cut
Definition: datamc.C:28
const SpillCut kStandardSpillCuts
Apply this unless you&#39;re doing something special.
Definition: SpillCuts.h:49
void containment_optimization(int axis=0)
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