52 for(
int i=0; i < 16; ++i)
57 for(
int i=0; i < 16; ++i)
62 for(
int i=0; i < 16; ++i)
83 _track = phtrk_iter->second;
87 std::cout << std::endl
89 <<
": Processing track itrack: " << phtrk_iter->first
90 <<
": nhits: " <<
_track-> size_cluster_keys()
97 std::vector<Acts::Vector3D> globalClusterPositions;
98 std::map<TrkrDefs::cluskey, Acts::Vector3D> tpc_clusters;
114 if( !cluster )
continue;
118 newclus->CopyFrom( cluster );
129 if(
Verbosity() > 2) std::cout <<
" layer " << layer <<
" distorted cluster position: " << global[0] <<
" " << global[1] <<
" " << global[2];
131 if(
Verbosity() > 2) std::cout <<
" corrected cluster position: " << global[0] <<
" " << global[1] <<
" " << global[2] << std::endl;
134 globalClusterPositions.push_back(global);
135 tpc_clusters.insert(std::make_pair(cluster_key, global));
140 if(globalClusterPositions.size() < 3)
142 if(
Verbosity() > 3) std::cout <<
PHWHERE <<
" -- skip this tpc track, not enough clusters: " << globalClusterPositions.size() << std::endl;
152 std::cout <<
" Fitted circle has R " << R <<
" X0 " << X0 <<
" Y0 " << Y0 << std::endl;
158 double A = 0;
double B = 0;
159 line_fit(globalClusterPositions, A, B);
161 std::cout <<
" Fitted line has A " << A <<
" B " << B << std::endl;
164 for (
auto clus_iter = tpc_clusters.begin();
165 clus_iter != tpc_clusters.end();
177 _z_proj = B + A * target_radius;
180 double cluster_radius = sqrt(global[0] * global[0] + global[1] * global[1]);
209 std::cout <<
PHWHERE <<
"Failed to find surface for cluster " << cluskey << std::endl;
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];
253 localPos(0) = rClusPhi - surfRphiCenter;
254 localPos(1) = znew - surfZCenter;
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;
282 _cluster_map = findNode::getClass<TrkrClusterContainer>(topNode,
"TRKR_CLUSTER");
285 std::cout <<
PHWHERE <<
" ERROR: Can't find node TRKR_CLUSTER" << std::endl;
289 _track_map = findNode::getClass<SvtxTrackMap>(topNode,
"SvtxTrackMap");
292 std::cout <<
PHWHERE <<
" ERROR: Can't find SvtxTrackMap: " << std::endl;
296 _surfmaps = findNode::getClass<ActsSurfaceMaps>(topNode,
"ActsSurfaceMaps");
299 std::cout <<
PHWHERE <<
"Error, can't find acts surface maps" << std::endl;
303 _tGeometry = findNode::getClass<ActsTrackingGeometry>(topNode,
"ActsTrackingGeometry");
306 std::cout <<
PHWHERE <<
"Error, can't find acts tracking geometry" << std::endl;
311 _dcc = findNode::getClass<TpcDistortionCorrectionContainer>(topNode,
"TpcDistortionCorrectionContainer");
314 std::cout <<
"PHTpcClusterMover: found TPC distortion correction container" << std::endl;
321 std::cout <<
"Creating node CORRECTED_TRKR_CLUSTER" << std::endl;
329 std::cout <<
PHWHERE <<
"DST Node missing, doing nothing." << std::endl;
344 DetNode->
addNode(TrkrClusterContainerNode);
367 if(std::isnan(xplus))
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;
378 if(fabs(xclus - xplus) < 5.0 && fabs(yclus - yplus) < 5.0)
406 double Mz,Mxy,Mxx,Myy,Mxz,Myz,Mzz,Cov_xy,Var_z;
407 double A0,A1,A2,
A22,A3,
A33;
409 double DET,Xcenter,Ycenter;
415 for(
unsigned int iclus = 0; iclus < clusters.size(); ++iclus)
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];
427 Mxx=Myy=Mxy=Mxz=Myz=Mzz=0.;
429 for (
unsigned int i=0; i<clusters.size(); i++)
431 double Xi = clusters[i][0] - meanX;
432 double Yi = clusters[i][1] - meanY;
433 double Zi = Xi*Xi + Yi*Yi;
452 Cov_xy = Mxx*Myy - Mxy*Mxy;
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;
465 for (x=0.,y=A0,iter=0; iter<IterMAX; iter++)
467 double Dy = A1 + x*(A22 + A33*
x);
468 double xnew = x - y/
Dy;
470 double ynew = A0 + xnew*(A1 + xnew*(A2 + xnew*A3));
471 if (fabs(ynew)>=fabs(y))
break;
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;
483 X0 = Xcenter + meanX;
484 Y0 = Ycenter + meanY;
485 R = sqrt(Xcenter*Xcenter + Ycenter*Ycenter + Mz);
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;
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);
514 yplus = - (2*x2*xplus -
D) / (2.0*y2);
515 yminus = -(2*x2*xminus -
D) / (2.0*y2);
524 double xsum=0,x2sum=0,ysum=0,xysum=0;
525 for (
unsigned int i=0; i<clusters.size(); ++i)
527 double z = clusters[i][2];
528 double r = sqrt(pow(clusters[i][0],2) + pow(clusters[i][1], 2));
532 x2sum=x2sum+pow(r,2);
535 a=(clusters.size()*xysum-xsum*ysum)/(clusters.size()*x2sum-xsum*xsum);
536 b=(x2sum*ysum-xsum*xysum)/(x2sum*clusters.size()-xsum*xsum);
548 std::map<unsigned int, std::vector<Surface>>::iterator mapIter;
554 <<
"Error: hitsetkey not found in clusterSurfaceMap, hitsetkey = "
555 << hitsetkey << std::endl;
559 double world_phi = atan2(world[1], world[0]);
560 double world_z = world[2];
562 std::vector<Surface> surf_vec = mapIter->second;
563 unsigned int surf_index = 999;
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;
571 nsurf += surf_vec.size()/2;
573 Surface this_surf = surf_vec[nsurf];
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};
577 double surf_z = surf_center[2];
578 double surf_phi = atan2(surf_center[1], surf_center[0]);
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) )
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
598 <<
"Error: TPC surface index not defined, skipping cluster!"
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;
609 return surf_vec[surf_index];