Project
Loading...
Searching...
No Matches
GPUTPCGMO2Output.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#include "GPUTPCDef.h"
16#include "GPUTPCGMO2Output.h"
17#include "GPUCommonAlgorithm.h"
21#include "TPCFastTransform.h"
23
24#ifndef GPUCA_GPUCODE
27#include "GPUQAHelper.h"
28#endif
29
30using namespace o2::gpu;
31using namespace o2::tpc;
32using namespace o2::tpc::constants;
33
34GPUdi() static constexpr uint8_t getFlagsReject() { return GPUTPCGMMergedTrackHit::flagReject | GPUTPCGMMergedTrackHit::flagNotFit; }
35GPUdi() static uint32_t getFlagsRequired(const GPUSettingsRec& rec) { return rec.tpc.dropSecondaryLegsInOutput ? gputpcgmmergertypes::attachGoodLeg : gputpcgmmergertypes::attachZero; }
36
37template <>
38GPUdii() void GPUTPCGMO2Output::Thread<GPUTPCGMO2Output::prepare>(int32_t nBlocks, int32_t nThreads, int32_t iBlock, int32_t iThread, GPUsharedref() GPUSharedMemory& smem, processorType& GPUrestrict() merger)
39{
40 const GPUTPCGMMergedTrack* tracks = merger.OutputTracks();
41 const uint32_t nTracks = merger.NOutputTracks();
42 const GPUTPCGMMergedTrackHit* trackClusters = merger.Clusters();
43 const GPUdEdxInfo* tracksdEdx = merger.OutputTracksdEdx();
44
45 constexpr uint8_t flagsReject = getFlagsReject();
46 const uint32_t flagsRequired = getFlagsRequired(merger.Param().rec);
47 bool cutOnTrackdEdx = merger.Param().par.dodEdx && merger.Param().dodEdxDownscaled && merger.Param().rec.tpc.minTrackdEdxMax2Tot > 0.f;
48
49 GPUTPCGMMerger::tmpSort* GPUrestrict() trackSort = merger.TrackSortO2();
50 uint2* GPUrestrict() tmpData = merger.ClusRefTmp();
51 for (uint32_t i = get_global_id(0); i < nTracks; i += get_global_size(0)) {
52 if (!tracks[i].OK()) {
53 continue;
54 }
55 uint32_t nCl = 0;
56 for (uint32_t j = 0; j < tracks[i].NClusters(); j++) {
57 if ((trackClusters[tracks[i].FirstClusterRef() + j].state & flagsReject) || (merger.ClusterAttachment()[trackClusters[tracks[i].FirstClusterRef() + j].num] & flagsRequired) != flagsRequired) {
58 continue;
59 }
60 if (merger.Param().rec.tpc.dropSecondaryLegsInOutput && trackClusters[tracks[i].FirstClusterRef() + j].leg != trackClusters[tracks[i].FirstClusterRef() + tracks[i].NClusters() - 1].leg) {
61 continue;
62 }
63 nCl++;
64 }
65 if (nCl == 0) {
66 continue;
67 }
68 if (merger.Param().rec.tpc.dropSecondaryLegsInOutput && nCl + 2 < GPUCA_TRACKLET_SELECTOR_MIN_HITS_B5(tracks[i].GetParam().GetQPt() * merger.Param().qptB5Scaler)) { // Give 2 hits tolerance in the primary leg, compared to the full fit of the looper
69 continue;
70 }
71 if (merger.Param().rec.tpc.minNClustersFinalTrack != -1 && nCl < (uint32_t)merger.Param().rec.tpc.minNClustersFinalTrack) {
72 continue;
73 }
74 if (cutOnTrackdEdx && (tracksdEdx[i].dEdxMaxTPC < merger.Param().rec.tpc.minTrackdEdxMax || tracksdEdx[i].dEdxMaxTPC < tracksdEdx[i].dEdxTotTPC * merger.Param().rec.tpc.minTrackdEdxMax2Tot) && !(tracksdEdx[i].dEdxMaxTPC == 0 && CAMath::Abs(tracks[i].GetParam().GetDzDs()) > 0.03f)) {
75 continue;
76 }
77 uint32_t myId = CAMath::AtomicAdd(&merger.Memory()->nO2Tracks, 1u);
78 tmpData[i] = {nCl, CAMath::AtomicAdd(&merger.Memory()->nO2ClusRefs, nCl + (nCl + 1) / 2)};
79 trackSort[myId] = {i, (merger.Param().par.earlyTpcTransform || tracks[i].CSide()) ? tracks[i].GetParam().GetTZOffset() : -tracks[i].GetParam().GetTZOffset()};
80 }
81}
82
83template <>
84GPUdii() void GPUTPCGMO2Output::Thread<GPUTPCGMO2Output::sort>(int32_t nBlocks, int32_t nThreads, int32_t iBlock, int32_t iThread, GPUsharedref() GPUSharedMemory& smem, processorType& GPUrestrict() merger)
85{
86#ifndef GPUCA_SPECIALIZE_THRUST_SORTS
87 if (iThread || iBlock) {
88 return;
89 }
90 GPUTPCGMMerger::tmpSort* GPUrestrict() trackSort = merger.TrackSortO2();
91 auto comp = [](const auto& a, const auto& b) { return (a.y > b.y); };
92 GPUCommonAlgorithm::sortDeviceDynamic(trackSort, trackSort + merger.Memory()->nO2Tracks, comp);
93#endif
94}
95
96template <>
97GPUdii() void GPUTPCGMO2Output::Thread<GPUTPCGMO2Output::output>(int32_t nBlocks, int32_t nThreads, int32_t iBlock, int32_t iThread, GPUsharedref() GPUSharedMemory& smem, processorType& GPUrestrict() merger)
98{
99 constexpr float MinDelta = 0.1f;
100 const GPUTPCGMMergedTrack* tracks = merger.OutputTracks();
101 GPUdEdxInfo* tracksdEdx = merger.OutputTracksdEdx();
102 GPUdEdxInfo* tracksdEdxAlt = merger.OutputTracksdEdxAlt();
103 const int32_t nTracks = merger.NOutputTracksTPCO2();
104 const GPUTPCGMMergedTrackHit* trackClusters = merger.Clusters();
105 constexpr uint8_t flagsReject = getFlagsReject();
106 const uint32_t flagsRequired = getFlagsRequired(merger.Param().rec);
107 TrackTPC* outputTracks = merger.OutputTracksTPCO2();
108 uint32_t* clusRefs = merger.OutputClusRefsTPCO2();
109
110 GPUTPCGMMerger::tmpSort* GPUrestrict() trackSort = merger.TrackSortO2();
111 uint2* GPUrestrict() tmpData = merger.ClusRefTmp();
112 float const SNPThresh = 0.999990f;
113
114 for (int32_t iTmp = get_global_id(0); iTmp < nTracks; iTmp += get_global_size(0)) {
115 TrackTPC oTrack;
116 const int32_t i = trackSort[iTmp].x;
117 auto snpIn = tracks[i].GetParam().GetSinPhi();
118 if (snpIn > SNPThresh) {
119 snpIn = SNPThresh;
120 } else if (snpIn < -SNPThresh) {
121 snpIn = -SNPThresh;
122 }
123 oTrack.set(tracks[i].GetParam().GetX(), tracks[i].GetAlpha(),
124 {tracks[i].GetParam().GetY(), tracks[i].GetParam().GetZ(), snpIn, tracks[i].GetParam().GetDzDs(), tracks[i].GetParam().GetQPt()},
125 {tracks[i].GetParam().GetCov(0),
126 tracks[i].GetParam().GetCov(1), tracks[i].GetParam().GetCov(2),
127 tracks[i].GetParam().GetCov(3), tracks[i].GetParam().GetCov(4), tracks[i].GetParam().GetCov(5),
128 tracks[i].GetParam().GetCov(6), tracks[i].GetParam().GetCov(7), tracks[i].GetParam().GetCov(8), tracks[i].GetParam().GetCov(9),
129 tracks[i].GetParam().GetCov(10), tracks[i].GetParam().GetCov(11), tracks[i].GetParam().GetCov(12), tracks[i].GetParam().GetCov(13), tracks[i].GetParam().GetCov(14)});
130
131 oTrack.setChi2(tracks[i].GetParam().GetChi2());
132 auto& outerPar = tracks[i].OuterParam();
133 if (merger.Param().par.dodEdx && merger.Param().dodEdxDownscaled) {
134 oTrack.setdEdx(tracksdEdx[i]);
135 oTrack.setdEdxAlt(tracksdEdxAlt[i]);
136 }
137
138 auto snpOut = outerPar.P[2];
139 if (snpOut > SNPThresh) {
140 snpOut = SNPThresh;
141 } else if (snpOut < -SNPThresh) {
142 snpOut = -SNPThresh;
143 }
144 oTrack.setOuterParam(o2::track::TrackParCov(
145 outerPar.X, outerPar.alpha,
146 {outerPar.P[0], outerPar.P[1], snpOut, outerPar.P[3], outerPar.P[4]},
147 {outerPar.C[0], outerPar.C[1], outerPar.C[2], outerPar.C[3], outerPar.C[4], outerPar.C[5],
148 outerPar.C[6], outerPar.C[7], outerPar.C[8], outerPar.C[9], outerPar.C[10], outerPar.C[11],
149 outerPar.C[12], outerPar.C[13], outerPar.C[14]}));
150
151 if (merger.Param().par.dodEdx && merger.Param().dodEdxDownscaled && merger.Param().rec.tpc.enablePID) {
152 PIDResponse pidResponse{};
153 auto pid = pidResponse.getMostProbablePID(oTrack, merger.Param().rec.tpc.PID_EKrangeMin, merger.Param().rec.tpc.PID_EKrangeMax, merger.Param().rec.tpc.PID_EPrangeMin, merger.Param().rec.tpc.PID_EPrangeMax, merger.Param().rec.tpc.PID_EDrangeMin, merger.Param().rec.tpc.PID_EDrangeMax, merger.Param().rec.tpc.PID_ETrangeMin, merger.Param().rec.tpc.PID_ETrangeMax, merger.Param().rec.tpc.PID_useNsigma, merger.Param().rec.tpc.PID_sigma);
154 auto pidRemap = merger.Param().rec.tpc.PID_remap[pid];
155 if (pidRemap >= 0) {
156 pid = pidRemap;
157 }
158 oTrack.setPID(pid, true);
159 oTrack.getParamOut().setPID(pid, true);
160 }
161
162 uint32_t nOutCl = tmpData[i].x;
163 uint32_t clBuff = tmpData[i].y;
164 oTrack.setClusterRef(clBuff, nOutCl);
165 uint32_t* clIndArr = &clusRefs[clBuff];
166 uint8_t* sectorIndexArr = reinterpret_cast<uint8_t*>(clIndArr + nOutCl);
167 uint8_t* rowIndexArr = sectorIndexArr + nOutCl;
168
169 uint32_t nOutCl2 = 0;
170 float t1 = 0, t2 = 0;
171 int32_t sector1 = 0, sector2 = 0;
172 const o2::tpc::ClusterNativeAccess* GPUrestrict() clusters = merger.GetConstantMem()->ioPtrs.clustersNative;
173 for (uint32_t j = 0; j < tracks[i].NClusters(); j++) {
174 if ((trackClusters[tracks[i].FirstClusterRef() + j].state & flagsReject) || (merger.ClusterAttachment()[trackClusters[tracks[i].FirstClusterRef() + j].num] & flagsRequired) != flagsRequired) {
175 continue;
176 }
177 if (merger.Param().rec.tpc.dropSecondaryLegsInOutput && trackClusters[tracks[i].FirstClusterRef() + j].leg != trackClusters[tracks[i].FirstClusterRef() + tracks[i].NClusters() - 1].leg) {
178 continue;
179 }
180 int32_t clusterIdGlobal = trackClusters[tracks[i].FirstClusterRef() + j].num;
181 int32_t sector = trackClusters[tracks[i].FirstClusterRef() + j].sector;
182 int32_t globalRow = trackClusters[tracks[i].FirstClusterRef() + j].row;
183 int32_t clusterIdInRow = clusterIdGlobal - clusters->clusterOffset[sector][globalRow];
184 clIndArr[nOutCl2] = clusterIdInRow;
185 sectorIndexArr[nOutCl2] = sector;
186 rowIndexArr[nOutCl2] = globalRow;
187 if (nOutCl2 == 0) {
188 t1 = clusters->clustersLinear[clusterIdGlobal].getTime();
189 sector1 = sector;
190 }
191 if (++nOutCl2 == nOutCl) {
192 t2 = clusters->clustersLinear[clusterIdGlobal].getTime();
193 sector2 = sector;
194 }
195 }
196
197 bool cce = tracks[i].CCE() && ((sector1 < MAXSECTOR / 2) ^ (sector2 < MAXSECTOR / 2));
198 float time0 = 0.f, tFwd = 0.f, tBwd = 0.f;
199 if (merger.Param().par.continuousTracking) {
200 time0 = tracks[i].GetParam().GetTZOffset();
201 if (cce) {
202 bool lastSide = trackClusters[tracks[i].FirstClusterRef()].sector < MAXSECTOR / 2;
203 float delta = 0.f;
204 for (uint32_t iCl = 1; iCl < tracks[i].NClusters(); iCl++) {
205 auto& cacl1 = trackClusters[tracks[i].FirstClusterRef() + iCl];
206 if (lastSide ^ (cacl1.sector < MAXSECTOR / 2)) {
207 auto& cl1 = clusters->clustersLinear[cacl1.num];
208 auto& cl2 = clusters->clustersLinear[trackClusters[tracks[i].FirstClusterRef() + iCl - 1].num];
209 delta = CAMath::Abs(cl1.getTime() - cl2.getTime()) * 0.5f;
210 if (delta < MinDelta) {
211 delta = MinDelta;
212 }
213 break;
214 }
215 }
216 tFwd = tBwd = delta;
217 } else {
218 // estimate max/min time increments which still keep track in the physical limits of the TPC
219 const float tmin = CAMath::Min(t1, t2);
220 const float maxDriftTime = merger.GetConstantMem()->calibObjects.fastTransformHelper->getCorrMap()->getMaxDriftTime(t1 > t2 ? sector1 : sector2);
221 const float clusterT0 = merger.GetConstantMem()->calibObjects.fastTransformHelper->getCorrMap()->getT0();
222 const float tmax = CAMath::Min(tmin + maxDriftTime, CAMath::Max(t1, t2));
223 float delta = 0.f;
224 if (time0 + maxDriftTime < tmax) {
225 delta = tmax - time0 - maxDriftTime;
226 }
227 if (tmin - clusterT0 < time0 + delta) {
228 delta = tmin - clusterT0 - time0;
229 }
230 if (delta != 0.f) {
231 time0 += delta;
232 const float deltaZ = merger.GetConstantMem()->calibObjects.fastTransformHelper->getCorrMap()->convDeltaTimeToDeltaZinTimeFrame(sector2, delta);
233 oTrack.setZ(oTrack.getZ() + deltaZ);
234 }
235 tFwd = tmin - clusterT0 - time0;
236 tBwd = time0 - tmax + maxDriftTime;
237 }
238 }
239 if (tBwd < 0.f) {
240 tBwd = 0.f;
241 }
242 oTrack.setTime0(time0);
243 oTrack.setDeltaTBwd(tBwd);
244 oTrack.setDeltaTFwd(tFwd);
245 if (cce) {
246 oTrack.setHasCSideClusters();
247 oTrack.setHasASideClusters();
248 } else if (tracks[i].CSide()) {
249 oTrack.setHasCSideClusters();
250 } else {
251 oTrack.setHasASideClusters();
252 }
253 outputTracks[iTmp] = oTrack;
254 }
255}
256
257template <>
258GPUdii() void GPUTPCGMO2Output::Thread<GPUTPCGMO2Output::mc>(int32_t nBlocks, int32_t nThreads, int32_t iBlock, int32_t iThread, GPUsharedref() GPUSharedMemory& smem, processorType& GPUrestrict() merger)
259{
260#ifndef GPUCA_GPUCODE
261 const o2::tpc::ClusterNativeAccess* GPUrestrict() clusters = merger.GetConstantMem()->ioPtrs.clustersNative;
262 if (clusters == nullptr || clusters->clustersMCTruth == nullptr) {
263 return;
264 }
265 if (merger.OutputTracksTPCO2MC() == nullptr) {
266 return;
267 }
268
269 auto labelAssigner = GPUTPCTrkLbl(clusters->clustersMCTruth, 0.1f);
270 uint32_t* clusRefs = merger.OutputClusRefsTPCO2();
271 for (uint32_t i = get_global_id(0); i < merger.NOutputTracksTPCO2(); i += get_global_size(0)) {
272 labelAssigner.reset();
273 const auto& trk = merger.OutputTracksTPCO2()[i];
274 for (int32_t j = 0; j < trk.getNClusters(); j++) {
275 uint8_t sectorIndex, rowIndex;
276 uint32_t clusterIndex;
277 trk.getClusterReference(clusRefs, j, sectorIndex, rowIndex, clusterIndex);
278 uint32_t clusterIdGlobal = clusters->clusterOffset[sectorIndex][rowIndex] + clusterIndex;
279 labelAssigner.addLabel(clusterIdGlobal);
280 }
281 merger.OutputTracksTPCO2MC()[i] = labelAssigner.computeLabel();
282 }
283#endif
284}
benchmark::State & state
A const (ready only) version of MCTruthContainer.
Helper class to access correction maps.
int32_t i
#define GPUsharedref()
#define get_global_size(dim)
#define GPUrestrict()
#define get_global_id(dim)
#define GPUCA_TRACKLET_SELECTOR_MIN_HITS_B5(QPTB5)
GPUdii() void GPUTPCGMO2Output
void output(const std::map< std::string, ChannelStat > &channels)
Definition rawdump.cxx:197
uint32_t j
Definition RawData.h:0
uint16_t pid
Definition RawData.h:2
Definition of TPCFastTransform class.
PID response class.
Definition PIDResponse.h:39
GLboolean GLboolean GLboolean b
Definition glcorearb.h:1233
typedef void(APIENTRYP PFNGLCULLFACEPROC)(GLenum mode)
GLboolean GLboolean GLboolean GLboolean a
Definition glcorearb.h:1233
GLuint GLfloat GLfloat GLfloat GLfloat GLfloat GLfloat GLfloat GLfloat GLfloat t1
Definition glcorearb.h:5034
constexpr int MAXSECTOR
Definition Constants.h:28
Global TPC definitions and constants.
Definition SimTraits.h:167
GPUdi() T BetheBlochAleph(T bg
TrackParCovF TrackParCov
Definition Track.h:33
GPUReconstruction * rec
uint32_t x
uint32_t y
std::vector< Cluster > clusters