Project
Loading...
Searching...
No Matches
GPUTRDTracker.cxx
Go to the documentation of this file.
1// Copyright 2019-2020 CERN and copyright holders of ALICE O2.
2// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders.
3// All rights not expressly granted are reserved.
4//
5// This software is distributed under the terms of the GNU General Public
6// License v3 (GPL Version 3), copied verbatim in the file "COPYING".
7//
8// In applying this license CERN does not waive the privileges and immunities
9// granted to it by virtue of its status as an Intergovernmental Organization
10// or submit itself to any jurisdiction.
11
14
15// #define ENABLE_GPUTRDDEBUG
16#define ENABLE_WARNING 0
17#define ENABLE_INFO 0
18
19#include "GPUTRDTracker.h"
20#include "GPUTRDTrackletWord.h"
21#include "GPUTRDGeometry.h"
22#include "GPUTRDTrackerDebug.h"
23#include "GPUCommonMath.h"
24#include "GPUCommonAlgorithm.h"
25#include "GPUConstantMem.h"
26
27using namespace o2::gpu;
28
30
31#ifndef GPUCA_GPUCODE
32#include "GPUMemoryResource.h"
33#include "GPUReconstruction.h"
34#include <chrono>
35#include <vector>
36
37#include "GPUChainTracking.h"
38
39template <class TRDTRK, class PROP>
41{
42 mNMaxTracks = std::max(std::max(io.nOutputTracksTPCO2, io.nTracksTPCITSO2), std::max(io.nMergedTracks, io.nOutputTracksTPCO2)); // TODO: This is a bit stupid, we should just take the correct number, not the max of all
43 mNMaxSpacePoints = io.nTRDTracklets;
44 mNMaxCollisions = io.nTRDTriggerRecords;
45}
46
47template <class TRDTRK, class PROP>
49{
50 AllocateAndInitializeLate();
51 mMemoryPermanent = mRec->RegisterMemoryAllocation(this, &GPUTRDTracker_t<TRDTRK, PROP>::SetPointersBase, GPUMemoryResource::MEMORY_PERMANENT, "TRDInitialize");
52 mMemoryTracklets = mRec->RegisterMemoryAllocation(this, &GPUTRDTracker_t<TRDTRK, PROP>::SetPointersTracklets, GPUMemoryResource::MEMORY_INPUT, "TRDTracklets");
53 mMemoryTracks = mRec->RegisterMemoryAllocation(this, &GPUTRDTracker_t<TRDTRK, PROP>::SetPointersTracks, GPUMemoryResource::MEMORY_INOUT, "TRDTracks");
56template <class TRDTRK, class PROP>
59 //--------------------------------------------------------------------
60 // Allocate memory for fixed size objects (needs to be done only once)
61 //--------------------------------------------------------------------
62 mMaxBackendThreads = mRec->GetMaxBackendThreads();
63 computePointerWithAlignment(base, mR, kNChambers);
64 computePointerWithAlignment(base, mHypothesis, mNCandidates * mMaxBackendThreads);
65 computePointerWithAlignment(base, mCandidates, mNCandidates * 2 * mMaxBackendThreads);
66 return base;
67}
68
69template <class TRDTRK, class PROP>
71{
72 //--------------------------------------------------------------------
73 // Allocate memory for tracklets and space points
74 // (size might change for different events)
75 //--------------------------------------------------------------------
76 if (mGenerateSpacePoints) {
77 computePointerWithAlignment(base, mSpacePoints, mNMaxSpacePoints);
78 }
79 computePointerWithAlignment(base, mTrackletIndexArray, (kNChambers + 1) * mNMaxCollisions);
80 return base;
81}
82
83template <class TRDTRK, class PROP>
85{
86 //--------------------------------------------------------------------
87 // Allocate memory for tracks (this is done once per event)
88 //--------------------------------------------------------------------
89 computePointerWithAlignment(base, mTracks, mNMaxTracks);
90 computePointerWithAlignment(base, mTrackAttribs, mNMaxTracks);
91 return base;
92}
93
94template <class TRDTRK, class PROP>
95GPUTRDTracker_t<TRDTRK, PROP>::GPUTRDTracker_t() : mR(nullptr), mIsInitialized(false), mGenerateSpacePoints(false), mProcessPerTimeFrame(false), mNAngleHistogramBins(25), mAngleHistogramRange(50), mMemoryPermanent(-1), mMemoryTracklets(-1), mMemoryTracks(-1), mNMaxCollisions(0), mNMaxTracks(0), mNMaxSpacePoints(0), mTracks(nullptr), mTrackAttribs(nullptr), mNCandidates(1), mNTracks(0), mNEvents(0), mMaxBackendThreads(100), mTrackletIndexArray(nullptr), mHypothesis(nullptr), mCandidates(nullptr), mSpacePoints(nullptr), mGeo(nullptr), mRPhiA2(0), mRPhiB(0), mRPhiC2(0), mDyA2(0), mDyB(0), mDyC2(0), mAngleToDyA(0), mAngleToDyB(0), mAngleToDyC(0), mDebugOutput(false), mMaxEta(0.84f), mRoadZ(18.f), mZCorrCoefNRC(1.4f), mTPCVdrift(2.58f), mTPCTDriftOffset(0.f), mDebug(new GPUTRDTrackerDebug<TRDTRK>())
96{
97 //--------------------------------------------------------------------
98 // Default constructor
99 //--------------------------------------------------------------------
100}
101
102template <class TRDTRK, class PROP>
104{
105 //--------------------------------------------------------------------
106 // Destructor
107 //--------------------------------------------------------------------
108 delete mDebug;
109}
110
111template <class TRDTRK, class PROP>
113{
114 //--------------------------------------------------------------------
115 // Initialise tracker
116 //--------------------------------------------------------------------
117
118 UpdateGeometry();
119
120 mDebug->ExpandVectors();
121 mIsInitialized = true;
122}
123
124template <class TRDTRK, class PROP>
126{
127 //--------------------------------------------------------------------
128 // Update Geometry of TRDTracker
129 //--------------------------------------------------------------------
130 mGeo = (const GPUTRDGeometry*)GetConstantMem()->calibObjects.trdGeometry;
131 if (!mGeo) {
132 GPUFatal("TRD geometry must be provided externally");
133 }
134 float Bz = Param().bzkG;
135 float resRPhiIdeal2 = Param().rec.trd.trkltResRPhiIdeal * Param().rec.trd.trkltResRPhiIdeal;
136 GPUInfo("Initializing with B-field: %f kG", Bz);
137 if (CAMath::Abs(CAMath::Abs(Bz) - 2) < 0.1f) {
138 // magnetic field +-0.2 T
139 if (Bz > 0) {
140 GPUInfo("Loading error parameterization for Bz = +2 kG");
141 mRPhiA2 = resRPhiIdeal2, mRPhiB = -1.43e-2f, mRPhiC2 = 4.55e-2f;
142 mDyA2 = 1.225e-3f, mDyB = -9.8e-3f, mDyC2 = 3.88e-2f;
143 mAngleToDyA = -0.1f, mAngleToDyB = 1.89f, mAngleToDyC = -0.4f;
144 } else {
145 GPUInfo("Loading error parameterization for Bz = -2 kG");
146 mRPhiA2 = resRPhiIdeal2, mRPhiB = 1.43e-2f, mRPhiC2 = 4.55e-2f;
147 mDyA2 = 1.225e-3f, mDyB = 9.8e-3f, mDyC2 = 3.88e-2f;
148 mAngleToDyA = 0.1f, mAngleToDyB = 1.89f, mAngleToDyC = 0.4f;
149 }
150 } else if (CAMath::Abs(CAMath::Abs(Bz) - 5) < 0.1f) {
151 // magnetic field +-0.5 T
152 if (Bz > 0) {
153 GPUInfo("Loading error parameterization for Bz = +5 kG");
154 mRPhiA2 = resRPhiIdeal2, mRPhiB = 0.125f, mRPhiC2 = 0.0961f;
155 mDyA2 = 1.681e-3f, mDyB = 0.15f, mDyC2 = 0.1849f;
156 mAngleToDyA = 0.13f, mAngleToDyB = 2.43f, mAngleToDyC = -0.58f;
157 } else {
158 GPUInfo("Loading error parameterization for Bz = -5 kG");
159 mRPhiA2 = resRPhiIdeal2, mRPhiB = -0.14f, mRPhiC2 = 0.1156f;
160 mDyA2 = 2.209e-3f, mDyB = -0.15f, mDyC2 = 0.2025f;
161 mAngleToDyA = -0.15f, mAngleToDyB = 2.34f, mAngleToDyC = 0.56f;
162 }
163 } else {
164 // magnetic field 0 T or another value which is not covered by the error parameterizations
165 // using default values instead
166 GPUWarning("No error parameterization available for Bz = %.2f kG. Keeping default value (sigma_y = const. = 1cm)", Bz);
167 mRPhiA2 = 1.f;
168 }
169
170 // obtain average radius of TRD chambers
171 float x0[kNLayers] = {300.2f, 312.8f, 325.4f, 338.0f, 350.6f, 363.2f}; // used as default value in case no transformation matrix can be obtained
172 auto* matrix = mGeo->GetClusterMatrix(0);
173 float loc[3] = {mGeo->AnodePos(), 0.f, 0.f};
174 float glb[3] = {0.f, 0.f, 0.f};
175 for (int32_t iDet = 0; iDet < kNChambers; ++iDet) {
176 matrix = mGeo->GetClusterMatrix(iDet);
177 if (!matrix) {
178 mR[iDet] = x0[mGeo->GetLayer(iDet)];
179 continue;
180 }
181 matrix->LocalToMaster(loc, glb);
182 mR[iDet] = glb[0];
183 }
184}
185
186template <class TRDTRK, class PROP>
188{
189 //--------------------------------------------------------------------
190 // Reset tracker
191 //--------------------------------------------------------------------
192 mNTracks = 0;
193}
194
195template <class TRDTRK, class PROP>
197{
198 //--------------------------------------------------------------------
199 // Prepare tracklet index array and if requested calculate space points
200 // in part duplicated from DoTracking() method to allow for calling
201 // this function on the host prior to GPU processing
202 //--------------------------------------------------------------------
203 for (uint32_t iColl = 0; iColl < GetConstantMem()->ioPtrs.nTRDTriggerRecords; ++iColl) {
204 if (GetConstantMem()->ioPtrs.trdTrigRecMask && GetConstantMem()->ioPtrs.trdTrigRecMask[iColl] == 0) {
205 // this trigger is masked as there is no ITS information available for it
206 continue;
207 }
208 int32_t nTrklts = 0;
209 int32_t idxOffset = 0;
210 if (mProcessPerTimeFrame) {
211 idxOffset = GetConstantMem()->ioPtrs.trdTrackletIdxFirst[iColl];
212 nTrklts = (iColl < GetConstantMem()->ioPtrs.nTRDTriggerRecords - 1) ? GetConstantMem()->ioPtrs.trdTrackletIdxFirst[iColl + 1] - GetConstantMem()->ioPtrs.trdTrackletIdxFirst[iColl] : GetConstantMem()->ioPtrs.nTRDTracklets - GetConstantMem()->ioPtrs.trdTrackletIdxFirst[iColl];
213 } else {
214 nTrklts = GetConstantMem()->ioPtrs.nTRDTracklets;
215 }
216 const GPUTRDTrackletWord* tracklets = &((GetConstantMem()->ioPtrs.trdTracklets)[idxOffset]);
217 int32_t* trkltIndexArray = &mTrackletIndexArray[iColl * (kNChambers + 1) + 1];
218 trkltIndexArray[-1] = 0;
219 int32_t currDet = 0;
220 int32_t nextDet = 0;
221 int32_t trkltCounter = 0;
222 for (int32_t iTrklt = 0; iTrklt < nTrklts; ++iTrklt) {
223 if (tracklets[iTrklt].GetDetector() > currDet) {
224 nextDet = tracklets[iTrklt].GetDetector();
225 for (int32_t iDet = currDet; iDet < nextDet; ++iDet) {
226 trkltIndexArray[iDet] = trkltCounter;
227 }
228 currDet = nextDet;
229 }
230 ++trkltCounter;
231 }
232 for (int32_t iDet = currDet; iDet <= kNChambers; ++iDet) {
233 trkltIndexArray[iDet] = trkltCounter;
234 }
235 if (mGenerateSpacePoints) {
236 if (!CalculateSpacePoints(iColl)) {
237 GPUError("Space points for at least one chamber could not be calculated (for interaction %i)", iColl);
238 break;
239 }
240 }
241 }
242 if (mGenerateSpacePoints) {
243 chainTracking->mIOPtrs.trdSpacePoints = mSpacePoints;
244 }
245 mNEvents++;
246}
247
248template <class TRDTRK, class PROP>
250{
251 //--------------------------------------------------------------------
252 // set the number of candidates to be used
253 //--------------------------------------------------------------------
254 if (!mIsInitialized) {
255 mNCandidates = n;
256 } else {
257 GPUError("Cannot change mNCandidates after initialization");
258 }
259}
260
261template <class TRDTRK, class PROP>
263{
264 //--------------------------------------------------------------------
265 // print current settings to screen
266 //--------------------------------------------------------------------
267 GPUInfo("##############################################################");
268 GPUInfo("Current settings for GPU TRD tracker:");
269 GPUInfo(" maxChi2(%.2f), chi2Penalty(%.2f), nCandidates(%i), maxMissingLayers(%i)", Param().rec.trd.maxChi2, Param().rec.trd.penaltyChi2, mNCandidates, Param().rec.trd.stopTrkAfterNMissLy);
270 GPUInfo(" ptCut = %.2f GeV, abs(eta) < %.2f", Param().rec.trd.minTrackPt, mMaxEta);
271 GPUInfo("##############################################################");
272}
273
274template <class TRDTRK, class PROP>
276{
277 mDebug->CreateStreamer();
278}
279
280#endif
281
282template <>
283GPUdi() const GPUTRDPropagatorGPU::propagatorParam* GPUTRDTracker_t<GPUTRDTrackGPU, GPUTRDPropagatorGPU>::getPropagatorParam()
284{
285 return &Param().polynomialField;
286}
287
288template <class TRDTRK, class PROP>
289GPUdi() const typename PROP::propagatorParam* GPUTRDTracker_t<TRDTRK, PROP>::getPropagatorParam()
290{
291 return GetConstantMem()->calibObjects.o2Propagator;
292}
293
294template <class TRDTRK, class PROP>
295GPUd() bool GPUTRDTracker_t<TRDTRK, PROP>::CheckTrackTRDCandidate(const TRDTRK& trk) const
296{
297 if (!trk.CheckNumericalQuality()) {
298 return false;
299 }
300 if (CAMath::Abs(trk.getEta()) > mMaxEta) {
301 return false;
302 }
303 if (trk.getPt() < Param().rec.trd.minTrackPt) {
304 return false;
305 }
306 return true;
307}
308
309template <class TRDTRK, class PROP>
310GPUd() int32_t GPUTRDTracker_t<TRDTRK, PROP>::LoadTrack(const TRDTRK& trk, uint32_t tpcTrackId, bool checkTrack, HelperTrackAttributes* attribs)
311{
312 if (mNTracks >= mNMaxTracks) {
313#ifndef GPUCA_GPUCODE
314 GPUError("Error: Track dropped (no memory available) -> must not happen");
315#endif
316 return (1);
317 }
318 if (checkTrack && !CheckTrackTRDCandidate(trk)) {
319 return 2;
320 }
321 mTracks[mNTracks] = trk;
322 mTracks[mNTracks].setRefGlobalTrackIdRaw(tpcTrackId);
323 if (attribs) {
324 mTrackAttribs[mNTracks] = *attribs;
325 }
326 mNTracks++;
327 return (0);
328}
329
330template <class TRDTRK, class PROP>
331GPUd() void GPUTRDTracker_t<TRDTRK, PROP>::DumpTracks()
332{
333 //--------------------------------------------------------------------
334 // helper function (only for debugging purposes)
335 //--------------------------------------------------------------------
336 GPUInfo("There are in total %i tracklets loaded", GetConstantMem()->ioPtrs.nTRDTracklets);
337 GPUInfo("There are %i tracks loaded. mNMaxTracks(%i)", mNTracks, mNMaxTracks);
338 for (int32_t i = 0; i < mNTracks; ++i) {
339 auto* trk = &(mTracks[i]);
340 GPUInfo("track %i: x=%f, alpha=%f, nTracklets=%i, pt=%f, time=%f", i, trk->getX(), trk->getAlpha(), trk->getNtracklets(), trk->getPt(), mTrackAttribs[i].mTime);
341 }
342}
343
344template <class TRDTRK, class PROP>
345GPUd() int32_t GPUTRDTracker_t<TRDTRK, PROP>::GetCollisionIDs(int32_t iTrk, int32_t* collisionIds) const
346{
347 //--------------------------------------------------------------------
348 // Check which TRD trigger times possibly match given input track.
349 // If ITS-TPC matches or CE-crossing TPC tracks the time is precisely
350 // known and max 1 trigger time can be assigned.
351 // For TPC-only tracks the collision IDs are stored in collisionIds array
352 // and the number of valid entries in the array is returned
353 //--------------------------------------------------------------------
354 int32_t nColls = 0;
355 for (uint32_t iColl = 0; iColl < GetConstantMem()->ioPtrs.nTRDTriggerRecords; ++iColl) {
356 if (GetConstantMem()->ioPtrs.trdTrigRecMask && GetConstantMem()->ioPtrs.trdTrigRecMask[iColl] == 0) {
357 continue;
358 }
359 if (GetConstantMem()->ioPtrs.trdTriggerTimes[iColl] > mTrackAttribs[iTrk].GetTimeMin() && GetConstantMem()->ioPtrs.trdTriggerTimes[iColl] < mTrackAttribs[iTrk].GetTimeMax()) {
360 if (nColls == 20) {
361 GPUError("Found too many collision candidates for track with tMin(%f) and tMax(%f)", mTrackAttribs[iTrk].GetTimeMin(), mTrackAttribs[iTrk].GetTimeMax());
362 return nColls;
363 }
364 collisionIds[nColls++] = iColl;
365 }
366 }
367 return nColls;
368}
369
370template <class TRDTRK, class PROP>
371GPUd() void GPUTRDTracker_t<TRDTRK, PROP>::DoTrackingThread(int32_t iTrk, int32_t threadId)
372{
373 //--------------------------------------------------------------------
374 // perform the tracking for one track (must be threadsafe)
375 //--------------------------------------------------------------------
376 int32_t collisionIds[20] = {0}; // due to the dead time there will never exist more possible TRD triggers for a single track
377 int32_t nCollisionIds = 1; // initialize with 1 for AliRoot compatibility
378 if (mProcessPerTimeFrame) {
379 nCollisionIds = GetCollisionIDs(iTrk, collisionIds);
380 if (nCollisionIds == 0) {
381 if (ENABLE_INFO) {
382 GPUInfo("Did not find TRD data for track %i with t=%f. tMin(%f), tMax(%f)", iTrk, mTrackAttribs[iTrk].mTime, mTrackAttribs[iTrk].GetTimeMin(), mTrackAttribs[iTrk].GetTimeMax());
383 }
384 // no TRD data available for the bunch crossing this track originates from
385 return;
386 }
387 }
388 PROP prop(getPropagatorParam());
389 mTracks[iTrk].setChi2(Param().rec.trd.penaltyChi2); // TODO check if this should not be higher
390 auto trkStart = mTracks[iTrk];
391 for (int32_t iColl = 0; iColl < nCollisionIds; ++iColl) {
392 // do track following for each collision candidate and keep best track
393 auto trkCopy = trkStart;
394 prop.setTrack(&trkCopy);
395 prop.setFitInProjections(true);
396 if (!FollowProlongation(&prop, &trkCopy, iTrk, threadId, collisionIds[iColl])) {
397 // track following failed
398 continue;
399 }
400 if (trkCopy.getReducedChi2() < mTracks[iTrk].getReducedChi2()) {
401 mTracks[iTrk] = trkCopy; // copy back the resulting track
402 }
403 }
404}
405
406template <class TRDTRK, class PROP>
407GPUd() bool GPUTRDTracker_t<TRDTRK, PROP>::CalculateSpacePoints(int32_t iCollision)
408{
409 //--------------------------------------------------------------------
410 // Calculates TRD space points in sector tracking coordinates
411 // from online tracklets
412 //--------------------------------------------------------------------
413
414 bool result = true;
415 int32_t idxOffset = iCollision * (kNChambers + 1); // offset for accessing mTrackletIndexArray for collision iCollision
416
417 const GPUTRDTrackletWord* tracklets = GetConstantMem()->ioPtrs.trdTracklets;
418
419 for (int32_t iDet = 0; iDet < kNChambers; ++iDet) {
420 int32_t iFirstTrackletInDet = mTrackletIndexArray[idxOffset + iDet];
421 int32_t iFirstTrackletInNextDet = mTrackletIndexArray[idxOffset + iDet + 1];
422 int32_t nTrackletsInDet = iFirstTrackletInNextDet - iFirstTrackletInDet;
423 if (nTrackletsInDet == 0) {
424 continue;
425 }
426 if (!mGeo->ChamberInGeometry(iDet)) {
427 GPUError("Found TRD tracklets in chamber %i which is not included in the geometry", iDet);
428 return false;
429 }
430 auto* matrix = mGeo->GetClusterMatrix(iDet);
431 if (!matrix) {
432 GPUError("No cluster matrix available for chamber %i. Skipping it...", iDet);
433 result = false;
434 continue;
435 }
436 const GPUTRDpadPlane* pp = mGeo->GetPadPlane(iDet);
437
438 int32_t trkltIdxOffset = (mProcessPerTimeFrame) ? GetConstantMem()->ioPtrs.trdTrackletIdxFirst[iCollision] : 0; // global index of first tracklet in iCollision
439 int32_t trkltIdxStart = trkltIdxOffset + iFirstTrackletInDet;
440 for (int32_t trkltIdx = trkltIdxStart; trkltIdx < trkltIdxStart + nTrackletsInDet; ++trkltIdx) {
441 int32_t trkltZbin = tracklets[trkltIdx].GetZbin();
442 float xTrkltDet[3] = {0.f}; // trklt position in chamber coordinates
443 float xTrkltSec[3] = {0.f}; // trklt position in sector coordinates
444 xTrkltDet[0] = mGeo->AnodePos() + sRadialOffset;
445 xTrkltDet[1] = tracklets[trkltIdx].GetY();
446 xTrkltDet[2] = pp->GetRowPos(trkltZbin) - pp->GetRowSize(trkltZbin) / 2.f - pp->GetRowPos(pp->GetNrows() / 2);
447 // GPUInfo("Space point local %i: x=%f, y=%f, z=%f", trkltIdx, xTrkltDet[0], xTrkltDet[1], xTrkltDet[2]);
448 matrix->LocalToMaster(xTrkltDet, xTrkltSec);
449 mSpacePoints[trkltIdx].setX(xTrkltSec[0]);
450 mSpacePoints[trkltIdx].setY(xTrkltSec[1]);
451 mSpacePoints[trkltIdx].setZ(xTrkltSec[2]);
452 mSpacePoints[trkltIdx].setDy(tracklets[trkltIdx].GetdY());
453
454 // GPUInfo("Space point global %i: x=%f, y=%f, z=%f", trkltIdx, mSpacePoints[trkltIdx].getX(), mSpacePoints[trkltIdx].getY(), mSpacePoints[trkltIdx].getZ());
455 }
456 }
457 return result;
458}
459
460template <class TRDTRK, class PROP>
461GPUd() bool GPUTRDTracker_t<TRDTRK, PROP>::FollowProlongation(PROP* prop, TRDTRK* t, int32_t iTrk, int32_t threadId, int32_t collisionId)
462{
463 //--------------------------------------------------------------------
464 // Propagate TPC track layerwise through TRD and pick up closest
465 // tracklet(s) on the way
466 // -> returns false if prolongation could not be executed fully
467 // or track does not fullfill threshold conditions
468 //--------------------------------------------------------------------
469
470 if (ENABLE_INFO) {
471 GPUInfo("Start track following for track %i at x=%f with pt=%f", t->getRefGlobalTrackIdRaw(), t->getX(), t->getPt());
472 }
473 mDebug->Reset();
474 t->setChi2(0.f);
475 float zShiftTrk = 0.f;
476 if (mProcessPerTimeFrame) {
477 zShiftTrk = (mTrackAttribs[iTrk].mTime - GetConstantMem()->ioPtrs.trdTriggerTimes[collisionId]) * mTPCVdrift * mTrackAttribs[iTrk].mSide;
478 // float addZerr = (mTrackAttribs[iTrk].mTimeAddMax + mTrackAttribs[iTrk].mTimeSubMax) * .5f * mTPCVdrift;
479 // increase Z error based on time window
480 // -> this is here since it was done before, but the efficiency seems to be better if the covariance is not updated (more tracklets are attached)
481 // t->updateCovZ2(addZerr * addZerr); // TODO check again once detailed performance study tools are available, maybe this can be tuned
482 }
483 const GPUTRDpadPlane* pad = nullptr;
484 const GPUTRDTrackletWord* tracklets = GetConstantMem()->ioPtrs.trdTracklets;
485 const GPUTRDSpacePoint* spacePoints = GetConstantMem()->ioPtrs.trdSpacePoints;
486
487#ifdef ENABLE_GPUTRDDEBUG
488 TRDTRK trackNoUp(*t);
489#endif
490
491 int32_t candidateIdxOffset = threadId * 2 * mNCandidates;
492 int32_t hypothesisIdxOffset = threadId * mNCandidates;
493 int32_t trkltIdxOffset = collisionId * (kNChambers + 1); // offset for accessing mTrackletIndexArray for given collision
494 int32_t glbTrkltIdxOffset = (mProcessPerTimeFrame) ? GetConstantMem()->ioPtrs.trdTrackletIdxFirst[collisionId] : 0; // offset of first tracklet in given collision in global tracklet array
495
496 auto trkWork = t;
497 if (mNCandidates > 1) {
498 // copy input track to first candidate
499 mCandidates[candidateIdxOffset] = *t;
500 }
501
502 int32_t nCandidates = 1; // we always start with one candidate
503 int32_t nCurrHypothesis = 0; // the number of track hypothesis in given iLayer
504
505 // search window
506 float roadY = 0.f;
507 float roadZ = 0.f;
508 const int32_t nMaxChambersToSearch = 4;
509
510 mDebug->SetGeneralInfo(mNEvents, mNTracks, iTrk, t->getPt());
511
512 for (int32_t iLayer = 0; iLayer < kNLayers; ++iLayer) {
513 nCurrHypothesis = 0;
514 bool isOK = false; // if at least one candidate could be propagated or the track was stopped this becomes true
515 int32_t currIdx = candidateIdxOffset + iLayer % 2;
516 int32_t nextIdx = candidateIdxOffset + (iLayer + 1) % 2;
517 pad = mGeo->GetPadPlane(iLayer, 0);
518 float tilt = CAMath::Tan(CAMath::Pi() / 180.f * pad->GetTiltingAngle()); // tilt is signed!
519 const float zMaxTRD = pad->GetRow0();
520
521 // --------------------------------------------------------------------------------
522 //
523 // for each candidate, propagate to layer radius and look for close tracklets
524 //
525 // --------------------------------------------------------------------------------
526 for (int32_t iCandidate = 0; iCandidate < nCandidates; iCandidate++) {
527
528 int32_t det[nMaxChambersToSearch] = {-1, -1, -1, -1}; // TRD chambers to be searched for tracklets
529
530 if (mNCandidates > 1) {
531 trkWork = &mCandidates[2 * iCandidate + currIdx];
532 prop->setTrack(trkWork);
533 }
534
535 if (trkWork->getIsStopped()) {
536 Hypothesis hypo(trkWork->getNlayersFindable(), iCandidate, -1, trkWork->getChi2());
537 InsertHypothesis(hypo, nCurrHypothesis, hypothesisIdxOffset);
538 isOK = true;
539 continue;
540 }
541
542 // propagate track to average radius of TRD layer iLayer (sector 0, stack 2 is chosen as a reference)
543 if (!prop->propagateToX(mR[2 * kNLayers + iLayer], .8f, 2.f)) {
544 if (ENABLE_INFO) {
545 GPUInfo("Track propagation failed for track %i candidate %i in layer %i (pt=%f, x=%f, mR[layer]=%f)", iTrk, iCandidate, iLayer, trkWork->getPt(), trkWork->getX(), mR[2 * kNLayers + iLayer]);
546 }
547 continue;
548 }
549
550 // rotate track in new sector in case of sector crossing
551 if (!AdjustSector(prop, trkWork)) {
552 if (ENABLE_INFO) {
553 GPUInfo("Adjusting sector failed for track %i candidate %i in layer %i", iTrk, iCandidate, iLayer);
554 }
555 continue;
556 }
557
558 // check if track is findable
559 if (IsGeoFindable(trkWork, iLayer, prop->getAlpha(), zShiftTrk)) {
560 trkWork->setIsFindable(iLayer);
561 }
562
563 // define search window
564 roadY = 7.f * CAMath::Sqrt(trkWork->getSigmaY2() + 0.1f * 0.1f) + Param().rec.trd.extraRoadY; // add constant to the road to account for uncertainty due to radial deviations (few mm)
565 // roadZ = 7.f * CAMath::Sqrt(trkWork->getSigmaZ2() + 9.f * 9.f / 12.f); // take longest pad length
566 roadZ = mRoadZ + Param().rec.trd.extraRoadZ; // simply twice the longest pad length -> efficiency 99.996%
567 //
568 if (CAMath::Abs(trkWork->getZ() + zShiftTrk) - roadZ >= zMaxTRD) {
569 if (ENABLE_INFO) {
570 GPUInfo("Track out of TRD acceptance with z=%f in layer %i (eta=%f)", trkWork->getZ() + zShiftTrk, iLayer, trkWork->getEta());
571 }
572 continue;
573 }
574
575 // determine chamber(s) to be searched for tracklets
576 FindChambersInRoad(trkWork, roadY, roadZ, iLayer, det, zMaxTRD, prop->getAlpha(), zShiftTrk);
577
578 // track debug information to be stored in case no matching tracklet can be found
579 mDebug->SetTrackParameter(*trkWork, iLayer);
580
581 // look for tracklets in chamber(s)
582 for (int32_t iDet = 0; iDet < nMaxChambersToSearch; iDet++) {
583 int32_t currDet = det[iDet];
584 if (currDet == -1) {
585 continue;
586 }
587 pad = mGeo->GetPadPlane(currDet);
588 int32_t currSec = mGeo->GetSector(currDet);
589 if (currSec != GetSector(prop->getAlpha())) {
590 if (!prop->rotate(GetAlphaOfSector(currSec))) {
591 if (ENABLE_WARNING) {
592 GPUWarning("Track could not be rotated in tracklet coordinate system");
593 }
594 break;
595 }
596 }
597 if (currSec != GetSector(prop->getAlpha())) {
598 GPUError("Track is in sector %i and sector %i is searched for tracklets", GetSector(prop->getAlpha()), currSec);
599 continue;
600 }
601 // propagate track to radius of chamber
602 if (!prop->propagateToX(mR[currDet], .8f, .2f)) {
603 if (ENABLE_WARNING) {
604 GPUWarning("Track parameter for track %i, x=%f at chamber %i x=%f in layer %i cannot be retrieved", iTrk, trkWork->getX(), currDet, mR[currDet], iLayer);
605 }
606 }
607 // first propagate track to x of tracklet
608 for (int32_t trkltIdx = glbTrkltIdxOffset + mTrackletIndexArray[trkltIdxOffset + currDet]; trkltIdx < glbTrkltIdxOffset + mTrackletIndexArray[trkltIdxOffset + currDet + 1]; ++trkltIdx) {
609 if (CAMath::Abs(trkWork->getY() - spacePoints[trkltIdx].getY()) > roadY || CAMath::Abs(trkWork->getZ() + zShiftTrk - spacePoints[trkltIdx].getZ()) > roadZ) {
610 // skip tracklets which are too far away
611 // although the radii of space points and tracks may differ by ~ few mm the roads are large enough to allow no efficiency loss by this cut
612 continue;
613 }
614 float projY, projZ;
615 prop->getPropagatedYZ(spacePoints[trkltIdx].getX(), projY, projZ);
616 // correction for tilted pads (only applied if deltaZ < lPad && track z err << lPad)
617 float tiltCorr = tilt * (spacePoints[trkltIdx].getZ() - projZ);
618 float lPad = pad->GetRowSize(tracklets[trkltIdx].GetZbin());
619 if (!((CAMath::Abs(spacePoints[trkltIdx].getZ() - projZ) < lPad) && (trkWork->getSigmaZ2() < (lPad * lPad / 12.f)))) {
620 tiltCorr = 0.f; // will be zero also for TPC tracks which are shifted in z
621 }
622 // correction for mean z position of tracklet (is not the center of the pad if track eta != 0)
623 float zPosCorr = spacePoints[trkltIdx].getZ() + mZCorrCoefNRC * trkWork->getTgl();
624 float yPosCorr = spacePoints[trkltIdx].getY() - tiltCorr;
625 zPosCorr -= zShiftTrk; // shift tracklet instead of track in order to avoid having to do a re-fit for each collision
626 float deltaY = yPosCorr - projY;
627 float deltaZ = zPosCorr - projZ;
628 float trkltPosTmpYZ[2] = {yPosCorr, zPosCorr};
629 float trkltCovTmp[3] = {0.f};
630 if ((CAMath::Abs(deltaY) < roadY) && (CAMath::Abs(deltaZ) < roadZ)) { // TODO: check if this is still necessary after the cut before propagation of track
631 // tracklet is in windwow: get predicted chi2 for update and store tracklet index if best guess
632 RecalcTrkltCov(tilt, trkWork->getSnp(), pad->GetRowSize(tracklets[trkltIdx].GetZbin()), trkltCovTmp);
633 float chi2 = prop->getPredictedChi2(trkltPosTmpYZ, trkltCovTmp);
634 // TODO cut on angular pull should be made stricter when proper v-drift calibration for the TRD tracklets is implemented
635 if ((chi2 > Param().rec.trd.maxChi2) || (Param().rec.trd.applyDeflectionCut && CAMath::Abs(GetAngularPull(spacePoints[trkltIdx].getDy(), trkWork->getSnp())) > 4)) {
636 continue;
637 }
638 Hypothesis hypo(trkWork->getNlayersFindable(), iCandidate, trkltIdx, trkWork->getChi2() + chi2);
639 InsertHypothesis(hypo, nCurrHypothesis, hypothesisIdxOffset);
640 } // end tracklet in window
641 } // tracklet loop
642 } // chamber loop
643
644 // add no update to hypothesis list
645 Hypothesis hypoNoUpdate(trkWork->getNlayersFindable(), iCandidate, -1, trkWork->getChi2() + Param().rec.trd.penaltyChi2);
646 InsertHypothesis(hypoNoUpdate, nCurrHypothesis, hypothesisIdxOffset);
647 isOK = true;
648 } // end candidate loop
649
650 mDebug->SetChi2Update(mHypothesis[0 + hypothesisIdxOffset].mChi2 - t->getChi2(), iLayer); // only meaningful for ONE candidate!!!
651 mDebug->SetRoad(roadY, roadZ, iLayer); // only meaningful for ONE candidate
652 bool wasTrackStored = false;
653 // --------------------------------------------------------------------------------
654 //
655 // loop over the best N_candidates hypothesis
656 //
657 // --------------------------------------------------------------------------------
658 // GPUInfo("nCurrHypothesis=%i, nCandidates=%i", nCurrHypothesis, nCandidates);
659 // for (int32_t idx=0; idx<10; ++idx) { GPUInfo("mHypothesis[%i]: candidateId=%i, nLayers=%i, trackletId=%i, chi2=%f", idx, mHypothesis[idx].mCandidateId, mHypothesis[idx].mLayers, mHypothesis[idx].mTrackletId, mHypothesis[idx].mChi2); }
660 for (int32_t iUpdate = 0; iUpdate < nCurrHypothesis && iUpdate < mNCandidates; iUpdate++) {
661 if (mHypothesis[iUpdate + hypothesisIdxOffset].mCandidateId == -1) {
662 // no more candidates
663 if (iUpdate == 0) {
664 return false; // no valid candidates for this track (probably propagation failed)
665 }
666 break; // go to next layer
667 }
668 nCandidates = iUpdate + 1;
669 if (mNCandidates > 1) {
670 mCandidates[2 * iUpdate + nextIdx] = mCandidates[2 * mHypothesis[iUpdate + hypothesisIdxOffset].mCandidateId + currIdx];
671 trkWork = &mCandidates[2 * iUpdate + nextIdx];
672 }
673 if (mHypothesis[iUpdate + hypothesisIdxOffset].mTrackletId == -1) {
674 // no matching tracklet found
675 if (trkWork->getIsFindable(iLayer)) {
676 if (trkWork->getNmissingConsecLayers(iLayer) > Param().rec.trd.stopTrkAfterNMissLy) {
677 trkWork->setIsStopped();
678 }
679 trkWork->setChi2(trkWork->getChi2() + Param().rec.trd.penaltyChi2);
680 }
681 if (iUpdate == 0 && mNCandidates > 1) { // TODO: is thie really necessary????? CHECK!
682 *t = mCandidates[2 * iUpdate + nextIdx];
683 }
684 continue;
685 }
686 // matching tracklet found
687 if (mNCandidates > 1) {
688 prop->setTrack(trkWork);
689 }
690 int32_t trkltSec = mGeo->GetSector(tracklets[mHypothesis[iUpdate + hypothesisIdxOffset].mTrackletId].GetDetector());
691 if (trkltSec != GetSector(prop->getAlpha())) {
692 // if after a matching tracklet was found another sector was searched for tracklets the track needs to be rotated back
693 prop->rotate(GetAlphaOfSector(trkltSec));
694 }
695 if (!prop->propagateToX(spacePoints[mHypothesis[iUpdate + hypothesisIdxOffset].mTrackletId].getX(), .8f, 2.f)) {
696 if (ENABLE_WARNING) {
697 GPUWarning("Final track propagation for track %i update %i in layer %i failed", iTrk, iUpdate, iLayer);
698 }
699 trkWork->setChi2(trkWork->getChi2() + Param().rec.trd.penaltyChi2);
700 if (trkWork->getIsFindable(iLayer)) {
701 if (trkWork->getNmissingConsecLayers(iLayer) >= Param().rec.trd.stopTrkAfterNMissLy) {
702 trkWork->setIsStopped();
703 }
704 }
705 if (iUpdate == 0 && mNCandidates > 1) {
706 *t = mCandidates[2 * iUpdate + nextIdx];
707 }
708 continue;
709 }
710
711 pad = mGeo->GetPadPlane(tracklets[mHypothesis[iUpdate + hypothesisIdxOffset].mTrackletId].GetDetector());
712 float tiltCorrUp = tilt * (spacePoints[mHypothesis[iUpdate + hypothesisIdxOffset].mTrackletId].getZ() - trkWork->getZ());
713 float zPosCorrUp = spacePoints[mHypothesis[iUpdate + hypothesisIdxOffset].mTrackletId].getZ() + mZCorrCoefNRC * trkWork->getTgl();
714 zPosCorrUp -= zShiftTrk;
715 float padLength = pad->GetRowSize(tracklets[mHypothesis[iUpdate + hypothesisIdxOffset].mTrackletId].GetZbin());
716 if (!((trkWork->getSigmaZ2() < (padLength * padLength / 12.f)) && (CAMath::Abs(spacePoints[mHypothesis[iUpdate + hypothesisIdxOffset].mTrackletId].getZ() - trkWork->getZ()) < padLength))) {
717 tiltCorrUp = 0.f;
718 }
719 float trkltPosUp[2] = {spacePoints[mHypothesis[iUpdate + hypothesisIdxOffset].mTrackletId].getY() - tiltCorrUp, zPosCorrUp};
720 float trkltCovUp[3] = {0.f};
721 RecalcTrkltCov(tilt, trkWork->getSnp(), pad->GetRowSize(tracklets[mHypothesis[iUpdate + hypothesisIdxOffset].mTrackletId].GetZbin()), trkltCovUp);
722
723#ifdef ENABLE_GPUTRDDEBUG
724 prop->setTrack(&trackNoUp);
725 prop->rotate(GetAlphaOfSector(trkltSec));
726 // prop->propagateToX(spacePoints[mHypothesis[iUpdate + hypothesisIdxOffset].mTrackletId].getX(), .8f, 2.f);
727 prop->propagateToX(mR[tracklets[mHypothesis[iUpdate + hypothesisIdxOffset].mTrackletId].GetDetector()], .8f, 2.f);
728 prop->setTrack(trkWork);
729#endif
730
731 if (!wasTrackStored) {
732#ifdef ENABLE_GPUTRDDEBUG
733 mDebug->SetTrackParameterNoUp(trackNoUp, iLayer);
734#endif
735 mDebug->SetTrackParameter(*trkWork, iLayer);
736 mDebug->SetRawTrackletPosition(spacePoints[mHypothesis[iUpdate + hypothesisIdxOffset].mTrackletId].getX(), spacePoints[mHypothesis[iUpdate + hypothesisIdxOffset].mTrackletId].getY(), spacePoints[mHypothesis[iUpdate + hypothesisIdxOffset].mTrackletId].getZ(), iLayer);
737 mDebug->SetCorrectedTrackletPosition(trkltPosUp, iLayer);
738 mDebug->SetTrackletCovariance(trkltCovUp, iLayer);
739 mDebug->SetTrackletProperties(spacePoints[mHypothesis[iUpdate + hypothesisIdxOffset].mTrackletId].getDy(), tracklets[mHypothesis[iUpdate + hypothesisIdxOffset].mTrackletId].GetDetector(), iLayer);
740 wasTrackStored = true;
741 }
742
743 if (!prop->update(trkltPosUp, trkltCovUp)) {
744 if (ENABLE_WARNING) {
745 GPUWarning("Failed to update track %i with space point in layer %i", iTrk, iLayer);
746 }
747 trkWork->setChi2(trkWork->getChi2() + Param().rec.trd.penaltyChi2);
748 if (trkWork->getIsFindable(iLayer)) {
749 if (trkWork->getNmissingConsecLayers(iLayer) >= Param().rec.trd.stopTrkAfterNMissLy) {
750 trkWork->setIsStopped();
751 }
752 }
753 if (iUpdate == 0 && mNCandidates > 1) {
754 *t = mCandidates[2 * iUpdate + nextIdx];
755 }
756 continue;
757 }
758 if (!trkWork->CheckNumericalQuality()) {
759 if (ENABLE_INFO) {
760 GPUInfo("Track %i has invalid covariance matrix. Aborting track following\n", iTrk);
761 }
762 return false;
763 }
764 trkWork->addTracklet(iLayer, mHypothesis[iUpdate + hypothesisIdxOffset].mTrackletId);
765 trkWork->setChi2(mHypothesis[iUpdate + hypothesisIdxOffset].mChi2);
766 trkWork->setIsFindable(iLayer);
767 trkWork->setCollisionId(collisionId);
768 // check if track crosses padrows
769 float projZEntry, projYEntry;
770 // Get Z for Entry of Track
771 prop->getPropagatedYZ(trkWork->getX() - mGeo->GetCdrHght(), projYEntry, projZEntry);
772 // Get Padrow number for Exit&Entry and compare. If is not equal mark
773 // as padrow crossing. While simple, this comes with the cost of not
774 // properly handling Tracklets that hit a different Pad but still get
775 // attached to the Track. It is estimated that the probability for
776 // this is low.
777 const auto padrowEntry = pad->GetPadRowNumber(projZEntry);
778 const auto padrowExit = pad->GetPadRowNumber(trkWork->getZ());
779 if (padrowEntry != padrowExit) {
780 trkWork->setIsCrossingNeighbor(iLayer);
781 trkWork->setHasPadrowCrossing();
782 }
783 const auto currDet = tracklets[mHypothesis[iUpdate + hypothesisIdxOffset].mTrackletId].GetDetector();
784 // Mark tracklets as Padrow crossing if they have a neighboring tracklet.
785 for (int32_t trkltIdx = glbTrkltIdxOffset + mTrackletIndexArray[trkltIdxOffset + currDet]; trkltIdx < glbTrkltIdxOffset + mTrackletIndexArray[trkltIdxOffset + currDet + 1]; ++trkltIdx) {
786 // skip orig tracklet
787 if (mHypothesis[iUpdate + hypothesisIdxOffset].mTrackletId == trkltIdx) {
788 continue;
789 }
790 if (CAMath::Abs(tracklets[mHypothesis[iUpdate + hypothesisIdxOffset].mTrackletId].GetZbin() - tracklets[trkltIdx].GetZbin()) == 1 &&
791 CAMath::Abs(tracklets[mHypothesis[iUpdate + hypothesisIdxOffset].mTrackletId].GetY() - tracklets[trkltIdx].GetY()) < 1) {
792 trkWork->setIsCrossingNeighbor(iLayer);
793 trkWork->setHasNeighbor();
794 break;
795 }
796 }
797 if (iUpdate == 0 && mNCandidates > 1) {
798 *t = mCandidates[2 * iUpdate + nextIdx];
799 }
800 } // end update loop
801
802 if (!isOK) {
803 if (ENABLE_INFO) {
804 GPUInfo("Track %i cannot be followed. Stopped in layer %i", iTrk, iLayer);
805 }
806 return false;
807 }
808 } // end layer loop
809
810 // --------------------------------------------------------------------------------
811 // add some debug information (compare labels of attached tracklets to track label)
812 // and store full track information
813 // --------------------------------------------------------------------------------
814 if (mDebugOutput) {
815 mDebug->SetTrack(*t);
816 mDebug->Output();
817 }
818 if (ENABLE_INFO) {
819 GPUInfo("Ended track following for track %i at x=%f with pt=%f. Attached %i tracklets", t->getRefGlobalTrackIdRaw(), t->getX(), t->getPt(), t->getNtracklets());
820 }
821 if (nCurrHypothesis > 1) {
822 if (CAMath::Abs(mHypothesis[hypothesisIdxOffset + 1].GetReducedChi2() - mHypothesis[hypothesisIdxOffset].GetReducedChi2()) < Param().rec.trd.chi2SeparationCut) {
823 t->setIsAmbiguous();
824 }
825 }
826 return true;
827}
828
829template <class TRDTRK, class PROP>
830GPUd() void GPUTRDTracker_t<TRDTRK, PROP>::InsertHypothesis(Hypothesis hypo, int32_t& nCurrHypothesis, int32_t idxOffset)
831{
832 // Insert hypothesis into the array. If the array is full and the reduced chi2 is worse
833 // than the worst hypothesis in the array it is dropped.
834 // The hypothesis array is always sorted.
835
836 if (nCurrHypothesis == 0) {
837 // this is the first hypothesis in the array
838 mHypothesis[idxOffset] = hypo;
839 ++nCurrHypothesis;
840 } else if (nCurrHypothesis > 0 && nCurrHypothesis < mNCandidates) {
841 // insert the hypothesis into the right position and shift all worse hypothesis to the right
842 for (int32_t i = idxOffset; i < nCurrHypothesis + idxOffset; ++i) {
843 if (hypo.GetReducedChi2() < mHypothesis[i].GetReducedChi2()) {
844 for (int32_t k = nCurrHypothesis + idxOffset; k > i; --k) {
845 mHypothesis[k] = mHypothesis[k - 1];
846 }
847 mHypothesis[i] = hypo;
848 ++nCurrHypothesis;
849 return;
850 }
851 }
852 mHypothesis[nCurrHypothesis + idxOffset] = hypo;
853 ++nCurrHypothesis;
854 return;
855 } else {
856 // array is already full, check if new hypothesis should be inserted
857 int32_t i = nCurrHypothesis + idxOffset - 1;
858 for (; i >= idxOffset; --i) {
859 if (mHypothesis[i].GetReducedChi2() < hypo.GetReducedChi2()) {
860 break;
861 }
862 }
863 if (i < (nCurrHypothesis + idxOffset - 1)) {
864 // new hypothesis should be inserted into the array
865 for (int32_t k = nCurrHypothesis + idxOffset - 1; k > i + 1; --k) {
866 mHypothesis[k] = mHypothesis[k - 1];
867 }
868 mHypothesis[i + 1] = hypo;
869 }
870 }
871}
872
873template <class TRDTRK, class PROP>
874GPUd() int32_t GPUTRDTracker_t<TRDTRK, PROP>::GetDetectorNumber(const float zPos, const float alpha, const int32_t layer) const
875{
876 //--------------------------------------------------------------------
877 // if track position is within chamber return the chamber number
878 // otherwise return -1
879 //--------------------------------------------------------------------
880 int32_t stack = mGeo->GetStack(zPos, layer);
881 if (stack < 0) {
882 return -1;
883 }
884 int32_t sector = GetSector(alpha);
885
886 return mGeo->GetDetector(layer, stack, sector);
887}
888
889template <class TRDTRK, class PROP>
890GPUd() bool GPUTRDTracker_t<TRDTRK, PROP>::AdjustSector(PROP* prop, TRDTRK* t) const
891{
892 //--------------------------------------------------------------------
893 // rotate track in new sector if necessary and
894 // propagate to previous x afterwards
895 // cancel if track crosses two sector boundaries
896 //--------------------------------------------------------------------
897 float alpha = mGeo->GetAlpha();
898 float xTmp = t->getX();
899 float y = t->getY();
900 float yMax = t->getX() * CAMath::Tan(0.5f * alpha);
901 float alphaCurr = t->getAlpha();
902
903 if (CAMath::Abs(y) > 2.f * yMax) {
904 if (ENABLE_INFO) {
905 GPUInfo("AdjustSector: Track %i with pT = %f crossing two sector boundaries at x = %f", t->getRefGlobalTrackIdRaw(), t->getPt(), t->getX());
906 }
907 return false;
908 }
909
910 int32_t nTries = 0;
911 while (CAMath::Abs(y) > yMax) {
912 if (nTries >= 2) {
913 return false;
914 }
915 int32_t sign = (y > 0) ? 1 : -1;
916 float alphaNew = alphaCurr + alpha * sign;
917 if (alphaNew > CAMath::Pi()) {
918 alphaNew -= 2 * CAMath::Pi();
919 } else if (alphaNew < -CAMath::Pi()) {
920 alphaNew += 2 * CAMath::Pi();
921 }
922 if (!prop->rotate(alphaNew)) {
923 return false;
924 }
925 if (!prop->propagateToX(xTmp, .8f, 2.f)) {
926 return false;
927 }
928 y = t->getY();
929 ++nTries;
930 }
931 return true;
932}
933
934template <class TRDTRK, class PROP>
935GPUd() int32_t GPUTRDTracker_t<TRDTRK, PROP>::GetSector(float alpha) const
936{
937 //--------------------------------------------------------------------
938 // TRD sector number for reference system alpha
939 //--------------------------------------------------------------------
940 if (alpha < 0) {
941 alpha += 2.f * CAMath::Pi();
942 } else if (alpha >= 2.f * CAMath::Pi()) {
943 alpha -= 2.f * CAMath::Pi();
944 }
945 return (int32_t)(alpha * (float)kNSectors / (2.f * CAMath::Pi()));
946}
947
948template <class TRDTRK, class PROP>
949GPUd() float GPUTRDTracker_t<TRDTRK, PROP>::GetAlphaOfSector(const int32_t sec) const
950{
951 //--------------------------------------------------------------------
952 // rotation angle for TRD sector sec
953 //--------------------------------------------------------------------
954 float alpha = 2.0f * CAMath::Pi() / (float)kNSectors * ((float)sec + 0.5f);
955 if (alpha > CAMath::Pi()) {
956 alpha -= 2 * CAMath::Pi();
957 }
958 return alpha;
959}
960
961template <class TRDTRK, class PROP>
962GPUd() void GPUTRDTracker_t<TRDTRK, PROP>::RecalcTrkltCov(const float tilt, const float snp, const float rowSize, float (&cov)[3])
963{
964 //--------------------------------------------------------------------
965 // recalculate tracklet covariance taking track phi angle into account
966 // store the new values in cov
967 //--------------------------------------------------------------------
968 float t2 = tilt * tilt; // tan^2 (tilt)
969 float c2 = 1.f / (1.f + t2); // cos^2 (tilt)
970 float sy2 = GetRPhiRes(snp);
971 float sz2 = rowSize * rowSize / 12.f;
972 cov[0] = c2 * (sy2 + t2 * sz2);
973 cov[1] = c2 * tilt * (sz2 - sy2);
974 cov[2] = c2 * (t2 * sy2 + sz2);
975}
976
977template <class TRDTRK, class PROP>
978GPUd() float GPUTRDTracker_t<TRDTRK, PROP>::GetAngularPull(float dYtracklet, float snp) const
979{
980 float dYtrack = ConvertAngleToDy(snp);
981 float dYresolution = GetAngularResolution(snp);
982 if (dYresolution < 1e-6f) {
983 return 999.f;
984 }
985 return (dYtracklet - dYtrack) / CAMath::Sqrt(dYresolution);
986}
987
988template <class TRDTRK, class PROP>
989GPUd() void GPUTRDTracker_t<TRDTRK, PROP>::FindChambersInRoad(const TRDTRK* t, const float roadY, const float roadZ, const int32_t iLayer, int32_t* det, const float zMax, const float alpha, const float zShiftTrk) const
990{
991 //--------------------------------------------------------------------
992 // determine initial chamber where the track ends up
993 // add more chambers of the same sector or (and) neighbouring
994 // stack if track is close the edge(s) of the chamber
995 //--------------------------------------------------------------------
996
997 const float yMax = CAMath::Abs(mGeo->GetCol0(iLayer));
998 float zTrk = t->getZ() + zShiftTrk;
999
1000 int32_t currStack = mGeo->GetStack(zTrk, iLayer);
1001 int32_t currSec = GetSector(alpha);
1002 int32_t currDet;
1003
1004 int32_t nDets = 0;
1005
1006 if (currStack > -1) {
1007 // chamber unambiguous
1008 currDet = mGeo->GetDetector(iLayer, currStack, currSec);
1009 det[nDets++] = currDet;
1010 const GPUTRDpadPlane* pp = mGeo->GetPadPlane(iLayer, currStack);
1011 int32_t lastPadRow = mGeo->GetRowMax(iLayer, currStack, 0);
1012 float zCenter = pp->GetRowPos(lastPadRow / 2);
1013 if ((zTrk + roadZ) > pp->GetRow0() || (zTrk - roadZ) < pp->GetRowEnd()) {
1014 int32_t addStack = zTrk > zCenter ? currStack - 1 : currStack + 1;
1015 if (addStack < kNStacks && addStack > -1) {
1016 det[nDets++] = mGeo->GetDetector(iLayer, addStack, currSec);
1017 }
1018 }
1019 } else {
1020 if (CAMath::Abs(zTrk) > zMax) {
1021 // shift track in z so it is in the TRD acceptance
1022 if (zTrk > 0) {
1023 currDet = mGeo->GetDetector(iLayer, 0, currSec);
1024 } else {
1025 currDet = mGeo->GetDetector(iLayer, kNStacks - 1, currSec);
1026 }
1027 det[nDets++] = currDet;
1028 currStack = mGeo->GetStack(currDet);
1029 } else {
1030 // track in between two stacks, add both surrounding chambers
1031 // gap between two stacks is 4 cm wide
1032 currDet = GetDetectorNumber(zTrk + 4.0f, alpha, iLayer);
1033 if (currDet != -1) {
1034 det[nDets++] = currDet;
1035 }
1036 currDet = GetDetectorNumber(zTrk - 4.0f, alpha, iLayer);
1037 if (currDet != -1) {
1038 det[nDets++] = currDet;
1039 }
1040 }
1041 }
1042 // add chamber(s) from neighbouring sector in case the track is close to the boundary
1043 if ((CAMath::Abs(t->getY()) + roadY) > yMax) {
1044 const int32_t nStacksToSearch = nDets;
1045 int32_t newSec;
1046 if (t->getY() > 0) {
1047 newSec = (currSec + 1) % kNSectors;
1048 } else {
1049 newSec = (currSec > 0) ? currSec - 1 : kNSectors - 1;
1050 }
1051 for (int32_t idx = 0; idx < nStacksToSearch; ++idx) {
1052 currStack = mGeo->GetStack(det[idx]);
1053 det[nDets++] = mGeo->GetDetector(iLayer, currStack, newSec);
1054 }
1055 }
1056 // skip PHOS hole and non-existing chamber 17_4_4
1057 for (int32_t iDet = 0; iDet < nDets; iDet++) {
1058 if (!mGeo->ChamberInGeometry(det[iDet])) {
1059 det[iDet] = -1;
1060 }
1061 }
1062}
1063
1064template <class TRDTRK, class PROP>
1065GPUd() bool GPUTRDTracker_t<TRDTRK, PROP>::IsGeoFindable(const TRDTRK* t, const int32_t layer, const float alpha, const float zShiftTrk) const
1066{
1067 //--------------------------------------------------------------------
1068 // returns true if track position inside active area of the TRD
1069 // and not too close to the boundaries
1070 //--------------------------------------------------------------------
1071
1072 float zTrk = t->getZ() + zShiftTrk;
1073
1074 int32_t det = GetDetectorNumber(zTrk, alpha, layer);
1075
1076 // reject tracks between stacks
1077 if (det < 0) {
1078 return false;
1079 }
1080
1081 // reject tracks in PHOS hole and for non existent chamber 17_4_4
1082 if (!mGeo->ChamberInGeometry(det)) {
1083 return false;
1084 }
1085
1086 const GPUTRDpadPlane* pp = mGeo->GetPadPlane(det);
1087 float yMax = pp->GetColEnd();
1088 float zMax = pp->GetRow0();
1089 float zMin = pp->GetRowEnd();
1090
1091 float epsY = 5.f;
1092 float epsZ = 5.f;
1093
1094 // reject tracks closer than epsY cm to pad plane boundary
1095 if (yMax - CAMath::Abs(t->getY()) < epsY) {
1096 return false;
1097 }
1098 // reject tracks closer than epsZ cm to stack boundary
1099 if (!((zTrk > zMin + epsZ) && (zTrk < zMax - epsZ))) {
1100 return false;
1101 }
1102
1103 return true;
1104}
1105
1106#ifndef GPUCA_GPUCODE
1107namespace o2::gpu
1108{
1111} // namespace o2::gpu
1112#endif
int32_t i
bool const GPUTPCGMMerger::trackCluster const clcomparestruct * c2
float float float & zMax
float float & zMin
float & yMax
For performance analysis + error parametrization of the TRD tracker.
#define ENABLE_INFO
#define ENABLE_WARNING
Online TRD tracker based on extrapolated TPC tracks.
TRD Tracklet word for GPU tracker - 32bit tracklet info + half chamber ID + index.
uint32_t stack
Definition RawData.h:1
GPUTrackingInOutPointers & mIOPtrs
void * SetPointersTracks(void *base)
void * SetPointersBase(void *base)
void SetNCandidates(int32_t n)
void SetMaxData(const GPUTrackingInOutPointers &io)
void PrepareTracking(GPUChainTracking *chainTracking)
void * SetPointersTracklets(void *base)
GLdouble n
Definition glcorearb.h:1982
GLfloat GLfloat GLfloat alpha
Definition glcorearb.h:279
GLuint64EXT * result
Definition glcorearb.h:5662
typedef void(APIENTRYP PFNGLCULLFACEPROC)(GLenum mode)
GLuint GLfloat x0
Definition glcorearb.h:5034
GLenum GLuint GLint GLint layer
Definition glcorearb.h:1310
GPUdi() o2
Definition TrackTRD.h:38
GPUd() const expr uint32_t MultivariatePolynomialHelper< Dim
double * getX(double *xyDxy, int N)
double * getY(double *xyDxy, int N)
bool isOK(const SanityError &error)
GPUReconstruction * rec
GPUChainTracking * chainTracking
const GPUTRDSpacePoint * trdSpacePoints
std::vector< Tracklet64 > tracklets