ECCE @ EIC 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"
9 #include "GPUTPCTrackLinearisation.h"
10 #include "GPUTPCTrackParam.h"
11 
12 #include "../PHTpcTracker/externals/kdfinder.hpp"
13 #include "../PHTpcTracker/PHTpcTrackerUtil.h"
14 
15 // trackbase_historic includes
16 #include <trackbase_historic/SvtxTrackMap.h>
17 #include <trackbase_historic/SvtxTrack_v2.h>
18 #include <trackbase_historic/SvtxTrack.h>
19 #include <trackbase_historic/SvtxVertex.h>
20 #include <trackbase_historic/SvtxVertexMap.h>
21 
22 #include <trackbase/TrkrCluster.h> // for TrkrCluster
23 #include <trackbase/TrkrClusterContainer.h>
24 #include <trackbase/TrkrDefs.h> // for getLayer, clu...
25 
26 // sPHENIX Geant4 includes
27 #include <g4detectors/PHG4CylinderCellGeom.h>
28 #include <g4detectors/PHG4CylinderCellGeomContainer.h>
29 #include <g4detectors/PHG4CylinderGeom.h>
30 #include <g4detectors/PHG4CylinderGeomContainer.h>
31 
32 // sPHENIX includes
33 #include <fun4all/Fun4AllReturnCodes.h>
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 
75 PHHybridSeeding::PHHybridSeeding(
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 
100 int PHHybridSeeding::InitializeGeometry(PHCompositeNode *topNode)
101 {
102  PHG4CylinderCellGeomContainer *cellgeos = findNode::getClass<
103  PHG4CylinderCellGeomContainer>(topNode, "CYLINDERCELLGEOM_SVTX");
104  PHG4CylinderGeomContainer *laddergeos = findNode::getClass<
105  PHG4CylinderGeomContainer>(topNode, "CYLINDERGEOM_INTT");
106  PHG4CylinderGeomContainer *mapsladdergeos = findNode::getClass<
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  {
118  PHG4CylinderCellGeomContainer::ConstRange layerrange =
119  cellgeos->get_begin_end();
120  for (PHG4CylinderCellGeomContainer::ConstIterator layeriter =
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  {
132  PHG4CylinderGeomContainer::ConstRange layerrange =
133  laddergeos->get_begin_end();
134  for (PHG4CylinderGeomContainer::ConstIterator layeriter =
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  {
146  PHG4CylinderGeomContainer::ConstRange layerrange =
147  mapsladdergeos->get_begin_end();
148  for (PHG4CylinderGeomContainer::ConstIterator layeriter =
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  {
171  PHG4CylinderCellGeomContainer::ConstRange begin_end =
172  cellgeos->get_begin_end();
173  PHG4CylinderCellGeomContainer::ConstIterator miter = begin_end.first;
174  for (; miter != begin_end.second; ++miter)
175  {
176  PHG4CylinderCellGeom *geo = miter->second;
177  _radii_all[_layer_ilayer_map_all[geo->get_layer()]] =
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  {
190  PHG4CylinderGeomContainer::ConstRange begin_end =
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;
196  _radii_all[_layer_ilayer_map_all[geo->get_layer()]] =
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  {
208  PHG4CylinderGeomContainer::ConstRange begin_end =
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
221  _radii_all[_layer_ilayer_map_all[geo->get_layer()]] =
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  }
230  return Fun4AllReturnCodes::EVENT_OK;
231 }
232 
233 int PHHybridSeeding::Process(PHCompositeNode */*topNode*/)
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
271  _search_radius2,
272  _search_angle2,
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_v2> seeds = fitter->ALICEKalmanFilter(clusterLists,false);
308  if(Verbosity()>0) cout << "nseeds: " << seeds.size() << "\n";
309  publishSeeds(seeds);
310  return Fun4AllReturnCodes::EVENT_OK;
311 }
312 
313 void PHHybridSeeding::publishSeeds(vector<SvtxTrack_v2> seeds)
314 {
315  for(size_t i=0;i<seeds.size();i++)
316  {
317  _track_map->insert(&(seeds[i]));
318  }
319 }
320 
321 int PHHybridSeeding::Setup(PHCompositeNode *topNode)
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;
331  return Fun4AllReturnCodes::ABORTEVENT;
332  }
333 
334  auto surfmaps = findNode::getClass<ActsSurfaceMaps>(topNode, "ActsSurfaceMaps");
335  if(!surfmaps)
336  {
337  std::cout << "No Acts surface maps, bailing." << std::endl;
338  return Fun4AllReturnCodes::ABORTEVENT;
339  }
340 
341  auto tGeometry = findNode::getClass<ActsTrackingGeometry>(topNode,"ActsTrackingGeometry");
342  if(!tGeometry)
343  {
344  std::cout << "No Acts tracking geometry, exiting." << std::endl;
345  return Fun4AllReturnCodes::ABORTEVENT;
346  }
347 
348  InitializeGeometry(topNode);
349  fitter = std::make_shared<ALICEKF>(topNode,_cluster_map,surfmaps, tGeometry,
350  _fieldDir,_min_fit_track_size,_max_sin_phi,Verbosity());
351  return Fun4AllReturnCodes::EVENT_OK;
352 }
353 
354 int PHHybridSeeding::End()
355 {
356  if(Verbosity()>0) cout << "Called End " << endl;
357  return Fun4AllReturnCodes::EVENT_OK;
358 }