35 #include "larevt/CalibrationDBI/Interface/ChannelStatusProvider.h" 36 #include "larevt/CalibrationDBI/Interface/ChannelStatusService.h" 47 return c1.length > c2.length;
56 fNumPass = pset.
get<
unsigned short>(
"NumPass", 0);
57 fMaxHitsFit = pset.
get<std::vector<unsigned short>>(
"MaxHitsFit");
58 fMinHits = pset.
get<std::vector<unsigned short>>(
"MinHits");
59 fNHitsAve = pset.
get<std::vector<unsigned short>>(
"NHitsAve");
60 fChgCut = pset.
get<std::vector<float>>(
"ChgCut");
61 fChiCut = pset.
get<std::vector<float>>(
"ChiCut");
62 fMaxWirSkip = pset.
get<std::vector<unsigned short>>(
"MaxWirSkip");
63 fMinWirAfterSkip = pset.
get<std::vector<unsigned short>>(
"MinWirAfterSkip");
64 fKinkChiRat = pset.
get<std::vector<float>>(
"KinkChiRat");
65 fKinkAngCut = pset.
get<std::vector<float>>(
"KinkAngCut");
66 fDoMerge = pset.
get<std::vector<bool>>(
"DoMerge");
67 fTimeDelta = pset.
get<std::vector<float>>(
"TimeDelta");
68 fMergeChgCut = pset.
get<std::vector<float>>(
"MergeChgCut");
69 fFindVertices = pset.
get<std::vector<bool>>(
"FindVertices");
70 fLACrawl = pset.
get<std::vector<bool>>(
"LACrawl");
71 fMinAmp = pset.
get<std::vector<float>>(
"MinAmp", {5, 5, 5});
72 fChgNearWindow = pset.
get<
float>(
"ChgNearWindow");
73 fChgNearCut = pset.
get<
float>(
"ChgNearCut");
75 fChkClusterDS = pset.
get<
bool>(
"ChkClusterDS",
false);
76 fVtxClusterSplit = pset.
get<
bool>(
"VtxClusterSplit",
false);
77 fFindStarVertices = pset.
get<
bool>(
"FindStarVertices",
false);
78 if (pset.
has_key(
"HammerCluster")) {
80 <<
"fcl setting HammerCluster is replaced by FindHammerClusters. Ignoring...";
82 fFindHammerClusters = pset.
get<
bool>(
"FindHammerClusters",
false);
83 fKillGarbageClusters = pset.
get<
float>(
"KillGarbageClusters", 0);
84 fRefineVertexClusters = pset.
get<
bool>(
"RefineVertexClusters",
false);
85 fHitErrFac = pset.
get<
float>(
"HitErrFac", 0.2);
86 fHitMinAmp = pset.
get<
float>(
"HitMinAmp", 0.2);
87 fClProjErrFac = pset.
get<
float>(
"ClProjErrFac", 4);
88 fMinHitFrac = pset.
get<
float>(
"MinHitFrac", 0.6);
90 fLAClusAngleCut = pset.
get<
float>(
"LAClusAngleCut", 45);
91 fLAClusMaxHitsFit = pset.
get<
unsigned short>(
"LAClusMaxHitsFit");
92 fMergeAllHits = pset.
get<
bool>(
"MergeAllHits",
false);
93 fHitMergeChiCut = pset.
get<
float>(
"HitMergeChiCut", 2.5);
94 fMergeOverlapAngCut = pset.
get<
float>(
"MergeOverlapAngCut");
95 fAllowNoHitWire = pset.
get<
unsigned short>(
"AllowNoHitWire", 0);
96 fVertex2DCut = pset.
get<
float>(
"Vertex2DCut", 5);
97 fVertex2DWireErrCut = pset.
get<
float>(
"Vertex2DWireErrCut", 5);
98 fVertex3DCut = pset.
get<
float>(
"Vertex3DCut", 5);
100 fDebugPlane = pset.
get<
int>(
"DebugPlane", -1);
101 fDebugWire = pset.
get<
int>(
"DebugWire", -1);
102 fDebugHit = pset.
get<
int>(
"DebugHit", -1);
105 bool badinput =
false;
106 if (fNumPass > fMaxHitsFit.size()) badinput =
true;
107 if (fNumPass > fMinHits.size()) badinput =
true;
108 if (fNumPass > fNHitsAve.size()) badinput =
true;
109 if (fNumPass > fChgCut.size()) badinput =
true;
110 if (fNumPass > fChiCut.size()) badinput =
true;
111 if (fNumPass > fMaxWirSkip.size()) badinput =
true;
112 if (fNumPass > fMinWirAfterSkip.size()) badinput =
true;
113 if (fNumPass > fKinkChiRat.size()) badinput =
true;
114 if (fNumPass > fKinkAngCut.size()) badinput =
true;
115 if (fNumPass > fDoMerge.size()) badinput =
true;
116 if (fNumPass > fTimeDelta.size()) badinput =
true;
117 if (fNumPass > fMergeChgCut.size()) badinput =
true;
118 if (fNumPass > fFindVertices.size()) badinput =
true;
119 if (fNumPass > fLACrawl.size()) badinput =
true;
123 <<
"ClusterCrawlerAlg: Bad input from fcl file";
137 if (cmp_res != 0)
return cmp_res < 0;
145 void ClusterCrawlerAlg::ClearResults()
155 void ClusterCrawlerAlg::CrawlInit()
181 WireHitRange.clear();
187 void ClusterCrawlerAlg::ClusterInit()
203 std::vector<recob::Hit>
const& srchits)
211 if (fHits.size() < 3)
return;
212 if (fHits.size() > UINT_MAX) {
213 mf::LogWarning(
"CC") <<
"Too many hits for ClusterCrawler " << fHits.size();
218 if (fNumPass == 0)
return;
223 std::sort(fHits.begin(), fHits.end(), &SortByMultiplet);
225 inClus.resize(fHits.size());
226 mergeAvailable.resize(fHits.size());
227 for (
unsigned int iht = 0; iht < inClus.size(); ++iht) {
229 mergeAvailable[iht] =
false;
234 cstat = tpcid.Cryostat;
237 plane = planeid.Plane;
238 WireHitRange.clear();
240 clCTP =
EncodeCTP(tpcid.Cryostat, tpcid.TPC, planeid.Plane);
241 cstat = tpcid.Cryostat;
248 if (WireHitRange.empty() || (fFirstWire == fLastWire))
continue;
252 float wirePitch = geom->WirePitch(geom->View(channel));
255 fScaleF = tickToDist / wirePitch;
257 if (fLAClusAngleCut > 0)
258 fLAClusSlopeCut = std::tan(3.142 * fLAClusAngleCut / 180.) / fScaleF;
260 fNumWires = geom->Nwires(planeid);
262 if (fNumPass > 0) ClusterLoop();
264 if (fVertex3DCut > 0) {
266 VtxMatch(det_prop, tpcid);
267 Vtx3ClusterMatch(det_prop, tpcid);
268 if (fFindHammerClusters) FindHammerClusters(det_prop);
270 Vtx3ClusterSplit(det_prop, tpcid);
272 if (fDebugPlane >= 0) {
279 WireHitRange.clear();
286 RemoveObsoleteHits();
291 void ClusterCrawlerAlg::ClusterLoop()
295 unsigned int nHitsUsed = 0, iwire, jwire, kwire;
296 bool AllDone =
false, SigOK =
false, HitOK =
false;
297 unsigned int ihit, jhit;
298 for (
unsigned short thispass = 0; thispass < fNumPass; ++thispass) {
301 unsigned int span = 3;
302 if (fMinHits[pass] < span) span = fMinHits[pass];
303 for (iwire = fLastWire; iwire > fFirstWire +
span; --iwire) {
305 if (WireHitRange[iwire].first < 0)
continue;
306 auto ifirsthit = (
unsigned int)WireHitRange[iwire].first;
307 auto ilasthit = (
unsigned int)WireHitRange[iwire].
second;
308 for (ihit = ifirsthit; ihit < ilasthit; ++ihit) {
309 bool ClusterAdded =
false;
312 if (ihit > fHits.size() - 1) {
313 mf::LogError(
"CC") <<
"ClusterLoop bad ihit " << ihit <<
" fHits size " << fHits.size();
317 if (inClus[ihit] != 0)
continue;
319 if (fHits[ihit].PeakAmplitude() < fHitMinAmp)
continue;
320 if ((iwire + 1) < span)
continue;
321 jwire = iwire - span + 1;
325 if (WireHitRange[jwire].first == -2)
continue;
326 if (WireHitRange[jwire].first == -1) {
328 unsigned int nmissed = 0;
329 while (WireHitRange[jwire].first == -1 && jwire > 1 && nmissed < fMaxWirSkip[pass]) {
335 <<
" new jwire " << jwire <<
" dead? " << WireHitRange[jwire].first;
336 if (WireHitRange[jwire].first < 0)
continue;
340 unsigned int useHit = 0;
341 bool doConstrain =
false;
342 VtxConstraint(iwire, ihit, jwire, useHit, doConstrain);
343 unsigned int jfirsthit = (
unsigned int)WireHitRange[jwire].first;
344 unsigned int jlasthit = (
unsigned int)WireHitRange[jwire].
second;
345 if (jfirsthit > fHits.size() - 1 || jfirsthit > fHits.size() - 1)
347 <<
"ClusterLoop jwire " << jwire <<
" bad firsthit " << jfirsthit <<
" lasthit " 348 << jlasthit <<
" fhits size " << fHits.size();
349 for (jhit = jfirsthit; jhit < jlasthit; ++jhit) {
350 if (jhit > fHits.size() - 1)
352 <<
"ClusterLoop bad jhit " << jhit <<
" firsthit " << jfirsthit <<
" lasthit " 353 << jlasthit <<
" fhits size" << fHits.size();
355 if (doConstrain && jhit != useHit)
continue;
358 if (inClus[jhit] != 0)
continue;
360 if (fHits[jhit].PeakAmplitude() < fHitMinAmp)
continue;
363 fcl2hits.push_back(ihit);
364 chifits.push_back(0.);
365 hitNear.push_back(0);
366 chgNear.push_back(0);
368 fcl2hits.push_back(jhit);
369 chifits.push_back(0.);
370 hitNear.push_back(0);
371 chgNear.push_back(0);
376 clparerr[1] = 0.2 *
std::abs(clpar[1]);
377 clpar[2] = fHits[jhit].WireID().Wire;
381 for (kwire = jwire + 1; kwire < iwire; ++kwire) {
383 if (CrawlVtxChk(kwire)) {
387 AddHit(kwire, HitOK, SigOK);
389 mf::LogVerbatim(
"CC") <<
" HitOK " << HitOK <<
" clChisq " << clChisq <<
" cut " 390 << fChiCut[pass] <<
" ClusterHitsOK " << ClusterHitsOK(-1);
394 if (clChisq > 2)
break;
396 if (!ClusterHitsOK(-1))
continue;
401 prt = (fDebugPlane == (int)plane && (
int)iwire == fDebugWire &&
405 <<
"ADD >>>>> Starting cluster with hits " <<
PrintHit(fcl2hits[0]) <<
" " 406 <<
PrintHit(fcl2hits[1]) <<
" " <<
PrintHit(fcl2hits[2]) <<
" on pass " << pass;
410 clBeginSlp = clpar[1];
413 if (!fLACrawl[pass] &&
std::abs(clBeginSlp) > fLAClusSlopeCut)
continue;
416 if (CrawlVtxChk2())
continue;
417 clBeginSlpErr = clparerr[1];
421 for (
unsigned short kk = 0; kk < fcl2hits.size(); ++kk) {
422 fAveHitWidth += fHits[fcl2hits[kk]].EndTick() - fHits[fcl2hits[kk]].StartTick();
424 fAveHitWidth /= (float)fcl2hits.size();
430 if (fLACrawl[pass] && fLAClusSlopeCut > 0) {
432 if (
std::abs(clBeginSlp) > fLAClusSlopeCut) { LACrawlUS(); }
441 if (fcl2hits.size() >= fMinHits[pass]) {
444 clEndSlpErr = clparerr[1];
447 mf::LogError(
"CC") <<
"Failed to store cluster in plane " << plane;
451 nHitsUsed += fcl2hits.size();
452 AllDone = (nHitsUsed == fHits.size());
459 if (ClusterAdded || AllDone)
break;
467 if (fDoMerge[pass]) ChkMerge();
469 if (fFindVertices[pass]) FindVertices();
476 if (fKillGarbageClusters > 0 && !tcl.empty()) KillGarbageClusters();
480 if (fChkClusterDS) ChkClusterDS();
482 if (fVtxClusterSplit) {
483 bool didSomething = VtxClusterSplit();
485 if (didSomething) VtxClusterSplit();
488 if (fFindStarVertices) FindStarVertices();
490 if (fDebugPlane == (
int)plane) {
495 CheckHitClusterAssociations();
500 void ClusterCrawlerAlg::KillGarbageClusters()
504 if (tcl.size() < 2)
return;
506 unsigned short icl, jcl;
509 std::vector<float> iHits, jHits;
513 float sepcut = 0, iFrac, jFrac;
515 unsigned short iLoIndx, jLoIndx, olapSize, iop, ii, jj;
516 unsigned short nclose;
519 std::vector<unsigned short> killMe;
521 for (icl = 0; icl < tcl.size() - 1; ++icl) {
522 if (tcl[icl].ID < 0)
continue;
523 if (tcl[icl].CTP != clCTP)
continue;
527 iHits.resize(tcl[icl].BeginWir - tcl[icl].EndWir + 1, 1000);
529 for (
auto iht : tcl[icl].tclhits) {
530 indx = fHits[iht].WireID().Wire - tcl[icl].EndWir;
531 if (indx > iHits.size() - 1) {
532 mf::LogWarning(
"CC") <<
"KillGarbageClusters: icl ID " << tcl[icl].ID <<
" Bad indx " 533 << indx <<
" " << iHits.size() <<
"\n";
536 iHits[indx] = fHits[iht].PeakTime();
537 iChg += fHits[iht].Integral();
538 if (first) sepcut += fHits[iht].RMS();
541 sepcut /= (float)tcl[icl].tclhits.size();
547 for (jcl = icl + 1; jcl < tcl.size(); ++jcl) {
548 if (tcl[jcl].ID < 0)
continue;
549 if (tcl[jcl].CTP != clCTP)
continue;
551 if (tcl[icl].BeginWir < tcl[jcl].EndWir)
continue;
552 if (tcl[icl].EndWir > tcl[jcl].BeginWir)
continue;
554 if (
std::abs(tcl[icl].BeginAng - tcl[jcl].BeginAng) > fKillGarbageClusters)
continue;
556 if (tcl[icl].EndWir < tcl[jcl].EndWir) {
560 iLoIndx = tcl[jcl].EndWir - tcl[icl].EndWir;
562 if (tcl[icl].BeginWir < tcl[jcl].BeginWir) {
566 olapSize = tcl[icl].BeginWir - tcl[jcl].EndWir + 1;
572 olapSize = tcl[jcl].BeginWir - tcl[jcl].EndWir + 1;
580 jLoIndx = tcl[icl].EndWir - tcl[icl].EndWir;
581 if (tcl[icl].BeginWir < tcl[jcl].BeginWir) {
585 olapSize = tcl[icl].BeginWir - tcl[icl].EndWir + 1;
591 olapSize = tcl[jcl].BeginWir - tcl[icl].EndWir + 1;
596 jHits.resize(tcl[jcl].BeginWir - tcl[jcl].EndWir + 1, -1000);
598 for (
auto jht : tcl[jcl].tclhits) {
599 indx = fHits[jht].WireID().Wire - tcl[jcl].EndWir;
600 if (indx > jHits.size() - 1) {
601 mf::LogWarning(
"CC") <<
"KillGarbageClusters: jcl ID " << tcl[jcl].ID <<
" Bad indx " 602 << indx <<
" " << jHits.size() <<
"\n";
605 jHits[indx] = fHits[jht].PeakTime();
606 jChg += fHits[jht].Integral();
610 for (iop = 0; iop < olapSize; ++iop) {
612 if (ii > iHits.size() - 1)
continue;
614 if (jj > jHits.size() - 1)
continue;
615 if (
std::abs(iHits[ii] - jHits[jj]) < sepcut) ++nclose;
617 iFrac = (float)nclose / (
float)iHits.size();
618 jFrac = (float)nclose / (
float)jHits.size();
619 if (iFrac < 0.5 && jFrac < 0.5)
continue;
620 doKill = (iFrac < jFrac && iChg > jChg);
621 if (doKill) killMe.push_back(jcl);
625 if (killMe.size() == 0)
return;
626 for (
auto icl : killMe) {
628 if (tcl[icl].ID < 0)
continue;
629 tcl[icl].ProcCode = 666;
630 MakeClusterObsolete(icl);
649 unsigned short icl, jcl;
651 bool chkprt = (fDebugWire == 666);
652 if (chkprt)
mf::LogVerbatim(
"CC") <<
"Inside MergeOverlap using clCTP " << clCTP;
654 unsigned short minLen = 6;
655 unsigned short minOvrLap = 2;
658 float maxDTick = fTimeDelta[fTimeDelta.size() - 1];
660 unsigned int overlapSize, ii, indx, bWire, eWire;
661 unsigned int iht, jht;
662 float dang, prtime, dTick;
663 for (icl = 0; icl < tcl.size(); ++icl) {
664 if (tcl[icl].ID < 0)
continue;
665 if (tcl[icl].CTP != clCTP)
continue;
666 prt = chkprt && fDebugPlane == (int)clCTP;
667 if (tcl[icl].BeginVtx >= 0)
continue;
668 if (tcl[icl].tclhits.size() < minLen)
continue;
669 for (jcl = 0; jcl < tcl.size(); ++jcl) {
670 if (icl == jcl)
continue;
671 if (tcl[jcl].ID < 0)
continue;
672 if (tcl[jcl].CTP != clCTP)
continue;
673 if (tcl[jcl].EndVtx >= 0)
continue;
674 if (tcl[jcl].tclhits.size() < minLen)
continue;
676 if (tcl[icl].BeginWir < tcl[jcl].EndWir + minOvrLap)
continue;
678 if (tcl[icl].BeginWir > tcl[jcl].BeginWir - minOvrLap)
continue;
680 if (tcl[jcl].EndWir < tcl[icl].EndWir + minOvrLap)
continue;
681 dang =
std::abs(tcl[icl].BeginAng - tcl[jcl].EndAng);
683 mf::LogVerbatim(
"CC") <<
"MergeOverlap icl ID " << tcl[icl].ID <<
" jcl ID " 684 << tcl[jcl].ID <<
" dang " << dang;
685 if (dang > 0.5)
continue;
686 overlapSize = tcl[icl].BeginWir - tcl[jcl].EndWir + 1;
687 eWire = tcl[jcl].EndWir;
688 bWire = tcl[icl].BeginWir;
690 mf::LogVerbatim(
"CC") <<
" Candidate icl ID " << tcl[icl].ID <<
" " << tcl[icl].EndWir
691 <<
"-" << tcl[icl].BeginWir <<
" jcl ID " << tcl[jcl].ID <<
" " 692 << tcl[jcl].EndWir <<
"-" << tcl[jcl].BeginWir <<
" overlapSize " 693 << overlapSize <<
" bWire " << bWire <<
" eWire " << eWire;
696 for (ii = 0; ii < tcl[icl].tclhits.size(); ++ii) {
697 iht = tcl[icl].tclhits[ii];
698 if (fHits[iht].
WireID().Wire < eWire)
break;
701 dTick =
std::abs(fHits[iht].PeakTime() - tcl[jcl].EndTim);
702 if (dTick > maxDTick)
continue;
705 << tcl[jcl].EndTim <<
" dTick " << dTick;
706 for (ii = 0; ii < tcl[jcl].tclhits.size(); ++ii) {
707 jht = tcl[jcl].tclhits[tcl[jcl].tclhits.size() - ii - 1];
708 if (fHits[jht].
WireID().Wire > bWire)
break;
710 dTick =
std::abs(fHits[jht].PeakTime() - tcl[icl].BeginTim);
711 if (dTick > maxDTick)
continue;
714 << tcl[icl].BeginTim <<
" dTick " << dTick;
716 clpar[0] = fHits[iht].PeakTime();
717 clpar[2] = fHits[iht].WireID().Wire;
718 clpar[1] = (fHits[jht].PeakTime() - fHits[iht].PeakTime()) /
719 ((
float)fHits[jht].WireID().Wire - clpar[2]);
721 std::vector<unsigned int> oWireHits(overlapSize, INT_MAX);
722 std::vector<float> delta(overlapSize, maxDTick);
723 for (ii = 0; ii < tcl[icl].tclhits.size(); ++ii) {
724 iht = tcl[icl].tclhits[ii];
725 if (fHits[iht].
WireID().Wire < eWire)
break;
726 prtime = clpar[0] + clpar[1] * ((float)fHits[iht].
WireID().Wire - clpar[2]);
727 dTick =
std::abs(fHits[iht].PeakTime() - prtime);
728 indx = fHits[iht].WireID().Wire - eWire;
729 if (dTick > delta[indx])
continue;
731 oWireHits[indx] = iht;
734 for (ii = 0; ii < tcl[jcl].tclhits.size(); ++ii) {
735 jht = tcl[jcl].tclhits[tcl[jcl].tclhits.size() - ii - 1];
736 if (fHits[jht].
WireID().Wire > bWire)
break;
737 prtime = clpar[0] + clpar[1] * ((float)fHits[jht].
WireID().Wire - clpar[2]);
738 dTick =
std::abs(fHits[jht].PeakTime() - prtime);
739 indx = fHits[jht].WireID().Wire - eWire;
740 if (dTick > delta[indx])
continue;
742 oWireHits[indx] = jht;
746 for (ii = 0; ii < oWireHits.size(); ++ii) {
747 if (oWireHits[ii] == INT_MAX)
continue;
749 fcl2hits.push_back(iht);
752 if (fcl2hits.size() < 0.5 * overlapSize)
continue;
753 if (fcl2hits.size() < 3)
continue;
754 std::sort(fcl2hits.begin(), fcl2hits.end(),
SortByLowHit);
757 mf::LogVerbatim(
"CC") <<
" Overlap size " << overlapSize <<
" fit chisq " << clChisq
758 <<
" nhits " << fcl2hits.size();
759 if (clChisq > 20)
continue;
761 std::vector<unsigned int> oHits = fcl2hits;
765 unsigned short jclNewSize;
766 for (jclNewSize = 0; jclNewSize < fcl2hits.size(); ++jclNewSize) {
767 iht = fcl2hits[jclNewSize];
768 if (fHits[iht].
WireID().Wire <= bWire)
break;
771 mf::LogVerbatim(
"CC") <<
"jcl old size " << fcl2hits.size() <<
" newSize " << jclNewSize;
772 iht = fcl2hits[fcl2hits.size() - 1];
773 unsigned int iiht = fcl2hits[jclNewSize - 1];
774 mf::LogVerbatim(
"CC") <<
"jcl old last wire " << fHits[iht].WireID().Wire
775 <<
" After resize last wire " << fHits[iiht].WireID().Wire;
777 fcl2hits.resize(jclNewSize);
779 fcl2hits.insert(fcl2hits.end(), oHits.begin(), oHits.end());
781 for (ii = 0; ii < tcl[icl].tclhits.size(); ++ii) {
782 iht = tcl[icl].tclhits[ii];
783 if ((
unsigned int)fHits[iht].WireID().Wire >= eWire)
continue;
784 fcl2hits.insert(fcl2hits.end(), tcl[icl].tclhits.begin() + ii, tcl[icl].tclhits.end());
787 clBeginSlp = tcl[jcl].BeginSlp;
788 clBeginSlpErr = tcl[jcl].BeginSlpErr;
789 clBeginAng = tcl[jcl].BeginAng;
790 clBeginWir = tcl[jcl].BeginWir;
791 clBeginTim = tcl[jcl].BeginTim;
792 clBeginChg = tcl[jcl].BeginChg;
793 clBeginChgNear = tcl[jcl].BeginChgNear;
795 clEndSlp = tcl[icl].EndSlp;
796 clEndSlpErr = tcl[icl].EndSlpErr;
797 clEndAng = tcl[icl].EndAng;
798 clEndWir = tcl[icl].EndWir;
799 clEndTim = tcl[icl].EndTim;
800 clEndChg = tcl[icl].EndChg;
801 clEndChgNear = tcl[icl].EndChgNear;
802 clStopCode = tcl[icl].StopCode;
803 clProcCode = tcl[icl].ProcCode + 500;
804 MakeClusterObsolete(icl);
805 MakeClusterObsolete(jcl);
808 RestoreObsoleteCluster(icl);
809 RestoreObsoleteCluster(jcl);
810 CheckHitClusterAssociations();
813 CheckHitClusterAssociations();
814 tcl[tcl.size() - 1].BeginVtx = tcl[jcl].BeginVtx;
815 tcl[tcl.size() - 1].EndVtx = tcl[icl].BeginVtx;
818 mf::LogVerbatim(
"CC") <<
"MergeOverlap new long cluster ID " << tcl[tcl.size() - 1].ID
819 <<
" in clCTP " << clCTP;
821 if (tcl[icl].ID < 0)
break;
828 void ClusterCrawlerAlg::MakeClusterObsolete(
unsigned short icl)
830 short& ID = tcl[icl].ID;
832 mf::LogError(
"CC") <<
"Trying to make already-obsolete cluster obsolete ID = " << ID;
838 for (
unsigned int iht = 0; iht < tcl[icl].tclhits.size(); ++iht)
839 inClus[tcl[icl].tclhits[iht]] = 0;
844 void ClusterCrawlerAlg::RestoreObsoleteCluster(
unsigned short icl)
846 short& ID = tcl[icl].ID;
848 mf::LogError(
"CC") <<
"Trying to restore non-obsolete cluster ID = " << ID;
853 for (
unsigned short iht = 0; iht < tcl[icl].tclhits.size(); ++iht)
854 inClus[tcl[icl].tclhits[iht]] = ID;
859 void ClusterCrawlerAlg::FclTrimUS(
unsigned short nTrim)
863 if (nTrim == 0)
return;
865 if (nTrim >= fcl2hits.size()) nTrim = fcl2hits.size();
868 for (
unsigned short ii = 0; ii < nTrim; ++ii) {
878 void ClusterCrawlerAlg::CheckHitClusterAssociations()
882 if (fHits.size() != inClus.size()) {
883 mf::LogError(
"CC") <<
"CHCA: Sizes wrong " << fHits.size() <<
" " << inClus.size();
887 unsigned int iht, nErr = 0;
891 for (
unsigned short icl = 0; icl < tcl.size(); ++icl) {
892 if (tcl[icl].ID < 0)
continue;
894 for (
unsigned short ii = 0; ii < tcl[icl].tclhits.size(); ++ii) {
895 iht = tcl[icl].tclhits[ii];
896 if (iht > fHits.size() - 1) {
897 mf::LogError(
"CC") <<
"CHCA: Bad tclhits index " << iht <<
" fHits size " << fHits.size();
900 if (inClus[iht] != clID) {
901 mf::LogError(
"CC") <<
"CHCA: Bad cluster -> hit association. clID " << clID
902 <<
" hit inClus " << inClus[iht] <<
" ProcCode " << tcl[icl].ProcCode
903 <<
" CTP " << tcl[icl].CTP;
905 if (nErr > 10)
return;
912 for (iht = 0; iht < fHits.size(); ++iht) {
913 if (inClus[iht] <= 0)
continue;
914 icl = inClus[iht] - 1;
916 if (tcl[icl].ID < 0) {
917 mf::LogError(
"CC") <<
"CHCA: Hit associated with an obsolete cluster. hit W:T " 918 << fHits[iht].WireID().Wire <<
":" << (int)fHits[iht].PeakTime()
919 <<
" tcl[icl].ID " << tcl[icl].ID;
921 if (nErr > 10)
return;
923 if (std::find(tcl[icl].tclhits.begin(), tcl[icl].tclhits.end(), iht) ==
924 tcl[icl].tclhits.end()) {
925 mf::LogError(
"CC") <<
"CHCA: Bad hit -> cluster association. hit index " << iht <<
" W:T " 926 << fHits[iht].WireID().Wire <<
":" << (int)fHits[iht].PeakTime()
927 <<
" inClus " << inClus[iht];
929 if (nErr > 10)
return;
936 void ClusterCrawlerAlg::RemoveObsoleteHits()
939 unsigned int destHit = 0;
941 if (fHits.size() != inClus.size()) {
942 mf::LogError(
"CC") <<
"RemoveObsoleteHits size mis-match " << fHits.size() <<
" " 948 for (
unsigned int srcHit = 0; srcHit < fHits.size(); ++srcHit) {
949 if (inClus[srcHit] < 0)
continue;
950 if (srcHit != destHit) {
951 fHits[destHit] = std::move(fHits[srcHit]);
952 inClus[destHit] = inClus[srcHit];
953 if (inClus[destHit] > 0) {
955 icl = inClus[destHit] - 1;
956 auto&
hits = tcl[icl].tclhits;
957 auto iHitIndex = std::find(
hits.begin(),
hits.end(), srcHit);
958 if (iHitIndex ==
hits.end()) {
959 mf::LogError(
"CC") <<
"RemoveObsoleteHits: Hit #" << srcHit
960 <<
" not found in cluster ID " << inClus[destHit];
963 *iHitIndex = destHit;
970 fHits.resize(destHit);
971 inClus.resize(destHit);
976 void ClusterCrawlerAlg::ChkClusterDS()
983 prt = (fDebugPlane == 3);
986 float dThCut = fKinkAngCut[fNumPass - 1];
989 mf::LogVerbatim(
"CC") <<
"ChkClusterDS clCTP " << clCTP <<
" kink angle cut " << dThCut;
991 const unsigned short tclsize = tcl.size();
992 bool didMerge, skipit;
993 unsigned short icl, ii, nhm, iv;
997 for (icl = 0; icl < tclsize; ++icl) {
998 if (tcl[icl].ID < 0)
continue;
999 if (tcl[icl].CTP != clCTP)
continue;
1001 if (tcl[icl].BeginVtx >= 0)
continue;
1004 for (
unsigned short iv = 0; iv < vtx.size(); ++iv) {
1005 if (vtx[iv].CTP != clCTP)
continue;
1006 if (
std::abs(vtx[iv].Wire - tcl[icl].BeginWir) > 4)
continue;
1007 if (
std::abs(vtx[iv].Time - tcl[icl].BeginTim) > 20)
continue;
1011 if (skipit)
continue;
1014 for (ii = 0; ii < 3; ++ii) {
1016 if (fHits[iht].Multiplicity() > 1) {
1017 MergeHits(iht, didMerge);
1018 if (didMerge) ++nhm;
1023 FitClusterMid(icl, 0, 3);
1024 tcl[icl].BeginTim = clpar[0];
1025 tcl[icl].BeginSlp = clpar[1];
1026 tcl[icl].BeginAng = atan(fScaleF * clpar[1]);
1027 tcl[icl].BeginSlpErr = clparerr[1];
1028 tcl[icl].BeginChg = fAveChg;
1029 tcl[icl].ProcCode += 5000;
1030 if (prt)
mf::LogVerbatim(
"CC") <<
"ChkClusterDS: Merge hits on cluster " << tcl[icl].ID;
1034 float thhits, prevth, hitrms, rmsrat;
1036 std::vector<unsigned int> dshits;
1037 for (icl = 0; icl < tclsize; ++icl) {
1038 if (tcl[icl].ID < 0)
continue;
1039 if (tcl[icl].CTP != clCTP)
continue;
1041 if (tcl[icl].BeginVtx >= 0)
continue;
1044 for (iv = 0; iv < vtx.size(); ++iv) {
1045 if (vtx[iv].CTP != clCTP)
continue;
1046 if (
std::abs(vtx[iv].Wire - tcl[icl].BeginWir) > 4)
continue;
1047 if (
std::abs(vtx[iv].Time - tcl[icl].BeginTim) > 20)
continue;
1051 if (skipit)
continue;
1053 unsigned int ih0 = tcl[icl].tclhits[1];
1054 unsigned int ih1 = tcl[icl].tclhits[0];
1055 const float slp = (fHits[ih1].PeakTime() - fHits[ih0].PeakTime()) /
1056 (fHits[ih1].
WireID().Wire - fHits[ih0].WireID().Wire);
1057 prevth = std::atan(fScaleF * slp);
1060 unsigned int wire = fHits[ih0].WireID().Wire;
1061 hitrms = fHits[ih0].RMS();
1062 float time0 = fHits[ih0].PeakTime();
1067 for (ii = 0; ii < 4; ++ii) {
1069 if (wire > fLastWire)
break;
1070 prtime = time0 + slp;
1072 mf::LogVerbatim(
"CC") <<
"ChkClusterDS: Try to extend " << tcl[icl].ID <<
" to W:T " 1073 << wire <<
" hitrms " << hitrms <<
" prevth " << prevth
1074 <<
" prtime " << (int)prtime;
1076 if (WireHitRange[wire].first == -2)
break;
1077 unsigned int firsthit = WireHitRange[wire].first;
1078 unsigned int lasthit = WireHitRange[wire].second;
1079 bool hitAdded =
false;
1080 for (ih1 = firsthit; ih1 < lasthit; ++ih1) {
1081 if (inClus[ih1] != 0)
continue;
1082 if (prtime < fHits[ih1].PeakTimeMinusRMS(5))
continue;
1083 if (prtime > fHits[ih1].PeakTimePlusRMS(5))
continue;
1084 const float slp = (fHits[ih1].PeakTime() - fHits[ih0].PeakTime()) /
1085 (fHits[ih1].
WireID().Wire - fHits[ih0].WireID().Wire);
1086 thhits = std::atan(fScaleF * slp);
1088 mf::LogVerbatim(
"CC") <<
" theta " << thhits <<
" prevth " << prevth <<
" cut " 1090 if (
std::abs(thhits - prevth) > dThCut)
continue;
1093 if (
std::abs(slp) < fLAClusSlopeCut) {
1094 rmsrat = fHits[ih1].RMS() / hitrms;
1095 ratOK = rmsrat > 0.3 && rmsrat < 3;
1100 MergeHits(ih1, didMerge);
1102 if (prt)
mf::LogVerbatim(
"CC") <<
" rmsrat " << rmsrat <<
" OK? " << ratOK;
1107 dshits.push_back(ih1);
1112 mf::LogVerbatim(
"CC") <<
" Add hit " << fHits[ih1].WireID().Wire <<
":" 1113 << (int)fHits[ih1].PeakTime() <<
" rmsrat " << rmsrat;
1118 if (!hitAdded)
break;
1121 if (dshits.size() > 0) {
1127 if (dshits.size() > 1) std::sort(dshits.begin(), dshits.end(),
SortByLowHit);
1131 for (ii = 0; ii < tcl[icl].tclhits.size(); ++ii) {
1133 iht = tcl[icl].tclhits[ii];
1135 fcl2hits.push_back(iht);
1138 pass = fNumPass - 1;
1140 clBeginChg = fAveChg;
1142 MakeClusterObsolete(icl);
1145 mf::LogError(
"CC") <<
"ChkClusterDS TmpStore failed while extending cluster ID " 1149 const size_t newcl = tcl.size() - 1;
1151 tcl[newcl].BeginVtx = tcl[icl].BeginVtx;
1152 tcl[newcl].EndVtx = tcl[icl].EndVtx;
1158 bool ClusterCrawlerAlg::CrawlVtxChk2()
1164 if (vtx.size() == 0)
return false;
1165 if (fcl2hits.size() == 0)
return false;
1167 unsigned int iht = fcl2hits.size() - 1;
1169 float wire0 = (0.5 + fHits[fcl2hits[iht]].WireID().Wire);
1171 float thclus = std::atan(fScaleF * clpar[1]);
1173 for (
unsigned short iv = 0; iv < vtx.size(); ++iv) {
1174 if (vtx[iv].CTP != clCTP)
continue;
1175 if (wire0 < vtx[iv].Wire)
continue;
1176 if (wire0 > vtx[iv].Wire + 10)
continue;
1177 dt = clpar[0] + (vtx[iv].Wire - wire0) * clpar[1] - vtx[iv].Time;
1181 for (icl = 0; icl < tcl.size(); ++icl) {
1182 if (tcl[icl].CTP != clCTP)
continue;
1183 if (tcl[icl].EndVtx != iv)
continue;
1184 if (
std::abs(tcl[icl].EndAng - thclus) < fKinkAngCut[pass])
return true;
1193 bool ClusterCrawlerAlg::CrawlVtxChk(
unsigned int kwire)
1197 if (vtx.size() == 0)
return false;
1198 unsigned int iht = fcl2hits.size() - 1;
1199 float wire0 = (0.5 + fHits[fcl2hits[iht]].WireID().Wire);
1200 float prtime = clpar[0] + (kwire - wire0) * clpar[1];
1201 for (
unsigned short iv = 0; iv < vtx.size(); ++iv) {
1202 if (vtx[iv].CTP != clCTP)
continue;
1203 if ((
unsigned int)(0.5 + vtx[iv].Wire) != kwire)
continue;
1204 if (
std::abs(prtime - vtx[iv].Time) < 10)
return true;
1210 void ClusterCrawlerAlg::VtxConstraint(
unsigned int iwire,
1213 unsigned int& useHit,
1219 doConstrain =
false;
1220 if (vtx.size() == 0)
return;
1222 if (pass == 0)
return;
1224 if (!fFindVertices[pass - 1])
return;
1226 if (jwire > WireHitRange.size() - 1) {
1227 mf::LogError(
"CC") <<
"VtxConstraint fed bad jwire " << jwire <<
" WireHitRange size " 1228 << WireHitRange.size();
1232 unsigned int jfirsthit = WireHitRange[jwire].first;
1233 unsigned int jlasthit = WireHitRange[jwire].second;
1234 for (
unsigned short iv = 0; iv < vtx.size(); ++iv) {
1235 if (vtx[iv].CTP != clCTP)
continue;
1237 if (vtx[iv].Wire > jwire)
continue;
1239 if (vtx[iv].Wire < jwire - 10)
continue;
1242 clpar[1] = (vtx[iv].Time - hit.
PeakTime()) / (vtx[iv].Wire - iwire);
1244 float prtime = clpar[0] + clpar[1] * (jwire - iwire);
1245 for (
unsigned int jhit = jfirsthit; jhit < jlasthit; ++jhit) {
1246 if (inClus[jhit] != 0)
continue;
1247 const float tdiff =
std::abs(fHits[jhit].TimeDistanceAsRMS(prtime));
1259 void ClusterCrawlerAlg::RefineVertexClusters(
unsigned short iv)
1264 std::vector<unsigned short> begClusters;
1265 std::vector<short> begdW;
1266 std::vector<unsigned short> endClusters;
1267 std::vector<short> enddW;
1269 unsigned int vWire = (
unsigned int)(vtx[iv].Wire + 0.5);
1270 unsigned int vWireErr = (
unsigned int)(2 * vtx[iv].WireErr);
1271 unsigned int vWireLo = vWire - vWireErr;
1272 unsigned int vWireHi = vWire + vWireErr;
1274 unsigned short icl, ii;
1276 bool needsWork =
false;
1279 for (icl = 0; icl < tcl.size(); ++icl) {
1280 if (tcl[icl].ID < 0)
continue;
1281 if (tcl[icl].CTP != vtx[iv].CTP)
continue;
1282 if (tcl[icl].BeginVtx == iv) {
1283 dW = vWire - tcl[icl].BeginWir;
1284 if (dW > maxdW) maxdW = dW;
1285 if (dW < mindW) mindW = dW;
1286 if (
std::abs(dW) > 1) needsWork =
true;
1288 begClusters.push_back(icl);
1289 begdW.push_back(dW);
1291 if (tcl[icl].EndVtx == iv) {
1292 dW = vWire - tcl[icl].EndWir;
1293 if (dW > maxdW) maxdW = dW;
1294 if (dW < mindW) mindW = dW;
1295 if (
std::abs(dW) > 1) needsWork =
true;
1296 endClusters.push_back(icl);
1297 enddW.push_back(dW);
1302 mf::LogVerbatim(
"CC") <<
"RefineVertexClusters: vertex " << iv <<
" needsWork " << needsWork
1303 <<
" mindW " << mindW <<
" maxdW " << maxdW <<
" vWireErr " << vWireErr;
1305 if (!needsWork)
return;
1309 if (((
unsigned int)
std::abs(mindW) < vWireErr) && ((
unsigned int)
std::abs(maxdW) < vWireErr) &&
1311 if (vtxprt)
mf::LogVerbatim(
"CC") <<
" Move vtx wire " << vtx[iv].Wire;
1312 vtx[iv].Wire -= (float)(maxdW + mindW) / 2;
1315 vtx[iv].Fixed =
true;
1322 unsigned short newSize;
1323 for (ii = 0; ii < endClusters.size(); ++ii) {
1324 icl = endClusters[ii];
1326 mf::LogVerbatim(
"CC") <<
" endCluster " << tcl[icl].ID <<
" dW " << enddW[ii] <<
" vWire " 1328 if (tcl[icl].EndWir < vWire) {
1331 newSize = fcl2hits.size();
1332 for (
auto hiter = fcl2hits.rbegin(); hiter < fcl2hits.rend(); ++hiter) {
1333 if (fHits[*hiter].
WireID().Wire > vWire)
break;
1337 for (
auto hiter = fcl2hits.begin(); hiter < fcl2hits.end(); ++hiter)
1340 fcl2hits.resize(newSize);
1341 MakeClusterObsolete(icl);
1347 tcl[tcl.size() - 1].EndVtx = iv;
1350 mf::LogVerbatim(
"CC") <<
" modified cluster " << tcl[icl].ID <<
" -> " 1351 << tcl[tcl.size() - 1].ID;
1353 else if (tcl[icl].EndWir > vWire) {
1354 mf::LogVerbatim(
"CC") <<
"RefineVertexClusters: Write some EndVtx code";
1358 if (begClusters.size() > 0)
1359 mf::LogVerbatim(
"CC") <<
"RefineVertexClusters: Write some BeginVtx code";
1361 if (mindW < 0 && maxdW > 0) {
1366 unsigned short clsBigChg = 0;
1371 for (ii = 0; ii < begClusters.size(); ++ii) {
1372 icl = begClusters[ii];
1373 for (iht = 0; iht < tcl[icl].tclhits.size(); ++iht) {
1374 ihit = tcl[icl].tclhits[iht];
1375 if (fHits[ihit].Integral() > bigChg) {
1376 bigChg = fHits[ihit].Integral();
1380 if (fHits[ihit].
WireID().Wire < vWireLo)
break;
1384 for (ii = 0; ii < endClusters.size(); ++ii) {
1385 icl = endClusters[ii];
1386 for (iht = 0; iht < tcl[icl].tclhits.size(); ++iht) {
1387 ihit = tcl[icl].tclhits[tcl[icl].tclhits.size() - 1 - iht];
1388 if (fHits[ihit].Integral() > bigChg) {
1389 bigChg = fHits[ihit].Integral();
1393 if (fHits[ihit].
WireID().Wire > vWireHi)
break;
1398 mf::LogVerbatim(
"CC") <<
" moving vertex location to hit " << fHits[vtxHit].WireID().Wire
1399 <<
":" << (int)fHits[vtxHit].PeakTime() <<
" on cluster " 1400 << tcl[clsBigChg].ID;
1401 vtx[iv].Wire = fHits[vtxHit].WireID().Wire;
1402 vtx[iv].Time = fHits[vtxHit].PeakTime();
1403 vtx[iv].Fixed =
true;
1412 bool ClusterCrawlerAlg::VtxClusterSplit()
1417 vtxprt = (fDebugPlane >= 0 && fDebugWire == 5555);
1421 if (vtx.size() == 0)
return false;
1422 unsigned short tclsize = tcl.size();
1423 if (tclsize < 2)
return false;
1426 bool didSomething =
false;
1428 for (
unsigned short icl = 0; icl < tclsize; ++icl) {
1429 if (tcl[icl].ID < 0)
continue;
1430 if (tcl[icl].CTP != clCTP)
continue;
1432 if (tcl[icl].tclhits.size() < 5)
continue;
1435 for (
unsigned short ivx = 0; ivx < vtx.size(); ++ivx) {
1436 if (vtx[ivx].CTP != clCTP)
continue;
1438 if (vtx[ivx].NClusters == 0)
continue;
1440 if (tcl[icl].BeginVtx == ivx)
continue;
1441 if (tcl[icl].EndVtx == ivx)
continue;
1443 if (vtx[ivx].Wire < tcl[icl].EndWir)
continue;
1444 if (vtx[ivx].Wire > tcl[icl].BeginWir)
continue;
1446 float hiTime = tcl[icl].BeginTim;
1447 if (tcl[icl].EndTim > hiTime) hiTime = tcl[icl].EndTim;
1448 if (vtx[ivx].Time > hiTime + 5)
continue;
1449 float loTime = tcl[icl].BeginTim;
1450 if (tcl[icl].EndTim < loTime) loTime = tcl[icl].EndTim;
1451 if (vtx[ivx].Time < loTime - 5)
continue;
1455 mf::LogVerbatim(
"CC") <<
" Chk cluster ID " << tcl[icl].ID <<
" with vertex " << ivx;
1459 unsigned short nSplit = 0;
1460 unsigned short nLop = 0;
1462 for (
unsigned short ii = tcl[icl].tclhits.size() - 1; ii > 0; --ii) {
1463 iht = tcl[icl].tclhits[ii];
1465 if (fHits[iht].
WireID().Wire >= vtx[ivx].Wire) {
1468 mf::LogVerbatim(
"CC") <<
"Split cluster " << tcl[icl].ID <<
" at wire " 1469 << fHits[iht].WireID().Wire <<
" nSplit " << nSplit;
1475 if (ihvx < 0)
continue;
1476 if (fabs(fHits[ihvx].PeakTime() - vtx[ivx].Time) > 10)
continue;
1481 if (nSplit > tcl[icl].tclhits.size() / 2) { iclAng = tcl[icl].EndAng; }
1483 iclAng = tcl[icl].BeginAng;
1488 for (
unsigned short jcl = 0; jcl < tclsize; ++jcl) {
1489 if (tcl[jcl].ID < 0)
continue;
1490 if (tcl[jcl].CTP != clCTP)
continue;
1491 if (tcl[jcl].BeginVtx == ivx) {
1492 if (fabs(tcl[jcl].BeginAng - iclAng) > 0.4) {
1498 if (tcl[jcl].EndVtx == ivx) {
1499 if (fabs(tcl[jcl].EndAng - iclAng) > 0.4) {
1513 for (
unsigned short ii = 0; ii < nLop; ++ii) {
1514 unsigned int iht = fcl2hits[fcl2hits.size() - 1];
1516 fcl2hits.pop_back();
1521 MakeClusterObsolete(icl);
1523 if (!TmpStore())
continue;
1524 unsigned short newcl = tcl.size() - 1;
1525 tcl[newcl].BeginVtx = tcl[icl].BeginVtx;
1526 tcl[newcl].EndVtx = ivx;
1532 if (SplitCluster(icl, nSplit, ivx)) {
1533 tcl[tcl.size() - 1].ProcCode += 1000;
1534 tcl[tcl.size() - 2].ProcCode += 1000;
1538 didSomething =
true;
1544 return didSomething;
1549 void ClusterCrawlerAlg::MergeHits(
const unsigned int theHit,
bool& didMerge)
1563 if (theHit > fHits.size() - 1) {
return; }
1568 if (fHits[theHit].GoodnessOfFit() == 6666) {
1570 mf::LogVerbatim(
"CC") <<
"MergeHits Trying to merge already merged hit " 1573 <<
" theHit " << theHit;
1578 if (!mergeAvailable[theHit]) {
return; }
1583 std::pair<size_t, size_t> MultipletRange = FindHitMultiplet(theHit);
1586 if (MultipletRange.second <= MultipletRange.first)
return;
1589 unsigned short nAvailable = 0;
1590 unsigned short nInClus = 0;
1591 for (
size_t jht = MultipletRange.first; jht < MultipletRange.second; ++jht) {
1592 if (jht == theHit)
continue;
1593 if (fHits[jht].GoodnessOfFit() == 6666)
continue;
1594 if (inClus[jht] != 0) {
1600 if (nAvailable == 0)
return;
1602 if (nInClus > 0)
return;
1608 short loTime = 9999;
1610 unsigned short nGaus = 1;
1613 unsigned short nclose = 0;
1615 for (
size_t jht = MultipletRange.first; jht < MultipletRange.second; ++jht) {
1616 if (inClus[jht] < 0)
continue;
1624 if (inClus[jht] != 0)
continue;
1626 if (arg < loTime) loTime = arg;
1628 if (arg > hiTime) hiTime = arg;
1629 if (jht != theHit) ++nGaus;
1631 if (jht != theHit && hitSep < 3) ++nclose;
1634 if (nGaus < 2)
return;
1637 const short int NewMultiplicity = hit.
Multiplicity() + 1 - nGaus;
1639 if (loTime < 0) loTime = 0;
1642 std::vector<double> signal(hiTime - loTime, 0.);
1645 for (
size_t jht = MultipletRange.first; jht < MultipletRange.second; ++jht) {
1647 if (jht != theHit) {
1649 if (inClus[jht] != 0)
continue;
1655 for (
unsigned short time = loTime; time < hiTime; ++time) {
1656 unsigned short indx = time - loTime;
1657 double arg = (other_hit.
PeakTime() - (double)time) / other_hit.
RMS();
1658 signal[indx] += other_hit.
PeakAmplitude() * exp(-0.5 * arg * arg);
1663 double sigsumt = 0.;
1664 for (
unsigned short time = loTime; time < hiTime; ++time) {
1665 sigsum += signal[time - loTime];
1666 sigsumt += signal[time - loTime] * time;
1672 double aveTime = sigsumt / sigsum;
1675 for (
unsigned short time = loTime; time < hiTime; ++time) {
1676 double dtime = time - aveTime;
1677 sigsumt += signal[time - loTime] * dtime * dtime;
1679 const float RMS = std::sqrt(sigsumt / sigsum);
1681 const float amplitude = chgsum * chgNorm / (2.507 * RMS);
1703 FixMultipletLocalIndices(MultipletRange.first, MultipletRange.second);
1709 void ClusterCrawlerAlg::FixMultipletLocalIndices(
size_t begin,
1711 short int multiplicity )
1723 if (multiplicity < 0) {
1725 for (
size_t iHit = begin; iHit <
end; ++iHit) {
1726 if (inClus[iHit] < 0)
continue;
1732 short int local_index = 0;
1733 for (
size_t iHit = begin; iHit <
end; ++iHit) {
1734 if (inClus[iHit] < 0)
continue;
1765 void ClusterCrawlerAlg::FindStarVertices()
1772 if (tcl.size() < 2)
return;
1777 vtxprt = (fDebugPlane == (int)plane && fDebugHit < 0);
1783 unsigned short vtxSizeIn = vtx.size();
1787 float dsl = 0, dth = 0;
1788 float es1 = 0, es2 = 0;
1789 float eth1 = 0, eth2 = 0;
1790 float bt1 = 0, bt2 = 0;
1791 float et1 = 0, et2 = 0;
1792 float bw1 = 0, bw2 = 0;
1793 float ew1 = 0, ew2 = 0;
1794 float lotime = 0, hitime = 0, nwirecut = 0;
1795 unsigned short tclsize = tcl.size();
1796 for (
unsigned short it1 = 0; it1 < tclsize - 1; ++it1) {
1798 if (tcl[it1].ID < 0)
continue;
1799 if (tcl[it1].CTP != clCTP)
continue;
1801 if (tcl[it1].tclhits.size() > 100) {
1802 dth = tcl[it1].BeginAng - tcl[it1].EndAng;
1805 es1 = tcl[it1].EndSlp;
1806 eth1 = tcl[it1].EndAng;
1807 ew1 = tcl[it1].EndWir;
1808 et1 = tcl[it1].EndTim;
1809 bw1 = tcl[it1].BeginWir;
1810 bt1 = tcl[it1].BeginTim;
1811 for (
unsigned short it2 = it1 + 1; it2 < tclsize; ++it2) {
1812 if (tcl[it2].ID < 0)
continue;
1813 if (tcl[it2].CTP != clCTP)
continue;
1815 if (tcl[it2].tclhits.size() > 100) {
1816 dth = tcl[it2].BeginAng - tcl[it2].EndAng;
1817 if (
std::abs(dth) < 0.05)
continue;
1819 es2 = tcl[it2].EndSlp;
1820 eth2 = tcl[it2].EndAng;
1821 ew2 = tcl[it2].EndWir;
1822 et2 = tcl[it2].EndTim;
1823 bw2 = tcl[it2].BeginWir;
1824 bt2 = tcl[it2].BeginTim;
1827 if (dth < 0.3)
continue;
1829 fvw = (et1 - ew1 * es1 - et2 + ew2 * es2) / dsl;
1831 if (fvw < ew1 || fvw > bw1)
continue;
1832 if (fvw < ew2 || fvw > bw2)
continue;
1834 mf::LogVerbatim(
"CC") <<
"Chk clusters " << tcl[it1].ID <<
" " << tcl[it2].ID
1835 <<
" topo5 vtx wire " << fvw;
1838 if (tcl[it1].tclhits.size() > tcl[it2].tclhits.size()) {
1839 if (tcl[it1].tclhits.size() > 20) {
1840 nwirecut = 0.5 * tcl[it1].tclhits.size();
1841 if ((fvw - ew1) > nwirecut)
continue;
1845 if (tcl[it2].tclhits.size() > 20) {
1846 nwirecut = 0.5 * tcl[it2].tclhits.size();
1847 if ((fvw - ew2) > nwirecut)
continue;
1850 fvt = et1 + (fvw - ew1) * es1;
1853 if (et1 < lotime) lotime = et1;
1854 if (bt1 < lotime) lotime = bt1;
1855 if (et2 < lotime) lotime = et2;
1856 if (bt2 < lotime) lotime = bt2;
1858 if (et1 > hitime) hitime = et1;
1859 if (bt1 > hitime) hitime = bt1;
1860 if (et2 > hitime) hitime = et2;
1861 if (bt2 > hitime) hitime = bt2;
1862 if (fvt < lotime || fvt > hitime)
continue;
1864 mf::LogVerbatim(
"CC") <<
" vertex time " << fvt <<
" lotime " << lotime <<
" hitime " 1866 unsigned int vw = (0.5 + fvw);
1868 unsigned int pos1 = 0;
1869 for (
unsigned short ii = 0; ii < tcl[it1].tclhits.size(); ++ii) {
1870 if (fHits[tcl[it1].tclhits[ii]].
WireID().Wire <= vw) {
1871 if (
std::abs(
int(fHits[tcl[it1].tclhits[ii]].PeakTime() - fvt)) < 10) pos1 = ii;
1876 if (pos1 == 0)
continue;
1877 unsigned short pos2 = 0;
1878 for (
unsigned short ii = 0; ii < tcl[it2].tclhits.size(); ++ii) {
1879 if (fHits[tcl[it2].tclhits[ii]].
WireID().Wire <= vw) {
1880 if (
std::abs(
int(fHits[tcl[it2].tclhits[ii]].PeakTime() - fvt)) < 10) pos2 = ii;
1885 if (pos2 == 0)
continue;
1887 for (
unsigned short iv = 0; iv < vtx.size(); ++iv) {
1888 if (
std::abs(fvw - vtx[iv].Wire) < 2 &&
std::abs(fvt - vtx[iv].Time) < 10)
continue;
1898 newvx.
Fixed =
false;
1899 vtx.push_back(newvx);
1900 unsigned short ivx = vtx.size() - 1;
1902 mf::LogVerbatim(
"CC") <<
" new vertex " << ivx <<
" cluster " << tcl[it1].ID
1903 <<
" split pos " << pos1;
1904 if (!SplitCluster(it1, pos1, ivx))
continue;
1905 tcl[tcl.size() - 1].ProcCode += 1000;
1906 tcl[tcl.size() - 2].ProcCode += 1000;
1908 mf::LogVerbatim(
"CC") <<
" new vertex " << ivx <<
" cluster " << tcl[it2].ID
1909 <<
" split pos " << pos2;
1910 if (!SplitCluster(it2, pos2, ivx))
continue;
1911 tcl[tcl.size() - 1].ProcCode += 1000;
1912 tcl[tcl.size() - 2].ProcCode += 1000;
1918 if (vtx.size() > vtxSizeIn) {
1922 VertexCluster(vtx.size() - 1);
1928 for (
unsigned short iv = 0; iv < vtx.size(); ++iv) {
1929 if (vtx[iv].CTP != clCTP)
continue;
1930 mf::LogVerbatim(
"CC") <<
"vtx " << iv <<
" wire " << vtx[iv].Wire <<
" time " 1931 << (int)vtx[iv].Time <<
" NClusters " << vtx[iv].NClusters <<
" topo " 1940 void ClusterCrawlerAlg::VertexCluster(
unsigned short iv)
1943 if (vtx[iv].NClusters == 0)
return;
1945 short dwb, dwe, dtb, dte;
1948 for (
unsigned short icl = 0; icl < tcl.size(); ++icl) {
1949 if (tcl[icl].ID < 0)
continue;
1950 if (tcl[icl].CTP != vtx[iv].CTP)
continue;
1952 dwb = vtx[iv].Wire - tcl[icl].BeginWir;
1953 dtb = vtx[iv].Time - tcl[icl].BeginTim;
1954 dwe = vtx[iv].Wire - tcl[icl].EndWir;
1955 dte = vtx[iv].Time - tcl[icl].EndTim;
1957 float drb = dwb * dwb + dtb * dtb;
1958 float dre = dwe * dwe + dte * dte;
1960 bool bCloser = (drb < dre);
1964 if (tcl[icl].BeginChgNear > fChgNearCut)
continue;
1967 if (tcl[icl].EndChgNear > fChgNearCut)
continue;
1971 mf::LogVerbatim(
"CC") <<
"VertexCluster: Try icl ID " << tcl[icl].ID <<
" w vtx " << iv
1972 <<
" dwb " << dwb <<
" dwe " << dwe <<
" drb " << drb <<
" dre " 1973 << dre <<
" Begin closer? " << bCloser;
1975 if (tcl[icl].BeginVtx < 0 && bCloser && dwb > -3 && dwb < 3 && tcl[icl].EndVtx != iv) {
1976 sigOK = ChkSignal(tcl[icl].BeginWir, tcl[icl].BeginTim, vtx[iv].Wire, vtx[iv].Time);
1978 mf::LogVerbatim(
"CC") <<
" Attach cluster Begin to vtx? " << iv <<
" sigOK " << sigOK;
1981 mf::LogVerbatim(
"CC") <<
" check ClusterVertexChi " << ClusterVertexChi(icl, 0, iv);
1982 if (ClusterVertexChi(icl, 0, iv) < fVertex2DCut) {
1984 tcl[icl].BeginVtx = iv;
1986 if (vtx[iv].ChiDOF > fVertex2DCut || vtx[iv].WireErr > fVertex2DWireErrCut) {
1987 tcl[icl].BeginVtx = -99;
1994 if (tcl[icl].EndVtx < 0 && !bCloser && dwe > -3 && dwe < 3 && tcl[icl].BeginVtx != iv) {
1995 sigOK = ChkSignal(tcl[icl].EndWir, tcl[icl].EndTim, vtx[iv].Wire, vtx[iv].Time);
1997 mf::LogVerbatim(
"CC") <<
" Attach cluster End to vtx? " << iv <<
" sigOK " << sigOK;
2000 mf::LogVerbatim(
"CC") <<
" check ClusterVertexChi " << ClusterVertexChi(icl, 1, iv);
2001 if (ClusterVertexChi(icl, 1, iv) < 3) {
2003 tcl[icl].EndVtx = iv;
2005 if (vtx[iv].ChiDOF > fVertex2DCut || vtx[iv].WireErr > fVertex2DWireErrCut) {
2006 tcl[icl].EndVtx = -99;
2016 void ClusterCrawlerAlg::FitAllVtx(
CTP_t inCTP)
2019 for (
unsigned short iv = 0; iv < vtx.size(); ++iv) {
2020 if (vtx[iv].CTP != inCTP)
continue;
2027 void ClusterCrawlerAlg::FindVertices()
2031 if (tcl.size() < 2)
return;
2034 std::vector<CluLen> clulens;
2036 for (
unsigned short icl = 0; icl < tcl.size(); ++icl) {
2037 if (tcl[icl].ID < 0)
continue;
2038 if (tcl[icl].CTP != clCTP)
continue;
2039 if (tcl[icl].BeginVtx >= 0 && tcl[icl].EndVtx >= 0)
continue;
2041 clulen.length = tcl[icl].tclhits.size();
2042 clulens.push_back(clulen);
2044 if (
empty(clulens))
return;
2045 std::sort(clulens.begin(), clulens.end(),
greaterThan);
2047 float nwires = fNumWires;
2048 float maxtime = fMaxTime;
2050 unsigned short vtxSizeIn = vtx.size();
2052 vtxprt = (fDebugPlane == (short)plane && fDebugHit < 0);
2054 mf::LogVerbatim(
"CC") <<
"FindVertices plane " << plane <<
" pass " << pass;
2058 float es1 = 0, eth1 = 0, ew1 = 0, et1 = 0;
2059 float bs1 = 0, bth1 = 0, bw1 = 0, bt1 = 0;
2060 float es2 = 0, eth2 = 0, ew2 = 0, et2 = 0;
2061 float bs2 = 0, bth2 = 0, bw2 = 0, bt2 = 0;
2062 float dth = 0, dsl = 0, fvw = 0, fvt = 0;
2064 unsigned int vw = 0, pass1, pass2, ii1, it1, ii2, it2;
2066 for (ii1 = 0; ii1 < clulens.size() - 1; ++ii1) {
2067 it1 = clulens[ii1].index;
2068 es1 = tcl[it1].EndSlp;
2069 eth1 = tcl[it1].EndAng;
2070 ew1 = tcl[it1].EndWir;
2071 et1 = tcl[it1].EndTim;
2072 bs1 = tcl[it1].BeginSlp;
2073 bth1 = tcl[it1].BeginAng;
2074 bw1 = tcl[it1].BeginWir;
2075 bt1 = tcl[it1].BeginTim;
2076 pass1 = tcl[it1].ProcCode - 10 * (tcl[it1].ProcCode / 10);
2077 for (ii2 = ii1 + 1; ii2 < clulens.size(); ++ii2) {
2078 it2 = clulens[ii2].index;
2082 if (tcl[it1].tclhits.size() < 5 && tcl[it2].tclhits.size() < 5)
continue;
2083 es2 = tcl[it2].EndSlp;
2084 eth2 = tcl[it2].EndAng;
2085 ew2 = tcl[it2].EndWir;
2086 et2 = tcl[it2].EndTim;
2087 bs2 = tcl[it2].BeginSlp;
2088 bth2 = tcl[it2].BeginAng;
2089 bw2 = tcl[it2].BeginWir;
2090 bt2 = tcl[it2].BeginTim;
2091 pass2 = tcl[it2].ProcCode - 10 * (tcl[it2].ProcCode / 10);
2092 if (pass1 < pass2) { angcut = fKinkAngCut[pass2]; }
2094 angcut = fKinkAngCut[pass1];
2097 dth = fabs(eth1 - eth2);
2098 if (tcl[it1].EndVtx < 0 && tcl[it2].EndVtx < 0 && dth > 0.1) {
2101 fvw = (et1 - ew1 * es1 - et2 + ew2 * es2) / dsl;
2103 if (fvw > 0. && fvw < nwires) {
2108 if (vw >= fFirstWire - 1 && fvw <= ew1 + 3 && fvw <= ew2 + 3 && fvw > (ew1 - 10) &&
2110 fvt = et1 + (fvw - ew1) * es1;
2113 <<
"Chk clusters topo1 " << tcl[it1].ID <<
" " << tcl[it2].ID <<
" vtx wire " 2114 << vw <<
" time " << (int)fvt <<
" dth " << dth;
2115 if (fvt > 0 && fvt < maxtime) {
2119 SigOK = ChkSignal(vw, fvt, vw, fvt);
2123 SigOK = ChkSignal(vw, fvt, vw, fvt);
2126 if (SigOK) ChkVertex(fvw, fvt, it1, it2, 1);
2133 if (tcl[it1].EndVtx < 0 && tcl[it2].BeginVtx < 0 && dth > angcut) {
2135 if (fabs(ew1 - bw2) < 3 && fabs(et1 - bt2) < 500) { fvw = 0.5 * (ew1 + bw2); }
2137 fvw = (et1 - ew1 * es1 - bt2 + bw2 * bs2) / dsl;
2139 if (fvw > 0 && fvw < nwires) {
2144 if (fvw <= ew1 + 2 && fvw >= bw2 - 2) {
2145 fvt = et1 + (vw - ew1) * es1;
2146 fvt += bt2 + (vw - bw2) * bs2;
2150 <<
"Chk clusters topo2 " << tcl[it1].ID <<
" " << tcl[it2].ID <<
" vtx wire " 2151 << vw <<
" time " << (int)fvt <<
" dth " << dth;
2152 if (fvt > 0. && fvt < maxtime) {
2153 ChkVertex(fvw, fvt, it1, it2, 2);
2160 if (tcl[it1].BeginVtx < 0 && tcl[it2].EndVtx < 0 && dth > angcut) {
2162 if (fabs(bw1 - ew2) < 3 && fabs(bt1 - et2) < 500) { fvw = 0.5 * (bw1 + ew2); }
2164 fvw = (et2 - ew2 * es2 - bt1 + bw1 * bs1) / dsl;
2166 if (fvw > 0 && fvw < nwires) {
2170 if (fvw <= ew2 + 2 && fvw >= bw1 - 2) {
2171 fvt = et2 + (fvw - ew2) * es2;
2172 fvt += bt1 + (fvw - bw1) * es1;
2176 <<
"Chk clusters topo3 " << tcl[it1].ID <<
" " << tcl[it2].ID <<
" vtx wire " 2177 << vw <<
" time " << (int)fvt <<
" dth " << dth;
2178 if (fvt > 0. && fvt < maxtime) {
2179 ChkVertex(fvw, fvt, it1, it2, 3);
2186 if (tcl[it1].BeginVtx < 0 && tcl[it2].BeginVtx < 0 && dth > 0.1) {
2190 fvw = (bt1 - bw1 * bs1 - bt2 + bw2 * bs2) / dsl;
2192 mf::LogVerbatim(
"CC") <<
"Chk clusters topo4 " << bw1 <<
":" << (int)bt1 <<
" " << bw2
2193 <<
":" << (
int)bt2 <<
" fvw " << fvw <<
" nwires " << nwires;
2194 if (fvw > 0 && fvw < nwires) {
2200 if (tcl[it1].tclhits.size() < 2 * dwcut) dwcut = tcl[it1].tclhits.size() / 2;
2201 if (tcl[it2].tclhits.size() < 2 * dwcut) dwcut = tcl[it2].tclhits.size() / 2;
2202 if (fvw <= fLastWire + 1 && fvw >= bw1 - dwcut && fvw <= bw1 + dwcut &&
2203 fvw >= bw2 - dwcut && fvw <= bw1 + dwcut) {
2204 fvt = bt1 + (fvw - bw1) * bs1;
2207 <<
" vtx wire " << vw <<
" time " << fvt <<
" dth " << dth <<
" dwcut " << dwcut;
2208 if (fvt > 0. && fvt < maxtime) {
2212 SigOK = ChkSignal(vw, fvt, vw, fvt);
2216 SigOK = ChkSignal(vw, fvt, vw, fvt);
2219 if (SigOK) ChkVertex(fvw, fvt, it1, it2, 4);
2229 for (
unsigned short it = 0; it < tcl.size(); ++it) {
2230 if (tcl[it].ID < 0)
continue;
2231 if (tcl[it].CTP != clCTP)
continue;
2232 if (tcl[it].BeginVtx > -1 && tcl[it].BeginVtx == tcl[it].EndVtx) {
2233 unsigned short iv = tcl[it].BeginVtx;
2234 float dwir = tcl[it].BeginWir - vtx[iv].Wire;
2235 float dtim = tcl[it].BeginTim - vtx[iv].Time;
2236 float rbeg = dwir * dwir + dtim * dtim;
2237 dwir = tcl[it].EndWir - vtx[iv].Wire;
2238 dtim = tcl[it].EndTim - vtx[iv].Time;
2239 float rend = dwir * dwir + dtim * dtim;
2240 if (rend < rbeg) { tcl[it].EndVtx = -99; }
2242 tcl[it].BeginVtx = -99;
2247 if (vtx.size() > vtxSizeIn) FitAllVtx(clCTP);
2250 for (
unsigned short ivx = 0; ivx < vtx.size(); ++ivx) {
2251 if (vtx[ivx].CTP != clCTP)
continue;
2252 if (vtx[ivx].NClusters == 1) {
2253 vtx[ivx].NClusters = 0;
2254 for (
unsigned short icl = 0; icl < tcl.size(); ++icl) {
2255 if (tcl[icl].CTP != clCTP)
continue;
2256 if (tcl[icl].BeginVtx == ivx) {
2257 tcl[icl].BeginVtx = -99;
2260 if (tcl[icl].EndVtx == ivx) {
2261 tcl[icl].EndVtx = -99;
2271 void ClusterCrawlerAlg::ClusterVertex(
unsigned short it)
2275 if (vtx.size() == 0)
return;
2277 unsigned short iv, jv;
2278 short dwib, dwjb, dwie, dwje;
2279 bool matchEnd, matchBegin;
2281 for (iv = 0; iv < vtx.size(); ++iv) {
2283 if (vtx[iv].CTP != clCTP)
continue;
2285 if (vtx[iv].NClusters == 0)
continue;
2289 if (tcl[it].tclhits.size() < 6) {
2291 dwib =
std::abs(vtx[iv].Wire - tcl[it].BeginWir);
2292 if (dwib > 2) dwib = 2;
2293 dwie =
std::abs(vtx[iv].Wire - tcl[it].EndWir);
2294 if (dwie > 2) dwie = 2;
2297 for (jv = 0; jv < vtx.size(); ++jv) {
2298 if (iv == jv)
continue;
2299 if (
std::abs(vtx[jv].Time - tcl[it].BeginTim) < 50) {
2300 if (
std::abs(vtx[jv].Wire - tcl[it].BeginWir) < dwjb)
2301 dwjb =
std::abs(vtx[jv].Wire - tcl[it].BeginWir);
2303 if (
std::abs(vtx[jv].Time - tcl[it].EndTim) < 50) {
2304 if (
std::abs(vtx[jv].Wire - tcl[it].EndWir) < dwje)
2305 dwje =
std::abs(vtx[jv].Wire - tcl[it].EndWir);
2308 matchEnd = tcl[it].EndVtx != iv && dwie < 3 && dwie < dwje && dwie < dwib;
2309 matchBegin = tcl[it].BeginVtx != iv && dwib < 3 && dwib < dwjb && dwib < dwie;
2312 matchEnd = tcl[it].EndVtx < 0 && vtx[iv].Wire <= tcl[it].EndWir + 2;
2313 matchBegin = tcl[it].BeginVtx < 0 && vtx[iv].Wire >= tcl[it].BeginWir - 2;
2317 mf::LogVerbatim(
"CC") <<
" Match End chi " << ClusterVertexChi(it, 1, iv) <<
" to vtx " 2318 << iv <<
" signalOK " 2320 vtx[iv].Wire, vtx[iv].Time, tcl[it].EndWir, tcl[it].EndTim);
2321 if (ClusterVertexChi(it, 1, iv) < fVertex2DCut &&
2322 ChkSignal(vtx[iv].Wire, vtx[iv].Time, tcl[it].EndWir, tcl[it].EndTim)) {
2324 tcl[it].EndVtx = iv;
2328 mf::LogVerbatim(
"CC") <<
" Add End " << tcl[it].ID <<
" to vtx " << iv <<
" NClusters " 2329 << vtx[iv].NClusters;
2330 if (vtx[iv].ChiDOF < fVertex2DCut && vtx[iv].WireErr < fVertex2DWireErrCut)
return;
2332 mf::LogVerbatim(
"CC") <<
" Bad fit. ChiDOF " << vtx[iv].ChiDOF <<
" WireErr " 2333 << vtx[iv].WireErr <<
" Undo it";
2334 tcl[it].EndVtx = -99;
2341 mf::LogVerbatim(
"CC") <<
" Match Begin chi " << ClusterVertexChi(it, 0, iv) <<
" to vtx " 2342 << iv <<
" signalOK " 2343 << ChkSignal(vtx[iv].Wire,
2347 if (ClusterVertexChi(it, 0, iv) < fVertex2DCut &&
2348 ChkSignal(vtx[iv].Wire, vtx[iv].Time, tcl[it].BeginWir, tcl[it].BeginTim)) {
2350 tcl[it].BeginVtx = iv;
2354 mf::LogVerbatim(
"CC") <<
" Add Begin " << tcl[it].ID <<
" to vtx " << iv
2355 <<
" NClusters " << vtx[iv].NClusters;
2356 if (vtx[iv].ChiDOF < fVertex2DCut && vtx[iv].WireErr < fVertex2DWireErrCut)
return;
2358 mf::LogVerbatim(
"CC") <<
" Bad fit. ChiDOF " << vtx[iv].ChiDOF <<
" WireErr " 2359 << vtx[iv].WireErr <<
" Undo it";
2360 tcl[it].BeginVtx = -99;
2368 void ClusterCrawlerAlg::ChkVertex(
float fvw,
2383 mf::LogVerbatim(
"CC") <<
" ChkVertex " << tcl[it1].EndWir <<
":" << (int)tcl[it1].EndTim
2384 <<
" - " << tcl[it1].BeginWir <<
":" << (
int)tcl[it1].BeginTim
2385 <<
" and " << tcl[it2].EndWir <<
":" << (int)tcl[it2].EndTim <<
" - " 2386 << tcl[it2].BeginWir <<
":" << (
int)tcl[it2].BeginTim <<
" topo " 2387 << topo <<
" fvw " << fvw <<
" fvt " << fvt;
2389 unsigned int vw = (
unsigned int)(0.5 + fvw);
2395 if (tcl[it1].tclhits.size() < 10 && tcl[it2].tclhits.size() < 10) {
2396 for (
unsigned short it = 0; it < tcl.size(); ++it) {
2397 if (it == it1 || it == it2)
continue;
2398 if (tcl[it].ID < 0)
continue;
2399 if (tcl[it].CTP != clCTP)
continue;
2400 if (tcl[it].tclhits.size() < 100)
continue;
2401 if (
std::abs(tcl[it].BeginAng - tcl[it].EndAng) > 0.1)
continue;
2403 if (vw < tcl[it].EndWir + 5)
continue;
2404 if (vw > tcl[it].BeginWir - 5)
continue;
2405 for (
unsigned short ii = 0; ii < tcl[it].tclhits.size(); ++ii) {
2406 iht = tcl[it].tclhits[ii];
2407 if (fHits[iht].
WireID().Wire <= vw) {
2408 if (
std::abs(fHits[iht].PeakTime() - fvt) < 60)
return;
2415 unsigned short nFitOK = 0;
2416 for (
unsigned short iv = 0; iv < vtx.size(); ++iv) {
2417 if (vtx[iv].CTP != clCTP)
continue;
2419 if (PointVertexChi(fvw, fvt, iv) > 300)
continue;
2422 << PointVertexChi(fvw, fvt, iv);
2424 if ((topo == 1 || topo == 2) && tcl[it1].EndVtx < 0) {
2425 if (ChkSignal(vw, fvt, tcl[it1].EndWir, tcl[it1].EndTim)) {
2427 tcl[it1].EndVtx = iv;
2430 mf::LogVerbatim(
"CC") <<
" FitVtx " << iv <<
" WireErr " << vtx[iv].WireErr
2431 <<
" ChiDOF " << vtx[iv].ChiDOF;
2432 if (vtx[iv].WireErr < fVertex2DWireErrCut && vtx[iv].ChiDOF < 5) { ++nFitOK; }
2435 tcl[it1].EndVtx = -99;
2440 else if ((topo == 3 || topo == 4) && tcl[it1].BeginVtx < 0) {
2441 if (ChkSignal(vw, fvt, tcl[it1].BeginWir, tcl[it1].BeginTim)) {
2442 tcl[it1].BeginVtx = iv;
2445 mf::LogVerbatim(
"CC") <<
" FitVtx " << iv <<
" WireErr " << vtx[iv].WireErr
2446 <<
" ChiDOF " << vtx[iv].ChiDOF;
2447 if (vtx[iv].WireErr < fVertex2DWireErrCut && vtx[iv].ChiDOF < 5) { ++nFitOK; }
2450 tcl[it1].BeginVtx = -99;
2455 if ((topo == 1 || topo == 3) && tcl[it2].EndVtx < 0) {
2456 if (ChkSignal(vw, fvt, tcl[it2].EndWir, tcl[it2].EndTim)) {
2457 tcl[it2].EndVtx = iv;
2460 mf::LogVerbatim(
"CC") <<
" FitVtx " << iv <<
" WireErr " << vtx[iv].WireErr
2461 <<
" ChiDOF " << vtx[iv].ChiDOF;
2462 if (vtx[iv].WireErr < fVertex2DWireErrCut && vtx[iv].ChiDOF < 5) { ++nFitOK; }
2465 tcl[it2].EndVtx = -99;
2470 else if ((topo == 2 || topo == 4) && tcl[it2].BeginVtx < 0) {
2471 if (ChkSignal(vw, fvt, tcl[it2].BeginWir, tcl[it2].BeginTim)) {
2472 tcl[it2].BeginVtx = iv;
2475 mf::LogVerbatim(
"CC") <<
" FitVtx " << iv <<
" WireErr " << vtx[iv].WireErr
2476 <<
" ChiDOF " << vtx[iv].ChiDOF;
2477 if (vtx[iv].WireErr < fVertex2DWireErrCut && vtx[iv].ChiDOF < 5) { ++nFitOK; }
2480 tcl[it2].BeginVtx = -99;
2486 if (vtxprt)
mf::LogVerbatim(
"CC") <<
" Attached " << nFitOK <<
" clusters to vertex " << iv;
2494 if (topo == 1 || topo == 2) { SigOK = ChkSignal(vw, fvt, tcl[it1].EndWir, tcl[it1].EndTim); }
2496 SigOK = ChkSignal(vw, fvt, tcl[it1].BeginWir, tcl[it1].BeginTim);
2500 if (topo == 1 || topo == 3) { SigOK = ChkSignal(vw, fvt, tcl[it2].EndWir, tcl[it2].EndTim); }
2502 SigOK = ChkSignal(vw, fvt, tcl[it2].BeginWir, tcl[it2].BeginTim);
2511 newvx.
Fixed =
false;
2512 vtx.push_back(newvx);
2513 unsigned short iv = vtx.size() - 1;
2514 if (topo == 1 || topo == 2) {
2515 if (tcl[it1].EndVtx >= 0) {
2519 tcl[it1].EndVtx = iv;
2522 if (tcl[it1].BeginVtx >= 0) {
2526 tcl[it1].BeginVtx = iv;
2528 if (topo == 1 || topo == 3) {
2529 if (tcl[it2].EndVtx >= 0) {
2533 tcl[it2].EndVtx = iv;
2536 if (tcl[it2].BeginVtx >= 0) {
2540 tcl[it2].BeginVtx = iv;
2545 if (vtx[iv].ChiDOF < 5 && vtx[iv].WireErr < fVertex2DWireErrCut && vtx[iv].TimeErr < 20) {
2547 mf::LogVerbatim(
"CC") <<
" New vtx " << iv <<
" in plane " << plane <<
" topo " << topo
2548 <<
" cls " << tcl[it1].ID <<
" " << tcl[it2].ID <<
" W:T " << (int)vw
2549 <<
":" << (
int)fvt <<
" NClusters " << vtx[iv].NClusters;
2551 if (fRefineVertexClusters) RefineVertexClusters(iv);
2555 mf::LogVerbatim(
"CC") <<
" Bad vtx fit " << vtx[iv].ChiDOF <<
" wire err " 2556 << vtx[iv].WireErr <<
" TimeErr " << vtx[iv].TimeErr;
2559 if (tcl[it1].BeginVtx == iv) tcl[it1].BeginVtx = -99;
2560 if (tcl[it1].EndVtx == iv) tcl[it1].EndVtx = -99;
2561 if (tcl[it2].BeginVtx == iv) tcl[it2].BeginVtx = -99;
2562 if (tcl[it2].EndVtx == iv) tcl[it2].EndVtx = -99;
2568 bool ClusterCrawlerAlg::ChkSignal(
unsigned int iht,
unsigned int jht)
2570 if (iht > fHits.size() - 1)
return false;
2571 if (jht > fHits.size() - 1)
return false;
2572 unsigned int wire1 = fHits[iht].WireID().Wire;
2573 float time1 = fHits[iht].PeakTime();
2574 unsigned int wire2 = fHits[jht].WireID().Wire;
2575 float time2 = fHits[jht].PeakTime();
2576 return ChkSignal(wire1, time1, wire2, time2);
2580 bool ClusterCrawlerAlg::ChkSignal(
unsigned int wire1,
2590 const float gausAmp[20] = {1, 0.99, 0.96, 0.90, 0.84, 0.75, 0.67, 0.58, 0.49, 0.40,
2591 0.32, 0.26, 0.20, 0.15, 0.11, 0.08, 0.06, 0.04, 0.03, 0.02};
2594 if (fMinAmp[plane] <= 0)
return true;
2597 unsigned int wireb = wire1;
2598 float timeb = time1;
2599 unsigned int wiree = wire2;
2600 float timee = time2;
2602 if (wiree > wireb) {
2608 if (wiree < fFirstWire || wiree > fLastWire)
return false;
2609 if (wireb < fFirstWire || wireb > fLastWire)
return false;
2614 bool oneWire =
false;
2615 float prTime, prTimeLo = 0, prTimeHi = 0;
2616 if (wireb == wiree) {
2618 if (time1 < time2) {
2628 slope = (timeb - timee) / (wireb - wiree);
2632 unsigned short nmissed = 0;
2633 for (
unsigned int wire = wiree; wire < wireb + 1; ++wire) {
2634 if (oneWire) { prTime = (prTimeLo + prTimeHi) / 2; }
2636 prTime = timee + (wire - wire0) * slope;
2639 if (WireHitRange[wire].first == -1)
continue;
2641 if (WireHitRange[wire].first == -2) ++nmissed;
2642 unsigned int firsthit = WireHitRange[wire].first;
2643 unsigned int lasthit = WireHitRange[wire].second;
2645 for (
unsigned int khit = firsthit; khit < lasthit; ++khit) {
2649 if (prTime < fHits[khit].StartTick())
continue;
2650 if (prTime > fHits[khit].EndTick())
continue;
2655 if (fHits[khit].PeakTime() - prTime > 500)
continue;
2656 bin =
std::abs(fHits[khit].PeakTime() - prTime) / fHits[khit].RMS();
2658 if (bin > 19)
continue;
2659 if (bin < 0)
continue;
2661 amp += fHits[khit].PeakAmplitude() * gausAmp[
bin];
2664 if (amp < fMinAmp[plane]) ++nmissed;
2666 if (nmissed > fAllowNoHitWire)
return false;
2672 bool ClusterCrawlerAlg::SplitCluster(
unsigned short icl,
unsigned short pos,
unsigned short ivx)
2676 if (tcl[icl].ID < 0)
return false;
2679 MakeClusterObsolete(icl);
2681 unsigned short ii, iclnew;
2687 clBeginSlp = tcl[icl].BeginSlp;
2688 clBeginSlpErr = tcl[icl].BeginSlpErr;
2689 clBeginAng = tcl[icl].BeginAng;
2690 clBeginWir = tcl[icl].BeginWir;
2691 clBeginTim = tcl[icl].BeginTim;
2692 clBeginChg = tcl[icl].BeginChg;
2694 clProcCode = tcl[icl].ProcCode;
2699 for (ii = 0; ii < pos; ++ii) {
2700 iht = tcl[icl].tclhits[ii];
2701 fcl2hits.push_back(iht);
2704 pass = tcl[icl].ProcCode - 10 * (tcl[icl].ProcCode / 10);
2705 if (pass > fNumPass - 1) pass = fNumPass - 1;
2708 clEndSlp = clpar[1];
2709 clEndSlpErr = clparerr[1];
2710 clEndAng = std::atan(fScaleF * clEndSlp);
2715 RestoreObsoleteCluster(icl);
2719 iclnew = tcl.size() - 1;
2720 tcl[iclnew].EndVtx = ivx;
2721 tcl[iclnew].BeginVtx = tcl[icl].BeginVtx;
2723 mf::LogVerbatim(
"CC") <<
"SplitCluster made cluster " << iclnew <<
" attached to Begin vtx " 2727 if (pos < tcl[icl].tclhits.size() - 3) {
2730 clEndSlp = tcl[icl].EndSlp;
2731 clEndSlpErr = tcl[icl].EndSlpErr;
2732 clEndAng = std::atan(fScaleF * clEndSlp);
2733 clEndWir = tcl[icl].EndWir;
2734 clEndTim = tcl[icl].EndTim;
2735 clEndChg = tcl[icl].EndChg;
2737 clProcCode = tcl[icl].ProcCode;
2742 bool didFit =
false;
2743 for (ii = pos; ii < tcl[icl].tclhits.size(); ++ii) {
2744 iht = tcl[icl].tclhits[ii];
2745 if (inClus[iht] != 0) {
2746 RestoreObsoleteCluster(icl);
2749 fcl2hits.push_back(iht);
2751 if (fcl2hits.size() == fMaxHitsFit[pass] || fcl2hits.size() == fMinHits[pass]) {
2753 clBeginSlp = clpar[1];
2754 clBeginAng = std::atan(fScaleF * clBeginSlp);
2755 clBeginSlpErr = clparerr[1];
2758 if ((
unsigned short)fcl2hits.size() == fNHitsAve[pass] + 1) {
2760 clBeginChg = fAveChg;
2768 clBeginChg = fAveChg;
2772 MakeClusterObsolete(tcl.size() - 1);
2773 RestoreObsoleteCluster(icl);
2777 iclnew = tcl.size() - 1;
2778 tcl[iclnew].BeginVtx = ivx;
2779 tcl[iclnew].EndVtx = tcl[icl].EndVtx;
2781 mf::LogVerbatim(
"CC") <<
"SplitCluster made cluster " << iclnew <<
" attached to End vtx " 2789 void ClusterCrawlerAlg::ChkMerge()
2794 if (tcl.size() < 2)
return;
2799 prt = (fDebugPlane == (short)plane && fDebugWire < 0);
2802 unsigned short it1, it2, nh1, pass1, pass2;
2803 float bs1, bth1, bt1, bc1, es1, eth1, et1, ec1;
2804 float bs2, bth2, bt2, bc2, es2, eth2, et2, ec2;
2805 int bw1, ew1, bw2, ew2, ndead;
2806 float dth, dw, angcut, chgrat, chgcut, dtim, timecut, bigslp;
2807 bool bothLong, NoVtx;
2809 int skipcut, tclsize = tcl.size();
2812 for (it1 = 0; it1 < tclsize - 1; ++it1) {
2814 if (tcl[it1].ID < 0)
continue;
2815 if (tcl[it1].CTP != clCTP)
continue;
2816 bs1 = tcl[it1].BeginSlp;
2818 bth1 = std::atan(fScaleF * bs1);
2820 bw1 = tcl[it1].BeginWir;
2821 bt1 = tcl[it1].BeginTim;
2822 bc1 = tcl[it1].BeginChg;
2823 es1 = tcl[it1].EndSlp;
2824 eth1 = tcl[it1].EndAng;
2825 ew1 = tcl[it1].EndWir;
2826 et1 = tcl[it1].EndTim;
2827 ec1 = tcl[it1].EndChg;
2828 nh1 = tcl[it1].tclhits.size();
2829 pass1 = tcl[it1].ProcCode - 10 * (tcl[it1].ProcCode / 10);
2830 if (pass1 > fNumPass) pass1 = fNumPass;
2831 for (it2 = it1 + 1; it2 < tclsize; ++it2) {
2833 if (tcl[it1].ID < 0)
continue;
2834 if (tcl[it2].ID < 0)
continue;
2836 if (tcl[it2].CTP != clCTP)
continue;
2839 if (!ChkMergedClusterHitFrac(it1, it2)) {
continue; }
2840 bs2 = tcl[it2].BeginSlp;
2841 bth2 = std::atan(fScaleF * bs2);
2842 bw2 = tcl[it2].BeginWir;
2843 bt2 = tcl[it2].BeginTim;
2844 bc2 = tcl[it2].BeginChg;
2845 es2 = tcl[it2].EndSlp;
2846 eth2 = tcl[it2].EndAng;
2847 ew2 = tcl[it2].EndWir;
2848 et2 = tcl[it2].EndTim;
2849 ec2 = tcl[it2].EndChg;
2850 pass2 = tcl[it2].ProcCode - 10 * (tcl[it2].ProcCode / 10);
2851 if (pass2 > fNumPass) pass2 = fNumPass;
2853 angcut = fKinkAngCut[pass1];
2854 if (fKinkAngCut[pass2] > angcut) angcut = fKinkAngCut[pass2];
2855 skipcut = fMaxWirSkip[pass1];
2856 if (fMaxWirSkip[pass2] > skipcut) skipcut = fMaxWirSkip[pass2];
2857 chgcut = fMergeChgCut[pass1];
2858 if (fMergeChgCut[pass2] > chgcut) chgcut = fMergeChgCut[pass2];
2860 bothLong = (nh1 > 100 && tcl[it2].tclhits.size() > 100);
2868 dw = ew1 - bw2 - ndead;
2870 NoVtx = (tcl[it1].EndVtx < 0) && (tcl[it2].BeginVtx < 0);
2871 if (prt && bw2 < (ew1 + maxOverlap))
2872 mf::LogVerbatim(
"CC") <<
"Chk1 ID1-2 " << tcl[it1].ID <<
"-" << tcl[it2].ID <<
" " << ew1
2873 <<
":" << (int)et1 <<
" " << bw2 <<
":" << (
int)bt2 <<
" dw " << dw
2874 <<
" ndead " << ndead <<
" skipcut " << skipcut <<
" dth " << dth
2875 <<
" angcut " << angcut;
2876 if (NoVtx && bw2 < (ew1 + maxOverlap) && dw < skipcut && dth < angcut) {
2877 chgrat = 2 * fabs(bc2 - ec1) / (bc2 + ec1);
2879 if (bothLong && dth < 0.05) chgrat = 0.;
2881 dtim = fabs(bt2 + (ew1 - bw2) * bs2 - et1);
2884 timecut = fTimeDelta[pass2] * AngleFactor(bigslp);
2886 mf::LogVerbatim(
"CC") <<
" dtim " << dtim <<
" timecut " << (int)timecut <<
" ec1 " 2887 << ec1 <<
" bc2 " << bc2 <<
" chgrat " << chgrat <<
" chgcut " 2888 << chgcut <<
" es1 " << es1 <<
" ChkSignal " 2889 << ChkSignal(ew1, et1, bw2, bt2);
2890 if (chgrat < chgcut && dtim < timecut) {
2892 if (ChkSignal(ew1, et1, bw2, bt2)) {
2893 DoMerge(it2, it1, 10);
2894 tclsize = tcl.size();
2902 dth = fabs(bth1 - eth2);
2904 dw = ew2 - bw1 - ndead;
2905 if (prt && bw1 < (ew2 + maxOverlap) && dw < skipcut)
2906 mf::LogVerbatim(
"CC") <<
"Chk2 ID1-2 " << tcl[it1].ID <<
"-" << tcl[it2].ID <<
" " << bw1
2907 <<
":" << (int)bt1 <<
" " << bw2 <<
":" << (
int)et2 <<
" dw " << dw
2908 <<
" ndead " << ndead <<
" skipcut " << skipcut <<
" dth " << dth
2909 <<
" angcut " << angcut;
2911 NoVtx = (tcl[it2].EndVtx < 0) && (tcl[it1].BeginVtx < 0);
2912 if (NoVtx && bw1 < (ew2 + maxOverlap) && dw < skipcut && dth < angcut) {
2913 chgrat = 2 * fabs((bc1 - ec2) / (bc1 + ec2));
2915 if (bothLong && dth < 0.05) chgrat = 0.;
2917 dtim =
std::abs(bt1 + (ew2 - bw1) * bs1 - et2);
2920 timecut = fTimeDelta[pass2] * AngleFactor(bigslp);
2922 mf::LogVerbatim(
"CC") <<
" dtim " << dtim <<
" err " << dtim <<
" timecut " 2923 << (int)timecut <<
" chgrat " << chgrat <<
" chgcut " << chgcut
2924 <<
" ChkSignal " << ChkSignal(bw1, bt1, ew2, et2);
2926 if (chgrat < chgcut && dtim < timecut) {
2927 if (ChkSignal(bw1, bt1, ew2, et2)) {
2928 DoMerge(it1, it2, 10);
2929 tclsize = tcl.size();
2935 if (bw2 < bw1 && ew2 > ew1) {
2938 dth = fabs(eth2 - eth1);
2939 dtim = fabs(et1 + (ew2 - ew1 - 1) * es1 - et2);
2942 timecut = fTimeDelta[pass2] * AngleFactor(bigslp);
2944 short nmiss1 = bw1 - ew1 + 1 - tcl[it1].tclhits.size();
2946 short nin2 = tcl[it2].tclhits.size();
2948 mf::LogVerbatim(
"CC") <<
"cl2: " << ew2 <<
":" << (int)et2 <<
" - " << bw2 <<
":" 2949 << (
int)bt2 <<
" within cl1 " << ew1 <<
":" << (int)et1 <<
" - " 2950 << bw1 <<
":" << (
int)bt1 <<
" ? dth " << dth <<
" dtim " << dtim
2951 <<
" nmissed " << nmiss1 <<
" timecut " << timecut
2952 <<
" FIX THIS CODE";
2957 if (dth < 1 && dtim < timecut && nmiss1 >= nin2) ChkMerge12(it1, it2, didit);
2959 tclsize = tcl.size();
2964 if (bw1 < bw2 && ew1 > ew2) {
2968 dtim =
std::abs(et2 + (ew1 - ew2 - 1) * es2 - et1);
2971 timecut = fTimeDelta[pass2] * AngleFactor(bigslp);
2973 short nmiss2 = bw2 - ew2 + 1 - tcl[it2].tclhits.size();
2975 short nin1 = tcl[it1].tclhits.size();
2977 mf::LogVerbatim(
"CC") <<
"cl1: " << ew1 <<
":" << (int)et1 <<
" - " << bw1 <<
":" 2978 << (
int)bt1 <<
" within cl2 " << ew2 <<
":" << (int)et2 <<
" - " 2979 << bw2 <<
":" << (
int)bt2 <<
" ? dth " << dth <<
" dtim " << dtim
2980 <<
" nmissed " << nmiss2 <<
" timecut " << timecut
2981 <<
" FIX THIS CODE";
2985 if (dth < 1 && dtim < timecut && nmiss2 >= nin1) ChkMerge12(it2, it1, didit);
2987 tclsize = tcl.size();
2992 if (tcl[it1].ID < 0)
break;
2994 if (tcl[it1].ID < 0)
continue;
2999 void ClusterCrawlerAlg::ChkMerge12(
unsigned short it1,
unsigned short it2,
bool& didit)
3009 if (prt)
mf::LogVerbatim(
"CC") <<
"ChkMerge12 " << tcl[it1].ID <<
" " << tcl[it2].ID;
3014 int ew1 = tcl[it1].
EndWir;
3015 int bw1 = tcl[it1].BeginWir;
3016 unsigned int iht,
hit;
3018 std::vector<unsigned int> cl1hits(bw1 + 1 - ew1);
3020 for (iht = 0; iht < cl1.
tclhits.size(); ++iht) {
3022 wire = fHits[hit].WireID().Wire;
3023 if (wire - ew1 < 0 || wire - ew1 > (
int)cl1hits.size()) {
return; }
3024 cl1hits[wire - ew1] = hit;
3026 unsigned int ew2 = tcl[it2].EndWir;
3027 float et2 = tcl[it2].EndTim;
3029 unsigned int wiron1 = 0;
3032 for (wire = ew2 - 1; wire > ew1; --wire) {
3033 if (cl1hits[wire - ew1] > 0) {
3039 if (prt)
mf::LogVerbatim(
"CC") <<
"chk next US wire " << wiron1 <<
" missed " << nmiss;
3040 if (wiron1 == 0)
return;
3041 if (nmiss > fMaxWirSkip[pass])
return;
3045 unsigned int hiton2;
3047 unsigned short nfit = 0;
3048 for (iht = 0; iht < tcl[it2].tclhits.size(); ++iht) {
3049 hiton2 = tcl[it2].tclhits[iht];
3050 wiron2 = fHits[hiton2].WireID().Wire;
3051 if (wiron2 < ew1 || wiron2 > bw1)
return;
3052 if (cl1hits[wiron2 - ew1] == 0) ++nfit;
3055 if (nfit < tcl[it2].tclhits.size())
return;
3058 unsigned short pass1 = tcl[it1].ProcCode - 10 * (tcl[it1].ProcCode / 10);
3059 unsigned short pass2 = tcl[it2].ProcCode - 10 * (tcl[it2].ProcCode / 10);
3060 unsigned short cpass = pass1;
3062 if (pass2 < pass1) cpass = pass2;
3065 if ((
int)wiron1 - ew1 < 0)
return;
3066 unsigned int hiton1 = cl1hits[wiron1 - ew1];
3067 if (hiton1 > fHits.size() - 1) {
return; }
3069 float timon1 = fHits[hiton1].PeakTime();
3070 float dtim =
std::abs(et2 + (wiron1 - ew2) * tcl[it2].EndSlp - timon1);
3071 if (dtim > fTimeDelta[cpass])
return;
3074 FitClusterMid(it1, hiton1, 3);
3075 if (clChisq > 20.)
return;
3078 float dth =
std::abs(std::atan(fScaleF * clpar[1]) - std::atan(fScaleF * tcl[it2].EndSlp));
3079 if (prt)
mf::LogVerbatim(
"CC") <<
"US dtheta " << dth <<
" cut " << fKinkAngCut[cpass];
3080 if (dth > fKinkAngCut[cpass])
return;
3082 float chgrat = 2 *
std::abs(fAveChg - tcl[it2].EndChg) / (fAveChg + tcl[it2].EndChg);
3083 if (prt)
mf::LogVerbatim(
"CC") <<
"US chgrat " << chgrat <<
" cut " << fMergeChgCut[pass];
3086 SigOK = ChkSignal(wiron1, timon1, ew2, et2);
3091 unsigned int bw2 = tcl[it2].BeginWir;
3092 float bt2 = tcl[it2].BeginTim;
3095 for (wire = bw2 + 1; wire < bw1; ++wire) {
3096 if (cl1hits[wire - ew1] > 0) {
3102 if (wiron1 == 0)
return;
3103 if (nmiss > fMaxWirSkip[pass])
return;
3106 hiton1 = cl1hits[wiron1 - ew1];
3107 if (hiton1 > fHits.size() - 1) {
return; }
3108 timon1 = fHits[hiton1].PeakTime();
3109 dtim =
std::abs(bt2 - (wiron1 - bw2) * tcl[it2].BeginSlp - timon1);
3110 if (dtim > fTimeDelta[cpass])
return;
3111 FitClusterMid(it1, hiton1, -3);
3112 if (clChisq > 20.)
return;
3114 dth =
std::abs(std::atan(fScaleF * clpar[1]) - std::atan(fScaleF * tcl[it2].BeginSlp));
3115 if (prt)
mf::LogVerbatim(
"CC") <<
"DS dtheta " << dth <<
" cut " << fKinkAngCut[cpass];
3116 if (dth > fKinkAngCut[cpass])
return;
3118 chgrat = 2 *
std::abs(fAveChg - tcl[it2].BeginChg) / (fAveChg + tcl[it2].BeginChg);
3119 if (prt)
mf::LogVerbatim(
"CC") <<
"DS chgrat " << chgrat <<
" cut " << fMergeChgCut[pass];
3121 SigOK = ChkSignal(wiron1, timon1, bw2, bt2);
3127 DoMerge(it1, it2, 100);
3132 void ClusterCrawlerAlg::DoMerge(
unsigned short it1,
unsigned short it2,
short inProcCode)
3139 if (cl1.
tclhits.size() < 3)
return;
3140 if (cl2.
tclhits.size() < 3)
return;
3142 unsigned int lowire, hiwire, whsize, ii, iht, indx;
3145 unsigned int fithit;
3160 whsize = hiwire + 2 - lowire;
3161 std::vector<int> wirehit(whsize, -1);
3163 for (ii = 0; ii < cl2.
tclhits.size(); ++ii) {
3165 indx = fHits[iht].WireID().Wire - lowire;
3166 wirehit[indx] = iht;
3169 for (ii = 0; ii < cl1.
tclhits.size(); ++ii) {
3171 indx = fHits[iht].WireID().Wire - lowire;
3172 wirehit[indx] = iht;
3176 for (std::vector<int>::reverse_iterator rit = wirehit.rbegin(); rit != wirehit.rend(); ++rit) {
3177 if (*rit < 0)
continue;
3178 fcl2hits.push_back(*rit);
3183 FitClusterMid(fcl2hits, fithit, nhitfit);
3184 if (clChisq > 5)
return;
3187 MakeClusterObsolete(it1);
3188 MakeClusterObsolete(it2);
3192 short del1Vtx = -99;
3193 short del2Vtx = -99;
3242 if (!TmpStore())
return;
3243 unsigned short itnew = tcl.size() - 1;
3245 mf::LogVerbatim(
"CC") <<
"DoMerge " << cl1.
ID <<
" " << cl2.
ID <<
" -> " << tcl[itnew].ID;
3247 tcl[itnew].ProcCode = inProcCode + pass;
3250 if (del1Vtx >= 0 && del1Vtx == del2Vtx) vtx[del1Vtx].NClusters = 0;
3252 tcl[itnew].BeginVtx = begVtx;
3253 tcl[itnew].EndVtx = endVtx;
3257 void ClusterCrawlerAlg::PrintVertices()
3262 if (vtx3.size() > 0) {
3265 <<
"****** 3D vertices ******************************************__2DVtx_Indx__*******\n";
3267 <<
"Vtx Cstat TPC Proc X Y Z XEr YEr ZEr pln0 pln1 pln2 Wire\n";
3268 for (
unsigned short iv = 0; iv < vtx3.size(); ++iv) {
3269 myprt <<
std::right << std::setw(3) << std::fixed << iv << std::setprecision(1);
3270 myprt <<
std::right << std::setw(7) << vtx3[iv].CStat;
3271 myprt <<
std::right << std::setw(5) << vtx3[iv].TPC;
3272 myprt <<
std::right << std::setw(5) << vtx3[iv].ProcCode;
3273 myprt <<
std::right << std::setw(8) << vtx3[iv].X;
3274 myprt <<
std::right << std::setw(8) << vtx3[iv].Y;
3275 myprt <<
std::right << std::setw(8) << vtx3[iv].Z;
3276 myprt <<
std::right << std::setw(5) << vtx3[iv].XErr;
3277 myprt <<
std::right << std::setw(5) << vtx3[iv].YErr;
3278 myprt <<
std::right << std::setw(5) << vtx3[iv].ZErr;
3279 myprt <<
std::right << std::setw(5) << vtx3[iv].Ptr2D[0];
3280 myprt <<
std::right << std::setw(5) << vtx3[iv].Ptr2D[1];
3281 myprt <<
std::right << std::setw(5) << vtx3[iv].Ptr2D[2];
3282 myprt <<
std::right << std::setw(5) << vtx3[iv].Wire;
3283 if (vtx3[iv].Wire < 0) { myprt <<
" Matched in all planes"; }
3285 myprt <<
" Incomplete";
3291 if (vtx.size() > 0) {
3293 myprt <<
"************ 2D vertices ************\n";
3294 myprt <<
"Vtx CTP wire error tick error ChiDOF NCl topo cluster IDs\n";
3295 for (
unsigned short iv = 0; iv < vtx.size(); ++iv) {
3296 if (fDebugPlane < 3 && fDebugPlane != (
int)vtx[iv].CTP)
continue;
3297 myprt <<
std::right << std::setw(3) << std::fixed << iv << std::setprecision(1);
3298 myprt <<
std::right << std::setw(6) << vtx[iv].CTP;
3299 myprt <<
std::right << std::setw(8) << vtx[iv].Wire <<
" +/- ";
3300 myprt <<
std::right << std::setw(4) << vtx[iv].WireErr;
3301 myprt <<
std::right << std::setw(8) << vtx[iv].Time <<
" +/- ";
3302 myprt <<
std::right << std::setw(4) << vtx[iv].TimeErr;
3303 myprt <<
std::right << std::setw(8) << vtx[iv].ChiDOF;
3304 myprt <<
std::right << std::setw(5) << vtx[iv].NClusters;
3305 myprt <<
std::right << std::setw(6) << vtx[iv].Topo;
3308 for (
unsigned short ii = 0; ii < tcl.size(); ++ii) {
3309 if (fDebugPlane < 3 && fDebugPlane != (
int)tcl[ii].CTP)
continue;
3310 if (tcl[ii].ID < 0)
continue;
3311 if (tcl[ii].BeginVtx == (
short)iv) myprt <<
std::right << std::setw(4) << tcl[ii].ID;
3312 if (tcl[ii].EndVtx == (
short)iv) myprt <<
std::right << std::setw(4) << tcl[ii].ID;
3329 float aveRMS, aveRes;
3330 myprt <<
"*************************************** Clusters " 3331 "*********************************************************************\n";
3332 myprt <<
" ID CTP nht Stop Proc beg_W:T bAng bSlp Err bChg end_W:T eAng eSlp " 3333 "Err eChg bVx eVx aveRMS Qual cnt\n";
3334 for (
unsigned short ii = 0; ii < tcl.size(); ++ii) {
3336 if (fDebugPlane < 3 && fDebugPlane != (
int)tcl[ii].CTP)
continue;
3337 myprt <<
std::right << std::setw(4) << tcl[ii].ID;
3338 myprt <<
std::right << std::setw(3) << tcl[ii].CTP;
3339 myprt <<
std::right << std::setw(5) << tcl[ii].tclhits.size();
3340 myprt <<
std::right << std::setw(4) << tcl[ii].StopCode;
3341 myprt <<
std::right << std::setw(6) << tcl[ii].ProcCode;
3342 unsigned int iTime = tcl[ii].BeginTim;
3343 myprt <<
std::right << std::setw(6) << tcl[ii].BeginWir <<
":" << iTime;
3344 if (iTime < 10) { myprt <<
" "; }
3345 else if (iTime < 100) {
3348 else if (iTime < 1000)
3350 myprt <<
std::right << std::setw(7) << std::fixed << std::setprecision(2) << tcl[ii].BeginAng;
3351 if (
std::abs(tcl[ii].BeginSlp) < 100) {
3352 myprt <<
std::right << std::setw(6) << std::fixed << std::setprecision(2)
3353 << tcl[ii].BeginSlp;
3356 myprt <<
std::right << std::setw(6) << (int)tcl[ii].BeginSlp;
3358 myprt <<
std::right << std::setw(6) << std::fixed << std::setprecision(2)
3359 << tcl[ii].BeginSlpErr;
3360 myprt <<
std::right << std::setw(5) << (int)tcl[ii].BeginChg;
3361 iTime = tcl[ii].EndTim;
3362 myprt <<
std::right << std::setw(6) << tcl[ii].EndWir <<
":" << iTime;
3363 if (iTime < 10) { myprt <<
" "; }
3364 else if (iTime < 100) {
3367 else if (iTime < 1000)
3369 myprt <<
std::right << std::setw(7) << std::fixed << std::setprecision(2) << tcl[ii].EndAng;
3370 if (
std::abs(tcl[ii].EndSlp) < 100) {
3371 myprt <<
std::right << std::setw(6) << std::fixed << std::setprecision(2) << tcl[ii].EndSlp;
3374 myprt <<
std::right << std::setw(6) << (int)tcl[ii].EndSlp;
3376 myprt <<
std::right << std::setw(6) << std::fixed << std::setprecision(2)
3377 << tcl[ii].EndSlpErr;
3378 myprt <<
std::right << std::setw(5) << (int)tcl[ii].EndChg;
3379 myprt <<
std::right << std::setw(5) << tcl[ii].BeginVtx;
3380 myprt <<
std::right << std::setw(5) << tcl[ii].EndVtx;
3382 unsigned int iht = 0;
3383 for (
unsigned short jj = 0; jj < tcl[ii].tclhits.size(); ++jj) {
3384 iht = tcl[ii].tclhits[jj];
3385 aveRMS += fHits[iht].RMS();
3387 aveRMS /= (float)tcl[ii].tclhits.size();
3388 myprt <<
std::right << std::setw(5) << std::fixed << std::setprecision(1) << aveRMS;
3391 unsigned int hit0, hit1, hit2, cnt = 0;
3393 for (
unsigned short iht = 1; iht < tcl[ii].tclhits.size() - 1; ++iht) {
3394 hit1 = tcl[ii].tclhits[iht];
3395 hit0 = tcl[ii].tclhits[iht - 1];
3396 hit2 = tcl[ii].tclhits[iht + 1];
3398 if (fHits[hit1].
WireID().Wire + 1 != fHits[hit0].WireID().Wire)
continue;
3399 if (fHits[hit2].
WireID().Wire + 1 != fHits[hit1].WireID().Wire)
continue;
3400 arg = (fHits[hit0].PeakTime() + fHits[hit2].PeakTime()) / 2 - fHits[hit1].PeakTime();
3401 aveRes += arg * arg;
3405 aveRes /= (float)cnt;
3406 aveRes = sqrt(aveRes);
3408 aveRes /= (aveRMS * fHitErrFac);
3409 myprt <<
std::right << std::setw(6) << std::fixed << std::setprecision(1) << aveRes;
3410 myprt <<
std::right << std::setw(5) << std::fixed << cnt;
3414 myprt <<
std::right << std::setw(5) << std::fixed << cnt;
3422 void ClusterCrawlerAlg::TmpGet(
unsigned short it1)
3427 if (it1 > tcl.size())
return;
3429 clBeginSlp = tcl[it1].BeginSlp;
3430 clBeginSlpErr = tcl[it1].BeginSlpErr;
3431 clBeginAng = tcl[it1].BeginAng;
3432 clBeginWir = tcl[it1].BeginWir;
3433 clBeginTim = tcl[it1].BeginTim;
3434 clBeginChg = tcl[it1].BeginChg;
3435 clBeginChgNear = tcl[it1].BeginChgNear;
3436 clEndSlp = tcl[it1].EndSlp;
3437 clEndSlpErr = tcl[it1].EndSlpErr;
3438 clEndAng = tcl[it1].EndAng;
3439 clEndWir = tcl[it1].EndWir;
3440 clEndTim = tcl[it1].EndTim;
3441 clEndChg = tcl[it1].EndChg;
3442 clEndChgNear = tcl[it1].EndChgNear;
3443 clStopCode = tcl[it1].StopCode;
3444 clProcCode = tcl[it1].ProcCode;
3445 clCTP = tcl[it1].CTP;
3446 fcl2hits = tcl[it1].tclhits;
3450 bool ClusterCrawlerAlg::TmpStore()
3453 if (fcl2hits.size() < 2)
return false;
3454 if (fcl2hits.size() > fHits.size())
return false;
3456 if (NClusters == SHRT_MAX)
return false;
3460 unsigned int hit0 = fcl2hits[0];
3462 unsigned int tCST = fHits[hit0].WireID().Cryostat;
3463 unsigned int tTPC = fHits[hit0].WireID().TPC;
3464 unsigned int tPlane = fHits[hit0].WireID().Plane;
3465 unsigned int lastWire = 0;
3467 for (
unsigned short it = 0; it < fcl2hits.size(); ++it) {
3469 if (inClus[hit] != 0) {
3474 if (fHits[hit].
WireID().Cryostat != tCST || fHits[hit].WireID().TPC != tTPC ||
3475 fHits[hit].WireID().Plane != tPlane) {
3480 if (clProcCode < 900 && it > 0 && fHits[hit].
WireID().Wire >= lastWire) {
3484 lastWire = fHits[hit].WireID().Wire;
3485 inClus[hit] = NClusters;
3491 if (clEndChg <= 0) {
3493 unsigned int ih0 = fcl2hits.size() - 2;
3494 hit = fcl2hits[ih0];
3495 clEndChg = fHits[hit].Integral();
3496 hit = fcl2hits[ih0 - 1];
3497 clEndChg += fHits[hit].Integral();
3498 clEndChg = clEndChg / 2.;
3500 if (clBeginChg <= 0) {
3503 clBeginChg = fHits[hit].Integral();
3505 clBeginChg += fHits[hit].Integral();
3506 clBeginChg = clBeginChg / 2.;
3510 hit = fcl2hits[fcl2hits.size() - 1];
3515 clstr.
ID = NClusters;
3518 clstr.
BeginAng = std::atan(fScaleF * clBeginSlp);
3519 clstr.
BeginWir = fHits[hit0].WireID().Wire;
3520 clstr.
BeginTim = fHits[hit0].PeakTime();
3525 clstr.
EndAng = std::atan(fScaleF * clEndSlp);
3526 clstr.
EndWir = fHits[hit].WireID().Wire;
3527 clstr.
EndTim = fHits[hit].PeakTime();
3536 tcl.push_back(clstr);
3542 void ClusterCrawlerAlg::LACrawlUS()
3547 unsigned int dhit = fcl2hits[0];
3548 short dwir = fHits[dhit].WireID().Wire;
3551 if (fDebugPlane == (
short)plane && dwir == fDebugWire && fDebugHit > 0)
3552 prt =
std::abs(fHits[dhit].PeakTime() - fDebugHit) < 40;
3556 myprt <<
"******************* LACrawlUS PASS " << pass <<
" Hits ";
3557 for (
unsigned short ii = 0; ii < fcl2hits.size(); ++ii) {
3558 unsigned int iht = fcl2hits[fcl2hits.size() - 1 - ii];
3559 myprt << fHits[iht].WireID().Wire <<
":" << (int)fHits[iht].PeakTime() <<
" ";
3568 unsigned short kinkOnWire = 0;
3569 unsigned int it = fcl2hits.size() - 1;
3570 unsigned int lasthit = fcl2hits[it];
3571 unsigned int lastwire = fHits[lasthit].WireID().Wire;
3572 unsigned int prevHit, prevWire;
3573 bool ChkCharge =
false;
3574 for (
unsigned int nextwire = lastwire - 1; nextwire >= fFirstWire; --nextwire) {
3576 mf::LogVerbatim(
"CC") <<
"LACrawlUS: next wire " << nextwire <<
" HitRange " 3577 << WireHitRange[nextwire].first;
3579 if (CrawlVtxChk(nextwire)) {
3585 AddLAHit(nextwire, ChkCharge, HitOK, SigOK);
3587 mf::LogVerbatim(
"CC") <<
"LACrawlUS: HitOK " << HitOK <<
" SigOK " << SigOK
3588 <<
" nmissed SigOK " << nmissed <<
" cut " << fAllowNoHitWire;
3589 if (SigOK) nmissed = 0;
3592 if (nmissed > fAllowNoHitWire) {
3601 prevHit = fcl2hits[fcl2hits.size() - 2];
3602 prevWire = fHits[prevHit].WireID().Wire;
3603 if (prevWire > nextwire + 2) {
3604 if (!ChkSignal(fcl2hits[fcl2hits.size() - 1], fcl2hits[fcl2hits.size() - 2])) {
3614 if (fcl2hits.size() == 4) {
3616 for (
unsigned short kk = 0; kk < fcl2hits.size() - 1; ++kk) {
3617 unsigned int hit = fcl2hits[kk];
3618 if (mergeAvailable[hit]) MergeHits(hit, didMerge);
3622 clBeginSlp = clpar[1];
3627 unsigned short chsiz = chifits.size() - 1;
3629 if (chsiz < 6)
continue;
3630 if (fKinkChiRat[pass] <= 0)
continue;
3631 if (chifits.size() != fcl2hits.size()) {
3632 mf::LogError(
"CC") <<
"LACrawlUS: chifits size error " << chifits.size() <<
" " 3637 mf::LogVerbatim(
"CC") <<
"Kink chk " << chifits[chsiz] <<
" " << chifits[chsiz - 1] <<
" " 3638 << chifits[chsiz - 2] <<
" " << chifits[chsiz - 3];
3639 if (chifits[chsiz - 1] > fKinkChiRat[pass] * chifits[chsiz - 2] &&
3640 chifits[chsiz] > fKinkChiRat[pass] * chifits[chsiz - 1]) {
3642 unsigned int ih0 = fcl2hits.size() - 1;
3643 unsigned int hit0 = fcl2hits[ih0];
3644 unsigned int ih2 = ih0 - 2;
3645 unsigned int hit2 = fcl2hits[ih2];
3646 float dt02 = fHits[hit2].PeakTime() - fHits[hit0].PeakTime();
3647 float dw02 = fHits[hit2].WireID().Wire - fHits[hit0].WireID().Wire;
3648 float th02 = std::atan(fScaleF * dt02 / dw02);
3650 unsigned int ih3 = ih2 - 1;
3651 unsigned int hit3 = fcl2hits[ih3];
3652 unsigned int ih5 = ih3 - 2;
3653 unsigned int hit5 = fcl2hits[ih5];
3654 float dt35 = fHits[hit5].PeakTime() - fHits[hit3].PeakTime();
3655 float dw35 = fHits[hit5].WireID().Wire - fHits[hit3].WireID().Wire;
3656 float th35 = std::atan(fScaleF * dt35 / dw35);
3659 mf::LogVerbatim(
"CC") <<
" Kink angle " << std::setprecision(3) << dth <<
" cut " 3660 << fKinkAngCut[pass];
3661 if (dth > fKinkAngCut[pass]) {
3667 if (kinkOnWire > 0) {
3668 if (kinkOnWire - nextwire < 4) {
3671 <<
"Hit a second kink. kinkOnWire = " << kinkOnWire <<
" Stopping";
3677 kinkOnWire = nextwire;
3683 CheckClusterHitFrac(prt);
3686 if (prt)
mf::LogVerbatim(
"CC") <<
"LACrawlUS done. Nhits = " << fcl2hits.size();
3691 void ClusterCrawlerAlg::CrawlUS()
3695 if (fcl2hits.size() < 2)
return;
3697 unsigned int dhit = fcl2hits[0];
3698 int dwir = fHits[dhit].WireID().Wire;
3702 if (fDebugPlane == (
short)plane && dwir == fDebugWire && fDebugHit > 0)
3703 prt =
std::abs(fHits[dhit].PeakTime() - fDebugHit) < 20;
3707 myprt <<
"******************* Start CrawlUS on pass " << pass <<
" Hits: ";
3708 for (
unsigned short ii = 0; ii < fcl2hits.size(); ++ii) {
3709 unsigned int iht = fcl2hits[fcl2hits.size() - 1 - ii];
3710 myprt << fHits[iht].WireID().Wire <<
":" << (int)fHits[iht].PeakTime() <<
" ";
3721 short nHitAfterSkip = 0;
3722 bool DidaSkip =
false;
3723 bool PostSkip =
false;
3724 unsigned int it = fcl2hits.size() - 1;
3725 unsigned int lasthit = fcl2hits[it];
3726 if (lasthit > fHits.size() - 1) {
3727 mf::LogError(
"CC") <<
"CrawlUS bad lasthit " << lasthit;
3730 unsigned int lastwire = fHits[lasthit].WireID().Wire;
3731 if (prt)
mf::LogVerbatim(
"CC") <<
"CrawlUS: last wire " << lastwire <<
" hit " << lasthit;
3733 unsigned int lastWireWithSignal = lastwire;
3734 unsigned short nDroppedHit = 0;
3736 for (
unsigned int nextwire = lastwire - 1; nextwire > fFirstWire; --nextwire) {
3738 mf::LogVerbatim(
"CC") <<
"CrawlUS: next wire " << nextwire <<
" HitRange " 3739 << WireHitRange[nextwire].first;
3741 if (CrawlVtxChk(nextwire)) {
3747 if (
std::abs(clpar[1]) > fLAClusSlopeCut) {
3748 if (prt)
mf::LogVerbatim(
"CC") <<
">>>>> CrawlUS: Switching to LACrawlUS";
3753 AddHit(nextwire, HitOK, SigOK);
3755 mf::LogVerbatim(
"CC") <<
"CrawlUS: HitOK " << HitOK <<
" SigOK " << SigOK <<
" nmissed " 3757 if (SigOK) lastWireWithSignal = nextwire;
3764 if (PostSkip && nmissed > fMinWirAfterSkip[pass]) {
3766 if ((
int)(fcl2hits.size() - nHitAfterSkip) < 4) {
3770 if (prt)
mf::LogVerbatim(
"CC") <<
" PostSkip && nmissed = " << nmissed;
3772 FclTrimUS(nHitAfterSkip);
3774 if (clChisq > 90.) {
3789 lasthit = fcl2hits[fcl2hits.size() - 1];
3790 if ((lastWireWithSignal - nextwire) > fAllowNoHitWire) {
3793 <<
"No hit or signal on wire " << nextwire <<
" last wire with signal " 3794 << lastWireWithSignal <<
" exceeding fAllowNoHitWire " << fAllowNoHitWire
3802 mf::LogVerbatim(
"CC") <<
" CrawlUS check clChisq " << clChisq <<
" cut " << fChiCut[pass];
3803 if (clChisq > fChiCut[pass]) {
3804 if (fcl2hits.size() < 3) {
3809 if (prt)
mf::LogVerbatim(
"CC") <<
"ADD- Bad clChisq " << clChisq <<
" dropping hit";
3813 if (nDroppedHit > 4) {
3816 <<
" Too many dropped hits: " << nDroppedHit <<
" Stop crawling";
3819 if (clChisq > fChiCut[pass]) {
3822 <<
" Bad clChisq " << clChisq <<
" after dropping hit. Stop crawling";
3830 if (chifits.size() > 5 && fKinkChiRat[pass] > 0) {
3831 if (chifits.size() != fcl2hits.size()) {
3832 mf::LogError(
"CC") <<
"CrawlUS: chifits size error " << chifits.size() <<
" " 3836 unsigned short chsiz = chifits.size() - 1;
3838 mf::LogVerbatim(
"CC") <<
"Kink chk " << chifits[chsiz] <<
" " << chifits[chsiz - 1]
3839 <<
" " << chifits[chsiz - 2] <<
" " << chifits[chsiz - 3];
3840 if (chifits[chsiz - 1] > fKinkChiRat[pass] * chifits[chsiz - 2] &&
3841 chifits[chsiz] > fKinkChiRat[pass] * chifits[chsiz - 1]) {
3842 if (fcl2hits.size() != chifits.size()) {
3843 mf::LogError(
"CC") <<
"bad kink check size " << chifits.size() <<
" " 3844 << fcl2hits.size() <<
" plane " << plane <<
" cluster " << dwir
3848 if (EndKinkAngle() > fKinkAngCut[pass]) {
3851 <<
"******************* Stopped tracking - kink angle " << EndKinkAngle() <<
" > " 3852 << fKinkAngCut[pass] <<
" Removing 3 hits";
3862 if (fcl2hits.size() == fMaxHitsFit[pass] || fcl2hits.size() == fMinHits[pass]) {
3863 clBeginSlp = clpar[1];
3864 clBeginSlpErr = clparerr[1];
3867 if (clBeginChg <= 0 && fAveChg > 0) {
3868 clBeginChg = fAveChg;
3884 if (nHitAfterSkip == fMinWirAfterSkip[pass]) PostSkip =
false;
3887 if (clChisq > fChiCut[pass]) {
3891 unsigned short lopped = 0;
3892 for (
unsigned short nlop = 0; nlop < 4; ++nlop) {
3893 unsigned short cfsize = chifits.size() - 1;
3894 chirat = chifits[cfsize] / chifits[cfsize - 1];
3897 <<
"chirat " << chirat <<
" last hit " << fcl2hits[fcl2hits.size() - 1];
3898 if (chirat < 1.2)
break;
3899 if (prt)
mf::LogVerbatim(
"CC") <<
"<<ADD- Bad chisq. Bad chirat " << chirat;
3902 if (fcl2hits.size() < 6)
break;
3903 if (chifits.size() < 6)
break;
3905 if (fcl2hits.size() < 6) {
3907 if (prt)
mf::LogVerbatim(
"CC") <<
"Bad fit chisq - short cluster. Break";
3910 if (lopped == 0 && fcl2hits.size() > 5) {
3918 mf::LogVerbatim(
"CC") <<
"Bad fit chisq - removed " << lopped <<
" hits. Continue";
3923 mf::LogVerbatim(
"CC") <<
"******************* CrawlUS normal stop. size " << fcl2hits.size();
3927 if (fcl2hits.size() > 5) {
3930 mf::LogVerbatim(
"CC") <<
"EndKinkAngle check " << EndKinkAngle() <<
" cut " 3931 << fKinkAngCut[pass];
3932 if (EndKinkAngle() > fKinkAngCut[pass]) {
3941 if ((
unsigned short)fcl2hits.size() > fMinWirAfterSkip[pass] + 3) {
3942 unsigned int ih0 = fcl2hits.size() - 1;
3943 unsigned int hit0 = fcl2hits[ih0];
3944 unsigned int uswir = fHits[hit0].WireID().Wire;
3945 unsigned int nxtwir;
3946 unsigned short nAdjHit = 0;
3947 for (
unsigned short ii = ih0 - 1; ii > 0; --ii) {
3948 nxtwir = fHits[fcl2hits[ii]].WireID().Wire;
3950 if (nxtwir != uswir + 1)
break;
3952 if (nAdjHit == fMinWirAfterSkip[pass])
break;
3956 if (nAdjHit < fMinWirAfterSkip[pass]) {
3957 if (prt)
mf::LogVerbatim(
"CC") <<
"fMinWirAfterSkip removes " << nAdjHit <<
" hits ";
3964 if (!reFit && fcl2hits.size() > 3) {
3965 float chirat = chifits[chifits.size() - 1] / chifits[chifits.size() - 2];
3967 mf::LogVerbatim(
"CC") <<
"Last hit chirat " << chirat <<
" cut " << fKinkChiRat[pass];
3969 mf::LogVerbatim(
"CC") <<
"Check " << clChisq <<
" " << chifits[chifits.size() - 1] <<
" " 3970 << chifits[chifits.size() - 2];
3971 if (chirat > fKinkChiRat[pass]) {
3982 CheckClusterHitFrac(prt);
3984 mf::LogVerbatim(
"CC") <<
"******************* CrawlUS done. Size " << fcl2hits.size()
3985 <<
" min length for this pass " << fMinHits[pass];
3991 float ClusterCrawlerAlg::EndKinkAngle()
3995 unsigned int ih0 = fcl2hits.size() - 1;
3996 unsigned int hit0 = fcl2hits[ih0];
3997 unsigned int ih2 = ih0 - 2;
3998 unsigned int hit2 = fcl2hits[ih2];
3999 float dt02 = fHits[hit2].PeakTime() - fHits[hit0].PeakTime();
4000 float dw02 = fHits[hit2].WireID().Wire - fHits[hit0].WireID().Wire;
4001 float th02 = std::atan(fScaleF * dt02 / dw02);
4003 unsigned int ih3 = ih2 - 1;
4004 unsigned int hit3 = fcl2hits[ih3];
4005 unsigned int ih5 = ih3 - 2;
4006 unsigned int hit5 = fcl2hits[ih5];
4007 float dt35 = fHits[hit5].PeakTime() - fHits[hit3].PeakTime();
4008 float dw35 = fHits[hit5].WireID().Wire - fHits[hit3].WireID().Wire;
4009 float th35 = std::atan(fScaleF * dt35 / dw35);
4014 bool ClusterCrawlerAlg::ChkMergedClusterHitFrac(
unsigned short it1,
unsigned short it2)
4018 if (it1 > tcl.size() - 1)
return false;
4019 if (it2 > tcl.size() - 1)
return false;
4020 unsigned int eWire = 99999;
4021 unsigned int bWire = 0, wire;
4022 if (tcl[it1].BeginWir > bWire) bWire = tcl[it1].BeginWir;
4023 if (tcl[it2].BeginWir > bWire) bWire = tcl[it2].BeginWir;
4024 if (tcl[it1].EndWir < eWire) eWire = tcl[it1].EndWir;
4025 if (tcl[it2].EndWir < eWire) eWire = tcl[it2].EndWir;
4026 unsigned int mergedClusterLength = bWire - eWire + 1;
4028 std::vector<bool> cHits(mergedClusterLength,
false);
4030 unsigned int ii, iht, indx;
4031 for (ii = 0; ii < tcl[it1].tclhits.size(); ++ii) {
4032 iht = tcl[it1].tclhits[ii];
4033 if (iht > fHits.size() - 1) {
4034 mf::LogError(
"CC") <<
"ChkMergedClusterHitFrac bad iht ";
4037 indx = fHits[iht].WireID().Wire - eWire;
4040 for (ii = 0; ii < tcl[it2].tclhits.size(); ++ii) {
4041 iht = tcl[it2].tclhits[ii];
4042 if (iht > fHits.size() - 1) {
4043 mf::LogError(
"CC") <<
"ChkMergedClusterHitFrac bad iht ";
4046 indx = fHits[iht].WireID().Wire - eWire;
4050 for (ii = 0; ii < cHits.size(); ++ii) {
4052 if (WireHitRange[wire].first == -1) cHits[ii] =
true;
4055 float nhits = std::count(cHits.begin(), cHits.end(),
true);
4056 float hitFrac = nhits / (float)cHits.size();
4057 return (hitFrac > fMinHitFrac);
4061 void ClusterCrawlerAlg::CheckClusterHitFrac(
bool prt)
4066 unsigned int iht = fcl2hits[fcl2hits.size() - 1];
4067 clEndWir = fHits[iht].WireID().Wire;
4068 clBeginWir = fHits[fcl2hits[0]].WireID().Wire;
4069 float hitFrac = (float)(fcl2hits.size() +
DeadWireCount()) / (
float)(clBeginWir - clEndWir + 1);
4071 if (hitFrac < fMinHitFrac) {
4073 mf::LogVerbatim(
"CC") <<
"CheckClusterHitFrac: Poor hit fraction " << hitFrac
4074 <<
" clBeginWir " << clBeginWir <<
" clEndWir " << clEndWir
4075 <<
" size " << fcl2hits.size() <<
" DeadWireCount " 4088 if (fcl2hits.size() < 5) {
4089 unsigned short nsing = 0;
4090 for (
unsigned short iht = 0; iht < fcl2hits.size(); ++iht)
4091 if (fHits[fcl2hits[iht]].Multiplicity() == 1) ++nsing;
4092 hitFrac = (float)nsing / (
float)fcl2hits.size();
4093 if (hitFrac < fMinHitFrac) {
4096 mf::LogVerbatim(
"CC") <<
"CheckClusterHitFrac: Poor short track hit fraction " << hitFrac;
4103 if (clBeginChg <= 0) {
4104 unsigned int iht, nht = 0;
4105 for (
unsigned short ii = 0; ii < fcl2hits.size(); ++ii) {
4107 clBeginChg += fHits[iht].Integral();
4109 if (nht == 8)
break;
4111 clBeginChg /= (float)nht;
4114 unsigned short cnt = chgNear.size() / 2;
4116 if (chgNear.size() > 60) cnt = 30;
4119 for (
unsigned short ids = 0; ids < cnt; ++ids) {
4120 clBeginChgNear += chgNear[ids];
4121 clEndChgNear += chgNear[chgNear.size() - 1 - ids];
4123 clBeginChgNear /= (float)cnt;
4124 clEndChgNear /= (float)cnt;
4129 void ClusterCrawlerAlg::FitClusterMid(
unsigned short it1,
unsigned int ihtin,
short nhit)
4131 FitClusterMid(tcl[it1].tclhits, ihtin, nhit);
4135 void ClusterCrawlerAlg::FitClusterMid(std::vector<unsigned int>& hitVec,
4148 if (hitVec.size() < 3)
return;
4150 std::vector<float> xwir;
4151 std::vector<float> ytim;
4152 std::vector<float> ytimerr2;
4154 unsigned short ii, hitcnt = 0, nht = 0, usnhit;
4164 for (ii = 0; ii < hitVec.size(); ++ii) {
4166 if (iht > fHits.size() - 1) {
4173 wire0 = fHits[iht].WireID().Wire;
4177 xwir.push_back((
float)fHits[iht].
WireID().Wire - wire0);
4178 ytim.push_back(fHits[iht].PeakTime());
4180 float terr = fHitErrFac * fHits[iht].RMS();
4181 ytimerr2.push_back(terr * terr);
4182 fAveChg += fHits[iht].Integral();
4184 if (hitcnt == usnhit)
break;
4192 for (
auto ii = hitVec.crbegin(); ii != hitVec.crend(); ++ii) {
4194 if (iht > fHits.size() - 1) {
4201 wire0 = fHits[iht].WireID().Wire;
4205 xwir.push_back((
float)fHits[iht].
WireID().Wire - wire0);
4206 ytim.push_back(fHits[iht].PeakTime());
4207 float terr = fHitErrFac * fHits[iht].RMS();
4208 ytimerr2.push_back(terr * terr);
4209 fAveChg += fHits[iht].Integral();
4211 if (hitcnt == usnhit)
break;
4217 if (nht < 2)
return;
4218 fAveChg = fAveChg / (float)nht;
4223 float intcpterr = 0.;
4224 float slopeerr = 0.;
4226 fLinFitAlg.LinFit(xwir, ytim, ytimerr2, intcpt, slope, intcpterr, slopeerr, chidof);
4228 if (clChisq > fChiCut[0])
return;
4232 clparerr[0] = intcpterr;
4233 clparerr[1] = slopeerr;
4237 void ClusterCrawlerAlg::FitCluster()
4246 if (pass > fNumPass - 1) {
4247 mf::LogError(
"CC") <<
"FitCluster called on invalid pass " << pass;
4251 unsigned short ii, nht = 0;
4253 nht = fcl2hits.size();
4256 if (nht > fLAClusMaxHitsFit) nht = fLAClusMaxHitsFit;
4259 if (nht > fMaxHitsFit[pass]) nht = fMaxHitsFit[pass];
4261 if (nht < 2)
return;
4263 std::vector<float> xwir;
4264 std::vector<float> ytim;
4265 std::vector<float> ytimerr2;
4267 float angfactor = AngleFactor(clpar[1]);
4270 unsigned int wire0 = fHits[fcl2hits[fcl2hits.size() - 1]].WireID().Wire;
4272 float terr, qave = 0, qwt;
4273 for (ii = 0; ii < nht; ++ii) {
4274 ihit = fcl2hits[fcl2hits.size() - 1 - ii];
4275 qave += fHits[ihit].Integral();
4278 for (ii = 0; ii < nht; ++ii) {
4279 ihit = fcl2hits[fcl2hits.size() - 1 - ii];
4280 wire = fHits[ihit].WireID().Wire;
4281 xwir.push_back(wire - wire0);
4282 ytim.push_back(fHits[ihit].PeakTime());
4285 terr = fHitErrFac * fHits[ihit].RMS() * fHits[ihit].Multiplicity();
4288 qwt = (fHits[ihit].Integral() / qave) - 1;
4289 if (qwt < 1) qwt = 1;
4292 if (terr < 0.01) terr = 0.01;
4293 ytimerr2.push_back(angfactor * terr * terr);
4295 CalculateAveHitWidth();
4298 myprt <<
"FitCluster W:T ";
4299 unsigned short cnt = 0;
4300 for (std::vector<unsigned int>::reverse_iterator it = fcl2hits.rbegin();
4301 it != fcl2hits.rend();
4303 unsigned int ihit = *it;
4304 unsigned short wire = fHits[ihit].WireID().Wire;
4305 myprt << wire <<
":" << (short)fHits[ihit].PeakTime() <<
" ";
4314 if (xwir.size() < 2)
return;
4318 float intcpterr = 0.;
4319 float slopeerr = 0.;
4321 fLinFitAlg.LinFit(xwir, ytim, ytimerr2, intcpt, slope, intcpterr, slopeerr, chidof);
4323 if (chidof > fChiCut[0])
return;
4327 clparerr[0] = intcpterr;
4328 clparerr[1] = slopeerr;
4331 mf::LogVerbatim(
"CC") <<
"nht " << nht <<
" fitpar " << (int)clpar[0] <<
"+/-" 4332 << (
int)intcpterr <<
" " << clpar[1] <<
"+/-" << slopeerr <<
" clChisq " 4336 float ClusterCrawlerAlg::AngleFactor(
float slope)
4342 if (slp > 15) slp = 15;
4344 float angfac = 1 + 0.03 * slp * slp;
4349 void ClusterCrawlerAlg::CalculateAveHitWidth()
4352 for (
unsigned short ii = 0; ii < fcl2hits.size(); ++ii)
4353 fAveHitWidth += fHits[fcl2hits[ii]].EndTick() - fHits[fcl2hits[ii]].StartTick();
4354 fAveHitWidth /= (float)fcl2hits.size();
4358 void ClusterCrawlerAlg::FitClusterChg()
4363 if (fcl2hits.size() == 0)
return;
4364 unsigned int ih0 = fcl2hits.size() - 1;
4366 if (pass >= fNumPass) {
4367 mf::LogError(
"CC") <<
"FitClusterChg bad pass " << pass;
4375 unsigned short nhave = fLAClusMaxHitsFit;
4376 if (nhave > fcl2hits.size()) nhave = fcl2hits.size();
4381 for (
unsigned short ii = 0; ii < nhave; ++ii) {
4382 iht = fcl2hits[fcl2hits.size() - 1 - ii];
4383 fAveChg += fHits[iht].Integral();
4384 fAveHitWidth += (fHits[iht].EndTick() - fHits[iht].StartTick());
4386 fAveChg /= (float)fcl2hits.size();
4387 fAveHitWidth /= (float)fcl2hits.size();
4392 unsigned short fitLen = fNHitsAve[pass];
4396 fcl2hits.size() > 5 &&
4397 fcl2hits.size() < fitLen)
4401 if (fNHitsAve[pass] < 1)
return;
4403 if (fNHitsAve[pass] == 1) {
4405 fAveChg = fHits[fcl2hits[ih0]].Integral();
4408 else if (fNHitsAve[pass] == 2) {
4410 fAveChg = (fHits[fcl2hits[ih0]].Integral() + fHits[fcl2hits[ih0 - 1]].Integral()) / 2.;
4413 else if ((
unsigned short)fcl2hits.size() > fitLen) {
4415 std::vector<float> xwir;
4416 std::vector<float> ychg;
4417 std::vector<float> ychgerr2;
4419 unsigned int wire0 = fHits[fcl2hits[fcl2hits.size() - 1]].WireID().Wire;
4421 unsigned short npt = 0;
4422 unsigned short imlast = 0;
4426 for (
unsigned int ii = fcl2hits.size() - 1; ii > 0; --ii) {
4428 float chg = fHits[fcl2hits[ii]].Integral();
4431 if (npt == fitLen) {
4438 rms = std::sqrt((rms - fnpt * ave * ave) / (fnpt - 1));
4439 float chgcut = ave + rms;
4442 for (
unsigned short ii = fcl2hits.size() - 1; ii > imlast; --ii) {
4443 wire = fHits[fcl2hits[ii]].WireID().Wire;
4444 chg = fHits[fcl2hits[ii]].Integral();
4445 if (chg > chgcut)
continue;
4446 xwir.push_back((
float)(wire - wire0));
4447 ychg.push_back(chg);
4448 ychgerr2.push_back(chg);
4450 if (ychg.size() < 3)
return;
4456 fLinFitAlg.LinFit(xwir, ychg, ychgerr2, intcpt, slope, intcpterr, slopeerr, chidof);
4458 mf::LogVerbatim(
"CC") <<
"FitClusterChg wire " << wire0 <<
" chidof " << (int)chidof
4459 <<
" npt " << xwir.size() <<
" charge = " << (int)intcpt
4460 <<
" slope = " << (
int)slope <<
" first ave " << (int)ave <<
" rms " 4462 if (chidof > 100.)
return;
4465 if (intcpt > ave)
return;
4468 ave = intcpt / fAveChg;
4469 if (ave > 1.3)
return;
4470 if (ave < 0.77)
return;
4478 void ClusterCrawlerAlg::AddLAHit(
unsigned int kwire,
bool& ChkCharge,
bool& HitOK,
bool& SigOK)
4486 if (kwire < fFirstWire || kwire > fLastWire)
return;
4488 if (fcl2hits.size() == 0)
return;
4491 if (WireHitRange[kwire].first == -1) {
4496 if (WireHitRange[kwire].first == -2)
return;
4498 unsigned int firsthit = WireHitRange[kwire].first;
4499 unsigned int lasthit = WireHitRange[kwire].second;
4502 float timeDiff = 40 * AngleFactor(clpar[1]);
4506 unsigned int lastClHit = UINT_MAX;
4507 if (fcl2hits.size() > 0) {
4508 lastClHit = fcl2hits[fcl2hits.size() - 1];
4509 if (lastClHit == 0) {
4510 fAveChg = fHits[lastClHit].Integral();
4511 fAveHitWidth = fHits[lastClHit].EndTick() - fHits[lastClHit].StartTick();
4516 float prtime = clpar[0] + ((float)kwire - clpar[2]) * clpar[1];
4517 float chgWinLo = prtime - fChgNearWindow;
4518 float chgWinHi = prtime + fChgNearWindow;
4519 float chgrat, hitWidth;
4520 float hitWidthCut = 0.5 * fAveHitWidth;
4524 mf::LogVerbatim(
"CC") <<
"AddLAHit: wire " << kwire <<
" prtime " << prtime
4525 <<
" max timeDiff " << timeDiff <<
" fAveChg " << fAveChg;
4526 unsigned int imbest = 0;
4528 for (khit = firsthit; khit < lasthit; ++khit) {
4530 if (inClus[khit] < 0)
continue;
4531 dtime =
std::abs(fHits[khit].PeakTime() - prtime);
4532 hitWidth = fHits[khit].EndTick() - fHits[khit].StartTick();
4534 if (ChkCharge && fAveChg > 0) { chgrat = fHits[khit].Integral() / fAveChg; }
4536 mf::LogVerbatim(
"CC") <<
" Chk W:T " << kwire <<
":" << (short)fHits[khit].PeakTime()
4537 <<
" dT " << std::fixed << std::setprecision(1) << dtime <<
" InClus " 4538 << inClus[khit] <<
" mult " << fHits[khit].Multiplicity() <<
" width " 4539 << (int)hitWidth <<
" MergeAvail " << mergeAvailable[khit] <<
" Chi2 " 4540 << std::fixed << std::setprecision(2) << fHits[khit].GoodnessOfFit()
4541 <<
" Charge " << (int)fHits[khit].Integral() <<
" chgrat " 4542 << std::fixed << std::setprecision(1) << chgrat <<
" index " << khit;
4544 if (fHits[khit].PeakTime() > chgWinLo && fHits[khit].PeakTime() < chgWinHi)
4545 cnear += fHits[khit].Integral();
4547 if (prtime < fHits[khit].StartTick() - timeDiff)
continue;
4548 if (prtime > fHits[khit].EndTick() + timeDiff)
continue;
4551 if (inClus[khit] > 0)
continue;
4553 if (hitWidth < hitWidthCut)
continue;
4555 if (chgrat < 0.1)
continue;
4556 if (dtime < timeDiff) {
4568 mf::LogVerbatim(
"CC") <<
" Pick hit time " << (int)fHits[imbest].PeakTime() <<
" hit index " 4573 if (lastClHit != UINT_MAX && fHits[imbest].Multiplicity() > 1) {
4574 bool doMerge =
true;
4577 for (
unsigned short ivx = 0; ivx < vtx.size(); ++ivx) {
4578 if (vtx[ivx].CTP != clCTP)
continue;
4580 mf::LogVerbatim(
"CC") <<
" close vtx chk W:T " << vtx[ivx].Wire <<
":" 4581 << (int)vtx[ivx].Time;
4582 if (
std::abs(kwire - vtx[ivx].Wire) < 5 &&
4583 std::abs(
int(fHits[imbest].PeakTime() - vtx[ivx].Time)) < 20) {
4584 if (prt)
mf::LogVerbatim(
"CC") <<
" Close to a vertex. Don't merge hits";
4591 unsigned short nused = 0;
4593 float multipletChg = 0.;
4594 float chicut = AngleFactor(clpar[1]) * fHitMergeChiCut * fHits[lastClHit].RMS();
4596 std::pair<size_t, size_t> MultipletRange = FindHitMultiplet(imbest);
4597 for (
size_t jht = MultipletRange.first; jht < MultipletRange.second; ++jht) {
4598 if (inClus[jht] < 0)
continue;
4599 if (inClus[jht] == 0)
4600 multipletChg += fHits[jht].Integral();
4604 if (jht > MultipletRange.first) {
4606 float hitRMS = fHits[jht].RMS();
4607 if (fHits[jht - 1].RMS() > hitRMS) hitRMS = fHits[jht - 1].RMS();
4609 std::abs(fHits[jht].PeakTime() - fHits[jht - 1].PeakTime()) / hitRMS;
4610 if (prt)
mf::LogVerbatim(
"CC") <<
" Hit RMS chisq " << tdiff <<
" chicut " << chicut;
4611 if (tdiff > chicut) doMerge =
false;
4615 if (!doMerge)
mf::LogVerbatim(
"CC") <<
" Hits are well separated. Don't merge them ";
4617 if (doMerge && nused == 0) {
4622 float chgrat = multipletChg / fHits[lastClHit].Integral();
4624 mf::LogVerbatim(
"CC") <<
" merge hits charge check " << (int)multipletChg
4625 <<
" Previous hit charge " << (
int)fHits[lastClHit].Integral();
4626 if (chgrat > 2) doMerge =
false;
4634 MergeHits(imbest, didMerge);
4639 fcl2hits.push_back(imbest);
4642 chifits.push_back(clChisq);
4643 hitNear.push_back(hnear);
4645 cnear -= fHits[imbest].Integral();
4646 if (cnear < 0) cnear = 0;
4648 cnear /= fHits[imbest].Integral();
4649 chgNear.push_back(cnear);
4651 hitWidth = fHits[imbest].EndTick() - fHits[imbest].StartTick();
4653 << timeDiff <<
" clChisq " << clChisq <<
" Chg " 4654 << (int)fHits[imbest].Integral() <<
" AveChg " << (int)fAveChg
4655 <<
" width " << (
int)hitWidth <<
" AveWidth " << (int)fAveHitWidth
4656 <<
" fcl2hits size " << fcl2hits.size();
4659 if (clChisq > fChiCut[pass]) {
4664 if (prt)
mf::LogVerbatim(
"CC") <<
" LADD- Removed bad hit. Stopped tracking";
4669 bool ClusterCrawlerAlg::ClusterHitsOK(
short nHitChk)
4682 if (fcl2hits.size() == 0)
return true;
4684 unsigned short nHitToChk = fcl2hits.size();
4685 if (nHitChk > 0) nHitToChk = nHitChk + 1;
4686 unsigned short ii, indx;
4692 if (plane < geom->TPC(
geo::TPCID(cstat, tpc)).Nplanes() - 1) tol = 40;
4695 (fHits[fcl2hits[0]].PeakTime() > fHits[fcl2hits[fcl2hits.size() - 1]].PeakTime());
4698 for (ii = 0; ii < nHitToChk; ++ii) {
4699 indx = fcl2hits.size() - 1 - ii;
4700 mf::LogVerbatim(
"CC") <<
"ClusterHitsOK chk " << fHits[fcl2hits[indx]].WireID().Wire
4701 <<
" start " << fHits[fcl2hits[indx]].StartTick() <<
" peak " 4702 << fHits[fcl2hits[indx]].PeakTime() <<
" end " 4703 << fHits[fcl2hits[indx]].EndTick() <<
" posSlope " << posSlope;
4708 for (ii = 0; ii < nHitToChk - 1; ++ii) {
4709 indx = fcl2hits.size() - 1 - ii;
4712 fHits[fcl2hits[indx - 1]].
WireID().Wire) > 1)
4715 std::max(fHits[fcl2hits[indx]].StartTick(), fHits[fcl2hits[indx - 1]].StartTick());
4716 loEndTick = std::min(fHits[fcl2hits[indx]].EndTick(), fHits[fcl2hits[indx - 1]].EndTick());
4718 if (loEndTick + tol < hiStartTick) {
return false; }
4721 if (loEndTick + tol < hiStartTick) {
return false; }
4728 void ClusterCrawlerAlg::AddHit(
unsigned int kwire,
bool& HitOK,
bool& SigOK)
4739 if (kwire < fFirstWire || kwire > fLastWire)
return;
4741 unsigned int lastClHit = UINT_MAX;
4742 if (fcl2hits.size() > 0) lastClHit = fcl2hits[fcl2hits.size() - 1];
4745 unsigned int wire0 = clpar[2];
4748 if (fAllowNoHitWire == 0) {
4749 if (WireHitRange[kwire].first == -2)
return;
4753 if (WireHitRange[kwire].first == -2 && (wire0 - kwire) > fAllowNoHitWire) {
4759 if (WireHitRange[kwire].first == -1) {
4764 unsigned int firsthit = WireHitRange[kwire].first;
4765 unsigned int lasthit = WireHitRange[kwire].second;
4768 float dw = (float)kwire - (
float)wire0;
4769 float prtime = clpar[0] + dw * clpar[1];
4770 if (prtime < 0 || (
unsigned int)prtime > fMaxTime)
return;
4773 float prtimerr2 =
std::abs(dw) * clparerr[1] * clparerr[1];
4777 if (lastClHit != UINT_MAX) hiterr = 3 * fHitErrFac * fHits[lastClHit].RMS();
4778 float err = std::sqrt(prtimerr2 + hiterr * hiterr);
4780 float hitWin = fClProjErrFac * err;
4782 float prtimeLo = prtime - hitWin;
4783 float prtimeHi = prtime + hitWin;
4784 float chgWinLo = prtime - fChgNearWindow;
4785 float chgWinHi = prtime + fChgNearWindow;
4787 mf::LogVerbatim(
"CC") <<
"AddHit: wire " << kwire <<
" prtime Lo " << (int)prtimeLo
4788 <<
" prtime " << (
int)prtime <<
" Hi " << (int)prtimeHi <<
" prtimerr " 4789 << sqrt(prtimerr2) <<
" hiterr " << hiterr <<
" fAveChg " 4790 << (int)fAveChg <<
" fAveHitWidth " << std::setprecision(3)
4794 unsigned int imbest = INT_MAX;
4795 float best = 9999., dtime;
4797 float hitTime, hitChg, hitStartTick, hitEndTick;
4798 for (
unsigned int khit = firsthit; khit < lasthit; ++khit) {
4800 if (inClus[khit] < 0)
continue;
4801 hitTime = fHits[khit].PeakTime();
4802 dtime =
std::abs(hitTime - prtime);
4803 if (dtime > 1000)
continue;
4804 hitStartTick = fHits[khit].StartTick();
4805 hitEndTick = fHits[khit].EndTick();
4807 if (fAveChg > 0) dtime *=
std::abs(fHits[khit].Integral() - fAveChg) / fAveChg;
4808 if (prt &&
std::abs(dtime) < 100) {
4810 << std::setprecision(1) << (hitTime - prtime) <<
" InClus " 4811 << inClus[khit] <<
" mult " << fHits[khit].Multiplicity() <<
" RMS " 4812 << std::fixed << std::setprecision(1) << fHits[khit].RMS() <<
" Chi2 " 4813 << std::fixed << std::setprecision(1) << fHits[khit].GoodnessOfFit()
4814 <<
" Charge " << (int)fHits[khit].Integral() <<
" Peak " << std::fixed
4815 << std::setprecision(1) << fHits[khit].PeakAmplitude() <<
" LoT " 4816 << (int)hitStartTick <<
" HiT " << (
int)hitEndTick <<
" index " 4820 if (fHits[khit].StartTick() > chgWinLo && fHits[khit].EndTick() < chgWinHi)
4821 cnear += fHits[khit].Integral();
4823 if (prtimeHi < hitStartTick)
continue;
4824 if (prtimeLo > hitEndTick)
continue;
4827 if (hitTime < prtimeLo)
continue;
4828 if (hitTime > prtimeHi)
continue;
4830 if (inClus[khit] > 0)
continue;
4838 if (fAllowNoHitWire == 0)
return;
4840 mf::LogVerbatim(
"CC") <<
" wire0 " << wire0 <<
" kwire " << kwire <<
" max " 4841 << fAllowNoHitWire <<
" imbest " << imbest;
4842 if ((wire0 - kwire) > fAllowNoHitWire)
return;
4846 if (imbest == INT_MAX)
return;
4855 bool didMerge =
false;
4856 if (lastClHit != UINT_MAX && fAveHitWidth > 0 && fHitMergeChiCut > 0 &&
4858 bool doMerge =
true;
4859 for (
unsigned short ivx = 0; ivx < vtx.size(); ++ivx) {
4860 if (
std::abs(kwire - vtx[ivx].Wire) < 10 &&
4867 if (hit.
LocalIndex() != 0 && imbest == 0) doMerge =
false;
4871 if (hit.
LocalIndex() == 0) { oht = imbest + 1; }
4878 hitSep /= hit.
RMS();
4880 float totChg = hitChg + other_hit.
Integral();
4881 float lastHitChg = fAveChg;
4882 if (lastHitChg < 0) lastHitChg = fHits[lastClHit].Integral();
4885 mf::LogVerbatim(
"CC") <<
" Chk hit merge hitsep " << hitSep <<
" dChg " 4886 <<
std::abs(totChg - lastHitChg) <<
" Cut " 4888 if (inClus[oht] == 0 && hitSep < fHitMergeChiCut) {
4890 MergeHits(imbest, didMerge);
4900 float chgrat = (hitChg - fAveChg) / fAveChg;
4901 if (prt)
mf::LogVerbatim(
"CC") <<
" Chgrat " << std::setprecision(2) << chgrat;
4904 if (chgrat > 3 * fChgCut[pass]) {
4906 mf::LogVerbatim(
"CC") <<
" fails 3 x high charge cut " << fChgCut[pass] <<
" on pass " 4914 float bigchgcut = 1.5 * fChgCut[pass];
4915 bool lasthitbig =
false;
4916 bool lasthitlow =
false;
4918 float lastchgrat = (fHits[lastClHit].Integral() - fAveChg) / fAveChg;
4919 lasthitbig = (lastchgrat > bigchgcut);
4920 lasthitlow = (lastchgrat < -fChgCut[pass]);
4924 if (lasthitlow && chgrat < -fChgCut[pass]) {
4925 if (prt)
mf::LogVerbatim(
"CC") <<
" fails low charge cut. Stop crawling.";
4931 if (lasthitbig && chgrat > fChgCut[pass]) {
4933 mf::LogVerbatim(
"CC") <<
" fails 2nd hit high charge cut. Last hit was high also. ";
4938 if (chgrat > fChgCut[pass]) {
4939 if (best > 2 * err) {
4940 if (prt)
mf::LogVerbatim(
"CC") <<
" high charge && bad dT= " << best <<
" err= " << err;
4946 fitChg = (chgrat <
std::abs(fChgCut[pass]));
4950 fcl2hits.push_back(imbest);
4952 if (fcl2hits.size() == 3) std::sort(fcl2hits.begin(), fcl2hits.end(),
SortByLowHit);
4954 chifits.push_back(clChisq);
4955 hitNear.push_back(hnear);
4957 cnear -= fHits[imbest].Integral();
4958 if (cnear < 0) cnear = 0;
4960 cnear /= fHits[imbest].Integral();
4961 chgNear.push_back(cnear);
4966 if (chgNear.size() != fcl2hits.size()) {
4973 <<
" clChisq " << clChisq <<
" Chg " << (int)fHits[imbest].Integral()
4974 <<
" AveChg " << (int)fAveChg <<
" fcl2hits size " << fcl2hits.size();
4976 if (!fitChg)
return;
4982 void ClusterCrawlerAlg::ChkClusterNearbyHits(
bool prt)
4989 if (fHitMergeChiCut <= 0)
return;
4991 if (hitNear.size() != fcl2hits.size()) {
4992 mf::LogWarning(
"CC") <<
"Coding error: hitNear size != fcl2hits";
4997 if (hitNear.size() < 12)
return;
5000 unsigned short ii, indx;
5001 unsigned short merged = 0;
5002 unsigned short notmerged = 0;
5003 for (ii = 0; ii < 6; ++ii) {
5004 indx = hitNear.size() - 1 - ii;
5005 if (hitNear[indx] > 0) ++notmerged;
5006 if (hitNear[indx] < 0) ++merged;
5010 mf::LogVerbatim(
"CC") <<
"ChkClusterNearbyHits: nearby hits merged " << merged
5011 <<
" not merged " << notmerged;
5013 if (notmerged < 2)
return;
5019 for (ii = 0; ii < 6; ++ii) {
5020 indx = fcl2hits.size() - 1 - ii;
5021 const unsigned int iht = fcl2hits[indx];
5025 if (hit.
LocalIndex() != 0 && iht == 0)
continue;
5028 if (hit.
LocalIndex() == 0) { oht = iht + 1; }
5035 hitSep /= hit.
RMS();
5036 if (hitSep < fHitMergeChiCut && inClus[oht] == 0) {
5039 << fHits[iht].WireID().Wire <<
":" << fHits[iht].PeakTime();
5040 MergeHits(iht, didMerge);
5041 if (didMerge) hitNear[indx] = -1;
5050 if (prt)
mf::LogVerbatim(
"CC") <<
"ChkClusterNearbyHits refit cluster. fAveChg= " << fAveChg;
5055 void ClusterCrawlerAlg::FitVtx(
unsigned short iv)
5057 std::vector<float>
x;
5058 std::vector<float>
y;
5059 std::vector<float> ey2;
5063 if (vtx[iv].Fixed)
return;
5066 vtx[iv].ChiDOF = 99;
5070 std::vector<unsigned short> vcl;
5071 for (icl = 0; icl < tcl.size(); ++icl) {
5072 if (tcl[icl].ID < 0)
continue;
5073 if (tcl[icl].CTP != vtx[iv].CTP)
continue;
5074 if (tcl[icl].EndVtx == iv) vcl.push_back(icl);
5075 if (tcl[icl].BeginVtx == iv) vcl.push_back(icl);
5078 vtx[iv].NClusters = vcl.size();
5080 if (vcl.size() == 0)
return;
5086 unsigned short indx = tcl[icl].tclhits.size() - 1;
5087 unsigned int hit = tcl[icl].tclhits[indx];
5088 float minTimeErr = fHitErrFac * fHits[hit].RMS() * fHits[hit].Multiplicity();
5090 if (vcl.size() == 1) {
5093 if (tcl[icl].EndVtx == iv) {
5094 vtx[iv].Wire = tcl[icl].EndWir;
5095 vtx[iv].WireErr = 1;
5096 vtx[iv].Time = tcl[icl].EndTim;
5098 indx = tcl[icl].tclhits.size() - 1;
5099 hit = tcl[icl].tclhits[indx];
5100 vtx[iv].TimeErr = fHitErrFac * fHits[hit].RMS() * fHits[hit].Multiplicity();
5103 if (tcl[icl].BeginVtx == iv) {
5104 vtx[iv].Wire = tcl[icl].BeginWir;
5105 vtx[iv].WireErr = 1;
5106 vtx[iv].Time = tcl[icl].BeginTim;
5108 hit = tcl[icl].tclhits[0];
5109 vtx[iv].TimeErr = fHitErrFac * fHits[hit].RMS() * fHits[hit].Multiplicity();
5115 std::vector<double> slps;
5116 std::vector<double> slperrs;
5117 for (
unsigned short ii = 0; ii < vcl.size(); ++ii) {
5119 if (tcl[icl].EndVtx == iv) {
5120 x.push_back(tcl[icl].EndSlp);
5121 slps.push_back(tcl[icl].EndSlp);
5122 slperrs.push_back(tcl[icl].EndSlpErr);
5123 arg = tcl[icl].EndSlp * tcl[icl].EndWir - tcl[icl].EndTim;
5125 if (tcl[icl].EndSlpErr > 0.) { arg = tcl[icl].EndSlpErr * tcl[icl].EndWir; }
5127 arg = .1 * tcl[icl].EndWir;
5129 ey2.push_back(arg * arg);
5131 else if (tcl[icl].BeginVtx == iv) {
5132 x.push_back(tcl[icl].BeginSlp);
5133 slps.push_back(tcl[icl].BeginSlp);
5134 slperrs.push_back(tcl[icl].BeginSlpErr);
5135 arg = tcl[icl].BeginSlp * tcl[icl].BeginWir - tcl[icl].BeginTim;
5137 if (tcl[icl].BeginSlpErr > 0.) { arg = tcl[icl].BeginSlpErr * tcl[icl].BeginWir; }
5139 arg = .1 * tcl[icl].BeginWir;
5141 ey2.push_back(arg * arg);
5144 if (x.size() < 2)
return;
5147 double sumerr = 0, cnt = 0;
5148 for (
unsigned short ii = 0; ii < slps.size() - 1; ++ii) {
5149 for (
unsigned short jj = ii + 1; jj < slps.size(); ++jj) {
5150 arg = std::min(slperrs[ii], slperrs[jj]);
5151 arg /= (slps[ii] - slps[jj]);
5152 sumerr += arg * arg;
5159 float vTimeErr = 0.;
5161 float vWireErr = 0.;
5163 fLinFitAlg.LinFit(x, y, ey2, vTime, vWire, vTimeErr, vWireErr, chiDOF);
5164 if (chiDOF > 900)
return;
5167 if (vTime < 0 || vTime > fMaxTime)
return;
5170 if (vWire < 0 || vWire > geom->Nwires(iplID))
return;
5171 vtx[iv].ChiDOF = chiDOF;
5172 vtx[iv].Wire = vWire;
5173 vtx[iv].Time = vTime;
5174 vtx[iv].WireErr = vWire * sqrt(sumerr);
5175 vtx[iv].TimeErr = vTime * fabs(sumerr);
5177 if (vtx[iv].WireErr < 1) vtx[iv].WireErr = 2;
5179 if (vtx[iv].TimeErr < minTimeErr) vtx[iv].TimeErr = minTimeErr;
5189 if (
empty(vtx3))
return;
5191 const unsigned int cstat = tpcid.
Cryostat;
5192 const unsigned int tpc = tpcid.
TPC;
5194 unsigned int thePlane, theWire;
5198 for (
unsigned short ivx = 0; ivx < vtx3.size(); ++ivx) {
5200 if (vtx3[ivx].Wire < 0)
continue;
5201 if (vtx3[ivx].CStat != cstat || vtx3[ivx].TPC != tpc)
continue;
5204 theWire = vtx3[ivx].Wire;
5205 for (plane = 0; plane < 3; ++plane) {
5206 if (vtx3[ivx].Ptr2D[plane] >= 0)
continue;
5210 if (thePlane > 2)
continue;
5212 clCTP =
EncodeCTP(cstat, tpc, thePlane);
5215 vnew.
Wire = theWire;
5216 vnew.
Time = theTime;
5220 vtx.push_back(vnew);
5221 unsigned short ivnew = vtx.size() - 1;
5222 std::vector<short> vclIndex;
5223 for (
unsigned short icl = 0; icl < tcl.size(); ++icl) {
5224 if (tcl[icl].ID < 0)
continue;
5225 if (tcl[icl].CTP != clCTP)
continue;
5229 if (dwb > 10 && dwe > 10)
continue;
5230 if (dwb < dwe && dwb < 10 && tcl[icl].BeginVtx < 0) {
5232 if (theWire < tcl[icl].BeginWir + 5)
continue;
5233 if (ClusterVertexChi(icl, 0, ivnew) > fVertex3DCut)
continue;
5234 tcl[icl].BeginVtx = ivnew;
5235 vclIndex.push_back(icl);
5237 else if (dwe < 10 && tcl[icl].EndVtx < 0) {
5239 if (theWire > tcl[icl].EndWir - 5)
continue;
5240 if (ClusterVertexChi(icl, 1, ivnew) > fVertex3DCut)
continue;
5241 tcl[icl].EndVtx = ivnew;
5242 vclIndex.push_back(icl);
5245 bool goodVtx =
false;
5246 if (vclIndex.size() > 0) {
5248 goodVtx = (vtx[ivnew].ChiDOF < fVertex3DCut);
5249 vtx3[ivx].Ptr2D[thePlane] = ivnew;
5252 vtx3[ivx].Ptr2D[thePlane] = ivnew;
5253 vtx3[ivx].Wire = -1;
5258 for (
unsigned short ii = 0; ii < vclIndex.size(); ++ii) {
5259 unsigned short icl = vclIndex[ii];
5260 if (tcl[icl].BeginVtx == ivnew) tcl[icl].BeginVtx = -99;
5261 if (tcl[icl].EndVtx == ivnew) tcl[icl].EndVtx = -99;
5273 if (
empty(vtx3))
return;
5275 const unsigned int cstat = tpcid.
Cryostat;
5276 const unsigned int tpc = tpcid.
TPC;
5278 vtxprt = (fDebugPlane >= 0) && (fDebugHit == 6666);
5280 unsigned int lastplane = 5, kcl, kclID;
5282 unsigned int thePlane, theWire, plane;
5283 unsigned int loWire, hiWire;
5285 for (
unsigned short ivx = 0; ivx < vtx3.size(); ++ivx) {
5286 if (vtx3[ivx].CStat != cstat || vtx3[ivx].TPC != tpc)
continue;
5289 mf::LogVerbatim(
"CC") <<
"Vtx3ClusterSplit ivx " << ivx <<
" Ptr2D " << vtx3[ivx].Ptr2D[0]
5290 <<
" " << vtx3[ivx].Ptr2D[1] <<
" " << vtx3[ivx].Ptr2D[2] <<
" wire " 5292 if (vtx3[ivx].Wire < 0)
continue;
5295 theWire = vtx3[ivx].Wire;
5296 for (plane = 0; plane < 3; ++plane) {
5297 if (vtx3[ivx].Ptr2D[plane] >= 0)
continue;
5301 if (thePlane > 2)
continue;
5303 (
double)vtx3[ivx].
X, (
int)thePlane, (
int)tpcid.
TPC, (
int)tpcid.
Cryostat);
5305 if (thePlane != lastplane) {
5308 lastplane = thePlane;
5311 std::vector<short> clIDs;
5312 if (theWire > fFirstWire + 5) { loWire = theWire - 5; }
5314 loWire = fFirstWire;
5316 if (theWire < fLastWire - 5) { hiWire = theWire + 5; }
5321 mf::LogVerbatim(
"CC") <<
"3DVtx " << ivx <<
" look for cluster hits near P:W:T " << thePlane
5322 <<
":" << theWire <<
":" << (int)theTime <<
" Wire range " << loWire
5323 <<
" to " << hiWire;
5324 for (
unsigned int wire = loWire; wire < hiWire; ++wire) {
5326 if (WireHitRange[wire].first < 0)
continue;
5327 unsigned int firsthit = WireHitRange[wire].first;
5328 unsigned int lasthit = WireHitRange[wire].second;
5329 for (
unsigned int khit = firsthit; khit < lasthit; ++khit) {
5331 if (inClus[khit] <= 0)
continue;
5332 if ((
unsigned int)inClus[khit] > tcl.size() + 1) {
5333 mf::LogError(
"CC") <<
"Invalid hit InClus. " << khit <<
" " << inClus[khit];
5337 if (theTime < fHits[khit].StartTick() - 20)
continue;
5338 if (theTime > fHits[khit].EndTick() + 20)
continue;
5339 kclID = inClus[khit];
5342 if (tcl[kcl].ID < 0)
continue;
5344 if (tcl[kcl].tclhits.size() < 6)
continue;
5346 if (tcl[kcl].tclhits.size() > 100 &&
std::abs(tcl[kcl].BeginAng - tcl[kcl].EndAng) < 0.1)
5350 mf::LogVerbatim(
"CC") <<
"Bingo " << ivx <<
" plane " << thePlane <<
" wire " << wire
5351 <<
" hit " << fHits[khit].WireID().Wire <<
":" 5352 << (int)fHits[khit].PeakTime() <<
" inClus " << inClus[khit]
5353 <<
" P:W:T " << thePlane <<
":" << tcl[kcl].BeginWir <<
":" 5354 << (int)tcl[kcl].BeginTim;
5355 if (std::find(clIDs.begin(), clIDs.end(), kclID) == clIDs.end()) {
5356 clIDs.push_back(kclID);
5360 if (clIDs.size() == 0)
continue;
5362 for (
unsigned int ii = 0; ii < clIDs.size(); ++ii)
5365 unsigned short ii, icl, jj;
5372 unsigned short i2Dvx = 0;
5373 for (ii = 0; ii < 3; ++ii) {
5374 if (ii == thePlane)
continue;
5375 i2Dvx = vtx3[ivx].Ptr2D[ii];
5376 if (i2Dvx > vtx.size() - 1) {
5377 mf::LogError(
"CC") <<
"Vtx3ClusterSplit: Coding error";
5380 if (vtx[i2Dvx].TimeErr > tErr) tErr = vtx[i2Dvx].TimeErr;
5384 for (ii = 0; ii < clIDs.size(); ++ii) {
5385 icl = clIDs[ii] - 1;
5387 for (jj = 0; jj < tcl[icl].tclhits.size(); ++jj) {
5388 iht = tcl[icl].tclhits[jj];
5389 if (fHits[iht].
WireID().Wire < theWire) {
5391 if (jj > 3) nhitfit = -3;
5392 FitClusterMid(icl, iht, nhitfit);
5393 float doca = DoCA(-1, 1, theWire, theTime);
5395 mf::LogVerbatim(
"CC") <<
" cls " << icl <<
" DoCA " << doca <<
" tErr " << tErr;
5396 if ((doca / tErr) > 2) clIDs[ii] = -1;
5406 for (ii = 0; ii < clIDs.size(); ++ii)
5411 unsigned short nok = 0;
5412 for (ii = 0; ii < clIDs.size(); ++ii)
5413 if (clIDs[ii] >= 0) ++nok;
5414 if (nok == 0)
continue;
5418 vnew.
Wire = theWire;
5420 vnew.
Time = theTime;
5425 vtx.push_back(vnew);
5427 unsigned short ivnew = vtx.size() - 1;
5429 mf::LogVerbatim(
"CC") <<
"Make new 2D vtx " << ivnew <<
" in plane " << thePlane
5430 <<
" from 3D vtx " << ivx;
5431 vtx3[ivx].Ptr2D[thePlane] = ivnew;
5433 for (ii = 0; ii < clIDs.size(); ++ii) {
5434 if (clIDs[ii] < 0)
continue;
5435 icl = clIDs[ii] - 1;
5438 unsigned short pos = 0;
5439 for (
unsigned short jj = 0; jj < tcl[icl].tclhits.size(); ++jj) {
5440 iht = tcl[icl].tclhits[jj];
5441 if (fHits[iht].
WireID().Wire < theWire) {
5448 tcl[icl].BeginVtx = ivnew;
5451 else if (pos > tcl[icl].tclhits.size()) {
5453 tcl[icl].EndVtx = ivnew;
5458 if (vtxprt)
mf::LogVerbatim(
"CC") <<
"Split cluster " << clIDs[ii] <<
" at pos " << pos;
5459 if (!SplitCluster(icl, pos, ivnew)) {
5463 tcl[icl].ProcCode += 10000;
5464 tcl[tcl.size() - 1].ProcCode += 10000;
5469 if (vtx[ivnew].ChiDOF < 5 && vtx[ivnew].WireErr < fVertex2DWireErrCut) {
5471 vtx3[ivx].Wire = -1;
5475 mf::LogVerbatim(
"CC") <<
"Bad vtx fit " << ivnew <<
" ChiDOF " << vtx[ivnew].ChiDOF
5476 <<
" WireErr " << vtx[ivnew].WireErr;
5479 vtx3[ivx].Ptr2D[thePlane] = -1;
5481 for (jj = 0; jj < tcl.size(); ++jj) {
5482 if (tcl[jj].BeginVtx == ivnew) tcl[jj].BeginVtx = -99;
5483 if (tcl[jj].EndVtx == ivnew) tcl[jj].EndVtx = -99;
5498 unsigned int nPln = geom->TPC(
geo::TPCID(cstat, tpc)).Nplanes();
5499 if (nPln != 3)
return;
5501 float ew1, ew2, bw2, fvw;
5508 unsigned short longClIndex;
5509 unsigned short shortClIndex;
5510 unsigned short splitPos;
5512 std::array<std::vector<Hammer>, 3> hamrVec;
5516 for (ipl = 0; ipl < 3; ++ipl) {
5518 for (
unsigned short ic1 = 0; ic1 < tcl.size(); ++ic1) {
5519 if (tcl[ic1].ID < 0)
continue;
5521 if (tcl[ic1].tclhits.size() < 20)
continue;
5522 if (tcl[ic1].CTP != clCTP)
continue;
5524 if (tcl[ic1].EndVtx >= 0)
continue;
5525 ew1 = tcl[ic1].EndWir;
5526 for (
unsigned short ic2 = 0; ic2 < tcl.size(); ++ic2) {
5527 if (tcl[ic2].ID < 0)
continue;
5529 if (tcl[ic2].tclhits.size() > 20)
continue;
5531 if (tcl[ic2].tclhits.size() < 6)
continue;
5532 if (tcl[ic2].CTP != clCTP)
continue;
5533 ew2 = tcl[ic2].EndWir;
5534 bw2 = tcl[ic2].BeginWir;
5537 if (ew1 < ew2 || ew1 > bw2)
continue;
5541 unsigned short spos = 0;
5542 for (
unsigned short ii = 0; ii < tcl[ic2].tclhits.size(); ++ii) {
5543 unsigned int iht = tcl[ic2].tclhits[ii];
5544 float dw = fHits[iht].WireID().Wire - tcl[ic1].EndWir;
5545 float dt = fabs(fHits[iht].PeakTime() - tcl[ic1].EndTim - tcl[ic1].EndSlp * dw);
5552 if (ibst < 0)
continue;
5553 fvw = 0.5 + fHits[ibst].WireID().Wire;
5556 aHam.Wire = (0.5 + fvw);
5557 aHam.Tick = fHits[ibst].PeakTime();
5558 aHam.X = det_prop.
ConvertTicksToX((
double)aHam.Tick, (
int)ipl, (
int)tpc, (
int)cstat);
5559 aHam.longClIndex = ic1;
5560 aHam.shortClIndex = ic2;
5561 aHam.splitPos = spos;
5562 unsigned short indx = hamrVec[ipl].size();
5563 hamrVec[ipl].resize(indx + 1);
5564 hamrVec[ipl][indx] = aHam;
5571 unsigned short noham = 0;
5572 for (ipl = 0; ipl < 3; ++ipl)
5573 if (hamrVec[ipl].
size() == 0) ++noham;
5574 if (noham > 1)
return;
5581 float YLo = world.Y() - thetpc.
HalfHeight() + 1;
5582 float YHi = world.Y() + thetpc.
HalfHeight() - 1;
5583 float ZLo = world.Z() - thetpc.
Length() / 2 + 1;
5584 float ZHi = world.Z() + thetpc.
Length() / 2 - 1;
5589 unsigned short icl, jpl, jcl, kpl, splitPos;
5590 for (ipl = 0; ipl < 3; ++ipl) {
5591 if (hamrVec[ipl].
size() == 0)
continue;
5592 jpl = (ipl + 1) % nPln;
5593 kpl = (jpl + 1) % nPln;
5594 for (
unsigned short ii = 0; ii < hamrVec[ipl].size(); ++ii) {
5595 if (hamrVec[ipl][ii].Used)
continue;
5596 for (
unsigned short jj = 0; jj < hamrVec[jpl].size(); ++jj) {
5597 if (hamrVec[jpl][jj].Used)
continue;
5598 dX = hamrVec[ipl][ii].X - hamrVec[jpl][jj].X;
5599 if (fabs(dX) > fVertex3DCut)
continue;
5602 geom->IntersectionPoint(
geo::WireID{plane_i, hamrVec[ipl][ii].Wire},
5606 if (y < YLo || y > YHi || z < ZLo || z > ZHi)
continue;
5608 hamrVec[ipl][ii].Used =
true;
5609 hamrVec[jpl][jj].Used =
true;
5613 newVtx3.
X = 0.5 * (hamrVec[ipl][ii].X + hamrVec[jpl][jj].X);
5615 newVtx3.
XErr = fabs(hamrVec[ipl][ii].
X - hamrVec[jpl][jj].
X);
5620 newVtx3.
CStat = cstat;
5625 newVtx2.
Wire = hamrVec[ipl][ii].Wire;
5627 newVtx2.
Time = hamrVec[ipl][ii].Tick;
5630 newVtx2.
Fixed =
false;
5631 icl = hamrVec[ipl][ii].longClIndex;
5632 newVtx2.
CTP = tcl[icl].CTP;
5633 vtx.push_back(newVtx2);
5634 unsigned short ivnew = vtx.size() - 1;
5636 tcl[icl].EndVtx = ivnew;
5639 newVtx3.
Ptr2D[ipl] = (short)ivnew;
5641 icl = hamrVec[ipl][ii].shortClIndex;
5642 splitPos = hamrVec[ipl][ii].splitPos;
5643 if (!SplitCluster(icl, splitPos, ivnew))
return;
5646 newVtx2.
Wire = hamrVec[jpl][jj].Wire;
5647 newVtx2.
Time = hamrVec[jpl][jj].Tick;
5649 jcl = hamrVec[jpl][jj].longClIndex;
5650 newVtx2.
CTP = tcl[jcl].CTP;
5651 vtx.push_back(newVtx2);
5652 ivnew = vtx.size() - 1;
5654 tcl[jcl].EndVtx = ivnew;
5656 newVtx3.
Ptr2D[jpl] = (short)(vtx.size() - 1);
5659 jcl = hamrVec[jpl][jj].shortClIndex;
5660 splitPos = hamrVec[jpl][jj].splitPos;
5661 if (!SplitCluster(jcl, splitPos, vtx.size() - 1))
return;
5665 newVtx3.
Ptr2D[kpl] = -1;
5668 newVtx3.
Wire = geom->NearestWireID(WPos,
geo::PlaneID{cstat, tpc, kpl}).Wire;
5673 vtx3.push_back(newVtx3);
5693 float YLo = world.Y() - TPC.
HalfHeight() + 1;
5694 float YHi = world.Y() + TPC.
HalfHeight() - 1;
5695 float ZLo = world.Z() - TPC.
Length() / 2 + 1;
5696 float ZHi = world.Z() + TPC.
Length() / 2 - 1;
5698 vtxprt = (fDebugPlane >= 0) && (fDebugHit == 6666);
5706 float wirePitch = geom->WirePitch(
geo::PlaneID{tpcid, 0});
5709 std::vector<float> vX(vtx.size());
5710 std::vector<float> vXsigma(vtx.size());
5712 for (
unsigned short ivx = 0; ivx < vtx.size(); ++ivx) {
5713 if (vtx[ivx].NClusters == 0)
continue;
5715 if (iplID.
TPC != tpc || iplID.
Cryostat != cstat)
continue;
5720 (
double)(vtx[ivx].Time + vtx[ivx].TimeErr), (
int)iplID.
Plane, (
int)tpc, (
int)cstat);
5721 vXsigma[ivx] = fabs(vXp - vX[ivx]);
5725 std::array<std::vector<unsigned short>, 3> vIndex;
5726 unsigned short indx, ipl;
5727 for (
unsigned short ivx = 0; ivx < vtx.size(); ++ivx) {
5728 if (vtx[ivx].NClusters == 0)
continue;
5730 if (iplID.
TPC != tpc || iplID.
Cryostat != cstat)
continue;
5732 if (ipl > 2)
continue;
5733 indx = vIndex[ipl].size();
5734 vIndex[ipl].resize(indx + 1);
5735 vIndex[ipl][indx] = ivx;
5739 std::vector<short> vPtr;
5740 for (
unsigned short ii = 0; ii < vtx.size(); ++ii)
5744 std::vector<Vtx3Store> v3temp;
5746 double y = 0,
z = 0;
5749 unsigned short ii, jpl, jj, kpl, kk, ivx, jvx, kvx;
5750 unsigned int iWire, jWire;
5751 unsigned short v3dBest = 0;
5752 float xbest = 0, ybest = 0, zbest = 0;
5756 for (ipl = 0; ipl < 2; ++ipl) {
5758 for (ii = 0; ii < vIndex[ipl].size(); ++ii) {
5759 ivx = vIndex[ipl][ii];
5760 if (ivx > vtx.size() - 1) {
5765 if (vPtr[ivx] >= 0)
continue;
5766 iWire = vtx[ivx].Wire;
5767 float best = fVertex3DCut;
5771 std::array<short, 3> t2dIndex = {{-1, -1, -1}};
5772 std::array<short, 3> tmpIndex = {{-1, -1, -1}};
5773 for (jpl = ipl + 1; jpl < 3; ++jpl) {
5775 for (jj = 0; jj < vIndex[jpl].size(); ++jj) {
5776 jvx = vIndex[jpl][jj];
5777 if (jvx > vtx.size() - 1) {
5782 if (vPtr[jvx] >= 0)
continue;
5783 jWire = vtx[jvx].Wire;
5785 float dX = fabs(vX[ivx] - vX[jvx]);
5786 float dXSigma = sqrt(vXsigma[ivx] * vXsigma[ivx] + vXsigma[jvx] * vXsigma[jvx]);
5787 float dXChi = dX / dXSigma;
5791 <<
"VtxMatch: ipl " << ipl <<
" ivx " << ivx <<
" ivX " << vX[ivx] <<
" jpl " << jpl
5792 <<
" jvx " << jvx <<
" jvX " << vX[jvx] <<
" W:T " << (int)vtx[jvx].Wire <<
":" 5793 << (
int)vtx[jvx].Time <<
" dXChi " << dXChi <<
" fVertex3DCut " << fVertex3DCut;
5795 if (dXChi > fVertex3DCut)
continue;
5797 if (y < YLo || y > YHi || z < ZLo || z > ZHi)
continue;
5800 kpl = 3 - ipl - jpl;
5801 kX = 0.5 * (vX[ivx] + vX[jvx]);
5805 kWire = geom->NearestWireID(WPos,
geo::PlaneID{cstat, tpc, kpl}).Wire;
5811 kpl = 3 - ipl - jpl;
5815 tmpIndex[ipl] = ivx;
5816 tmpIndex[jpl] = jvx;
5818 v3d.
Ptr2D = tmpIndex;
5822 float yzSigma = wirePitch * sqrt(vtx[ivx].WireErr * vtx[ivx].WireErr +
5823 vtx[jvx].WireErr * vtx[jvx].WireErr);
5830 v3temp.push_back(v3d);
5834 <<
"VtxMatch: 2 Plane match ivx " << ivx <<
" P:W:T " << ipl <<
":" 5835 << (int)vtx[ivx].Wire <<
":" << (
int)vtx[ivx].Time <<
" jvx " << jvx <<
" P:W:T " 5836 << jpl <<
":" << (int)vtx[jvx].Wire <<
":" << (
int)vtx[jvx].Time <<
" dXChi " 5837 << dXChi <<
" yzSigma " << yzSigma;
5839 if (TPC.
Nplanes() == 2)
continue;
5841 best = fVertex3DCut;
5842 for (kk = 0; kk < vIndex[kpl].size(); ++kk) {
5843 kvx = vIndex[kpl][kk];
5844 if (vPtr[kvx] >= 0)
continue;
5846 float dW = wirePitch * (vtx[kvx].Wire -
kWire) / yzSigma;
5848 float dX = (vX[kvx] -
kX) / dXSigma;
5849 float kChi = 0.5 * sqrt(dW * dW + dX * dX);
5852 xbest = (vX[kvx] + 2 *
kX) / 3;
5855 t2dIndex[ipl] = ivx;
5856 t2dIndex[jpl] = jvx;
5857 t2dIndex[kpl] = kvx;
5858 v3dBest = v3temp.size() - 1;
5863 <<
" kvx " << kvx <<
" kpl " << kpl <<
" wire " << (int)vtx[kvx].Wire <<
" kTime " 5864 << (
int)vtx[kvx].Time <<
" kChi " << kChi <<
" best " << best <<
" dW " 5865 << vtx[kvx].Wire -
kWire;
5869 mf::LogVerbatim(
"CC") <<
" done best = " << best <<
" fVertex3DCut " << fVertex3DCut;
5870 if (TPC.
Nplanes() > 2 && best < fVertex3DCut) {
5872 if (v3dBest > v3temp.size() - 1) {
5873 mf::LogError(
"CC") <<
"VtxMatch: bad v3dBest " << v3dBest;
5877 v3d.
Ptr2D = t2dIndex;
5883 vtx3.push_back(v3d);
5886 for (
unsigned short jj = 0; jj < 3; ++jj)
5887 if (t2dIndex[jj] >= 0) vPtr[t2dIndex[jj]] = vtx3.size() - 1;
5891 <<
"New 3D vtx " << vtx3.size() <<
" X " << v3d.
X <<
" Y " << v3d.
Y <<
" Z " 5892 << v3d.
Z <<
" t2dIndex " << t2dIndex[0] <<
" " << t2dIndex[1] <<
" " 5893 << t2dIndex[2] <<
" best Chi " << best;
5905 unsigned short vsize = vtx3.size();
5906 for (
unsigned short it = 0; it < v3temp.size(); ++it) {
5908 for (
unsigned short i3d = 0; i3d < vsize; ++i3d) {
5909 for (
unsigned short plane = 0; plane < 3; ++plane) {
5910 if (v3temp[it].Ptr2D[plane] == vtx3[i3d].Ptr2D[plane]) {
5918 if (keepit) vtx3.push_back(v3temp[it]);
5923 for (
unsigned short iv3 = 0; iv3 < vtx3.size(); ++iv3) {
5924 vtx3[iv3].Ptr2D[2] = 666;
5929 for (
unsigned short it = 0; it < vtx3.size(); ++it) {
5930 mf::LogVerbatim(
"CC") <<
"vtx3 " << it <<
" Ptr2D " << vtx3[it].Ptr2D[0] <<
" " 5931 << vtx3[it].Ptr2D[1] <<
" " << vtx3[it].Ptr2D[2] <<
" wire " 5939 void ClusterCrawlerAlg::GetHitRange(
CTP_t CTP)
5945 unsigned int nwires = geom->Nwires(planeID);
5946 WireHitRange.resize(nwires + 1);
5952 unsigned int wire, iht;
5953 unsigned int nHitInPlane;
5954 std::pair<int, int> flag;
5959 for (
auto& apair : WireHitRange)
5964 std::vector<bool> firsthit;
5965 firsthit.resize(nwires + 1,
true);
5966 bool firstwire =
true;
5967 for (iht = 0; iht < fHits.size(); ++iht) {
5968 if (fHits[iht].
WireID().asPlaneID() != planeID)
continue;
5969 wire = fHits[iht].WireID().Wire;
5971 if (firsthit[wire]) {
5972 WireHitRange[wire].first = iht;
5973 firsthit[wire] =
false;
5979 WireHitRange[wire].second = iht + 1;
5980 fLastWire = wire + 1;
5984 lariov::ChannelStatusProvider
const& channelStatus =
5989 unsigned int nbad = 0;
5990 for (wire = 0; wire < nwires; ++wire) {
5992 if (!channelStatus.IsGood(chan)) {
5993 WireHitRange[wire] = flag;
5998 if (mergeAvailable.size() < fHits.size())
6000 <<
"GetHitRange: Invalid mergeAvailable vector size " << mergeAvailable.size()
6002 unsigned int firstHit, lastHit;
6005 float maxRMS, chiSep, peakCut;
6006 for (wire = 0; wire < nwires; ++wire) {
6008 if (WireHitRange[wire].first < 0)
continue;
6009 firstHit = WireHitRange[wire].first;
6010 lastHit = WireHitRange[wire].second;
6011 for (iht = firstHit; iht < lastHit; ++iht) {
6012 if (fHits[iht].
WireID().Wire != wire)
6014 <<
"Bad WireHitRange on wire " << wire <<
"\n";
6016 if (fHits[iht].Multiplicity() > 1) {
6017 peakCut = 0.6 * fHits[iht].PeakAmplitude();
6018 std::pair<size_t, size_t> MultipletRange = FindHitMultiplet(iht);
6019 for (
size_t jht = MultipletRange.first; jht < MultipletRange.second; ++jht) {
6020 if (jht == iht)
continue;
6022 if (fHits[jht].PeakAmplitude() < peakCut)
continue;
6023 maxRMS = std::max(fHits[iht].RMS(), fHits[jht].RMS());
6024 chiSep =
std::abs(fHits[iht].PeakTime() - fHits[jht].PeakTime()) / maxRMS;
6025 if (chiSep < fHitMergeChiCut) {
6026 mergeAvailable[iht] =
true;
6033 if (cnt != nHitInPlane)
6034 mf::LogWarning(
"CC") <<
"Bad WireHitRange count " << cnt <<
" " << nHitInPlane <<
"\n";
6036 if (!fMergeAllHits)
return;
6040 for (wire = 0; wire < nwires; ++wire) {
6041 if (WireHitRange[wire].first < 0)
continue;
6042 firstHit = WireHitRange[wire].first;
6043 lastHit = WireHitRange[wire].second;
6044 for (iht = firstHit; iht < lastHit; ++iht) {
6045 if (!mergeAvailable[iht])
continue;
6047 if (fHits[iht].GoodnessOfFit() == 6666)
continue;
6048 MergeHits(iht, didMerge);
6049 mergeAvailable[iht] =
false;
6058 if (inWire1 > inWire2) {
6060 unsigned int tmp = inWire1;
6065 unsigned int wire, ndead = 0;
6066 for (wire = inWire1; wire < inWire2; ++wire)
6067 if (WireHitRange[wire].first == -1) ++ndead;
6075 if (fcl2hits.size() < 2)
return 0;
6076 unsigned int wire, ndead = 0;
6077 unsigned int iht = fcl2hits[fcl2hits.size() - 1];
6078 unsigned int eWire = fHits[iht].WireID().Wire;
6080 unsigned int bWire = fHits[iht].WireID().Wire + 1;
6081 for (wire = eWire; wire < bWire; ++wire)
6082 if (WireHitRange[wire].first == -1) ++ndead;
6087 bool ClusterCrawlerAlg::areInSameMultiplet(
recob::Hit const& first_hit,
6095 std::pair<size_t, size_t> ClusterCrawlerAlg::FindHitMultiplet(
size_t iHit)
const 6097 std::pair<size_t, size_t> range{iHit, iHit};
6099 range.first = iHit - fHits[iHit].LocalIndex();
6100 range.second = range.first + fHits[iHit].Multiplicity();
6106 bool ClusterCrawlerAlg::CheckHitDuplicates(std::string location,
6107 std::string marker )
const 6110 unsigned int nDuplicates = 0;
6111 std::set<unsigned int>
hits;
6112 for (
unsigned int hit : fcl2hits) {
6113 if (hits.count(
hit)) {
6116 log <<
"Hit #" <<
hit 6117 <<
" being included twice in the future cluster (ID=" << (tcl.size() + 1)
6118 <<
"?) at location: " << location;
6119 if (!marker.empty()) log <<
" (marker: '" << marker <<
"')";
6123 return nDuplicates > 0;
6127 float ClusterCrawlerAlg::DoCA(
short icl,
unsigned short end,
float vwire,
float vtick)
6132 if (icl > (
short)tcl.size())
return 9999;
6134 float cwire, cslp, ctick;
6137 if (fcl2hits.size() == 0)
return 9999;
6153 cwire = tcl[icl].BeginWir;
6154 cslp = tcl[icl].BeginSlp;
6155 ctick = tcl[icl].BeginTim;
6158 cwire = tcl[icl].EndWir;
6159 cslp = tcl[icl].EndSlp;
6160 ctick = tcl[icl].EndTim;
6165 float docaW = (vwire + cslp * (vtick - ctick) + cwire * cslp * cslp) / (1 + cslp * cslp);
6166 float dW = docaW - vwire;
6167 float dT = ctick + (vwire - cwire) * cslp - vtick;
6168 return sqrt(dW * dW + dT * dT);
6173 float ClusterCrawlerAlg::ClusterVertexChi(
short icl,
unsigned short end,
unsigned short ivx)
6177 if (icl > (
short)tcl.size())
return 9999;
6178 if (ivx > vtx.size())
return 9999;
6180 float cwire, cslp, cslpErr, ctick;
6183 if (fcl2hits.size() == 0)
return 9999;
6188 cslpErr = clBeginSlpErr;
6194 cslpErr = clparerr[1];
6201 cwire = tcl[icl].BeginWir;
6202 cslp = tcl[icl].BeginSlp;
6203 cslpErr = tcl[icl].BeginSlpErr;
6204 ctick = tcl[icl].BeginTim;
6207 cwire = tcl[icl].EndWir;
6208 cslp = tcl[icl].EndSlp;
6209 cslpErr = tcl[icl].EndSlpErr;
6210 ctick = tcl[icl].EndTim;
6216 (vtx[ivx].Wire + cslp * (vtx[ivx].Time - ctick) + cwire * cslp * cslp) / (1 + cslp * cslp);
6217 float dW = docaW - vtx[ivx].Wire;
6218 float chi = dW / vtx[ivx].WireErr;
6219 float totChi = chi * chi;
6220 dW = vtx[ivx].Wire - cwire;
6221 float dT = ctick + dW * cslp - vtx[ivx].Time;
6222 if (cslpErr < 1
E-3) cslpErr = 1
E-3;
6224 cslpErr *= AngleFactor(cslp);
6226 float dTErr = dW * cslpErr;
6230 dTErr += vtx[ivx].TimeErr * vtx[ivx].TimeErr;
6231 if (dTErr < 1
E-3) dTErr = 1
E-3;
6232 totChi += dT * dT / dTErr;
6240 float ClusterCrawlerAlg::PointVertexChi(
float wire,
float tick,
unsigned short ivx)
6244 if (ivx > vtx.size())
return 9999;
6246 float dW = wire - vtx[ivx].Wire;
6247 float chi = dW / vtx[ivx].WireErr;
6248 float totChi = chi * chi;
6249 float dT = tick - vtx[ivx].Time;
6250 chi = dT / vtx[ivx].TimeErr;
6251 totChi += chi * chi;
6261 if (iht > fHits.size() - 1)
return "Bad Hit";
short int LocalIndex() const
How well do we believe we know this hit?
MaybeLogger_< ELseverityLevel::ELsev_info, true > LogVerbatim
Functions to help with numbers.
geo::SigType_t SignalType() const
Signal type for the plane of the hit.
constexpr auto const & right(const_AssnsIter< L, R, D, Dir > const &a, const_AssnsIter< L, R, D, Dir > const &b)
Utilities related to art service access.
static unsigned int kWire
Encapsulate the construction of a single cyostat.
struct of temporary 2D vertices (end points)
unsigned int Nplanes() const
Number of planes in this tpc.
float RMS() const
RMS of the hit shape, in tick units.
Declaration of signal hit object.
Planes which measure X direction.
The data type to uniquely identify a Plane.
Geometry information for a single TPC.
double Temperature() const
In kelvin.
float SigmaPeakAmplitude() const
Uncertainty on estimated amplitude of the hit at its peak, in ADC units.
float SigmaIntegral() const
Initial tdc tick for hit.
constexpr auto abs(T v)
Returns the absolute value of the argument.
int DegreesOfFreedom() const
Initial tdc tick for hit.
CryostatID_t Cryostat
Index of cryostat.
float Integral() const
Integral under the calibrated signal waveform of the hit, in tick x ADC units.
struct of temporary clusters
geo::View_t View() const
View for the plane of the hit.
WireID_t Wire
Index of the wire within its plane.
geo::WireID const & WireID() const
Initial tdc tick for hit.
MaybeLogger_< ELseverityLevel::ELsev_error, false > LogError
Cluster finding and building.
float GoodnessOfFit() const
Degrees of freedom in the determination of the hit signal shape (-1 by default)
short int Multiplicity() const
How many hits could this one be shared with.
struct of temporary 3D vertices
float PeakAmplitude() const
The estimated amplitude of the hit at its peak, in ADC units.
double Length() const
Length is associated with z coordinate [cm].
constexpr int cmp(WireID const &other) const
Returns < 0 if this is smaller than tpcid, 0 if equal, > 0 if larger.
double Efield(unsigned int planegap=0) const
kV/cm
int TDCtick_t
Type representing a TDC tick.
decltype(auto) constexpr end(T &&obj)
ADL-aware version of std::end.
decltype(auto) constexpr size(T &&obj)
ADL-aware version of std::size.
Access the description of detector geometry.
Point_t GetCenter() const
Returns the center of the TPC volume in world coordinates [cm].
Collection of exceptions for Geometry system.
IDparameter< geo::WireID > WireID
Member type of validated geo::WireID parameter.
float DeadWireCount(const TCSlice &slc, const TrajPoint &tp1, const TrajPoint &tp2)
decltype(auto) constexpr to_string(T &&obj)
ADL-aware version of std::to_string.
T get(std::string const &key) const
tick_as<> tick
Tick number, represented by std::ptrdiff_t.
double ConvertXToTicks(double X, int p, int t, int c) const
void MergeOverlap(std::string inFcnLabel, TCSlice &slc, const CTP_t &inCTP, bool prt)
unsigned int NumberTimeSamples() const
std::string PrintHit(const TCHit &tch)
bool SortByLowHit(unsigned int i, unsigned int j)
bool has_key(std::string const &key) const
double DriftVelocity(double efield=0., double temperature=0.) const
cm/us
The data type to uniquely identify a TPC.
PlaneID_t Plane
Index of the plane within its TPC.
raw::TDCtick_t StartTick() const
Initial tdc tick for hit.
float PeakTimeMinusRMS(float sigmas=+1.) const
Returns a time sigmas RMS away from the peak time.
std::array< short, 3 > Ptr2D
float HitSummedADC() const
The sum of calibrated ADC counts of the hit (0. by default)
Detector simulation of raw signals on wires.
double ConvertTicksToX(double ticks, int p, int t, int c) const
raw::TDCtick_t EndTick() const
Final tdc tick for hit.
double HalfHeight() const
Height is associated with y coordinate [cm].
bool greaterThan(CluLen c1, CluLen c2)
cet::coded_exception< errors::ErrorCodes, ExceptionDetail::translate > Exception
float ROISummedADC() const
The sum of calibrated ADC counts of the ROI (0. by default)
ROOT::Math::PositionVector3D< ROOT::Math::Cartesian3D< double >, ROOT::Math::GlobalCoordinateSystemTag > Point_t
Type for representation of position in physical 3D space.
float PeakTime() const
Time of the signal peak, in tick units.
Contains all timing reference information for the detector.
std::vector< unsigned int > tclhits
geo::PlaneID DecodeCTP(CTP_t CTP)
constexpr auto absDiff(A const &a, B const &b)
Returns the absolute value of the difference between two values.
MaybeLogger_< ELseverityLevel::ELsev_warning, false > LogWarning
CTP_t EncodeCTP(unsigned int cryo, unsigned int tpc, unsigned int plane)
span(IterB &&b, IterE &&e, Adaptor &&adaptor) -> span< decltype(adaptor(std::forward< IterB >(b))), decltype(adaptor(std::forward< IterE >(e))) >
Exception thrown on invalid wire number.
decltype(auto) constexpr begin(T &&obj)
ADL-aware version of std::begin.
float SigmaPeakTime() const
Uncertainty for the signal peak, in tick units.
2D representation of charge deposited in the TDC/wire plane
float PeakTimePlusRMS(float sigmas=+1.) const
Returns a time sigmas RMS away from the peak time.
unsigned int ChannelID_t
Type representing the ID of a readout channel.
TPCID_t TPC
Index of the TPC within its cryostat.
second_as<> second
Type of time stored in seconds, in double precision.
geo::WireID suggestedWireID() const
Returns a better wire ID.
double sampling_rate(DetectorClocksData const &data)
Returns the period of the TPC readout electronics clock.
recob::tracking::Plane Plane
raw::ChannelID_t Channel() const
ID of the readout channel the hit was extracted from.
art framework interface to geometry description
decltype(auto) constexpr empty(T &&obj)
ADL-aware version of std::empty.
Encapsulate the construction of a single detector plane.