50 return (c1.length > c2.length);
54 return (c1.length < c2.length);
127 std::array<float, 2>
X;
144 std::array<std::vector<clPar>, 3>
cls;
149 std::array<float, 2>
X;
180 std::array<std::vector<unsigned short>, 3>
vxCls;
185 std::array<std::vector<art::Ptr<recob::Hit>>, 3>
TrkHits;
201 std::array<std::vector<art::Ptr<recob::Hit>>, 3>
trkHits;
203 std::array<std::vector<art::Ptr<recob::Hit>>, 3>
seedHits;
215 std::array<unsigned short, 3>
End;
235 void PrintTracks()
const;
243 unsigned short end1);
249 void FindMaybeVertices(
geo::TPCID const& tpcid);
272 float ChargeAsym(std::array<float, 3>& mChg);
274 bool FindMissingCluster(
unsigned short kpl,
276 unsigned short& kend,
282 bool DupMatch(
MatchPars& match,
unsigned short nplanes);
286 unsigned short procCode,
292 unsigned short nplanes);
298 unsigned short procCode,
303 unsigned short wire1,
305 unsigned short wire2,
309 float AngleFactor(
float slope);
320 fMatchAlgs = pset.get<std::vector<short>>(
"MatchAlgs");
321 fXMatchErr = pset.get<std::vector<float>>(
"XMatchErr");
324 fMatchMinLen = pset.get<std::vector<float>>(
"MatchMinLen");
327 fMaxDAng = pset.get<
float>(
"MaxDAng");
340 fuBCode = pset.get<
bool>(
"uBCode");
348 if (fMatchAlgs.size() > fXMatchErr.size() || fMatchAlgs.size() > fAngleMatchErr.size() ||
349 fMatchAlgs.size() > fChgAsymFactor.size() || fMatchAlgs.size() > fMatchMinLen.size() ||
350 fMatchAlgs.size() > fMakeAlgTracks.size()) {
351 mf::LogError(
"CCTM") <<
"Incompatible fcl input vector sizes";
355 for (
unsigned short ii = 0; ii < fMatchAlgs.size(); ++ii) {
356 if (fAngleMatchErr[ii] <= 0 || fXMatchErr[ii] <= 0) {
357 mf::LogError(
"CCTM") <<
"Invalid matching parameters " << fAngleMatchErr[ii] <<
" " 363 produces<std::vector<recob::PFParticle>>();
364 produces<art::Assns<recob::PFParticle, recob::Track>>();
365 produces<art::Assns<recob::PFParticle, recob::Cluster>>();
366 produces<art::Assns<recob::PFParticle, recob::Seed>>();
367 produces<art::Assns<recob::PFParticle, recob::Vertex>>();
368 produces<std::vector<recob::Vertex>>();
369 produces<std::vector<recob::Track>>();
370 produces<art::Assns<recob::Track, recob::Hit>>();
371 produces<std::vector<recob::Seed>>();
381 auto tcol = std::make_unique<std::vector<recob::Track>>();
382 auto thassn = std::make_unique<art::Assns<recob::Track, recob::Hit>>();
383 auto vcol = std::make_unique<std::vector<recob::Vertex>>();
384 auto pcol = std::make_unique<std::vector<recob::PFParticle>>();
385 auto ptassn = std::make_unique<art::Assns<recob::PFParticle, recob::Track>>();
386 auto pcassn = std::make_unique<art::Assns<recob::PFParticle, recob::Cluster>>();
387 auto psassn = std::make_unique<art::Assns<recob::PFParticle, recob::Seed>>();
388 auto pvassn = std::make_unique<art::Assns<recob::PFParticle, recob::Vertex>>();
391 auto scol = std::make_unique<std::vector<recob::Seed>>();
392 auto shassn = std::make_unique<art::Assns<recob::Seed, recob::Hit>>();
398 std::vector<art::Ptr<recob::Cluster>> clusterlist;
401 std::vector<art::Ptr<recob::Vertex>> vtxlist;
411 if (clusterlist.size() == 0)
return;
420 std::vector<CluLen> clulens;
422 unsigned short ipl, icl,
end, itr, tID, tIndex;
430 std::vector<art::Ptr<recob::Hit>> tmpHits;
431 std::vector<art::Ptr<recob::Cluster>> tmpCls;
432 std::vector<art::Ptr<recob::Vertex>> tmpVtx;
435 std::vector<size_t> dtrIndices;
438 double sPos[3], sDir[3];
439 double sErr[3] = {0, 0, 0};
445 std::vector<art::Ptr<recob::Hit>> clusterhits;
446 for (icl = 0; icl < clusterlist.size(); ++icl) {
447 ipl = clusterlist[icl]->Plane().Plane;
448 clusterhits = fmCluHits.at(icl);
449 if (clusterhits[0]->
WireID().Wire != std::nearbyint(clusterlist[icl]->EndWire())) {
450 std::cout <<
"CCTM Cluster-Hit End wire mis-match " << clusterhits[0]->WireID().Wire
451 <<
" vs " << std::nearbyint(clusterlist[icl]->EndWire()) <<
" Bail out! \n";
454 for (
unsigned short iht = 0; iht < clusterhits.size(); ++iht) {
455 if (clusterhits[iht]->
WireID().Plane != ipl) {
456 std::cout <<
"CCTM Cluster-Hit plane mis-match " << ipl <<
" vs " 457 << clusterhits[iht]->WireID().Plane <<
" on hit " << iht <<
" Bail out! \n";
467 unsigned int const nplanes = tpcgeom.Nplanes();
468 if (nplanes > 3)
continue;
469 for (ipl = 0; ipl < 3; ++ipl) {
479 for (icl = 0; icl < clusterlist.size(); ++icl) {
480 if (clusterlist[icl]->
Plane() != planeid)
continue;
483 clulen.length = clusterlist[icl]->EndWire();
484 clulens.push_back(clulen);
486 if (clulens.empty())
continue;
488 std::sort(clulens.begin(), clulens.end(),
lessThan);
489 if (clulens.empty())
continue;
490 for (
unsigned short ii = 0; ii < clulens.size(); ++ii) {
491 const unsigned short icl = clulens[ii].index;
498 clstr.
X[1] = (float)detProp.ConvertTicksToX(cluster.
StartTick(), planeid);
512 clstr.
X[0] = (float)detProp.ConvertTicksToX(cluster.
EndTick(), planeid);
516 if (clstr.
Time[1] > clstr.
Time[0]) {
533 clstr.
Length = (
unsigned short)(0.5 + clstr.
Wire[1] - clstr.
Wire[0]);
536 clusterhits = fmCluHits.at(icl);
537 if (clusterhits.size() == 0) {
538 mf::LogError(
"CCTM") <<
"No associated cluster hits " << icl;
542 clstr.
TotChg *= clstr.
Length / (float)clusterhits.size();
543 cls[ipl].push_back(clstr);
551 for (
unsigned short ivx = 0; ivx < vtxlist.size(); ++ivx) {
555 vtxlist[ivx]->XYZ(xyz);
562 std::vector<art::Ptr<recob::Cluster>>
const& vtxCls = fmVtxCls.at(ivx);
563 std::vector<const unsigned short*>
const& vtxClsEnd = fmVtxCls.
data(ivx);
564 for (
unsigned short vcass = 0; vcass < vtxCls.size(); ++vcass) {
565 icl = vtxCls[vcass].key();
567 ipl = vtxCls[vcass]->Plane().Plane;
568 end = *vtxClsEnd[vcass];
571 <<
"Invalid end data from vertex - cluster association" <<
end;
575 for (
unsigned short jcl = 0; jcl <
cls[ipl].size(); ++jcl) {
576 if (
cls[ipl][jcl].EvtIndex == icl) {
577 cls[ipl][jcl].VtxIndex[
end] = ivx;
584 throw cet::exception(
"CCTM") <<
"Bad index from vertex - cluster association" << icl;
597 VtxMatch(detProp, fmCluHits, tpcgeom.ID());
617 for (
unsigned short ipf = 0; ipf <
pfpToTrkID.size(); ++ipf) {
620 if (tID >
trk.size() + 1) {
627 for (
unsigned short jpf = 0; jpf <
pfpToTrkID.size(); ++jpf) {
629 if (
trk[itr].MomID == tID) dtrIndices.push_back(jpf);
630 if (
trk[itr].MomID == tID)
633 unsigned short parentIndex = USHRT_MAX;
637 pcol->emplace_back(std::move(pfp));
638 for (
unsigned short ivx = 0; ivx <
vtx.size(); ++ivx) {
639 if (!
vtx[ivx].Neutrino)
continue;
644 size_t vStart = vcol->size();
646 vcol->emplace_back(std::move(vertex));
647 size_t vEnd = vcol->size();
650 vtx[ivx].ID = -
vtx[ivx].ID;
659 for (
unsigned short ii = 0; ii <
pfpToTrkID.size(); ++ii) {
667 <<
" ipf " << ipf <<
" parentIndex " << parentIndex
668 <<
" dtr size " << dtrIndices.size();
670 pcol->emplace_back(std::move(pfp));
673 size_t tStart = tcol->size();
685 tcol->emplace_back(std::move(
track));
686 size_t tEnd = tcol->size();
690 trk[tIndex].ID = -
trk[tIndex].ID;
693 for (ipl = 0; ipl < nplanes; ++ipl)
695 tmpHits.end(),
trk[tIndex].TrkHits[ipl].begin(),
trk[tIndex].TrkHits[ipl].end());
699 unsigned short itj = 0;
700 if (end > 0) itj =
trk[tIndex].TrjPos.size() - 1;
701 for (
unsigned short ii = 0; ii < 3; ++ii) {
702 sPos[ii] =
trk[tIndex].TrjPos[itj](ii);
703 sDir[ii] =
trk[tIndex].TrjDir[itj](ii);
705 size_t sStart = scol->size();
707 scol->emplace_back(std::move(seed));
708 size_t sEnd = scol->size();
713 for (ipl = 0; ipl < nplanes; ++ipl)
719 for (
unsigned short ii = 0; ii <
trk[tIndex].ClsEvtIndices.size(); ++ii) {
720 icl =
trk[tIndex].ClsEvtIndices[ii];
721 tmpCls.push_back(clusterlist[icl]);
727 for (
unsigned short itr = 0; itr <
trk.size(); ++itr) {
729 if (
trk[itr].ID < 0)
continue;
741 tcol->emplace_back(std::move(
track));
743 for (ipl = 0; ipl < nplanes; ++ipl)
744 tmpHits.insert(tmpHits.end(),
trk[itr].TrkHits[ipl].begin(),
trk[itr].TrkHits[ipl].end());
748 for (
unsigned short ivx = 0; ivx <
vtx.size(); ++ivx) {
749 if (
vtx[ivx].ID < 0)
continue;
754 vcol->emplace_back(std::move(vertex));
758 double orphanLen = 0;
759 for (ipl = 0; ipl < nplanes; ++ipl) {
760 for (icl = 0; icl <
cls[ipl].size(); ++icl) {
761 if (
cls[ipl][icl].
Length > 40 &&
cls[ipl][icl].InTrack < 0) {
762 orphanLen +=
cls[ipl][icl].Length;
765 <<
"Orphan long cluster " << ipl <<
":" << icl <<
":" <<
cls[ipl][icl].Wire[0] <<
":" 766 << (int)
cls[ipl][icl].Time[0] <<
" length " <<
cls[ipl][icl].
Length;
773 std::cout <<
"Total orphan length " << orphanLen <<
"\n";
779 evt.
put(std::move(pcol));
780 evt.
put(std::move(ptassn));
781 evt.
put(std::move(pcassn));
782 evt.
put(std::move(pvassn));
783 evt.
put(std::move(psassn));
784 evt.
put(std::move(tcol));
785 evt.
put(std::move(thassn));
786 evt.
put(std::move(scol));
787 evt.
put(std::move(vcol));
802 std::vector<std::vector<geo::WireID>> hitWID;
803 std::vector<std::vector<double>> hitX;
804 std::vector<std::vector<double>> hitXErr;
805 TVector3 pos, posErr;
806 std::vector<TVector3> trkDir;
807 std::vector<TVector3> trkDirErr;
812 unsigned short indx, indx2, iht, nHitsFit;
814 for (
unsigned short ivx = 0; ivx <
vtx.size(); ++ivx) {
815 if (!
vtx[ivx].Neutrino)
continue;
821 unsigned int thePln, theTPC, theCst;
822 for (
unsigned short itk = 0; itk <
trk.size(); ++itk) {
823 for (
unsigned short end = 0;
end < 2; ++
end) {
824 if (
trk[itk].VtxIndex[
end] != ivx)
continue;
825 unsigned short itj = 0;
826 if (
end == 1) itj =
trk[itk].TrjPos.size() - 1;
829 hitWID.resize(indx + 1);
830 hitX.resize(indx + 1);
831 hitXErr.resize(indx + 1);
832 trkDir.resize(indx + 1);
833 trkDir[indx] =
trk[itk].TrjDir[itj];
835 for (
unsigned short ipl = 0; ipl < nplanes; ++ipl) {
836 if (
trk[itk].TrkHits[ipl].
size() < 2)
continue;
838 nHitsFit =
trk[itk].TrkHits[ipl].size();
840 indx2 = hitWID[indx].size();
841 hitWID[indx].resize(indx2 + nHitsFit);
842 hitX[indx].resize(indx2 + nHitsFit);
843 hitXErr[indx].resize(indx2 + nHitsFit);
844 for (
unsigned short ii = 0; ii < nHitsFit; ++ii) {
845 if (
end == 0) { iht = ii; }
847 iht =
trk[itk].TrkHits[ipl].size() - ii - 1;
849 hitWID[indx][indx2 + ii] =
trk[itk].TrkHits[ipl][iht]->WireID();
850 thePln =
trk[itk].TrkHits[ipl][iht]->WireID().Plane;
851 theTPC =
trk[itk].TrkHits[ipl][iht]->WireID().TPC;
852 theCst =
trk[itk].TrkHits[ipl][iht]->WireID().Cryostat;
854 trk[itk].TrkHits[ipl][iht]->PeakTime(), thePln, theTPC, theCst);
855 hitXErr[indx][indx2 + ii] =
fHitFitErrFac *
trk[itk].TrkHits[ipl][iht]->RMS();
860 if (hitX.size() < 2) {
868 if (ChiDOF > 3000)
continue;
874 unsigned short fitTrk = 0;
875 for (
unsigned short itk = 0; itk <
trk.size(); ++itk) {
876 for (
unsigned short end = 0;
end < 2; ++
end) {
877 if (
trk[itk].VtxIndex[
end] != ivx)
continue;
878 unsigned short itj = 0;
879 if (
end == 1) itj =
trk[itk].TrjPos.size() - 1;
880 trk[itk].TrjDir[itj] = trkDir[fitTrk];
893 unsigned short wire, nwires, indx;
894 float dir, ctime, cx, chgWinLo, chgWinHi;
898 auto const ipl = planeid.Plane;
899 for (
unsigned short icl = 0; icl <
cls[ipl].size(); ++icl) {
901 nwires =
cls[ipl][icl].Length / 2;
902 if (nwires < 2)
continue;
904 if (nwires > 30) nwires = 30;
905 for (
unsigned short end = 0;
end < 2; ++
end) {
909 for (
unsigned short w = 0;
w < nwires; ++
w) {
910 wire =
cls[ipl][icl].Wire[
end] + dir *
w;
919 for (
unsigned int hit = firhit;
hit < lashit; ++
hit) {
922 <<
"FillChgNear bad lashit " << lashit <<
" size " <<
allhits.size() <<
"\n";
925 if (
allhits[
hit]->PeakTime() < chgWinLo)
continue;
926 if (
allhits[
hit]->PeakTime() > chgWinHi)
continue;
930 cnear /= (float)(nwires - 1);
931 if (cnear >
cls[ipl][icl].Charge[
end]) {
935 cls[ipl][icl].ChgNear[
end] = 0;
948 unsigned short ivx, itr, ipl, ii, jtr;
949 unsigned short nus, nds, nuhs, ndhs;
950 float longUSTrk, longDSTrk, qual;
954 float tgMuonCut2 = 9;
958 unsigned short VtxIndex;
959 unsigned short nUSTk;
960 unsigned short nDSTk;
961 unsigned short nUSHit;
962 unsigned short nDSHit;
967 std::vector<NuVtx> nuVtxCand;
972 float best = 999, dx, dy, dz, dr;
977 for (ivx = 0; ivx <
vtx.size(); ++ivx) {
978 vtx[ivx].Neutrino =
false;
987 for (itr = 0; itr <
trk.size(); ++itr) {
988 if (
trk[itr].ID < 0)
continue;
989 if (
trk[itr].PDGCode != 13)
continue;
990 for (itj = 0; itj <
trk[itr].TrjPos.size(); ++itj) {
991 dx =
trk[itr].TrjPos[itj](0) -
vtx[ivx].
X;
992 dy =
trk[itr].TrjPos[itj](1) -
vtx[ivx].
Y;
993 dz =
trk[itr].TrjPos[itj](2) -
vtx[ivx].
Z;
994 dr = dx * dx + dy * dy + dz * dz;
995 if (dr < tgMuonCut2) {
1003 if (skipVtx)
continue;
1004 for (itr = 0; itr <
trk.size(); ++itr) {
1005 if (
trk[itr].ID < 0)
continue;
1006 if (
trk[itr].VtxIndex[0] == ivx) {
1009 if (
trk[itr].
Length > longDSTrk) longDSTrk =
trk[itr].Length;
1010 for (ipl = 0; ipl < nplanes; ++ipl)
1011 ndhs +=
trk[itr].TrkHits[ipl].
size();
1015 if (
trk[itr].VtxIndex[1] == ivx &&
trk[itr].VtxIndex[0] >= 0) {
1019 if (
trk[itr].VtxIndex[1] == ivx &&
trk[itr].VtxIndex[0] < 0) {
1022 if (
trk[itr].
Length > longUSTrk) longUSTrk =
trk[itr].Length;
1023 for (ipl = 0; ipl < nplanes; ++ipl)
1024 nuhs +=
trk[itr].TrkHits[ipl].
size();
1027 if (skipVtx)
continue;
1028 if (nds == 0)
continue;
1029 qual = 1 / (float)nds;
1030 qual /= (float)ndhs;
1031 if (nus > 0) qual *= (float)nuhs / (
float)ndhs;
1036 if (nds > 0 && longDSTrk > 5) {
1038 aNuVtx.VtxIndex = ivx;
1041 aNuVtx.nUSHit = nuhs;
1042 aNuVtx.nDSHit = ndhs;
1043 aNuVtx.longUSTrk = longUSTrk;
1044 aNuVtx.longDSTrk = longDSTrk;
1046 nuVtxCand.push_back(aNuVtx);
1049 if (imbest < 0)
return;
1053 vtx[ivx].Neutrino =
true;
1061 std::vector<unsigned short> dtrGen;
1062 std::vector<unsigned short> dtrNextGen;
1063 for (itr = 0; itr <
trk.size(); ++itr) {
1064 if (
trk[itr].ID < 0)
continue;
1065 if (
trk[itr].VtxIndex[0] == ivx) {
1069 trk[itr].PDGCode = 2212;
1071 dtrGen.push_back(itr);
1073 if (
trk[itr].VtxIndex[1] == ivx) {
1077 trk[itr].PDGCode = 2212;
1080 std::reverse(
trk[itr].TrjPos.begin(),
trk[itr].TrjPos.end());
1081 for (ii = 0; ii <
trk[itr].TrjDir.size(); ++ii)
1082 trk[itr].TrjDir[ii] = -
trk[itr].TrjDir[ii];
1086 if (dtrGen.size() == 0)
return;
1088 unsigned short tmp, indx;
1089 unsigned short nit = 0;
1097 for (ii = 0; ii < dtrGen.size(); ++ii) {
1099 if (
trk[itr].VtxIndex[1] >= 0) {
1101 ivx =
trk[itr].VtxIndex[1];
1103 for (jtr = 0; jtr <
trk.size(); ++jtr) {
1104 if (jtr == itr)
continue;
1105 if (
trk[jtr].VtxIndex[0] == ivx) {
1107 indx =
trk[itr].DtrID.size();
1108 trk[itr].DtrID.resize(indx + 1);
1109 trk[itr].DtrID[indx] = jtr;
1110 trk[jtr].MomID =
trk[itr].ID;
1112 trk[jtr].PDGCode = 211;
1113 dtrNextGen.push_back(jtr);
1116 if (
trk[jtr].VtxIndex[1] == ivx) {
1118 indx =
trk[itr].DtrID.size();
1119 trk[itr].DtrID.resize(indx + 1);
1120 trk[itr].DtrID[indx] = jtr;
1121 trk[jtr].MomID =
trk[itr].ID;
1123 trk[jtr].PDGCode = 211;
1124 dtrNextGen.push_back(jtr);
1127 std::reverse(
trk[jtr].TrjPos.begin(),
trk[jtr].TrjPos.end());
1128 for (
unsigned short jj = 0; jj <
trk[jtr].TrjDir.size(); ++jj)
1129 trk[jtr].TrjDir[jj] = -
trk[jtr].TrjDir[jj];
1131 tmp =
trk[jtr].VtxIndex[0];
1132 trk[jtr].VtxIndex[0] =
trk[jtr].VtxIndex[1];
1133 trk[jtr].VtxIndex[1] =
tmp;
1139 if (dtrNextGen.size() == 0)
break;
1140 dtrGen = dtrNextGen;
1149 unsigned short ipf, itj;
1163 bool startsIn, endsIn;
1165 for (
unsigned short itk = 0; itk <
trk.size(); ++itk) {
1167 if (
trk[itk].ID < 0)
continue;
1172 for (ipf = 0; ipf <
pfpToTrkID.size(); ++ipf) {
1178 if (skipit)
continue;
1180 if (
trk[itk].TrjPos[0](0) < XLo ||
trk[itk].TrjPos[0](0) > XHi) startsIn =
false;
1181 if (
trk[itk].TrjPos[0](1) < YLo ||
trk[itk].TrjPos[0](1) > YHi) startsIn =
false;
1182 if (
trk[itk].TrjPos[0](2) < ZLo ||
trk[itk].TrjPos[0](2) > ZHi) startsIn =
false;
1183 if (startsIn)
continue;
1185 itj =
trk[itk].TrjPos.size() - 1;
1186 if (
trk[itk].TrjPos[itj](0) < XLo ||
trk[itk].TrjPos[itj](0) > XHi) endsIn =
false;
1187 if (
trk[itk].TrjPos[itj](1) < YLo ||
trk[itk].TrjPos[itj](1) > YHi) endsIn =
false;
1188 if (
trk[itk].TrjPos[itj](2) < ZLo ||
trk[itk].TrjPos[itj](2) > ZHi) endsIn =
false;
1189 if (endsIn)
continue;
1191 trk[itk].PDGCode = 13;
1197 for (
unsigned short itk = 0; itk <
trk.size(); ++itk) {
1199 if (
trk[itk].PDGCode != 13)
continue;
1211 unsigned short ivx, ii, ipl, icl, jj, jpl, jcl, kk, kpl, kcl;
1212 short idir, iend, jdir, jend, kdir, kend, ioend;
1215 for (ivx = 0; ivx <
vtx.size(); ++ivx) {
1218 for (ipl = 0; ipl < nplanes; ++ipl) {
1220 for (icl = 0; icl <
clsChain[ipl].size(); ++icl) {
1221 if (
clsChain[ipl][icl].InTrack >= 0)
continue;
1222 for (iend = 0; iend < 2; ++iend) {
1223 if (
clsChain[ipl][icl].VtxIndex[iend] ==
vtx[ivx].EvtIndex)
vxCls[ipl].push_back(icl);
1230 myprt <<
"VtxMatch: Vertex ID " <<
vtx[ivx].EvtIndex <<
"\n";
1231 for (ipl = 0; ipl < nplanes; ++ipl) {
1232 myprt <<
"ipl " << ipl <<
" cls";
1233 for (
unsigned short ii = 0; ii <
vxCls[ipl].size(); ++ii)
1234 myprt <<
" " <<
vxCls[ipl][ii];
1243 for (ipl = 0; ipl < nplanes; ++ipl) {
1244 if (nplanes == 2 && ipl > 0)
continue;
1245 for (ii = 0; ii <
vxCls[ipl].size(); ++ii) {
1246 icl =
vxCls[ipl][ii];
1248 if (
clsChain[ipl][icl].InTrack >= 0)
continue;
1249 jpl = (ipl + 1) % nplanes;
1250 kpl = (jpl + 1) % nplanes;
1251 for (jj = 0; jj <
vxCls[jpl].size(); ++jj) {
1252 jcl =
vxCls[jpl][jj];
1253 if (
clsChain[jpl][jcl].InTrack >= 0)
continue;
1254 for (iend = 0; iend < 2; ++iend) {
1255 if (
clsChain[ipl][icl].VtxIndex[iend] !=
vtx[ivx].EvtIndex)
continue;
1257 idir =
clsChain[ipl][icl].Dir[iend];
1258 for (jend = 0; jend < 2; ++jend) {
1259 if (
clsChain[jpl][jcl].VtxIndex[jend] !=
vtx[ivx].EvtIndex)
continue;
1260 jdir =
clsChain[jpl][jcl].Dir[jend];
1261 if (idir != 0 && jdir != 0 && idir != jdir)
continue;
1266 match.
Cls[ipl] = icl;
1267 match.
End[ipl] = iend;
1268 match.
Cls[jpl] = jcl;
1269 match.
End[jpl] = jend;
1279 <<
"FillEndMatch2: Err " << match.
Err <<
" oErr " << match.
oErr;
1280 if (match.
Err + match.
oErr > 100)
continue;
1281 if (
DupMatch(match, nplanes))
continue;
1285 match.
Cls[kpl] = -1;
1289 <<
"VtxMatch: check " << ipl <<
":" << icl <<
":" << iend <<
" and " << jpl
1290 <<
":" << jcl <<
":" << jend <<
" for cluster in kpl " << kpl;
1292 for (kk = 0; kk <
vxCls[kpl].size(); ++kk) {
1293 kcl =
vxCls[kpl][kk];
1294 if (
clsChain[kpl][kcl].InTrack >= 0)
continue;
1295 for (kend = 0; kend < 2; ++kend) {
1296 kdir =
clsChain[kpl][kcl].Dir[kend];
1297 if (idir != 0 && kdir != 0 && idir != kdir)
continue;
1298 if (
clsChain[kpl][kcl].VtxIndex[kend] !=
vtx[ivx].EvtIndex)
continue;
1301 match.
Cls[kpl] = kcl;
1302 match.
End[kpl] = kend;
1304 if (
DupMatch(match, nplanes))
continue;
1307 if (match.
Chg[kpl] <= 0)
continue;
1308 if (match.
Err + match.
oErr > 100)
continue;
1310 if (
DupMatch(match, nplanes))
continue;
1316 if (gotkcl)
continue;
1320 unsigned short kbend = 0;
1322 mf::LogVerbatim(
"CCTM") <<
"VtxMatch: look for missed cluster chain in kpl";
1323 for (kcl = 0; kcl <
clsChain[kpl].size(); ++kcl) {
1324 if (
clsChain[kpl][kcl].InTrack >= 0)
continue;
1325 for (kend = 0; kend < 2; ++kend) {
1326 kdir =
clsChain[kpl][kcl].Dir[kend];
1327 if (idir != 0 && kdir != 0 && idir != kdir)
continue;
1328 if (
clsChain[kpl][kcl].VtxIndex[kend] >= 0)
continue;
1330 if (fabs(
clsChain[kpl][kcl].
X[kend] -
vtx[ivx].
X) > 5)
continue;
1332 if (fabs(
clsChain[kpl][kcl].X[1 - kend] -
clsChain[ipl][icl].X[ioend]) > 50)
1335 match.
Cls[kpl] = kcl;
1336 match.
End[kpl] = kend;
1337 if (
DupMatch(match, nplanes))
continue;
1339 totErr = match.
Err + match.
oErr;
1342 myprt <<
"VtxMatch: Chk missing cluster match ";
1343 for (
unsigned short ii = 0; ii < nplanes; ++ii)
1344 myprt <<
" " << ii <<
":" << match.
Cls[ii] <<
":" << match.
End[ii];
1345 myprt <<
" Err " << match.
Err <<
"\n";
1347 if (totErr > 100)
continue;
1348 if (totErr < best) {
1357 match.
Cls[kpl] = kbst;
1358 match.
End[kpl] = kbend;
1362 clsChain[kpl][kbst].VtxIndex[kbend] = ivx;
1364 vxCls[kpl].push_back(kbst);
1368 match.
Cls[kpl] = -1;
1370 if (
DupMatch(match, nplanes))
continue;
1380 if (
matcomb.size() == 0)
continue;
1385 for (ipl = 0; ipl < 3; ++ipl)
1397 if (
vtx.empty())
return;
1400 unsigned int const ipl = planeID.Cryostat;
1401 for (
unsigned int icl = 0; icl <
cls[ipl].size(); ++icl) {
1402 for (
unsigned int end = 0;
end < 2u; ++
end) {
1404 if (
cls[ipl][icl].VtxIndex[
end] >= 0)
continue;
1408 unsigned int oend = 1 -
end;
1409 for (std::size_t ivx = 0; ivx <
vtx.size(); ++ivx) {
1411 if (
cls[ipl][icl].VtxIndex[oend] == static_cast<short>(ivx))
continue;
1416 if (dWire > 30 || dWire < -2)
continue;
1419 if (dWire < -30 || dWire > 2)
continue;
1422 float const dX = fabs(
cls[ipl][icl].
X[
end] +
1431 cls[ipl][icl].VtxIndex[
end] = ibstvx;
1432 cls[ipl][icl].mVtxIndex[
end] = ibstvx;
1445 unsigned short icl, icl1, icl2;
1446 float dw, dx, dWCut, dw1Max, dw2Max;
1447 float dA, dA2, dACut =
fMaxDAng, chgAsymCut;
1448 float dXCut, chgasym, mrgErr;
1451 bool gotprt =
false;
1453 for (
unsigned short ipl = 0; ipl < nplanes; ++ipl) {
1455 for (icl1 = 0; icl1 <
cls[ipl].size() - 1; ++icl1) {
1457 if (
prt) gotprt =
true;
1459 dw1Max = 0.6 *
cls[ipl][icl1].Length;
1460 ls1 = (
cls[ipl][icl1].Length > 100 &&
1461 fabs(
cls[ipl][icl1].Angle[0] -
cls[ipl][icl1].Angle[1]) < 0.04);
1462 for (icl2 = icl1 + 1; icl2 <
cls[ipl].size(); ++icl2) {
1463 ls2 = (
cls[ipl][icl2].Length > 100 &&
1464 fabs(
cls[ipl][icl2].Angle[0] -
cls[ipl][icl2].Angle[1]) < 0.04);
1465 dw2Max = 0.6 *
cls[ipl][icl2].Length;
1468 if (dw2Max < dWCut) dWCut = dw2Max;
1470 if (dWCut > 100) dWCut = 100;
1471 if (dWCut < 2) dWCut = 2;
1477 <<
"MCC P:C:W icl1 " << ipl <<
":" << icl1 <<
":" <<
cls[ipl][icl1].Wire[1]
1478 <<
" vtx " <<
cls[ipl][icl1].VtxIndex[1] <<
" ls1 " << ls1 <<
" icl2 " << ipl <<
":" 1479 << icl2 <<
":" <<
cls[ipl][icl2].Wire[0] <<
" vtx " <<
cls[ipl][icl2].VtxIndex[0]
1480 <<
" ls2 " << ls2 <<
" dWCut " << dWCut;
1481 if (
std::abs(
cls[ipl][icl1].Wire[1] -
cls[ipl][icl2].Wire[0]) > dWCut)
continue;
1488 dA = fabs(
cls[ipl][icl1].Angle[1] -
cls[ipl][icl2].Angle[0]);
1491 dA2 = fabs(
cls[ipl][icl1].Angle[0] -
cls[ipl][icl2].Angle[1]);
1495 <<
" dA " << dA <<
" dA2 " << dA2 <<
" DACut " << dACut <<
" dXCut " << dXCut;
1497 if (dA2 < dA) dA = dA2;
1503 cls[ipl][icl2].X[0];
1504 unsigned short ivx =
cls[ipl][icl1].VtxIndex[1];
1505 if (
vtx[ivx].nClusInPln[ipl] == 2 && fabs(dx) < 1) {
1506 cls[ipl][icl1].VtxIndex[1] = -2;
1507 cls[ipl][icl2].VtxIndex[0] = -2;
1508 vtx[ivx].nClusInPln[ipl] = 0;
1514 if (
cls[ipl][icl1].VtxIndex[1] >= 0)
continue;
1515 if (
cls[ipl][icl2].VtxIndex[0] >= 0)
continue;
1518 if (
cls[ipl][icl2].Wire[0] -
cls[ipl][icl1].Wire[1] < 3 &&
1529 dx =
cls[ipl][icl2].X[0] -
cls[ipl][icl1].X[1];
1530 float dA3 =
std::abs(atan(dx / dw) -
cls[ipl][icl1].Angle[1]);
1532 if (dA3 > dA) dA = dA3;
1536 if (dA > dACut)
continue;
1540 <<
" rough dX " << fabs(
cls[ipl][icl1].
X[1] -
cls[ipl][icl2].
X[0]) <<
" cut = 20";
1543 if (fabs(
cls[ipl][icl1].X[1] -
cls[ipl][icl2].X[0]) > 20)
continue;
1551 chgasym = fabs(
cls[ipl][icl1].Charge[1] -
cls[ipl][icl2].Charge[0]);
1552 chgasym /=
cls[ipl][icl1].Charge[1] +
cls[ipl][icl2].Charge[0];
1554 if (chgasym > chgAsymCut)
continue;
1560 cls[ipl][icl2].X[0];
1565 cls[ipl][icl1].X[1];
1569 if (dA2 < 0.01 &&
abs(dx) > dXCut && dx < -1) {
1570 dx =
dXClTraj(fmCluHits, ipl, icl1, 1);
1576 <<
" X0 " <<
cls[ipl][icl1].X[1] <<
" slp " <<
cls[ipl][icl1].Slope[1] <<
" dw " 1577 << dw <<
" oX " <<
cls[ipl][icl2].X[0] <<
" dx " << dx <<
" cut " << dXCut;
1579 if (fabs(dx) > dXCut)
continue;
1582 float xerr = dx / dXCut;
1583 float aerr = dA / dACut;
1584 mrgErr = xerr * xerr + aerr * aerr;
1588 <<
"icl1 mrgErr " << mrgErr <<
" MergeError " <<
cls[ipl][icl1].MergeError[1]
1589 <<
" icl2 MergeError " <<
cls[ipl][icl2].MergeError[0];
1592 if (mrgErr >
cls[ipl][icl1].MergeError[1])
continue;
1593 if (mrgErr >
cls[ipl][icl2].MergeError[0])
continue;
1596 if (
cls[ipl][icl1].BrkIndex[1] >= 0) {
1597 unsigned short ocl =
cls[ipl][icl1].BrkIndex[1];
1599 if (
cls[ipl][ocl].BrkIndex[0] == icl1) {
1600 cls[ipl][ocl].BrkIndex[0] = -1;
1603 if (
cls[ipl][ocl].BrkIndex[1] == icl1) {
1604 cls[ipl][ocl].BrkIndex[1] = -1;
1608 cls[ipl][icl1].BrkIndex[1] = icl2;
1609 cls[ipl][icl1].MergeError[1] = mrgErr;
1612 if (
cls[ipl][icl2].BrkIndex[0] >= 0) {
1613 unsigned short ocl =
cls[ipl][icl2].BrkIndex[0];
1615 if (
cls[ipl][ocl].BrkIndex[0] == icl1) {
1616 cls[ipl][ocl].BrkIndex[0] = -1;
1619 if (
cls[ipl][ocl].BrkIndex[1] == icl1) {
1620 cls[ipl][ocl].BrkIndex[1] = -1;
1624 cls[ipl][icl2].BrkIndex[0] = icl1;
1625 cls[ipl][icl2].MergeError[0] = mrgErr;
1634 for (icl1 = 0; icl1 <
cls[ipl].size() - 1; ++icl1) {
1636 for (icl2 = icl1 + 1; icl2 <
cls[ipl].size(); ++icl2) {
1638 for (
unsigned short end = 0;
end < 2; ++
end) {
1640 if (
cls[ipl][icl1].BrkIndex[
end] >= 0)
continue;
1641 if (
cls[ipl][icl2].BrkIndex[
end] >= 0)
continue;
1643 if (fabs(
cls[ipl][icl1].Angle[
end]) < 1)
continue;
1645 if (fabs(
cls[ipl][icl2].Angle[end]) < 1)
continue;
1648 <<
"BrokenC: clusters " <<
cls[ipl][icl1].Wire[
end] <<
":" 1649 << (int)
cls[ipl][icl1].Time[end] <<
" " <<
cls[ipl][icl2].Wire[end] <<
":" 1650 << (
int)
cls[ipl][icl2].Time[
end] <<
" angles " <<
cls[ipl][icl1].Angle[
end] <<
" " 1651 <<
cls[ipl][icl2].Angle[
end] <<
" dWire " 1652 << fabs(
cls[ipl][icl1].Wire[end] -
cls[ipl][icl2].Wire[end]);
1653 if (fabs(
cls[ipl][icl1].Wire[end] -
cls[ipl][icl2].Wire[end]) > 5)
continue;
1656 float dsl =
cls[ipl][icl2].Slope[
end] -
cls[ipl][icl1].Slope[
end];
1662 if (fabs(
cls[ipl][icl1].Wire[end] - fvw) > 4)
continue;
1663 if (fabs(
cls[ipl][icl2].Wire[end] - fvw) > 4)
continue;
1664 cls[ipl][icl1].BrkIndex[
end] = icl2;
1666 cls[ipl][icl1].MergeError[
end] = 1;
1667 cls[ipl][icl2].BrkIndex[
end] = icl1;
1668 cls[ipl][icl2].MergeError[
end] = 1;
1670 dx = fabs(
cls[ipl][icl1].
X[end] -
cls[ipl][icl2].
X[end]);
1673 <<
"BrokenC: icl1:W " << icl1 <<
":" <<
cls[ipl][icl1].Wire[
end] <<
" icl2:W " 1674 << icl2 <<
":" <<
cls[ipl][icl2].Wire[
end] <<
" end " << end <<
" dx " << dx;
1683 unsigned short end, mom, momBrkEnd, dtrBrkEnd, nit;
1686 std::vector<bool> gotcl(
cls[ipl].
size());
1687 for (icl = 0; icl <
cls[ipl].size(); ++icl)
1690 mf::LogVerbatim(
"CCTM") <<
"ipl " << ipl <<
" cls.size() " <<
cls[ipl].size() <<
"\n";
1692 std::vector<unsigned short> sCluster;
1693 std::vector<unsigned short> sOrder;
1694 for (icl = 0; icl <
cls[ipl].size(); ++icl) {
1697 if (gotcl[icl])
continue;
1699 if (
cls[ipl][icl].BrkIndex[0] >= 0 &&
cls[ipl][icl].BrkIndex[1] >= 0)
continue;
1700 for (end = 0; end < 2; ++
end) {
1701 if (
cls[ipl][icl].BrkIndex[end] < 0)
continue;
1707 sCluster.push_back(mom);
1708 if (momBrkEnd == 1) {
1710 sOrder.push_back(0);
1714 sOrder.push_back(1);
1716 dtr =
cls[ipl][icl].BrkIndex[
end];
1718 while (dtr >= 0 && dtr < (
short)
cls[ipl].
size() && nit <
cls[ipl].
size()) {
1720 for (dtrBrkEnd = 0; dtrBrkEnd < 2; ++dtrBrkEnd)
1721 if (
cls[ipl][dtr].BrkIndex[dtrBrkEnd] == mom)
break;
1722 if (dtrBrkEnd == 2) {
1728 sCluster.push_back(dtr);
1729 sOrder.push_back(dtrBrkEnd);
1736 momBrkEnd = 1 - dtrBrkEnd;
1738 dtr =
cls[ipl][mom].BrkIndex[momBrkEnd];
1740 if (dtrBrkEnd == 2)
continue;
1745 sCluster.push_back(icl);
1746 sOrder.push_back(0);
1749 if (sCluster.size() == 0) {
1750 mf::LogError(
"CCTM") <<
"MakeClusterChains error in plane " << ipl <<
" cluster " << icl;
1756 unsigned short jcl = sCluster[0];
1757 if (jcl >
cls[ipl].
size()) std::cout <<
"oops MCC\n";
1758 unsigned short oend;
1759 for (end = 0; end < 2; ++
end) {
1761 if (sOrder[0] > 0) oend = 1 -
end;
1764 ccp.
X[
end] =
cls[ipl][jcl].X[oend];
1776 for (
unsigned short ii = 1; ii < sCluster.size(); ++ii) {
1778 if (jcl >
cls[ipl].
size()) std::cout <<
"oops MCC\n";
1781 if (end > 1) std::cout <<
"oops2 MCC\n";
1784 ccp.
Wire[1] =
cls[ipl][jcl].Wire[oend];
1785 ccp.
Time[1] =
cls[ipl][jcl].Time[oend];
1786 ccp.
X[1] =
cls[ipl][jcl].X[oend];
1787 ccp.
Slope[1] =
cls[ipl][jcl].Slope[oend];
1788 ccp.
Angle[1] =
cls[ipl][jcl].Angle[oend];
1789 ccp.
Dir[1] =
cls[ipl][jcl].Dir[oend];
1791 ccp.
ChgNear[1] =
cls[ipl][jcl].ChgNear[oend];
1812 unsigned short brkCls;
1814 for (
unsigned short ccl = 0; ccl <
clsChain[ipl].size(); ++ccl) {
1815 for (
unsigned short end = 0; end < 2; ++
end) {
1816 if (
clsChain[ipl][ccl].mBrkIndex[end] < 0)
continue;
1820 for (
unsigned short ccl2 = 0; ccl2 <
clsChain[ipl].size(); ++ccl2) {
1821 if (ccl2 == ccl)
continue;
1822 if (std::find(
clsChain[ipl][ccl2].ClsIndex.begin(),
1823 clsChain[ipl][ccl2].ClsIndex.end(),
1824 brkCls) ==
clsChain[ipl][ccl2].ClsIndex.end())
1832 mf::LogError(
"CCTM") <<
"MCC: Cluster chain " << ccl <<
" end " << end
1833 <<
" Failed to find brkCls " << brkCls <<
" in plane " << ipl;
1847 unsigned short icl1,
1848 unsigned short end1)
1851 float dw, dx, best = 999;
1852 std::vector<art::Ptr<recob::Hit>> clusterhits = fmCluHits.at(
cls[ipl][icl1].EvtIndex);
1853 for (
unsigned short hit = 0;
hit < clusterhits.size(); ++
hit) {
1854 dw = clusterhits[
hit]->WireID().Wire -
cls[ipl][icl1].Wire[end1];
1855 dx = fabs(
cls[ipl][icl1].Time[end1] + dw *
fWirePitch *
cls[ipl][icl1].Slope[end1] -
1856 clusterhits[
hit]->PeakTime());
1857 if (dx < best) best = dx;
1858 if (dx < 0.01)
break;
1866 unsigned short imat,
1867 unsigned short procCode,
1874 if (imat >
matcomb.size() - 1) {
1880 unsigned short nhitinpl = 0;
1882 for (
unsigned short ipl = 0; ipl < nplanes; ++ipl)
1885 mf::LogError(
"CCTM") <<
"StoreTrack: Not enough hits in each plane\n";
1889 mf::LogVerbatim(
"CCTM") <<
"In StoreTrack: matcomb " << imat <<
" cluster chains " 1894 std::array<std::vector<geo::WireID>, 3> trkWID;
1895 std::array<std::vector<double>, 3> trkX;
1896 std::array<std::vector<double>, 3> trkXErr;
1899 std::vector<TVector3> trkPos;
1900 std::vector<TVector3> trkDir;
1902 newtrk.
ID =
trk.size() + 1;
1903 newtrk.
Proc = procCode;
1912 newtrk.
EndInTPC = {{
false,
false}};
1913 newtrk.
GoodEnd = {{
false,
false}};
1917 unsigned short icl, iht;
1920 mf::LogVerbatim(
"CCTM") <<
"CCTM: Make traj for track " << newtrk.
ID <<
" procCode " 1921 << procCode <<
" nhits in planes " <<
trkHits[0].size() <<
" " 1925 trkWID[2].resize(0);
1927 trkXErr[2].resize(0);
1930 auto const ipl = planeid.Plane;
1934 for (iht = 0; iht <
trkHits[ipl].size(); ++iht) {
1935 trkWID[ipl][iht] =
trkHits[ipl][iht]->WireID();
1942 if (trkPos.size() < 2) {
1943 mf::LogError(
"CCTM") <<
"StoreTrack: No trajectory points on failed track " << newtrk.
ID 1944 <<
" in StoreTrack: matcomb " << imat <<
" cluster chains " 1960 unsigned short end, nClose, indx, jndx;
1962 for (end = 0; end < 2; ++
end) {
1964 for (
unsigned short ipl = 0; ipl < nplanes - 1; ++ipl) {
1965 if (trkX[ipl].
size() == 0)
continue;
1966 for (
unsigned short jpl = ipl + 1; jpl < nplanes; ++jpl) {
1967 if (trkX[jpl].
size() == 0)
continue;
1973 indx = trkXErr[ipl].size() - 1;
1974 jndx = trkXErr[jpl].size() - 1;
1976 xErr = 3 * (trkXErr[ipl][indx] + trkXErr[jpl][jndx]);
1977 if (
std::abs(trkX[ipl][indx] - trkX[jpl][jndx]) < xErr) ++nClose;
1980 if (nClose == nplanes) newtrk.
GoodEnd[
end] =
true;
1984 unsigned short ivx, itj, ccl;
1985 float dx, dy, dz, dr0, dr1;
1986 unsigned short attachEnd;
1987 for (end = 0; end < 2; ++
end) {
1990 if (end == 1 &&
matcomb[imat].oVtx >= 0) ivx =
matcomb[imat].oVtx;
1991 if (ivx == USHRT_MAX)
continue;
1994 dx =
vtx[ivx].X - newtrk.
TrjPos[itj](0);
1995 dy =
vtx[ivx].Y - newtrk.
TrjPos[itj](1);
1996 dz =
vtx[ivx].Z - newtrk.
TrjPos[itj](2);
1997 dr0 = dx * dx + dy * dy + dz * dz;
1998 itj = newtrk.
TrjPos.size() - 1;
1999 dx =
vtx[ivx].X - newtrk.
TrjPos[itj](0);
2000 dy =
vtx[ivx].Y - newtrk.
TrjPos[itj](1);
2001 dz =
vtx[ivx].Z - newtrk.
TrjPos[itj](2);
2002 dr1 = dx * dx + dy * dy + dz * dz;
2008 if (dr0 > 5)
return;
2012 if (dr1 > 5)
return;
2022 newtrk.
TrjDir[0] = dir.Unit();
2026 newtrk.
TrjDir[itj] = dir.Unit();
2032 <<
"StoreTrack: Trying to attach a vertex to both ends of a track. imat = " << imat;
2040 for (
unsigned short itj = 1; itj < newtrk.
TrjPos.size(); ++itj) {
2044 norm = sqrt(X * X + Y * Y + Z * Z);
2050 for (
unsigned short ipl = 0; ipl < nplanes; ++ipl) {
2051 if (
matcomb[imat].Cls[ipl] < 0)
continue;
2053 if (ccl >
clsChain[ipl].
size()) std::cout <<
"oops StoreTrack\n";
2055 for (
unsigned short icc = 0; icc <
clsChain[ipl][ccl].ClsIndex.size(); ++icc) {
2056 icl =
clsChain[ipl][ccl].ClsIndex[icc];
2057 if (icl >
cls[ipl].
size()) std::cout <<
"oops StoreTrack\n";
2058 cls[ipl][icl].InTrack = newtrk.
ID;
2059 if (
cls[ipl][icl].EvtIndex > fmCluHits.size() - 1) {
2060 std::cout <<
"ooops2 store track EvtIndex " <<
cls[ipl][icl].EvtIndex <<
" size " 2061 << fmCluHits.size() <<
" icl " << icl <<
"\n";
2070 trk.push_back(newtrk);
2079 float kSlp, kAng,
kX, kWir, okWir;
2080 short idir, ioend, jdir, joend, kdir;
2095 std::array<float, 3> mchg;
2097 for (
unsigned short ipl = 0; ipl < nplanes; ++ipl) {
2099 for (
unsigned short icl = 0; icl <
clsChain[ipl].size(); ++icl) {
2100 if (
clsChain[ipl][icl].InTrack >= 0)
continue;
2103 unsigned short jpl = (ipl + 1) % nplanes;
2104 unsigned short kpl = (jpl + 1) % nplanes;
2106 for (
unsigned short jcl = 0; jcl <
clsChain[jpl].size(); ++jcl) {
2107 if (
clsChain[jpl][jcl].InTrack >= 0)
continue;
2111 mchg[0] =
clsChain[ipl][icl].TotChg;
2112 mchg[1] =
clsChain[jpl][jcl].TotChg;
2115 for (
unsigned short iend = 0; iend < 2; ++iend) {
2116 idir =
clsChain[ipl][icl].Dir[iend];
2117 for (
unsigned short jend = 0; jend < 2; ++jend) {
2118 jdir =
clsChain[jpl][jcl].Dir[jend];
2119 if (idir != 0 && jdir != 0 && idir != jdir)
continue;
2126 ipl,
clsChain[ipl][icl].Slope[iend], jpl,
clsChain[jpl][jcl].Slope[jend], tpcid);
2134 if (yp > tpcSizeY || yp < -tpcSizeY)
continue;
2135 if (zp < 0 || zp > tpcSizeZ)
continue;
2144 if (yp > tpcSizeY || yp < -tpcSizeY)
continue;
2145 if (zp < 0 || zp > tpcSizeZ)
continue;
2149 <<
"PlnMatch: chk i " << ipl <<
":" << icl <<
":" << iend <<
" idir " << idir
2150 <<
" X " <<
clsChain[ipl][icl].X[iend] <<
" j " << jpl <<
":" << jcl <<
":" 2151 << jend <<
" jdir " << jdir <<
" X " <<
clsChain[jpl][jcl].X[jend];
2155 <<
"PlnMatch: chk j " << ipl <<
":" << icl <<
":" << iend <<
" " << jpl <<
":" 2156 << jcl <<
":" << jend <<
" iSlp " << std::setprecision(2)
2157 <<
clsChain[ipl][icl].Slope[iend] <<
" jSlp " << std::setprecision(2)
2158 <<
clsChain[jpl][jcl].Slope[jend] <<
" kWir " << (int)kWir <<
" okWir " 2159 << (
int)okWir <<
" kSlp " << std::setprecision(2) << kSlp <<
" kAng " 2160 << std::setprecision(2) << kAng <<
" kX " << std::setprecision(1) <<
kX;
2164 ignoreSign = (fabs(kSlp) > 1.5);
2165 if (ignoreSign) kAng = fabs(kAng);
2167 bool gotkcl =
false;
2168 for (
unsigned short kcl = 0; kcl <
clsChain[kpl].size(); ++kcl) {
2169 if (
clsChain[kpl][kcl].InTrack >= 0)
continue;
2171 mchg[0] =
clsChain[ipl][icl].TotChg;
2172 mchg[1] =
clsChain[jpl][jcl].TotChg;
2173 mchg[2] =
clsChain[kpl][kcl].TotChg;
2175 for (
unsigned short kend = 0; kend < 2; ++kend) {
2176 kdir =
clsChain[kpl][kcl].Dir[kend];
2177 if (idir != 0 && kdir != 0 && idir != kdir)
continue;
2180 <<
" kcl " << kcl <<
" kend " << kend <<
" dx " 2186 <<
" kcl " << kcl <<
" kend " << kend <<
" dw " 2187 << (
clsChain[kpl][kcl].Wire[kend] - kWir) <<
" ignoreSign " << ignoreSign;
2188 if (fabs(
clsChain[kpl][kcl].Wire[kend] - kWir) > dwcut)
continue;
2191 match.
Cls[ipl] = icl;
2192 match.
End[ipl] = iend;
2193 match.
Cls[jpl] = jcl;
2194 match.
End[jpl] = jend;
2195 match.
Cls[kpl] = kcl;
2196 match.
End[kpl] = kend;
2198 if (
DupMatch(match, nplanes))
continue;
2207 <<
":" << match.
End[kpl] <<
" oChg " << match.
Chg[kpl]
2208 <<
" mErr " << match.
Err <<
" oErr " << match.
oErr;
2209 if (match.
Chg[kpl] == 0)
continue;
2210 if (match.
Err > 10 || match.
oErr > 10)
continue;
2212 if (
DupMatch(match, nplanes))
continue;
2221 match.
Cls[ipl] = icl;
2222 match.
End[ipl] = iend;
2223 match.
Cls[jpl] = jcl;
2224 match.
End[jpl] = jend;
2225 match.
Cls[kpl] = -1;
2228 if (
DupMatch(match, nplanes))
continue;
2237 <<
" Tried 2-plane match" 2238 <<
" k " << kpl <<
":" << match.
Cls[kpl] <<
":" << match.
End[kpl] <<
" Chg " 2239 << match.
Chg[kpl] <<
" Err " << match.
Err <<
" match.oErr " << match.
oErr;
2240 if (match.
Chg[kpl] <= 0)
continue;
2241 if (match.
Err > 10 || match.
oErr > 10)
continue;
2254 unsigned short nMatCl, nMiss;
2255 float toterr = match.
Err + match.
oErr;
2256 for (
unsigned int imat = 0; imat <
matcomb.size(); ++imat) {
2283 for (
unsigned short ipl = 0; ipl < nplanes; ++ipl) {
2284 if (match.
Cls[ipl] >= 0) {
2285 if (match.
Cls[ipl] ==
matcomb[imat].Cls[ipl] &&
2293 if (nMatCl == 2 && nMiss == 1)
return true;
2301 unsigned short const procCode,
2307 std::vector<CluLen> materr;
2308 unsigned int ii, im;
2310 if (
matcomb.size() == 0)
return;
2313 for (ii = 0; ii <
matcomb.size(); ++ii) {
2316 materr.push_back(merr);
2318 std::sort(materr.begin(), materr.end(),
lessThan);
2323 myprt <<
"SortMatches\n";
2324 myprt <<
" ii im Vx Err dW dA dX oVx oErr odW odA odX Asym " 2326 for (ii = 0; ii < materr.size(); ++ii) {
2327 im = materr[ii].index;
2332 myprt << std::fixed <<
std::right << std::setw(5) << ii << std::setw(5) << im
2333 << std::setw(4) <<
matcomb[im].Vtx << std::setw(7) << std::setprecision(2)
2334 <<
matcomb[im].Err << std::setw(7) << std::setprecision(1) <<
matcomb[im].dWir
2335 << std::setw(7) << std::setprecision(2) <<
matcomb[im].dAng << std::setw(7)
2336 << std::setprecision(2) <<
matcomb[im].dX << std::setw(4) <<
matcomb[im].oVtx
2337 << std::setw(7) << std::setprecision(2) <<
matcomb[im].oErr << std::setw(7)
2338 << std::setprecision(1) <<
matcomb[im].odWir << std::setw(7) << std::setprecision(2)
2339 <<
matcomb[im].odAng << std::setw(7) << std::setprecision(2) <<
matcomb[im].odX
2340 << std::setw(7) << std::setprecision(3) << asym <<
" 0:" <<
matcomb[im].Cls[0] <<
":" 2342 if (nplanes > 2) myprt <<
" 2:" <<
matcomb[im].Cls[2] <<
":" <<
matcomb[im].End[2];
2348 std::array<std::vector<bool>, 3> pclUsed;
2350 for (ipl = 0; ipl < nplanes; ++ipl) {
2355 unsigned short matcombTotCl = 0;
2356 float matcombTotLen = 0;
2358 for (ii = 0; ii <
matcomb.size(); ++ii) {
2359 for (ipl = 0; ipl < nplanes; ++ipl) {
2360 if (
matcomb[ii].Cls[ipl] < 0)
continue;
2363 matcombTotLen +=
clsChain[ipl][icl].Length;
2368 mf::LogVerbatim(
"CCTM") <<
"Number of clusters to match " << matcombTotCl <<
" total length " 2371 if (matcombTotLen <= 0) {
2372 mf::LogError(
"CCTM") <<
"SortMatches: bad matcomb total length " << matcombTotLen;
2377 std::vector<unsigned short> matIndex;
2379 std::vector<unsigned short> bestMatIndex;
2380 float totLen, totErr, bestTotErr = 9999;
2382 unsigned short jj, jm, nused, jcl;
2386 for (ii = 0; ii < materr.size(); ++ii) {
2387 im = materr[ii].index;
2390 if (
matcomb[im].Err > bestTotErr)
continue;
2393 for (ipl = 0; ipl < nplanes; ++ipl) {
2397 if (
matcomb[im].Cls[ipl] < 0)
continue;
2399 pclUsed[ipl][icl] =
true;
2400 totLen +=
clsChain[ipl][icl].Length;
2405 matIndex.push_back(im);
2407 for (jj = 0; jj < materr.size(); ++jj) {
2408 if (jj == ii)
continue;
2409 jm = materr[jj].index;
2411 if (
matcomb[jm].Err > bestTotErr)
continue;
2414 for (ipl = 0; ipl < nplanes; ++ipl) {
2415 if (
matcomb[jm].Cls[ipl] < 0)
continue;
2417 if (pclUsed[ipl][jcl]) ++nused;
2419 if (nused > 0)
break;
2420 totLen +=
clsChain[ipl][jcl].Length;
2423 if (nused != 0)
continue;
2426 matIndex.push_back(jm);
2428 for (ipl = 0; ipl < nplanes; ++ipl) {
2429 if (
matcomb[jm].Cls[ipl] < 0)
continue;
2431 pclUsed[ipl][jcl] =
true;
2434 if (totLen == 0)
continue;
2436 for (ipl = 0; ipl < nplanes; ++ipl) {
2437 for (
unsigned short indx = 0; indx < pclUsed[ipl].size(); ++indx)
2438 if (pclUsed[ipl][indx]) ++nused;
2440 if (totLen > matcombTotLen) std::cout <<
"Oops " << totLen <<
" " << matcombTotLen <<
"\n";
2442 fracLen = totLen / matcombTotLen;
2446 myprt <<
"match " << im <<
" totErr " << totErr <<
" nused " << nused <<
" fracLen " 2447 << fracLen <<
" totLen " << totLen <<
" mat: ";
2448 for (
unsigned short indx = 0; indx < matIndex.size(); ++indx)
2449 myprt <<
" " << matIndex[indx];
2452 if (totErr < bestTotErr) {
2453 bestTotErr = totErr;
2454 bestMatIndex = matIndex;
2455 if (nused == matcombTotCl)
break;
2458 myprt <<
"bestTotErr " << bestTotErr <<
" nused " << nused <<
" matcombTotCl " 2459 << matcombTotCl <<
" mat: ";
2460 for (
unsigned short indx = 0; indx < bestMatIndex.size(); ++indx)
2461 myprt <<
" " << bestMatIndex[indx];
2464 if (fracLen > 0.999)
break;
2468 if (bestTotErr > 9000)
return;
2470 for (ii = 0; ii < bestMatIndex.size(); ++ii) {
2471 im = bestMatIndex[ii];
2475 StoreTrack(detProp, fmCluHits, im, procCode, tpcid);
2493 unsigned short ipl = 0;
2494 unsigned short jpl = 1;
2496 if (match.
Cls[0] < 0 || match.
Cls[1] < 0)
return;
2498 unsigned short icl = match.
Cls[ipl];
2499 unsigned short iend = match.
End[ipl];
2502 float miX =
clsChain[ipl][icl].X[iend];
2504 unsigned short oiend = 1 - iend;
2505 float oiX =
clsChain[ipl][icl].X[oiend];
2507 unsigned short jcl = match.
Cls[jpl];
2508 unsigned short jend = match.
End[jpl];
2511 float mjX =
clsChain[jpl][jcl].X[jend];
2513 unsigned short ojend = 1 - jend;
2514 float ojX =
clsChain[jpl][jcl].X[ojend];
2518 if (
clsChain[ipl][icl].VtxIndex[iend] >= 0 &&
2527 if (
clsChain[ipl][icl].VtxIndex[oiend] >= 0 &&
2528 clsChain[ipl][icl].VtxIndex[oiend] ==
clsChain[jpl][jcl].VtxIndex[ojend]) {
2537 chgAsym = fabs(match.
Chg[ipl] - match.
Chg[jpl]) / (match.
Chg[ipl] + match.
Chg[jpl]);
2538 if (chgAsym > 0.5)
return;
2543 float maxSlp = fabs(
clsChain[ipl][icl].Slope[iend]);
2544 if (fabs(
clsChain[jpl][jcl].Slope[jend]) > maxSlp)
2545 maxSlp = fabs(
clsChain[jpl][jcl].Slope[jend]);
2547 match.
dX = fabs(miX - mjX);
2548 match.
Err = chgAsym * match.
dX / sigmaX;
2551 maxSlp = fabs(
clsChain[ipl][icl].Slope[oiend]);
2552 if (fabs(
clsChain[jpl][jcl].Slope[ojend]) > maxSlp)
2553 maxSlp = fabs(
clsChain[jpl][jcl].Slope[ojend]);
2555 match.
odX = fabs(oiX - ojX);
2556 match.
oErr = chgAsym * match.
odX / sigmaX;
2559 mf::LogVerbatim(
"CCTM") <<
"FEM2: m " << ipl <<
":" << icl <<
":" << iend <<
" miX " << miX
2560 <<
" - " << jpl <<
":" << jcl <<
":" << jend <<
" mjX " << mjX
2561 <<
" match.dX " << match.
dX <<
" match.Err " << match.
Err 2562 <<
" chgAsym " << chgAsym <<
" o " 2563 <<
" oiX " << oiX <<
" ojX " << ojX <<
" match.odX " << match.
odX 2564 <<
" match.oErr " << match.
oErr <<
"\n";
2586 std::array<short, 3> mVtx;
2587 std::array<short, 3> oVtx;
2588 std::array<float, 3> oWir;
2589 std::array<float, 3> oSlp;
2590 std::array<float, 3> oAng;
2591 std::array<float, 3> oX;
2593 std::array<float, 3> mChg;
2595 unsigned short ii, ipl, iend, jpl, jend, kpl, kend, oend;
2596 short icl, jcl, kcl;
2598 for (ipl = 0; ipl < 3; ++ipl) {
2622 for (ipl = 0; ipl < nplanes; ++ipl) {
2623 myprt <<
" " << ipl <<
":" << match.
Cls[ipl] <<
":" << match.
End[ipl];
2627 short missingPlane = -1;
2628 unsigned short nClInPln = 0;
2631 unsigned short novxmat = 0;
2633 unsigned short nvxmat = 0;
2634 unsigned short nShortCl = 0;
2636 for (ipl = 0; ipl < nplanes; ++ipl) {
2637 if (match.
Cls[ipl] < 0) {
2642 icl = match.
Cls[ipl];
2644 mChg[ipl] =
clsChain[ipl][icl].TotChg;
2645 iend = match.
End[ipl];
2646 mVtx[ipl] =
clsChain[ipl][icl].VtxIndex[iend];
2648 if (mVtx[ipl] >= 0) {
2649 if (aVtx < 0) aVtx = mVtx[ipl];
2650 if (mVtx[ipl] == aVtx) ++nvxmat;
2653 mf::LogVerbatim(
"CCTM") <<
"FEM: m " << ipl <<
":" << icl <<
":" << iend <<
" Vtx " 2654 << mVtx[ipl] <<
" Wir " <<
clsChain[ipl][icl].Wire[iend]
2655 << std::fixed << std::setprecision(3) <<
" Slp " 2656 <<
clsChain[ipl][icl].Slope[iend] << std::fixed
2657 << std::setprecision(1) <<
" X " <<
clsChain[ipl][icl].X[iend];
2660 oWir[ipl] =
clsChain[ipl][icl].Wire[oend];
2661 oAng[ipl] =
clsChain[ipl][icl].Angle[oend];
2662 oSlp[ipl] =
clsChain[ipl][icl].Slope[oend];
2663 oX[ipl] =
clsChain[ipl][icl].X[oend];
2664 oVtx[ipl] =
clsChain[ipl][icl].VtxIndex[oend];
2665 if (oVtx[ipl] >= 0) {
2666 if (aoVtx < 0) aoVtx = oVtx[ipl];
2667 if (oVtx[ipl] == aoVtx) ++novxmat;
2671 mf::LogVerbatim(
"CCTM") <<
" o " << ipl <<
":" << icl <<
":" << oend <<
" oVtx " 2672 << oVtx[ipl] <<
" oWir " << oWir[ipl] << std::fixed
2673 << std::setprecision(3) <<
" oSlp " << oSlp[ipl] << std::fixed
2674 << std::setprecision(1) <<
" oX " << oX[ipl] <<
" Chg " 2679 bool isShort = (nShortCl > 1);
2687 mf::LogVerbatim(
"CCTM") <<
"FEM: Vtx m " << aVtx <<
" count " << nvxmat <<
" o " << aoVtx
2688 <<
" count " << novxmat <<
" missingPlane " << missingPlane
2689 <<
" nClInPln " << nClInPln;
2692 if (nvxmat == 3 && novxmat == 3) {
2703 float ovxFactor = 1;
2704 if (nClInPln == 3) {
2737 float kSlp, okSlp, kAng, okAng, okX,
kX, kTim, okTim;
2739 if (nClInPln == 3) {
2751 else if (kpl == 1) {
2760 iend = match.
End[ipl];
2761 jend = match.
End[jpl];
2762 icl = match.
Cls[ipl];
2763 jcl = match.
Cls[jpl];
2765 kcl = match.
Cls[kpl];
2766 kend = match.
End[kpl];
2774 okAng = atan(okSlp);
2777 bool ignoreSign = (fabs(okSlp) > 10);
2778 if (ignoreSign) okAng = fabs(okAng);
2779 if (match.
oVtx >= 0) {
2789 okX = 0.5 * (oX[ipl] + oX[jpl]);
2794 <<
" kpl " << kpl <<
" okSlp " << okSlp <<
" okAng " << okAng
2795 <<
" okWir " << (int)okWir <<
" okX " << okX;
2800 ipl,
clsChain[ipl][icl].Slope[iend], jpl,
clsChain[jpl][jcl].Slope[jend], kplaneid);
2802 if (ignoreSign) kAng = fabs(kAng);
2803 if (match.
Vtx >= 0) {
2804 if (
vtx.size() == 0 || (
unsigned int)match.
Vtx >
vtx.size() - 1) {
2805 mf::LogError(
"CCTM") <<
"FEM: Bad match.Vtx " << match.
Vtx <<
" vtx size " <<
vtx.size();
2824 <<
" kpl " << kpl <<
" kSlp " << kSlp <<
" kAng " << kAng <<
" kX " 2831 match.
Cls[kpl] = kcl;
2832 match.
End[kpl] = kend;
2834 mChg[kpl] =
clsChain[kpl][kcl].TotChg;
2836 oWir[kpl] =
clsChain[kpl][kcl].Wire[oend];
2837 oX[kpl] =
clsChain[kpl][kcl].X[oend];
2838 oAng[kpl] =
clsChain[kpl][kcl].Angle[oend];
2839 oSlp[kpl] =
clsChain[kpl][kcl].Slope[oend];
2844 if (nClInPln == 2 && fabs(okWir - kWir) > 3)
return;
2850 if (nClInPln < 3 && mChg[missingPlane] <= 0) {
2851 if (missingPlane != kpl)
2852 mf::LogError(
"CCTM") <<
"FEM bad missingPlane " << missingPlane <<
" " << kpl <<
"\n";
2853 mChg[kpl] =
ChargeNear(kplaneid, (
unsigned short)kWir, kTim, (
unsigned short)okWir, okTim);
2854 match.
Chg[kpl] = mChg[kpl];
2856 mf::LogVerbatim(
"CCTM") <<
"FEM: Missing cluster in " << kpl <<
" ChargeNear " << (int)kWir
2857 <<
":" << (
int)kTim <<
" " << (int)okWir <<
":" << (
int)okTim
2858 <<
" chg " << mChg[kpl];
2859 if (mChg[kpl] <= 0)
return;
2864 if (chgAsym > 0.5)
return;
2869 float sigmaX, sigmaA;
2875 bool allPlnVtxMatch =
false;
2876 if (nClInPln == 3) {
2877 unsigned short nmvtx = 0;
2878 for (ii = 0; ii < nplanes; ++ii) {
2879 if (mVtx[ii] >= 0) {
2880 if (aVtx < 0) aVtx = mVtx[ii];
2885 if (nmvtx) allPlnVtxMatch =
true;
2895 if (nClInPln == 3) {
2896 kcl = match.
Cls[kpl];
2897 kend = match.
End[kpl];
2898 dw = kWir -
clsChain[kpl][kcl].Wire[kend];
2900 if (fabs(match.
dWir) > 100)
return;
2901 if (match.
Vtx >= 0) { match.
dX = kX -
clsChain[kpl][kcl].X[kend]; }
2909 if (ignoreSign) { match.
dAng = kAng - fabs(
clsChain[kpl][kcl].Angle[kend]); }
2914 da = fabs(match.
dAng) / sigmaA;
2915 dx = fabs(match.
dX) / sigmaX;
2916 if (allPlnVtxMatch) {
2918 match.
Err = vxFactor * chgAsym * da / 3;
2924 match.
Err = vxFactor * chgAsym * sqrt(dx * dx + da * da + dw * dw) / 9;
2934 match.
Err = 3 + vxFactor * chgAsym * fabs(match.
dX) / sigmaX;
2939 if (nClInPln == 3) {
2941 dw = okWir - oWir[kpl];
2943 if (match.
oVtx >= 0) { match.
odX = okX - oX[kpl]; }
2952 if (ignoreSign) { match.
odAng = okAng - fabs(oAng[kpl]); }
2954 match.
odAng = okAng - oAng[kpl];
2957 da = match.
odAng / sigmaA;
2958 dx = fabs(match.
odX) / sigmaX;
2962 match.
oErr = ovxFactor * chgAsym * sqrt(dx * dx + da * da + dw * dw) / 9;
2967 match.
odX = (oX[ipl] - oX[jpl]) / sigmaX;
2968 match.
oErr = 3 + ovxFactor * chgAsym * fabs(match.
odX);
2976 unsigned short& kend,
2984 unsigned short okend;
2987 if (kcl >= 0)
return false;
2990 float kslp = fabs((okX - kX) / (okWir - kWir));
2991 if (kslp > 20) kslp = 20;
2994 unsigned short nfound = 0;
2995 unsigned short foundCl = 0, foundEnd = 0;
2996 for (
unsigned short ccl = 0; ccl <
clsChain[kpl].size(); ++ccl) {
2997 if (
clsChain[kpl][ccl].InTrack >= 0)
continue;
2999 for (
unsigned short end = 0;
end < 2; ++
end) {
3001 if (fabs(
clsChain[kpl][ccl].Wire[
end] - kWir) > 4)
continue;
3002 if (fabs(
clsChain[kpl][ccl].Wire[okend] - okWir) > 4)
continue;
3005 fabs(
clsChain[kpl][ccl].
X[okend] - okX) > dxcut)
3012 if (nfound == 0)
return false;
3014 mf::LogVerbatim(
"CCTM") <<
"FindMissingCluster: Found too many matches. Write some code " 3029 float big = 0, small = 1.E9;
3030 for (
unsigned short ii = 0; ii < 3; ++ii) {
3031 if (mChg[ii] < small) small = mChg[ii];
3032 if (mChg[ii] > big) big = mChg[ii];
3035 return (big - small) / (big + small);
3040 unsigned short imat,
3041 unsigned short nplanes)
3047 for (ipl = 0; ipl < 3; ++ipl)
3050 if (imat >
matcomb.size())
return;
3052 unsigned short indx;
3053 std::vector<art::Ptr<recob::Hit>> clusterhits;
3054 unsigned short icc, ccl, icl, ecl, iht, ii;
3055 short endOrder, fillOrder;
3058 mf::LogVerbatim(
"CCTM") <<
"In FillTrkHits: matcomb " << imat <<
" cluster chains " 3062 for (ipl = 0; ipl < nplanes; ++ipl) {
3063 if (
matcomb[imat].Cls[ipl] < 0)
continue;
3067 endOrder = 1 - 2 *
matcomb[imat].End[ipl];
3070 std::reverse(
clsChain[ipl][ccl].ClsIndex.begin(),
clsChain[ipl][ccl].ClsIndex.end());
3071 std::reverse(
clsChain[ipl][ccl].Order.begin(),
clsChain[ipl][ccl].Order.end());
3072 for (ii = 0; ii <
clsChain[ipl][ccl].Order.size(); ++ii)
3076 mf::LogError(
"CCTM") <<
"Bad cluster chain index " << ccl <<
" in plane " << ipl;
3080 for (icc = 0; icc <
clsChain[ipl][ccl].ClsIndex.size(); ++icc) {
3081 icl =
clsChain[ipl][ccl].ClsIndex[icc];
3082 if (icl > fmCluHits.size() - 1) {
3083 std::cout <<
"oops in FTH " << icl <<
" clsChain size " <<
clsChain[ipl].size() <<
"\n";
3086 ecl =
cls[ipl][icl].EvtIndex;
3087 if (ecl > fmCluHits.size() - 1) {
3088 std::cout <<
"FTH bad EvtIndex " << ecl <<
" fmCluHits size " << fmCluHits.size() <<
"\n";
3091 clusterhits = fmCluHits.at(ecl);
3092 if (clusterhits.size() == 0) {
3093 std::cout <<
"FTH no cluster hits for EvtIndex " <<
cls[ipl][icl].EvtIndex <<
"\n";
3097 trkHits[ipl].resize(indx + clusterhits.size());
3099 fillOrder = 1 - 2 *
clsChain[ipl][ccl].Order[icc];
3100 if (fillOrder == 1) {
3101 for (iht = 0; iht < clusterhits.size(); ++iht) {
3102 if (indx + iht >
trkHits[ipl].
size() - 1) std::cout <<
"FTH oops3\n";
3103 trkHits[ipl][indx + iht] = clusterhits[iht];
3107 for (ii = 0; ii < clusterhits.size(); ++ii) {
3108 iht = clusterhits.size() - 1 - ii;
3109 if (indx + ii >
trkHits[ipl].
size() - 1) std::cout <<
"FTH oops4\n";
3110 trkHits[ipl][indx + ii] = clusterhits[iht];
3117 <<
" w " <<
trkHits[ipl][0]->WireID().Wire <<
":" 3118 << (int)
trkHits[ipl][0]->PeakTime() <<
" last p " 3119 <<
trkHits[ipl][ii]->WireID().Plane <<
" w " 3120 <<
trkHits[ipl][ii]->WireID().Wire <<
":" 3121 << (int)
trkHits[ipl][ii]->PeakTime();
3133 myprt <<
"********* PrintTracks \n";
3134 myprt <<
"vtx Index X Y Z\n";
3135 for (
unsigned short ivx = 0; ivx <
vtx.size(); ++ivx) {
3136 myprt <<
std::right << std::setw(4) << ivx << std::setw(4) <<
vtx[ivx].EvtIndex;
3137 myprt << std::fixed;
3138 myprt <<
std::right << std::setw(10) << std::setprecision(1) <<
vtx[ivx].X;
3139 myprt <<
std::right << std::setw(7) << std::setprecision(1) <<
vtx[ivx].Y;
3140 myprt <<
std::right << std::setw(7) << std::setprecision(1) <<
vtx[ivx].Z;
3141 if (
vtx[ivx].Neutrino) myprt <<
" Neutrino vertex";
3145 myprt <<
">>>>>>>>>> Tracks \n";
3146 myprt <<
"trk ID Proc nht nTrj sX sY sZ eX eY eZ sVx eVx sGd eGd " 3147 "ChgOrd dirZ Mom PDG ClsIndices\n";
3148 for (
unsigned short itr = 0; itr <
trk.size(); ++itr) {
3149 myprt <<
std::right << std::setw(3) << itr << std::setw(4) <<
trk[itr].ID;
3150 myprt <<
std::right << std::setw(5) << std::setw(4) <<
trk[itr].Proc;
3151 unsigned short nht = 0;
3152 for (
unsigned short ii = 0; ii < 3; ++ii)
3153 nht +=
trk[itr].TrkHits[ii].
size();
3155 myprt << std::setw(4) <<
trk[itr].TrjPos.size();
3156 myprt << std::fixed;
3157 myprt <<
std::right << std::setw(7) << std::setprecision(1) <<
trk[itr].TrjPos[0](0);
3158 myprt <<
std::right << std::setw(7) << std::setprecision(1) <<
trk[itr].TrjPos[0](1);
3159 myprt <<
std::right << std::setw(7) << std::setprecision(1) <<
trk[itr].TrjPos[0](2);
3160 unsigned short itj =
trk[itr].TrjPos.size() - 1;
3161 myprt <<
std::right << std::setw(7) << std::setprecision(1) <<
trk[itr].TrjPos[itj](0);
3162 myprt <<
std::right << std::setw(7) << std::setprecision(1) <<
trk[itr].TrjPos[itj](1);
3163 myprt <<
std::right << std::setw(7) << std::setprecision(1) <<
trk[itr].TrjPos[itj](2);
3164 myprt << std::setw(4) <<
trk[itr].VtxIndex[0] << std::setw(4) <<
trk[itr].VtxIndex[1];
3165 myprt << std::setw(4) <<
trk[itr].GoodEnd[0];
3166 myprt << std::setw(4) <<
trk[itr].GoodEnd[1];
3167 myprt << std::setw(4) <<
trk[itr].ChgOrder;
3168 myprt <<
std::right << std::setw(10) << std::setprecision(3) <<
trk[itr].TrjDir[itj](2);
3170 myprt <<
std::right << std::setw(5) <<
trk[itr].PDGCode <<
" ";
3171 for (
unsigned short ii = 0; ii <
trk[itr].ClsEvtIndices.size(); ++ii)
3172 myprt <<
" " <<
trk[itr].ClsEvtIndices[ii];
3183 unsigned short iTime;
3185 myprt <<
"******* PrintClusters ********* Num_Clusters_in Wire:Time\n";
3186 myprt <<
"vtx Index X Y Z Pln0 Pln1 Pln2 Pln0 Pln1 Pln2\n";
3187 for (
unsigned short ivx = 0; ivx <
vtx.size(); ++ivx) {
3188 myprt <<
std::right << std::setw(3) << ivx << std::setw(7) << ivx;
3189 myprt << std::fixed;
3190 myprt <<
std::right << std::setw(7) << std::setprecision(1) <<
vtx[ivx].X;
3191 myprt <<
std::right << std::setw(7) << std::setprecision(1) <<
vtx[ivx].Y;
3192 myprt <<
std::right << std::setw(7) << std::setprecision(1) <<
vtx[ivx].Z;
3193 myprt <<
std::right << std::setw(5) <<
vtx[ivx].nClusInPln[0];
3194 myprt <<
std::right << std::setw(5) <<
vtx[ivx].nClusInPln[1];
3195 myprt <<
std::right << std::setw(5) <<
vtx[ivx].nClusInPln[2];
3200 myprt <<
std::right << std::setw(7) << wire <<
":" << time;
3207 for (
unsigned short ipl = 0; ipl < nplanes; ++ipl) {
3208 myprt <<
">>>>>>>>>> Cluster chains in Plane " << ipl <<
"\n";
3209 myprt <<
"ipl ccl Len Chg W0:T0 Ang0 Dir0 Vx0 mBk0 W1:T1 Ang1 Dir1 Vx1 " 3210 " mBk1 InTk cls:Order \n";
3211 for (
unsigned short ccl = 0; ccl <
clsChain[ipl].size(); ++ccl) {
3216 for (
unsigned short end = 0;
end < 2; ++
end) {
3219 << std::setprecision(1) << iTime;
3220 if (iTime < 10) { myprt <<
" "; }
3221 else if (iTime < 100) {
3224 else if (iTime < 1000)
3226 myprt <<
std::right << std::setw(7) << std::setprecision(2)
3230 myprt << std::fixed <<
std::right << std::setw(6) << std::setprecision(1)
3235 for (
unsigned short ii = 0; ii <
clsChain[ipl][ccl].ClsIndex.size(); ++ii)
3236 myprt <<
" " <<
clsChain[ipl][ccl].ClsIndex[ii] <<
":" <<
clsChain[ipl][ccl].Order[ii];
3240 myprt <<
">>>>>>>>>> Clusters in Plane " << ipl <<
"\n";
3241 myprt <<
"ipl icl Evt Len Chg W0:T0 Ang0 Dir0 Vx0 CN0 W1:T1 Ang1 " 3242 "Dir1 Vx1 CN1 InTk Brk0 MrgEr0 Brk1 MrgEr1\n";
3243 for (
unsigned short icl = 0; icl <
cls[ipl].size(); ++icl) {
3246 myprt <<
std::right << std::setw(5) <<
cls[ipl][icl].EvtIndex;
3247 myprt <<
std::right << std::setw(6) <<
cls[ipl][icl].Length;
3248 myprt <<
std::right << std::setw(8) << (int)
cls[ipl][icl].TotChg;
3249 for (
unsigned short end = 0;
end < 2; ++
end) {
3250 iTime =
cls[ipl][icl].Time[
end];
3251 myprt <<
std::right << std::setw(5) << (int)
cls[ipl][icl].Wire[
end] <<
":" << iTime;
3252 if (iTime < 10) { myprt <<
" "; }
3253 else if (iTime < 100) {
3256 else if (iTime < 1000)
3258 myprt <<
std::right << std::setw(7) << std::setprecision(2) <<
cls[ipl][icl].Angle[
end];
3261 myprt << std::fixed <<
std::right << std::setw(5) << std::setprecision(1)
3262 <<
cls[ipl][icl].ChgNear[
end];
3264 myprt << std::fixed;
3265 myprt <<
std::right << std::setw(5) <<
cls[ipl][icl].InTrack;
3266 myprt <<
std::right << std::setw(5) << (int)
cls[ipl][icl].BrkIndex[0];
3267 myprt <<
std::right << std::setw(7) << std::setprecision(1)
3268 <<
cls[ipl][icl].MergeError[0];
3269 myprt <<
std::right << std::setw(5) << (int)
cls[ipl][icl].BrkIndex[1];
3270 myprt <<
std::right << std::setw(7) << std::setprecision(1)
3271 <<
cls[ipl][icl].MergeError[1];
3281 float slp = fabs(slope);
3282 if (slp > 10.) slp = 30.;
3284 return 1 + 0.05 * slp * slp;
3289 unsigned short wire1,
3291 unsigned short wire2,
3298 unsigned short w1 = wire1;
3299 unsigned short w2 = wire2;
3303 if (w1 == w2) { slp = 0; }
3311 slp = (t2 -
t1) / (w2 - w1);
3314 unsigned short wire;
3320 if (wire < w1)
continue;
3321 if (wire > w2)
continue;
3322 prtime = t1 + (wire - w1) * slp;
3323 if (prtime >
allhits[
hit]->PeakTimePlusRMS(3))
continue;
3324 if (prtime <
allhits[
hit]->PeakTimeMinusRMS(3))
continue;
3338 for (ipl = 0; ipl < 3; ++ipl) {
3348 unsigned short oldipl = 0;
3350 if (static_cast<geo::TPCID const&>(
allhits[
hit]->
WireID()) != tpcid)
continue;
3360 mf::LogError(
"CCTM") <<
"Hits are not in increasing-plane order\n";
3374 for (ipl = 0; ipl < nplanes; ++ipl) {
3376 mf::LogError(
"CCTM") <<
"Invalid first/last wire in plane " << ipl;
3383 for (ipl = 0; ipl < nplanes; ++ipl)
3391 int sflag, nwires, wire;
3392 unsigned int indx, thisWire, thisHit, lastFirstHit;
3394 for (ipl = 0; ipl < nplanes; ++ipl) {
3399 for (wire = 0; wire < nwires; ++wire)
3400 WireHitRange[ipl][wire] = std::make_pair(sflag, sflag);
3410 indx =
lastWire[ipl] - firstWire[ipl];
3411 int tmp1 = lastFirstHit;
3415 lastFirstHit = thisHit;
3417 else if (thisWire <
lastWire[ipl]) {
3418 mf::LogError(
"CCTM") <<
"Hit not in proper order in plane " << ipl;
3424 indx =
lastWire[ipl] - firstWire[ipl];
3425 int tmp1 = lastFirstHit;
3434 for (ipl = 0; ipl < nplanes; ++ipl) {
3437 if (
firstHit[ipl] < INT_MAX)
continue;
3438 if (
lastHit[ipl] > 0)
continue;
3439 std::cout <<
"FWHR problem\n";
3443 unsigned int nht = 0;
3444 std::vector<bool> hchk(
allhits.size());
3445 for (
unsigned int ii = 0; ii < hchk.size(); ++ii)
3447 for (ipl = 0; ipl < nplanes; ++ipl) {
3450 if (indx > lastWire[ipl]) {
3451 std::cout <<
"FWHR bad index " << indx <<
"\n";
3458 for (
unsigned int hit = firhit;
hit < lashit; ++
hit) {
3461 std::cout <<
"FWHR: Bad hchk " <<
hit <<
" size " <<
allhits.size() <<
"\n";
3466 std::cout <<
"FWHR bad plane " <<
allhits[
hit]->WireID().Plane <<
" " << ipl
3467 <<
" or wire " <<
allhits[
hit]->WireID().Wire <<
" " <<
w <<
"\n";
3475 std::cout <<
"FWHR hit count problem " << nht <<
" " <<
allhits.size() <<
"\n";
3476 for (
unsigned int ii = 0; ii < hchk.size(); ++ii)
3478 std::cout <<
" " << ii <<
" " <<
allhits[ii]->WireID().Plane <<
" " 3479 <<
allhits[ii]->WireID().Wire <<
" " << (int)
allhits[ii]->PeakTime() <<
"\n";
std::string fHitModuleLabel
details::range_type< T > Iterate() const
Initializes the specified ID with the ID of the first cryostat.
MaybeLogger_< ELseverityLevel::ELsev_info, true > LogVerbatim
float Length(const PFPStruct &pfp)
std::array< short, 2 > VtxIndex
Length_t WireCoordinate(Point_t const &pos, PlaneID const &planeid) const
Returns the index of the nearest wire to the specified position.
std::vector< MatchPars > matcomb
float ChargeNear(geo::PlaneID const &planeid, unsigned short wire1, float time1, unsigned short wire2, float time2)
constexpr auto const & right(const_AssnsIter< L, R, D, Dir > const &a, const_AssnsIter< L, R, D, Dir > const &b)
art::ServiceHandle< geo::Geometry const > geom
std::string fVertexModuleLabel
void StoreTrack(detinfo::DetectorPropertiesData const &detProp, art::FindManyP< recob::Hit > const &fmCluHits, unsigned short imat, unsigned short procCode, geo::TPCID const &tpcid)
std::array< float, 3 > ChgNorm
std::vector< Vector_t > convertCollToVector(std::vector< Vector > const &coll)
std::vector< float > fChgAsymFactor
std::array< unsigned short, 3 > End
Declaration of signal hit object.
std::array< unsigned int, 3 > firstWire
void MakeClusterChains(detinfo::DetectorPropertiesData const &detProp, art::FindManyP< recob::Hit > const &fmCluHits, geo::TPCID const &tpcid)
Length_t DetHalfWidth(TPCID const &tpcid=tpc_zero) const
Returns the half width of the active volume of the specified TPC.
Planes which measure X direction.
The data type to uniquely identify a Plane.
Geometry information for a single TPC.
std::vector< art::Ptr< recob::Hit > > allhits
void FindMaybeVertices(geo::TPCID const &tpcid)
std::array< short, 2 > Dir
TrackTrajectory::Flags_t Flags_t
constexpr auto abs(T v)
Returns the absolute value of the argument.
float StartWire() const
Returns the wire coordinate of the start of the cluster.
float dXClTraj(art::FindManyP< recob::Hit > const &fmCluHits, unsigned short ipl, unsigned short icl1, unsigned short end1)
std::array< std::vector< art::Ptr< recob::Hit > >, 3 > seedHits
std::string fClusterModuleLabel
void FitVertices(detinfo::DetectorPropertiesData const &detProp, geo::TPCID const &tpcid)
TrackTrajectoryAlg fTrackTrajectoryAlg
ROOT::Math::SMatrix< Double32_t, 5, 5, ROOT::Math::MatRepSym< Double32_t, 5 >> SMatrixSym55
std::array< float, 2 > ChgNear
void FillWireHitRange(geo::TPCID inTPCID)
Set of hits with a 2D structure.
std::array< short, 2 > BrkIndex
float EndTick() const
Returns the tick coordinate of the end of the cluster.
std::vector< unsigned short > pfpToTrkID
void FillTrkHits(art::FindManyP< recob::Hit > const &fmCluHits, unsigned short imat, unsigned short nplanes)
data_const_reference data(size_type i) const
std::vector< TVector3 > TrjPos
std::vector< TrkPar > trk
bool FindMissingCluster(unsigned short kpl, short &kcl, unsigned short &kend, float kWir, float kX, float okWir, float okX)
MaybeLogger_< ELseverityLevel::ELsev_error, false > LogError
Cluster finding and building.
std::vector< float > fAngleMatchErr
std::array< float, 2 > ChgNear
std::array< std::vector< clPar >, 3 > cls
float StartAngle() const
Returns the starting angle of the cluster.
Definition of vertex object for LArSoft.
std::array< float, 2 > Slope
Length_t DetLength(TPCID const &tpcid=tpc_zero) const
Returns the length of the active volume of the specified TPC.
std::array< bool, 2 > EndInTPC
double Length() const
Length is associated with z coordinate [cm].
TPCGeo const & TPC(TPCID const &tpcid=tpc_zero) const
Returns the specified TPC.
PutHandle< PROD > put(std::unique_ptr< PROD > &&edp, std::string const &instance={})
decltype(auto) constexpr end(T &&obj)
ADL-aware version of std::end.
decltype(auto) constexpr size(T &&obj)
ADL-aware version of std::size.
std::array< std::vector< art::Ptr< recob::Hit > >, 3 > trkHits
std::array< std::vector< ClsChainPar >, 3 > clsChain
float EndCharge() const
Returns the charge on the last wire of the cluster.
Point_t GetCenter() const
Returns the center of the TPC volume in world coordinates [cm].
void TagCosmics(geo::TPCID const &tpcid)
#define DEFINE_ART_MODULE(klass)
IDparameter< geo::WireID > WireID
Member type of validated geo::WireID parameter.
A trajectory in space reconstructed from hits.
void FillEndMatch(detinfo::DetectorPropertiesData const &detProp, geo::TPCID const &tpcid, MatchPars &match)
std::vector< float > fMatchMinLen
std::vector< unsigned short > Order
float AngleFactor(float slope)
std::array< short, 2 > VtxIndex
std::array< short, 2 > mVtxIndex
std::vector< short > DtrID
double ConvertXToTicks(double X, int p, int t, int c) const
std::vector< vtxPar > vtx
std::array< float, 2 > MergeError
Provides recob::Track data product.
void produce(art::Event &evt) override
void fill(const art::PtrVector< recob::Hit > &hits, int only_plane)
std::array< bool, 2 > GoodEnd
std::vector< Point_t > convertCollToPoint(std::vector< Point > const &coll)
std::array< float, 2 > Wire
The data type to uniquely identify a TPC.
PlaneID_t Plane
Index of the plane within its TPC.
Declaration of cluster object.
void TrackTrajectory(std::array< std::vector< geo::WireID >, 3 > trkWID, std::array< std::vector< double >, 3 > trkX, std::array< std::vector< double >, 3 > trkXErr, std::vector< TVector3 > &TrajPos, std::vector< TVector3 > &TrajDir)
bool DupMatch(MatchPars &match, unsigned short nplanes)
std::array< unsigned int, 3 > firstHit
std::array< float, 2 > Charge
void VtxMatch(detinfo::DetectorPropertiesData const &detProp, art::FindManyP< recob::Hit > const &fmCluHits, geo::TPCID const &tpcid)
std::vector< unsigned short > ClsIndex
std::vector< unsigned short > ClsEvtIndices
Detector simulation of raw signals on wires.
std::array< unsigned int, 3 > lastWire
Class defining a plane for tracking. It provides various functionalities to convert track parameters ...
void MakeFamily(geo::TPCID const &tpcid)
double ConvertTicksToX(double ticks, int p, int t, int c) const
double HalfHeight() const
Height is associated with y coordinate [cm].
bool greaterThan(CluLen c1, CluLen c2)
bool CreateAssn(art::Event &evt, std::vector< T > const &a, art::Ptr< U > const &b, art::Assns< U, T > &assn, std::string a_instance, size_t index=UINT_MAX)
Creates a single one-to-one association.
ROOT::Math::PositionVector3D< ROOT::Math::Cartesian3D< double >, ROOT::Math::GlobalCoordinateSystemTag > Point_t
Type for representation of position in physical 3D space.
void VertexFit(std::vector< std::vector< geo::WireID >> const &hitWID, std::vector< std::vector< double >> const &hitX, std::vector< std::vector< double >> const &hitXErr, TVector3 &VtxPos, TVector3 &VtxPosErr, std::vector< TVector3 > &TrkDir, std::vector< TVector3 > &TrkDirErr, float &ChiDOF) const
Hierarchical representation of particle flow.
bool getByLabel(std::string const &label, std::string const &instance, Handle< PROD > &result) const
void FillEndMatch2(MatchPars &match)
Utility object to perform functions of association.
std::array< float, 2 > Time
float StartCharge() const
Returns the charge on the first wire of the cluster.
std::array< float, 2 > Slope
std::array< short, 2 > VtxIndex
void PrintClusters(detinfo::DetectorPropertiesData const &detProp, geo::TPCID const &tpcid) const
std::array< unsigned int, 3 > lastHit
MaybeLogger_< ELseverityLevel::ELsev_warning, false > LogWarning
double ThirdPlaneSlope(PlaneID const &pid1, double slope1, PlaneID const &pid2, double slope2, PlaneID const &output_plane) const
Returns the slope on the third plane, given it in the other two.
void PlnMatch(detinfo::DetectorPropertiesData const &detProp, geo::TPCID const &tpcid)
unsigned int Nwires(PlaneID const &planeid) const
Returns the total number of wires in the specified plane.
bool lessThan(CluLen c1, CluLen c2)
std::array< short, 2 > Dir
unsigned int Nplanes(TPCID const &tpcid=tpc_zero) const
Returns the total number of planes in the specified TPC.
decltype(auto) constexpr begin(T &&obj)
ADL-aware version of std::begin.
std::vector< TVector3 > TrjDir
void FillWireHitRange(geo::TPCID const &tpcid)
bool IntersectionPoint(WireID const &wid1, WireID const &wid2, double &y, double &z) const
Returns the intersection point of two wires.
std::array< float, 2 > Angle
std::vector< bool > fMakeAlgTracks
std::vector< short > fMatchAlgs
std::vector< float > fXMatchErr
void fill_ptr_vector(std::vector< Ptr< T >> &ptrs, H const &h)
float EndAngle() const
Returns the ending angle of the cluster.
recob::tracking::Plane Plane
float ChargeAsym(std::array< float, 3 > &mChg)
std::array< float, 2 > Wire
std::array< short, 3 > Cls
float StartTick() const
Returns the tick coordinate of the start of the cluster.
std::array< std::vector< art::Ptr< recob::Hit > >, 3 > TrkHits
Length_t WirePitch(PlaneID const &planeid=plane_zero) const
Returns the distance between two consecutive wires.
void FillChgNear(detinfo::DetectorPropertiesData const &detProp, geo::TPCID const &tpcid)
VertexFitAlg fVertexFitAlg
std::array< short, 2 > Time
std::array< float, 3 > Chg
art framework interface to geometry description
double HalfWidth() const
Width is associated with x coordinate [cm].
Track from a non-cascading particle.A recob::Track consists of a recob::TrackTrajectory, plus additional members relevant for a "fitted" track:
std::array< float, 2 > Angle
cet::coded_exception< error, detail::translate > exception
Encapsulate the construction of a single detector plane.
std::array< std::vector< std::pair< int, int > >, 3 > WireHitRange
std::array< std::vector< unsigned short >, 3 > vxCls
std::array< unsigned short, 3 > nClusInPln
unsigned short fNVtxTrkHitsFit
float Integral() const
Returns the total charge of the cluster from hit shape.
float EndWire() const
Returns the wire coordinate of the end of the cluster.
void SortMatches(detinfo::DetectorPropertiesData const &detProp, art::FindManyP< recob::Hit > const &fmCluHits, unsigned short procCode, geo::TPCID const &tpcid)
std::array< short, 2 > mBrkIndex