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