29 if (trk == t.first.Track())
return true;
36 if (!
Has(t.first.Track()))
return false;
43 if (!rootTrk)
throw cet::exception(
"pma::VtxCandidate") <<
"Broken track.";
47 if (!rootAssn)
throw cet::exception(
"pma::VtxCandidate") <<
"Broken track.";
63 for (
size_t t = 0; t <
fAssigned.size(); t++) {
65 if (!trk_t)
throw cet::exception(
"pma::VtxCandidate") <<
"Broken track.";
67 for (
size_t u = 0; u <
fAssigned.size(); u++)
70 if (!trk_u)
throw cet::exception(
"pma::VtxCandidate") <<
"Broken track.";
82 if (c.first.Track()->Length() > minLength) n++;
97 for (
size_t n = 0;
n < trk.
Track()->
Nodes().size() - 1;
n++) {
114 if (d_best < kMaxDistToTrack) {
130 size_t n_best = 0, m_best = 0;
133 double lm, ln, l_best = 0;
134 for (
size_t m = 0; m < p0->
Nodes().size() - 1; m++) {
141 for (
size_t n = 0;
n < trk.
Track()->
Nodes().size() - 1;
n++) {
153 double d_dist = (d_best -
d) / d_best;
154 if (lm + ln > 0.8 * d_dist * l_best)
166 if (d_best < kMaxDistToTrack) {
184 for (
size_t n = 0;
n < trk.
Track()->
Nodes().size() - 1;
n++) {
206 int tpc = trk->
Nodes()[t.second]->TPC();
207 int cryo = trk->
Nodes()[t.second]->Cryo();
212 if (tpcgeom.HasPlane(
geo::kU)) {
217 if (tpcgeom.HasPlane(
geo::kV)) {
222 if (tpcgeom.HasPlane(
geo::kZ)) {
227 mse += m / (double)k;
230 return mse / fAssigned.size();
238 double dw =
fErr[0] * other.
fErr[0] * dx * dx +
fErr[1] * other.
fErr[1] * dy * dy +
246 size_t max_l_idx = 0;
247 double l, max_l = 0.0;
248 for (
size_t i = 0; i <
fAssigned.size() - 1; i++) {
249 l =
fAssigned[i].first.Track()->Length();
257 dir_i *= 1.0 / dir_i.Mag();
262 for (
size_t j = 0; j <
fAssigned.size(); j++)
263 if ((j != max_l_idx) && (
fAssigned[j].first.Track()->Length() > minLength)) {
268 dir_j *= 1.0 / dir_j.Mag();
269 a = fabs(dir_i * dir_j);
270 if (a < min) min = a;
273 return 180.0 * acos(min) / TMath::Pi();
284 double dw =
Test(other);
292 if (!
Has(t.first.Track())) {
298 double mse0 =
fMse, mse1 = other.
fMse;
300 <<
"try: " << d <<
" mse0:" << sqrt(mse0) <<
" mse1:" << sqrt(mse1);
305 <<
"out: " <<
Size() <<
" mse:" << sqrt(mse) <<
" dw:" << dw;
332 std::vector<pma::Segment3D*> segments;
333 std::vector<std::pair<TVector3, TVector3>> lines;
334 std::vector<double> weights;
341 double segLength = seg->
Length();
345 std::pair<TVector3, TVector3> endpoints(vtx1->
Point3D(), vtx2->
Point3D());
346 double dy = endpoints.first.Y() - endpoints.second.Y();
347 double fy_norm = asin(fabs(dy) / segLength) / (0.5 * TMath::Pi());
348 double w = 1.0 - pow(fy_norm - 1.0, 12);
349 if (w < 0.3) w = 0.3;
351 lines.push_back(endpoints);
352 segments.push_back(seg);
353 weights.push_back(w);
358 fErr.SetXYZ(0., 0., 0.);
362 if (resultMse < 0.0) {
363 mf::LogWarning(
"pma::VtxCandidate") <<
"Cannot compute crossing point.";
370 for (
size_t s = 0; s < segments.size(); s++) {
380 fErr[0] += weights[s] * weights[s];
384 fCenter[0] += weights[s] * pproj.X();
393 fErr *= 1.0 / segments.size();
406 mf::LogError(
"pma::VtxCandidate") <<
"Tracks already attached to the vertex.";
412 <<
"JoinTracks (" <<
fAssigned.size() <<
") at:" 417 while (t < src.
size()) {
418 if (c.first.Track() == src[t].Track()) {
429 for (
auto& c : fAssigned)
430 for (
auto const& t : tracks.
tracks())
431 if (c.first.Track() == t.Track()) {
432 c.first.SetTreeId(t.TreeId());
437 std::vector<int> treeIds;
441 bool hasInnerCenter =
false;
443 for (
size_t i = 0; i < fAssigned.size(); i++) {
447 int key = fAssigned[i].first.Key();
448 int tid = fAssigned[i].first.TreeId();
449 size_t idx = fAssigned[i].second;
452 <<
" (nodes:" << trk->
Nodes().size() <<
")";
454 if (!
has(treeIds, tid))
456 treeIds.push_back(tid);
461 mf::LogError(
"pma::VtxCandidate") <<
"Root of the tree not found in tracks collection.";
464 TVector3 p0(trk->
Nodes()[idx]->Point3D());
465 TVector3 p1(trk->
Nodes()[idx + 1]->Point3D());
467 int tpc0 = trk->
Nodes()[idx]->TPC();
468 int tpc1 = trk->
Nodes()[idx + 1]->TPC();
470 int cryo0 = trk->
Nodes()[idx]->Cryo();
471 int cryo1 = trk->
Nodes()[idx + 1]->Cryo();
482 vtxCenter = trk->
Nodes().front();
488 if (trk->
AttachTo(vtxCenter)) nOK++;
494 mf::LogVerbatim(
"pma::VtxCandidate") <<
" flip trk to make new center";
496 vtxCenter = trk->
Nodes().front();
500 vtxCenter = trk->
Nodes().back();
507 mf::LogVerbatim(
"pma::VtxCandidate") <<
" flip trk to attach to inner";
509 if (trk->
AttachTo(vtxCenter)) nOK++;
518 bool canFlipPrev =
true;
519 if (vtxCenter && vtxCenter->
Prev()) {
527 if (hasInnerCenter || !canFlipPrev) {
548 mf::LogVerbatim(
"pma::VtxCandidate") <<
" add center at end of segment";
552 mf::LogVerbatim(
"pma::VtxCandidate") <<
" center at start of segment - no action";
560 <<
" trk size:" << trk->
size() <<
" (nodes:" << trk->
Nodes().size() <<
")";
565 <<
" t0 size:" << t0->
size() <<
" (nodes:" << t0->
Nodes().size() <<
")";
569 tracks.
tracks().emplace_back(t0, key, tid);
572 vtxCenter = trk->
Nodes().front();
577 if (trk->
AttachTo(vtxCenter)) nOK += 2;
584 hasInnerCenter =
true;
604 mf::LogVerbatim(
"pma::VtxCandidate") <<
" add center at end of segment";
608 mf::LogVerbatim(
"pma::VtxCandidate") <<
" center at start of segment - no action";
619 rootBranch = seg->
Parent();
626 for (
size_t j = 0; j < branches.size(); ++j) {
627 if (branches[j]->AttachTo(innerCenter,
true)) {}
629 vtxCenter = innerCenter;
632 vtxCenter = innerCenter;
646 rootSeg = static_cast<pma::Segment3D*>(vtxCenter->
Next(0));
647 else if (vtxCenter->
Prev())
648 rootSeg = static_cast<pma::Segment3D*>(vtxCenter->
Prev());
650 throw cet::exception(
"pma::VtxCandidate") <<
"Vertex with no segments attached.";
653 if (!rootTrk) rootTrk = rootSeg->
Parent();
655 std::vector<pma::Track3D const*> branchesToRemove;
656 bool noLoops = rootTrk->
GetBranches(branchesToRemove);
659 if (noLoops && (nOK > 1)) {
676 if (noLoops && tuneOK) {
678 for (
auto& c : backupTracks.
tracks())
682 mf::LogVerbatim(
"pma::VtxCandidate") <<
"restore tracks from backup....";
683 for (
int tid : treeIds) {
685 while (t < tracks.
size()) {
686 if (tracks[t].TreeId() == tid) {
687 tracks[t].DeleteTrack();
694 for (
const auto& c : backupTracks.
tracks())
700 mf::LogError(
"pma::VtxCandidate") <<
"Cannot create common vertex";
Vertex finding helper for the Projection Matching Algorithm.
code to link reconstructed objects back to the MC truth information
MaybeLogger_< ELseverityLevel::ELsev_info, true > LogVerbatim
TVector3 const & Point3D() const
bool Has(pma::Track3D *trk) const
double SolveLeastSquares3D(const std::vector< std::pair< TVector3, TVector3 >> &lines, TVector3 &result)
pma::Track3D * getTreeCopy(pma::TrkCandidateColl &dst, size_t trkIdx, bool isRoot=true)
bool AttachTo(pma::Node3D *vStart, bool noFlip=false)
Implementation of the Projection Matching Algorithm.
bool IsAttachedTo(pma::Track3D const *trk) const
double Dist2(const TVector2 &v1, const TVector2 &v2)
double Test(const VtxCandidate &other) const
bool JoinTracks(detinfo::DetectorPropertiesData const &detProp, pma::TrkCandidateColl &tracks, pma::TrkCandidateColl &src)
void erase_at(size_t pos)
Implementation of the Projection Matching Algorithm.
bool GetBranches(std::vector< pma::Track3D const * > &branches, bool skipFirst=false) const
std::vector< std::pair< pma::TrkCandidate, size_t > > fAssigned
virtual unsigned int NextCount(void) const
std::vector< pma::Track3D * > GetBranches() const
Planes which measure Z direction.
Implementation of the Projection Matching Algorithm.
virtual pma::SortedObjectBase * Next(unsigned int=0) const
MaybeLogger_< ELseverityLevel::ELsev_error, false > LogError
double TuneFullTree(double eps=0.001, double gmax=50.0)
TPCGeo const & TPC(TPCID const &tpcid=tpc_zero) const
Returns the specified TPC.
pma::Track3D * Split(detinfo::DetectorPropertiesData const &detProp, size_t idx, bool try_start_at_idx=true)
bool CanFlip() const
Check if the track can be flipped without breaking any other track.
TVector2 GetProjectionToSegment(const TVector2 &p, const TVector2 &p0, const TVector2 &p1)
static constexpr double kMinDistToNode
int getCandidateIndex(pma::Track3D const *candidate) const
bool SetPoint3D(const TVector3 &p3d)
static constexpr double kMaxDistToTrack
The data type to uniquely identify a TPC.
Definition of data types for geometry description.
double MaxAngle(double minLength=0.0) const
void InsertNode(detinfo::DetectorPropertiesData const &detProp, TVector3 const &p3d, size_t at_idx, unsigned int tpc, unsigned int cryo)
bool MergeWith(const VtxCandidate &other)
Implementation of the Projection Matching Algorithm.
virtual pma::SortedObjectBase * Next(unsigned int index=0) const
bool Add(const pma::TrkCandidate &trk)
double GetSegmentProjVector(const TVector2 &p, const TVector2 &p0, const TVector2 &p1)
TVector2 GetProjectionToPlane(const TVector3 &p, unsigned int plane, unsigned int tpc, unsigned int cryo)
MaybeLogger_< ELseverityLevel::ELsev_warning, false > LogWarning
bool IsAttached(pma::Track3D *trk) const
double GetDistance2To(const TVector3 &p3d) const override
Distance [cm] from the 3D segment to the point 3D.
std::vector< pma::Node3D * > const & Nodes() const noexcept
virtual pma::SortedObjectBase * Prev(void) const
bool Flip(const detinfo::DetectorPropertiesData &detProp, std::vector< pma::Track3D * > &allTracks)
pma::Track3D * Track() const
pma::Track3D * Parent(void) const
bool AttachBackTo(pma::Node3D *vStart)
pma::Segment3D * NextSegment(pma::Node3D *vtx) const
double Length(void) const
void push_back(const TrkCandidate &trk)
art framework interface to geometry description
std::vector< TrkCandidate > const & tracks() const
cet::coded_exception< error, detail::translate > exception
bool has(const std::vector< int > &v, int id) const