ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
PHRTreeSeeding.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file PHRTreeSeeding.cc
1 
8 //begin
9 
10 #include "PHRTreeSeeding.h"
11 
12 // trackbase_historic includes
17 
18 #include <trackbase/TrkrCluster.h> // for TrkrCluster
20 #include <trackbase/TrkrDefs.h> // for getLayer, clu...
21 #include <trackbase/TrkrHitSet.h>
23 
24 // sPHENIX Geant4 includes
29 
30 // sPHENIX includes
32 
33 #include <phool/PHTimer.h> // for PHTimer
34 #include <phool/getClass.h>
35 #include <phool/phool.h> // for PHWHERE
36 
37 //ROOT includes for debugging
38 #include <TFile.h>
39 #include <TNtuple.h>
40 #include <TVector3.h> // for TVector3
41 
42 //BOOST for combi seeding
43 #include <boost/geometry.hpp>
44 #include <boost/geometry/geometries/box.hpp>
45 #include <boost/geometry/geometries/point.hpp>
46 #include <boost/geometry/index/rtree.hpp>
47 
48 #include <algorithm>
49 #include <cmath>
50 #include <iostream>
51 #include <numeric>
52 #include <utility> // for pair, make_pair
53 #include <vector>
54 
55 // forward declarations
57 
58 //end
59 
60 typedef bg::model::point<float, 3, bg::cs::cartesian> point;
61 typedef bg::model::box<point> box;
62 typedef std::pair<point, TrkrDefs::cluskey> pointKey;
63 
64 #define LogDebug(exp) std::cout << "DEBUG: " << __FILE__ << ": " << __LINE__ << ": " << exp
65 #define LogError(exp) std::cout << "ERROR: " << __FILE__ << ": " << __LINE__ << ": " << exp
66 #define LogWarning(exp) std::cout << "WARNING: " << __FILE__ << ": " << __LINE__ << ": " << exp
67 
68 using namespace std;
69 //using namespace ROOT::Minuit2;
70 namespace bg = boost::geometry;
71 namespace bgi = boost::geometry::index;
72 
73 vector<TrkrCluster *> clusterpoints;
76  const string &name,
77  unsigned int nlayers_maps,
78  unsigned int nlayers_intt,
79  unsigned int nlayers_tpc,
80  unsigned int start_layer)
81  : PHTrackSeeding(name)
82  , _g4clusters(nullptr)
83  , _g4tracks(nullptr)
84  , _g4vertexes(nullptr)
85  , _cluster_map(nullptr)
86  , _svtxhitsmap(nullptr)
87  , _hit_used_map(nullptr)
88  , _hit_used_map_size(0)
89  , phisr(0.005)
90  , etasr(0.0035)
91  , phist(0.001)
92  , etast(0.003)
93  , phixt(0.008)
94  , etaxt(0.005)
95  , _nlayers_maps(nlayers_maps)
96  , _nlayers_intt(nlayers_intt)
97  , _nlayers_tpc(nlayers_tpc)
98  , _start_layer(start_layer)
99  , _phi_scale(2)
100  , _z_scale(2)
101 {
102 }
103 
105 {
107  PHG4CylinderCellGeomContainer>(topNode, "CYLINDERCELLGEOM_SVTX");
109  PHG4CylinderGeomContainer>(topNode, "CYLINDERGEOM_INTT");
111  PHG4CylinderGeomContainer>(topNode, "CYLINDERGEOM_MVTX");
112 
113  //_nlayers_seeding = _seeding_layer.size();
114  //_radii.assign(_nlayers_seeding, 0.0);
115  map<float, int> radius_layer_map;
116 
117  _radii_all.assign(60, 0.0);
118  _layer_ilayer_map.clear();
119  _layer_ilayer_map_all.clear();
120  if (cellgeos)
121  {
123  cellgeos->get_begin_end();
125  layerrange.first;
126  layeriter != layerrange.second; ++layeriter)
127  {
128  radius_layer_map.insert(
129  make_pair(layeriter->second->get_radius(),
130  layeriter->second->get_layer()));
131  }
132  }
133 
134  if (laddergeos)
135  {
137  laddergeos->get_begin_end();
139  layerrange.first;
140  layeriter != layerrange.second; ++layeriter)
141  {
142  radius_layer_map.insert(
143  make_pair(layeriter->second->get_radius(),
144  layeriter->second->get_layer()));
145  }
146  }
147 
148  if (mapsladdergeos)
149  {
151  mapsladdergeos->get_begin_end();
153  layerrange.first;
154  layeriter != layerrange.second; ++layeriter)
155  {
156  radius_layer_map.insert(
157  make_pair(layeriter->second->get_radius(),
158  layeriter->second->get_layer()));
159  }
160  }
161  for (map<float, int>::iterator iter = radius_layer_map.begin();
162  iter != radius_layer_map.end(); ++iter)
163  {
164  _layer_ilayer_map_all.insert(make_pair(iter->second, _layer_ilayer_map_all.size()));
165 
166  /*if (std::find(_seeding_layer.begin(), _seeding_layer.end(),
167  iter->second) != _seeding_layer.end())
168  {
169  _layer_ilayer_map.insert(make_pair(iter->second, ilayer));
170  ++ilayer;
171  }*/
172  }
173  if (cellgeos)
174  {
176  cellgeos->get_begin_end();
177  PHG4CylinderCellGeomContainer::ConstIterator miter = begin_end.first;
178  for (; miter != begin_end.second; miter++)
179  {
180  PHG4CylinderCellGeom *geo = miter->second;
182  geo->get_radius() + 0.5 * geo->get_thickness();
183 
184  /*if (_layer_ilayer_map.find(geo->get_layer()) != _layer_ilayer_map.end())
185  {
186  _radii[_layer_ilayer_map[geo->get_layer()]] =
187  geo->get_radius();
188  }*/
189  }
190  }
191 
192  if (laddergeos)
193  {
195  laddergeos->get_begin_end();
196  PHG4CylinderGeomContainer::ConstIterator miter = begin_end.first;
197  for (; miter != begin_end.second; miter++)
198  {
199  PHG4CylinderGeom *geo = miter->second;
201  geo->get_radius() + 0.5 * geo->get_thickness();
202 
203  /*if (_layer_ilayer_map.find(geo->get_layer()) != _layer_ilayer_map.end())
204  {
205  _radii[_layer_ilayer_map[geo->get_layer()]] = geo->get_radius();
206  }*/
207  }
208  }
209 
210  if (mapsladdergeos)
211  {
213  mapsladdergeos->get_begin_end();
214  PHG4CylinderGeomContainer::ConstIterator miter = begin_end.first;
215  for (; miter != begin_end.second; miter++)
216  {
217  PHG4CylinderGeom *geo = miter->second;
218 
219  //if(geo->get_layer() > (int) _radii.size() ) continue;
220 
221  // if (Verbosity() >= 2)
222  // geo->identify();
223 
224  //TODO
226  geo->get_radius();
227 
228  /*if (_layer_ilayer_map.find(geo->get_layer()) != _layer_ilayer_map.end())
229  {
230  _radii[_layer_ilayer_map[geo->get_layer()]] = geo->get_radius();
231  }*/
232  }
233  }
235 }
236 
238 {
239  //TrkrCluster tca = *tc1;
240  //TrkrCluster tcb = *tc2;
241  TrkrDefs::cluskey ck1 = tc1->getClusKey();
242  TrkrDefs::cluskey ck2 = tc2->getClusKey();
243  if (TrkrDefs::getLayer(ck1) != TrkrDefs::getLayer(ck2)) return TrkrDefs::getLayer(ck1) < TrkrDefs::getLayer(ck2);
244  return ck1 < ck2;
245  //return ck1/16<ck2/16;
246 }
247 
248 bool largeintersection(const vector<TrkrCluster *> &v1,
249  const vector<TrkrCluster *> &v2)
250 {
251  std::vector<TrkrCluster *> v3;
252  vector<TrkrCluster *> va(v1);
253  vector<TrkrCluster *> vb(v2);
254  //std::sort(v1.begin(), v1.end());
255  //std::sort(v2.begin(), v2.end());
256  std::set_intersection(va.begin(), va.end(),
257  vb.begin(), vb.end(),
258  back_inserter(v3), comparing);
259  //cout << "Got an intersection of size " << v3.size() << endl;
260  if (v3.size() < 3) return false;
261  //TrkrCluster *temp1 = new TrkrCluster;
262  //*temp1 = *va[0];
263 
264  /*TrkrCluster *temp2 = vb[0];
265  unsigned int ii = va.size();
266  unsigned int jj = vb.size();
267  TrkrCluster *temp3 = va[ii-1];
268  TrkrCluster *temp4 = vb[jj-1];
269  TrkrDefs::cluskey ck1 = (va[0])->getClusKey();
270  TrkrDefs::cluskey ck2 = temp2->getClusKey();
271  TrkrDefs::cluskey ck3 = temp3->getClusKey();
272  TrkrDefs::cluskey ck4 = temp4->getClusKey();*/
273  /*bool to_ret = (TrkrDefs::getLayer(ck1)-TrkrDefs::getLayer(ck2))*(TrkrDefs::getLayer(ck3)-TrkrDefs::getLayer(ck4))>0;
274  cout << to_ret << endl;
275  return to_ret;*/
276  return (TrkrDefs::getLayer(va[0]->getClusKey()) - TrkrDefs::getLayer(vb[0]->getClusKey())) * (TrkrDefs::getLayer(va[va.size() - 1]->getClusKey()) - TrkrDefs::getLayer(vb[vb.size() - 1]->getClusKey())) > 0 && max(TrkrDefs::getLayer(va[va.size() - 1]->getClusKey()) - TrkrDefs::getLayer(vb[0]->getClusKey()), TrkrDefs::getLayer(vb[vb.size() - 1]->getClusKey()) - TrkrDefs::getLayer(va[0]->getClusKey())) + 1 == (int) (va.size() + vb.size() - v3.size());
277 }
278 
279 vector<TrkrCluster *> yunion(const vector<TrkrCluster *> &v1, const vector<TrkrCluster *> &v2)
280 {
281  std::vector<TrkrCluster *> v3;
282 
283  //std::sort(v1.begin(), v1.end());
284  //std::sort(v2.begin(), v2.end());
285 
286  std::set_union(v1.begin(), v1.end(),
287  v2.begin(), v2.end(),
288  back_inserter(v3), comparing);
289  return v3;
290 }
291 
292 bool issuperset(const vector<TrkrCluster *> &larger,
293  const vector<TrkrCluster *> &smaller)
294 {
295  for (unsigned int i = 0; i < smaller.size(); i++)
296  {
297  if (!binary_search(larger.begin(), larger.end(), smaller.at(i), comparing)) return false;
298  }
299  return true;
300 }
301 
302 void changetounion(vector<TrkrCluster *> &v1,
303  const vector<TrkrCluster *> &v2)
304 {
305  vector<TrkrCluster *> temp(yunion(v1, v2));
306  v1 = temp;
307 }
308 
309 double RosenBrock(const double *xx)
310 {
311  /* const Double_t x = xx[0];
312  const Double_t y = xx[1];
313  const Double_t tmp1 = y-x*x;
314  const Double_t tmp2 = 1-x;
315  return 100*tmp1*tmp1+tmp2*tmp2;*/
316  //printf("arguments: %f %f %f %f %f", xx[0],xx[1],xx[2],xx[3],xx[4]);
317  double chi2 = 0.;
318  for (unsigned int i = 0; i < clusterpoints.size(); i++)
319  {
320  TrkrCluster *tc = clusterpoints.at(i);
321  double qz = xx[i + 5] - clusterpoints.at(0)->getZ();
322  double predictedx = xx[0] - sin(xx[2]) / xx[4] + cos(xx[4] / tan(xx[3]) * qz + xx[2] - M_PI / 2) / xx[4];
323  double predictedy = xx[1] + cos(xx[2]) / xx[4] + sin(xx[4] / tan(xx[3]) * qz + xx[2] - M_PI / 2) / xx[4];
324  double xresid = tc->getX() - predictedx;
325  double yresid = tc->getY() - predictedy;
326  double zresid = tc->getZ() - xx[i + 5];
327  //cout << "zresid=" << zresid << endl;
328  double contribution = xresid * xresid / tc->getError(0, 0) + yresid * yresid / tc->getError(1, 1) + zresid * zresid / tc->getError(2, 2);
329  //cout << "contribution of " << contribution << endl;
330  chi2 += contribution;
331  }
332  double qz = xx[clusterpoints.size() + 5] - _vertex->get_z();
333  double predictedx = xx[0] - sin(xx[2]) / xx[4] + cos(xx[4] / tan(xx[3]) * qz + xx[2] - M_PI / 2) / xx[4];
334  double predictedy = xx[1] + cos(xx[2]) / xx[4] + sin(xx[4] / tan(xx[3]) * qz + xx[2] - M_PI / 2) / xx[4];
335  double xresid = _vertex->get_x() - predictedx;
336  double yresid = _vertex->get_y() - predictedy;
337  double zresid = _vertex->get_z() - xx[clusterpoints.size() + 5];
338  //cout << "zresid=" << zresid << endl;
339  double contribution = xresid * xresid / _vertex->get_error(0, 0) + yresid * yresid / _vertex->get_error(1, 1) + zresid * zresid / _vertex->get_error(2, 2);
340  //cout << "contribution of " << contribution << endl;
341  chi2 += contribution;
342  return chi2;
343 }
344 
346 {
347  //---------------------------------
348  // Get Objects off of the Node Tree
349  //---------------------------------
350 
351  _cluster_map = findNode::getClass<TrkrClusterContainer>(topNode, "TRKR_CLUSTER");
352  if (!_cluster_map)
353  {
354  cerr << PHWHERE << " ERROR: Can't find node TRKR_CLUSTER" << endl;
356  }
358 }
359 
360 double PHRTreeSeeding::phiadd(double phi1, double phi2)
361 {
362  double s = phi1 + phi2;
363  if (s > 2 * M_PI)
364  return s - 2 * M_PI;
365  else if (s < 0)
366  return s + 2 * M_PI;
367  else
368  return s;
369 }
370 
371 double PHRTreeSeeding::phidiff(double phi1, double phi2)
372 {
373  double d = phi1 - phi2;
374  if (d > M_PI)
375  return d - 2 * M_PI;
376  else if (d < -M_PI)
377  return d + 2 * M_PI;
378  else
379  return d;
380 }
381 
382 void PHRTreeSeeding::QueryTree(const bgi::rtree<pointKey, bgi::quadratic<16>> &rtree, double phimin, double etamin, double lmin, double phimax, double etamax, double lmax, std::vector<pointKey> &returned_values)
383 {
384  rtree.query(bgi::intersects(box(point(phimin, etamin, lmin), point(phimax, etamax, lmax))), std::back_inserter(returned_values));
385  if (phimin < 0) rtree.query(bgi::intersects(box(point(2 * M_PI + phimin, etamin, lmin), point(2 * M_PI, etamax, lmax))), std::back_inserter(returned_values));
386  if (phimax > 2 * M_PI) rtree.query(bgi::intersects(box(point(0, etamin, lmin), point(phimax - 2 * M_PI, etamax, lmax))), std::back_inserter(returned_values));
387 }
388 
389 /*double PHRTreeSeeding::pointKeyToTuple(pointKey *pK)
390 {
391  TVector3 *v;
392  v->SetPtEtaPhi(_radii_all[(int)(pK->first.get<2>+0.4)], pK->first.get<1>, pK->first.get<0>);
393  return make_tuple(v.X(), v.Y(), v.Z());
394  }*/
395 
396 double /*PHRTreeSeeding::*/ costfunction(const double *xx)
397 // xcenter, ycenter, radius, frequency, phase
398 {
399  /*double cost = 0;
400  //should add a term for the vertex
401  for(vector<TrkrCluster*>::iterator i = clusterpoints.begin(); i != clusterpoints.end(); i++)
402  {
403  TrkrCluster* tc = *i;
404  double qa = tc->getError(0,0);
405  double qb = tc->getError(0,1);
406  double qc = tc->getError(1,0);
407  double qd = tc->getError(1,1);
408  double qz = tc->getZ();
409  double qx = tc->getX()-xx[0]-xx[2]*cos(xx[3]*qz+xx[4]);
410  double qy = tc->getY()-xx[1]-xx[2]*sin(xx[3]*qz+xx[4]);
411 
412  cost += (qd*qx*qx-(qb+qc)*qx*qy+qa*qy*qy)/(qa*qd-qb*qc);
413  //cost += pow(get<0>(clusterpoints.at(i))-xx[0]-xx[2]*cos(xx[3]*get<2>(clusterpoints.at(i))+xx[4]),2)+pow(get<1>(clusterpoints.at(i))-xx[1]-xx[2]*sin(xx[3]*get<2>(clusterpoints.at(i))+xx[4]),2)
414  }
415  return cost;*/
416  return xx[0] + xx[1] + xx[2] + xx[3] + xx[4];
417 }
418 
419 double PHRTreeSeeding::chisq(const double *xx)
420 {
421  double chi2 = 0.;
422  for (vector<TrkrCluster *>::iterator i = clusterpoints.begin(); i != clusterpoints.end(); ++i)
423  {
424  TrkrCluster *tc = *i;
425  double qa = tc->getError(0, 0);
426  double qb = tc->getError(0, 1);
427  double qc = tc->getError(1, 0);
428  double qd = tc->getError(1, 1);
429  double qz = tc->getZ() - clusterpoints.at(0)->getZ();
430  double predictedx = xx[0] - sin(xx[2]) / xx[4] + cos(xx[4] / tan(xx[3]) * qz + xx[2] - M_PI / 2) / xx[4];
431  double predictedy = xx[1] + cos(xx[2]) / xx[4] + sin(xx[4] / tan(xx[3]) * qz + xx[2] - M_PI / 2) / xx[4];
432  double qx = tc->getX() - predictedx;
433  double qy = tc->getY() - predictedy;
434 
435  chi2 += (qd * qx * qx - (qb + qc) * qx * qy + qa * qy * qy) / (qa * qd - qb * qc);
436  }
437 
438  return chi2;
439 }
440 
442 {
443  // bgi::rtree<pointKey, bgi::quadratic<16> > rtree;
444  PHTimer *t_fill = new PHTimer("t_fill");
445  t_fill->stop();
446  int n_dupli = 0;
447  int nlayer[60];
448  for (int j = 0; j < 60; j++) nlayer[j] = 0;
449 
450  auto hitsetrange = _hitsets->getHitSets(TrkrDefs::TrkrId::tpcId);
451  for (auto hitsetitr = hitsetrange.first;
452  hitsetitr != hitsetrange.second;
453  ++hitsetitr){
454  auto range = _cluster_map->getClusters(hitsetitr->first);
455  for( auto clusIter = range.first; clusIter != range.second; ++clusIter ){
456  TrkrCluster *cluster = clusIter->second;
457  TrkrDefs::cluskey ckey = clusIter->first;
458  unsigned int layer = TrkrDefs::getLayer(ckey);
459  if (layer < 39) continue;
460 
461  TVector3 vec(cluster->getPosition(0) - _vertex->get_x(), cluster->getPosition(1) - _vertex->get_y(), cluster->getPosition(2) - _vertex->get_z());
462 
463  double clus_phi = vec.Phi();
464  clus_phi -= 2 * M_PI * floor(clus_phi / (2 * M_PI));
465  double clus_eta = vec.Eta();
466  double clus_l = layer; // _radii_all[layer];
467 
468  vector<pointKey> testduplicate;
469  QueryTree(_rtree, clus_phi - 0.00001, clus_eta - 0.00001, layer - 0.5, clus_phi + 0.00001, clus_eta + 0.00001, layer + 0.5, testduplicate);
470  if (!testduplicate.empty())
471  {
472  n_dupli++;
473  continue;
474  }
475  nlayer[layer]++;
476  t_fill->restart();
477  _rtree.insert(std::make_pair(point(clus_phi, clus_eta, clus_l), ckey));
478  t_fill->stop();
479  }
480  }
481  std::cout << "fill time: " << t_fill->get_accumulated_time() / 1000. << " sec" << std::endl;
482  std::cout << "number of duplicates : " << n_dupli << std::endl;
483 }
484 
486 {
487  TFile fpara("rtree_para.root", "RECREATE");
488  TNtuple *NT = new TNtuple("NT", "NT", "pt:dpt:z:dz:phi:dphi:c:dc:nhit");
489  _vertex = _vertex_map->get(0);
490 
491  //for different purpose
492  phisr = 0.005 * _phi_scale; // *2
493  etasr = 0.0035 * _z_scale; // *2;
494  /* 0.7 version
495  phist = 0.001*1;
496  etast = 0.003*1;
497  0.9 version
498  phist = 0.001*2;
499  etast = 0.003*2;
500  */
501  phist = 0.001 * _phi_scale; // *5 *7;
502  etast = 0.003 * _z_scale; // *5/ *7;
503 
504  PHTimer *t_seed = new PHTimer("t_seed");
505  t_seed->stop();
506  t_seed->restart();
507 
508  _rtree.clear();
509  FillTree();
510 
511  int numberofseeds = 0;
512  cout << " entries in tree: " << _rtree.size() << endl;
513 
514  for (unsigned int iteration = 0; iteration < 1; ++iteration)
515  {
516  if (iteration == 1) _start_layer -= 7;
517  vector<pointKey> StartLayerClusters;
518  _rtree.query(bgi::intersects(
519  box(point(0, -3, ((double) _start_layer - 0.5)),
520  point(2 * M_PI, 3, ((double) _start_layer + 0.5)))),
521  std::back_inserter(StartLayerClusters));
522 
523  for (vector<pointKey>::iterator StartCluster = StartLayerClusters.begin(); StartCluster != StartLayerClusters.end(); ++StartCluster)
524  {
525  double StartPhi = StartCluster->first.get<0>();
526  double StartEta = StartCluster->first.get<1>();
527 
528  vector<pointKey> SecondLayerClusters;
530  StartPhi - phisr,
531  StartEta - etasr,
532  (double) _start_layer - 1.5,
533  StartPhi + phisr,
534  StartEta + etasr,
535  (double) _start_layer - 0.5,
536  SecondLayerClusters);
537  cout << " entries in second layer: " << SecondLayerClusters.size() << endl;
538 
539  for (vector<pointKey>::iterator SecondCluster = SecondLayerClusters.begin(); SecondCluster != SecondLayerClusters.end(); ++SecondCluster)
540  {
541  double currentphi = SecondCluster->first.get<0>();
542  double currenteta = SecondCluster->first.get<1>();
543  int lastgoodlayer = _start_layer - 1;
544 
545  int failures = 0;
546 
547  double dphidr = phidiff(StartPhi, currentphi) / (_radii_all[_start_layer] - _radii_all[_start_layer - 1]);
548  double detadr = (StartEta - currenteta) / (_radii_all[_start_layer] - _radii_all[_start_layer - 1]);
549 
550  double ther = (_radii_all[_start_layer] + _radii_all[_start_layer - 1]) / 2;
551 
552  vector<double> curvatureestimates;
553  curvatureestimates.push_back(copysign(2 / sqrt(ther * ther + 1 / dphidr / dphidr), dphidr));
554  vector<double> phi_zigzag;
555  phi_zigzag.push_back(dphidr);
556  vector<double> z_zigzag;
557  z_zigzag.push_back(detadr);
558 
559  vector<TrkrDefs::cluskey> cluskeys;
560  cluskeys.push_back(StartCluster->second);
561  cluskeys.push_back(SecondCluster->second);
562  cout << " phi 1: " << StartPhi
563  << " phi 2: " << currentphi
564  << " dphi: " << dphidr
565  << " eta 1: " << StartEta
566  << " eta 2: " << currenteta
567  << " deta: " << detadr
568  << endl;
569 
570  for (unsigned int newlayer = _start_layer - 2; newlayer >= (_start_layer - 7); newlayer--)
571  {
572  vector<pointKey> newlayer_clusters;
573  cout << " window - "
574  << " phimin " << currentphi - dphidr * (_radii_all[lastgoodlayer] - _radii_all[newlayer]) - phist
575  << " phimax " << currentphi - dphidr * (_radii_all[lastgoodlayer] - _radii_all[newlayer]) + phist
576  << " etamin " << currenteta - etast
577  << " etamax " << currenteta + etast
578  << endl;
580  currentphi - dphidr * (_radii_all[lastgoodlayer] - _radii_all[newlayer]) - phist,
581  currenteta - etast,
582  newlayer - 0.5,
583  currentphi - dphidr * (_radii_all[lastgoodlayer] - _radii_all[newlayer]) + phist,
584  currenteta + etast,
585  newlayer + 0.5,
586  newlayer_clusters);
587 
588  if (newlayer_clusters.empty())
589  {
590  failures += 1;
591  if (failures > 2) break; //0.7 2 0.9 3
592  }
593  else
594  {
595  double xinrecord = 100.0;
596  pointKey *xinkey = &*newlayer_clusters.begin();
597  for (std::vector<pointKey>::iterator it = newlayer_clusters.begin(); it != newlayer_clusters.end(); ++it)
598  {
599  double dist = abs(phidiff(it->first.get<0>(), currentphi - dphidr * (_radii_all[lastgoodlayer] - _radii_all[newlayer]))) + abs(it->first.get<1>() - currenteta);
600 
601  cout << " nuphi: " << it->first.get<0>()
602  << " nueta: " << it->first.get<1>()
603  << " dist: " << dist
604  << " lay: " << newlayer
605  << " dl: " << lastgoodlayer - newlayer
606  << " r: " << _radii_all[newlayer]
607  << " dr: " << _radii_all[lastgoodlayer] - _radii_all[newlayer]
608  << endl;
609  if (dist < xinrecord)
610  {
611  *xinkey = *it;
612  xinrecord = dist;
613  }
614  }
615 
616  dphidr = phidiff(currentphi, xinkey->first.get<0>()) / (_radii_all[lastgoodlayer] - _radii_all[newlayer]);
617  detadr = (currenteta - xinkey->first.get<1>()) / (_radii_all[lastgoodlayer] - _radii_all[newlayer]);
618  ther = (_radii_all[lastgoodlayer] - _radii_all[newlayer]) / 2;
619 
620  curvatureestimates.push_back(copysign(2 / sqrt(ther * ther + 1 / dphidr / dphidr), dphidr));
621 
622  phi_zigzag.push_back(dphidr);
623  z_zigzag.push_back(detadr);
624 
625  cluskeys.push_back(xinkey->second);
626 
627  currentphi = xinkey->first.get<0>();
628  currenteta = (currenteta + xinkey->first.get<1>()) / 2;
629 
630  lastgoodlayer = newlayer;
631  }
632  }
633  if (failures > 2) continue; //0.7 2 0.9 3
634 
635  double phi_sum = std::accumulate(phi_zigzag.begin(), phi_zigzag.end(), 0.0);
636  double phi_mean = phi_sum / phi_zigzag.size();
637 
638  std::vector<double> phi_diff(phi_zigzag.size());
639  std::transform(phi_zigzag.begin(), phi_zigzag.end(), phi_diff.begin(),
640  std::bind2nd(std::minus<double>(), phi_mean));
641  double phi_sq_sum = std::inner_product(phi_diff.begin(), phi_diff.end(), phi_diff.begin(), 0.0);
642  double phi_stdev = std::sqrt(phi_sq_sum / (phi_zigzag.size() - 1));
643 
644  double z_sum = std::accumulate(z_zigzag.begin(), z_zigzag.end(), 0.0);
645  double z_mean = z_sum / z_zigzag.size();
646 
647  std::vector<double> z_diff(z_zigzag.size());
648  std::transform(z_zigzag.begin(), z_zigzag.end(), z_diff.begin(),
649  std::bind2nd(std::minus<double>(), z_mean));
650  double z_sq_sum = std::inner_product(z_diff.begin(), z_diff.end(), z_diff.begin(), 0.0);
651  double z_stdev = std::sqrt(z_sq_sum / (z_zigzag.size() - 1));
652 
653  double curv_sum = std::accumulate(curvatureestimates.begin(), curvatureestimates.end(), 0.0);
654  double curv_mean = curv_sum / curvatureestimates.size();
655 
656  std::vector<double> curv_diff(curvatureestimates.size());
657  std::transform(curvatureestimates.begin(), curvatureestimates.end(), curv_diff.begin(),
658  std::bind2nd(std::minus<double>(), curv_mean));
659  double curv_sq_sum = std::inner_product(curv_diff.begin(), curv_diff.end(), curv_diff.begin(), 0.0);
660  double curv_stdev = std::sqrt(curv_sq_sum / (curvatureestimates.size() - 1));
661 
662  const double BQ = 0.01 * 1.4 * 0.299792458;
663  double pt = BQ / abs(curv_mean);
664  double pterror = BQ * curv_stdev / (curv_mean * curv_mean);
665 
666  // pt:z:dz:phi:dphi:c:dc
667  NT->Fill(pt, pterror, z_mean, z_stdev, phi_mean, phi_stdev, curv_mean, curv_stdev, cluskeys.size());
669  track.set_id(numberofseeds);
670 
671  for (unsigned int j = 0; j < cluskeys.size(); j++)
672  {
673  track.insert_cluster_key(cluskeys.at(j));
674  }
675  track.set_ndf(2 * cluskeys.size() - 5);
676  short int helicity = 1;
677  if (StartPhi * curv_mean < 0) helicity -= 2;
678  track.set_charge(-helicity);
679 
680  TrkrCluster *cl = _cluster_map->findCluster(StartCluster->second);
681  track.set_x(_vertex->get_x()); //track.set_x(cl->getX());
682  track.set_y(_vertex->get_y()); //track.set_y(cl->getY());
683  track.set_z(_vertex->get_z()); //track.set_z(cl->getZ());
684  track.set_px(pt * cos(StartPhi));
685  track.set_py(pt * sin(StartPhi));
686  track.set_pz(pt / tan(2 * atan(exp(-StartEta))));
687  track.set_error(0, 0, cl->getError(0, 0));
688  track.set_error(0, 1, cl->getError(0, 1));
689  track.set_error(0, 2, cl->getError(0, 2));
690  track.set_error(1, 1, cl->getError(1, 1));
691  track.set_error(1, 2, cl->getError(1, 2));
692  track.set_error(2, 2, cl->getError(2, 2));
693  track.set_error(3, 3, pterror * pterror * cos(StartPhi) * cos(StartPhi));
694  track.set_error(4, 4, pterror * pterror * sin(StartPhi) * sin(StartPhi));
695  track.set_error(5, 5, pterror * pterror / tan(2 * atan(exp(-StartEta))) / tan(2 * atan(exp(-StartEta))));
696  _track_map->insert(&track);
697 
698  numberofseeds++;
699  }
700  }
701  }
702  t_seed->stop();
703  cout << "number of seeds " << numberofseeds << endl;
704  cout << "seeding time: " << t_seed->get_accumulated_time() / 1000 << " s" << endl;
705 
706  /*
707 
708  */
709 
710  fpara.cd();
711  NT->Write();
712  fpara.Close();
713 
715 }
716 
718 {
719  cout << "Called Setup" << endl;
720  cout << "topNode:" << topNode << endl;
721  PHTrackSeeding::Setup(topNode);
722  // int ret = GetNodes(topNode);
723  //return ret;
724  GetNodes(topNode);
725  InitializeGeometry(topNode);
727 }
728 
730 {
731  cout << "Called End " << endl;
733 }