ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
PHMicromegasTpcTrackMatching.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file PHMicromegasTpcTrackMatching.cc
2 
3 #include "AssocInfoContainer.h"
4 //#include "PHTrackPropagating.h" // for PHTrackPropagating
5 
7 
10 
12 #include <trackbase/TrkrCluster.h> // for TrkrCluster
13 #include <trackbase/TrkrDefs.h> // for cluskey, getLayer, TrkrId
15 #include <trackbase/TrkrHitSet.h>
18 
19 #include <trackbase_historic/SvtxTrack.h> // for SvtxTrack, SvtxTrack::C...
22 
24 #include <fun4all/SubsysReco.h> // for SubsysReco
25 
26 #include <phool/getClass.h>
27 #include <phool/phool.h>
28 
29 #include <TF1.h>
30 #include <TVector3.h>
31 
32 #include <array>
33 #include <cmath> // for sqrt, std::abs, atan2, cos
34 #include <iostream> // for operator<<, basic_ostream
35 #include <map> // for map
36 #include <set> // for _Rb_tree_const_iterator
37 #include <utility> // for pair, make_pair
38 
39 namespace
40 {
41 
43  template<class T>
44  inline constexpr T square( const T& x ) { return x*x; }
45 
47  template<class T>
48  inline constexpr T get_r( const T& x, const T& y ) { return std::sqrt( square(x) + square(y) ); }
49 
62  void CircleFitByTaubin (std::vector<Acts::Vector3D>& clusters, double &R, double &X0, double &Y0)
63  {
64  int iter,IterMAX=99;
65 
66  double Mz,Mxy,Mxx,Myy,Mxz,Myz,Mzz,Cov_xy,Var_z;
67  double A0,A1,A2,A22,A3,A33;
68  double x,y;
69  double DET,Xcenter,Ycenter;
70 
71  // Compute x- and y- sample means
72  double meanX = 0;
73  double meanY = 0;
74  double weight = 0;
75 
76  for(auto& global : clusters)
77  {
78  meanX += global(0);
79  meanY += global(1);
80  weight++;
81  }
82  meanX /= weight;
83  meanY /= weight;
84 
85  // computing moments
86 
87  Mxx=Myy=Mxy=Mxz=Myz=Mzz=0.;
88 
89  for (auto& global : clusters)
90  {
91  double Xi = global(0) - meanX; // centered x-coordinates
92  double Yi = global(1) - meanY; // centered y-coordinates
93  double Zi = Xi*Xi + Yi*Yi;
94 
95  Mxy += Xi*Yi;
96  Mxx += Xi*Xi;
97  Myy += Yi*Yi;
98  Mxz += Xi*Zi;
99  Myz += Yi*Zi;
100  Mzz += Zi*Zi;
101  }
102  Mxx /= weight;
103  Myy /= weight;
104  Mxy /= weight;
105  Mxz /= weight;
106  Myz /= weight;
107  Mzz /= weight;
108 
109  // computing coefficients of the characteristic polynomial
110 
111  Mz = Mxx + Myy;
112  Cov_xy = Mxx*Myy - Mxy*Mxy;
113  Var_z = Mzz - Mz*Mz;
114  A3 = 4*Mz;
115  A2 = -3*Mz*Mz - Mzz;
116  A1 = Var_z*Mz + 4*Cov_xy*Mz - Mxz*Mxz - Myz*Myz;
117  A0 = Mxz*(Mxz*Myy - Myz*Mxy) + Myz*(Myz*Mxx - Mxz*Mxy) - Var_z*Cov_xy;
118  A22 = A2 + A2;
119  A33 = A3 + A3 + A3;
120 
121  // finding the root of the characteristic polynomial
122  // using Newton's method starting at x=0
123  // (it is guaranteed to converge to the right root)
124 
125  for (x=0.,y=A0,iter=0; iter<IterMAX; iter++) // usually, 4-6 iterations are enough
126  {
127  double Dy = A1 + x*(A22 + A33*x);
128  double xnew = x - y/Dy;
129  if ((xnew == x)||(!isfinite(xnew))) break;
130  double ynew = A0 + xnew*(A1 + xnew*(A2 + xnew*A3));
131  if (std::abs(ynew)>=std::abs(y)) break;
132  x = xnew; y = ynew;
133  }
134 
135  // computing parameters of the fitting circle
136 
137  DET = x*x - x*Mz + Cov_xy;
138  Xcenter = (Mxz*(Myy - x) - Myz*Mxy)/DET/2;
139  Ycenter = (Myz*(Mxx - x) - Mxz*Mxy)/DET/2;
140 
141  // assembling the output
142 
143  X0 = Xcenter + meanX;
144  Y0 = Ycenter + meanY;
145  R = sqrt(Xcenter*Xcenter + Ycenter*Ycenter + Mz);
146  }
147 
148  // 2D linear fit
149 void line_fit(std::vector<Acts::Vector3D>& clusters, double &a, double &b)
150  {
151  // copied from: https://www.bragitoff.com
152  // we want to fit z vs radius
153 
154  //variables for sums/sigma of xi,yi,xi^2,xiyi etc
155  double xsum=0,x2sum=0,ysum=0,xysum=0;
156 
157  for (auto& global : clusters)
158  {
159  const double z = global(2);
160  const double r = get_r( global(0), global(1) );
161 
162  xsum=xsum+r; //calculate sigma(xi)
163  ysum=ysum+z; //calculate sigma(yi)
164  x2sum=x2sum+ square(r); //calculate sigma(x^2i)
165  xysum=xysum+r*z; //calculate sigma(xi*yi)
166  }
167 
168  a=(clusters.size()*xysum-xsum*ysum)/(clusters.size()*x2sum-xsum*xsum); //calculate slope
169  b=(x2sum*ysum-xsum*xysum)/(x2sum*clusters.size()-xsum*xsum); //calculate intercept
170  return;
171  }
172 
173 
187  bool circle_circle_intersection(double r1, double r2, double x2, double y2, double &xplus, double &yplus, double &xminus, double &yminus)
188  {
189 
190  const double D = square(r1) - square(r2) + square(x2) + square(y2);
191  const double a = 1.0 + square(x2/y2);
192  const double b = - D * x2/square(y2);
193  const double c = square(D/(2.*y2)) - square(r1);
194  const double delta = square(b)-4.*a*c;
195  if( delta < 0 ) return false;
196 
197  const double sqdelta = std::sqrt( delta );
198 
199  xplus = (-b + sqdelta ) / (2. * a);
200  xminus = (-b - sqdelta ) / (2. * a);
201 
202  // both values of x are valid
203  // but for each of those values, there are two possible y values on circle 1
204  // but only one of those falls on the radical line:
205 
206  yplus = -(2*x2*xplus - D) / (2.*y2);
207  yminus = -(2*x2*xminus - D) / (2.*y2);
208  return true;
209 
210  }
211 
213 
220  bool circle_line_intersection(
221  double r, double xc, double yc,
222  double x0, double y0, double nx, double ny,
223  double &xplus, double &yplus, double &xminus, double &yminus)
224  {
225  if( ny == 0 )
226  {
227  // vertical lines are defined by ny=0 and x = x0
228  xplus = xminus = x0;
229 
230  // calculate y accordingly
231  const double delta = square(r) - square(x0-xc);
232  if( delta < 0 ) return false;
233 
234  const double sqdelta = std::sqrt( delta );
235  yplus = yc + sqdelta;
236  yminus = yc - sqdelta;
237 
238  } else {
239 
240  const double a = square(nx) + square(ny);
241  const double b = -2.*( square(ny)*xc + square(nx)*x0 + nx*ny*(y0-yc) );
242  const double c = square(ny)*(square(xc)-square(r)) + square(ny*(y0-yc)+nx*x0);
243  const double delta = square(b) - 4.*a*c;
244  if( delta < 0 ) return false;
245 
246  const double sqdelta = std::sqrt( delta );
247  xplus = (-b + sqdelta)/(2.*a);
248  xminus = (-b - sqdelta)/(2.*a);
249 
250  yplus = y0 -(nx/ny)*(xplus-x0);
251  yminus = y0 - (nx/ny)*(xminus-x0);
252 
253  }
254 
255  return true;
256 
257  }
258 
259  // streamer of TVector3
260  [[maybe_unused]] inline std::ostream& operator << (std::ostream& out, const TVector3& vector )
261  {
262  out << "( " << vector.x() << "," << vector.y() << "," << vector.z() << ")";
263  return out;
264  }
265 
266 }
267 
268 //____________________________________________________________________________..
270  SubsysReco(name)
271 {}
272 
273 //____________________________________________________________________________..
275 {
276 
277  std::cout << std::endl << PHWHERE
278  << " rphi_search_win inner layer " << _rphi_search_win[0]
279  << " z_search_win inner layer " << _z_search_win[0]
280  << " rphi_search_win outer layer " << _rphi_search_win[1]
281  << " z_search_win outer layer " << _z_search_win[1]
282  << std::endl;
283 
284  // load micromegas geometry
285  _geomContainerMicromegas = findNode::getClass<PHG4CylinderGeomContainer>(topNode, "CYLINDERGEOM_MICROMEGAS_FULL" );
287  {
288  std::cout << PHWHERE << "Could not find CYLINDERGEOM_MICROMEGAS_FULL." << std::endl;
290  }
291 
292  // ensures there are at least two micromegas layers
294  {
295  std::cout << PHWHERE << "Inconsistent number of micromegas layers." << std::endl;
297  }
298 
299  // get first micromegas layer
301 
302  fdrphi = new TF1("fdrphi", "[0] + [1]*std::abs(x)");
304  fdrphi->SetParameter(1, _par1);
305 
306  // base class setup
307 
308  int ret = GetNodes(topNode);
309  if (ret != Fun4AllReturnCodes::EVENT_OK) return ret;
310 
311  return ret;
312 }
313 
314 //____________________________________________________________________________..
316 {
317 
318  if(_n_iteration >0)
319  {
320  _iteration_map = findNode::getClass<TrkrClusterIterationMapv1>(topNode, "CLUSTER_ITERATION_MAP");
321  if (!_iteration_map)
322  {
323  std::cerr << PHWHERE << "Cluster Iteration Map missing, aborting." << std::endl;
325  }
326  }
327  // _track_map contains the TPC seed track stubs
328  // We will add the micromegas cluster to the TPC tracks already on the node tree
329  // We will have to expand the number of tracks whenever we find multiple matches to the silicon
330 
331  _event++;
332 
333  if(Verbosity() > 0)
334  std::cout << PHWHERE << " Event " << _event << " TPC track map size " << _track_map->size() << std::endl;
335 
336  // We remember the original size of the TPC track map here
337  const unsigned int original_track_map_lastkey = _track_map->empty() ? 0:std::prev(_track_map->end())->first;
338 
339  // loop over the original TPC tracks
340  for (auto phtrk_iter = _track_map->begin(); phtrk_iter != _track_map->end(); ++phtrk_iter)
341  {
342 
343  // we may add tracks to the map, so we stop at the last original track
344  if(phtrk_iter->first > original_track_map_lastkey) break;
345 
346  auto tracklet_tpc = phtrk_iter->second;
347 
348  if (Verbosity() >= 1)
349  {
350  std::cout << std::endl
351  << __LINE__
352  << ": Processing seed itrack: " << phtrk_iter->first
353  << ": nhits: " << tracklet_tpc-> size_cluster_keys()
354  << ": Total tracks: " << _track_map->size()
355  << ": phi: " << tracklet_tpc->get_phi()
356  << std::endl;
357  }
358 
359  // Get the outermost TPC clusters for this tracklet
360  std::map<unsigned int, TrkrCluster*> outer_clusters;
361  std::vector<TrkrCluster*> clusters;
362  std::vector<Acts::Vector3D> clusGlobPos;
363  ActsTransformations transformer;
364 
365  for (SvtxTrack::ConstClusterKeyIter key_iter = tracklet_tpc->begin_cluster_keys(); key_iter != tracklet_tpc->end_cluster_keys(); ++key_iter)
366  {
367  TrkrDefs::cluskey cluster_key = *key_iter;
368  unsigned int layer = TrkrDefs::getLayer(cluster_key);
369 
370  if(layer < _min_tpc_layer) continue;
371  if(layer >= _min_mm_layer) continue;
372 
373  // get the cluster
374  TrkrCluster *tpc_clus = _cluster_map->findCluster(cluster_key);
375  if(!tpc_clus) continue;
376 
377  outer_clusters.insert(std::make_pair(layer, tpc_clus));
378  clusters.push_back(tpc_clus);
379  clusGlobPos.push_back(transformer.getGlobalPosition(tpc_clus,_surfmaps,
380  _tGeometry));
381  if(Verbosity() > 10)
382  {
383  std::cout
384  << " TPC cluster in layer " << layer << " with position " << tpc_clus->getLocalX()
385  << " " << tpc_clus->getLocalY() << " " << " outer_clusters.size() " << outer_clusters.size() << std::endl;
386  }
387  }
388 
389  // need at least 3 clusters to fit a circle
390  if(outer_clusters.size() < 3)
391  {
392  if(Verbosity() > 3) std::cout << PHWHERE << " -- skip this tpc tracklet, not enough outer clusters " << std::endl;
393  continue; // skip to the next TPC tracklet
394  }
395 
396  // fit a circle to the clusters
397  double R = 0;
398  double X0 = 0;
399  double Y0 = 0;
400  CircleFitByTaubin(clusGlobPos, R, X0, Y0);
401  if(Verbosity() > 10) std::cout << " Fitted circle has R " << R << " X0 " << X0 << " Y0 " << Y0 << std::endl;
402 
403  // toss tracks for which the fitted circle could not have come from the vertex
404  if(R < 40.0) continue;
405 
406  // get the straight line representing the z trajectory in the form of z vs radius
407  double A = 0; double B = 0;
408  line_fit(clusGlobPos, A, B);
409  if(Verbosity() > 10) std::cout << " Fitted line has A " << A << " B " << B << std::endl;
410 
411  // loop over micromegas layer
412  for(unsigned int imm = 0; imm < _n_mm_layers; ++imm)
413  {
414 
415  // get micromegas geometry object
416  const unsigned int layer = _min_mm_layer + imm;
417  const auto layergeom = static_cast<CylinderGeomMicromegas*>(_geomContainerMicromegas->GetLayerGeom(layer));
418  const auto layer_radius = layergeom->get_radius();
419 
420  // method to find where fitted circle intersects this layer
421  double xplus = 0;
422  double xminus = 0;
423  double yplus = 0;
424  double yminus = 0;
425 
426  // finds the intersection of the fitted circle with the micromegas layer
427  if( !circle_circle_intersection( layer_radius, R, X0, Y0, xplus, yplus, xminus, yminus) )
428  {
429  if(Verbosity() > 10)
430  {
431  std::cout << PHWHERE << " circle/circle intersection calculation failed, skip this case" << std::endl;
432  std::cout << PHWHERE << " mm_radius " << layer_radius << " fitted R " << R << " fitted X0 " << X0 << " fitted Y0 " << Y0 << std::endl;
433  }
434 
435  continue;
436  }
437 
438  // we can figure out which solution is correct based on the last cluster position in the TPC
439  const double last_clus_phi = std::atan2(clusGlobPos.back()(1), clusGlobPos.back()(0));
440  double phi_plus = std::atan2(yplus, xplus);
441  double phi_minus = std::atan2(yminus, xminus);
442 
443  // calculate z
444  double r = layer_radius;
445  double z = B + A*r;
446 
447  // select the angle that is the closest to last cluster
448  // store phi, apply coarse space charge corrections in calibration mode
449  double phi = std::abs(last_clus_phi - phi_plus) < std::abs(last_clus_phi - phi_minus) ? phi_plus:phi_minus;
450  if(_sc_calib_mode) phi -= fdrphi->Eval(z)/r;
451 
452  // create cylinder intersection point in world coordinates
453  const TVector3 world_intersection_cylindrical( r*std::cos(phi), r*std::sin(phi), z );
454 
455  // find matching tile
456  int tileid = layergeom->find_tile_cylindrical( world_intersection_cylindrical );
457  if( tileid < 0 ) continue;
458 
459  // get tile coordinates
460  const auto tile_center_phi = layergeom->get_center_phi( tileid );
461  const double x0 = r*std::cos( tile_center_phi );
462  const double y0 = r*std::sin( tile_center_phi );
463 
464  const double nx = x0;
465  const double ny = y0;
466 
467  // calculate intersection to tile
468  if( !circle_line_intersection( R, X0, Y0, x0, y0, nx, ny, xplus, yplus, xminus, yminus) )
469  {
470  if(Verbosity() > 10)
471  { std::cout << PHWHERE << "circle_line_intersection - failed" << std::endl; }
472  continue;
473  }
474 
475  // select again angle closest to last cluster
476  phi_plus = std::atan2(yplus, xplus);
477  phi_minus = std::atan2(yminus, xminus);
478  const bool is_plus = (std::abs(last_clus_phi - phi_plus) < std::abs(last_clus_phi - phi_minus));
479 
480  // calculate x, y and z
481  const double x = (is_plus ? xplus:xminus);
482  const double y = (is_plus ? yplus:yminus);
483  r = get_r( x, y );
484  z = B + A*r;
485 
486  /*
487  * create planar intersection point in world coordinates
488  * this is the position to be compared to the clusters
489  */
490  const TVector3 world_intersection_planar( x, y, z );
491 
492  // convert to tile local reference frame, apply SC correction
493  TVector3 local_intersection_planar = layergeom->get_local_from_world_coords( tileid, world_intersection_planar );
494  if(_sc_calib_mode)
495  {
496  /*
497  * apply SC correction to the local intersection,
498  * to make sure that the corrected intersection is still in the micromegas plane
499  * in local tile coordinates, the rphi direction, to which the correction is applied, corresponds to the x direction
500  */
501  local_intersection_planar.SetX( local_intersection_planar.x() - fdrphi->Eval(z) );
502  }
503 
504  // store segmentation type
505  const auto segmentation_type = layergeom->get_segmentation_type();
506 
507  // generate tilesetid and get corresponding clusters
508  const auto tilesetid = MicromegasDefs::genHitSetKey(layer, segmentation_type, tileid);
509  const auto mm_clusrange = _cluster_map->getClusters(tilesetid);
510 
511  // convert to tile local coordinate and compare
512  for(auto clusiter = mm_clusrange.first; clusiter != mm_clusrange.second; ++clusiter)
513  {
514  TrkrDefs::cluskey ckey = clusiter->first;
515  if(_iteration_map != NULL ){
516  if( _iteration_map->getIteration(ckey) > 0)
517  continue; // skip hits used in a previous iteration
518  }
519  // store cluster and key
520  const auto& [key, cluster] = *clusiter;
521  const auto glob = transformer.getGlobalPosition(cluster,_surfmaps,_tGeometry);
522  const TVector3 world_cluster(glob(0), glob(1), glob(2));
523  const TVector3 local_cluster = layergeom->get_local_from_world_coords( tileid, world_cluster );
524 
525  // compute residuals and store
526  /* in local tile coordinate, x is along rphi, and z is along z) */
527  const double drphi = local_intersection_planar.x() - local_cluster.x();
528  const double dz = local_intersection_planar.z() - local_cluster.z();
529 
530  // compare to cuts and add to track if matching
531  if( std::abs(drphi) < _rphi_search_win[imm] && std::abs(dz) < _z_search_win[imm] )
532  {
533  tracklet_tpc->insert_cluster_key(key);
534  _assoc_container->SetClusterTrackAssoc(key, tracklet_tpc->get_id());
535 
536  // prints out a line that can be grep-ed from the output file to feed to a display macro
537  if( _test_windows )
538  {
539  // cluster rphi and z
540  const double mm_clus_rphi = get_r( glob(0), glob(1) ) * std::atan2( glob(1), glob(0) );
541  const double mm_clus_z = glob(2);
542 
543  // projection phi and z, without correction
544  const double rphi_proj = get_r( world_intersection_planar.x(), world_intersection_planar.y() ) * std::atan2( world_intersection_planar.y(), world_intersection_planar.x() );
545  const double z_proj = world_intersection_planar.z();
546 
547  /*
548  * Note: drphi and dz might not match the difference of the rphi and z quoted values. This is because
549  * 1/ drphi and dz are actually calculated in Tile's local reference frame, not in world coordinates
550  * 2/ drphi also includes SC distortion correction, which the world coordinates don't
551  */
552  std::cout
553  << " Try_mms: " << (int) layer
554  << " drphi " << drphi
555  << " dz " << dz
556  << " mm_clus_rphi " << mm_clus_rphi << " mm_clus_z " << mm_clus_z
557  << " rphi_proj " << rphi_proj << " z_proj " << z_proj
558  << std::endl;
559  }
560  }
561 
562  } // end loop over clusters
563 
564  } // end loop over Micromegas layers
565 
566  if(Verbosity() > 3)
567  { tracklet_tpc->identify(); }
568 
569  }
570 
571  // loop over all tracks and copy the silicon clusters to the corrected cluster map
574 
575  if(Verbosity() > 0)
576  { std::cout << " Final track map size " << _track_map->size() << std::endl; }
577 
579 }
580 
581 //_________________________________________________________________________________________________
584 
585 //_________________________________________________________________________________________________
587 {
588  // tpc-distortion-corrected clusters
589  _corrected_cluster_map = findNode::getClass<TrkrClusterContainer>(topNode,"CORRECTED_TRKR_CLUSTER");
591  {
592  std::cout << " Found CORRECTED_TRKR_CLUSTER node " << std::endl;
593  }
594 
595  // all clusters
597  _cluster_map = findNode::getClass<TrkrClusterContainer>(topNode, "TRKR_CLUSTER_TRUTH");
598  else
599  {
600  _cluster_map = findNode::getClass<TrkrClusterContainer>(topNode, "TRKR_CLUSTER");
601  }
602 
603  if (!_cluster_map)
604  {
605  std:: cerr << PHWHERE << " ERROR: Can't find node TRKR_CLUSTER" << std::endl;
607  }
608 
609 
610  _tGeometry = findNode::getClass<ActsTrackingGeometry>(topNode, "ActsTrackingGeometry");
611  if(!_tGeometry)
612  {
613  std::cerr << PHWHERE << "No acts tracking geometry, can't continue."
614  << std::endl;
616  }
617 
618  _surfmaps = findNode::getClass<ActsSurfaceMaps>(topNode, "ActsSurfaceMaps");
619  if(!_surfmaps)
620  {
621  std::cerr << PHWHERE << "No acts surface son node tree, exiting."
622  << std::endl;
624  }
625 
626  _track_map = findNode::getClass<SvtxTrackMap>(topNode, "SvtxTrackMap");
627  if (!_track_map)
628  {
629  std::cerr << PHWHERE << " ERROR: Can't find SvtxTrackMap" << std::endl;
631  }
632 
633  _assoc_container = findNode::getClass<AssocInfoContainer>(topNode, "AssocInfoContainer");
634  if (!_assoc_container)
635  {
636  std::cerr << PHWHERE << " ERROR: Can't find AssocInfoContainer." << std::endl;
638  }
639 
640  // micromegas geometry
641  _geomContainerMicromegas = findNode::getClass<PHG4CylinderGeomContainer>(topNode, "CYLINDERGEOM_MICROMEGAS_FULL" );
643  {
644  std::cout << PHWHERE << "Could not find CYLINDERGEOM_MICROMEGAS_FULL." << std::endl;
646  }
647 
649 }
650 
652 {
653  // loop over final track map, copy silicon clusters to corrected cluster map
654  for (auto phtrk_iter = _track_map->begin();
655  phtrk_iter != _track_map->end();
656  ++phtrk_iter)
657  {
658  SvtxTrack *track = phtrk_iter->second;
659 
660  // loop over associated clusters to get keys for silicon cluster
662  iter != track->end_cluster_keys();
663  ++iter)
664  {
665  TrkrDefs::cluskey cluster_key = *iter;
666  const unsigned int trkrid = TrkrDefs::getTrkrId(cluster_key);
667  if(trkrid == TrkrDefs::micromegasId)
668  {
669  TrkrCluster *cluster = _cluster_map->findCluster(cluster_key);
670  if( !cluster ) continue;
671 
672  TrkrCluster *newclus = _corrected_cluster_map->findOrAddCluster(cluster_key)->second;
673  newclus->CopyFrom( cluster );
674  }
675  }
676  }
677 }
678