ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
PHTpcTrackSeedCircleFit.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file PHTpcTrackSeedCircleFit.cc
2 
3 #include "AssocInfoContainer.h"
4 
5 #include <trackbase/TrkrDefs.h> // for cluskey, getTrkrId, tpcId
11 
13 
14 #include <g4main/PHG4Hit.h> // for PHG4Hit
15 #include <g4main/PHG4Particle.h> // for PHG4Particle
16 #include <g4main/PHG4HitDefs.h> // for keytype
17 
19 
20 #include <phool/getClass.h>
21 #include <phool/phool.h>
22 
23 #if __cplusplus < 201402L
24 #include <boost/make_unique.hpp>
25 #endif
26 
27 #include <TF1.h>
28 
29 #include <climits> // for UINT_MAX
30 #include <iostream> // for operator<<, basic_ostream
31 #include <cmath> // for fabs, std::sqrt
32 #include <set> // for _Rb_tree_const_iterator
33 #include <utility> // for pair
34 #include <memory>
35 
36 namespace
37 {
38 
39  // square
40  template<class T> inline constexpr T square( const T& x ) { return x*x; }
41 
42  // radius
43  template<class T> T get_r( const T& x, const T& y ) { return std::sqrt( square(x) + square(y) ); }
44 
45  void line_fit(std::vector<std::pair<double,double>> points, double &a, double &b)
46  {
47  // copied from: https://www.bragitoff.com
48  // we want to fit z vs radius
49 
50  //variables for sums/sigma of xi,yi,xi^2,xiyi etc
51  double xsum=0,x2sum=0,ysum=0,xysum=0;
52  for( const auto& [r,z]:points )
53  {
54  xsum+=r; //calculate sigma(xi)
55  ysum+=z; //calculate sigma(yi)
56  x2sum+=square(r); //calculate sigma(x^2i)
57  xysum+=r*z; //calculate sigma(xi*yi)
58  }
59  a=(points.size()*xysum-xsum*ysum)/(points.size()*x2sum-square(xsum)); //calculate slope
60  b=(x2sum*ysum-xsum*xysum)/(points.size()*x2sum-square(xsum)); //calculate intercept
61 
62  return;
63  }
64 
65  void line_fit_clusters(std::vector<Acts::Vector3D>& globPos, double &a, double &b)
66  {
67  std::vector<std::pair<double,double>> points;
68  std::transform( globPos.begin(), globPos.end(), std::back_inserter( points ), []( const Acts::Vector3D& pos )
69  { return std::make_pair( get_r( pos(0), pos(1) ), pos(2)); } );
70 
71  line_fit(points, a, b);
72  }
73 
74  void CircleFitByTaubin (std::vector<Acts::Vector3D> points, double &R, double &X0, double &Y0)
75  /*
76  Circle fit to a given set of data points (in 2D)
77  This is an algebraic fit, due to Taubin, based on the journal article
78  G. Taubin, "Estimation Of Planar Curves, Surfaces And Nonplanar
79  Space Curves Defined By Implicit Equations, With
80  Applications To Edge And Range Image Segmentation",
81  IEEE Trans. PAMI, Vol. 13, pages 1115-1138, (1991)
82  */
83  {
84  int iter,IterMAX=99;
85 
86  double Mz,Mxy,Mxx,Myy,Mxz,Myz,Mzz,Cov_xy,Var_z;
87  double A0,A1,A2,A22,A3,A33;
88  double x,y;
89  double DET,Xcenter,Ycenter;
90 
91  // Compute x- and y- sample means
92  double meanX = 0;
93  double meanY = 0;
94  double weight = 0;
95  for(const auto& position : points)
96  {
97  meanX += position(0);
98  meanY += position(1);
99  weight++;
100  }
101  meanX /= weight;
102  meanY /= weight;
103 
104  // computing moments
105 
106  Mxx=Myy=Mxy=Mxz=Myz=Mzz=0.;
107 
108  for (const auto& position : points)
109  {
110  double Xi = position(0) - meanX; // centered x-coordinates
111  double Yi = position(1) - meanY; // centered y-coordinates
112  double Zi = Xi*Xi + Yi*Yi;
113 
114  Mxy += Xi*Yi;
115  Mxx += Xi*Xi;
116  Myy += Yi*Yi;
117  Mxz += Xi*Zi;
118  Myz += Yi*Zi;
119  Mzz += Zi*Zi;
120  }
121  Mxx /= weight;
122  Myy /= weight;
123  Mxy /= weight;
124  Mxz /= weight;
125  Myz /= weight;
126  Mzz /= weight;
127 
128  // computing coefficients of the characteristic polynomial
129 
130  Mz = Mxx + Myy;
131  Cov_xy = Mxx*Myy - Mxy*Mxy;
132  Var_z = Mzz - Mz*Mz;
133  A3 = 4*Mz;
134  A2 = -3*Mz*Mz - Mzz;
135  A1 = Var_z*Mz + 4*Cov_xy*Mz - Mxz*Mxz - Myz*Myz;
136  A0 = Mxz*(Mxz*Myy - Myz*Mxy) + Myz*(Myz*Mxx - Mxz*Mxy) - Var_z*Cov_xy;
137  A22 = A2 + A2;
138  A33 = A3 + A3 + A3;
139 
140  // finding the root of the characteristic polynomial
141  // using Newton's method starting at x=0
142  // (it is guaranteed to converge to the right root)
143 
144  for (x=0.,y=A0,iter=0; iter<IterMAX; ++iter) // usually, 4-6 iterations are enough
145  {
146  const double Dy = A1 + x*(A22 + A33*x);
147  const double xnew = x - y/Dy;
148  if ((xnew == x)||(!isfinite(xnew))) break;
149  const double ynew = A0 + xnew*(A1 + xnew*(A2 + xnew*A3));
150  if (std::abs(ynew)>=std::abs(y)) break;
151  x = xnew; y = ynew;
152  }
153 
154  // computing parameters of the fitting circle
155 
156  DET = x*x - x*Mz + Cov_xy;
157  Xcenter = (Mxz*(Myy - x) - Myz*Mxy)/DET/2;
158  Ycenter = (Myz*(Mxx - x) - Mxz*Mxy)/DET/2;
159 
160  // assembling the output
161 
162  X0 = Xcenter + meanX;
163  Y0 = Ycenter + meanY;
164  R = std::sqrt(square(Xcenter) + square(Ycenter) + Mz);
165  }
166 
167 // std::vector<double> GetCircleClusterResiduals(std::vector<std::pair<double,double>> points, double R, double X0, double Y0)
168 // {
169 // std::vector<double> residues;
170 // std::transform( points.begin(), points.end(), std::back_inserter( residues ),
171 // [R,X0,Y0]( const std::pair<double,double>& point )
172 // {
173 // // The shortest distance of a point from a circle is along the radial; line from the circle center to the point
174 // const auto& x = point.first;
175 // const auto& y = point.second;
176 // return std::sqrt( square(x-X0) + square(y-Y0))-R;
177 // } );
178 // return residues;
179 // }
180 
181 // std::vector<double> GetLineClusterResiduals(std::vector<std::pair<double,double>> points, double A, double B)
182 // {
183 // std::vector<double> residues;
184 // std::transform( points.begin(), points.end(), std::back_inserter( residues ),
185 // [A,B]( const std::pair<double,double>& point )
186 // {
187 // // The shortest distance of a point from a circle is along the radial; line from the circle center to the point
188 // const auto& r = point.first;
189 // const auto& z = point.second;
190 // const double a = -A;
191 // const double b = 1.0;
192 // const double c = -B;
193 // return std::abs(a*r+b*z+c)/std::sqrt(square(a)+square(b));
194 // });
195 // return residues;
196 // }
197 
198  void findRoot(
199  const double R, const double X0,
200  const double Y0, double& x,
201  double& y)
202  {
216  const double miny = (std::sqrt(square(X0)*square(R)*square(Y0) + square(R)
217  * pow(Y0,4)) + square(X0)*Y0 + pow(Y0, 3))
218  / (square(X0)+square(Y0));
219 
220  const double miny2 = (-std::sqrt(square(X0)*square(R)*square(Y0) + square(R)
221  * pow(Y0,4)) + square(X0) * Y0 + pow(Y0, 3))
222  / (square(X0)+square(Y0));
223 
224  const double minx = std::sqrt(square(R) - square(miny-Y0)) + X0;
225  const double minx2 = -std::sqrt(square(R) - square(miny2-Y0)) + X0;
226 
228  if(std::abs(minx) < std::abs(minx2))
229  x = minx;
230  else
231  x = minx2;
232 
233  if(std::abs(miny) < std::abs(miny2))
234  y = miny;
235  else
236  y = miny2;
237  }
238 
239 }
240 
241 //____________________________________________________________________________..
243  SubsysReco(name)
244 {}
245 
246 //____________________________________________________________________________..
248 {
249  // get relevant nodes
250  int ret = GetNodes( topnode );
251  if( ret != Fun4AllReturnCodes::EVENT_OK ) return ret;
252 
254 
255 //____________________________________________________________________________..
257 {
258 
259  // _track_map contains the TPC seed track stubs
260  // We want to associate these TPC track seeds with a collision vertex
261  // All we need is to project the TPC clusters in Z to the beam line.
262  // The TPC track seeds are given a position that is the PCA of the line and circle fit to the beam line
263 
264  if(Verbosity() > 0)
265  std::cout << PHWHERE << " TPC track map size " << _track_map->size() << std::endl;
266 
267 
268  // loop over the TPC track seeds
269  for (auto phtrk_iter = _track_map->begin();
270  phtrk_iter != _track_map->end();
271  ++phtrk_iter)
272  {
273  const auto& track_key = phtrk_iter->first;
274  SvtxTrack* tracklet_tpc = phtrk_iter->second;
275 
276  if (Verbosity() > 1)
277  {
278  std::cout
279  << __LINE__
280  << ": Processing seed itrack: " << phtrk_iter->first
281  << ": nhits: " << tracklet_tpc-> size_cluster_keys()
282  << ": pT: " << tracklet_tpc->get_pt()
283  << ": phi: " << tracklet_tpc->get_phi()
284  << ": eta: " << tracklet_tpc->get_eta()
285  << std::endl;
286  }
287 
288  // get the tpc track seed cluster positions in z and r
289 
290  // Get the TPC clusters for this tracklet
291  std::vector<TrkrCluster*> clusters = getTrackClusters(tracklet_tpc);
292  if(Verbosity() > 3) std::cout << " TPC tracklet " << tracklet_tpc->get_id() << " clusters.size " << clusters.size() << std::endl;
293 
294  // count TPC layers for this track
295  bool reject_track = false;
296  std::set<unsigned int> layers;
297  for (unsigned int i=0; i<clusters.size(); ++i)
298  {
299  if(!clusters[i])
300  {
301  if(Verbosity() > 0) std::cout << " trackid " << phtrk_iter->first << " no cluster found, skip track" << std::endl;
302  reject_track = true;
303  break;
304  }
305  unsigned int layer = TrkrDefs::getLayer(clusters[i]->getClusKey());
306  layers.insert(layer);
307  }
308 
309  if(reject_track) continue;
310 
311  unsigned int nlayers = layers.size();
312  if(Verbosity() > 2) std::cout << " TPC layers this track: " << nlayers << std::endl;
313 
314 
315  std::vector<Acts::Vector3D> globalClusterPositions;
316  for (unsigned int i=0; i<clusters.size(); ++i)
317  {
318  const Acts::Vector3D global = getGlobalPosition(clusters.at(i));
319  globalClusterPositions.push_back(global);
320 
321  if(Verbosity() > 3)
322  {
323  ActsTransformations transformer;
324  auto global_before = transformer.getGlobalPosition(clusters.at(i),
325  _surfmaps,
326  _tGeometry);
327  TrkrDefs::cluskey key = clusters.at(i)->getClusKey();
328  std::cout << "CircleFit: Cluster: " << key << " _corrected_clusters " << _are_clusters_corrected << std::endl;
329  std::cout << " Global before: " << global_before[0] << " " << global_before[1] << " " << global_before[2] << std::endl;
330  std::cout << " Global after : " << global[0] << " " << global[1] << " " << global[2] << std::endl;
331  }
332  }
333 
334  if(clusters.size() < 3)
335  {
336  if(Verbosity() > 3) std::cout << PHWHERE << " -- skip this tpc tracklet, not enough TPC clusters " << std::endl;
337  continue; // skip to the next TPC tracklet
338  }
339 
340  // get the straight line representing the z trajectory in the form of z vs radius
341  double A = 0;
342  double B = 0;
343  line_fit_clusters(globalClusterPositions, A, B);
344  if(Verbosity() > 2)
345  { std::cout << "PHTpcTrackSeedCircleFit::process_event - track: " << track_key << " A=" << A << " B=" << B << std::endl; }
346 
347  // Project this TPC tracklet to the beam line and store the projections
348  const double z_proj = B;
349 
350  // set the track Z position to the Z dca
351  tracklet_tpc->set_z(z_proj);
352 
353  // Finished association of track with vertex, now we modify the track parameters
354 
355  // extract the track theta
356  double track_angle = atan(A); // referenced to 90 degrees
357 
358 
359  // make circle fit
360  double R, X0, Y0;
361  CircleFitByTaubin(globalClusterPositions, R, X0, Y0);
362  if(Verbosity() > 2)
363  { std::cout << "PHTpcTrackSeedCircleFit::process_event - track: " << track_key << " R=" << R << " X0=" << X0 << " Y0=" << Y0 << std::endl; }
364 
365  // set the track x and y positions to the circle PCA
366  double dcax, dcay;
367  findRoot(R, X0, Y0, dcax, dcay);
368  if(std::isnan(dcax) or std::isnan(dcay))
369  { continue; }
370 
371  tracklet_tpc->set_x(dcax);
372  tracklet_tpc->set_y(dcay);
373 
374  double pt_track = tracklet_tpc->get_pt();
375  if(Verbosity() > 5)
376  {
377  // get the pT from radius of circle as a check
378  static constexpr double Bfield = 1.4; // T
379  double pt_circle = 0.3 * Bfield * R * 0.01; // convert cm to m
380  std::cout << "PHTpcTrackSeedCircleFit::process_event-"
381  << " R: " << R
382  << " Bfield: " << Bfield
383  << " pt_circle: " << pt_circle << " pt_track: " << pt_track
384  << std::endl;
385  }
386 
387  // We want the angle of the tangent relative to the positive x axis
388  // start with the angle of the radial line from vertex to circle center
389 
390  double dx = X0 - dcax;
391  double dy = Y0 - dcay;
392  double phi= atan2(-dx,dy);
393 
394  // convert to the angle of the tangent to the circle
395  // we need to know if the track proceeds clockwise or CCW around the circle
396  double dx0 = globalClusterPositions.at(0)(0) - X0;
397  double dy0 = globalClusterPositions.at(0)(1) - Y0;
398  double phi0 = atan2(dy0, dx0);
399  double dx1 = globalClusterPositions.at(1)(0) - X0;
400  double dy1 = globalClusterPositions.at(1)(1) - Y0;
401  double phi1 = atan2(dy1, dx1);
402  double dphi = phi1 - phi0;
403 
404  // need to deal with the switch from -pi to +pi at phi = 180 degrees
405  // final phi - initial phi must be < 180 degrees for it to be a valid track
406  if(dphi > M_PI) dphi -= 2.0 * M_PI;
407  if(dphi < - M_PI) dphi += M_PI;
408 
409  if(Verbosity() > 5)
410  {
411  std::cout << " charge " << tracklet_tpc->get_charge() << " phi0 " << phi0*180.0 / M_PI << " phi1 " << phi1*180.0 / M_PI
412  << " dphi " << dphi*180.0 / M_PI << std::endl;
413  }
414 
415  // whether we add 180 degrees depends on the angle of the bend
416  if(dphi < 0)
417  {
418  phi += M_PI;
419  if(phi > M_PI)
420  { phi -= 2. * M_PI; }
421  }
422 
423  if(Verbosity() > 5)
424  std::cout << " input track phi " << tracklet_tpc->get_phi() << " new phi " << phi << std::endl;
425 
426  // get the updated values of px, py, pz from the pT and the angles found here
427  double px_new = pt_track * cos(phi);
428  double py_new = pt_track * sin(phi);
429  double ptrack_new = pt_track / cos(track_angle);
430  double pz_new = ptrack_new * sin(track_angle);
431 
432  if(Verbosity() > 4)
433  std::cout << " input track id " << tracklet_tpc->get_id()
434  << " input track mom " << tracklet_tpc->get_p() << " new mom " << ptrack_new
435  << " px in " << tracklet_tpc->get_px() << " px " << px_new
436  << " py in " << tracklet_tpc->get_py() << " py " << py_new
437  << " pz in " << tracklet_tpc->get_pz() << " pz " << pz_new
438  << " eta in " << tracklet_tpc->get_eta() << " phi in " << tracklet_tpc->get_phi() * 180.0 / M_PI
439  << " track angle " << track_angle * 180.0 / M_PI
440  << std::endl;
441 
442  // update track on node tree
443  tracklet_tpc->set_px(px_new);
444  tracklet_tpc->set_py(py_new);
445  tracklet_tpc->set_pz(pz_new);
446 
447  if(Verbosity() > 5)
448  {
449  std::cout << " new mom " << tracklet_tpc->get_p() << " new eta " << tracklet_tpc->get_eta()
450  << " new phi " << tracklet_tpc->get_phi() * 180.0 / M_PI << std::endl;
451  }
452  } // end loop over TPC track seeds
453 
454  if(Verbosity() > 0)
455  std::cout << " Final track map size " << _track_map->size() << std::endl;
456 
457  if (Verbosity() > 0)
458  std::cout << "PHTpcTrackSeedCircleFit::process_event(PHCompositeNode *topNode) Leaving process_event" << std::endl;
459 
461 }
462 
465 
467 {
468  _surfmaps = findNode::getClass<ActsSurfaceMaps>(topNode,"ActsSurfaceMaps");
469  if(!_surfmaps)
470  {
471  std::cout << PHWHERE << "Error, can't find acts surface maps" << std::endl;
473  }
474 
475  _tGeometry = findNode::getClass<ActsTrackingGeometry>(topNode,"ActsTrackingGeometry");
476  if(!_tGeometry)
477  {
478  std::cout << PHWHERE << "Error, can't find acts tracking geometry" << std::endl;
480  }
481 
483  _cluster_map = findNode::getClass<TrkrClusterContainer>(topNode, "TRKR_CLUSTER_TRUTH");
484  else
485  {
486  _cluster_map = findNode::getClass<TrkrClusterContainer>(topNode, "CORRECTED_TRKR_CLUSTER");
487  if(_cluster_map)
488  {
489  std::cout << " using CORRECTED_TRKR_CLUSTER node "<< std::endl;
490  }
491  else
492  {
493  std::cout << " CORRECTED_TRKR_CLUSTER node not found, using TRKR_CLUSTER " << std::endl;
494  _cluster_map = findNode::getClass<TrkrClusterContainer>(topNode, "TRKR_CLUSTER");
495  _are_clusters_corrected = false;
496  }
497  }
498  if (!_cluster_map)
499  {
500  std::cerr << PHWHERE << " ERROR: Can't find node TRKR_CLUSTER" << std::endl;
502  }
503 
504  // tpc distortion correction
505  _dcc = findNode::getClass<TpcDistortionCorrectionContainer>(topNode,"TpcDistortionCorrectionContainer");
506  if( _dcc )
507  { std::cout << "PHTpcTrackSeedCircleFit::get_Nodes - found TPC distortion correction container" << std::endl; }
508 
509  _track_map = findNode::getClass<SvtxTrackMap>(topNode, _track_map_name);
510  if (!_track_map)
511  {
512  std::cerr << PHWHERE << " ERROR: Can't find SvtxTrackMap" << std::endl;
514  }
515 
517 }
518 
519 std::vector<TrkrCluster*> PHTpcTrackSeedCircleFit::getTrackClusters(SvtxTrack *tracklet_tpc)
520 {
521  std::vector<TrkrCluster*> clusters;
522 
523  for(SvtxTrack::ConstClusterKeyIter key_iter = tracklet_tpc->begin_cluster_keys();
524  key_iter != tracklet_tpc->end_cluster_keys();
525  ++key_iter)
526  {
527  const TrkrDefs::cluskey& cluster_key = *key_iter;
528 
529  // only consider TPC clusters
530  if( TrkrDefs::getTrkrId(cluster_key) != TrkrDefs::tpcId ) continue;
531 
532  // get the cluster
533  TrkrCluster *tpc_clus = _cluster_map->findCluster(cluster_key);
534 
535  // insert in list
536  clusters.push_back(tpc_clus);
537 
538  if(Verbosity() > 5)
539  {
540  unsigned int layer = TrkrDefs::getLayer(cluster_key);
541  std::cout << " TPC cluster in layer " << layer << " with local position " << tpc_clus->getLocalX()
542  << " " << tpc_clus->getLocalY() << " clusters.size() " << clusters.size() << std::endl;
543  }
544  }
545  return clusters;
546  }
547 
549 {
550  // get global position from Acts transform
551  ActsTransformations transformer;
552  auto globalpos = transformer.getGlobalPosition(cluster,
553  _surfmaps,
554  _tGeometry);
555 
556  // check if TPC distortion correction are in place and apply if clusters are not from the corrected node
558  if(_dcc) { globalpos = _distortionCorrection.get_corrected_position( globalpos, _dcc ); }
559 
560  return globalpos;
561 }
562