LArSoft  v07_13_02
Liquid Argon Software toolkit - http://larsoft.org/
Public Member Functions | Static Public Attributes | Private Member Functions | Private Attributes | List of all members
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 (void) const
 
bool Add (const pma::TrkCandidate &trk)
 
double ComputeMse2D (void)
 
double Test (const VtxCandidate &other) const
 
double MaxAngle (double minLength=0.0) const
 
size_t Size (void) const
 
size_t Size (double minLength) const
 
bool MergeWith (const VtxCandidate &other)
 
double Compute (void)
 
bool JoinTracks (pma::TrkCandidateColl &tracks, pma::TrkCandidateColl &src)
 
const TVector3 & Center (void) const
 
double Mse (void) const
 
double Mse2D (void) const
 
std::pair< pma::Track3D *, size_t > Track (size_t i) const
 

Static Public Attributes

static const double kMaxDistToTrack = 4.0
 
static const 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 22 of file PmaVtxCandidate.h.

Constructor & Destructor Documentation

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

Definition at line 28 of file PmaVtxCandidate.h.

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

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

Member Function Documentation

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

Definition at line 87 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().

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

Definition at line 63 of file PmaVtxCandidate.h.

References fCenter.

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

Definition at line 350 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(), s, pma::SolveLeastSquares3D(), and w.

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

351 {
352  std::vector< pma::Segment3D* > segments;
353  std::vector< std::pair<TVector3, TVector3> > lines;
354  std::vector< double > weights;
355  for (const auto & v : fAssigned)
356  {
357  pma::Track3D* trk = v.first.Track();
358  int vIdx = v.second;
359 
360  pma::Node3D* vtx1 = trk->Nodes()[vIdx];
361  pma::Segment3D* seg = trk->NextSegment(vtx1);
362  double segLength = seg->Length();
363  if (segLength >= fSegMinLength)
364  {
365  pma::Node3D* vtx2 = static_cast< pma::Node3D* >(seg->Next(0));
366 
367  std::pair<TVector3, TVector3> endpoints(vtx1->Point3D(), vtx2->Point3D());
368  double dy = endpoints.first.Y() - endpoints.second.Y();
369  double fy_norm = asin(fabs(dy) / segLength) / (0.5 * TMath::Pi());
370  double w = 1.0 - pow(fy_norm - 1.0, 12);
371  if (w < 0.3) w = 0.3;
372 
373  lines.push_back(endpoints);
374  segments.push_back(seg);
375  weights.push_back(w);
376  }
377  }
378 
379  fCenter.SetXYZ(0., 0., 0.); fErr.SetXYZ(0., 0., 0.);
380 
381  TVector3 result;
382  double resultMse = pma::SolveLeastSquares3D(lines, result);
383  if (resultMse < 0.0)
384  {
385  mf::LogWarning("pma::VtxCandidate") << "Cannot compute crossing point.";
386  return 1.0E+6;
387  }
388 
389  TVector3 pproj;
390  //double dx, dy, dz
391  double wsum = 0.0;
392  for (size_t s = 0; s < segments.size(); s++)
393  {
394  pma::Node3D* vprev = static_cast< pma::Node3D* >(segments[s]->Prev());
395  pma::Node3D* vnext = static_cast< pma::Node3D* >(segments[s]->Next(0));
396 
397  pproj = pma::GetProjectionToSegment(result, vprev->Point3D(), vnext->Point3D());
398 
399  //dx = weights[s] * (result.X() - pproj.X());
400  //dy = result.Y() - pproj.Y();
401  //dz = result.Z() - pproj.Z();
402 
403  fErr[0] += weights[s] * weights[s];
404  fErr[1] += 1.0;
405  fErr[2] += 1.0;
406 
407  fCenter[0] += weights[s] * pproj.X();
408  fCenter[1] += pproj.Y();
409  fCenter[2] += pproj.Z();
410  wsum += weights[s];
411  }
412  fCenter[0] /= wsum;
413  fCenter[1] /= segments.size();
414  fCenter[2] /= segments.size();
415 
416  fErr *= 1.0 / segments.size();
417  fErr[0] = sqrt(fErr[0]);
418  fErr[1] = sqrt(fErr[1]);
419  fErr[2] = sqrt(fErr[2]);
420 
421  return resultMse;
422 }
s
Float_t s
Definition: plot.C:23
pma::Node3D::Point3D
TVector3 const & Point3D(void) const
Definition: PmaNode3D.h:33
pma::Track3D::Nodes
std::vector< pma::Node3D * > const & Nodes(void) const
Definition: PmaTrack3D.h:232
pma::GetProjectionToSegment
TVector2 GetProjectionToSegment(const TVector2 &p, const TVector2 &p0, const TVector2 &p1)
Definition: Utilities.cxx:156
pma::VtxCandidate::fAssigned
std::vector< std::pair< pma::TrkCandidate, size_t > > fAssigned
Definition: PmaVtxCandidate.h:81
pma::SolveLeastSquares3D
double SolveLeastSquares3D(const std::vector< std::pair< TVector3, TVector3 > > &lines, TVector3 &result)
Definition: Utilities.cxx:176
mf::LogWarning
MaybeLogger_< ELseverityLevel::ELsev_warning, false > LogWarning
Definition: MessageLogger.h:213
w
Float_t w
Definition: plot.C:23
pma::SortedObjectBase::Next
virtual pma::SortedObjectBase * Next(unsigned int index=0) const
Definition: SortedObjects.h:45
pma::Track3D::NextSegment
pma::Segment3D * NextSegment(pma::Node3D *vtx) const
Definition: PmaTrack3D.cxx:966
pma::Element3D::Length
double Length(void) const
Definition: PmaElement3D.h:49
double pma::VtxCandidate::ComputeMse2D ( void  )

Definition at line 208 of file PmaVtxCandidate.cxx.

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

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

209 {
210  art::ServiceHandle< geo::Geometry > geom;
211 
212  double mse = 0.0;
213  TVector2 center2d;
214  for (const auto & t : fAssigned)
215  {
216  pma::Track3D* trk = t.first.Track();
217  pma::Segment3D* seg = trk->NextSegment(trk->Nodes()[t.second]);
218 
219  int tpc = trk->Nodes()[t.second]->TPC();
220  int cryo = trk->Nodes()[t.second]->Cryo();
221 
222  size_t k = 0;
223  double m = 0.0;
224  if (geom->TPC(tpc, cryo).HasPlane(geo::kU))
225  {
226  center2d = GetProjectionToPlane(fCenter, geo::kU, tpc, cryo);
227  m += seg->GetDistance2To(center2d, geo::kU); k++;
228  }
229  if (geom->TPC(tpc, cryo).HasPlane(geo::kV))
230  {
231  center2d = GetProjectionToPlane(fCenter, geo::kV, tpc, cryo);
232  m += seg->GetDistance2To(center2d, geo::kV); k++;
233  }
234  if (geom->TPC(tpc, cryo).HasPlane(geo::kZ))
235  {
236  center2d = GetProjectionToPlane(fCenter, geo::kZ, tpc, cryo);
237  m += seg->GetDistance2To(center2d, geo::kZ); k++;
238  }
239  mse += m / (double)k;
240  }
241 
242  return mse / fAssigned.size();
243 }
geo::TPCGeo::HasPlane
bool HasPlane(unsigned int iplane) const
Returns whether a plane with index iplane is present in this TPC.
Definition: TPCGeo.h:155
geo::kV
Planes which measure V.
Definition: geo_types.h:77
pma::Track3D::Nodes
std::vector< pma::Node3D * > const & Nodes(void) const
Definition: PmaTrack3D.h:232
geo::kZ
Planes which measure Z direction.
Definition: geo_types.h:79
geo::kU
Planes which measure U.
Definition: geo_types.h:76
pma::VtxCandidate::fAssigned
std::vector< std::pair< pma::TrkCandidate, size_t > > fAssigned
Definition: PmaVtxCandidate.h:81
pma::GetProjectionToPlane
TVector2 GetProjectionToPlane(const TVector3 &p, unsigned int plane, unsigned int tpc, unsigned int cryo)
Definition: Utilities.cxx:280
pma::Segment3D::GetDistance2To
double GetDistance2To(const TVector3 &p3d) const override
Distance [cm] from the 3D segment to the point 3D.
Definition: PmaSegment3D.cxx:24
geo::GeometryCore::TPC
TPCGeo const & TPC(unsigned int const tpc=0, unsigned int const cstat=0) const
Returns the specified TPC.
Definition: GeometryCore.h:2393
pma::Track3D::NextSegment
pma::Segment3D * NextSegment(pma::Node3D *vtx) const
Definition: PmaTrack3D.cxx:966
bool pma::VtxCandidate::Has ( pma::Track3D *  trk) const

Definition at line 24 of file PmaVtxCandidate.cxx.

References fAssigned.

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

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

Definition at line 31 of file PmaVtxCandidate.cxx.

References fAssigned, and Has().

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

Definition at line 73 of file PmaVtxCandidate.h.

Referenced by JoinTracks().

74  {
75  for (auto c : v) if (c == id) return true;
76  return false;
77  }
bool pma::VtxCandidate::HasLoops ( void  ) const

Definition at line 60 of file PmaVtxCandidate.cxx.

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

Referenced by VtxCandidate().

61 {
62  for (size_t t = 0; t < fAssigned.size(); t++)
63  {
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  {
70  pma::Track3D const * trk_u = fAssigned[u].first.Track()->GetRoot();
71  if (!trk_u) throw cet::exception("pma::VtxCandidate") << "Broken track.";
72 
73  if (trk_t->IsAttachedTo(trk_u)) return true;
74  }
75  }
76  return false;
77 }
pma::Track3D::IsAttachedTo
bool IsAttachedTo(pma::Track3D const *trk) const
Definition: PmaTrack3D.cxx:1752
pma::VtxCandidate::fAssigned
std::vector< std::pair< pma::TrkCandidate, size_t > > fAssigned
Definition: PmaVtxCandidate.h:81
fhicl::exception
cet::coded_exception< error, detail::translate > exception
Definition: exception.h:33
bool pma::VtxCandidate::IsAttached ( pma::Track3D *  trk) const

Definition at line 38 of file PmaVtxCandidate.cxx.

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

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

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

Definition at line 53 of file PmaVtxCandidate.cxx.

References fAssigned, and IsAttached().

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

Definition at line 424 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::Track3D::size(), pma::TrkCandidateColl::size(), pma::Track3D::Split(), pma::TrkCandidateColl::tracks(), tracksJoined, and pma::Track3D::TuneFullTree().

Referenced by Size().

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

Definition at line 256 of file PmaVtxCandidate.cxx.

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

Referenced by VtxCandidate().

257 {
258  TVector3 dir_i;
259  size_t max_l_idx = 0;
260  double l, max_l = 0.0;
261  for (size_t i = 0; i < fAssigned.size() - 1; i++)
262  {
263  l = fAssigned[i].first.Track()->Length();
264  if (l > max_l)
265  {
266  max_l = l; max_l_idx = i;
267  pma::Track3D* trk_i = fAssigned[i].first.Track();
268  pma::Node3D* vtx_i0 = trk_i->Nodes()[fAssigned[i].second];
269  pma::Node3D* vtx_i1 = trk_i->Nodes()[fAssigned[i].second + 1];
270  dir_i = vtx_i1->Point3D() - vtx_i0->Point3D();
271  dir_i *= 1.0 / dir_i.Mag();
272  }
273  }
274 
275  double a, min = 1.0;
276  for (size_t j = 0; j < fAssigned.size(); j++)
277  if ((j != max_l_idx) && (fAssigned[j].first.Track()->Length() > minLength))
278  {
279  pma::Track3D* trk_j = fAssigned[j].first.Track();
280  pma::Node3D* vtx_j0 = trk_j->Nodes()[fAssigned[j].second];
281  pma::Node3D* vtx_j1 = trk_j->Nodes()[fAssigned[j].second + 1];
282  TVector3 dir_j = vtx_j1->Point3D() - vtx_j0->Point3D();
283  dir_j *= 1.0 / dir_j.Mag();
284  a = fabs(dir_i * dir_j);
285  if (a < min) min = a;
286  }
287 
288  return 180.0 * acos(min) / TMath::Pi();
289 }
pma::Node3D::Point3D
TVector3 const & Point3D(void) const
Definition: PmaNode3D.h:33
pma::Track3D::Nodes
std::vector< pma::Node3D * > const & Nodes(void) const
Definition: PmaTrack3D.h:232
pma::VtxCandidate::fAssigned
std::vector< std::pair< pma::TrkCandidate, size_t > > fAssigned
Definition: PmaVtxCandidate.h:81
min
Int_t min
Definition: plot.C:26
bool pma::VtxCandidate::MergeWith ( const VtxCandidate &  other)

Definition at line 291 of file PmaVtxCandidate.cxx.

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

Referenced by Size().

292 {
293  double d = sqrt( pma::Dist2(fCenter, other.fCenter) );
294  if (d > 10.0)
295  {
296  mf::LogVerbatim("pma::VtxCandidate") << "too far..";
297  return false;
298  }
299 
300  double dw = Test(other);
301 
302  size_t ntrk = 0;
303  for (const auto & t : other.fAssigned)
304  {
305  if (IsAttached(t.first.Track()))
306  {
307  mf::LogVerbatim("pma::VtxCandidate") << "already attached..";
308  return false;
309  }
310  if (!Has(t.first.Track()))
311  {
312  fAssigned.push_back(t);
313  ntrk++;
314  }
315  }
316  if (ntrk)
317  {
318  double mse0 = fMse, mse1 = other.fMse;
319  mf::LogVerbatim("pma::VtxCandidate")
320  << "try: " << d << " mse0:" << sqrt(mse0) << " mse1:" << sqrt(mse1);
321  //std::cout << "try: " << d << " mse0:" << sqrt(mse0) << " mse1:" << sqrt(mse1) << std::endl;
322 
323  double mse = Compute();
324  mf::LogVerbatim("pma::VtxCandidate")
325  << "out: " << Size() << " mse:" << sqrt(mse) << " dw:" << dw;
326  //std::cout << "out: " << Size() << " mse:" << sqrt(mse) << " dw:" << dw << std::endl;
327 
328  if (mse < 1.0) // kMaxDistToTrack * kMaxDistToTrack)
329  {
330  fMse = mse;
331  fMse2D = ComputeMse2D();
332  return true;
333  }
334  else
335  {
336  mf::LogVerbatim("pma::VtxCandidate") << "high mse..";
337  while (ntrk--) { fAssigned.pop_back(); }
338  fMse = Compute();
339  fMse2D = ComputeMse2D();
340  return false;
341  }
342  }
343  else
344  {
345  mf::LogVerbatim("pma::VtxCandidate") << "no tracks..";
346  return false;
347  }
348 }
mf::LogVerbatim
MaybeLogger_< ELseverityLevel::ELsev_info, true > LogVerbatim
Definition: MessageLogger.h:219
pma::VtxCandidate::ComputeMse2D
double ComputeMse2D(void)
Definition: PmaVtxCandidate.cxx:208
pma::VtxCandidate::Has
bool Has(pma::Track3D *trk) const
Definition: PmaVtxCandidate.cxx:24
pma::Dist2
double Dist2(const TVector2 &v1, const TVector2 &v2)
Definition: Utilities.cxx:19
pma::VtxCandidate::Test
double Test(const VtxCandidate &other) const
Definition: PmaVtxCandidate.cxx:245
d
Float_t d
Definition: plot.C:237
pma::VtxCandidate::fAssigned
std::vector< std::pair< pma::TrkCandidate, size_t > > fAssigned
Definition: PmaVtxCandidate.h:81
pma::VtxCandidate::IsAttached
bool IsAttached(pma::Track3D *trk) const
Definition: PmaVtxCandidate.cxx:38
pma::VtxCandidate::Size
size_t Size(void) const
Definition: PmaVtxCandidate.h:54
double pma::VtxCandidate::Mse ( void  ) const
inline

Definition at line 64 of file PmaVtxCandidate.h.

References fMse.

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

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

Definition at line 65 of file PmaVtxCandidate.h.

References fMse2D.

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

Definition at line 54 of file PmaVtxCandidate.h.

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

Referenced by MergeWith().

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

Definition at line 79 of file PmaVtxCandidate.cxx.

References fAssigned, and n.

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

Definition at line 245 of file PmaVtxCandidate.cxx.

References fCenter, and fErr.

Referenced by MergeWith(), and VtxCandidate().

246 {
247  double dx = fCenter[0] - other.fCenter[0];
248  double dy = fCenter[1] - other.fCenter[1];
249  double dz = fCenter[2] - other.fCenter[2];
250  double dw = fErr[0] * other.fErr[0] * dx * dx
251  + fErr[1] * other.fErr[1] * dy * dy
252  + fErr[2] * other.fErr[2] * dz * dz;
253  return sqrt(dw);
254 }
std::pair< pma::Track3D*, size_t > pma::VtxCandidate::Track ( size_t  i) const
inline

Definition at line 67 of file PmaVtxCandidate.h.

References fAssigned.

68  {
69  return std::pair< pma::Track3D*, size_t >(fAssigned[i].first.Track(), fAssigned[i].second);
70  }
pma::VtxCandidate::fAssigned
std::vector< std::pair< pma::TrkCandidate, size_t > > fAssigned
Definition: PmaVtxCandidate.h:81

Member Data Documentation

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

Definition at line 81 of file PmaVtxCandidate.h.

Referenced by Add(), Compute(), ComputeMse2D(), Has(), HasLoops(), IsAttached(), JoinTracks(), MaxAngle(), MergeWith(), Size(), and Track().

TVector3 pma::VtxCandidate::fCenter
private

Definition at line 82 of file PmaVtxCandidate.h.

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

TVector3 pma::VtxCandidate::fErr
private

Definition at line 82 of file PmaVtxCandidate.h.

Referenced by Compute(), and Test().

double pma::VtxCandidate::fMse
private

Definition at line 80 of file PmaVtxCandidate.h.

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

double pma::VtxCandidate::fMse2D
private

Definition at line 80 of file PmaVtxCandidate.h.

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

double pma::VtxCandidate::fSegMinLength
private

Definition at line 80 of file PmaVtxCandidate.h.

Referenced by Add(), and Compute().

const double pma::VtxCandidate::kMaxDistToTrack = 4.0
static

Definition at line 25 of file PmaVtxCandidate.h.

Referenced by Add().

const double pma::VtxCandidate::kMinDistToNode = 2.0
static

Definition at line 26 of file PmaVtxCandidate.h.

Referenced by JoinTracks().

bool pma::VtxCandidate::tracksJoined
private

Definition at line 79 of file PmaVtxCandidate.h.

Referenced by JoinTracks().

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