ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
PHSiliconTpcTrackMatching.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file PHSiliconTpcTrackMatching.cc
2 
3 #include "AssocInfoContainer.h"
4 
6 #include <trackbase/TrkrDefs.h> // for cluskey, getTrkrId, tpcId
12 #include <trackbase_historic/SvtxVertex.h> // for SvtxVertex
15 
16 #include <g4main/PHG4Hit.h> // for PHG4Hit
17 #include <g4main/PHG4Particle.h> // for PHG4Particle
18 #include <g4main/PHG4HitDefs.h> // for keytype
19 
21 
22 #include <phool/phool.h>
23 #include <phool/PHCompositeNode.h>
24 #include <phool/getClass.h>
25 
26 
27 #if __cplusplus < 201402L
28 #include <boost/make_unique.hpp>
29 #endif
30 
31 #include <TF1.h>
32 
33 #include <climits> // for UINT_MAX
34 #include <iostream> // for operator<<, basic_ostream
35 #include <cmath> // for fabs, sqrt
36 #include <set> // for _Rb_tree_const_iterator
37 #include <utility> // for pair
38 #include <memory>
39 
40 using namespace std;
41 
42 //____________________________________________________________________________..
44  SubsysReco(name)
45  , _track_map_name_silicon("SvtxSiliconTrackMap")
46 {
47  //cout << "PHSiliconTpcTrackMatching::PHSiliconTpcTrackMatching(const std::string &name) Calling ctor" << endl;
48 }
49 
50 //____________________________________________________________________________..
52 {
53 
54 }
55 
56 //____________________________________________________________________________..
58 {
59  // put these in the output file
60  cout << PHWHERE "_is_ca_seeder " << _is_ca_seeder << " Search windows: phi " << _phi_search_win << " eta "
61  << _eta_search_win << " pp_mode " << _pp_mode << endl;
62 
63  // corrects the PHTpcTracker phi bias
64  fdphi = new TF1("f1", "[0] + [1]/x^[2]");
65  fdphi->SetParameter(0, _par0);
66  fdphi->SetParameter(1, _par1);
67  fdphi->SetParameter(2, _par2);
68 
69  // corrects the space charge distortion phi bias
70  if(!_is_ca_seeder)
71  {
72  // PHTpcTracker correction is opposite in sign
73  // and different in magnitude - why?
74  _parsc0 *= -1.0 * 0.7;
75  _parsc1 *= -1.0 * 0.7;
76  }
77  fscdphi = new TF1("f2","[0] + [1]*x^2");
80 
81  int ret = GetNodes(topNode);
82  if (ret != Fun4AllReturnCodes::EVENT_OK) return ret;
83 
84  return ret;
85 }
86 
87 //____________________________________________________________________________..
89 {
90  // _track_map contains the TPC seed track stubs
91  // _track_map_silicon contains the silicon seed track stubs
92  // We will add the silicon clusters to the TPC tracks already on the node tree
93  // We will have to expand the number of tracks whenever we find multiple matches to the silicon
94 
96  _z_mismatch_map.clear();
97 
98  if(Verbosity() > 0)
99  cout << PHWHERE << " TPC track map size " << _track_map->size() << " Silicon track map size " << _track_map_silicon->size() << endl;
100 
101  if(_track_map->size() == 0)
103 
104  if(Verbosity() > 2)
105  {
106  // list silicon tracks
107  for (auto phtrk_iter_si = _track_map_silicon->begin();
108  phtrk_iter_si != _track_map_silicon->end();
109  ++phtrk_iter_si)
110  {
111  _tracklet_si = phtrk_iter_si->second;
112 
113  double si_phi = atan2(_tracklet_si->get_py(), _tracklet_si->get_px());
114  double si_eta = _tracklet_si->get_eta();
115 
116  cout << " Si track " << _tracklet_si->get_id() << " si_phi " << si_phi << " si_eta " << si_eta << endl;
117  }
118  }
119 
120  // Find all matches of tpc and si tracklets in eta and phi, x and y
121  // If _pp_mode is not set, a match in z is also required - gives same behavior as old code
122  // In any case, multiple matches are handled by duplicating the tpc tracklet into a new track
123  // "_seed_track_map" records (original id,duplicate id) so that the track cleaner can choose the one with the best Acts fit later
124  std::multimap<unsigned int, unsigned int> tpc_matches;
125  std::set<unsigned int> tpc_matched_set;
126  findEtaPhiMatches(tpc_matched_set, tpc_matches);
127 
128  // We have a complete list of all eta/phi matched tracks in the map "tpc_matches"
129  // In all cases, the crossing value is only rough at this point (it is only a tag)
130  std::multimap<int, std::pair<unsigned int, unsigned int>> crossing_matches;
131  std::map<unsigned int, int> tpc_crossing_map;
132  std::set<int> crossing_set;
133  if(_pp_mode)
134  {
135  // This section is to correct the TPC z positions of tracks for all bunch crossings.
136  //=================================================================================
137 
138  // All tracks treated as if we do not know the bunch crossing
139  // The crossing estimate here is crude, for now
140  tagMatchCrossing(tpc_matches, crossing_set, crossing_matches, tpc_crossing_map);
141 
142  // Sort candidates by the silicon tracklet Z position by putting them in si_sorted map
143  // -- captures crossing, tpc_id and si_id
144  std::multimap<double, std::pair<unsigned int, unsigned int>> si_sorted_map;
145  for( auto ncross : crossing_set)
146  {
147  if(Verbosity() > 1) std::cout << " ncross = " << ncross << std::endl;
148 
149  auto ret = crossing_matches. equal_range(ncross);
150  for(auto it = ret.first; it != ret.second; ++it)
151  {
152  if(Verbosity() > 1)
153  std::cout << " crossing " << it->first << " tpc_id " << it->second.first << " si_id " << it->second.second
154  << " pT " << _track_map->get(it->second.first)->get_pt()
155  << " eta " << _track_map->get(it->second.first)->get_eta()
156  << " tpc_z " << _track_map->get(it->second.first)->get_z()
157  << " si_z " << _track_map_silicon->get(it->second.second)->get_z()
158  << std::endl;
159 
160  double z_si = _track_map_silicon->get(it->second.second)->get_z();
161  si_sorted_map.insert(std::make_pair(z_si, std::make_pair(it->second.first, it->second.second)));
162  }
163  }
164 
165  // make a list of silicon vertices, and a multimap with all associated tpc/si pairs
166  std::multimap<unsigned int, std::pair<unsigned int, unsigned int>> vertex_map;
167  std::vector<double> vertex_list;
168  getSiVertexList(si_sorted_map, vertex_list, vertex_map);
169 
170  // find the crossing number for every track from the mismatch with the associated silicon vertex z position
171  std::map<unsigned int, double> vertex_crossings_map;
172  getCrossingNumber(vertex_list, vertex_map, vertex_crossings_map);
173 
174  // Remove tracks from vertex_map where the vertex crossing is badly inconsistent with the initial crossing estimate
175  cleanVertexMap( vertex_crossings_map, vertex_map, tpc_crossing_map );
176 
177  // correct the TPC cluster z values for the bunch crossing offset
178  correctTpcClusterZ(vertex_crossings_map, vertex_map);
179 
180  // add silicon clusters to the surviving tracks
181  addSiliconClusters(vertex_map);
182 
183  //===================================================
184  }
185  else
186  {
187  // only crossing zero has been added to the tpc_matches map, just add silicon clusters
188  addSiliconClusters(tpc_matches);
189  }
190 
191  // loop over all tracks and copy the silicon clusters to the corrected cluster map
194 
195  if(Verbosity() > 0)
196  cout << " Final track map size " << _track_map->size()
197  << " seed-track map size " << _seed_track_map->size() << endl;
198 
199  if (Verbosity() > 0)
200  cout << "PHSiliconTpcTrackMatching::process_event(PHCompositeNode *topNode) Leaving process_event" << endl;
201 
203  }
204 
206 {
208 }
209 
211 {
212  //---------------------------------
213  // Get additional objects off the Node Tree
214  //---------------------------------
215 
216  _assoc_container = findNode::getClass<AssocInfoContainer>(topNode, "AssocInfoContainer");
217  if (!_assoc_container)
218  {
219  cerr << PHWHERE << " ERROR: Can't find AssocInfoContainer." << endl;
221  }
222 
223  _track_map_silicon = findNode::getClass<SvtxTrackMap>(topNode, _silicon_track_map_name);
224  if (!_track_map_silicon)
225  {
226  cerr << PHWHERE << " ERROR: Can't find SvtxSiliconTrackMap: " << endl;
228  }
229 
230  _track_map = findNode::getClass<SvtxTrackMap>(topNode, _track_map_name);
231  if (!_track_map)
232  {
233  cerr << PHWHERE << " ERROR: Can't find SvtxTrackMap " << endl;
235  }
236 
237  _seed_track_map = findNode::getClass<TpcSeedTrackMap>(topNode, _tpcseed_track_map_name);
238  if(!_seed_track_map)
239  {
240  std::cout << "Creating node TpcSeedTrackMap" << std::endl;
241 
243  PHNodeIterator iter(topNode);
244  PHCompositeNode *dstNode = dynamic_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", "DST"));
245 
247  if (!dstNode)
248  {
249  std::cerr << "DST Node missing, quitting" << std::endl;
250  throw std::runtime_error("failed to find DST node in PHActsSourceLinks::createNodes");
251  }
252 
254  PHCompositeNode *svtxNode = dynamic_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", "SVTX"));
255 
257  if (!svtxNode)
258  {
259  svtxNode = new PHCompositeNode("SVTX");
260  dstNode->addNode(svtxNode);
261  }
262 
266  svtxNode->addNode(node);
267  }
268 
269  _corrected_cluster_map = findNode::getClass<TrkrClusterContainer>(topNode,"CORRECTED_TRKR_CLUSTER");
271  {
272  std::cout << " Found CORRECTED_TRKR_CLUSTER node " << std::endl;
273  }
274 
275  _cluster_map = findNode::getClass<TrkrClusterContainer>(topNode, "TRKR_CLUSTER");
276  if (!_cluster_map)
277  {
278  std::cout << PHWHERE << " ERROR: Can't find node TRKR_CLUSTER" << std::endl;
280  }
281 
282  _tGeometry = findNode::getClass<ActsTrackingGeometry>(topNode,"ActsTrackingGeometry");
283  if(!_tGeometry)
284  {
285  std::cout << PHWHERE << "Error, can't find acts tracking geometry" << std::endl;
287  }
288 
289  _surfmaps = findNode::getClass<ActsSurfaceMaps>(topNode,"ActsSurfaceMaps");
290  if(!_surfmaps)
291  {
292  std::cout << PHWHERE << "Error, can't find acts surface maps" << std::endl;
294  }
295 
297 }
298 
299 double PHSiliconTpcTrackMatching::getBunchCrossing(unsigned int trid, double z_mismatch )
300 {
301  // the bunch crossing separation is 110 nsec
302  double vdrift = 8.00; // cm /microsecond
303  double z_bunch_separation = 0.107 * vdrift; // 107 ns bunch crossing interval
304 
305  // The sign of z_mismatch will depend on which side of the TPC the tracklet is in
306  SvtxTrack *track = _track_map->get(trid);
307 
308  double crossings = z_mismatch / z_bunch_separation;
309 
310  // Check the sign of z for the first cluster in the track
311  // loop over associated clusters to get hits for TPC only, add to new track copy
312  double clus_z = 0.0;
313  ActsTransformations transformer;
315  iter != track->end_cluster_keys();
316  ++iter)
317  {
318  TrkrDefs::cluskey cluster_key = *iter;
319  unsigned int trkrid = TrkrDefs::getTrkrId(cluster_key);
320  if(trkrid == TrkrDefs::tpcId)
321  {
322 
323  TrkrCluster *tpc_clus;
325  tpc_clus = _corrected_cluster_map->findCluster(cluster_key);
326  else
327  tpc_clus = _cluster_map->findCluster(cluster_key);
328 
329  auto global = transformer.getGlobalPosition(tpc_clus,
330  _surfmaps,
331  _tGeometry);
332  clus_z = global[2];
333  break; // we only need the first one
334  }
335  }
336 
337  // if true z > 0, the z offset will be negative from a positive t0, so z_tpc - z_si will be negative for a positive crossing offset
338  // -- reverse the sign of crossings to get a positive crossing for a positive t0
339  // if true z < 0, the z offset will be positive from a positive t0, z_tpc - z_si will be positive for a positive crossing offset
340  // -- the sign of crossings is OK
341  // have to correct clus_z for the z_mismatch to see if it was really positive or negative
342  if(clus_z -z_mismatch > 0)
343  crossings *= -1.0;
344 
345  if(Verbosity() > 1) std::cout << " trackid " << trid << " clus_z " << clus_z << " z_mismatch " << z_mismatch << " crossings " << crossings << std::endl;
346 
347  return crossings;
348 }
349 
350 double PHSiliconTpcTrackMatching::getMedian(std::vector<double> &v)
351 {
352  if(v.size() == 0) return NAN;
353 
354  double median = 0.0;
355 
356  if( (v.size() % 2) == 0)
357  {
358  // even number of entries
359  // we want the average of the middle two numbers, v.size()/2 and v.size()/2-1
360  auto m1 = v.begin() + v.size()/2;
361  std::nth_element(v.begin(), m1, v.end());
362  double median1 = v[v.size()/2];
363 
364  auto m2 = v.begin() + v.size()/2 - 1;
365  std::nth_element(v.begin(), m2, v.end());
366  double median2 = v[v.size()/2 - 1];
367 
368  median = (median1 + median2) / 2.0;
369  if(Verbosity() > 2) std::cout << "The vector size is " << v.size()
370  << " element m is " << v.size() / 2 << " = " << v[v.size()/2]
371  << " element m-1 is " << v.size() / 2 -1 << " = " << v[v.size()/2-1]
372  << std::endl;
373  }
374  else
375  {
376  // odd number of entries
377  auto m = v.begin() + v.size()/2;
378  std::nth_element(v.begin(), m, v.end());
379  median = v[v.size()/2];
380  if(Verbosity() > 2) std::cout << "The vector size is " << v.size() << " element m is " << v.size() / 2 << " = " << v[v.size()/2] << std::endl;
381  }
382 
383  return median ;
384 }
385 
386 void PHSiliconTpcTrackMatching::addSiliconClusters( std::multimap<unsigned int, unsigned int> &tpc_matches )
387 {
388 
389  for(auto it = tpc_matches.begin(); it != tpc_matches.end(); ++it)
390  {
391  unsigned int tpcid = it->first;
392  SvtxTrack *tpc_track = _track_map->get(tpcid);
393  if(Verbosity() > 1) std::cout << " tpcid " << tpcid << " original z " << tpc_track->get_z() << std::endl;
394 
395  // add the silicon cluster keys to the track
396  unsigned int si_id = it->second;
397  SvtxTrack *si_track = _track_map_silicon->get(si_id);
398  if(Verbosity() > 1) std::cout << " si track id " << si_id << std::endl;
399  for (SvtxTrack::ConstClusterKeyIter si_iter = si_track->begin_cluster_keys();
400  si_iter != si_track->end_cluster_keys();
401  ++si_iter)
402  {
403  TrkrDefs::cluskey si_cluster_key = *si_iter;
404 
405  if(Verbosity() > 1)
406  cout << " inserting si cluster key " << si_cluster_key << " into existing TPC track " << tpc_track->get_id() << endl;
407 
408  tpc_track->insert_cluster_key(si_cluster_key);
409  _assoc_container->SetClusterTrackAssoc(si_cluster_key, tpc_track->get_id());
410 
411 
412  }
413 
414  // update the track position to the si one
415  tpc_track->set_x(si_track->get_x());
416  tpc_track->set_y(si_track->get_y());
417  tpc_track->set_z(si_track->get_z());
418 
419  if(Verbosity() > 2)
420  std::cout << " TPC seed track ID " << tpc_track->get_id() << " si track id " << si_track->get_id()
421  << " new nclus " << tpc_track->size_cluster_keys() << std::endl;
422 
423  if(Verbosity() > 2)
424  tpc_track->identify();
425  }
426 
427  return;
428 }
429 
430 // not used
431 void PHSiliconTpcTrackMatching::addSiliconClusters( std::multimap<int, std::pair<unsigned int, unsigned int>> &crossing_matches )
432 {
433 
434  for(auto it = crossing_matches.begin(); it != crossing_matches.end(); ++it)
435  {
436  unsigned int tpcid = it->second.first;
437  SvtxTrack *tpc_track = _track_map->get(tpcid);
438  if(Verbosity() > 1) std::cout << " tpcid " << tpcid << " original z " << tpc_track->get_z() << std::endl;
439 
440  // add the silicon cluster keys to the track
441  unsigned int si_id = it->second.second;
442  SvtxTrack *si_track = _track_map_silicon->get(si_id);
443  if(Verbosity() > 1) std::cout << " si track id " << si_id << std::endl;
444  for (SvtxTrack::ConstClusterKeyIter si_iter = si_track->begin_cluster_keys();
445  si_iter != si_track->end_cluster_keys();
446  ++si_iter)
447  {
448  TrkrDefs::cluskey si_cluster_key = *si_iter;
449 
450  if(Verbosity() > 1)
451  cout << " inserting si cluster key " << si_cluster_key << " into existing TPC track " << tpc_track->get_id() << endl;
452 
453  tpc_track->insert_cluster_key(si_cluster_key);
454  _assoc_container->SetClusterTrackAssoc(si_cluster_key, tpc_track->get_id());
455 
456  }
457 
458  // update the track position to the si one
459  tpc_track->set_x(si_track->get_x());
460  tpc_track->set_y(si_track->get_y());
461  tpc_track->set_z(si_track->get_z());
462 
463  if(Verbosity() > 2)
464  std::cout << " TPC seed track ID " << tpc_track->get_id() << " si track id " << si_track->get_id()
465  << " new nclus " << tpc_track->size_cluster_keys() << std::endl;
466 
467  if(Verbosity() > 2)
468  tpc_track->identify();
469  }
470 
471  return;
472 }
473 
474 void PHSiliconTpcTrackMatching::addSiliconClusters( std::multimap<unsigned int, std::pair<unsigned int, unsigned int>> &vertex_map)
475 {
476 
477  for(auto it = vertex_map.begin(); it != vertex_map.end(); ++it)
478  {
479  unsigned int tpcid = it->second.first;
480  SvtxTrack *tpc_track = _track_map->get(tpcid);
481  if(Verbosity() > 1) std::cout << " tpcid " << tpcid << " original z " << tpc_track->get_z() << std::endl;
482 
483  // add the silicon cluster keys to the track
484  unsigned int si_id = it->second.second;
485  SvtxTrack *si_track = _track_map_silicon->get(si_id);
486  if(Verbosity() > 1) std::cout << " si track id " << si_id << std::endl;
487  for (SvtxTrack::ConstClusterKeyIter si_iter = si_track->begin_cluster_keys();
488  si_iter != si_track->end_cluster_keys();
489  ++si_iter)
490  {
491  TrkrDefs::cluskey si_cluster_key = *si_iter;
492 
493  if(Verbosity() > 1)
494  cout << " inserting si cluster key " << si_cluster_key << " into existing TPC track " << tpc_track->get_id() << endl;
495 
496  tpc_track->insert_cluster_key(si_cluster_key);
497  _assoc_container->SetClusterTrackAssoc(si_cluster_key, tpc_track->get_id());
498 
499  }
500 
501  // update the track position to the si one
502  tpc_track->set_x(si_track->get_x());
503  tpc_track->set_y(si_track->get_y());
504  tpc_track->set_z(si_track->get_z());
505 
506  if(Verbosity() > 2)
507  {
508  std::cout << " TPC seed track ID " << tpc_track->get_id() << " si track id " << si_track->get_id()
509  << " new nclus " << tpc_track->size_cluster_keys() << std::endl;
510  }
511 
512  if(Verbosity() > 2)
513  tpc_track->identify();
514  }
515 
516  return;
517 }
518 
520  std::map<unsigned int, double> &vertex_crossings_map,
521  std::multimap<unsigned int, std::pair<unsigned int, unsigned int>> &vertex_map )
522 {
523  if(vertex_crossings_map.size() == 0) return;
524 
525  // Correct the z positions of the clusters
526  double vdrift = 8.00; // cm /microsecond
527  double z_bunch_separation = 0.107 * vdrift; // 107 ns bunch crossing interval
528  ActsTransformations transformer;
529  for(auto [ivert, crossing] : vertex_crossings_map)
530  {
531  // the direction of the z shift depends on
532  // a) the sign of crossings (ie. t0 earlier or later), b) the side of the TPC
533  // positive (corrected) z values: z values decrease for positive crossings, correction is positive = crossings * z_bunch_separation
534  // negative z values: z values increase for positive crossings, correction is negative = - crossings * z_bunch_separation
535 
536  if(Verbosity() > 1) std::cout << " ivert " << ivert << " crossing " << crossing << std::endl;
537 
538  if(crossing == 0) continue;
539 
540  // loop over all tracks associated with this vertex
541  // remember that any tracks not associated with one of the vertices are unusable
542  auto ret = vertex_map.equal_range(ivert);
543  for(auto it = ret.first; it != ret.second; ++it)
544  {
545  unsigned int tpcid = it->second.first;
546  SvtxTrack *tpc_track = _track_map->get(tpcid);
547  if(Verbosity() > 1) std::cout << " tpcid " << tpcid << " original z " << tpc_track->get_z() << std::endl;
548 
549  // loop over associated clusters to get hits for TPC only, shift the z positions by z_shift
550  for (SvtxTrack::ConstClusterKeyIter iter = tpc_track->begin_cluster_keys();
551  iter != tpc_track->end_cluster_keys();
552  ++iter)
553  {
554  TrkrDefs::cluskey cluster_key = *iter;
555  unsigned int trkrid = TrkrDefs::getTrkrId(cluster_key);
556  if(trkrid == TrkrDefs::tpcId)
557  {
558  // get the cluster z
559  TrkrCluster *tpc_clus;
561  tpc_clus = _corrected_cluster_map->findCluster(cluster_key);
562  else
563  tpc_clus = _cluster_map->findCluster(cluster_key);
564 
565  if(Verbosity() > 2) std::cout << " original local cluster z " << tpc_clus->getLocalY() << std::endl;
566  auto global = transformer.getGlobalPosition(tpc_clus,
567  _surfmaps,
568  _tGeometry);
569  double clus_z = global[2];
570  if(Verbosity() > 2) std::cout << " global: x " << global[0] << " y " << global[1] << " z " << global[2] << std::endl;
571 
572  double corrected_clus_z;
573  if(clus_z > 0)
574  corrected_clus_z = clus_z + crossing * z_bunch_separation;
575  else
576  corrected_clus_z = clus_z - crossing * z_bunch_separation;
577 
578  auto surface = transformer.getSurface(tpc_clus, _surfmaps);
580  double surfZCenter = center[2];
581  double local_z = corrected_clus_z - surfZCenter;
582  tpc_clus->setLocalY(local_z);
583  if(Verbosity() > 2) std::cout << " original clus_z " << clus_z << " corrected_clus_z " << corrected_clus_z << " local z " << local_z << std::endl;
584 
585 
586  if(Verbosity() > 2)
587  {
588  TrkrCluster *tpc_clus_check;
590  tpc_clus_check = _corrected_cluster_map->findCluster(cluster_key);
591  else
592  tpc_clus_check = _cluster_map->findCluster(cluster_key);
593 
594  auto global_check = transformer.getGlobalPosition(tpc_clus_check,
595  _surfmaps,
596  _tGeometry);
597  std::cout << " global check: x " << global_check[0] << " y " << global_check[1] << " z " << global_check[2] << std::endl;
598  }
599  }
600  }
601  }
602  }
603  return;
604 }
605 
607  std::map<unsigned int, double> &vertex_crossings_map,
608  std::multimap<unsigned int, std::pair<unsigned int, unsigned int>> &vertex_map,
609  std::map<unsigned int, int> &tpc_crossing_map )
610 {
611  if(vertex_crossings_map.size() == 0) return;
612 
613  // Sanity check: is the initial rough crossing estimate consistent with the crossing number for this si vertex?
614 
615  std::multimap<unsigned int, std::pair<unsigned int, unsigned int>> bad_map;
616  for(auto [ivert, crossing] : vertex_crossings_map)
617  {
618  if(Verbosity() > 1) std::cout << " CleanVertexMap: ivert " << ivert << " crossing " << crossing << std::endl;
619 
620  // loop over all tracks associated with this vertex
621  // remember that any tracks not associated with one of the vertices are unusable, they are ignored
622  auto ret = vertex_map.equal_range(ivert);
623  for(auto it = ret.first; it != ret.second; ++it)
624  {
625  unsigned int tpcid = it->second.first;
626  unsigned int si_id = it->second.second;
627 
628  // need the initial crossing estimate for comparison - stored in tpc_ crossing_map
629  auto sit = tpc_crossing_map.find(tpcid);
630  if(abs(sit->second - crossing) > 2)
631  {
632  // Fails the sanity check, mark for removal from the vertex map
633  if(Verbosity() > 1)
634  std::cout << " Crossing mismatch: ivert " << ivert << " tpc id " << tpcid << " vert crossing " << crossing << " track crossing " << sit->second << std::endl;
635 
636  bad_map.insert(std::make_pair(ivert, std::make_pair(tpcid, si_id)));
637  }
638  }
639  }
640 
641  // If we delete an entry from vertex map, the TPC track is not associated with a bunch crossing,
642  // -- the cluster z values are not changed, and the silicon clusters are not added.
643  // The track remains in the track map as a TPC-only track, by default assumed to be in bunch crossing zero
644  // If multiple entries are deleted from vertex map, the TPC-only track copies are left in the track map, and also in _seed-track_map
645  // The track cleaner will remove all but one of them, based on chisq from the Acts fitter.
646 
647  // remove bad entries from vertex_map so the wrong silicon is not associated
648  for(auto [ivert, id_pair] : bad_map)
649  {
650  unsigned int tpc_id = id_pair.first;
651  unsigned int si_id = id_pair.second;
652 
653  // Have to iterate over vertex_map and examine each pair to find the one matching bad_map
654  // this logic works because we call the equal range on vertex_map for every id_pair
655  // so we only delete one entry per equal range call
656  auto ret = vertex_map.equal_range(ivert);
657  for(auto it = ret.first; it != ret.second; ++it)
658  {
659  if(it->second.first == tpc_id && it->second.second == si_id)
660  {
661  if(Verbosity() > 1) std::cout << " erasing entry for ivert " << ivert << " tpc_id " << tpc_id << " si_id " << si_id << std::endl;
662  vertex_map.erase(it);
663  break; // the iterator is no longer valid
664  }
665  }
666  }
667 
668  return;
669 }
670 
672  std::vector<double> &vertex_list,
673  std::multimap<unsigned int, std::pair<unsigned int, unsigned int>> &vertex_map,
674  std::map<unsigned int, double> &vertex_crossings_map)
675 {
676  if(vertex_list.size() == 0) return;
677 
678  for(unsigned int ivert = 0; ivert<vertex_list.size(); ++ivert)
679  {
680  std::vector<double> crossing_vec;
681 
682  double vert_z = vertex_list[ivert];
683  if(Verbosity() > 1) std::cout << "Vertex " << ivert << " vertex z " << vert_z << std::endl;
684 
685  auto ret = vertex_map.equal_range(ivert);
686  for(auto it = ret.first; it != ret.second; ++it)
687  {
688  unsigned int tpcid = it->second.first;
689  double tpc_z = _track_map->get(tpcid)->get_z();
690  double z_mismatch = tpc_z - vert_z;
691 
692  double crossings = getBunchCrossing(tpcid, z_mismatch);
693  crossing_vec.push_back(crossings);
694 
695  if(Verbosity() > 1) std::cout << " tpc_id " << tpcid << " tpc pT " << _track_map->get(tpcid)->get_pt()
696  << " tpc eta " << _track_map->get(tpcid)->get_eta() << " si_id " << it->second.second
697  << " tpc z " << tpc_z << " crossings " << crossings << std::endl;
698  }
699 
700  if(crossing_vec.size() == 0) continue;
701  double crossing_median = getMedian(crossing_vec);
702 
703  double crossing_avge = 0.0;
704  double crossing_wt = 0.0;
705  for(auto cross : crossing_vec)
706  {
707  if(fabs(cross-crossing_median) < 1.0)
708  {
709  crossing_avge += cross;
710  crossing_wt++;
711  }
712  }
713  crossing_avge /= crossing_wt;
714  double crossing_avge_rounded = round(crossing_avge);
715  vertex_crossings_map.insert(std::make_pair(ivert, crossing_avge_rounded));
716 
717  if(Verbosity() > 1) std::cout << " crossing_median " << crossing_median << " crossing average = " << crossing_avge << " crossing_wt " << crossing_wt
718  << " crossing integer " << crossing_avge_rounded << std::endl;
719  }
720  return;
721 }
722 
724  std::multimap<double, std::pair<unsigned int, unsigned int>> &si_sorted_map,
725  std::vector<double> &vertex_list,
726  std::multimap<unsigned int, std::pair<unsigned int, unsigned int>> &vertex_map)
727 {
728  if(si_sorted_map.size() == 0) return;
729 
730  // process si_sorted_map to get the vertex list, and the average silicon z value for each vertex
731  double zkeep = 999;
732  double zavge = 0.0;
733  double zwt = 0.0;
734  unsigned int nvert = 0;
735  for(auto it = si_sorted_map.begin(); it != si_sorted_map.end(); ++it)
736  {
737  double z_si = it->first;
738  std::pair<unsigned int, unsigned int> id_pair = it->second;
739  if(Verbosity() > 1) std::cout << " z_si " << z_si << " tpc_id " << it->second.first << " si_id " << it->second.second << std::endl;
740 
741  if( (zkeep == 999) || (fabs(z_si - zkeep) < _si_vertex_dzmax) )
742  {
743  vertex_map.insert(std::make_pair(nvert, id_pair));
744  zavge+= z_si;
745  zwt++;
746  zkeep = z_si;
747  }
748  else
749  {
750  zavge /= zwt;
751  vertex_list.push_back(zavge);
752 
753  // start a new vertex
754  nvert++;
755  zavge = z_si;
756  zwt = 1.0;
757  zkeep = z_si;
758  vertex_map.insert(std::make_pair(nvert, id_pair));
759  }
760  }
761  // get the last one
762  zavge /= zwt;
763  vertex_list.push_back(zavge);
764 
765  return;
766 }
767 
769  std::set<unsigned int> &tpc_matched_set,
770  std::multimap<unsigned int, unsigned int> &tpc_matches )
771 {
772  // loop over the original TPC tracks
773  for (auto phtrk_iter = _track_map->begin();
774  phtrk_iter != _track_map->end();
775  ++phtrk_iter)
776  {
777  // we may add tracks to the map, so we stop at the last original track
778 
779  _tracklet_tpc = phtrk_iter->second;
780 
781  if (Verbosity() > 1)
782  {
783  std::cout
784  << __LINE__
785  << ": Processing seed itrack: " << phtrk_iter->first
786  << ": nhits: " << _tracklet_tpc-> size_cluster_keys()
787  << ": Total tracks: " << _track_map->size()
788  << ": phi: " << _tracklet_tpc->get_phi()
789  << endl;
790  }
791 
792  // This will always end up as a final track, no matter what - add it to the seed-track-map
793  if(Verbosity() > 1)
794  std::cout << " TPC seed ID " << _tracklet_tpc->get_id() << " original nclus " << _tracklet_tpc->size_cluster_keys() << std::endl;
796 
797  double tpc_phi = atan2(_tracklet_tpc->get_py(), _tracklet_tpc->get_px());
798  double tpc_eta = _tracklet_tpc->get_eta();
799  double tpc_pt = sqrt( pow(_tracklet_tpc->get_px(),2) + pow(_tracklet_tpc->get_py(),2) );
800 
801  // phi correction for PHTpcTracker tracklets is charge dependent
802  double sign_phi_correction = _tracklet_tpc->get_charge();
803 
807  if(_field.find("2d") != std::string::npos)
808  {
809  sign_phi_correction *= -1;
810  if(_fieldDir > 0)
811  sign_phi_correction *= -1;
812  }
813 
814  // this factor will increase the window size at low pT
815  // otherwise the matching efficiency drops off at low pT
816  // it would be better if this was a smooth function
817  double mag = 1.0;
818  if(tpc_pt < 6.0) mag = 2;
819  if(tpc_pt < 3.0) mag = 4.0;
820 
821  if(Verbosity() > 3)
822  {
823  cout << "Original TPC tracklet:" << endl;
825  }
826 
827  // This is no longer applicable I believe, and now wrong - look into it
828  //=============================
829  // correct the TPC tracklet phi for the space charge offset, if this is the calib pass
830  // this is done just to let us tighten up the matching window
831  if(_sc_calib_flag)
832  {
833  tpc_phi -= fscdphi->Eval(tpc_eta);
834  }
835  // the distortion correction can push tpc_phi outside +/- M_PI
836  if(tpc_phi < - M_PI) tpc_phi += 2.0*M_PI;
837  if(tpc_phi > M_PI) tpc_phi -= 2.0*M_PI;
838  //=============================
839 
840  double tpc_x = _tracklet_tpc->get_x();
841  double tpc_y = _tracklet_tpc->get_y();
842  double tpc_z = _tracklet_tpc->get_z();
843 
844  // Now search the silicon track list for a match in eta and phi
845  for (auto phtrk_iter_si = _track_map_silicon->begin();
846  phtrk_iter_si != _track_map_silicon->end();
847  ++phtrk_iter_si)
848  {
849  _tracklet_si = phtrk_iter_si->second;
850 
851  double si_phi = atan2(_tracklet_si->get_py(), _tracklet_si->get_px());
852  double si_eta = _tracklet_si->get_eta();
853  double si_x = _tracklet_si->get_x();
854  double si_y = _tracklet_si->get_y();
855  double si_z = _tracklet_si->get_z();
856 
857  if(Verbosity() > 2)
858  {
859  cout << " testing for a match for TPC track " << _tracklet_tpc->get_id() << " with pT " << _tracklet_tpc->get_pt()
860  << " with Si track " << _tracklet_si->get_id() << endl;
861  cout << " tpc_phi " << tpc_phi << " si_phi " << si_phi << " dphi " << tpc_phi-si_phi << " phi search " << _phi_search_win*mag << " tpc_eta " << tpc_eta
862  << " si_eta " << si_eta << " deta " << tpc_eta-si_eta << " eta search " << _eta_search_win*mag << endl;
863  std::cout << " tpc x " << tpc_x << " si x " << si_x << " tpc y " << tpc_y << " si y " << si_y << " tpc_z " << tpc_z << " si z " << si_z << std::endl;
864  std::cout << " x search " << _x_search_win*mag << " y search " << _y_search_win*mag << " z search " << _z_search_win*mag << std::endl;
865  }
866 
867  bool eta_match = false;
868  bool phi_match = false;
869  bool position_match = false;
870  if( fabs(tpc_eta - si_eta) < _eta_search_win * mag) eta_match = true;
871 
872  // PHTpcTracker has a bias in the tracklet phi that depends on charge sign, PHCASeeding does not
873  if(_is_ca_seeder)
874  {
875  if( fabs(tpc_phi - si_phi) < _phi_search_win * mag) phi_match = true;
876  }
877  else
878  {
879  // PHTpcTracker
880  //double si_pt = sqrt( pow(_tracklet_si->get_px(),2) + pow(_tracklet_si->get_py(),2) );
881  double phi_search_win_lo = fdphi->Eval(tpc_pt) * sign_phi_correction - _phi_search_win * mag;
882  double phi_search_win_hi = fdphi->Eval(tpc_pt) * sign_phi_correction + _phi_search_win * mag;
883 
884  if(Verbosity() > 10)
885  cout << " phi_search_win_lo " << phi_search_win_lo << " phi_search_win_hi " << phi_search_win_hi << endl;
886 
887  if( (tpc_phi - si_phi) > phi_search_win_lo && (tpc_phi - si_phi) < phi_search_win_hi) phi_match = true;
888  }
889 
890  if(_pp_mode)
891  {
892  if(
893  fabs(tpc_x - si_x) < _x_search_win * mag
894  && fabs(tpc_y - si_y) < _y_search_win * mag
895  )
896  position_match = true;
897  }
898  else
899  {
900  if(
901  fabs(tpc_x - si_x) < _x_search_win * mag
902  && fabs(tpc_y - si_y) < _y_search_win * mag
903  && fabs(tpc_z - si_z) < _z_search_win * mag
904  )
905  position_match = true;
906  }
907 
908  if(eta_match && phi_match && position_match)
909  {
910  // got a match, add to the list
911  // These stubs are matched in eta, phi, x and y already
912  tpc_matches.insert(std::make_pair(_tracklet_tpc->get_id(), _tracklet_si->get_id()));
913  tpc_matched_set.insert(_tracklet_tpc->get_id());
914 
915  if(Verbosity() > 1)
916  {
917  cout << " found a match for TPC track " << _tracklet_tpc->get_id() << " with Si track " << _tracklet_si->get_id() << endl;
918  cout << " tpc_phi " << tpc_phi << " si_phi " << si_phi << " phi_match " << phi_match
919  << " tpc_eta " << tpc_eta << " si_eta " << si_eta << " eta_match " << eta_match << endl;
920  std::cout << " tpc x " << tpc_x << " si x " << si_x << " tpc y " << tpc_y << " si y " << si_y << " tpc_z " << tpc_z << " si z " << si_z << std::endl;
921  }
922 
923  // temporary!
924  if(_test_windows)
925  cout << " Try_silicon: pt " << tpc_pt << " tpc_phi " << tpc_phi << " si_phi " << si_phi << " dphi " << tpc_phi-si_phi
926  << " tpc_eta " << tpc_eta << " si_eta " << si_eta << " deta " << tpc_eta-si_eta << " tpc_x " << tpc_x << " tpc_y " << tpc_y << " tpc_z " << tpc_z
927  << " dx " << tpc_x - si_x << " dy " << tpc_y - si_y << " dz " << tpc_z - si_z
928  << endl;
929  }
930  }
931  }
932 
933  // We have all of the candidates for association, but we have to deal with cases of multiple matches for a tpc tracklet
934  // We cannot move the TPC clusters independently for multiple matches, so we duplicate the tpc tracklet
935  std::multimap<unsigned int, unsigned int> additional_tpc_matches;
936  std::multimap<unsigned int, unsigned int> remove_tpc_matches;
937  std::set<unsigned int> additional_tpc_matched_set;
938  for(unsigned int tpcid : tpc_matched_set)
939  {
940  auto ret = tpc_matches.equal_range(tpcid);
941 
942  unsigned int size = std::distance(ret.first, ret.second);
943  if(Verbosity() > 1) std::cout << " tpcid " << tpcid << " number of matches " << size << std::endl;
944 
945  if(size == 1) continue;
946 
947  // make copy of all but the first TPC tracklet
948  int counter = -1;
949  for(auto it = ret.first; it != ret.second; ++it)
950  {
951  counter++;
952  if(counter == 0) continue;
953 
954  unsigned int si_id = it->second;
955  auto this_track = _track_map->get(tpcid);
956  auto newTrack = std::make_unique<SvtxTrack_v2>();
957  const unsigned int lastTrackKey = _track_map->empty() ? 0:std::prev(_track_map->end())->first;
958  if(Verbosity() > 1) cout << "Extra match, add a new track to node tree with key " << lastTrackKey + 1 << endl;
959 
960  newTrack->set_id(lastTrackKey+1);
961 
962  newTrack->set_x(this_track->get_x());
963  newTrack->set_y(this_track->get_y());
964  newTrack->set_z(this_track->get_z());
965 
966  newTrack->set_charge(this_track->get_charge());
967  newTrack->set_px(this_track->get_px());
968  newTrack->set_py(this_track->get_py());
969  newTrack->set_pz(this_track->get_pz());
970 
971  for(int i = 0; i < 6; ++i)
972  {
973  for(int j = 0; j < 6; ++j)
974  {
975  newTrack->set_error(i,j, this_track->get_error(i,j));
976  }
977  }
978 
979  _seed_track_map->addAssoc(this_track->get_id(), newTrack->get_id());
980 
981  // loop over associated clusters to get hits for TPC only, add to new track copy
982  for (SvtxTrack::ConstClusterKeyIter iter = this_track->begin_cluster_keys();
983  iter != this_track->end_cluster_keys();
984  ++iter)
985  {
986  TrkrDefs::cluskey cluster_key = *iter;
987  unsigned int trkrid = TrkrDefs::getTrkrId(cluster_key);
988  if(trkrid == TrkrDefs::tpcId)
989  {
990  newTrack->insert_cluster_key(cluster_key);
991  _assoc_container->SetClusterTrackAssoc(cluster_key, newTrack->get_id());
992  }
993  }
994 
995  if(Verbosity() > 1) cout << " -- inserting new track with id " << newTrack->get_id() << " from TPC tracklet " << tpcid << " into trackmap " << endl;
996  _track_map->insert(newTrack.get());
997 
998  // add a map remove_tpc_matches so we can remove old matches with the same tpc id
999  remove_tpc_matches.insert(std::make_pair(tpcid, si_id));
1000 
1001  additional_tpc_matches.insert(std::make_pair(newTrack->get_id(), si_id));
1002  additional_tpc_matched_set.insert(newTrack->get_id());
1003  }
1004  }
1005  for(auto [tpcid, si_id] : additional_tpc_matches)
1006  {
1007  tpc_matched_set.insert(tpcid);
1008  tpc_matches.insert(std::make_pair(tpcid, si_id));
1009  if(Verbosity() > 1) std::cout << "add to tpc_matches tpc_id " << tpcid << " si_id " << si_id << std::endl;
1010  }
1011  for(auto [tpcid, si_id] : remove_tpc_matches)
1012  {
1013  if(Verbosity() > 1) std::cout << "remove from tpc_matches tpc_id " << tpcid << " si_id " << si_id << std::endl;
1014 
1015  // Have to iterate over the multimap and examine each si_id to find the one matching remove_tpc_matches
1016  auto ret = tpc_matches.equal_range(tpcid);
1017  for(auto it = ret.first; it != ret.second; ++it)
1018  {
1019  if(it->first == tpcid && it->second == si_id)
1020  {
1021  tpc_matches.erase(it);
1022  break; // the iterator is no longer valid
1023  }
1024  }
1025  }
1026 
1027  return;
1028 }
1029 
1031  std::multimap<unsigned int, unsigned int> &tpc_matches,
1032  std::set<int> &crossing_set,
1033  std::multimap<int, std::pair<unsigned int, unsigned int>> &crossing_matches,
1034  std::map<unsigned int, int> &tpc_crossing_map )
1035 {
1036 
1037  for(auto [tpcid, si_id] : tpc_matches)
1038  {
1039  SvtxTrack *tpc_track = _track_map->get(tpcid);
1040  double tpc_z = tpc_track->get_z();
1041  double tpc_pt = tpc_track->get_pt();
1042 
1043  SvtxTrack *si_track =_track_map_silicon->get(si_id);
1044  double si_z = si_track->get_z();
1045 
1046  // this factor will increase the window size at low pT
1047  // otherwise the matching efficiency drops off at low pT
1048  // it would be better if this was a smooth function
1049  double mag = 1.0;
1050  if(tpc_pt < 6.0) mag = 2;
1051  if(tpc_pt < 3.0) mag = 4.0;
1052 
1053  // Check for a match in z
1054  if(fabs(tpc_z - si_z) < _z_search_win * mag)
1055  {
1056  //track from triggered event
1057  int crossing = 0;
1058  crossing_matches.insert(std::make_pair(crossing,std::make_pair(tpc_track->get_id(), si_track->get_id())));
1059  crossing_set.insert(crossing);
1060  tpc_crossing_map.insert(std::make_pair(tpc_track->get_id(), crossing));
1061  if(Verbosity() > 1)
1062  std::cout << " triggered: tpc_trackid " << tpc_track->get_id()
1063  << " eta " << tpc_track->get_eta()
1064  << " pT " << tpc_track->get_pt()
1065  << " si_trackid " << si_track->get_id()
1066  << " z_tpc " << tpc_track->get_z()
1067  << " z_si " << si_track->get_z()
1068  << " crossing " << crossing
1069  << std::endl;
1070  }
1071  }
1072  return;
1073 }
1074 
1076  std::multimap<unsigned int, unsigned int> &tpc_matches,
1077  std::set<int> &crossing_set,
1078  std::multimap<int, std::pair<unsigned int, unsigned int>> &crossing_matches,
1079  std::map<unsigned int, int> &tpc_crossing_map )
1080 {
1081 
1082  for(auto [tpcid, si_id] : tpc_matches)
1083  {
1084  SvtxTrack *tpc_track = _track_map->get(tpcid);
1085  double tpc_z = tpc_track->get_z();
1086  //double tpc_pt = tpc_track->get_pt();
1087 
1088  SvtxTrack *si_track =_track_map_silicon->get(si_id);
1089  double si_z = si_track->get_z();
1090 
1091  // this is an initial estimate of the bunch crossing based on the z-mismatch for this track
1092  int crossing = (int) getBunchCrossing(tpc_track->get_id(), tpc_z - si_z);
1093  crossing_matches.insert(std::make_pair(crossing,std::make_pair(tpc_track->get_id(), si_track->get_id())));
1094  crossing_set.insert(crossing);
1095  tpc_crossing_map.insert(std::make_pair(tpc_track->get_id(), crossing));
1096  if(Verbosity() > 1)
1097  std::cout << " pileup: tpc_trackid " << tpc_track->get_id()
1098  << " eta " << tpc_track->get_eta()
1099  << " pT " << tpc_track->get_pt()
1100  << " si_trackid " << si_track->get_id()
1101  << " z_tpc " << tpc_track->get_z()
1102  << " z_si " << si_track->get_z()
1103  << " crossing " << crossing
1104  << std::endl;
1105  }
1106 
1107 return;
1108 }
1109 
1111 {
1112  // loop over final track map, copy silicon clusters to corrected cluster map
1113  for (auto phtrk_iter = _track_map->begin();
1114  phtrk_iter != _track_map->end();
1115  ++phtrk_iter)
1116  {
1117  SvtxTrack *track = phtrk_iter->second;
1118 
1119  // loop over associated clusters to get keys for silicon cluster
1121  iter != track->end_cluster_keys();
1122  ++iter)
1123  {
1124  TrkrDefs::cluskey cluster_key = *iter;
1125  const unsigned int trkrid = TrkrDefs::getTrkrId(cluster_key);
1126  if(trkrid == TrkrDefs::mvtxId || trkrid == TrkrDefs::inttId)
1127  {
1128  TrkrCluster *cluster = _cluster_map->findCluster(cluster_key);
1129  if( !cluster ) continue;
1130 
1131  TrkrCluster *newclus = _corrected_cluster_map->findOrAddCluster(cluster_key)->second;
1132  newclus->CopyFrom( cluster );
1133  }
1134  }
1135  }
1136 }
1137 
1138