ATLAS Offline Software
Loading...
Searching...
No Matches
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
18TrigTrackSeedGenerator::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
51
52void 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
119void 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
209void 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
299bool 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
434bool 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
475int 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
555int 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
610void 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
786void 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
1030void 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
1256void 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
1269void 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}
M_PI
#define M_PI
Definition ActiveFraction.h:11
L_PHI_SECTOR
struct LPhiSector L_PHI_SECTOR
INDEXED_SP
struct IndexedSP INDEXED_SP
L_PHI_STORAGE
struct LPhi_Storage L_PHI_STORAGE
SP_RANGE
std::pair< std::vector< const INDEXED_SP * >::const_iterator, std::vector< const INDEXED_SP * >::const_iterator > SP_RANGE
Definition TrigTrackSeedGenerator.h:308
IRoiDescriptor
Describes the API of the Region of Ineterest geometry.
Definition IRoiDescriptor.h:23
IRoiDescriptor::isFullscan
virtual bool isFullscan() const =0
is this a full detector RoI?
IRoiDescriptor::phiPlus
virtual double phiPlus() const =0
extreme phi values
IRoiDescriptor::zedPlus
virtual double zedPlus() const =0
the zed and eta values at the most forward and most rear ends of the RoI
IRoiDescriptor::phiMinus
virtual double phiMinus() const =0
IRoiDescriptor::zedMinus
virtual double zedMinus() const =0
InDet::SiCluster::width
const InDet::SiWidth & width() const
return width class reference
InDet::SiWidth::widthPhiRZ
const Amg::Vector2D & widthPhiRZ() const
Definition SiWidth.h:121
TrigSiSpacePointBase::r
void r(const double r)
Definition TrigSiSpacePointBase.h:51
TrigSiSpacePointBase::dr
void dr(const double dr)
Definition TrigSiSpacePointBase.h:56
TrigSiSpacePointBase::y
void y(const double y)
Definition TrigSiSpacePointBase.h:55
TrigSiSpacePointBase::dz
void dz(const double dz)
Definition TrigSiSpacePointBase.h:57
TrigSiSpacePointBase::x
void x(const double x)
Definition TrigSiSpacePointBase.h:54
TrigSiSpacePointBase::z
void z(const double z)
Definition TrigSiSpacePointBase.h:53
TrigTrackSeedGenerator::m_settings
const TrigCombinatorialSettings & m_settings
Definition TrigTrackSeedGenerator.h:340
TrigTrackSeedGenerator::loadSpacePoints
void loadSpacePoints(const std::vector< TrigSiSpacePointBase > &)
Definition TrigTrackSeedGenerator.cxx:52
TrigTrackSeedGenerator::m_maxTau
std::vector< float > m_maxTau
Definition TrigTrackSeedGenerator.h:329
TrigTrackSeedGenerator::createTripletsNew
void createTripletsNew(const TrigSiSpacePointBase *, int, int, std::vector< TrigInDetTriplet > &, const IRoiDescriptor *)
Definition TrigTrackSeedGenerator.cxx:786
TrigTrackSeedGenerator::storeTriplets
void storeTriplets(std::vector< TrigInDetTriplet > &)
Definition TrigTrackSeedGenerator.cxx:1256
TrigTrackSeedGenerator::m_minTau
std::vector< float > m_minTau
Definition TrigTrackSeedGenerator.h:328
TrigTrackSeedGenerator::getSpacepointRange
bool getSpacepointRange(int, const std::vector< const INDEXED_SP * > &, SP_RANGE &)
Definition TrigTrackSeedGenerator.cxx:434
TrigTrackSeedGenerator::createSeeds
void createSeeds(const IRoiDescriptor *)
Definition TrigTrackSeedGenerator.cxx:119
TrigTrackSeedGenerator::getSeeds
void getSeeds(std::vector< TrigInDetTriplet > &)
Definition TrigTrackSeedGenerator.cxx:1269
TrigTrackSeedGenerator::m_spStorage
std::vector< INDEXED_SP > m_spStorage
Definition TrigTrackSeedGenerator.h:327
TrigTrackSeedGenerator::createTriplets
void createTriplets(const TrigSiSpacePointBase *, int, int, std::vector< TrigInDetTriplet > &, const IRoiDescriptor *)
Definition TrigTrackSeedGenerator.cxx:610
TrigTrackSeedGenerator::processSpacepointRangeZv
int processSpacepointRangeZv(const INDEXED_SP *, bool, const SP_RANGE &, bool, const float &, const float &)
Definition TrigTrackSeedGenerator.cxx:555
TrigTrackSeedGenerator::m_innerMarkers
std::vector< int > m_innerMarkers
Definition TrigTrackSeedGenerator.h:355
TrigTrackSeedGenerator::m_triplets
std::vector< TrigInDetTriplet > m_triplets
Definition TrigTrackSeedGenerator.h:350
TrigTrackSeedGenerator::createConfirmedTriplets
void createConfirmedTriplets(const TrigSiSpacePointBase *, int, int, std::vector< TrigInDetTriplet > &, const IRoiDescriptor *)
Definition TrigTrackSeedGenerator.cxx:1030
TrigTrackSeedGenerator::TrigTrackSeedGenerator
TrigTrackSeedGenerator(const TrigCombinatorialSettings &)
Definition TrigTrackSeedGenerator.cxx:18
TrigTrackSeedGenerator::validateLayerPairNew
bool validateLayerPairNew(int, int, float, float)
Definition TrigTrackSeedGenerator.cxx:299
TrigTrackSeedGenerator::m_outerMarkers
std::vector< int > m_outerMarkers
Definition TrigTrackSeedGenerator.h:355
TrigTrackSeedGenerator::processSpacepointRange
int processSpacepointRange(int, const INDEXED_SP *, bool, const SP_RANGE &, const IRoiDescriptor *)
Definition TrigTrackSeedGenerator.cxx:475
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
dirs
std::map< std::string, int > dirs
list of directories to be explicitly included, together with corresponding depths of subdirectories
Definition hcg.cxx:102
RoiUtil::containsPhi
bool containsPhi(const IRoiDescriptor &roi, double phi)
test whether a stub is contained within the roi
Definition RoiUtil.cxx:79
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
std::sort
void sort(typename DataModel_detail::iterator< DVL > beg, typename DataModel_detail::iterator< DVL > end)
Specialization of sort for DataVector/List.
Definition DVL_algorithms.h:554
A
hold the test vectors and ease the comparison
type
IndexedSP::m_pSP
const TrigSiSpacePointBase * m_pSP
Definition TrigTrackSeedGenerator.h:27
LPhiSector::greaterThanR_i
Definition TrigTrackSeedGenerator.h:72
LPhiSector::smallerThanR_i
Definition TrigTrackSeedGenerator.h:84
TrigSeedML_LUT
Definition TrigInDetPattRecoTools/TrigInDetPattRecoTools/TrigSeedML_LUT.h:10
TrigSeedML_LUT::getValidRange
bool getValidRange(float fX, float &min, float &max) const
Definition TrigInDetPattRecoTools/TrigInDetPattRecoTools/TrigSeedML_LUT.h:37
Generated on for ATLAS Offline Software by doxygen 1.14.0