13 #include "CAFAna/Extrap/ExtrapSterile.h" 71 const HistAxis kNCHighEAxis(
"Calorimetric Energy (GeV)",
73 const HistAxis kHighENumuCCAxis(
"Reconstructed Neutrino Energy (GeV)",
76 std::string labelRecoE =
"Calorimetric Energy (GeV)";
79 "sel.cosrej.numucontpid",
"sel.nuecosrej.partptp",
80 "vtx.elastic.fuzzyk.nshwlid",
"vtx.elastic.fuzzyk.png.shwlid.lid.ismuon",
81 "slc.calE",
"slc.nhit",
"sel.nuecosrej.starttop",
82 "sel.nuecosrej.stoptop" 89 if(numucontpid > 0.5)
return false;
90 if(numucontpid <= 0.42)
return false;
91 if(partptp >= 0.8)
return false;
94 if(nhit <= 0.)
return false;
95 if(
sr->
slc.
calE/nhit <= 0.018)
return false;
100 std::map<std::string, Sample> sidebands;
102 kNCHighEAxis, kHighENumuCCAxis,
117 std::map<std::string, IDecomp*> ncdecomps;
118 std::map<std::string, IDecomp*> numudecomps;
119 std::map<std::string, ModularExtrapSterile*> extrap123bs;
120 std::map<std::string, ModularExtrapSterile*> extrap3c3ds;
121 std::map<std::string, PredictionSterile*> predst123bs;
122 std::map<std::string, PredictionSterile*> predst3c3ds;
123 std::map<std::string, PredictionCombinePeriods*> preds;
124 std::map<std::string, Spectrum*>
cosmics;
125 std::map<std::string, Spectrum*> dataspecs;
127 for(
const auto& sideband : sidebands) {
129 Sample sample = sideband.second;
132 loaderNDMC, loaderNDdata,
133 *sample.axnc, *sample.nd,
138 loaderNDMC, loaderNDdata,
145 loaderNDMC, loaderFDMC_swp, loaderFDMC_non, loaderFDMC_tau,
146 *ncdecomps[samplelabel], *numudecomps[samplelabel],
147 *sample.axnc, *sample.axnumu,
148 *sample.fd, *sample.nd,
kNumuND,
155 loaderNDMC, loaderFDMC_swp_3c, loaderFDMC_non_3c, loaderFDMC_tau_3c,
156 *ncdecomps[samplelabel], *numudecomps[samplelabel],
157 *sample.axnc, *sample.axnumu,
158 *sample.fd, *sample.nd,
kNumuND,
166 cosmics[samplelabel] =
new Spectrum(loaderFDdata, *sample.axnc,
179 loaderFDMC_non_3c.
Go();
180 loaderFDMC_swp_3c.
Go();
181 loaderFDMC_tau_3c.
Go();
188 for(
const auto& sideband : sidebands) {
197 TFile* rootF =
new TFile(outfile.c_str(),
"RECREATE");
200 TDirectory*
tmp = gDirectory;
201 TDirectory* saveDir = gDirectory;
205 for(
const auto& cafanaobj : ncdecomps) {
206 dir =
"ncdecomp" + sep + cafanaobj.first;
207 saveDir = rootF->mkdir(dir.c_str());
208 cafanaobj.second->SaveTo(saveDir);
211 for(
const auto& cafanaobj : numudecomps) {
212 dir =
"numudecomp" + sep + cafanaobj.first;
213 saveDir = rootF->mkdir(dir.c_str());
214 cafanaobj.second->SaveTo(saveDir);
217 for(
const auto& cafanaobj : extrap123bs) {
218 dir =
"extrap123b" + sep + cafanaobj.first;
219 saveDir = rootF->mkdir(dir.c_str());
220 cafanaobj.second->SaveTo(saveDir);
223 for(
const auto& cafanaobj : extrap3c3ds) {
224 dir =
"extrap3c3d" + sep + cafanaobj.first;
225 saveDir = rootF->mkdir(dir.c_str());
226 cafanaobj.second->SaveTo(saveDir);
229 for(
const auto& cafanaobj : predst123bs) {
230 dir =
"predst123b" + sep + cafanaobj.first;
231 saveDir = rootF->mkdir(dir.c_str());
232 cafanaobj.second->SaveTo(saveDir);
235 for(
const auto& cafanaobj : predst3c3ds) {
236 dir =
"predst3c3d" + sep + cafanaobj.first;
237 saveDir = rootF->mkdir(dir.c_str());
238 cafanaobj.second->SaveTo(saveDir);
241 for(
const auto& cafanaobj : preds) {
242 dir =
"pred" + sep + cafanaobj.first;
243 saveDir = rootF->mkdir(dir.c_str());
244 cafanaobj.second->SaveTo(saveDir);
247 for(
const auto& cafanaobj : cosmics) {
248 dir =
"cosmic" + sep + cafanaobj.first;
249 saveDir = rootF->mkdir(dir.c_str());
250 cafanaobj.second->SaveTo(saveDir);
253 for(
const auto& cafanaobj : dataspecs) {
254 dir =
"dataspec" + sep + cafanaobj.first;
255 saveDir = rootF->mkdir(dir.c_str());
256 cafanaobj.second->SaveTo(saveDir);
267 ret.axnc =
new HistAxis(axnc.label, axnc.bins, axnc.var);
268 ret.axnumu =
new HistAxis(axnm.label, axnm.bins, axnm.var);
269 ret.nd =
new Cut(nd);
270 ret.fd =
new Cut(fd);
277 std::cout <<
"!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!" <<
std::endl;
279 std::cout <<
" Be very cautious about what distributions you are making," <<
std::endl;
281 std::cout <<
"!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!" <<
std::endl;
283 std::vector<std::string> fname1 =
Wildcard(
284 diskdirFD +
"by_period/prod_restricteddecaf_S16-05-20_fd_numi_fhc_period1_nus_contain_v1_goodruns_prod2-snapshot/" 285 +
"prod_restricteddecaf_S16-05-20_fd_numi_fhc_period1_nus_contain_v1_goodruns_prod2-snapshot*.root" 287 std::vector<std::string> fname2 =
Wildcard(
288 diskdirFD +
"by_period/prod_restricteddecaf_S16-05-20_fd_numi_fhc_period2_nus_contain_v1_goodruns_prod2-snapshot/" 289 +
"prod_restricteddecaf_S16-05-20_fd_numi_fhc_period2_nus_contain_v1_goodruns_prod2-snapshot*.root" 291 std::vector<std::string> fname3b =
Wildcard(
292 diskdirFD +
"by_period/prod_restricteddecaf_S16-05-20_fd_numi_fhc_epoch3*_nus_contain_v1_goodruns_prod2-snapshot.root" 294 std::vector<std::string> fname3c =
Wildcard(
295 diskdirFD +
"by_period/prod_restricteddecaf_S16-05-20_fd_numi_fhc_epoch3c_nus_contain_v1_goodruns_prod2-snapshot/" 296 +
"prod_restricteddecaf_S16-05-20_fd_numi_fhc_epoch3c_nus_contain_v1_goodruns_prod2-snapshot*.root" 298 std::vector<std::string> fname3d =
Wildcard(
299 diskdirFD +
"by_period/prod_restricteddecaf_S16-05-20_fd_numi_fhc_epoch3d_nus_contain_v1_goodruns_prod2-snapshot/" 300 +
"prod_restricteddecaf_S16-05-20_fd_numi_fhc_epoch3d_nus_contain_v1_goodruns_prod2-snapshot*.root" 302 std::vector<std::string>
ret;
303 for(
const auto&
file : fname1) { ret.push_back(
file); }
304 for(
const auto&
file : fname2) { ret.push_back(
file); }
305 for(
const auto&
file : fname3b) { ret.push_back(
file); }
306 for(
const auto&
file : fname3c) { ret.push_back(
file); }
307 for(
const auto&
file : fname3d) { ret.push_back(
file); }
_HistAxis< Var > HistAxis
Represent the binning of a Spectrum's x-axis.
Cuts and Vars for the 2020 FD DiF Study.
const double kSecondAnaPeriod2POT
Sample(std::string shortName, std::string longName, Cut cut)
SRCosRej cosrej
Output from CosRej (Cosmic Rejection)
std::vector< std::string > MakeUnblindList()
std::vector< SRFuzzyKProng > png
Vector of 3D prong objects.
Proxy for caf::StandardRecord.
const double kSecondAnaEpoch3dPOT
const double kSecondAnaEpoch3bPOT
const std::vector< std::string > fFDMC_non
void SetSpillCut(const SpillCut &cut)
std::pair< Spectrum *, CheatDecomp * > make_pair(SpectrumLoaderBase &loader_data, SpectrumLoaderBase &loader_mc, HistAxis *axis, Cut *cut, const SystShifts &shift, const Var &wei)
const double kSecondAnaPeriod1POT
const Cut kNusNCSel([](const caf::SRProxy *sr){if(sr->slc.nhit >=200) return false;if(sr->slc.nhit< 20) return false;if(sr->sel.cvn.ncid< 0.2) return false;return true;})
Cut that is more CC rejection than NC selection from docdb 15285.
const Cut kInTimingSideband([](const caf::SRProxy *sr){if(sr->spill.run > util::kLastBadTimingRun) return(kInTimingSideband_before(sr)|| kInTimingSideband_after(sr));else return(kInTimingSideband_before(sr)|| kInTimingSideband_afterA(sr)|| kInTimingSideband_afterB(sr));}, [](const caf::SRSpillProxy *spill){if(spill->run > util::kLastBadTimingRun) return(kInTimingSideband_before.Livetime(spill)+ kInTimingSideband_after.Livetime(spill));else return(kInTimingSideband_before.Livetime(spill)+ kInTimingSideband_afterA.Livetime(spill)+ kInTimingSideband_afterB.Livetime(spill));}, [](const caf::SRSpillProxy *spill){return 0;})
const std::string fnamenear_concat
std::vector< std::string > Wildcard(const std::string &wildcardString)
Find files matching a UNIX glob, plus expand environment variables.
unsigned int nshwlid
number of shwlid showers - either 0 or number of 3d prongs
const std::string fnameneardata_concat
const Var kCaloE
Summed calorimetric energy of all hits in slice, uncorrected.
const std::vector< std::string > fFDMC_swp
SRNueCosRej nuecosrej
Output from NueCosRej (Nue Cosmic Rejection)
const std::string fFDMC_non_3c
float calE
Calorimetric energy of the cluster [GeV].
Sample MakeSample(const HistAxis &axnc, const HistAxis &axnm, const Cut &nd, const Cut &fd)
_Cut< caf::SRProxy > Cut
Representation of a cut (selection) to be applied to a caf::StandardRecord object.
const std::string diskdirFD
virtual void Go() override
Load all the registered spectra.
const HistAxis kNCAxis("Calorimetric Energy (GeV)", kNCDisappearanceEnergyBinning, kCaloE)
Axes used in Ana01 analysis by nus group.
std::vector< float > Spectrum
unsigned int nhit
number of hits
const SystShifts kNoShift
const std::string fFDMC_tau_3c
const HistAxis kNCBinsNumuCCAxis("Reconstructed Neutrino Energy (GeV)", kNCDisappearanceEnergyBinning, kCCE)
Splits Data proportionally according to MC.
Collaborates with Spectrum and OscillatableSpectrum to fill spectra from CAF files.
void SideBandLoad(std::string outfile)
const std::vector< std::string > fFDMC_tau
SRIDBranch sel
Selector (PID) branch.
SRElastic elastic
Single vertex found by Elastic Arms.
const std::string fFDMC_swp_3c
const Cut kInBeamSpill([](const caf::SRProxy *sr){if(sr->spill.run > util::kLastBadTimingRun) return kInBeamSpill_main(sr);else return kInBeamSpill_main(sr)||kInBeamSpill_shifted(sr);}, [](const caf::SRSpillProxy *spill){if(spill->run > util::kLastBadTimingRun) return kInBeamSpill_main.Livetime(spill);else return kInBeamSpill_main.Livetime(spill)+kInBeamSpill_shifted.Livetime(spill);}, [](const caf::SRSpillProxy *spill) -> double{return spill->spillpot;})
Does the event fall inside the window we call the beam spill?
const Cut kNusCosRejMod([](const caf::SRProxy *sr){double partptp=sr->sel.nuecosrej.partptp;if(sr->vtx.elastic.fuzzyk.nshwlid==0) return false;if(partptp >=0.8) return false;if(sr->vtx.elastic.fuzzyk.png[0].shwlid.lid.ismuon==1) return false;double nhit=(double) sr->slc.nhit;if(nhit<=0.) return false;if(sr->slc.calE/nhit<=0.009) return false;return true;})
SRSlice slc
Slice branch: nhit, extents, time, etc.
A prediction object compatible with sterile oscillations.
const double kSecondAnaEpoch3cPOT
SRFuzzyK fuzzyk
Primary 3D prong object.
const SpillCut kStandardSpillCuts
Apply this unless you're doing something special.
Sum MC predictions from different periods scaled according to data POT targets.
T min(const caf::Proxy< T > &a, T b)
static Binning Simple(int n, double lo, double hi, const std::vector< std::string > &labels={})
const Cut kNusCosRej([](const caf::SRProxy *sr){double numucontpid2019=sr->sel.cosrej.numucontpid2019;double partptp=sr->sel.nuecosrej.partptp;if(sr->vtx.elastic.fuzzyk.nshwlid==0) return false;if(numucontpid2019<=0.5) return false;if(partptp >=0.8) return false;if(sr->vtx.elastic.fuzzyk.png[0].shwlid.lid.ismuon==1) return false;double nhit=(double) sr->slc.nhit;if(nhit<=0.) return false;if(sr->slc.calE/nhit<=0.018) return false;if(std::min(sr->sel.nuecosrej.starttop, sr->sel.nuecosrej.stoptop)< 480) return false;return true;})
Cosmic rejection for the NC sample from docdb 15241.
SRVertexBranch vtx
Vertex branch: location, time, etc.