ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
PHHoughSeeding.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file PHHoughSeeding.cc
1 
8 #include "PHHoughSeeding.h"
9 
10 #include "AssocInfoContainer.h" // for AssocInfoCont...
11 
12 // Helix Hough includes
13 #include <HelixHough/HelixKalmanState.h> // for HelixKalmanState
14 #include <HelixHough/HelixRange.h>
15 #include <HelixHough/SimpleHit3D.h>
16 #include <HelixHough/SimpleTrack3D.h>
17 #include <HelixHough/sPHENIXSeedFinder.h> // for sPHENIXSeedFi...
18 
19 
20 // trackbase_historic includes
21 #include <trackbase_historic/SvtxTrackMap.h>
22 #include <trackbase_historic/SvtxTrack_v1.h>
23 #include <trackbase_historic/SvtxVertex.h>
24 #include <trackbase_historic/SvtxVertexMap.h>
25 #include <trackbase_historic/SvtxVertex_v1.h>
26 
27 #include <trackbase/TrkrCluster.h> // for TrkrCluster
28 #include <trackbase/TrkrDefs.h> // for getLayer, clu...
29 #include <trackbase/TrkrClusterContainer.h>
30 #include <trackbase/TrkrHitSet.h>
31 #include <trackbase/TrkrHitSetContainer.h>
32 
33 // sPHENIX Geant4 includes
34 #include <g4detectors/PHG4CylinderCellGeom.h>
35 #include <g4detectors/PHG4CylinderCellGeomContainer.h>
36 #include <g4detectors/PHG4CylinderGeom.h>
37 #include <g4detectors/PHG4CylinderGeomContainer.h>
38 
39 #include <g4bbc/BbcVertex.h>
40 #include <g4bbc/BbcVertexMap.h>
41 
42 // sPHENIX includes
43 #include <fun4all/Fun4AllReturnCodes.h>
44 
45 #include <phool/PHTimer.h> // for PHTimer
46 #include <phool/getClass.h>
47 #include <phool/phool.h> // for PHWHERE
48 
49 #include <Eigen/Core> // for Matrix
50 #include <Eigen/Dense>
51 
52 //ROOT includes for debugging
53 #include <TFile.h>
54 #include <TNtuple.h>
55 
56 // standard includes
57 #include <algorithm>
58 #include <climits> // for UINT_MAX
59 #include <cmath>
60 #include <iostream>
61 #include <memory>
62 #include <set> // for set
63 #include <tuple>
64 #include <utility> // for pair, make_pair
65 
66 
67 #define LogDebug(exp) std::cout << "DEBUG: " << __FILE__ << ": " << __LINE__ << ": " << exp
68 #define LogError(exp) std::cout << "ERROR: " << __FILE__ << ": " << __LINE__ << ": " << exp
69 #define LogWarning(exp) std::cout << "WARNING: " << __FILE__ << ": " << __LINE__ << ": " << exp
70 
71 //#define _DEBUG_
72 
73 //#define _USE_ALAN_FULL_VERTEXING_
74 #define _USE_ALAN_TRACK_REFITTING_
75 
76 //#define _MEARGE_SEED_CLUSTER_
77 //#define _USE_ZERO_SEED_
78 
79 //#define _USE_CONSTANT_SEARCH_WIN_
80 
81 //#define _DO_FULL_FITTING_
82 
83 using namespace std;
84 
85 PHHoughSeeding::PHHoughSeeding(
86  const string& name,
87  unsigned int nlayers_maps,
88  unsigned int nlayers_intt,
89  unsigned int nlayers_tpc,
90  unsigned int nlayers_seeding,
91  unsigned int min_nlayers_seeding)
92  : PHTrackSeeding(name)
93  , _t_seeding(nullptr)
94  , _t_seed_init1(nullptr)
95  , _t_seed_init2(nullptr)
96  , _t_seed_init3(nullptr)
97  , _t_seeds_cleanup(nullptr)
98  , _t_translate_to_PHGenFitTrack(nullptr)
99  , _t_translate1(nullptr)
100  , _t_translate2(nullptr)
101  , _t_translate3(nullptr)
102  , _t_kalman_pat_rec(nullptr)
103  , _t_search_clusters(nullptr)
104  , _t_search_clusters_encoding(nullptr)
105  , _t_search_clusters_map_iter(nullptr)
106  , _t_track_propagation(nullptr)
107  , _t_full_fitting(nullptr)
108  , _t_output_io(nullptr)
109  , _seeding_layer()
110  , _nlayers_seeding(nlayers_seeding)
111  , _min_nlayers_seeding(min_nlayers_seeding)
112  , _radii()
113  , _material()
114  , _user_material()
115  , _magField(1.4)
116  , _reject_ghosts(true)
117  , _remove_hits(false)
118  , _min_pt(0.2)
119  , _min_z0(-14.0)
120  , _max_z0(+14.0)
121  , _max_r(1.0)
122  , _cut_on_dca(true)
123  , _dcaxy_cut(0.2)
124  , _dcaz_cut(0.2)
125  , _chi2_cut_fast_par0(10.0)
126  , _chi2_cut_fast_par1(50.0)
127  , _chi2_cut_fast_max(75.0)
128  , _chi2_cut_full(5.0)
129  , _ca_chi2_cut(5.0)
130  , _cos_angle_cut(0.985)
131  , _bin_scale(0.8)
132  , _z_bin_scale(0.8)
133  , _min_combo_hits(min_nlayers_seeding)
134  , _max_combo_hits(nlayers_seeding * 4)
135  , _pt_rescale(0.9972 / 1.00117)
136  , // 1.0
137  _fit_error_scale(_nlayers_seeding, 1.0 / sqrt(12.0))
138  , _vote_error_scale(_nlayers_seeding, 1.0)
139  , _layer_ilayer_map()
140  , _clusters()
141  , _tracks()
142  , _track_errors()
143  , _track_covars()
144  , _vertex()
145  , _tracker(nullptr)
146  , _tracker_vertex(nullptr)
147  , _tracker_etap_seed(nullptr)
148  , _tracker_etam_seed(nullptr)
149  , _vertexFinder()
150  , _bbc_vertexes(nullptr)
151  , _hit_used_map(nullptr)
152  , _hit_used_map_size(0)
153  , _geom_container_intt(nullptr)
154  , _geom_container_maps(nullptr)
155  , _analyzing_mode(false)
156  , _analyzing_file(nullptr)
157  , _analyzing_ntuple(nullptr)
158  , _max_merging_dphi(0.1)
159  , _max_merging_deta(0.1)
160  , _max_merging_dr(0.1)
161  , _max_merging_dz(0.1)
162  , _max_share_hits(3)
163  , _cut_min_pT(0.2)
164  , _nlayers_maps(nlayers_maps)
165  , _nlayers_intt(nlayers_intt)
166  , _nlayers_tpc(nlayers_tpc)
167  , _nlayers_all(_nlayers_maps + _nlayers_intt + _nlayers_tpc)
168  , _layer_ilayer_map_all()
169  , _radii_all()
170 {
171  _event = 0;
172 
173  _user_material.clear();
174  for (unsigned int i = 0; i < _nlayers_maps; ++i)
175  _user_material[i] = 0.003;
176  for (unsigned int i = _nlayers_maps; i < _nlayers_maps + _nlayers_intt; ++i)
177  _user_material[i] = 0.008;
178 
179  //int seeding_layers[] = {7,15,25,35,45,55,66};
180  int ninner_layer = _nlayers_maps + _nlayers_intt;
181  int incr_layer = floor(_nlayers_tpc / 6.);
182  int seeding_layers[] = {
183  ninner_layer,
184  ninner_layer + incr_layer * 1,
185  ninner_layer + incr_layer * 2,
186  ninner_layer + incr_layer * 3,
187  ninner_layer + incr_layer * 4,
188  ninner_layer + incr_layer * 5,
189  _nlayers_all - 1};
190  this->set_seeding_layer(seeding_layers, 7);
191 
192  _vertex_error.clear();
193  //_vertex_error.assign(3, 0.0100);
194 
195  if (_analyzing_mode)
196  {
197  cout << "Ana Mode, creating ntuples! " << endl;
198  _analyzing_file = new TFile("./PHHoughSeedingAna.root", "RECREATE");
199  _analyzing_ntuple = new TNtuple("ana_nt", "ana_nt", "pt:kappa:d:phi:dzdl:z0:nhit:ml:rec:dt");
200  cout << "Done" << endl;
201  }
202 }
203 
204 int PHHoughSeeding::Setup(PHCompositeNode* topNode)
205 {
206  //cout << PHWHERE << "Entering Setup" << endl;
207  // Start new interface ----
208  int ret = PHTrackSeeding::Setup(topNode);
209  if (ret != Fun4AllReturnCodes::EVENT_OK) return ret;
210 
211  ret = GetNodes(topNode);
212  if (ret != Fun4AllReturnCodes::EVENT_OK) return ret;
213  // End new interface ----
214  // , _nlayers_seeding(nlayers_seeding)
215  // , _min_nlayers_seeding(min_nlayers_seeding)
216  if(_nlayers_seeding == 12){
217  int seeding_layers[] = {
218  (int) (_nlayers_maps + _nlayers_intt),
219  (int) (_nlayers_maps + _nlayers_intt + 1),
220  (int) (_nlayers_maps + _nlayers_intt + 2),
221 
222  (int) (_nlayers_maps + _nlayers_intt + 9),
223  (int) (_nlayers_maps + _nlayers_intt + 10),
224  (int) (_nlayers_maps + _nlayers_intt + 11),
225 
226  (int) (_nlayers_maps + _nlayers_intt + 20),
227  (int) (_nlayers_maps + _nlayers_intt + 21),
228 
229  (int) (_nlayers_maps + _nlayers_intt + 31),
230  (int) (_nlayers_maps + _nlayers_intt + 32),
231 
232  (int) (_nlayers_maps + _nlayers_intt + 38),
233 
234  (int) (_nlayers_maps + _nlayers_intt + 45)
235  };
236  set_seeding_layer(seeding_layers, _nlayers_seeding);
237  set_min_nlayers_seeding(_min_nlayers_seeding);
238  }else{
239  int seeding_layers[] = {
240  // (int) (_nlayers_maps + _nlayers_intt),
241  // (int) (_nlayers_maps + _nlayers_intt + 6),
242  // (int) (_nlayers_maps + _nlayers_intt + 12),
243  // (int) (_nlayers_maps + _nlayers_intt + 18),
244  // (int) (_nlayers_maps + _nlayers_intt + 24),
245  // (int) (_nlayers_maps + _nlayers_intt + 30),
246  // (int) (_nlayers_maps + _nlayers_intt + 39)
247  //7,13,19,25,31,37,46
248  (int) (_nlayers_maps + _nlayers_intt + 1),
249  (int) (_nlayers_maps + _nlayers_intt + 7),
250  (int) (_nlayers_maps + _nlayers_intt + 15),
251  (int) (_nlayers_maps + _nlayers_intt + 23),
252  (int) (_nlayers_maps + _nlayers_intt + 31),
253  (int) (_nlayers_maps + _nlayers_intt + 39),
254  (int) (_nlayers_maps + _nlayers_intt + 46)};
255  set_seeding_layer(seeding_layers, _nlayers_seeding);
256  set_min_nlayers_seeding(_min_nlayers_seeding);
257  }
258 
259 
260 
261  ret = InitializeGeometry(topNode);
262  if (ret != Fun4AllReturnCodes::EVENT_OK)
263  return ret;
264 
269  _t_seeding = new PHTimer("_t_seeding");
270  _t_seeding->stop();
271 
272  _t_seed_init1 = new PHTimer("_t_seed_init1");
273  _t_seed_init1->stop();
274 
275  _t_seed_init2 = new PHTimer("_t_seed_init2");
276  _t_seed_init2->stop();
277 
278  _t_seed_init3 = new PHTimer("_t_seed_init3");
279  _t_seed_init3->stop();
280 
281  _t_seeds_cleanup = new PHTimer("_t_seeds_cleanup");
282  _t_seeds_cleanup->stop();
283 
284  _t_translate_to_PHGenFitTrack = new PHTimer("_t_translate_to_PHGenFitTrack");
285  _t_translate_to_PHGenFitTrack->stop();
286 
287  _t_translate1 = new PHTimer("_t_translate1");
288  _t_translate1->stop();
289  _t_translate2 = new PHTimer("_t_translate2");
290  _t_translate2->stop();
291  _t_translate3 = new PHTimer("_t_translate3");
292  _t_translate3->stop();
293 
294  _t_kalman_pat_rec = new PHTimer("_t_kalman_pat_rec");
295  _t_kalman_pat_rec->stop();
296 
297  _t_search_clusters = new PHTimer("_t_search_clusters");
298  _t_search_clusters->stop();
299 
300  _t_search_clusters_encoding = new PHTimer("_t_search_clusters_encoding");
301  _t_search_clusters_encoding->stop();
302 
303  _t_search_clusters_map_iter = new PHTimer("_t_search_clusters_map_iter");
304  _t_search_clusters_map_iter->stop();
305 
306  _t_track_propagation = new PHTimer("_t_track_propergation");
307  _t_track_propagation->stop();
308 
309  _t_full_fitting = new PHTimer("_t_full_fitting");
310  _t_full_fitting->stop();
311 
312  _t_output_io = new PHTimer("_t_output_io");
313  _t_output_io->stop();
314 
315  if (Verbosity() > 0)
316  {
317  cout
318  << "====================== PHHoughSeeding::InitRun() ======================"
319  << endl;
320  cout << " Magnetic field set to: " << _magField << " Tesla" << endl;
321  cout << " Number of tracking layers: " << _nlayers_seeding << endl;
322  for (unsigned int i = 0; i < _nlayers_seeding; ++i)
323  {
324  cout << " Tracking layer #" << i << " "
325  << "radius = "
326  << _radii[i] << " cm, "
327  << "material = " << _material[i]
328  << endl;
329  cout << " Tracking layer #" << i << " "
330  << "vote error scale = "
331  << _vote_error_scale[i] << ", "
332  << "fit error scale = "
333  << _fit_error_scale[i] << endl;
334  }
335  cout << " Required hits: " << _min_nlayers_seeding << endl;
336  cout << " Minimum pT: " << _min_pt << endl;
337  cout << " Fast fit chisq cut min(par0+par1/pt,max): min( "
338  << _chi2_cut_fast_par0 << " + " << _chi2_cut_fast_par1
339  << " / pt, " << _chi2_cut_fast_max << " )" << endl;
340  cout << " Maximum chisq (kalman fit): " << _chi2_cut_full << endl;
341  cout << " Cell automaton chisq: " << _ca_chi2_cut << endl;
342  cout << " Cos Angle Cut: " << _cos_angle_cut << endl;
343  cout << " Ghost rejection: " << boolalpha << _reject_ghosts
344  << noboolalpha << endl;
345  cout << " Hit removal: " << boolalpha << _remove_hits << noboolalpha
346  << endl;
347  cout << " Maximum DCA: " << boolalpha << _cut_on_dca << noboolalpha
348  << endl;
349  if (_cut_on_dca)
350  {
351  cout << " Maximum DCA cut: " << _dcaxy_cut << endl;
352  }
353  cout << " Maximum DCAZ cut: " << _dcaz_cut << endl;
354  cout << " Phi bin scale: " << _bin_scale << endl;
355  cout << " Z bin scale: " << _z_bin_scale << endl;
356  cout << " Momentum rescale factor: " << _pt_rescale << endl;
357  cout
358  << "==========================================================================="
359  << endl;
360  }
361 
362  return ret;
363 }
364 
365 int PHHoughSeeding::Process(PHCompositeNode */*topNode*/)
366 {
367  if (Verbosity() > 1)
368  {
369  cout << "PHHoughSeeding::process_event -- entered" << endl;
370  cout << "nMapsLayers = " << _nlayers_maps << endl;
371  cout << "nInttLayers = " << _nlayers_intt << endl;
372  cout << "nTPCLayers = " << _nlayers_tpc << endl;
373  cout << "n seeding layers = " << _nlayers_seeding << endl;
374  cout << "n min layers = " << _min_nlayers_seeding << endl;
375  }
376  // start fresh
377  int code;
378 
379  _clusters.clear();
380  _all_tracks.clear();
381  _all_track_errors.clear();
382  _all_track_covars.clear();
383 
384  _vertex.clear();
385  _vertex_error.clear();
386  _vertex_tracks.clear();
387 
388  /*iteration from Christof*/ {
389  // currently only one pass is made
390  _tracks.clear();
391  _track_errors.clear();
392  _track_covars.clear();
393 
394  if (Verbosity() >= 1) _t_seeding->restart();
395 
396  //-----------------------------------
397  // Translate into Helix_Hough objects
398  //-----------------------------------
399 
400  code = translate_input(); //Check if cluster is already used in here
401  if (code != Fun4AllReturnCodes::EVENT_OK) return code;
402 
403  //-----------------------------------
404  // Obtain initial vertex from SvtxVtxMap[0]
405  //-----------------------------------
406  {
407  code = vertexing();
408  if (code != Fun4AllReturnCodes::EVENT_OK) return code;
409  // here expect vertex to be better than +/- 500 um
410  }
411 
412  if(Verbosity() > 1) cout << PHWHERE << " _vertex.size() = " << _vertex.size() << endl;
413  for(unsigned int ivert = 0; ivert < _vertex.size(); ++ivert)
414  {
415  //-----------------------------------
416  // Seeding - Alan's Hough Tracking with selected layers
417  //-----------------------------------
418 
419  // We have a list of vertices, now we want to seed tracks for each vertex.
420  // loop over vertices and call full_track_seeding for each one
421  if(Verbosity() > 1) cout << "Call full_track_seeding for ivert = " << ivert << " at Z = " << _vertex[ivert][2] << endl;
422 
423  code = full_track_seeding(ivert);
424  if (code != Fun4AllReturnCodes::EVENT_OK)
425  return code;
426 
427  if (Verbosity() >= 1) _t_seeding->stop();
428 
429  _t_seed_init1->stop();
430 
431  if (Verbosity() > 1) print_timers();
432  }
433 
434  } /*end of iteration*/
435 
436  // CleanupTracksByHitPattern();
437 
438  //-----------------------------------
439  // Alan's exportation
440  //-----------------------------------
441 
442  code = export_output();
443  if (code != Fun4AllReturnCodes::EVENT_OK)
444  return code;
445 
446  ++_event;
447 
448  return Fun4AllReturnCodes::EVENT_OK;
449 }
450 
451 void PHHoughSeeding::print_timers()
452 {
453  std::cout << "=============== PHHoughSeeding::print_timers: ===============" << std::endl;
454  std::cout << "CPUSCALE Seeding time: " << _t_seeding->get_accumulated_time() / 1000. << " sec" << std::endl;
455  std::cout << "CPUSCALE Init Seed1 time: " << _t_seed_init1->get_accumulated_time() / 1000. << " sec" << std::endl;
456  std::cout << "CPUSCALE Init Seed2 time: " << _t_seed_init2->get_accumulated_time() / 1000. << " sec" << std::endl;
457  std::cout << "CPUSCALE Init Seed3 time: " << _t_seed_init3->get_accumulated_time() / 1000. << " sec" << std::endl;
458  std::cout << "\t - Seeds Cleanup: " << _t_seeds_cleanup->get_accumulated_time() / 1000. << " sec" << std::endl;
459  std::cout << "CPUSCALE Pattern recognition time: " << _t_kalman_pat_rec->get_accumulated_time() / 1000. << " sec" << std::endl;
460  std::cout << "\t - Track Translation time: " << _t_translate_to_PHGenFitTrack->get_accumulated_time() / 1000. << " sec" << std::endl;
461  std::cout << "\t - - Translation1 time: " << _t_translate1->get_accumulated_time() / 1000. << " sec" << std::endl;
462  std::cout << "\t - - Translation2 time: " << _t_translate2->get_accumulated_time() / 1000. << " sec" << std::endl;
463  std::cout << "\t - - Translation3 time: " << _t_translate3->get_accumulated_time() / 1000. << " sec" << std::endl;
464  std::cout << "\t - Cluster searching time: " << _t_search_clusters->get_accumulated_time() / 1000. << " sec" << std::endl;
465  std::cout << "\t\t - Encoding time: " << _t_search_clusters_encoding->get_accumulated_time() / 1000. << " sec" << std::endl;
466  std::cout << "\t\t - Map iteration: " << _t_search_clusters_map_iter->get_accumulated_time() / 1000. << " sec" << std::endl;
467  std::cout << "\t - Kalman updater time: " << _t_track_propagation->get_accumulated_time() / 1000. << " sec" << std::endl;
468  std::cout << "Full fitting time: " << _t_full_fitting->get_accumulated_time() / 1000. << " sec" << std::endl;
469  std::cout << "Output IO time: " << _t_output_io->get_accumulated_time() / 1000. << " sec" << std::endl;
470  std::cout << "=======================================" << std::endl;
471 }
472 
473 int PHHoughSeeding::End()
474 {
475  delete _t_seeding;
476  delete _t_seeds_cleanup;
477  delete _t_translate_to_PHGenFitTrack;
478  delete _t_translate1;
479  delete _t_translate2;
480  delete _t_translate3;
481  delete _t_kalman_pat_rec;
482  delete _t_full_fitting;
483  delete _t_search_clusters;
484  delete _t_track_propagation;
485  delete _t_output_io;
486 
487  delete _tracker_etap_seed;
488  _tracker_etap_seed = nullptr;
489  delete _tracker_etam_seed;
490  _tracker_etam_seed = nullptr;
491  delete _tracker_vertex;
492  _tracker_vertex = nullptr;
493  delete _tracker;
494  _tracker = nullptr;
495 
496 #ifdef _DEBUG_
497  LogDebug("Leaving End \n");
498 #endif
499 
500  if (_analyzing_mode)
501  {
502  cout << " cleaning up " << endl;
503  _analyzing_file->cd();
504  _analyzing_ntuple->Write();
505  _analyzing_file->Close();
506  // delete _analyzing_ntuple;
507  // delete _analyzing_file;
508  }
509 
510  return Fun4AllReturnCodes::EVENT_OK;
511 }
512 
513 void PHHoughSeeding::set_material(int layer, float value)
514 {
515  _user_material[layer] = value;
516 }
517 
518 int PHHoughSeeding::InitializeGeometry(PHCompositeNode* topNode)
519 {
520  //---------------------------------------------------------
521  // Grab Run-Dependent Detector Geometry and Configure Hough
522  //---------------------------------------------------------
523 
524  PHG4CylinderCellGeomContainer* cellgeos = findNode::getClass<
525  PHG4CylinderCellGeomContainer>(topNode, "CYLINDERCELLGEOM_SVTX");
526  PHG4CylinderGeomContainer* laddergeos = findNode::getClass<
527  PHG4CylinderGeomContainer>(topNode, "CYLINDERGEOM_INTT");
528  PHG4CylinderGeomContainer* mapsladdergeos = findNode::getClass<
529  PHG4CylinderGeomContainer>(topNode, "CYLINDERGEOM_MVTX");
530 
531  cout << "layer array size: " << _seeding_layer.size() << endl;
532  cout << "init _nlayers_seeding " << _nlayers_seeding << endl;
533  cout << "init _min_layers_seeding " << _nlayers_seeding << endl;
534  cout << "init n min layers = " << _min_nlayers_seeding << endl;
535  _nlayers_seeding = _seeding_layer.size();
536 
537  //=================================================//
538  // Initializing HelixHough objects //
539  //=================================================//
540 
541  // Since the G4 layers don't necessarily correspond to the
542  // silicon layers, and don't necessarily start from zero (argh),
543  // we create our own layers numbers that are consecutive
544  // starting from zero.
545 
546  // Now that we have two kinds of layers, I won't know in principle
547  // which type is in what order, so I figure that out now...
548 
549  _radii.assign(_nlayers_seeding, 0.0);
550  map<float, int> radius_layer_map;
551 
552  _radii_all.assign(_nlayers_all, 0.0);
553  _layer_ilayer_map.clear();
554  _layer_ilayer_map_all.clear();
555  if (cellgeos)
556  {
557  PHG4CylinderCellGeomContainer::ConstRange layerrange =
558  cellgeos->get_begin_end();
559  for (PHG4CylinderCellGeomContainer::ConstIterator layeriter =
560  layerrange.first;
561  layeriter != layerrange.second; ++layeriter)
562  {
563  radius_layer_map.insert(
564  make_pair(layeriter->second->get_radius(),
565  layeriter->second->get_layer()));
566  //cout << " making layer map for TPC " << endl;
567  }
568  }
569 
570  if (laddergeos)
571  {
572  PHG4CylinderGeomContainer::ConstRange layerrange =
573  laddergeos->get_begin_end();
574  for (PHG4CylinderGeomContainer::ConstIterator layeriter =
575  layerrange.first;
576  layeriter != layerrange.second; ++layeriter)
577  {
578  radius_layer_map.insert(
579  make_pair(layeriter->second->get_radius(),
580  layeriter->second->get_layer()));
581  //cout << " making layer map for intt " << endl;
582  }
583  }
584 
585  if (mapsladdergeos)
586  {
587  PHG4CylinderGeomContainer::ConstRange layerrange =
588  mapsladdergeos->get_begin_end();
589  for (PHG4CylinderGeomContainer::ConstIterator layeriter =
590  layerrange.first;
591  layeriter != layerrange.second; ++layeriter)
592  {
593  radius_layer_map.insert(
594  make_pair(layeriter->second->get_radius(),
595  layeriter->second->get_layer()));
596  //cout << " making layer map for mvtx " << endl;
597  }
598  }
599 
600  if (Verbosity() > 1) {
601  for (map<float, int>::const_iterator iter = radius_layer_map.begin();
602  iter != radius_layer_map.end(); ++iter) {
603  //cout << "radius_layer_map: first: " << iter->first << "; second: "
604  // << iter->second << endl;
605  }
606  }
607 
608  // now that the layer ids are sorted by radius, I can create a storage
609  // index, ilayer, that is 0..N-1 and sorted by radius
610 
611  int ilayer = 0;
612  for (map<float, int>::iterator iter = radius_layer_map.begin();
613  iter != radius_layer_map.end(); ++iter)
614  {
615  _layer_ilayer_map_all.insert(make_pair(iter->second, _layer_ilayer_map_all.size()));
616 
617  if (std::find(_seeding_layer.begin(), _seeding_layer.end(),
618  iter->second) != _seeding_layer.end())
619  {
620  _layer_ilayer_map.insert(make_pair(iter->second, ilayer));
621  ++ilayer;
622  }
623  //if(ilayer >= (int) _radii.size()) break; //yuhw
624  }
625 
626  // if (Verbosity() >= 10) {
627  // for (map<int, unsigned int>::const_iterator iter = _layer_ilayer_map_all.begin();
628  // iter != _layer_ilayer_map_all.end(); iter++) {
629  // cout << "_layer_ilayer_map_all: first: " << iter->first << "; second: "
630  // << iter->second << endl;
631  // }
632  // }
633 
634  // now we extract the information from the cellgeos first
635  if (cellgeos)
636  {
637  PHG4CylinderCellGeomContainer::ConstRange begin_end =
638  cellgeos->get_begin_end();
639  PHG4CylinderCellGeomContainer::ConstIterator miter = begin_end.first;
640  for (; miter != begin_end.second; miter++)
641  {
642  PHG4CylinderCellGeom* geo = miter->second;
643 
644  //if(cellgeo->get_layer() > (int) _radii.size() ) continue;
645 
646  // if (Verbosity() >= 2)
647  // cellgeo->identify();
648 
649  //TODO
650  _radii_all[_layer_ilayer_map_all[geo->get_layer()]] =
651  geo->get_radius() + 0.5 * geo->get_thickness();
652 
653  if (_layer_ilayer_map.find(geo->get_layer()) != _layer_ilayer_map.end())
654  {
655  _radii[_layer_ilayer_map[geo->get_layer()]] =
656  geo->get_radius();
657  }
658  }
659  }
660 
661  if (laddergeos)
662  {
663  PHG4CylinderGeomContainer::ConstRange begin_end =
664  laddergeos->get_begin_end();
665  PHG4CylinderGeomContainer::ConstIterator miter = begin_end.first;
666  for (; miter != begin_end.second; miter++)
667  {
668  PHG4CylinderGeom* geo = miter->second;
669 
670  //if(geo->get_layer() > (int) _radii.size() ) continue;
671 
672  // if (Verbosity() >= 2)
673  // geo->identify();
674 
675  _radii_all[_layer_ilayer_map_all[geo->get_layer()]] =
676  geo->get_radius() + 0.5 * geo->get_thickness();
677 
678  if (_layer_ilayer_map.find(geo->get_layer()) != _layer_ilayer_map.end())
679  {
680  _radii[_layer_ilayer_map[geo->get_layer()]] = geo->get_radius();
681  }
682  }
683  }
684 
685  if (mapsladdergeos)
686  {
687  PHG4CylinderGeomContainer::ConstRange begin_end =
688  mapsladdergeos->get_begin_end();
689  PHG4CylinderGeomContainer::ConstIterator miter = begin_end.first;
690  for (; miter != begin_end.second; miter++)
691  {
692  PHG4CylinderGeom* geo = miter->second;
693 
694  //if(geo->get_layer() > (int) _radii.size() ) continue;
695 
696  // if (Verbosity() >= 2)
697  // geo->identify();
698 
699  //TODO
700  _radii_all[_layer_ilayer_map_all[geo->get_layer()]] =
701  geo->get_radius();
702 
703  if (_layer_ilayer_map.find(geo->get_layer()) != _layer_ilayer_map.end())
704  {
705  _radii[_layer_ilayer_map[geo->get_layer()]] = geo->get_radius();
706  }
707  }
708  }
709  // set material on each layer
710 
711  _material.assign(_radii.size(), 0.03);
712 
713  map<int, float>::iterator mat_it;
714  for (map<int, float>::iterator iter = _user_material.begin();
715  iter != _user_material.end(); ++iter)
716  {
717  if (_layer_ilayer_map.find(iter->first) != _layer_ilayer_map.end())
718  {
719  _material[_layer_ilayer_map[iter->first]] = iter->second;
720  }
721  }
722  if (_tracker) delete _tracker;
723  if (_tracker_vertex) delete _tracker_vertex;
724  if (_tracker_etap_seed) delete _tracker_etap_seed;
725 
726  // initialize the pattern recogition tools
727  setup_tracker_object();
728  setup_initial_tracker_object();
729  setup_seed_tracker_objects();
730 
735  /*
736  _cells_svtx = findNode::getClass<PHG4CellContainer>(topNode,
737  "G4CELL_TPC");
738 
739  _cells_intt = findNode::getClass<PHG4CellContainer>(
740  topNode, "G4CELL_INTT");
741 
742  _cells_maps = findNode::getClass<PHG4CellContainer>(
743  topNode, "G4CELL_MVTX");
744 
745  if (!_cells_svtx and !_cells_intt and !_cells_maps)
746  {
747  if (Verbosity() >= 0)
748  {
749  LogError("No PHG4CellContainer found!");
750  }
751  return Fun4AllReturnCodes::ABORTRUN;
752  }
753  */
754 
755  _geom_container_intt = findNode::getClass<
756  PHG4CylinderGeomContainer>(topNode, "CYLINDERGEOM_INTT");
757 
758  _geom_container_maps = findNode::getClass<
759  PHG4CylinderGeomContainer>(topNode, "CYLINDERGEOM_MVTX");
760 
761  /*
762  if (!_cells_svtx && !_cells_maps && !_cells_intt)
763  {
764  cout << PHWHERE << "ERROR: Can't find any cell node!" << endl;
765  return Fun4AllReturnCodes::ABORTRUN;
766  }
767  */
768 
769  return Fun4AllReturnCodes::EVENT_OK;
770 }
771 
772 int PHHoughSeeding::setup_seed_tracker_objects()
773 {
774  float kappa_max = ptToKappa(_min_pt);
775 
776  // for the initial tracker we may not have the best guess on the vertex yet
777  // so I've doubled the search range on dca and dcaz
778 
779  std::vector<unsigned int> onezoom(5, 0);
780  std::vector<vector<unsigned int> > zoomprofile;
781  zoomprofile.assign(5, onezoom);
782  zoomprofile[0][0] = 16;
783  zoomprofile[0][1] = 1;
784  zoomprofile[0][2] = 4;
785  zoomprofile[0][3] = 8;
786  zoomprofile[0][4] = 1;
787 
788  zoomprofile[1][0] = 16;
789  zoomprofile[1][1] = 1;
790  zoomprofile[1][2] = 4;
791  zoomprofile[1][3] = 4;
792  zoomprofile[1][4] = 2;
793 
794  zoomprofile[2][0] = 4;
795  zoomprofile[2][1] = 3;
796  zoomprofile[2][2] = 2;
797  zoomprofile[2][3] = 1;
798  zoomprofile[2][4] = 3;
799 
800  for (unsigned int i = 2; i <= 3; ++i)
801  {
802  zoomprofile[i][0] = 3;
803  zoomprofile[i][1] = 3;
804  zoomprofile[i][2] = 3;
805  zoomprofile[i][3] = 3;
806  zoomprofile[i][4] = 3;
807  }
808 
809  HelixRange pos_range(0.0, 2. * M_PI, // center of rotation azimuthal angles
810  -_max_r, _max_r, // 2d dca range
811  0.0, kappa_max, // curvature range
812  0.0, 0.9, // dzdl range
813  _min_z0, _max_z0); // dca_z range
814 
815  _tracker_etap_seed = new sPHENIXSeedFinder(zoomprofile, 1, pos_range,
816  _material, _radii, _magField);
817  _tracker_etap_seed->setNLayers(_nlayers_seeding);
818  _tracker_etap_seed->requireLayers(_min_nlayers_seeding);
819  _tracker_etap_seed->setClusterStartBin(1);
820  _tracker_etap_seed->setRejectGhosts(_reject_ghosts);
821  _tracker_etap_seed->setFastChi2Cut(_chi2_cut_fast_par0, _chi2_cut_fast_par1,
822  _chi2_cut_fast_max);
823  _tracker_etap_seed->setChi2Cut(_chi2_cut_full);
824  _tracker_etap_seed->setChi2RemovalCut(_chi2_cut_full * 0.5);
825  _tracker_etap_seed->setCellularAutomatonChi2Cut(_ca_chi2_cut);
826  _tracker_etap_seed->setPrintTimings(false);
827  _tracker_etap_seed->setCutOnDca(false);
828  _tracker_etap_seed->setSmoothBack(true);
829  _tracker_etap_seed->setBinScale(_bin_scale);
830  _tracker_etap_seed->setZBinScale(_z_bin_scale);
831  _tracker_etap_seed->setRemoveHits(_remove_hits);
832  _tracker_etap_seed->setSeparateByHelicity(true);
833  _tracker_etap_seed->setMaxHitsPairs(0);
834  _tracker_etap_seed->setCosAngleCut(_cos_angle_cut);
835 
836  for (unsigned int ilayer = 0; ilayer < _fit_error_scale.size(); ++ilayer)
837  {
838  float scale1 = _fit_error_scale[ilayer];
839  float scale2 = _vote_error_scale[ilayer];
840  float scale = scale1 / scale2;
841  _tracker_etap_seed->setHitErrorScale(ilayer, scale);
842  }
843 
844  // for the initial tracker we may not have the best guess on the vertex yet
845  // so I've doubled the search range on dca and dcaz
846 
847  HelixRange neg_range(0.0, 2. * M_PI, // center of rotation azimuthal angles
848  -_max_r, _max_r, // 2d dca range
849  0.0, kappa_max, // curvature range
850  -0.9, 0.0, // dzdl range
851  _min_z0, _max_z0); // dca_z range
852 
853  _tracker_etam_seed = new sPHENIXSeedFinder(zoomprofile, 1, neg_range,
854  _material, _radii, _magField);
855  _tracker_etam_seed->setNLayers(_nlayers_seeding);
856  _tracker_etam_seed->requireLayers(_min_nlayers_seeding);
857  _tracker_etam_seed->setClusterStartBin(1);
858  _tracker_etam_seed->setRejectGhosts(_reject_ghosts);
859  _tracker_etam_seed->setFastChi2Cut(_chi2_cut_fast_par0, _chi2_cut_fast_par1,
860  _chi2_cut_fast_max);
861  _tracker_etam_seed->setChi2Cut(_chi2_cut_full);
862  _tracker_etam_seed->setChi2RemovalCut(_chi2_cut_full * 0.5);
863  _tracker_etam_seed->setCellularAutomatonChi2Cut(_ca_chi2_cut);
864  _tracker_etam_seed->setPrintTimings(false);
865  _tracker_etam_seed->setCutOnDca(false);
866  _tracker_etam_seed->setSmoothBack(true);
867  _tracker_etam_seed->setBinScale(_bin_scale);
868  _tracker_etam_seed->setZBinScale(_z_bin_scale);
869  _tracker_etam_seed->setRemoveHits(_remove_hits);
870  _tracker_etam_seed->setSeparateByHelicity(true);
871  _tracker_etam_seed->setMaxHitsPairs(0);
872  _tracker_etam_seed->setCosAngleCut(_cos_angle_cut);
873 
874  for (unsigned int ilayer = 0; ilayer < _fit_error_scale.size(); ++ilayer)
875  {
876  float scale1 = _fit_error_scale[ilayer];
877  float scale2 = _vote_error_scale[ilayer];
878  float scale = scale1 / scale2;
879  _tracker_etam_seed->setHitErrorScale(ilayer, scale);
880  }
881 
882  return Fun4AllReturnCodes::EVENT_OK;
883 }
884 
885 int PHHoughSeeding::setup_initial_tracker_object()
886 {
887  // copy of the final tracker modified to:
888  // expand the DCA search regions (2.0 cm z search > 3 sigma of BBC z vertex
889  // remove the DCA cut on the track output
890 
891  // tell the initial tracker object the phase space extent of the search region
892  // and the recursive zoom factors to utilize
893 
894  float kappa_max = ptToKappa(_min_pt);
895 
896  // for the initial tracker we may not have the best guess on the vertex yet
897  // so I've doubled the search range on dca and dcaz
898 
899  HelixRange top_range(0.0, 2. * M_PI, // center of rotation azimuthal angles
900  -1.0, +1.0, // 2d dca range
901  0.0, kappa_max, // curvature range
902  -0.9, 0.9, // dzdl range
903  -2.0, +2.0); // dca_z range
904 
905  vector<unsigned int> onezoom(5, 0);
906  vector<vector<unsigned int> > zoomprofile;
907  zoomprofile.assign(5, onezoom);
908  zoomprofile[0][0] = 16;
909  zoomprofile[0][1] = 1;
910  zoomprofile[0][2] = 4;
911  zoomprofile[0][3] = 8;
912  zoomprofile[0][4] = 1;
913 
914  zoomprofile[1][0] = 16;
915  zoomprofile[1][1] = 1;
916  zoomprofile[1][2] = 4;
917  zoomprofile[1][3] = 4;
918  zoomprofile[1][4] = 2;
919 
920  zoomprofile[2][0] = 4;
921  zoomprofile[2][1] = 3;
922  zoomprofile[2][2] = 2;
923  zoomprofile[2][3] = 1;
924  zoomprofile[2][4] = 3;
925 
926  for (unsigned int i = 2; i <= 3; ++i)
927  {
928  zoomprofile[i][0] = 3;
929  zoomprofile[i][1] = 3;
930  zoomprofile[i][2] = 3;
931  zoomprofile[i][3] = 3;
932  zoomprofile[i][4] = 3;
933  }
934 
935  _tracker_vertex = new sPHENIXSeedFinder(zoomprofile, 1, top_range,
936  _material, _radii, _magField);
937  _tracker_vertex->setNLayers(_nlayers_seeding);
938  _tracker_vertex->requireLayers(_min_nlayers_seeding);
939  _tracker_vertex->setClusterStartBin(1);
940  _tracker_vertex->setRejectGhosts(_reject_ghosts);
941  _tracker_vertex->setFastChi2Cut(_chi2_cut_fast_par0, _chi2_cut_fast_par1,
942  _chi2_cut_fast_max);
943  _tracker_vertex->setChi2Cut(_chi2_cut_full);
944  _tracker_vertex->setChi2RemovalCut(_chi2_cut_full * 0.5);
945  _tracker_vertex->setCellularAutomatonChi2Cut(_ca_chi2_cut);
946  _tracker_vertex->setPrintTimings(false);
947  _tracker_vertex->setCutOnDca(false);
948  _tracker_vertex->setSmoothBack(true);
949  _tracker_vertex->setBinScale(_bin_scale);
950  _tracker_vertex->setZBinScale(_z_bin_scale);
951  _tracker_vertex->setRemoveHits(_remove_hits);
952  _tracker_vertex->setSeparateByHelicity(true);
953  _tracker_vertex->setMaxHitsPairs(0);
954  _tracker_vertex->setCosAngleCut(_cos_angle_cut);
955 
956  for (unsigned int ilayer = 0; ilayer < _fit_error_scale.size(); ++ilayer)
957  {
958  float scale1 = _fit_error_scale[ilayer];
959  float scale2 = _vote_error_scale[ilayer];
960  float scale = scale1 / scale2;
961  _tracker_vertex->setHitErrorScale(ilayer, scale);
962  }
963 
964  return Fun4AllReturnCodes::EVENT_OK;
965 }
966 
967 int PHHoughSeeding::setup_tracker_object()
968 {
969  // input vertex must be within 500 um of final
970 
971  // tell the tracker object the phase space extent of the search region
972  // and the recursive zoom factors to utilize
973 
974  float kappa_max = ptToKappa(_min_pt);
975 
976  HelixRange top_range(0.0, 2. * M_PI, // center of rotation azimuthal angles
977  -_dcaxy_cut, _dcaxy_cut, // 2d dca range
978  0.0, kappa_max, // curvature range
979  -0.9, 0.9, // dzdl range
980  -_dcaz_cut, _dcaz_cut); // dca_z range
981 
982  vector<unsigned int> onezoom(5, 0);
983  vector<vector<unsigned int> > zoomprofile;
984  zoomprofile.assign(5, onezoom);
985  zoomprofile[0][0] = 16;
986  zoomprofile[0][1] = 1;
987  zoomprofile[0][2] = 4;
988  zoomprofile[0][3] = 8;
989  zoomprofile[0][4] = 1;
990 
991  zoomprofile[1][0] = 16;
992  zoomprofile[1][1] = 1;
993  zoomprofile[1][2] = 4;
994  zoomprofile[1][3] = 4;
995  zoomprofile[1][4] = 2;
996 
997  zoomprofile[2][0] = 4;
998  zoomprofile[2][1] = 3;
999  zoomprofile[2][2] = 2;
1000  zoomprofile[2][3] = 1;
1001  zoomprofile[2][4] = 3;
1002 
1003  for (unsigned int i = 2; i <= 3; ++i)
1004  {
1005  zoomprofile[i][0] = 3;
1006  zoomprofile[i][1] = 3;
1007  zoomprofile[i][2] = 3;
1008  zoomprofile[i][3] = 3;
1009  zoomprofile[i][4] = 3;
1010  }
1011 
1012  _tracker = new sPHENIXSeedFinder(zoomprofile, 1, top_range, _material,
1013  _radii, _magField);
1014  _tracker->setNLayers(_nlayers_seeding);
1015  _tracker->requireLayers(_min_nlayers_seeding);
1016  _tracker->setClusterStartBin(1);
1017  _tracker->setRejectGhosts(_reject_ghosts);
1018  _tracker->setFastChi2Cut(_chi2_cut_fast_par0, _chi2_cut_fast_par1,
1019  _chi2_cut_fast_max);
1020  _tracker->setChi2Cut(_chi2_cut_full);
1021  _tracker->setChi2RemovalCut(_chi2_cut_full * 0.5);
1022  _tracker->setCellularAutomatonChi2Cut(_ca_chi2_cut);
1023  _tracker->setPrintTimings(false);
1024  if (Verbosity() >= 2)
1025  _tracker->setPrintTimings(true);
1026  _tracker->setVerbosity(Verbosity());
1027  _tracker->setCutOnDca(_cut_on_dca);
1028  _tracker->setDcaCut(_dcaxy_cut);
1029  _tracker->setSmoothBack(true);
1030  _tracker->setBinScale(_bin_scale);
1031  _tracker->setZBinScale(_z_bin_scale);
1032  _tracker->setRemoveHits(_remove_hits);
1033  _tracker->setSeparateByHelicity(true);
1034  _tracker->setMaxHitsPairs(0);
1035  _tracker->setCosAngleCut(_cos_angle_cut);
1036 
1037  for (unsigned int ilayer = 0; ilayer < _fit_error_scale.size(); ++ilayer)
1038  {
1039  float scale1 = _fit_error_scale[ilayer];
1040  float scale2 = _vote_error_scale[ilayer];
1041  float scale = scale1 / scale2;
1042  _tracker->setHitErrorScale(ilayer, scale);
1043  }
1044 
1045  return Fun4AllReturnCodes::EVENT_OK;
1046 }
1047 
1048 int PHHoughSeeding::GetNodes(PHCompositeNode* topNode)
1049 {
1050  //---------------------------------
1051  // Get Objects off of the Node Tree
1052  //---------------------------------
1053 
1054  // used in fast vertexing from BBC
1055  _bbc_vertexes = findNode::getClass<BbcVertexMap>(topNode, "BbcVertexMap");
1056 
1057  /*
1058  // get node containing the digitized hits
1059  _svtxhitsmap = findNode::getClass<SvtxHitMap>(topNode, "SvtxHitMap");
1060  if (!_svtxhitsmap)
1061  {
1062  cout << PHWHERE << "ERROR: Can't find node SvtxHitMap" << endl;
1063  return Fun4AllReturnCodes::ABORTRUN;
1064  }
1065  */
1066 
1067  // if(_hit_used_map_size!=0) delete[] _hit_used_map;
1068  // _hit_used_map_size = static_cast<int>(_cluster_map->size());
1069  // _hit_used_map = new int[_hit_used_map_size];
1070  // for (Int_t i=0;i<_hit_used_map_size;i++){
1071  // _hit_used_map[i] = 0;
1072  // }
1073 
1074  return Fun4AllReturnCodes::EVENT_OK;
1075 }
1076 
1077 int PHHoughSeeding::translate_input()
1078 {
1079  _clusters.clear();
1080  int count = 0;
1081  // int count7 = 0;
1082  // int count46 = 0;
1083  int nhits[60];
1084  int nhits_all[60];
1085  for (int i = 0; i < 60; i++)
1086  {
1087  nhits[i] = 0;
1088  nhits_all[i] = 0;
1089  }
1090 
1091  // loop over all clusters
1092  //cout << "_clusters size " << _clusters.size() << endl;
1093  int nhits3d = -1;
1094  unsigned int clusid = -1;
1095 
1096  auto hitsetrange = _hitsets->getHitSets();
1097  for (auto hitsetitr = hitsetrange.first;
1098  hitsetitr != hitsetrange.second;
1099  ++hitsetitr){
1100  auto range = _cluster_map->getClusters(hitsetitr->first);
1101  for( auto clusIter = range.first; clusIter != range.second; ++clusIter ){
1102  clusid += 1;
1103  TrkrCluster *cluster = clusIter->second;
1104  TrkrDefs::cluskey cluskey = clusIter->first;
1105  unsigned int layer = TrkrDefs::getLayer(cluskey);
1106 
1107  count++;
1108 
1109  nhits_all[layer]++;
1110 
1111  unsigned int ilayer = UINT_MAX;
1112  std::map<int, unsigned int>::const_iterator it = _layer_ilayer_map.find(layer);
1113  if (it != _layer_ilayer_map.end())
1114  ilayer = it->second;
1115  if (ilayer >= _nlayers_seeding) continue;
1116 
1117  SimpleHit3D hit3d;
1118  nhits3d++;
1119  // capture the cluster keys so the cluster can be found in the cluster container
1120  hit3d.set_cluskey( cluskey );
1121 
1122  // this is just a cluster index - should be from 0 to number of clusters for seeding to work
1123  hit3d.set_id(clusid);
1124 
1125  if(Verbosity() > 40)
1126  {
1127  unsigned int layer_tmp = TrkrDefs::getLayer(cluster->getClusKey());
1128  cout << " found in seeding layer # " << ilayer << " layer " << layer_tmp << " cluskey " << cluster->getClusKey() << " clusid " << clusid << endl;
1129  }
1130 
1131  hit3d.set_layer(ilayer);
1132 
1133  hit3d.set_x(cluster->getPosition(0));
1134  hit3d.set_y(cluster->getPosition(1));
1135  hit3d.set_z(cluster->getPosition(2));
1136 
1137  // copy covariance over
1138  for (int i = 0; i < 3; ++i)
1139  {
1140  for (int j = i; j < 3; ++j)
1141  {
1142  hit3d.set_error(i, j, cluster->getError(i, j));
1143 
1144  //FIXME
1145  //hit3d.set_size(i, j, cluster->get_size(i, j)); // original
1146  hit3d.set_size(i, j, cluster->getError(i, j) * sqrt(12.)); // yuhw 2017-05-08
1147  //cout << " i " << i << " j " << j << " error " << cluster->getError(i,j) << endl;
1148  }
1149  }
1150 
1151  nhits[ilayer]++;
1152  //cout << " adding cluster " << clusid << " with key " << cluskey << endl;
1153  //hit3d.print();
1154  _clusters.push_back(hit3d);
1155  //cout << " ilayer " << ilayer << " nhits " << nhits[ilayer] << " _clusters size now " << _clusters.size() << endl;
1156  }
1157  }
1158  //cout << "_clusters size " << _clusters.size() << endl;
1159 
1160  if (Verbosity() > 10)
1161  {
1162  cout
1163  << "-------------------------------------------------------------------"
1164  << endl;
1165  cout
1166  << "PHHoughSeeding::process_event has the following input clusters:"
1167  << endl;
1168 
1169  cout << " _clusters.size = " << _clusters.size() << endl;
1170 
1171  for (unsigned int i = 0; i < _clusters.size(); ++i)
1172  {
1173  cout << "n init clusters = " << _clusters.size() << endl;
1174  _clusters[i].print();
1175  }
1176 
1177  cout
1178  << "-------------------------------------------------------------------"
1179  << endl;
1180  }
1181 
1182  if (Verbosity() >= 10)
1183  {
1184  cout << "CPUSCALE hits: " << count << endl;
1185  }
1186  if (Verbosity() >= 10)
1187  {
1188  for (int i = 0; i < 60; i++)
1189  {
1190  cout << "layer: " << i << " << hits: " << nhits[i] << " | " << nhits_all[i] << endl;
1191  }
1192  }
1193  return Fun4AllReturnCodes::EVENT_OK;
1194 }
1195 
1196 int PHHoughSeeding::fast_vertex_from_bbc()
1197 {
1198  // fail over to bbc vertex if no tracks were found...
1199  if (_bbc_vertexes)
1200  {
1201  BbcVertex* vertex = _bbc_vertexes->begin()->second;
1202 
1203  if (vertex)
1204  {
1205  _vertex[0][0] = 0.0;
1206  _vertex[0][1] = 0.0;
1207  _vertex[0][2] = vertex->get_z();
1208 
1209  if (Verbosity() > 1)
1210  cout << " initial bbc vertex guess: " << _vertex[0][0] << " "
1211  << _vertex[0][1] << " " << _vertex[0][2] << endl;
1212  }
1213  }
1214 
1215  return Fun4AllReturnCodes::EVENT_OK;
1216 }
1217 /*
1218 int PHHoughSeeding::fast_vertex_guessing()
1219 {
1220  // fast vertex guessing uses two tracker objects
1221  // one looks for postive going eta tracks, the other for negative going tracks
1222  // it asks for just a handful of each, but searches
1223  // over the broadest possible phase space for the vertex origin
1224 
1225  // limit each window to no more than N tracks
1226  unsigned int maxtracks = 10;
1227 
1228  _tracks.clear();
1229  _track_errors.clear();
1230  _track_covars.clear();
1231 
1232  // loop over initial tracking windows
1233  std::vector<SimpleTrack3D> newtracks;
1234 
1235  _tracker_etap_seed->clear();
1236  _tracker_etap_seed->findHelices(_clusters, _min_combo_hits, _max_combo_hits,
1237  newtracks, maxtracks);
1238 
1239  for (unsigned int t = 0; t < newtracks.size(); ++t)
1240  {
1241  _tracks.push_back(newtracks[t]);
1242  _track_covars.push_back((_tracker_etap_seed->getKalmanStates())[t].C);
1243  }
1244 
1245  _tracker_etap_seed->clear();
1246  newtracks.clear();
1247 
1248  _tracker_etam_seed->clear();
1249  _tracker_etam_seed->findHelices(_clusters, _min_combo_hits, _max_combo_hits,
1250  newtracks, maxtracks);
1251 
1252  for (unsigned int t = 0; t < newtracks.size(); ++t)
1253  {
1254  _tracks.push_back(newtracks[t]);
1255  _track_covars.push_back((_tracker_etam_seed->getKalmanStates())[t].C);
1256  }
1257 
1258  _tracker_etam_seed->clear();
1259  newtracks.clear();
1260 
1261  _vertex.clear();
1262  std::vector<float> vert;
1263  vert.assign(3, 0.0);
1264  _vertex.pushback(vert);
1265 
1266 
1267  if (Verbosity() > 1)
1268  cout << " seed track finding count: " << _tracks.size() << endl;
1269 
1270  if (!_tracks.empty())
1271  {
1272  // --- compute seed vertex from initial track finding --------------------
1273 
1274  double zsum = 0.0;
1275 
1276  for (unsigned int i = 0; i < _tracks.size(); ++i)
1277  {
1278  zsum += _tracks[i].z0;
1279  }
1280 
1281  _vertex[0][2] = zsum / _tracks.size();
1282 
1283  if (Verbosity() > 1)
1284  {
1285  cout << " seed track vertex pre-fit: " << _vertex[0][0] << " "
1286  << _vertex[0][1] << " " << _vertex[0][2] << endl;
1287  }
1288 
1289  // start with the average position and converge from there
1290  _vertexFinder.findVertex(_tracks, _track_covars, _vertex[0], 3.00, true);
1291  _vertexFinder.findVertex(_tracks, _track_covars, _vertex[0], 0.10, true);
1292  _vertexFinder.findVertex(_tracks, _track_covars, _vertex[0], 0.02, true);
1293  }
1294 
1295  // we don't need the tracks anymore
1296  _tracks.clear();
1297  _track_errors.clear();
1298  _track_covars.clear();
1299 
1300  if (Verbosity() > 1)
1301  {
1302  cout << " seed track vertex post-fit: " << _vertex[0][0] << " "
1303  << _vertex[0][1] << " " << _vertex[0][2] << endl;
1304  }
1305 
1306  return Fun4AllReturnCodes::EVENT_OK;
1307 }
1308 */
1309 
1310 /*
1311 int PHHoughSeeding::initial_vertex_finding()
1312 {
1313  // shift to the guess vertex position
1314  // run the tracking pattern recognition, stop after some number of tracks
1315  // have been found, fit those tracks to a vertex
1316  // nuke out tracks, leave vertex info, shift back
1317 
1318  float shift_dx = -_vertex[0];
1319  float shift_dy = -_vertex[1];
1320  float shift_dz = -_vertex[2];
1321 
1322  // shift to vertex guess position
1323  shift_coordinate_system(shift_dx, shift_dy, shift_dz);
1324 
1325  // limit each window to no more than N tracks
1326  unsigned int maxtracks = 40;
1327 
1328  // reset track storage and tracker
1329  _tracks.clear();
1330  _track_errors.clear();
1331  _track_covars.clear();
1332 
1333  _tracker_vertex->clear();
1334 
1335  // initial track finding
1336  _tracker_vertex->findHelices(_clusters, _min_combo_hits, _max_combo_hits,
1337  _tracks, maxtracks);
1338 
1339  for (unsigned int t = 0; t < _tracks.size(); ++t)
1340  {
1341  _track_covars.push_back((_tracker_vertex->getKalmanStates())[t].C);
1342  }
1343 
1344  // don't need the tracker object anymore
1345  _tracker_vertex->clear();
1346 
1347  if (Verbosity() > 1)
1348  cout << " initial track finding count: " << _tracks.size() << endl;
1349 
1350  if (!_tracks.empty())
1351  {
1352  // --- compute seed vertex from initial track finding --------------------
1353 
1354  double xsum = 0.0;
1355  double ysum = 0.0;
1356  double zsum = 0.0;
1357 
1358  for (unsigned int i = 0; i < _tracks.size(); ++i)
1359  {
1360  xsum += _tracks[i].d * cos(_tracks[i].phi);
1361  ysum += _tracks[i].d * sin(_tracks[i].phi);
1362  zsum += _tracks[i].z0;
1363  }
1364 
1365  _vertex[0] = xsum / _tracks.size();
1366  _vertex[1] = ysum / _tracks.size();
1367  _vertex[2] = zsum / _tracks.size();
1368 
1369  if (Verbosity() > 1)
1370  {
1371  cout << " initial track vertex pre-fit: " << _vertex[0] - shift_dx
1372  << " " << _vertex[1] - shift_dy << " "
1373  << _vertex[2] - shift_dz << endl;
1374  }
1375 
1376  // start with the average position and converge from there
1377  _vertexFinder.findVertex(_tracks, _track_covars, _vertex, 3.00, true);
1378  _vertexFinder.findVertex(_tracks, _track_covars, _vertex, 0.10, true);
1379  _vertexFinder.findVertex(_tracks, _track_covars, _vertex, 0.02, false);
1380  }
1381 
1382  // don't need the tracks anymore
1383  _tracks.clear();
1384  _track_errors.clear();
1385  _track_covars.clear();
1386 
1387  // shift back to the global coordinates
1388  shift_coordinate_system(-shift_dx, -shift_dy, -shift_dz);
1389 
1390  if (Verbosity() > 1)
1391  {
1392  cout << " initial track vertex post-fit: " << _vertex[0] << " "
1393  << _vertex[1] << " " << _vertex[2] << endl;
1394  }
1395 
1396  return Fun4AllReturnCodes::EVENT_OK;
1397 }
1398 */
1399 
1400 int PHHoughSeeding::vertexing()
1401 {
1402  _vertex.clear();
1403  _vertex_error.clear();
1404  _vertex_tracks.clear();
1405 
1406  // make a vector of vectors containing all of the vertex locations from the node tree
1407  for(unsigned int ivert=0;ivert<_vertex_map->size(); ++ivert)
1408  {
1409  SvtxVertex* svtx_vtx = _vertex_map->get(ivert);
1410 
1411  if(svtx_vtx->get_z() < -30.0 || svtx_vtx->get_z() > 30.0)
1412  continue;
1413 
1414  std::vector<float> this_vertex_pos;
1415  this_vertex_pos.assign(3,0.0);
1416  this_vertex_pos[0] = svtx_vtx->get_x();
1417  this_vertex_pos[1] = svtx_vtx->get_y();
1418  this_vertex_pos[2] = svtx_vtx->get_z();
1419 
1420  std::vector<float> this_vertex_error;
1421  this_vertex_error.assign(3,0.0);
1422  this_vertex_error[0] = sqrt(svtx_vtx->get_error(0,0));
1423  this_vertex_error[1] = sqrt(svtx_vtx->get_error(1,1));
1424  this_vertex_error[2] = sqrt(svtx_vtx->get_error(2,2));
1425 
1426  _vertex.push_back(this_vertex_pos);
1427  _vertex_error.push_back(this_vertex_error);
1428  _vertex_tracks.push_back( svtx_vtx->size_tracks());
1429  }
1430 
1431  if(Verbosity() > 10) cout << " vertex list has " << _vertex.size() << " Svtxmap has " << _vertex_map->size() << endl;
1432 
1433  if(_vertex.size() == 0)
1434  {
1435  cout << endl << PHWHERE << "Do not have a valid vertex, skip track seeding for this event " << endl << endl;
1436  }
1437 
1438  return Fun4AllReturnCodes::EVENT_OK;
1439 }
1440 
1441 int PHHoughSeeding::full_track_seeding(int ivert)
1442 {
1443  if(Verbosity() > 10) cout << "Entering full_track_seeding with ivert = " << ivert << endl;
1444 
1445  float shift_dx = -_vertex[ivert][0];
1446  float shift_dy = -_vertex[ivert][1];
1447  float shift_dz = -_vertex[ivert][2];
1448  // shift to initial vertex position
1449  //cout << "shifting coord system to vertex " << ivert << " dx " << shift_dx << " dy " << shift_dy << " dz " << shift_dz << endl;
1450  shift_coordinate_system(shift_dx, shift_dy, shift_dz, ivert);
1451 
1452  // _tracks is a vector of SimpleTrack3D objects
1453  // it may already contain entries from previous vertices
1454  // We do not want to overwrite those, so save them for later
1455  std::vector<SimpleTrack3D> previous_tracks = _tracks;
1456  std::vector<double> previous_track_errors = _track_errors;
1457  std::vector<Eigen::Matrix<float, 5, 5> > previous_track_covars = _track_covars;
1458 
1459  // reset track storage and tracker to use it for this vertex
1460  _tracks.clear();
1461  _track_errors.clear();
1462  _track_covars.clear();
1463 
1464  _tracker->clear();
1465  // final track finding
1466  _tracker->findHelices(_clusters, _min_combo_hits, _max_combo_hits, _tracks);
1467  if (Verbosity() >= 1)
1468  cout << "SEEDSTUDY nbefore clean (" << _min_nlayers_seeding << "): " << _tracks.size() << endl;
1469  // Cleanup Seeds
1470 #ifdef _USE_ALAN_TRACK_REFITTING_
1471 #else
1472  if (Verbosity() >= 1) _t_seeds_cleanup->restart();
1473  CleanupSeeds();
1474  if (Verbosity() >= 1) _t_seeds_cleanup->stop();
1475 #endif
1476  if (Verbosity() >= 1)
1477  cout << "SEEDSTUDY nafter clean: " << _tracks.size() << endl;
1478  for (unsigned int tt = 0; tt < _tracks.size(); ++tt)
1479  {
1480  _tracks[tt].set_vertex_id(ivert);
1481  _track_covars.push_back((_tracker->getKalmanStates())[tt].C);
1482  _track_errors.push_back(_tracker->getKalmanStates()[tt].chi2);
1483  }
1484 
1485  // we will need the tracker object below to refit the tracks... so we won't
1486  // reset it here
1487 
1488  if (Verbosity() > 1)
1489  cout << " final track count: " << _tracks.size() << endl;
1490 #ifdef _USE_ALAN_FULL_VERTEXING_
1491  if (!_tracks.empty())
1492  {
1493  if (Verbosity() > 1)
1494  {
1495  cout << " final vertex pre-fit: " << _vertex[0] - shift_dx << " "
1496  << _vertex[1] - shift_dy << " " << _vertex[2] - shift_dz
1497  << endl;
1498  }
1499 
1500  _vertexFinder.findVertex(_tracks, _track_covars, _vertex[ivert], 0.300, false);
1501  _vertexFinder.findVertex(_tracks, _track_covars, _vertex[ivert], 0.100, false);
1502  _vertexFinder.findVertex(_tracks, _track_covars, _vertex[ivert], 0.020, false);
1503  _vertexFinder.findVertex(_tracks, _track_covars, _vertex[ivert], 0.005, false);
1504 
1505  if (Verbosity() > 1)
1506  {
1507  cout << " final vertex post-fit: " << _vertex[ivert][0] - shift_dx << " "
1508  << _vertex[ivert][1] - shift_dy << " " << _vertex[ivert][2] - shift_dz
1509  << endl;
1510  }
1511  }
1512 #endif
1513  // shift back to global coordinates
1514  shift_coordinate_system(-shift_dx, -shift_dy, -shift_dz, ivert);
1515 
1516 #ifdef _USE_ALAN_TRACK_REFITTING_
1517  if (Verbosity() >= 1) _t_seeds_cleanup->restart();
1518  // we still need to update the track fields for DCA and PCA
1519  // we can do that from the final vertex position
1520 
1521  shift_dx = -_vertex[ivert][0];
1522  shift_dy = -_vertex[ivert][1];
1523  shift_dz = -_vertex[ivert][2];
1524 
1525  // shift to precision final vertex
1526  shift_coordinate_system(shift_dx, shift_dy, shift_dz, ivert);
1527 
1528  // recompute track fits to fill dca and pca + error fields
1529  std::vector<SimpleTrack3D> refit_tracks;
1530  std::vector<double> refit_errors;
1531  std::vector<Eigen::Matrix<float, 5, 5> > refit_covars;
1532 
1533  if (Verbosity() >= 1)
1534  {
1535  cout << __LINE__ << ": Event: " << _event << ": # tracks before cleanup: " << _tracks.size() << endl;
1536  }
1537 
1538  _tracker->finalize(_tracks, refit_tracks);
1539 
1540  if (Verbosity() >= 1)
1541  {
1542  cout << __LINE__ << ": Event: " << _event << ": # tracks after cleanup: " << _tracks.size() << endl;
1543  }
1544 
1545  for (unsigned int tt = 0; tt < refit_tracks.size(); ++tt)
1546  {
1547  refit_tracks[tt].set_vertex_id(ivert);
1548  refit_errors.push_back(_tracker->getKalmanStates()[tt].chi2);
1549  refit_covars.push_back(_tracker->getKalmanStates()[tt].C);
1550  }
1551 
1552  _tracks = refit_tracks;
1553  _track_errors = refit_errors;
1554  _track_covars = refit_covars;
1555 
1556  // shift back to global coordinates
1557  shift_coordinate_system(-shift_dx, -shift_dy, -shift_dz, ivert);
1558  if (Verbosity() >= 1) _t_seeds_cleanup->stop();
1559 #endif
1560 
1561  // okay now we are done with the tracker
1562  _tracker->clear();
1563 
1564  // Now we add back the tracks from previous vertices at the start of the track list
1565  previous_tracks.insert( previous_tracks.end(), _tracks.begin(), _tracks.end() );
1566  previous_track_errors.insert( previous_track_errors.end(), _track_errors.begin(), _track_errors.end() );
1567  previous_track_covars.insert( previous_track_covars.end(), _track_covars.begin(), _track_covars.end() );
1568  _tracks = previous_tracks;
1569  _track_errors = previous_track_errors;
1570  _track_covars = previous_track_covars;
1571 
1572  // is this necessary before going out of scope?
1573  previous_tracks.clear();
1574  previous_track_errors.clear();
1575  previous_track_covars.clear();
1576 
1577  if(Verbosity() > 2)
1578  {
1579  cout << "Leaving full_track_seeding with ivert = " << ivert << " _tracks.size() = " << _tracks.size() << " list of tracks:" << endl;
1580  for(unsigned int itrack = 0; itrack < _tracks.size(); ++itrack)
1581  {
1582  cout << " trackid " << itrack << " vertexid " << _tracks[itrack].vertex_id << endl;
1583  }
1584  }
1585 
1586 
1587  return Fun4AllReturnCodes::EVENT_OK;
1588 }
1589 
1590 int PHHoughSeeding::export_output()
1591 {
1592 
1593  // _tracks etc. contain the SimpleTrack3D tracks accumulated for all vertices
1594  // Each SimpleTrack3D knows which vertex ID it was seeded from
1595 
1596  _all_tracks = _tracks;
1597  _all_track_errors = _track_errors;
1598  _all_track_covars = _track_covars;
1599 
1600  if (_all_tracks.empty())
1601  return Fun4AllReturnCodes::EVENT_OK;
1602 
1603  // We have possibly multiple vertices in the event
1604  // These are stored as a vector of vertex position vectors in _vertex
1605 
1606  // clear the SvtxMap on the node tree, will be rewritten
1607  _vertex_map->clear();
1608 
1609  // loop over all collision vertices
1610  for(unsigned int ivert = 0; ivert < _vertex.size(); ++ivert)
1611  {
1612  if(Verbosity() > 10) cout << PHWHERE << "processing vertex " << ivert << endl;
1613 
1614  SvtxVertex_v1 vertex;
1615  vertex.set_t0(0.0);
1616  for (int i = 0; i < 3; ++i)
1617  vertex.set_position(i, _vertex[ivert][i]);
1618  vertex.set_chisq(0.0);
1619  vertex.set_ndof(0);
1620  vertex.set_error(0, 0, 0.0);
1621  vertex.set_error(0, 1, 0.0);
1622  vertex.set_error(0, 2, 0.0);
1623  vertex.set_error(1, 0, 0.0);
1624  vertex.set_error(1, 1, 0.0);
1625  vertex.set_error(1, 2, 0.0);
1626  vertex.set_error(2, 0, 0.0);
1627  vertex.set_error(2, 1, 0.0);
1628  vertex.set_error(2, 2, 0.0);
1629 
1630  // at this point we should already have an initial pt and pz guess...
1631  // need to translate this into the SvtxTrack object...
1632 
1633  vector<SimpleHit3D> track_hits;
1634  TrkrDefs::cluskey clusterkey;
1635 
1636  // _all_tracks is a SimpleTrack3D object filled in helixhough, track_hits is a SimpleHit3D object set equal to the simple hits associated with the track in helixhough
1637  // The SimpleHit3D knows the cluskey for its clusters
1638  // Here we make SvtxTrack objects
1639  for (unsigned int itrack = 0; itrack < _all_tracks.size(); itrack++)
1640  {
1641  // only tracks from this vertex
1642  if(_all_tracks.at(itrack).vertex_id != ivert)
1643  continue;
1644 
1645  if(Verbosity() > 10) cout << PHWHERE << " processing track " << itrack << endl;
1646  SvtxTrack_v1 track;
1647  track.set_id(itrack);
1648  track.set_vertex_id(ivert);
1649  track_hits.clear();
1650  track_hits = _all_tracks.at(itrack).hits;
1651 
1652  //cout << " insert cluster in svtxtrack " << endl;
1653  for (unsigned int ihit = 0; ihit < track_hits.size(); ihit++)
1654  {
1655  //cout << " cluster id " << track_hits.at(ihit).get_id() << " _cluster_map->size " << _cluster_map->size() << " cluskey " << track_hits.at(ihit).get_cluskey() << endl;
1656  if (track_hits.at(ihit).get_id() >= _cluster_map->size())
1657  {
1658  continue;
1659  }
1660  // Note: clusters can be accessd only by clusterkey
1661  clusterkey = track_hits.at(ihit).get_cluskey();
1662 
1663  //mark hits as used by iteration number n
1664  //_hit_used_map[track_hits.at(ihit).get_id()] = _n_iteration;
1665  _assoc_container->SetClusterTrackAssoc(clusterkey, track.get_id());
1666 
1667 #ifdef _DEBUG_
1668  TrkrCluster* cluster = _cluster_map->findCluster(clusterkey);
1669  cout
1670  << __LINE__
1671  << ": itrack: " << itrack
1672  << ": nhits: " << track_hits.size()
1673  << ": cluster key: " << clusterkey
1674  << endl;
1675  cluster->identify();
1676 #endif
1677 
1678  //TODO verify this change
1679  //if ((clusterLayer < (int) _nlayers_seeding) && (clusterLayer >= 0)) {
1680  track.insert_cluster_key(clusterkey);
1681  //}
1682  }
1683 
1684  float kappa = _all_tracks.at(itrack).kappa;
1685  float d = _all_tracks.at(itrack).d;
1686  float phi = _all_tracks.at(itrack).phi;
1687  float dzdl = _all_tracks.at(itrack).dzdl;
1688  float z0 = _all_tracks.at(itrack).z0;
1689 
1690  // track.set_helix_phi(phi);
1691  // track.set_helix_kappa(kappa);
1692  // track.set_helix_d(d);
1693  // track.set_helix_z0(z0);
1694  // track.set_helix_dzdl(dzdl);
1695 
1696  float pT = kappaToPt(kappa);
1697 
1698  float x_center = cos(phi) * (d + 1 / kappa); // x coordinate of circle center
1699  float y_center = sin(phi) * (d + 1 / kappa); // y " " " "
1700 
1701  // find helicity from cross product sign
1702  short int helicity;
1703  if ((track_hits[0].get_x() - x_center) * (track_hits[track_hits.size() - 1].get_y() - y_center) - (track_hits[0].get_y() - y_center) * (track_hits[track_hits.size() - 1].get_x() - x_center) > 0)
1704  {
1705  helicity = 1;
1706  }
1707  else
1708  {
1709  helicity = -1;
1710  }
1711  float pZ = 0;
1712  if (dzdl != 1)
1713  {
1714  pZ = pT * dzdl / sqrt(1.0 - dzdl * dzdl);
1715  }
1716  int ndf = 2 * _all_tracks.at(itrack).hits.size() - 5;
1717  track.set_chisq(_all_track_errors[itrack]);
1718  track.set_ndf(ndf);
1719  track.set_px(pT * cos(phi - helicity * M_PI / 2));
1720  track.set_py(pT * sin(phi - helicity * M_PI / 2));
1721  track.set_pz(pZ);
1722 
1723  track.set_dca2d(d);
1724  track.set_dca2d_error(sqrt(_all_track_covars[itrack](1, 1)));
1725 
1726  if (_magField > 0)
1727  {
1728  track.set_charge(helicity);
1729  }
1730  else
1731  {
1732  track.set_charge(-1.0 * helicity);
1733  }
1734 
1735  Eigen::Matrix<float, 6, 6> euclidean_cov =
1736  Eigen::Matrix<float, 6, 6>::Zero(6, 6);
1737  convertHelixCovarianceToEuclideanCovariance(_magField, phi, d, kappa,
1738  z0, dzdl, _all_track_covars[itrack], euclidean_cov);
1739 
1740  for (unsigned int row = 0; row < 6; ++row)
1741  {
1742  for (unsigned int col = 0; col < 6; ++col)
1743  {
1744  track.set_error(row, col, euclidean_cov(row, col));
1745  }
1746  }
1747 
1748  track.set_x(vertex.get_x() + d * cos(phi));
1749  track.set_y(vertex.get_y() + d * sin(phi));
1750  track.set_z(vertex.get_z() + z0);
1751 
1752 #ifdef _DEBUG_
1753  cout
1754  << __LINE__
1755  << ": itrack: " << itrack
1756  << ": nhits: " << track_hits.size()
1757  << endl;
1758 #endif
1759  //cout << " insert track in trackmap " << endl;
1760  _track_map->insert(&track);
1761  vertex.insert_track(track.get_id());
1762 
1763  if (Verbosity() > 5)
1764  {
1765  cout << "track " << itrack << " quality = " << track.get_quality()
1766  << endl;
1767  cout << "px = " << track.get_px() << " py = " << track.get_py()
1768  << " pz = " << track.get_pz() << endl;
1769  cout << " cluster keys size " << track.size_cluster_keys() << endl;
1770  }
1771  } // track loop
1772 
1773  if(Verbosity() > 2) cout << " adding vertex " << ivert << " to node tree" << endl;
1774  SvtxVertex* vtxptr = _vertex_map->insert_clone(&vertex);
1775  if (Verbosity() > 5)
1776  vtxptr->identify();
1777 
1778  } // vertex loop
1779 
1780 #ifdef _DEBUG_
1781  _track_map->identify();
1782  for(unsigned int i=0;i<_track_map->size();i++)
1783  {
1784  SvtxTrack *tr = _track_map->get(i);
1785  tr->identify();
1786  }
1787 #endif
1788 
1789  if (Verbosity() > 2)
1790  {
1791  cout << "PHHoughSeeding::process_event -- leaving process_event"
1792  << endl;
1793  }
1794 
1795  // we are done with these now...
1796  _clusters.clear();
1797  _all_tracks.clear();
1798  _all_track_errors.clear();
1799  _all_track_covars.clear();
1800  _vertex.clear();
1801  //_vertex.assign(3, 0.0);
1802 
1803  return Fun4AllReturnCodes::EVENT_OK;
1804 }
1805 
1806 float PHHoughSeeding::kappaToPt(float kappa)
1807 {
1808  return _pt_rescale * _magField / 333.6 / kappa;
1809 }
1810 
1811 float PHHoughSeeding::ptToKappa(float pt)
1812 {
1813  return _pt_rescale * _magField / 333.6 / pt;
1814 }
1815 
1816 void PHHoughSeeding::convertHelixCovarianceToEuclideanCovariance(float B,
1817  float phi, float d, float kappa, float /*z0*/, float dzdl,
1818  Eigen::Matrix<float, 5, 5> const& input,
1819  Eigen::Matrix<float, 6, 6>& output)
1820 {
1821  Eigen::Matrix<float, 6, 5> J = Eigen::Matrix<float, 6, 5>::Zero(6, 5);
1822 
1823  // phi,d,nu,z0,dzdl
1824  // -->
1825  // x,y,z,px,py,pz
1826 
1827  float nu = sqrt(kappa);
1828  float dk_dnu = 2 * nu;
1829 
1830  float cosphi = cos(phi);
1831  float sinphi = sin(phi);
1832 
1833  J(0, 0) = -sinphi * d;
1834  J(0, 1) = cosphi;
1835  J(1, 0) = d * cosphi;
1836  J(1, 1) = sinphi;
1837  J(2, 3) = 1;
1838 
1839  float pt = 0.003 * B / kappa;
1840  float dpt_dk = -0.003 * B / (kappa * kappa);
1841 
1842  J(3, 0) = -cosphi * pt;
1843  J(3, 2) = -sinphi * dpt_dk * dk_dnu;
1844  J(4, 0) = -sinphi * pt;
1845  J(4, 2) = cosphi * dpt_dk * dk_dnu;
1846 
1847  float alpha = 1. / (1. - dzdl * dzdl);
1848  float alpha_half = pow(alpha, 0.5);
1849  float alpha_3_half = alpha * alpha_half;
1850 
1851  J(5, 2) = dpt_dk * dzdl * alpha_half * dk_dnu;
1852  J(5, 4) = pt * (alpha_half + dzdl * dzdl * alpha_3_half) * dk_dnu;
1853 
1854  output = J * input * (J.transpose());
1855 }
1856 
1857 void PHHoughSeeding::shift_coordinate_system(double dx, double dy,
1858  double dz, int ivertex)
1859 {
1860 
1861  for (unsigned int ht = 0; ht < _clusters.size(); ++ht)
1862  {
1863  _clusters[ht].set_x(_clusters[ht].get_x() + dx);
1864  _clusters[ht].set_y(_clusters[ht].get_y() + dy);
1865  _clusters[ht].set_z(_clusters[ht].get_z() + dz);
1866  }
1867 
1868  for (unsigned int tt = 0; tt < _tracks.size(); ++tt)
1869  {
1870  for (unsigned int hh = 0; hh < _tracks[tt].hits.size(); ++hh)
1871  {
1872  _tracks[tt].hits[hh].set_x(_tracks[tt].hits[hh].get_x() + dx);
1873  _tracks[tt].hits[hh].set_y(_tracks[tt].hits[hh].get_y() + dy);
1874  _tracks[tt].hits[hh].set_z(_tracks[tt].hits[hh].get_z() + dz);
1875  }
1876  }
1877 
1878  // This has to be modified to work on a specific vertex
1879 
1880  _vertex[ivertex][0] += dx;
1881  _vertex[ivertex][1] += dy;
1882  _vertex[ivertex][2] += dz;
1883 
1884  return;
1885 }
1886 
1887 int PHHoughSeeding::CleanupSeedsByHitPattern()
1888 {
1889  std::vector<SimpleTrack3D> _tracks_cleanup;
1890  _tracks_cleanup.clear();
1891 
1892  if (Verbosity() >= 1)
1893  {
1894  cout << __LINE__ << ": Event: " << _event << ": # tracks before cleanup: " << _tracks.size() << endl;
1895  }
1896 
1897  /*
1898  std::vector<double> _track_errors_cleanup;
1899  _track_errors_cleanup.clear();
1900  std::vector<Eigen::Matrix<float, 5, 5> > _track_covars_cleanup;
1901  _track_covars_cleanup.clear();
1902 
1903  std::vector<HelixKalmanState> _kalman_states_cleanup;
1904  _kalman_states_cleanup.clear();
1905  */
1906 
1907  typedef std::tuple<int, int, int, int> KeyType;
1908  typedef std::multimap<KeyType, unsigned int> MapKeyTrkID;
1909 
1910  std::set<KeyType> keys;
1911  std::vector<bool> v_track_used;
1912  MapKeyTrkID m_key_itrack;
1913 
1914  typedef std::set<unsigned int> TrackList;
1915 
1916  std::set<unsigned int> OutputList;
1917  std::multimap<int, unsigned int> hitIdTrackList;
1918 
1919  unsigned int max_hit_id = 0;
1920  //For each hit make list of all associated tracks
1921 
1922  std::vector<bool> good_track;
1923  // printf("build hit track map\n");
1924  for (unsigned int itrack = 0; itrack < _tracks.size(); ++itrack)
1925  {
1926  good_track.push_back(true);
1927  SimpleTrack3D track = _tracks[itrack];
1928  for (SimpleHit3D hit : track.hits)
1929  {
1930  hitIdTrackList.insert(std::make_pair(hit.get_id(), itrack));
1931  if (hit.get_id() > max_hit_id) max_hit_id = hit.get_id();
1932  }
1933  }
1934  // printf("build track duplicate map\n");
1935  //Check Tracks for duplicates by looking for hits shared
1936  for (unsigned int itrack = 0; itrack < _tracks.size(); ++itrack)
1937  {
1938  if (good_track[itrack] == false) continue; //already checked this one
1939  if (OutputList.count(itrack) > 0) continue; //already got this one
1940 
1941  SimpleTrack3D track = _tracks[itrack];
1942 
1943  int trackid_max_nhit = itrack;
1944  unsigned int max_nhit = track.hits.size();
1945  int onhit = track.hits.size();
1946 
1947  TrackList tList;
1948  for (SimpleHit3D hit : track.hits)
1949  {
1950  int nmatch = hitIdTrackList.count(hit.get_id());
1951  if (nmatch > 1)
1952  {
1953  multimap<int, unsigned int>::iterator it = hitIdTrackList.find(hit.get_id());
1954  //Loop over track matches and add them to the list, select longest in the process
1955  for (; it != hitIdTrackList.end(); ++it)
1956  {
1957  unsigned int match_trackid = (*it).second;
1958  if (match_trackid == itrack) continue; //original track
1959  if (good_track[match_trackid] == false) continue;
1960  tList.insert(match_trackid);
1961  SimpleTrack3D mtrack = _tracks[match_trackid];
1962  }
1963  }
1964  }
1965  // int tlsize = tList.size();
1966 
1967  // cout << "remove bad matches " << tList.size() << "itrk: " << itrack << endl;
1968  //loop over matches and remove matches with too few shared hits
1969  TrackList mergeList;
1970  for (unsigned int match : tList)
1971  {
1972  // cout << "processing " << match << " of " << tList.size() << " itrk " << itrack << endl;
1973  if (match == itrack) continue;
1974  if (good_track[match] == false) continue;
1975 
1976  SimpleTrack3D mtrack = _tracks[match]; //matched track
1977  int mnhit = mtrack.hits.size();
1978  std::set<unsigned int> HitList;
1979  //put hits from both tracks in a set
1980  for (SimpleHit3D hit : track.hits) HitList.insert(hit.get_id());
1981  for (SimpleHit3D hit : mtrack.hits) HitList.insert(hit.get_id());
1982  //set stores only unique hits, tracks overlap if:
1983  int sumnhit = HitList.size();
1984  if (sumnhit < (onhit + mnhit - 3))
1985  { // more than 3 overlaps
1986  //not enough overlap, drop track from list
1987  //tList.erase(match);
1988  //good_track[match] = false;
1989  if (sumnhit != onhit)
1990  { //no subset
1991  mergeList.insert(match);
1992  }
1993  }
1994  }
1995 
1996  tList.clear();
1997  // cout << "flag bad matches done " << mergeList.size() << " itrk " << itrack << endl;
1998  //loop over matches and flag all tracks bad except the longest
1999  std::set<unsigned int> MergedHitList;
2000  if (mergeList.size() == 0)
2001  {
2002  for (SimpleHit3D hit : track.hits) MergedHitList.insert(hit.get_id());
2003  }
2004  // cout << "merge good matches itrk " << itrack << " #" << mergeList.size() << endl;
2005  for (unsigned int match : mergeList)
2006  {
2007  if (match == itrack) continue;
2008  if (good_track[match] == false) continue;
2009  // cout << " adding " << match << endl;
2010  //check number of shared hits
2011  //get tracks
2012 
2013  SimpleTrack3D mtrack = _tracks[match]; //matched track
2014  if (mtrack.hits.size() > max_nhit)
2015  {
2016  max_nhit = mtrack.hits.size();
2017  trackid_max_nhit = match;
2018  good_track[itrack] = false;
2019  }
2020  else
2021  {
2022  good_track[match] = false;
2023  }
2024  for (SimpleHit3D hit : track.hits) MergedHitList.insert(hit.get_id());
2025  for (SimpleHit3D hit : mtrack.hits) MergedHitList.insert(hit.get_id());
2026  }
2027 
2028  // int ntracks = _tracks.size();
2029  //int outtracks = OutputList.size();
2030  // printf("CLEANUP: itrack: %5d(%d) => %5d matches max %d(%d) tracks kept: %d\n",
2031  // itrack, ntracks,tlsize, max_nhit, trackid_max_nhit, outtracks);
2032 
2033  // printf("keep track %d\n",trackid_max_nhit);
2034  //add merged hit list to merged track
2035  if (OutputList.count(trackid_max_nhit) == 0)
2036  {
2037  _tracks_cleanup.push_back(_tracks[trackid_max_nhit]);
2038 
2039  /* _kalman_states_cleanup.push_back((_tracker->getKalmanStates())[trackid_max_nhit]);
2040  _track_covars_cleanup.push_back((_tracker->getKalmanStates())[trackid_max_nhit].C);
2041  _track_errors_cleanup.push_back(_tracker->getKalmanStates()[trackid_max_nhit].chi2);
2042  */
2043  }
2044  OutputList.insert(trackid_max_nhit);
2045 
2046  _tracks_cleanup.back().hits.clear();
2047 
2048  for (unsigned int hitID : MergedHitList)
2049  {
2050  SimpleHit3D hit;
2051  hit.set_id(hitID);
2052  _tracks_cleanup.back().hits.push_back(hit);
2053  }
2054  }
2055 
2056  _tracks.clear();
2057  _tracks = _tracks_cleanup;
2058 
2059  /* _track_errors.clear();
2060  _track_errors = _track_errors_cleanup;
2061  _track_covars.clear();
2062  _track_covars = _track_covars_cleanup;
2063  _tracker->getKalmanStates().clear();
2064  for(auto &kstate : _kalman_states_cleanup){
2065  _tracker->getKalmanStates().push_back(kstate);
2066  }
2067  */
2068 
2069  if (Verbosity() >= 1)
2070  {
2071  cout << __LINE__ << ": Event: " << _event << endl;
2072  cout << ": # tracks after cleanup: " << _tracks.size() << " ol:" << OutputList.size() << endl;
2073  }
2074 
2075  return Fun4AllReturnCodes::EVENT_OK;
2076 }
2077 
2078 int PHHoughSeeding::CleanupTracksByHitPattern()
2079 {
2080  std::vector<SimpleTrack3D> _tracks_cleanup;
2081  _tracks_cleanup.clear();
2082 
2083  // if(Verbosity() >= 1)
2084  {
2085  cout << __LINE__ << ": Event: " << _event << ": # tracks before cleanup: " << _tracks.size() << endl;
2086  }
2087 
2088  std::vector<double> _track_errors_cleanup;
2089  _track_errors_cleanup.clear();
2090  std::vector<Eigen::Matrix<float, 5, 5> > _track_covars_cleanup;
2091  _track_covars_cleanup.clear();
2092 
2093  std::vector<HelixKalmanState> _kalman_states_cleanup;
2094  _kalman_states_cleanup.clear();
2095 
2096  typedef std::tuple<int, int, int, int> KeyType;
2097  typedef std::multimap<KeyType, unsigned int> MapKeyTrkID;
2098 
2099  std::set<KeyType> keys;
2100  std::vector<bool> v_track_used;
2101  MapKeyTrkID m_key_itrack;
2102 
2103  typedef std::set<unsigned int> TrackList;
2104 
2105  std::set<unsigned int> OutputList;
2106  std::multimap<int, unsigned int> hitIdTrackList;
2107 
2108  unsigned int max_hit_id = 0;
2109  //For each hit make list of all associated tracks
2110 
2111  std::vector<bool> good_track;
2112  // printf("build hit track map\n");
2113  for (unsigned int itrack = 0; itrack < _tracks.size(); ++itrack)
2114  {
2115  good_track.push_back(true);
2116  SimpleTrack3D track = _tracks[itrack];
2117  for (SimpleHit3D hit : track.hits)
2118  {
2119  hitIdTrackList.insert(std::make_pair(hit.get_id(), itrack));
2120  if (hit.get_id() > max_hit_id) max_hit_id = hit.get_id();
2121  }
2122  }
2123  // printf("build track duplicate map\n");
2124  //Check Tracks for duplicates by looking for hits shared
2125  for (unsigned int itrack = 0; itrack < _tracks.size(); ++itrack)
2126  {
2127  if (good_track[itrack] == false) continue; //already checked this one
2128  if (OutputList.count(itrack) > 0) continue; //already got this one
2129 
2130  SimpleTrack3D track = _tracks[itrack];
2131 
2132  int trackid_max_nhit = itrack;
2133  unsigned int max_nhit = track.hits.size();
2134  int onhit = track.hits.size();
2135 
2136  TrackList tList;
2137  for (SimpleHit3D hit : track.hits)
2138  {
2139  int nmatch = hitIdTrackList.count(hit.get_id());
2140  if (nmatch > 1)
2141  {
2142  multimap<int, unsigned int>::iterator it = hitIdTrackList.find(hit.get_id());
2143  //Loop over track matches and add them to the list, select longest in the process
2144  for (; it != hitIdTrackList.end(); ++it)
2145  {
2146  unsigned int match_trackid = (*it).second;
2147  if (match_trackid == itrack) continue; //original track
2148  if (good_track[match_trackid] == false) continue;
2149  tList.insert(match_trackid);
2150  SimpleTrack3D mtrack = _tracks[match_trackid];
2151  }
2152  }
2153  }
2154  // int tlsize = tList.size();
2155 
2156  // cout << "remove bad matches " << tList.size() << "itrk: " << itrack << endl;
2157  //loop over matches and remove matches with too few shared hits
2158  TrackList mergeList;
2159  for (unsigned int match : tList)
2160  {
2161  // cout << "processing " << match << " of " << tList.size() << " itrk " << itrack << endl;
2162  if (match == itrack) continue;
2163  if (good_track[match] == false) continue;
2164 
2165  SimpleTrack3D mtrack = _tracks[match]; //matched track
2166  int mnhit = mtrack.hits.size();
2167  std::set<unsigned int> HitList;
2168  //put hits from both tracks in a set
2169  for (SimpleHit3D hit : track.hits) HitList.insert(hit.get_id());
2170  for (SimpleHit3D hit : mtrack.hits) HitList.insert(hit.get_id());
2171  //set stores only unique hits, tracks overlap if:
2172  int sumnhit = HitList.size();
2173  if (sumnhit < (onhit + mnhit - 3))
2174  { // more than 3 overlaps
2175  //not enough overlap, drop track from list
2176  //tList.erase(match);
2177  //good_track[match] = false;
2178  if (sumnhit != onhit)
2179  { //no subset
2180  mergeList.insert(match);
2181  }
2182  }
2183  }
2184 
2185  tList.clear();
2186  // cout << "flag bad matches done " << mergeList.size() << " itrk " << itrack << endl;
2187  //loop over matches and flag all tracks bad except the longest
2188  std::set<unsigned int> MergedHitList;
2189  if (mergeList.size() == 0)
2190  {
2191  for (SimpleHit3D hit : track.hits) MergedHitList.insert(hit.get_id());
2192  }
2193  // cout << "merge good matches itrk " << itrack << " #" << mergeList.size() << endl;
2194  for (unsigned int match : mergeList)
2195  {
2196  if (match == itrack) continue;
2197  if (good_track[match] == false) continue;
2198  // cout << " adding " << match << endl;
2199  //check number of shared hits
2200  //get tracks
2201 
2202  SimpleTrack3D mtrack = _tracks[match]; //matched track
2203  if (mtrack.hits.size() > max_nhit)
2204  {
2205  max_nhit = mtrack.hits.size();
2206  trackid_max_nhit = match;
2207  good_track[itrack] = false;
2208  }
2209  else
2210  {
2211  good_track[match] = false;
2212  }
2213  for (SimpleHit3D hit : track.hits) MergedHitList.insert(hit.get_id());
2214  for (SimpleHit3D hit : mtrack.hits) MergedHitList.insert(hit.get_id());
2215  }
2216 
2217  // int ntracks = _tracks.size();
2218  //int outtracks = OutputList.size();
2219  // printf("CLEANUP: itrack: %5d(%d) => %5d matches max %d(%d) tracks kept: %d\n",
2220  // itrack, ntracks,tlsize, max_nhit, trackid_max_nhit, outtracks);
2221 
2222  // printf("keep track %d\n",trackid_max_nhit);
2223  //add merged hit list to merged track
2224  if (OutputList.count(trackid_max_nhit) == 0)
2225  {
2226  _tracks_cleanup.push_back(_tracks[trackid_max_nhit]);
2227 
2228  _kalman_states_cleanup.push_back((_tracker->getKalmanStates())[trackid_max_nhit]);
2229  _track_covars_cleanup.push_back((_tracker->getKalmanStates())[trackid_max_nhit].C);
2230  _track_errors_cleanup.push_back(_tracker->getKalmanStates()[trackid_max_nhit].chi2);
2231  }
2232  OutputList.insert(trackid_max_nhit);
2233 
2234  _tracks_cleanup.back().hits.clear();
2235 
2236  for (unsigned int hitID : MergedHitList)
2237  {
2238  SimpleHit3D hit;
2239  hit.set_id(hitID);
2240  _tracks_cleanup.back().hits.push_back(hit);
2241  }
2242  }
2243 
2244  _tracks.clear();
2245  _tracks = _tracks_cleanup;
2246 
2247  _track_errors.clear();
2248  _track_errors = _track_errors_cleanup;
2249  _track_covars.clear();
2250  _track_covars = _track_covars_cleanup;
2251  _tracker->getKalmanStates().clear();
2252  for (auto& kstate : _kalman_states_cleanup)
2253  {
2254  _tracker->getKalmanStates().push_back(kstate);
2255  }
2256 
2257  if (Verbosity() >= 1)
2258  {
2259  cout << __LINE__ << ": Event: " << _event << endl;
2260  cout << ": # tracks after cleanup: " << _tracks.size() << " ol:" << OutputList.size() << endl;
2261  }
2262 
2263  return Fun4AllReturnCodes::EVENT_OK;
2264 }
2265 
2266 int PHHoughSeeding::CleanupSeeds()
2267 {
2268  std::vector<SimpleTrack3D> _tracks_cleanup;
2269  _tracks_cleanup.clear();
2270 
2271  typedef std::tuple<int, int, int, int> KeyType;
2272  typedef std::multimap<KeyType, unsigned int> MapKeyTrkID;
2273 
2274  std::set<KeyType> keys;
2275  std::vector<bool> v_track_used;
2276  MapKeyTrkID m_key_itrack;
2277 
2278 #ifdef _DEBUG_
2279  cout << __LINE__ << ": CleanupSeeds: Event: " << _event << endl;
2280 #endif
2281 
2282  for (unsigned int itrack = 0; itrack < _tracks.size(); ++itrack)
2283  {
2284  SimpleTrack3D track = _tracks[itrack];
2285 
2286  int id = track.d / _max_merging_dr;
2287  int iz = track.z0 / _max_merging_dz;
2288  int iphi = track.phi / _max_merging_dphi;
2289  int idzdl = track.dzdl / _max_merging_deta;
2290 
2291 #ifdef _DEBUG_
2292 // printf("itrack: %d: \t%e, \t%e, \t%e, \t%e, %5d, %5d, %5d, %5d \n",
2293 // itrack,
2294 // track.d, track.z0, track.phi, track.dzdl,
2295 // id, iz, iphi, idzdl
2296 // );
2297 #endif
2298 
2299  KeyType key = std::make_tuple(id, iz, iphi, idzdl);
2300 
2301  keys.insert(key);
2302  m_key_itrack.insert(
2303  std::make_pair(key,
2304  itrack));
2305 
2306  v_track_used.push_back(false);
2307  }
2308 
2309 #ifdef _DEBUG_
2310  for (auto it = m_key_itrack.begin();
2311  it != m_key_itrack.end();
2312  ++it)
2313  {
2314  KeyType key = it->first;
2315  unsigned int itrack = it->second;
2316 
2317  int id = std::get<0>(key);
2318  int iz = std::get<1>(key);
2319  int iphi = std::get<2>(key);
2320  int idzdl = std::get<3>(key);
2321 
2322  SimpleTrack3D track = _tracks[itrack];
2323 
2324  cout << __LINE__ << endl;
2325  printf("itrack: %5u => {%5d, %5d, %5d, %5d} \n",
2326  itrack,
2327  id, iz, iphi, idzdl);
2328  }
2329 #endif
2330 
2331  for (KeyType key : keys)
2332  {
2333  unsigned int itrack = m_key_itrack.equal_range(key).first->second;
2334 
2335 #ifdef _DEBUG_
2336  cout << "---------------------------------------------------\n";
2337  cout << __LINE__ << ": processing: " << itrack << endl;
2338  cout << "---------------------------------------------------\n";
2339 #endif
2340 
2341  if (v_track_used[itrack] == true)
2342  continue;
2343 #ifdef _DEBUG_
2344  cout << __LINE__ << ": itrack: " << itrack << ": {";
2345 #endif
2346  std::set<unsigned int> hitIDs;
2347  for (SimpleHit3D hit : _tracks[itrack].hits)
2348  {
2349  hitIDs.insert(hit.get_id());
2350 #ifdef _DEBUG_
2351  cout << hit.get_id() << ", ";
2352 #endif
2353  }
2354 #ifdef _DEBUG_
2355  cout << "}" << endl;
2356 #endif
2357 
2359  std::vector<unsigned int> v_related_tracks;
2360  for (int id = std::get<0>(key) - 1; id <= std::get<0>(key) + 1;
2361  ++id)
2362  {
2363  for (int iz = std::get<1>(key) - 1;
2364  iz <= std::get<1>(key) + 1; ++iz)
2365  {
2366  for (int iphi = std::get<2>(key) - 1;
2367  iphi <= std::get<2>(key) + 1; ++iphi)
2368  {
2369  for (int idzdl = std::get<3>(key) - 1;
2370  idzdl <= std::get<3>(key) + 1; ++idzdl)
2371  {
2372  KeyType key_temp = std::make_tuple(id, iz, iphi, idzdl);
2373 
2374  if (m_key_itrack.find(key_temp) != m_key_itrack.end())
2375  {
2376  for (auto it =
2377  m_key_itrack.equal_range(key_temp).first;
2378  it != m_key_itrack.equal_range(
2379  key_temp)
2380  .second;
2381  ++it)
2382  {
2383  if (it->second == itrack)
2384  continue;
2385 
2386  unsigned int share_hits = 0;
2387  for (SimpleHit3D hit : _tracks[it->second].hits)
2388  {
2389  unsigned int hitID = hit.get_id();
2390  if (std::find(
2391  hitIDs.begin(),
2392  hitIDs.end(),
2393  hitID) != hitIDs.end())
2394  {
2395  ++share_hits;
2396  if (share_hits > _max_share_hits)
2397  {
2398  v_related_tracks.push_back(it->second);
2399 #ifdef _DEBUG_
2400  cout << __LINE__ << ": rel track: " << it->second << ": {";
2401  for (SimpleHit3D hit_3d : _tracks[it->second].hits)
2402  {
2403  cout << hit_3d.get_id() << ", ";
2404  }
2405  cout << "}" << endl;
2406 #endif
2407  break;
2408  }
2409  }
2410  } //loop to find common hits
2411  }
2412  //#ifdef _DEBUG_
2413  // cout << __LINE__ << ": ";
2414  // printf("{%5d, %5d, %5d, %5d} => {%d, %d} \n", id,
2415  // iz, iphi, idzdl,
2416  // m_key_itrack.equal_range(key_temp).first->second,
2417  // m_key_itrack.equal_range(key_temp).second->second);
2418  //#endif
2419  }
2420  }
2421  }
2422  }
2423  }
2424 
2425  if (v_related_tracks.size() == 0)
2426  {
2427  _tracks_cleanup.push_back(_tracks[itrack]);
2428  }
2429  else
2430  {
2431  _tracks_cleanup.push_back(_tracks[itrack]);
2432 
2433  _tracks_cleanup.back().hits.clear();
2434 
2435 #ifdef _DEBUG_
2436  int n_merge_track = 1;
2437  cout << __LINE__ << ": nclusters before merge: " << hitIDs.size() << endl;
2438 #endif
2439 
2441  //std::set<unsigned int> hitIDs;
2442  for (unsigned int irel : v_related_tracks)
2443  {
2444  if (v_track_used[irel] == true) continue;
2445 
2447  //if(irel == itrack) continue;
2448 
2449 #ifdef _DEBUG_
2450  ++n_merge_track;
2451 #endif
2452 
2453 #ifdef _MEARGE_SEED_CLUSTER_
2454  SimpleTrack3D track = _tracks[irel];
2455  for (SimpleHit3D hit : track.hits)
2456  {
2457  hitIDs.insert(hit.get_id());
2458  }
2459 #endif
2460  v_track_used[irel] = true;
2461  }
2462 
2463 #ifdef _DEBUG_
2464  cout << __LINE__ << ": # tracks merged: " << n_merge_track << endl;
2465  cout << "{ ";
2466 #endif
2467  for (unsigned int hitID : hitIDs)
2468  {
2469  SimpleHit3D hit;
2470  hit.set_id(hitID);
2471 #ifdef _DEBUG_
2472  cout << hitID << ", ";
2473 #endif
2474  _tracks_cleanup.back().hits.push_back(hit);
2475  }
2476 #ifdef _DEBUG_
2477  cout << "}" << endl;
2478  cout << __LINE__ << ": nclusters after merge: " << hitIDs.size() << endl;
2479  cout << __LINE__ << ": nclusters after merge: " << _tracks_cleanup.back().hits.size() << endl;
2480 #endif
2481  }
2482 
2483  v_track_used[itrack] = true;
2484  }
2485 
2486 #ifdef _DEBUG_
2487  cout << __LINE__ << ": Event: " << _event << endl;
2488  cout << ": # tracks before cleanup: " << _tracks.size() << endl;
2489  cout << ": # tracks after cleanup: " << _tracks_cleanup.size() << endl;
2490 #endif
2491  _tracks.clear();
2492  _tracks = _tracks_cleanup;
2493 
2494  return Fun4AllReturnCodes::EVENT_OK;
2495 }