33 #include <TDatabasePDG.h>
34 #include <TParticlePDG.h>
44 #define LogDebug(exp) std::cout << "DEBUG: " << __FILE__ << ": " << __LINE__ << ": " << exp << std::endl
45 #define LogError(exp) std::cout << "ERROR: " << __FILE__ << ": " << __LINE__ << ": " << exp << std::endl
46 #define LogWarning(exp) std::cout << "WARNING: " << __FILE__ << ": " << __LINE__ << ": " << exp << std::endl
54 , _clustereval(nullptr)
59 cout <<
"Enter PHTruthTrackSeeding:: Setup" << endl;
78 typedef std::map<int, std::set<TrkrCluster*> > TrkClustersMap;
79 TrkClustersMap m_trackID_clusters;
81 vector<TrkrDefs::cluskey> ClusterKeyList;
88 ClusterKeyList.clear();
91 if (g4particle==NULL){
92 cout <<__PRETTY_FUNCTION__<<
" - validity check failed: missing truth particle" <<endl;
100 const double monentum2 =
107 cout <<__PRETTY_FUNCTION__<<
" ignore low momentum particle "<< gtrackID <<endl;
117 ClusterKeyList.push_back(cluskey);
122 auto svtx_track = std::make_unique<SvtxTrack_FastSim_v2>();
124 svtx_track->set_truth_track_id(gtrackID);
127 svtx_track->set_vertex_id(vertexID-1);
137 const auto particle = TDatabasePDG::Instance()->GetParticle(g4particle->
get_pid());
148 svtx_track->set_px(g4particle->
get_px() * (1 + random));
149 svtx_track->set_py(g4particle->
get_py() * (1 + random));
150 svtx_track->set_pz(g4particle->
get_pz() * (1 + random));
154 svtx_track->
set_x(g4vertex->get_x() * (1 + random));
155 svtx_track->set_y(g4vertex->get_y() * (1 + random));
156 svtx_track->set_z(g4vertex->get_z() * (1 + random));
159 svtx_track->insert_cluster_key(
cluskey);
165 double x,
y,
z, px, py, pz;
167 px, py, pz, svtx_track->get_charge());
168 svtx_track->set_x(x);
169 svtx_track->set_y(y);
170 svtx_track->set_z(z);
171 svtx_track->set_px(px);
172 svtx_track->set_py(py);
173 svtx_track->set_pz(pz);
176 svtx_track->set_ndf(ClusterKeyList.size()*3-5);
177 svtx_track->set_chisq(1.5*ClusterKeyList.size()*3-5);
183 cout <<
"Loop over TrackMap " <<
_track_map_name <<
" entries " << endl;
192 cout <<
"Track ID: " << svtx_track->
get_id() <<
", Dummy Track pT: "
193 << svtx_track->
get_pt() <<
", Truth Track/Particle ID: "
195 <<
" (X,Y,Z) " << svtx_track->
get_x() <<
", " << svtx_track->
get_y() <<
", " << svtx_track->
get_z()
204 cout <<
"Key: " << cluster_key<< endl;
208 cout <<
" cluster ID: "
209 << cluster->
getClusKey() <<
", cluster radius: " << radius
223 m_clusterMap = findNode::getClass<TrkrClusterContainer>(topNode,
228 cerr <<
PHWHERE <<
"Error: Can't find node TRKR_CLUSTER" << endl;
232 _g4truth_container = findNode::getClass<PHG4TruthInfoContainer>(topNode,
"G4TruthInfo");
235 cerr <<
PHWHERE <<
" ERROR: Can't find node G4TruthInfo" << endl;
239 hittruthassoc = findNode::getClass<TrkrHitTruthAssoc>(topNode,
"TRKR_HITTRUTHASSOC");
242 cout <<
PHWHERE <<
"Failed to find TRKR_HITTRUTHASSOC node, quit!" << endl;
246 using nodePair = std::pair<std::string, PHG4HitContainer*&>;
247 std::initializer_list<nodePair> nodes =
255 for(
auto&& node: nodes )
257 if( !( node.second = findNode::getClass<PHG4HitContainer>( topNode, node.first ) ) )
258 { std::cerr <<
PHWHERE <<
" PHG4HitContainer " << node.first <<
" not found" << std::endl; }
272 double&
x,
double&
y,
double&
z,
273 double& px,
double& py,
double& pz,
281 double phi = atan2(-1 * (X0-x), Y0-y);
294 double theta = atan(1./m);
300 float pt = 0.3 * 1.4 * R / 100.;
301 float eta = -log(tan(theta/2.));
302 float p = pt * cosh(eta);
306 px = p * sin(theta) * cos(phi);
307 py = p * sin(theta) * sin(phi);
313 double&
A,
double&
B)
316 double xsum = 0,x2sum = 0,ysum = 0,xysum = 0;
317 for(
auto& clusterkey : clusters)
320 double z = cluster->getZ();
321 double r = sqrt(pow(cluster->getX(),2) + pow(cluster->getY(), 2));
325 x2sum=x2sum+pow(r,2);
330 A = (clusters.size()*xysum-xsum*ysum) / (clusters.size()*x2sum-xsum*xsum);
333 B = (x2sum*ysum-xsum*xysum) / (x2sum*clusters.size()-xsum*xsum);
339 double miny = (sqrt(pow(X0, 2) * pow(R, 2) * pow(Y0, 2) + pow(R, 2)
340 * pow(Y0,4)) + pow(X0,2) * Y0 + pow(Y0, 3))
341 / (pow(X0, 2) + pow(Y0, 2));
343 double miny2 = (-sqrt(pow(X0, 2) * pow(R, 2) * pow(Y0, 2) + pow(R, 2)
344 * pow(Y0,4)) + pow(X0,2) * Y0 + pow(Y0, 3))
345 / (pow(X0, 2) + pow(Y0, 2));
347 double minx = sqrt(pow(R, 2) - pow(miny - Y0, 2)) + X0;
348 double minx2 = -sqrt(pow(R, 2) - pow(miny2 - Y0, 2)) + X0;
351 if(fabs(minx) < fabs(minx2))
356 if(fabs(miny) < fabs(miny2))
364 double&
R,
double& X0,
double& Y0)
381 int iter, IterMAX=99;
383 double Mz, Mxy, Mxx, Myy, Mxz, Myz, Mzz, Cov_xy, Var_z;
384 double A0, A1, A2,
A22, A3,
A33;
386 double DET, Xcenter, Ycenter;
396 meanX += clus->
getX();
397 meanY += clus->getY();
403 Mxx=Myy=Mxy=Mxz=Myz=Mzz=0.;
409 double Xi = clus->
getX() - meanX;
410 double Yi = clus->getY() - meanY;
411 double Zi = Xi * Xi + Yi * Yi;
429 Cov_xy = Mxx * Myy - Mxy * Mxy;
430 Var_z = Mzz - Mz * Mz;
432 A2 = -3 * Mz * Mz - Mzz;
433 A1 = Var_z * Mz + 4 * Cov_xy * Mz - Mxz * Mxz - Myz * Myz;
434 A0 = Mxz * (Mxz * Myy - Myz * Mxy) +
435 Myz * (Myz * Mxx - Mxz * Mxy) - Var_z * Cov_xy;
439 for (x=0., y=A0, iter=0; iter<IterMAX; iter++)
441 double Dy = A1 + x * (A22 + A33 *
x);
442 double xnew = x - y /
Dy;
444 double ynew = A0 + xnew * (A1 + xnew * (A2 + xnew * A3));
445 if (fabs(ynew)>=fabs(y))
break;
451 DET = x*x - x*Mz + Cov_xy;
452 Xcenter = (Mxz*(Myy -
x) - Myz*Mxy)/DET/2;
453 Ycenter = (Myz*(Mxx -
x) - Mxz*Mxy)/DET/2;
457 X0 = Xcenter + meanX;
458 Y0 = Ycenter + meanY;
459 R = sqrt(Xcenter*Xcenter + Ycenter*Ycenter + Mz);