ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
AtlasSeedfinder.ipp
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file AtlasSeedfinder.ipp
1 // This file is part of the Acts project.
2 //
3 // Copyright (C) 2018 CERN for the benefit of the Acts project
4 //
5 // This Source Code Form is subject to the terms of the Mozilla Public
6 // License, v. 2.0. If a copy of the MPL was not distributed with this
7 // file, You can obtain one at http://mozilla.org/MPL/2.0/.
8 
10 // AtlasSeedfinder.ipp Acts project
12 
13 #include <algorithm>
14 
16 // Constructor
18 
19 template <typename SpacePoint>
21  m_checketa = false;
22  m_maxsize = 50000;
23  m_ptmin = 400.;
24  m_etamin = 0.;
25  m_etamax = 2.7;
26  // delta R minimum & maximum within a seed
27  m_drmin = 5.;
28  m_drminv = 20.;
29  m_drmax = 270.;
30  // restrict z coordinate of particle origin
31  m_zmin = -250.;
32  m_zmax = +250.;
33  // radius of detector in mm
34  r_rmax = 600.;
35  r_rstep = 2.;
36 
37  // checking if Z is compatible:
38  // m_dzver is related to delta-Z
39  // m_dzdrver is related to delta-Z divided by delta-R
40  m_dzver = 5.;
41  m_dzdrver = .02;
42 
43  // shift all spacepoints by this offset such that the beam can be assumed to
44  // be at 0,0
45  // z shift should not matter as beam is assumed to be parallel to central
46  // detector axis,
47  // but spacepoints will be shifted by z as well anyway.
48  m_xbeam = 0.;
49  m_ybeam = 0.;
50  m_zbeam = 0.;
51 
52  // config
53  // max impact parameters
54  // m_diver = max ppp impact params
55  m_diver = 10.;
56  m_diverpps = 1.7;
57  m_diversss = 50;
58  m_divermax = 20.;
59  // delta azimuth (phi)
60  m_dazmax = .02;
61  // limit for sp compatible with 1 middle space point
62  // ridiculously large to be EXTRA safe
63  // only actually keep 5 of these max 5000 (m_maxOneSize of m_maxsizeSP)
64  m_maxsizeSP = 5000;
65  m_maxOneSize = 5;
66 
67  // cache: counting if ran already
68  m_state = 0;
69 
70  m_nlist = 0;
71  m_endlist = true;
72  r_Sorted = 0;
73  r_index = 0;
74  r_map = 0;
75  m_SP = 0;
76  m_R = 0;
77  m_Tz = 0;
78  m_Er = 0;
79  m_U = 0;
80  m_V = 0;
81  m_Zo = 0;
82  m_OneSeeds = 0;
83  m_seedOutput = 0;
84 
85  // Build framework
86  //
87  buildFrameWork();
88  m_CmSp.reserve(500);
89 }
90 
92 // Destructor
94 template <class SpacePoint>
96  if (r_index) {
97  delete[] r_index;
98  }
99  if (r_map) {
100  delete[] r_map;
101  }
102  if (r_Sorted) {
103  delete[] r_Sorted;
104  }
105 
106  // Delete seeds
107  //
108  for (i_seed = l_seeds.begin(); i_seed != l_seeds.end(); ++i_seed) {
109  delete *i_seed;
110  }
111  // Delete space points for reconstruction
112  //
113  i_spforseed = l_spforseed.begin();
114  for (; i_spforseed != l_spforseed.end(); ++i_spforseed) {
115  delete *i_spforseed;
116  }
117  if (m_seedOutput) {
118  delete m_seedOutput;
119  }
120 
121  if (m_SP) {
122  delete[] m_SP;
123  }
124  if (m_R) {
125  delete[] m_R;
126  }
127  if (m_Tz) {
128  delete[] m_Tz;
129  }
130  if (m_Er) {
131  delete[] m_Er;
132  }
133  if (m_U) {
134  delete[] m_U;
135  }
136  if (m_V) {
137  delete[] m_V;
138  }
139  if (m_Zo) {
140  delete[] m_Zo;
141  }
142  if (m_OneSeeds) {
143  delete[] m_OneSeeds;
144  }
145 }
146 
148 // Initialize tool for new event
150 template <typename SpacePoint>
151 template <class RandIter>
153  RandIter spBegin,
154  RandIter spEnd) {
155  iteration <= 0 ? m_iteration = 0 : m_iteration = iteration;
156  erase();
157  m_dzdrmin = m_dzdrmin0;
158  m_dzdrmax = m_dzdrmax0;
159  m_umax = 100.;
160  // if first event
161  r_first = 0;
162  if (!m_iteration) {
163  buildBeamFrameWork();
164 
165  // curvature depending on bfield
166  m_K = 2. / (300. * m_config.bFieldInZ);
167  // curvature of minimum pT track
168  m_ipt2K = m_ipt2 / (m_K * m_K);
169  // scattering of min pT track
170  m_ipt2C = m_ipt2 * m_COF;
171  // scattering times curvature (missing: div by pT)
172  m_COFK = m_COF * (m_K * m_K);
173  i_spforseed = l_spforseed.begin();
174  }
175  // only if not first event
176  else {
177  fillLists();
178  return;
179  }
180 
181  float irstep = 1. / r_rstep;
182  int irmax = r_size - 1;
183  // TODO using 3 member vars to clear r_Sorted
184  for (int i = 0; i != m_nr; ++i) {
185  int n = r_index[i];
186  r_map[n] = 0;
187  r_Sorted[n].clear();
188  }
189  m_nr = 0;
190 
191  // convert space-points to SPForSeed and sort into radius-binned array
192  // r_Sorted
193  // store number of SP per bin in r_map
194  RandIter sp = spBegin;
195  for (; sp != spEnd; ++sp) {
196  Acts::Legacy::SPForSeed<SpacePoint>* sps = newSpacePoint((*sp));
197  if (!sps) {
198  continue;
199  }
200  int ir = int(sps->radius() * irstep);
201  if (ir > irmax) {
202  ir = irmax;
203  }
204  r_Sorted[ir].push_back(sps);
205  // if there is exactly one SP in current bin, add bin nr in r_index (for
206  // deletion later)
207  // TODO overly complicated
208  ++r_map[ir];
209  if (r_map[ir] == 1) {
210  r_index[m_nr++] = ir;
211  }
212  }
213 
214  fillLists();
215 }
216 
218 // Methods to initialize different strategies of seeds production
219 // with three space points with or without vertex constraint
221 
222 template <class SpacePoint>
224  m_zminU = m_zmin;
225  m_zmaxU = m_zmax;
226 
227  if ((m_state == 0) || m_nlist) {
228  i_seede = l_seeds.begin();
229  m_state = 1;
230  m_nlist = 0;
231  m_endlist = true;
232  m_fvNmin = 0;
233  m_fNmin = 0;
234  m_zMin = 0;
235  production3Sp();
236  }
237  i_seed = l_seeds.begin();
238  m_seed = m_seeds.begin();
239 }
240 
242 // Find next set space points
244 template <class SpacePoint>
246  if (m_endlist) {
247  return;
248  }
249 
250  i_seede = l_seeds.begin();
251 
252  production3Sp();
253 
254  i_seed = l_seeds.begin();
255  m_seed = m_seeds.begin();
256  ++m_nlist;
257 }
258 
260 // Initiate frame work for seed generator
262 template <class SpacePoint>
264  m_ptmin = fabs(m_ptmin);
265 
266  if (m_ptmin < 100.) {
267  m_ptmin = 100.;
268  }
269 
270  if (m_diversss < m_diver) {
271  m_diversss = m_diver;
272  }
273  if (m_divermax < m_diversss) {
274  m_divermax = m_diversss;
275  }
276 
277  if (fabs(m_etamin) < .1) {
278  m_etamin = -m_etamax;
279  }
280  m_dzdrmax0 = 1. / tan(2. * atan(exp(-m_etamax)));
281  m_dzdrmin0 = 1. / tan(2. * atan(exp(-m_etamin)));
282 
283  // scattering factor. depends on error, forward direction and distance between
284  // SP
285  m_COF = 134 * .05 * 9.;
286  m_ipt = 1. / fabs(.9 * m_ptmin);
287  m_ipt2 = m_ipt * m_ipt;
288  m_K = 0.;
289 
290  // set all counters zero
291  m_nsaz = m_nsazv = m_nr = m_nrfz = 0;
292 
293  // Build radius sorted containers
294  //
295  r_size = int((r_rmax + .1) / r_rstep);
296  r_Sorted = new std::list<Acts::Legacy::SPForSeed<SpacePoint>*>[r_size];
297  r_index = new int[r_size];
298  r_map = new int[r_size];
299  for (int i = 0; i != r_size; ++i) {
300  r_index[i] = 0;
301  r_map[i] = 0;
302  }
303 
304  // Build radius-azimuthal sorted containers
305  //
306  const float pi2 = 2. * M_PI;
307  const int NFmax = 53;
308  const float sFmax = float(NFmax) / pi2;
309  const float m_sFmin = 100. / 60.;
310  // make phi-slices for 400MeV tracks, unless ptMin is even smaller
311  float ptm = 400.;
312  if (m_ptmin < ptm) {
313  ptm = m_ptmin;
314  }
315 
316  m_sF = ptm / 60.;
317  if (m_sF > sFmax) {
318  m_sF = sFmax;
319  } else if (m_sF < m_sFmin) {
320  m_sF = m_sFmin;
321  }
322  m_fNmax = int(pi2 * m_sF);
323  if (m_fNmax >= NFmax) {
324  m_fNmax = NFmax - 1;
325  }
326 
327  // Build radius-azimuthal-Z sorted containers
328  //
329  m_nrfz = 0;
330  for (int i = 0; i != 583; ++i) {
331  rfz_index[i] = 0;
332  rfz_map[i] = 0;
333  }
334 
335  // Build maps for radius-azimuthal-Z sorted collections
336  //
337  for (int f = 0; f <= m_fNmax; ++f) {
338  int fb = f - 1;
339  if (fb < 0) {
340  fb = m_fNmax;
341  }
342  int ft = f + 1;
343  if (ft > m_fNmax) {
344  ft = 0;
345  }
346 
347  // For each azimuthal region loop through all Z regions
348  //
349  for (int z = 0; z != 11; ++z) {
350  int a = f * 11 + z;
351  int b = fb * 11 + z;
352  int c = ft * 11 + z;
353  rfz_b[a] = 3;
354  rfz_t[a] = 3;
355  rfz_ib[a][0] = a;
356  rfz_it[a][0] = a;
357  rfz_ib[a][1] = b;
358  rfz_it[a][1] = b;
359  rfz_ib[a][2] = c;
360  rfz_it[a][2] = c;
361  if (z == 5) {
362  rfz_t[a] = 9;
363  rfz_it[a][3] = a + 1;
364  rfz_it[a][4] = b + 1;
365  rfz_it[a][5] = c + 1;
366  rfz_it[a][6] = a - 1;
367  rfz_it[a][7] = b - 1;
368  rfz_it[a][8] = c - 1;
369  } else if (z > 5) {
370  rfz_b[a] = 6;
371  rfz_ib[a][3] = a - 1;
372  rfz_ib[a][4] = b - 1;
373  rfz_ib[a][5] = c - 1;
374 
375  if (z < 10) {
376  rfz_t[a] = 6;
377  rfz_it[a][3] = a + 1;
378  rfz_it[a][4] = b + 1;
379  rfz_it[a][5] = c + 1;
380  }
381  } else {
382  rfz_b[a] = 6;
383  rfz_ib[a][3] = a + 1;
384  rfz_ib[a][4] = b + 1;
385  rfz_ib[a][5] = c + 1;
386 
387  if (z > 0) {
388  rfz_t[a] = 6;
389  rfz_it[a][3] = a - 1;
390  rfz_it[a][4] = b - 1;
391  rfz_it[a][5] = c - 1;
392  }
393  }
394 
395  if (z == 3) {
396  rfz_b[a] = 9;
397  rfz_ib[a][6] = a + 2;
398  rfz_ib[a][7] = b + 2;
399  rfz_ib[a][8] = c + 2;
400  } else if (z == 7) {
401  rfz_b[a] = 9;
402  rfz_ib[a][6] = a - 2;
403  rfz_ib[a][7] = b - 2;
404  rfz_ib[a][8] = c - 2;
405  }
406  }
407  }
408 
409  if (!m_SP) {
410  m_SP = new Acts::Legacy::SPForSeed<SpacePoint>*[m_maxsizeSP];
411  }
412  if (m_R == nullptr) {
413  m_R = new float[m_maxsizeSP];
414  }
415  if (m_Tz == nullptr) {
416  m_Tz = new float[m_maxsizeSP];
417  }
418  if (m_Er == nullptr) {
419  m_Er = new float[m_maxsizeSP];
420  }
421  if (m_U == nullptr) {
422  m_U = new float[m_maxsizeSP];
423  }
424  if (m_V == nullptr) {
425  m_V = new float[m_maxsizeSP];
426  }
427  if (m_Zo == nullptr) {
428  m_Zo = new float[m_maxsizeSP];
429  }
430  if (!m_OneSeeds) {
431  m_OneSeeds = new Acts::Legacy::InternalSeed<SpacePoint>[m_maxOneSize];
432  }
433 
434  if (!m_seedOutput) {
435  m_seedOutput = new Acts::Legacy::Seed<SpacePoint>();
436  }
437 
438  i_seed = l_seeds.begin();
439  i_seede = l_seeds.end();
440 }
441 
443 // Initiate beam frame work for seed generator
445 template <class SpacePoint>
447  double bx = m_config.beamPosX;
448  double by = m_config.beamPosY;
449  double bz = m_config.beamPosZ;
450 
451  m_xbeam = float(bx);
452  m_ybeam = float(by);
453  m_zbeam = float(bz);
454 }
455 
457 // Initiate beam frame work for seed generator
459 template <class SpacePoint>
461  SpacePoint* const& sp, float* r) {
462  r[0] = float(sp->x) - m_xbeam;
463  r[1] = float(sp->y) - m_ybeam;
464  r[2] = float(sp->z) - m_zbeam;
465 }
466 
468 // Initiate space points seed maker
470 template <class SpacePoint>
472  const float pi2 = 2. * M_PI;
473  typename std::list<Acts::Legacy::SPForSeed<SpacePoint>*>::iterator r, re;
474 
475  int ir0 = 0;
476  bool ibl = false;
477 
478  r_first = 0;
479  if (m_iteration) {
480  r_first = m_config.SCT_rMin / r_rstep;
481  }
482  for (int i = r_first; i != r_size; ++i) {
483  if (!r_map[i]) {
484  continue;
485  }
486 
487  r = r_Sorted[i].begin();
488  re = r_Sorted[i].end();
489 
490  if (ir0 == 0) {
491  ir0 = i;
492  }
493  // if not 1st event
494  if (m_iteration) {
495  //
496  if (!(*r)->spacepoint->clusterList().second) {
497  if (i < 20) {
498  ibl = true;
499  }
500  } else if (ibl) {
501  break;
502  } else if (i > 175) {
503  break;
504  }
505  }
506 
507  for (; r != re; ++r) {
508  // Azimuthal (Phi) angle sort
509  // bin nr "f" is phi * phi-slices
510  float F = (*r)->phi();
511  if (F < 0.) {
512  F += pi2;
513  }
514 
515  int f = int(F * m_sF);
516  if (f < 0) {
517  f = m_fNmax;
518  } else if (f > m_fNmax) {
519  f = 0;
520  }
521 
522  int z;
523  float Z = (*r)->z();
524 
525  // Azimuthal angle and Z-coordinate sort
526  // assign z-bin a value between 0 and 10 identifying the z-slice of a
527  // space-point
528  if (Z > 0.) {
529  Z < 250. ? z = 5
530  : Z < 450. ? z = 6
531  : Z < 925. ? z = 7
532  : Z < 1400. ? z = 8
533  : Z < 2500. ? z = 9 : z = 10;
534  } else {
535  Z > -250.
536  ? z = 5
537  : Z > -450.
538  ? z = 4
539  : Z > -925. ? z = 3
540  : Z > -1400. ? z = 2 : Z > -2500. ? z = 1 : z = 0;
541  }
542  // calculate bin nr "n" for self made r-phi-z sorted 3D array "rfz_Sorted"
543  // record number of sp in m_nsaz
544  int n = f * 11 + z;
545  ++m_nsaz;
546  // push back sp into rfz_Sorted, record new number of entries in bin in
547  // rfz_map,
548  // if 1st entry record non-empty bin in "rfz_index"
549  rfz_Sorted[n].push_back(*r);
550  if (!rfz_map[n]++) {
551  rfz_index[m_nrfz++] = n;
552  }
553  }
554  }
555  m_state = 0;
556 }
557 
559 // Erase space point information
561 template <class SpacePoint>
563  for (int i = 0; i != m_nrfz; ++i) {
564  int n = rfz_index[i];
565  rfz_map[n] = 0;
566  rfz_Sorted[n].clear();
567  }
568 
569  m_state = 0;
570  m_nsaz = 0;
571  m_nsazv = 0;
572  m_nrfz = 0;
573 }
574 
576 // Production 3 space points seeds
578 template <class SpacePoint>
580  // if less than 3 sp in total
581  if (m_nsaz < 3) {
582  return;
583  }
584  m_seeds.clear();
585 
586  // indices for the z-regions array in weird order.
587  // ensures creating seeds first for barrel, then left EC, then right EC
588  const int ZI[11] = {5, 6, 7, 8, 9, 10, 4, 3, 2, 1, 0};
589  typename std::list<Acts::Legacy::SPForSeed<SpacePoint>*>::iterator rt[9],
590  rte[9], rb[9], rbe[9];
591  int nseed = 0;
592 
593  // Loop through all azimuthal regions
594  //
595  for (int f = m_fNmin; f <= m_fNmax; ++f) {
596  // For each azimuthal region loop through all Z regions
597  // first for barrel, then left EC, then right EC
598  int z = 0;
599  if (!m_endlist) {
600  z = m_zMin;
601  }
602  for (; z != 11; ++z) {
603  int a = f * 11 + ZI[z];
604  if (!rfz_map[a]) {
605  continue;
606  }
607  int NB = 0, NT = 0;
608  for (int i = 0; i != rfz_b[a]; ++i) {
609  int an = rfz_ib[a][i];
610  // if bin has no entry: continue
611  if (!rfz_map[an]) {
612  continue;
613  }
614  // assign begin-pointer and end-pointer of current bin to rb and rbe
615  rb[NB] = rfz_Sorted[an].begin();
616  rbe[NB++] = rfz_Sorted[an].end();
617  }
618  for (int i = 0; i != rfz_t[a]; ++i) {
619  int an = rfz_it[a][i];
620  // if bin has no entry: continue
621  if (!rfz_map[an]) {
622  continue;
623  }
624  // assign begin-pointer and end-pointer of current bin to rt and rte
625  rt[NT] = rfz_Sorted[an].begin();
626  rte[NT++] = rfz_Sorted[an].end();
627  }
628  production3Sp(rb, rbe, rt, rte, NB, NT, nseed);
629  if (!m_endlist) {
630  m_fNmin = f;
631  m_zMin = z;
632  return;
633  }
634  }
635  }
636  m_endlist = true;
637 }
638 
640 // Production 3 space points seeds for full scan
642 template <class SpacePoint>
644  typename std::list<Acts::Legacy::SPForSeed<SpacePoint>*>::iterator* rb,
645  typename std::list<Acts::Legacy::SPForSeed<SpacePoint>*>::iterator* rbe,
646  typename std::list<Acts::Legacy::SPForSeed<SpacePoint>*>::iterator* rt,
647  typename std::list<Acts::Legacy::SPForSeed<SpacePoint>*>::iterator* rte,
648  int NB, int NT, int& nseed) {
649  typename std::list<Acts::Legacy::SPForSeed<SpacePoint>*>::iterator r0 = rb[0],
650  r;
651  if (!m_endlist) {
652  r0 = m_rMin;
653  m_endlist = true;
654  }
655 
656  float ipt2K = m_ipt2K;
657  float ipt2C = m_ipt2C;
658  float COFK = m_COFK;
659  float imaxp = m_diver;
660  float imaxs = m_divermax;
661 
662  m_CmSp.clear();
663 
664  // Loop through all space points
665  // first bottom bin used as "current bin" for middle spacepoints
666  for (; r0 != rbe[0]; ++r0) {
667  m_nOneSeeds = 0;
668  m_mapOneSeeds.clear();
669 
670  float R = (*r0)->radius();
671 
672  const int sur0 = (*r0)->surface();
673  float X = (*r0)->x();
674  float Y = (*r0)->y();
675  float Z = (*r0)->z();
676  int Nb = 0;
677 
678  // Bottom links production
679  //
680  for (int i = 0; i != NB; ++i) {
681  for (r = rb[i]; r != rbe[i]; ++r) {
682  float Rb = (*r)->radius();
683  float dR = R - Rb;
684  // if deltaR larger than deltaRMax, store spacepoint counter position in
685  // rb[i]
686  // this is not necessary at all because r is the loop counter (rb never
687  // read again)
688  if (dR > m_drmax) {
689  rb[i] = r;
690  continue;
691  }
692  // if dR is smaller than drmin, stop processing this bottombin
693  // because it has no more sp with lower radii
694  // OR break because processing PPP and encountered SCT SP
695  if (dR < m_drmin ||
696  (m_iteration && (*r)->spacepoint->clusterList().second)) {
697  break;
698  }
699  if ((*r)->surface() == sur0) {
700  continue;
701  }
702  // forward direction of SP duplet
703  float Tz = (Z - (*r)->z()) / dR;
704  float aTz = fabs(Tz);
705  // why also exclude seeds with small pseudorapidity??
706  if (aTz < m_dzdrmin || aTz > m_dzdrmax) {
707  continue;
708  }
709 
710  // Comparison with vertices Z coordinates
711  // continue if duplet origin not within collision region on z axis
712  float Zo = Z - R * Tz;
713  if (!isZCompatible(Zo)) {
714  continue;
715  }
716  m_SP[Nb] = (*r);
717  if (++Nb == m_maxsizeSP) {
718  goto breakb;
719  }
720  }
721  }
722  breakb:
723  if ((Nb == 0) || Nb == m_maxsizeSP) {
724  continue;
725  }
726  int Nt = Nb;
727 
728  // Top links production
729  //
730  for (int i = 0; i != NT; ++i) {
731  for (r = rt[i]; r != rte[i]; ++r) {
732  float Rt = (*r)->radius();
733  float dR = Rt - R;
734 
735  if (dR < m_drmin) {
736  rt[i] = r;
737  continue;
738  }
739  if (dR > m_drmax) {
740  break;
741  }
742  // TODO: is check if in same surface necessary?
743  if ((*r)->surface() == sur0) {
744  continue;
745  }
746 
747  float Tz = ((*r)->z() - Z) / dR;
748  float aTz = fabs(Tz);
749 
750  if (aTz < m_dzdrmin || aTz > m_dzdrmax) {
751  continue;
752  }
753 
754  // Comparison with vertices Z coordinates
755  //
756  float Zo = Z - R * Tz;
757  if (!isZCompatible(Zo)) {
758  continue;
759  }
760  m_SP[Nt] = (*r);
761  if (++Nt == m_maxsizeSP) {
762  goto breakt;
763  }
764  }
765  }
766 
767  breakt:
768  if ((Nt - Nb) == 0) {
769  continue;
770  }
771  float covr0 = (*r0)->covr();
772  float covz0 = (*r0)->covz();
773  float ax = X / R;
774  float ay = Y / R;
775 
776  for (int i = 0; i != Nt; ++i) {
778 
779  float dx = sp->x() - X;
780  float dy = sp->y() - Y;
781  float dz = sp->z() - Z;
782  // calculate projection fraction of spM->spT vector pointing in same
783  // direction as
784  // vector origin->spM (x) and projection fraction of spM->spT vector
785  // pointing
786  // orthogonal to origin->spM (y)
787  float x = dx * ax + dy * ay;
788  float y = dy * ax - dx * ay;
789  // 1/(r*r) = 1/dx*dx+dy*dy = 1/(distance between the two points)^2
790  float r2 = 1. / (x * x + y * y);
791  // 1/r
792  float dr = sqrt(r2);
793  float tz = dz * dr;
794  if (i < Nb) {
795  tz = -tz;
796  }
797 
798  m_Tz[i] = tz;
799  m_Zo[i] = Z - R * tz;
800  m_R[i] = dr;
801  m_U[i] = x * r2;
802  m_V[i] = y * r2;
803  m_Er[i] = ((covz0 + sp->covz()) + (tz * tz) * (covr0 + sp->covr())) * r2;
804  }
805  covr0 *= .5;
806  covz0 *= 2.;
807 
808  // Three space points comparison
809  //
810  for (int b = 0; b != Nb; ++b) {
811  float Zob = m_Zo[b];
812  float Tzb = m_Tz[b];
813  float Rb2r = m_R[b] * covr0;
814  float Rb2z = m_R[b] * covz0;
815  float Erb = m_Er[b];
816  float Vb = m_V[b];
817  float Ub = m_U[b];
818  // Tzb2 = 1/sin^2(theta)
819  float Tzb2 = (1. + Tzb * Tzb);
820  float sTzb2 = sqrt(Tzb2);
821  // CSA = curvature^2/(pT^2 * sin^2(theta)) * scatteringFactor
822  float CSA = Tzb2 * COFK;
823  // ICSA = (1/(pT^2 * sin^2(theta)) * scatteringFactor
824  float ICSA = Tzb2 * ipt2C;
825  float imax = imaxp;
826  if (m_SP[b]->spacepoint->clusterList().second) {
827  imax = imaxs;
828  }
829 
830  for (int t = Nb; t != Nt; ++t) {
831  float dT = ((Tzb - m_Tz[t]) * (Tzb - m_Tz[t]) - m_R[t] * Rb2z -
832  (Erb + m_Er[t])) -
833  (m_R[t] * Rb2r) * ((Tzb + m_Tz[t]) * (Tzb + m_Tz[t]));
834  if (dT > ICSA) {
835  continue;
836  }
837 
838  float dU = m_U[t] - Ub;
839  if (dU == 0.) {
840  continue;
841  }
842  float A = (m_V[t] - Vb) / dU;
843  float S2 = 1. + A * A;
844  float B = Vb - A * Ub;
845  float B2 = B * B;
846  // B2/S2=1/helixradius^2
847  if (B2 > ipt2K * S2 || dT * S2 > B2 * CSA) {
848  continue;
849  }
850 
851  float Im = fabs((A - B * R) * R);
852 
853  if (Im <= imax) {
854  // Add penalty factor dependent on difference between cot(theta) to
855  // the quality Im (previously Impact)
856  float dr;
857  m_R[t] < m_R[b] ? dr = m_R[t] : dr = m_R[b];
858  Im += fabs((Tzb - m_Tz[t]) / (dr * sTzb2));
859  // B/sqrt(S2) = 1/helixradius
860  m_CmSp.push_back(std::make_pair(B / sqrt(S2), m_SP[t]));
861  m_SP[t]->setParam(Im);
862  }
863  }
864  if (!m_CmSp.empty()) {
865  newOneSeedWithCurvaturesComparison(m_SP[b], (*r0), Zob);
866  }
867  }
868  fillSeeds();
869  nseed += m_fillOneSeeds;
870  if (nseed >= m_maxsize) {
871  m_endlist = false;
872  ++r0;
873  m_rMin = r0;
874  return;
875  }
876  }
877 }
878 
880 // New 3 space points pro seeds
882 template <class SpacePoint>
886  Acts::Legacy::SPForSeed<SpacePoint>*& p3, float z, float q) {
887  // if the number of seeds already in m_OneSeeds does not exceed m_maxOneSize
888  // then insert the current SP into m_mapOneSeeds and m_OneSeeds.
889  if (m_nOneSeeds < m_maxOneSize) {
890  m_OneSeeds[m_nOneSeeds].set(p1, p2, p3, z);
891  m_mapOneSeeds.insert(std::make_pair(q, m_OneSeeds + m_nOneSeeds));
892  ++m_nOneSeeds;
893  }
894  // if not, check the q(uality) of the LAST seed of m_mapOneSeeds
895  // if the quality of the new seed is worse, replace the value this
896  // (better) quality-key pointed to with the SP from the new seed.
897  // Then, add the new seed a SECOND time to the map with the worse quality as
898  // key.
899  // Then remove all entries after the newly inserted entry equal to the new
900  // seed.
901  else {
902  typename std::multimap<
903  float, Acts::Legacy::InternalSeed<SpacePoint>*>::reverse_iterator l =
904  m_mapOneSeeds.rbegin();
905 
906  if ((*l).first <= q) {
907  return;
908  }
909 
911  s->set(p1, p2, p3, z);
912 
913  typename std::multimap<
914  float, Acts::Legacy::InternalSeed<SpacePoint>*>::iterator i =
915  m_mapOneSeeds.insert(std::make_pair(q, s));
916 
917  for (++i; i != m_mapOneSeeds.end(); ++i) {
918  if ((*i).second == s) {
919  m_mapOneSeeds.erase(i);
920  return;
921  }
922  }
923  }
924 }
925 
927 // New 3 space points pro seeds production
929 template <class SpacePoint>
933  Acts::Legacy::SPForSeed<SpacePoint>*& SP0, float Zob) {
934  // allowed (1/helixradius)-delta between 2 seeds
935  const float dC = .00003;
936 
937  bool pixb = !SPb->spacepoint->clusterList().second;
938  float ub = SPb->quality();
939  float u0 = SP0->quality();
940 
941  std::sort(m_CmSp.begin(), m_CmSp.end(), Acts::Legacy::comCurvature());
942  typename std::vector<
943  std::pair<float, Acts::Legacy::SPForSeed<SpacePoint>*>>::iterator j,
944  jn, i = m_CmSp.begin(), ie = m_CmSp.end();
945  jn = i;
946 
947  for (; i != ie; ++i) {
948  float u = (*i).second->param();
949  float Im = (*i).second->param();
950 
951  bool pixt = !(*i).second->spacepoint->clusterList().second;
952 
953  const int Sui = (*i).second->surface();
954  float Ri = (*i).second->radius();
955  float Ci1 = (*i).first - dC;
956  float Ci2 = (*i).first + dC;
957  float Rmi = 0.;
958  float Rma = 0.;
959  bool in = false;
960 
961  if (!pixb) {
962  u -= 400.;
963  } else if (pixt) {
964  u -= 200.;
965  }
966 
967  for (j = jn; j != ie; ++j) {
968  if (j == i) {
969  continue;
970  }
971  if ((*j).first < Ci1) {
972  jn = j;
973  ++jn;
974  continue;
975  }
976  if ((*j).first > Ci2) {
977  break;
978  }
979  if ((*j).second->surface() == Sui) {
980  continue;
981  }
982  // Compared seeds should have at least deltaRMin distance
983  float Rj = (*j).second->radius();
984  if (fabs(Rj - Ri) < m_drmin) {
985  continue;
986  }
987 
988  if (in) {
989  if (Rj > Rma) {
990  Rma = Rj;
991  } else if (Rj < Rmi) {
992  Rmi = Rj;
993  } else {
994  continue;
995  }
996  // add 200 to quality if 2 compatible seeds with high deltaR of their
997  // spT have been found
998  // (i.e. potential track spans 5 surfaces)
999  if ((Rma - Rmi) > 20.) {
1000  u -= 200.;
1001  break;
1002  }
1003  } else {
1004  // first found compatible seed: add 200 to quality
1005  in = true;
1006  Rma = Rmi = Rj;
1007  u -= 200.;
1008  }
1009  }
1010  // if quality is below threshold, discard seed
1011  if (u > m_umax) {
1012  continue;
1013  }
1014  // if mixed seed and no compatible seed was found, discard seed
1015  if (pixb != pixt) {
1016  if (u > 0. || (u > ub && u > u0 && u > (*i).second->quality())) {
1017  continue;
1018  }
1019  }
1020  // if sct seed and at least 1 comp seed was found, accept even seeds with
1021  // larger impact params + cot theta penalty
1022  // m_diversss < Impact parameters + penalty for cot(theta) difference <
1023  // m_divermax
1024  // (if m_diversss set smaller than m_divermax, else m_diversss=m_divermax)
1025  if (!pixb && Im > m_diversss && u > Im - 500.) {
1026  continue;
1027  }
1028 
1029  newOneSeed(SPb, SP0, (*i).second, Zob, u);
1030  }
1031  m_CmSp.clear();
1032 }
1033 
1035 // Fill seeds
1037 template <class SpacePoint>
1039  m_fillOneSeeds = 0;
1040 
1041  typename std::multimap<float, Acts::Legacy::InternalSeed<SpacePoint>
1042  *>::iterator lf = m_mapOneSeeds.begin(),
1043  l = m_mapOneSeeds.begin(),
1044  le = m_mapOneSeeds.end();
1045 
1046  if (l == le) {
1047  return;
1048  }
1049 
1051 
1052  for (; l != le; ++l) {
1053  float w = (*l).first;
1054  s = (*l).second;
1055  if (l != lf && s->spacepoint0()->radius() < 43. && w > -200.) {
1056  continue;
1057  }
1058  if (!s->setQuality(w)) {
1059  continue;
1060  }
1061 
1062  if (i_seede != l_seeds.end()) {
1063  s = (*i_seede++);
1064  *s = *(*l).second;
1065  } else {
1066  s = new Acts::Legacy::InternalSeed<SpacePoint>(*(*l).second);
1067  l_seeds.push_back(s);
1068  i_seede = l_seeds.end();
1069  }
1070 
1071  if (s->spacepoint0()->spacepoint->clusterList().second) {
1072  w -= 3000.;
1073  } else if (s->spacepoint1()->spacepoint->clusterList().second) {
1074  w -= 2000.;
1075  } else if (s->spacepoint2()->spacepoint->clusterList().second) {
1076  w -= 1000.;
1077  }
1078 
1079  m_seeds.insert(std::make_pair(w, s));
1080  ++m_fillOneSeeds;
1081  }
1082 }