12 #include "CAFAna/Extrap/ExtrapSterile.h" 75 std::string labelTop =
"Distance to Top of Detector (cm)";
116 if(nhit <= 0.) {
return 0.; }
125 const HistAxis kNCAllEAxis(
"Calorimetric Energy (GeV)", kNCAllEBins,
kCaloE);
126 const HistAxis kNumuCCAllEAxis(
"Reconstructed Neutrino Energy (GeV)", kNCAllEBins,
kCCE);
132 const HistAxis kAxisNCP( labelNCP, kFatNCPBins, kNCP);
133 const HistAxis kAxisPTP( labelPTP, kFatPIDBins, kPartPTP);
142 const Cut kFidMinusX(
148 if(vtx.Y() < -720.0)
return false;
149 if(vtx.Y() > 600.0)
return false;
150 if(vtx.Z() < 50.0)
return false;
151 if(vtx.Z() > 5450.0)
return false;
154 const Cut kFidMinusY(
160 if(vtx.X() < -680.0)
return false;
161 if(vtx.X() > 650.0)
return false;
162 if(vtx.Z() < 50.0)
return false;
163 if(vtx.Z() > 5450.0)
return false;
166 const Cut kFidMinusZ(
172 if(vtx.X() < -680.0)
return false;
173 if(vtx.X() > 650.0)
return false;
174 if(vtx.Y() < -720.0)
return false;
175 if(vtx.Y() > 600.0)
return false;
179 const Cut kNCSelMinusHit(
185 const Cut kNCSelMinusCVN(
188 if(sr->
slc.
nhit >= 200)
return false;
189 if(sr->
slc.
nhit < 20)
return false;
193 const Cut kCosRejMinusNCP(
198 if(partptp >= 0.8)
return false;
201 if(
nhit <= 0.)
return false;
206 const Cut kCosRejMinusPTP(
209 double numucontpid = sr->
sel.
cosrej.numucontpid;
211 if(numucontpid <= 0.5)
return false;
214 if(
nhit <= 0.)
return false;
219 const Cut kCosRejMinusEpH(
222 double numucontpid = sr->
sel.
cosrej.numucontpid;
225 if(numucontpid <= 0.5)
return false;
226 if(partptp >= 0.8)
return false;
231 const Cut kCosRejMinusTop(
234 double numucontpid = sr->
sel.
cosrej.numucontpid;
237 if(numucontpid <= 0.5)
return false;
238 if(partptp >= 0.8)
return false;
241 if(
nhit <= 0.)
return false;
249 const Cut kCosRejModMinusPTP(
255 if(
nhit <= 0.)
return false;
259 const Cut kCosRejModMinusEpH(
264 if(partptp >= 0.8)
return false;
267 if(
nhit <= 0.)
return false;
271 std::map<std::string, Sample> samples;
273 kNCAllEAxis, kNumuCCAllEAxis,
278 kNusND, kNM1Fid && kFidMinusX
282 kNusND, kNM1Fid && kFidMinusY
286 kNusND, kNM1Fid && kFidMinusZ
290 kNM1NCSelND && kNCSelMinusCVN,
291 kNM1NCSel && kNCSelMinusCVN
295 kNusND, kNM1CosRej && kCosRejMinusNCP
299 kNM1CosRejMod && kCosRejModMinusPTP,
300 kNM1CosRej && kCosRejMinusPTP
304 kNM1CosRejMod && kCosRejModMinusEpH,
305 kNM1CosRej && kCosRejMinusEpH
309 kNusND, kNM1CosRej && kCosRejMinusTop
312 std::map<std::string, IDecomp*> ncdecomps;
313 std::map<std::string, IDecomp*> numudecomps;
314 std::map<std::string, ModularExtrapSterile*> extrap123bs;
315 std::map<std::string, ModularExtrapSterile*> extrap3c3ds;
316 std::map<std::string, PredictionSterile*> predst123bs;
317 std::map<std::string, PredictionSterile*> predst3c3ds;
318 std::map<std::string, PredictionCombinePeriods*> preds;
319 std::map<std::string, Spectrum*>
cosmics;
320 std::map<std::string, Spectrum*> dataspecs;
322 for(
const auto& sample : samples) {
326 loaderNDMC, loaderNDdata,
327 *sample.second.axnc, *sample.second.nd,
332 loaderNDMC, loaderNDdata,
333 *sample.second.axnumu,
kNumuND,
339 loaderNDMC, loaderFDMC_swp, loaderFDMC_non, loaderFDMC_tau,
340 *ncdecomps[axlabel], *numudecomps[axlabel],
341 *sample.second.axnc, *sample.second.axnumu,
342 *sample.second.fd, *sample.second.nd,
kNumuND,
349 loaderNDMC, loaderFDMC_swp_3c, loaderFDMC_non_3c, loaderFDMC_tau_3c,
350 *ncdecomps[axlabel], *numudecomps[axlabel],
351 *sample.second.axnc, *sample.second.axnumu,
352 *sample.second.fd, *sample.second.nd,
kNumuND,
360 cosmics[axlabel] =
new Spectrum(loaderFDdata, *sample.second.axnc,
364 dataspecs[axlabel] =
new Spectrum(loaderFDdata, *sample.second.axnc,
374 loaderFDMC_non_3c.
Go();
375 loaderFDMC_swp_3c.
Go();
376 loaderFDMC_tau_3c.
Go();
383 for(
const auto& sample : samples) {
393 TFile* rootF =
new TFile(outfile.c_str(),
"RECREATE");
396 TDirectory*
tmp = gDirectory;
397 TDirectory* saveDir = gDirectory;
401 for(
const auto& sample : samples) {
404 dir =
"ncdecomp" + sep + axlabel;
405 saveDir = rootF->mkdir(dir.c_str());
406 ncdecomps[axlabel]->SaveTo(saveDir);
408 dir =
"numudecomp" + sep + axlabel;
409 saveDir = rootF->mkdir(dir.c_str());
410 numudecomps[axlabel]->SaveTo(saveDir);
412 dir =
"extrap123b" + sep + axlabel;
413 saveDir = rootF->mkdir(dir.c_str());
414 extrap123bs[axlabel]->SaveTo(saveDir);
416 dir =
"extrap3c3d" + sep + axlabel;
417 saveDir = rootF->mkdir(dir.c_str());
418 extrap3c3ds[axlabel]->SaveTo(saveDir);
420 dir =
"predst123b" + sep + axlabel;
421 saveDir = rootF->mkdir(dir.c_str());
422 predst123bs[axlabel]->SaveTo(saveDir);
424 dir =
"predst3c3d" + sep + axlabel;
425 saveDir = rootF->mkdir(dir.c_str());
426 predst3c3ds[axlabel]->SaveTo(saveDir);
428 dir =
"pred" + sep + axlabel;
429 saveDir = rootF->mkdir(dir.c_str());
430 preds[axlabel]->SaveTo(saveDir);
432 dir =
"cosmic" + sep + axlabel;
433 saveDir = rootF->mkdir(dir.c_str());
434 cosmics[axlabel]->SaveTo(saveDir);
436 dir =
"dataspec" + sep + axlabel;
437 saveDir = rootF->mkdir(dir.c_str());
438 dataspecs[axlabel]->SaveTo(saveDir);
449 ret.axnc =
new HistAxis(axnc.label, axnc.bins, axnc.var);
450 ret.axnumu =
new HistAxis(axnm.label, axnm.bins, axnm.var);
451 ret.nd =
new Cut(nd);
452 ret.fd =
new Cut(fd);
caf::Proxy< unsigned int > nshwlid
const Cut kNusFDFiducial([](const caf::SRProxy *sr){ if(!sr->vtx.elastic.IsValid) return false;const TVector3 vtx=sr->vtx.elastic.vtx;if(vtx.X()< -680.0) return false;if(vtx.X() > 650.0) return false;if(vtx.Y()< -720.0) return false;if(vtx.Y() > 500.0) return false;if(vtx.Z()< 50.0) return false;if(vtx.Z() > 5450.0) return false;return true;})
FD Fiducial volume from docdb 15285.
caf::Proxy< caf::SRFuzzyK > fuzzyk
_HistAxis< Var > HistAxis
Represent the binning of a Spectrum's x-axis.
Cuts and Vars for the 2020 FD DiF Study.
caf::Proxy< float > stoptop
void ExtendedAxesLoad(std::string outfile)
const double kSecondAnaPeriod2POT
Sample(std::string shortName, std::string longName, Cut cut)
const Cut kNusEventQuality([](const caf::SRProxy *sr){if(!sr->vtx.elastic.IsValid) return false;if(sr->sel.nuecosrej.hitsperplane >=8) return false;if(sr->vtx.elastic.fuzzyk.nshwlid==0) return false;if(sr->vtx.elastic.fuzzyk.png[0].shwlid.gap >=100.) return false;if(sr->slc.ncontplanes<=2) return false;return true;})
Data Quality cuts from docdb 14241.
const Var kDistTop([](const caf::SRProxy *sr){return std::min(sr->sel.nuecosrej.starttop, sr->sel.nuecosrej.stoptop);})
Proxy for caf::StandardRecord.
const double kSecondAnaEpoch3dPOT
const double kSecondAnaEpoch3bPOT
const std::vector< std::string > fFDMC_non
const Cut kNusFDContain([](const caf::SRProxy *sr){const caf::SRNueCosRejProxy &cr=sr->sel.nuecosrej;if(std::min(cr.starteast, cr.stopeast)< 10) return false;if(std::min(cr.startwest, cr.stopwest)< 10) return false;if(std::min(cr.starttop, cr.stoptop) < 10) return false;if(std::min(cr.startbottom, cr.stopbottom)< 10) return false;if(std::min(cr.startfront, cr.stopfront)< 10) return false;if(std::min(cr.startback, cr.stopback)< 10) return false;return true;})
Containment variable for NC events from docdb 14241.
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 Var kEPerHit([](const caf::SRProxy *sr){if(sr->slc.nhit >0) return 1000.0 *(sr->slc.calE/sr->slc.nhit);else return-5.;})
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.
caf::Proxy< float > starttop
caf::Proxy< caf::SRCosRej > cosrej
int isnan(const stan::math::var &a)
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;})
caf::Proxy< caf::SRElastic > elastic
caf::Proxy< caf::SRNueCosRej > nuecosrej
caf::Proxy< std::vector< caf::SRFuzzyKProng > > png
const std::string fnamenear_concat
caf::Proxy< unsigned int > nhit
const std::string fnameneardata_concat
const Var kCaloE
Summed calorimetric energy of all hits in slice, uncorrected.
const std::vector< std::string > fFDMC_swp
caf::Proxy< caf::SRCVNResult > cvn
const std::string fFDMC_non_3c
_Cut< caf::SRProxy > Cut
Representation of a cut (selection) to be applied to a caf::StandardRecord object.
caf::Proxy< float > partptp
virtual void Go() override
Load all the registered spectra.
std::vector< float > Spectrum
const SystShifts kNoShift
const std::string fFDMC_tau_3c
const HistAxis kNCBinsNumuCCAxis("Reconstructed Neutrino Energy (GeV)", kNCDisappearanceEnergyBinning, kCCE)
caf::Proxy< bool > IsValid
caf::Proxy< caf::SRSlice > slc
Splits Data proportionally according to MC.
Collaborates with Spectrum and OscillatableSpectrum to fill spectra from CAF files.
const std::vector< std::string > fFDMC_tau
#define SIMPLEVAR(CAFNAME)
For Vars where literally all you need is a single CAF variable.
caf::Proxy< caf::SRVector3D > vtx
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?
A prediction object compatible with sterile oscillations.
const double kSecondAnaEpoch3cPOT
const SpillCut kStandardSpillCuts
Apply this unless you're doing something special.
Sum MC predictions from different periods scaled according to data POT targets.
caf::Proxy< caf::SRIDBranch > sel
T min(const caf::Proxy< T > &a, T b)
caf::Proxy< caf::SRVertexBranch > vtx
static Binning Simple(int n, double lo, double hi, const std::vector< std::string > &labels={})
Sample MakeSample(const HistAxis &axnc, const HistAxis &axnm, const Cut &nd, const Cut &fd)
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.
const std::vector< std::string > fnamefardata_unblind(MakeUnblindList())