EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
PHGhostRejection.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file PHGhostRejection.cc
1 #include "PHGhostRejection.h"
2 
3 #include "PHGhostRejection.h"
4 
6 
7 #include <trackbase/TrkrCluster.h> // for TrkrCluster
8 #include <trackbase/TrkrDefs.h> // for cluskey, getLayer, TrkrId
11 #include <trackbase_historic/SvtxTrack.h> // for SvtxTrack, SvtxTrack::C...
13 
15 
16 #include <phool/getClass.h>
17 #include <phool/phool.h>
18 
19 #include <cmath> // for sqrt, fabs, atan2, cos
20 #include <iostream> // for operator<<, basic_ostream
21 #include <map> // for map
22 #include <set> // for _Rb_tree_const_iterator
23 #include <utility> // for pair, make_pair
24 
25 //____________________________________________________________________________..
27  : SubsysReco(name)
28 {
29 
30 }
31 
32 //____________________________________________________________________________..
34 {
35 
36 }
37 
38 //____________________________________________________________________________..
40 {
41  int ret = GetNodes(topNode);
42  if (ret != Fun4AllReturnCodes::EVENT_OK) return ret;
43 
44  return ret;
45 }
46 
47 //____________________________________________________________________________..
49 {
50 
51  if(Verbosity() > 0)
52  std::cout << PHWHERE << " Beginning track map size " << _track_map->size() << std::endl;
53 
54  // Try to eliminate repeated tracks
55 
57 
58  if(Verbosity() > 0)
59  std::cout << "Track map size after deleting ghost tracks: " << _track_map->size() << std::endl;
60 
61 
63 }
64 
66 {
68 }
69 
71 {
72 
73  _track_map = findNode::getClass<SvtxTrackMap>(topNode, _track_map_name);
74  if (!_track_map)
75  {
76  std::cout << PHWHERE << " ERROR: Can't find SvtxTrackMap: " << std::endl;
78  }
79 
81 }
82 
84 {
85  std::set<unsigned int> matches_set;
86  std::multimap<unsigned int, unsigned int> matches;
87  for (auto tr1_iter = _track_map->begin();
88  tr1_iter != _track_map->end();
89  ++tr1_iter)
90  {
91  auto track1 = (tr1_iter)->second;
92 
93  for (auto tr2_iter = tr1_iter;
94  tr2_iter != _track_map->end();
95  ++tr2_iter)
96  {
97  if((tr2_iter)->first == (tr1_iter)->first) continue;
98 
99  auto track2 = (tr2_iter)->second;
100  if(
101  fabs( track1->get_phi() - track2->get_phi() ) < _phi_cut &&
102  fabs( track1->get_eta() - track2->get_eta() ) < _eta_cut &&
103  fabs( track1->get_x() - track2->get_x() ) < _x_cut &&
104  fabs( track1->get_y() - track2->get_y() ) < _y_cut &&
105  fabs( track1->get_z() - track2->get_z() ) < _z_cut
106  )
107  {
108  matches_set.insert(tr1_iter->first);
109  matches.insert( std::pair( (tr1_iter)->first, (tr2_iter)->first) );
110 
111  if(Verbosity() > 1)
112  std::cout << "Found match for tracks " << (tr1_iter)->first << " and " << (tr2_iter)->first << std::endl;
113  }
114  }
115  }
116 
117  std::set<unsigned int> ghost_reject_list;
118 
119  for(auto set_it : matches_set)
120  {
121  if(ghost_reject_list.find(set_it) != ghost_reject_list.end()) continue; // already rejected
122 
123  auto match_list = matches.equal_range(set_it);
124 
125  auto tr1 = _track_map->get(set_it);
126  double best_qual = tr1->get_chisq() / tr1->get_ndf();
127  unsigned int best_track = set_it;
128 
129  if(Verbosity() > 1)
130  std::cout << " ****** start checking track " << set_it << " with best quality " << best_qual << " best_track " << best_track << std::endl;
131 
132  for (auto it=match_list.first; it!=match_list.second; ++it)
133  {
134  if(Verbosity() > 1)
135  std::cout << " match of track " << it->first << " to track " << it->second << std::endl;
136 
137  auto tr2 = _track_map->get(it->second);
138 
139  // Check that these two tracks actually share the same clusters, if not skip this pair
140 
141  bool is_same_track = checkClusterSharing(tr1, tr2);
142  if(!is_same_track) continue;
143 
144  // which one has the best quality?
145  double tr2_qual = tr2->get_chisq() / tr2->get_ndf();
146  if(Verbosity() > 1)
147  {
148  std::cout << " Compare: best quality " << best_qual << " track 2 quality " << tr2_qual << std::endl;
149  std::cout << " tr1: phi " << tr1->get_phi() << " eta " << tr1->get_eta()
150  << " x " << tr1->get_x() << " y " << tr1->get_y() << " z " << tr1->get_z() << std::endl;
151  std::cout << " tr2: phi " << tr2->get_phi() << " eta " << tr2->get_eta()
152  << " x " << tr2->get_x() << " y " << tr2->get_y() << " z " << tr2->get_z() << std::endl;
153  }
154 
155  if(tr2_qual < best_qual)
156  {
157  if(Verbosity() > 1)
158  std::cout << " --------- Track " << it->second << " has better quality, erase track " << best_track << std::endl;
159  ghost_reject_list.insert(best_track);
160  best_qual = tr2_qual;
161  best_track = it->second;
162  }
163  else
164  {
165  if(Verbosity() > 1)
166  std::cout << " --------- Track " << best_track << " has better quality, erase track " << it->second << std::endl;
167  ghost_reject_list.insert(it->second);
168  }
169 
170  }
171  if(Verbosity() > 1)
172  std::cout << " best track " << best_track << " best_qual " << best_qual << std::endl;
173 
174  }
175 
176  // delete ghost tracks
177  for(auto it : ghost_reject_list)
178  {
179  if(Verbosity() > 1)
180  std::cout << " erasing track ID " << it << std::endl;
181 
182  _track_map->erase(it);
183  }
184 
185  return;
186 }
187 
189 {
190  // Check that tr1 and tr2 share many clusters
191 
192  bool is_same_track = false;
193 
194  std::multimap<TrkrDefs::cluskey, unsigned int> cluskey_map;
195  std::vector<TrkrDefs::cluskey> clusterkeys;
196 
197  for (SvtxTrack::ConstClusterKeyIter key_iter = tr1->begin_cluster_keys();
198  key_iter != tr1->end_cluster_keys();
199  ++key_iter)
200  {
201  TrkrDefs::cluskey cluster_key = *key_iter;
202  if(Verbosity() > 2) std::cout << " track id: " << tr1->get_id() << " adding clusterkey " << cluster_key << std::endl;
203  cluskey_map.insert( std::make_pair(cluster_key, tr1->get_id()) );
204  clusterkeys.push_back(cluster_key);
205  }
206 
207  for (SvtxTrack::ConstClusterKeyIter key_iter = tr2->begin_cluster_keys();
208  key_iter != tr2->end_cluster_keys();
209  ++key_iter)
210  {
211  TrkrDefs::cluskey cluster_key = *key_iter;
212  if(Verbosity() > 2) std::cout << " track id: " << tr2->get_id() << " adding clusterkey " << cluster_key << std::endl;
213  cluskey_map.insert( std::make_pair(cluster_key, tr2->get_id()) );
214  }
215 
216  unsigned int nclus = clusterkeys.size();
217 
218  unsigned int nclus_used = 0;
219  for (TrkrDefs::cluskey cluskey : clusterkeys)
220  {
221  if(cluskey_map.count(cluskey)>0) nclus_used++;
222  }
223 
224  if( (float) nclus_used / (float) nclus > 0.5)
225  is_same_track = true;
226 
227  if(Verbosity() > 1)
228  std::cout << " tr1 " << tr1->get_id() << " tr2 " << tr2->get_id() << " nclus_used " << nclus_used << " nclus " << nclus << std::endl;
229  if(Verbosity() > 0)
230  if(!is_same_track) std::cout << " ***** not the same track! ********" << " tr1 " << tr1->get_id() << " tr2 "
231  << tr2->get_id() << " nclus_used " << nclus_used << " nclus " << nclus << std::endl;
232 
233  return is_same_track;
234 }