12 if(rawEM < 0. || rawHA < 0.)
return -5.;
13 const double p0 = 0.0;
14 const double p1 = 9.54388e-01;
15 const double p2 = 1.13311;
16 const double p3 = 0.0;
17 const double p4 = 3.39617e-06;
18 const double p5 = 3.49874e-06;
19 const double recoEn = (1./(1+0.0089272))*(rawHA*rawHA*p5 + rawEM*rawEM*p4 + rawEM*rawHA*p3 + rawHA*p2 + rawEM*p1 + p0);
26 if(rawEM < 0. || rawHA < 0.)
return -5.;
27 const double p0 = 0.0;
28 const double p1 = 9.88258e-01;
29 const double p2 = 1.20084;
30 const double p3 = 0.0;
31 const double p4 = 1.92904e-07;
32 const double p5 = 1.20704e-07;
33 const double recoEn = (1./(1+0.0111873))*(rawHA*rawHA*p5 + rawEM*rawEM*p4 + rawEM*rawHA*p3 + rawHA*p2 + rawEM*p1 + p0);
40 if(rawEM < 0. || rawHA < 0.)
return -5.;
41 const double p0 = 0.0;
42 const double p1 = 9.24605e-01;
43 const double p2 = 1.15426;
44 const double p3 = 0.0;
45 const double p4 = 1.58171e-02;
46 const double p5 = 6.81363e-03;
47 const double recoEn = (1./(1+0.0217326))*(rawHA*rawHA*p5 + rawEM*rawEM*p4 + rawEM*rawHA*p3 + rawHA*p2 + rawEM*p1 + p0);
53 if(rawEM < 0. || rawHA < 0.)
return -5.;
54 const double p0 = 0.0;
55 const double p1 = 1.01777;
56 const double p2 = 1.10868;
57 const double p3 = 0.0;
58 const double p4 = 1.43541e-03;
59 const double p5 = 1.09628e-01;
60 const double recoEn = (1./(1+0.0355622))*(rawHA*rawHA*p5 + rawEM*rawEM*p4 + rawEM*rawHA*p3 + rawHA*p2 + rawEM*p1 + p0);
69 double CVNem_CalE = 0.0;
71 double png_CalE = png.shwlid.calE;
72 double emPID = ( (double)png.cvnpart.photonid +(double)png.cvnpart.electronid);
74 if ( emPID <= 0 ) continue;
75 if ( emPID >= 0.5 ) CVNem_CalE += png_CalE;
79 double png_CalE = png.calE;
80 double emPID = ( (double)png.cvnpart.photonid + (double)png.cvnpart.electronid);
82 if ( emPID <= 0 ) continue;
83 if ( emPID >= 0.7 ) CVNem_CalE += png_CalE;
85 if(CVNem_CalE <=0 ||
kLongestProng(sr) >= 500) CVNem_CalE = prim_png.shwlid.calE;
94 double CVNhad_CalE = 0.0;
97 double png_CalE = png.shwlid.calE;
98 double emPID = ( (double)png.cvnpart.photonid +(double)png.cvnpart.electronid);
100 if ( emPID < 0.5 ) CVNhad_CalE += png_CalE;
105 double png_CalE = png.calE;
106 double emPID = ( (double)png.cvnpart.photonid + (double)png.cvnpart.electronid);
108 if ( emPID < 0.7 ) CVNhad_CalE += png_CalE;
113 if (CVNhad_CalE<0) abort();
124 double CVNem_CalE = 0.0;
127 double png_CalE = png.shwlid.calE;
129 double emPID = ( (double)png.cvnpart.photonid +
130 (double)png.cvnpart.pizeroid +
131 (double)png.cvnpart.electronid);
132 double haPID = ( (double)png.cvnpart.protonid +
133 (double)png.cvnpart.pionid +
134 (double)png.cvnpart.neutronid +
135 (double)png.cvnpart.otherid +
136 (double)png.cvnpart.muonid);
138 if ( emPID < 0 ) continue;
139 if ( emPID >= haPID ) CVNem_CalE += png_CalE;
152 const double CVNha_CalE = sr->
slc.
calE;
223 if ( sr->spill.isRHC )
268 float ct = elecDir.Dot( beamDir );
270 if ( p < 0 ||
fabs( ct ) > 1 )
return -1000.f;
caf::Proxy< size_t > npng
T max(const caf::Proxy< T > &a, T b)
caf::Proxy< caf::SRFuzzyK > fuzzyk
Far Detector at Ash River.
Cuts and Vars for the 2020 FD DiF Study.
fvar< T > fabs(const fvar< T > &x)
const Var kNueEnergyRHC2020([](const caf::SRProxy *sr){if(!sr->vtx.elastic.IsValid) return-1.0;if(sr->vtx.elastic.fuzzyk.npng< 1) return-1.0;return NueRecoE_2020RHCFit(kEME_2020(sr), kHADE_2020(sr));})
caf::Proxy< std::vector< caf::SRProng > > png2d
caf::Proxy< caf::SRHeader > hdr
Proxy for caf::StandardRecord.
const Var kNueEnergy2020_2D3D([](const caf::SRProxy *sr){if(!sr->vtx.elastic.IsValid) return-1.0;if(sr->vtx.elastic.fuzzyk.npng< 1) return-1.0;if(kIsRHC(sr)) return kNueEnergyRHC2020_2D3D(sr);else return kNueEnergyFHC2020_2D3D(sr);})
Nue energy with 3d prong and 2d prong info. (new version of vars)
Proxy for caf::SRFuzzyKProng.
const Var kEME2020_2D3D([](const caf::SRProxy *sr){if(!sr->vtx.elastic.IsValid) return-1.0;if(sr->vtx.elastic.fuzzyk.npng< 1) return-1.0;const caf::SRFuzzyKProngProxy &prim_png=sr->vtx.elastic.fuzzyk.png[0];double CVNem_CalE=0.0;for(const caf::SRFuzzyKProngProxy &png:sr->vtx.elastic.fuzzyk.png){double png_CalE=png.shwlid.calE;double emPID=((double) png.cvnpart.photonid+(double) png.cvnpart.electronid);if(emPID<=0) continue;if(emPID >=0.5 ) CVNem_CalE+=png_CalE;}for(const caf::SRProngProxy &png:sr->vtx.elastic.fuzzyk.png2d){double png_CalE=png.calE;double emPID=((double) png.cvnpart.photonid+(double) png.cvnpart.electronid);if(emPID<=0) continue;if(emPID >=0.7 ) CVNem_CalE+=png_CalE;}if(CVNem_CalE<=0||kLongestProng(sr) >=500) CVNem_CalE=prim_png.shwlid.calE;return CVNem_CalE;})
double NueRecoE_2020RHCFit_2D3D(double rawEM, double rawHA)
double NueRecoE_2020FHCFit(double rawEM, double rawHA)
const Var kHADE2020_2D3D([](const caf::SRProxy *sr){if(!sr->vtx.elastic.IsValid) return-1.0;if(sr->vtx.elastic.fuzzyk.npng< 1) return-1.0;double CVNhad_CalE=0.0;double em=0.0;for(const caf::SRFuzzyKProngProxy &png:sr->vtx.elastic.fuzzyk.png){double png_CalE=png.shwlid.calE;double emPID=((double) png.cvnpart.photonid+(double) png.cvnpart.electronid);if(emPID< 0.5 ) CVNhad_CalE+=png_CalE;else em+=png_CalE;}for(const caf::SRProngProxy &png:sr->vtx.elastic.fuzzyk.png2d){double png_CalE=png.calE;double emPID=((double) png.cvnpart.photonid+(double) png.cvnpart.electronid);if(emPID< 0.7 ) CVNhad_CalE+=png_CalE;else em+=png_CalE;}CVNhad_CalE+=sr->vtx.elastic.fuzzyk.orphCalE;if(em<=0||kLongestProng(sr) >=500) CVNhad_CalE-=sr->vtx.elastic.fuzzyk.png[0].shwlid.calE;if(CVNhad_CalE< 0) abort();return CVNhad_CalE;})
caf::Proxy< caf::SRElastic > elastic
caf::Proxy< std::vector< caf::SRFuzzyKProng > > png
const Var kLongestProng([](const caf::SRProxy *sr){if(!sr->vtx.elastic.IsValid) return 0.f;if(sr->vtx.elastic.fuzzyk.npng==0) return 0.f;auto idx=sr->vtx.elastic.fuzzyk.longestidx;return float(sr->vtx.elastic.fuzzyk.png[idx].len);})
const Var kNueEnergyRHC2020_2D3D([](const caf::SRProxy *sr){if(!sr->vtx.elastic.IsValid) return-1.0;if(sr->vtx.elastic.fuzzyk.npng< 1) return-1.0;return NueRecoE_2020RHCFit_2D3D(kEME2020_2D3D(sr), kHADE2020_2D3D(sr));})
double NueRecoE_2020RHCFit(double rawEM, double rawHA)
const Var kNueElectronE2020([](const caf::SRProxy *sr){if(!sr->vtx.elastic.IsValid) return-1000.f;if(sr->vtx.elastic.fuzzyk.npng< 1) return-1000.f;float x=sr->vtx.elastic.fuzzyk.png[0].calE;float energy=x;double a0, a1, a2, a3, a4, a5, a6;if(sr->hdr.det==caf::kFARDET){if(sr->spill.isRHC){a0=-0.727211;a1=2.319804;a2=-2.697100;a3=1.604802;a4=-0.493659;a5=0.075183;a6=-0.004478;}else{a0=-0.428873;a1=1.486957;a2=-1.772140;a3=1.114639;a4=-0.358208;a5=0.056065;a6=-0.003376;}}else{return energy;}if(x > 0) energy=x-(a0+a1 *x+a2 *pow(x, 2)+ a3 *pow(x, 3)+a4 *pow(x, 4)+ a5 *pow(x, 5)+a6 *pow(x, 6));if(energy< 0) energy=x;return energy;})
const Var kHADE_2020([](const caf::SRProxy *sr){if(!sr->vtx.elastic.IsValid) return-1.0;if(sr->vtx.elastic.fuzzyk.npng< 1) return-1.0;const double CVNha_CalE=sr->slc.calE;return std::max(CVNha_CalE-kEME_2020(sr), 0.);})
const Cut kIsRHC([](const caf::SRProxy *sr){return sr->spill.isRHC;})
caf::Proxy< float > orphCalE
const Var kNueEnergy2020([](const caf::SRProxy *sr){if(!sr->vtx.elastic.IsValid) return-1.0;if(sr->vtx.elastic.fuzzyk.npng< 1) return-1.0;if(kIsRHC(sr)) return kNueEnergyRHC2020(sr);else return kNueEnergyFHC2020(sr);})
Nue energy with 3d prong info. only (old version of vars)
caf::Proxy< bool > IsValid
caf::Proxy< caf::SRSlice > slc
const Var kNueEnergyFHC2020_2D3D([](const caf::SRProxy *sr){if(!sr->vtx.elastic.IsValid) return-1.0;if(sr->vtx.elastic.fuzzyk.npng< 1) return-1.0;return NueRecoE_2020FHCFit_2D3D(kEME2020_2D3D(sr), kHADE2020_2D3D(sr));})
const Var kNueElectronPt2020([](const caf::SRProxy *sr){if(!sr->vtx.elastic.IsValid) return-1000.f;if(sr->vtx.elastic.fuzzyk.npng< 1) return-1000.f;TVector3 elecDir=(TVector3) sr->vtx.elastic.fuzzyk.png[0].dir;TVector3 beamDir=ana::NuMIBeamDirection(sr->hdr.det);float ct=elecDir.Dot(beamDir);float p=kNueElectronE2020(sr);if(p< 0||fabs(ct) > 1) return-1000.f;return float(p *sin(acos(ct)));})
TVector3 NuMIBeamDirection(caf::Det_t det)
Average direction of NuMI neutrinos in a given detector This function is a copy of geo::GeometryBase:...
double NueRecoE_2020FHCFit_2D3D(double rawEM, double rawHA)
const Var kNueEnergyFHC2020([](const caf::SRProxy *sr){if(!sr->vtx.elastic.IsValid) return-1.0;if(sr->vtx.elastic.fuzzyk.npng< 1) return-1.0;return NueRecoE_2020FHCFit(kEME_2020(sr), kHADE_2020(sr));})
caf::Proxy< caf::SRVertexBranch > vtx
const Var kEME_2020([](const caf::SRProxy *sr){if(!sr->vtx.elastic.IsValid) return-1.0;if(sr->vtx.elastic.fuzzyk.npng< 1) return-1.0;double CVNem_CalE=0.0;const caf::SRFuzzyKProngProxy &prim_png=sr->vtx.elastic.fuzzyk.png[0];for(const caf::SRFuzzyKProngProxy &png:sr->vtx.elastic.fuzzyk.png){double png_CalE=png.shwlid.calE;double emPID=((double) png.cvnpart.photonid+ (double) png.cvnpart.pizeroid+ (double) png.cvnpart.electronid);double haPID=((double) png.cvnpart.protonid+ (double) png.cvnpart.pionid+ (double) png.cvnpart.neutronid+ (double) png.cvnpart.otherid+ (double) png.cvnpart.muonid);if(emPID< 0) continue;if(emPID >=haPID ) CVNem_CalE+=png_CalE;else continue;}if(kLongestProng(sr) >=500) CVNem_CalE=prim_png.shwlid.calE;return CVNem_CalE;})