4 #include <boost/algorithm/string/classification.hpp> 5 #include <boost/algorithm/string/split.hpp> 16 #include "larevt/CalibrationDBI/Interface/ChannelStatusProvider.h" 17 #include "larevt/CalibrationDBI/Interface/ChannelStatusService.h" 49 using namespace detail;
65 std::vector<int> dtrs;
66 for (
auto& dtj : slc.
tjs) {
67 if (dtj.AlgMod[
kKilled])
continue;
68 if (dtj.ParentID != muTj.
ID)
continue;
69 dtrs.push_back(dtj.ID);
71 if (prt)
mf::LogVerbatim(
"TC") <<
"MakeHaloTj: Killing delta-ray T" << dtj.ID;
75 if (pfpIndex == USHRT_MAX) {
76 if (prt)
mf::LogVerbatim(
"TC") <<
" No PFP found for 3D-matched delta-ray";
79 auto& pfp = slc.
pfps[pfpIndex];
80 if (prt)
mf::LogVerbatim(
"TC") <<
" Killing delta-ray PFParticle P" << pfp.UID;
83 if (pfp.ParentUID > 0) {
85 if (parentIndx.first != USHRT_MAX) {
86 auto& parent =
slices[parentIndx.first].pfps[parentIndx.second];
87 std::vector<int> newDtrUIDs;
88 for (
auto uid : parent.DtrUIDs)
89 if (uid != dtj.UID) newDtrUIDs.push_back(uid);
90 parent.DtrUIDs = newDtrUIDs;
102 tj.
ID = slc.
tjs.size() + 1;
118 std::vector<int> closeTjs;
119 for (
unsigned short ipt = muTj.
EndPt[0]; ipt <= muTj.
EndPt[1]; ++ipt) {
120 auto tp = muTj.
Pts[ipt];
131 if (tp.Dir[0] != 0) window *=
std::abs(1 / tp.Dir[0]);
134 bool hitsAdded =
false;
135 for (
unsigned short ii = 0; ii < tp.Hits.size(); ++ii) {
136 unsigned int iht = tp.Hits[ii];
137 auto inTraj = slc.
slHits[iht].InTraj;
138 if (inTraj < 0)
continue;
140 tp.UseHit[ii] =
true;
146 if (inTraj != muTj.
ID &&
147 std::find(closeTjs.begin(), closeTjs.end(), inTraj) == closeTjs.end())
148 closeTjs.push_back(inTraj);
155 tj.
Pts.push_back(tp);
158 if (tj.
Pts.empty())
return;
162 myprt <<
"MHTj: T" << muTj.
ID <<
" npts " << tj.
Pts.size() <<
" close";
163 for (
auto tid : closeTjs)
164 myprt <<
" T" << tid;
168 slc.
tjs.push_back(tj);
205 for (
auto& tj : slc.
tjs) {
206 if (tj.AlgMod[
kKilled])
continue;
213 std::vector<int> temp;
214 for (
auto& vx3 : slc.
vtx3s) {
215 if (vx3.ID == 0)
continue;
218 temp.push_back(vx3.ID);
220 if (temp.empty())
return;
223 std::vector<int> masterlist;
224 for (
auto vx3id : temp) {
225 auto& vx3 = slc.
vtx3s[vx3id - 1];
228 for (
auto tjid : tjlist) {
229 auto& tj = slc.
tjs[tjid - 1];
230 if (tj.ParentID != 0) tj.ParentID = 0;
231 if (std::find(masterlist.begin(), masterlist.end(), tjid) == masterlist.end())
232 masterlist.push_back(tjid);
237 myprt <<
"DTP: masterlist Tjs";
238 for (
auto tjid : masterlist)
239 myprt <<
" " << tjid;
243 std::vector<SortEntry> sortVec(temp.size());
244 for (
unsigned short indx = 0; indx < temp.size(); ++indx) {
245 auto& vx3 = slc.
vtx3s[temp[indx] - 1];
246 sortVec[indx].index = indx;
247 sortVec[indx].val = vx3.Score;
249 if (sortVec.size() > 1) std::sort(sortVec.begin(), sortVec.end(),
valsDecreasing);
252 for (
unsigned short indx = 0; indx < temp.size(); ++indx)
253 vlist[indx] = temp[sortVec[indx].index];
257 auto& vx3 = slc.
vtx3s[vlist[0] - 1];
263 auto& sf = neutrinoPFP.SectionFits[0];
269 neutrinoPFP.PDGCode = 14;
270 neutrinoPFP.Vx3ID[1] = vx3.ID;
271 neutrinoPFP.Vx3ID[0] = vx3.ID;
274 if (!
StorePFP(slc, neutrinoPFP))
return;
278 std::vector<bool> lookedAt3(slc.
vtx3s.size() + 1,
false);
279 std::vector<bool> lookedAt2(slc.
vtxs.size() + 1,
false);
281 std::vector<std::pair<int, int>> pardtr;
283 for (
unsigned short indx = 0; indx < vlist.size(); ++indx) {
284 auto& vx3 = slc.
vtx3s[vlist[indx] - 1];
285 if (lookedAt3[vx3.ID])
continue;
287 lookedAt3[vx3.ID] =
true;
291 if (primTjList.empty())
continue;
293 for (
auto primTjID : primTjList) {
294 auto& primTj = slc.
tjs[primTjID - 1];
296 if (primTj.ParentID != -1)
continue;
297 if (prt)
mf::LogVerbatim(
"TC") <<
"Vx3 " << vx3.ID <<
" Primary tj " << primTj.ID;
302 for (
unsigned short end = 0;
end < 2; ++
end) {
303 if (primTj.VtxID[
end] == 0)
continue;
304 auto& vx2 = slc.
vtxs[primTj.VtxID[
end] - 1];
305 if (vx2.Vx3ID == vx3.ID)
continue;
308 for (
auto dtrID : dtrList) {
310 if (dtrID == primTjID)
continue;
311 auto& dtj = slc.
tjs[dtrID - 1];
312 if (dtj.ParentID != -1)
continue;
313 pardtr.push_back(std::make_pair(primTjID, dtrID));
314 if (prt)
mf::LogVerbatim(
"TC") <<
" primTj " << primTjID <<
" dtrID " << dtrID;
318 for (
unsigned short end = 0;
end < 2; ++
end) {
319 if (primTj.VtxID[
end] == 0)
continue;
320 auto& vx2 = slc.
vtxs[primTj.VtxID[
end] - 1];
324 if (pardtr.empty())
continue;
328 for (
auto pdtr : pardtr)
329 myprt <<
" " << pdtr.first <<
"_" << pdtr.second;
333 for (
unsigned short nit = 0; nit < 100; ++nit) {
334 auto lastPair = pardtr[pardtr.size() - 1];
335 auto& dtj = slc.
tjs[lastPair.second - 1];
336 dtj.ParentID = lastPair.first;
339 unsigned short dpt = 0, ppt = 0;
340 auto& ptj = slc.
tjs[lastPair.first - 1];
344 if (prt)
mf::LogVerbatim(
"TC") <<
"Set parent " << ptj.ID <<
" dtr " << dtj.ID;
348 for (
unsigned short end = 0;
end < 2; ++
end) {
349 if (dtj.VtxID[
end] == 0)
continue;
350 auto& vx2 = slc.
vtxs[dtj.VtxID[
end] - 1];
351 if (lookedAt2[vx2.ID])
continue;
352 lookedAt2[vx2.ID] =
true;
354 for (
auto tjid : tjlist) {
355 if (tjid == dtj.ID || tjid == ptj.ID)
continue;
356 pardtr.push_back(std::make_pair(dtj.ID, tjid));
359 myprt <<
" add par_dtr";
360 for (
auto pdtr : pardtr)
361 myprt <<
" " << pdtr.first <<
"_" << pdtr.second;
365 if (pardtr.empty())
break;
369 for (
auto tjid : masterlist) {
370 auto& tj = slc.
tjs[tjid - 1];
371 if (tj.ParentID < 0) tj.ParentID = tj.ID;
380 if (tjIDs.size() < 2)
return 1;
381 std::vector<float> plnchg(slc.
nPlanes);
382 for (
auto tjid : tjIDs) {
383 if (tjid <= 0 || tjid > (
int)slc.
tjs.size())
return 1;
384 auto& tj = slc.
tjs[tjid - 1];
387 plnchg[plane] += tj.TotChg;
391 for (
unsigned short plane = 0; plane < slc.
nPlanes; ++plane) {
392 if (plnchg[plane] == 0)
continue;
393 aveChg += plnchg[plane];
396 if (cnt < 2)
return 1;
399 for (
unsigned short plane = 0; plane < slc.
nPlanes; ++plane) {
401 if (plnchg[plane] == 0)
continue;
402 float asym =
std::abs(plnchg[plane] - aveChg) / (plnchg[plane] + aveChg);
403 if (asym > maxAsym) maxAsym = asym;
418 std::array<int, 5> codeList = {{0, 11, 13, 111, 211}};
419 unsigned short codeIndex = 0;
420 if (tjIDs.empty())
return codeList[codeIndex];
422 std::array<unsigned short, 5> cnts;
425 for (
auto tjid : tjIDs) {
426 if (tjid <= 0 || tjid > (
int)slc.
tjs.size())
continue;
427 auto& tj = slc.
tjs[tjid - 1];
428 for (
unsigned short ii = 0; ii < 5; ++ii)
429 if (tj.PDGCode == codeList[ii]) ++cnts[ii];
431 if (len > maxLen) maxLen = len;
435 for (
unsigned short ii = 1; ii < 5; ++ii) {
436 if (cnts[ii] > maxCnt) {
441 return codeList[codeIndex];
452 if (primID <= 0 || primID > (
int)slc.
tjs.size())
return -1;
455 auto& ptj = slc.
tjs[primID - 1];
456 for (
unsigned short end = 0;
end < 2; ++
end) {
457 if (ptj.VtxID[
end] == 0)
continue;
458 auto& vx2 = slc.
vtxs[ptj.VtxID[
end] - 1];
459 if (vx2.Vx3ID == 0)
continue;
460 auto& vx3 = slc.
vtx3s[vx2.Vx3ID - 1];
461 if (vx3.Neutrino)
return primID;
475 for (
unsigned short nit = 0; nit < 10; ++nit) {
476 if (parid < 1 || parid > (
int)slc.
tjs.size())
break;
477 auto& tj = slc.
tjs[parid - 1];
492 int dtruid = pfp.
UID;
493 unsigned short nit = 0;
496 auto& parent =
slices[slcIndx.first].pfps[slcIndx.second];
498 if (parent.PDGCode == 14 || parent.PDGCode == 12)
return dtruid;
500 if (parent.ParentUID == 0)
return parent.UID;
501 if (
int(parent.ParentUID) == parent.UID)
return parent.UID;
503 paruid = parent.ParentUID;
504 if (paruid < 0)
return 0;
506 if (nit == 10)
return 0;
514 if (mtjid > (
int)slc.
tjs.size())
return false;
515 auto const& mtj = slc.
tjs[mtjid - 1];
518 for (
auto tjid : pfp.
TjIDs) {
519 auto const& otj = slc.
tjs[tjid - 1];
520 if (otj.CTP == mtj.CTP) {
525 if (otjid == 0)
return false;
527 int newtjid = slc.
tjs.size();
529 mf::LogVerbatim(
"TC") <<
"MergeTjIntoPFP: merged T" << otjid <<
" with T" << mtjid
530 <<
" -> T" << newtjid;
531 std::replace(pfp.
TjIDs.begin(), pfp.
TjIDs.begin(), otjid, newtjid);
536 mf::LogVerbatim(
"TC") <<
"MergeTjIntoPFP: merge T" << otjid <<
" with T" << mtjid
548 if (tj.
AveChg <= 0)
return 100;
550 unsigned short closePt = USHRT_MAX;
552 for (
unsigned short ipt = tj.
EndPt[0]; ipt <= tj.
EndPt[1]; ++ipt) {
553 auto& tp = tj.
Pts[ipt];
554 float sep2 =
PosSep2(pos, tp.Pos);
555 if (sep2 > close)
continue;
559 if (closePt == USHRT_MAX)
return 100;
561 auto& tp = tj.
Pts[closePt];
564 float posErr = tp.DeltaRMS;
565 if (tp.AngErr > 0 && close > 10) posErr += sqrt(tp.AngErr * sqrt(close));
566 if (posErr < 0.1) posErr = 0.1;
567 float posPull = delta / posErr;
569 if (chgErr < 0.15) chgErr = 0.15;
572 return 0.5 * (posPull + chgPull);
590 if (tjIDs.size() < 2)
return false;
591 unsigned short lasttj = tjIDs[tjIDs.size() - 1] - 1;
592 auto const& mtj = slc.
tjs[lasttj];
593 bool mtjIsShort = (mtj.Pts.size() < 5);
595 std::array<float, 2> minsep2{{1000, 1000}};
597 std::array<int, 2> minsepTj{{0, 0}};
599 std::array<unsigned short, 2> minsepPt;
602 std::array<unsigned short, 2> minsepEnd;
603 for (
auto tjid : tjIDs) {
604 auto& tj = slc.
tjs[tjid - 1];
605 if (tj.CTP != mtj.CTP)
continue;
606 if (tj.ID == mtj.ID)
continue;
607 for (
unsigned short mend = 0; mend < 2; ++mend) {
608 Point2_t mendPos = mtj.Pts[mtj.EndPt[mend]].Pos;
609 float sep2 = minsep2[mend];
610 unsigned short closePt = 0;
612 minsep2[mend] = sep2;
613 minsepTj[mend] = tjid;
614 minsepPt[mend] = closePt;
617 short dend0 =
abs((
short)closePt - tj.EndPt[0]);
618 short dend1 =
abs((
short)closePt - tj.EndPt[1]);
619 if (dend0 < dend1 && dend0 < 3) minsepEnd[mend] = 0;
620 if (dend1 < dend0 && dend1 < 3) minsepEnd[mend] = 1;
626 bool isCompatible = (minsepEnd[0] != 2 && minsepEnd[1] != 2);
628 if (isCompatible && mtjIsShort) {
629 float minminsep = minsep2[0];
630 if (minsep2[1] < minminsep) minminsep = minsep2[1];
632 isCompatible = minminsep < 5;
636 myprt <<
"CompatibleMerge: T" << mtj.ID <<
" end";
637 for (
unsigned short end = 0;
end < 2; ++
end)
638 myprt <<
" T" << minsepTj[
end] <<
"_I" << minsepPt[
end] <<
"_E" << minsepEnd[
end]
639 <<
" minsep " << sqrt(minsep2[
end]);
640 myprt <<
" Compatible? " << isCompatible;
653 if (tj1.
CTP != tj2.
CTP)
return false;
654 unsigned short end1 = -1, end2 = 0;
657 if (len2 < minLen) minLen = len2;
659 if (minLen > 10) minLen = 10;
660 for (
unsigned short e1 = 0; e1 < 2; ++e1) {
662 for (
unsigned short e2 = 0; e2 < 2; ++e2) {
664 float sep =
PosSep(tp1.Pos, tp2.Pos);
673 if (end2 != 1 - end1)
return false;
676 if (overlapFraction > 0.25) {
679 << overlapFraction <<
" > 0.25 ";
683 auto& tp1 = tj1.
Pts[tj1.
EndPt[end1]];
684 auto& tp2 = tj2.
Pts[tj2.
EndPt[end2]];
687 if (doca1 > 2 && doca2 > 2) {
710 float maxWire = -1E6;
713 for (
auto& tp : tj1.
Pts) {
714 if (tp.Chg == 0)
continue;
715 if (tp.Pos[0] < 0)
continue;
716 if (tp.Pos[0] < minWire) minWire = tp.Pos[0];
717 if (tp.Pos[0] > maxWire) maxWire = tp.Pos[0];
720 if (cnt1 == 0)
return 0;
722 for (
auto& tp : tj2.
Pts) {
723 if (tp.Chg == 0)
continue;
724 if (tp.Pos[0] < 0)
continue;
725 if (tp.Pos[0] < minWire) minWire = tp.Pos[0];
726 if (tp.Pos[0] > maxWire) maxWire = tp.Pos[0];
729 if (cnt2 == 0)
return 0;
730 int span = maxWire - minWire;
731 if (span <= 0)
return 0;
732 std::vector<unsigned short> wcnt(span);
733 for (
auto& tp : tj1.
Pts) {
734 if (tp.Chg == 0)
continue;
735 if (tp.Pos[0] < -0.4)
continue;
736 int indx = std::nearbyint(tp.Pos[0] - minWire);
737 if (indx < 0 || indx > span - 1)
continue;
740 for (
auto& tp : tj2.
Pts) {
741 if (tp.Chg == 0)
continue;
742 if (tp.Pos[0] < -0.4)
continue;
743 int indx = std::nearbyint(tp.Pos[0] - minWire);
744 if (indx < 0 || indx > span - 1)
continue;
747 float cntOverlap = 0;
748 for (
auto cnt : wcnt)
749 if (cnt > 1) ++cntOverlap;
750 if (cnt1 < cnt2) {
return cntOverlap / cnt1; }
752 return cntOverlap / cnt2;
786 if (angle > M_PI) angle = M_PI;
787 if (angle < -M_PI) angle = M_PI;
788 if (angle < 0) angle = -angle;
789 if (angle > M_PI / 2) angle = M_PI - angle;
800 unsigned short originPt = tj.
EndPt[1];
801 unsigned short npts = tj.
Pts[originPt].NTPsFit;
803 unsigned short fitDir = -1;
804 FitTraj(slc, tj, originPt, npts, fitDir, tpFit);
805 tj.
Pts[originPt] = tpFit;
812 unsigned short originPt,
832 if (originPt > tj.
Pts.size() - 1) {
833 mf::LogWarning(
"TC") <<
"FitTraj: Requesting fit of invalid TP " << originPt;
838 tpFit = tj.
Pts[originPt];
841 if (fitDir < -1 || fitDir > 1)
return;
843 std::vector<double>
x,
y;
846 if (tj.
Pts[originPt].Chg == 0) origin = tj.
Pts[originPt].Pos;
850 for (
unsigned short ipt = tj.
EndPt[0]; ipt < tj.
EndPt[1]; ++ipt) {
851 if (tj.
Pts[ipt].Chg <= 0)
continue;
852 double xx = tj.
Pts[ipt].HitPos[0] - origin[0];
853 double yy = tj.
Pts[ipt].HitPos[1] - origin[1];
857 if (x.size() != 2)
return;
860 tpFit.
Ang = M_PI / 2;
861 if (y[1] < y[0]) tpFit.
Ang = -tpFit.
Ang;
864 double dx = x[1] - x[0];
865 double dy = y[1] - y[0];
866 tpFit.
Ang = atan2(dy, dx);
868 tpFit.
Dir[0] = cos(tpFit.
Ang);
869 tpFit.
Dir[1] = sin(tpFit.
Ang);
870 tpFit.
Pos[0] += origin[0];
871 tpFit.
Pos[1] += origin[1];
878 std::vector<double>
w, q;
879 std::array<double, 2>
dir;
880 double xx, yy, xr, yr;
885 double rotAngle = tj.
Pts[originPt].Ang;
886 double cs = cos(-rotAngle);
887 double sn = sin(-rotAngle);
890 if (tj.
Pts[originPt].Chg > 0) {
891 xx = tj.
Pts[originPt].HitPos[0] - origin[0];
892 yy = tj.
Pts[originPt].HitPos[1] - origin[1];
893 xr = cs * xx - sn * yy;
894 yr = sn * xx + cs * yy;
897 chgWt = tj.
Pts[originPt].ChgPull;
898 if (chgWt < 1) chgWt = 1;
900 w.push_back(chgWt * tj.
Pts[originPt].HitPosErr2);
904 if (fitDir != 0) --npts;
908 unsigned short cnt = 0;
909 for (
unsigned short ipt = originPt + 1; ipt < tj.
Pts.size(); ++ipt) {
910 if (tj.
Pts[ipt].Chg <= 0)
continue;
911 xx = tj.
Pts[ipt].HitPos[0] - origin[0];
912 yy = tj.
Pts[ipt].HitPos[1] - origin[1];
913 xr = cs * xx - sn * yy;
914 yr = sn * xx + cs * yy;
917 chgWt = tj.
Pts[ipt].ChgPull;
918 if (chgWt < 1) chgWt = 1;
920 w.push_back(chgWt * tj.
Pts[ipt].HitPosErr2);
922 if (cnt == npts)
break;
927 if (fitDir != 1 && originPt > 0) {
928 unsigned short cnt = 0;
929 for (
unsigned short ii = 1; ii < tj.
Pts.size(); ++ii) {
930 unsigned short ipt = originPt - ii;
931 if (ipt > tj.
Pts.size() - 1)
continue;
932 if (tj.
Pts[ipt].Chg == 0)
continue;
933 xx = tj.
Pts[ipt].HitPos[0] - origin[0];
934 yy = tj.
Pts[ipt].HitPos[1] - origin[1];
935 xr = cs * xx - sn * yy;
936 yr = sn * xx + cs * yy;
939 chgWt = tj.
Pts[ipt].ChgPull;
940 if (chgWt < 1) chgWt = 1;
942 w.push_back(chgWt * tj.
Pts[ipt].HitPosErr2);
944 if (cnt == npts)
break;
950 if (x.size() < 2)
return;
961 for (
unsigned short ipt = 0; ipt < x.size(); ++ipt) {
962 if (w[ipt] < 0.00001) w[ipt] = 0.00001;
965 sumx += wght * x[ipt];
966 sumy += wght * y[ipt];
967 sumx2 += wght * x[ipt] * x[ipt];
968 sumy2 += wght * y[ipt] * y[ipt];
969 sumxy += wght * x[ipt] * y[ipt];
972 double delta = sum * sumx2 - sumx * sumx;
973 if (delta == 0)
return;
975 double A = (sumx2 * sumy - sumx * sumxy) / delta;
977 double B = (sumxy * sum - sumx * sumy) / delta;
982 double newang = atan(B);
983 dir[0] = cos(newang);
984 dir[1] = sin(newang);
988 tpFit.
Dir[0] = cs * dir[0] - sn * dir[1];
989 tpFit.
Dir[1] = sn * dir[0] + cs * dir[1];
991 bool flipDir =
false;
993 flipDir = std::signbit(tpFit.
Dir[1]) != std::signbit(tj.
Pts[originPt].Dir[1]);
996 flipDir = std::signbit(tpFit.
Dir[0]) != std::signbit(tj.
Pts[originPt].Dir[0]);
999 tpFit.
Dir[0] = -tpFit.
Dir[0];
1000 tpFit.
Dir[1] = -tpFit.
Dir[1];
1002 tpFit.
Ang = atan2(tpFit.
Dir[1], tpFit.
Dir[0]);
1006 tpFit.
Pos[0] = -sn * A + origin[0];
1007 tpFit.
Pos[1] = cs * A + origin[1];
1011 if (x.size() < 3)
return;
1014 double ndof = x.size() - 2;
1016 (sumy2 + A * A * sum + B * B * sumx2 - 2 * (A * sumy + B * sumxy - A * B * sumx)) / ndof;
1020 double slopeError = sqrt(varnce * sum / delta);
1029 for (
unsigned short ii = 0; ii < y.size(); ++ii) {
1030 arg = y[ii] - A - B * x[ii];
1031 sum += arg * arg / w[ii];
1033 tpFit.
FitChi = sum / ndof;
1039 if (slc.
pfps.empty())
return USHRT_MAX;
1040 for (
unsigned int ipfp = 0; ipfp < slc.
pfps.size(); ++ipfp) {
1041 const auto& pfp = slc.
pfps[ipfp];
1042 if (std::find(pfp.TjIDs.begin(), pfp.TjIDs.end(), tjID) != pfp.TjIDs.end())
return ipfp;
1051 for (
auto const& tp : tj.
Pts) {
1052 for (
auto iht : tp.Hits) {
1063 for (
unsigned short ii = 0; ii < tp.
Hits.size(); ++ii) {
1076 for (
auto& tp : tj.
Pts) {
1077 if (tp.Hits.size() > 16)
return false;
1090 unsigned short npts = tj.
EndPt[1] - tj.
EndPt[0] + 1;
1091 if (npts < 2)
return false;
1093 auto& endTp0 = tj.
Pts[tj.
EndPt[0]];
1094 auto& endTp1 = tj.
Pts[tj.
EndPt[1]];
1098 if (endTp0.AngErr == 0.1 && endTp1.AngErr != 0.1) { endTp0.AngErr = endTp1.AngErr; }
1099 else if (endTp0.AngErr != 0.1 && endTp1.AngErr == 0.1) {
1100 endTp1.AngErr = endTp0.AngErr;
1106 if (endTp0.AveChg <= 0) {
1107 unsigned short cnt = 0;
1109 for (
unsigned short ipt = tj.
EndPt[0]; ipt <= tj.
EndPt[1]; ++ipt) {
1110 if (tj.
Pts[ipt].Chg == 0)
continue;
1111 sum += tj.
Pts[ipt].Chg;
1113 if (cnt == 4)
break;
1115 tj.
Pts[tj.
EndPt[0]].AveChg = sum / (float)cnt;
1117 if (endTp1.AveChg <= 0 && npts < 5) endTp1.AveChg = endTp0.AveChg;
1118 if (endTp1.AveChg <= 0) {
1120 unsigned short cnt = 0;
1121 for (
unsigned short ii = 0; ii < tj.
Pts.size(); ++ii) {
1122 short ipt = tj.
EndPt[1] - ii;
1124 if (tj.
Pts[ipt].Chg == 0)
continue;
1125 sum += tj.
Pts[ipt].Chg;
1127 if (cnt == 4)
break;
1128 if (ipt == 0)
break;
1130 tj.
Pts[tj.
EndPt[1]].AveChg = sum / (float)cnt;
1137 if (npts > 2 * nPtsFit) {
1138 for (
unsigned short ipt = tj.
EndPt[0] + nPtsFit; ipt < tj.
EndPt[1] - nPtsFit; ++ipt) {
1139 auto& tp = tj.
Pts[ipt];
1140 if (tp.KinkSig < 0) tp.KinkSig =
KinkSignificance(slc, tj, ipt, nPtsFit, useChg,
false);
1147 int trID = slc.
tjs.size() + 1;
1149 for (
unsigned short ipt = tj.
EndPt[0]; ipt <= tj.
EndPt[1]; ++ipt) {
1150 for (
unsigned short ii = 0; ii < tj.
Pts[ipt].Hits.size(); ++ii) {
1151 if (tj.
Pts[ipt].UseHit[ii]) {
1152 unsigned int iht = tj.
Pts[ipt].Hits[ii];
1153 if (iht > slc.
slHits.size() - 1) {
1157 if (slc.
slHits[iht].InTraj > 0) {
1161 slc.
slHits[iht].InTraj = trID;
1167 for (
unsigned int iht = 0; iht < slc.
slHits.size(); ++iht) {
1168 if (slc.
slHits[iht].InTraj == tj.
ID) {
1170 <<
" thinks it belongs to T" << tj.
ID <<
" but it isn't in the Tj\n";
1182 slc.
tjs.push_back(tj);
1185 for (
unsigned short ipt = 0; ipt < tj.
Pts.size(); ++ipt) {
1186 for (
unsigned short ii = 0; ii < tj.
Pts[ipt].Hits.size(); ++ii) {
1187 unsigned int iht = tj.
Pts[ipt].Hits[ii];
1189 std::cout <<
"Debug hit appears in trajectory w WorkID " << tj.
WorkID <<
" UseHit " 1190 << tj.
Pts[ipt].UseHit[ii] <<
"\n";
1202 unsigned short originPt,
1203 unsigned short npts,
1206 unsigned short usePar)
1213 if (originPt > tj.
Pts.size() - 1)
return;
1214 if (fitDir != 1 && fitDir != -1)
return;
1218 Fit2D(0, inPt, pErr, outVec, outVecErr, chiDOF);
1219 unsigned short cnt = 0;
1220 for (
unsigned short ii = 0; ii < tj.
Pts.size(); ++ii) {
1221 unsigned short ipt = originPt + ii * fitDir;
1222 if (ipt < tj.
EndPt[0] || ipt > tj.
EndPt[1])
break;
1223 auto& tp = tj.
Pts[ipt];
1224 if (tp.Chg <= 0)
continue;
1226 inPt[0] =
std::abs(tp.Pos[0] - tj.
Pts[originPt].Pos[0]);
1227 float parVal = tp.Chg;
1229 pErr = 0.1 * parVal;
1233 pErr = sqrt(tp.HitPosErr2);
1237 if (!
Fit2D(2, inPt, pErr, outVec, outVecErr, chiDOF))
break;
1239 if (cnt == npts)
break;
1241 if (cnt < npts)
return;
1243 if (!
Fit2D(-1, inPt, pErr, outVec, outVecErr, chiDOF))
return;
1244 pFit.
Pos = tj.
Pts[originPt].Pos;
1245 pFit.
Par0 = outVec[0];
1246 pFit.
AvePar /= (float)cnt;
1247 pFit.
ParErr = outVecErr[0];
1248 pFit.
Pos = tj.
Pts[originPt].Pos;
1262 unsigned short itj = 0;
1263 std::vector<unsigned int> tHits;
1264 std::vector<unsigned int> atHits;
1265 for (
auto& tj : slc.
tjs) {
1267 if (tj.AlgMod[
kKilled])
continue;
1270 if (tHits.size() < 2)
continue;
1271 std::sort(tHits.begin(), tHits.end());
1273 for (iht = 0; iht < slc.
slHits.size(); ++iht) {
1274 if (slc.
slHits[iht].InTraj == tID) atHits.push_back(iht);
1276 if (atHits.size() < 2)
continue;
1277 if (!std::equal(tHits.begin(), tHits.end(), atHits.begin())) {
1279 myprt << someText <<
" ChkInTraj failed: inTraj - UseHit mis-match for T" << tID
1280 <<
" tj.WorkID " << tj.WorkID <<
" atHits size " << atHits.size() <<
" tHits size " 1281 << tHits.size() <<
" in CTP " << tj.CTP <<
"\n";
1282 myprt <<
"AlgMods: ";
1283 for (
unsigned short ib = 0; ib <
AlgBitNames.size(); ++ib)
1284 if (tj.AlgMod[ib]) myprt <<
" " <<
AlgBitNames[ib];
1286 myprt <<
"index inTraj UseHit \n";
1287 for (iht = 0; iht < atHits.size(); ++iht) {
1288 myprt <<
"iht " << iht <<
" " <<
PrintHit(slc.
slHits[atHits[iht]]);
1289 if (iht < tHits.size()) myprt <<
" " <<
PrintHit(slc.
slHits[tHits[iht]]);
1290 if (atHits[iht] != tHits[iht]) myprt <<
" <<< " << atHits[iht] <<
" != " << tHits[iht];
1293 if (tHits.size() > atHits.size()) {
1294 for (iht = atHits.size(); iht < atHits.size(); ++iht) {
1295 myprt <<
"atHits " << iht <<
" " <<
PrintHit(slc.
slHits[atHits[iht]]) <<
"\n";
1302 for (
unsigned short end = 0;
end < 2; ++
end) {
1303 if (tj.VtxID[
end] > slc.
vtxs.size()) {
1304 mf::LogVerbatim(
"TC") << someText <<
" ChkInTraj: Bad VtxID " << tj.ID;
1328 if (itj > slc.
tjs.size() - 1)
return;
1329 auto& tj = slc.
tjs[itj];
1332 if (tj.EndFlag[0][
kBragg])
return;
1335 if (tj.Pts.size() < 20)
return;
1340 float chg2 = tj.Pts[tj.EndPt[0] + 2].AveChg;
1342 float chg15 = tj.Pts[tj.EndPt[0] + 15].AveChg;
1343 if (chg2 < 3 * chg15)
return;
1346 float midChg = 0.5 * (chg2 + chg15);
1348 unsigned short breakPt = USHRT_MAX;
1349 for (
unsigned short ipt = tj.EndPt[0] + 3; ipt < 15; ++ipt) {
1350 float chgm2 = tj.Pts[ipt - 2].Chg;
1351 if (chgm2 == 0)
continue;
1352 float chgm1 = tj.Pts[ipt - 1].Chg;
1353 if (chgm1 == 0)
continue;
1354 float chgp1 = tj.Pts[ipt + 1].Chg;
1355 if (chgp1 == 0)
continue;
1356 float chgp2 = tj.Pts[ipt + 2].Chg;
1357 if (chgp2 == 0)
continue;
1358 if (chgm2 > midChg && chgm1 > midChg && chgp1 < midChg && chgp2 < midChg) {
1363 if (breakPt == USHRT_MAX)
return;
1365 std::array<double, 2> cnt,
sum, sum2;
1366 for (
unsigned short ipt = tj.EndPt[0]; ipt <= tj.EndPt[1]; ++ipt) {
1367 auto& tp = tj.Pts[ipt];
1368 if (tp.Chg <= 0)
continue;
1369 unsigned short end = 0;
1370 if (ipt > breakPt) end = 1;
1373 sum2[
end] += tp.Chg * tp.Chg;
1375 for (
unsigned short end = 0;
end < 2; ++
end) {
1376 if (cnt[
end] < 3)
return;
1377 double ave = sum[
end] / cnt[
end];
1378 double arg = sum2[
end] - cnt[
end] * ave * ave;
1379 if (arg <= 0)
return;
1380 sum2[
end] = sqrt(arg / (cnt[
end] - 1));
1384 bool doSplit =
true;
1387 if (tj.ChgRMS > 0.5 && sum2[0] > 0.3 && sum2[1] > 0.3) doSplit =
false;
1390 myprt <<
"CTBC: T" << tj.ID <<
" chgRMS " << tj.ChgRMS;
1391 myprt <<
" AveChg before split point " << (int)sum[0] <<
" rms " << sum2[0];
1392 myprt <<
" after " << (int)sum[1] <<
" rms " << sum2[1] <<
" doSplit? " << doSplit;
1394 if (!doSplit)
return;
1397 aVtx.
Pos = tj.Pts[breakPt].Pos;
1399 aVtx.
Pass = tj.Pass;
1403 aVtx.
ID = slc.
vtxs.size() + 1;
1405 unsigned short ivx = slc.
vtxs.size();
1407 if (!
SplitTraj(slc, itj, breakPt, ivx, prt)) {
1426 if (itj > slc.
tjs.size() - 1)
return false;
1430 auto& tj = slc.
tjs[itj];
1432 if (npwc < 4)
return false;
1433 if (npwc < nPtsToCheck) nPtsToCheck = npwc;
1436 unsigned short maxPullPt = USHRT_MAX;
1437 for (
unsigned short ipt = tj.EndPt[0]; ipt < tj.EndPt[1]; ++ipt) {
1438 auto& tp = tj.Pts[ipt];
1439 if (tp.ChgPull < maxPull)
continue;
1440 maxPull = tp.ChgPull;
1443 if (maxPullPt == USHRT_MAX)
return false;
1445 if (maxPullPt < 0.5 * (tj.EndPt[0] + tj.EndPt[1])) { dpt = maxPullPt - tj.EndPt[0]; }
1447 dpt = tj.EndPt[1] - maxPullPt;
1449 if (dpt < 3)
return false;
1452 mf::LogVerbatim(
"TC") <<
"BS: T" << tj.ID <<
" maxPull " << maxPull <<
" at " 1453 <<
PrintPos(tj.Pts[maxPullPt]) <<
" dpt " << dpt;
1454 unsigned short breakPt = USHRT_MAX;
1456 unsigned short bestBragg = 0;
1460 for (
unsigned short ipt = maxPullPt - 2; ipt <= maxPullPt + 2; ++ipt) {
1461 FitTraj(slc, tj, ipt - 1, nPtsFit, -1, tp1);
1462 if (tp1.
FitChi > 10)
continue;
1463 FitTraj(slc, tj, ipt + 1, nPtsFit, 1, tp2);
1464 if (tp2.
FitChi > 10)
continue;
1466 FitPar(tj, ipt - 1, nPtsToCheck, -1, chgFit1, 1);
1467 if (chgFit1.
ChiDOF > 100)
continue;
1469 FitPar(tj, ipt + 1, nPtsToCheck, 1, chgFit2, 1);
1470 if (chgFit2.
ChiDOF > 100)
continue;
1475 unsigned short bragg = 1;
1476 float bchi = chgFit1.
ChiDOF;
1483 if (bchi < 1) bchi = 1;
1484 float fom = 10 * dang * chgAsym * slpAsym / bchi;
1487 myprt <<
"pt " <<
PrintPos(tj.Pts[ipt]) <<
" " << std::setprecision(2) << dang;
1488 myprt <<
" chg1 " << (int)chgFit1.
Par0 <<
" slp " << chgFit1.
ParSlp <<
" chi " 1490 myprt <<
" chg2 " << (
int)chgFit2.
Par0 <<
" slp " << chgFit2.
ParSlp <<
" chi " 1492 myprt <<
" chgAsym " << chgAsym;
1493 myprt <<
" slpAsym " << slpAsym;
1494 myprt <<
" fom " << fom;
1495 myprt <<
" bragg " << bragg;
1497 if (fom < bestFOM)
continue;
1502 if (breakPt == USHRT_MAX)
return false;
1507 aVtx.
Pos = tj.Pts[breakPt].Pos;
1509 aVtx.
Pass = tj.Pass;
1513 aVtx.
ID = slc.
vtxs.size() + 1;
1515 unsigned short ivx = slc.
vtxs.size();
1517 if (!
SplitTraj(slc, itj, breakPt, ivx, prt)) {
1524 unsigned short otj = slc.
tjs.size() - 1;
1525 if (bestBragg == 2) std::swap(itj, otj);
1526 slc.
tjs[itj].PDGCode = 211;
1528 slc.
tjs[otj].PDGCode = 13;
1538 if (tj.
PDGCode == 111)
return;
1541 if (npwc < 50)
return;
1547 if (nPtsMax > 8) nPtsMax = 8;
1552 unsigned short firstBad = USHRT_MAX;
1553 for (
unsigned short ii = 0; ii < nPtsMax; ++ii) {
1554 unsigned short ipt = tj.
EndPt[1] - nPtsMax + ii;
1555 auto& tp = tj.
Pts[ipt];
1556 if (tp.Chg <= 0)
continue;
1557 if (tp.ChgPull < 3)
continue;
1559 if (firstBad == USHRT_MAX) firstBad = ipt;
1561 if (firstBad == USHRT_MAX)
return;
1563 float cntTot = tj.
EndPt[1] - firstBad;
1565 float fracBad = cntBad / cntTot;
1566 if (fracBad < 0.5)
return;
1569 for (
unsigned short ipt = firstBad; ipt <= tj.
EndPt[1]; ++ipt)
1578 const std::vector<float>& fQualityCuts,
1591 if (tj.
PDGCode == 111)
return;
1595 short minPts = fQualityCuts[1];
1596 if (minPts < 1)
return;
1597 if (npwc < minPts)
return;
1599 if (npwc < 8)
return;
1602 if (npwc == minPts + 1) {
1603 unsigned short endPt1 = tj.
EndPt[1];
1604 auto& tp = tj.
Pts[endPt1];
1605 auto& ptp = tj.
Pts[endPt1 - 1];
1608 float dwire =
std::abs(ptp.Pos[0] - tp.Pos[0]);
1609 if (ptp.Chg == 0 || dwire > 1.1) {
1619 for (lastPt = tj.
EndPt[1]; lastPt >= minPts; --lastPt) {
1621 if (lastPt == 1)
break;
1622 if (tj.
Pts[lastPt].Chg == 0)
continue;
1624 unsigned short nadj = 0;
1625 unsigned short npwc = 0;
1626 for (
short ipt = lastPt - minPts; ipt < lastPt; ++ipt) {
1629 auto& tp = tj.
Pts[ipt];
1631 auto& ptp = tj.
Pts[ipt - 1];
1632 if (tp.Chg > 0 && ptp.Chg > 0) {
1634 if (
std::abs(tp.Pos[0] - ptp.Pos[0]) < 1.5) ++nadj;
1639 float hitFrac = ntpwc / nwires;
1641 mf::LogVerbatim(
"TC") << fcnLabel <<
"-TEP: T" << tj.
ID <<
" lastPt " << lastPt <<
" npwc " 1642 << npwc <<
" ntpwc " << ntpwc <<
" nadj " << nadj <<
" hitFrac " 1644 if (hitFrac > fQualityCuts[0] && npwc == minPts && nadj >= minPts - 1)
break;
1649 if (tj.
Pts[lastPt].Pos[0] > -0.4) {
1650 unsigned int prevWire = std::nearbyint(tj.
Pts[lastPt].Pos[0]);
1651 if (tj.
StepDir > 0) { --prevWire; }
1656 mf::LogVerbatim(
"TC") << fcnLabel <<
"-TEP: is prevWire " << prevWire <<
" dead? ";
1663 if (lastPt == tj.
EndPt[1]) {
1669 for (
unsigned short ipt = lastPt + 1; ipt <= tj.
EndPt[1]; ++ipt)
1674 fcnLabel +=
"-TEPo";
1685 if (tj.
PDGCode == 111)
return;
1691 unsigned short withNptsFit = 0;
1694 for (
unsigned short nptsf = 3; nptsf < nPtsFit; ++nptsf) {
1695 unsigned short ipt = tj.
EndPt[1] - nptsf;
1699 withNptsFit = nptsf;
1702 if (withNptsFit > 0) {
1703 unsigned short ipt = tj.
EndPt[1] - withNptsFit;
1704 std::cout <<
"CEK: T" << tj.
ID <<
" ipt " << ipt;
1706 auto& tp = tj.
Pts[ipt];
1707 std::cout <<
" " <<
PrintPos(tp) <<
" withNptsFit " << withNptsFit <<
" ks " << ks <<
"\n";
1720 if (tj.
PDGCode == 111)
return;
1721 unsigned short npts = tj.
EndPt[1] - tj.
EndPt[0];
1724 if (npts > 50)
return;
1726 if (npts < 8)
return;
1729 unsigned short atPt = 0;
1731 for (
unsigned short ipt = tj.
EndPt[0] + 2; ipt <= tj.
EndPt[1] - 2; ++ipt) {
1732 auto& tp = tj.
Pts[ipt];
1733 if (tp.ChgPull > bigPull) {
1734 bigPull = tp.ChgPull;
1738 if (atPt == 0)
return;
1740 if ((atPt - tj.
EndPt[0]) < 0.5 * npts)
return;
1743 <<
". Check charge asymmetry around it.";
1744 unsigned short nchk = 0;
1745 unsigned short npos = 0;
1746 unsigned short nneg = 0;
1747 for (
short ii = 1; ii < 5; ++ii) {
1748 short iplu = atPt + ii;
1749 if (iplu > tj.
EndPt[1])
break;
1750 short ineg = atPt - ii;
1751 if (ineg < tj.
EndPt[0])
break;
1752 if (tj.
Pts[iplu].Chg == 0)
continue;
1753 if (tj.
Pts[ineg].Chg == 0)
continue;
1754 float asym = (tj.
Pts[iplu].Chg - tj.
Pts[ineg].Chg) / (tj.
Pts[iplu].Chg + tj.
Pts[ineg].Chg);
1756 if (asym > 0.5) ++npos;
1757 if (asym < -0.5) ++nneg;
1759 mf::LogVerbatim(
"TC") <<
" ineg " << ineg <<
" iplu " << iplu <<
" asym " << asym
1760 <<
" nchk " << nchk;
1762 if (nchk < 3)
return;
1765 bool doTrim = (nneg > nchk) || (npos > nchk);
1766 if (!doTrim)
return;
1768 auto& prevTP = tj.
Pts[atPt - 1];
1769 if (
std::abs(prevTP.ChgPull) > 2) --atPt;
1770 for (
unsigned short ipt = atPt; ipt <= tj.
EndPt[1]; ++ipt)
1781 if (MinWireSignalFraction == 0)
return true;
1783 if (tp1.
Pos[0] < -0.4 || tp2.
Pos[0] < -0.4)
return false;
1784 int fromWire = std::nearbyint(tp1.
Pos[0]);
1785 int toWire = std::nearbyint(tp2.
Pos[0]);
1787 if (fromWire == toWire) {
1790 tp.
Pos[1] = 0.5 * (tp1.
Pos[1] + tp2.
Pos[1]);
1812 if (tp.
Pos[0] < -0.4 || toPos0 < -0.4)
return 0;
1813 int fromWire = std::nearbyint(tp.
Pos[0]);
1814 int toWire = std::nearbyint(toPos0);
1816 if (fromWire == toWire)
return SignalAtTp(tp);
1818 int nWires =
abs(toWire - fromWire) + 1;
1823 if (toWire > fromWire && tp.
Dir[0] < 0) stepSize = -stepSize;
1824 if (toWire < fromWire && tp.
Dir[0] > 0) stepSize = -stepSize;
1827 for (
unsigned short cnt = 0; cnt < nWires; ++cnt) {
1830 tp.
Pos[0] += tp.
Dir[0] * stepSize;
1831 tp.
Pos[1] += tp.
Dir[1] * stepSize;
1833 float sigFrac = nsig / num;
1839 const std::vector<unsigned int>& iHitsInMultiplet,
1840 const std::vector<unsigned int>& jHitsInMultiplet)
1844 if (iHitsInMultiplet.empty() || jHitsInMultiplet.empty())
return false;
1848 if (cvI < 0)
return false;
1851 for (
auto& iht : iHitsInMultiplet) {
1853 float cv =
hit.PeakTime();
1854 float rms =
hit.RMS();
1855 float arg = cv - 3.1 * rms;
1856 if (arg < minI) minI = arg;
1857 arg = cv + 3.1 * rms;
1858 if (arg > maxI) maxI = arg;
1862 if (cvJ < 0)
return false;
1865 for (
auto& jht : jHitsInMultiplet) {
1867 float cv =
hit.PeakTime();
1868 float rms =
hit.RMS();
1869 float arg = cv - 3.1 * rms;
1870 if (arg < minJ) minJ = arg;
1871 arg = cv + 3.1 * rms;
1872 if (arg > maxJ) maxJ = arg;
1876 if (maxI > minJ)
return true;
1879 if (minI < maxJ)
return true;
1888 if (iht > slc.
slHits.size() - 1)
return false;
1889 if (jht > slc.
slHits.size() - 1)
return false;
1893 int iwire = ihit.WireID().Wire;
1894 int jwire = jhit.WireID().Wire;
1895 if (
std::abs(iwire - jwire) > 1)
return false;
1896 if (ihit.PeakTime() > jhit.PeakTime()) {
1897 float minISignal = ihit.PeakTime() - 3 * ihit.RMS();
1898 float maxJSignal = jhit.PeakTime() + 3 * ihit.RMS();
1899 if (maxJSignal > minISignal)
return true;
1902 float maxISignal = ihit.PeakTime() + 3 * ihit.RMS();
1903 float minJSignal = jhit.PeakTime() - 3 * ihit.RMS();
1904 if (minJSignal > maxISignal)
return true;
1926 if (tp.
Pos[0] < -0.4)
return false;
1928 unsigned short pln = planeID.
Plane;
1929 unsigned int wire = std::nearbyint(tp.
Pos[0]);
1930 if (wire >
evt.
goodWire[pln].size() - 1)
return false;
1934 if (slc.
wireHitRange[pln][wire].first == UINT_MAX)
return false;
1937 float tickRange = 0;
1941 if (tickRange > 40) tickRange = 40;
1943 float loTpTick = projTick - tickRange;
1944 float hiTpTick = projTick + tickRange;
1945 for (
unsigned int iht = slc.
wireHitRange[pln][wire].first;
1948 unsigned int ahi = slc.
slHits[iht].allHitsIndex;
1950 if (projTick <
hit.PeakTime()) {
1951 float loHitTick =
hit.PeakTime() - 3 *
hit.RMS();
1952 if (hiTpTick > loHitTick)
return true;
1955 float hiHitTick =
hit.PeakTime() + 3 *
hit.RMS();
1956 if (loTpTick < hiHitTick)
return true;
1976 if (tp.
Pos[0] < -0.4)
return false;
1978 unsigned short pln = planeID.
Plane;
1979 unsigned int wire = std::nearbyint(tp.
Pos[0]);
1980 if (wire >
evt.
goodWire[pln].size() - 1)
return false;
1986 float tickRange = 0;
1990 if (tickRange > 40) tickRange = 40;
1992 float loTpTick = projTick - tickRange;
1993 float hiTpTick = projTick + tickRange;
2003 const auto& wid =
hit.WireID();
2004 if (wid.Cryostat != planeID.
Cryostat)
continue;
2005 if (wid.TPC != planeID.
TPC)
continue;
2006 if (wid.Plane != planeID.
Plane)
continue;
2007 if (projTick <
hit.PeakTime()) {
2008 float loHitTick =
hit.PeakTime() - 3 *
hit.RMS();
2009 if (hiTpTick > loHitTick)
return true;
2012 float hiHitTick =
hit.PeakTime() + 3 *
hit.RMS();
2013 if (loTpTick < hiHitTick)
return true;
2034 unsigned int pln = plnID.
Plane;
2035 if (pln == 2)
return false;
2037 unsigned int tpc = plnID.
TPC;
2042 float atTick = 0.5 * (loTick + hiTick);
2046 if (
hit.Channel() != chan)
continue;
2047 if (atTick <
hit.PeakTime()) {
2048 float loHitTick =
hit.PeakTime() - 3 *
hit.RMS();
2049 if (hiTick > loHitTick)
return true;
2052 float hiHitTick =
hit.PeakTime() + 3 *
hit.RMS();
2053 if (loTick < hiHitTick)
return true;
2063 for (
size_t i = 0; i < tp.
Hits.size(); ++i) {
2064 if (!tp.
UseHit[i])
continue;
2073 unsigned short firstPt = tj.
EndPt[0];
2074 unsigned short lastPt = tj.
EndPt[1];
2081 bool includeDeadWires,
2082 unsigned short firstPt,
2083 unsigned short lastPt)
2085 unsigned short ntp = 0;
2086 for (
unsigned short ipt = firstPt; ipt <= lastPt; ++ipt)
2087 if (tj.
Pts[ipt].Chg > 0) ++ntp;
2101 const float& inWirePos1,
2102 const float& inWirePos2,
2105 if (inWirePos1 < -0.4 || inWirePos2 < -0.4)
return 0;
2106 unsigned int inWire1 = std::nearbyint(inWirePos1);
2107 unsigned int inWire2 = std::nearbyint(inWirePos2);
2109 unsigned short plane = planeID.
Plane;
2110 if (inWire1 > slc.
nWires[plane] || inWire2 > slc.
nWires[plane])
return 0;
2111 if (inWire1 > inWire2) {
2113 unsigned int tmp = inWire1;
2118 unsigned int wire, ndead = 0;
2119 for (wire = inWire1; wire < inWire2; ++wire)
2127 unsigned short pdg =
abs(PDGCode);
2128 if (pdg == 11)
return 0;
2129 if (pdg == 13)
return 1;
2130 if (pdg == 211)
return 2;
2131 if (pdg == 321)
return 3;
2132 if (pdg == 2212)
return 4;
2141 if (itj > slc.
tjs.size() - 1)
return;
2142 int killTjID = slc.
tjs[itj].ID;
2144 if (
hit.InTraj == killTjID)
hit.InTraj = 0;
2151 if (itj > slc.
tjs.size() - 1)
return;
2154 <<
"RestoreObsoleteTrajectory: Trying to restore not-obsolete trajectory " 2159 for (
auto& tp : slc.
tjs[itj].Pts) {
2160 for (
unsigned short ii = 0; ii < tp.Hits.size(); ++ii) {
2161 if (tp.UseHit[ii]) {
2163 if (slc.
slHits[iht].InTraj == 0) { slc.
slHits[iht].InTraj = slc.
tjs[itj].ID; }
2176 for (
auto& shortTj : slc.
tjs) {
2177 if (shortTj.AlgMod[
kKilled] || shortTj.AlgMod[
kHaloTj])
continue;
2178 if (shortTj.CTP != inCTP)
continue;
2179 unsigned short spts = shortTj.EndPt[1] - shortTj.EndPt[0];
2180 if (spts > 20)
continue;
2182 if (shortTj.PDGCode == 11)
continue;
2184 if (shortTj.SSID > 0)
continue;
2186 if (tjhits.empty())
continue;
2187 std::vector<int> tids;
2188 std::vector<unsigned short> tcnt;
2189 for (
auto iht : tjhits) {
2191 if (
hit.InTraj <= 0)
continue;
2192 if ((
unsigned int)
hit.InTraj > slc.
tjs.size())
continue;
2193 if (
hit.InTraj == shortTj.ID)
continue;
2194 unsigned short indx = 0;
2195 for (indx = 0; indx < tids.size(); ++indx)
2196 if (
hit.InTraj == tids[indx])
break;
2197 if (indx == tids.size()) {
2198 tids.push_back(
hit.InTraj);
2205 if (tids.empty())
continue;
2207 unsigned short maxcnt = 0;
2208 for (
unsigned short indx = 0; indx < tids.size(); ++indx) {
2209 if (tcnt[indx] > maxcnt) {
2210 auto& ltj = slc.
tjs[tids[indx] - 1];
2211 unsigned short lpts = ltj.EndPt[1] - ltj.EndPt[0];
2212 if (lpts < spts)
continue;
2213 maxcnt = tcnt[indx];
2216 float hitFrac = (float)maxcnt / (
float)tjhits.size();
2217 if (hitFrac < 0.1)
continue;
2231 if (itj > slc.
tjs.size() - 1)
return false;
2233 auto& tj = slc.
tjs[itj];
2236 unsigned short atPt = USHRT_MAX;
2237 for (
unsigned short ipt = tj.EndPt[0] + 1; ipt <= tj.EndPt[1]; ++ipt) {
2238 if (tj.Pts[ipt].Pos[1] > tj.Pts[ipt - 1].Pos[1]) {
2240 if (tj.Pts[ipt - 1].Pos[1] < atPos1 && tj.Pts[ipt].Pos[1] >= atPos1) {
2247 if (tj.Pts[ipt - 1].Pos[1] >= atPos1 && tj.Pts[ipt].Pos[1] < atPos1) {
2253 if (atPt == USHRT_MAX)
return false;
2254 unsigned short vx2Index = USHRT_MAX;
2257 newVx2.
CTP = tj.CTP;
2258 newVx2.
Pos[0] = 0.5 * (tj.Pts[atPt - 1].Pos[0] + tj.Pts[atPt].Pos[0]);
2259 newVx2.
Pos[1] = 0.5 * (tj.Pts[atPt - 1].Pos[1] + tj.Pts[atPt].Pos[1]);
2264 return SplitTraj(slc, itj, atPt, vx2Index, prt);
2268 bool SplitTraj(
TCSlice& slc,
unsigned short itj,
unsigned short pos,
unsigned short ivx,
bool prt)
2275 if (itj > slc.
tjs.size() - 1)
return false;
2276 if (pos < slc.
tjs[itj].EndPt[0] + 1 || pos > slc.
tjs[itj].EndPt[1] - 1)
return false;
2277 if (ivx != USHRT_MAX && ivx > slc.
vtxs.size() - 1)
return false;
2282 bool splittingMuon = (tj.
PDGCode == 13);
2283 if (splittingMuon) tj.
PDGCode = 0;
2287 myprt <<
"SplitTraj: Split T" << tj.
ID <<
" at point " <<
PrintPos(tj.
Pts[pos]);
2288 if (ivx < slc.
vtxs.size()) myprt <<
" with Vtx 2V" << slc.
vtxs[ivx].ID;
2292 unsigned short ntp = 0;
2293 for (
unsigned short ipt = 0; ipt <= pos; ++ipt) {
2294 if (tj.
Pts[ipt].Chg > 0) ++ntp;
2298 if (prt)
mf::LogVerbatim(
"TC") <<
" Split point to small at begin " << ntp <<
" pos " << pos;
2302 for (
unsigned short ipt = pos + 1; ipt <= tj.
EndPt[1]; ++ipt) {
2303 if (tj.
Pts[ipt].Chg > 0) ++ntp;
2308 mf::LogVerbatim(
"TC") <<
" Split point too small at end " << ntp <<
" pos " << pos
2309 <<
" EndPt " << tj.
EndPt[1];
2315 newTj.
ID = slc.
tjs.size() + 1;
2324 for (
unsigned short ipt = pos + 1; ipt <= tj.
EndPt[1]; ++ipt) {
2325 tj.
Pts[ipt].Chg = 0;
2326 for (
unsigned short ii = 0; ii < tj.
Pts[ipt].Hits.size(); ++ii) {
2327 if (!tj.
Pts[ipt].UseHit[ii])
continue;
2328 iht = tj.
Pts[ipt].Hits[ii];
2330 if (slc.
slHits[iht].InTraj != tj.
ID)
continue;
2332 tj.
Pts[ipt].UseHit[ii] =
false;
2344 unsigned short eraseSize = pos - 2;
2345 if (eraseSize > newTj.
Pts.size() - 1) {
2358 newTj.
Pts.erase(newTj.
Pts.begin(), newTj.
Pts.begin() + eraseSize);
2360 for (
unsigned short ipt = 0; ipt < 3; ++ipt) {
2361 for (
unsigned short ii = 0; ii < newTj.
Pts[ipt].Hits.size(); ++ii)
2362 newTj.
Pts[ipt].UseHit[ii] =
false;
2363 newTj.
Pts[ipt].Chg = 0;
2369 if (ivx < slc.
vtxs.size()) newTj.
VtxID[0] = slc.
vtxs[ivx].ID;
2372 slc.
tjs.push_back(newTj);
2385 unsigned short& closePt,
2389 float best = minSep * minSep;
2390 closePt = USHRT_MAX;
2393 for (ipt = tj.
EndPt[0]; ipt <= tj.
EndPt[1]; ++ipt) {
2394 dw = tj.
Pts[ipt].Pos[0] - tp.
Pos[0];
2395 dt = tj.
Pts[ipt].Pos[1] - tp.
Pos[1];
2396 dp2 = dw * dw + dt * dt;
2402 minSep = sqrt(best);
2409 unsigned short& ipt1,
2410 unsigned short& ipt2,
2413 return TrajTrajDOCA(slc, tj1, tj2, ipt1, ipt2, minSep,
false);
2420 unsigned short& ipt1,
2421 unsigned short& ipt2,
2423 bool considerDeadWires)
2428 for (
unsigned short iwt = 0; iwt < 2; ++iwt) {
2431 float wt0 = tj1.
Pts[tj1.
EndPt[0]].Pos[iwt];
2432 float wt1 = tj1.
Pts[tj1.
EndPt[1]].Pos[iwt];
2440 wt0 = tj2.
Pts[tj2.
EndPt[0]].Pos[iwt];
2441 wt1 = tj2.
Pts[tj2.
EndPt[1]].Pos[iwt];
2451 if (lowt2 > hiwt1 + minSep)
return false;
2453 if (lowt1 > hiwt2 + minSep)
return false;
2456 float best = minSep * minSep;
2460 bool isClose =
false;
2461 for (
unsigned short i1 = tj1.
EndPt[0]; i1 < tj1.
EndPt[1] + 1; ++i1) {
2462 for (
unsigned short i2 = tj2.
EndPt[0]; i2 < tj2.
EndPt[1] + 1; ++i2) {
2464 float dw = tj1.
Pts[i1].Pos[0] - tj2.
Pts[i2].Pos[0] - dwc;
2465 if (
std::abs(dw) > minSep)
continue;
2466 float dt = tj1.
Pts[i1].Pos[1] - tj2.
Pts[i2].Pos[1];
2467 if (
std::abs(dt) > minSep)
continue;
2468 float dp2 = dw * dw + dt * dt;
2477 minSep = sqrt(best);
2485 if (iht > slc.
slHits.size() - 1 || jht > slc.
slHits.size() - 1)
return 1E6;
2488 float dw = (float)ihit.WireID().Wire - (float)jhit.WireID().Wire;
2490 return dw * dw + dt * dt;
2496 unsigned short endPt = tj.
EndPt[0];
2497 auto& tp0 = tj.
Pts[endPt];
2498 endPt = tj.
EndPt[1];
2499 auto& tp1 = tj.
Pts[endPt];
2507 float dw = wire - tp.
Pos[0];
2508 float dt = time - tp.
Pos[1];
2509 return dw * dw + dt * dt;
2515 if (iht > slc.
slHits.size() - 1)
return 1E6;
2517 float wire =
hit.WireID().Wire;
2534 double t = (double)(wire - tp.
Pos[0]) * tp.
Dir[0] + (double)(time - tp.
Pos[1]) * tp.
Dir[1];
2535 double dw = tp.
Pos[0] + t * tp.
Dir[0] - wire;
2536 double dt = tp.
Pos[1] + t * tp.
Dir[1] - time;
2537 return (
float)(dw * dw + dt * dt);
2554 double arg1 = tp1.
Pos[0] * tp1.
Dir[1] - tp1.
Pos[1] * tp1.
Dir[0];
2555 double arg2 = tp2.
Pos[0] * tp1.
Dir[1] - tp2.
Pos[1] * tp1.
Dir[0];
2556 double arg3 = tp2.
Dir[0] * tp1.
Dir[1] - tp2.
Dir[1] * tp1.
Dir[0];
2557 if (arg3 == 0)
return;
2558 double s = (arg1 - arg2) / arg3;
2560 x = (float)(tp2.
Pos[0] + s * tp2.
Dir[0]);
2561 y = (float)(tp2.
Pos[1] + s * tp2.
Dir[1]);
2569 if (tjIDs.empty())
return 0;
2571 for (
auto tjid : tjIDs) {
2572 if (tjid < 1 || tjid > (
int)slc.
tjs.size())
continue;
2573 auto& tj = slc.
tjs[tjid - 1];
2574 float sep2 =
PosSep2(tj.Pts[tj.EndPt[0]].Pos, tj.Pts[tj.EndPt[1]].Pos);
2575 if (sep2 > maxLen) maxLen = sep2;
2577 return sqrt(maxLen);
2583 float len = 0, dx, dy;
2585 unsigned short prevPt = tj.
EndPt[0];
2586 for (ipt = tj.
EndPt[0] + 1; ipt < tj.
EndPt[1] + 1; ++ipt) {
2587 if (tj.
Pts[ipt].Chg == 0)
continue;
2588 dx = tj.
Pts[ipt].Pos[0] - tj.
Pts[prevPt].Pos[0];
2589 dy = tj.
Pts[ipt].Pos[1] - tj.
Pts[prevPt].Pos[1];
2590 len += sqrt(dx * dx + dy * dy);
2599 return sqrt(
PosSep2(pos1, pos2));
2606 float d0 = pos1[0] - pos2[0];
2607 float d1 = pos1[1] - pos2[1];
2608 return d0 * d0 + d1 * d1;
2615 float dx = tp1.
Pos[0] - tp2.
Pos[0];
2616 float dy = tp1.
Pos[1] - tp2.
Pos[1];
2617 return sqrt(dx * dx + dy * dy);
2624 unsigned short& closePt,
2631 float close2 = DOCA * DOCA;
2633 bool foundClose =
false;
2635 for (
unsigned short ipt = tj.
EndPt[0]; ipt < tj.
EndPt[1] + 1; ++ipt) {
2636 if (tj.
Pts[ipt].Chg == 0)
continue;
2637 float dx = tj.
Pts[ipt].Pos[0] -
x;
2639 float dy = tj.
Pts[ipt].Pos[1] -
y;
2641 float sep2 = dx * dx + dy * dy;
2642 if (sep2 < close2) {
2649 DOCA = sqrt(close2);
2658 float dw = tp2.
Pos[0] - tp1.
Pos[0];
2659 float dt = tp2.
Pos[1] - tp1.
Pos[1];
2660 return atan2(dw, dt);
2669 std::vector<unsigned int> hitVec;
2670 if (pfp.
TP3Ds.empty())
return hitVec;
2672 for (
auto& tp3d : pfp.
TP3Ds) {
2673 if (tp3d.Flags[
kTP3DBad])
continue;
2674 if (tp3d.TjID <= 0)
continue;
2675 auto& tp = slc.
tjs[tp3d.TjID - 1].Pts[tp3d.TPIndex];
2676 for (
unsigned short ii = 0; ii < tp.Hits.size(); ++ii) {
2677 unsigned int iht = tp.Hits[ii];
2678 bool useit = (hitRequest ==
kAllHits);
2679 if (tp.UseHit[ii] && hitRequest ==
kUsedHits) useit =
true;
2680 if (!tp.UseHit[ii] && hitRequest ==
kUnusedHits) useit =
true;
2681 if (useit) hitVec.push_back(iht);
2691 std::vector<unsigned int> hitVec;
2695 for (
auto& tp : tj.
Pts)
2696 hitVec.insert(hitVec.end(), tp.Hits.begin(), tp.Hits.end());
2701 hitVec.reserve(tj.
Pts.size());
2702 for (
unsigned short ipt = 0; ipt < tj.
Pts.size(); ++ipt) {
2703 for (
unsigned short ii = 0; ii < tj.
Pts[ipt].Hits.size(); ++ii) {
2704 unsigned int iht = tj.
Pts[ipt].Hits[ii];
2705 bool useit = (hitRequest ==
kAllHits);
2706 if (tj.
Pts[ipt].UseHit[ii] && hitRequest ==
kUsedHits) useit =
true;
2707 if (!tj.
Pts[ipt].UseHit[ii] && hitRequest ==
kUnusedHits) useit =
true;
2708 if (useit) hitVec.push_back(iht);
2722 if (tj.
Pts.size() > 10)
return;
2723 if (tj.
PDGCode == 111)
return;
2725 unsigned short nhm = 0;
2726 unsigned short npwc = 0;
2727 for (
auto& tp : tj.
Pts) {
2728 if (tp.Chg == 0)
continue;
2730 unsigned short nused = 0;
2731 for (
unsigned short ii = 0; ii < tp.Hits.size(); ++ii) {
2732 if (tp.UseHit[ii]) ++nused;
2734 if (nused > 3) ++nhm;
2739 mf::LogVerbatim(
"TC") <<
"TGT: T" << tj.
ID <<
" npwc " << npwc <<
" nhm " << nhm <<
" junk? " 2748 for (
unsigned short ii = 0; ii < tjHits.size() - 1; ++ii) {
2749 for (
unsigned short jj = ii + 1; jj < tjHits.size(); ++jj) {
2750 if (tjHits[ii] == tjHits[jj]) {
2765 if (tp.
Dir[0] == 0)
return;
2766 float dw = wire - tp.
Pos[0];
2769 tp.
Pos[1] += dw * tp.
Dir[1] / tp.
Dir[0];
2774 std::array<int, 2>
const& wireWindow,
2776 const unsigned short plane,
2795 std::vector<unsigned int> closeHits;
2796 if (plane > slc.
firstWire.size() - 1)
return closeHits;
2798 int loWire = wireWindow[0];
2800 int hiWire = wireWindow[1];
2805 for (
int wire = loWire; wire <= hiWire; ++wire) {
2808 if (slc.
wireHitRange[plane][wire].first == UINT_MAX)
continue;
2809 unsigned int firstHit = slc.
wireHitRange[plane][wire].first;
2810 unsigned int lastHit = slc.
wireHitRange[plane][wire].second;
2811 for (
unsigned int iht = firstHit; iht <= lastHit; ++iht) {
2814 if (
hit.PeakTime() < minTick)
continue;
2815 if (
hit.PeakTime() > maxTick)
break;
2819 if (
hit.StartTick() > hiLo) hiLo =
hit.StartTick();
2821 if (
hit.EndTick() < loHi) loHi =
hit.EndTick();
2822 if (loHi < hiLo)
continue;
2823 if (hiLo > loHi)
break;
2826 bool takeit = (hitRequest ==
kAllHits);
2827 if (hitRequest ==
kUsedHits && slc.
slHits[iht].InTraj > 0) takeit =
true;
2829 if (takeit) closeHits.push_back(iht);
2838 float const maxDelta,
2849 if (tp.
Pos[0] < -0.4)
return false;
2851 unsigned int wire = std::nearbyint(tp.
Pos[0]);
2852 if (wire < slc.
firstWire[plane])
return false;
2853 if (wire > slc.
lastWire[plane] - 1)
return false;
2862 if (slc.
wireHitRange[plane][wire].first == UINT_MAX)
return false;
2864 unsigned int firstHit = slc.
wireHitRange[plane][wire].first;
2865 unsigned int lastHit = slc.
wireHitRange[plane][wire].second;
2868 for (
unsigned int iht = firstHit; iht <= lastHit; ++iht) {
2869 if ((
unsigned int)slc.
slHits[iht].InTraj > slc.
tjs.size())
continue;
2870 bool useit = (hitRequest ==
kAllHits);
2871 if (hitRequest ==
kUsedHits && slc.
slHits[iht].InTraj > 0) useit =
true;
2873 if (!useit)
continue;
2877 if (delta < maxDelta) tp.
Hits.push_back(iht);
2879 if (tp.
Hits.size() > 16) { tp.
Hits.resize(16); }
2882 return (!tp.
Hits.empty());
2894 if (end > 1)
return USHRT_MAX;
2896 if (end == 1) dir = -1;
2897 for (
short ii = 0; ii < (short)tj.
Pts.size(); ++ii) {
2898 short ipt = tj.
EndPt[
end] + dir * ii;
2899 if (ipt < 0 || ipt >= (
short)tj.
Pts.size())
return USHRT_MAX;
2900 auto& tp = tj.
Pts[ipt];
2910 const float& maxDelta)
2928 std::vector<int>
tmp;
2929 if (fromTp.
Pos[0] < -0.4 || toTp.
Pos[0] < -0.4)
return tmp;
2933 unsigned int firstWire, lastWire;
2934 if (toTp.
Pos[0] > fromTp.
Pos[0]) {
2936 firstWire = std::nearbyint(fromTp.
Pos[0]);
2937 lastWire = std::nearbyint(toTp.
Pos[0]);
2939 else if (toTp.
Pos[0] < fromTp.
Pos[0]) {
2941 firstWire = std::nearbyint(toTp.
Pos[0]);
2942 lastWire = std::nearbyint(fromTp.
Pos[0]);
2946 float tmp = fromTp.
Pos[0] - maxDelta;
2947 if (tmp < 0) tmp = 0;
2948 firstWire = std::nearbyint(tmp);
2949 tmp = fromTp.
Pos[0] + maxDelta;
2950 lastWire = std::nearbyint(tmp);
2956 if (firstWire > slc.
lastWire[plane] - 1)
return tmp;
2958 if (lastWire > slc.
lastWire[plane] - 1) lastWire = slc.
lastWire[plane] - 1;
2960 for (
unsigned int wire = firstWire; wire <= lastWire; ++wire) {
2961 if (slc.
wireHitRange[plane][wire].first == UINT_MAX)
continue;
2966 unsigned int firstHit = slc.
wireHitRange[plane][wire].first;
2967 unsigned int lastHit = slc.
wireHitRange[plane][wire].second;
2968 for (
unsigned int iht = firstHit; iht <= lastHit; ++iht) {
2969 if (slc.
slHits[iht].InTraj <= 0)
continue;
2970 if ((
unsigned int)slc.
slHits[iht].InTraj > slc.
tjs.size())
continue;
2972 if (
hit.PeakTime() < minTick)
continue;
2974 if (
hit.PeakTime() > maxTick)
break;
2975 if (std::find(tmp.begin(), tmp.end(), slc.
slHits[iht].InTraj) != tmp.end())
continue;
2976 tmp.push_back(slc.
slHits[iht].InTraj);
2987 unsigned short end1,
2989 unsigned short end2,
2990 unsigned short nPtsFit,
2997 if (tj1.
CTP != tj2.
CTP)
return -1;
2998 if (end1 > 1 || end2 > 1)
return -1;
3006 if (end1 == 1) dir = -1;
3007 unsigned short cnt = 0;
3009 for (
short ii = 0; ii < (short)tj1.
Pts.size(); ++ii) {
3010 short ipt = tj1.
EndPt[end1] + dir * ii;
3012 if (ipt >= (
short)tj1.
Pts.size())
break;
3013 auto& tp = tj1.
Pts[ipt];
3014 if (tp.Chg <= 0)
continue;
3015 tj.
Pts.push_back(tp);
3017 if (cnt == nPtsFit + 1)
break;
3019 if (cnt < nPtsFit)
return -1;
3022 if (end2 == 1) dir = -1;
3024 for (
short ii = 0; ii < (short)tj2.
Pts.size(); ++ii) {
3025 short ipt = tj2.
EndPt[end2] + dir * ii;
3027 if (ipt >= (
short)tj2.
Pts.size())
break;
3028 auto& tp = tj2.
Pts[ipt];
3029 if (tp.Chg <= 0)
continue;
3030 tj.
Pts.push_back(tp);
3032 if (cnt == nPtsFit + 1)
break;
3042 unsigned short kinkPt,
3043 unsigned short nPtsFit,
3052 if (kinkPt < tj.
EndPt[0] + 2)
return -1;
3053 if (kinkPt > tj.
EndPt[1] - 2)
return -1;
3056 if (nPtsFit < 3)
return -1;
3059 if (npwc < 2 * nPtsFit + 1)
return -1;
3064 double chgRMS = 0.07;
3068 double tFactor = 1 + 0.3 / double(nPtsFit - 2);
3074 FitTraj(slc, tj, kinkPt, nPtsFit, fitDir, tpPos);
3075 if (tpPos.
FitChi > 900)
return -1;
3079 FitTraj(slc, tj, kinkPt, nPtsFit, fitDir, tpNeg);
3080 if (tpNeg.
FitChi > 900)
return -1;
3081 double angErr = tpNeg.
AngErr;
3085 double dangSig = dang / angErr;
3092 unsigned short cntNeg = 0;
3093 for (
unsigned short ipt = kinkPt - 1; ipt >= tj.
EndPt[0]; --ipt) {
3094 auto& tp = tj.
Pts[ipt];
3095 if (tp.Chg <= 0)
continue;
3098 if (cntNeg == nPtsFit)
break;
3099 if (ipt == 0)
break;
3101 if (cntNeg != nPtsFit) {
3102 if (prt)
mf::LogVerbatim(
"TC") <<
" KL: Bad cntNeg " << cntNeg <<
" != " << nPtsFit;
3107 unsigned short cntPos = 0;
3108 for (
unsigned short ipt = kinkPt + 1; ipt <= tj.
EndPt[1]; ++ipt) {
3109 auto& tp = tj.
Pts[ipt];
3110 if (tp.Chg <= 0)
continue;
3113 if (cntPos == nPtsFit)
break;
3115 if (cntPos != nPtsFit) {
3116 if (prt)
mf::LogVerbatim(
"TC") <<
" KL: Bad cntPos " << cntPos <<
" != " << nPtsFit;
3119 chgNeg /= (float)nPtsFit;
3120 chgPos /= (float)nPtsFit;
3122 chgAsym =
std::abs(chgPos - chgNeg) / (chgPos + chgNeg);
3124 chgSig = chgAsym / chgRMS;
3126 double kinkSig = sqrt(dangSig * dangSig + chgSig * chgSig);
3130 myprt <<
"KL: T" << tj.
ID <<
" kinkPt " <<
PrintPos(tj.
Pts[kinkPt]);
3131 myprt <<
" nPtsFit " << nPtsFit;
3132 myprt <<
" dang " << std::fixed << std::setprecision(3) << dang;
3133 myprt << std::fixed << std::setprecision(3) <<
" angErr " << angErr;
3134 myprt << std::setprecision(2) <<
" sig " << dangSig;
3135 myprt <<
" chgAsym " << chgAsym;
3136 myprt <<
" chgSig " << chgSig;
3137 myprt <<
" kinkSig " << kinkSig;
3139 return (
float)kinkSig;
3149 unsigned short midPt = 0.5 * (tj.
EndPt[0] + tj.
EndPt[1]);
3150 double rms0 = 0, rms1 = 0;
3154 float asym =
std::abs(rms0 - rms1) / (rms0 + rms1);
3155 float chgFact = (tj.
ChgRMS - 0.1) * 5;
3156 float elh = 5 * asym * chgFact;
3157 if (elh > 1) elh = 1;
3166 if (tjIDs.empty())
return 0;
3167 std::array<int, 2> wireWindow;
3170 constexpr
float NNDelta = 5;
3171 wireWindow[0] = pos[0] - NNDelta;
3172 wireWindow[1] = pos[0] + NNDelta;
3173 timeWindow[0] = pos[1] - NNDelta;
3174 timeWindow[1] = pos[1] + NNDelta;
3176 for (
auto& tjID : tjIDs)
3177 if (tjID <= 0 || tjID > (
int)slc.
tjs.size())
return 0;
3182 std::vector<unsigned int> closeHits =
3184 if (closeHits.empty())
return 0;
3189 for (
auto& iht : closeHits) {
3191 chg +=
hit.Integral();
3192 if (slc.
slHits[iht].InTraj == 0)
continue;
3193 if (std::find(tjIDs.begin(), tjIDs.end(), slc.
slHits[iht].InTraj) != tjIDs.end())
3194 tchg +=
hit.Integral();
3196 if (chg == 0)
return 0;
3203 float delta, md = 0;
3206 for (
auto& tp : tj.
Pts) {
3207 for (ii = 0; ii < tp.Hits.size(); ++ii) {
3208 if (!tp.UseHit[ii])
continue;
3211 if (delta > md) md = delta;
3221 if (tj.
Pts.empty())
return;
3227 std::reverse(tj.
Pts.begin(), tj.
Pts.end());
3232 for (
unsigned short ipt = 0; ipt < tj.
Pts.size(); ++ipt) {
3233 if (tj.
Pts[ipt].Dir[0] != 0) tj.
Pts[ipt].Dir[0] = -tj.
Pts[ipt].Dir[0];
3234 if (tj.
Pts[ipt].Dir[1] != 0) tj.
Pts[ipt].Dir[1] = -tj.
Pts[ipt].Dir[1];
3235 if (tj.
Pts[ipt].Ang > 0) { tj.
Pts[ipt].Ang -= M_PI; }
3237 tj.
Pts[ipt].Ang += M_PI;
3253 unsigned short nvx = Envelope.size();
3254 double angleSum = 0;
3255 for (
unsigned short ii = 0; ii < Envelope.size(); ++ii) {
3256 p1[0] = Envelope[ii][0] - Point[0];
3257 p1[1] = Envelope[ii][1] - Point[1];
3258 p2[0] = Envelope[(ii + 1) % nvx][0] - Point[0];
3259 p2[1] = Envelope[(ii + 1) % nvx][1] - Point[1];
3262 if (
abs(angleSum) < M_PI)
return false;
3270 double den = v1[0] * v1[0] + v1[1] * v1[1];
3271 if (den == 0)
return false;
3285 if (pos1[0] == pos2[0] && pos1[1] == pos2[1])
return;
3286 pos1[0] = pos2[0] - pos1[0];
3287 pos1[1] = pos2[1] - pos1[1];
3288 double sep = sqrt(pos1[0] * pos1[0] + pos1[1] * pos1[1]);
3289 if (sep < 1
E-6)
return;
3291 ptDir[0] = pos1[0] / sep;
3292 ptDir[1] = pos1[1] / sep;
3294 double costh =
DotProd(dir1, ptDir);
3295 if (costh > 1.0 || costh < -1.0)
return;
3296 alongTrans[0] = costh * sep;
3297 double sinth = sqrt(1 - costh * costh);
3298 alongTrans[1] = sinth * sep;
3305 double ang1 = atan2(p1[1], p1[0]);
3306 double ang2 = atan2(p2[1], p2[0]);
3313 constexpr
double twopi = 2 * M_PI;
3314 double dang = Ang1 - Ang2;
3317 while (dang < -M_PI)
3325 return std::abs(std::remainder(Ang1 - Ang2, M_PI));
3338 if (tj.
Pts.size() == 0)
return;
3341 for (
unsigned short ipt = 0; ipt < tj.
Pts.size(); ++ipt) {
3342 if (tj.
Pts[ipt].Chg != 0) {
3347 for (
unsigned short ii = 0; ii < tj.
Pts.size(); ++ii) {
3348 unsigned short ipt = tj.
Pts.size() - 1 - ii;
3349 if (tj.
Pts[ipt].Chg != 0) {
3361 unsigned short nUsed = 0;
3362 unsigned short nTotHits = 0;
3363 for (
unsigned short ipt = tj.
EndPt[0]; ipt <= tj.
EndPt[1]; ++ipt) {
3365 nTotHits += tp.
Hits.size();
3366 for (
unsigned short ii = 0; ii < tp.
Hits.size(); ++ii) {
3367 if (tp.
UseHit[ii]) ++nUsed;
3370 if (nTotHits == 0)
return false;
3371 float fracUsed = (float)nUsed / (
float)nTotHits;
3373 mf::LogVerbatim(
"TC") <<
"TrajIsClean: nTotHits " << nTotHits <<
" nUsed " << nUsed
3374 <<
" fracUsed " << fracUsed;
3376 return fracUsed > 0.9;
3383 if (tjIDs.empty())
return 0;
3386 for (
auto tjid : tjIDs) {
3387 auto& tj = slc.
tjs[tjid - 1];
3388 float npts = tj.EndPt[1] - tj.EndPt[0] + 1;
3389 summ += npts * tj.MCSMom;
3392 return (
short)(summ / suml);
3404 unsigned short firstPt,
3405 unsigned short lastPt)
3409 if (firstPt == lastPt)
return 0;
3410 if (firstPt > lastPt) std::swap(firstPt, lastPt);
3414 if (firstPt >= lastPt)
return 0;
3416 if (firstPt < tj.
EndPt[0])
return 0;
3417 if (lastPt > tj.
EndPt[1])
return 0;
3424 if (tjLen < 1)
return 0;
3426 double thetaRMS =
MCSThetaRMS(tj, firstPt, lastPt);
3427 if (thetaRMS < 0.001)
return 999;
3428 double mom = 13.8 * sqrt(tjLen / 14) / thetaRMS;
3429 if (mom > 999) mom = 999;
3437 if (thePt > tj.
EndPt[1])
return thePt;
3438 if (tj.
Pts[thePt].Chg > 0)
return thePt;
3440 short endPt0 = tj.
EndPt[0];
3441 short endPt1 = tj.
EndPt[1];
3442 for (
short off = 1; off < 10; ++off) {
3443 short ipt = thePt + off;
3444 if (ipt <= endPt1 && tj.
Pts[ipt].Chg > 0)
return (
unsigned short)ipt;
3446 if (ipt >= endPt0 && tj.
Pts[ipt].Chg > 0)
return (
unsigned short)ipt;
3459 if (tps < 1)
return 1;
3471 if (firstPt < tj.
EndPt[0])
return 1;
3472 if (lastPt > tj.
EndPt[1])
return 1;
3476 if (firstPt >= lastPt)
return 1;
3480 TjDeltaRMS(tj, firstPt, lastPt, sigmaS, cnt);
3481 if (sigmaS < 0)
return 1;
3483 if (tjLen < 1)
return 1;
3485 return (6.8 * sigmaS / tjLen);
3491 unsigned short firstPt,
3492 unsigned short lastPt,
3494 unsigned short& cnt)
3499 if (firstPt < tj.
EndPt[0])
return;
3500 if (lastPt > tj.
EndPt[1])
return;
3504 if (firstPt >= lastPt)
return;
3517 for (
unsigned short ipt = firstPt + 1; ipt < lastPt; ++ipt) {
3518 if (tj.
Pts[ipt].Chg == 0)
continue;
3520 if (tj.
Pts[ipt].HitPosErr2 > 4)
continue;
3524 if (cnt < 2)
return;
3525 rms = sqrt(dsum / (
double)cnt);
3536 std::array<int, 2> wireWindow;
3542 for (
auto& mutj : slc.
tjs) {
3543 if (mutj.AlgMod[
kKilled])
continue;
3544 if (mutj.CTP != inCTP)
continue;
3545 if (mutj.PDGCode != 13)
continue;
3546 unsigned short nnear = 0;
3547 for (
unsigned short ipt = mutj.EndPt[0]; ipt <= mutj.EndPt[1]; ++ipt) {
3548 auto& tp = mutj.Pts[ipt];
3549 wireWindow[0] = tp.Pos[0];
3550 wireWindow[1] = tp.Pos[0];
3551 timeWindow[0] = tp.Pos[1] - delta;
3552 timeWindow[1] = tp.Pos[1] + delta;
3557 if (closeHits.empty())
continue;
3558 for (
auto iht : closeHits) {
3559 auto inTraj = slc.
slHits[iht].InTraj;
3560 if (inTraj <= 0)
continue;
3561 if (inTraj == mutj.ID)
continue;
3562 auto& dtj = slc.
tjs[inTraj - 1];
3563 if (dtj.PDGCode == 13)
continue;
3564 for (
unsigned short jpt = dtj.EndPt[0]; jpt <= dtj.EndPt[1]; ++jpt) {
3565 auto& dtp = dtj.Pts[jpt];
3566 if (std::find(dtp.Hits.begin(), dtp.Hits.end(), iht) == dtp.Hits.end())
continue;
3585 for (
auto& tp : tj.
Pts) {
3586 if (tp.Chg <= 0)
continue;
3588 for (
unsigned short ii = 0; ii < tp.Hits.size(); ++ii) {
3589 if (tp.UseHit[ii])
continue;
3590 unsigned int iht = tp.Hits[ii];
3616 for (
unsigned short ipt = tj.
EndPt[0] + 1; ipt < tj.
EndPt[1]; ++ipt) {
3617 auto& tp = tj.
Pts[ipt];
3618 if (tp.Chg > bigChg) bigChg = tp.Chg;
3625 for (
unsigned short ipt = tj.
EndPt[0] + 1; ipt < tj.
EndPt[1]; ++ipt) {
3626 auto& tp = tj.
Pts[ipt];
3627 if (tp.Chg <= 0)
continue;
3629 if (tp.Chg == bigChg)
continue;
3634 bsum2 += tp.Chg * tp.Chg;
3635 if (tp.Chg > bigChg) bigChg = tp.Chg;
3641 for (
unsigned short ii = 0; ii < tj.
Pts[ipt].Hits.size(); ++ii) {
3642 if (!tp.UseHit[ii])
continue;
3643 unsigned int iht = tp.Hits[ii];
3647 vsum2 += tpchg * tpchg;
3650 if (bcnt == 0)
return;
3658 if (arg > 0) tj.
ChgRMS = sqrt(arg / (bcnt - 1));
3660 for (
auto& tp : tj.
Pts)
3663 mf::LogVerbatim(
"TC") << inFcnLabel <<
".UpdateTjChgProperties: backup sum Set tj.AveChg " 3668 double nWires = tj.
EndPt[1] - tj.
EndPt[0] + 1;
3669 if (nWires < 2)
return;
3674 for (
unsigned short end = 0;
end < 2; ++
end) {
3678 int dw =
std::abs(tp.Pos[0] - vx2.Pos[0]);
3690 if (arg > 0) rms = sqrt(arg / (vcnt - 1));
3693 if (rms < 0.1) rms = 0.1;
3697 double defFrac = 1 / vcnt;
3698 rms = defFrac * 0.5 + (1 - defFrac) * rms;
3702 mf::LogVerbatim(
"TC") << inFcnLabel <<
".UpdateTjChgProperties: Set tj.AveChg " 3709 for (
unsigned short ipt = tj.
EndPt[0]; ipt <= tj.
EndPt[1]; ++ipt) {
3710 auto& tp = tj.
Pts[ipt];
3711 if (tp.Chg <= 0)
continue;
3718 for (
auto& tp : tj.
Pts)
3724 for (
auto& tp : tj.
Pts)
3728 unsigned short minPt = tj.
EndPt[0] + nptsave;
3730 for (
unsigned short ii = 0; ii < tj.
Pts.size(); ++ii) {
3731 unsigned short ipt = tj.
EndPt[1] - ii;
3732 if (ipt < minPt)
break;
3735 for (
unsigned short iii = 0; iii < nptsave; ++iii) {
3736 unsigned short iipt = ipt - iii;
3738 if (iipt == tj.
EndPt[0])
break;
3739 auto& tp = tj.
Pts[iipt];
3740 if (tp.Chg <= 0)
continue;
3745 tj.
Pts[ipt].AveChg = sum / cnt;
3746 lastAve = tj.
Pts[ipt].AveChg;
3750 for (
unsigned short ii = tj.
EndPt[0]; ii <= tj.
EndPt[1]; ++ii) {
3751 unsigned short ipt = tj.
EndPt[1] - ii;
3752 auto& tp = tj.
Pts[ipt];
3753 if (tp.AveChg == 0) { tp.AveChg = lastAve; }
3755 lastAve = tp.AveChg;
3772 for (
auto& tj : slc.
tjs) {
3773 if (tj.AlgMod[
kKilled])
continue;
3774 for (
auto& tp : tj.Pts)
3778 for (
auto& vx : slc.
vtxs) {
3779 if (vx.ID <= 0)
continue;
3790 if (vx2.
ID == 0)
return;
3793 if (prt)
mf::LogVerbatim(
"TC") <<
"UpdateVxEnvironment check Tjs attached to vx2 " << vx2.
ID;
3795 std::vector<int> tjlist;
3796 std::vector<unsigned short> tjends;
3797 if (vx2.
Pos[0] < -0.4)
return;
3798 unsigned int vxWire = std::nearbyint(vx2.
Pos[0]);
3799 unsigned int loWire = vxWire;
3800 unsigned int hiWire = vxWire;
3801 for (
auto& tj : slc.
tjs) {
3803 if (tj.CTP != vx2.
CTP)
continue;
3805 if (tj.AlgMod[
kPhoton])
continue;
3806 for (
unsigned short end = 0;
end < 2; ++
end) {
3807 if (tj.VtxID[
end] != vx2.
ID)
continue;
3808 tjlist.push_back(tj.ID);
3809 tjends.push_back(
end);
3810 if (tj.Pts[tj.EndPt[
end]].Pos[0] < -0.4)
return;
3811 unsigned int endWire = std::nearbyint(tj.Pts[tj.EndPt[
end]].Pos[0]);
3812 if (endWire < loWire) loWire = endWire;
3813 if (endWire > hiWire) hiWire = endWire;
3816 if (tjlist.size() < 2)
return;
3817 if (hiWire < loWire + 1)
return;
3819 mf::LogVerbatim(
"TC") <<
" check Tjs on wires in the range " << loWire <<
" to " << hiWire;
3823 std::vector<std::vector<TrajPoint>> wire_tjpt;
3825 std::vector<int> tjids;
3827 unsigned short nwires = hiWire - loWire + 1;
3828 for (
unsigned short itj = 0; itj < tjlist.size(); ++itj) {
3829 auto& tj = slc.
tjs[tjlist[itj] - 1];
3830 unsigned short end = tjends[itj];
3831 std::vector<TrajPoint> tjpt(nwires);
3833 for (
unsigned short ii = 0; ii < tj.Pts.size(); ++ii) {
3835 if (end == 0) { ipt = tj.EndPt[0] + ii; }
3837 ipt = tj.EndPt[1] - ii;
3839 if (ipt > tj.Pts.size() - 1)
break;
3841 auto tp = tj.Pts[ipt];
3842 if (tp.Chg <= 0)
continue;
3845 if (tp.Pos[0] < -0.4)
continue;
3846 unsigned int wire = std::nearbyint(tp.Pos[0]);
3847 unsigned short indx = wire - loWire;
3848 if (indx > nwires - 1)
break;
3858 if (ltp.
Dir[0] == 0)
continue;
3859 if (ltp.
Pos[0] < -0.4)
continue;
3860 unsigned int wire = std::nearbyint(ltp.
Pos[0]);
3862 unsigned short indx = wire - loWire;
3864 if (tjpt[indx].Chg == 0) tjpt[indx] = ltp;
3866 for (
unsigned short ii = 0; ii < nwires; ++ii) {
3868 for (
unsigned short iwt = 0; iwt < 2; ++iwt)
3869 ltp.
Pos[iwt] += ltp.
Dir[iwt] * stepSize;
3870 if (ltp.
Pos[0] < -0.4)
break;
3871 wire = std::nearbyint(ltp.
Pos[0]);
3872 if (wire < loWire || wire > hiWire)
break;
3873 indx = wire - loWire;
3874 if (tjpt[indx].Chg > 0)
continue;
3879 myprt <<
" T" << tj.ID;
3880 for (
auto& tp : tjpt)
3881 myprt <<
" " <<
PrintPos(tp.Pos) <<
"_" << tp.Step <<
"_" << (int)tp.Chg;
3883 wire_tjpt.push_back(tjpt);
3884 tjids.push_back(tj.ID);
3888 for (
unsigned short indx = 0; indx < nwires; ++indx) {
3890 unsigned short npts = 0;
3892 unsigned short npwc = 0;
3893 for (
unsigned short itj = 0; itj < wire_tjpt.size(); ++itj) {
3894 if (wire_tjpt[itj][indx].Pos[0] == 0)
continue;
3897 if (wire_tjpt[itj][indx].Chg > 0) ++npwc;
3900 if (npts == 0)
continue;
3902 if (npwc == npts)
continue;
3904 for (
unsigned short itj = 0; itj < wire_tjpt.size(); ++itj) {
3905 if (wire_tjpt[itj][indx].Pos[0] == 0)
continue;
3906 if (wire_tjpt[itj][indx].Chg == 0)
continue;
3907 auto& tj = slc.
tjs[tjids[itj] - 1];
3908 unsigned short ipt = wire_tjpt[itj][indx].Step;
3910 tj.NeedsUpdate =
true;
3911 if (prt)
mf::LogVerbatim(
"TC") <<
" Set kEnvOverlap bit on T" << tj.ID <<
" ipt " << ipt;
3917 for (
auto tjid : tjids) {
3918 auto& tj = slc.
tjs[tjid - 1];
3919 if (!tj.NeedsUpdate)
continue;
3920 if (tj.CTP != vx2.
CTP)
continue;
3964 if (dir[0] == 0 && dir[1] == 0 && dir[2] == 0)
return tp;
3968 Point3_t pos3 = {{100 * dir[0], 100 * dir[1], 100 * dir[2]}};
3970 std::array<double, 2> ori2;
3971 std::array<double, 2> pos2;
3972 std::array<double, 2> dir2;
3980 dir2[0] = pos2[0] - ori2[0];
3981 dir2[1] = pos2[1] - ori2[1];
3983 double norm = sqrt(dir2[0] * dir2[0] + dir2[1] * dir2[1]);
3986 tp.
Ang = atan2(dir2[1], dir2[0]);
3987 tp.
Delta = norm / 100;
3995 norm = sqrt(cs * cs + sn * sn);
4000 tp.
DeltaRMS = 100 / (pos2[0] - ori2[0]);
4007 unsigned int fromHit,
4011 if (fromHit > slc.
slHits.size() - 1)
return false;
4012 if (toHit > slc.
slHits.size() - 1)
return false;
4018 (float)thit.WireID().Wire,
4034 tp.
Pos[0] = fromWire;
4036 tp.
Dir[0] = toWire - fromWire;
4039 if (norm == 0)
return false;
4049 tpOut.
Pos = fromPos;
4051 tpOut.
Ang = atan2(tpOut.
Dir[1], tpOut.
Dir[0]);
4062 tpOut.
Ang = atan2(tpOut.
Dir[1], tpOut.
Dir[0]);
4070 if (tj.
ID == 0)
return 0;
4080 for (
unsigned short xyz = 0; xyz < 2; ++xyz)
4081 dir[xyz] = p2[xyz] - p1[xyz];
4082 if (dir[0] == 0 && dir[1] == 0)
return dir;
4083 double norm = sqrt(dir[0] * dir[0] + dir[1] * dir[1]);
4101 if (tp.
Hits.empty())
return 0;
4102 float minVal = 9999;
4104 for (
unsigned short ii = 0; ii < tp.
Hits.size(); ++ii) {
4105 bool useit = (hitRequest ==
kAllHits);
4108 if (!useit)
continue;
4109 unsigned int iht = tp.
Hits[ii];
4111 float cv =
hit.PeakTime();
4112 float rms =
hit.RMS();
4113 float arg = cv - rms;
4114 if (arg < minVal) minVal = arg;
4116 if (arg > maxVal) maxVal = arg;
4118 if (maxVal == 0)
return 0;
4119 return (maxVal - minVal) / 2;
4124 const std::vector<unsigned int>& hitsInMultiplet,
4132 const std::vector<unsigned int>& hitsInMultiplet,
4135 if (hitsInMultiplet.empty())
return 0;
4137 if (hitsInMultiplet.size() == 1) {
4142 float minVal = 9999;
4144 for (
unsigned short ii = 0; ii < hitsInMultiplet.size(); ++ii) {
4145 unsigned int iht = hitsInMultiplet[ii];
4146 bool useit = (hitRequest ==
kAllHits);
4147 if (hitRequest ==
kUsedHits && slc.
slHits[iht].InTraj > 0) useit =
true;
4149 if (!useit)
continue;
4151 float cv =
hit.PeakTime();
4152 float rms =
hit.RMS();
4153 float arg = cv - rms;
4154 if (arg < minVal) minVal = arg;
4156 if (arg > maxVal) maxVal = arg;
4158 if (maxVal == 0)
return 0;
4159 return (maxVal - minVal) / 2;
4164 const std::vector<unsigned int>& hitsInMultiplet,
4173 const std::vector<unsigned int>& hitsInMultiplet,
4180 for (
unsigned short ii = 0; ii < hitsInMultiplet.size(); ++ii) {
4181 unsigned int iht = hitsInMultiplet[ii];
4182 bool useit = (hitRequest ==
kAllHits);
4183 if (hitRequest ==
kUsedHits && slc.
slHits[iht].InTraj > 0) useit =
true;
4185 if (!useit)
continue;
4187 float chg =
hit.Integral();
4188 pos += chg *
hit.PeakTime();
4191 if (sum <= 0)
return -1;
4199 if (tp.
Hits.empty())
return 0;
4203 unsigned short nhits = 0;
4204 for (
unsigned short ii = 0; ii < tp.
Hits.size(); ++ii) {
4206 if (tp.
UseHit[ii]) ++nhits;
4210 if (!tp.
UseHit[ii]) ++nhits;
4219 if (itj > slc.
tjs.size() - 1)
return;
4253 if (npwc > 500) isAMuon =
true;
4272 std::vector<float> cnt(nplanes, 0);
4273 for (
unsigned short iht = 0; iht < (*
evt.
allHits).
size(); ++iht) {
4275 unsigned short plane =
hit.WireID().Plane;
4276 if (plane > nplanes - 1)
return false;
4277 if (cnt[plane] > 200)
continue;
4279 if (
hit.Multiplicity() != 1)
continue;
4281 if (
hit.GoodnessOfFit() < 0 ||
hit.GoodnessOfFit() > 500)
continue;
4283 if (
hit.PeakAmplitude() < 1)
continue;
4287 bool allDone =
true;
4288 for (
unsigned short plane = 0; plane < nplanes; ++plane)
4289 if (cnt[plane] < 200) allDone =
false;
4295 for (
unsigned short plane = 0; plane < nplanes; ++plane) {
4296 if (cnt[plane] > 4) {
evt.
aveHitRMS[plane] /= cnt[plane]; }
4304 std::cout <<
"Analyze hits aveHitRMS";
4305 std::cout << std::fixed << std::setprecision(1);
4307 std::cout <<
" " << rms;
4336 lariov::ChannelStatusProvider
const& channelStatus =
4343 for (
unsigned int wire = 0; wire < nwires; ++wire) {
4345 evt.
goodWire[
id.Plane][wire] = channelStatus.IsGood(chan);
4367 for (
unsigned int wire = 0; wire < nwires; ++wire)
4372 unsigned int nBadWireFix = 0;
4375 auto const& wid =
hit.WireID();
4377 if (static_cast<geo::TPCID const&>(wid) != inTPCID)
continue;
4378 unsigned short pln = wid.Plane;
4379 unsigned int wire = wid.Wire;
4389 std::cout <<
"FillWireHitRange found hits on " << nBadWireFix
4390 <<
" wires that were declared not-good by the ChannelStatus service. Fixed it...\n";
4406 unsigned short nplanes = tpcgeom.
Nplanes();
4408 if (nplanes > 3)
return false;
4411 auto const world = tpcgeom.GetCenter();
4414 slc.
xLo = world.X() - tpcgeom.HalfWidth() + 1;
4415 slc.
xHi = world.X() + tpcgeom.HalfWidth() - 1;
4416 slc.
yLo = world.Y() - tpcgeom.HalfHeight() + 1;
4417 slc.
yHi = world.Y() + tpcgeom.HalfHeight() - 1;
4418 slc.
zLo = world.Z() - tpcgeom.Length() / 2 + 1;
4419 slc.
zHi = world.Z() + tpcgeom.Length() / 2 - 1;
4425 slc.
nWires.resize(nplanes);
4430 std::pair<unsigned int, unsigned int> flag{UINT_MAX, UINT_MAX};
4444 for (
unsigned short plane = 0; plane < nplanes; ++plane) {
4453 unsigned int lastWire = 0, lastPlane = 0;
4454 for (
unsigned int iht = 0; iht < slc.
slHits.size(); ++iht) {
4455 unsigned int ahi = slc.
slHits[iht].allHitsIndex;
4456 if (ahi > (*
evt.
allHits).size() - 1)
return false;
4458 if (
hit.WireID().Cryostat != cstat)
continue;
4459 if (
hit.WireID().TPC != tpc)
continue;
4460 unsigned short plane =
hit.WireID().Plane;
4461 unsigned int wire =
hit.WireID().Wire;
4462 if (wire > slc.
nWires[plane] - 1) {
4463 mf::LogWarning(
"TC") <<
"FillWireHitRange: Invalid wire number " << wire <<
" > " 4464 << slc.
nWires[plane] - 1 <<
" in plane " << plane <<
" Quitting";
4467 if (plane == lastPlane && wire < lastWire) {
4469 <<
"FillWireHitRange: Hits are not in increasing wire order. Quitting ";
4481 unsigned int slhitsSize = slc.
slHits.size();
4482 for (
unsigned short plane = 0; plane < nplanes; ++plane) {
4483 for (
unsigned int wire = slc.
firstWire[plane]; wire < slc.
lastWire[plane]; ++wire) {
4484 if (slc.
wireHitRange[plane][wire].first == UINT_MAX)
continue;
4485 if (slc.
wireHitRange[plane][wire].first > slhitsSize - 1 &&
4495 std::cout <<
"Slice ID/Index " << slc.
ID <<
"/" <<
slices.size() <<
" tpc " << tpc
4497 std::cout <<
" Active volume (";
4498 std::cout << std::fixed << std::setprecision(1) << slc.
xLo <<
" < X < " << slc.
xHi <<
") (";
4499 std::cout << std::fixed << std::setprecision(1) << slc.
yLo <<
" < Y < " << slc.
yHi <<
") (";
4500 std::cout << std::fixed << std::setprecision(1) << slc.
zLo <<
" < Z < " << slc.
zHi <<
")\n";
4525 if (itj1 > slc.
tjs.size() - 1)
return false;
4526 if (itj2 > slc.
tjs.size() - 1)
return false;
4537 if (pfp1 != USHRT_MAX || pfp2 != USHRT_MAX) {
4538 if (pfp1 != USHRT_MAX && pfp2 != USHRT_MAX)
return false;
4540 if (pfp1 == USHRT_MAX) std::swap(itj1, itj2);
4562 std::swap(tj1, tj2);
4563 std::swap(tp1e0, tp2e0);
4564 std::swap(tp1e1, tp2e1);
4582 if (tp2e0[0] > tp1e0[0] && tp2e1[0] < tp1e1[0])
return false;
4586 if (tp1e0[0] > tp2e0[0] && tp1e1[0] < tp2e1[0])
return false;
4590 if (tp2e1[0] > tp1e1[0] && tp2e0[0] < tp1e0[0])
return false;
4591 if (tp1e1[0] > tp2e1[0] && tp1e0[0] < tp2e0[0])
return false;
4606 mf::LogVerbatim(
"TC") <<
"MergeAndStore: You are merging the end of trajectory T" << tj1.
ID 4607 <<
" with a Bragg peak. Not merging\n";
4618 float minSep = 1000;
4619 unsigned short tj2ClosePt = 0;
4624 << tj2ClosePt <<
" Pos " <<
PrintPos(tj2.
Pts[tj2ClosePt]);
4626 if (tj2ClosePt > tj2.
EndPt[1])
return false;
4635 std::vector<unsigned int> tj1Hits;
4636 for (
unsigned short ii = 0; ii < tj1.
Pts.size(); ++ii) {
4639 unsigned short ipt = tj1.
Pts.size() - 1 - ii;
4640 tj1Hits.insert(tj1Hits.end(), tj1.
Pts[ipt].Hits.begin(), tj1.
Pts[ipt].Hits.end());
4641 if (ipt == 0)
break;
4644 bool bumpedPt =
true;
4647 for (
unsigned short ii = 0; ii < tj2.
Pts[tj2ClosePt].Hits.size(); ++ii) {
4648 unsigned int iht = tj2.
Pts[tj2ClosePt].Hits[ii];
4649 if (std::find(tj1Hits.begin(), tj1Hits.end(), iht) != tj1Hits.end()) bumpedPt =
true;
4651 if (bumpedPt && tj2ClosePt < tj2.
EndPt[1]) { ++tj2ClosePt; }
4656 if (doPrt)
mf::LogVerbatim(
"TC") <<
" revised tj2ClosePt " << tj2ClosePt;
4659 tj1.
Pts.insert(tj1.
Pts.end(), tj2.
Pts.begin() + tj2ClosePt, tj2.
Pts.end());
4666 if (tj2.
VtxID[1] > 0) {
4687 int newTjID = slc.
tjs.size();
4691 if (doPrt)
mf::LogVerbatim(
"TC") <<
" MAS success. Created T" << newTjID;
4693 for (
auto& tj : slc.
tjs)
4694 if (tj.ParentID == tj1ID || tj.ParentID == tj2ID) tj.ParentID = newTjID;
4702 std::string type1Name,
4704 std::string type2Name)
4709 std::vector<int>
tmp;
4710 if (
id <= 0)
return tmp;
4711 unsigned int uid = id;
4713 if (type1Name ==
"T" && uid <= slc.
tjs.size() && type2Name ==
"P") {
4715 for (
auto& pfp : slc.
pfps) {
4716 if (pfp.ID <= 0)
continue;
4717 if (std::find(pfp.TjIDs.begin(), pfp.TjIDs.end(), id) != pfp.TjIDs.end())
4718 tmp.push_back(pfp.ID);
4723 if (type1Name ==
"P" && uid <= slc.
pfps.size() && (type2Name ==
"2S" || type2Name ==
"3S")) {
4725 auto& pfp = slc.
pfps[uid - 1];
4727 std::vector<int> ssid;
4728 for (
auto& ss : slc.
cots) {
4729 if (ss.ID <= 0)
continue;
4731 if (!shared.empty() && std::find(ssid.begin(), ssid.end(), ss.ID) == ssid.end())
4732 ssid.push_back(ss.ID);
4734 if (type2Name ==
"2S")
return ssid;
4735 for (
auto& ss3 : slc.
showers) {
4736 if (ss3.ID <= 0)
continue;
4738 if (!shared.empty() && std::find(tmp.begin(), tmp.end(), ss3.ID) == tmp.end())
4739 tmp.push_back(ss3.ID);
4744 if (type1Name ==
"2V" && uid <= slc.
vtxs.size() && type2Name ==
"T") {
4746 for (
auto& tj : slc.
tjs) {
4748 for (
unsigned short end = 0;
end < 2; ++
end) {
4749 if (tj.VtxID[
end] !=
id)
continue;
4750 if (std::find(tmp.begin(), tmp.end(), tj.ID) == tmp.end()) tmp.push_back(tj.ID);
4756 if (type1Name ==
"3V" && uid <= slc.
vtx3s.size() && type2Name ==
"P") {
4757 for (
auto& pfp : slc.
pfps) {
4758 if (pfp.ID == 0)
continue;
4759 for (
unsigned short end = 0;
end < 2; ++
end) {
4760 if (pfp.Vx3ID[
end] !=
id)
continue;
4762 if (std::find(tmp.begin(), tmp.end(), pfp.ID) == tmp.end()) tmp.push_back(pfp.ID);
4768 if (type1Name ==
"3V" && uid <= slc.
vtx3s.size() && type2Name ==
"T") {
4770 for (
auto& tj : slc.
tjs) {
4772 for (
unsigned short end = 0;
end < 2; ++
end) {
4773 if (tj.VtxID[
end] > 0 && tj.VtxID[
end] <= slc.
vtxs.size()) {
4774 auto& vx2 = slc.
vtxs[tj.VtxID[
end] - 1];
4775 if (vx2.Vx3ID !=
id)
continue;
4776 if (std::find(tmp.begin(), tmp.end(), tj.ID) == tmp.end()) tmp.push_back(tj.ID);
4783 if (type1Name ==
"3V" && uid <= slc.
vtx3s.size() && type2Name ==
"2V") {
4785 for (
auto& vx2 : slc.
vtxs) {
4786 if (vx2.ID == 0)
continue;
4787 if (vx2.Vx3ID ==
id) tmp.push_back(vx2.ID);
4792 if (type1Name ==
"3S" && uid <= slc.
showers.size() && type2Name ==
"T") {
4794 auto& ss3 = slc.
showers[uid - 1];
4795 if (ss3.ID == 0)
return tmp;
4796 for (
auto cid : ss3.CotIDs) {
4797 auto& ss = slc.
cots[cid - 1];
4798 if (ss.ID == 0)
continue;
4799 tmp.insert(tmp.end(), ss.TjIDs.begin(), ss.TjIDs.end());
4805 if (type1Name ==
"2S" && uid <= slc.
cots.size() && type2Name ==
"T") {
4807 auto& ss = slc.
cots[uid - 1];
4811 if (type1Name ==
"3S" && uid <= slc.
showers.size() && type2Name ==
"P") {
4813 auto& ss3 = slc.
showers[uid - 1];
4814 if (ss3.ID == 0)
return tmp;
4815 for (
auto cid : ss3.CotIDs) {
4816 auto& ss = slc.
cots[cid - 1];
4817 if (ss.ID == 0)
continue;
4818 for (
auto tid : ss.TjIDs) {
4819 auto& tj = slc.
tjs[tid - 1];
4821 if (!tj.AlgMod[
kMat3D])
continue;
4822 for (
auto& pfp : slc.
pfps) {
4823 if (pfp.ID <= 0)
continue;
4824 if (std::find(pfp.TjIDs.begin(), pfp.TjIDs.end(), tj.ID) == pfp.TjIDs.end())
continue;
4825 if (std::find(tmp.begin(), tmp.end(), pfp.ID) == tmp.end()) tmp.push_back(pfp.ID);
4832 if (type1Name ==
"T" && uid <= slc.
tjs.size() && type2Name ==
"2S") {
4834 for (
auto& ss : slc.
cots) {
4835 if (ss.ID == 0)
continue;
4836 if (std::find(ss.TjIDs.begin(), ss.TjIDs.end(), id) != ss.TjIDs.end()) tmp.push_back(ss.ID);
4841 if (type1Name ==
"T" && uid <= slc.
tjs.size() && type2Name ==
"3S") {
4843 for (
auto& ss : slc.
cots) {
4844 if (ss.ID == 0)
continue;
4845 if (std::find(ss.TjIDs.begin(), ss.TjIDs.end(), id) == ss.TjIDs.end())
continue;
4846 if (ss.SS3ID > 0) tmp.push_back(ss.SS3ID);
4857 unsigned int fromhit,
4859 unsigned short pass)
4865 float fromWire = fromHit.WireID().Wire;
4866 float fromTick = fromHit.PeakTime();
4867 float toWire = toHit.WireID().Wire;
4868 float toTick = toHit.PeakTime();
4870 bool success =
StartTraj(tj, fromWire, fromTick, toWire, toTick, tCTP, pass);
4871 if (!success)
return false;
4876 auto& tp = tj.
Pts[0];
4877 mf::LogVerbatim(
"TC") <<
"StartTraj T" << tj.
ID <<
" from " << (int)fromWire <<
":" 4878 << (
int)fromTick <<
" -> " << (int)toWire <<
":" << (
int)toTick
4879 <<
" StepDir " << tj.
StepDir <<
" dir " << tp.Dir[0] <<
" " << tp.Dir[1]
4880 <<
" ang " << tp.Ang <<
" AngleCode " << tp.AngleCode <<
" angErr " 4893 unsigned short pass)
4904 int fWire = std::nearbyint(fromWire);
4905 int tWire = std::nearbyint(toWire);
4906 if (tWire < fWire) { stepdir = -1; }
4907 else if (tWire == fWire) {
4909 if (toTick < fromTick) stepdir = -1;
4919 if (!
MakeBareTrajPoint(fromWire, fromTick, toWire, toTick, tCTP, tp))
return false;
4922 tj.
Pts.push_back(tp);
4927 auto& tp = tj.
Pts[0];
4928 mf::LogVerbatim(
"TC") <<
"StartTraj T" << tj.
ID <<
" from " << (int)fromWire <<
":" 4929 << (
int)fromTick <<
" -> " << (int)toWire <<
":" << (
int)toTick
4930 <<
" StepDir " << tj.
StepDir <<
" dir " << tp.
Dir[0] <<
" " << tp.
Dir[1]
4931 <<
" ang " << tp.
Ang <<
" AngleCode " << tp.
AngleCode <<
" angErr " 4939 std::pair<unsigned short, unsigned short>
GetSliceIndex(std::string typeName,
int uID)
4942 for (
unsigned short isl = 0; isl <
slices.size(); ++isl) {
4944 if (typeName ==
"T") {
4945 for (
unsigned short indx = 0; indx < slc.tjs.size(); ++indx) {
4946 if (slc.tjs[indx].UID == uID) {
return std::make_pair(isl, indx); }
4949 if (typeName ==
"P") {
4950 for (
unsigned short indx = 0; indx < slc.pfps.size(); ++indx) {
4951 if (slc.pfps[indx].UID == uID) {
return std::make_pair(isl, indx); }
4954 if (typeName ==
"2V") {
4955 for (
unsigned short indx = 0; indx < slc.vtxs.size(); ++indx) {
4956 if (slc.vtxs[indx].UID == uID) {
return std::make_pair(isl, indx); }
4959 if (typeName ==
"3V") {
4960 for (
unsigned short indx = 0; indx < slc.vtx3s.size(); ++indx) {
4961 if (slc.vtx3s[indx].UID == uID) {
return std::make_pair(isl, indx); }
4964 if (typeName ==
"2S") {
4965 for (
unsigned short indx = 0; indx < slc.cots.size(); ++indx) {
4966 if (slc.cots[indx].UID == uID) {
return std::make_pair(isl, indx); }
4969 if (typeName ==
"3S") {
4970 for (
unsigned short indx = 0; indx < slc.showers.size(); ++indx) {
4971 if (slc.showers[indx].UID == uID) {
return std::make_pair(isl, indx); }
4975 return std::make_pair(USHRT_MAX, USHRT_MAX);
4992 static double sum, sumx, sumy, sumx2, sumy2, sumxy;
4993 static unsigned short cnt;
4994 static std::vector<Point2_t> fitPts;
4995 static std::vector<double> fitWghts;
5012 if (inPtErr <= 0.)
return false;
5014 double wght = 1 / (inPtErr * inPtErr);
5016 sumx += wght * inPt[0];
5017 sumx2 += wght * inPt[0] * inPt[0];
5018 sumy += wght * inPt[1];
5019 sumy2 += wght * inPt[1] * inPt[1];
5020 sumxy += wght * inPt[0] * inPt[1];
5021 if (mode == 1)
return true;
5022 fitPts.push_back(inPt);
5023 fitWghts.push_back(wght);
5027 if (cnt < 2)
return false;
5029 double delta = sum * sumx2 - sumx * sumx;
5030 if (delta == 0.)
return false;
5031 double A = (sumx2 * sumy - sumx * sumxy) / delta;
5032 double B = (sumxy * sum - sumx * sumy) / delta;
5036 if (cnt == 2 || fitPts.empty())
return true;
5039 if (fitPts.size() != cnt)
return false;
5040 double ndof = cnt - 2;
5042 (sumy2 + A * A * sum + B * B * sumx2 - 2 * (A * sumy + B * sumxy - A * B * sumx)) / ndof;
5044 outVecErr[0] = sqrt(varnce * sumx2 / delta);
5045 outVecErr[1] = sqrt(varnce * sum / delta);
5053 for (
unsigned short ii = 0; ii < fitPts.size(); ++ii) {
5054 double arg = fitPts[ii][1] - A - B * fitPts[ii][0];
5055 sum += fitWghts[ii] * arg * arg;
5057 chiDOF = sum / ndof;
5070 if (strng ==
"instruct") {
5071 std::cout <<
"****** Unrecognized DebugConfig. Here are your options\n";
5072 std::cout <<
" 'C:T:P:W:Tick' where C = cryostat, T = TPC, W = wire, Tick (+/-5) to debug " 5073 "stepping (DUNE)\n";
5074 std::cout <<
" 'P:W:Tick' for single cryostat/TPC detectors (uB, LArIAT, etc)\n";
5075 std::cout <<
" 'WorkID <id> <slice index>' where <id> is a tj work ID (< 0) in slice <slice " 5076 "index> (default = 0)\n";
5077 std::cout <<
" 'Merge <CTP>' to debug trajectory merging\n";
5078 std::cout <<
" '2V <CTP>' to debug 2D vertex finding\n";
5079 std::cout <<
" '3V' to debug 3D vertex finding\n";
5080 std::cout <<
" 'VxMerge' to debug 2D vertex merging\n";
5081 std::cout <<
" 'JunkVx' to debug 2D junk vertex finder\n";
5082 std::cout <<
" 'PFP' to debug 3D matching and PFParticles\n";
5083 std::cout <<
" 'MVI <MVI> <MVI Iteration>' for detailed debugging of one PFP MatchVecIndex\n";
5084 std::cout <<
" 'DeltaRay' to debug delta ray tagging\n";
5085 std::cout <<
" 'Muon' to debug muon tagging\n";
5086 std::cout <<
" '2S <CTP>' to debug a 2D shower in CTP\n";
5087 std::cout <<
" 'Reco TPC <TPC>' to only reconstruct hits in the specified TPC\n";
5088 std::cout <<
" 'Reco Slice <ID>' to reconstruct all sub-slices in the recob::Slice with the " 5090 std::cout <<
" 'SubSlice <sub-slice index>' where <slice index> restricts output to the " 5091 "specified sub-slice index\n";
5092 std::cout <<
" 'Stitch' to debug PFParticle stitching between TPCs\n";
5093 std::cout <<
" 'Sum' or 'Summary' to print a debug summary report\n";
5094 std::cout <<
" 'Dump <WorkID>' or 'Dump <UID>' to print all TPs in the trajectory to " 5095 "tcdump<UID>.csv\n";
5096 std::cout <<
" Note: Algs with debug printing include HamVx, HamVx2, SplitTjCVx, Comp3DVx, " 5097 "Comp3DVxIG, VtxHitsSwap\n";
5098 std::cout <<
" Set SkipAlgs: [\"bogusText\"] to print a list of algorithm names\n";
5103 if (strng.find(
"3V") != std::string::npos) {
5108 if (strng.find(
"3S") != std::string::npos) {
5113 if (strng.find(
"VxMerge") != std::string::npos) {
5118 if (strng.find(
"JunkVx") != std::string::npos) {
5123 if (strng.find(
"DeltaRay") != std::string::npos) {
5128 if (strng.find(
"Muon") != std::string::npos) {
5133 if (strng.find(
"Stitch") != std::string::npos) {
5138 if (strng.find(
"HamVx") != std::string::npos) {
5143 if (strng.find(
"HamVx2") != std::string::npos) {
5148 if (strng.find(
"Sum") != std::string::npos) {
5154 std::vector<std::string> words;
5155 boost::split(words, strng, boost::is_any_of(
" :"), boost::token_compress_on);
5156 if (words.size() == 5) {
5169 if (words[0] ==
"PFP" || words[0] ==
"MVI") {
5173 if (words.size() > 2) {
5175 if (words.size() == 3)
debug.
MVI_Iter = std::stoi(words[2]);
5179 if (words.size() == 2 && words[0] ==
"Dump") {
5186 if (words.size() > 1 && words[0] ==
"WorkID") {
5191 if (words.size() > 2)
debug.
Slice = std::stoi(words[2]);
5197 if (words.size() == 3 && words[0] ==
"Reco" && words[1] ==
"TPC") {
5200 std::cout <<
"Reconstructing only in TPC " <<
tcc.
recoTPC <<
"\n";
5203 if (words.size() == 3 && words[0] ==
"Reco" && words[1] ==
"Slice") {
5205 std::cout <<
"Reconstructing Slice " <<
tcc.
recoSlice <<
"\n";
5208 if (words.size() == 3) {
5220 if (words.size() == 2 && words[0] ==
"Merge") {
5226 if (words.size() == 2 && words[0] ==
"2V") {
5232 if (words.size() == 2 && words[0] ==
"2S") {
5239 if (words.size() == 2 && words[0] ==
"SubSlice") {
5254 for (
auto const& slc :
slices) {
5255 for (
auto& tj : slc.tjs) {
5258 std::ofstream outfile;
5260 outfile.open(fname, std::ios::out | std::ios::trunc);
5261 outfile <<
"Dump trajectory T" << tj.UID <<
" WorkID " << tj.WorkID;
5262 outfile <<
" ChgRMS " << std::setprecision(2) << tj.ChgRMS;
5264 outfile <<
"Wire, Chg T" << tj.UID
5265 <<
", totChg, Tick, Delta, NTPsFit, Ang, ChiDOF, KinkSig, HitPosErr\n";
5266 for (
unsigned short ipt = tj.EndPt[0]; ipt <= tj.EndPt[1]; ++ipt) {
5267 auto& tp = tj.Pts[ipt];
5268 outfile << std::fixed;
5269 outfile << std::setprecision(0) << std::nearbyint(tp.Pos[0]);
5270 outfile <<
"," << (int)tp.Chg;
5273 for (
auto iht : tp.Hits) {
5275 totChg +=
hit.Integral();
5277 outfile <<
"," << (int)totChg;
5278 outfile <<
"," << std::setprecision(0) << std::nearbyint(tp.Pos[1] /
tcc.
unitsPerTick);
5279 outfile <<
"," << std::setprecision(2) << tp.Delta;
5280 outfile <<
"," << tp.NTPsFit;
5281 outfile <<
"," << std::setprecision(3) << tp.Ang;
5282 outfile <<
"," << std::setprecision(2) << tp.FitChi;
5283 outfile <<
"," << std::setprecision(2) << tp.KinkSig;
5284 outfile <<
"," << std::setprecision(2) << sqrt(tp.HitPosErr2);
5288 std::cout <<
"Points on T" << tj.UID <<
" dumped to " << fname <<
"\n";
5300 std::cout <<
"*** TrajCluster debug mode configuration in";
5301 std::cout <<
" CTP=";
5302 if (
debug.
CTP == UINT_MAX) { std::cout <<
"NA"; }
5311 std::cout <<
" Hit=";
5312 if (
debug.
Hit == UINT_MAX) { std::cout <<
"NA"; }
5316 std::cout <<
" WorkID=";
5321 std::cout <<
" Slice=";
5322 if (
debug.
Slice == -1) { std::cout <<
"All"; }
5327 std::cout <<
"*** tcc.dbg modes:";
5331 if (
tcc.
dbg2V) std::cout <<
" dbg2V";
5332 if (
tcc.
dbg2S) std::cout <<
" dbg2S";
5336 if (
tcc.
dbg3V) std::cout <<
" dbg3V";
5344 std::cout <<
"*** Using algs:";
5345 unsigned short cnt = 0;
5346 for (
unsigned short ib = 0; ib <
AlgBitNames.size(); ++ib) {
5349 if (cnt % 10 == 0) std::cout <<
"\n ";
5354 std::cout <<
"*** Skipping algs:";
5356 for (
unsigned short ib = 0; ib <
AlgBitNames.size(); ++ib) {
5359 if (cnt % 10 == 0) std::cout <<
"\n ";
5375 for (
size_t isl = 0; isl <
slices.size(); ++isl) {
5377 auto const& slc =
slices[isl];
5378 if (!slc.vtx3s.empty()) prt3V =
true;
5379 if (!slc.vtxs.empty()) prt2V =
true;
5380 if (!slc.tjs.empty()) prtT =
true;
5381 if (!slc.pfps.empty()) prtP =
true;
5382 if (!slc.showers.empty()) prtS3 =
true;
5385 myprt <<
"Debug report from caller " << someText <<
"\n";
5386 myprt <<
" 'prodID' = <sliceID>:<subSliceIndex>:<productID>/<productUID>\n";
5388 myprt <<
"************ Showers ************\n";
5389 myprt <<
" prodID Vtx parUID ___ChgPos____ ______Dir_____ ____posInPln____ " 5390 "___projInPln____ 2D shower UIDs\n";
5391 for (
size_t isl = 0; isl <
slices.size(); ++isl) {
5393 auto const& slc =
slices[isl];
5394 if (slc.showers.empty())
continue;
5395 for (
auto& ss3 : slc.showers)
5400 bool printHeader =
true;
5401 for (
size_t isl = 0; isl <
slices.size(); ++isl) {
5403 auto const& slc =
slices[isl];
5404 if (slc.pfps.empty())
continue;
5405 for (
auto& pfp : slc.pfps)
5406 PrintP(myprt, pfp, printHeader);
5410 bool printHeader =
true;
5411 myprt <<
"****** 3D vertices " 5412 "******************************************__2DVtx_UID__*******\n";
5413 myprt <<
" prodID Cstat TPC X Y Z XEr YEr " 5414 "ZEr pln0 pln1 pln2 Wire score Prim? Nu? nTru";
5415 myprt <<
" ___________2D_Pos____________ _____Tj UIDs________\n";
5416 for (
size_t isl = 0; isl <
slices.size(); ++isl) {
5418 auto const& slc =
slices[isl];
5419 if (slc.vtx3s.empty())
continue;
5420 for (
auto& vx3 : slc.vtx3s)
5421 Print3V(detProp, myprt, vx3, printHeader);
5425 bool printHeader =
true;
5426 myprt <<
"************ 2D vertices ************\n";
5427 myprt <<
" prodID CTP wire err tick err ChiDOF NTj Pass " 5428 " Topo ChgFrac Score v3D Tj UIDs\n";
5429 for (
size_t isl = 0; isl <
slices.size(); ++isl) {
5431 auto const& slc =
slices[isl];
5432 if (slc.vtxs.empty())
continue;
5433 for (
auto& vx2 : slc.vtxs)
5434 Print2V(myprt, vx2, printHeader);
5438 bool printHeader =
true;
5439 for (
size_t isl = 0; isl <
slices.size(); ++isl) {
5441 auto const& slc =
slices[isl];
5442 if (slc.tjs.empty())
continue;
5443 for (
auto& tj : slc.tjs)
5444 PrintT(myprt, tj, printHeader);
5452 if (pfp.
ID <= 0)
return;
5454 myprt <<
"************ PFParticles ************\n";
5455 myprt <<
" prodID sVx _____sPos____ CS _______sDir______ ____sdEdx_____ eVx " 5456 "_____ePos____ CS ____edEdx_____ MVI MCSMom Len nTP3 nSec SLk? PDG Par \n";
5457 printHeader =
false;
5460 if (sIndx.first == USHRT_MAX)
return;
5461 auto const& slc =
slices[sIndx.first];
5465 myprt << std::setw(12) << str;
5467 for (
unsigned short end = 0;
end < 2; ++
end) {
5470 myprt << std::setw(6) << str;
5471 myprt << std::fixed <<
std::right << std::setprecision(0);
5473 myprt << std::setw(5) << pos[0];
5474 myprt << std::setw(5) << pos[1];
5475 myprt << std::setw(5) << pos[2];
5483 myprt << std::fixed <<
std::right << std::setprecision(2);
5485 myprt << std::setw(6) <<
dir[0];
5486 myprt << std::setw(6) << dir[1];
5487 myprt << std::setw(6) << dir[2];
5489 for (
auto& dedx : pfp.
dEdx[
end]) {
5490 if (dedx < 50) { myprt << std::setw(5) << std::setprecision(1) << dedx; }
5492 myprt << std::setw(5) << std::setprecision(0) << dedx;
5495 if (pfp.
dEdx[end].size() < 3) {
5496 for (
size_t i = 0; i < 3 - pfp.
dEdx[
end].size(); ++i) {
5497 myprt << std::setw(6) <<
' ';
5501 myprt << std::setw(6) << pfp.
MVI;
5504 float length =
Length(pfp);
5505 if (length < 100) { myprt << std::setw(5) << std::setprecision(1) << length; }
5507 myprt << std::setw(5) << std::setprecision(0) << length;
5509 myprt << std::setw(5) << pfp.
TP3Ds.size();
5512 myprt << std::setw(5) << pfp.
PDGCode;
5514 if (!pfp.
TjIDs.empty()) {
5515 if (pfp.
TjUIDs.empty()) {
5517 for (
auto tjid : pfp.
TjIDs)
5518 myprt <<
" TU" << slc.tjs[tjid - 1].UID;
5522 for (
auto tjuid : pfp.
TjUIDs)
5523 myprt <<
" TU" << tjuid;
5528 for (
auto dtruid : pfp.
DtrUIDs)
5529 myprt <<
" PU" << dtruid;
5541 if (vx3.
ID <= 0)
return;
5543 if (sIndx.first == USHRT_MAX)
return;
5544 auto const& slc =
slices[sIndx.first];
5547 <<
"****** 3D vertices ******************************************__2DVtx_UID__*******\n";
5548 myprt <<
" prodID Cstat TPC X Y Z pln0 pln1 pln2 Wire score " 5550 myprt <<
" ___________2D_Pos____________ _____Tj UIDs________\n";
5551 printHeader =
false;
5554 myprt <<
std::right << std::setw(12) << std::fixed << str;
5555 myprt << std::setprecision(0);
5561 for (
auto vx2id : vx3.
Vx2ID) {
5571 unsigned short nTruMatch = 0;
5572 for (
unsigned short ipl = 0; ipl < slc.nPlanes; ++ipl) {
5573 if (vx3.
Vx2ID[ipl] == 0)
continue;
5574 unsigned short iv2 = vx3.
Vx2ID[ipl] - 1;
5577 myprt <<
std::right << std::setw(6) << std::setprecision(1) << vx3.
Score;
5578 myprt << std::setw(6) << vx3.
Primary;
5579 myprt << std::setw(4) << vx3.
Neutrino;
5580 myprt <<
std::right << std::setw(5) << nTruMatch;
5582 for (
unsigned short plane = 0; plane < slc.nPlanes; ++plane) {
5586 if (vx3.
Wire == -2) {
5588 for (
unsigned short end = 0;
end < 2; ++
end) {
5589 for (
auto& pfp : slc.pfps) {
5590 if (pfp.Vx3ID[
end] == vx3.
ID) {
5591 for (
auto tjID : pfp.TjIDs) {
5592 auto& tj = slc.tjs[tjID - 1];
5593 myprt <<
" T" << tj.UID;
5600 auto vxtjs =
GetAssns(slc,
"3V", vx3.
ID,
"T");
5601 for (
auto tjid : vxtjs) {
5602 auto& tj = slc.tjs[tjid - 1];
5603 myprt <<
" TU" << tj.UID;
5613 if (vx2.
ID <= 0)
return;
5616 if (sIndx.first == USHRT_MAX)
return;
5617 auto const& slc =
slices[sIndx.first];
5619 myprt <<
"************ 2D vertices ************\n";
5620 myprt <<
" prodID CTP wire err tick err ChiDOF NTj Pass Topo ChgFrac Score " 5622 printHeader =
false;
5625 myprt <<
std::right << std::setw(12) << std::fixed << str;
5627 myprt <<
std::right << std::setw(8) << std::setprecision(0) << std::nearbyint(vx2.
Pos[0]);
5628 myprt <<
std::right << std::setw(5) << std::setprecision(1) << vx2.
PosErr[0];
5629 myprt <<
std::right << std::setw(8) << std::setprecision(0)
5637 myprt <<
std::right << std::setw(6) << std::setprecision(1) << vx2.
Score;
5639 if (vx2.
Vx3ID > 0) v3id = slc.vtx3s[vx2.
Vx3ID - 1].UID;
5643 for (
unsigned short ii = 0; ii < slc.tjs.size(); ++ii) {
5644 auto const& tj = slc.tjs[ii];
5645 if (tj.AlgMod[
kKilled])
continue;
5646 for (
unsigned short end = 0;
end < 2; ++
end) {
5647 if (tj.VtxID[
end] != (
short)vx2.
ID)
continue;
5654 for (
unsigned short ib = 1; ib <
VtxBitNames.size(); ++ib)
5664 if (ss3.
ID <= 0)
return;
5666 if (sIndx.first == USHRT_MAX)
return;
5667 auto const& slc =
slices[sIndx.first];
5671 myprt << std::fixed << std::setw(12) << str;
5674 myprt << std::setw(6) << str;
5675 for (
unsigned short xyz = 0; xyz < 3; ++xyz)
5676 myprt << std::setprecision(0) << std::setw(5) << ss3.
ChgPos[xyz];
5677 for (
unsigned short xyz = 0; xyz < 3; ++xyz)
5678 myprt << std::setprecision(2) << std::setw(5) << ss3.
Dir[xyz];
5679 std::vector<float> projInPlane(slc.nPlanes);
5680 for (
unsigned short plane = 0; plane < slc.nPlanes; ++plane) {
5684 projInPlane[plane] = tp.Delta;
5686 for (
unsigned short plane = 0; plane < slc.nPlanes; ++plane) {
5687 myprt << std::setprecision(2) << std::setw(5) << projInPlane[plane];
5689 for (
auto cid : ss3.
CotIDs) {
5690 auto& ss = slc.cots[cid - 1];
5692 myprt << std::setw(5) << str;
5694 if (ss3.
NeedsUpdate) myprt <<
" *** Needs update";
5702 if (tj.
ID <= 0)
return;
5705 myprt <<
"************ Trajectories ************\n";
5706 myprt <<
"Tj AngleCode-EndFlag decoder (EF): <AngleCode> + <reason for stopping>";
5707 myprt <<
" (B=Bragg Peak, V=Vertex, A=AngleKink, C=ChargeKink, T=Trajectory)\n";
5708 myprt <<
" prodID CTP Pass Pts W:T Ang EF AveQ W:T Ang EF AveQ " 5709 "Chg(k) chgRMS Mom __Vtx__ PDG eLike Par Pri NuPar WorkID \n";
5710 printHeader =
false;
5713 if (sIndx.first == USHRT_MAX)
return;
5714 auto const& slc =
slices[sIndx.first];
5716 myprt << std::fixed << std::setw(12) << str;
5717 myprt << std::setw(6) << tj.
CTP;
5718 myprt << std::setw(5) << tj.
Pass;
5719 myprt << std::setw(5) << tj.
EndPt[1] - tj.
EndPt[0] + 1;
5720 unsigned short endPt0 = tj.
EndPt[0];
5721 auto& tp0 = tj.
Pts[endPt0];
5723 if (itick < 0) itick = 0;
5724 myprt << std::setw(6) << (int)(tp0.Pos[0] + 0.5) <<
":" << itick;
5725 if (itick < 10) { myprt <<
" "; }
5726 if (itick < 100) { myprt <<
" "; }
5727 if (itick < 1000) { myprt <<
" "; }
5728 myprt << std::setw(6) << std::setprecision(2) << tp0.Ang;
5729 myprt << std::setw(2) << tp0.AngleCode;
5743 myprt << std::setw(5) << (int)tp0.AveChg;
5744 unsigned short endPt1 = tj.
EndPt[1];
5745 auto& tp1 = tj.
Pts[endPt1];
5747 myprt << std::setw(6) << (int)(tp1.Pos[0] + 0.5) <<
":" << itick;
5748 if (itick < 10) { myprt <<
" "; }
5749 if (itick < 100) { myprt <<
" "; }
5750 if (itick < 1000) { myprt <<
" "; }
5751 myprt << std::setw(6) << std::setprecision(2) << tp1.Ang;
5752 myprt << std::setw(2) << tp1.AngleCode;
5766 myprt << std::setw(5) << (int)tp1.AveChg;
5767 myprt << std::setw(7) << std::setprecision(1) << tj.
TotChg / 1000;
5768 myprt << std::setw(7) << std::setprecision(2) << tj.
ChgRMS;
5769 myprt << std::setw(5) << tj.
MCSMom;
5771 if (tj.
VtxID[0] > 0) vxid = slc.vtxs[tj.
VtxID[0] - 1].UID;
5772 myprt << std::setw(4) << vxid;
5774 if (tj.
VtxID[1] > 0) vxid = slc.vtxs[tj.
VtxID[1] - 1].UID;
5775 myprt << std::setw(4) << vxid;
5776 myprt << std::setw(5) << tj.
PDGCode;
5778 myprt << std::setw(5) << tj.
ParentID;
5779 myprt << std::setw(5) <<
PrimaryID(slc, tj);
5781 myprt << std::setw(7) << tj.
WorkID;
5782 for (
unsigned short ib = 0; ib <
AlgBitNames.size(); ++ib)
5791 std::string someText,
5801 if (!slc.
vtx3s.empty()) {
5805 <<
"****** 3D vertices ******************************************__2DVtx_ID__*******\n";
5807 <<
" Vtx Cstat TPC X Y Z XEr YEr ZEr pln0 pln1 pln2 Wire " 5808 "score Prim? Nu? nTru";
5809 myprt <<
" ___________2D_Pos____________ _____Tjs________\n";
5810 for (
unsigned short iv = 0; iv < slc.
vtx3s.size(); ++iv) {
5811 if (slc.
vtx3s[iv].ID == 0)
continue;
5815 myprt <<
std::right << std::setw(5) << std::fixed << vid;
5816 myprt << std::setprecision(1);
5829 unsigned short nTruMatch = 0;
5830 for (
unsigned short ipl = 0; ipl < slc.
nPlanes; ++ipl) {
5831 if (vx3.
Vx2ID[ipl] == 0)
continue;
5832 unsigned short iv2 = vx3.
Vx2ID[ipl] - 1;
5835 myprt <<
std::right << std::setw(6) << std::setprecision(1) << vx3.
Score;
5836 myprt << std::setw(6) << vx3.
Primary;
5837 myprt << std::setw(4) << vx3.
Neutrino;
5838 myprt <<
std::right << std::setw(5) << nTruMatch;
5840 for (
unsigned short plane = 0; plane < slc.
nPlanes; ++plane) {
5844 if (vx3.
Wire == -2) {
5846 for (
auto& pfp : slc.
pfps) {
5847 if (pfp.Vx3ID[0] == slc.
vtx3s[iv].ID) {
5848 for (
auto& tjID : pfp.TjIDs)
5849 myprt <<
" t" << tjID;
5851 if (pfp.Vx3ID[1] == slc.
vtx3s[iv].ID) {
5852 for (
auto& tjID : pfp.TjIDs)
5853 myprt <<
" t" << tjID;
5858 auto vxtjs =
GetAssns(slc,
"3V", vx3.
ID,
"T");
5859 for (
auto tjid : vxtjs)
5860 myprt <<
" t" << tjid;
5865 if (!slc.
vtxs.empty()) {
5866 bool foundOne =
false;
5867 for (
unsigned short iv = 0; iv < slc.
vtxs.size(); ++iv) {
5868 auto& vx2 = slc.
vtxs[iv];
5870 if (vx2.NTraj == 0)
continue;
5875 myprt << someText <<
"************ 2D vertices ************\n";
5877 <<
" ID CTP wire err tick err ChiDOF NTj Pass Topo ChgFrac Score v3D " 5879 for (
auto& vx2 : slc.
vtxs) {
5880 if (vx2.ID == 0)
continue;
5884 myprt <<
std::right << std::setw(5) << std::fixed << vid;
5885 myprt <<
std::right << std::setw(6) << vx2.CTP;
5886 myprt <<
std::right << std::setw(8) << std::setprecision(0)
5887 << std::nearbyint(vx2.Pos[0]);
5888 myprt <<
std::right << std::setw(5) << std::setprecision(1) << vx2.PosErr[0];
5889 myprt <<
std::right << std::setw(8) << std::setprecision(0)
5891 myprt <<
std::right << std::setw(5) << std::setprecision(1)
5893 myprt <<
std::right << std::setw(7) << vx2.ChiDOF;
5894 myprt <<
std::right << std::setw(5) << vx2.NTraj;
5895 myprt <<
std::right << std::setw(5) << vx2.Pass;
5896 myprt <<
std::right << std::setw(6) << vx2.Topo;
5897 myprt <<
std::right << std::setw(9) << std::setprecision(2) << vx2.TjChgFrac;
5898 myprt <<
std::right << std::setw(6) << std::setprecision(1) << vx2.Score;
5899 myprt <<
std::right << std::setw(5) << vx2.Vx3ID;
5902 for (
unsigned short ii = 0; ii < slc.
tjs.size(); ++ii) {
5903 auto const& aTj = slc.
tjs[ii];
5905 if (aTj.AlgMod[
kKilled])
continue;
5906 for (
unsigned short end = 0;
end < 2; ++
end) {
5907 if (aTj.VtxID[
end] != (
short)vx2.ID)
continue;
5913 for (
unsigned short ib = 1; ib <
VtxBitNames.size(); ++ib)
5914 if (vx2.Stat[ib]) myprt <<
" " <<
VtxBitNames[ib];
5921 if (slc.
tjs.empty()) {
5922 mf::LogVerbatim(
"TC") << someText <<
" No allTraj trajectories to print";
5928 if (itj == USHRT_MAX) {
5930 myprt <<
"Tj AngleCode-EndFlag (EF) decoder: <AngleCode> + <reason for stopping>";
5931 myprt <<
" (B=Bragg Peak, V=Vertex, A=AngleKink, C=ChargeKink, T=Trajectory)\n";
5932 std::vector<unsigned int>
tmp;
5934 <<
" UID CTP Pass Pts W:T Ang EF AveQ W:T Ang EF AveQ Chg(k) " 5935 "chgRMS Mom SDr __Vtx__ PDG Par Pri NuPar WorkID \n";
5936 for (
unsigned short ii = 0; ii < slc.
tjs.size(); ++ii) {
5937 auto& aTj = slc.
tjs[ii];
5939 myprt << someText <<
" ";
5945 myprt << std::fixed << std::setw(5) << tid;
5946 myprt << std::setw(6) << aTj.CTP;
5947 myprt << std::setw(5) << aTj.Pass;
5948 myprt << std::setw(5) << aTj.EndPt[1] - aTj.EndPt[0] + 1;
5949 unsigned short endPt0 = aTj.EndPt[0];
5950 auto& tp0 = aTj.Pts[endPt0];
5952 if (itick < 0) itick = 0;
5953 myprt << std::setw(6) << (int)(tp0.Pos[0] + 0.5) <<
":" << itick;
5954 if (itick < 10) { myprt <<
" "; }
5955 if (itick < 100) { myprt <<
" "; }
5956 if (itick < 1000) { myprt <<
" "; }
5957 myprt << std::setw(6) << std::setprecision(2) << tp0.Ang;
5958 myprt << std::setw(2) << tp0.AngleCode;
5959 if (aTj.EndFlag[0][
kBragg]) { myprt <<
"B"; }
5960 else if (aTj.EndFlag[0][
kAtVtx]) {
5963 else if (aTj.EndFlag[0][
kAtKink]) {
5966 else if (aTj.EndFlag[0][
kAtTj]) {
5972 myprt << std::setw(5) << (int)tp0.AveChg;
5973 unsigned short endPt1 = aTj.EndPt[1];
5974 auto& tp1 = aTj.Pts[endPt1];
5976 myprt << std::setw(6) << (int)(tp1.Pos[0] + 0.5) <<
":" << itick;
5977 if (itick < 10) { myprt <<
" "; }
5978 if (itick < 100) { myprt <<
" "; }
5979 if (itick < 1000) { myprt <<
" "; }
5980 myprt << std::setw(6) << std::setprecision(2) << tp1.Ang;
5981 myprt << std::setw(2) << tp1.AngleCode;
5982 if (aTj.EndFlag[1][
kBragg]) { myprt <<
"B"; }
5983 else if (aTj.EndFlag[1][
kAtVtx]) {
5989 myprt << std::setw(5) << (int)tp1.AveChg;
5990 myprt << std::setw(7) << std::setprecision(1) << aTj.TotChg / 1000;
5991 myprt << std::setw(7) << std::setprecision(2) << aTj.ChgRMS;
5992 myprt << std::setw(5) << aTj.MCSMom;
5993 myprt << std::setw(4) << aTj.StepDir;
5994 myprt << std::setw(4) << aTj.VtxID[0];
5995 myprt << std::setw(4) << aTj.VtxID[1];
5996 myprt << std::setw(5) << aTj.PDGCode;
5997 myprt << std::setw(5) << aTj.ParentID;
5998 myprt << std::setw(5) <<
PrimaryID(slc, aTj);
6000 myprt << std::setw(7) << aTj.WorkID;
6001 for (
unsigned short ib = 0; ib <
AlgBitNames.size(); ++ib)
6002 if (aTj.AlgMod[ib]) myprt <<
" " <<
AlgBitNames[ib];
6008 if (itj > slc.
tjs.size() - 1)
return;
6010 auto const& aTj = slc.
tjs[itj];
6012 mf::LogVerbatim(
"TC") <<
"Print slc.tjs[" << itj <<
"] Vtx[0] " << aTj.VtxID[0] <<
" Vtx[1] " 6015 for (
unsigned short ib = 0; ib <
AlgBitNames.size(); ++ib)
6016 if (aTj.AlgMod[ib]) myprt <<
" " <<
AlgBitNames[ib];
6020 if (ipt == USHRT_MAX) {
6022 for (
unsigned short ii = 0; ii < aTj.Pts.size(); ++ii)
6023 PrintTP(someText, slc, ii, aTj.StepDir, aTj.Pass, aTj.Pts[ii]);
6027 PrintTP(someText, slc, ipt, aTj.StepDir, aTj.Pass, aTj.Pts[ipt]);
6035 unsigned short tPoint)
6039 if (tPoint == USHRT_MAX) {
6042 myprt << someText <<
" ";
6043 myprt <<
"Work: UID " << tj.
UID <<
" CTP " << tj.
CTP <<
" StepDir " << tj.
StepDir 6045 <<
" nPts " << tj.
Pts.size() <<
" EndPts " << tj.
EndPt[0] <<
" " << tj.
EndPt[1];
6046 myprt <<
" MCSMom " << tj.
MCSMom;
6048 myprt <<
" AlgMods:";
6049 for (
unsigned short ib = 0; ib <
AlgBitNames.size(); ++ib)
6054 myprt << someText <<
" ";
6055 myprt <<
"slcID " << slc.
ID <<
" T" << tj.
ID <<
" uT" << tj.
UID <<
" WorkID " << tj.
WorkID 6057 <<
" " << tj.
VtxID[1] <<
" nPts " << tj.
Pts.size() <<
" EndPts " << tj.
EndPt[0] <<
" " 6059 myprt <<
" MCSMom " << tj.
MCSMom;
6061 myprt <<
" AlgMods:";
6062 for (
unsigned short ib = 0; ib <
AlgBitNames.size(); ++ib)
6066 for (
unsigned short ipt = 0; ipt < tj.
Pts.size(); ++ipt)
6070 for (
unsigned short ic = 0; ic < slc.
cots.size(); ++ic) {
6071 if (slc.
cots[ic].TjIDs.empty())
continue;
6073 if (slc.
cots[ic].ShowerTjID != tj.
ID)
continue;
6076 myprt <<
"cots index " << ic <<
" ";
6077 myprt << someText <<
" Envelope";
6078 if (ss.
Envelope.empty()) { myprt <<
" NA"; }
6081 myprt <<
" " << (int)vtx[0] <<
":" << (
int)(vtx[1] /
tcc.
unitsPerTick);
6083 myprt <<
" Energy " << (int)ss.
Energy;
6084 myprt <<
" Area " << std::fixed << std::setprecision(1) << (int)ss.
EnvelopeArea 6086 myprt <<
"\nInShower TjIDs";
6088 myprt <<
" " << tjID;
6092 myprt <<
"NearTjIDs";
6094 myprt <<
" " << tjID;
6098 myprt <<
"Angle " << std::fixed << std::setprecision(2) << ss.
Angle <<
" +/- " 6100 myprt <<
" AspectRatio " << std::fixed << std::setprecision(2) << ss.
AspectRatio;
6101 myprt <<
" DirectionFOM " << std::fixed << std::setprecision(2) << ss.
DirectionFOM;
6104 myprt <<
" No parent";
6106 myprt <<
" TruParentID " << ss.
TruParentID <<
" SS3ID " << ss.
SS3ID <<
"\n";
6107 if (ss.
NeedsUpdate) myprt <<
"*********** This shower needs to be updated ***********";
6108 myprt <<
"................................................";
6114 if (tPoint > tj.
Pts.size() - 1) {
6115 mf::LogVerbatim(
"TC") <<
"Can't print non-existent traj point " << tPoint;
6126 <<
" TRP CTP Ind Stp Delta RMS Ang C Err Dir0 Dir1 Q " 6127 " AveQ Pull FitChi NTPF KinkSig Hits ";
6135 unsigned short pass,
6139 myprt << someText <<
" TRP" << std::fixed;
6141 if (dir > 0) { myprt <<
"+"; }
6145 myprt << std::setw(6) << tp.
CTP;
6146 myprt << std::setw(5) << ipt;
6147 myprt << std::setw(5) << tp.
Step;
6148 myprt << std::setw(6) << std::setprecision(2) << tp.
Delta;
6149 myprt << std::setw(6) << std::setprecision(2) << tp.
DeltaRMS;
6150 myprt << std::setw(6) << std::setprecision(2) << tp.
Ang;
6152 myprt << std::setw(6) << std::setprecision(2) << tp.
AngErr;
6153 myprt << std::setw(6) << std::setprecision(2) << tp.
Dir[0];
6154 myprt << std::setw(6) << std::setprecision(2) << tp.
Dir[1];
6155 myprt << std::setw(7) << (int)tp.
Chg;
6156 myprt << std::setw(8) << (int)tp.
AveChg;
6157 myprt << std::setw(6) << std::setprecision(1) << tp.
ChgPull;
6158 myprt << std::setw(7) << tp.
FitChi;
6159 myprt << std::setw(6) << tp.
NTPsFit;
6160 myprt << std::setw(7) << std::setprecision(3) << tp.
KinkSig;
6162 if (tp.
Hits.size() > 16) {
6164 myprt <<
" " << tp.
Hits.size() <<
" shower hits";
6167 for (
unsigned short ii = 0; ii < tp.
Hits.size(); ++ii) {
6168 unsigned int iht = tp.
Hits[ii];
6170 myprt <<
" " <<
hit.WireID().Wire <<
":" << (int)
hit.PeakTime();
6178 myprt <<
"T" << slc.
slHits[iht].InTraj;
6180 if (tp.
InPFP > 0) myprt <<
" inP" << tp.
InPFP;
6190 std::string str =
"";
6191 for (
unsigned short ib = 0; ib < 8; ++ib) {
6200 if (ib ==
kEnvFlag) str +=
" Flag";
6211 myprt <<
" PFP sVx ________sPos_______ EF _______sDir______ ____sdEdx_____ eVx " 6212 "________ePos_______ EF _______eDir______ ____edEdx____ Len nTp3 MCSMom ShLike? " 6217 myprt << std::setw(5) << pid;
6219 for (
unsigned short end = 0;
end < 2; ++
end) {
6220 myprt << std::setw(4) << pfp.
Vx3ID[
end];
6221 myprt << std::fixed <<
std::right << std::setprecision(1);
6223 myprt << std::setw(7) << pos[0];
6224 myprt << std::setw(7) << pos[1];
6225 myprt << std::setw(7) << pos[2];
6233 myprt << std::setw(6) << ef;
6234 myprt << std::fixed <<
std::right << std::setprecision(2);
6236 myprt << std::setw(6) <<
dir[0];
6237 myprt << std::setw(6) << dir[1];
6238 myprt << std::setw(6) << dir[2];
6239 for (
auto& dedx : pfp.
dEdx[
end]) {
6240 if (dedx < 50) { myprt << std::setw(5) << std::setprecision(1) << dedx; }
6242 myprt << std::setw(5) << std::setprecision(0) << dedx;
6245 if (pfp.
dEdx[end].size() < 3) {
6246 for (
size_t i = 0; i < 3 - pfp.
dEdx[
end].size(); ++i) {
6247 myprt << std::setw(6) <<
' ';
6252 float length =
Length(pfp);
6253 if (length < 100) { myprt << std::setw(5) << std::setprecision(1) << length; }
6255 myprt << std::setw(5) << std::setprecision(0) << length;
6257 myprt << std::setw(5) << std::setprecision(2) << pfp.
TP3Ds.size();
6260 myprt << std::setw(5) << pfp.
PDGCode;
6264 if (!pfp.
TjIDs.empty()) {
6265 for (
auto& tjID : pfp.
TjIDs)
6266 myprt <<
" T" << tjID;
6270 for (
auto& dtrUID : pfp.
DtrUIDs)
6271 myprt <<
" P" << dtrUID;
6278 if (slc.
pfps.empty())
return;
6283 <<
" PFP sVx ________sPos_______ ______sDir______ ______sdEdx_____ eVx " 6284 "________ePos_______ ______eDir______ ______edEdx_____ BstPln PDG TruPDG Par Prim E*P\n";
6285 bool printHeader =
true;
6286 for (
auto& pfp : slc.
pfps) {
6287 PrintPFP(someText, slc, pfp, printHeader);
6288 printHeader =
false;
6296 if (end > 1)
return "Invalid end";
6299 for (
unsigned short ib = 0; ib <
EndFlagNames.size(); ++ib) {
6310 if (first) tmp =
" none";
6317 if (end > 1)
return "Invalid end";
6320 for (
unsigned short ib = 0; ib <
EndFlagNames.size(); ++ib) {
6361 unsigned int wire = 0;
6362 if (pos[0] > -0.4) wire = std::nearbyint(pos[0]);
Expect tracks entering from the front face. Don't create neutrino PFParticles.
void PrintAll(detinfo::DetectorPropertiesData const &detProp, std::string someText)
float HitsPosTime(const TCSlice &slc, const std::vector< unsigned int > &hitsInMultiplet, float &sum, HitStatus_t hitRequest)
void CheckTrajBeginChg(TCSlice &slc, unsigned short itj)
float AveChg
Calculated using ALL hits.
details::range_type< T > Iterate() const
Initializes the specified ID with the ID of the first cryostat.
short MCSMom(const TCSlice &slc, const std::vector< int > &tjIDs)
MaybeLogger_< ELseverityLevel::ELsev_info, true > LogVerbatim
std::vector< Trajectory > tjs
vector of all trajectories in each plane
void FitTraj(TCSlice const &slc, Trajectory &tj)
Point2_t dEdx
dE/dx for 3D matched trajectories
float Length(const PFPStruct &pfp)
bool dbgStitch
debug PFParticle stitching
Length_t WireCoordinate(Point_t const &pos, PlaneID const &planeid) const
Returns the index of the nearest wire to the specified position.
void SetEndPoints(Trajectory &tj)
bool TrajHitsOK(TCSlice const &slc, const std::vector< unsigned int > &iHitsInMultiplet, const std::vector< unsigned int > &jHitsInMultiplet)
bool MakeVertexObsolete(std::string fcnLabel, TCSlice &slc, VtxStore &vx2, bool forceKill)
std::vector< float > kinkCuts
kink finder algorithm
constexpr auto const & right(const_AssnsIter< L, R, D, Dir > const &a, const_AssnsIter< L, R, D, Dir > const &b)
bool InTrajOK(TCSlice &slc, std::string someText)
bool InsideFV(const TCSlice &slc, const PFPStruct &pfp, unsigned short end)
struct of temporary 2D vertices (end points)
std::vector< unsigned int > PutTrajHitsInVector(const Trajectory &tj, HitStatus_t hitRequest)
const std::vector< std::string > AlgBitNames
std::vector< ShowerStruct > cots
bool AttachAnyVertexToTraj(TCSlice &slc, int tjID, bool prt)
CTP_t CTP
Cryostat, TPC, Plane code.
std::vector< float > maxPos0
short recoTPC
only reconstruct in the seleted TPC
std::vector< int > NearTjIDs
std::array< double, 3 > Point3_t
std::vector< ShowerStruct3D > showers
bool SignalAtTp(TrajPoint &tp)
void SetPDGCode(TCSlice &slc, unsigned short itj)
std::vector< Point2_t > Envelope
tagged as a vertex between Tjs that are matched to MC truth neutrino interaction particles ...
unsigned int Nplanes() const
Number of planes in this tpc.
vertex position fixed manually - no fitting done
Declaration of signal hit object.
void FindAlongTrans(Point3_t pos1, Vector3_t dir1, Point3_t pos2, Point2_t &alongTrans)
void PrintTrajectory(std::string someText, const TCSlice &slc, const Trajectory &tj, unsigned short tPoint)
bool StartTraj(TCSlice const &slc, Trajectory &tj, unsigned int fromhit, unsigned int tohit, unsigned short pass)
void PrintTP(std::string someText, const TCSlice &slc, unsigned short ipt, short dir, unsigned short pass, const TrajPoint &tp)
std::vector< std::vector< std::pair< unsigned int, unsigned int > > > wireHitRange
The data type to uniquely identify a Plane.
void UnsetUsedHits(TCSlice &slc, TrajPoint &tp)
int ParentID
ID of the parent, or the ID of the Tj this one was merged with if it is killed.
std::vector< unsigned int > Hits
std::string PrintEndFlag(const PFPStruct &pfp, unsigned short end)
bool StoreTraj(TCSlice &slc, Trajectory &tj)
float TotChg
Total including an estimate for dead wires.
double Temperature() const
In kelvin.
short MCSMom
Normalized RMS using ALL hits. Assume it is 50% to start.
bool IsShowerLike(TCSlice const &slc, std::vector< int > const &TjIDs)
float TpSumHitChg(const TCSlice &slc, TrajPoint const &tp)
constexpr auto abs(T v)
Returns the absolute value of the argument.
void ChkChgAsymmetry(TCSlice &slc, Trajectory &tj, bool prt)
CryostatID_t Cryostat
Index of cryostat.
for(Int_t i=0;i< nentries;i++)
bool TrajClosestApproach(Trajectory const &tj, float x, float y, unsigned short &closePt, float &DOCA)
void FillWireHitRange(geo::TPCID inTPCID)
float TrajPointSeparation(const TrajPoint &tp1, const TrajPoint &tp2)
std::vector< std::pair< unsigned int, unsigned int > > tpcSrcHitRange
bool CompatibleMerge(const TCSlice &slc, std::vector< int > const &tjIDs, bool prt)
void SetAngleCode(TrajPoint &tp)
PFPStruct CreatePFP(const TCSlice &slc)
Point3_t PosAtEnd(const PFPStruct &pfp, unsigned short end)
float PointTrajDOCA(const TCSlice &slc, unsigned int iht, TrajPoint const &tp)
double DeltaAngle(const Vector3_t v1, const Vector3_t v2)
a general purpose flag bit used in 3D matching
bool dbgSlc
debug only in the user-defined slice? default is all slices
bool MakeBareTrajPoint(const TCSlice &slc, unsigned int fromHit, unsigned int toHit, TrajPoint &tp)
void Print2V(mf::LogVerbatim &myprt, VtxStore const &vx2, bool &printHeader)
std::vector< unsigned int > lastWire
the last wire with a hit
float GoodnessOfFit() const
Degrees of freedom in the determination of the hit signal shape (-1 by default)
bool LongPulseHit(const recob::Hit &hit)
std::array< int, 2 > Vx3ID
bool expectSlicedHits
info passed from the module - used to (not) define wireHitRange
void PrintPFPs(std::string someText, TCSlice const &slc)
short int Multiplicity() const
How many hits could this one be shared with.
float TPHitsRMSTime(const TCSlice &slc, const TrajPoint &tp, HitStatus_t hitRequest)
unsigned int MVI
MatchVec Index for detailed 3D matching.
std::string TPEnvString(const TrajPoint &tp)
pure virtual base interface for detector clocks
std::vector< int > CotIDs
std::string PrintPos(const TrajPoint &tp)
bool StoreVertex(TCSlice &slc, VtxStore &vx)
float MaxTjLen(const TCSlice &slc, std::vector< int > &tjIDs)
bool MergeTjIntoPFP(TCSlice &slc, int mtjid, PFPStruct &pfp, bool prt)
void Print3S(detinfo::DetectorPropertiesData const &detProp, mf::LogVerbatim &myprt, ShowerStruct3D const &ss3)
double DeltaAngle2(double Ang1, double Ang2)
std::vector< T > SetIntersection(const std::vector< T > &set1, const std::vector< T > &set2)
std::vector< float > angleRanges
list of max angles for each angle range
const std::vector< std::string > EndFlagNames
std::vector< float > showerTag
shower-like trajectory tagging + shower reconstruction
unsigned short Pass
the pass on which it was created
bool dbg3V
debug 3D vertex finding
double Efield(unsigned int planegap=0) const
kV/cm
bool SignalBetween(const TrajPoint &tp1, const TrajPoint &tp2, const float MinWireSignalFraction)
float MaxChargeAsymmetry(TCSlice &slc, std::vector< int > const &tjIDs)
TPCGeo const & TPC(TPCID const &tpcid=tpc_zero) const
Returns the specified TPC.
float OverlapFraction(const Trajectory &tj1, const Trajectory &tj2)
std::string PrintHitShort(const TCHit &tch)
decltype(auto) constexpr end(T &&obj)
ADL-aware version of std::end.
decltype(auto) constexpr size(T &&obj)
ADL-aware version of std::size.
float HitsRMSTime(const TCSlice &slc, const std::vector< unsigned int > &hitsInMultiplet, HitStatus_t hitRequest)
Access the description of detector geometry.
unsigned int Hit
set to the hit index in evt.allHits if a Plane:Wire:Tick match is found
std::vector< int > TjUIDs
void FitPar(const Trajectory &tj, unsigned short originPt, unsigned short npts, short fitDir, ParFit &pFit, unsigned short usePar)
std::vector< unsigned int > PutHitsInVector(const TCSlice &slc, PFPStruct const &pfp, HitStatus_t hitRequest)
const std::vector< std::string > StrategyBitNames
float HitSep2(const TCSlice &slc, unsigned int iht, unsigned int jht)
int Cryostat
Select Cryostat.
int NeutrinoPrimaryTjID(const TCSlice &slc, const Trajectory &tj)
float HitsRMSTick(const TCSlice &slc, const std::vector< unsigned int > &hitsInMultiplet, HitStatus_t hitRequest)
struct of temporary 3D vertices
short nPtsAve
dump trajectory points
int PrimaryID(const TCSlice &slc, const Trajectory &tj)
unsigned short NearestPtWithChg(const Trajectory &tj, unsigned short thePt)
bool dbgStp
debug stepping using debug.Cryostat, debug.TPC, etc
int Wire
Select hit Wire for debugging.
std::array< float, 2 > Point2_t
std::vector< float > maxPos1
int PDGCodeVote(detinfo::DetectorClocksData const &clockData, detinfo::DetectorPropertiesData const &detProp, const TCSlice &slc, PFPStruct &pfp)
float unitsPerTick
scale factor from Tick to WSE equivalent units
IDparameter< geo::WireID > WireID
Member type of validated geo::WireID parameter.
int WorkID
Select the StartWorkID for debugging.
float DeadWireCount(const TCSlice &slc, const TrajPoint &tp1, const TrajPoint &tp2)
bool BraggSplit(TCSlice &slc, unsigned short itj)
void DefineHitPos(TCSlice &slc, TrajPoint &tp)
Vector3_t PointDirection(const Point3_t p1, const Point3_t p2)
float PointTrajDOCA2(float wire, float time, TrajPoint const &tp)
TP is near a hit in the srcHit collection but no allHit hit exists (DUNE disambiguation error) ...
CTP_t CTP
Cryostat, TPC, Plane code.
decltype(auto) constexpr to_string(T &&obj)
ADL-aware version of std::to_string.
bool MergeShowerTjsAndStore(TCSlice &slc, unsigned short istj, unsigned short jstj, bool prt)
bool dbg2V
debug 2D vertex finding
std::vector< float > aveHitRMS
average RMS of an isolated hit
std::vector< TrajPoint > Pts
Trajectory points.
float ElectronLikelihood(const TCSlice &slc, const Trajectory &tj)
float TwoTPAngle(const TrajPoint &tp1, const TrajPoint &tp2)
void TjDeltaRMS(const Trajectory &tj, unsigned short firstPt, unsigned short lastPt, double &rms, unsigned short &cnt)
std::vector< std::vector< bool > > goodWire
float ChgFracNearPos(const TCSlice &slc, const Point2_t &pos, const std::vector< int > &tjIDs)
std::vector< unsigned int > FindCloseHits(const TCSlice &slc, std::array< int, 2 > const &wireWindow, Point2_t const &timeWindow, const unsigned short plane, HitStatus_t hitRequest, bool usePeakTime, bool &hitsNear)
float ChgFracBetween(detinfo::DetectorPropertiesData const &detProp, const TCSlice &slc, Point3_t pos1, Point3_t pos2)
void MakeHaloTj(TCSlice &slc, Trajectory &muTj, bool prt)
TrajPoint MakeBareTP(detinfo::DetectorPropertiesData const &detProp, const Point3_t &pos, CTP_t inCTP)
void PrintT(mf::LogVerbatim &myprt, Trajectory const &tj, bool &printHeader)
std::vector< VtxStore > vtxs
2D vertices
std::array< unsigned short, 2 > EndPt
First and last point in the trajectory that has charge.
unsigned short PDGCode
shower-like or track-like {default is track-like}
double ConvertXToTicks(double X, int p, int t, int c) const
bool PointInsideEnvelope(const Point2_t &Point, const std::vector< Point2_t > &Envelope)
std::vector< float > match3DCuts
3D matching cuts
std::vector< SectionFit > SectionFits
unsigned int NumberTimeSamples() const
void MakeTrajectoryObsolete(TCSlice &slc, unsigned int itj)
short StartEnd
The starting end (-1 = don't know)
double PosSep2(const Point3_t &pos1, const Point3_t &pos2)
unsigned short GetPFPIndex(const TCSlice &slc, int tjID)
std::string PrintHit(const TCHit &tch)
bool Fit2D(short mode, Point2_t inPt, float &inPtErr, Vector2_t &outVec, Vector2_t &outVecErr, float &chiDOF)
bool SplitTraj(detinfo::DetectorPropertiesData const &detProp, TCSlice &slc, unsigned short itj, float XPos, bool makeVx2, bool prt)
const geo::GeometryCore * geom
double DriftVelocity(double efield=0., double temperature=0.) const
cm/us
int UID
a unique ID for all slices
The data type to uniquely identify a TPC.
PlaneID_t Plane
Index of the plane within its TPC.
void PrintP(mf::LogVerbatim &myprt, PFPStruct const &pfp, bool &printHeader)
float ExpectedHitsRMS(const TrajPoint &tp)
bool valsIncreasing(const SortEntry &c1, const SortEntry &c2)
void ReleaseHits(TCSlice &slc, Trajectory const &tj)
bool SetMag(Vector3_t &v1, double mag)
std::array< double, 2 > Vector2_t
Definition of data types for geometry description.
unsigned short NumHitsInTP(const TrajPoint &tp, HitStatus_t hitRequest)
void TrimEndPts(std::string fcnLabel, TCSlice &slc, Trajectory &tj, const std::vector< float > &fQualityCuts, bool prt)
unsigned short FarEnd(const PFPStruct &pfp, const Point3_t &pos)
void DefineTjParents(TCSlice &slc, bool prt)
float HitsPosTick(const TCSlice &slc, const std::vector< unsigned int > &hitsInMultiplet, float &sum, HitStatus_t hitRequest)
std::vector< unsigned int > firstWire
the first wire with a hit
int ID
ID that is local to one slice.
std::vector< TCSlice > slices
std::array< unsigned short, 2 > VtxID
ID of 2D vertex.
Detector simulation of raw signals on wires.
float MCSThetaRMS(const Trajectory &tj)
std::bitset< 16 > modes
number of points to find AveChg
std::vector< TCHit > slHits
unsigned short NearbyCleanPt(const Trajectory &tj, unsigned short end)
bool MergeAndStore(TCSlice &slc, unsigned int itj1, unsigned int itj2, bool doPrt)
void RestoreObsoleteTrajectory(TCSlice &slc, unsigned int itj)
bool aveHitRMSValid
set true when the average hit RMS is well-known
std::bitset< 16 > Stat
Vertex status bits using kVtxBit_t.
void PrintPFP(std::string someText, TCSlice const &slc, const PFPStruct &pfp, bool printHeader)
ROOT::Math::PositionVector3D< ROOT::Math::Cartesian3D< double >, ROOT::Math::GlobalCoordinateSystemTag > Point_t
Type for representation of position in physical 3D space.
std::vector< int > GetAssns(TCSlice const &slc, std::string type1Name, int id, std::string type2Name)
int Tick
Select hit PeakTime for debugging (< 0 for vertex finding)
void TrajIntersection(TrajPoint const &tp1, TrajPoint const &tp2, Point2_t &pos)
raw::ChannelID_t PlaneWireToChannel(WireID const &wireid) const
Returns the ID of the TPC channel connected to the specified wire.
const std::vector< std::string > VtxBitNames
int PrimaryUID(const PFPStruct &pfp)
double DotProd(const Vector3_t &v1, const Vector3_t &v2)
void UpdateTjChgProperties(std::string inFcnLabel, TCSlice const &slc, Trajectory &tj, bool prt)
std::vector< Vtx3Store > vtx3s
3D vertices
float KinkSignificance(TCSlice const &slc, Trajectory const &tj1, unsigned short end1, Trajectory const &tj2, unsigned short end2, unsigned short nPtsFit, bool useChg, bool prt)
void Print3V(detinfo::DetectorPropertiesData const &detProp, mf::LogVerbatim &myprt, Vtx3Store const &vx3, bool &printHeader)
bool DecodeDebugString(std::string strng)
std::vector< recob::Hit > const * srcHits
std::bitset< 128 > useAlg
Max hit separation for making junk trajectories. < 0 to turn off.
Contains all timing reference information for the detector.
short StepDir
-1 = going US (-> small wire#), 1 = going DS (-> large wire#)
float MaxHitDelta(TCSlice &slc, Trajectory &tj)
float TrajLength(const Trajectory &tj)
std::bitset< 128 > AlgMod
Bit set if algorithm AlgBit_t modifed the trajectory.
std::vector< short > muonTag
void UpdateVxEnvironment(TCSlice &slc)
std::tuple< double, double, const reco::ClusterHit3D * > Point
Definitions used by the VoronoiDiagram algorithm.
geo::PlaneID DecodeCTP(CTP_t CTP)
void TagJunkTj(Trajectory &tj, bool prt)
View_t View(PlaneID const &pid) const
Returns the view (wire orientation) on the channels of specified TPC plane.
std::vector< int > GetVtxTjIDs(const TCSlice &slc, const VtxStore &vx2)
double PosSep(const Point3_t &pos1, const Point3_t &pos2)
MaybeLogger_< ELseverityLevel::ELsev_warning, false > LogWarning
void TrimHiChgEndPts(TCSlice &slc, Trajectory &tj, bool prt)
std::bitset< 128 > dbgAlg
Allow user to turn on debug printing in algorithms (that print...)
unsigned short PDGCodeIndex(int PDGCode)
float PointPull(const TP3D &tp3d)
std::bitset< 8 > Environment
std::vector< recob::Hit > const * allHits
bool valsDecreasing(const SortEntry &c1, const SortEntry &c2)
CTP_t EncodeCTP(unsigned int cryo, unsigned int tpc, unsigned int plane)
unsigned int Nwires(PlaneID const &planeid) const
Returns the total number of wires in the specified plane.
std::pair< unsigned short, unsigned short > GetSliceIndex(std::string typeName, int uID)
unsigned short MVI_Iter
MVI iteration - see FindPFParticles.
span(IterB &&b, IterE &&e, Adaptor &&adaptor) -> span< decltype(adaptor(std::forward< IterB >(b))), decltype(adaptor(std::forward< IterE >(e))) >
std::vector< unsigned int > nWires
use the stiff electron strategy
std::array< std::bitset< 8 >, 2 > EndFlag
float PointTrajSep2(float wire, float time, TrajPoint const &tp)
std::array< double, 3 > Vector3_t
std::vector< float > chkStopCuts
Bragg peak finder configuration.
bool SignalAtTpInSlc(const TCSlice &slc, const TrajPoint &tp)
std::vector< TP3D > TP3Ds
int ID
ID of the recob::Slice (not the sub-slice)
std::array< std::vector< float >, 2 > dEdx
void SetVx2Score(TCSlice &slc)
unsigned int Nplanes(TPCID const &tpcid=tpc_zero) const
Returns the total number of planes in the specified TPC.
std::vector< float > vtx2DCuts
Max position pull, max Position error rms.
std::bitset< 8 > Strategy
bool NearbySrcHit(geo::PlaneID plnID, unsigned int wire, float loTick, float hiTick)
std::vector< PFPStruct > pfps
std::vector< int > FindCloseTjs(const TCSlice &slc, const TrajPoint &fromTp, const TrajPoint &toTp, const float &maxDelta)
2D representation of charge deposited in the TDC/wire plane
void PrintAllTraj(detinfo::DetectorPropertiesData const &detProp, std::string someText, TCSlice const &slc, unsigned short itj, unsigned short ipt, bool prtVtx)
void MoveTPToWire(TrajPoint &tp, float wire)
unsigned short NumPtsWithCharge(const TCSlice &slc, const Trajectory &tj, bool includeDeadWires)
unsigned short CloseEnd(const Trajectory &tj, const Point2_t &pos)
void MergeGhostTjs(TCSlice &slc, CTP_t inCTP)
bool TrajIsClean(Trajectory const &tj, bool prt)
unsigned int ChannelID_t
Type representing the ID of a readout channel.
unsigned int allHitsIndex
TPCID_t TPC
Index of the TPC within its cryostat.
bool StorePFP(TCSlice &slc, PFPStruct &pfp)
double sampling_rate(DetectorClocksData const &data)
Returns the period of the TPC readout electronics clock.
bool TrajTrajDOCA(const TCSlice &slc, const Trajectory &tj1, const Trajectory &tj2, unsigned short &ipt1, unsigned short &ipt2, float &minSep)
std::vector< std::vector< std::pair< unsigned int, unsigned int > > > wireHitRange
Vector3_t DirAtEnd(const PFPStruct &pfp, unsigned short end)
short recoSlice
only reconstruct the slice with ID (0 = all)
float TjChgFrac
Fraction of charge near the vertex that is from hits on the vertex Tjs.
bool dbgSummary
print a summary report
bool NeedsUpdate
Set true when the Tj needs to be updated.
bool WireHitRangeOK(TCSlice const &slc, const CTP_t &inCTP)
void TrajPointTrajDOCA(TrajPoint const &tp, Trajectory const &tj, unsigned short &closePt, float &minSep)
master switch for turning on debug mode
Length_t WirePitch(PlaneID const &planeid=plane_zero) const
Returns the distance between two consecutive wires.
float TPHitsRMSTick(const TCSlice &slc, const TrajPoint &tp, HitStatus_t hitRequest)
void PrintTPHeader(std::string someText)
void ChkEndKink(TCSlice const &slc, Trajectory &tj, bool prt)
CTP_t CTP
set to an invalid CTP
std::vector< int > DtrUIDs
constexpr Point origin()
Returns a origin position with a point of the specified type.
unsigned short AngleRange(TrajPoint const &tp)
use the stiff muon strategy
void PosInPlane(detinfo::DetectorPropertiesData const &detProp, const Vtx3Store &vx3, unsigned short plane, Point2_t &pos)
Encapsulate the construction of a single detector plane.
std::array< std::bitset< 8 >, 2 > EndFlag
void SetTPEnvironment(TCSlice &slc, CTP_t inCTP)
bool HasDuplicateHits(const TCSlice &slc, Trajectory const &tj, bool prt)
void ReverseTraj(Trajectory &tj)