ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
PHSimpleVertexFinder.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file PHSimpleVertexFinder.cc
1 #include "PHSimpleVertexFinder.h"
2 
4 
5 #include <trackbase/TrkrCluster.h> // for TrkrCluster
6 #include <trackbase/TrkrDefs.h> // for cluskey, getLayer, TrkrId
9 #include <trackbase_historic/SvtxTrack.h> // for SvtxTrack, SvtxTrack::C...
14 
15 #include <phool/PHCompositeNode.h>
16 #include <phool/PHNode.h>
17 #include <phool/PHNodeIterator.h>
19 
20 #include <phool/getClass.h>
21 #include <phool/phool.h>
22 
23 #include <cmath> // for sqrt, fabs, atan2, cos
24 #include <iostream> // for operator<<, basic_ostream
25 #include <map> // for map
26 #include <set> // for _Rb_tree_const_iterator
27 #include <utility> // for pair, make_pair
28 #include <iomanip>
29 
30 #include <vector>
31 #include <cassert>
32 #include <numeric>
33 #include <iostream>
34 #include <algorithm>
35 #include <functional>
36 
37 #include <Eigen/Dense>
38 
39 //____________________________________________________________________________..
41  : SubsysReco(name)
42 {
43 
44 }
45 
46 //____________________________________________________________________________..
48 {
49 
50 }
51 
52 //____________________________________________________________________________..
54 {
55  int ret = GetNodes(topNode);
56  if (ret != Fun4AllReturnCodes::EVENT_OK) return ret;
57  ret = CreateNodes(topNode);
58  if (ret != Fun4AllReturnCodes::EVENT_OK) return ret;
59  return ret;
60 }
61 
62 //____________________________________________________________________________..
64 {
65 
66  if(Verbosity() > 0)
67  std::cout << PHWHERE << " track map size " << _track_map->size() << std::endl;
68 
69  // reset maps
70  _vertex_track_map.clear();
71  _track_pair_map.clear();
72  _track_pair_pca_map.clear();
73  _vertex_position_map.clear();
74  _vertex_covariance_map.clear();
75  _vertex_set.clear();
76 
77  // Find all instances where two tracks have a dca of < _dcacut, and capture the pair details
78  // Fills _track_pair_map and _track_pair_pca_map
79  checkDCAs();
80 
82  if(_track_pair_map.size() == 0)
83  {
84  _dcacut *= 1.5;
85  checkDCAs();
86  }
87 
88  // get all connected pairs of tracks by looping over the track_pair map
89  std::vector<std::set<unsigned int>> connected_tracks = findConnectedTracks();
90 
91  // make vertices - each set of connected tracks is a vertex
92  for(unsigned int ivtx = 0; ivtx < connected_tracks.size(); ++ivtx)
93  {
94  if(Verbosity() > 1) std::cout << "process vertex " << ivtx << std::endl;
95 
96  for(auto it : connected_tracks[ivtx])
97  {
98  unsigned int id = it;
99  _vertex_track_map.insert(std::make_pair(ivtx, id));
100  if(Verbosity() > 2) std::cout << " adding track " << id << " to vertex " << ivtx << std::endl;
101  }
102  }
103 
104  // make a list of vertices
105  for(auto it : _vertex_track_map)
106  {
107  if(Verbosity() > 1) std::cout << " vertex " << it.first << " track " << it.second << std::endl;
108  _vertex_set.insert(it.first);
109  }
110 
111  // this finds average vertex positions after removal of outlying track pairs
113 
114  // average covariance for accepted tracks
115  for(auto it : _vertex_set)
116  {
117  matrix_t avgCov = matrix_t::Zero();
118  double cov_wt = 0.0;
119 
120  auto ret = _vertex_track_map.equal_range(it);
121  for (auto cit=ret.first; cit!=ret.second; ++cit)
122  {
123  unsigned int trid = cit->second;
124  matrix_t cov;
125  auto track = _track_map->get(trid);
126  for(int i = 0; i < 3; ++i)
127  for(int j = 0; j < 3; ++j)
128  cov(i,j) = track->get_error(i,j);
129 
130  avgCov += cov;
131  cov_wt++;
132  }
133 
134  avgCov /= sqrt(cov_wt);
135  if(Verbosity() > 2)
136  {
137  std::cout << "Average covariance for vertex " << it << " is:" << std::endl;
138  std::cout << std::setprecision(8) << avgCov << std::endl;
139  }
140  _vertex_covariance_map.insert(std::make_pair(it, avgCov));
141  }
142 
143  // Write the vertices to the vertex map on the node tree
144  //==============================================
145  if(_vertex_track_map.size() > 0)
147 
148  for(auto it : _vertex_set)
149  {
150  auto svtxVertex = std::make_unique<SvtxVertex_v1>();
151 
152  svtxVertex->set_chisq(0.0);
153  svtxVertex->set_ndof(0);
154  svtxVertex->set_t0(0);
155  svtxVertex->set_id(it);
156 
157  auto ret = _vertex_track_map.equal_range(it);
158  for (auto cit=ret.first; cit!=ret.second; ++cit)
159  {
160  unsigned int trid = cit->second;
161 
162  if(Verbosity() > 1) std::cout << " vertex " << it << " insert track " << trid << std::endl;
163  svtxVertex->insert_track(trid);
164  _track_map->get(trid)->set_vertex_id(it);
165  }
166 
167  Eigen::Vector3d pos = _vertex_position_map.find(it)->second;
168  svtxVertex->set_x(pos.x());
169  svtxVertex->set_y(pos.y());
170  svtxVertex->set_z(pos.z());
171  if(Verbosity() > 1) std::cout << " vertex " << it << " insert pos.x " << pos.x() << " pos.y " << pos.y() << " pos.z " << pos.z() << std::endl;
172 
173  auto vtxCov = _vertex_covariance_map.find(it)->second;
174  for(int i = 0; i < 3; ++i)
175  {
176  for(int j = 0; j < 3; ++j)
177  {
178  svtxVertex->set_error(i, j, vtxCov(i,j));
179  }
180  }
181 
182  _svtx_vertex_map->insert(svtxVertex.release());
183  }
184 
189  //=================================================
190  for(const auto& [trackkey, track] : *_track_map)
191  {
192  auto vtxid = track->get_vertex_id();
193 
195  if(_svtx_vertex_map->get(vtxid))
196  { continue; }
197 
198  float maxdz = std::numeric_limits<float>::max();
199  unsigned int newvtxid = std::numeric_limits<unsigned int>::max();
200 
201  for(const auto& [vtxkey, vertex] : *_svtx_vertex_map)
202  {
203  float dz = track->get_z() - vertex->get_z();
204  if(fabs(dz) < maxdz)
205  {
206  maxdz = dz;
207  newvtxid = vtxkey;
208  }
209  }
210 
211  track->set_vertex_id(newvtxid);
212  }
213 
215 }
216 
217 
218 
220 {
222 }
223 
225 {
226  PHNodeIterator iter(topNode);
227 
228  PHCompositeNode *dstNode = dynamic_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", "DST"));
229 
231  if (!dstNode)
232  {
233  std::cerr << "DST Node missing, quitting" << std::endl;
234  throw std::runtime_error("failed to find DST node in PHActsInitialVertexFinder::createNodes");
235  }
236 
238  PHCompositeNode *svtxNode = dynamic_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", "SVTX"));
239 
240  if (!svtxNode)
241  {
242  svtxNode = new PHCompositeNode("SVTX");
243  dstNode->addNode(svtxNode);
244  }
245 
246  _svtx_vertex_map = findNode::getClass<SvtxVertexMap>(topNode,"SvtxVertexMap");
247 
248  if(!_svtx_vertex_map)
249  {
252  _svtx_vertex_map, "SvtxVertexMap","PHObject");
253 
254  svtxNode->addNode(vertexNode);
255 
256  }
257 
259 }
261 {
262 
263  _track_map = findNode::getClass<SvtxTrackMap>(topNode, "SvtxTrackMap");
264  if (!_track_map)
265  {
266  std::cout << PHWHERE << " ERROR: Can't find SvtxTrackMap: " << std::endl;
268  }
269 
270 
271 
273 }
274 
276 {
277  // Loop over tracks and check for close DCA match with all other tracks
278  for(auto tr1_it = _track_map->begin(); tr1_it != _track_map->end(); ++tr1_it)
279  {
280  auto id1 = tr1_it->first;
281  auto tr1 = tr1_it->second;
282  if(tr1->get_quality() > _qual_cut) continue;
283  if(_require_mvtx)
284  {
285  unsigned int nmvtx = 0;
286  for(auto clusit = tr1->begin_cluster_keys(); clusit != tr1->end_cluster_keys(); ++clusit)
287  {
288  if(TrkrDefs::getTrkrId(*clusit) == TrkrDefs::mvtxId )
289  {
290  nmvtx++;
291  }
292  if(nmvtx >= _nmvtx_required) break;
293  }
294  if(nmvtx < _nmvtx_required) continue;
295  if(Verbosity() > 3) std::cout << " tr1 has nmvtx at least " << nmvtx << std::endl;
296  }
297 
298  // look for close DCA matches with all other such tracks
299  for(auto tr2_it = std::next(tr1_it); tr2_it != _track_map->end(); ++tr2_it)
300  {
301  auto id2 = tr2_it->first;
302  auto tr2 = tr2_it->second;
303  if(tr2->get_quality() > _qual_cut) continue;
304  if(_require_mvtx)
305  {
306  unsigned int nmvtx = 0;
307  for(auto clusit = tr2->begin_cluster_keys(); clusit != tr2->end_cluster_keys(); ++clusit)
308  {
309  if(TrkrDefs::getTrkrId(*clusit) == TrkrDefs::mvtxId)
310  {
311  nmvtx++;
312  }
313  if(nmvtx >= _nmvtx_required) break;
314  }
315  if(nmvtx < _nmvtx_required) continue;
316  if(Verbosity() > 3) std::cout << " tr2 has nmvtx at least " << nmvtx << std::endl;
317  }
318 
319  // find DCA of these two tracks
320  if(Verbosity() > 3) std::cout << "Check DCA for tracks " << id1 << " and " << id2 << std::endl;
321 
322  findDcaTwoTracks(tr1, tr2);
323 
324  }
325  }
326 
327 
328 }
329 
331 {
332  if(tr1->get_pt() < _track_pt_cut) return;
333  if(tr2->get_pt() < _track_pt_cut) return;
334 
335  unsigned int id1 = tr1->get_id();
336  unsigned int id2 = tr2->get_id();
337 
338  // get the line equations for the tracks
339 
340  Eigen::Vector3d a1(tr1->get_x(), tr1->get_y(), tr1->get_z());
341  Eigen::Vector3d b1(tr1->get_px() / tr1->get_p(), tr1->get_py() / tr1->get_p(), tr1->get_pz() / tr1->get_p());
342  Eigen::Vector3d a2(tr2->get_x(), tr2->get_y(), tr2->get_z());
343  Eigen::Vector3d b2(tr2->get_px() / tr2->get_p(), tr2->get_py() / tr2->get_p(), tr2->get_pz() / tr2->get_p());
344 
345  Eigen::Vector3d PCA1(0,0,0);
346  Eigen::Vector3d PCA2(0,0,0);
347  double dca = dcaTwoLines(a1, b1, a2, b2, PCA1, PCA2);
348 
349  // check dca cut is satisfied, and that PCA is close to beam line
350 if( fabs(dca) < _dcacut && (fabs(PCA1.x()) < _beamline_xy_cut && fabs(PCA1.y()) < _beamline_xy_cut) )
351  {
352  if(Verbosity() > 3)
353  {
354  std::cout << " good match for tracks " << tr1->get_id() << " and " << tr2->get_id() << " with pT " << tr1->get_pt() << " and " << tr2->get_pt() << std::endl;
355  std::cout << " a1.x " << a1.x() << " a1.y " << a1.y() << " a1.z " << a1.z() << std::endl;
356  std::cout << " a2.x " << a2.x() << " a2.y " << a2.y() << " a2.z " << a2.z() << std::endl;
357  std::cout << " PCA1.x() " << PCA1.x() << " PCA1.y " << PCA1.y() << " PCA1.z " << PCA1.z() << std::endl;
358  std::cout << " PCA2.x() " << PCA2.x() << " PCA2.y " << PCA2.y() << " PCA2.z " << PCA2.z() << std::endl;
359  std::cout << " dca " << dca << std::endl;
360  }
361 
362  // capture the results for successful matches
363  _track_pair_map.insert(std::make_pair(id1,std::make_pair(id2, dca)));
364  _track_pair_pca_map.insert( std::make_pair(id1, std::make_pair(id2, std::make_pair(PCA1, PCA2))) );
365  }
366 
367  return;
368 }
369 
370 double PHSimpleVertexFinder::dcaTwoLines(const Eigen::Vector3d &a1,const Eigen::Vector3d &b1,
371  const Eigen::Vector3d &a2,const Eigen::Vector3d &b2,
372  Eigen::Vector3d &PCA1, Eigen::Vector3d &PCA2)
373 {
374  // The shortest distance between two skew lines described by
375  // a1 + c * b1
376  // a2 + d * b2
377  // where a1, a2, are vectors representing points on the lines, b1, b2 are direction vectors, and c and d are scalars
378  // is:
379  // dca = (b1 x b2) .(a2-a1) / |b1 x b2|
380 
381  // bcrossb/mag_bcrossb is a unit vector perpendicular to both direction vectors b1 and b2
382  auto bcrossb = b1.cross(b2);
383  auto mag_bcrossb = bcrossb.norm();
384  // a2-a1 is the vector joining any arbitrary points on the two lines
385  auto aminusa = a2-a1;
386 
387  // The DCA of these two lines is the projection of a2-a1 along the direction of the perpendicular to both
388  // remember that a2-a1 is longer than (or equal to) the dca by definition
389  double dca = 999;
390  if( mag_bcrossb != 0)
391  dca = bcrossb.dot(aminusa) / mag_bcrossb;
392  else
393  return dca; // same track, skip combination
394 
395  // get the points at which the normal to the lines intersect the lines, where the lines are perpendicular
396  // Assume the shortest distance occurs at points A on line 1, and B on line 2, call the line AB
397  // AB = a2+d*b2 - (a1+c*b1)
398  // we need to find c and d where AB is perpendicular to both lines. so AB.b1 = 0 and AB.b2 = 0
399  // ( (a2 -a1) + d*b2 -c*b1 ).b1 = 0
400  // ( (a2 -a1) + d*b2 -c*b1 ).b2 = 0
401  // so we have two simultaneous equations in 2 unknowns
402  // (a2-a1).b1 + d*b2.b1 -c*b1.b1 = 0 => d*b2.b1 = c*b1.b1 - (a2-a1).b1 => d = (1/b2.b1) * (c*b1.b1 - (a2-a1).b1)
403  // (a2-a1).b2 + d*b2.b2 - c*b1.b2 = 0 => c*b1.b2 = (a2-a1).b2 + [(1/b2.b1) * (c*b1*b1 -(a2-a1).b1)}*b2.b2
404  // c*b1.b2 - (1/b2.b1) * c*b1.b1*b2.b2 = (a2-a1).b2 - (1/b2.b1)*(a2-a1).b1*b2.b2
405  // c *(b1.b2 - (1/b2.b1)*b1.b1*b2.b2) = (a2-a1).b2 - (1/b2.b1)*(a2-a1).b1*b2.b2
406  // call this: c*X = Y
407  // plug back into the d equation
408  // d = c*b1.b1 / b2.b1 - (a2-a1).b1 / b2.b1
409  // and call the d equation: d = c * F - G
410 
411  double X = b1.dot(b2) - b1.dot(b1) * b2.dot(b2) / b2.dot(b1);
412  double Y = (a2.dot(b2) - a1.dot(b2)) - (a2.dot(b1) - a1.dot(b1)) * b2.dot(b2) / b2.dot(b1) ;
413  double c = Y/X;
414 
415  double F = b1.dot(b1) / b2.dot(b1);
416  double G = - (a2.dot(b1) - a1.dot(b1)) / b2.dot(b1);
417  double d = c * F + G;
418 
419  // then the points of closest approach are:
420  PCA1 = a1+c*b1;
421  PCA2 = a2+d*b2;
422 
423  return dca;
424 
425 
426 }
427 
428 std::vector<std::set<unsigned int>> PHSimpleVertexFinder::findConnectedTracks()
429 {
430  std::vector<std::set<unsigned int>> connected_tracks;
431  std::set<unsigned int> connected;
432  std::set<unsigned int> used;
433  for(auto it : _track_pair_map)
434  {
435  unsigned int id1 = it.first;
436  unsigned int id2 = it.second.first;
437 
438  if( (used.find(id1) != used.end()) && (used.find(id2) != used.end()) )
439  {
440  if(Verbosity() > 3) std::cout << " tracks " << id1 << " and " << id2 << " are both in used , skip them" << std::endl;
441  continue;
442  }
443  else if( (used.find(id1) == used.end()) && (used.find(id2) == used.end()))
444  {
445  if(Verbosity() > 3) std::cout << " tracks " << id1 << " and " << id2 << " are both not in used , start a new connected set" << std::endl;
446  // close out and start a new connections set
447  if(connected.size() > 0)
448  {
449  connected_tracks.push_back(connected);
450  connected.clear();
451  if(Verbosity() > 3) std::cout << " closing out set " << std::endl;
452  }
453  }
454 
455  // get everything connected to id1 and id2
456  connected.insert(id1);
457  used.insert(id1);
458  connected.insert(id2);
459  used.insert(id2);
460  for(auto cit : _track_pair_map)
461  {
462  unsigned int id3 = cit.first;
463  unsigned int id4 = cit.second.first;
464  if( (connected.find(id3) != connected.end()) || (connected.find(id4) != connected.end()) )
465  {
466  if(Verbosity() > 3) std::cout << " found connection to " << id3 << " and " << id4 << std::endl;
467  connected.insert(id3);
468  used.insert(id3);
469  connected.insert(id4);
470  used.insert(id4);
471  }
472  }
473  }
474 
475  // close out the last set
476  if(connected.size() > 0)
477  {
478  connected_tracks.push_back(connected);
479  connected.clear();
480  if(Verbosity() > 3) std::cout << " closing out last set " << std::endl;
481  }
482 
483  if(Verbosity() > 3) std::cout << "connected_tracks size " << connected_tracks.size() << std::endl;
484 
485  return connected_tracks;
486 }
487 
489 {
490  // Note: std::multimap<unsigned int, std::pair<unsigned int, std::pair<Eigen::Vector3d, Eigen::Vector3d>>> _track_pair_pca_map
491 
492  for(auto it : _vertex_set)
493  {
494  unsigned int vtxid = it;
495  if(Verbosity() > 1) std::cout << "calculate average position for vertex " << vtxid << std::endl;
496 
497  // we need the median values of the x and y positions
498  std::vector<double> vx;
499  std::vector<double> vy;
500  std::vector<double> vz;
501 
502  double pca_median_x = 0.;
503  double pca_median_y = 0.;
504  double pca_median_z = 0.;
505 
506  Eigen::Vector3d new_pca_avge(0.,0.,0.);
507  double new_wt = 0.0;
508 
509  auto ret = _vertex_track_map.equal_range(vtxid);
510 
511  // Start by getting the positions for this vertex into vectors for the median calculation
512  for (auto cit=ret.first; cit!=ret.second; ++cit)
513  {
514  unsigned int tr1id = cit->second;
515  if(Verbosity() > 2) std::cout << " vectors: get entries for track " << tr1id << " for vertex " << vtxid << std::endl;
516 
517  // find all pairs for this vertex with tr1id
518  auto pca_range = _track_pair_pca_map.equal_range(tr1id);
519  for (auto pit=pca_range.first; pit!=pca_range.second; ++pit)
520  {
521  unsigned int tr2id = pit->second.first;
522 
523  Eigen::Vector3d PCA1 = pit->second.second.first;
524  Eigen::Vector3d PCA2 = pit->second.second.second;
525 
526  if(Verbosity() > 2)
527  std::cout << " vectors: tr1id " << tr1id << " tr2id " << tr2id
528  << " PCA1 " << PCA1.x() << " " << PCA1.y() << " " << PCA1.z()
529  << " PCA2 " << PCA2.x() << " " << PCA2.y() << " " << PCA2.z()
530  << std::endl;
531 
532  vx.push_back(PCA1.x());
533  vx.push_back(PCA2.x());
534  vy.push_back(PCA1.y());
535  vy.push_back(PCA2.y());
536  vz.push_back(PCA1.z());
537  vz.push_back(PCA2.z());
538  }
539  }
540 
541  // Get the medians for this vertex
542  // Using the median as a reference for rejecting outliers only makes sense for more than 2 tracks
543  if(vx.size() < 3)
544  {
545  new_pca_avge.x() = getAverage(vx);
546  new_pca_avge.y() = getAverage(vy);
547  new_pca_avge.z() = getAverage(vz);
548  _vertex_position_map.insert(std::make_pair(vtxid, new_pca_avge));
549  if(Verbosity() > 1)
550  std::cout << " Vertex has only 2 tracks, use average for PCA: " << new_pca_avge.x() << " " << new_pca_avge.y() << " " << new_pca_avge.z() << std::endl;
551 
552  // done with this vertex
553  continue;
554  }
555 
556  pca_median_x = getMedian(vx);
557  pca_median_y = getMedian(vy);
558  pca_median_z = getMedian(vz);
559  if(Verbosity() > 1) std::cout << "Median values: x " << pca_median_x << " y " << pca_median_y << std::endl;
560 
561  // Make the average vertex position with outlier rejection wrt the median
562  for (auto cit=ret.first; cit!=ret.second; ++cit)
563  {
564  unsigned int tr1id = cit->second;
565  if(Verbosity() > 2) std::cout << " average: get entries for track " << tr1id << " for vertex " << vtxid << std::endl;
566 
567  // find all pairs for this vertex with tr1id
568  auto pca_range = _track_pair_pca_map.equal_range(tr1id);
569  for (auto pit=pca_range.first; pit!=pca_range.second; ++pit)
570  {
571  unsigned int tr2id = pit->second.first;
572 
573  Eigen::Vector3d PCA1 = pit->second.second.first;
574  Eigen::Vector3d PCA2 = pit->second.second.second;
575 
576  if(
577  fabs(PCA1.x() - pca_median_x) < _outlier_cut &&
578  fabs(PCA1.y() - pca_median_y) < _outlier_cut &&
579  fabs(PCA2.x() - pca_median_x) < _outlier_cut &&
580  fabs(PCA2.y() - pca_median_y) < _outlier_cut
581  )
582  {
583  // good track pair, add to new average
584 
585  new_pca_avge += PCA1;
586  new_wt++;
587  new_pca_avge += PCA2;
588  new_wt++;
589  }
590  else
591  {
592  if(Verbosity() > 1) std::cout << "Reject pair with tr1id " << tr1id << " tr2id " << tr2id << std::endl;
593  }
594  }
595  }
596  if(new_wt > 0.0)
597  new_pca_avge = new_pca_avge / new_wt;
598  else
599  {
600  // There were no pairs that survived the track cuts, use the median values
601  new_pca_avge.x() = pca_median_x;
602  new_pca_avge.y() = pca_median_y;
603  new_pca_avge.z() = pca_median_z;
604  }
605 
606  _vertex_position_map.insert(std::make_pair(vtxid, new_pca_avge));
607 
608  }
609 
610  return;
611  }
612 
613 double PHSimpleVertexFinder::getMedian(std::vector<double> &v)
614 {
615  double median = 0.0;
616 
617  if( (v.size() % 2) == 0)
618  {
619  // even number of entries
620  // we want the average of the middle two numbers, v.size()/2 and v.size()/2-1
621  auto m1 = v.begin() + v.size()/2;
622  std::nth_element(v.begin(), m1, v.end());
623  double median1 = v[v.size()/2];
624 
625  auto m2 = v.begin() + v.size()/2 - 1;
626  std::nth_element(v.begin(), m2, v.end());
627  double median2 = v[v.size()/2 - 1];
628 
629  median = (median1 + median2) / 2.0;
630  if(Verbosity() > 2) std::cout << "The vector size is " << v.size()
631  << " element m is " << v.size() / 2 << " = " << v[v.size()/2]
632  << " element m-1 is " << v.size() / 2 -1 << " = " << v[v.size()/2-1]
633  << std::endl;
634  }
635  else
636  {
637  // odd number of entries
638  auto m = v.begin() + v.size()/2;
639  std::nth_element(v.begin(), m, v.end());
640  median = v[v.size()/2];
641  if(Verbosity() > 2) std::cout << "The vector size is " << v.size() << " element m is " << v.size() / 2 << " = " << v[v.size()/2] << std::endl;
642  }
643 
644  return median ;
645 }
646 double PHSimpleVertexFinder::getAverage(std::vector<double> &v)
647 {
648  double avge = 0.0;
649  double wt = 0.0;
650  for(auto it : v)
651  {
652  avge += it;
653  wt++;
654  }
655 
656  avge /= wt;
657  if(Verbosity() > 2)
658  {
659  std::cout << " average = " << avge << std::endl;
660  }
661 
662  return avge;
663 }
664 
665