4 #include "CAFAna/Core/HistAxis.h" 19 std::pair<Spectrum*, CheatDecomp*>
make_pair(
28 std::map<std::string, std::pair<Spectrum*, CheatDecomp*> >
MakeMap(
31 std::map<std::string, HistAxis* >
axes,
32 std::map<std::string, Cut* >
cuts,
37 std::map<std::string, std::pair<Spectrum*, CheatDecomp*> >
MakeMap(
39 std::map<std::string, HistAxis*>
axes,
40 std::map<std::string, Cut*>
cuts,
52 const std::string fneardatafull =
"defname: prod_caf_R16-11-12-feature_caf_size.a_nd_numi_fhc_full_v1_goodruns with limit 1";
55 const std::string fnearcaffull =
"defname: prod_caf_R16-11-12-feature_caf_size.a_nd_genie_nonswap_genierw_fhc_nova_v08_full_v1 with limit 1";
65 std::string labelGap =
"Primary Shower Gap to Vertex (cm)";
69 std::string labelCont =
"Minimum Distance to Any Detector Face (cm)";
78 std::string labelTop =
"Distance to Top of Detector (cm)";
93 return sr->vtx.elastic.fuzzyk.gap;
123 float mindist = 10000.;
130 return (mindist > 9999. ? (
float)-1. : mindist);
149 if(nhit <= 0.) {
return 0.; }
158 const HistAxis kNCAllEAxis(
"Calorimetric Energy (GeV)", kNCAllEBins,
kCaloE);
159 const HistAxis kAxisHpP( labelHpP, kHpPBins, kHpP);
160 const HistAxis kAxisGap( labelGap, kGapBins, kGap);
164 const HistAxis kAxisCont(labelCont, kNearEdgeBins, kMinDistP);
166 const HistAxis kAxisLID( labelLID, kFinePIDBins, kElecID);
169 const HistAxis kAxisNCP( labelNCP, kFinePIDBins, kNCP);
170 const HistAxis kAxisPTP( labelPTP, kFinePIDBins, kPartPTP);
184 const Cut kFidMinusX(
190 if(vtx.Y() < -100.0)
return false;
191 if(vtx.Y() > 100.0)
return false;
192 if(vtx.Z() < 200.0)
return false;
193 if(vtx.Z() > 1000.0)
return false;
196 const Cut kFidMinusY(
202 if(vtx.X() < -100.0)
return false;
203 if(vtx.X() > 100.0)
return false;
204 if(vtx.Z() < 200.0)
return false;
205 if(vtx.Z() > 1000.0)
return false;
208 const Cut kFidMinusZ(
214 if(vtx.X() < -100.0)
return false;
215 if(vtx.X() > 100.0)
return false;
216 if(vtx.Y() < -100.0)
return false;
217 if(vtx.Y() > 100.0)
return false;
221 const Cut kNCSelMinusHit(
227 const Cut kNCSelMinusCVN(
230 if(sr->
slc.
nhit >= 200)
return false;
231 if(sr->
slc.
nhit < 20)
return false;
235 const Cut kCosRejMinusPTP(
241 if(
nhit <= 0.)
return false;
245 const Cut kCosRejMinusIsM(
249 if(partptp >= 0.8)
return false;
251 if(nhit <= 0.)
return false;
252 if(sr->
slc.
calE/nhit <= 0.009)
return false;
255 const Cut kCosRejMinusEpH(
260 if(partptp >= 0.8)
return false;
266 std::map<std::string, HistAxis*>
axes;
268 axes[
"ECalLong"] =
new HistAxis(kNCAllEAxis);
269 axes[
"HitpPlane"] =
new HistAxis(kAxisHpP);
270 axes[
"Gap"] =
new HistAxis(kAxisGap);
271 axes[
"VtxX"] =
new HistAxis(kAxisVtxX);
272 axes[
"VtxY"] =
new HistAxis(kAxisVtxY);
273 axes[
"VtxZ"] =
new HistAxis(kAxisVtxZ);
274 axes[
"Cont"] =
new HistAxis(kAxisCont);
275 axes[
"CVN"] =
new HistAxis(kAxisCVN);
276 axes[
"RemID"] =
new HistAxis(kAxisRem);
277 axes[
"LID"] =
new HistAxis(kAxisLID);
278 axes[
"NHit"] =
new HistAxis(kAxisNHit);
279 axes[
"EpH"] =
new HistAxis(kAxisEpH);
280 axes[
"PTP"] =
new HistAxis(kAxisPTP);
281 axes[
"IsMu"] =
new HistAxis(kAxisIsMu);
282 axes[
"DistTop"] =
new HistAxis(kAxisTop);
285 std::map<std::string, Cut*>
cuts;
288 cuts[
"EQFid"] =
new Cut(kEQFid);
290 cuts[
"PreselNCSel"] =
new Cut(kPreselNCSel);
295 std::map<std::string, std::pair<ana::Spectrum*, ana::CheatDecomp*>> loaded_specs =
MakeMap(
299 std::map<std::string, Cut*> nm1s;
300 nm1s[
"VtxX"] =
new Cut(kNM1Fid && kFidMinusX);
301 nm1s[
"VtxY"] =
new Cut(kNM1Fid && kFidMinusY);
302 nm1s[
"VtxZ"] =
new Cut(kNM1Fid && kFidMinusZ);
303 nm1s[
"Cont"] =
new Cut(kNM1Cont);
304 nm1s[
"CVN"] =
new Cut(kNM1NCSel && kNCSelMinusCVN);
305 nm1s[
"NHit"] =
new Cut(kNM1NCSel && kNCSelMinusHit);
306 nm1s[
"PTP"] =
new Cut(kNM1CosRej && kCosRejMinusPTP);
307 nm1s[
"EpH"] =
new Cut(kNM1CosRej && kCosRejMinusEpH);
308 nm1s[
"IsMu"] =
new Cut(kNM1CosRej && kCosRejMinusIsM);
309 for(
const auto& nm1 : nm1s) {
313 loaded_specs[strcut + strNM1] =
make_pair(
325 TFile* rootOut =
new TFile(
"/nova/ana/steriles/Ana01/DataMC.root",
"RECREATE");
327 for(
const auto& loaded_spec : loaded_specs) {
328 loaded_spec.second.first->SaveTo(rootOut, (loaded_spec.first +
"Data").c_str());
329 loaded_spec.second.second->SaveTo(rootOut, (loaded_spec.first +
"MC").c_str());
343 std::pair<Spectrum*, CheatDecomp*>
output;
359 std::map<std::string, std::pair<Spectrum*, CheatDecomp*> >
MakeMap(
362 std::map<std::string, HistAxis* >
axes,
363 std::map<std::string, Cut* >
cuts,
368 std::map<std::string, std::pair<Spectrum*, CheatDecomp*> >
output;
369 for(
const auto&
cut : cuts) {
370 for(
const auto&
axis : axes) {
372 loader_data, loader_mc,
373 axis.second,
cut.second, shift, wei
382 std::map<std::string, std::pair<Spectrum*, CheatDecomp*> >
MakeMap(
384 std::map<std::string, HistAxis*>
axes,
385 std::map<std::string, Cut*>
cuts,
caf::Proxy< unsigned int > nshwlid
Near Detector underground.
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
std::map< std::string, std::pair< Spectrum *, CheatDecomp * > > MakeMap(SpectrumLoaderBase &loader_data, SpectrumLoaderBase &loader_mc, std::map< std::string, HistAxis * > axes, std::map< std::string, Cut * > cuts, const SystShifts &shift, const Var &wei)
const Binning kRemidBinning
Binning for plotting remid attractively.
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);})
Simple record of shifts applied to systematic parameters.
Proxy for caf::StandardRecord.
Collection of SpectrumLoaders for many configurations.
const Cut kNusNDFiducial([](const caf::SRProxy *sr){assert(sr->vtx.elastic.IsValid &&"Must apply DQ cuts");if(sr->vtx.elastic.vtx.X()< -100.0) return false;if(sr->vtx.elastic.vtx.X() > 100.0) return false;if(sr->vtx.elastic.vtx.Y()< -100.0) return false;if(sr->vtx.elastic.vtx.Y() > 100.0) return false;if(sr->vtx.elastic.vtx.Z()< 200.0) return false;if(sr->vtx.elastic.vtx.Z() > 1000.0) return false;return true;})
ND Fiducial volume from docdb 15242.
std::pair< Spectrum *, CheatDecomp * > make_pair(SpectrumLoaderBase &loader_data, SpectrumLoaderBase &loader_mc, HistAxis *axis, Cut *cut, const SystShifts &shift, const Var &wei)
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
int isnan(const stan::math::var &a)
void Go()
Call Go() on all the loaders.
caf::Proxy< caf::SRElastic > elastic
caf::Proxy< caf::SRNueCosRej > nuecosrej
caf::Proxy< std::vector< caf::SRFuzzyKProng > > png
void SetSpillCut(const SpillCut &cut)
caf::Proxy< unsigned int > nhit
const Var kCaloE
Summed calorimetric energy of all hits in slice, uncorrected.
caf::Proxy< caf::SRCVNResult > cvn
SpectrumLoaderBase & GetLoader(caf::Det_t det, DataMC datamc, DataSource src=kBeam, SwappingConfig swap=kNonSwap)
Retrieve a specific loader.
const Cut kIsMuon([](const caf::SRProxy *sr){return(abs(sr->mc.nu[0].beam.ptype)==13);})
_Cut< caf::SRProxy > Cut
Representation of a cut (selection) to be applied to a caf::StandardRecord object.
const Binning kBooleanBins
caf::Proxy< float > partptp
const HistAxis kNCAxis("Calorimetric Energy (GeV)", kNCDisappearanceEnergyBinning, kCaloE)
Axes used in Ana01 analysis by nus group.
std::vector< float > Spectrum
Output from Cosmic Rejection (Nuecosrej) module.
const SystShifts kNoShift
Base class for the various types of spectrum loader.
caf::Proxy< bool > IsValid
caf::Proxy< caf::SRSlice > slc
const Cut kNusNDContain([](const caf::SRProxy *sr){const caf::SRNueCosRejProxy &cr=sr->sel.nuecosrej;if(std::min(cr.starteast, cr.stopeast)< 25) return false;if(std::min(cr.startwest, cr.stopwest)< 25) return false;if(std::min(cr.starttop, cr.stoptop) < 25) return false;if(std::min(cr.startbottom, cr.stopbottom)< 25) return false;if(std::min(cr.startfront, cr.stopfront)< 25) return false;if(std::min(cr.startback, cr.stopback)< 25) return false;return true;})
Containment variable for NC events from docdb 15242.
#define SIMPLEVAR(CAFNAME)
For Vars where literally all you need is a single CAF variable.
caf::Proxy< caf::SRVector3D > vtx
const Binning kFinePIDBins
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;})
const SpillCut kStandardSpillCuts
Apply this unless you're doing something special.
caf::Proxy< caf::SRIDBranch > sel
T min(const caf::Proxy< T > &a, T b)
caf::Proxy< caf::SRVertexBranch > vtx
void SetLoaderPath(const std::string &path, caf::Det_t det, DataMC datamc, DataSource src=kBeam, SwappingConfig swap=kNonSwap)
Configure loader via wildcard path.
static Binning Simple(int n, double lo, double hi, const std::vector< std::string > &labels={})
Just return the ND truth spectra as the decomposition.
const Cut kNoCut
The simplest possible cut: pass everything, used as a default.