LArSoft  v09_90_00
Liquid Argon Software toolkit - https://larsoft.org/
pma::VtxCandidate Class Reference

#include "PmaVtxCandidate.h"

Public Member Functions

 VtxCandidate (double segMinLength=0.5)
 
bool Has (pma::Track3D *trk) const
 
bool Has (const VtxCandidate &other) const
 
bool IsAttached (pma::Track3D *trk) const
 
bool IsAttached (const VtxCandidate &other) const
 
bool HasLoops () const
 
bool Add (const pma::TrkCandidate &trk)
 
double ComputeMse2D ()
 
double Test (const VtxCandidate &other) const
 
double MaxAngle (double minLength=0.0) const
 
size_t Size () const
 
size_t Size (double minLength) const
 
bool MergeWith (const VtxCandidate &other)
 
double Compute ()
 
bool JoinTracks (detinfo::DetectorPropertiesData const &detProp, pma::TrkCandidateColl &tracks, pma::TrkCandidateColl &src)
 
const TVector3 & Center () const
 
double Mse () const
 
double Mse2D () const
 
std::pair< pma::Track3D *, size_t > Track (size_t i) const
 

Static Public Attributes

static constexpr double kMaxDistToTrack {4.0}
 
static constexpr double kMinDistToNode {2.0}
 

Private Member Functions

bool has (const std::vector< int > &v, int id) const
 

Private Attributes

bool tracksJoined
 
double fSegMinLength
 
double fMse
 
double fMse2D
 
std::vector< std::pair< pma::TrkCandidate, size_t > > fAssigned
 
TVector3 fCenter
 
TVector3 fErr
 

Detailed Description

Definition at line 24 of file PmaVtxCandidate.h.

Constructor & Destructor Documentation

pma::VtxCandidate::VtxCandidate ( double  segMinLength = 0.5)
inline

Definition at line 29 of file PmaVtxCandidate.h.

References Add(), ComputeMse2D(), Has(), HasLoops(), IsAttached(), MaxAngle(), fhicl::other, and Test().

30  : tracksJoined(false)
31  , fSegMinLength(segMinLength)
32  , fMse(0.0)
33  , fMse2D(0.0)
34  , fCenter(0., 0., 0.)
35  , fErr(0., 0., 0.)
36  {}

Member Function Documentation

bool pma::VtxCandidate::Add ( const pma::TrkCandidate trk)

Definition at line 86 of file PmaVtxCandidate.cxx.

References Compute(), ComputeMse2D(), d, fAssigned, fCenter, fMse, fMse2D, fSegMinLength, pma::Segment3D::GetDistance2To(), IsAttached(), kMaxDistToTrack, pma::Element3D::Length(), n, pma::Track3D::NextSegment(), pma::Track3D::Nodes(), and pma::TrkCandidate::Track().

Referenced by pma::PMAlgVertexing::firstPassCandidates(), pma::PMAlgVertexing::secondPassCandidates(), and VtxCandidate().

87 {
88  if (IsAttached(trk.Track())) return false;
89 
90  fAssigned.emplace_back(trk, 0);
91 
92  double d, d_best;
93  double mse, min_mse = kMaxDistToTrack * kMaxDistToTrack;
94  if (fAssigned.size() > 2) {
95  size_t n_best = 0;
96  d_best = kMaxDistToTrack;
97  for (size_t n = 0; n < trk.Track()->Nodes().size() - 1; n++) {
98  pma::Segment3D* seg = trk.Track()->NextSegment(trk.Track()->Nodes()[n]);
99  if (seg->Length() < fSegMinLength) continue;
100 
101  fAssigned.back().second = n;
102 
103  mse = Compute();
104  if (mse < min_mse) {
105  d = sqrt(seg->GetDistance2To(fCenter));
106  if (d < d_best) {
107  min_mse = mse;
108  n_best = n;
109  d_best = d;
110  }
111  }
112  }
113 
114  if (d_best < kMaxDistToTrack) {
115  fAssigned.back().second = n_best;
116  fMse = Compute();
117  fMse2D = ComputeMse2D();
118  return true;
119  }
120  else {
121  fAssigned.pop_back();
122  fMse = Compute();
123  fMse2D = ComputeMse2D();
124  return false;
125  }
126  }
127  else if (fAssigned.size() == 2) {
128  pma::Track3D* p0 = fAssigned.front().first.Track();
129 
130  size_t n_best = 0, m_best = 0;
131  d_best = kMaxDistToTrack;
132 
133  double lm, ln, l_best = 0;
134  for (size_t m = 0; m < p0->Nodes().size() - 1; m++) {
135  pma::Segment3D* seg_m = p0->NextSegment(p0->Nodes()[m]);
136  lm = seg_m->Length();
137  if (lm < fSegMinLength) continue;
138 
139  fAssigned.front().second = m;
140 
141  for (size_t n = 0; n < trk.Track()->Nodes().size() - 1; n++) {
142  pma::Segment3D* seg_n = trk.Track()->NextSegment(trk.Track()->Nodes()[n]);
143  ln = seg_n->Length();
144  if (ln < fSegMinLength) continue;
145 
146  fAssigned.back().second = n;
147 
148  mse = Compute(); // std::cout << mse << std::endl;
149 
150  d = sqrt(ComputeMse2D());
151 
152  if (d < d_best) {
153  double d_dist = (d_best - d) / d_best;
154  if (lm + ln > 0.8 * d_dist * l_best) // take closer if not much shorter
155  {
156  min_mse = mse;
157  n_best = n;
158  m_best = m;
159  d_best = d;
160  l_best = lm + ln;
161  }
162  }
163  }
164  }
165 
166  if (d_best < kMaxDistToTrack) {
167  fAssigned.front().second = m_best;
168  fAssigned.back().second = n_best;
169  fMse = Compute();
170 
171  fMse2D = ComputeMse2D();
172 
173  return true;
174  }
175  else {
176  fAssigned.pop_back();
177  fCenter.SetXYZ(0., 0., 0.);
178  fMse = 0;
179  fMse2D = 0;
180  return false;
181  }
182  }
183  else {
184  for (size_t n = 0; n < trk.Track()->Nodes().size() - 1; n++) {
185  pma::Segment3D* seg = trk.Track()->NextSegment(trk.Track()->Nodes()[n]);
186  if (seg->Length() >= fSegMinLength) { return true; }
187  }
188  fAssigned.pop_back();
189  fCenter.SetXYZ(0., 0., 0.);
190  fMse = 0;
191  fMse2D = 0;
192  return false;
193  }
194 }
std::vector< std::pair< pma::TrkCandidate, size_t > > fAssigned
Float_t d
Definition: plot.C:235
static constexpr double kMaxDistToTrack
bool IsAttached(pma::Track3D *trk) const
double GetDistance2To(const TVector3 &p3d) const override
Distance [cm] from the 3D segment to the point 3D.
Char_t n[5]
std::vector< pma::Node3D * > const & Nodes() const noexcept
Definition: PmaTrack3D.h:274
pma::Track3D * Track() const
pma::Segment3D * NextSegment(pma::Node3D *vtx) const
double Length(void) const
Definition: PmaElement3D.h:53
const TVector3& pma::VtxCandidate::Center ( ) const
inline

Definition at line 67 of file PmaVtxCandidate.h.

References fCenter.

67 { return fCenter; }
double pma::VtxCandidate::Compute ( )

Definition at line 330 of file PmaVtxCandidate.cxx.

References fAssigned, fCenter, fErr, fSegMinLength, pma::GetProjectionToSegment(), pma::Element3D::Length(), pma::SortedObjectBase::Next(), pma::Track3D::NextSegment(), pma::Track3D::Nodes(), pma::Node3D::Point3D(), pma::SolveLeastSquares3D(), and w.

Referenced by Add(), MergeWith(), and Size().

331 {
332  std::vector<pma::Segment3D*> segments;
333  std::vector<std::pair<TVector3, TVector3>> lines;
334  std::vector<double> weights;
335  for (const auto& v : fAssigned) {
336  pma::Track3D* trk = v.first.Track();
337  int vIdx = v.second;
338 
339  pma::Node3D* vtx1 = trk->Nodes()[vIdx];
340  pma::Segment3D* seg = trk->NextSegment(vtx1);
341  double segLength = seg->Length();
342  if (segLength >= fSegMinLength) {
343  pma::Node3D* vtx2 = static_cast<pma::Node3D*>(seg->Next());
344 
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;
350 
351  lines.push_back(endpoints);
352  segments.push_back(seg);
353  weights.push_back(w);
354  }
355  }
356 
357  fCenter.SetXYZ(0., 0., 0.);
358  fErr.SetXYZ(0., 0., 0.);
359 
360  TVector3 result;
361  double resultMse = pma::SolveLeastSquares3D(lines, result);
362  if (resultMse < 0.0) {
363  mf::LogWarning("pma::VtxCandidate") << "Cannot compute crossing point.";
364  return 1.0E+6;
365  }
366 
367  TVector3 pproj;
368  //double dx, dy, dz
369  double wsum = 0.0;
370  for (size_t s = 0; s < segments.size(); s++) {
371  pma::Node3D* vprev = static_cast<pma::Node3D*>(segments[s]->Prev());
372  pma::Node3D* vnext = static_cast<pma::Node3D*>(segments[s]->Next());
373 
374  pproj = pma::GetProjectionToSegment(result, vprev->Point3D(), vnext->Point3D());
375 
376  //dx = weights[s] * (result.X() - pproj.X());
377  //dy = result.Y() - pproj.Y();
378  //dz = result.Z() - pproj.Z();
379 
380  fErr[0] += weights[s] * weights[s];
381  fErr[1] += 1.0;
382  fErr[2] += 1.0;
383 
384  fCenter[0] += weights[s] * pproj.X();
385  fCenter[1] += pproj.Y();
386  fCenter[2] += pproj.Z();
387  wsum += weights[s];
388  }
389  fCenter[0] /= wsum;
390  fCenter[1] /= segments.size();
391  fCenter[2] /= segments.size();
392 
393  fErr *= 1.0 / segments.size();
394  fErr[0] = sqrt(fErr[0]);
395  fErr[1] = sqrt(fErr[1]);
396  fErr[2] = sqrt(fErr[2]);
397 
398  return resultMse;
399 }
TVector3 const & Point3D() const
Definition: PmaNode3D.h:44
double SolveLeastSquares3D(const std::vector< std::pair< TVector3, TVector3 >> &lines, TVector3 &result)
Definition: Utilities.cxx:158
std::vector< std::pair< pma::TrkCandidate, size_t > > fAssigned
virtual pma::SortedObjectBase * Next(unsigned int=0) const
Definition: SortedObjects.h:43
TVector2 GetProjectionToSegment(const TVector2 &p, const TVector2 &p0, const TVector2 &p1)
Definition: Utilities.cxx:144
MaybeLogger_< ELseverityLevel::ELsev_warning, false > LogWarning
std::vector< pma::Node3D * > const & Nodes() const noexcept
Definition: PmaTrack3D.h:274
Float_t w
Definition: plot.C:20
pma::Segment3D * NextSegment(pma::Node3D *vtx) const
double Length(void) const
Definition: PmaElement3D.h:53
double pma::VtxCandidate::ComputeMse2D ( )

Definition at line 196 of file PmaVtxCandidate.cxx.

References fAssigned, fCenter, pma::Segment3D::GetDistance2To(), pma::GetProjectionToPlane(), geo::kU, geo::kV, geo::kZ, pma::Track3D::NextSegment(), pma::Track3D::Nodes(), and geo::GeometryCore::TPC().

Referenced by Add(), MergeWith(), and VtxCandidate().

197 {
199 
200  double mse = 0.0;
201  TVector2 center2d;
202  for (const auto& t : fAssigned) {
203  pma::Track3D* trk = t.first.Track();
204  pma::Segment3D* seg = trk->NextSegment(trk->Nodes()[t.second]);
205 
206  int tpc = trk->Nodes()[t.second]->TPC();
207  int cryo = trk->Nodes()[t.second]->Cryo();
208 
209  size_t k = 0;
210  double m = 0.0;
211  auto const& tpcgeom = geom->TPC(geo::TPCID(cryo, tpc));
212  if (tpcgeom.HasPlane(geo::kU)) {
213  center2d = GetProjectionToPlane(fCenter, geo::kU, tpc, cryo);
214  m += seg->GetDistance2To(center2d, geo::kU);
215  k++;
216  }
217  if (tpcgeom.HasPlane(geo::kV)) {
218  center2d = GetProjectionToPlane(fCenter, geo::kV, tpc, cryo);
219  m += seg->GetDistance2To(center2d, geo::kV);
220  k++;
221  }
222  if (tpcgeom.HasPlane(geo::kZ)) {
223  center2d = GetProjectionToPlane(fCenter, geo::kZ, tpc, cryo);
224  m += seg->GetDistance2To(center2d, geo::kZ);
225  k++;
226  }
227  mse += m / (double)k;
228  }
229 
230  return mse / fAssigned.size();
231 }
Planes which measure V.
Definition: geo_types.h:136
std::vector< std::pair< pma::TrkCandidate, size_t > > fAssigned
Planes which measure Z direction.
Definition: geo_types.h:138
TPCGeo const & TPC(TPCID const &tpcid=tpc_zero) const
Returns the specified TPC.
Definition: GeometryCore.h:722
Planes which measure U.
Definition: geo_types.h:135
The data type to uniquely identify a TPC.
Definition: geo_types.h:381
TVector2 GetProjectionToPlane(const TVector3 &p, unsigned int plane, unsigned int tpc, unsigned int cryo)
Definition: Utilities.cxx:263
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
Definition: PmaTrack3D.h:274
pma::Segment3D * NextSegment(pma::Node3D *vtx) const
bool pma::VtxCandidate::Has ( pma::Track3D trk) const

Definition at line 26 of file PmaVtxCandidate.cxx.

References fAssigned.

Referenced by Has(), MergeWith(), and VtxCandidate().

27 {
28  for (const auto& t : fAssigned)
29  if (trk == t.first.Track()) return true;
30  return false;
31 }
std::vector< std::pair< pma::TrkCandidate, size_t > > fAssigned
bool pma::VtxCandidate::Has ( const VtxCandidate other) const

Definition at line 33 of file PmaVtxCandidate.cxx.

References fAssigned, and Has().

34 {
35  for (const auto& t : other.fAssigned)
36  if (!Has(t.first.Track())) return false;
37  return true;
38 }
bool Has(pma::Track3D *trk) const
bool pma::VtxCandidate::has ( const std::vector< int > &  v,
int  id 
) const
inlineprivate

Definition at line 77 of file PmaVtxCandidate.h.

Referenced by JoinTracks().

78  {
79  for (auto c : v)
80  if (c == id) return true;
81  return false;
82  }
bool pma::VtxCandidate::HasLoops ( ) const

Definition at line 61 of file PmaVtxCandidate.cxx.

References fAssigned, and pma::Track3D::IsAttachedTo().

Referenced by VtxCandidate().

62 {
63  for (size_t t = 0; t < fAssigned.size(); t++) {
64  pma::Track3D const* trk_t = fAssigned[t].first.Track()->GetRoot();
65  if (!trk_t) throw cet::exception("pma::VtxCandidate") << "Broken track.";
66 
67  for (size_t u = 0; u < fAssigned.size(); u++)
68  if (t != u) {
69  pma::Track3D const* trk_u = fAssigned[u].first.Track()->GetRoot();
70  if (!trk_u) throw cet::exception("pma::VtxCandidate") << "Broken track.";
71 
72  if (trk_t->IsAttachedTo(trk_u)) return true;
73  }
74  }
75  return false;
76 }
bool IsAttachedTo(pma::Track3D const *trk) const
std::vector< std::pair< pma::TrkCandidate, size_t > > fAssigned
cet::coded_exception< error, detail::translate > exception
Definition: exception.h:33
bool pma::VtxCandidate::IsAttached ( pma::Track3D trk) const

Definition at line 40 of file PmaVtxCandidate.cxx.

References fAssigned, pma::Track3D::GetRoot(), and pma::Track3D::IsAttachedTo().

Referenced by Add(), IsAttached(), MergeWith(), and VtxCandidate().

41 {
42  pma::Track3D const* rootTrk = trk->GetRoot();
43  if (!rootTrk) throw cet::exception("pma::VtxCandidate") << "Broken track.";
44 
45  for (const auto& t : fAssigned) {
46  pma::Track3D const* rootAssn = t.first.Track()->GetRoot();
47  if (!rootAssn) throw cet::exception("pma::VtxCandidate") << "Broken track.";
48 
49  if (rootTrk->IsAttachedTo(rootAssn)) return true;
50  }
51  return false;
52 }
bool IsAttachedTo(pma::Track3D const *trk) const
std::vector< std::pair< pma::TrkCandidate, size_t > > fAssigned
pma::Track3D * GetRoot()
cet::coded_exception< error, detail::translate > exception
Definition: exception.h:33
bool pma::VtxCandidate::IsAttached ( const VtxCandidate other) const

Definition at line 54 of file PmaVtxCandidate.cxx.

References fAssigned, and IsAttached().

55 {
56  for (const auto& t : other.fAssigned)
57  if (IsAttached(t.first.Track())) return true;
58  return false;
59 }
bool IsAttached(pma::Track3D *trk) const
bool pma::VtxCandidate::JoinTracks ( detinfo::DetectorPropertiesData const &  detProp,
pma::TrkCandidateColl tracks,
pma::TrkCandidateColl src 
)

Definition at line 401 of file PmaVtxCandidate.cxx.

References pma::Track3D::AttachBackTo(), pma::Track3D::AttachTo(), pma::Track3D::CanFlip(), pma::Dist2(), pma::TrkCandidateColl::erase_at(), f, fAssigned, fCenter, pma::Track3D::Flip(), fMse, fMse2D, pma::Node3D::GetBranches(), pma::Track3D::GetBranches(), pma::TrkCandidateColl::getCandidateIndex(), pma::GetProjectionToSegment(), pma::Track3D::GetRoot(), pma::GetSegmentProjVector(), pma::TrkCandidateColl::getTreeCopy(), has(), pma::Track3D::InsertNode(), kMinDistToNode, pma::Track3D::MakeProjection(), pma::SortedBranchBase::Next(), pma::SortedBranchBase::NextCount(), pma::Track3D::NextSegment(), pma::Track3D::Nodes(), pma::Segment3D::Parent(), pma::Node3D::Point3D(), pma::SortedObjectBase::Prev(), proj, pma::TrkCandidateColl::push_back(), pma::Node3D::SetPoint3D(), pma::TrkCandidateColl::setTreeIds(), pma::TrkCandidateColl::size(), pma::Track3D::size(), pma::Track3D::Split(), pma::TrkCandidateColl::tracks(), tracksJoined, and pma::Track3D::TuneFullTree().

Referenced by Size().

404 {
405  if (tracksJoined) {
406  mf::LogError("pma::VtxCandidate") << "Tracks already attached to the vertex.";
407  return false;
408  }
409  tracksJoined = true;
410 
411  mf::LogVerbatim("pma::VtxCandidate")
412  << "JoinTracks (" << fAssigned.size() << ") at:"
413  << " vx:" << fCenter.X() << " vy:" << fCenter.Y() << " vz:" << fCenter.Z();
414 
415  for (auto const& c : fAssigned) {
416  size_t t = 0;
417  while (t < src.size()) {
418  if (c.first.Track() == src[t].Track()) { // move from "src" to "tracks"
419  tracks.push_back(src[t]);
420  src.erase_at(t);
421  break;
422  }
423  else
424  t++;
425  }
426  }
427  // all involved tracks are in the same container, so:
428  tracks.setTreeIds();
429  for (auto& c : fAssigned) // update ids
430  for (auto const& t : tracks.tracks())
431  if (c.first.Track() == t.Track()) {
432  c.first.SetTreeId(t.TreeId());
433  break;
434  }
435 
436  // backup in case of fitting problems
437  std::vector<int> treeIds;
438  pma::TrkCandidateColl backupTracks;
439 
440  pma::Node3D* vtxCenter = 0;
441  bool hasInnerCenter = false;
442  size_t nOK = 0;
443  for (size_t i = 0; i < fAssigned.size(); i++) {
444  mf::LogVerbatim("pma::VtxCandidate") << "----------> track #" << i;
445 
446  pma::Track3D* trk = fAssigned[i].first.Track();
447  int key = fAssigned[i].first.Key();
448  int tid = fAssigned[i].first.TreeId();
449  size_t idx = fAssigned[i].second;
450 
451  mf::LogVerbatim("pma::VtxCandidate") << " track tid:" << tid << ", size:" << trk->size()
452  << " (nodes:" << trk->Nodes().size() << ")";
453 
454  if (!has(treeIds, tid)) // backup in case of fitting problems
455  {
456  treeIds.push_back(tid);
457  int rootIdx = tracks.getCandidateIndex(trk->GetRoot());
458  if (rootIdx >= 0)
459  tracks.getTreeCopy(backupTracks, rootIdx);
460  else
461  mf::LogError("pma::VtxCandidate") << "Root of the tree not found in tracks collection.";
462  }
463 
464  TVector3 p0(trk->Nodes()[idx]->Point3D());
465  TVector3 p1(trk->Nodes()[idx + 1]->Point3D());
466 
467  int tpc0 = trk->Nodes()[idx]->TPC();
468  int tpc1 = trk->Nodes()[idx + 1]->TPC();
469 
470  int cryo0 = trk->Nodes()[idx]->Cryo();
471  int cryo1 = trk->Nodes()[idx + 1]->Cryo();
472 
473  double d0 = sqrt(pma::Dist2(p0, fCenter));
474  double d1 = sqrt(pma::Dist2(p1, fCenter));
475  double ds = sqrt(pma::Dist2(p0, p1));
476  double f = pma::GetSegmentProjVector(fCenter, p0, p1);
477  TVector3 proj = pma::GetProjectionToSegment(fCenter, p0, p1);
478 
479  if ((idx == 0) && (f * ds <= kMinDistToNode)) {
480  if (i == 0) {
481  mf::LogVerbatim("pma::VtxCandidate") << " new at front";
482  vtxCenter = trk->Nodes().front();
483  vtxCenter->SetPoint3D(fCenter);
484  nOK++;
485  }
486  else {
487  mf::LogVerbatim("pma::VtxCandidate") << " front to center";
488  if (trk->AttachTo(vtxCenter)) nOK++;
489  }
490  }
491  else if ((idx + 2 == trk->Nodes().size()) && ((1.0 - f) * ds <= kMinDistToNode)) {
492  if (i == 0) {
493  if (trk->CanFlip()) {
494  mf::LogVerbatim("pma::VtxCandidate") << " flip trk to make new center";
495  trk->Flip();
496  vtxCenter = trk->Nodes().front();
497  }
498  else {
499  mf::LogVerbatim("pma::VtxCandidate") << " new center at the endpoint";
500  vtxCenter = trk->Nodes().back();
501  }
502  vtxCenter->SetPoint3D(fCenter);
503  nOK++;
504  }
505  else {
506  if (vtxCenter->Prev() && trk->CanFlip()) {
507  mf::LogVerbatim("pma::VtxCandidate") << " flip trk to attach to inner";
508  trk->Flip();
509  if (trk->AttachTo(vtxCenter)) nOK++;
510  }
511  else {
512  mf::LogVerbatim("pma::VtxCandidate") << " endpoint to center";
513  if (trk->AttachBackTo(vtxCenter)) nOK++;
514  }
515  }
516  }
517  else {
518  bool canFlipPrev = true;
519  if (vtxCenter && vtxCenter->Prev()) {
520  pma::Segment3D* seg = static_cast<pma::Segment3D*>(vtxCenter->Prev());
521  if (seg->Parent()->NextSegment(vtxCenter))
522  canFlipPrev = false;
523  else
524  canFlipPrev = seg->Parent()->CanFlip();
525  }
526 
527  if (hasInnerCenter || !canFlipPrev) {
528  mf::LogVerbatim("pma::VtxCandidate") << " split track";
529 
530  if ((f >= 0.0F) && (f <= 1.0) && (f * ds > kMinDistToNode) &&
531  ((1.0 - f) * ds > kMinDistToNode)) {
532  mf::LogVerbatim("pma::VtxCandidate") << " add center inside segment";
533 
534  int tpc, cryo;
535  if (f < 0.5) {
536  tpc = tpc0;
537  cryo = cryo0;
538  }
539  else {
540  tpc = tpc1;
541  cryo = cryo1;
542  }
543 
544  trk->InsertNode(detProp, fCenter, ++idx, tpc, cryo);
545  }
546  else {
547  if (d1 < d0) {
548  mf::LogVerbatim("pma::VtxCandidate") << " add center at end of segment";
549  ++idx;
550  }
551  else {
552  mf::LogVerbatim("pma::VtxCandidate") << " center at start of segment - no action";
553  }
554  }
555 
556  pma::Track3D* t0 = trk->Split(detProp, idx); // makes both tracks attached to each other
557 
558  if (t0) {
559  mf::LogVerbatim("pma::VtxCandidate")
560  << " trk size:" << trk->size() << " (nodes:" << trk->Nodes().size() << ")";
561 
562  trk->MakeProjection();
563 
564  mf::LogVerbatim("pma::VtxCandidate")
565  << " t0 size:" << t0->size() << " (nodes:" << t0->Nodes().size() << ")";
566 
567  t0->MakeProjection();
568 
569  tracks.tracks().emplace_back(t0, key, tid);
570  if (i == 0) {
571  mf::LogVerbatim("pma::VtxCandidate") << " center at trk0 back";
572  vtxCenter = trk->Nodes().front();
573  nOK += 2;
574  }
575  else {
576  mf::LogVerbatim("pma::VtxCandidate") << " attach trk to trk0";
577  if (trk->AttachTo(vtxCenter)) nOK += 2;
578  }
579  }
580  mf::LogVerbatim("pma::VtxCandidate") << " done";
581  }
582  else {
583  mf::LogVerbatim("pma::VtxCandidate") << " inner center";
584  hasInnerCenter = true;
585 
586  if ((f >= 0.0F) && (f <= 1.0) && (f * ds > kMinDistToNode) &&
587  ((1.0 - f) * ds > kMinDistToNode)) {
588  mf::LogVerbatim("pma::VtxCandidate") << " add center inside segment";
589 
590  int tpc, cryo;
591  if (f < 0.5) {
592  tpc = tpc0;
593  cryo = cryo0;
594  }
595  else {
596  tpc = tpc1;
597  cryo = cryo1;
598  }
599 
600  trk->InsertNode(detProp, fCenter, ++idx, tpc, cryo);
601  }
602  else {
603  if (d1 < d0) {
604  mf::LogVerbatim("pma::VtxCandidate") << " add center at end of segment";
605  ++idx;
606  }
607  else {
608  mf::LogVerbatim("pma::VtxCandidate") << " center at start of segment - no action";
609  }
610  }
611 
612  pma::Node3D* innerCenter = trk->Nodes()[idx];
613  if (i > 0) // already has vtxCenter
614  {
615  // prepare for prev...
616  pma::Track3D* rootBranch = 0;
617  pma::Segment3D* seg = static_cast<pma::Segment3D*>(vtxCenter->Prev());
618  if (seg) {
619  rootBranch = seg->Parent();
620  rootBranch->Flip();
621  }
622 
623  // ...but nexts reattached first, then...
624  auto branches = vtxCenter->GetBranches();
625 
626  for (size_t j = 0; j < branches.size(); ++j) {
627  if (branches[j]->AttachTo(innerCenter, true)) {}
628  }
629  vtxCenter = innerCenter; // vtxCenter is deleted after full reattach
630  }
631  else {
632  vtxCenter = innerCenter; // vtx from inner center
633  } // set vtxCenter on i == 0
634 
635  nOK++;
636 
637  mf::LogVerbatim("pma::VtxCandidate") << " done";
638  }
639  }
640  }
641 
642  bool result = false;
643  if (vtxCenter) {
644  pma::Segment3D* rootSeg = 0;
645  if (vtxCenter->NextCount())
646  rootSeg = static_cast<pma::Segment3D*>(vtxCenter->Next(0));
647  else if (vtxCenter->Prev())
648  rootSeg = static_cast<pma::Segment3D*>(vtxCenter->Prev());
649  else
650  throw cet::exception("pma::VtxCandidate") << "Vertex with no segments attached.";
651 
652  pma::Track3D* rootTrk = rootSeg->Parent()->GetRoot();
653  if (!rootTrk) rootTrk = rootSeg->Parent();
654 
655  std::vector<pma::Track3D const*> branchesToRemove; // change to a more simple check
656  bool noLoops = rootTrk->GetBranches(branchesToRemove); // if there are loops
657 
658  bool tuneOK = true;
659  if (noLoops && (nOK > 1)) {
660  fAssigned.clear();
661  fCenter = vtxCenter->Point3D();
662  fMse = 0.0;
663  fMse2D = 0.0;
664 
665  double g = rootTrk->TuneFullTree();
666 
667  if (g > -2.0) // -1.0: high value of g; -2.0: inf. value of g.
668  {
669  result = true; // all OK, new vertex added
670  }
671  else {
672  tuneOK = false; // inf. g, remove tracks
673  }
674  }
675 
676  if (noLoops && tuneOK) {
677  mf::LogVerbatim("pma::VtxCandidate") << "remove backup tracks";
678  for (auto& c : backupTracks.tracks())
679  c.DeleteTrack();
680  }
681  else {
682  mf::LogVerbatim("pma::VtxCandidate") << "restore tracks from backup....";
683  for (int tid : treeIds) {
684  size_t t = 0;
685  while (t < tracks.size()) {
686  if (tracks[t].TreeId() == tid) {
687  tracks[t].DeleteTrack();
688  tracks.erase_at(t);
689  }
690  else
691  t++;
692  }
693  }
694  for (const auto& c : backupTracks.tracks())
695  tracks.push_back(c);
696  mf::LogVerbatim("pma::VtxCandidate") << " done";
697  }
698  }
699  else
700  mf::LogError("pma::VtxCandidate") << "Cannot create common vertex";
701 
702  return result;
703 }
void MakeProjection()
code to link reconstructed objects back to the MC truth information
MaybeLogger_< ELseverityLevel::ELsev_info, true > LogVerbatim
TVector3 const & Point3D() const
Definition: PmaNode3D.h:44
size_t size() const
pma::Track3D * getTreeCopy(pma::TrkCandidateColl &dst, size_t trkIdx, bool isRoot=true)
bool AttachTo(pma::Node3D *vStart, bool noFlip=false)
double Dist2(const TVector2 &v1, const TVector2 &v2)
Definition: Utilities.cxx:39
void erase_at(size_t pos)
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
Definition: SortedObjects.h:84
std::vector< pma::Track3D * > GetBranches() const
Definition: PmaNode3D.cxx:443
MaybeLogger_< ELseverityLevel::ELsev_error, false > LogError
double TuneFullTree(double eps=0.001, double gmax=50.0)
TFile f
Definition: plotHisto.C:6
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.
Definition: PmaTrack3D.cxx:644
TVector2 GetProjectionToSegment(const TVector2 &p, const TVector2 &p0, const TVector2 &p1)
Definition: Utilities.cxx:144
pma::Track3D * GetRoot()
static constexpr double kMinDistToNode
int getCandidateIndex(pma::Track3D const *candidate) const
bool SetPoint3D(const TVector3 &p3d)
Definition: PmaNode3D.cxx:162
void InsertNode(detinfo::DetectorPropertiesData const &detProp, TVector3 const &p3d, size_t at_idx, unsigned int tpc, unsigned int cryo)
virtual pma::SortedObjectBase * Next(unsigned int index=0) const
Definition: SortedObjects.h:77
double GetSegmentProjVector(const TVector2 &p, const TVector2 &p0, const TVector2 &p1)
Definition: Utilities.cxx:112
Float_t proj
Definition: plot.C:35
std::vector< pma::Node3D * > const & Nodes() const noexcept
Definition: PmaTrack3D.h:274
virtual pma::SortedObjectBase * Prev(void) const
Definition: SortedObjects.h:42
bool Flip(const detinfo::DetectorPropertiesData &detProp, std::vector< pma::Track3D * > &allTracks)
Definition: PmaTrack3D.cxx:527
size_t size() const
Definition: PmaTrack3D.h:89
pma::Track3D * Parent(void) const
Definition: PmaSegment3D.h:66
bool AttachBackTo(pma::Node3D *vStart)
pma::Segment3D * NextSegment(pma::Node3D *vtx) const
void push_back(const TrkCandidate &trk)
std::vector< TrkCandidate > const & tracks() const
cet::coded_exception< error, detail::translate > exception
Definition: exception.h:33
bool has(const std::vector< int > &v, int id) const
double pma::VtxCandidate::MaxAngle ( double  minLength = 0.0) const

Definition at line 243 of file PmaVtxCandidate.cxx.

References fAssigned, pma::Track3D::Nodes(), and pma::Node3D::Point3D().

Referenced by VtxCandidate().

244 {
245  TVector3 dir_i;
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();
250  if (l > max_l) {
251  max_l = l;
252  max_l_idx = i;
253  pma::Track3D* trk_i = fAssigned[i].first.Track();
254  pma::Node3D* vtx_i0 = trk_i->Nodes()[fAssigned[i].second];
255  pma::Node3D* vtx_i1 = trk_i->Nodes()[fAssigned[i].second + 1];
256  dir_i = vtx_i1->Point3D() - vtx_i0->Point3D();
257  dir_i *= 1.0 / dir_i.Mag();
258  }
259  }
260 
261  double a, min = 1.0;
262  for (size_t j = 0; j < fAssigned.size(); j++)
263  if ((j != max_l_idx) && (fAssigned[j].first.Track()->Length() > minLength)) {
264  pma::Track3D* trk_j = fAssigned[j].first.Track();
265  pma::Node3D* vtx_j0 = trk_j->Nodes()[fAssigned[j].second];
266  pma::Node3D* vtx_j1 = trk_j->Nodes()[fAssigned[j].second + 1];
267  TVector3 dir_j = vtx_j1->Point3D() - vtx_j0->Point3D();
268  dir_j *= 1.0 / dir_j.Mag();
269  a = fabs(dir_i * dir_j);
270  if (a < min) min = a;
271  }
272 
273  return 180.0 * acos(min) / TMath::Pi();
274 }
TVector3 const & Point3D() const
Definition: PmaNode3D.h:44
std::vector< std::pair< pma::TrkCandidate, size_t > > fAssigned
std::vector< pma::Node3D * > const & Nodes() const noexcept
Definition: PmaTrack3D.h:274
bool pma::VtxCandidate::MergeWith ( const VtxCandidate other)

Definition at line 276 of file PmaVtxCandidate.cxx.

References Compute(), ComputeMse2D(), d, pma::Dist2(), fAssigned, fCenter, fMse, fMse2D, Has(), IsAttached(), Size(), and Test().

Referenced by Size().

277 {
278  double d = sqrt(pma::Dist2(fCenter, other.fCenter));
279  if (d > 10.0) {
280  mf::LogVerbatim("pma::VtxCandidate") << "too far..";
281  return false;
282  }
283 
284  double dw = Test(other);
285 
286  size_t ntrk = 0;
287  for (const auto& t : other.fAssigned) {
288  if (IsAttached(t.first.Track())) {
289  mf::LogVerbatim("pma::VtxCandidate") << "already attached..";
290  return false;
291  }
292  if (!Has(t.first.Track())) {
293  fAssigned.push_back(t);
294  ntrk++;
295  }
296  }
297  if (ntrk) {
298  double mse0 = fMse, mse1 = other.fMse;
299  mf::LogVerbatim("pma::VtxCandidate")
300  << "try: " << d << " mse0:" << sqrt(mse0) << " mse1:" << sqrt(mse1);
301  //std::cout << "try: " << d << " mse0:" << sqrt(mse0) << " mse1:" << sqrt(mse1) << std::endl;
302 
303  double mse = Compute();
304  mf::LogVerbatim("pma::VtxCandidate")
305  << "out: " << Size() << " mse:" << sqrt(mse) << " dw:" << dw;
306  //std::cout << "out: " << Size() << " mse:" << sqrt(mse) << " dw:" << dw << std::endl;
307 
308  if (mse < 1.0) // kMaxDistToTrack * kMaxDistToTrack)
309  {
310  fMse = mse;
311  fMse2D = ComputeMse2D();
312  return true;
313  }
314  else {
315  mf::LogVerbatim("pma::VtxCandidate") << "high mse..";
316  while (ntrk--) {
317  fAssigned.pop_back();
318  }
319  fMse = Compute();
320  fMse2D = ComputeMse2D();
321  return false;
322  }
323  }
324  else {
325  mf::LogVerbatim("pma::VtxCandidate") << "no tracks..";
326  return false;
327  }
328 }
MaybeLogger_< ELseverityLevel::ELsev_info, true > LogVerbatim
bool Has(pma::Track3D *trk) const
double Dist2(const TVector2 &v1, const TVector2 &v2)
Definition: Utilities.cxx:39
double Test(const VtxCandidate &other) const
std::vector< std::pair< pma::TrkCandidate, size_t > > fAssigned
size_t Size() const
Float_t d
Definition: plot.C:235
bool IsAttached(pma::Track3D *trk) const
double pma::VtxCandidate::Mse ( ) const
inline

Definition at line 68 of file PmaVtxCandidate.h.

References fMse.

Referenced by pma::PMAlgVertexing::firstPassCandidates(), and pma::PMAlgVertexing::secondPassCandidates().

68 { return fMse; }
double pma::VtxCandidate::Mse2D ( ) const
inline

Definition at line 69 of file PmaVtxCandidate.h.

References fMse2D.

69 { return fMse2D; }
size_t pma::VtxCandidate::Size ( ) const
inline

Definition at line 56 of file PmaVtxCandidate.h.

References Compute(), fAssigned, JoinTracks(), and MergeWith().

Referenced by MergeWith().

56 { return fAssigned.size(); }
std::vector< std::pair< pma::TrkCandidate, size_t > > fAssigned
size_t pma::VtxCandidate::Size ( double  minLength) const

Definition at line 78 of file PmaVtxCandidate.cxx.

References fAssigned, and n.

79 {
80  size_t n = 0;
81  for (auto const& c : fAssigned)
82  if (c.first.Track()->Length() > minLength) n++;
83  return n;
84 }
std::vector< std::pair< pma::TrkCandidate, size_t > > fAssigned
Char_t n[5]
double pma::VtxCandidate::Test ( const VtxCandidate other) const

Definition at line 233 of file PmaVtxCandidate.cxx.

References fCenter, and fErr.

Referenced by MergeWith(), and VtxCandidate().

234 {
235  double dx = fCenter[0] - other.fCenter[0];
236  double dy = fCenter[1] - other.fCenter[1];
237  double dz = fCenter[2] - other.fCenter[2];
238  double dw = fErr[0] * other.fErr[0] * dx * dx + fErr[1] * other.fErr[1] * dy * dy +
239  fErr[2] * other.fErr[2] * dz * dz;
240  return sqrt(dw);
241 }
std::pair<pma::Track3D*, size_t> pma::VtxCandidate::Track ( size_t  i) const
inline

Definition at line 71 of file PmaVtxCandidate.h.

References fAssigned.

72  {
73  return std::pair<pma::Track3D*, size_t>(fAssigned[i].first.Track(), fAssigned[i].second);
74  }
std::vector< std::pair< pma::TrkCandidate, size_t > > fAssigned

Member Data Documentation

std::vector<std::pair<pma::TrkCandidate, size_t> > pma::VtxCandidate::fAssigned
private
TVector3 pma::VtxCandidate::fCenter
private

Definition at line 87 of file PmaVtxCandidate.h.

Referenced by Add(), Center(), Compute(), ComputeMse2D(), JoinTracks(), MergeWith(), and Test().

TVector3 pma::VtxCandidate::fErr
private

Definition at line 87 of file PmaVtxCandidate.h.

Referenced by Compute(), and Test().

double pma::VtxCandidate::fMse
private

Definition at line 85 of file PmaVtxCandidate.h.

Referenced by Add(), JoinTracks(), MergeWith(), and Mse().

double pma::VtxCandidate::fMse2D
private

Definition at line 85 of file PmaVtxCandidate.h.

Referenced by Add(), JoinTracks(), MergeWith(), and Mse2D().

double pma::VtxCandidate::fSegMinLength
private

Definition at line 85 of file PmaVtxCandidate.h.

Referenced by Add(), and Compute().

constexpr double pma::VtxCandidate::kMaxDistToTrack {4.0}
static

Definition at line 26 of file PmaVtxCandidate.h.

Referenced by Add().

constexpr double pma::VtxCandidate::kMinDistToNode {2.0}
static

Definition at line 27 of file PmaVtxCandidate.h.

Referenced by JoinTracks().

bool pma::VtxCandidate::tracksJoined
private

Definition at line 84 of file PmaVtxCandidate.h.

Referenced by JoinTracks().


The documentation for this class was generated from the following files: