EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
PHPatternReco.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file PHPatternReco.cc
1 
2 #include "PHPatternReco.h"
3 #include "CellularAutomaton_v1.h"
4 #include "CellularAutomaton.h" // for CellularAutom...
5 #include "HelixHoughBin.h" // for HelixHoughBin
6 #include "HelixHoughBin_v1.h" // for HelixHoughBin_v1
7 #include "HelixHoughFuncs.h" // for HelixHoughFuncs
8 #include "HelixHoughFuncs_v1.h" // for HelixHoughFun...
9 #include "HelixHoughSpace.h" // for HelixHoughSpace
10 #include "HelixHoughSpace_v1.h" // for HelixHoughSpa...
11 
12 
13 // g4hough includes
14 #include <trackbase_historic/SvtxVertexMap.h>
15 #include <trackbase_historic/SvtxVertexMap_v1.h>
16 #include <trackbase_historic/SvtxVertex.h>
17 #include <trackbase_historic/SvtxVertex_v1.h>
18 #include <trackbase_historic/SvtxTrackMap.h>
19 #include <trackbase_historic/SvtxTrackMap_v1.h>
20 #include <trackbase_historic/SvtxTrack.h>
21 #include <trackbase_historic/SvtxTrack_v1.h>
22 
23 #include <trackbase/TrkrClusterContainer.h>
24 #include <trackbase/TrkrCluster.h>
25 #include <trackbase/TrkrHitSetContainer.h>
26 
27 // sPHENIX Geant4 includes
28 #include <g4detectors/PHG4CylinderGeom.h>
29 #include <g4detectors/PHG4CylinderCellGeomContainer.h>
30 #include <g4detectors/PHG4CylinderCellGeom.h>
31 #include <g4detectors/PHG4CylinderGeomContainer.h>
32 
33 #include <g4bbc/BbcVertexMap.h>
34 
35 // sPHENIX includes
36 #include <fun4all/Fun4AllReturnCodes.h>
37 #include <fun4all/SubsysReco.h> // for SubsysReco
38 
39 #include <phool/PHCompositeNode.h>
40 #include <phool/PHIODataNode.h>
41 #include "phool/PHNode.h" // for PHNode
42 #include <phool/PHNodeIterator.h>
43 #include <phool/PHTimer.h> // for PHTimer
44 #include <phool/getClass.h>
45 #include <phool/phool.h> // for PHWHERE
46 
47 // standard includes
48 #include <TFile.h> // for TFile
49 #include <TNtuple.h>
50 
51 #include <algorithm> // for find
52 #include <climits> // for UINT_MAX
53 #include <cmath>
54 #include <cstdlib> // for NULL, exit
55 #include <iostream>
56 #include <map>
57 #include <memory> // for allocator_tra...
58 #include <utility> // for pair, make_pair
59 #include <limits>
60 
61 class PHObject;
62 
63 #define LogDebug(exp) std::cout<<"DEBUG: " <<__FILE__<<": "<<__LINE__<<": "<< exp
64 #define LogError(exp) std::cout<<"ERROR: " <<__FILE__<<": "<<__LINE__<<": "<< exp
65 #define LogWarning(exp) std::cout<<"WARNING: "<<__FILE__<<": "<<__LINE__<<": "<< exp
66 
67 //#define _DEBUG_
68 //#define _HOUGHTRANSF_
69 #define _MULTIVTX_
70 #define _TRIPLETS_
71 
72 using namespace std;
73 using namespace Eigen;
74 
75 
76 PHPatternReco::PHPatternReco(unsigned int nlayers,
77  unsigned int min_nlayers,
78  const string& name)
79  : SubsysReco(name),
80  _t_output_io(nullptr),
81  _seeding_layer(),
82  _nlayers(nlayers),
83  _min_nlayers(min_nlayers),
84  _ca_nlayers(nlayers),
85  _ca_min_nlayers(min_nlayers),
86  _ca_chi2(2.),
87  _radii(),
88  _material(),
89  _user_material(),
90  _mag_field(1.4),
91  _remove_hits(true),
92  _use_max_kappa(false),
93  fill_multi_zvtx(true),
94  _min_pt(0.2),
95  _max_kappa(numeric_limits<float>::max()),
96  _min_d(-1.0),
97  _max_d(1.0),
98  _min_dzdl(-0.9),
99  _max_dzdl(0.9),
100  _min_z0(-14.0),
101  _max_z0(+14.0),
102  _cut_on_dca(false),
103  _dcaxy_cut(0.2),
104  _dcaz_cut(0.05),
105  _pt_rescale(1.0),
106  _ca_phi_cut(M_PI/36.),
107  _ca_z_cut(20.),
108  _z_cut(1.),
109  _mult_onebin(1.),
110  _mult_twobins(1.),
111  _mult_threebins(1.),
112  _min_zvtx_tracks(1),
113  bin(0),
114  ik(0),
115  ip(0),
116  id(0),
117  il(0),
118  iz(0),
119  nkappa(0),
120  nphi(0),
121  nd(0),
122  ndzdl(0),
123  nz0(0),
124  cluster_id(0),
125  _layer_ilayer_map(),
126  _temp_tracks(),
127  _tracks(),
128  _track_states(),
129  _track_errors(),
130  _track_covars(),
131  _vertex(),
132  _vertex_list(),
133  _bbc_vertexes(nullptr),
134  _clustermap(nullptr),
135  _hitsets(nullptr),
136  _trackmap(nullptr),
137  _vertexmap(nullptr),
138  _vertex_finder(),
139  _hough_space(nullptr),
140  _hough_funcs(nullptr),
141  _mode(0),
142  _ntp_zvtx_by_event(nullptr),
143  _ntp_zvtx_by_track(nullptr),
144  _z0_dzdl(nullptr),
145  _kappa_phi(nullptr),
146  _d_phi(nullptr),
147  _kappa_d_phi(nullptr),
148  _ofile(nullptr),
149  _ofile2(nullptr),
150  _fname("test.root"),
151  _nlayers_all(67),
152  _layer_ilayer_map_all(),
153  _radii_all(),
154  zooms_vec(),
155  hits_used(),
156  hits_map(),
157  bins_map(),
158  bins_map_prev(),
159  bins_map_cur(),
160  bins_map_sel(),
161  kappa_phi_map(),
162  d_phi_map(),
163  kappa_d_phi_map(),
164  hit(),
165  hitpos3d(),
166  nzooms(0),
167  separate_helicity(false),
168  helicity(0),
169  n_vtx_tracks(0)
170  {
171  _event = 0;
172  zooms_vec.clear();
173 }
174 
175 PHPatternReco::~PHPatternReco()
176 {
177  delete _t_output_io;
178  if (_hough_space != nullptr) delete _hough_space;
179  if (_hough_funcs != nullptr) delete _hough_funcs;
180 // if (ca != nullptr) delete ca;
181 }
182 
183 int PHPatternReco::Init(PHCompositeNode* topNode) {
184 
185 #ifdef _DEBUG_
186 // _ofile = new TFile("z0_dzdl_kappa_phi_d.root","recreate");
187 #endif
188 
189  if (!_use_max_kappa){
190  _max_kappa = pt_to_kappa(_min_pt);
191 #ifdef _DEBUG_
192  cout<<"kappa max "<<_max_kappa<<endl;
193 #endif
194  }
195 
196  set_nzooms();
197  _hough_space = new HelixHoughSpace_v1();
198  for (unsigned int izoom =0; izoom<nzooms ; ++izoom)
199  _hough_space->add_one_zoom(zooms_vec[izoom]);
200 
201  _hough_space->set_kappa_min(0);
202  _hough_space->set_kappa_max(_max_kappa);
203  _hough_space->set_phi_min(0);
204  _hough_space->set_phi_max(2.*M_PI);
205  _hough_space->set_d_min(_min_d);
206  _hough_space->set_d_max(_max_d);
207  _hough_space->set_dzdl_min(_min_dzdl);// 0 if separate helicity/chargge
208  _hough_space->set_dzdl_max(_max_dzdl);// 0.9 default
209  _hough_space->set_z0_min(_min_z0);// -14 for non-zero vertex
210  _hough_space->set_z0_max(_max_z0);
211 
212  _hough_funcs = new HelixHoughFuncs_v1();
213  _hough_funcs->set_hough_space(_hough_space);
214 
215 #ifdef _DEBUG_
216 #ifdef _HOUGHTRANSF_
217 // _hough_space->print_zoom_profile();
218 // _hough_space->print_para_range();
219  unsigned int n_z0_bins= _hough_space->get_n_z0_bins(0);
220  unsigned int n_dzdl_bins= _hough_space->get_n_dzdl_bins(0);
221  unsigned int n_kappa_bins = _hough_space->get_n_kappa_bins(0);
222  unsigned int n_d_bins = _hough_space->get_n_d_bins(0);
223  unsigned int n_phi_bins = _hough_space->get_n_phi_bins(0);
224 
225  _z0_dzdl = new TH2D("z0_dzdl","z0_dzdl",n_z0_bins,_hough_space->get_z0_min(),_hough_space->get_z0_max(),n_dzdl_bins,_hough_space->get_dzdl_min(),_hough_space->get_dzdl_max());
226 // _z0_dzdl = new TH2D("z0_dzdl","z0_dzdl",n_z0_bins*_hough_space->get_n_z0_bins(1),_hough_space->get_z0_min(),_hough_space->get_z0_max(),n_dzdl_bins*_hough_space->get_n_dzdl_bins(1),_hough_space->get_dzdl_min(),_hough_space->get_dzdl_max());
227  _kappa_phi = new TH2D("kappa_phi","kappa_phi",n_kappa_bins,_hough_space->get_kappa_min(),_hough_space->get_kappa_max(),n_phi_bins,_hough_space->get_phi_min(),_hough_space->get_phi_max());
228 
229 // _kappa_phi = new TH2D("kappa_phi","kappa_phi",n_kappa_bins*_hough_space->get_n_kappa_bins(1),_hough_space->get_kappa_min(),_hough_space->get_kappa_max(),n_phi_bins*_hough_space->get_n_phi_bins(1),_hough_space->get_phi_min(),_hough_space->get_phi_max());
230  _d_phi = new TH2D("d_phi","d_phi",n_d_bins,_hough_space->get_d_min(),_hough_space->get_d_max(),n_phi_bins,_hough_space->get_phi_min(),_hough_space->get_phi_max());
231 // _d_phi = new TH2D("d_phi","d_phi",n_d_bins*_hough_space->get_n_d_bins(1),_hough_space->get_d_min(),_hough_space->get_d_max(),n_phi_bins*_hough_space->get_n_phi_bins(1),_hough_space->get_phi_min(),_hough_space->get_phi_max());
232 
233  _kappa_d_phi = new TH3D("kappa_d_phi","kappa_d_phi",n_kappa_bins,_hough_space->get_kappa_min(),_hough_space->get_kappa_max(),n_d_bins,_hough_space->get_d_min(),_hough_space->get_d_max(),n_phi_bins,_hough_space->get_phi_min(),_hough_space->get_phi_max());
234 #endif
235 #endif
236  return Fun4AllReturnCodes::EVENT_OK;
237 }
238 
239 int PHPatternReco::InitRun(PHCompositeNode* topNode) {
240 
241  int code = Fun4AllReturnCodes::ABORTRUN;
242 
243  code = create_nodes(topNode);
244  if(code != Fun4AllReturnCodes::EVENT_OK)
245  return code;
246 
247  code = initialize_geometry(topNode);
248  if(code != Fun4AllReturnCodes::EVENT_OK)
249  return code;
250 
251 #ifdef _MULTIVTX_
252  _ofile2 = new TFile(_fname.c_str(),"recreate");
253  _ntp_zvtx_by_event = new TNtuple("ntp_zvtx_by_event","all vertices found event-by-event","event:zvtx");
254  _ntp_zvtx_by_track = new TNtuple("ntp_zvtx_by_track","track-by-track (zvtx + z0) distribution", "event:zvtx");
255 #endif
256 
257  _t_output_io = new PHTimer("_t_output_io");
258  _t_output_io->stop();
259 
260  if (Verbosity() > 0) {
261  cout
262  << "====================== PHPatternReco::InitRun() ======================"
263  << endl;
264  cout << " Magnetic field set to: " << _mag_field << " Tesla" << endl;
265  cout << " Number of tracking layers: " << _nlayers << endl;
266  for (unsigned int i = 0; i < _nlayers; ++i) {
267  cout << " Tracking layer #" << i << " " << "radius = "
268  << _radii[i] << " cm, " << "material = " << _material[i]
269  << endl;
270  }
271  cout << " Required hits: " << _min_nlayers << endl;
272  cout << " Minimum pT: " << _min_pt << endl;
273  cout << " Maximum DCA: " << boolalpha << _cut_on_dca << noboolalpha
274  << endl;
275  if (_cut_on_dca) {
276  cout << " Maximum DCA cut: " << _dcaxy_cut << endl;
277  }
278  cout << " Maximum DCAZ cut: " << _dcaz_cut << endl;
279  cout << " Momentum rescale factor: " << _pt_rescale << endl;
280  cout
281  << "==========================================================================="
282  << endl;
283  }
284 
285 
286  return code;
287 }
288 
289 int PHPatternReco::process_event(PHCompositeNode *topNode)
290 {
291 
292  if (Verbosity() > 0) cout << "PHPatternReco::process_event -- entered" << endl;
293 
294  // start fresh
295  for(unsigned int i=0; i<_tracks.size(); ++i) _tracks[i].reset();
296  _tracks.clear();
297  _track_states.clear();
298  _track_errors.clear();
299  _track_covars.clear();
300  _vertex.clear();
301  _vertex.assign(3, 0.0);
302  _vertex_list.clear();
303 
304  hits_map.clear();
305  hits_used.clear();
306 
307  //-----------------------------------
308  // Get Objects off of the Node Tree
309  //-----------------------------------
310 
311  get_nodes(topNode);
312 
313  int code = translate_input(topNode);
314  if (code != Fun4AllReturnCodes::EVENT_OK)
315  return code;
316 
317  // First iteration
318 
319 #ifdef _HOUGHTRANSF_
320  reset_zooms();
321  add_zoom(1,18,1,16,16);// zoom 0
322  add_zoom(1,3,1,3,2);// zoom 1
323  add_zoom(1,3,1,3,2);// zoom 2
324  separate_helicity = true;
325 #endif
326 
327  Init(topNode);
328 
329  turnoff_hits_used_in_triplets();
330 
331 // bool zvtx_found = false;
332  _ca_nlayers = 4; // 3(vertexing)/4(track seeding) for p+p
333 
334  // p+p
335 // unsigned int nseq =4; // 65(0.5) 103(0.5)
336 // unsigned int nattempt =4;
337  // N+N
338  unsigned int nseq = 1;
339  unsigned int nattempt = 1; // 59(1.5) 63(1.5)
340 
341  unsigned int iseq=0;
342  while(1)
343  {// iseq
344 
345  if (iseq==nseq) break;
346  for(unsigned int i=0; i<_tracks.size(); ++i) _tracks[i].reset();
347  _tracks.clear();
348  _track_states.clear();
349  _track_errors.clear();
350  _track_covars.clear();
351 // reset_hits_used();
352 
353  float shift_dx = -_vertex[0];
354  float shift_dy = -_vertex[1];
355  float shift_dz = -_vertex[2];
356  shift_coordinate_system(shift_dx,shift_dy,shift_dz);
357 
358  for (unsigned int iattempt =0; iattempt<nattempt; ++iattempt ){
359 #ifdef _DEBUG_
360  cout<<iattempt << " th attempt "<<endl;
361 #endif
362  helicity = 1;
363 
364 #ifdef _TRIPLETS_
365  _temp_tracks.clear();
366 
367  build_triplets_to_SimpleTrack3D(_temp_tracks, true); // false : backward
368 #endif
369 
370 #ifdef _HOUGHTRANSF_
371 
372  unsigned int zoomlevel = 0;
373  zoomlevel = 0;
374  set_nbins(zoomlevel);
375 
376 #ifdef _DEBUG_
377 // cout<<"PHPatternReco:: nkappa " <<nkappa<<" nphi "<<nphi<<" nd "<<nd<<" ndzdl "<<ndzdl<<" nz0 " <<nz0<<endl;
378 #endif
379 
380  for(unsigned int i=0; i<_temp_tracks.size(); ++i) _temp_tracks[i].reset();
381  _temp_tracks.clear();
382  auto it_begin = bins_map.begin();
383  auto it_end = bins_map.end();
384  while(it_begin != it_end)
385  {
386  delete it_begin->second;
387  bins_map.erase(it_begin++);
388  }
389  it_begin = bins_map_prev.begin();
390  it_end = bins_map_prev.end();
391  while(it_begin != it_end)
392  {
393  delete it_begin->second;
394  bins_map_prev.erase(it_begin++);
395  }
396  it_begin = bins_map_cur.begin();
397  it_end = bins_map_cur.end();
398  while(it_begin != it_end)
399  {
400  delete it_begin->second;
401  bins_map_cur.erase(it_begin++);
402  }
403  it_begin = bins_map_sel.begin();
404  it_end = bins_map_sel.end();
405  while(it_begin != it_end)
406  {
407  delete it_begin->second;
408  bins_map_sel.erase(it_begin++);
409  }
410 
411  bins_map.clear();
412  bins_map_prev.clear();
413  bins_map_cur.clear();
414  bins_map_sel.clear();
415 
416  initialize_houghbin();
417  vote_z_init(zoomlevel);
418  find_track_candidates_z_init(zoomlevel);
419 
420  vote_xy(zoomlevel);
421  find_track_candidates_xy(zoomlevel);
422  prune_xy(zoomlevel);// on prev
423 
424  bins_map_prev.swap(bins_map_sel);
425  prune_z(zoomlevel); // on sel
426  bins_map_sel.swap(bins_map_prev);
427 
428  for (zoomlevel =1; zoomlevel<2; ++zoomlevel){ // depends on the total number of tracks thrown in
429  vote_z(zoomlevel); // on prev
430  find_track_candidates_z(zoomlevel);
431  prune_z(zoomlevel);
432  vote_xy(zoomlevel);
433  find_track_candidates_xy(zoomlevel);
434  prune_xy(zoomlevel);
435  } //zoomlevel
436 
437 
438  // high pt, primary vertex originating tracks
439  bins_to_SimpleTrack3D(_temp_tracks,1,zoomlevel);// 0(z) 1(xy) : create SimpleTrack3D objects, save cluster_id
440 #endif
441 
442  code = 0;
443  code = cellular_automaton_zvtx_init(_temp_tracks); // remove bad hits, save fit parameters && turn off used hits
444  if (code != Fun4AllReturnCodes::EVENT_OK)
445  {
446  cout<<"CellularAutomaton failed. "<<endl;
447  exit(1);
448  }
449 
450 
451  }// iattempt
452 
453  shift_coordinate_system(-shift_dx,-shift_dy,-shift_dz);
454  ++iseq;
455 
456  }//iseq
457 
458 #ifdef _DEBUG_
459  cout<<"export output"<<endl;
460 #endif
461  code = export_output();
462  if (code != Fun4AllReturnCodes::EVENT_OK)
463  return code;
464 
465  ++_event;
466 
467 // BbcVertex* vertex = _bbc_vertexes->begin()->second;
468 // cout<< "true z-vertex " << vertex->get_z()<<endl;
469 
470  return Fun4AllReturnCodes::EVENT_OK;
471 }
472 
473 int PHPatternReco::End(PHCompositeNode *topNode) {
474 #ifdef _DEBUG_
475  _ofile->cd();
476  _z0_dzdl->Write();
477  _kappa_phi->Write();
478  _d_phi->Write();
479  _kappa_d_phi->Write();
480  _ofile->Close();
481 #endif
482 
483 #ifdef _MULTIVTX_
484  _ofile2->cd();
485  _ntp_zvtx_by_event->Write();
486  _ntp_zvtx_by_track->Write();
487  _ofile2->Close();
488 #endif
489 
490  _temp_tracks.clear();
491 
492  return Fun4AllReturnCodes::EVENT_OK;
493 }
494 
495 
496 void PHPatternReco::set_material(int layer, float value) {
497  _user_material[layer] = value;
498 }
499 
500 void PHPatternReco::add_zoom(unsigned int n_kappa, unsigned int n_phi, unsigned int n_d,unsigned int n_dzdl, unsigned int n_z0) {
501  std::vector<unsigned int> zoom {n_kappa, n_phi, n_d, n_dzdl, n_z0};
502  zooms_vec.push_back(zoom);
503 }
504 
505 int PHPatternReco::create_nodes(PHCompositeNode* topNode) {
506  // create nodes...
507  PHNodeIterator iter(topNode);
508 
509  PHCompositeNode* dstNode = static_cast<PHCompositeNode*>(iter.findFirst(
510  "PHCompositeNode", "DST"));
511  if (!dstNode) {
512  cerr << PHWHERE << "DST Node missing, doing nothing." << endl;
513  return Fun4AllReturnCodes::ABORTEVENT;
514  }
515  PHNodeIterator iter_dst(dstNode);
516 
517  // Create the SVTX node
518  PHCompositeNode* tb_node =
519  dynamic_cast<PHCompositeNode*>(iter_dst.findFirst("PHCompositeNode",
520  "SVTX"));
521  if (!tb_node) {
522  tb_node = new PHCompositeNode("SVTX");
523  dstNode->addNode(tb_node);
524  if (Verbosity() > 0)
525  cout << "SVTX node added" << endl;
526  }
527 
528  _trackmap = new SvtxTrackMap_v1;
529  PHIODataNode<PHObject>* tracks_node = new PHIODataNode<PHObject>(_trackmap,
530  "SvtxTrackMap", "PHObject");
531  tb_node->addNode(tracks_node);
532  if (Verbosity() > 0)
533  cout << "Svtx/SvtxTrackMap node added" << endl;
534 
535  _vertexmap = new SvtxVertexMap_v1;
536  PHIODataNode<PHObject>* vertexes_node = new PHIODataNode<PHObject>(
537  _vertexmap, "SvtxVertexMap", "PHObject");
538  tb_node->addNode(vertexes_node);
539  if (Verbosity() > 0)
540  cout << "Svtx/SvtxVertexMap node added" << endl;
541 
542 
543  return Fun4AllReturnCodes::EVENT_OK;
544 }
545 
546 int PHPatternReco::initialize_geometry(PHCompositeNode *topNode) {
547 
548  //---------------------------------------------------------
549  // Grab Run-Dependent Detector Geometry and Configure Hough
550  //---------------------------------------------------------
551 
552  PHG4CylinderCellGeomContainer* cellgeos = findNode::getClass<
553  PHG4CylinderCellGeomContainer>(topNode, "CYLINDERCELLGEOM_SVTX");
554  PHG4CylinderGeomContainer* laddergeos = findNode::getClass<
555  PHG4CylinderGeomContainer>(topNode, "CYLINDERGEOM_INTT");
556  PHG4CylinderGeomContainer* mvtxladdergeos = findNode::getClass<
557  PHG4CylinderGeomContainer>(topNode, "CYLINDERGEOM_MVTX");
558 
559 
560  _nlayers = _seeding_layer.size();
561 
562 
563  _radii.assign(_nlayers, 0.0);
564  map<float, int> radius_layer_map;
565 
566  _radii_all.assign(_nlayers_all, 0.0);
567 
568  if (cellgeos) {
569  PHG4CylinderCellGeomContainer::ConstRange layerrange =
570  cellgeos->get_begin_end();
571  for (PHG4CylinderCellGeomContainer::ConstIterator layeriter =
572  layerrange.first; layeriter != layerrange.second; ++layeriter) {
573  radius_layer_map.insert(
574  make_pair(layeriter->second->get_radius(),
575  layeriter->second->get_layer()));
576  }
577  }
578 
579  if (laddergeos) {
580  PHG4CylinderGeomContainer::ConstRange layerrange =
581  laddergeos->get_begin_end();
582  for (PHG4CylinderGeomContainer::ConstIterator layeriter =
583  layerrange.first; layeriter != layerrange.second; ++layeriter) {
584  radius_layer_map.insert(
585  make_pair(layeriter->second->get_radius(),
586  layeriter->second->get_layer()));
587  }
588  }
589 
590  if (mvtxladdergeos) {
591  PHG4CylinderGeomContainer::ConstRange layerrange =
592  mvtxladdergeos->get_begin_end();
593  for (PHG4CylinderGeomContainer::ConstIterator layeriter =
594  layerrange.first; layeriter != layerrange.second; ++layeriter) {
595  radius_layer_map.insert(
596  make_pair(layeriter->second->get_radius(),
597  layeriter->second->get_layer()));
598  }
599  }
600 
601 
602  // now that the layer ids are sorted by radius, I can create a storage
603  // index, ilayer, that is 0..N-1 and sorted by radius
604 
605  int ilayer = 0;
606  for (map<float, int>::iterator iter = radius_layer_map.begin();
607  iter != radius_layer_map.end(); ++iter) {
608 
609  _layer_ilayer_map_all.insert(make_pair(iter->second, _layer_ilayer_map_all.size()));
610 
611  if (std::find(_seeding_layer.begin(), _seeding_layer.end(),
612  iter->second) != _seeding_layer.end()) {
613  _layer_ilayer_map.insert(make_pair(iter->second, ilayer));
614  ++ilayer;
615  }
616  }
617 
618 
619  // now we extract the information from the cellgeos first
620  if (cellgeos) {
621  PHG4CylinderCellGeomContainer::ConstRange begin_end =
622  cellgeos->get_begin_end();
623  PHG4CylinderCellGeomContainer::ConstIterator miter = begin_end.first;
624  for (; miter != begin_end.second; ++miter) {
625  PHG4CylinderCellGeom *geo = miter->second;
626 
627  //if(cellgeo->get_layer() > (int) _radii.size() ) continue;
628 
629 // if (Verbosity() >= 2)
630 // cellgeo->identify();
631 
632  //TODO
633  _radii_all[_layer_ilayer_map_all[geo->get_layer()]] =
634  geo->get_radius() + 0.5*geo->get_thickness();
635 
636  if (_layer_ilayer_map.find(geo->get_layer())
637  != _layer_ilayer_map.end()) {
638  _radii[_layer_ilayer_map[geo->get_layer()]] =
639  geo->get_radius();
640  }
641  }
642  }
643 
644  if (laddergeos) {
645  PHG4CylinderGeomContainer::ConstRange begin_end =
646  laddergeos->get_begin_end();
647  PHG4CylinderGeomContainer::ConstIterator miter = begin_end.first;
648  for (; miter != begin_end.second; ++miter) {
649  PHG4CylinderGeom *geo = miter->second;
650 
651  //if(geo->get_layer() > (int) _radii.size() ) continue;
652 
653 // if (Verbosity() >= 2)
654 // geo->identify();
655 
656  _radii_all[_layer_ilayer_map_all[geo->get_layer()]] =
657  geo->get_radius() + 0.5*geo->get_thickness();
658 
659  if (_layer_ilayer_map.find(geo->get_layer())
660  != _layer_ilayer_map.end()) {
661  _radii[_layer_ilayer_map[geo->get_layer()]] = geo->get_radius();
662  }
663  }
664  }
665 
666  if (mvtxladdergeos) {
667  PHG4CylinderGeomContainer::ConstRange begin_end =
668  mvtxladdergeos->get_begin_end();
669  PHG4CylinderGeomContainer::ConstIterator miter = begin_end.first;
670  for (; miter != begin_end.second; ++miter) {
671  PHG4CylinderGeom *geo = miter->second;
672 
673  //if(geo->get_layer() > (int) _radii.size() ) continue;
674 
675 // if (Verbosity() >= 2)
676 // geo->identify();
677 
678  _radii_all[_layer_ilayer_map_all[geo->get_layer()]] =
679  geo->get_radius() + 0.5*geo->get_thickness();
680 
681  if (_layer_ilayer_map.find(geo->get_layer())
682  != _layer_ilayer_map.end()) {
683  _radii[_layer_ilayer_map[geo->get_layer()]] = geo->get_radius();
684  }
685  }
686  }
687 
688  // set material on each layer
689 
690  _material.assign(_radii.size(), 0.03);
691 
692  map<int, float>::iterator mat_it;
693  for (map<int, float>::iterator iter = _user_material.begin();
694  iter != _user_material.end(); ++iter) {
695  if (_layer_ilayer_map.find(iter->first) != _layer_ilayer_map.end()) {
696  _material[_layer_ilayer_map[iter->first]] = iter->second;
697  }
698  }
699 
700  // initialize the pattern recogition tools
701 
702 
703  return Fun4AllReturnCodes::EVENT_OK;
704 }
705 
706 
707 int PHPatternReco::get_nodes(PHCompositeNode* topNode) {
708 
709  //---------------------------------
710  // Get Objects off of the Node Tree
711  //---------------------------------
712 
713  // Pull the reconstructed track information off the node tree...
714  _bbc_vertexes = findNode::getClass<BbcVertexMap>(topNode, "BbcVertexMap");
715 
716  /*
717  _clustermap = findNode::getClass<SvtxClusterMap>(topNode, "SvtxClusterMap");
718  if (!_clustermap) {
719  cerr << PHWHERE << " ERROR: Can't find node SvtxClusterMap" << endl;
720  return Fun4AllReturnCodes::ABORTEVENT;
721  }
722  */
723 
724  _clustermap = findNode::getClass<TrkrClusterContainer>(topNode, "TRKR_CLUSTER");
725  if (!_clustermap)
726  {
727  cerr << PHWHERE << " ERROR: Can't find node TrkrClusterContainer" << endl;
728  return Fun4AllReturnCodes::ABORTEVENT;
729  }
730  _hitsets = findNode::getClass<TrkrHitSetContainer>(topNode, "TRKR_HITSET");
731  if(!_hitsets)
732  {
733  std::cout << PHWHERE << "No hitset container on node tree. Bailing."
734  << std::endl;
735  return Fun4AllReturnCodes::ABORTEVENT;
736  }
737  // Pull the reconstructed track information off the node tree...
738  _trackmap = findNode::getClass<SvtxTrackMap>(topNode, "SvtxTrackMap");
739  if (!_trackmap) {
740  cerr << PHWHERE << " ERROR: Can't find SvtxTrackMap." << endl;
741  return Fun4AllReturnCodes::ABORTEVENT;
742  }
743 
744  // Pull the reconstructed track information off the node tree...
745  _vertexmap = findNode::getClass<SvtxVertexMap>(topNode, "SvtxVertexMap");
746  if (!_vertexmap) {
747  cerr << PHWHERE << " ERROR: Can't find SvtxVertexMap." << endl;
748  return Fun4AllReturnCodes::ABORTEVENT;
749  }
750 
751  return Fun4AllReturnCodes::EVENT_OK;
752 }
753 
754 int PHPatternReco::translate_input(PHCompositeNode* topNode) {
755 
756  unsigned int clusid = 0;
757  unsigned int ilayer = 0;
758 
759  auto hitsetrange = _hitsets->getHitSets();
760  for (auto hitsetitr = hitsetrange.first;
761  hitsetitr != hitsetrange.second;
762  ++hitsetitr){
763  auto range = _clustermap->getClusters(hitsetitr->first);
764  for( auto clusIter = range.first; clusIter != range.second; ++clusIter ){
765  TrkrCluster *cluster = clusIter->second;
766  TrkrDefs::cluskey cluskey = clusIter->first;
767  unsigned int layer = TrkrDefs::getLayer(cluskey);
768  //cout << " cluster with layer " << layer << " _nlayers = " << _nlayers << endl;
769  std::map<int, unsigned int>::const_iterator it = _layer_ilayer_map.find(layer);
770  if(it != _layer_ilayer_map.end())
771  {
772  ilayer = it->second;
773  //cout << " set ilayer to " << ilayer << endl;
774 
775  if(ilayer >= _nlayers) continue;
776 
777  SimpleHit3D hit3d;
778 
779  hit3d.set_cluskey(cluskey); // will be added to SvtxTrack later, needed to find the cluster in TrkrClusterContainer
780  hit3d.set_id(clusid);
781  hit3d.set_layer(ilayer);
782 
783  // cout<<" adding cluster, id = "<< cluster->get_id()<<", layer = "<<cluster->get_layer()<<endl;
784 
785  hit3d.set_x(cluster->getPosition(0));
786  hit3d.set_y(cluster->getPosition(1));
787  hit3d.set_z(cluster->getPosition(2));
788  for (int i = 0; i < 3; ++i) {
789  for (int j = i; j < 3; ++j) {
790  hit3d.set_error(i, j, cluster->getError(i, j));
791 
792  //hit3d.set_size(i, j, cluster->get_size(i, j)); // original
793  hit3d.set_size(i, j, cluster->getError(i, j)*sqrt(12.)); // yuhw 2017-05-08
794  }
795  }
796  cout << " inserting hit with layer " << ilayer << " into hits_map" << endl;
797  hits_map.insert(std::pair<unsigned int, SimpleHit3D>(hit3d.get_id(),hit3d));
798  hits_used.insert(std::pair<unsigned int, bool>(hit3d.get_id(),false));
799 
800  clusid += 1;
801  }
802  }
803  }
804 
805  if (Verbosity() > 10) {
806  cout << "-------------------------------------------------------------------"
807  << endl;
808  cout << "PHPatternReco::process_event has the following input clusters:"
809  << endl;
810 
811  cout << "n init clusters = " << hits_map.size() << endl;
812 
813  for (std::map<unsigned int,SimpleHit3D>::iterator it=hits_map.begin();
814  it!=hits_map.end();
815  ++it)
816  {
817  SimpleHit3D hit = it->second;
818  hit.print();
819  }
820 
821 
822  cout << "-------------------------------------------------------------------"
823  << endl;
824  }
825 
826  return Fun4AllReturnCodes::EVENT_OK;
827 }
828 
829 int PHPatternReco::export_output(){
830 
831  // clear mvtx at the begining of an event
832  // _vertexmap->clear();
833  // _trackmap->clear();
834  // if (_tracks.empty())
835  // return Fun4AllReturnCodes::EVENT_OK;
836  /*
837  std::vector<SvtxVertex_v1> svtx_vertex_list;
838  unsigned int nvertex = _vertex_list.size();
839  for (unsigned int vid = 0; vid < nvertex; ++vid ){
840 
841  SvtxVertex_v1 vertex;
842  vertex.set_t0(0.0);
843  for (int i = 0; i < 3; ++i)
844  vertex.set_position(i, _vertex_list[vid][i]);
845  vertex.set_chisq(0.0);
846  vertex.set_ndof(0);
847  vertex.set_error(0, 0, 0.0);
848  vertex.set_error(0, 1, 0.0);
849  vertex.set_error(0, 2, 0.0);
850  vertex.set_error(1, 0, 0.0);
851  vertex.set_error(1, 1, 0.0);
852  vertex.set_error(1, 2, 0.0);
853  vertex.set_error(2, 0, 0.0);
854  vertex.set_error(2, 1, 0.0);
855  vertex.set_error(2, 2, 0.0);
856 
857  svtx_vertex_list.push_back(vertex);
858  }
859 
860  cout<<"vertex list "<<endl;
861  // at this point we should already have an initial pt and pz guess...
862  // need to translate this into the PHG4Track object...
863  */
864  // read vertices from SvtxVertexMap
865  // store vertex object in a map (so that it can have an vertex_id = -1)
866  // and assign a vertex_id to each track
867  std::map<int, SvtxVertex*> svtx_vertex_list;
868 
869  // loop over all vertexes on node
870  for (SvtxVertexMap::Iter iter = _vertexmap->begin();
871  iter != _vertexmap->end();
872  ++iter) {
873  SvtxVertex* vertex = iter->second;
874  int vertex_id = (int) vertex->get_id();
875  cout<<"vertex id "<< vertex->get_id()<<", z " <<vertex->get_position(2)<<endl;
876  svtx_vertex_list[vertex_id] = vertex;
877  }
878  // add a vertex with a vertex_id = 9999 for triplets not tied to a reco vertex
879  SvtxVertex_v1 fake_vertex;
880  fake_vertex.set_t0(0.0);
881  for (int i = 0; i < 3; ++i)
882  fake_vertex.set_position(i,999.0);
883  fake_vertex.set_chisq(0.0);
884  fake_vertex.set_ndof(0);
885  fake_vertex.set_error(0, 0, 0.0);
886  fake_vertex.set_error(0, 1, 0.0);
887  fake_vertex.set_error(0, 2, 0.0);
888  fake_vertex.set_error(1, 0, 0.0);
889  fake_vertex.set_error(1, 1, 0.0);
890  fake_vertex.set_error(1, 2, 0.0);
891  fake_vertex.set_error(2, 0, 0.0);
892  fake_vertex.set_error(2, 1, 0.0);
893  fake_vertex.set_error(2, 2, 0.0);
894  // svtx_vertex_list[9999] = fake_vertex;
895 
896  std::vector<SimpleHit3D> track_hits;
897  track_hits.clear();
898  TrkrDefs::cluskey clusterkey;
899 
900  for (unsigned int itrack = 0; itrack < _tracks.size(); ++itrack) {
901  // if (_tracks.at(itrack).vertex_id >= nvertex) continue;
902  // a triplet is allowed to have no reco vertex tied to it, assign -1
903  SvtxTrack_v1 track;
904  track.set_id(itrack);
905  track_hits.clear();
906  track_hits = _tracks.at(itrack).hits;
907  // cout<<"hits vector assigned " <<endl;
908  for (unsigned int ihit = 0; ihit < track_hits.size(); ++ihit) {
909  if ((track_hits.at(ihit).get_id()) >= _clustermap->size()) {
910  continue;
911  }
912 
913  clusterkey = track_hits.at(ihit).get_cluskey();
914 
915  if (Verbosity() > 5) {
916  int layer = TrkrDefs::getLayer(clusterkey);
917  cout
918  <<__LINE__
919  <<": itrack: " << itrack
920  <<": nhits: " << track_hits.size()
921  <<": clusterkey: " << clusterkey
922  <<": layer: " << layer
923  <<endl;
924  }
925  track.insert_cluster_key(clusterkey);
926  }
927 
928  // cout<<"hits added to a track "<<endl;
929  float kappa = _tracks.at(itrack).kappa;
930  float d = _tracks.at(itrack).d;
931  float phi = _tracks.at(itrack).phi;
932  float dzdl = _tracks.at(itrack).dzdl;
933  float z0 = _tracks.at(itrack).z0;
934 
935  float pT = kappa_to_pt(kappa);
936 
937  float x_center = cos(phi) * (d + 1 / kappa); // x coordinate of circle center
938  float y_center = sin(phi) * (d + 1 / kappa); // y " " " "
939 
940  // find helicity from cross product sign
941  short int helicity;
942  if ((track_hits[0].get_x()- x_center)*(track_hits[track_hits.size()-1].get_y()- y_center)
943  - (track_hits[0].get_y()- y_center)*(track_hits[track_hits.size()-1].get_x()- x_center)> 0)
944  {
945  helicity = 1;
946  } else {
947  helicity = -1;
948  }
949 
950  float pZ = 0;
951  if (dzdl != 1) {
952  pZ = pT * dzdl / sqrt(1.0 - dzdl * dzdl);
953  }
954  int ndf = 2 * _tracks.at(itrack).hits.size() - 5;
955  // track.set_chisq(_track_errors[itrack]);
956  track.set_ndf(ndf);
957  track.set_px(pT * cos(phi - helicity * M_PI / 2));
958  track.set_py(pT * sin(phi - helicity * M_PI / 2));
959  track.set_pz(pZ);
960 
961  track.set_dca2d(d);
962  // track.set_dca2d_error(sqrt(_track_covars[itrack](1, 1)));
963 
964  // cout<<"track parameters set "<<endl;
965  if (_mag_field > 0) {
966  track.set_charge(helicity);
967  } else {
968  track.set_charge(-1.0 * helicity);
969  }
970 
971 
972  /* Eigen::Matrix<float, 6, 6> euclidean_cov =
973  Eigen::Matrix<float, 6, 6>::Zero(6, 6);
974  convertHelixCovarianceToEuclideanCovariance(_mag_field, phi, d, kappa,
975  z0, dzdl, _track_covars[itrack], euclidean_cov);
976  */
977  for (unsigned int row = 0; row < 6; ++row) {
978  for (unsigned int col = 0; col < 6; ++col) {
979  // track.set_error(row, col, euclidean_cov(row, col));
980  }
981  }
982  // cout<<"set track "<<endl;
983  unsigned int vid =9999;// _tracks.at(itrack).vertex_id
984  float distance = 9999;
985  for (std::map<int,SvtxVertex*>::iterator it=svtx_vertex_list.begin();
986  it!=svtx_vertex_list.end();
987  ++it){
988  unsigned int iv = it->second->get_id();
989  float zvtx = it->second->get_position(2);
990  if (fabs(z0-zvtx) <distance){
991  vid = iv;
992  distance = fabs(z0-zvtx);
993  }
994  }
995 
996  //if a triplet doesn't have a reco vertex tied to it, assign 9999
997  if (fabs(z0 - distance)>0.05) vid = 9999;
998  if (Verbosity() > 5) cout<<"vertex_id "<<vid<<endl;
999  // pca
1000  if (vid==9999){
1001  track.set_x(0.);
1002  track.set_y(0.);
1003  track.set_z(0.);
1004  }else {
1005  track.set_x(svtx_vertex_list[vid]->get_x() + d * cos(phi));
1006  track.set_y(svtx_vertex_list[vid]->get_y() + d * sin(phi));
1007  track.set_z(svtx_vertex_list[vid]->get_z() + z0);
1008  }
1009 
1010  _trackmap->insert(&track);
1011  if (vid==9999) fake_vertex.insert_track(track.get_id());
1012  else svtx_vertex_list[vid]->insert_track(track.get_id());
1013 
1014  if (Verbosity() > 5) {
1015  cout << "track " << itrack << " quality = " << track.get_quality() << endl;
1016  cout << "px = " << track.get_px() << " py = " << track.get_py() << " pz = " << track.get_pz() << endl;
1017  }
1018  } // track loop
1019 
1020  // do not repeat saving vertices this time, add only fake vertex
1021  SvtxVertex *vtxptr = _vertexmap->insert_clone(&fake_vertex);
1022  if (Verbosity() > 5) vtxptr->identify();
1023 
1024  /*
1025  for (unsigned int vid = 0; vid < _vertex_list.size(); ++vid ){
1026  SvtxVertex *vtxptr = _vertexmap->insert_clone(&svtx_vertex_list[vid]);
1027  if (Verbosity() > 5) vtxptr->identify();
1028  }
1029  */
1030  hits_map.clear();
1031 
1032  for(unsigned int i=0; i<_tracks.size(); ++i) _tracks[i].reset();
1033  _tracks.clear();
1034  _track_errors.clear();
1035  _track_covars.clear();
1036  _vertex.clear();
1037  _vertex.assign(3, 0.0);
1038  _vertex_list.clear();
1039 
1040  return Fun4AllReturnCodes::EVENT_OK;
1041 
1042 }
1043 
1044 float PHPatternReco::kappa_to_pt(float kappa) {
1045  return _pt_rescale * _mag_field / 333.6 / kappa;
1046 }
1047 
1048 float PHPatternReco::pt_to_kappa(float pt) {
1049  return _pt_rescale * _mag_field / 333.6 / pt;
1050 }
1051 
1052 void PHPatternReco::set_max_kappa(float kappa_max){
1053  _max_kappa = kappa_max;
1054  _use_max_kappa = true;
1055 }
1056 
1057 void PHPatternReco::set_nbins(unsigned int izoom){
1058  nkappa = _hough_space->get_n_kappa_bins(izoom);
1059  nphi = _hough_space->get_n_phi_bins(izoom);
1060  nd = _hough_space->get_n_d_bins(izoom);
1061  ndzdl = _hough_space->get_n_dzdl_bins(izoom);
1062  nz0 = _hough_space->get_n_z0_bins(izoom);
1063 }
1064 
1065 void PHPatternReco::initialize_houghbin(){
1066  for (ik=0; ik<nkappa; ++ik){
1067  for (ip=0; ip<nphi; ++ip) {
1068  for (id=0;id<nd; ++id) {
1069  for (il=0; il<ndzdl; ++il) {
1070  for (iz=0; iz<nz0; ++iz) {
1071 
1072  unsigned int bins[5]={ik,ip,id,il,iz};
1073 // cout<<"InitZVertexing:: izoom "<<izoom<<" il " <<il<<" iz "<<iz<<endl;
1074  bin = _hough_space->get_bin(0,bins);
1075 // cout<<"InitZVertexing:: bin " << bin<<endl;
1076  HelixHoughBin* _bin = new HelixHoughBin_v1(bin);
1077  _bin->set_zoomlevel(0);
1078  _bin->set_hough_space(_hough_space);
1079  _bin->init();
1080  if (bin == 0) _bin->identify();
1081  bins_map.insert(std::pair<unsigned int, HelixHoughBin*>(bin,_bin));
1082 // cout<<"InitZVertexing:: lbin "<<_bin->get_dzdl_bin(0)<<" zbin "<<_bin->get_z0_bin(0)<<endl;
1083  }
1084  }
1085  }
1086  }
1087  }
1088 
1089  cout<<"bins_map.size " <<bins_map.size()<<endl;
1090 }
1091 
1092 void PHPatternReco::vote_z_init(unsigned int zoomlevel){
1093 
1094  float hitpos3d[3];
1095  unsigned int cluster_id=-999;
1096  std::vector<float> kappa_phi_d_ranges;
1097  float kappa_min = _hough_space->get_kappa_min();
1098  float kappa_max = _hough_space->get_kappa_max();
1099  float phi_min = _hough_space->get_phi_min();
1100  float phi_max = _hough_space->get_phi_max();
1101  float d_min = _hough_space->get_d_min();
1102  float d_max = _hough_space->get_d_max();
1103  kappa_phi_d_ranges.push_back(kappa_min);
1104  kappa_phi_d_ranges.push_back(kappa_max);
1105  kappa_phi_d_ranges.push_back(phi_min);
1106  kappa_phi_d_ranges.push_back(phi_max);
1107  kappa_phi_d_ranges.push_back(d_min);
1108  kappa_phi_d_ranges.push_back(d_max);
1109 
1110  float dzdl_min = _hough_space->get_dzdl_min();
1111  float dzdl_max = _hough_space->get_dzdl_max();
1112  int icluster=0;//test
1113 
1114 // cout<<"kappa min " <<kappa_min<<" kappa max "<<kappa_max<<" phi_min " <<phi_min<<endl;
1115  for (std::map<unsigned int,SimpleHit3D>::iterator it=hits_map.begin();
1116  it!=hits_map.end();
1117  ++it)
1118  {
1119 
1120  cluster_id = it->first;
1121 // ?????????????? What is This ???????????????????
1122 // if (cluster_id != cluster_id) continue;
1123  bool used = false;
1124  auto hitused = hits_used.find(cluster_id);
1125  if (hitused != hits_used.end())
1126  used = hits_used.find(cluster_id)->second;
1127  if(used) continue;
1128  SimpleHit3D hit = it->second;
1129  hitpos3d[0] = hit.get_x();
1130  hitpos3d[1] = hit.get_y();
1131  hitpos3d[2] = hit.get_z();
1132 
1133 // cout<<"x "<< hitpos3d[0]<<" y "<<hitpos3d[1]<<" z "<<hitpos3d[2]<<endl;
1134 
1135  for (iz=0; iz < _hough_space->get_n_z0_bins(zoomlevel); ++iz ) {
1136  std::vector<float> z0_range;
1137  float z0_min = _hough_space->get_z0_min()+iz*_hough_space->get_z0_bin_size(zoomlevel);
1138  float z0_max = _hough_space->get_z0_min()+(iz+1)*_hough_space->get_z0_bin_size(zoomlevel);
1139  z0_range.push_back(z0_min);
1140  z0_range.push_back(z0_max);
1141  float dzdl_range[2] = {4,5};//test
1142  //compute min_dzdl & max_dzdl for given iz
1143  _hough_funcs->set_current_zoom(zoomlevel);
1144 
1145 #ifdef _DEBUG_
1146 // cout<<"set zoom "<<endl;
1147 // cout<<"i "<<i<<", iz "<<iz<<endl;
1148 #endif
1149 
1150  _hough_funcs->calculate_dzdl_range(hitpos3d,z0_range, kappa_phi_d_ranges, dzdl_range);
1151  unsigned int dzdl_bin_range[2];
1152  dzdl_bin_range[0] = _hough_space->get_dzdl_bin(0,dzdl_range[0]);
1153  dzdl_bin_range[1] = _hough_space->get_dzdl_bin(0,dzdl_range[1]);
1154  if (dzdl_range[0] > dzdl_max) continue;
1155  else if (dzdl_range[0] <dzdl_min) dzdl_bin_range[0]=0;
1156  if (dzdl_range[1] < dzdl_min) continue;
1157  else if (dzdl_range[1] >dzdl_max) dzdl_bin_range[1]=_hough_space->get_n_dzdl_bins(zoomlevel)-1;
1158 #ifdef _DEBUG_
1159  cout<<"dzdl_min " <<dzdl_range[0]<<" dzdl_max "<<dzdl_range[1]<<endl;
1160 #endif
1161  for (il =dzdl_bin_range[0]; il<dzdl_bin_range[1]+1; ++il) {
1162  unsigned int zbins[5] = {0,0,0,il,iz};
1163  bin = _hough_space->get_bin(zoomlevel,zbins);
1164 #ifdef _DEBUG_
1165  cout<<"bin "<<bin<<endl;
1166 #endif
1167  auto search = bins_map.find(bin);
1168  if (search != bins_map.end())
1169  {
1170  bins_map.find(bin)->second->add_cluster_ID(cluster_id);
1171  }
1172 #ifdef _DEBUG_
1173  cout<<"count "<<bins_map.find(bin)->second->get_count()<<endl;
1174 #endif
1175  } // bins_map_prev.clear(); // for zoomlevel >0
1176  z0_range.clear();
1177  } //iz
1178 
1179  kappa_phi_d_ranges.clear();
1180  hitpos3d[0]=-999; hitpos3d[1]= -999; hitpos3d[2]=-999;
1181  ++icluster;//test
1182  }
1183  cout<<"total number of clusters "<<icluster<<endl;
1184 
1185 }
1186 
1187 
1188 void PHPatternReco::find_track_candidates_z_init(unsigned int zoomlevel){
1189  bins_map_sel.erase(bins_map_sel.begin(), bins_map_sel.end());
1190  bins_map_sel.clear();
1191 
1192  unsigned int max_i=-1;
1193  unsigned int max_count = 0;
1194 
1195  for(unsigned int i =0; i<_hough_space->get_n_z0_bins(0); ++i){
1196  for(unsigned int j =0; j<_hough_space->get_n_dzdl_bins(0); ++j ){
1197  unsigned int bins[5]={0,0,0,j,i};
1198  bin = _hough_space->get_bin(0,bins);
1199  unsigned int count = bins_map.find(bin)->second->get_count();
1200  if (count <1) continue;
1201  if (count>max_count) {
1202  max_i=i;
1203  max_count=count;
1204  }
1205 #ifdef _DEBUG_
1206  cout<<"bin "<<bin<<endl;
1207  cout<<"z0 bin "<<i<<" dzdl bin "<<j<< " count "<< count<<endl;
1208  cout<<"z0 bin "<< bins_map.find(bin)->second->get_z0_bin(0)<<" dzdl bin "<<bins_map.find(bin)->second->get_dzdl_bin(0)<<endl;
1209 
1210  _z0_dzdl->Fill(bins_map.find(bin)->second->get_z0_center(zoomlevel),bins_map.find(bin)->second->get_dzdl_center(zoomlevel), count);
1211 #endif
1212 // if (count>= _min_nlayers)
1213 // {
1214 // bins_map_sel.insert(std::pair<unsigned int, HelixHoughBin*>(bin,bins_map.find(bin)->second) );
1215 // cout<<"inserting selected bin.. "<<endl;
1216 
1217 // print cluster ids
1218 // for(HelixHoughBin::ClusterIter iter = bins_map.find(bin)->second->begin_clusters();
1219 // iter != bins_map.find(bin)->second->end_clusters();
1220 // iter++){
1221 // cout<<"z0 bin "<< bins_map.find(bin)->second->get_z0_bin(0)<<" dzdl bin "<<bins_map.find(bin)->second->get_dzdl_bin(0)<<endl;
1222 // cout<< "cluster_id "<<*iter<<endl;
1223 // }
1224 // }
1225  }
1226  }
1227 
1228  for (unsigned int k =0; k<_hough_space->get_n_dzdl_bins(0); ++k){
1229  unsigned int bins[5]={0,0,0,k,max_i};
1230  bin = _hough_space->get_bin(0,bins);
1231  unsigned int count = bins_map.find(bin)->second->get_count();
1232  if (count <1) continue;
1233 
1234  if (count>= _min_nlayers)
1235  {
1236  bins_map_sel.insert(std::pair<unsigned int, HelixHoughBin*>(bin,bins_map.find(bin)->second) );
1237  }
1238 
1239  }
1240  cout<<"bins_map_sel.size at zoom "<<zoomlevel << " : (find_track_candidates_z_init) " <<bins_map_sel.size()<<endl;
1241 
1242 }
1243 
1244 
1245 void PHPatternReco::vote_z(unsigned int zoomlevel){
1246 
1247  set_nbins(zoomlevel); // at next zoom, bin number changes
1248 
1249  float hitpos3d[3];
1250  unsigned int cluster_id=-999;
1251  for (std::map<unsigned int,HelixHoughBin*>::iterator it=bins_map_prev.begin(); it!=bins_map_prev.end(); ++it)
1252  {
1253  bin = it->first;
1254  HelixHoughBin* houghbin = it->second;
1255 
1256 #ifdef _DEBUG_
1257 // unsigned int izprev = houghbin->get_z0_bin(zoomlevel-1);
1258 // unsigned int ilprev = houghbin->get_dzdl_bin(zoomlevel-1);
1259 // unsigned int ipprev = houghbin->get_phi_bin(zoomlevel-1);
1260 
1261 // unsigned int ikprev = houghbin->get_kappa_bin(zoomlevel-1);
1262 // unsigned int idprev = houghbin->get_d_bin(0);
1263 // bool fillhisto = (izprev==5) && (ilprev==11)&& (ipprev==21);
1264 // cout<<"bin "<<bin<<" ik " <<ikprev<<" ip "<<ipprev<<" id "<<idprev<<endl;
1265 #endif
1266  for(HelixHoughBin::ClusterIter iter = bins_map_prev.find(bin)->second->begin_clusters();
1267  iter != bins_map_prev.find(bin)->second->end_clusters();
1268  ++iter){
1269  cluster_id = *iter;
1270 
1271  SimpleHit3D hit = hits_map.find(cluster_id)->second;
1272  hitpos3d[0] = hit.get_x();
1273  hitpos3d[1] = hit.get_y();
1274  hitpos3d[2] = hit.get_z();
1275 #ifdef _DEBUG_
1276 // cout<<"cluster_id "<<cluster_id<<endl;
1277 #endif
1278 
1279  std::vector<float> kappa_phi_d_ranges;
1280  float kappa_min = houghbin->get_kappa_center(zoomlevel-1)-0.5*_hough_space->get_kappa_bin_size(zoomlevel-1);
1281  float kappa_max = houghbin->get_kappa_center(zoomlevel-1)+0.5*_hough_space->get_kappa_bin_size(zoomlevel-1);
1282  float phi_min = houghbin->get_phi_center(zoomlevel-1)-0.5*_hough_space->get_phi_bin_size(zoomlevel-1);
1283  float phi_max = houghbin->get_phi_center(zoomlevel-1)+0.5*_hough_space->get_phi_bin_size(zoomlevel-1);
1284  float d_min = houghbin->get_d_center(zoomlevel-1)-0.5*_hough_space->get_d_bin_size(zoomlevel-1);
1285  float d_max = houghbin->get_d_center(zoomlevel-1)+0.5*_hough_space->get_d_bin_size(zoomlevel-1);
1286  kappa_phi_d_ranges.push_back(kappa_min);
1287  kappa_phi_d_ranges.push_back(kappa_max);
1288  kappa_phi_d_ranges.push_back(phi_min);
1289  kappa_phi_d_ranges.push_back(phi_max);
1290  kappa_phi_d_ranges.push_back(d_min);
1291  kappa_phi_d_ranges.push_back(d_max);
1292  float dzdl_min = houghbin->get_dzdl_center(zoomlevel-1)-0.5*_hough_space->get_dzdl_bin_size(zoomlevel-1);
1293  float dzdl_max = houghbin->get_dzdl_center(zoomlevel-1)+0.5*_hough_space->get_dzdl_bin_size(zoomlevel-1);
1294 #ifdef _DEBUG_
1295 // if (fillhisto){
1296 // cout<<"kappa_min "<<kappa_min <<" kappa_max " <<kappa_max<<" phi_min "<<phi_min<<" phi_max "<<phi_max<<endl;
1297 // cout<<"d_min "<<d_min<< " d_max "<<d_max<<" dzdl_min "<<dzdl_min<<" dzdl_max "<<dzdl_max<<endl;
1298 // }
1299 #endif
1300 
1301  for (iz=0; iz < _hough_space->get_n_z0_bins(zoomlevel); ++iz ) {
1302  std::vector<float> z0_range;
1303  float z0_min = houghbin->get_z0_center(zoomlevel-1)-0.5*_hough_space->get_z0_bin_size(zoomlevel-1)+iz*_hough_space->get_z0_bin_size(zoomlevel);
1304  float z0_max = z0_min + _hough_space->get_z0_bin_size(zoomlevel);;
1305  z0_range.push_back(z0_min);
1306  z0_range.push_back(z0_max);
1307  float dzdl_range[2];
1308 #ifdef _DEBUG_
1309 // if (fillhisto) cout<<"z0 min "<<z0_min<<" z0 max "<<z0_max<<endl;
1310 #endif
1311  unsigned int dzdl_bin_range[2];
1312  //compute min_dzdl & max_dzdl for given iz
1313  _hough_funcs->set_current_zoom(zoomlevel);
1314  _hough_funcs->calculate_dzdl_range(hitpos3d,z0_range, kappa_phi_d_ranges, dzdl_range);
1315  dzdl_bin_range[0] = _hough_space->get_dzdl_bin(zoomlevel,dzdl_range[0]);
1316  dzdl_bin_range[1] = _hough_space->get_dzdl_bin(zoomlevel,dzdl_range[1]);
1317 
1318  if (dzdl_range[0] > dzdl_max) continue;
1319  else if (dzdl_range[0] <dzdl_min) dzdl_bin_range[0]=0;
1320  if (dzdl_range[1] < dzdl_min) continue;
1321  else if (dzdl_range[1] >dzdl_max) dzdl_bin_range[1]=_hough_space->get_n_dzdl_bins(zoomlevel)-1;
1322 #ifdef _DEBUG_
1323 // if (fillhisto)cout<<" dzdl_bin_range[0] " <<dzdl_bin_range[0]<<" dzdl_bin_range[1] "<<dzdl_bin_range[1]<<endl;
1324 #endif
1325  for (il =dzdl_bin_range[0]; il<dzdl_bin_range[1]+1; ++il) {
1326  unsigned int zbins[5] = {0,0,0,il,iz};
1327  unsigned int curbin = _hough_space->get_bin(zoomlevel,zbins);
1328  houghbin->set_bin(zoomlevel,curbin);
1329  houghbin->set_bins(zoomlevel,curbin);
1330  unsigned int curgbin = houghbin->get_global_bin(zoomlevel);
1331 #ifdef _DEBUG_
1332 // if (fillhisto) cout<<"global bin "<<curgbin<<" z0 bin "<<houghbin->get_z0_bin(zoomlevel)<<" dzdl bin "<<houghbin->get_dzdl_bin(zoomlevel) <<endl;
1333 #endif
1334  auto search = bins_map_cur.find(curgbin);
1335  if(search != bins_map_cur.end())
1336  {
1337  bins_map_cur.find(curgbin)->second->add_cluster_ID(cluster_id);
1338  }
1339  else
1340  {
1341  HelixHoughBin* cur_bin = houghbin->Clone();
1342  cur_bin->clear_clusters();
1343  cur_bin->set_zoomlevel(zoomlevel);
1344  cur_bin->set_bin(zoomlevel,curbin);
1345  cur_bin->set_bins(zoomlevel,curbin);
1346  cur_bin->add_cluster_ID(cluster_id);
1347  bins_map_cur.insert(std::pair<unsigned int,HelixHoughBin*>(curgbin,cur_bin));
1348 #ifdef _DEBUG_
1349 // if (fillhisto) cout<<"global bin "<<cur_bin->get_global_bin(zoomlevel)<<" z0 bin "<<cur_bin->get_z0_bin(zoomlevel)<<" dzdl bin "<<cur_bin->get_dzdl_bin(zoomlevel) <<endl;
1350 #endif
1351  }
1352 #ifdef _DEBUG_
1353 // cout<<" il "<<il<<" iz "<<iz<<" il from bin"<< bins_map_cur.find(curgbin)->second->get_dzdl_bin(zoomlevel)<<" iz from bin "<< bins_map_cur.find(curgbin)->second->get_z0_bin(zoomlevel)<<endl;
1354 
1355  unsigned int count = bins_map_cur.find(curgbin)->second->get_count();
1356 // cout<<"count "<<count<<endl;
1357 // cout<<"z0 "<<bins_map_cur.find(curgbin)->second->get_z0_center(zoomlevel)<<" dzdl "<< bins_map_cur.find(curgbin)->second->get_dzdl_center(zoomlevel)<<endl;
1358  _z0_dzdl->Fill(bins_map_cur.find(curgbin)->second->get_z0_center(zoomlevel),bins_map_cur.find(curgbin)->second->get_dzdl_center(zoomlevel), count);
1359 #endif
1360  } // il
1361  } //iz
1362  hitpos3d[0]=-999; hitpos3d[1]= -999; hitpos3d[2]=-999;
1363  kappa_phi_d_ranges.clear();
1364  }//clusters
1365  } //bins_map_prev
1366  bins_map_prev.erase(bins_map_prev.begin(),bins_map_prev.end());
1367  bins_map_prev.clear();
1368 
1369  cout<<"bins_map_cur.size at zoom "<<zoomlevel << " (vote_z) : " <<bins_map_cur.size()<<endl;
1370 }
1371 
1372 
1373 void PHPatternReco::find_track_candidates_z(unsigned int zoomlevel){
1374  bins_map_sel.erase(bins_map_sel.begin(),bins_map_sel.end());
1375  bins_map_sel.clear();
1376 
1377 
1378  for (std::map<unsigned int,HelixHoughBin*>::iterator it = bins_map_cur.begin(); it != bins_map_cur.end(); ++it)
1379  {
1380  bin = it->first;
1381  HelixHoughBin* houghbin = it->second;
1382  unsigned int count = houghbin->get_count();
1383  if (count <1) continue;
1384 
1385  il = houghbin->get_dzdl_bin(zoomlevel);
1386  iz = houghbin->get_z0_bin(zoomlevel);
1387 
1388 #ifdef _DEBUG_
1389 // ik = houghbin->get_kappa_bin(zoomlevel);
1390 // ip = houghbin->get_phi_bin(zoomlevel);
1391 // id = houghbin->get_d_bin(zoomlevel);
1392 
1393 // cout<<"bin "<<bin<<endl;
1394 // cout<<"kappa bin "<<ik<<" phi bin "<<ip<< " d bin "<<id <<" il "<<il<<" iz "<<iz<< " count "<< count<<endl;
1395 #endif
1396  if (count>= _min_nlayers)
1397  {
1398  bins_map_sel.insert(std::pair<unsigned int, HelixHoughBin*>(bin,bins_map_cur.find(bin)->second) );
1399 #ifdef _DEBUG_
1400  cout<<"il "<<il<<" iz "<<iz<< " count "<< count<<endl;
1401 #endif
1402  }
1403 
1404  }
1405  bins_map_cur.erase(bins_map_cur.begin(), bins_map_cur.end());
1406  bins_map_cur.clear();
1407 
1408  cout<<"bins_map_sel.size at zoom "<< zoomlevel<<" (find_track_candidates_z) : " <<bins_map_sel.size()<<endl;
1409 }
1410 
1411 void PHPatternReco::vote_xy(unsigned int zoomlevel){
1412 
1413  int count = -1;
1414 #ifdef _DEBUG_
1415  unsigned int ilfill= -1;
1416  unsigned int izfill= -1;
1417 #endif
1418  set_nbins(zoomlevel);
1419  // hand select bins over to bins_map_cur
1420  kappa_phi_map.clear();
1421  d_phi_map.clear();
1422  unsigned int phi_bin_range[2];
1423  for (std::map<unsigned int,HelixHoughBin*>::iterator it=bins_map_sel.begin(); it!=bins_map_sel.end(); ++it)
1424  {
1425  bin = it->first;
1426  HelixHoughBin* houghbin = it->second;
1427  il = houghbin->get_dzdl_bin(zoomlevel);
1428  iz = houghbin->get_z0_bin(zoomlevel);
1429  ++count;
1430 
1431 #ifdef _DEBUG_
1432  if (count==0){ ilfill = il; izfill = iz;}
1433 
1434  unsigned int ikprev =0;
1435  unsigned int ipprev =0;
1436  unsigned int idprev =0;
1437  if (zoomlevel>0){
1438  ikprev = houghbin->get_kappa_bin(zoomlevel-1);
1439  ipprev = houghbin->get_phi_bin(zoomlevel-1);
1440  idprev = houghbin->get_d_bin(zoomlevel-1);
1441  }
1442 
1443  cout<<"bin "<< bin<<endl;
1444  cout<<"il " <<il<<" iz "<<iz<<" ikprev "<<ikprev<<" ipprev "<<ipprev<<" idprev "<<idprev<<" count"<<houghbin->get_count()<<endl;
1445 #endif
1446 
1447  float phi_min;
1448  if (zoomlevel==0) phi_min = _hough_space->get_phi_min();
1449  else
1450  phi_min = houghbin->get_phi_center(zoomlevel-1)-0.5*_hough_space->get_phi_bin_size(zoomlevel-1);
1451  float phi_max;
1452  if (zoomlevel==0) phi_max = _hough_space->get_phi_max();
1453  else
1454  phi_max = phi_min + _hough_space->get_phi_bin_size(zoomlevel-1);
1455 
1456 #ifdef _DEBUG_
1457 // cout<<"phi_min "<<phi_min <<" phi_max "<<phi_max<<endl;
1458 #endif
1459  for (HelixHoughBin::ClusterIter iter = houghbin->begin_clusters();
1460  iter != houghbin->end_clusters(); ++iter){
1461  cluster_id = *iter;
1462 
1463  SimpleHit3D hit = hits_map.find(cluster_id)->second;
1464  hitpos3d[0] = hit.get_x();
1465  hitpos3d[1] = hit.get_y();
1466 #ifdef _DEBUG_
1467  cout<<"cluster id "<<cluster_id<<" x "<< hitpos3d[0]<<" y "<<hitpos3d[1]<<endl;
1468 #endif
1469  for (ik=0; ik<nkappa; ++ik){
1470  for (id=0; id<nd; ++id){
1471  std::vector<float> kappa_d_ranges;
1472 
1473  float kappa_min;
1474  if (zoomlevel == 0) kappa_min = _hough_space->get_kappa_min();
1475  else
1476  kappa_min= houghbin->get_kappa_center(zoomlevel-1)-0.5*_hough_space->get_kappa_bin_size(zoomlevel-1);
1477  kappa_min += ik*_hough_space->get_kappa_bin_size(zoomlevel);
1478  float kappa_max = kappa_min + _hough_space->get_kappa_bin_size(zoomlevel);
1479  float d_min;
1480  if (zoomlevel==0) d_min = _hough_space->get_d_min();
1481  else d_min = houghbin->get_d_center(zoomlevel-1)-0.5*_hough_space->get_d_bin_size(zoomlevel-1);
1482  d_min += id*_hough_space->get_d_bin_size(zoomlevel);
1483  float d_max = d_min + _hough_space->get_d_bin_size(zoomlevel);
1484 
1485  kappa_d_ranges.push_back(kappa_min);
1486  kappa_d_ranges.push_back(kappa_max);
1487  kappa_d_ranges.push_back(d_min);
1488  kappa_d_ranges.push_back(d_max);
1489  _hough_funcs->set_current_zoom(zoomlevel);
1490 
1491  if (!separate_helicity){
1492  float phi_r_range[2];
1493  float phi_l_range[2];
1494  unsigned int phi_r_bin_range[2];
1495  unsigned int phi_l_bin_range[2];
1496  //cout<<"before calc phi "<<endl;
1497  _hough_funcs->calculate_phi_range(hitpos3d,kappa_d_ranges, phi_r_range,phi_l_range);
1498 // cout<<"phi_l_min "<<phi_l_range[0]<<" phi_l_max "<<phi_l_range[1]<<endl;
1499 // cout<<"phi_r_min "<<phi_r_range[1]<<" phi_r_max "<<phi_r_range[1]<<endl;
1500  phi_r_bin_range[0] = _hough_space->get_phi_bin(zoomlevel,phi_r_range[0]);
1501  phi_r_bin_range[1] = _hough_space->get_phi_bin(zoomlevel,phi_r_range[1]);
1502  phi_l_bin_range[0] = _hough_space->get_phi_bin(zoomlevel,phi_l_range[0]);
1503  phi_l_bin_range[1] = _hough_space->get_phi_bin(zoomlevel,phi_l_range[1]);
1504 // cout<<"phi_r_bin_min "<<phi_r_bin_range[0]<<" phi_r_bin_max "<<phi_r_bin_range[1]<<endl;
1505 
1506  bool phi_r_out_of_range=false;
1507  bool phi_l_out_of_range=false;
1508  if (phi_r_range[0] > phi_max) phi_r_out_of_range=true;
1509  else if (phi_r_range[0] <phi_min) phi_r_bin_range[0]=0;
1510  if (phi_r_range[1] < phi_min) phi_r_out_of_range=true;
1511  else if (phi_r_range[1] >phi_max) phi_r_bin_range[1]=_hough_space->get_n_phi_bins(zoomlevel)-1;
1512 
1513  if (phi_l_range[0] > phi_max) phi_r_out_of_range=true;
1514  else if (phi_r_range[0] <phi_min) phi_r_bin_range[0]=0;
1515  if (phi_r_range[1] < phi_min) phi_r_out_of_range=true;
1516  else if (phi_r_range[1] >phi_max) phi_r_bin_range[1]=_hough_space->get_n_phi_bins(zoomlevel)-1;
1517 
1518  for (int i = 0; i<2;++i){// helicity
1519  //cout<<"i "<<i<<endl;
1520  if ((i==0&&phi_r_out_of_range) || (i==1&&phi_l_out_of_range)) continue;
1521  switch (i){
1522  case 0:
1523  phi_bin_range[0]=phi_r_bin_range[0];
1524  phi_bin_range[1]=phi_r_bin_range[1];
1525  break;
1526  case 1:
1527  phi_bin_range[0]=phi_l_bin_range[0];
1528  phi_bin_range[1]=phi_l_bin_range[1];
1529  break;
1530  }
1531 // cout<<"phi_bin_range[0] " <<phi_bin_range[0]<<" phi_bin_range[1] "<<phi_bin_range[1]<<endl;
1532 
1533  for (ip= phi_bin_range[0]; ip<phi_bin_range[1]+1;++ip){
1534  unsigned int xyzbins[5]={ik,ip,id,il,iz};
1535  unsigned int curbin = _hough_space->get_bin(zoomlevel,xyzbins);
1536  houghbin->set_bin(zoomlevel,curbin);
1537  houghbin->set_bins(zoomlevel,curbin);
1538  unsigned int curgbin = houghbin->get_global_bin(zoomlevel);
1539 // cout<<" curgbin "<<curgbin<<endl;
1540  auto search = bins_map_cur.find(curgbin);
1541  if(search != bins_map_cur.end())
1542  {
1543  bins_map_cur.find(curgbin)->second->add_cluster_ID(cluster_id);
1544  }
1545  else
1546  {
1547  HelixHoughBin* cur_bin = houghbin->Clone();
1548  cur_bin->clear_clusters();
1549  cur_bin->set_zoomlevel(zoomlevel);
1550  cur_bin->set_bin(zoomlevel,curbin);
1551  cur_bin->set_bins(zoomlevel,curbin);
1552  cur_bin->add_cluster_ID(cluster_id);
1553  bins_map_cur.insert(std::pair<unsigned int,HelixHoughBin*>(curgbin,cur_bin));
1554 // cout<<" il "<<il<<" iz "<<iz<<" il from bin"<< bins_map_cur.find(curbin)->second->get_dzdl_bin(0)<<" iz from bin "<< bins_map_cur.find(curbin)->second->get_z0_bin(0)<<endl;
1555  }
1556 
1557  auto search0 = kappa_d_phi_map.find(curgbin);
1558  if(search0 != kappa_d_phi_map.end())
1559  {
1560  kappa_d_phi_map.find(curgbin)->second = kappa_d_phi_map.find(curgbin)->second + 1;
1561  }
1562  else
1563  {
1564  kappa_d_phi_map.insert(std::pair<unsigned int,unsigned int>(curgbin,1));
1565  }
1566 
1567  unsigned int kpbins[5]={ik,ip,0,il,iz};
1568  unsigned int kpbin = _hough_space->get_bin(zoomlevel,kpbins);
1569  auto search1 = kappa_phi_map.find(kpbin);
1570  if(search1 != kappa_phi_map.end())
1571  {
1572  kappa_phi_map.find(kpbin)->second = kappa_phi_map.find(kpbin)->second +1;
1573  }
1574  else
1575  {
1576  kappa_phi_map.insert(std::pair<unsigned int,unsigned int>(kpbin,1));
1577  }
1578 
1579  unsigned int dpbins[5]={0,ip,id,il,iz};
1580  unsigned int dpbin = _hough_space->get_bin(zoomlevel,dpbins);
1581  auto search2 = d_phi_map.find(dpbin);
1582  if (search2 != d_phi_map.end())
1583  {
1584  d_phi_map.find(dpbin)->second = d_phi_map.find(dpbin)->second + 1;
1585  }
1586  else
1587  {
1588  d_phi_map.insert(std::pair<unsigned int,unsigned int>(dpbin,1));
1589  }
1590 
1591 #ifdef _DEBUG_
1592  if (il ==ilfill && iz == izfill){
1593 // for(HelixHoughBin::ClusterIter iter = bins_map.find(bin)->second->begin_clusters();
1594 // iter != bins_map.find(bin)->second->end_clusters();
1595 // ++iter){
1596 // cout<< "cluster_id "<<*iter<<endl;
1597 // }
1598 
1599  cout << "il "<<il<<" iz "<<iz<<endl;
1600  _kappa_phi->Fill(bins_map_cur.find(curgbin)->second->get_kappa_center(zoomlevel),bins_map_cur.find(curgbin)->second->get_phi_center(zoomlevel),1);
1601  _d_phi->Fill(bins_map_cur.find(curgbin)->second->get_d_center(zoomlevel),bins_map_cur.find(curgbin)->second->get_phi_center(zoomlevel),1);
1602  _kappa_d_phi->Fill(bins_map_cur.find(curgbin)->second->get_kappa_center(zoomlevel),bins_map_cur.find(curgbin)->second->get_d_center(zoomlevel),bins_map_cur.find(curgbin)->second->get_phi_center(zoomlevel),1);
1603  }
1604 #endif
1605  }//ip
1606  }// helicity
1607 
1608  }// if not separating helicities
1609  else
1610  {
1611  float phi_range[2];
1612  float phi_prev_range[2];
1613  float phi_next_range[2];
1614  unsigned int phi_bin_range[2];
1615  //cout<<"before calc phi "<<endl;
1616 
1617  if (id==0)
1618  _hough_funcs->calculate_phi_range(hitpos3d,kappa_d_ranges, helicity, phi_range,phi_next_range);
1619  else
1620  _hough_funcs->calculate_phi_range(hitpos3d,kappa_d_ranges, helicity, phi_range,phi_prev_range,phi_next_range);
1621 
1622  phi_prev_range[0]= phi_next_range[0];
1623  phi_prev_range[1]= phi_next_range[1];
1624 
1625  //cout<<"phi_min "<<phi_range[0]<<" phi_max "<<phi_range[1]<<endl;
1626  phi_bin_range[0] = _hough_space->get_phi_bin(zoomlevel,phi_range[0]);
1627  phi_bin_range[1] = _hough_space->get_phi_bin(zoomlevel,phi_range[1]);
1628  //cout<<"phi_bin_min "<<phi_bin_range[0]<<" phi_bin_max "<<phi_bin_range[1]<<endl;
1629 
1630  bool phi_out_of_range=false;
1631  if (phi_range[0] > phi_max) phi_out_of_range=true;
1632  else if (phi_range[0] <phi_min) phi_bin_range[0]=0;
1633  if (phi_range[1] < phi_min) phi_out_of_range=true;
1634  else if (phi_range[1] >phi_max) phi_bin_range[1]=_hough_space->get_n_phi_bins(zoomlevel)-1;
1635 
1636  if (phi_out_of_range) continue;
1637 
1638  for (ip= phi_bin_range[0]; ip<phi_bin_range[1]+1;++ip){
1639  unsigned int xyzbins[5]={ik,ip,id,il,iz};
1640  unsigned int curbin = _hough_space->get_bin(zoomlevel,xyzbins);
1641  houghbin->set_bin(zoomlevel,curbin);
1642  houghbin->set_bins(zoomlevel,curbin);
1643  unsigned int curgbin = houghbin->get_global_bin(zoomlevel);
1644 #ifdef _DEBUG_
1645  cout<<" curgbin "<<curgbin<<endl;
1646 #endif
1647  auto search = bins_map_cur.find(curgbin);
1648  if(search != bins_map_cur.end())
1649  {
1650  bins_map_cur.find(curgbin)->second->add_cluster_ID(cluster_id);
1651  }
1652  else
1653  {
1654  HelixHoughBin* cur_bin = houghbin->Clone();
1655  cur_bin->clear_clusters();
1656  cur_bin->set_zoomlevel(zoomlevel);
1657  cur_bin->set_bin(zoomlevel,curbin);
1658  cur_bin->set_bins(zoomlevel,curbin);
1659  cur_bin->add_cluster_ID(cluster_id);
1660  bins_map_cur.insert(std::pair<unsigned int,HelixHoughBin*>(curgbin,cur_bin));
1661 #ifdef _DEBUG_
1662  cout<<" il "<<il<<" iz "<<iz<<" il from bin"<< bins_map_cur.find(curbin)->second->get_dzdl_bin(0)<<" iz from bin "<< bins_map_cur.find(curbin)->second->get_z0_bin(0)<<endl;
1663 #endif
1664  }
1665 
1666 #ifdef _DEBUG_
1667 // if (il ==ilfill && iz == izfill){
1668 
1669 // for(HelixHoughBin::ClusterIter iter = bins_map.find(bin)->second->begin_clusters();
1670 // iter != bins_map.find(bin)->second->end_clusters();
1671 // ++iter){
1672 // cout<< "cluster_id "<<*iter<<endl;
1673 // }
1674 // cout << "il "<<il<<" iz "<<iz<<endl;
1675 
1676 // _kappa_phi->Fill(bins_map_cur.find(curgbin)->second->get_kappa_center(zoomlevel),bins_map_cur.find(curgbin)->second->get_phi_center(zoomlevel),1);
1677 // _d_phi->Fill(bins_map_cur.find(curgbin)->second->get_d_center(zoomlevel),bins_map_cur.find(curgbin)->second->get_phi_center(zoomlevel),1);
1678 // _kappa_d_phi->Fill(bins_map_cur.find(curgbin)->second->get_kappa_center(zoomlevel),bins_map_cur.find(curgbin)->second->get_d_center(zoomlevel),bins_map_cur.find(curgbin)->second->get_phi_center(zoomlevel),1);
1679 // }
1680 #endif
1681  }//ip
1682 
1683  }// if separating helicities
1684 
1685  kappa_d_ranges.clear();
1686  }//id
1687  }//ik bins
1688  }//cluster
1689 
1690  }//binsmap_sel
1691  cout<<"bins_map_cur.size at zoom "<<zoomlevel <<" (vote_xy) : " <<bins_map_cur.size()<<endl;
1692  bins_map_sel.erase(bins_map_sel.begin(), bins_map_sel.end());
1693  bins_map_sel.clear();
1694 
1695 }
1696 
1697 void PHPatternReco::find_track_candidates_xy(unsigned int zoomlevel){
1698 
1699  bins_map_prev.erase(bins_map_prev.begin(), bins_map_prev.end());
1700  bins_map_prev.clear();
1701 
1702  for (std::map<unsigned int,HelixHoughBin*>::iterator it=bins_map_cur.begin(); it!=bins_map_cur.end(); ++it)
1703  {
1704  bin = it->first;
1705  HelixHoughBin* houghbin = it->second;
1706  unsigned int count = houghbin->get_count();
1707  if (count <1) continue;
1708 
1709  il = houghbin->get_dzdl_bin(zoomlevel);
1710  iz = houghbin->get_z0_bin(zoomlevel);
1711  ik = houghbin->get_kappa_bin(zoomlevel);
1712  ip = houghbin->get_phi_bin(zoomlevel);
1713  id = houghbin->get_d_bin(zoomlevel);
1714 
1715 #ifdef _DEBUG_
1716  cout<<"bin "<<bin<<endl;
1717  cout<<"kappa bin "<<ik<<" phi bin "<<ip<< " d bin "<<id <<" il "<<il<<" iz "<<iz<< " count "<< count<<endl;
1718 #endif
1719  if (count >= _min_nlayers
1720 //&& count >= cnt_kappa_up_d_phi && count >= cnt_kappa_dw_d_phi && count >= cnt_kappa_d_dw_phi && count >= cnt_kappa_d_up_phi && cnt_kappa_d_phi >= cnt_kappa_d_phi_up && count >= cnt_kappa_d_phi_dw
1721  )
1722  {
1723  bins_map_prev.insert(std::pair<unsigned int, HelixHoughBin*>(bin,bins_map_cur.find(bin)->second) );
1724 #ifdef _DEBUG_
1725  cout<<"inserting selected bin.. "<<endl;
1726  cout<<"bin "<<bin<<endl;
1727  cout<<"kappa bin "<<ik<<" phi bin "<<ip<< " d bin "<<id <<" il "<<il<<" iz "<<iz<< " count "<< count<<endl;
1728 #endif
1729  }
1730 
1731 
1732  }
1733 
1734  cout<<"bins_map_prev.size at zoom "<<zoomlevel <<" (find_track_candidates_xy) " <<bins_map_prev.size()<<endl;
1735  bins_map_cur.erase(bins_map_cur.begin(), bins_map_cur.end());
1736  bins_map_cur.clear();
1737 
1738 }
1739 
1740 void PHPatternReco::prune_z(unsigned int zoomlevel){
1741  set_nbins(zoomlevel);
1742  for (std::map<unsigned int,HelixHoughBin*>::iterator it=bins_map_sel.begin(); it!=bins_map_sel.end(); ++it){
1743  HelixHoughBin* houghbin = it->second;
1744  unsigned int nclust = houghbin->get_count();
1745  unsigned int globalbin = houghbin->get_global_bin(zoomlevel);
1746 #ifdef _DEBUG_
1747  cout<<"global bin "<<houghbin->get_global_bin(zoomlevel)<<endl;
1748 
1749  if (zoomlevel>0)
1750  cout<<"ikprev "<<houghbin->get_kappa_bin(zoomlevel-1)<<" ipprev "<<houghbin->get_phi_bin(zoomlevel-1)<<" idprev "<< houghbin->get_d_bin(zoomlevel-1) <<" ilprev "<< houghbin->get_dzdl_bin(zoomlevel-1)<<" izprev "<<houghbin->get_z0_bin(zoomlevel-1)<<endl;
1751  cout<<"il "<<houghbin->get_dzdl_bin(zoomlevel)<<" iz "<<houghbin->get_z0_bin(zoomlevel)<<" count "<<houghbin->get_count()<<endl;
1752 #endif
1753  unsigned int var[3] = {1<<3,1<<4,(1<<3)+(1<<4)};
1754  unsigned int sign[3][4] ={{0<<3,1<<3,0,0},
1755  {0<<4,1<<4,0,0},
1756  {(0<<3)+(0<<4),(0<<3)+(1<<4),(1<<3)+(0<<4),(1<<3)+(1<<4)}};
1757  for (unsigned int lz =0; lz<3; ++lz){
1758  for (unsigned int ss = 0; ss<4; ++ss ){
1759  if (lz<2 && ss>1) continue;
1760  unsigned int neighborbin = houghbin->get_neighbors_global_bin(zoomlevel,var[lz],sign[lz][ss]);
1761 #ifdef _DEBUG_
1762  cout<<"var "<<var[lz]<<" ss "<<ss<<" neighbor_global bin "<<neighborbin<<endl;
1763 #endif
1764  if (globalbin==neighborbin) continue;
1765  auto search = bins_map_sel.find(neighborbin);
1766  if (search != bins_map_sel.end()){
1767 #ifdef _DEBUG_
1768  cout<<"neighbor bin found"<<endl;
1769 #endif
1770  unsigned int nclust_neighbor = bins_map_sel.find(neighborbin)->second->get_count();
1771  if(nclust >= nclust_neighbor && nclust <= 1.1*_nlayers) {
1772 #ifdef _DEBUG_
1773  cout<<"merging "<<endl;
1774 #endif
1775  for(HelixHoughBin::ClusterIter iter = bins_map_sel.find(neighborbin)->second->begin_clusters();
1776  iter != bins_map_sel.find(neighborbin)->second->end_clusters();
1777  ++iter){
1778  houghbin->add_cluster_ID(*iter);
1779  }
1780  }//nclust
1781  }//search
1782  bins_map_sel.erase(neighborbin);
1783  }//ss
1784  }//lz
1785 
1786 // }//i
1787 
1788  }
1789 
1790  cout<<"bins_map_sel.size at zoom " <<zoomlevel<<" (prune_z) : " <<bins_map_sel.size()<<endl;
1791 }
1792 
1793 void PHPatternReco::prune_xy(unsigned int zoomlevel){
1794 
1795  set_nbins(zoomlevel);
1796  cout<<"bins_map_prev.size at zoom " <<zoomlevel<<" (pre prune_xy) : " <<bins_map_prev.size()<<endl;
1797  for (std::map<unsigned int,HelixHoughBin*>::iterator it=bins_map_prev.begin(); it!=bins_map_prev.end(); ++it){
1798  HelixHoughBin* houghbin = it->second;
1799  unsigned int nclust = houghbin->get_count();
1800  unsigned int globalbin = houghbin->get_global_bin(zoomlevel);
1801 #ifdef _DEBUG_
1802  cout<<"global bin "<<globalbin<<endl;
1803 
1804  if (zoomlevel>0)
1805  cout<<"ikprev "<<houghbin->get_kappa_bin(zoomlevel-1)<<" ipprev "<<houghbin->get_phi_bin(zoomlevel-1)<<" idprev "<< houghbin->get_d_bin(zoomlevel-1) <<" ilprev "<< houghbin->get_dzdl_bin(zoomlevel-1)<<" izprev "<<houghbin->get_z0_bin(zoomlevel-1)<<endl;
1806  cout<<"ik "<<houghbin->get_kappa_bin(zoomlevel)<<" ip "<<houghbin->get_phi_bin(zoomlevel)<<" id "<<houghbin->get_d_bin(zoomlevel)<<" il "<<houghbin->get_dzdl_bin(zoomlevel)<<" iz "<<houghbin->get_z0_bin(zoomlevel)<<" count "<<houghbin->get_count()<<endl;
1807 #endif
1808  // k, p, d, kd, kp, dp
1809  unsigned int var[6] = {1,1<<1,1<<2,(1)+(1<<1),(1)+(1<<2),(1<<1)+(1<<2)};
1810  unsigned int sign[6][4]={{0,1,0,0},
1811  {0<<1,1<<1,0,0},
1812  {0<<2,1<<2,0,0},
1813  {0,1<<1,1,(1)+(1<<1)},
1814  {0,1<<2,1, (1)+(1<<2)},
1815  {0,1<<2,1<<1,(1<<1)+(1<<2)}};
1816  for (unsigned int kpd =0; kpd<6; ++kpd){
1817  for (unsigned int ss = 0; ss<4; ++ss ){
1818  if (kpd<3 && ss>1) continue;
1819  unsigned int neighborbin = houghbin->get_neighbors_global_bin(zoomlevel,var[kpd],sign[kpd][ss]);
1820 // cout<<"var "<<var[kpd]<<" ss "<<ss<<" neighbor_global bin "<<neighborbin<<endl;
1821  if (globalbin==neighborbin) continue;
1822  auto search = bins_map_prev.find(neighborbin);
1823  if (search != bins_map_prev.end()){
1824 // cout<<"neighbor bin found"<<endl;
1825  unsigned int nclust_neighbor = bins_map_prev.find(neighborbin)->second->get_count();
1826  if(nclust >= nclust_neighbor && nclust <= 1.1*_nlayers) {
1827 #ifdef _DEBUG_
1828  cout<<"merging "<<endl;
1829 #endif
1830  for(HelixHoughBin::ClusterIter iter = bins_map_prev.find(neighborbin)->second->begin_clusters();
1831  iter != bins_map_prev.find(neighborbin)->second->end_clusters();
1832  ++iter){
1833  houghbin->add_cluster_ID(*iter);
1834  }
1835  }//nclust
1836  }//search
1837  bins_map_prev.erase(neighborbin);
1838  }//ss
1839  }//lz
1840 
1841 // }//i
1842 
1843  }
1844 
1845  cout<<"bins_map_prev.size at zoom " <<zoomlevel<<" (prune_xy) : " <<bins_map_prev.size()<<endl;
1846 
1847 }
1848 
1849 void PHPatternReco::reset_hits_used(){
1850 
1851  for (std::map<unsigned int,bool>::iterator it = hits_used.begin();
1852  it != hits_used.end(); ++it){
1853  it->second = false;
1854  }
1855 }
1856 
1857 void PHPatternReco::bins_to_SimpleTrack3D(std::vector<SimpleTrack3D>& new_tracks, int imap, unsigned int zoomlevel){
1858 
1859  unsigned int icluster = 0;
1860  switch (imap){
1861  case 0: // not implemented yet
1862  for (std::map<unsigned int,HelixHoughBin*>::iterator it = bins_map_sel.begin(); it != bins_map_sel.end(); ++it){
1863  HelixHoughBin *houghbin = it->second;
1864  SimpleTrack3D track;
1865  icluster = 0;
1866  for(HelixHoughBin::ClusterIter iter = houghbin->begin_clusters();
1867  iter != houghbin->end_clusters();
1868  ++iter){
1869  //cout<< "cluster_id "<<*iter<<endl;
1870  SimpleHit3D hit = hits_map.find(*iter)->second;
1871  track.hits.push_back(hit);
1872  track.hits[icluster].set_id(*iter);
1873  ++icluster;
1874  }
1875  _tracks.push_back(track);
1876  }
1877  break;
1878 
1879  case 1:
1880  cout<<"bins_map_prev selected, total bins : "<< bins_map_prev.size()<<endl;
1881  for (std::map<unsigned int,HelixHoughBin*>::iterator it = bins_map_prev.begin(); it != bins_map_prev.end(); ++it){
1882  HelixHoughBin *houghbin = it->second;
1883  SimpleTrack3D track;
1884  icluster = 0;
1885  track.hits.assign(houghbin->get_count(),SimpleHit3D());
1886 #ifdef _DEBUG_
1887 // cout<<"bin "<< houghbin->get_global_bin(zoomlevel);
1888 // cout<<" : start loop over clusters "<<endl;
1889 #endif
1890  for(HelixHoughBin::ClusterIter iter = houghbin->begin_clusters();
1891  iter != houghbin->end_clusters();
1892  ++iter){
1893 #ifdef _DEBUG_
1894 // cout <<"cluster_id "<< *iter<<endl;
1895 #endif
1896  auto search = hits_map.find(*iter);
1897  if (search != hits_map.end())
1898  {
1899  SimpleHit3D hit = hits_map.find(*iter)->second;
1900  track.hits[icluster] = hit;
1901  track.hits[icluster].set_id(*iter);
1902  track.kappa = houghbin->get_kappa_center(zoomlevel);
1903  track.phi = houghbin->get_phi_center(zoomlevel);
1904  track.d = 0.;
1905  track.dzdl = houghbin->get_dzdl_center(zoomlevel);
1906  track.z0 = houghbin->get_z0_center(zoomlevel);
1907  ++icluster;
1908  }
1909  }
1910 // _tracks.push_back(track);// test
1911  new_tracks.push_back(track);
1912  }
1913  break;
1914  }
1915  cout<< "Number of tracks added : (to be used for initial vertexing or track seeding)"<<new_tracks.size()<<endl;
1916 
1917 }
1918 
1919 int PHPatternReco::build_triplets_to_SimpleTrack3D(std::vector<SimpleTrack3D>& new_tracks, bool forward){
1920 
1921  unsigned int layer0 = 0;
1922  unsigned int layer1 = 1;
1923  unsigned int layer2 = 2;
1924  if (!forward) {
1925  layer0 = _nlayers -1;
1926  layer1 = _nlayers -2;
1927  layer2 = _nlayers -3;
1928  }
1929 
1930  std::vector<unsigned int> layer_sorted_0;
1931  std::vector<unsigned int> layer_sorted_1;
1932  std::vector<unsigned int> layer_sorted_2;
1933 
1934  for (std::map<unsigned int,SimpleHit3D>::iterator it=hits_map.begin();
1935  it!=hits_map.end();
1936  ++it)
1937  {
1938  cluster_id = it->first;
1939  if (cluster_id != cluster_id) continue;
1940  bool used = false;
1941  auto hitused = hits_used.find(cluster_id);
1942  if (hitused != hits_used.end())
1943  used = hits_used.find(cluster_id)->second;
1944  if(used) continue;
1945 
1946  SimpleHit3D hit = it->second;
1947  unsigned int layer = hit.get_layer();
1948  unsigned int hitid = hit.get_id();
1949  if (hitid != hitid) continue;
1950  //cout<<"cluster_id "<< it->first <<" layer "<< layer<<" hitid "<<hitid<<endl;
1951  if (layer == layer0 ) layer_sorted_0.push_back(hitid);
1952  else if (layer == layer1) layer_sorted_1.push_back(hitid);
1953  else if (layer == layer2) layer_sorted_2.push_back(hitid);
1954  else continue;
1955  }
1956 
1957  cout<<"layer 0 " <<layer_sorted_0.size()<<endl;
1958  cout<<"layer 1 " <<layer_sorted_1.size()<<endl;
1959  cout<<"layer 2 " <<layer_sorted_2.size()<<endl;
1960 
1961 
1962  std::set<unsigned int> clusters;
1963  bool fill_track = false;
1964  unsigned int cluster_layer0 = 0;
1965  for (std::vector<unsigned int>::iterator it0 = layer_sorted_0.begin() ; it0 != layer_sorted_0.end(); ++it0)
1966  {
1967 #ifdef _DEBUG_
1968 // cout<<"icluster layer 0 "<<cluster_layer0<<endl;
1969 #endif
1970  ++cluster_layer0;
1971  clusters.clear();
1972  fill_track = false;
1973 // cout<<"cluster_id for hit 0 "<< *it0<<endl;
1974  auto search0 = hits_map.find(*it0);
1975  if (search0 == hits_map.end()) continue;
1976  SimpleHit3D hit0 = hits_map.find(*it0)->second;
1977  float phi_h0 = atan2(hit0.get_y(),hit0.get_x());
1978  phi_h0 = shift_phi_range(phi_h0);
1979  float z_h0 = hit0.get_z();
1980 
1981  // one track candidate for each hit on the first layer
1982  SimpleTrack3D track;
1983  clusters.insert(*it0);
1984 
1985  for (std::vector<unsigned int>::iterator it1 = layer_sorted_1.begin(); it1 != layer_sorted_1.end(); ++it1 )
1986  {
1987  SimpleHit3D hit1 = hits_map.find(*it1)->second;
1988  float phi_h1 = atan2(hit1.get_y(),hit1.get_x());
1989  phi_h1 = shift_phi_range(phi_h1);
1990  float z_h1 = hit1.get_z();
1991 
1992  float phi_h0h1 = phi_h1-phi_h0;
1993  if (phi_h1< M_PI/2. && phi_h0 > 3*M_PI/2.) phi_h0h1 += 2.*M_PI;
1994  else if (phi_h1> 3*M_PI/2. && phi_h0 <M_PI/2.) phi_h0h1 -= 2.*M_PI;
1995 
1996  if (fabs(phi_h0h1)> _ca_phi_cut || fabs(z_h1-z_h0) > _z_cut) continue;
1997 
1998 
1999  for(std::vector<unsigned int>::iterator it2 = layer_sorted_2.begin(); it2 != layer_sorted_2.end(); ++it2)
2000  {
2001  SimpleHit3D hit2 = hits_map.find(*it2)->second;
2002  float phi_h2 = atan2(hit2.get_y(),hit2.get_x());
2003  phi_h2 = shift_phi_range(phi_h2);
2004  float phi_h1h2 = phi_h2-phi_h1;
2005  if (phi_h2< M_PI/2. && phi_h1 > 3*M_PI/2.) phi_h1h2 += 2.*M_PI;
2006  else if (phi_h2> 3*M_PI/2. && phi_h1 <M_PI/2.) phi_h1h2 -= 2.*M_PI;
2007 
2008 
2009  float z_h2 = hit2.get_z();
2010  if (fabs(phi_h1h2) < _ca_phi_cut && (z_h1-z_h0)*(z_h2-z_h1)>0 && fabs(z_h2-z_h1)<_z_cut ){
2011 
2012  clusters.insert(*it1);
2013  clusters.insert(*it2);
2014  fill_track = true;
2015  }
2016  } // layer 2
2017  } // layer 1
2018 
2019  if (fill_track) {
2020 // cout<<" fill_track "<<endl;
2021  unsigned int nclusters =0 ;
2022  track.hits.assign(clusters.size(),SimpleHit3D());
2023  for (std::set<unsigned int>::iterator it=clusters.begin(); it!=clusters.end(); ++it ){
2024  //cout<<"cluster id "<<*it<<endl;
2025  SimpleHit3D hit = hits_map.find(*it)->second;
2026  track.hits[nclusters] = hit;
2027  track.hits[nclusters].set_id(*it);
2028  ++nclusters;
2029  }
2030  new_tracks.push_back(track);
2031  }
2032  }// layer 0
2033 
2034 
2035  cout<<"number of triplets "<<new_tracks.size()<<endl;
2036  return 1;
2037 }
2038 
2039 int PHPatternReco::turnoff_hits_used_in_triplets(){
2040 
2041  // loop over all triplets
2042  for (SvtxTrackMap::Iter iter = _trackmap->begin();
2043  iter != _trackmap->end();
2044  ++iter) {
2045  SvtxTrack* track = iter->second;
2046  for (SvtxTrack::ConstClusterIter iter = track->begin_clusters();
2047  iter != track->end_clusters();
2048  ++iter) {
2049  unsigned int cluster_id = *iter;
2050  // SvtxCluster* cluster = _clustermap->get(cluster_id);
2051  // Turn off hits used in triplets ; hits_used
2052  auto hitused = hits_used.find(cluster_id);
2053  if (hitused != hits_used.end())
2054  hits_used.find(cluster_id)->second = true;
2055  }
2056  }
2057  return 1;
2058 
2059 }
2060 
2061 int PHPatternReco::cellular_automaton_zvtx_init(std::vector<SimpleTrack3D>& candidate_tracks){
2062 
2063  cout<<"Entering cellular autumaton : processing "<< candidate_tracks.size()<<" tracks. "<<endl;
2064  std::unique_ptr<CellularAutomaton_v1> ca(new CellularAutomaton_v1(candidate_tracks,_radii,_material));
2065  ca->set_hough_space(_hough_space);
2066  ca->set_mag_field(_mag_field);
2067  ca->set_pt_rescale(_pt_rescale);
2068  ca->set_remove_hits(true);
2069 // ca->set_propagate_forward(false);// need to implement triplet in forward propagation
2070  ca->set_propagate_forward(true);
2071 
2072 // ca->set_mode(0);
2073  ca->set_triplet_mode(true); // triplet
2074  ca->set_seeding_mode(true);
2075  ca->set_hits_map(hits_map);
2076 
2077  ca->set_remove_inner_hits(true);
2078  ca->set_n_layers(_ca_nlayers);
2079  ca->set_required_layers(_ca_nlayers);
2080  ca->set_ca_chi2_layer(2.0);// 1.0
2081  ca->set_ca_chi2(_ca_chi2);
2082  ca->set_ca_phi_cut(_ca_phi_cut);
2083  ca->set_ca_z_cut(_ca_z_cut);
2084  ca->set_ca_dcaxy_cut(_dcaxy_cut);
2085  int code = ca->run(_tracks, _track_states, hits_used); // push_back output tracks into _tracks
2086  ca->Reset();
2087  for(unsigned int i=0; i<candidate_tracks.size(); ++i) candidate_tracks[i].reset();
2088  candidate_tracks.clear();
2089  if (!code)
2090  return 0;
2091 
2092  return Fun4AllReturnCodes::EVENT_OK;
2093 
2094 }
2095 
2096 /*
2097 int PHPatternReco::cellular_automaton_zvtx_second(std::vector<SimpleTrack3D>& candidate_tracks){
2098 
2099 
2100  cout<<"Entering cellular autumaton zvtx_second : processing "<< candidate_tracks.size()<<" tracks. "<<endl;
2101  ca = new CellularAutomaton_v1(candidate_tracks,_radii,_material);
2102  ca->set_hough_space(_hough_space);
2103  ca->set_mag_field(_mag_field);
2104  ca->set_pt_rescale(_pt_rescale);
2105  ca->set_n_layers(_ca_nlayers);
2106  ca->set_required_layers(_ca_min_nlayers);
2107  ca->set_ca_chi2(_ca_chi2);
2108  ca->set_ca_chi2_layer(2.);
2109  ca->set_propagate_forward(true);
2110  ca->set_mode(0);
2111  ca->set_remove_hits(true);
2112  ca->set_remove_inner_hits(true);
2113  ca->set_require_inner_hits(false);//>1
2114  int code = ca->run(_tracks,_track_states, hits_used); // push_back output tracks into _tracks
2115  ca->Reset();
2116  for(unsigned int i=0; i<candidate_tracks.size(); ++i) candidate_tracks[i].reset();
2117  candidate_tracks.clear();
2118  if (!code) return 0;
2119  return Fun4AllReturnCodes::EVENT_OK;
2120 
2121 }
2122 
2123 int PHPatternReco::cellular_automaton_zvtx_third(std::vector<SimpleTrack3D>& candidate_tracks){
2124 
2125  cout<<"Entering cellular autumaton zvtx_third : processing "<< candidate_tracks.size()<<" tracks. "<<endl;
2126  ca = new CellularAutomaton_v1(candidate_tracks,_radii,_material);
2127  ca->set_hough_space(_hough_space);
2128  ca->set_mag_field(_mag_field);
2129  ca->set_pt_rescale(_pt_rescale);
2130  ca->set_n_layers(_ca_nlayers);
2131  ca->set_required_layers(_ca_min_nlayers);
2132  ca->set_ca_chi2(_ca_chi2);
2133  ca->set_ca_chi2_layer(2.);
2134  ca->set_propagate_forward(true);
2135  ca->set_remove_hits(true);
2136  ca->set_remove_inner_hits(true);
2137  int code = ca->run(_tracks,_track_states, hits_used); // push_back output tracks into _tracks
2138  ca->Reset();
2139  for(unsigned int i=0; i<candidate_tracks.size(); ++i) candidate_tracks[i].reset();
2140  candidate_tracks.clear();
2141  if (!code) return 0;
2142  return Fun4AllReturnCodes::EVENT_OK;
2143 
2144 }
2145 
2146 */
2147 
2148 int PHPatternReco::fit_vertex(){
2149  cout<<"tracks size " << _tracks.size()<<endl;
2150 
2151  if (_tracks.empty()) return -1;
2152  std::vector<SimpleTrack3D> vtx_tracks;
2153  vtx_tracks.clear();
2154  _tracks.swap(vtx_tracks);
2155 
2156  unsigned int nzvtx = vtx_tracks.size();
2157 
2158  // range (-32, 32)
2159  unsigned int nzbins= 2*1280;
2160  float binsize = (_max_z0-_min_z0)/nzbins;
2161  float zvalues[2*1280];
2162  unsigned int zcounts[2*1280];
2163  for (unsigned int i = 0; i<nzbins; ++i)
2164  {
2165  zvalues[i] = _min_z0 + (i+0.5)*binsize;
2166  zcounts[i] = 0;
2167  }
2168 
2169  for (unsigned int i = 0; i < nzvtx; ++i)
2170  {
2171  double zvtx = vtx_tracks[i].z0;
2172  if (zvtx != zvtx || zvtx<_min_z0 || zvtx>_max_z0) continue;
2173  unsigned int zbin = (zvtx-_min_z0)/binsize;
2174  cout<<"zvtx " <<zvtx <<" zbin "<< zbin <<endl;
2175  ++zcounts[zbin];
2176  }
2177 
2178  std::vector<float> zvertices; // peak bin position
2179  std::vector<float> zmeans; // averaged z for tracks within 500um window
2180  std::vector<float> zsums;
2181  std::vector<float> zsigmas;
2182  std::vector<unsigned int> nzvtxes;
2183 
2184  for (unsigned int j=2; j<(nzbins-3); ++j)
2185  {
2186  cout<<"z countz "<<j<<" "<<zcounts[j]<<endl;
2187 // bool threebinspeak = _mult_threebins*(zcounts[j-2]) < (zcounts[j-1]+zcounts[j]+zcounts[j+1])
2188 // && _mult_threebins*(zcounts[j+2]) < (zcounts[j-1]+zcounts[j]+zcounts[j+1]);
2189  bool twobinspeak = _mult_twobins*(zcounts[j-1]+zcounts[j-2])<(zcounts[j]+zcounts[j+1])
2190  && _mult_twobins*(zcounts[j+2]+zcounts[j+3])< (zcounts[j]+zcounts[j+1]);
2191  bool onebinpeak = _mult_onebin*(zcounts[j-1])<zcounts[j] && _mult_onebin*(zcounts[j+1]< zcounts[j]);
2192 
2193 
2194  if ((zcounts[j]>=_min_zvtx_tracks && onebinpeak) || ( (zcounts[j]+zcounts[j+1])>= _min_zvtx_tracks && twobinspeak)
2195  /*|| ((zcounts[j-1]+zcounts[j]+zcounts[j+1]) > _min_zvtx_tracks && threebinspeak)*/ ){
2196  float zvertex=-999.;
2197 // if (threebinspeak) zvertex = (zvalues[j-1]*zcounts[j-1] + zvalues[j]*zcounts[j] + zvalues[j+1]*zcounts[j+1])
2198 // /(zcounts[j-1]+zcounts[j]+zcounts[j+1]);
2199 
2200  if (onebinpeak) zvertex = zvalues[j];
2201  if (twobinspeak) zvertex = (zvalues[j]*zcounts[j] + zvalues[j+1]*zcounts[j+1])/(zcounts[j]+zcounts[j+1]);
2202  zvertices.push_back(zvertex);
2203  nzvtxes.push_back(0);
2204  zmeans.push_back(0.0);
2205  zsums.push_back(0.0);
2206  zsigmas.push_back(0.0);
2207  cout<<"vertex "<< zvertex<<endl;
2208  }
2209  }
2210 
2211  cout<<"number of candidate z-vertices : "<< zvertices.size()<<endl;
2212  for (unsigned int j = 0; j <zvertices.size(); ++j)
2213  cout<<"primary vertex position : "<<zvertices[j]<<endl;
2214 
2215  unsigned int izvtx= 999;
2216  for (unsigned int i = 0; i < nzvtx; ++i) {
2217  double zvtx = vtx_tracks[i].z0;
2218  if (zvtx != zvtx) continue;
2219  bool pass_dcaz= false;
2220  for (unsigned int k = 0; k<zvertices.size(); ++k) {
2221  bool new_vtx = fabs(zvtx-zvertices[k]) < _dcaz_cut;
2222  if (new_vtx) izvtx=k;
2223  pass_dcaz = pass_dcaz || new_vtx;
2224  }
2225 
2226  if (!pass_dcaz || izvtx>99) continue;
2227  ++nzvtxes[izvtx];
2228  zsums[izvtx] += zvtx;
2229  }
2230 
2231  for (unsigned int iz = 0; iz<zvertices.size();++iz ){
2232  if (nzvtxes[iz]==0) continue;
2233  zmeans[iz] = zsums[iz]/nzvtxes[iz];
2234  cout<<"zmean for vertex "<< iz <<" = "<<zmeans[iz]<<endl;
2235  zvertices[iz] = zmeans[iz];
2236  zsums[iz]=0.;
2237  if (fill_multi_zvtx){
2238  float ntp_data[2];
2239  ntp_data[0] = _event;
2240  ntp_data[1] = zmeans[iz];
2241  _ntp_zvtx_by_event->Fill(ntp_data);
2242  }
2243  }
2244 
2245  for (unsigned int j = 0; j < nzvtx; ++j){
2246  double zvtx = vtx_tracks[j].z0;
2247  if (zvtx != zvtx) continue;
2248 
2249 #ifdef _MULTIVTX_
2250  if (fill_multi_zvtx){
2251  float ntp_data[2];
2252  ntp_data[0] = _event;
2253  ntp_data[1] = vtx_tracks[j].z0;
2254 // ntp_data[1] = vtx_tracks[j].d;
2255 // cout<<"event "<<ntp_data[0]<<" "<< "vtx "<<ntp_data[1]<<endl;
2256  _ntp_zvtx_by_track->Fill(ntp_data);
2257  }
2258 #endif
2259  bool pass_dcaz= false;
2260  for (unsigned int k = 0; k<zvertices.size(); ++k) {
2261  bool new_vtx = fabs(zvtx-zvertices[k]) < _dcaz_cut;
2262  if (new_vtx) izvtx=k;
2263  pass_dcaz = pass_dcaz || new_vtx;
2264  }
2265 
2266  if (!pass_dcaz || izvtx>99) continue;
2267  zsums[izvtx] += pow(fabs(zmeans[izvtx] - vtx_tracks[j].z0),2);
2268  }
2269 
2270  std::vector<float> _multi_vtx;
2271  std::vector<std::vector<SimpleTrack3D> > _multi_vtx_tracks;
2272  std::vector<std::vector<double> > _multi_vtx_track_errors;
2273  std::vector<std::vector<Eigen::Matrix<float, 5, 5> > > _multi_vtx_track_covars;
2274 
2275 
2276  for (unsigned int i = 0; i < zvertices.size(); ++i)
2277  {
2278  if (nzvtxes[i]==0) continue;
2279  _multi_vtx.push_back(zvertices[i]);
2280  std::vector<SimpleTrack3D> one_vtx_tracks;
2281  _multi_vtx_tracks.push_back(one_vtx_tracks);
2282  std::vector<double> one_track_errors;
2283  _multi_vtx_track_errors.push_back(one_track_errors);
2284  std::vector<Eigen::Matrix<float, 5, 5> > one_track_covars;
2285  _multi_vtx_track_covars.push_back(one_track_covars);
2286  zsigmas[i] = sqrt(zsums[i]/(nzvtxes[i]-1));
2287  cout<<"zsigma for vertex "<<i <<" = "<<zsigmas[i]<<endl;
2288  }
2289 
2290  unsigned int nzvtx_final = 0;
2291  for (unsigned int k = 0; k < nzvtx; ++k){
2292 
2293  float zvtx = vtx_tracks[k].z0;
2294  bool pass_dcaz= false;
2295  for (unsigned int iz = 0; iz<_multi_vtx.size(); ++iz) {
2296  bool new_vtx = fabs(zvtx-_multi_vtx[iz]) < _dcaz_cut;
2297  if (new_vtx) izvtx=iz;
2298  pass_dcaz = pass_dcaz || new_vtx;
2299  }
2300 
2301  if (!pass_dcaz || izvtx >= _multi_vtx.size()) continue;
2302 // cout<<"adding a track to vtx "<<izvtx<<endl;
2303 // if (fabs(vtx_tracks[k].z0-zmeans[k])> _dcaz_cut/*10*zsigmz*/)continue;
2304  ++nzvtx_final;
2305 
2306  _multi_vtx_tracks[izvtx].push_back(vtx_tracks[k]);
2307  _multi_vtx_track_covars[izvtx].push_back(_track_states[k].C);
2308  _multi_vtx_track_errors[izvtx].push_back(_track_states[k].chi2);
2309 
2310  }
2311 
2312  cout<<"start fitting vertex.. "<<endl;
2313  for (unsigned int i = 0; i<_multi_vtx.size(); ++i)
2314  {
2315  if (_multi_vtx_tracks[i].size()==0) continue;
2316  _vertex[2] = _multi_vtx[i];
2317  if (Verbosity() > 0) {
2318  cout << " seed track vertex pre-fit: "
2319  << _vertex[0] << " "
2320  << _vertex[1] << " "
2321  << _vertex[2] << endl;
2322  }
2323 
2324  _vertex_finder.findVertex( _multi_vtx_tracks[i], _multi_vtx_track_covars[i], _vertex, 0.10, true);
2325  _vertex_finder.findVertex( _multi_vtx_tracks[i], _multi_vtx_track_covars[i], _vertex, 0.02, true);
2326  _vertex_finder.findVertex( _multi_vtx_tracks[i], _multi_vtx_track_covars[i], _vertex, 0.005, true);
2327 
2328  n_vtx_tracks = _multi_vtx_tracks[i].size();
2329  cout<<"number of fitted tracks for vertex "<<i<< " : "<<n_vtx_tracks <<endl;
2330 
2331  if (Verbosity() > 0) {
2332  cout << " seed track vertex post-fit: "
2333  << _vertex[0] << " " << _vertex[1] << " " << _vertex[2] << endl;
2334  }
2335 
2336  // add vertex_id to SimpleTrack3D
2337 // std::vector<float> pre_vtx(3,0.0);
2338 
2339  _vertex_list.push_back(_vertex);
2340 
2341  }// loop over vertices
2342 
2343  cout<<"final vertices : "<<endl;
2344  for (unsigned int i = 0; i<_vertex_list.size(); ++i){
2345  cout<<i<<" "<<_vertex_list[i][2]<<endl;
2346  }
2347 
2348 
2349  for (unsigned int k = 0; k < nzvtx; ++k){
2350 
2351  float zvtx = vtx_tracks[k].z0;
2352  bool pass_dcaz= false;
2353  for (unsigned int iz = 0; iz < _vertex_list.size(); ++iz) {
2354  bool new_vtx = fabs(zvtx - _vertex_list[iz][2]) < _dcaz_cut;
2355  if (new_vtx) izvtx=iz;
2356  pass_dcaz = pass_dcaz || new_vtx;
2357  }
2358 
2359 
2360  if (!pass_dcaz || izvtx >= _vertex_list.size()) continue;
2361 // if (pass_dcaz && izvtx <9999)
2362  ++nzvtx_final;
2363 
2364  SimpleTrack3D vtx_track = vtx_tracks[k];
2365  vtx_track.set_vertex_id(izvtx);
2366 // cout<<"vertex_id "<< izvtx<<endl;
2367  _tracks.push_back(vtx_track);
2368  _track_covars.push_back(_track_states[k].C);
2369  _track_errors.push_back(_track_states[k].chi2 );
2370 
2371  izvtx=0;
2372  }
2373 
2374  for(unsigned int i=0; i<vtx_tracks.size(); ++i) vtx_tracks[i].reset();
2375  vtx_tracks.clear();
2376  _multi_vtx_tracks.clear();
2377  _multi_vtx_track_covars.clear();
2378  _multi_vtx_track_errors.clear();
2379 
2380  return Fun4AllReturnCodes::EVENT_OK;
2381 
2382 }
2383 
2384 void PHPatternReco::convertHelixCovarianceToEuclideanCovariance(float B,
2385  float phi, float d, float kappa, float z0, float dzdl,
2386  Eigen::Matrix<float, 5, 5> const& input,
2387  Eigen::Matrix<float, 6, 6>& output) {
2388 
2389  Eigen::Matrix<float, 6, 5> J = Eigen::Matrix<float, 6, 5>::Zero(6, 5);
2390 
2391  // phi,d,nu,z0,dzdl
2392  // -->
2393  // x,y,z,px,py,pz
2394 
2395  float nu = sqrt(kappa);
2396  float dk_dnu = 2 * nu;
2397 
2398  float cosphi = cos(phi);
2399  float sinphi = sin(phi);
2400 
2401  J(0, 0) = -sinphi * d;
2402  J(0, 1) = cosphi;
2403  J(1, 0) = d * cosphi;
2404  J(1, 1) = sinphi;
2405  J(2, 3) = 1;
2406 
2407  float pt = 0.003 * B / kappa;
2408  float dpt_dk = -0.003 * B / (kappa * kappa);
2409 
2410  J(3, 0) = -cosphi * pt;
2411  J(3, 2) = -sinphi * dpt_dk * dk_dnu;
2412  J(4, 0) = -sinphi * pt;
2413  J(4, 2) = cosphi * dpt_dk * dk_dnu;
2414 
2415  float alpha = 1. / (1. - dzdl * dzdl);
2416  float alpha_half = pow(alpha, 0.5);
2417  float alpha_3_half = alpha * alpha_half;
2418 
2419  J(5, 2) = dpt_dk * dzdl * alpha_half * dk_dnu;
2420  J(5, 4) = pt * (alpha_half + dzdl * dzdl * alpha_3_half) * dk_dnu;
2421 
2422  output = J * input * (J.transpose());
2423 }
2424 
2425 void PHPatternReco::shift_coordinate_system(double dx, double dy, double dz) {
2426 
2427  for (std::map<unsigned int,SimpleHit3D>::iterator it=hits_map.begin();
2428  it!=hits_map.end();
2429  ++it)
2430  {
2431  SimpleHit3D hit = it->second;
2432  hit.set_x(hit.get_x() + dx);
2433  hit.set_y(hit.get_y() + dy);
2434  hit.set_z(hit.get_z() + dz);
2435  it->second = hit;
2436 
2437  }
2438 
2439  for (unsigned int itrack = 0; itrack < _tracks.size(); ++itrack) {
2440  for (unsigned int ihit = 0; ihit < _tracks[itrack].hits.size(); ++ihit) {
2441  _tracks[itrack].hits[ihit].set_x(_tracks[itrack].hits[ihit].get_x() + dx);
2442  _tracks[itrack].hits[ihit].set_y(_tracks[itrack].hits[ihit].get_y() + dy);
2443  _tracks[itrack].hits[ihit].set_z(_tracks[itrack].hits[ihit].get_z() + dz);
2444  }
2445  }
2446 
2447  _vertex[0] += dx;
2448  _vertex[1] += dy;
2449  _vertex[2] += dz;
2450 
2451  return;
2452 }
2453 
2454 float PHPatternReco::shift_phi_range(float _phi){
2455 
2456  if (_phi < 0.) _phi += 2.*M_PI;
2457  return _phi;
2458 }
2459