ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
PHTruthTrackSeeding.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file PHTruthTrackSeeding.cc
1 #include "PHTruthTrackSeeding.h"
2 
3 #include "AssocInfoContainer.h"
4 
5 #include <trackbase_historic/SvtxTrack.h> // for SvtxTrack, SvtxTra...
6 #include <trackbase_historic/SvtxTrackMap.h> // for SvtxTrackMap, Svtx...
8 
11 
15 
16 #include <g4eval/SvtxClusterEval.h>
17 #include <g4eval/SvtxEvalStack.h>
18 #include <g4eval/SvtxEvaluator.h>
19 #include <g4eval/SvtxTrackEval.h>
20 
21 #include <g4main/PHG4Hit.h> // for PHG4Hit
22 #include <g4main/PHG4Particle.h> // for PHG4Particle
24 #include <g4main/PHG4HitDefs.h> // for keytype
26 #include <g4main/PHG4VtxPoint.h>
27 
29 
30 #include <phool/getClass.h>
31 #include <phool/phool.h>
32 
33 #include <TDatabasePDG.h>
34 #include <TParticlePDG.h>
35 
36 #include <cstdlib> // for exit
37 #include <iostream> // for operator<<, endl
38 #include <map> // for multimap, map<>::c...
39 #include <memory>
40 #include <utility> // for pair
41 #include <cassert>
42 #include <set>
43 
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
47 
48 class PHCompositeNode;
49 
50 using namespace std;
51 
53  : PHTrackSeeding(name)
54  , _clustereval(nullptr)
55 {}
56 
58 {
59  cout << "Enter PHTruthTrackSeeding:: Setup" << endl;
60 
61  int ret = PHTrackSeeding::Setup(topNode);
62  if (ret != Fun4AllReturnCodes::EVENT_OK) return ret;
63 
64  ret = GetNodes(topNode);
65  if (ret != Fun4AllReturnCodes::EVENT_OK) return ret;
66  _clustereval = new SvtxClusterEval(topNode);
67  _clustereval->do_caching(true);
68  // _clustereval.set_strict(strict);
69  // _clustereval.set_verbosity(verbosity);
71 }
72 
74 {
75  _clustereval->next_event(topNode);
76  _track_map->Clear();
77 
78  typedef std::map<int, std::set<TrkrCluster*> > TrkClustersMap;
79  TrkClustersMap m_trackID_clusters;
80 
81  vector<TrkrDefs::cluskey> ClusterKeyList;
82 
84  // Float_t gntracks = (Float_t) truthinfo->GetNumPrimaryVertexParticles();
85  for (PHG4TruthInfoContainer::ConstIterator iter = range.first;
86  iter != range.second;
87  ++iter){
88  ClusterKeyList.clear();
89  PHG4Particle* g4particle = iter->second;
90 
91  if (g4particle==NULL){
92  cout <<__PRETTY_FUNCTION__<<" - validity check failed: missing truth particle" <<endl;
93  exit(1);
94  }
95 
96  const float gtrackID = g4particle->get_track_id();
97  const int vertexID = g4particle->get_vtx_id();
98  // monentum cut-off
99  if (_min_momentum>0){
100  const double monentum2 =
101  g4particle->get_px() * g4particle->get_px()+
102  g4particle->get_py() * g4particle->get_py()+
103  g4particle->get_pz() * g4particle->get_pz();
104 
105  if (monentum2 < _min_momentum * _min_momentum){
106  if (Verbosity() >= 3){
107  cout <<__PRETTY_FUNCTION__<<" ignore low momentum particle "<< gtrackID <<endl;
108  g4particle->identify();
109  }
110  continue;
111  }
112  }
113 
114  for(unsigned int layer = _min_layer;layer < _max_layer;layer++){
116  if(cluskey!=0)
117  ClusterKeyList.push_back(cluskey);
118  }
119  if(ClusterKeyList.size()< _min_clusters_per_track)
120  continue;
121 
122  auto svtx_track = std::make_unique<SvtxTrack_FastSim_v2>();
123  svtx_track->set_id(_track_map->size());
124  svtx_track->set_truth_track_id(gtrackID);
126  // This is te truth particle vertex, and does not necessarily give the ID in the vertex map
127  svtx_track->set_vertex_id(vertexID-1);
128 
129  // set track charge
130  /*
131  * having the proper track charge is necessary for the ACTS fit to work properly
132  * with normal tracking, it is set by the track seeding. Here we get it from the G4Particle
133  * unfortunately, particle charge is not stored into PHG4Particle.
134  * we need to retrieve it from TParticlePDG instead
135  */
136  {
137  const auto particle = TDatabasePDG::Instance()->GetParticle(g4particle->get_pid());
138  if( particle ) svtx_track->set_charge(particle->Charge());
139  }
140 
141  // Smear the truth values out by 5% so that the seed momentum and
142  // position aren't completely exact
143 
144  double random = 0;//((double) rand() / (RAND_MAX)) * 0.05;
145  // make it negative sometimes
146  if(rand() % 2)
147  random *= -1;
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));
151 
152  // assign the track position using the truth vertex for this track
153  auto g4vertex = _g4truth_container->GetVtx(vertexID);
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));
157 
158  for (TrkrDefs::cluskey cluskey : ClusterKeyList){
159  svtx_track->insert_cluster_key(cluskey);
160  _assoc_container->SetClusterTrackAssoc(cluskey,svtx_track->get_id());
161  }
162 
164  {
165  double x, y, z, px, py, pz;
166  circleFitSeed(ClusterKeyList, x, y, z,
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);
174  }
175 
176  svtx_track->set_ndf(ClusterKeyList.size()*3-5);
177  svtx_track->set_chisq(1.5*ClusterKeyList.size()*3-5);
178  _track_map->insert(svtx_track.get());
179  }
180 
181  if (Verbosity() >= 5)
182  {
183  cout << "Loop over TrackMap " << _track_map_name << " entries " << endl;
184  for (SvtxTrackMap::Iter iter = _track_map->begin();
185  iter != _track_map->end(); ++iter)
186  {
187  SvtxTrack* svtx_track = iter->second;
188 
189  //svtx_track->identify();
190  //continue;
191 
192  cout << "Track ID: " << svtx_track->get_id() << ", Dummy Track pT: "
193  << svtx_track->get_pt() << ", Truth Track/Particle ID: "
194  << svtx_track->get_truth_track_id()
195  << " (X,Y,Z) " << svtx_track->get_x() << ", " << svtx_track->get_y() << ", " << svtx_track->get_z()
196  << endl;
197  cout << " nhits: " << svtx_track->size_cluster_keys()<< endl;
198  //Print associated clusters;
199  for (SvtxTrack::ConstClusterKeyIter iter_clus =
200  svtx_track->begin_cluster_keys();
201  iter_clus != svtx_track->end_cluster_keys(); ++iter_clus)
202  {
203  TrkrDefs::cluskey cluster_key = *iter_clus;
204  cout << "Key: " << cluster_key<< endl;
205  TrkrCluster* cluster = _cluster_map->findCluster(cluster_key);
206  float radius = sqrt(
207  cluster->getX() * cluster->getX() + cluster->getY() * cluster->getY());
208  cout << " cluster ID: "
209  << cluster->getClusKey() << ", cluster radius: " << radius
210  << endl;
211  }
212  }
213  }
214 
215  //==================================
216 
218 }
219 
221 {
222 
223  m_clusterMap = findNode::getClass<TrkrClusterContainer>(topNode,
224  "TRKR_CLUSTER");
225 
226  if(!m_clusterMap)
227  {
228  cerr << PHWHERE << "Error: Can't find node TRKR_CLUSTER" << endl;
230  }
231 
232  _g4truth_container = findNode::getClass<PHG4TruthInfoContainer>(topNode, "G4TruthInfo");
233  if (!_g4truth_container)
234  {
235  cerr << PHWHERE << " ERROR: Can't find node G4TruthInfo" << endl;
237  }
238 
239  hittruthassoc = findNode::getClass<TrkrHitTruthAssoc>(topNode, "TRKR_HITTRUTHASSOC");
240  if (!hittruthassoc)
241  {
242  cout << PHWHERE << "Failed to find TRKR_HITTRUTHASSOC node, quit!" << endl;
243  exit(1);
244  }
245 
246  using nodePair = std::pair<std::string, PHG4HitContainer*&>;
247  std::initializer_list<nodePair> nodes =
248  {
249  { "G4HIT_TPC", phg4hits_tpc },
250  { "G4HIT_INTT", phg4hits_intt },
251  { "G4HIT_MVTX", phg4hits_mvtx },
252  { "G4HIT_MICROMEGAS", phg4hits_micromegas }
253  };
254 
255  for( auto&& node: nodes )
256  {
257  if( !( node.second = findNode::getClass<PHG4HitContainer>( topNode, node.first ) ) )
258  { std::cerr << PHWHERE << " PHG4HitContainer " << node.first << " not found" << std::endl; }
259  }
260 
262 }
263 
265 {
266  return 0;
267 }
268 
269 
270 
271 void PHTruthTrackSeeding::circleFitSeed(std::vector<TrkrDefs::cluskey> clusters,
272  double& x, double& y, double& z,
273  double& px, double& py, double& pz,
274  int charge)
275 {
276  double R, X0, Y0;
277  circleFitByTaubin(clusters, R, X0, Y0);
278 
279  findRoot(R, X0, Y0, x, y);
280 
281  double phi = atan2(-1 * (X0-x), Y0-y);
282 
283  if(charge > 0)
284  {
285  phi += M_PI;
286  if(phi > M_PI)
287  phi -= 2. * M_PI;
288  }
289 
290  double m, B;
291  lineFit(clusters, m, B);
292  z = B;
293 
294  double theta = atan(1./m);
295  if(theta < 0)
296  theta =+ M_PI;
297 
300  float pt = 0.3 * 1.4 * R / 100.;
301  float eta = -log(tan(theta/2.));
302  float p = pt * cosh(eta);
303 
306  px = p * sin(theta) * cos(phi);
307  py = p * sin(theta) * sin(phi);
308  pz = p * cos(theta);
309 
310 }
311 
312 void PHTruthTrackSeeding::lineFit(std::vector<TrkrDefs::cluskey>& clusters,
313  double& A, double& B)
314 {
315 
316  double xsum = 0,x2sum = 0,ysum = 0,xysum = 0;
317  for(auto& clusterkey : clusters)
318  {
319  auto cluster = m_clusterMap->findCluster(clusterkey);
320  double z = cluster->getZ();
321  double r = sqrt(pow(cluster->getX(),2) + pow(cluster->getY(), 2));
322 
323  xsum=xsum+r; // calculate sigma(xi)
324  ysum=ysum+z; // calculate sigma(yi)
325  x2sum=x2sum+pow(r,2); // calculate sigma(x^2i)
326  xysum=xysum+r*z; // calculate sigma(xi*yi)
327  }
328 
330  A = (clusters.size()*xysum-xsum*ysum) / (clusters.size()*x2sum-xsum*xsum);
331 
333  B = (x2sum*ysum-xsum*xysum) / (x2sum*clusters.size()-xsum*xsum);
334 
335 }
336 void PHTruthTrackSeeding::findRoot(const double& R, const double& X0, const double& Y0, double& x, double& y)
337 {
338 
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));
342 
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));
346 
347  double minx = sqrt(pow(R, 2) - pow(miny - Y0, 2)) + X0;
348  double minx2 = -sqrt(pow(R, 2) - pow(miny2 - Y0, 2)) + X0;
349 
351  if(fabs(minx) < fabs(minx2))
352  x = minx;
353  else
354  x = minx2;
355 
356  if(fabs(miny) < fabs(miny2))
357  y = miny;
358  else
359  y = miny2;
360 
361 }
362 
363 void PHTruthTrackSeeding::circleFitByTaubin(std::vector<TrkrDefs::cluskey>& clusters,
364  double& R, double& X0, double& Y0)
365 {
381  int iter, IterMAX=99;
382 
383  double Mz, Mxy, Mxx, Myy, Mxz, Myz, Mzz, Cov_xy, Var_z;
384  double A0, A1, A2, A22, A3, A33;
385  double x, y;
386  double DET, Xcenter, Ycenter;
387 
388  // Compute x- and y- sample means
389  double meanX = 0;
390  double meanY = 0;
391  double weight = 0;
392 
393  for(auto cluskey : clusters)
394  {
395  auto clus = m_clusterMap->findCluster(cluskey);
396  meanX += clus->getX();
397  meanY += clus->getY();
398  weight++;
399  }
400  meanX /= weight;
401  meanY /= weight;
402 
403  Mxx=Myy=Mxy=Mxz=Myz=Mzz=0.;
404 
405  for(auto cluskey : clusters)
406  {
407  auto clus = m_clusterMap->findCluster(cluskey);
408 
409  double Xi = clus->getX() - meanX;
410  double Yi = clus->getY() - meanY;
411  double Zi = Xi * Xi + Yi * Yi;
412 
413  Mxy += Xi*Yi;
414  Mxx += Xi*Xi;
415  Myy += Yi*Yi;
416  Mxz += Xi*Zi;
417  Myz += Yi*Zi;
418  Mzz += Zi*Zi;
419  }
420 
421  Mxx /= weight;
422  Myy /= weight;
423  Mxy /= weight;
424  Mxz /= weight;
425  Myz /= weight;
426  Mzz /= weight;
427 
428  Mz = Mxx + Myy;
429  Cov_xy = Mxx * Myy - Mxy * Mxy;
430  Var_z = Mzz - Mz * Mz;
431  A3 = 4 * 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;
436  A22 = A2 + A2;
437  A33 = A3 + A3 + A3;
438 
439  for (x=0., y=A0, iter=0; iter<IterMAX; iter++) // usually, 4-6 iterations are enough
440  {
441  double Dy = A1 + x * (A22 + A33 * x);
442  double xnew = x - y / Dy;
443  if ((xnew == x)||(!std::isfinite(xnew))) break;
444  double ynew = A0 + xnew * (A1 + xnew * (A2 + xnew * A3));
445  if (fabs(ynew)>=fabs(y)) break;
446  x = xnew; y = ynew;
447  }
448 
449  // computing parameters of the fitting circle
450 
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;
454 
455  // assembling the output
456 
457  X0 = Xcenter + meanX;
458  Y0 = Ycenter + meanY;
459  R = sqrt(Xcenter*Xcenter + Ycenter*Ycenter + Mz);
460 
461 }