ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
PHTpcClusterMover.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file PHTpcClusterMover.cc
1 #include "PHTpcClusterMover.h"
2 
3 #include "PHTpcClusterMover.h"
4 
6 
7 #include <trackbase/TrkrClusterv3.h> // for TrkrCluster
8 #include <trackbase/TrkrDefs.h> // for cluskey, getLayer, TrkrId
10 #include <trackbase_historic/SvtxTrack.h> // for SvtxTrack, SvtxTrack::C...
13 
14 
16 
17 #include <phool/getClass.h>
18 #include <phool/phool.h>
19 #include <phool/PHCompositeNode.h>
20 
21 #include <TF1.h>
22 
23 #include <cmath> // for sqrt, fabs, atan2, cos
24 #include <iostream> // for operator<<, basic_ostream
25 #include <map> // for map
26 #include <set> // for _Rb_tree_const_iterator
27 #include <utility> // for pair, make_pair
28 
29 //____________________________________________________________________________..
31  : SubsysReco(name)
32 {
33 
34 }
35 
36 //____________________________________________________________________________..
38 {
39 
40 }
41 
42 //____________________________________________________________________________..
44 {
45  int ret = GetNodes(topNode);
46  if (ret != Fun4AllReturnCodes::EVENT_OK) return ret;
47 
48  // initialize layer radii
52  for(int i=0; i < 16; ++i)
53  {
55  if(Verbosity() > 4) std::cout << " i " << i << " layer_radius " << layer_radius[i] << std::endl;
56  }
57  for(int i=0; i < 16; ++i)
58  {
59  layer_radius[i+16] = mid_tpc_min_radius + (double) i * mid_tpc_spacing + 0.5 * mid_tpc_spacing;
60  if(Verbosity() > 4) std::cout << " i " << i << " layer_radius " << layer_radius[i+16] << std::endl;
61  }
62  for(int i=0; i < 16; ++i)
63  {
65  if(Verbosity() > 4) std::cout << " i " << i << " layer_radius " << layer_radius[i+32] << std::endl;
66  }
67 
68  return ret;
69 }
70 
71 //____________________________________________________________________________..
73 {
74 
75  if(Verbosity() > 0)
76  std::cout << PHWHERE << " track map size " << _track_map->size() << std::endl;
77 
78  // loop over the tracks
79  for (auto phtrk_iter = _track_map->begin();
80  phtrk_iter != _track_map->end();
81  ++phtrk_iter)
82  {
83  _track = phtrk_iter->second;
84 
85  if (Verbosity() >= 1)
86  {
87  std::cout << std::endl
88  << __LINE__
89  << ": Processing track itrack: " << phtrk_iter->first
90  << ": nhits: " << _track-> size_cluster_keys()
91  << ": Total tracks: " << _track_map->size()
92  << ": phi: " << _track->get_phi()
93  << std::endl;
94  }
95 
96  // Get the TPC clusters for this track and correct them for distortions
97  std::vector<Acts::Vector3D> globalClusterPositions;
98  std::map<TrkrDefs::cluskey, Acts::Vector3D> tpc_clusters;
99 
101  key_iter != _track->end_cluster_keys();
102  ++key_iter)
103  {
104  TrkrDefs::cluskey cluster_key = *key_iter;
105  unsigned int trkrId = TrkrDefs::getTrkrId(cluster_key);
106  unsigned int layer = TrkrDefs::getLayer(cluster_key);
107 
108  // non Tpc clusters are copied unchanged to the new map if they are present
109  // This is needed for truth seeding case, where the tracks already have silicon clusters
110  if(trkrId != TrkrDefs::tpcId)
111  {
112  // get cluster from original map
113  auto cluster = _cluster_map->findCluster(cluster_key);
114  if( !cluster ) continue;
115 
116  // create in corrected map and copy content
117  auto newclus = _corrected_cluster_map->findOrAddCluster(cluster_key)->second;
118  newclus->CopyFrom( cluster );
119  continue;
120  }
121 
122  // get the cluster in 3D coordinates
123  TrkrCluster *tpc_clus = _cluster_map->findCluster(cluster_key);
124  auto global = _transformer.getGlobalPosition(tpc_clus,
125  _surfmaps,
126  _tGeometry);
127 
128  // check if TPC distortion correction are in place and apply
129  if(Verbosity() > 2) std::cout << " layer " << layer << " distorted cluster position: " << global[0] << " " << global[1] << " " << global[2];
130  if( _dcc ) global = _distortionCorrection.get_corrected_position( global, _dcc );
131  if(Verbosity() > 2) std::cout << " corrected cluster position: " << global[0] << " " << global[1] << " " << global[2] << std::endl;
132 
133  // Store the corrected 3D cluster positions
134  globalClusterPositions.push_back(global);
135  tpc_clusters.insert(std::make_pair(cluster_key, global));
136 
137  }
138 
139  // need at least 3 clusters to fit a circle
140  if(globalClusterPositions.size() < 3)
141  {
142  if(Verbosity() > 3) std::cout << PHWHERE << " -- skip this tpc track, not enough clusters: " << globalClusterPositions.size() << std::endl;
143  continue; // skip to the next TPC track
144  }
145 
146  // fit a circle to the clusters
147  double R = 0;
148  double X0 = 0;
149  double Y0 = 0;
150  CircleFitByTaubin(globalClusterPositions, R, X0, Y0);
151  if(Verbosity() > 10)
152  std::cout << " Fitted circle has R " << R << " X0 " << X0 << " Y0 " << Y0 << std::endl;
153 
154  // toss tracks for which the fitted circle could not have come from the vertex
155  //if(R < 30.0) continue;
156 
157  // get the straight line representing the z trajectory in the form of z vs radius
158  double A = 0; double B = 0;
159  line_fit(globalClusterPositions, A, B);
160  if(Verbosity() > 10)
161  std::cout << " Fitted line has A " << A << " B " << B << std::endl;
162 
163  // Now we need to move each cluster associated with this track to the readout layer radius
164  for (auto clus_iter = tpc_clusters.begin();
165  clus_iter != tpc_clusters.end();
166  ++clus_iter)
167  {
168  TrkrDefs::cluskey cluskey = clus_iter->first;
169  unsigned int layer = TrkrDefs::getLayer(cluskey);
170  Acts::Vector3D global = clus_iter->second;
171 
172  // get circle position at target surface radius
173  double target_radius = layer_radius[layer-7];
174  int ret = get_circle_circle_intersection(target_radius, R, X0, Y0, global[0], global[1], _x_proj, _y_proj);
175  if(ret == Fun4AllReturnCodes::ABORTEVENT) continue; // skip to next cluster
176  // z projection is unique
177  _z_proj = B + A * target_radius;
178 
179  // get circle position at cluster radius
180  double cluster_radius = sqrt(global[0] * global[0] + global[1] * global[1]);
181  ret = get_circle_circle_intersection(cluster_radius, R, X0, Y0, global[0], global[1], _x_start, _y_start);
182  if(ret == Fun4AllReturnCodes::ABORTEVENT) continue; // skip to next cluster
183  // z projection is unique
184  _z_start = B + A * cluster_radius;
185 
186  // calculate dx, dy, dz along circle trajectory from cluster radius to surface radius
187  double xnew = global[0] - (_x_start - _x_proj);
188  double ynew = global[1] - (_y_start - _y_proj);
189  double znew = global[2] - (_z_start - _z_proj);
190 
191  // now move the cluster to the surface radius
192  // we keep the cluster key fixed, change the surface if necessary
193  // write the new cluster position local coordinates on the surface
194 
195  Acts::Vector3D global_new(xnew, ynew, znew);
196 
200  global_new,
201  _surfmaps,
202  _tGeometry,
203  subsurfkey);
204 
205  if(!surface)
206  {
209  std::cout << PHWHERE << "Failed to find surface for cluster " << cluskey << std::endl;
210  continue;
211  }
212 
213  // get the original cluster
214  TrkrCluster *cluster = _cluster_map->findCluster(cluskey);
215 
216  // put the corrected cluster in the new cluster map
217 
218  // ghost tracks can have repeat clusters, so check if cluster already moved
219  if(_corrected_cluster_map->findCluster(cluskey)) continue;
220 
221  TrkrCluster *newclus = _corrected_cluster_map->findOrAddCluster(cluskey)->second;
222  newclus->setSubSurfKey(subsurfkey);
223  newclus->setAdc(cluster->getAdc());
224 
225  newclus->setActsLocalError(0,0,cluster->getActsLocalError(0,0));
226  newclus->setActsLocalError(1,0,cluster->getActsLocalError(1,0));
227  newclus->setActsLocalError(0,1,cluster->getActsLocalError(0,1));
228  newclus->setActsLocalError(1,1,cluster->getActsLocalError(1,1));
229 
230  // get local coordinates
231  Acts::Vector3D normal = surface->normal(_tGeometry->geoContext);
232  auto local = surface->globalToLocal(_tGeometry->geoContext,
233  global * Acts::UnitConstants::cm,
234  normal);
235 
236  Acts::Vector2D localPos;
237  if(local.ok())
238  {
239  localPos = local.value() / Acts::UnitConstants::cm;
240  }
241  else
242  {
244  Acts::Vector3D center = surface->center(_tGeometry->geoContext)/Acts::UnitConstants::cm;
245  double clusRadius = sqrt(xnew * xnew + ynew * ynew);
246  double clusphi = atan2(ynew, xnew);
247  double rClusPhi = clusRadius * clusphi;
248  double surfRadius = sqrt(center(0)*center(0) + center(1)*center(1));
249  double surfPhiCenter = atan2(center[1], center[0]);
250  double surfRphiCenter = surfPhiCenter * surfRadius;
251  double surfZCenter = center[2];
252 
253  localPos(0) = rClusPhi - surfRphiCenter;
254  localPos(1) = znew - surfZCenter;
255  }
256 
257  if(Verbosity() > 4)
258  {
259  std::cout << "*** cluster_radius " << cluster_radius << " cluster x,y,z: " << global[0] << " " << global[1] << " " << global[2] << std::endl;
260  std::cout << " projection_radius " << target_radius << " proj x,y,z: " << _x_proj << " " << _y_proj << " " << _z_proj << std::endl;
261  std::cout << " traj_start_radius " << cluster_radius << " start x,y,z: "<< _x_start << " " << _y_start << " " << _z_start << std::endl;
262  std::cout << " moved_clus_radius " << target_radius << " final x,y,z: "<< xnew << " " << ynew << " " << znew << std::endl;
263  }
264 
265  newclus->setLocalX(localPos(0));
266  newclus->setLocalY(localPos(1));
267  }
268 
269  // For normal reconstruction, the silicon clusters for this track will be copied over after the matching is done
270  }
271 
273 }
274 
276 {
278 }
279 
281 {
282  _cluster_map = findNode::getClass<TrkrClusterContainer>(topNode, "TRKR_CLUSTER");
283  if (!_cluster_map)
284  {
285  std::cout << PHWHERE << " ERROR: Can't find node TRKR_CLUSTER" << std::endl;
287  }
288 
289  _track_map = findNode::getClass<SvtxTrackMap>(topNode, "SvtxTrackMap");
290  if (!_track_map)
291  {
292  std::cout << PHWHERE << " ERROR: Can't find SvtxTrackMap: " << std::endl;
294  }
295 
296  _surfmaps = findNode::getClass<ActsSurfaceMaps>(topNode,"ActsSurfaceMaps");
297  if(!_surfmaps)
298  {
299  std::cout << PHWHERE << "Error, can't find acts surface maps" << std::endl;
301  }
302 
303  _tGeometry = findNode::getClass<ActsTrackingGeometry>(topNode,"ActsTrackingGeometry");
304  if(!_tGeometry)
305  {
306  std::cout << PHWHERE << "Error, can't find acts tracking geometry" << std::endl;
308  }
309 
310  // tpc distortion correction
311  _dcc = findNode::getClass<TpcDistortionCorrectionContainer>(topNode,"TpcDistortionCorrectionContainer");
312  if( _dcc )
313  {
314  std::cout << "PHTpcClusterMover: found TPC distortion correction container" << std::endl;
315  }
316 
317  // create the node for distortion corrected clusters, if it does not already exist
318  _corrected_cluster_map = findNode::getClass<TrkrClusterContainer>(topNode, "CORRECTED_TRKR_CLUSTER");
320  {
321  std::cout << "Creating node CORRECTED_TRKR_CLUSTER" << std::endl;
322 
323  PHNodeIterator iter(topNode);
324 
325  // Looking for the DST node
326  PHCompositeNode *dstNode = dynamic_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", "DST"));
327  if (!dstNode)
328  {
329  std::cout << PHWHERE << "DST Node missing, doing nothing." << std::endl;
331  }
332  PHNodeIterator dstiter(dstNode);
333  PHCompositeNode *DetNode =
334  dynamic_cast<PHCompositeNode *>(dstiter.findFirst("PHCompositeNode", "TRKR"));
335  if (!DetNode)
336  {
337  DetNode = new PHCompositeNode("TRKR");
338  dstNode->addNode(DetNode);
339  }
340 
342  PHIODataNode<PHObject> *TrkrClusterContainerNode =
343  new PHIODataNode<PHObject>(_corrected_cluster_map, "CORRECTED_TRKR_CLUSTER", "PHObject");
344  DetNode->addNode(TrkrClusterContainerNode);
345  }
346  else
347  {
349  }
350 
351 
352 
354 }
355 
356 int PHTpcClusterMover::get_circle_circle_intersection(double target_radius, double R, double X0, double Y0, double xclus, double yclus, double &x, double &y)
357 {
358  // finds the intersection of the fitted circle with the cylinder having radius = target_radius
359  double xplus = 0;
360  double yplus = 0;
361  double xminus = 0;
362  double yminus = 0;
363 
364  circle_circle_intersection(target_radius, R, X0, Y0, xplus, yplus, xminus, yminus);
365 
366  // We only need to check xplus for failure, skip this TPC cluster in that case
367  if(std::isnan(xplus))
368  {
369  if(Verbosity() > 0)
370  {
371  std::cout << " circle/circle intersection calculation failed, skip this cluster" << std::endl;
372  std::cout << " target_radius " << target_radius << " fitted R " << R << " fitted X0 " << X0 << " fitted Y0 " << Y0 << std::endl;
373  }
374  return Fun4AllReturnCodes::ABORTEVENT; // skip to next cluster
375  }
376 
377  // we can figure out which solution is correct based on the cluster position in the TPC
378  if(fabs(xclus - xplus) < 5.0 && fabs(yclus - yplus) < 5.0) // 5 cm, large and arbitrary
379  {
380  x = xplus;
381  y = yplus;
382  }
383  else
384  {
385  x = xminus;
386  y = yminus;
387  }
389  }
390 void PHTpcClusterMover::CircleFitByTaubin (std::vector<Acts::Vector3D> clusters, double &R, double &X0, double &Y0)
391 /*
392  Circle fit to a given set of data points (in 2D)
393  This is an algebraic fit, due to Taubin, based on the journal article
394  G. Taubin, "Estimation Of Planar Curves, Surfaces And Nonplanar
395  Space Curves Defined By Implicit Equations, With
396  Applications To Edge And Range Image Segmentation",
397  IEEE Trans. PAMI, Vol. 13, pages 1115-1138, (1991)
398  It works well whether data points are sampled along an entire circle or along a small arc.
399  It still has a small bias and its statistical accuracy is slightly lower than that of the geometric fit (minimizing geometric distances),
400  It provides a very good initial guess for a subsequent geometric fit.
401  Nikolai Chernov (September 2012)
402 */
403 {
404  int iter,IterMAX=99;
405 
406  double Mz,Mxy,Mxx,Myy,Mxz,Myz,Mzz,Cov_xy,Var_z;
407  double A0,A1,A2,A22,A3,A33;
408  double x,y;
409  double DET,Xcenter,Ycenter;
410 
411  // Compute x- and y- sample means
412  double meanX = 0;
413  double meanY = 0;
414  double weight = 0;
415  for(unsigned int iclus = 0; iclus < clusters.size(); ++iclus)
416  {
417  if(Verbosity() > 3) std::cout << " add cluster with x " << clusters[iclus][0] << " and y " << clusters[iclus][1] << std::endl;
418  meanX += clusters[iclus][0];
419  meanY += clusters[iclus][1];
420  weight++;
421  }
422  meanX /= weight;
423  meanY /= weight;
424 
425  // computing moments
426 
427  Mxx=Myy=Mxy=Mxz=Myz=Mzz=0.;
428 
429  for (unsigned int i=0; i<clusters.size(); i++)
430  {
431  double Xi = clusters[i][0] - meanX; // centered x-coordinates
432  double Yi = clusters[i][1] - meanY; // centered y-coordinates
433  double Zi = Xi*Xi + Yi*Yi;
434 
435  Mxy += Xi*Yi;
436  Mxx += Xi*Xi;
437  Myy += Yi*Yi;
438  Mxz += Xi*Zi;
439  Myz += Yi*Zi;
440  Mzz += Zi*Zi;
441  }
442  Mxx /= weight;
443  Myy /= weight;
444  Mxy /= weight;
445  Mxz /= weight;
446  Myz /= weight;
447  Mzz /= weight;
448 
449  // computing coefficients of the characteristic polynomial
450 
451  Mz = Mxx + Myy;
452  Cov_xy = Mxx*Myy - Mxy*Mxy;
453  Var_z = Mzz - Mz*Mz;
454  A3 = 4*Mz;
455  A2 = -3*Mz*Mz - Mzz;
456  A1 = Var_z*Mz + 4*Cov_xy*Mz - Mxz*Mxz - Myz*Myz;
457  A0 = Mxz*(Mxz*Myy - Myz*Mxy) + Myz*(Myz*Mxx - Mxz*Mxy) - Var_z*Cov_xy;
458  A22 = A2 + A2;
459  A33 = A3 + A3 + A3;
460 
461  // finding the root of the characteristic polynomial
462  // using Newton's method starting at x=0
463  // (it is guaranteed to converge to the right root)
464 
465  for (x=0.,y=A0,iter=0; iter<IterMAX; iter++) // usually, 4-6 iterations are enough
466  {
467  double Dy = A1 + x*(A22 + A33*x);
468  double xnew = x - y/Dy;
469  if ((xnew == x)||(!std::isfinite(xnew))) break;
470  double ynew = A0 + xnew*(A1 + xnew*(A2 + xnew*A3));
471  if (fabs(ynew)>=fabs(y)) break;
472  x = xnew; y = ynew;
473  }
474 
475  // computing parameters of the fitting circle
476 
477  DET = x*x - x*Mz + Cov_xy;
478  Xcenter = (Mxz*(Myy - x) - Myz*Mxy)/DET/2;
479  Ycenter = (Myz*(Mxx - x) - Mxz*Mxy)/DET/2;
480 
481  // assembling the output
482 
483  X0 = Xcenter + meanX;
484  Y0 = Ycenter + meanY;
485  R = sqrt(Xcenter*Xcenter + Ycenter*Ycenter + Mz);
486 }
487 
488 void PHTpcClusterMover::circle_circle_intersection(double r1, double r2, double x2, double y2, double &xplus, double &yplus, double &xminus, double &yminus)
489 {
490  // r1 is radius of sPHENIX layer
491  // r2, x2 and y2 are parameters of circle fitted to TPC clusters
492  // the solutions are xplus, xminus, yplus, yminus
493 
494  // The intersection of two circles occurs when
495  // (x-x1)^2 + (y-y1)^2 = r1^2, / (x-x2)^2 + (y-y2)^2 = r2^2
496  // Here we assume that circle 1 is an sPHENIX layer centered on x1=y1=0, and circle 2 is arbitrary
497  // x^2 +y^2 = r1^2, (x-x2)^2 + (y-y2)^2 = r2^2
498  // expand the equations and subtract to eliminate the x^2 and y^2 terms, gives the radical line connecting the intersection points
499  // iy = - (2*x2*x - D) / 2*y2,
500  // then substitute for y in equation of circle 1
501 
502  double D = r1*r1 - r2*r2 + x2*x2 + y2*y2;
503  double a = 1.0 + (x2*x2) / (y2*y2);
504  double b = - D * x2/( y2*y2);
505  double c = D*D / (4.0*y2*y2) - r1*r1;
506 
507  xplus = (-b + sqrt(b*b - 4.0* a * c) ) / (2.0 * a);
508  xminus = (-b - sqrt(b*b - 4.0* a * c) ) / (2.0 * a);
509 
510  // both values of x are valid
511  // but for each of those values, there are two possible y values on circle 1
512  // but only one of those falls on the radical line:
513 
514  yplus = - (2*x2*xplus - D) / (2.0*y2);
515  yminus = -(2*x2*xminus - D) / (2.0*y2);
516 
517 }
518 
519 void PHTpcClusterMover::line_fit(std::vector<Acts::Vector3D> clusters, double &a, double &b)
520 {
521  // copied from: https://www.bragitoff.com
522  // we want to fit z vs radius
523 
524  double xsum=0,x2sum=0,ysum=0,xysum=0; //variables for sums/sigma of xi,yi,xi^2,xiyi etc
525  for (unsigned int i=0; i<clusters.size(); ++i)
526  {
527  double z = clusters[i][2];
528  double r = sqrt(pow(clusters[i][0],2) + pow(clusters[i][1], 2));
529 
530  xsum=xsum+r; //calculate sigma(xi)
531  ysum=ysum+z; //calculate sigma(yi)
532  x2sum=x2sum+pow(r,2); //calculate sigma(x^2i)
533  xysum=xysum+r*z; //calculate sigma(xi*yi)
534  }
535  a=(clusters.size()*xysum-xsum*ysum)/(clusters.size()*x2sum-xsum*xsum); //calculate slope
536  b=(x2sum*ysum-xsum*xysum)/(x2sum*clusters.size()-xsum*xsum); //calculate intercept
537 
538  return;
539 }
540 
543  ActsSurfaceMaps *surfMaps,
546 {
547  unsigned int layer = TrkrDefs::getLayer(hitsetkey);
548  std::map<unsigned int, std::vector<Surface>>::iterator mapIter;
549  mapIter = surfMaps->tpcSurfaceMap.find(layer);
550 
551  if(mapIter == surfMaps->tpcSurfaceMap.end())
552  {
553  std::cout << PHWHERE
554  << "Error: hitsetkey not found in clusterSurfaceMap, hitsetkey = "
555  << hitsetkey << std::endl;
556  return nullptr;
557  }
558 
559  double world_phi = atan2(world[1], world[0]);
560  double world_z = world[2];
561 
562  std::vector<Surface> surf_vec = mapIter->second;
563  unsigned int surf_index = 999;
564 
565  // Predict which surface index this phi and z will correspond to
566  // assumes that the vector elements are ordered positive z, -pi to pi, then negative z, -pi to pi
567  double fraction = (world_phi + M_PI) / (2.0 * M_PI);
568  double rounded_nsurf = round( (double) (surf_vec.size()/2) * fraction - 0.5);
569  unsigned int nsurf = (unsigned int) rounded_nsurf;
570  if(world_z < 0)
571  nsurf += surf_vec.size()/2;
572 
573  Surface this_surf = surf_vec[nsurf];
574 
575  auto vec3d = this_surf->center(tGeometry->geoContext);
576  std::vector<double> surf_center = {vec3d(0) / 10.0, vec3d(1) / 10.0, vec3d(2) / 10.0}; // convert from mm to cm
577  double surf_z = surf_center[2];
578  double surf_phi = atan2(surf_center[1], surf_center[0]);
579  double surfStepPhi = tGeometry->tpcSurfStepPhi;
580  double surfStepZ = tGeometry->tpcSurfStepZ;
581 
582  if( (world_phi > surf_phi - surfStepPhi / 2.0 && world_phi < surf_phi + surfStepPhi / 2.0 ) &&
583  (world_z > surf_z - surfStepZ / 2.0 && world_z < surf_z + surfStepZ / 2.0) )
584  {
585  if(Verbosity() > 2)
586  std::cout << " got it: surf_phi " << surf_phi << " surf_z " << surf_z
587  << " surfStepPhi/2 " << surfStepPhi/2.0 << " surfStepZ/2 " << surfStepZ/2.0
588  << " world_phi " << world_phi << " world_z " << world_z
589  << " rounded_nsurf "<< rounded_nsurf << " surf_index " << nsurf
590  << std::endl;
591 
592  surf_index = nsurf;
593  subsurfkey = nsurf;
594  }
595  else
596  {
597  std::cout << PHWHERE
598  << "Error: TPC surface index not defined, skipping cluster!"
599  << std::endl;
600  std::cout << " coordinates: " << world[0] << " " << world[1] << " " << world[2]
601  << " radius " << sqrt(world[0]*world[0]+world[1]*world[1]) << std::endl;
602  std::cout << " world_phi " << world_phi << " world_z " << world_z << std::endl;
603  std::cout << " surf coords: " << surf_center[0] << " " << surf_center[1] << " " << surf_center[2] << std::endl;
604  std::cout << " surf_phi " << surf_phi << " surf_z " << surf_z << std::endl;
605  std::cout << " number of surfaces " << surf_vec.size() << std::endl;
606  return nullptr;
607  }
608 
609  return surf_vec[surf_index];
610 
611 }