Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
PHHybridSeeding.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file PHHybridSeeding.cc
1 
8 #include "PHHybridSeeding.h"
10 #include "GPUTPCTrackParam.h"
11 
12 #include "../PHTpcTracker/externals/kdfinder.hpp"
13 #include "../PHTpcTracker/PHTpcTrackerUtil.h"
14 
15 // trackbase_historic includes
21 
22 #include <trackbase/TrkrCluster.h> // for TrkrCluster
24 #include <trackbase/TrkrDefs.h> // for getLayer, clu...
25 
26 // sPHENIX Geant4 includes
31 
32 // sPHENIX includes
34 
35 #include <phool/PHTimer.h> // for PHTimer
36 #include <phool/getClass.h>
37 #include <phool/phool.h> // for PHWHERE
38 
39 //ROOT includes
40 #include <TVector3.h> // for TVector3
41 #include <TFile.h>
42 #include <TNtuple.h>
43 
44 #include <Eigen/Core>
45 #include <Eigen/Dense>
46 
47 #include <algorithm>
48 #include <cmath>
49 #include <iostream>
50 #include <numeric>
51 #include <utility> // for pair, make_pair
52 #include <vector>
53 #include <algorithm> // for find
54 #include <unordered_set>
55 
56 // forward declarations
57 class PHCompositeNode;
58 
59 #define _DEBUG_
60 
61 #if defined(_DEBUG_)
62 #define LogDebug(exp) std::cout << "DEBUG: " << __FILE__ << ": " << __LINE__ << ": " << exp
63 #else
64 #define LogDebug(exp) (void)0
65 #endif
66 
67 #define LogError(exp) std::cout << "ERROR: " << __FILE__ << ": " << __LINE__ << ": " << exp
68 #define LogWarning(exp) std::cout << "WARNING: " << __FILE__ << ": " << __LINE__ << ": " << exp
69 
70 //end
71 
72 typedef std::vector<TrkrDefs::cluskey> keylist;
73 using namespace std;
74 
76  const string &name,
77  double maxSinPhi,
78  double fieldDir,
79  double search_radius1,
80  double search_angle1,
81  size_t min_track_size1,
82  double search_radius2,
83  double search_angle2,
84  size_t min_track_size2,
85  size_t nthreads
86  )
87  : PHTrackSeeding(name)
88  , _max_sin_phi(maxSinPhi)
89  , _fieldDir(fieldDir)
90  , _search_radius1(search_radius1)
91  , _search_angle1(search_angle1)
92  , _min_track_size1(min_track_size1)
93  , _search_radius2(search_radius2)
94  , _search_angle2(search_angle2)
95  , _min_track_size2(min_track_size2)
96  , _nthreads(nthreads)
97 {
98 }
99 
101 {
103  PHG4TpcCylinderGeomContainer>(topNode, "CYLINDERCELLGEOM_SVTX");
105  PHG4CylinderGeomContainer>(topNode, "CYLINDERGEOM_INTT");
107  PHG4CylinderGeomContainer>(topNode, "CYLINDERGEOM_MVTX");
108 
109  //_nlayers_seeding = _seeding_layer.size();
110  //_radii.assign(_nlayers_seeding, 0.0);
111  map<double, int> radius_layer_map;
112 
113  _radii_all.assign(60, 0.0);
114  _layer_ilayer_map.clear();
115  _layer_ilayer_map_all.clear();
116  if (cellgeos)
117  {
119  cellgeos->get_begin_end();
121  layerrange.first;
122  layeriter != layerrange.second; ++layeriter)
123  {
124  radius_layer_map.insert(
125  make_pair(layeriter->second->get_radius(),
126  layeriter->second->get_layer()));
127  }
128  }
129 
130  if (laddergeos)
131  {
133  laddergeos->get_begin_end();
135  layerrange.first;
136  layeriter != layerrange.second; ++layeriter)
137  {
138  radius_layer_map.insert(
139  make_pair(layeriter->second->get_radius(),
140  layeriter->second->get_layer()));
141  }
142  }
143 
144  if (mapsladdergeos)
145  {
147  mapsladdergeos->get_begin_end();
149  layerrange.first;
150  layeriter != layerrange.second; ++layeriter)
151  {
152  radius_layer_map.insert(
153  make_pair(layeriter->second->get_radius(),
154  layeriter->second->get_layer()));
155  }
156  }
157  for (map<double, int>::iterator iter = radius_layer_map.begin();
158  iter != radius_layer_map.end(); ++iter)
159  {
160  _layer_ilayer_map_all.insert(make_pair(iter->second, _layer_ilayer_map_all.size()));
161 
162  /*if (std::find(_seeding_layer.begin(), _seeding_layer.end(),
163  iter->second) != _seeding_layer.end())
164  {
165  _layer_ilayer_map.insert(make_pair(iter->second, ilayer));
166  ++ilayer;
167  }*/
168  }
169  if (cellgeos)
170  {
172  cellgeos->get_begin_end();
173  PHG4TpcCylinderGeomContainer::ConstIterator miter = begin_end.first;
174  for (; miter != begin_end.second; ++miter)
175  {
176  PHG4TpcCylinderGeom *geo = miter->second;
178  geo->get_radius() + 0.5 * geo->get_thickness();
179 
180  /*if (_layer_ilayer_map.find(geo->get_layer()) != _layer_ilayer_map.end())
181  {
182  _radii[_layer_ilayer_map[geo->get_layer()]] =
183  geo->get_radius();
184  }*/
185  }
186  }
187 
188  if (laddergeos)
189  {
191  laddergeos->get_begin_end();
192  PHG4CylinderGeomContainer::ConstIterator miter = begin_end.first;
193  for (; miter != begin_end.second; ++miter)
194  {
195  PHG4CylinderGeom *geo = miter->second;
197  geo->get_radius() + 0.5 * geo->get_thickness();
198 
199  /*if (_layer_ilayer_map.find(geo->get_layer()) != _layer_ilayer_map.end())
200  {
201  _radii[_layer_ilayer_map[geo->get_layer()]] = geo->get_radius();
202  }*/
203  }
204  }
205 
206  if (mapsladdergeos)
207  {
209  mapsladdergeos->get_begin_end();
210  PHG4CylinderGeomContainer::ConstIterator miter = begin_end.first;
211  for (; miter != begin_end.second; ++miter)
212  {
213  PHG4CylinderGeom *geo = miter->second;
214 
215  //if(geo->get_layer() > (int) _radii.size() ) continue;
216 
217  // if (Verbosity() >= 2)
218  // geo->identify();
219 
220  //TODO
222  geo->get_radius();
223 
224  /*if (_layer_ilayer_map.find(geo->get_layer()) != _layer_ilayer_map.end())
225  {
226  _radii[_layer_ilayer_map[geo->get_layer()]] = geo->get_radius();
227  }*/
228  }
229  }
231 }
232 
234 {
235  // wipe previous vertex coordinates
236  _vertex_x.clear();
237  _vertex_y.clear();
238  _vertex_z.clear();
239  _vertex_xerr.clear();
240  _vertex_yerr.clear();
241  _vertex_zerr.clear();
242  _vertex_ids.clear();
243  // fill new vertex coordinates
244  for(map<unsigned int, SvtxVertex*>::iterator iter = _vertex_map->begin(); iter != _vertex_map->end(); ++iter)
245  {
246  SvtxVertex* v = dynamic_cast<SvtxVertex*>(iter->second->CloneMe());
247  _vertex_x.push_back(v->get_x());
248  _vertex_y.push_back(v->get_y());
249  _vertex_z.push_back(v->get_z());
250  _vertex_xerr.push_back(sqrt(v->get_error(0,0)));
251  _vertex_yerr.push_back(sqrt(v->get_error(1,1)));
252  _vertex_zerr.push_back(sqrt(v->get_error(2,2)));
253  _vertex_ids.push_back(v->get_id());
254  }
255  if(Verbosity()>1) cout << "vertices:\n";
256  for(size_t i=0;i<_vertex_x.size();i++)
257  {
258  if(Verbosity()>1) cout << "(" << _vertex_x[i] << "," << _vertex_y[i] << "," << _vertex_z[i] << ")\n";
259  }
260  vector<vector<double>> kdhits(PHTpcTrackerUtil::convert_clusters_to_hits(_cluster_map,_hitsets));
261  vector<vector<double> > unused_hits;
262  vector<vector<vector<double> > > kdtracks;
263 
264  bool print_stats = (Verbosity()>0);
265 
266  kdtracks = kdfinder::find_tracks_iterative<double>(kdhits, unused_hits,
267  _search_radius1, /* max distance in cm*/
268  _search_angle1, /* triplet angle */
269  _min_track_size1, /* min hits to keep track */
270  // first iteration
273  _min_track_size2, // second iteration params
274  _nthreads,
275  print_stats);
276  if(Verbosity()>0) cout << "n_kdtracks: " << kdtracks.size() << "\n";
277  vector<keylist> clusterLists;
278  for(auto track : kdtracks)
279  {
280  keylist k;
281  for(auto kdhit : track)
282  {
283  // see PHTpcTracker/PHTpcTrackerUtil.cc; this recovers the cluster key, apparently
284  k.push_back(*((int64_t*)&kdhit[3]));
285  }
286  clusterLists.push_back(k);
287  }
288  for(auto& clist : clusterLists)
289  {
290  if(Verbosity()>1) cout << "front layer: " << (int)TrkrDefs::getLayer(clist.front()) << " back layer: " << (int)TrkrDefs::getLayer(clist.back());
291  if(clist.size()>1 && ((int)TrkrDefs::getLayer(clist.front()))<((int)TrkrDefs::getLayer(clist.back())))
292  {
293  if(Verbosity()>1) cout << "reversing\n";
294  std::reverse(clist.begin(),clist.end());
295  }
296  if(Verbosity()>1) cout << "final layer order:\n";
297  for(auto c : clist) if(Verbosity()>1) cout << (int)TrkrDefs::getLayer(c) << endl;
298  }
299  for(auto clist : clusterLists)
300  {
301  if(Verbosity()>1) cout << "layers:\n";
302  for(auto c : clist)
303  {
304  if(Verbosity()>1) cout << (int)TrkrDefs::getLayer(c) << endl;
305  }
306  }
307  vector<SvtxTrack_v3> seeds = fitter->ALICEKalmanFilter(clusterLists,false);
308  if(Verbosity()>0) cout << "nseeds: " << seeds.size() << "\n";
309  publishSeeds(seeds);
311 }
312 
313 void PHHybridSeeding::publishSeeds(vector<SvtxTrack_v3> seeds)
314 {
315  for(size_t i=0;i<seeds.size();i++)
316  {
317  _track_map->insert(&(seeds[i]));
318  }
319 }
320 
322 {
323  if(Verbosity()>0) cout << "Called Setup" << endl;
324  if(Verbosity()>0) cout << "topNode:" << topNode << endl;
325  PHTrackSeeding::Setup(topNode);
326 
327  _vertex_map = findNode::getClass<SvtxVertexMap>(topNode, "SvtxVertexMap");
328  if (!_vertex_map)
329  {
330  cerr << PHWHERE << " ERROR: Can't find SvtxVertexMap." << endl;
332  }
333 
334  auto surfmaps = findNode::getClass<ActsSurfaceMaps>(topNode, "ActsSurfaceMaps");
335  if(!surfmaps)
336  {
337  std::cout << "No Acts surface maps, bailing." << std::endl;
339  }
340 
341  auto tGeometry = findNode::getClass<ActsTrackingGeometry>(topNode,"ActsTrackingGeometry");
342  if(!tGeometry)
343  {
344  std::cout << "No Acts tracking geometry, exiting." << std::endl;
346  }
347 
348  InitializeGeometry(topNode);
349  fitter = std::make_shared<ALICEKF>(topNode,_cluster_map,surfmaps, tGeometry,
352 }
353 
355 {
356  if(Verbosity()>0) cout << "Called End " << endl;
358 }