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