23 #if __cplusplus < 201402L
24 #include <boost/make_unique.hpp>
40 template<
class T>
inline constexpr
T square(
const T&
x ) {
return x*
x; }
45 void line_fit(std::vector<std::pair<double,double>> points,
double &
a,
double &
b)
51 double xsum=0,x2sum=0,ysum=0,xysum=0;
52 for(
const auto& [
r,
z]:points )
59 a=(points.size()*xysum-xsum*ysum)/(points.size()*x2sum-
square(xsum));
60 b=(x2sum*ysum-xsum*xysum)/(points.size()*x2sum-
square(xsum));
65 void line_fit_clusters(std::vector<Acts::Vector3D>& globPos,
double &a,
double &b)
67 std::vector<std::pair<double,double>> points;
71 line_fit(points, a, b);
74 void CircleFitByTaubin (std::vector<Acts::Vector3D> points,
double &
R,
double &X0,
double &Y0)
86 double Mz,Mxy,Mxx,Myy,Mxz,Myz,Mzz,Cov_xy,Var_z;
87 double A0,A1,A2,
A22,A3,
A33;
89 double DET,Xcenter,Ycenter;
106 Mxx=Myy=Mxy=Mxz=Myz=Mzz=0.;
112 double Zi = Xi*Xi + Yi*Yi;
131 Cov_xy = Mxx*Myy - Mxy*Mxy;
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;
144 for (x=0.,y=A0,iter=0; iter<IterMAX; ++iter)
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));
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;
162 X0 = Xcenter + meanX;
163 Y0 = Ycenter + meanY;
199 const double R,
const double X0,
200 const double Y0,
double& x,
217 * pow(Y0,4)) +
square(X0)*Y0 + pow(Y0, 3))
221 * pow(Y0,4)) +
square(X0) * Y0 + pow(Y0, 3))
224 const double minx = std::sqrt(
square(R) -
square(miny-Y0)) + X0;
225 const double minx2 = -std::sqrt(
square(R) -
square(miny2-Y0)) + X0;
273 const auto& track_key = phtrk_iter->first;
274 SvtxTrack* tracklet_tpc = phtrk_iter->second;
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()
292 if(
Verbosity() > 3) std::cout <<
" TPC tracklet " << tracklet_tpc->
get_id() <<
" clusters.size " << clusters.size() << std::endl;
295 bool reject_track =
false;
296 std::set<unsigned int> layers;
297 for (
unsigned int i=0; i<clusters.size(); ++i)
301 if(
Verbosity() > 0) std::cout <<
" trackid " << phtrk_iter->first <<
" no cluster found, skip track" << std::endl;
306 layers.insert(layer);
309 if(reject_track)
continue;
311 unsigned int nlayers = layers.size();
312 if(
Verbosity() > 2) std::cout <<
" TPC layers this track: " << nlayers << std::endl;
315 std::vector<Acts::Vector3D> globalClusterPositions;
316 for (
unsigned int i=0; i<clusters.size(); ++i)
319 globalClusterPositions.push_back(global);
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;
334 if(clusters.size() < 3)
336 if(
Verbosity() > 3) std::cout <<
PHWHERE <<
" -- skip this tpc tracklet, not enough TPC clusters " << std::endl;
343 line_fit_clusters(globalClusterPositions, A, B);
345 { std::cout <<
"PHTpcTrackSeedCircleFit::process_event - track: " << track_key <<
" A=" << A <<
" B=" << B << std::endl; }
348 const double z_proj =
B;
351 tracklet_tpc->
set_z(z_proj);
356 double track_angle = atan(A);
361 CircleFitByTaubin(globalClusterPositions, R, X0, Y0);
363 { std::cout <<
"PHTpcTrackSeedCircleFit::process_event - track: " << track_key <<
" R=" << R <<
" X0=" << X0 <<
" Y0=" << Y0 << std::endl; }
367 findRoot(R, X0, Y0, dcax, dcay);
368 if(std::isnan(dcax) or std::isnan(dcay))
371 tracklet_tpc->
set_x(dcax);
372 tracklet_tpc->
set_y(dcay);
374 double pt_track = tracklet_tpc->
get_pt();
378 static constexpr
double Bfield = 1.4;
379 double pt_circle = 0.3 * Bfield * R * 0.01;
380 std::cout <<
"PHTpcTrackSeedCircleFit::process_event-"
382 <<
" Bfield: " << Bfield
383 <<
" pt_circle: " << pt_circle <<
" pt_track: " << pt_track
390 double dx = X0 - dcax;
391 double dy = Y0 - dcay;
392 double phi= atan2(-dx,dy);
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;
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;
420 { phi -= 2. *
M_PI; }
424 std::cout <<
" input track phi " << tracklet_tpc->
get_phi() <<
" new phi " << phi << std::endl;
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);
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
443 tracklet_tpc->
set_px(px_new);
444 tracklet_tpc->
set_py(py_new);
445 tracklet_tpc->
set_pz(pz_new);
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;
455 std::cout <<
" Final track map size " <<
_track_map->
size() << std::endl;
458 std::cout <<
"PHTpcTrackSeedCircleFit::process_event(PHCompositeNode *topNode) Leaving process_event" << std::endl;
468 _surfmaps = findNode::getClass<ActsSurfaceMaps>(topNode,
"ActsSurfaceMaps");
471 std::cout <<
PHWHERE <<
"Error, can't find acts surface maps" << std::endl;
475 _tGeometry = findNode::getClass<ActsTrackingGeometry>(topNode,
"ActsTrackingGeometry");
478 std::cout <<
PHWHERE <<
"Error, can't find acts tracking geometry" << std::endl;
483 _cluster_map = findNode::getClass<TrkrClusterContainer>(topNode,
"TRKR_CLUSTER_TRUTH");
486 _cluster_map = findNode::getClass<TrkrClusterContainer>(topNode,
"CORRECTED_TRKR_CLUSTER");
489 std::cout <<
" using CORRECTED_TRKR_CLUSTER node "<< std::endl;
493 std::cout <<
" CORRECTED_TRKR_CLUSTER node not found, using TRKR_CLUSTER " << std::endl;
494 _cluster_map = findNode::getClass<TrkrClusterContainer>(topNode,
"TRKR_CLUSTER");
500 std::cerr <<
PHWHERE <<
" ERROR: Can't find node TRKR_CLUSTER" << std::endl;
505 _dcc = findNode::getClass<TpcDistortionCorrectionContainer>(topNode,
"TpcDistortionCorrectionContainer");
507 { std::cout <<
"PHTpcTrackSeedCircleFit::get_Nodes - found TPC distortion correction container" << std::endl; }
512 std::cerr <<
PHWHERE <<
" ERROR: Can't find SvtxTrackMap" << std::endl;
521 std::vector<TrkrCluster*> clusters;
536 clusters.push_back(tpc_clus);
541 std::cout <<
" TPC cluster in layer " << layer <<
" with local position " << tpc_clus->
getLocalX()
542 <<
" " << tpc_clus->
getLocalY() <<
" clusters.size() " << clusters.size() << std::endl;