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