ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
PHGenFitTrkProp.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file PHGenFitTrkProp.cc
1 
8 #include "PHGenFitTrkProp.h"
9 
10 #include "AssocInfoContainer.h"
11 
12 #include <trackbase/TrkrCluster.h> // for TrkrCluster
14 #include <trackbase/TrkrDefs.h> // for cluskey, getL...
15 #include <trackbase/TrkrHitSet.h>
17 
22 
23 // sPHENIX Geant4 includes
28 //
29 
30 #include <intt/CylinderGeomIntt.h>
31 #include <intt/InttDefs.h>
32 
33 #include <mvtx/CylinderGeom_Mvtx.h>
34 #include <mvtx/MvtxDefs.h>
35 
36 #include <g4bbc/BbcVertexMap.h>
37 
38 #include <phfield/PHFieldUtility.h>
39 
40 #include <phgeom/PHGeomUtility.h>
41 
42 // sPHENIX includes
44 
45 #include <phool/getClass.h>
46 #include <phool/PHTimer.h> // for PHTimer
47 #include <phool/phool.h> // for PHWHERE
48 
49 // GenFit
50 #include <GenFit/EventDisplay.h> // for EventDisplay
51 #include <GenFit/MeasuredStateOnPlane.h>
52 #include <GenFit/RKTrackRep.h>
53 #include <GenFit/Track.h>
54 #include <phgenfit/SpacepointMeasurement.h>
55 
56 #include <phgenfit/Fitter.h>
57 #include <phgenfit/Measurement.h> // for Measurement
58 #include <phgenfit/PlanarMeasurement.h>
59 #include <phgenfit/Track.h>
60 
61 //ROOT includes for debugging
62 #include <TFile.h>
63 #include <TMatrixDSymfwd.h> // for TMatrixDSym
64 #include <TMatrixTSym.h> // for TMatrixTSym
65 #include <TMatrixTUtils.h> // for TMatrixTRow
66 #include <TNtuple.h>
67 #include <TVector3.h> // for TVector3
68 #include <TVectorDfwd.h> // for TVectorD
69 #include <TVectorT.h> // for TVectorT
70 
71 // standard includes
72 #include <cassert> // for assert
73 #include <climits> // for UINT_MAX
74 #include <cmath>
75 #include <cstdlib> // for exit
76 #include <iostream>
77 #include <memory>
78 
79 class PHField;
80 class TGeoManager;
81 namespace genfit { class AbsTrackRep; }
82 
83 #define LogDebug(exp) std::cout << "DEBUG: " << __FILE__ << ": " << __LINE__ << ": " << exp
84 #define LogError(exp) std::cout << "ERROR: " << __FILE__ << ": " << __LINE__ << ": " << exp
85 #define LogWarning(exp) std::cout << "WARNING: " << __FILE__ << ": " << __LINE__ << ": " << exp
86 
87 //#define _DEBUG_
88 
89 //#define _USE_ALAN_FULL_VERTEXING_
90 #define _USE_ALAN_TRACK_REFITTING_
91 
92 //#define _MEARGE_SEED_CLUSTER_
93 //#define _USE_ZERO_SEED_
94 
95 //#define _USE_CONSTANT_SEARCH_WIN_
96 
97 //#define _DO_FULL_FITTING_
98 
99 using namespace std;
100 
101 #ifdef _DEBUG_
102 ofstream fout_kalman_pull("kalman_pull.txt");
103 ofstream fout_chi2("chi2.txt");
104 #endif
105 
107  const string& name,
108  unsigned int nlayers_maps,
109  unsigned int nlayers_intt,
110  unsigned int nlayers_tpc,
111  unsigned int nlayers_micromegas)
112  : PHTrackPropagating(name)
113  , _nlayers_maps(nlayers_maps)
114  , _nlayers_intt(nlayers_intt)
115  , _nlayers_tpc(nlayers_tpc)
116  , _nlayers_micromegas(nlayers_micromegas)
117  , _nlayers_all(_nlayers_maps + _nlayers_intt + _nlayers_tpc + _nlayers_micromegas)
118  , _firstlayer_maps(0)
119  , _firstlayer_intt(_firstlayer_maps + _nlayers_maps)
120  , _firstlayer_tpc(_firstlayer_intt + _nlayers_intt)
121  , _firstlayer_micromegas(_firstlayer_tpc + _nlayers_tpc)
122  , _use_beam(false)
123 {
124  _event = 0;
125 }
126 
128 {
129  // Start new interface ----
130 
131  int ret = PHTrackPropagating::Setup(topNode);
132  if (ret != Fun4AllReturnCodes::EVENT_OK) return ret;
133 
134  ret = GetNodes(topNode);
135  if (ret != Fun4AllReturnCodes::EVENT_OK) return ret;
136  // End new interface ----
137 
138  ret = InitializeGeometry(topNode);
139  if (ret != Fun4AllReturnCodes::EVENT_OK)
140  return ret;
141 
142  ret = InitializePHGenFit(topNode);
143  if (ret != Fun4AllReturnCodes::EVENT_OK)
144  return ret;
145 
146  if (_analyzing_mode)
147  {
148  cout << "Ana Mode, creating ntuples! " << endl;
149  _analyzing_file = new TFile("./PHGenFitTrkProp.root", "RECREATE");
150  // _analyzing_ntuple = new TNtuple("ana_nt","ana_nt","spt:seta:sphi:pt:eta:phi:layer:ncand:nmeas");
151  _analyzing_ntuple = new TNtuple("ana_nt", "ana_nt", "pt:kappa:d:phi:dzdl:z0:nhit:ml:rec:dt");
152  cout << "Done" << endl;
153  }
154 
155  // ret = InitializeGeometry(topNode);
156  // if(ret != Fun4AllReturnCodes::EVENT_OK)
157  // return ret;
158 
162  for (int layer = 0; layer < _nlayers_all; ++layer)
163  {
164  _search_wins_phi.insert(std::make_pair(layer, _search_win_phi));
165  _search_wins_theta.insert(std::make_pair(layer, _search_win_theta));
166  _max_incr_chi2s.insert(std::make_pair(layer, _max_incr_chi2));
167  }
168 
169 #ifdef _DEBUG_
170  for (int layer = 0; layer < _nlayers_all; ++layer)
171  {
172  cout
173  << __LINE__
174  << ": layer: " << layer
175  << ": search_wins_rphi: " << _search_wins_phi[layer]
176  << ": search_wins_z: " << _search_wins_theta[layer]
177  << ": max_incr_chi2: " << _max_incr_chi2s[layer]
178  << endl;
179  }
180 #endif
181 
182  _t_seeds_cleanup = new PHTimer("_t_seeds_cleanup");
184 
185  _t_translate_to_PHGenFitTrack = new PHTimer("_t_translate_to_PHGenFitTrack");
187 
188  _t_translate1 = new PHTimer("_t_translate1");
189  _t_translate1->stop();
190  _t_translate2 = new PHTimer("_t_translate2");
191  _t_translate2->stop();
192  _t_translate3 = new PHTimer("_t_translate3");
193  _t_translate3->stop();
194 
195  _t_kalman_pat_rec = new PHTimer("_t_kalman_pat_rec");
197 
198  _t_search_clusters = new PHTimer("_t_search_clusters");
200 
201  _t_search_clusters_encoding = new PHTimer("_t_search_clusters_encoding");
203 
204  _t_search_clusters_map_iter = new PHTimer("_t_search_clusters_map_iter");
206 
207  _t_track_propagation = new PHTimer("_t_track_propergation");
209 
210  _t_full_fitting = new PHTimer("_t_full_fitting");
212 
213  _t_output_io = new PHTimer("_t_output_io");
214  _t_output_io->stop();
215 
216  return ret;
217 }
218 
220 {
221  if (Verbosity() > 10)
222  {
223  cout << "PHGenFitTrkProp::process_event -- entered" << endl;
224  cout << "nMapsLayers = " << _nlayers_maps << endl;
225  cout << "nInttLayers = " << _nlayers_intt << endl;
226  cout << "nTPCLayers = " << _nlayers_tpc << endl;
227  }
228  // start fresh
229  _gftrk_hitkey_map.clear();
230 
231  _vertex.clear();
232  _vertex_error.clear();
233 
234  if (_vertex_map)
235  {
236  for(unsigned int ivert=0; ivert<_vertex_map->size(); ++ivert)
237  {
238  SvtxVertex* vertex = _vertex_map->get(ivert);
239 
240  std::vector<float> v = { vertex->get_x(), vertex->get_y(), vertex->get_z() };
241  _vertex.push_back(v);
242 
243  std::vector<float> v_err = { sqrt(vertex->get_error(0,0)), sqrt(vertex->get_error(1,1)), sqrt(vertex->get_error(2,2)) };
244  _vertex_error.push_back(v_err);
245  }
246  }
247 
248  {
249  //-----------------------------------
250  // Kalman track propagating
251  //-----------------------------------
252  if (Verbosity() >= 1) _t_kalman_pat_rec->restart();
253  int ret = KalmanTrkProp();
254  if (ret != Fun4AllReturnCodes::EVENT_OK)
255  return ret;
256  if (Verbosity() >= 1) _t_kalman_pat_rec->stop();
257  }
258  if (Verbosity() > 1) print_timers();
259 
260  ++_event;
261 
263 }
264 
266 {
267  std::cout << "=============== PHGenFitTrkProp::print_timers: ===============" << std::endl;
268  std::cout << "\t - Seeds Cleanup: " << _t_seeds_cleanup->get_accumulated_time() / 1000. << " sec" << std::endl;
269  std::cout << "CPUSCALE Pattern recognition time: " << _t_kalman_pat_rec->get_accumulated_time() / 1000. << " sec" << std::endl;
270  std::cout << "\t - Track Translation time: " << _t_translate_to_PHGenFitTrack->get_accumulated_time() / 1000. << " sec" << std::endl;
271  std::cout << "\t - - Translation1 time: " << _t_translate1->get_accumulated_time() / 1000. << " sec" << std::endl;
272  std::cout << "\t - - Translation2 time: " << _t_translate2->get_accumulated_time() / 1000. << " sec" << std::endl;
273  std::cout << "\t - - Translation3 time: " << _t_translate3->get_accumulated_time() / 1000. << " sec" << std::endl;
274  std::cout << "\t - Cluster searching time: " << _t_search_clusters->get_accumulated_time() / 1000. << " sec" << std::endl;
275  std::cout << "\t\t - Encoding time: " << _t_search_clusters_encoding->get_accumulated_time() / 1000. << " sec" << std::endl;
276  std::cout << "\t\t - Map iteration: " << _t_search_clusters_map_iter->get_accumulated_time() / 1000. << " sec" << std::endl;
277  std::cout << "\t - Kalman updater time: " << _t_track_propagation->get_accumulated_time() / 1000. << " sec" << std::endl;
278  std::cout << "Full fitting time: " << _t_full_fitting->get_accumulated_time() / 1000. << " sec" << std::endl;
279  std::cout << "Output IO time: " << _t_output_io->get_accumulated_time() / 1000. << " sec" << std::endl;
280  std::cout << "=======================================" << std::endl;
281 }
282 
284 {
285  if (_do_evt_display)
286  _fitter->displayEvent();
287 
288  delete _t_seeds_cleanup;
290  delete _t_translate1;
291  delete _t_translate2;
292  delete _t_translate3;
293  delete _t_kalman_pat_rec;
294  delete _t_full_fitting;
295  delete _t_search_clusters;
296  delete _t_track_propagation;
297  delete _t_output_io;
298 
299 #ifdef _DEBUG_
300  LogDebug("Leaving End \n");
301 #endif
302 
303 #ifdef _DEBUG_
304  fout_kalman_pull.close();
305  fout_chi2.close();
306 #endif
307 
308  if (_analyzing_mode)
309  {
310  cout << " cleaning up " << endl;
311  _analyzing_file->cd();
312  _analyzing_ntuple->Write();
313  _analyzing_file->Close();
314  // delete _analyzing_ntuple;
315  // delete _analyzing_file;
316  }
317 
319 }
320 
322 {
323  TGeoManager* tgeo_manager = PHGeomUtility::GetTGeoManager(topNode);
324 
325  PHField* field = PHFieldUtility::GetFieldMapNode(nullptr, topNode);
326 
327  _fitter.reset( PHGenFit::Fitter::getInstance(tgeo_manager, field, _track_fitting_alg_name, "RKTrackRep", _do_evt_display) );
328 
329  if (!_fitter)
330  {
331  cerr << PHWHERE << endl;
333  }
334 
335 #ifdef _DEBUG_
336  _fitter->set_verbosity(10);
337 #endif
338 
340 }
341 
343 {
344  //---------------------------------
345  // Get Objects off of the Node Tree
346  //---------------------------------
347  _topNode = topNode;
348 
349  // used in fast vertexing from BBC
350  _bbc_vertexes = findNode::getClass<BbcVertexMap>(topNode, "BbcVertexMap");
351 
352  // _cluster_map is the TrkrClusterMap container, and is found in PHTrackPropagating
353 
355  PHG4CylinderGeomContainer>(topNode, "CYLINDERGEOM_INTT");
356 
358  PHG4CylinderGeomContainer>(topNode, "CYLINDERGEOM_MVTX");
359 
361 }
362 
363 int PHGenFitTrkProp::check_track_exists(MapPHGenFitTrack::iterator iter, SvtxTrackMap::Iter /*phtrk_iter*/)
364 {
365  //Loop over hitIDs on current track and check if they have been used
366  unsigned int n_clu = iter->second->get_cluster_keys().size();
367 
368  unsigned int n_clu_used = 0;
369 
370  const std::vector<TrkrDefs::cluskey>& clusterkeys = iter->second->get_cluster_keys();
371 
372  int n = 0;
373  for (TrkrDefs::cluskey iCluId = 0; iCluId < clusterkeys.size(); ++iCluId)
374  {
375  TrkrDefs::cluskey cluster_ID = clusterkeys[iCluId];
376 
377  if(_gftrk_hitkey_map.count(iCluId)>0) n_clu_used++;
378  if (Verbosity() >= 10){
379  cout << " trk map size: " << _gftrk_hitkey_map.count(iCluId) << endl;
380  cout << "n: " << n << "#Clu_g = " << iCluId << " ntrack match: " << _assoc_container->GetTracksFromCluster(cluster_ID).size()
381  << " layer: " << (float)TrkrDefs::getLayer(cluster_ID)
382  << " r: " << sqrt(_cluster_map->findCluster(cluster_ID)->getX()*_cluster_map->findCluster(cluster_ID)->getX() +_cluster_map->findCluster(cluster_ID)->getY()*_cluster_map->findCluster(cluster_ID)->getY() )
383  << " used: " << n_clu_used
384  << endl;
385  n++;
386  }
387  }
388 
389  int code = 0;
390  if (((float) n_clu_used / n_clu) > 0.3)
391  {
392  if (Verbosity() >= 1)
393  cout << "Found duplicate track. n_clu: " << n_clu << " c_clu_used: " << n_clu_used << endl;
394  /*
395  for(TrkrDefs::cluskey iCluId = 0; iCluId < clusterkeys.size(); ++iCluId){
396  TrkrDefs::cluskey cluster_ID = clusterkeys[iCluId];
397 
398  }
399  */
400 
401  return code;
402  }
403  code = 1;
404  return code;
405 }
406 
408 {
409  // _track_map is created in PHHoughSeeding::ExportOutput at line 1714 and written to SvtxTrack node.
410  // PHTrackPropagating (base class of this one) gets it from the node tree and calls this "Process" method.
411  // The process method then calls this method.
412 
413 #ifdef _DEBUG_
414  std::cout << "=========================" << std::endl;
415  std::cout << "PHGenFitTrkProp::KalmanTrkProp: Start: Event: " << _event << std::endl;
416  std::cout << "Total Raw Tracks: " << _track_map->size() << std::endl;
417  std::cout << "=========================" << std::endl;
418 #endif
419 
423  // A map is needed for each vertex location. This is a vector of maps
425  /* for(unsigned int ivert=0; ivert<_vertex.size(); ++ivert)
426  {
427  BuildLayerZPhiHitMap(ivert);
428  }
429  */
431 
432  vector<genfit::Track*> evt_disp_copy;
433 
434  _PHGenFitTracks.clear();
435 
436  //_track_map->identify();
437 
438  if (Verbosity() > 1){
439  cout << " found " << _track_map->size() << " track seeds " << endl;
440  }
441 
442  for (auto phtrk_iter = _track_map->begin();
443  phtrk_iter != _track_map->end();)
444  {
445  SvtxTrack* tracklet = phtrk_iter->second;
446  if (Verbosity() >= 1){
447  std::cout
448  << __LINE__
449  << ": Processing seed itrack: " << phtrk_iter->first
450  << ": nhits: " << tracklet-> size_cluster_keys()
451  << ": Total tracks: " << _track_map->size()
452  << endl;
453  }
454  _ntrack = phtrk_iter->first;
459  /*
460  std::shared_ptr<PHGenFit::Track> rf_phgf_track = ReFitTrack(tracklet);
461  if(!rf_phgf_track)
462  {
463  cout << "FAIL ReFitting input track# " << phtrk_iter->first << endl;
464  }
465  */
466  //if (Verbosity() > 1) _t_translate_to_PHGenFitTrack->stop();
467 
471  bool is_splitting_track = false;
472 #ifdef _DEBUG_
473  int i = 0;
474 #endif
475  SvtxTrackToPHGenFitTracks(tracklet);
476  if (_PHGenFitTracks.empty())
477  {
478  cout << "FAIL SvtxTrackToGenFit " << phtrk_iter->first << endl;
479  //cout << "Warning: Conversion of SvtxTrack tracklet " << phtrk_iter->first << " to PHGenFitTrack failed, moving to next tracklet " << endl;
480  ++phtrk_iter;
481  continue;
482  }
483 
484  for (auto gftrk_iter = _PHGenFitTracks.begin();
485  gftrk_iter != _PHGenFitTracks.end(); ++gftrk_iter)
486  {
487  if(Verbosity() > 10)
488  {
489  cout
490  << __LINE__
491  << ": propagating Genfit track: " << phtrk_iter->first << endl;
492  }
493 
494  // associate this track with the same vertex as the seed track
495  // unsigned int ivert = tracklet->get_vertex_id();
496  unsigned int ivert = 0;
497  gftrk_iter->second->set_vertex_id(ivert);
498  //cout << PHWHERE << " read back track vertex ID from Genfit track " << gftrk_iter->second->get_vertex_id() << endl;
499  if(ivert > _vertex.size())
500  {
501  cout << PHWHERE << " Track vertex is screwed up, have to quit! #" << ivert << " v_size: " << _vertex.size() << endl;
503  }
504 
505  std::vector<TrkrDefs::cluskey> clusterkeys = gftrk_iter->second->get_cluster_keys();
506 
507  unsigned int init_layer = UINT_MAX;
508  unsigned int min_layer = UINT_MAX;
509  unsigned int max_layer = 0;
510 
511  for (auto iter = clusterkeys.begin();
512  iter != clusterkeys.end(); ++iter){
513  TrkrDefs::cluskey cluster_key = *iter;
514  if(cluster_key!=0){
515  unsigned int layer = TrkrDefs::getLayer(cluster_key);
516  if(layer < min_layer) min_layer = layer;
517  if(layer > max_layer) max_layer = layer;
518  }
519  }
520  if (!is_splitting_track)
521  {
522  if (_init_direction == 1)
523  {
524  //init_layer = _cluster_map->get(clusterkeys.front())->get_layer();
525  init_layer = min_layer;//TrkrDefs::getLayer(clusterkeys.front());
526  TrackPropPatRec(ivert, gftrk_iter, init_layer, _nlayers_all, true);
527  TrackPropPatRec(ivert, gftrk_iter, init_layer, 0, false);
528  }
529  else
530  {
531  //init_layer = _cluster_map->get(clusterkeys.back())->get_layer();
532  init_layer = max_layer;//TrkrDefs::getLayer(clusterkeys.back());
533  TrackPropPatRec(ivert, gftrk_iter, init_layer, 0, true);
534  TrackPropPatRec(ivert, gftrk_iter, init_layer, _nlayers_all, false);
535  }
536  is_splitting_track = true;
537  }
538  else
539  {
540  if (_init_direction == 1)
541  {
542  //init_layer = _cluster_map->get(clusterkeys.front())->get_layer();
543  init_layer = min_layer;//TrkrDefs::getLayer(clusterkeys.front());
544  TrackPropPatRec(ivert, gftrk_iter, init_layer, _nlayers_all, false);
545  }
546  else
547  {
548  //init_layer = _cluster_map->get(clusterkeys.back())->get_layer();
549  init_layer = max_layer;//TrkrDefs::getLayer(clusterkeys.back());
550  TrackPropPatRec(ivert, gftrk_iter, init_layer, 0, false);
551  }
552  }
553 
554 #ifdef _DEBUG_
555  cout
556  << __LINE__
557  << ": tracki: " << i
558  << ": clusterkeys size: " << gftrk_iter->second->get_cluster_keys().size()
559  << ": quality: " << gftrk_iter->first
560  << endl;
561  ++i;
562 #endif
563 
564  //_trackID_PHGenFitTrack.erase(iter);
565  } // loop _PHGenFitTracks
566 
567 
568 #ifdef _DEBUG_
569  i = 0;
570  for (auto iter = _PHGenFitTracks.begin();
571  iter != _PHGenFitTracks.end(); ++iter)
572  {
573  cout
574  << __LINE__
575  << ": track: " << i++
576  << ": clusterkeys size: " << iter->second->get_cluster_keys().size()
577  << ": quality: " << iter->first
578  << endl;
579  }
580 #endif
581 
582  //std::sort(_PHGenFitTracks.begin(), _PHGenFitTracks.end());
583  _PHGenFitTracks.sort();
584 
585 #ifdef _DEBUG_
586  for (auto iter = _PHGenFitTracks.begin();
587  iter != _PHGenFitTracks.end(); ++iter)
588  {
589  cout
590  << __LINE__
591  << ": clusterkeys size: " << iter->second->get_cluster_keys().size()
592  << ": quality: " << iter->first
593  << endl;
594  }
595 #endif
596 
597  auto gftrk_iter_best = _PHGenFitTracks.begin();
598 
599  int track_exists = check_track_exists(gftrk_iter_best,phtrk_iter);
600 
601  if (gftrk_iter_best->second->get_cluster_keys().size() >= _min_good_track_hits && track_exists)
602  {
603  OutputPHGenFitTrack(gftrk_iter_best, phtrk_iter);
604 #ifdef _DEBUG_
605  cout << __LINE__ << endl;
606 #endif
607  if (_do_evt_display)
608  {
609  evt_disp_copy.push_back(
610  new genfit::Track(*gftrk_iter_best->second->getGenFitTrack()));
611  }
612  }
613  else
614  {
615  auto key = phtrk_iter->first;
616  ++phtrk_iter;
617  _track_map->erase(key);
618  continue;
619  }
620 
621  _PHGenFitTracks.clear();
622 
623  ++phtrk_iter;
624  }
625 
626  if(Verbosity() > 1)
627  {
628  std::cout << "=========================" << std::endl;
629  std::cout << "PHGenFitTrkProp::KalmanTrkProp: End: Event: " << _event << std::endl;
630  std::cout << "PHGenFitTrkProp::KalmanTrkProp: End: Event: " << _event << std::endl;
631  std::cout << "Total Final Tracks: " << _track_map->size() << std::endl;
632  }
633 
634  if (_do_evt_display)
635  {
636  _fitter->getEventDisplay()->addEvent(evt_disp_copy);
637  }
638  else
639  {
640  evt_disp_copy.clear();
641  }
642 
644 }
645 
647  MapPHGenFitTrack::iterator gftrk_iter,
648  SvtxTrackMap::Iter phtrk_iter)
649 {
650  if(Verbosity() > 10)
651  {
652  std::cout << "=========================" << std::endl;
653  std::cout << __LINE__ << ": iPHGenFitTrack: " << phtrk_iter->first << std::endl;
654  std::cout << __LINE__ << ": _track_map->size(): " << _track_map->size() << std::endl;
655  std::cout << "Contains: " << gftrk_iter->second->get_cluster_keys().size() << " clusters." << std::endl;
656  std::cout << "=========================" << std::endl;
657  }
658 
659  // get the track from the node tree and extract the track and vertex id
660  SvtxTrack* track = phtrk_iter->second;
661  auto track_id = track->get_id();
662  auto vertex_id = track->get_vertex_id();
663  auto track_charge = track->get_charge();
664 
665  // now reset the track info and rewrite it using info from the Genfit track
666  track->Reset();
667  track->set_id(track_id);
668  track->set_vertex_id(vertex_id);
669  track->set_charge(track_charge);
670 
671 #ifdef _DO_FULL_FITTING_
672  if (Verbosity() >= 1) _t_full_fitting->restart();
673  if (_fitter->processTrack(gftrk_iter->second.get(), false) != 0)
674  {
675  if (Verbosity() >= 1)
676  LogWarning("Track fitting failed\n");
677  //delete track;
678  return -1;
679  }
680  if (Verbosity() >= 1) _t_full_fitting->stop();
681 
682  if (Verbosity() >= 1) _t_output_io->restart();
683  // iter->second->getGenFitTrack()->Print();
684 
685  track.set_chisq(gftrk_iter->second->get_chi2());
686  track.set_ndf(gftrk_iter->second->get_ndf());
687 
688  // Use fitted vertex
689  TVector3 vertex_position(_vertex[vertex_id][0], _vertex[vertex_id][1], _vertex[vertex_id][2]);
690  std::unique_ptr<genfit::MeasuredStateOnPlane> gf_state_vertex_ca = nullptr;
691  try
692  {
693  gf_state_vertex_ca.reset(gftrk_iter->second->extrapolateToPoint(vertex_position));
694  }
695  catch (...)
696  {
697  if (Verbosity() >= 2)
698  LogWarning("extrapolateToPoint failed!");
699  }
700  if (!gf_state_vertex_ca)
701  {
702  //delete out_track;
703  return -1;
704  }
705 
706  TVector3 mom = gf_state_vertex_ca->getMom();
707  TVector3 pos = gf_state_vertex_ca->getPos();
708  TMatrixDSym cov = gf_state_vertex_ca->get6DCov();
709 
710 #else
711  TVectorD state = gftrk_iter->second->getGenFitTrack()->getStateSeed();
712  TVector3 pos(state(0), state(1), state(2));
713  TVector3 mom(state(3), state(4), state(5));
714  TMatrixDSym cov = gftrk_iter->second->getGenFitTrack()->getCovSeed();
715 #endif
716 
717  for(int i=0; i<6; i++){
718  for(int j=0; j<6; j++){
719  track->set_error(i, j, cov[i][j]);
720  }
721  }
722 
723  track->set_px(mom.Px());
724  track->set_py(mom.Py());
725  track->set_pz(mom.Pz());
726 
727  track->set_x(pos.X());
728  track->set_y(pos.Y());
729  track->set_z(pos.Z());
730 
731  for (TrkrDefs::cluskey cluster_key : gftrk_iter->second->get_cluster_keys())
732  {
733  if(cluster_key<=0)continue;
734  if(Verbosity() > 10) cout << " track id: " << phtrk_iter->first << " adding clusterkey " << cluster_key << endl;
735  _gftrk_hitkey_map.insert(std::make_pair(cluster_key, phtrk_iter->first));
736  track->insert_cluster_key(cluster_key);
737  }
738 
739  //Check track quality
740  // bool is_good_track = true;
741 
742  Int_t n_maps = 0;
743  Int_t n_intt = 0;
744  Int_t n_tpc = 0;
745  Int_t n_micromegas = 0;
746 
747  for (auto iter = track->begin_cluster_keys(); iter != track->end_cluster_keys(); ++iter)
748  {
749  TrkrDefs::cluskey cluster_key = *iter;
750  unsigned int layer = TrkrDefs::getLayer(cluster_key);
751  if( is_maps_layer( layer ) )
752  {
753  n_maps++;
754  } else if( is_intt_layer( layer ) ) {
755  n_intt++;
756  if (n_intt > 4)
757  {
758  cout << PHWHERE << " Can not have more than 4 INTT layers, quit!" << endl;
759  exit(1);
760  }
761  } else if( is_tpc_layer( layer ) ) {
762  n_tpc++;
763  } else if( is_micromegas_layer( layer ) ) {
764  n_micromegas++;
765  }
766  }
767 
768  // Add the cluster-track association to the association table for later use
769  for (TrkrDefs::cluskey cluster_key : gftrk_iter->second->get_cluster_keys())
770  {
771  //cout << "Add cluster key " << cluster_key << " to ClusterTrackAssoc " << track->get_id() << endl;
772  _assoc_container->SetClusterTrackAssoc(cluster_key, track->get_id());
773  }
774 
775  if (Verbosity() > 10)
776  {
777  cout << "Output propagared track " << track->get_id() << " vertex " << track->get_vertex_id()<< " quality = " << track->get_quality()
778  << " clusters mvtx: " << n_maps << " intt: " << n_intt << " tpc: " << n_tpc << " micromegas: " << n_micromegas
779  << endl;
780  cout << "px = " << track->get_px() << " py = " << track->get_py()
781  << " pz = " << track->get_pz() << endl;
782  }
783 
784  if (Verbosity() >= 1) _t_output_io->stop();
785 
786  return 0;
787 }
788 
789 std::shared_ptr<PHGenFit::Track> PHGenFitTrkProp::ReFitTrack(SvtxTrack* intrack)
790 {
791  if (!intrack)
792  {
793  cerr << PHWHERE << " Input SvtxTrack is nullptr!" << endl;
794  return nullptr;
795  }
796 
797  PHG4CylinderGeomContainer* geom_container_intt = findNode::getClass<
798  PHG4CylinderGeomContainer>(_topNode, "CYLINDERGEOM_INTT");
799 
800  PHG4CylinderGeomContainer* geom_container_mvtx = findNode::getClass<
801  PHG4CylinderGeomContainer>(_topNode, "CYLINDERGEOM_MVTX");
802 
803  // prepare seed
804  TVector3 seed_mom(100, 0, 0);
805  TVector3 seed_pos(0, 0, 0);
806  TMatrixDSym seed_cov(6);
807  for (int i = 0; i < 6; i++)
808  {
809  for (int j = 0; j < 6; j++)
810  {
811  seed_cov[i][j] = 100.;
812  }
813  }
814 
815  //Sort hits in r
816  std::map<float, TrkrDefs::cluskey> m_r_cluster_id;
817  for (auto iter = intrack->begin_cluster_keys();
818  iter != intrack->end_cluster_keys(); ++iter){
819  TrkrDefs::cluskey cluster_key = *iter;
820  TrkrCluster* cluster = _cluster_map->findCluster(cluster_key);
821  float x = cluster->getPosition(0);
822  float y = cluster->getPosition(1);
823  float r = sqrt(x*x + y*y);
824  m_r_cluster_id.insert(std::pair<float, TrkrDefs::cluskey>(r, cluster_key));
825  int layer_out = TrkrDefs::getLayer(cluster_key);
826  if(Verbosity() >= 2) cout << " Layer " << layer_out << " cluster " << cluster_key << " radius " << r << endl;
827  }
828 
829  // Create measurements
830  std::vector<PHGenFit::Measurement*> measurements;
831  for (auto iter = m_r_cluster_id.begin();
832  iter != m_r_cluster_id.end();
833  ++iter)
834  {
835  TrkrDefs::cluskey cluster_key = iter->second;
836  const int layer = TrkrDefs::getLayer(cluster_key);
837 
838  TrkrCluster* cluster = _cluster_map->findCluster(cluster_key);
839  if (!cluster){
840  LogError("No cluster Found!");
841  continue;
842  }
843  // PHGenFit::Measurement* meas = TrkrClusterToPHGenFitMeasurement(cluster);
844 
845  TVector3 pos(cluster->getPosition(0), cluster->getPosition(1), cluster->getPosition(2));
846 
847  seed_mom.SetPhi(pos.Phi());
848  seed_mom.SetTheta(pos.Theta());
849 
850  //TODO use u, v explicitly?
851  TVector3 n(cluster->getPosition(0), cluster->getPosition(1), 0);
852 
853  // get the trkrid
854  unsigned int trkrid = TrkrDefs::getTrkrId(cluster_key);
855 
856  if(trkrid == TrkrDefs::mvtxId)
857  {
858  int stave_index = MvtxDefs::getStaveId(cluster_key);
859  int chip_index = MvtxDefs::getChipId(cluster_key);
860 
861  double ladder_location[3] = {0.0, 0.0, 0.0};
862  CylinderGeom_Mvtx* geom =
863  dynamic_cast<CylinderGeom_Mvtx*>(geom_container_mvtx->GetLayerGeom(layer));
864  // returns the center of the sensor in world coordinates - used to get the ladder phi location
865  geom->find_sensor_center(stave_index, 0,
866  0, chip_index, ladder_location);
867 
868  //cout << " MVTX stave phi tilt = " << geom->get_stave_phi_tilt()
869  // << " seg.X " << ladder_location[0] << " seg.Y " << ladder_location[1] << " seg.Z " << ladder_location[2] << endl;
870  n.SetXYZ(ladder_location[0], ladder_location[1], 0);
871  n.RotateZ(geom->get_stave_phi_tilt());
872  }
873  else if(trkrid == TrkrDefs::inttId)
874  {
875  CylinderGeomIntt* geom =
876  dynamic_cast<CylinderGeomIntt*>(geom_container_intt->GetLayerGeom(layer));
877  double hit_location[3] = {0.0, 0.0, 0.0};
878  geom->find_segment_center(InttDefs::getLadderZId(cluster_key),
879  InttDefs::getLadderPhiId(cluster_key), hit_location);
880 
881  //cout << " Intt strip phi tilt = " << geom->get_strip_phi_tilt()
882  // << " seg.X " << hit_location[0] << " seg.Y " << hit_location[1] << " seg.Z " << hit_location[2] << endl;
883  n.SetXYZ(hit_location[0], hit_location[1], 0);
884  n.RotateZ(geom->get_strip_phi_tilt());
885  }
886  // end new
887  //-----------------
888 
890  cluster->getRPhiError(), cluster->getZError());
891  cout << "Add meas layer " << layer << " cluskey " << cluster_key
892  << endl
893  << " pos.X " << pos.X() << " pos.Y " << pos.Y() << " pos.Z " << pos.Z()
894  << " RPhiErr " << cluster->getRPhiError()
895  << " ZErr " << cluster->getZError()
896  << endl;
897  /*
898  if(Verbosity() > 10)
899  {
900  cout << "Add meas layer " << layer << " cluskey " << cluster_key
901  << endl
902  << " pos.X " << pos.X() << " pos.Y " << pos.Y() << " pos.Z " << pos.Z()
903  << " n.X " << n.X() << " n.Y " << n.Y()
904  << " RPhiErr " << cluster->getRPhiError()
905  << " ZErr " << cluster->getZError()
906  << endl;
907  }
908  */
909  measurements.push_back(meas);
910  }
911 
920  //TODO Add multiple TrackRep choices.
921  //int pid = 211;
922 
923  genfit::AbsTrackRep* rep = new genfit::RKTrackRep(_primary_pid_guess);
924  std::shared_ptr<PHGenFit::Track> track(new PHGenFit::Track(rep, seed_pos, seed_mom,
925  seed_cov));
926 
927  //TODO unsorted measurements, should use sorted ones?
928  track->addMeasurements(measurements);
929  if (Verbosity() >= 1)
930  cout << "#trk: " << _ntrack << " ReFit nhits: " << measurements.size() << endl;
936  if (_fitter->processTrack(track.get(), false) != 0)
937  {
938  if (Verbosity() >= 1)
939  {
940  cout << "#trk: " << _ntrack << "Track fitting failed (ReFit) nhits: " << measurements.size() << endl;
941  // LogWarning("Track fitting failed ");
942  cout << " " << endl;
943  /*
944  cout << " track->getChisq() " << track->get_chi2() << " get_ndf " << track->get_ndf()
945  << " mom.X " << track->get_mom().X()
946  << " mom.Y " << track->get_mom().Y()
947  << " mom.Z " << track->get_mom().Z()
948  << endl;
949  */
950  }
951  //delete track;
952  return nullptr;
953  }
954  return track;
955 }
956 
958 {
959  // clean up working array for each event
960  _PHGenFitTracks.clear();
961 
962  if(Verbosity() > 10) cout << PHWHERE << "Converting SvtxTrack to PHGenFit track: track ID " << svtxtrack->get_id()
963  << " track z " << svtxtrack->get_z()
964  << " vertex ID " << svtxtrack->get_vertex_id()
965  << " vertex z " << _vertex[svtxtrack->get_vertex_id()][2]
966  << endl;
967  //svtxtrack->identify();
968 
969  std::vector<PHGenFit::Measurement*> measurements;
970 
971  if(_use_beam){
972  TVector3 posvtx(0.0, 0.0,0.0);
973  TVector3 resvtx(0.015, 0.015,30.0);
975  beam->set_cluster_key(0);
976  measurements.push_back(beam);
977  }
978 
979  //Sort hits in r
980  // std::multimap<float, TrkrDefs::cluskey> m_r_clusterID;
981  std::map<float, TrkrDefs::cluskey> m_r_clusterID;
982 
983  for (auto hit_iter = svtxtrack->begin_cluster_keys();
984  hit_iter != svtxtrack->end_cluster_keys(); ++hit_iter)
985  {
986  TrkrDefs::cluskey clusterkey = *hit_iter;
987  TrkrCluster *cluster = _cluster_map->findCluster(clusterkey);
988  float r = sqrt(
989  cluster->getPosition(0) * cluster->getPosition(0) +
990  cluster->getPosition(1) * cluster->getPosition(1));
991  m_r_clusterID.insert(std::pair<float, TrkrDefs::cluskey>(r, clusterkey));
992  if (Verbosity() >= 10){
993  cout << PHWHERE << " inserted r " << r << " clusterkey " << clusterkey
994  << " layer: " << (float)TrkrDefs::getLayer(clusterkey)
995  << " max: " << svtxtrack->size_cluster_keys()
996  << endl;
997  }
998  }
999 
1000  /*
1001  if (_vertex_map)
1002  {
1003  TVector3 v(_vertex[0], _vertex[1], _vertex[2]);
1004  TMatrixDSym cov(3);
1005  cov.Zero();
1006  cov(0, 0) = _vertex_error[0] * _vertex_error[0];
1007  cov(1, 1) = _vertex_error[1] * _vertex_error[1];
1008  cov(2, 2) = _vertex_error[2] * _vertex_error[2];
1009  PHGenFit::Measurement* meas = new PHGenFit::SpacepointMeasurement(v, cov);
1010  //FIXME re-use the first cluster id
1011  // TrkrDefs::cluskey id = m_r_clusterID.begin()->second;
1012  // meas->set_cluster_key(id);
1013  meas->set_cluster_key(0);
1014  measurements.push_back(meas);
1015  }
1016  */
1017  TrkrDefs::cluskey last_cluster_key = 0;
1018 
1019  for (auto iter = m_r_clusterID.begin();
1020  iter != m_r_clusterID.end();
1021  ++iter)
1022  {
1023  TrkrDefs::cluskey cluster_key = iter->second;
1024  TrkrCluster *cluster = _cluster_map->findCluster(cluster_key);
1025  if (!cluster){
1026  LogError("No cluster Found!\n");
1027  continue;
1028  }
1029  // const int layer = TrkrDefs::getLayer(cluster_key);
1030  last_cluster_key = cluster_key;
1032  if (meas){
1033  measurements.push_back(meas);
1034  }
1035  }
1036 
1037  TVector3 seed_pos;
1038  TVector3 seed_mom;
1039  const float blowup_factor = 100.;
1040  TMatrixDSym seed_cov(6);
1041 
1042  if(svtxtrack->get_py()==0&&svtxtrack->get_pz()==0){
1043  TrkrCluster *last_cluster = _cluster_map->findCluster(last_cluster_key);
1044  TVector3 pos(last_cluster->getPosition(0), last_cluster->getPosition(1), last_cluster->getPosition(2));
1045  seed_mom.SetXYZ(100,0,0);
1046  seed_mom.SetPhi(pos.Phi());
1047  seed_mom.SetTheta(pos.Theta());
1048  // seed_mom.SetPtThetaPhi(pos.Perp(),pos.Theta(),pos.Phi());
1049  for (int i = 0; i < 6; i++){
1050  for (int j = 0; j < 6; j++){
1051  seed_cov[i][j] = 100.;
1052  }
1053  }
1054  }else{
1055  seed_pos.SetXYZ(svtxtrack->get_x(),svtxtrack->get_y(),svtxtrack->get_z());
1056  seed_mom.SetXYZ(svtxtrack->get_px(),svtxtrack->get_py(),svtxtrack->get_pz());
1057 
1058  for (int i = 0; i < 6; i++){
1059  for (int j = 0; j < 6; j++){
1060  seed_cov[i][j] = blowup_factor * svtxtrack->get_error(i, j);
1061  }
1062  }
1063  }
1064 
1065  genfit::AbsTrackRep* rep = new genfit::RKTrackRep(_primary_pid_guess);
1066  std::shared_ptr<PHGenFit::Track> track(
1067  new PHGenFit::Track(rep, seed_pos, seed_mom, seed_cov));
1068 
1069 
1070  track->addMeasurements(measurements);
1071  if (Verbosity() >= 1)
1072  cout << "#trk: " << _ntrack << " ToPHGen size:" << measurements.size() <<endl;
1073 
1074  if (_fitter->processTrack(track.get(), false) != 0)
1075  {
1076  if (Verbosity() >= 1){
1077  cout << "#trk: " << _ntrack << " Seed fitting failed (ToPHGen) size:" << measurements.size() <<endl;
1078  }
1079  return -1;
1080  }
1081 
1082  int nhits = track->get_cluster_keys().size();
1083  float chi2 = track->get_chi2();
1084  float ndf = track->get_ndf();
1085 
1086  if (nhits > 0 and chi2 > 0 and ndf > 0)
1087  {
1088  _PHGenFitTracks.push_back(
1089  MapPHGenFitTrack::value_type(
1090  TrackQuality(nhits, chi2, ndf, 0, nhits, 0, 0), track));
1091  }
1092 
1094 }
1095 
1097  const unsigned int ivert,
1098  MapPHGenFitTrack::iterator& track_iter,
1099  unsigned int init_layer, unsigned int end_layer,
1100  const bool use_fitted_state_once)
1101 {
1102 #ifdef _DEBUG_
1103  cout
1104  << __LINE__
1105  << " TrackPropPatRec"
1106  << " : init_layer: " << init_layer
1107  << " : end_layer: " << end_layer
1108  << " : use_fitted_state_once: " << use_fitted_state_once
1109  << endl;
1110  #endif
1111 
1112  std::shared_ptr<PHGenFit::Track>& track = track_iter->second;
1113  if(Verbosity() > 10) cout << endl << PHWHERE << " Entering TrackPropPatRec for track with vertex ID " << ivert << endl;
1114  if(ivert >= _vertex.size())
1115  {
1116  cout << PHWHERE << " WARNING: vertex ID out of range, something wrong, quit! " << endl;
1118  }
1119 
1120  int direction = end_layer >= init_layer ? 1 : -1;
1121  assert(direction == 1 or direction == -1);
1122 
1123  int first_extrapolate_base_TP_id = -1;
1124 
1125  bool use_fitted_state = use_fitted_state_once;
1126  float blowup_factor = use_fitted_state ? _blowup_factor : 1.;
1127  unsigned int layer_occupied[_nlayers_all];
1128  for(int i = 0;i<_nlayers_all;i++) layer_occupied[i] = 0;
1135  std::vector<TrkrDefs::cluskey> clusterkeys = track->get_cluster_keys();
1136 
1137  for (unsigned int i = 0; i < clusterkeys.size(); ++i)
1138  {
1139  if(clusterkeys[i]==0)continue;
1140  unsigned int layer = TrkrDefs::getLayer(clusterkeys[i]);
1141  layer_occupied[layer] = 1;
1142  if(layer == init_layer)
1143  {
1144  first_extrapolate_base_TP_id = i;
1145  break;
1146  }
1147  }
1148 
1149 
1150  if (first_extrapolate_base_TP_id < 0)
1151  {
1152  if (Verbosity() > 0)
1153  LogError("first_extrapolate_base_TP_id < 0");
1154  return -1;
1155  }
1156 
1157  int extrapolate_base_TP_id = first_extrapolate_base_TP_id;
1158 
1159  unsigned int consecutive_missing_layer = 0;
1160 
1161  // unsigned int layer = init_layer + direction;
1162  // while (layer>=0 and layer < (unsigned int)_nlayers_all and layer!=end_layer) {
1163 
1164  for (unsigned int layer = init_layer + direction;
1165  layer != end_layer + direction;
1166  layer += direction)
1167  {
1168  // layer is unsigned int, check for >=0 is meaningless
1169  if (layer >= (unsigned int) _nlayers_all) break;
1170  if (layer_occupied[layer]) continue;
1171 
1175  if (consecutive_missing_layer > _max_consecutive_missing_layer)
1176  {
1177  if (Verbosity() > 1)
1178  {
1179  LogWarning("consecutive_missing_layer > ") << _max_consecutive_missing_layer << endl;
1180  }
1181  if (track->get_cluster_keys().size() >= _min_good_track_hits)
1182  return 0;
1183  else
1184  return -1;
1185  }
1186 
1187  bool layer_updated = false;
1188 
1189  float layer_r = _radii_all[_layer_ilayer_map_all[layer]];
1190 
1191  if(Verbosity() > 10)
1192  {
1193  std::cout << "=========================" << std::endl;
1194  std::cout << __LINE__ << ": Event: " << _event << ": _PHGenFitTracks.size(): " << _PHGenFitTracks.size() << ": layer: " << layer << std::endl;
1195  std::cout << "=========================" << std::endl;
1196  }
1197 
1198 #ifdef _DEBUG_
1199  {
1200  unsigned int tempIdx = extrapolate_base_TP_id >= 0 ? extrapolate_base_TP_id : extrapolate_base_TP_id + track->get_cluster_keys().size();
1201  cout
1202  << __LINE__
1203  << " tempIdx: " << tempIdx
1204  << endl;
1205  // tempIdx is unsigned int, checking for >=0 is meaningless
1206  if (tempIdx < track->get_cluster_keys().size())
1207  {
1208  TrkrDefs::cluskey extrapolate_base_cluster_id = track->get_cluster_keys()[tempIdx];
1209  //SvtxCluster* extrapolate_base_cluster = _cluster_map->get(extrapolate_base_cluster_id);
1210  TrkrCluster* extrapolate_base_cluster = _cluster_map->findCluster(extrapolate_base_cluster_id);
1211  int from_layer = TrkrDefs::getLayer(extrapolate_base_cluster_id);
1212  cout
1213  << __LINE__
1214  << ": Target layer: { " << layer
1215  << ", " << layer_r
1216  << "} : From layer: { " << from_layer
1217  << ", " << sqrt(extrapolate_base_cluster->getX() * extrapolate_base_cluster->getX() + extrapolate_base_cluster->getY() * extrapolate_base_cluster->getY())
1218  << "} : ID: " << extrapolate_base_cluster_id
1219  << endl;
1220  }
1221  }
1222 #endif
1223 
1224  // bool have_tp_with_fit_info = false;
1225  // std::vector<unsigned int> clusterkeys = track->get_cluster_keys();
1226  // for (unsigned int i = clusterkeys.size() - 1; i >= 0; --i) {
1227  // std::unique_ptr<genfit::MeasuredStateOnPlane> kfsop = nullptr;
1228  // genfit::Track genfit_track = track->getGenFitTrack();
1229  // if (genfit_track->getNumPointsWithMeasurement() > 0) {
1230  //
1231  // genfit::TrackPoint* tp = genfit_track->getPointWithMeasurement(tr_point_id);
1232  // if (dynamic_cast<genfit::KalmanFitterInfo*>(tp->getFitterInfo(rep))) {
1233  //
1234  // if (static_cast<genfit::KalmanFitterInfo*>(tp->getFitterInfo(rep))->getForwardUpdate()) {
1235  // have_tp_with_fit_info = true;
1236  // break;
1237  // }
1238  // }
1239  // }
1240  // }
1241 
1242  std::unique_ptr<genfit::MeasuredStateOnPlane> state = nullptr;
1243  try
1244  {
1245  state.reset(track->extrapolateToCylinder(layer_r, TVector3(0, 0, 0),
1246  TVector3(0, 0, 1), extrapolate_base_TP_id, direction));
1247  // genfit::MeasuredStateOnPlane *state = track->extrapolateToCylinder(
1248  // layer_r, TVector3(0, 0, 0), TVector3(0, 0, 1), 0);
1249  }
1250  catch (...)
1251  {
1252  if (Verbosity() > 1)
1253  {
1254  LogWarning("Can not extrapolate to Cylinder!") << std::endl;
1255  }
1256  continue;
1257  }
1258 
1259  if (!state)
1260  {
1261  if (Verbosity() > 1)
1262  {
1263  LogWarning("Can not extrapolate to Cylinder!") << std::endl;
1264  }
1265  continue;
1266  }
1267 
1268  // The vertex position is used here to caculate theta and phi for the
1269  // extrapolated track position. These values of theta and phi are compared
1270  // with the cluster theta and phi positions in the map "_layer_thetaID_phiID_clusterID"
1271 
1272  // Need to use the correct vertex position for this track as well as for the theta-phi map.
1273 
1274  TVector3 pos = state->getPos();
1275  if(Verbosity() > 10) cout << PHWHERE << " pos z " << pos.Z() << " _vertex z " << _vertex[ivert][2] << endl;
1276  pos.SetXYZ(
1277  pos.X() - _vertex[ivert][0],
1278  pos.Y() - _vertex[ivert][1],
1279  pos.Z() - _vertex[ivert][2]);
1280 
1281  float phi_center = pos.Phi();
1282  float theta_center = pos.Theta();
1283 
1284 #ifdef _USE_CONSTANT_SEARCH_WIN_
1285 
1286  float phi_window = 25e-4;
1287  float theta_window = 25e-4;
1288 
1289  if (layer >= 3 and layer <= 6)
1290  {
1291  phi_window = 300e-4;
1292  theta_window = 0.2;
1293  }
1294 
1295  if (layer <= 2)
1296  {
1297  phi_window = 3000e-4;
1298  theta_window = 3000e-4;
1299  }
1300 #else
1301  TMatrixDSym cov = state->get6DCov();
1302 
1303  float phi_window = _search_wins_phi[layer] * sqrt(cov[0][0] + cov[1][1] + cov[0][1] + cov[1][0]) / pos.Perp();
1304  float theta_window = _search_wins_theta[layer] * sqrt(cov[2][2]) / pos.Perp();
1305 
1306  if( is_maps_layer( layer ) )
1307  {
1308  if (phi_window > _max_search_win_phi_maps) phi_window = _max_search_win_phi_maps;
1309  if (phi_window < _min_search_win_phi_maps) phi_window = _min_search_win_phi_maps;
1310  if (theta_window > _max_search_win_theta_maps) theta_window = _max_search_win_theta_maps;
1311  if (theta_window < _min_search_win_theta_maps) theta_window = _min_search_win_theta_maps;
1312  }
1313  else if( is_intt_layer( layer ) )
1314  {
1315  const auto layer_intt = layer - _firstlayer_intt;
1316  if (phi_window > _max_search_win_phi_intt[layer_intt]) phi_window = _max_search_win_phi_intt[layer_intt];
1317  if (phi_window < _min_search_win_phi_intt[layer_intt]) phi_window = _min_search_win_phi_intt[layer_intt];
1318  if (theta_window > _max_search_win_theta_intt[layer_intt]) theta_window = _max_search_win_theta_intt[layer_intt];
1319  if (theta_window < _min_search_win_theta_intt[layer_intt]) theta_window = _min_search_win_theta_intt[layer_intt];
1320  }
1321  else if( is_tpc_layer( layer ) )
1322  {
1323  if (phi_window > _max_search_win_phi_tpc) phi_window = _max_search_win_phi_tpc;
1324  if (phi_window < _min_search_win_phi_tpc) phi_window = _min_search_win_phi_tpc;
1325  if (theta_window > _max_search_win_theta_tpc) theta_window = _max_search_win_theta_tpc;
1326  if (theta_window < _min_search_win_theta_tpc) theta_window = _min_search_win_theta_tpc;
1327  }
1328  else if( is_micromegas_layer( layer ) )
1329  {
1330  const auto layer_micromegas = layer - _firstlayer_micromegas;
1331  if (phi_window > _max_search_win_phi_micromegas[layer_micromegas]) phi_window = _max_search_win_phi_micromegas[layer_micromegas];
1332  if (phi_window < _min_search_win_phi_micromegas[layer_micromegas]) phi_window = _min_search_win_phi_micromegas[layer_micromegas];
1333  if (theta_window > _max_search_win_theta_micromegas[layer_micromegas]) theta_window = _max_search_win_theta_micromegas[layer_micromegas];
1334  if (theta_window < _min_search_win_theta_micromegas[layer_micromegas]) theta_window = _min_search_win_theta_micromegas[layer_micromegas];
1335  }
1336 #endif
1337 
1338 #ifdef _DEBUG_
1339  cout << __LINE__ << ": ";
1340  printf("layer: %u: r: %f: phi: %f +- %f; theta: %f +- %f\n",
1341  layer, pos.Perp(),
1342  phi_center, phi_window,
1343  theta_center, theta_window);
1344 
1345 // cout<<__LINE__<<": ";
1346 // printf("layer: %d: phi: %f +- %f\n",
1347 // layer,
1348 // pos.Phi(), phi_window
1349 // );
1350 #endif
1351 
1352  if (Verbosity() >= 1) _t_search_clusters->restart();
1353  std::vector<TrkrDefs::cluskey> new_cluster_keys = SearchHitsNearBy(ivert, layer,
1354  theta_center, phi_center, theta_window, phi_window);
1355  if (Verbosity() >= 1) _t_search_clusters->stop();
1356 
1357 #ifdef _DEBUG_
1358  cout << __LINE__ << ": new_cluster_keys size: " << new_cluster_keys.size() << std::endl;
1359 #endif
1360 
1361  std::vector<PHGenFit::Measurement*> measurements;
1362  for (TrkrDefs::cluskey cluster_key : new_cluster_keys)
1363  {
1364  TrkrCluster* cluster = _cluster_map->findCluster(cluster_key);
1365  if (!cluster)
1366  {
1367  LogError("No cluster Found!\n");
1368  continue;
1369  }
1370 
1372 
1373  if (meas)
1374  measurements.push_back(meas);
1375  }
1376  std::map<double, shared_ptr<PHGenFit::Track> > incr_chi2s_new_tracks;
1377 
1378 #ifdef _DEBUG_
1379  cout << __LINE__ << ": measurements.size(): " << measurements.size() << endl;
1380 #endif
1381 
1382  if (Verbosity() >= 1) _t_track_propagation->restart();
1383  track->updateOneMeasurementKalman(measurements, incr_chi2s_new_tracks, extrapolate_base_TP_id, direction, blowup_factor, use_fitted_state);
1384  use_fitted_state = false;
1385  blowup_factor = 1.;
1386  if (Verbosity() >= 1) _t_track_propagation->stop();
1387 
1388 #ifdef _DEBUG_
1389  cout << __LINE__ << ": incr_chi2s_new_tracks.size(): " << incr_chi2s_new_tracks.size() << endl;
1390 #endif
1391 
1392  PHGenFitTrkProp::TrackQuality tq(track_iter->first);
1393 
1394  // Update first track candidate
1395  if (incr_chi2s_new_tracks.size() > 0)
1396  {
1397  auto iter = incr_chi2s_new_tracks.begin();
1398 
1399  if (iter->first < _max_incr_chi2s[layer] and iter->first > 0)
1400  {
1401 #ifdef _DEBUG_
1402  cout
1403  << __LINE__
1404  << ": iPHGenFitTrack: " << iPHGenFitTrack << endl
1405  << ": First accepted IncrChi2: " << iter->first << endl
1406  << "; before update: " << track->get_cluster_keys().back()
1407  << endl;
1408 #endif
1409  // _PHGenFitTracks[iPHGenFitTrack] = std::shared_ptr
1410  // < PHGenFit::Track > (iter->second);
1411  // track = _PHGenFitTracks[iPHGenFitTrack];
1412 
1413  // track_iter->first += iter->first;
1414 
1415  track_iter->first.nhits = tq.nhits + 1;
1416  track_iter->first.chi2 = tq.chi2 + iter->first;
1417  track_iter->first.ndf = tq.ndf + 2;
1418  track_iter->first.nmicromegas = tq.nmicromegas + (is_micromegas_layer( layer ) ? 1 : 0);
1419  track_iter->first.ntpc = tq.ntpc + (is_tpc_layer(layer) ? 1 : 0);
1420  track_iter->first.nintt = tq.nintt + (is_intt_layer(layer) ? 1 : 0);
1421  track_iter->first.nmaps = tq.nmaps + (is_maps_layer(layer) ? 1 : 0);
1422 
1423  track_iter->second = std::shared_ptr<PHGenFit::Track>(iter->second);
1424 
1425  consecutive_missing_layer = 0;
1426  layer_updated = true;
1427  extrapolate_base_TP_id = -1;
1428 #ifdef _DEBUG_
1429  cout
1430  << __LINE__
1431  << ": after update: " << track->get_cluster_keys().back()
1432  << endl;
1433 
1434  fout_chi2
1435  << _event << "\t"
1436  << iPHGenFitTrack << "\t"
1437  << layer << "\t "
1438  << iter->first
1439  << endl;
1440 #endif
1441  }
1442  }
1443 
1444  // Update other candidates
1445  if (incr_chi2s_new_tracks.size() > 1 and
1446  (layer >= _nlayers_maps and layer < _nlayers_maps + _nlayers_intt))
1447  {
1448  for (auto iter = (++incr_chi2s_new_tracks.begin());
1449  iter != incr_chi2s_new_tracks.end(); ++iter)
1450  {
1451  if (!(iter->first < _max_splitting_chi2 and iter->first > 0))
1452  break;
1453 
1454 #ifdef _DEBUG_
1455  std::cout << __LINE__ << ": "
1456  << "Track Spliting with "
1457  << "IncrChi2: " << iter->first << std::endl;
1458 #endif
1459 
1460  // _PHGenFitTracks.insert(
1461  // MapPHGenFitTrack::value_type(track_iter->first + iter->first,
1462  // std::shared_ptr < PHGenFit::Track> (iter->second)));
1463 
1464  // this is a new track, have to associate it with the vertex
1465  iter->second->set_vertex_id(ivert);
1466  _PHGenFitTracks.push_back(
1467  MapPHGenFitTrack::value_type(
1468  TrackQuality(
1469  tq.nhits + 1,
1470  tq.chi2 + iter->first,
1471  tq.ndf + 2,
1472  tq.nmicromegas + (is_micromegas_layer( layer ) ? 1 : 0),
1473  tq.ntpc + (is_tpc_layer(layer) ? 1 : 0),
1474  tq.nintt + (is_intt_layer(layer) ? 1 : 0),
1475  tq.nmaps + (is_maps_layer(layer) ? 1 : 0)),
1476  iter->second));
1477  }
1478 
1479 #ifdef _DEBUG_
1480  std::cout << __LINE__ << ": "
1481  << "_PHGenFitTracksSize: " << _PHGenFitTracks.size() << std::endl;
1482  std::cout << __LINE__ << ": " << track_iter->second->get_cluster_keys().back() << std::endl;
1483 #endif
1484  }
1485 
1486 #ifdef _DEBUG_
1487  cout<<__LINE__<<": updateOneMeasurementKalman:"<<endl;
1488  std::cout<<"iPHGenFitTrack: "<<iPHGenFitTrack
1489  <<", layer: "<<layer
1490  <<", #meas: "<<measurements.size()
1491  <<", #tracks: "<<incr_chi2s_new_tracks.size()
1492  //<<", #totoal tracks: "<<_trackID_PHGenFitTrack.size()
1493  <<std::endl;
1494 
1495  for (auto iter =
1496  incr_chi2s_new_tracks.begin();
1497  iter != incr_chi2s_new_tracks.end(); iter++)
1498  {
1499  std::cout << __LINE__ << ": IncrChi2: " << iter->first << std::endl;
1500  }
1501 #endif
1502  // if (_analyzing_mode)
1503  // {
1504  // int ncand = 0;
1505  // for (auto iter =
1506  // incr_chi2s_new_tracks.begin();
1507  // iter != incr_chi2s_new_tracks.end(); iter++)
1508  // {
1509  // if (iter->first < _max_incr_chi2s[layer] and iter->first > 0) ncand++;
1510  // }
1511  // /*
1512  // float this_pt = 0.0;//track->getGenFitTrack()->getCardinalRep()->getMom(state).Pt();
1513  // float this_phi = 0.0;//track->getGenFitTrack()->getCardinalRep()->getMom(state).Phi();
1514  // float this_eta = 0.0;//track->getGenFitTrack()->getCardinalRep()->getMom(state).Eta();
1515  // */
1516  // //"spt:seta:sphi:pt:eta:phi:layer:ncand:nmeas"
1517  // //_analyzing_ntuple->Fill(init_pt,init_eta,init_phi,this_pt,this_eta,this_phi,layer,ncand,measurements.size());
1518  // }
1519  if (!layer_updated)
1520  ++consecutive_missing_layer;
1521  } // layer loop
1522 
1523 #ifdef _DEBUG_
1524  cout
1525  << __LINE__
1526  << ": cluster keys size: " << track->get_cluster_keys().size()
1527  << endl;
1528 #endif
1529 
1531  return 0;
1532 }
1533 
1535  const TrkrCluster* cluster)
1536 {
1537 
1538  if (!cluster) return nullptr;
1539 
1540  TVector3 pos(cluster->getPosition(0), cluster->getPosition(1), cluster->getPosition(2));
1541  TVector3 n(cluster->getPosition(0), cluster->getPosition(1), 0);
1542 
1543  // get the trkrid
1544  TrkrDefs::cluskey cluster_id = cluster->getClusKey();
1545  unsigned int trkrid = TrkrDefs::getTrkrId(cluster_id);
1546  int layer = TrkrDefs::getLayer(cluster_id);
1547 
1548  if(trkrid == TrkrDefs::mvtxId)
1549  {
1550  int stave_index = MvtxDefs::getStaveId(cluster_id);
1551  int chip_index = MvtxDefs::getChipId(cluster_id);
1552 
1553  double ladder_location[3] = {0.0, 0.0, 0.0};
1554  auto geom = dynamic_cast<CylinderGeom_Mvtx*>(_geom_container_maps->GetLayerGeom(layer));
1555  // returns the center of the sensor in world coordinates - used to get the ladder phi location
1556  geom->find_sensor_center(stave_index, 0,
1557  0, chip_index, ladder_location);
1558  //n.Print();
1559  n.SetXYZ(ladder_location[0], ladder_location[1], 0);
1560  n.RotateZ(geom->get_stave_phi_tilt());
1561  }
1562  else if(trkrid == TrkrDefs::inttId)
1563  {
1564  auto geom = dynamic_cast<CylinderGeomIntt*>(_geom_container_intt->GetLayerGeom(layer));
1565  double hit_location[3] = {0.0, 0.0, 0.0};
1566  geom->find_segment_center(InttDefs::getLadderZId(cluster_id),
1567  InttDefs::getLadderPhiId(cluster_id), hit_location);
1568 
1569  n.SetXYZ(hit_location[0], hit_location[1], 0);
1570  n.RotateZ(geom->get_strip_phi_tilt());
1571  }
1572 
1573  //double radius = sqrt(cluster->getPosition(0)*cluster->getPosition(0) + cluster->getPosition(1)*cluster->getPosition(1));
1574 
1576  cluster->getRPhiError(), cluster->getZError());
1577 
1578  meas->set_cluster_key(cluster->getClusKey());
1579 
1580 #ifdef _DEBUG_
1581  int layer_out = TrkrDefs::getLayer(cluster->getClusKey());
1582  cout
1583  << __LINE__
1584  << ": ID: " << cluster->getClusKey()
1585  << ": layer: " << layer_out
1586  << ": pos: {" << pos.X() << ", " << pos.Y() << ", " << pos.Z() << "}"
1587  << ": n: {" << n.X() << ", " << n.Y() << ", " << n.Z() << "}"
1588  << ": r*phi_error: " << cluster->getRPhiError()
1589  << ": z error: " << cluster->getZError()
1590  << endl;
1591 #endif
1592 
1593  return meas;
1594 }
1595 
1597 {
1598  // make this map for each collision vertex
1599  std::multimap<unsigned int, TrkrDefs::cluskey> this_layer_thetaID_phiID_clusterID;
1600 
1601  auto hitsetrange = _hitsets->getHitSets();
1602  for (auto hitsetitr = hitsetrange.first;
1603  hitsetitr != hitsetrange.second;
1604  ++hitsetitr){
1605  auto range = _cluster_map->getClusters(hitsetitr->first);
1606  for( auto clusIter = range.first; clusIter != range.second; ++clusIter ){
1607  TrkrCluster *cluster = clusIter->second;
1608  TrkrDefs::cluskey cluskey = clusIter->first;
1609 
1610  // This z-phi map relative to the collision vertex includes all clusters,
1611  // we do not know whch clusters belong to which vertices yet
1612  // make a map for every vertex?
1613 
1614  unsigned int layer = TrkrDefs::getLayer(cluskey);
1615  float x = cluster->getPosition(0) - _vertex[ivert][0];
1616  float y = cluster->getPosition(1) - _vertex[ivert][1];
1617  float z = cluster->getPosition(2) - _vertex[ivert][2];
1618 
1619  float phi = atan2(y, x);
1620  float r = sqrt(x * x + y * y);
1621  float theta = atan2(r, z);
1622 
1623 #ifdef _DEBUG_
1624  // //float rphi = r*phi;
1625  // std::cout
1626  // << __LINE__
1627  // <<": ID: " << cluster->get_id()
1628  // <<": layer: "<<cluster->get_layer()
1629  // <<", r: "<<r
1630  // <<", z: "<<z
1631  // <<", phi: "<<phi
1632  // <<", theta: "<<theta
1633  // <<endl;
1634 #endif
1635 
1636  unsigned int idx = encode_cluster_index(layer, theta, phi);
1637 
1638  // #ifdef _DEBUG_
1639  // cout
1640  // <<__LINE__<<": "
1641  // <<"{ "
1642  // <<layer <<", "
1643  // <<z <<", "
1644  // <<phi << "} =>"
1645  // <<idx << ": size: "
1646  // <<_layer_thetaID_phiID_clusterID.count(idx)
1647  // <<endl;
1648  // #endif
1649 
1650  this_layer_thetaID_phiID_clusterID.insert(std::make_pair(idx, cluster->getClusKey()));
1651  //cout << "BuildLayerPhiZHitmap for ivert "
1652  // << ivert << " : Inserted pair with x" << x << " y " << y << " z " << z << " idx " << idx << " cluskey " << cluster->getClusKey() << endl;
1653  }
1654  }
1655  _layer_thetaID_phiID_clusterID.push_back(this_layer_thetaID_phiID_clusterID);
1656 
1658 }
1659 
1660 std::vector<TrkrDefs::cluskey> PHGenFitTrkProp::SearchHitsNearBy(const unsigned int ivert, const unsigned int layer,
1661  const float theta_center, const float phi_center, const float theta_window,
1662  const float phi_window)
1663 {
1664  if(Verbosity() > 10) cout << "SearchHitsNearBy for ivert = " << ivert << endl;
1665 
1666  std::vector<TrkrDefs::cluskey> cluster_keys;
1667 
1668  const unsigned int max_phi_bin = 16383; //2^14 - 1
1669  const unsigned int max_z_bin = 2047; // 2^11 - 1
1670 
1671  float lower_phi = phi_center - phi_window + _half_max_phi;
1672  float upper_phi = phi_center + phi_window + _half_max_phi;
1673  float lower_z = theta_center - theta_window + _half_max_theta;
1674  float upper_z = theta_center + theta_window + _half_max_theta;
1675 
1676  unsigned int lower_phi_bin = (unsigned int) ((lower_phi) / _layer_thetaID_phiID_clusterID_phiSize);
1677  unsigned int upper_phi_bin = (unsigned int) ((upper_phi) / _layer_thetaID_phiID_clusterID_phiSize);
1678  unsigned int lower_z_bin = (unsigned int) ((lower_z) / _layer_thetaID_phiID_clusterID_zSize);
1679  unsigned int upper_z_bin = (unsigned int) ((upper_z) / _layer_thetaID_phiID_clusterID_zSize);
1680 
1681  if (lower_phi < 0) lower_phi_bin = 0;
1682  if (upper_phi_bin > max_phi_bin) upper_phi_bin = max_phi_bin;
1683 
1684  if (lower_z < 0) lower_z_bin = 0;
1685  if (upper_z_bin > max_z_bin) upper_z_bin = max_z_bin;
1686 
1687  for (unsigned int iz = lower_z_bin; iz <= upper_z_bin; ++iz)
1688  {
1689  for (unsigned int irphi = lower_phi_bin; irphi <= upper_phi_bin;
1690  ++irphi)
1691  {
1693  unsigned int idx = encode_cluster_index(layer, iz, irphi);
1695 
1696  if(Verbosity() > 10)
1697  {
1698  if(_layer_thetaID_phiID_clusterID[ivert].count(idx) > 0)
1699  {
1700  cout
1701  <<__LINE__<<": "
1702  <<"{ "
1703  <<layer <<", "
1704  <<iz <<", "
1705  <<irphi << "} =>"
1706  <<idx << ": size: "
1707  <<_layer_thetaID_phiID_clusterID[ivert].count(idx)
1708  <<endl;
1709  }
1710  }
1711 
1713  for (auto iter = _layer_thetaID_phiID_clusterID[ivert].lower_bound(idx);
1714  iter != _layer_thetaID_phiID_clusterID[ivert].upper_bound(idx);
1715  ++iter)
1716  {
1717  if(Verbosity() > 10) cout << " adding cluster with key " << iter->second << endl;
1718  cluster_keys.push_back(iter->second);
1719  }
1721  }
1722  }
1723 
1724  if(Verbosity() > 10)
1725  {
1726  // std::cout
1727  // <<__LINE__<<": "
1728  // <<"layer: "<<layer
1729  // <<": "<<phi_center <<" +- "<<phi_window
1730  // <<": "<<theta_center <<" +- "<<theta_window
1731  // <<endl;
1732  std::cout
1733  << __LINE__ << ": "
1734  << "layer: " << layer
1735  << ", rphi: {" << lower_phi_bin
1736  << ", " << upper_phi_bin
1737  << "}, z: {" << lower_z_bin
1738  << ", " << upper_z_bin
1739  << "}, found #clusters: " << cluster_keys.size()
1740  << endl;
1741  }
1742 
1743  return cluster_keys;
1744 }
1745 
1746 unsigned int PHGenFitTrkProp::encode_cluster_index(const unsigned int layer,
1747  const float theta, const float phi)
1748 {
1749  unsigned int idx = UINT_MAX;
1750 
1751  if (layer >= 128)
1752  {
1753  LogError("layer >= 128\n");
1754  return idx;
1755  }
1756 
1757  if ((theta + _half_max_theta) < 0)
1758  {
1759  LogError("(theta + _half_max_theta) < 0 \n");
1760  return idx;
1761  }
1762  unsigned int itheta = (theta + _half_max_theta) / _layer_thetaID_phiID_clusterID_zSize;
1763 
1764  if ((phi + _half_max_phi) < 0)
1765  {
1766  LogError("(rphi + _half_max_rphi) < 0 \n");
1767  return idx;
1768  }
1769  unsigned int irphi = (phi + _half_max_phi) / _layer_thetaID_phiID_clusterID_phiSize;
1770 
1771 #ifdef _DEBUG_
1772 // std::cout<<__LINE__<<": "
1773 // <<": layer: "<<layer
1774 // <<", irphi: "<<irphi
1775 // <<", itheta: "<<itheta
1776 // <<endl;
1777 #endif
1778 
1779  return encode_cluster_index(layer, itheta, irphi);
1780 }
1781 
1782 unsigned int PHGenFitTrkProp::encode_cluster_index(const unsigned int layer,
1783  const unsigned int iz, const unsigned int irphi)
1784 {
1785  if (layer >= 128)
1786  {
1787  LogError("layer >= 128: ") << layer << endl;
1788  ;
1789  return UINT_MAX;
1790  }
1791 
1792  if (iz >= 2048)
1793  {
1794  LogError("iz >= 2048: ") << iz << endl;
1795  return UINT_MAX;
1796  }
1797 
1798  if (irphi >= 16384)
1799  {
1800  LogError("irphi >= 16384: ") << irphi << endl;
1801  return UINT_MAX;
1802  }
1803 
1804  unsigned int index = 0;
1805 
1806  index |= (layer << 25);
1807  index |= (iz << 14);
1808  index |= (irphi);
1809 
1810  return index;
1811 }
1812 
1814 {
1815 
1816  auto micromegasgeos = findNode::getClass<PHG4CylinderGeomContainer>(topNode, "CYLINDERGEOM_MICROMEGAS");
1817  auto cellgeos = findNode::getClass<PHG4CylinderCellGeomContainer>(topNode, "CYLINDERCELLGEOM_SVTX");
1818  auto laddergeos = findNode::getClass<PHG4CylinderGeomContainer>(topNode, "CYLINDERGEOM_INTT");
1819  auto mapsladdergeos = findNode::getClass<PHG4CylinderGeomContainer>(topNode, "CYLINDERGEOM_MVTX");
1820 
1821  map<float, int> radius_layer_map;
1822 
1823  _radii_all.assign(_nlayers_all, 0.0);
1824  _layer_ilayer_map_all.clear();
1825 
1826  // micromegas
1827  if (micromegasgeos)
1828  {
1829  const auto range = micromegasgeos->get_begin_end();
1830  for (auto layeriter = range.first; layeriter != range.second; ++layeriter)
1831  { radius_layer_map.insert( std::make_pair(layeriter->second->get_radius(), layeriter->second->get_layer())); }
1832  }
1833 
1834  // tpc
1835  if (cellgeos)
1836  {
1838  cellgeos->get_begin_end();
1840  layerrange.first;
1841  layeriter != layerrange.second; ++layeriter)
1842  {
1843  radius_layer_map.insert(
1844  make_pair(layeriter->second->get_radius(),
1845  layeriter->second->get_layer()));
1846  }
1847  }
1848 
1849  // intt
1850  if (laddergeos)
1851  {
1853  laddergeos->get_begin_end();
1855  layerrange.first;
1856  layeriter != layerrange.second; ++layeriter)
1857  {
1858  radius_layer_map.insert(
1859  make_pair(layeriter->second->get_radius(),
1860  layeriter->second->get_layer()));
1861  }
1862  }
1863 
1864  // mvtx
1865  if (mapsladdergeos)
1866  {
1868  mapsladdergeos->get_begin_end();
1870  layerrange.first;
1871  layeriter != layerrange.second; ++layeriter)
1872  {
1873  radius_layer_map.insert(
1874  make_pair(layeriter->second->get_radius(),
1875  layeriter->second->get_layer()));
1876  }
1877  }
1878 
1879  for (map<float, int>::iterator iter = radius_layer_map.begin();
1880  iter != radius_layer_map.end(); ++iter)
1881  {
1882  _layer_ilayer_map_all.insert(make_pair(iter->second, _layer_ilayer_map_all.size()));
1883  }
1884 
1885  // now we extract the information from the cellgeos first
1886  // micromegas
1887  if (micromegasgeos)
1888  {
1889  const auto range = micromegasgeos->get_begin_end();
1890  for( auto iter = range.first; iter != range.second; iter++)
1891  {
1892  auto geo = iter->second;
1893  _radii_all[_layer_ilayer_map_all[geo->get_layer()]] = geo->get_radius() + 0.5 * geo->get_thickness();
1894  }
1895  }
1896 
1897  // tpc
1898  if (cellgeos)
1899  {
1901  cellgeos->get_begin_end();
1902  PHG4CylinderCellGeomContainer::ConstIterator miter = begin_end.first;
1903  for (; miter != begin_end.second; miter++)
1904  {
1905  PHG4CylinderCellGeom* geo = miter->second;
1906 
1908  geo->get_radius() + 0.5 * geo->get_thickness();
1909  }
1910  }
1911 
1912  // intt
1913  if (laddergeos)
1914  {
1916  laddergeos->get_begin_end();
1917  PHG4CylinderGeomContainer::ConstIterator miter = begin_end.first;
1918  for (; miter != begin_end.second; miter++)
1919  {
1920  PHG4CylinderGeom* geo = miter->second;
1921 
1923  geo->get_radius() + 0.5 * geo->get_thickness();
1924  }
1925  }
1926 
1927  // mvtx
1928  if (mapsladdergeos)
1929  {
1931  mapsladdergeos->get_begin_end();
1932  PHG4CylinderGeomContainer::ConstIterator miter = begin_end.first;
1933  for (; miter != begin_end.second; miter++)
1934  {
1935  PHG4CylinderGeom* geo = miter->second;
1936 
1938  geo->get_radius();
1939  }
1940  }
1942 }