27 #include <boost/math/special_functions.hpp> 28 #include <boost/math/special_functions/gamma.hpp> 44 #include "Utilities/func/MathUtil.h" 55 inline double getErr(
double PE) {
56 return 231594.0/(2267.16+
pow(PE,2.01011)) + 9.55689;
59 inline double getDist(
double x1,
double y1,
double z1,
60 double x2,
double y2,
double z2) {
61 return sqrt((x1-x2)*(x1-x2)+(y1-y2)*(y1-y2)+(z1-z2)*(z1-z2));
65 if (lhs.
Y() == rhs.
Y()) {
66 if (lhs.
Z() == rhs.
Z()) {
67 return lhs.
T() < rhs.
T();
69 else return lhs.
Z() < rhs.
Z();
71 else return lhs.
Y() < rhs.
Y();
74 bool pairHitComp (std::pair<rb::CellHit, rb::RecoHit> lhp,
75 std::pair<rb::CellHit, rb::RecoHit> rhp) {
125 double getLLR(std::set< std::pair<rb::CellHit,rb::RecoHit>,
126 bool(*)(std::pair<rb::CellHit,rb::RecoHit>,
127 std::pair<rb::CellHit,rb::RecoHit>)> ,
128 std::vector<rb::RecoHit> &outliers,
129 double &P_up,
double &P_dn,
double &
chi2,
double &slope);
130 void LLR(std::vector<double>& eT,
131 std::vector<double>& mT,
132 std::vector<double>& mTErr,
double& slope,
double&
chi2,
133 double& P_up,
double& P_dn, std::vector<rb::RecoHit> &
in,
134 std::vector<rb::RecoHit> &outliers);
136 void LinFit(
const std::vector<double>&
x,
137 const std::vector<double>&
y,
double *fitpar);
138 void LinFit(
const std::vector<double>&
x,
139 const std::vector<double>&
y,
140 const std::vector<double>& ye,
double *fitpar);
145 deltaT,
eT,
PE,
PECorr,
ADC,
Plane,
Cell,
View,
GoodTiming,
Reco,
170 , fHitsInstance (p.
get<
std::
string>(
"HitInstanceLabel"))
171 , fSliceModuleLabel (p.
get<
std::
string>(
"SliceModuleLabel"))
172 , fSliceInstance (p.
get<
std::
string>(
"SliceInstanceLabel"))
173 , fTrackModuleLabel (p.
get<
std::
string>(
"TrackModuleLabel"))
174 , fTrackToSliceInstance(p.
get<
std::
string>(
"TrackToSliceInstance"))
175 , fFuzzyKVertexLabel (p.
get<
std::
string>(
"FuzzyKVertexLabel"))
176 , fFuzzyKVertexInstance(p.
get<
std::
string>(
"FuzzyKVertexInstance"))
178 , fMinY (p.
get<double>(
"MinY"))
179 , fMinX (p.
get<double>(
"MinX"))
180 , fMinZ (p.
get<double>(
"MinZ"))
181 , fMaxY (p.
get<double>(
"MaxY"))
182 , fMaxX (p.
get<double>(
"MaxX"))
183 , fMaxZ (p.
get<double>(
"MaxZ"))
184 , _TrackLen (p.
get<double>(
"TrackLen"))
185 , _TrackHitsXY (p.
get<unsigned>(
"TrackHitsXY"))
186 , _TrackHitsX (p.
get<unsigned>(
"TrackHitsX"))
187 , _TrackHitsY (p.
get<unsigned>(
"TrackHitsY"))
191 , _Chi2 (p.
get<double>(
"Chi2"))
192 , _LLR (p.
get<double>(
"LLR"))
193 , _R2X (p.
get<double>(
"R2X"))
194 , _R2Y (p.
get<double>(
"R2Y"))
195 , _MinSlope (p.
get<double>(
"MinSlope"))
196 , _MaxSlope (p.
get<double>(
"MaxSlope"))
197 , fIsSim (p.
get<bool>(
"IsSim"))
198 , fParticleModuleLabel (p.
get<
std::
string>(
"ParticleModuleLabel"))
199 , fParticleInstance (p.
get<
std::
string>(
"ParticleInstance"))
214 "Run:SubRun:Event:SliceID:TrackID:Nhits:NRecohits:NOutliers:ProbUp:ProbDn:LLR:Chi2:Slope:LLRX:Chi2X:SlopeX:LLRY:Chi2Y:SlopeY:R2X:R2Y:StartX:StartY:StartZ:StartT:EndX:EndY:EndZ:EndT:TrackHitsX:TrackHitsY:Length:dirX:dirY:dirZ:eleAngle:totalE:containment:avgTX:avgTY:meRecoGeV:meRecoNHits:meRecoX:meRecoY:meRecoZ:meRecoDist:meRecoDeltaT:meRecoAtTrackEnd");
216 "Run:SubRun:Event:TrackID:HitID:X:Y:Z:T:deltaT:eT:PE:PECorr:ADC:Plane:Cell:View:GoodTiming:Reco:Outlier:dirX:dirY:dirZ:eleAngle:GeV:Diblock");
219 "Run:SubRun:Event:Ntracks:totalHits:totalE:containment:X:Y:Z:T");
222 "Run:SubRun:Event:Ntracks:NcontainedT:Nvertices:NcontainedV:avgNhits:avgNRecohits:avgLen:hasUpMu:hasFCUpMu:hasTGUpMu:nDiblocks");
225 "Run:SubRun:Event:SliceID:Ntracks:containment");
229 "Run:SubRun:Event:ParticleID:PDG:Length:StartX:StartY:StartZ:StartT:StartE:EndX:EndY:EndZ:EndT:EndE:StartPX:StartPY:StartPZ:StartEleAngle:containment:MotherID:TrackID:EDep:nFLS");
244 for(
int i = 0;
i < (
int)slicecol->size();
i++){
252 int nSlices = slces.
size();
255 for(
int s = 0;
s < nSlices;
s++){
256 for(
int c = 0;
c<(
int)slces[
s]->NCell();
c++)
266 e.
getByLabel(fHitsLabel, fHitsInstance, hits);
280 e.
getByLabel(fFuzzyKVertexLabel, fFuzzyKVertexInstance, fuzzyKprongs);
283 std::vector<art::Ptr<rb::Vertex> > vertices;
286 float hasUpMu=0, hasFCUpMu=0, hasTGUpMu=0;
288 size_t nContainedT = 0, nVertices = 0, nContainedV = 0;
289 float avgNhits = 0, avgNRecohits = 0, avgLen = 0;
291 for (
size_t i_prong=0; i_prong<fuzzyKprongs->size(); ++i_prong) {
292 auto theVertexVec = fmVert.at(i_prong);
294 for(
size_t i_vecvert=0; i_vecvert < theVertexVec.size(); ++i_vecvert) {
300 for (
size_t i_vert=0; i_vert<vertices.size(); ++i_vert) {
301 if (&(*vertices.at(i_vert)) == &(*theVertex)) {
305 vertices.push_back(theVertex);
312 vertTotalE=0, containment=4;
335 if (theVertex->
GetX() < fMaxX && theVertex->
GetX() > fMinX &&
336 theVertex->
GetY() < fMaxY && theVertex->
GetY() > fMinY &&
339 else containment = 1;
341 if (containment == 4) ++nContainedV;
344 (
float)e.
id().
event(), 0, totalHits, vertTotalE, containment,
351 unsigned int nGoodTiming = 0;
352 for (
size_t i_hit=0; i_hit < hits->size(); ++i_hit) {
353 if (hits->at(i_hit).GoodTiming()) ++nGoodTiming;
361 std::vector<const rb::Cluster*> slices;
362 std::vector<unsigned> tracks_per_slice;
363 std::vector<float> slice_containment;
365 std::vector<hit_info> allHitInfo;
366 std::vector<time_window> goodTimeWindows;
374 std::vector<float> activeDiblocks;
375 for (
size_t i_track=0; i_track < tracks->size(); ++i_track) {
376 bool keepHits =
true;
377 rb::Track theTrack = tracks->at(i_track);
381 bool found_slice =
false;
382 const rb::Cluster *theSlice = &(*fmSlice.at(i_track));
383 for (
size_t i_slice=0; i_slice<slices.size(); ++i_slice) {
385 if (theSlice == slices.at(i_slice))
386 { slice = (
float)i_slice; found_slice =
true;
387 ++tracks_per_slice.at(i_slice);
388 if(containment != 4) slice_containment.at(i_slice) = 1;
392 slice = (
float)slices.size();
393 slices.push_back(theSlice);
394 tracks_per_slice.push_back(1);
396 slice_containment.push_back(1);
398 slice_containment.push_back(4);
401 if (!theTrack.
Is3D())
continue;
407 std::set< std::pair<rb::CellHit,rb::RecoHit>,
408 bool(*)(std::pair<rb::CellHit,rb::RecoHit>,
409 std::pair<rb::CellHit,rb::RecoHit>)>
411 std::set< std::pair<rb::CellHit,rb::RecoHit>,
412 bool(*)(std::pair<rb::CellHit,rb::RecoHit>,
413 std::pair<rb::CellHit,rb::RecoHit>)>
415 std::set< std::pair<rb::CellHit,rb::RecoHit>,
416 bool(*)(std::pair<rb::CellHit,rb::RecoHit>,
417 std::pair<rb::CellHit,rb::RecoHit>)>
420 std::vector<double> x_hit;
421 std::vector<double> y_hit;
422 std::vector<double> zy_hit;
423 std::vector<double> zx_hit;
425 std::vector<art::Ptr<me::TrkME> > michelE = fo.at(i_track);
428 float meRecoGeV = -5;
429 float meRecoHits = -5;
430 float meRecoPos[3] = {-5, -5, -5};
431 float meRecoDist = -5;
432 float meRecoDeltaT = -5;
434 if(!michelE.empty()){
435 meRecoGeV = michelE[0]->TotalGeV();
436 meRecoHits = michelE[0]->NCell();
437 meRecoPos[0] = michelE[0]->MeanX();
438 meRecoPos[1] = michelE[0]->MeanY();
439 meRecoPos[2] = michelE[0]->MeanZ();
440 meRecoDist = michelE[0]->DistToTrk();
441 meRecoDeltaT = michelE[0]->DeltaT();
445 for (
size_t i_hit=0; i_hit < trackHits.
size(); ++i_hit) {
451 unsigned short iC = theHit->
Cell();
452 unsigned short iP = theHit->
Plane();
456 zx_hit.push_back(iP);
461 zy_hit.push_back(iP);
473 LinFit(x_hit, zx_hit, fitpar);
474 double r2x = fitpar[2];
475 LinFit(y_hit, zy_hit, fitpar);
476 double r2y = fitpar[2];
477 unsigned nxhit = x_hit.size();
478 unsigned nyhit = y_hit.size();
480 std::vector<rb::RecoHit> outliers, outliersX, outliersY;
481 double P_up, P_dn,
chi2, slope;
483 double llr =
getLLR(sortedTrackHits, outliers, P_up, P_dn, chi2,
486 double chi2X, slopeX, P_upX, P_dnX, chi2Y, slopeY, P_upY, P_dnY,
487 llrX =
getLLR(sortedTrackHitsX, outliersX, P_upX, P_dnX,
489 llrY =
getLLR(sortedTrackHitsY, outliersY, P_upY, P_dnY,
494 if (sortedTrackHits.size() >= 1) start = (sortedTrackHits.begin()->second);
496 if (sortedTrackHits.size() >= 2) end = (--sortedTrackHits.end())->
second;
505 float tDirX=0, tDirY=0, tDirZ=0, tEleAngle, trackTotalE=0;
508 if(chi2 > _Chi2) keepHits =
false;
509 if(dist < _TrackLen) keepHits =
false;
510 if(TMath::Abs(start.
X() - end.
X()) < _dX*cw) keepHits =
false;
511 if(TMath::Abs(start.
Y() - end.
Y()) < _dY*cw) keepHits =
false;
512 if(TMath::Abs(start.
Z() - end.
Z()) < _dZ*pw) keepHits =
false;
513 if(r2x < _R2X) keepHits =
false;
514 if(r2y < _R2Y) keepHits =
false;
515 if(nxhit < _TrackHitsX) keepHits =
false;
516 if(nyhit < _TrackHitsY) keepHits =
false;
521 window.
startT = start.
T() - 2000;
522 window.
endT = end.
T() + 2000;
523 goodTimeWindows.push_back(window);
526 float avgTX=0, avgTY=0;
529 for (
size_t i_hit=0; i_hit < trackHits.
size(); ++i_hit) {
531 float goodtime = theHit->
GoodTiming() ? 1 : 0;
533 float PE = theHit->
PE(),
539 float eleAngle =
atan(theDir.y()/
sqrt(theDir.x()*theDir.x() +
540 theDir.z()*theDir.z()));
542 tDirX += theDir.x()/trackHits.
size();
543 tDirY += theDir.y()/trackHits.
size();
544 tDirZ += theDir.z()/trackHits.
size();
550 GeV = theRecoHit.
GeV();
551 for (
size_t i_out=0; i_out < outliers.size(); ++i_out) {
552 if (!outliers.at(i_out).IsCalibrated())
continue;
553 if (outliers.at(i_out).PECorr() == theRecoHit.
PECorr() &&
554 outliers.at(i_out).X() == theRecoHit.
X() &&
555 outliers.at(i_out).Y() == theRecoHit.
Y() &&
556 outliers.at(i_out).Z() == theRecoHit.
Z() &&
557 outliers.at(i_out).T() == theRecoHit.
T())
558 { outlier = 1;
break; }
564 avgTX += theRecoHit.
T()/sortedTrackHitsX.size();
566 avgTY += theRecoHit.
T()/sortedTrackHitsY.size();
569 float eT = (
float)(
getDist(theRecoHit.
X(), theRecoHit.
Y(), theRecoHit.
Z(),
570 start.
X(), start.
Y(), start.
Z())/29.97),
579 info.
X = (
float)theRecoHit.
X();
580 info.
Y = (
float)theRecoHit.
Y();
581 info.
Z = (
float)theRecoHit.
Z();
582 info.
T = (
float)theRecoHit.
T() - sortedTrackHits.begin()->second.T();
599 info.
absT = theRecoHit.
T();
601 allHitInfo.push_back(info);
603 if(std::find(activeDiblocks.begin(),activeDiblocks.end(),
604 info.
Diblock)==activeDiblocks.end())
605 activeDiblocks.push_back(info.
Diblock);
609 tEleAngle =
atan(tDirY/
sqrt(tDirX*tDirX + tDirZ*tDirZ));
612 if (containment == 4) ++nContainedT;
613 avgNhits += trackHits.
size()/tracks->size();
614 avgNRecohits += sortedTrackHits.size()/tracks->size();
615 avgLen += dist/tracks->size();
617 float track_entries[48]=
620 (
float)trackHits.
size(), (
float)sortedTrackHits.size(),
622 (
float)
chi2, (
float)slope, (
float)llrX, (
float)chi2X, (
float)slopeX,
624 start.
X(), start.
Y(), start.
Z(), start.
T(), end.
X(), end.
Y(),
625 end.
Z(), end.
T(), (
float)nxhit,(
float)nyhit,
dist, tDirX, tDirY,
626 tDirZ, tEleAngle, trackTotalE, containment, avgTX, avgTY,
627 meRecoGeV, meRecoHits, meRecoPos[0], meRecoPos[1],
628 meRecoPos[2], meRecoDist, meRecoDeltaT, meAtEnd
635 for (
size_t i_slice = 0; i_slice < slices.size(); ++i_slice) {
638 (
float)tracks_per_slice.at(i_slice),
639 (
float)slice_containment.at(i_slice)};
643 for(
size_t i_hitInfo = 0; i_hitInfo < allHitInfo.size(); ++i_hitInfo){
646 for(
size_t i_window = 0; i_window < goodTimeWindows.size(); i_window++){
666 e.
getByLabel(fParticleModuleLabel, fParticleInstance, particles);
667 std::vector<sim::Particle> michelElectrons;
668 for (
size_t i_part=0; i_part<particles->size(); ++i_part) {
670 float startEleAngle =
atan(thePart.
Py() /
sqrt(thePart.
Px()*thePart.
Px() +
671 thePart.
Pz()*thePart.
Pz())),
675 thePart.
Py() > 0 &&
fabs(startEleAngle)*180/
M_PI > 10.0 &&
678 if (containment == 4) hasFCUpMu=1;
679 if (containment == 1) hasTGUpMu=1;
684 +(thePart.
Vy()-thePart.
EndY())*(thePart.
Vy()-thePart.
EndY())
685 +(thePart.
Vz()-thePart.
EndZ())*(thePart.
Vz()-thePart.
EndZ()));
696 for(
size_t i_fls = 0; i_fls < flsHits.size(); i_fls++){
759 float particle_entries[25] =
767 startEleAngle, containment,(
float)thePart.
Mother(),
775 (
float)nVertices, (
float)nContainedV, avgNhits, avgNRecohits,
776 avgLen, hasUpMu, hasFCUpMu, hasTGUpMu,
777 (
float) activeDiblocks.size() };
785 bool(*)(std::pair<rb::CellHit,rb::RecoHit>,
786 std::pair<rb::CellHit,rb::RecoHit>)>
788 std::vector<rb::RecoHit> &outliers,
789 double &P_up,
double &P_dn,
double &
chi2,
792 if (sortedHits.size() < 2) {
800 std::vector<double> measuredTimes;
801 std::vector<double> expectedTimes;
802 std::vector<double> wgts;
805 measuredTimes.push_back(0.0);
806 expectedTimes.push_back(0.0);
807 double err =
getErr(sortedHits.begin()->first.PE());
808 wgts.push_back(1.0/err/err);
810 double startX = sortedHits.begin()->second.X(),
811 startY = sortedHits.begin()->second.Y(),
812 startZ = sortedHits.begin()->second.Z(),
813 startT = sortedHits.begin()->second.T();
815 std::vector<rb::RecoHit>
in;
817 for (
auto i_hit = ++sortedHits.begin();
818 i_hit != sortedHits.end();
820 in.push_back(i_hit->second);
821 measuredTimes.push_back(i_hit->second.T() - startT);
822 double dist =
getDist(i_hit->second.X(), i_hit->second.Y(),
823 i_hit->second.Z(), startX, startY, startZ);
824 expectedTimes.push_back(dist/29.97);
825 err =
getErr(i_hit->first.PE());
826 wgts.push_back(1.0/err/err);
829 LLR(expectedTimes, measuredTimes, wgts, slope, chi2, P_up, P_dn,
832 return log(P_up/P_dn);
836 std::vector<double>& mT,
837 std::vector<double>& wgts,
double& slope,
838 double&
chi2,
double& P_up,
double& P_dn,
839 std::vector<rb::RecoHit>&
in,
840 std::vector<rb::RecoHit>& outliers)
843 size_t n = mT.size();
859 size_t totAllowOutlier = n/10;
860 size_t nOutliers = 0;
861 std::vector<double> x_filt, y_filt, w_filt;
862 for (
size_t i=0;
i <
n;
i++) {
863 double y_fit = a + b*eT.at(
i);
864 if (
fabs(mT.at(
i) - y_fit) < 5*(1.0/
sqrt(wgts.at(
i)))
865 || nOutliers>totAllowOutlier)
867 x_filt.push_back(eT.at(
i));
868 y_filt.push_back(mT.at(
i));
869 w_filt.push_back(wgts.at(
i));
873 outliers.push_back(in[
i]);
877 if (x_filt.size() >= 2)
902 chi2 /= (double)(n-2);
904 double one_over_ee = 0.0,
908 for (
size_t i=0;
i<
n; ++
i) {
910 one_over_ee += w_filt.at(
i);
911 x_over_ee += x_filt.at(
i)*w_filt.at(
i);
912 y_over_ee += y_filt.at(
i)*w_filt.at(
i);
915 double up_inter = (y_over_ee-x_over_ee)/one_over_ee;
916 double dn_inter = (y_over_ee+x_over_ee)/one_over_ee;
918 double up_chi2 = 0.0, dn_chi2 = 0.0;
919 for (
size_t i=0;
i<
n; ++
i) {
920 double e = 1.0/
sqrt(w_filt.at(
i));
921 double up_expected = up_inter + x_filt.at(
i);
922 double dn_expected = dn_inter - x_filt.at(
i);
923 up_chi2 +=
pow((y_filt.at(
i)-up_expected) / e, 2.0);
924 dn_chi2 +=
pow((y_filt.at(
i)-dn_expected) / e, 2.0);
930 if (prob_up < 1e-30) prob_up = 1e-30;
931 if (prob_dn < 1e-30) prob_dn = 1e-30;
967 TLorentzVector stop = TLorentzVector(thePart.
EndPosition());
969 if (start.Y() > stop.Y())
970 { TLorentzVector hold =
start; start = stop; stop = hold; }
972 if (start.Y() < fMinY || start.X() < fMinX || start.X() > fMaxX ||
973 start.Z() < fMinZ || start.Z() >
fMaxZ) {
974 if (stop.Y() > fMaxY ||
978 stop.Z() >
fMaxZ) {
return 1; }
981 if (stop.Y() > fMaxY || stop.X() < fMinX || stop.X() > fMaxX ||
982 stop.Z() < fMinZ || stop.Z() >
fMaxZ)
992 TVector3 stop = theTrack.
Stop();
997 if (start.y() > stop.y())
998 { TVector3 hold =
start; start = stop; stop = hold; }
1003 if (start.y() < fMinY || start.x() < fMinX || start.x() > fMaxX ||
1004 start.z() < fMinZ || start.z() >
fMaxZ) {
1005 if (stop.y() > fMaxY ||
1009 stop.z() >
fMaxZ) {
return 1; }
1012 if (stop.y() > fMaxY || stop.x() < fMinX || stop.x() > fMaxX ||
1013 stop.z() < fMinZ || stop.z() >
fMaxZ)
1023 const auto n = x_val.size();
1024 const auto x = (std::accumulate(x_val.begin(), x_val.end(), 0.0))/
n;
1025 const auto y = (std::accumulate(y_val.begin(), y_val.end(), 0.0))/
n;
1026 const auto xx = (std::inner_product(x_val.begin(), x_val.end(), x_val.begin(), 0.0))/
n;
1027 const auto yy = (std::inner_product(y_val.begin(), y_val.end(), y_val.begin(), 0.0))/
n;
1028 const auto xy = (std::inner_product(x_val.begin(), x_val.end(), y_val.begin(), 0.0))/
n;
1030 const auto b = (xy -
x*
y)/(
xx -
x*
x);
1031 const auto a =
y -
b*
x;
1032 const auto r = (xy - x*
y)/
sqrt((
xx - x*x)*(yy -
y*
y));
1041 int n = x_val.size();
1049 for (
int i=0;
i<
n;
i++ ){
1051 x = x + x_val[
i]/y_err[
i]/y_err[
i];
1052 y = y + y_val[
i]/y_err[
i]/y_err[
i];
1054 xx = xx + x_val[
i]*x_val[
i]/y_err[
i]/y_err[
i];
1055 yy = yy + y_val[
i]*y_val[
i]/y_err[
i]/y_err[
i];
1056 xy = xy + x_val[
i]*y_val[
i]/y_err[
i]/y_err[
i];
1057 ee = ee + 1./y_err[
i]/y_err[
i];
1060 const auto b = (xy*ee - x*
y)/(xx*ee - x*x);
1061 const auto a = (xy -
b*
xx)/x;
1062 const auto r = (xy - x*
y)/
sqrt((xx - x*x)*(yy - y*
y));
double E(const int i=0) const
void LinFit(const std::vector< double > &x, const std::vector< double > &y, double *fitpar)
const XML_Char XML_Encoding * info
back track the reconstruction to the simulation
const TLorentzVector & Position(const int i=0) const
double Py(const int i=0) const
std::string fSliceModuleLabel
float containmentType(rb::Track)
fvar< T > fabs(const fvar< T > &x)
const TLorentzVector & EndPosition() const
int getDiblock(int plane)
Float_t y1[n_points_granero]
enum geo::_plane_proj View_t
Enumerate the possible plane projections.
unsigned short Plane() const
Float_t x1[n_points_granero]
std::vector< sim::FLSHit > ParticleToFLSHit(const int &trackID) const
All the FLSHits that were created by the track id trackID, sorted from most to least light...
Vertical planes which measure X.
double Px(const int i=0) const
UpMuRecoAna(fhicl::ParameterSet const &p)
A collection of associated CellHits.
void LLR(std::vector< double > &eT, std::vector< double > &mT, std::vector< double > &mTErr, double &slope, double &chi2, double &P_up, double &P_dn, std::vector< rb::RecoHit > &in, std::vector< rb::RecoHit > &outliers)
A rb::Prong with full reconstructed trajectory.
A single unit of energy deposition in the liquid scintillator.
std::pair< Spectrum *, CheatDecomp * > make_pair(SpectrumLoaderBase &loader_data, SpectrumLoaderBase &loader_mc, HistAxis *axis, Cut *cut, const SystShifts &shift, const Var &wei)
DEFINE_ART_MODULE(TestTMapFile)
virtual TVector3 Start() const
Horizontal planes which measure Y.
art::PtrVector< rb::CellHit > AllCells() const
Get all cells from both views.
Calibrated quantities relying on position in the orthogonal view. To generate a rb::CellHit from a rb...
std::string fTrackToSliceInstance
virtual double TotalLength() const
Length (cm) of all the track segments.
unsigned short Cell() const
std::string fFuzzyKVertexInstance
std::string fParticleModuleLabel
void push_back(Ptr< U > const &p)
double getLLR(std::set< std::pair< rb::CellHit, rb::RecoHit >, bool(*)(std::pair< rb::CellHit, rb::RecoHit >, std::pair< rb::CellHit, rb::RecoHit >)>, std::vector< rb::RecoHit > &outliers, double &P_up, double &P_dn, double &chi2, double &slope)
bool recoHitComp(rb::RecoHit lhs, rb::RecoHit rhs)
std::string fHitsInstance
double T(const int i=0) const
rb::RecoHit RecoHit(const art::Ptr< rb::CellHit > &chit) const
Return calibrated hit based on assumed W coordinate.
static constexpr Double_t GeV
bool IsCalibrated() const
You MUST check here before accessing PECorr, MIP or GeV.
reference at(size_type n)
EDAnalyzer(Table< Config > const &config)
Vertex location in position and time.
EventNumber_t event() const
TNtuple * fParticleNtuple
void analyze(art::Event const &e) override
bool pairHitComp(std::pair< rb::CellHit, rb::RecoHit > lhp, std::pair< rb::CellHit, rb::RecoHit > rhp)
double Vx(const int i=0) const
SubRunNumber_t subRun() const
T * make(ARGS...args) const
std::string fTrackModuleLabel
std::string fSliceInstance
int16_t ADC(uint32_t i) const
bool getByLabel(std::string const &label, std::string const &productInstanceName, Handle< PROD > &result) const
double Pz(const int i=0) const
double LinFit(const std::vector< double > &x, const std::vector< double > &y, const std::vector< double > &w, double &m, double &c)
Find the best-fit line to a collection of points in 2-D by minimizing the squared vertical distance f...
double Vz(const int i=0) const
virtual TVector3 InterpolateDir(double z) const
std::string fMCTruthLabel
float GetEdep() const
Get total Energy deposited into the cell for the whole FLSHit.
std::string fMichelELabel
std::string fFuzzyKVertexLabel
TVector3 Stop() const
Position of the final trajectory point.
double getDist(double x1, double y1, double z1, double x2, double y2, double z2)
fvar< T > gamma_q(const fvar< T > &x1, const fvar< T > &x2)
double Vy(const int i=0) const
Encapsulate the geometry of one entire detector (near, far, ndos)
virtual bool Is3D() const
std::string fParticleInstance