Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
MvtxStandaloneTracking.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file MvtxStandaloneTracking.cc
2 
5 #include <trackbase/TrkrDefs.h>
6 
8 #include <phool/PHIODataNode.h>
9 #include <phool/PHNodeIterator.h>
10 #include <phool/getClass.h>
11 
12 #include <TDecompSVD.h>
13 #include <TMath.h>
14 
15 #include <iostream>
16 #include <utility>
17 #include <stdio.h>
18 #include <map>
19 
21  : window_x_(10)
22  , window_z_(10)
23  , ghostrejection_(false)
24  , verbosity_(0)
25 {
26 
27 }
28 
30 {
31 }
32 
33 void
34 MvtxStandaloneTracking::RunTracking(PHCompositeNode* topNode, MvtxTrackList &trklst, std::vector<int> &lyrs)
35 {
36  // check for appropriate layers
37  if ( lyrs.size() < 3 || lyrs.size() > 4 )
38  {
39  std::cout << PHWHERE << "ERROR: Inappropriate number of input layers - " << lyrs.size() << std::endl;
40  return;
41  }
42 
43  trklst.clear();
44 
45  //------
46  //--- get cluster container
47  //------
48  clusters_ = findNode::getClass<TrkrClusterContainer>(topNode, "TRKR_CLUSTER");
49  if (!clusters_)
50  {
51  std::cout << PHWHERE << "ERROR: Can't find node TRKR_CLUSTER" << std::endl;
52  return;
53  }
54 
55 
56  //------
57  //--- associate clusters
58  //------
59  AssociateClusters(trklst, lyrs);
60 
61  if ( verbosity_ > 0 )
62  std::cout << PHWHERE << " Finished associating clusters. N candidates:" << trklst.size() << std::endl;
63  //------
64  //--- fit tracks
65  //------
66  for ( unsigned int itrk = 0; itrk < trklst.size(); itrk++)
67  {
68  TrackFitXY(trklst.at(itrk));
69  TrackFitZY(trklst.at(itrk));
70  } // itrk
71 
72  if ( verbosity_ > 1 )
73  PrintTrackCandidates(trklst);
74 
75 
76  //------
77  // --- choose best tracks
78  //------
79  if ( ghostrejection_ && trklst.size() > 1 )
80  RunGhostRejection(trklst);
81 
82  // --- done
83  return;
84 }
85 
86 void
87 MvtxStandaloneTracking::AssociateClusters(MvtxTrackList &trklst, std::vector<int> &lyrs)
88 {
89 
90  // --- utility class
91  //TrkrDefUtil util;
92 
93 
94  // --- loop over all clusters in the first desired layer
97  for ( TrkrClusterContainer::ConstIterator iter0 = clusrange0.first;
98  iter0 != clusrange0.second;
99  ++iter0)
100  {
101 
102  // -- loop over all clusters in the second desired layer
105  for ( TrkrClusterContainer::ConstIterator iter1 = clusrange1.first;
106  iter1 != clusrange1.second;
107  ++iter1)
108  {
109 
110  // get clusters
111  TrkrCluster* clus0 = iter0->second;
112  TrkrCluster* clus1 = iter1->second;
113 
114  // calculate slope & interecept in xy plane
115  double mxy = CalcSlope(clus0->getY(), clus0->getX(), clus1->getY(), clus1->getX());
116  double bxy = CalcIntecept(clus0->getY(), clus0->getX(), mxy);
117 
118  // calculate slope & interecept in zy plane
119  double mzy = CalcSlope(clus0->getY(), clus0->getZ(), clus1->getY(), clus1->getZ());
120  double bzy = CalcIntecept(clus0->getY(), clus0->getZ(), mzy);
121 
122  // -- loop over all clusters in the third desired layer
125  for ( TrkrClusterContainer::ConstIterator iter2 = clusrange2.first;
126  iter2 != clusrange2.second;
127  ++iter2)
128  {
129 
130  // check that the projection is within our window
131  if ( fabs((iter2->second)->getX() - CalcProjection((iter2->second)->getY(), mxy, bxy)) < window_x_ &&
132  fabs((iter2->second)->getZ() - CalcProjection((iter2->second)->getY(), mzy, bzy)) < window_z_ )
133  {
134 
135  // make the candidate
136  MvtxTrack trk;
137  (trk.ClusterList).push_back(iter0->second);
138  (trk.ClusterList).push_back(iter1->second);
139  (trk.ClusterList).push_back(iter2->second);
140 
141  // if there's another layer, require it
142  if ( lyrs.size() > 3 )
143  {
144  bool found_clu3 = false;
147  for ( TrkrClusterContainer::ConstIterator iter3 = clusrange3.first;
148  iter3 != clusrange3.second;
149  ++iter3)
150  {
151 
152  // check that the projection is within our window
153  if ( fabs((iter3->second)->getX() - CalcProjection((iter3->second)->getY(), mxy, bxy)) < window_x_ &&
154  fabs((iter3->second)->getZ() - CalcProjection((iter3->second)->getY(), mzy, bzy)) < window_z_ )
155  {
156  found_clu3 = true;
157  MvtxTrack tmp_trk = trk;
158  (tmp_trk.ClusterList).push_back(iter3->second);
159  trklst.push_back(tmp_trk);
160  }
161  }
162  if ( !found_clu3 ) trklst.push_back(trk);
163  }
164  else
165  {
166  // else we're done
167  trklst.push_back(trk);
168  }
169  }
170  } // clusrange 2
171  } // clusrange 1
172  } // clusrange 0
173 
174 }
175 
176 
177 void
178 MvtxStandaloneTracking::RunGhostRejection(MvtxTrackList &trklst)
179 {
180  if ( verbosity_ > 0 )
181  std::cout << PHWHERE << " Running Ghost Rejection on "
182  << trklst.size() << " tracks" << std::endl;
183 
184  // --- First, make a map of all cluster keys & track index
185  std::multimap<TrkrDefs::cluskey, unsigned int> key_trk_map;
186 
187  for (unsigned int itrk = 0; itrk < trklst.size(); itrk++)
188  {
189  for (unsigned int iclus = 0; iclus < trklst.at(itrk).ClusterList.size(); iclus++)
190  {
191  TrkrDefs::cluskey ckey = trklst.at(itrk).ClusterList.at(iclus)->getClusKey();
192  key_trk_map.insert(std::make_pair(ckey, itrk));
193  } // iclus
194  } // itrk
195 
196  // --- find clusters associated with more than one track and pick the best track
197  std::set<unsigned int> remove_set;
198  for ( auto iter = key_trk_map.begin(); iter != key_trk_map.end(); ++iter)
199  {
200  int ntrk = key_trk_map.count(iter->first);
201  if ( ntrk > 1 ) //cluster belong to more than one track
202  {
203  if ( verbosity_ > 1 )
204  {
205  std::cout << PHWHERE << " Tracks sharing cluster:" << ntrk << std::endl;
206  }
207 
208  // get the upper bound for this key
209  auto upiter = key_trk_map.upper_bound(iter->first);
210 
211  // iterate over common clusters and get the best track
212  double chi2_best = 9999.;
213  unsigned int idx_best = 0;
214  int ntrk = 0;
215  for ( auto jter = iter; jter != upiter; ++jter)
216  {
217  ntrk++;
218  double chi2 = trklst.at(jter->second).chi2_xy + trklst.at(jter->second).chi2_zy;
219  if ( chi2 < chi2_best && chi2 > 0 )
220  {
221  chi2_best = chi2;
222  idx_best = jter->second;
223  }
224  }
225 
226  for ( auto jter = iter; jter != upiter; ++jter)
227  {
228  if ( verbosity_ > 1 )
229  {
230  double chi2 = trklst.at(jter->second).chi2_xy + trklst.at(jter->second).chi2_zy;
231  std::cout << " "
232  << " trk idx: " << jter->second
233  << " chi2:" << chi2
234  << " m_xy:" << trklst.at(jter->second).m_xy;
235 
236  if ( jter->second == idx_best)
237  {
238  std::cout << " <--- BEST" << std::endl;
239  }
240  else
241  {
242  std::cout << std::endl;
243  }
244  }
245  // remove all pairs that aren't the best
246  if ( jter->second != idx_best )
247  {
248  remove_set.insert(jter->second);
249  }
250  }
251  }
252 
253  } // iter
254 
255 
256  // --- Now remove tracks from the track list
257  // reverse iterate so we don't have to keep track of indeces changed by removal
258  if ( verbosity_ > 0 )
259  {
260  std::cout << PHWHERE << " List of idx to remove:";
261  for ( auto it = remove_set.begin(); it != remove_set.end(); ++it)
262  {
263  std::cout << " " << *it;
264  }
265  std::cout << std::endl;
266  }
267 
268 
269  for ( auto rit = remove_set.rbegin(); rit != remove_set.rend(); ++rit)
270  trklst.erase(trklst.begin() + *rit);
271 
272  // --- done
273 
274  if ( verbosity_ > 1 )
275  {
276  std::cout << PHWHERE << " Remaining tracks:" << std::endl;
277  for ( unsigned int i = 0; i < trklst.size(); i++)
278  {
279  double chi2 = trklst.at(i).chi2_xy + trklst.at(i).chi2_zy;
280  std::cout << " chi2:" << chi2 << " m_xy:" << trklst.at(i).m_xy << std::endl;
281  }
282  }
283 
284  return;
285 }
286 
287 
288 TVectorD
289 MvtxStandaloneTracking::SolveGLS(TMatrixD &X, TVectorD &y, TMatrixD &L)
290 {
291  // Simple generalized linear least-squares fitter.
292  // Solve y(X) = beta'*X + error(L) for beta (vector of regression coefs.)
293  // L is the inverse covariance matrix for y.
294  // Least-squares solution involves solving X' L X beta = X' L y
295 
296  TMatrixD XT(X); XT.T();
297  TMatrixD A = XT * L * X;
298  TVectorD b = XT * L * y;
299 
300  // Now solve A*beta = b using SVD. Decompose A = U S V'.
301  TDecompSVD svd(A);
302  TMatrixD UT = svd.GetU(); UT.T();
303 
304  // Construct Moore-Penrose pseudoinverse of S
305  TVectorD s = svd.GetSig();
306  TMatrixD Sd(s.GetNrows(), s.GetNrows());
307  for (int i = 0; i < s.GetNrows(); i++)
308  Sd(i, i) = s(i) > 0 ? 1. / s(i) : 0.;
309 
310  TVectorD beta = svd.GetV() * Sd * UT * b;
311 
312  return beta;
313 }
314 
315 void
317 {
318  // Longitudinal/polar component
319  // Fit track using a straight line x(y') = x0 + c*y'.
320  // Assigns residuals, z-intercept, and polar angle.
321 
322  // m = # measurements; n = # parameters.
323  int m = (trk.ClusterList).size(), n = 2;
324 
325  TMatrixD X(m, n);
326  TMatrixD Cinv(m, m);
327  TVectorD y(m);
328  TVectorD x(m);
329 
330  for (int iclus = 0; iclus < m; iclus++)
331  {
332  TrkrCluster *clus = (trk.ClusterList).at(iclus);
333 
334  x(iclus) = clus->getX();
335  y(iclus) = clus->getY();
336 
337  X(iclus, 0) = 1;
338  X(iclus, 1) = y(iclus);
339 
340  Cinv(iclus, iclus) = 2.0 * sqrt(clus->getSize(0, 0));
341  }
342 
343  TVectorD beta = SolveGLS(X, x, Cinv);
344  double x0 = beta(0), c = beta(1);
345 
346  trk.m_xy = c;
347  trk.b_xy = x0;
348 
349  double chi2 = 0;
350  for (int iclus = 0; iclus < m; iclus++)
351  {
352  TrkrCluster *clus = (trk.ClusterList).at(iclus);
353 
354  double cx = clus->getX();
355  double cy = clus->getY();
356 
357  double px = cy * c + x0;
358 
359  double dx = cx - px;
360  trk.dx.push_back(dx);
361 
362  chi2 += pow(dx, 2) / clus->getError(0, 0);
363  }
364  chi2 /= double(m - 2);
365  trk.chi2_xy = chi2;
366 
367  return;
368 }
369 
370 void
372 {
373  // Longitudinal/polar component
374  // Fit track using a straight line z(y') = z0 + c*y'.
375  // Assigns residuals, z-intercept, and polar angle.
376 
377  // m = # measurements; n = # parameters.
378  int m = (trk.ClusterList).size(), n = 2;
379 
380  TMatrixD X(m, n);
381  TMatrixD Cinv(m, m);
382  TVectorD y(m);
383  TVectorD z(m);
384 
385  for (int iclus = 0; iclus < m; iclus++)
386  {
387  TrkrCluster *clus = (trk.ClusterList).at(iclus);
388 
389  z(iclus) = clus->getZ();
390  y(iclus) = clus->getY();
391 
392  X(iclus, 0) = 1;
393  X(iclus, 1) = y(iclus);
394 
395  Cinv(iclus, iclus) = 2.0 * sqrt(clus->getSize(2, 2));
396  }
397 
398  TVectorD beta = SolveGLS(X, z, Cinv);
399  double z0 = beta(0), c = beta(1);
400 
401  trk.m_zy = c;
402  trk.b_zy = z0;
403 
404  double chi2 = 0;
405  for (int iclus = 0; iclus < m; iclus++)
406  {
407  TrkrCluster *clus = (trk.ClusterList).at(iclus);
408 
409  double cz = clus->getZ();
410  double cy = clus->getY();
411 
412  double pz = cy * c + z0;
413 
414  double dz = cz - pz;
415  trk.dz.push_back(dz);
416 
417  chi2 += pow(dz, 2) / clus->getError(2, 2);
418  }
419  chi2 /= double(m - 2);
420  trk.chi2_zy = chi2;
421 
422  return;
423 }
424 
425 void
427 {
428  std::cout << "===================================================" << std::endl;
429  std::cout << "== " << PHWHERE << " Found " << trklst.size() << " Track Candidates" << std::endl;
430  std::cout << "===================================================" << std::endl;
431  for ( unsigned int i = 0; i < trklst.size(); i++)
432  {
433  std::cout << "== " << i << std::endl;
434  for ( unsigned int j = 0; j < trklst.at(i).ClusterList.size(); j++)
435  {
436  std::cout << " clus " << j
437  << " key:0x" << std::hex << trklst.at(i).ClusterList.at(j)->getClusKey() << std::dec
438  << " (" << trklst.at(i).ClusterList.at(j)->getX()
439  << ", " << trklst.at(i).ClusterList.at(j)->getY()
440  << ", " << trklst.at(i).ClusterList.at(j)->getZ()
441  << ")"
442  << " dx:" << trklst.at(i).dx.at(j)
443  << " dz:" << trklst.at(i).dz.at(j)
444  << std::endl;
445  }
446  std::cout << " xy"
447  << " m:" << trklst.at(i).m_xy
448  << " b:" << trklst.at(i).b_xy
449  << " chi2:" << trklst.at(i).chi2_xy
450  << std::endl;
451  std::cout << " zy"
452  << " m:" << trklst.at(i).m_zy
453  << " b:" << trklst.at(i).b_zy
454  << " chi2:" << trklst.at(i).chi2_zy
455  << std::endl;
456  } // i
457  std::cout << "===================================================" << std::endl;
458 }
459 
460 double
461 MvtxStandaloneTracking::CalcSlope(double x0, double y0, double x1, double y1)
462 {
463  return (y1 - y0) / (x1 - x0);
464 }
465 
466 double
467 MvtxStandaloneTracking::CalcIntecept(double x0, double y0, double m)
468 {
469  return y0 - x0 * m;
470 }
471 
472 double
473 MvtxStandaloneTracking::CalcProjection(double x, double m, double b)
474 {
475  return m * x + b;
476 }
477