ATLAS Offline Software
TrigTrackSeedGenerator.cxx
Go to the documentation of this file.
1 /*
2  Copyright (C) 2002-2023 CERN for the benefit of the ATLAS collaboration
3 */
4 
5 #include <cmath>
6 #include <iostream>
7 #include <list>
8 #include <algorithm>
9 #include <iterator>
10 #include "TrigInDetEvent/TrigSiSpacePointBase.h"
11 #include "TrigInDetPattRecoEvent/TrigInDetTriplet.h"
12 #include "TrigInDetPattRecoEvent/TrigInDetSiLayer.h"
13 #include "TrigInDetPattRecoTools/TrigTrackSeedGenerator.h"
14 #include "IRegionSelector/IRoiDescriptor.h"
15 #include "IRegionSelector/RoiUtil.h"
16 #include "InDetPrepRawData/PixelCluster.h"
17 
18 TrigTrackSeedGenerator::TrigTrackSeedGenerator(const TrigCombinatorialSettings& tcs)
19  : m_settings(tcs),
20  m_minDeltaRadius(10.0),
21  m_zTol(3.0),
22  m_pStore(NULL)
23 {
24  m_maxDeltaRadius = m_settings.m_doublet_dR_Max;
25  m_maxDeltaRadiusConf = m_settings.m_doublet_dR_Max_Confirm;
26  m_phiSliceWidth = 2*M_PI/m_settings.m_nMaxPhiSlice;
27  m_pStore = new L_PHI_STORAGE(m_settings.m_nMaxPhiSlice, (int)m_settings.m_layerGeometry.size());
28 
29  //mult scatt. variance for doublet matching
30  const double radLen = 0.036;
31  m_CovMS = std::pow((13.6/m_settings.m_tripletPtMin),2)*radLen;
32  const double ptCoeff = 0.29997*1.9972/2.0;// ~0.3*B/2 - assumes nominal field of 2*T
33  m_minR_squ = m_settings.m_tripletPtMin*m_settings.m_tripletPtMin/std::pow(ptCoeff,2);
34  m_dtPreCut = std::tan(2.0*m_settings.m_tripletDtCut*std::sqrt(m_CovMS));
35 }
36 
37 TrigTrackSeedGenerator::~TrigTrackSeedGenerator() {
38  delete m_pStore;
39 
40  m_SoA.clear();
41 
42 }
43 
44 void TrigTrackSeedGenerator::loadSpacePoints(const std::vector<TrigSiSpacePointBase>& vSP) {
45 
46  m_pStore->reset();
47 
48  m_SoA.clear();
49 
50  m_spStorage.clear();
51 
52  m_spStorage.resize(vSP.size());
53 
54  m_minTau.resize(vSP.size(), 0.0);
55  m_maxTau.resize(vSP.size(), 100.0);
56 
57  int spIndex = 0;
58 
59  for(std::vector<TrigSiSpacePointBase>::const_iterator it = vSP.begin();it != vSP.end();++it) {
60  int layerId = (*it).layer();
61 
62  bool isPixel = (m_settings.m_layerGeometry[layerId].m_subdet == 1);
63  bool isEndcap = (m_settings.m_layerGeometry[layerId].m_type != 0);
64 
65  bool updateTauRange = false;
66  float minTau = 0.0;
67  float maxTau = 100.0;
68 
69  if((m_settings.m_useTrigSeedML > 0) && isPixel) {
70  const Trk::SpacePoint* osp = (*it).offlineSpacePoint();
71  const InDet::PixelCluster* pCL = dynamic_cast<const InDet::PixelCluster*>(osp->clusterList().first);
72 
73  float cluster_width = pCL->width().widthPhiRZ().y();
74  if(isEndcap) {
75  if(cluster_width > m_settings.m_maxEC_len) continue;
76  }
77  else {//Barrel
78  if(!m_settings.m_vLUT.empty()) {
79  const TrigSeedML_LUT& LUT = (*m_settings.m_vLUT.begin());
80  updateTauRange = LUT.getValidRange(cluster_width, minTau, maxTau);
81  }
82  }
83  }
84 
85  int phiIdx = ((*it).phi()+M_PI)/m_phiSliceWidth;
86  if (phiIdx >= m_settings.m_nMaxPhiSlice) {
87  phiIdx %= m_settings.m_nMaxPhiSlice;
88  }
89  else if (phiIdx < 0) {
90  phiIdx += m_settings.m_nMaxPhiSlice;
91  phiIdx %= m_settings.m_nMaxPhiSlice;
92  }
93  m_spStorage[spIndex].set(&(*it), spIndex);
94  m_pStore->addSpacePoint(phiIdx, layerId, &m_spStorage[spIndex]);
95 
96  if(updateTauRange) {
97  m_minTau[spIndex] = minTau;
98  m_maxTau[spIndex] = maxTau;
99  }
100  spIndex++;
101 
102  }
103 
104  m_pStore->sortSpacePoints2(m_settings.m_layerGeometry);
105 
106  m_SoA.resize(vSP.size());
107 
108 
109 }
110 
111 void TrigTrackSeedGenerator::createSeeds(const IRoiDescriptor* roiDescriptor) {
112 
113  m_zMinus = roiDescriptor->zedMinus() - m_zTol;
114  m_zPlus = roiDescriptor->zedPlus() + m_zTol;
115  m_zMinusEndcap = m_zMinus;
116  m_zPlusEndcap = m_zPlus;
117 
118  m_triplets.clear();
119 
120  int nLayers = (int) m_settings.m_layerGeometry.size();
121 
122 
123  //TO-DO: create vector of MIDDLE layers and for each middle layer create a list of
124  //ALLOWED layers for doublet making
125  //for example coaxial barrel layers should be separated by more than DR_Max
126  // no simultaneous +/- or -/+ endcaps
127 
128 
129  for(int layerI=1;layerI<nLayers;layerI++) {//skip layer 0 for middle spacepoint search
130 
131  const L_PHI_SECTOR& S = m_pStore->m_layers[layerI];
132  if(S.m_nSP==0) continue;
133 
134  bool isSct = (m_settings.m_layerGeometry[layerI].m_subdet == 2);
135 
136  for(int phiI=0;phiI<m_settings.m_nMaxPhiSlice;phiI++) {
137 
138  for(auto spm : S.m_phiSlices[phiI]) {//loop over middle spacepoints
139 
140  float zm = spm->m_pSP->z();
141  float rm = spm->m_pSP->r();
142 
143  m_nInner = 0;
144  m_nOuter = 0;
145 
146  m_innerMarkers.clear();
147  m_outerMarkers.clear();
148 
149  for(int layerJ=0;layerJ<nLayers;layerJ++) {
150 
151  bool isPixel2 = (m_settings.m_layerGeometry[layerJ].m_subdet == 1);
152 
153  bool checkPSS = (!m_settings.m_tripletDoPSS) && (isSct && isPixel2);
154  if(checkPSS) continue;//no mixed PSS seeds allowed
155 
156  if((!isSct) && (!isPixel2)) {// i.e. xPS (or SPx)
157  if((!m_settings.m_tripletDoConfirm) && (!m_settings.m_tripletDoPPS)) {
158  continue;//no mixed PPS seeds allowed
159  }
160  //but if m_tripletDoConfirm is true we will use PPS seeds to confirm PPP seeds
161  }
162 
163  if(!validateLayerPairNew(layerI, layerJ, rm, zm)) continue;
164 
165  const std::vector<const INDEXED_SP*>& spVec = m_pStore->m_layers[layerJ].m_phiThreeSlices.at(phiI);
166 
167  if(spVec.empty()) continue;
168 
169  SP_RANGE delta;
170 
171  if(!getSpacepointRange(layerJ, spVec, delta)) continue;
172 
173  int nI = m_nInner;
174  int nO = m_nOuter;
175 
176  processSpacepointRange(layerJ, spm, checkPSS, delta, roiDescriptor);
177 
178  if(m_nInner > nI) m_innerMarkers.push_back(m_nInner);
179  if(m_nOuter > nO) m_outerMarkers.push_back(m_nOuter);
180 
181 
182  }//loop over inner/outer layers
183 
184  if(m_nInner != 0 && m_nOuter != 0) {
185  std::vector<TrigInDetTriplet> tripletVec;
186 
187  if(m_settings.m_tripletDoConfirm && !isSct) {
188  createConfirmedTriplets(spm->m_pSP, m_nInner, m_nOuter, tripletVec, roiDescriptor);
189  }
190  else createTripletsNew(spm->m_pSP, m_nInner, m_nOuter, tripletVec, roiDescriptor);
191 
192  if(!tripletVec.empty()) storeTriplets(tripletVec);
193  tripletVec.clear();
194  }
195  }
196  }
197  }
198 
199 }
200 
201 void TrigTrackSeedGenerator::createSeeds(const IRoiDescriptor* roiDescriptor, const std::vector<float>& vZv) {
202 
203  m_triplets.clear();
204 
205  if(vZv.empty()) return;
206 
207  // With the vertex z information, there is no need to use the same z restriction as the ROI,
208  // but reasonable minimum and maximum bounds which should not be exceeded regardless of the
209  // z position of the vertex should still be set.
210  float roiZMinus = -225;
211  float roiZPlus = 225;
212 
213  int nLayers = (int) m_settings.m_layerGeometry.size();
214 
215 
216  for(int layerI=1;layerI<nLayers;layerI++) {//skip layer 0 for middle spacepoint search
217 
218  const L_PHI_SECTOR& S = m_pStore->m_layers[layerI];
219  if(S.m_nSP==0) continue;
220 
221  bool isSct = (m_settings.m_layerGeometry[layerI].m_subdet == 2);
222  bool isBarrel = (m_settings.m_layerGeometry[layerI].m_type == 0);
223 
224  bool checkWidth = isBarrel && (!isSct) && (m_settings.m_useTrigSeedML > 0);
225 
226  for(int phiI=0;phiI<m_settings.m_nMaxPhiSlice;phiI++) {
227 
228  for(auto spm : S.m_phiSlices[phiI]) {//loop over middle spacepoints
229 
230  float zm = spm->m_pSP->z();
231  float rm = spm->m_pSP->r();
232 
233  float minTau = 0.0;
234  float maxTau = 100.0;
235 
236  if(checkWidth) {
237  minTau = m_minTau[spm->m_idx];
238  maxTau = m_maxTau[spm->m_idx];
239  }
240 
241  std::vector<TrigInDetTriplet> tripletVec;
242 
243  for(const auto zVertex : vZv) {//loop over zvertices
244 
245  m_nInner = 0;
246  m_nOuter = 0;
247 
248  m_zMinus = std::max(zVertex - m_settings.m_zvError, roiZMinus);
249  m_zPlus = std::min(zVertex + m_settings.m_zvError, roiZPlus);
250  m_zPlusEndcap = std::min(zVertex + m_settings.m_zvErrorEndcap, roiZPlus);
251  m_zMinusEndcap = std::max(zVertex - m_settings.m_zvErrorEndcap, roiZMinus);
252 
253  for(int layerJ=0;layerJ<nLayers;layerJ++) {//loop over other layers
254 
255  bool isPixel2 = (m_settings.m_layerGeometry[layerJ].m_subdet == 1);
256 
257  if((!m_settings.m_tripletDoPSS) && (!m_settings.m_tripletDoPPS)) {//no mixed seeds allowed
258  if(isSct && isPixel2) continue;//no PSx
259  if((!isSct) && (!isPixel2)) continue;//no xPS
260  }
261 
262  if(!validateLayerPairNew(layerI, layerJ, rm, zm)) continue;
263 
264  bool checkPSS = (!m_settings.m_tripletDoPSS) && (isSct && isPixel2);
265 
266  // Get spacepoints from current and adjacent phi slices
267  const std::vector<const INDEXED_SP*>& spVec = m_pStore->m_layers[layerJ].m_phiThreeSlices.at(phiI);
268 
269  if(spVec.empty()) continue;
270 
271  SP_RANGE delta;
272 
273  if(!getSpacepointRange(layerJ, spVec, delta)) continue;
274 
275  processSpacepointRangeZv(spm, checkPSS, delta, checkWidth, minTau, maxTau);
276 
277  }//loop over inner/outer layers
278  if(m_nInner != 0 && m_nOuter != 0) {
279  createTriplets(spm->m_pSP, m_nInner, m_nOuter, tripletVec, roiDescriptor);
280  }
281  }//loop over zvertices
282  if(!tripletVec.empty()) storeTriplets(tripletVec);
283  tripletVec.clear();
284  }
285  }
286  }
287 
288 
289 }
290 
291 bool TrigTrackSeedGenerator::validateLayerPairNew(int layerI, int layerJ, float rm, float zm) {
292 
293  const float deltaRefCoord = 5.0;
294 
295  if(layerJ==layerI) return false;//skip the same layer ???
296 
297  if(m_pStore->m_layers[layerJ].m_nSP==0) return false;
298 
299  int typeI = m_settings.m_layerGeometry[layerI].m_type;
300 
301  float refCoordI = m_settings.m_layerGeometry[layerI].m_refCoord;
302 
303  int typeJ = m_settings.m_layerGeometry[layerJ].m_type;
304 
305  float refCoordJ = m_settings.m_layerGeometry[layerJ].m_refCoord;
306 
307  if((typeI!=0) && (typeJ!=0) && refCoordI*refCoordJ<0.0) return false;//only ++ or -- EC combinations are allowed
308 
309  //project beamline interval to the ref. coord of the layer
310 
311  bool isBarrel = (typeJ == 0);
312 
313  if(isBarrel && std::fabs(refCoordJ-rm)>m_maxDeltaRadius) return false;
314 
315  //boundaries for nextLayer
316 
317  m_minCoord = 10000.0;
318  m_maxCoord =-10000.0;
319 
320  if(isBarrel) {
321 
322  if(refCoordJ<rm){//inner layer
323 
324  m_minCoord = m_zMinus + refCoordJ*(zm-m_zMinus)/rm;//+deltaRefCoord
325  m_maxCoord = m_zPlus + refCoordJ*(zm-m_zPlus)/rm;//-deltaRefCoord
326  m_minCoord -= deltaRefCoord*std::fabs(zm-m_zMinus)/rm;//corrections due to the layer width
327  m_maxCoord += deltaRefCoord*std::fabs(zm-m_zPlus)/rm;
328  }
329  else {//outer layer
330 
331  m_minCoord = m_zPlus + refCoordJ*(zm-m_zPlus)/rm;//+deltaRefCoord
332  m_maxCoord = m_zMinus + refCoordJ*(zm-m_zMinus)/rm;//-deltaRefCoord
333  m_minCoord -= deltaRefCoord*std::fabs(zm-m_zPlus)/rm;
334  m_maxCoord += deltaRefCoord*std::fabs(zm-m_zMinus)/rm;
335  }
336  }
337  else {//endcap
338  float maxB =m_settings.m_layerGeometry[layerJ].m_maxBound;
339  float minB =m_settings.m_layerGeometry[layerJ].m_minBound;
340  if(maxB<=rm) return false;// This currently rejects SP type doublets
341  // Could correct this by retrieving if layers are pix or SCT, and not performing this check for SCt->pix doublets
342  if(refCoordJ>0) {//positive EC
343  if(refCoordJ > zm) {//outer layer
344 
345  if(zm < m_zMinusEndcap) return false;
346  if(zm == m_zPlusEndcap) return false;
347  float zMax = (zm*maxB-rm*refCoordJ)/(maxB-rm);
348  if( m_zMinusEndcap > zMax) return false;
349  if (rm < minB) {
350  float zMin = (zm*minB-rm*refCoordJ)/(minB-rm);
351  if(m_zPlusEndcap<zMin) return false;
352  }
353 
354  m_minCoord = (refCoordJ-m_zMinusEndcap)*rm/(zm-m_zMinusEndcap);
355  m_maxCoord = (refCoordJ-m_zPlusEndcap)*rm/(zm-m_zPlusEndcap);
356  m_minCoord -= deltaRefCoord*rm/std::fabs(zm-m_zMinusEndcap);
357  m_maxCoord += deltaRefCoord*rm/std::fabs(zm-m_zPlusEndcap);
358 
359  if(zm <= m_zPlusEndcap) m_maxCoord = maxB;
360 
361  if(m_minCoord > maxB) return false;
362  if(m_maxCoord < minB) return false;
363 
364  }
365  else {//inner layer
366  if(minB == rm) return false;
367  float zMax = (zm*minB-rm*refCoordJ)/(minB-rm);
368  if( m_zMinusEndcap > zMax) return false;
369  if (rm>maxB) {// otherwise, intersect of line from maxB through middle sp will be on the wrong side of the layer
370  float zMin = (zm*maxB-rm*refCoordJ)/(maxB-rm);
371  if(m_zPlusEndcap<zMin) return false;
372  }
373  if(zm == m_zPlusEndcap || zm == m_zMinusEndcap) return false;
374  m_minCoord = (refCoordJ-m_zPlusEndcap)*rm/(zm-m_zPlusEndcap);
375  m_maxCoord = (refCoordJ-m_zMinusEndcap)*rm/(zm-m_zMinusEndcap);
376  m_minCoord -= deltaRefCoord*rm/std::fabs(zm-m_zPlusEndcap);
377  m_maxCoord += deltaRefCoord*rm/std::fabs(zm-m_zMinusEndcap);
378  }
379  }
380  else {//negative EC
381  if(refCoordJ < zm) {//outer layer
382 
383  if(zm > m_zPlusEndcap) return false;
384  if(zm == m_zMinusEndcap) return false;
385  float zMin = (zm*maxB-rm*refCoordJ)/(maxB-rm);
386  if( m_zPlusEndcap < zMin) return false;
387  if (rm<minB) {// otherwise, intersect of line from minB through middle sp will be on the wrong side of the layer
388  float zMax = (zm*minB-rm*refCoordJ)/(minB-rm);
389  if(m_zMinusEndcap>zMax) return false;
390  }
391 
392  m_minCoord = (refCoordJ-m_zPlusEndcap)*rm/(zm-m_zPlusEndcap);
393  m_maxCoord = (refCoordJ-m_zMinusEndcap)*rm/(zm-m_zMinusEndcap);
394  m_minCoord -= deltaRefCoord*rm/std::fabs(zm-m_zPlusEndcap);
395  m_maxCoord += deltaRefCoord*rm/std::fabs(zm-m_zMinusEndcap);
396  if(zm > m_zMinusEndcap) m_maxCoord = maxB;
397  if(m_minCoord > maxB) return false;
398  if(m_maxCoord < minB) return false;
399 
400  }
401  else {//inner layer
402  if(minB == rm) return false;
403  float zMin = (zm*minB-rm*refCoordJ)/(minB-rm);
404  if( m_zPlusEndcap < zMin) return false;
405  if (rm>maxB) {// otherwise, intersect of line from maxB through middle sp will be on the wrong side of the layer
406  float zMax = (zm*maxB-rm*refCoordJ)/(maxB-rm);
407  if(m_zMinusEndcap>zMax) return false;
408  }
409  if(zm == m_zPlusEndcap || zm == m_zMinusEndcap) return false;
410  m_minCoord = (refCoordJ-m_zMinusEndcap)*rm/(zm-m_zMinusEndcap);
411  m_maxCoord = (refCoordJ-m_zPlusEndcap)*rm/(zm-m_zPlusEndcap);
412  m_minCoord -= deltaRefCoord*rm/std::fabs(zm-m_zMinusEndcap);
413  m_maxCoord += deltaRefCoord*rm/std::fabs(zm-m_zPlusEndcap);
414 
415  }
416  }
417  }
418  float minBoundJ = m_settings.m_layerGeometry[layerJ].m_minBound;
419  float maxBoundJ = m_settings.m_layerGeometry[layerJ].m_maxBound;
420  if(maxBoundJ<m_minCoord || minBoundJ>m_maxCoord) return false;
421  return true;
422 }
423 
424 
425 
426 bool TrigTrackSeedGenerator::getSpacepointRange(int lJ, const std::vector<const INDEXED_SP*>& spVec, SP_RANGE& delta) {
427 
428  int typeJ = m_settings.m_layerGeometry[lJ].m_type;
429  bool isBarrel = (typeJ == 0);
430  bool isPositiveEC = m_settings.m_layerGeometry[lJ].m_refCoord > 0;
431  float minSpCoord = (isBarrel) ? (*spVec.begin())->m_pSP->z() : (*spVec.begin())->m_pSP->r();
432  float maxSpCoord = (isBarrel) ? (*spVec.rbegin())->m_pSP->z() : (*spVec.rbegin())->m_pSP->r();
433 
434  if(!isBarrel && isPositiveEC) {//forward endcap SPs are sorted differently
435  float tmp = minSpCoord;minSpCoord = maxSpCoord;maxSpCoord = tmp;
436  }
437 
438  if(minSpCoord > m_maxCoord) return false;
439  if(maxSpCoord < m_minCoord) return false;
440 
441  //identify the range of spacepoints in the layer
442 
443  std::vector<const INDEXED_SP*>::const_iterator it1 = spVec.end();
444  std::vector<const INDEXED_SP*>::const_iterator it2 = spVec.end();
445 
446  if(isBarrel) {
447  it1 = std::lower_bound(spVec.begin(), spVec.end(), m_minCoord, L_PHI_SECTOR::smallerThanZ());
448  it2 = std::upper_bound(spVec.begin(), spVec.end(), m_maxCoord, L_PHI_SECTOR::greaterThanZ());
449  }
450  else {
451  if(isPositiveEC) {
452  it1 = std::lower_bound(spVec.begin(), spVec.end(), m_maxCoord, L_PHI_SECTOR::greaterThanR_i());
453  it2 = std::upper_bound(spVec.begin(), spVec.end(), m_minCoord, L_PHI_SECTOR::smallerThanR_i());
454  }
455  else {
456  it1 = std::lower_bound(spVec.begin(), spVec.end(), m_minCoord, L_PHI_SECTOR::smallerThanR());
457  it2 = std::upper_bound(spVec.begin(), spVec.end(), m_maxCoord, L_PHI_SECTOR::greaterThanR());
458  }
459  }
460 
461  if(std::distance(it1, it2)==0) return false;
462 
463  delta.first = it1;delta.second = it2;
464  return true;
465 }
466 
467 int TrigTrackSeedGenerator::processSpacepointRange(int lJ, const INDEXED_SP* spm, bool checkPSS, const SP_RANGE& delta, const IRoiDescriptor* roiDescriptor) {
468 
469  int nSP=0;
470 
471  bool isBarrel = (m_settings.m_layerGeometry[lJ].m_type==0);
472  float refCoord = m_settings.m_layerGeometry[lJ].m_refCoord;
473  bool isPixel = spm->m_pSP->isPixel();
474 
475  float rm = spm->m_pSP->r();
476  float zm = spm->m_pSP->z();
477 
478  float dZ = refCoord-zm;
479  float dR_i = isBarrel ? 1.0/(refCoord-rm) : 1.0;
480 
481  int SeedML = m_settings.m_useTrigSeedML;
482 
483  for(std::vector<const INDEXED_SP*>::const_iterator spIt=delta.first; spIt!=delta.second; ++spIt) {
484  // if(deIds == (*spIt)->offlineSpacePoint()->elementIdList()) continue;
485  float zsp = (*spIt)->m_pSP->z();
486  float rsp = (*spIt)->m_pSP->r();
487 
488  float dr = rsp - rm;
489 
490  if(std::fabs(dr)>m_maxDeltaRadius ) continue;
491 
492  if(m_settings.m_tripletDoConfirm) {
493  if(isPixel && !(*spIt)->m_pSP->isPixel()) {
494  if(std::fabs(dr)>m_maxDeltaRadiusConf ) continue;
495  }
496  }
497 
498  if(std::fabs(dr)<m_minDeltaRadius ) continue;
499 
500  //inner doublet check
501 
502  if (dr<0 && checkPSS) continue;//Allow PSS seeds ?
503 
504  float dz = zsp - zm;
505  float tau = dz/dr;
506  float ftau = std::fabs(tau);
507  if (ftau > 7.41) continue;
508 
509  if(isPixel && SeedML > 0) {
510  if(ftau < m_minTau[spm->m_idx] || ftau > m_maxTau[spm->m_idx]) continue;
511  }
512 
513 
514  float z0 = zm - rm*tau;
515  if (m_settings.m_doubletFilterRZ) {
516  if (!RoiUtil::contains( *roiDescriptor, z0, tau)) continue;
517  }
518 
519  float t = isBarrel ? dz*dR_i : dZ/dr;
520 
521  if(dr<0) {
522  if(isBarrel && (*spIt)->m_pSP->isPixel()) {
523  if(SeedML == 3 || SeedML == 4) {
524  if(ftau < m_minTau[(*spIt)->m_idx] || ftau > m_maxTau[(*spIt)->m_idx]) continue;
525  }
526  }
527  m_SoA.m_spi[m_nInner] = (*spIt)->m_pSP;
528  m_SoA.m_ti[m_nInner] = t;
529  m_nInner++;
530  }
531  else {
532  if(isBarrel && (*spIt)->m_pSP->isPixel()) {
533  if(SeedML == 2 || SeedML == 4) {
534  if(ftau < m_minTau[(*spIt)->m_idx] || ftau > m_maxTau[(*spIt)->m_idx]) continue;
535  }
536  }
537  m_SoA.m_spo[m_nOuter] = (*spIt)->m_pSP;
538  m_SoA.m_to[m_nOuter] = t;
539  m_nOuter++;
540  }
541  nSP++;
542  }
543 
544  return nSP;
545 }
546 
547 int TrigTrackSeedGenerator::processSpacepointRangeZv(const INDEXED_SP* spm, bool checkPSS, const SP_RANGE& delta, bool checkWidth, const float& minTau, const float& maxTau) {
548 
549  int nSP=0;
550 
551  float rm = spm->m_pSP->r();
552  float zm = spm->m_pSP->z();
553 
554  int SeedML = m_settings.m_useTrigSeedML;
555 
556  for(std::vector<const INDEXED_SP*>::const_iterator spIt=delta.first; spIt!=delta.second; ++spIt) {
557 
558  float rsp = (*spIt)->m_pSP->r();
559  float zsp = (*spIt)->m_pSP->z();
560  float dr = rsp - rm;
561 
562  if(std::fabs(dr)>m_maxDeltaRadius ) continue;
563  if(std::fabs(dr)<m_minDeltaRadius ) continue;
564 
565  //inner doublet check
566 
567  if (dr<0 && checkPSS) continue;//Allow PSS seeds ?
568 
569  float dz = zsp - zm;
570  float tau = dz/dr;
571  float ftau = std::fabs(tau);
572 
573  if (ftau > 7.41) continue;//|eta|<2.7
574 
575  if(checkWidth) {
576  if(ftau < minTau || ftau > maxTau) continue;
577  }
578 
579  if(dr<0) {
580  if((*spIt)->m_pSP->isPixel()) {
581  if(SeedML == 3 || SeedML == 4) {
582  if(ftau < m_minTau[(*spIt)->m_idx] || ftau > m_maxTau[(*spIt)->m_idx]) continue;
583  }
584  }
585  m_SoA.m_spi[m_nInner++] = (*spIt)->m_pSP;
586  }
587  else {
588  if((*spIt)->m_pSP->isPixel()) {
589  if(SeedML == 2 || SeedML == 4) {
590  if(ftau < m_minTau[(*spIt)->m_idx] || ftau > m_maxTau[(*spIt)->m_idx]) continue;
591  }
592  }
593  m_SoA.m_spo[m_nOuter++] = (*spIt)->m_pSP;
594  }
595 
596  nSP++;
597  }
598 
599  return nSP;
600 }
601 
602 void TrigTrackSeedGenerator::createTriplets(const TrigSiSpacePointBase* pS, int nInner, int nOuter,
603  std::vector<TrigInDetTriplet>& output, const IRoiDescriptor* roiDescriptor) {
604 
605  if(nInner==0 || nOuter==0) return;
606 
611  bool fullPhi = ( roiDescriptor->isFullscan() || ( std::fabs( roiDescriptor->phiPlus() - roiDescriptor->phiMinus() - 2*M_PI ) < 1e-6 ) );
612 
613  int nSP = nInner + nOuter;
614 
615  const double pS_r = pS->r();
616  const double pS_x = pS->x();
617  const double pS_y = pS->y();
618  const double pS_dr = pS->dr();
619  const double pS_dz = pS->dz();
620  const double cosA = pS_x/pS_r;
621  const double sinA = pS_y/pS_r;
622  const double covZ = pS_dz*pS_dz;
623  const double covR = pS_dr*pS_dr;
624  const bool isPixel = pS->isPixel();
625 
626  int idx = 0;
627  for(;idx<nInner;idx++) {
628  const TrigSiSpacePointBase* pSP = m_SoA.m_spi[idx];
629 
630  //1. transform in the pS-centric c.s
631 
632  const double dx = pSP->x() - pS_x;
633  const double dy = pSP->y() - pS_y;
634  const double R2inv = 1.0/(dx*dx+dy*dy);
635  const double Rinv = std::sqrt(R2inv);
636  const double xn = dx*cosA + dy*sinA;
637  const double yn =-dx*sinA + dy*cosA;
638  const double dz = pSP->z() - pS->z();
639  const double t = Rinv*dz;
640 
641  //2. Conformal mapping
642 
643  m_SoA.m_r[idx] = Rinv;
644  m_SoA.m_u[idx] = xn*R2inv;
645  m_SoA.m_v[idx] = yn*R2inv;
646 
647  //3. cot(theta) estimation for the doublet
648 
649  const double covZP = pSP->dz()*pSP->dz();
650  const double covRP = pSP->dr()*pSP->dr();
651 
652  m_SoA.m_t[idx] = t;
653  m_SoA.m_tCov[idx] = R2inv*(covZ + covZP + t*t*(covR+covRP));
654  }
655 
656  for(int k=0;k<nOuter;k++,idx++) {
657  const TrigSiSpacePointBase* pSP = m_SoA.m_spo[k];
658 
659  //1. transform in the pS-centric c.s
660 
661  const double dx = pSP->x() - pS_x;
662  const double dy = pSP->y() - pS_y;
663  const double R2inv = 1.0/(dx*dx+dy*dy);
664  const double Rinv = std::sqrt(R2inv);
665  const double xn = dx*cosA + dy*sinA;
666  const double yn =-dx*sinA + dy*cosA;
667  const double dz = -pSP->z() + pS->z();
668  const double t = Rinv*dz;
669 
670  //2. Conformal mapping
671 
672  m_SoA.m_r[idx] = Rinv;
673  m_SoA.m_u[idx] = xn*R2inv;
674  m_SoA.m_v[idx] = yn*R2inv;
675 
676  //3. cot(theta) estimation for the doublet
677 
678  const double covZP = pSP->dz()*pSP->dz();
679  const double covRP = pSP->dr()*pSP->dr();
680 
681  m_SoA.m_t[idx] = t;
682  m_SoA.m_tCov[idx] = R2inv*(covZ + covZP + t*t*(covR+covRP));
683  }
684 
685  //double loop
686 
687  for(int innIdx=0;innIdx<nInner;innIdx++) {
688 
689  //mult. scatt contribution due to the layer with middle SP
690 
691  const double r_inn = m_SoA.m_r[innIdx];
692  const double t_inn = m_SoA.m_t[innIdx];
693  const double v_inn = m_SoA.m_v[innIdx];
694  const double u_inn = m_SoA.m_u[innIdx];
695  const double tCov_inn = m_SoA.m_tCov[innIdx];
696  const double dCov = m_CovMS*(1+t_inn*t_inn);
697 
698 
699  for(int outIdx=nInner;outIdx<nSP;outIdx++) {
700 
701  //1. rz doublet matching
702  const double t_out = m_SoA.m_t[outIdx];
703 
704 
705 
706  const double dt2 = std::pow((t_inn - t_out), 2)*(1.0/9.0);
707 
708 
709  double covdt = (t_inn*t_out*covR + covZ);
710  covdt *= 2*r_inn*m_SoA.m_r[outIdx];
711  covdt += tCov_inn + m_SoA.m_tCov[outIdx];
712 
713  if(dt2 > covdt+dCov) continue;
714 
715  //2. pT estimate
716 
717  const double du = m_SoA.m_u[outIdx] - u_inn;
718  if(du==0.0) continue;
719  const double A = (m_SoA.m_v[outIdx] - v_inn)/du;
720  const double B = v_inn - A*u_inn;
721  const double R_squ = (1 + A*A)/(B*B);
722 
723  if(R_squ < m_minR_squ) continue;
724 
725  //3. the 3-sigma cut with estimated pT
726 
727  const double frac = m_minR_squ/R_squ;
728  if(dt2 > covdt+frac*dCov) continue;
729 
730  //4. d0 cut
731 
732  const double fabs_d0 = std::fabs(pS_r*(B*pS_r - A));
733 
734  if(fabs_d0 > m_settings.m_tripletD0Max) continue;
735 
736  if (m_SoA.m_spo[outIdx-nInner]->isSCT() && isPixel) {
737  if(fabs_d0 > m_settings.m_tripletD0_PPS_Max) continue;
738  }
739 
740  //5. phi0 cut
741 
742  if ( !fullPhi ) {
743  const double uc = 2*B*pS_r - A;
744  const double phi0 = atan2(sinA - uc*cosA, cosA + uc*sinA);
745 
746  if ( !RoiUtil::containsPhi( *roiDescriptor, phi0 ) ) {
747  continue;
748  }
749  }
750 
751  //6. add new triplet
752 
753  // Calculate triplet weight: No weighting in LRT mode, but d0**2 weighting for normal (non LRT) mode. Triplets are then sorted by lowest weight.
754  const double Q= (m_settings.m_LRTmode ? 0 : fabs_d0*fabs_d0);
755  if(output.size()>=m_settings.m_maxTripletBufferLength) {
756  if (m_settings.m_LRTmode) { // take the first m_maxTripletBufferLength triplets
757  continue;
758  } else { // choose smallest d0
759 
760  std::sort(output.begin(), output.end(),
761  [](const TrigInDetTriplet& A, const TrigInDetTriplet& B) {
762  return A.Q() > B.Q();
763  }
764  );
765 
766  std::vector<TrigInDetTriplet>::iterator it = output.begin();
767  if( Q >= (*it).Q()) {
768  continue;
769  }
770  output.erase(it);
771  }
772  }
773  output.emplace_back(*m_SoA.m_spi[innIdx], *pS, *m_SoA.m_spo[outIdx-nInner], Q);
774  }
775  }
776 }
777 
778 void TrigTrackSeedGenerator::createTripletsNew(const TrigSiSpacePointBase* pS, int nInner, int nOuter,
779  std::vector<TrigInDetTriplet>& output, const IRoiDescriptor* roiDescriptor) {
780 
781  if(nInner==0 || nOuter==0) return;
782 
784  bool fullPhi = ( roiDescriptor->isFullscan() || ( std::fabs( roiDescriptor->phiPlus() - roiDescriptor->phiMinus() - 2*M_PI ) < 1e-6 ) );
785 
786 
787  int nSP = nInner + nOuter;
788  output.reserve(m_settings.m_maxTripletBufferLength);
789 
790  int nL = m_outerMarkers.size() + m_innerMarkers.size();
791 
792  // There are 32 pixel + SCT layers altogether, and
793  // pairs of spacepoints from the same layer are not considered,
794  // so these arrays should need to contain at most 31 elements
795  int nleft[31];
796  int iter[31];
797  int dirs[31];
798  int type[31];
799  double values[31];
800 
801  iter[0] = m_innerMarkers[0]-1;
802  nleft[0] = m_innerMarkers[0];
803  dirs[0] = -1;
804  type[0] = 0;//inner sp
805  unsigned int kL=1;
806  for(;kL<m_innerMarkers.size();kL++) {
807  iter[kL] = m_innerMarkers[kL]-1;
808  nleft[kL] = m_innerMarkers[kL]-m_innerMarkers[kL-1];
809  dirs[kL] = -1;
810  type[kL] = 0;
811  }
812  iter[kL] = 0;
813  nleft[kL] = m_outerMarkers[0];
814  dirs[kL] = 1;
815  type[kL] = 1;//outer sp
816  kL++;
817  for(unsigned int k=1;k<m_outerMarkers.size();k++) {
818  iter[kL] = m_outerMarkers[k-1];
819  nleft[kL] = m_outerMarkers[k] - m_outerMarkers[k-1];
820  dirs[kL] = 1;
821  type[kL] = 1;
822  kL++;
823  }
824 
825  for(int k=0;k<nL;k++) {
826  values[k] = (type[k]==0) ? m_SoA.m_ti[iter[k]] : m_SoA.m_to[iter[k]];//initialization
827  }
828 
829  //merged sort
830 
831  int nActive=nL;
832  int iSP=0;
833  while(nActive!=0) {
834  nActive = 0;
835  //find min element
836  double min_el=1000.0;
837  int k_min=-1;
838  for(int k=0;k<nL;k++) {
839  if(nleft[k]==0) continue;
840  nActive++;
841  if(values[k]<min_el) {
842  min_el = values[k];
843  k_min = k;
844  }
845  }
846  if(k_min==-1) break;
847  //output
848 
849  int i_min = iter[k_min];
850 
851  if(type[k_min]==0){//inner
852  m_SoA.m_sorted_sp[iSP] = m_SoA.m_spi[i_min];
853  } else {//outer
854  m_SoA.m_sorted_sp[iSP] = m_SoA.m_spo[i_min];
855  }
856 
857  m_SoA.m_sorted_sp_type[iSP] = type[k_min];
858  m_SoA.m_sorted_sp_t[iSP] = min_el;
859  iSP++;
860 
861  iter[k_min] += dirs[k_min];
862  nleft[k_min]--;
863  if(nleft[k_min]==0) {
864  nActive--;
865  continue;
866  }
867  values[k_min] = (type[k_min]==0) ? m_SoA.m_ti[iter[k_min]] : m_SoA.m_to[iter[k_min]];
868  }
869 
870  const double pS_r = pS->r();
871  const double pS_x = pS->x();
872  const double pS_y = pS->y();
873  //const double pS_z = pS->z();
874  const double pS_dr = pS->dr();
875  const double pS_dz = pS->dz();
876  const double cosA = pS_x/pS_r;
877  const double sinA = pS_y/pS_r;
878  const double covZ = pS_dz*pS_dz;
879  const double covR = pS_dr*pS_dr;
880  const bool isPixel = pS->isPixel();
881 
882  for(int idx=0;idx<nSP;idx++) {
883 
884  const TrigSiSpacePointBase* pSP = m_SoA.m_sorted_sp[idx];
885 
886  //1. transform in the pS-centric c.s
887 
888  const double dx = pSP->x() - pS_x;
889  const double dy = pSP->y() - pS_y;
890  const double R2inv = 1.0/(dx*dx+dy*dy);
891  const double Rinv = std::sqrt(R2inv);
892  const double xn = dx*cosA + dy*sinA;
893  const double yn =-dx*sinA + dy*cosA;
894  const double dz = (m_SoA.m_sorted_sp_type[idx]==0) ? pSP->z() - pS->z() : -pSP->z() + pS->z();
895 
896  const double t = Rinv*dz;
897 
898  //2. Conformal mapping
899 
900  m_SoA.m_r[idx] = Rinv;
901  m_SoA.m_u[idx] = xn*R2inv;
902  m_SoA.m_v[idx] = yn*R2inv;
903 
904  //3. cot(theta) estimation for the doublet
905 
906  const double covZP = pSP->dz()*pSP->dz();
907  const double covRP = pSP->dr()*pSP->dr();
908 
909  m_SoA.m_t[idx] = t;
910  m_SoA.m_tCov[idx] = R2inv*(covZ + covZP + t*t*(covR+covRP));
911  }
912 
913  //double loop
914 
915  int iter1 = 0;
916 
917  while(iter1<nSP-1) {//outer loop
918 
919  int type1 = m_SoA.m_sorted_sp_type[iter1];//0 - inner, 1 - outer
920  double t1 = m_SoA.m_sorted_sp_t[iter1];
921  double tcut = (1 + t1*t1)*m_dtPreCut;
922 
923  //mult. scatt contribution due to the layer with middle SP
924 
925  const double r_inn = m_SoA.m_r[iter1];
926  const double t_inn = m_SoA.m_t[iter1];
927  const double v_inn = m_SoA.m_v[iter1];
928  const double u_inn = m_SoA.m_u[iter1];
929  const double tCov_inn = m_SoA.m_tCov[iter1];
930  const double dCov = m_CovMS*(1+t_inn*t_inn);
931 
932  for(int iter2=iter1+1;iter2<nSP;iter2++) {//inner loop
933 
934  if(type1==m_SoA.m_sorted_sp_type[iter2]) continue;
935 
936  // 1. rz doublet matching
937  if(m_SoA.m_sorted_sp_t[iter2]-t1>tcut) break;
938 
939  const double t_out = m_SoA.m_t[iter2];
940 
941  const double dt2 = std::pow((t_inn - t_out), 2)*(1.0/9.0);
942 
943  double covdt = (t_inn*t_out*covR + covZ);
944  covdt *= 2*r_inn*m_SoA.m_r[iter2];
945  covdt += tCov_inn + m_SoA.m_tCov[iter2];
946 
947  if(dt2 > covdt+dCov) continue;
948 
949  //2. pT estimate
950 
951  const double du = m_SoA.m_u[iter2] - u_inn;
952  if(du==0.0) continue;
953 
954  const double A = (m_SoA.m_v[iter2] - v_inn)/du;
955  //Branchless version of (type1==0) ? v_inn - A*u_inn : m_SoA.m_v[iter2] - A*m_SoA.m_u[iter2];
956  const double B = (1-type1)*(v_inn - A*u_inn) + type1*(m_SoA.m_v[iter2] - A*m_SoA.m_u[iter2]);
957  const double R_squ = (1 + A*A)/(B*B);
958 
959  if(R_squ < m_minR_squ) continue;
960 
961  //3. the 3-sigma cut with estimated pT
962 
963  const double frac = m_minR_squ/R_squ;
964  if(dt2 > covdt+frac*dCov) continue;
965 
966  //4. d0 cut
967  const double B_pS_r = B*pS_r;//Pre-calculate for use in phi check
968  const double fabs_d0 = std::fabs(pS_r*(B_pS_r - A));
969 
970  if(fabs_d0 > m_settings.m_tripletD0Max) continue;
971 
972  //Branchless version of (type1 == 1) ? m_SoA.m_sorted_sp[iter1]->isSCT() : m_SoA.m_sorted_sp[iter2]->isSCT();
973  bool isSCT = (1-type1)*m_SoA.m_sorted_sp[iter1]->isSCT() + type1*m_SoA.m_sorted_sp[iter2]->isSCT();
974 
975  if (isSCT && isPixel) {
976  if(fabs_d0 > m_settings.m_tripletD0_PPS_Max) continue;
977  }
978 
979  //5. phi0 cut
980 
981  if ( !fullPhi ) {
982  //previously an incorrect calculation was used uc = 2*(B*pS_r - A); see ATR-28202 for details
983  double uc = 2*B_pS_r - A;
984 
985  const double phi0 = atan2(sinA - uc*cosA, cosA + uc*sinA);
986 
987  if ( !RoiUtil::containsPhi( *roiDescriptor, phi0 ) ) {
988  continue;
989  }
990  }
991 
992  //6. add new triplet
993 
994  const double Q = fabs_d0*fabs_d0;
995 
996  if(output.size()>=m_settings.m_maxTripletBufferLength) {
997  std::sort(output.begin(), output.end(),
998  [](const TrigInDetTriplet& A, const TrigInDetTriplet& B) {
999  return A.Q() > B.Q();
1000  }
1001  );
1002 
1003  std::vector<TrigInDetTriplet>::iterator it = output.begin();
1004  if( Q >= (*it).Q()) {
1005  continue;
1006  }
1007  output.erase(it);
1008 
1009  }
1010 
1011  const TrigSiSpacePointBase* pSPI = (type1==0) ? m_SoA.m_sorted_sp[iter1] : m_SoA.m_sorted_sp[iter2];
1012  const TrigSiSpacePointBase* pSPO = (type1==0) ? m_SoA.m_sorted_sp[iter2] : m_SoA.m_sorted_sp[iter1];
1013 
1014  output.emplace_back(*pSPI, *pS, *pSPO, Q);
1015  }
1016 
1017  iter1++;
1018  }
1019 
1020 }
1021 
1022 void TrigTrackSeedGenerator::createConfirmedTriplets(const TrigSiSpacePointBase* pS, int nInner, int nOuter,
1023  std::vector<TrigInDetTriplet>& output, const IRoiDescriptor* roiDescriptor) {
1024 
1025 
1026  struct ProtoSeed {
1027  public:
1028  ProtoSeed(const TrigSiSpacePointBase* s, double c, double Q, bool conf = false) : m_sp(s), m_curv(c), m_Q(Q), m_confirmed(conf) {};
1029  ProtoSeed(const ProtoSeed& ps) : m_sp(ps.m_sp), m_curv(ps.m_curv), m_Q(ps.m_Q), m_confirmed(ps.m_confirmed) {};
1030  const TrigSiSpacePointBase* m_sp;
1031  double m_curv, m_Q;
1032  bool m_confirmed;
1033  };
1034 
1035  if(nInner==0 || nOuter==0) return;
1036 
1038  bool fullPhi = ( roiDescriptor->isFullscan() || ( std::fabs( roiDescriptor->phiPlus() - roiDescriptor->phiMinus() - 2*M_PI ) < 1e-6 ) );
1039 
1040 
1041  int nSP = nInner + nOuter;
1042 
1043  const double pS_r = pS->r();
1044  const double pS_x = pS->x();
1045  const double pS_y = pS->y();
1046  const double pS_dr = pS->dr();
1047  const double pS_dz = pS->dz();
1048  const double cosA = pS_x/pS_r;
1049  const double sinA = pS_y/pS_r;
1050  const double covZ = pS_dz*pS_dz;
1051  const double covR = pS_dr*pS_dr;
1052  const bool isPixel = pS->isPixel();
1053 
1054 
1055  double bestQ = 1e8;
1056 
1057 
1058  int idx = 0;
1059  for(;idx<nInner;idx++) {
1060  const TrigSiSpacePointBase* pSP = m_SoA.m_spi[idx];
1061 
1062  //1. transform in the pS-centric c.s
1063 
1064  const double dx = pSP->x() - pS_x;
1065  const double dy = pSP->y() - pS_y;
1066  const double R2inv = 1.0/(dx*dx+dy*dy);
1067  const double Rinv = std::sqrt(R2inv);
1068  const double xn = dx*cosA + dy*sinA;
1069  const double yn =-dx*sinA + dy*cosA;
1070  const double dz = pSP->z() - pS->z();
1071  const double t = Rinv*dz;
1072 
1073  //2. Conformal mapping
1074 
1075  m_SoA.m_r[idx] = Rinv;
1076  m_SoA.m_u[idx] = xn*R2inv;
1077  m_SoA.m_v[idx] = yn*R2inv;
1078 
1079  //3. cot(theta) estimation for the doublet
1080 
1081  const double covZP = pSP->dz()*pSP->dz();
1082  const double covRP = pSP->dr()*pSP->dr();
1083 
1084  m_SoA.m_t[idx] = t;
1085  m_SoA.m_tCov[idx] = R2inv*(covZ + covZP + t*t*(covR+covRP));
1086  }
1087 
1088  for(int k=0;k<nOuter;k++,idx++) {
1089  const TrigSiSpacePointBase* pSP = m_SoA.m_spo[k];
1090 
1091  //1. transform in the pS-centric c.s
1092 
1093  const double dx = pSP->x() - pS_x;
1094  const double dy = pSP->y() - pS_y;
1095  const double R2inv = 1.0/(dx*dx+dy*dy);
1096  const double Rinv = std::sqrt(R2inv);
1097  const double xn = dx*cosA + dy*sinA;
1098  const double yn =-dx*sinA + dy*cosA;
1099  const double dz = -pSP->z() + pS->z();
1100  const double t = Rinv*dz;
1101 
1102  //2. Conformal mapping
1103 
1104  m_SoA.m_r[idx] = Rinv;
1105  m_SoA.m_u[idx] = xn*R2inv;
1106  m_SoA.m_v[idx] = yn*R2inv;
1107 
1108  //3. cot(theta) estimation for the doublet
1109 
1110  const double covZP = pSP->dz()*pSP->dz();
1111  const double covRP = pSP->dr()*pSP->dr();
1112 
1113  m_SoA.m_t[idx] = t;
1114  m_SoA.m_tCov[idx] = R2inv*(covZ + covZP + t*t*(covR+covRP));
1115  }
1116 
1117  //double loop
1118 
1119  for(int innIdx=0;innIdx<nInner;innIdx++) {//loop over inner spacepoints
1120 
1121  //mult. scatt contribution due to the layer with middle SP
1122 
1123  const double r_inn = m_SoA.m_r[innIdx];
1124  const double t_inn = m_SoA.m_t[innIdx];
1125  const double v_inn = m_SoA.m_v[innIdx];
1126  const double u_inn = m_SoA.m_u[innIdx];
1127  const double tCov_inn = m_SoA.m_tCov[innIdx];
1128  const double dCov = m_CovMS*(1+t_inn*t_inn);
1129 
1130  std::vector<ProtoSeed> vPro;
1131  vPro.reserve(100);
1132 
1133  for(int outIdx=nInner;outIdx<nSP;outIdx++) {
1134 
1135  //1. rz doublet matching
1136  const double t_out = m_SoA.m_t[outIdx];
1137 
1138  const double dt2 = std::pow((t_inn - t_out), 2)*(1.0/9.0);//i.e. 3-sigma cut
1139 
1140 
1141  double covdt = (t_inn*t_out*covR + covZ);
1142  covdt *= 2*r_inn*m_SoA.m_r[outIdx];
1143  covdt += tCov_inn + m_SoA.m_tCov[outIdx];
1144 
1145  if(dt2 > covdt+dCov) continue;
1146 
1147  //2. pT estimate
1148 
1149  const double du = m_SoA.m_u[outIdx] - u_inn;
1150  if(du==0.0) continue;
1151  const double A = (m_SoA.m_v[outIdx] - v_inn)/du;
1152  const double B = v_inn - A*u_inn;
1153  const double R_squ = (1 + A*A)/(B*B);
1154 
1155  if(R_squ < m_minR_squ) continue;
1156 
1157 
1158  //3. the 3-sigma cut with estimated pT
1159 
1160  const double frac = m_minR_squ/R_squ;
1161  if(dt2 > covdt+frac*dCov) continue;
1162 
1163  //4. d0 cut
1164 
1165  const double fabs_d0 = std::fabs(pS_r*(B*pS_r - A));
1166 
1167  if(fabs_d0 > m_settings.m_tripletD0Max) continue;
1168 
1169  if (m_SoA.m_spo[outIdx-nInner]->isSCT() && isPixel) {
1170  if(fabs_d0 > m_settings.m_tripletD0_PPS_Max) continue;
1171  }
1172 
1173  // Calculate triplet weight: No weighting in LRT mode, but d0**2 weighting for normal (non LRT) mode. Triplets are then sorted by lowest weight.
1174 
1175  const double Q= (m_settings.m_LRTmode ? 0 : fabs_d0*fabs_d0);
1176 
1177  bool isOuterPixel = m_SoA.m_spo[outIdx-nInner]->isPixel();
1178 
1179  if(!m_settings.m_LRTmode) {
1180  if (isOuterPixel && (bestQ < Q)) continue;//this PPP seed has worse quality
1181  }
1182 
1183  //5. phi0 cut
1184 
1185  if ( !fullPhi ) {
1186  const double uc = 2*B*pS_r - A;
1187  const double phi0 = atan2(sinA - uc*cosA, cosA + uc*sinA);
1188 
1189  if ( !RoiUtil::containsPhi( *roiDescriptor, phi0 ) ) {
1190  continue;
1191  }
1192  }
1193 
1194  //6. add new triplet
1195 
1196  const double Curv = B/std::sqrt(1+A*A);
1197 
1198 
1199  bool isConfirmed = false;
1200 
1201  for(auto& ps : vPro) {
1202  double diffC = 1 - Curv/ps.m_curv;
1203  if(std::fabs(diffC) < m_settings.m_curv_delta) {
1204  ps.m_confirmed = true;
1205  isConfirmed = true;
1206  }
1207  }
1208 
1209  if(!isOuterPixel) {
1210  continue;//we use PPS seeds only to confirm PPP seeds
1211  }
1212 
1213  vPro.emplace_back(ProtoSeed(m_SoA.m_spo[outIdx-nInner], Curv, Q, isConfirmed));
1214 
1215  }//loop over outer spacepoints
1216 
1217  for(const auto& ps : vPro) {
1218  if(!ps.m_confirmed) continue;
1219  if(output.size()>=m_settings.m_maxTripletBufferLength) {
1220  if (m_settings.m_LRTmode) { // take the first m_maxTripletBufferLength triplets
1221  continue;
1222  }
1223  else { // choose largest d0
1224 
1225  std::sort(output.begin(), output.end(),
1226  [](const TrigInDetTriplet& A, const TrigInDetTriplet& B) {
1227  return A.Q() < B.Q();
1228  }
1229  );
1230  double QL = output.back().Q();
1231  if(bestQ > QL) {
1232  bestQ = QL;
1233  }
1234  if( ps.m_Q >= QL) {//reject this seed
1235  continue;
1236  }
1237  output.pop_back();
1238  }
1239  }
1240 
1241  output.emplace_back(*m_SoA.m_spi[innIdx], *pS, *ps.m_sp, ps.m_Q);
1242  }
1243  }
1244 }
1245 
1246 
1247 
1248 void TrigTrackSeedGenerator::storeTriplets(std::vector<TrigInDetTriplet>& tripletVec) {
1249  for(std::vector<TrigInDetTriplet>::iterator it=tripletVec.begin();it!=tripletVec.end();++it) {
1250  float newQ = (*it).Q();
1251 
1252  if((!m_settings.m_LRTmode) && (*it).s3().isSCT()) {
1253  // In normal (non LRT) mode penalise SSS by 1000, PSS (if enabled) and PPS by 10000
1254  newQ += (*it).s1().isSCT() ? 1000.0 : 10000.0;
1255  }
1256  (*it).Q(newQ);
1257  m_triplets.push_back((*it));
1258  }
1259 }
1260 
1261 void TrigTrackSeedGenerator::getSeeds(std::vector<TrigInDetTriplet>& vs) {
1262  vs.clear();
1263  std::sort(m_triplets.begin(), m_triplets.end(),
1264  [](const TrigInDetTriplet& A, const TrigInDetTriplet& B) {
1265  return A.Q() < B.Q();
1266  }
1267  );
1268  vs = std::move(m_triplets);
1269 }
TrigCombinatorialSettings::m_nMaxPhiSlice
int m_nMaxPhiSlice
Definition: TrigCombinatorialSettings.h:67
TrigInDetTriplet
Definition: TrigInDetTriplet.h:13
Trk::SpacePoint::clusterList
const std::pair< const PrepRawData *, const PrepRawData * > & clusterList() const
return the pair of cluster pointers by reference
Definition: Tracking/TrkEvent/TrkSpacePoint/TrkSpacePoint/SpacePoint.h:127
xAOD::iterator
JetConstituentVector::iterator iterator
Definition: JetConstituentVector.cxx:68
TrigTrackSeedGenerator::m_CovMS
double m_CovMS
Definition: TrigTrackSeedGenerator.h:348
Trk::SpacePoint
Definition: Tracking/TrkEvent/TrkSpacePoint/TrkSpacePoint/SpacePoint.h:35
LPhi_Storage::m_layers
std::vector< L_PHI_SECTOR > m_layers
Definition: TrigTrackSeedGenerator.h:234
InternalSoA::m_t
double * m_t
Definition: TrigTrackSeedGenerator.h:297
TrigCombinatorialSettings::m_zvErrorEndcap
float m_zvErrorEndcap
Definition: TrigCombinatorialSettings.h:71
python.SystemOfUnits.s
int s
Definition: SystemOfUnits.py:131
TrigTrackSeedGenerator::m_triplets
std::vector< TrigInDetTriplet > m_triplets
Definition: TrigTrackSeedGenerator.h:350
IndexedSP::m_idx
int m_idx
Definition: TrigTrackSeedGenerator.h:28
max
#define max(a, b)
Definition: cfImp.cxx:41
TrigTrackSeedGenerator::validateLayerPairNew
bool validateLayerPairNew(int, int, float, float)
Definition: TrigTrackSeedGenerator.cxx:291
TrigTrackSeedGenerator::m_maxDeltaRadius
double m_maxDeltaRadius
Definition: TrigTrackSeedGenerator.h:342
TrigCombinatorialSettings::m_tripletD0_PPS_Max
float m_tripletD0_PPS_Max
Definition: TrigCombinatorialSettings.h:58
CaloCellPos2Ntuple.int
int
Definition: CaloCellPos2Ntuple.py:24
InternalSoA::m_r
double * m_r
Definition: TrigTrackSeedGenerator.h:294
InternalSoA::clear
void clear()
Definition: TrigTrackSeedGenerator.h:249
TrigTrackSeedGenerator::createTripletsNew
void createTripletsNew(const TrigSiSpacePointBase *, int, int, std::vector< TrigInDetTriplet > &, const IRoiDescriptor *)
Definition: TrigTrackSeedGenerator.cxx:778
TrigCombinatorialSettings::m_tripletDtCut
float m_tripletDtCut
Definition: TrigCombinatorialSettings.h:66
TrigSiSpacePointBase::isPixel
bool isPixel() const
Definition: TrigSiSpacePointBase.h:70
SP_RANGE
std::pair< std::vector< const INDEXED_SP * >::const_iterator, std::vector< const INDEXED_SP * >::const_iterator > SP_RANGE
Definition: TrigTrackSeedGenerator.h:308
TrigTrackSeedGenerator::m_maxCoord
float m_maxCoord
Definition: TrigTrackSeedGenerator.h:352
PixelCluster.h
InternalSoA::m_sorted_sp_t
double * m_sorted_sp_t
Definition: TrigTrackSeedGenerator.h:303
TrigTrackSeedGenerator.h
InternalSoA::m_ti
double * m_ti
Definition: TrigTrackSeedGenerator.h:298
InDet::SiWidth::widthPhiRZ
const Amg::Vector2D & widthPhiRZ() const
Definition: SiWidth.h:121
LPhi_Storage::reset
void reset()
Definition: TrigTrackSeedGenerator.h:211
LPhiSector::greaterThanR_i
Definition: TrigTrackSeedGenerator.h:72
InternalSoA::m_v
double * m_v
Definition: TrigTrackSeedGenerator.h:296
InDetAccessor::phi0
@ phi0
Definition: InDetAccessor.h:33
conifer::pow
constexpr int pow(int x)
Definition: conifer.h:20
TrigTrackSeedGenerator::m_minR_squ
double m_minR_squ
Definition: TrigTrackSeedGenerator.h:348
RoiUtil::containsPhi
bool containsPhi(const IRoiDescriptor &roi, double phi)
test whether a stub is contained within the roi
Definition: RoiUtil.cxx:79
ALFA_EventTPCnv_Dict::t1
std::vector< ALFA_RawDataCollection_p1 > t1
Definition: ALFA_EventTPCnvDict.h:43
skel.it
it
Definition: skel.GENtoEVGEN.py:423
TrigCombinatorialSettings::m_layerGeometry
std::vector< TrigInDetSiLayer > m_layerGeometry
Definition: TrigCombinatorialSettings.h:78
TrigTrackSeedGenerator::m_zTol
double m_zTol
Definition: TrigTrackSeedGenerator.h:342
M_PI
#define M_PI
Definition: ActiveFraction.h:11
TrigCombinatorialSettings::m_curv_delta
float m_curv_delta
Definition: TrigCombinatorialSettings.h:64
LPhiSector::smallerThanR_i
Definition: TrigTrackSeedGenerator.h:84
TrigTrackSeedGenerator::m_zMinusEndcap
float m_zMinusEndcap
Definition: TrigTrackSeedGenerator.h:352
TrigCombinatorialSettings::m_LRTmode
bool m_LRTmode
Definition: TrigCombinatorialSettings.h:72
TrigCombinatorialSettings::m_tripletPtMin
float m_tripletPtMin
Definition: TrigCombinatorialSettings.h:59
TrigSiSpacePointBase::x
void x(const double x)
Definition: TrigSiSpacePointBase.h:54
TrigTrackSeedGenerator::m_spStorage
std::vector< INDEXED_SP > m_spStorage
Definition: TrigTrackSeedGenerator.h:327
IndexedSP::m_pSP
const TrigSiSpacePointBase * m_pSP
Definition: TrigTrackSeedGenerator.h:27
python.TurnDataReader.dr
dr
Definition: TurnDataReader.py:112
read_hist_ntuple.t
t
Definition: read_hist_ntuple.py:5
LPhiSector
Definition: TrigTrackSeedGenerator.h:32
TrigTrackSeedGenerator::TrigTrackSeedGenerator
TrigTrackSeedGenerator(const TrigCombinatorialSettings &)
Definition: TrigTrackSeedGenerator.cxx:18
TrigCombinatorialSettings::m_tripletDoConfirm
bool m_tripletDoConfirm
Definition: TrigCombinatorialSettings.h:63
TrigTrackSeedGenerator::getSpacepointRange
bool getSpacepointRange(int, const std::vector< const INDEXED_SP * > &, SP_RANGE &)
Definition: TrigTrackSeedGenerator.cxx:426
InternalSoA::m_to
double * m_to
Definition: TrigTrackSeedGenerator.h:299
JetTiledMap::S
@ S
Definition: TiledEtaPhiMap.h:44
TrigSeedML_LUT
Definition: TrigSeedML_LUT.h:10
TrigTrackSeedGenerator::m_settings
const TrigCombinatorialSettings & m_settings
Definition: TrigTrackSeedGenerator.h:340
TrigTrackSeedGenerator::processSpacepointRange
int processSpacepointRange(int, const INDEXED_SP *, bool, const SP_RANGE &, const IRoiDescriptor *)
Definition: TrigTrackSeedGenerator.cxx:467
L_PHI_STORAGE
struct LPhi_Storage L_PHI_STORAGE
TrigTrackSeedGenerator::loadSpacePoints
void loadSpacePoints(const std::vector< TrigSiSpacePointBase > &)
Definition: TrigTrackSeedGenerator.cxx:44
python.Bindings.values
values
Definition: Control/AthenaPython/python/Bindings.py:797
InternalSoA::m_spi
const TrigSiSpacePointBase ** m_spi
Definition: TrigTrackSeedGenerator.h:292
TrigTrackSeedGenerator::m_minCoord
float m_minCoord
Definition: TrigTrackSeedGenerator.h:352
TrigSiSpacePointBase::z
void z(const double z)
Definition: TrigSiSpacePointBase.h:53
TrigTrackSeedGenerator::m_outerMarkers
std::vector< int > m_outerMarkers
Definition: TrigTrackSeedGenerator.h:355
dqt_zlumi_alleff_HIST.A
A
Definition: dqt_zlumi_alleff_HIST.py:110
TrigCombinatorialSettings::m_maxEC_len
float m_maxEC_len
Definition: TrigCombinatorialSettings.h:82
TrigInDetSiLayer.h
TrigSiSpacePointBase::r
void r(const double r)
Definition: TrigSiSpacePointBase.h:51
InternalSoA::m_u
double * m_u
Definition: TrigTrackSeedGenerator.h:295
TrigCombinatorialSettings::m_zvError
float m_zvError
Definition: TrigCombinatorialSettings.h:70
InternalSoA::m_sorted_sp
const TrigSiSpacePointBase ** m_sorted_sp
Definition: TrigTrackSeedGenerator.h:301
TrigSiSpacePointBase::dr
void dr(const double dr)
Definition: TrigSiSpacePointBase.h:56
InternalSoA::m_tCov
double * m_tCov
Definition: TrigTrackSeedGenerator.h:300
python.ConfigurableDb.conf
def conf
Definition: ConfigurableDb.py:282
TrigSiSpacePointBase::isSCT
bool isSCT() const
Definition: TrigSiSpacePointBase.h:71
TrigCombinatorialSettings::m_doublet_dR_Max
float m_doublet_dR_Max
Definition: TrigCombinatorialSettings.h:54
TrigTrackSeedGenerator::m_maxDeltaRadiusConf
double m_maxDeltaRadiusConf
Definition: TrigTrackSeedGenerator.h:342
TrigTrackSeedGenerator::processSpacepointRangeZv
int processSpacepointRangeZv(const INDEXED_SP *, bool, const SP_RANGE &, bool, const float &, const float &)
Definition: TrigTrackSeedGenerator.cxx:547
TrigSiSpacePointBase.h
DetDescrDictionaryDict::it1
std::vector< HWIdentifier >::iterator it1
Definition: DetDescrDictionaryDict.h:17
TrigTrackSeedGenerator::m_nOuter
int m_nOuter
Definition: TrigTrackSeedGenerator.h:354
IRoiDescriptor
Describes the API of the Region of Ineterest geometry.
Definition: IRoiDescriptor.h:23
checkxAOD.frac
frac
Definition: Tools/PyUtils/bin/checkxAOD.py:256
TrigTrackSeedGenerator::m_dtPreCut
double m_dtPreCut
Definition: TrigTrackSeedGenerator.h:348
TrigCombinatorialSettings::m_tripletDoPPS
bool m_tripletDoPPS
Definition: TrigCombinatorialSettings.h:62
InternalSoA::m_spo
const TrigSiSpacePointBase ** m_spo
Definition: TrigTrackSeedGenerator.h:293
TrigTrackSeedGenerator::m_phiSliceWidth
double m_phiSliceWidth
Definition: TrigTrackSeedGenerator.h:341
TrigTrackSeedGenerator::m_zMinus
float m_zMinus
Definition: TrigTrackSeedGenerator.h:352
TrigTrackSeedGenerator::~TrigTrackSeedGenerator
~TrigTrackSeedGenerator()
Definition: TrigTrackSeedGenerator.cxx:37
TrigSiSpacePointBase::y
void y(const double y)
Definition: TrigSiSpacePointBase.h:55
LPhiSector::greaterThanR
Definition: TrigTrackSeedGenerator.h:66
InternalSoA::m_sorted_sp_type
int * m_sorted_sp_type
Definition: TrigTrackSeedGenerator.h:302
LPhiSector::smallerThanZ
Definition: TrigTrackSeedGenerator.h:60
TrigTrackSeedGenerator::getSeeds
void getSeeds(std::vector< TrigInDetTriplet > &)
Definition: TrigTrackSeedGenerator.cxx:1261
IRoiDescriptor::phiMinus
virtual double phiMinus() const =0
drawFromPickle.tan
tan
Definition: drawFromPickle.py:36
TRT::Track::z0
@ z0
Definition: InnerDetector/InDetCalibEvent/TRT_CalibData/TRT_CalibData/TrackInfo.h:63
DeMoUpdate.tmp
string tmp
Definition: DeMoUpdate.py:1167
TrigCombinatorialSettings::m_maxTripletBufferLength
unsigned int m_maxTripletBufferLength
Definition: TrigCombinatorialSettings.h:68
TrigCombinatorialSettings::m_tripletDoPSS
bool m_tripletDoPSS
Definition: TrigCombinatorialSettings.h:61
min
#define min(a, b)
Definition: cfImp.cxx:40
IRoiDescriptor::phiPlus
virtual double phiPlus() const =0
extreme phi values
TrigTrackSeedGenerator::m_minDeltaRadius
double m_minDeltaRadius
Definition: TrigTrackSeedGenerator.h:342
merge.output
output
Definition: merge.py:17
TrigTrackSeedGenerator::createTriplets
void createTriplets(const TrigSiSpacePointBase *, int, int, std::vector< TrigInDetTriplet > &, const IRoiDescriptor *)
Definition: TrigTrackSeedGenerator.cxx:602
TrigTrackSeedGenerator::m_maxTau
std::vector< float > m_maxTau
Definition: TrigTrackSeedGenerator.h:329
TrigCombinatorialSettings::m_useTrigSeedML
int m_useTrigSeedML
Definition: TrigCombinatorialSettings.h:80
RoiUtil.h
TrigSiSpacePointBase::dz
void dz(const double dz)
Definition: TrigSiSpacePointBase.h:57
TrigTrackSeedGenerator::m_SoA
INTERNAL_SOA m_SoA
Definition: TrigTrackSeedGenerator.h:346
dqt_zlumi_alleff_HIST.B
B
Definition: dqt_zlumi_alleff_HIST.py:110
dirs
std::map< std::string, int > dirs
list of directories to be explicitly included, together with corresponding depths of subdirectories
Definition: hcg.cxx:99
TrigTrackSeedGenerator::createSeeds
void createSeeds(const IRoiDescriptor *)
Definition: TrigTrackSeedGenerator.cxx:111
InternalSoA::resize
void resize(const int spSize)
Definition: TrigTrackSeedGenerator.h:276
TrigTrackSeedGenerator::m_innerMarkers
std::vector< int > m_innerMarkers
Definition: TrigTrackSeedGenerator.h:355
InDet::PixelCluster
Definition: InnerDetector/InDetRecEvent/InDetPrepRawData/InDetPrepRawData/PixelCluster.h:49
LPhi_Storage::addSpacePoint
void addSpacePoint(int phiIdx, int layerId, const INDEXED_SP *p)
Definition: TrigTrackSeedGenerator.h:207
makeTRTBarrelCans.dy
tuple dy
Definition: makeTRTBarrelCans.py:21
TrigCombinatorialSettings::m_doubletFilterRZ
bool m_doubletFilterRZ
Definition: TrigCombinatorialSettings.h:65
InDet::SiCluster::width
const InDet::SiWidth & width() const
return width class reference
TrigCombinatorialSettings::m_vLUT
std::vector< TrigSeedML_LUT > m_vLUT
Definition: TrigCombinatorialSettings.h:81
DiTauMassTools::MaxHistStrategyV2::e
e
Definition: PhysicsAnalysis/TauID/DiTauMassTools/DiTauMassTools/HelperFunctions.h:26
LPhiSector::greaterThanZ
Definition: TrigTrackSeedGenerator.h:54
IRoiDescriptor::zedPlus
virtual double zedPlus() const =0
the zed and eta values at the most forward and most rear ends of the RoI
TrigTrackSeedGenerator::createConfirmedTriplets
void createConfirmedTriplets(const TrigSiSpacePointBase *, int, int, std::vector< TrigInDetTriplet > &, const IRoiDescriptor *)
Definition: TrigTrackSeedGenerator.cxx:1022
LPhiSector::smallerThanR
Definition: TrigTrackSeedGenerator.h:78
TrigCombinatorialSettings::m_doublet_dR_Max_Confirm
float m_doublet_dR_Max_Confirm
Definition: TrigCombinatorialSettings.h:55
python.CaloScaleNoiseConfig.type
type
Definition: CaloScaleNoiseConfig.py:78
TrigTrackSeedGenerator::m_zPlus
float m_zPlus
Definition: TrigTrackSeedGenerator.h:352
makeTRTBarrelCans.dx
tuple dx
Definition: makeTRTBarrelCans.py:20
IndexedSP
Definition: TrigTrackSeedGenerator.h:16
IRoiDescriptor::isFullscan
virtual bool isFullscan() const =0
is this a full detector RoI?
LArNewCalib_DelayDump_OFC_Cali.idx
idx
Definition: LArNewCalib_DelayDump_OFC_Cali.py:69
python.LArCondContChannels.isBarrel
isBarrel
Definition: LArCondContChannels.py:659
IRoiDescriptor::zedMinus
virtual double zedMinus() const =0
TrigTrackSeedGenerator::m_nInner
int m_nInner
Definition: TrigTrackSeedGenerator.h:354
IRoiDescriptor.h
TrigCombinatorialSettings::m_tripletD0Max
float m_tripletD0Max
Definition: TrigCombinatorialSettings.h:57
TrigTrackSeedGenerator::storeTriplets
void storeTriplets(std::vector< TrigInDetTriplet > &)
Definition: TrigTrackSeedGenerator.cxx:1248
RoiUtil::contains
bool contains(const IRoiDescriptor &roi, double z0, double dzdr)
see whether a segment is contained within the roi in r-z
Definition: RoiUtil.cxx:42
TrigSiSpacePointBase
Definition: TrigSiSpacePointBase.h:23
TrigInDetTriplet.h
Amg::distance
float distance(const Amg::Vector3D &p1, const Amg::Vector3D &p2)
calculates the distance between two point in 3D space
Definition: GeoPrimitivesHelpers.h:54
TrigTrackSeedGenerator::m_zPlusEndcap
float m_zPlusEndcap
Definition: TrigTrackSeedGenerator.h:352
TrigTrackSeedGenerator::m_pStore
L_PHI_STORAGE * m_pStore
Definition: TrigTrackSeedGenerator.h:344
python.compressB64.c
def c
Definition: compressB64.py:93
TrigCombinatorialSettings
Definition: TrigCombinatorialSettings.h:15
LPhi_Storage::sortSpacePoints2
void sortSpacePoints2(const std::vector< TrigInDetSiLayer > &layerGeometry)
Definition: TrigTrackSeedGenerator.h:226
fitman.k
k
Definition: fitman.py:528
TrigTrackSeedGenerator::m_minTau
std::vector< float > m_minTau
Definition: TrigTrackSeedGenerator.h:328
TrigSeedML_LUT::getValidRange
bool getValidRange(float fX, float &min, float &max) const
Definition: TrigSeedML_LUT.h:37
Generated on Fri Jul 5 2024 21:31:15 for ATLAS Offline Software by   doxygen 1.8.18