Classes | Public Member Functions | Private Member Functions | Private Attributes | List of all members
pma::PMAlgVertexing Class Reference

#include <PMAlgVertexing.h>

Classes

struct  Config
 

Public Member Functions

 PMAlgVertexing (const Config &config)
 
 PMAlgVertexing (const fhicl::ParameterSet &pset)
 
 ~PMAlgVertexing ()
 
void reset ()
 
size_t run (const detinfo::DetectorPropertiesData &detProp, pma::TrkCandidateColl &trk_input)
 
size_t run (pma::TrkCandidateColl &trk_input, const std::vector< TVector3 > &vtx_input)
 
std::vector< std::pair< TVector3, std::vector< std::pair< size_t, bool > > > > getVertices (const pma::TrkCandidateColl &tracks, bool onlyBranching=false) const
 
std::vector< std::pair< TVector3, size_t > > getKinks (const pma::TrkCandidateColl &tracks) const
 

Private Member Functions

bool has (const std::vector< size_t > &v, size_t idx) const
 
std::vector< pma::VtxCandidatefirstPassCandidates () const
 
std::vector< pma::VtxCandidatesecondPassCandidates () const
 
size_t makeVertices (detinfo::DetectorPropertiesData const &detProp, std::vector< pma::VtxCandidate > &candidates)
 
std::vector< std::pair< double, double > > getdQdx (const pma::Track3D &trk) const
 Get dQ/dx sequence to detect various features. More...
 
double convolute (size_t idx, size_t len, double *adc, double const *shape) const
 Get convolution value. More...
 
bool isSingleParticle (pma::Track3D *trk1, pma::Track3D *trk2) const
 Check if colinear in 3D and dQ/dx with no significant step. More...
 
void mergeBrokenTracks (pma::TrkCandidateColl &trk_input) const
 
void splitMergedTracks (pma::TrkCandidateColl &trk_input) const
 Split track and add vertex and reoptimize when dQ/dx step detected. More...
 
void findKinksOnTracks (const detinfo::DetectorPropertiesData &detProp, pma::TrkCandidateColl &trk_input) const
 Remove penalty on the angle if kink detected and reopt track. More...
 
void cleanTracks ()
 
void sortTracks (const pma::TrkCandidateColl &trk_input)
 
void collectTracks (pma::TrkCandidateColl &result)
 

Private Attributes

pma::TrkCandidateColl fOutTracks
 
pma::TrkCandidateColl fShortTracks
 
pma::TrkCandidateColl fEmTracks
 
pma::TrkCandidateColl fExcludedTracks
 
double fMinTrackLength
 
bool fFindKinks
 
double fKinkMinDeg
 
double fKinkMinStd
 

Detailed Description

Definition at line 32 of file PMAlgVertexing.h.

Constructor & Destructor Documentation

pma::PMAlgVertexing::PMAlgVertexing ( const Config config)

Definition at line 15 of file PMAlgVertexing.cxx.

16 {
18 
22 }
fhicl::Atom< double > MinTrackLength
static Config * config
Definition: config.cpp:1054
fhicl::Atom< bool > FindKinks
fhicl::Atom< double > KinkMinStd
fhicl::Atom< double > KinkMinDeg
pma::PMAlgVertexing::PMAlgVertexing ( const fhicl::ParameterSet pset)
inline

Definition at line 54 of file PMAlgVertexing.h.

55  {}
PMAlgVertexing(const Config &config)
pma::PMAlgVertexing::~PMAlgVertexing ( )

Definition at line 25 of file PMAlgVertexing.cxx.

26 {
27  cleanTracks();
28 }

Member Function Documentation

void pma::PMAlgVertexing::cleanTracks ( )
private

Definition at line 32 of file PMAlgVertexing.cxx.

33 {
34  for (auto& t : fOutTracks.tracks())
35  t.DeleteTrack();
36  fOutTracks.clear();
37 
38  for (auto& t : fShortTracks.tracks())
39  t.DeleteTrack();
41 
42  for (auto& t : fEmTracks.tracks())
43  t.DeleteTrack();
44  fEmTracks.clear();
45 
46  for (auto& t : fExcludedTracks.tracks())
47  t.DeleteTrack();
49 }
pma::TrkCandidateColl fEmTracks
pma::TrkCandidateColl fShortTracks
pma::TrkCandidateColl fOutTracks
pma::TrkCandidateColl fExcludedTracks
std::vector< TrkCandidate > const & tracks() const
void pma::PMAlgVertexing::collectTracks ( pma::TrkCandidateColl result)
private

Definition at line 53 of file PMAlgVertexing.cxx.

54 {
55  mf::LogVerbatim("pma::PMAlgVertexing") << "clean input: " << result.size() << std::endl;
56  for (auto& t : result.tracks())
57  t.DeleteTrack();
58  result.clear();
59 
60  mf::LogVerbatim("pma::PMAlgVertexing")
61  << "fill input from out: " << fOutTracks.size() << std::endl;
62  for (auto const& t : fOutTracks.tracks())
63  result.push_back(t);
64  fOutTracks.clear();
65 
66  mf::LogVerbatim("pma::PMAlgVertexing")
67  << "fill input from short: " << fShortTracks.size() << std::endl;
68  for (auto const& t : fShortTracks.tracks())
69  result.push_back(t);
71 
72  mf::LogVerbatim("pma::PMAlgVertexing") << "fill input from em: " << fEmTracks.size() << std::endl;
73  for (auto const& t : fEmTracks.tracks())
74  result.push_back(t);
75  fEmTracks.clear();
76 
77  mf::LogVerbatim("pma::PMAlgVertexing")
78  << "copy back excluded tracks: " << fExcludedTracks.size() << std::endl;
79  for (auto const& t : fExcludedTracks.tracks())
80  result.push_back(t);
82 }
pma::TrkCandidateColl fEmTracks
MaybeLogger_< ELseverityLevel::ELsev_info, true > LogVerbatim
size_t size() const
pma::TrkCandidateColl fShortTracks
pma::TrkCandidateColl fOutTracks
pma::TrkCandidateColl fExcludedTracks
void push_back(const TrkCandidate &trk)
std::vector< TrkCandidate > const & tracks() const
QTextStream & endl(QTextStream &s)
double pma::PMAlgVertexing::convolute ( size_t  idx,
size_t  len,
double *  adc,
double const *  shape 
) const
private

Get convolution value.

Definition at line 390 of file PMAlgVertexing.cxx.

391 {
392  size_t half = len >> 1;
393  double v, mean = 0.0, stdev = 0.0;
394  for (size_t i = 0; i < len; i++) {
395  v = adc[idx - half + i];
396  mean += v;
397  stdev += v * v;
398  }
399  mean /= len;
400  stdev /= len;
401  stdev -= mean;
402 
403  double sum = 0.0;
404  for (size_t i = 0; i < len; i++)
405  sum += (adc[idx - half + i] - mean) * shape[i];
406 
407  return sum / sqrt(stdev);
408 }
def stdev(lst)
Definition: HandyFuncs.py:269
int16_t adc
Definition: CRTFragment.hh:202
double mean(sqlite3 *db, std::string const &table_name, std::string const &column_name)
Definition: statistics.cc:16
void pma::PMAlgVertexing::findKinksOnTracks ( const detinfo::DetectorPropertiesData detProp,
pma::TrkCandidateColl trk_input 
) const
private

Remove penalty on the angle if kink detected and reopt track.

Definition at line 542 of file PMAlgVertexing.cxx.

544 {
545  if (trk_input.size() < 1) return;
546 
547  mf::LogVerbatim("pma::PMAlgVertexing")
548  << "Find kinks on tracks, reopt with no penalty on angle where kinks.";
549  for (size_t t = 0; t < trk_input.size(); ++t) {
550  pma::Track3D* trk = trk_input[t].Track();
551  if (trk->Nodes().size() < 5) continue;
552 
553  trk->Optimize(detProp, 0, 1.0e-5, false);
554 
555  std::vector<size_t> tested_nodes;
556  bool kinkFound = true;
557  while (kinkFound) {
558  int kinkIdx = -1, nnodes = 0;
559  double mean = 0.0, stdev = 0.0, min = 1.0, max_a = 0.0;
560  for (size_t n = 1; n < trk->Nodes().size() - 1; ++n) {
561  auto const& node = *(trk->Nodes()[n]);
562 
563  if (node.IsVertex() || node.IsTPCEdge()) continue;
564  nnodes++;
565 
566  double c = -node.SegmentCosTransverse();
567  double a = 180.0 * (1 - std::acos(node.SegmentCosTransverse()) / TMath::Pi());
568  mean += c;
569  stdev += c * c;
570  if ((c < min) && !has(tested_nodes, n)) {
571  if ((n > 1) && (n < trk->Nodes().size() - 2)) kinkIdx = n;
572  min = c;
573  max_a = a;
574  }
575  }
576 
577  kinkFound = false;
578  if ((nnodes > 2) && (kinkIdx > 0) && (max_a > fKinkMinDeg)) {
579  mean /= nnodes;
580  stdev /= nnodes;
581  stdev -= mean * mean;
582 
583  double thr = 1.0 - fKinkMinStd * stdev;
584  if (min < thr) {
585  mf::LogVerbatim("pma::PMAlgVertexing") << " kink a: " << max_a << "deg";
586  trk->Nodes()[kinkIdx]->SetVertex(true);
587  tested_nodes.push_back(kinkIdx);
588  kinkFound = true;
589 
590  trk->Optimize(detProp, 0, 1.0e-5, false, false);
591  double c = -trk->Nodes()[kinkIdx]->SegmentCosTransverse();
592  double a =
593  180.0 * (1 - std::acos(trk->Nodes()[kinkIdx]->SegmentCosTransverse()) / TMath::Pi());
594 
595  if ((a <= fKinkMinDeg) || (c >= thr)) // not a significant kink after precise optimization
596  {
597  mf::LogVerbatim("pma::PMAlgVertexing") << " -> tag removed after reopt";
598  trk->Nodes()[kinkIdx]->SetVertex(
599  false); // kinkIdx is saved in tested_nodes to avoid inf loop
600  }
601  }
602  }
603  }
604 
605  mf::LogVerbatim("pma::PMAlgVertexing") << "-------- done --------";
606  }
607 }
MaybeLogger_< ELseverityLevel::ELsev_info, true > LogVerbatim
size_t size() const
def stdev(lst)
Definition: HandyFuncs.py:269
G4double thr[100]
Definition: G4S2Light.cc:59
double Optimize(const detinfo::DetectorPropertiesData &detProp, int nNodes=-1, double eps=0.01, bool selAllHits=true, bool setAllNodes=true, size_t selSegHits=0, size_t selVtxHits=0)
Main optimization method.
const double e
std::void_t< T > n
const double a
T min(sqlite3 *const db, std::string const &table_name, std::string const &column_name)
Definition: statistics.h:55
std::vector< pma::Node3D * > const & Nodes() const noexcept
Definition: PmaTrack3D.h:335
double mean(sqlite3 *db, std::string const &table_name, std::string const &column_name)
Definition: statistics.cc:16
bool has(const std::vector< size_t > &v, size_t idx) const
std::vector< pma::VtxCandidate > pma::PMAlgVertexing::firstPassCandidates ( ) const
private

Definition at line 127 of file PMAlgVertexing.cxx.

128 {
129  std::vector<pma::VtxCandidate> candidates;
130  for (size_t t = 0; t < fOutTracks.size() - 1; t++) {
131  for (size_t u = t + 1; u < fOutTracks.size(); u++) {
132  pma::VtxCandidate candidate;
133  if (!candidate.Add(fOutTracks[t])) break; // no segments with length > thr
134 
135  // **************************** try Mse2D / or only Mse ************************************
136  if (candidate.Add(fOutTracks[u]) && (sqrt(candidate.Mse()) < 1.0))
137  //if (candidate.Add(fOutTracks[u]) && (sqrt(candidate.Mse()) < 2.0) && (candidate.Mse2D() < 1.0))
138  {
139  candidates.push_back(candidate);
140  }
141  }
142  }
143  return candidates;
144 }
size_t size() const
double Mse() const
pma::TrkCandidateColl fOutTracks
bool Add(const pma::TrkCandidate &trk)
std::vector< std::pair< double, double > > pma::PMAlgVertexing::getdQdx ( const pma::Track3D trk) const
private

Get dQ/dx sequence to detect various features.

Definition at line 349 of file PMAlgVertexing.cxx.

350 {
351  std::vector<std::pair<double, double>> result;
352 
353  unsigned int view = geo::kZ;
354  unsigned int nhits = trk.NHits(view);
355  unsigned int max = nhits;
356 
357  nhits = trk.NHits(geo::kV);
358  if (nhits > max) {
359  max = nhits;
360  view = geo::kV;
361  }
362 
363  nhits = trk.NHits(geo::kU);
364  if (nhits > max) {
365  max = nhits;
366  view = geo::kU;
367  }
368 
369  if (max >= 16) {
370  std::map<size_t, std::vector<double>> dqdx;
371  trk.GetRawdEdxSequence(dqdx, view);
372 
373  for (size_t i = 0; i < trk.size(); i++) {
374  auto it = dqdx.find(i);
375  if (it != dqdx.end()) {
376  if (it->second[6] > 0.0) // dx > 0
377  {
378  double dvalue = it->second[5] / it->second[6];
379  result.emplace_back(std::pair<double, double>(dvalue, it->second[7]));
380  }
381  }
382  }
383  }
384 
385  return result;
386 }
static QCString result
Planes which measure V.
Definition: geo_types.h:130
Planes which measure Z direction.
Definition: geo_types.h:132
Planes which measure U.
Definition: geo_types.h:129
unsigned int NHits(unsigned int view) const
Definition: PmaTrack3D.cxx:422
static int max(int a, int b)
double GetRawdEdxSequence(std::map< size_t, std::vector< double >> &dedx, unsigned int view=geo::kZ, unsigned int skip=0, bool inclDisabled=false) const
size_t size() const
Definition: PmaTrack3D.h:108
std::vector< std::pair< TVector3, size_t > > pma::PMAlgVertexing::getKinks ( const pma::TrkCandidateColl tracks) const

Definition at line 652 of file PMAlgVertexing.cxx.

653 {
654  std::vector<std::pair<TVector3, size_t>> ksel;
655  for (size_t t = 0; t < tracks.size(); ++t) {
656  pma::Track3D const* trk = tracks[t].Track();
657  for (size_t n = 1; n < trk->Nodes().size() - 1; ++n) {
658  pma::Node3D const* node = trk->Nodes()[n];
659  if (!node->IsBranching() && node->IsVertex()) {
660  ksel.emplace_back(std::pair<TVector3, size_t>(node->Point3D(), t));
661  }
662  }
663  }
664  return ksel;
665 }
TVector3 const & Point3D() const
Definition: PmaNode3D.h:45
size_t size() const
std::void_t< T > n
bool IsVertex() const
Check fIsVertex flag.
Definition: PmaNode3D.h:74
bool IsBranching() const
Belongs to more than one track?
Definition: PmaNode3D.cxx:436
std::vector< pma::Node3D * > const & Nodes() const noexcept
Definition: PmaTrack3D.h:335
std::vector< std::pair< TVector3, std::vector< std::pair< size_t, bool > > > > pma::PMAlgVertexing::getVertices ( const pma::TrkCandidateColl tracks,
bool  onlyBranching = false 
) const

Definition at line 611 of file PMAlgVertexing.cxx.

612 {
613  std::vector<std::pair<TVector3, std::vector<std::pair<size_t, bool>>>> vsel;
614  std::vector<pma::Node3D const*> bnodes;
615 
616  for (size_t t = 0; t < tracks.size(); ++t) {
617  pma::Track3D const* trk = tracks[t].Track();
618  pma::Node3D const* firstNode = trk->Nodes().front();
619  if (!(onlyBranching || firstNode->IsBranching())) {
620  std::vector<std::pair<size_t, bool>> tidx;
621  tidx.emplace_back(std::pair<size_t, bool>(t, true));
622  vsel.emplace_back(
623  std::pair<TVector3, std::vector<std::pair<size_t, bool>>>(trk->front()->Point3D(), tidx));
624  }
625 
626  bool pri = true;
627  for (auto node : trk->Nodes())
628  if (node->IsBranching()) {
629  bool found = false;
630  for (size_t n = 0; n < bnodes.size(); n++)
631  if (node == bnodes[n]) {
632  vsel[n].second.emplace_back(std::pair<size_t, bool>(t, pri));
633  found = true;
634  break;
635  }
636  if (!found) {
637  std::vector<std::pair<size_t, bool>> tidx;
638  tidx.emplace_back(std::pair<size_t, bool>(t, pri));
639  vsel.emplace_back(
640  std::pair<TVector3, std::vector<std::pair<size_t, bool>>>(node->Point3D(), tidx));
641  bnodes.push_back(node);
642  }
643  pri = false;
644  }
645  }
646 
647  return vsel;
648 }
size_t size() const
pma::Hit3D const * front() const
Definition: PmaTrack3D.h:98
struct vector vector
TVector3 const & Point3D() const
Definition: PmaHit3D.h:55
std::void_t< T > n
bool IsBranching() const
Belongs to more than one track?
Definition: PmaNode3D.cxx:436
std::vector< pma::Node3D * > const & Nodes() const noexcept
Definition: PmaTrack3D.h:335
bool pma::PMAlgVertexing::has ( const std::vector< size_t > &  v,
size_t  idx 
) const
inlineprivate

Definition at line 86 of file PMAlgVertexing.h.

87  {
88  for (auto c : v)
89  if (c == idx) return true;
90  return false;
91  }
bool pma::PMAlgVertexing::isSingleParticle ( pma::Track3D trk1,
pma::Track3D trk2 
) const
private

Check if colinear in 3D and dQ/dx with no significant step.

Definition at line 411 of file PMAlgVertexing.cxx.

412 {
413  const double minCos = 0.996194698; // 5 deg (is it ok?)
414  double segCos =
415  trk1->Segments().back()->GetDirection3D().Dot(trk2->Segments().front()->GetDirection3D());
416  if (segCos < minCos) {
417  mf::LogVerbatim("pma::PMAlgVertexing") << " has large angle, cos: " << segCos;
418  return false;
419  }
420 
421  const size_t stepShapeLen = 16;
422  const size_t stepShapeHalf = stepShapeLen >> 1;
423  const double stepShape[stepShapeLen] = {
424  -1., -1., -1., -1., -1., -1., -1., -1., 1., 1., 1., 1., 1., 1., 1., 1.};
425 
426  auto dqdx1 = getdQdx(*trk1);
427  if (dqdx1.size() < stepShapeHalf) return false;
428  auto dqdx2 = getdQdx(*trk2);
429  if (dqdx2.size() < stepShapeHalf) return false;
430 
431  const size_t adcLen =
432  stepShapeLen + 2; // 1 sample before/after to check convolution at 3 points in total
433  const size_t adcHalf = adcLen >> 1;
434 
435  double dqdx[adcLen];
436  for (size_t i = 0; i < adcLen; i++)
437  dqdx[i] = 0.0;
438 
439  bool has_m = true;
440  for (int i = adcHalf - 1, j = dqdx1.size() - 1; i >= 0; i--, j--) {
441  if (j >= 0)
442  dqdx[i] = dqdx1[j].first;
443  else {
444  dqdx[i] = dqdx[i + 1];
445  has_m = false;
446  }
447  }
448  bool has_p = true;
449  for (size_t i = adcHalf, j = 0; i < adcLen; i++, j++) {
450  if (j < dqdx2.size())
451  dqdx[i] = dqdx2[j].first;
452  else {
453  dqdx[i] = dqdx[i - 1];
454  has_p = false;
455  }
456  }
457 
458  double sum_m = 0.0;
459  if (has_m) sum_m = convolute(adcHalf - 1, stepShapeLen, dqdx, stepShape);
460  double sum_0 = convolute(adcHalf, stepShapeLen, dqdx, stepShape);
461  double sum_p = 0.0;
462  if (has_p) sum_p = convolute(adcHalf + 1, stepShapeLen, dqdx, stepShape);
463 
464  const double convMin = 0.8;
465  if ((fabs(sum_m) >= convMin) || (fabs(sum_0) >= convMin) || (fabs(sum_p) >= convMin)) {
466  mf::LogVerbatim("pma::PMAlgVertexing")
467  << " has step in conv.values: " << sum_m << ", " << sum_0 << ", " << sum_p;
468  return false;
469  }
470  else {
471  mf::LogVerbatim("pma::PMAlgVertexing")
472  << " single particle, conv.values: " << sum_m << ", " << sum_0 << ", " << sum_p;
473  return true;
474  }
475 }
MaybeLogger_< ELseverityLevel::ELsev_info, true > LogVerbatim
std::vector< pma::Segment3D * > const & Segments() const noexcept
Definition: PmaTrack3D.h:326
std::vector< std::pair< double, double > > getdQdx(const pma::Track3D &trk) const
Get dQ/dx sequence to detect various features.
double convolute(size_t idx, size_t len, double *adc, double const *shape) const
Get convolution value.
size_t pma::PMAlgVertexing::makeVertices ( detinfo::DetectorPropertiesData const &  detProp,
std::vector< pma::VtxCandidate > &  candidates 
)
private

Definition at line 167 of file PMAlgVertexing.cxx.

169 {
170  bool merged = true;
171  while (merged && (candidates.size() > 1)) {
172  size_t k_best, l_best, k = 0;
173  while (k < candidates.size() - 1) {
174  size_t l = k + 1;
175  while (l < candidates.size()) {
176  if (candidates[l].Has(candidates[k])) {
177  candidates[k] = candidates[l];
178  candidates.erase(candidates.begin() + l);
179  }
180  else if (candidates[k].Has(candidates[l])) {
181  candidates.erase(candidates.begin() + l);
182  }
183  else
184  l++;
185  }
186  k++;
187  }
188 
189  merged = false;
190  double d_thr = 1.0; // 1.0 = max weighted dist. threshold
191  double d, dmin = d_thr;
192 
193  k = 0;
194  while (k < candidates.size() - 1) {
195  size_t l = k + 1;
196  while (l < candidates.size()) {
197  d = candidates[k].Test(candidates[l]);
198  if (d < dmin) {
199  dmin = d;
200  k_best = k;
201  l_best = l;
202  }
203  l++;
204  }
205  k++;
206  }
207  if ((dmin < d_thr) && candidates[k_best].MergeWith(candidates[l_best])) {
208  candidates.erase(candidates.begin() + l_best);
209  merged = true;
210  }
211  }
212 
213  mf::LogVerbatim("pma::PMAlgVertexing") << "*** Vtx candidates: " << candidates.size();
214  std::vector<pma::VtxCandidate> toJoin;
215  bool select = true;
216  while (select) {
217  int s, nmax = 0, c_best = -1;
218  double a, amax = 0.0;
219 
220  for (size_t v = 0; v < candidates.size(); v++) {
221  if (candidates[v].HasLoops()) continue;
222 
223  bool maybeCorrelated = false;
224  for (size_t u = 0; u < toJoin.size(); u++)
225  if (toJoin[u].IsAttached(candidates[v]) || // connected with tracks or close centers
226  (pma::Dist2(toJoin[u].Center(), candidates[v].Center()) < 15.0 * 15.0)) {
227  maybeCorrelated = true;
228  break;
229  }
230  if (maybeCorrelated) continue;
231 
232  s = (int)candidates[v].Size(2 * fMinTrackLength);
233  a = candidates[v].MaxAngle(1.0);
234 
235  if ((s > nmax) || ((s == nmax) && (a > amax))) {
236  nmax = s;
237  amax = a;
238  c_best = v;
239  }
240  /*
241  mf::LogVerbatim("pma::PMAlgVertexing")
242  << "center x:" << candidates[v].Center().X()
243  << " y:" << candidates[v].Center().Y()
244  << " z:" << candidates[v].Center().Z();
245 
246  for (size_t i = 0; i < candidates[v].Size(); i++)
247  mf::LogVerbatim("pma::PMAlgVertexing")
248  << " trk:" << i << " "
249  << candidates[v].Track(i).first->size();
250 
251  mf::LogVerbatim("pma::PMAlgVertexing")
252  << " dist 3D:" << sqrt(candidates[v].Mse())
253  << " 2D:" << sqrt(candidates[v].Mse2D())
254  << " max ang:" << a;
255 */
256  }
257  if (c_best >= 0) {
258  toJoin.push_back(candidates[c_best]);
259  candidates.erase(candidates.begin() + c_best);
260  }
261  else
262  select = false;
263  }
264  mf::LogVerbatim("pma::PMAlgVertexing") << "*** Vtx selected to join: " << toJoin.size();
265 
266  size_t njoined = 0;
267  for (auto& c : toJoin) {
268  if (c.JoinTracks(detProp, fOutTracks, fEmTracks)) njoined++;
269  }
270 
271  return njoined;
272 }
pma::TrkCandidateColl fEmTracks
MaybeLogger_< ELseverityLevel::ELsev_info, true > LogVerbatim
double Dist2(const TVector2 &v1, const TVector2 &v2)
Definition: Utilities.cxx:37
static QStrList * l
Definition: config.cpp:1044
const double a
pma::TrkCandidateColl fOutTracks
auto select(T const &...t)
Definition: select.h:146
static QCString * s
Definition: config.cpp:1042
void pma::PMAlgVertexing::mergeBrokenTracks ( pma::TrkCandidateColl trk_input) const
private

Find elastic scattering vertices on tracks, merge back tracks that were split during vertex finding. 3D angle between two tracks and dQ/dx is checked.

Definition at line 478 of file PMAlgVertexing.cxx.

479 {
480  if (trk_input.size() < 2) return;
481 
482  mf::LogVerbatim("pma::PMAlgVertexing") << "Find and merge tracks broken by vertices.";
483  bool merged = true;
484  while (merged) {
485  merged = false;
486  for (size_t t = 0; t < trk_input.size(); t++) {
487  pma::Track3D* trk1 = trk_input[t].Track();
488  pma::Track3D* trk2 = 0;
489 
490  pma::Node3D* node = trk1->Nodes().front();
491  if (node->Prev()) {
492  pma::Segment3D* seg = static_cast<pma::Segment3D*>(node->Prev());
493  trk2 = seg->Parent();
494  if ((trk1 != trk2) && isSingleParticle(trk2, trk1)) // note: reverse order
495  {
496  //merged = true;
497  break;
498  }
499  }
500 
501  trk2 = 0;
502  double c, maxc = 0.0;
503  pma::Vector3D dir1 = trk1->Segments().back()->GetDirection3D();
504  node = trk1->Nodes().back();
505  for (size_t n = 0; n < node->NextCount(); n++) {
506  pma::Segment3D* seg = static_cast<pma::Segment3D*>(node->Next(n));
507  pma::Track3D* tst = seg->Parent();
508  if (tst != trk1) // should always be true: the last node of trk1 is tested
509  {
510  c = dir1.Dot(tst->Segments().front()->GetDirection3D());
511  if (c > maxc) {
512  maxc = c;
513  trk2 = tst;
514  }
515  }
516  }
517  if ((trk2) && isSingleParticle(trk1, trk2)) {
518  //merged = true;
519  break;
520  }
521  }
522  }
523  mf::LogVerbatim("pma::PMAlgVertexing") << "-------- done --------";
524 }
MaybeLogger_< ELseverityLevel::ELsev_info, true > LogVerbatim
size_t size() const
virtual unsigned int NextCount(void) const
Definition: SortedObjects.h:84
recob::tracking::Vector_t Vector3D
Definition: Utilities.h:31
std::void_t< T > n
std::vector< pma::Segment3D * > const & Segments() const noexcept
Definition: PmaTrack3D.h:326
virtual pma::SortedObjectBase * Next(unsigned int index=0) const
Definition: SortedObjects.h:79
bool isSingleParticle(pma::Track3D *trk1, pma::Track3D *trk2) const
Check if colinear in 3D and dQ/dx with no significant step.
std::vector< pma::Node3D * > const & Nodes() const noexcept
Definition: PmaTrack3D.h:335
virtual pma::SortedObjectBase * Prev(void) const
Definition: SortedObjects.h:44
pma::Track3D * Parent(void) const
Definition: PmaSegment3D.h:68
void pma::PMAlgVertexing::reset ( )
inline

Definition at line 60 of file PMAlgVertexing.h.

61  {
62  cleanTracks();
63  }
size_t pma::PMAlgVertexing::run ( const detinfo::DetectorPropertiesData detProp,
pma::TrkCandidateColl trk_input 
)

Copy input tracks, find 3D vertices, connect tracks, break them or flip if needed, reoptimize track structures. Result is returned as a collection of new tracks, that replaces content of trk_input (old tracks are deleted). Vertices can be accessed with getVertices function.

Definition at line 276 of file PMAlgVertexing.cxx.

278 {
279  if (fFindKinks) findKinksOnTracks(detProp, trk_input);
280 
281  if (trk_input.size() < 2) {
282  mf::LogWarning("pma::PMAlgVertexing") << "need min two source tracks!";
283  return 0;
284  }
285 
286  sortTracks(trk_input); // copy input and split by tag/size
287 
288  size_t nvtx = 0;
289  mf::LogVerbatim("pma::PMAlgVertexing") << "Pass #1:";
290  //std::cout << "Pass #1:" << std::endl;
291  if (fOutTracks.size() > 1) {
292  size_t nfound = 0;
293  do {
294  auto candidates = firstPassCandidates();
295  if (candidates.size()) {
296  nfound = makeVertices(detProp, candidates);
297  nvtx += nfound;
298  }
299  else
300  nfound = 0;
301  } while (nfound > 0);
302  mf::LogVerbatim("pma::PMAlgVertexing") << " " << nvtx << " vertices.";
303  //std::cout << " " << nvtx << " vertices." << std::endl;
304  }
305  else
306  mf::LogVerbatim("pma::PMAlgVertexing") << " ...short tracks only.";
307 
308  mf::LogVerbatim("pma::PMAlgVertexing") << "Pass #2:";
309  //std::cout << "Pass #2:" << std::endl;
310  if (fOutTracks.size() && fEmTracks.size()) {
311  size_t nfound = 1; // just to start
312  while (nfound && fEmTracks.size()) {
313  auto candidates = secondPassCandidates();
314  if (candidates.size()) {
315  nfound = makeVertices(detProp, candidates);
316  nvtx += nfound;
317  }
318  else
319  nfound = 0;
320  }
321  mf::LogVerbatim("pma::PMAlgVertexing") << " " << nvtx << " vertices.";
322  //std::cout << " " << nvtx << " vertices." << std::endl;
323  }
324  else
325  mf::LogVerbatim("pma::PMAlgVertexing") << " ...no tracks.";
326 
327  collectTracks(trk_input);
328 
329  mergeBrokenTracks(trk_input);
330 
331  return nvtx;
332 }
pma::TrkCandidateColl fEmTracks
MaybeLogger_< ELseverityLevel::ELsev_info, true > LogVerbatim
size_t size() const
std::vector< pma::VtxCandidate > firstPassCandidates() const
std::vector< pma::VtxCandidate > secondPassCandidates() const
void collectTracks(pma::TrkCandidateColl &result)
pma::TrkCandidateColl fOutTracks
size_t makeVertices(detinfo::DetectorPropertiesData const &detProp, std::vector< pma::VtxCandidate > &candidates)
MaybeLogger_< ELseverityLevel::ELsev_warning, false > LogWarning
void sortTracks(const pma::TrkCandidateColl &trk_input)
void findKinksOnTracks(const detinfo::DetectorPropertiesData &detProp, pma::TrkCandidateColl &trk_input) const
Remove penalty on the angle if kink detected and reopt track.
void mergeBrokenTracks(pma::TrkCandidateColl &trk_input) const
size_t pma::PMAlgVertexing::run ( pma::TrkCandidateColl trk_input,
const std::vector< TVector3 > &  vtx_input 
)

Copy input tracks, use provided 3D vertices to connect tracks, break tracks or flip if needed, reoptimize track structures. Result is returned as a collection of new tracks, that replaces content of trk_input (old tracks are deleted). Input vertices that were actually associated to tracks are copied to the output collection (use getVertices function).

Definition at line 336 of file PMAlgVertexing.cxx.

337 {
338  sortTracks(trk_input); // copy input and split by tag/size
339 
340  // ....
341 
342  //collectTracks(trk_input); // return output in place of (deleted) input
343 
344  return 0;
345 }
void sortTracks(const pma::TrkCandidateColl &trk_input)
std::vector< pma::VtxCandidate > pma::PMAlgVertexing::secondPassCandidates ( ) const
private

Definition at line 147 of file PMAlgVertexing.cxx.

148 {
149  std::vector<pma::VtxCandidate> candidates;
150  for (size_t t = 0; t < fOutTracks.size(); t++)
151  if (fOutTracks[t].Track()->Length() > fMinTrackLength) {
152  for (size_t u = 0; u < fEmTracks.size(); u++) {
153  pma::VtxCandidate candidate;
154  if (!candidate.Add(fOutTracks[t])) break; // no segments with length > thr
155 
156  if (fOutTracks[t].Track() == fEmTracks[u].Track()) continue;
157 
158  if (candidate.Add(fEmTracks[u]) && (sqrt(candidate.Mse()) < 1.0)) {
159  candidates.push_back(candidate);
160  }
161  }
162  }
163  return candidates;
164 }
pma::TrkCandidateColl fEmTracks
size_t size() const
double Mse() const
pma::TrkCandidateColl fOutTracks
bool Add(const pma::TrkCandidate &trk)
TrackCollectionProxyElement< TrackCollProxy > Track
Proxy to an element of a proxy collection of recob::Track objects.
Definition: Track.h:1036
void pma::PMAlgVertexing::sortTracks ( const pma::TrkCandidateColl trk_input)
private

Definition at line 86 of file PMAlgVertexing.cxx.

87 {
88  cleanTracks();
89 
90  for (auto const& t : trk_input.tracks()) {
91  pma::Track3D* copy = new pma::Track3D(*(t.Track()));
92  int key = t.Key();
93 
94  if ((t.Track()->GetT0() != 0) || // do not create vertices on cosmic muons, now track with any
95  // non-zero t0 is a muon,
96  (t.Track()->HasTagFlag(pma::Track3D::kCosmic))) // tag for track recognized as a cosmic ray
97  // is starting to be used
98  {
99  fExcludedTracks.tracks().emplace_back(copy, key);
100  continue;
101  }
102 
103  double l = t.Track()->Length();
104 
105  if (t.Track()->GetTag() == pma::Track3D::kEmLike) {
106  if (l > 2 * fMinTrackLength)
107  fOutTracks.tracks().emplace_back(copy, key);
108  else
109  fEmTracks.tracks().emplace_back(copy, key);
110  }
111  else {
112  if (l > fMinTrackLength)
113  fOutTracks.tracks().emplace_back(copy, key);
114  else
115  fEmTracks.tracks().emplace_back(copy, key);
116  }
117  }
118  mf::LogVerbatim("pma::PMAlgVertexing") << "long tracks: " << fOutTracks.size() << std::endl;
119  mf::LogVerbatim("pma::PMAlgVertexing")
120  << "em and short tracks: " << fEmTracks.size() << std::endl;
121  mf::LogVerbatim("pma::PMAlgVertexing")
122  << "excluded cosmic tracks: " << fExcludedTracks.size() << std::endl;
123 }
pma::TrkCandidateColl fEmTracks
MaybeLogger_< ELseverityLevel::ELsev_info, true > LogVerbatim
size_t size() const
static QStrList * l
Definition: config.cpp:1044
def key(type, name=None)
Definition: graph.py:13
pma::TrkCandidateColl fOutTracks
pma::TrkCandidateColl fExcludedTracks
T copy(T const &v)
std::vector< TrkCandidate > const & tracks() const
QTextStream & endl(QTextStream &s)
void pma::PMAlgVertexing::splitMergedTracks ( pma::TrkCandidateColl trk_input) const
private

Split track and add vertex and reoptimize when dQ/dx step detected.

Definition at line 528 of file PMAlgVertexing.cxx.

529 {
530  if (trk_input.size() < 1) return;
531 
532  mf::LogVerbatim("pma::PMAlgVertexing") << "Find missed vtx by dQ/dx steps along merged tracks.";
533  size_t t = 0;
534  while (t < trk_input.size()) {
535  t++;
536  }
537  mf::LogVerbatim("pma::PMAlgVertexing") << "-------- done --------";
538 }
MaybeLogger_< ELseverityLevel::ELsev_info, true > LogVerbatim
size_t size() const

Member Data Documentation

pma::TrkCandidateColl pma::PMAlgVertexing::fEmTracks
private

Definition at line 118 of file PMAlgVertexing.h.

pma::TrkCandidateColl pma::PMAlgVertexing::fExcludedTracks
private

Definition at line 119 of file PMAlgVertexing.h.

bool pma::PMAlgVertexing::fFindKinks
private

Definition at line 130 of file PMAlgVertexing.h.

double pma::PMAlgVertexing::fKinkMinDeg
private

Definition at line 132 of file PMAlgVertexing.h.

double pma::PMAlgVertexing::fKinkMinStd
private

Definition at line 133 of file PMAlgVertexing.h.

double pma::PMAlgVertexing::fMinTrackLength
private

Definition at line 127 of file PMAlgVertexing.h.

pma::TrkCandidateColl pma::PMAlgVertexing::fOutTracks
private

Definition at line 118 of file PMAlgVertexing.h.

pma::TrkCandidateColl pma::PMAlgVertexing::fShortTracks
private

Definition at line 118 of file PMAlgVertexing.h.


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