Project
Loading...
Searching...
No Matches
TrackInterpolation.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
17
20#include "ITStracking/IOUtils.h"
23#include "TOFBase/Geo.h"
25#include "DataFormatsTPC/Defs.h"
30#include "MathUtils/Tsallis.h"
31#include "TRDBase/PadPlane.h"
32#include "TMath.h"
34#include "Framework/Logger.h"
36#include <set>
37#include <algorithm>
38#include <random>
39
40using namespace o2::tpc;
43
44bool UnbinnedResid::gInitDone = false;
45
47{
48 if (!isITS()) {
50 }
51 // ITS alpha repends on the chip ID
54}
55
57{
58 if (isTPC()) {
59 return param::RowX[row];
60 }
62 if (isITS()) {
63 return o2::its::GeometryTGeo::Instance()->getSensorRefX(channel); // ITS X repends on the chip ID
64 }
65 if (isTRD()) {
66 auto geo = o2::trd::Geometry::instance();
67 ROOT::Math::Impl::Transform3D<double>::Point local{geo->cdrHght() - 0.5 - 0.279, 0., 0.}; // see TrackletTransformer::transformTracklet
68 return (geo->getMatrixT2L(channel) ^ local).X();
69 }
70 if (isTOF()) {
71 int det[5];
73 float pos[3] = {0.f, 0.f, 0.f};
75 float posl[3] = {pos[0], pos[1], pos[2]};
77 return pos[2]; // coordinates in sector frame: note that the rotation above puts z in pos[1], the radial coordinate in pos[2], and the tangent coordinate in pos[0] (this is to match the TOF residual system, where we don't use the radial component), so we swap their positions.
78 }
79 LOGP(fatal, "Did not recognize detector type: row:{}, sec:{}, channel:{}", row, sec, channel);
80 return 0.;
81}
82
84{
85 if (!gInitDone) {
86 LOGP(warn, "geometry initialization was not done, doing this for the current timestamp");
87 init();
88 if (!gInitDone) {
89 LOGP(fatal, "geometry initialization failed");
90 }
91 }
92}
93
94void UnbinnedResid::init(long timestamp)
95{
96 if (gInitDone) {
97 LOGP(warn, "Initialization was already done");
98 return;
99 }
100 if (!gGeoManager) {
101 o2::ccdb::BasicCCDBManager::instance().getSpecific<TGeoManager>("GLO/Config/GeometryAligned", timestamp);
102 }
103 auto geoTRD = o2::trd::Geometry::instance();
104 geoTRD->createPadPlaneArray();
105 geoTRD->createClusterMatrixArray();
106 gInitDone = true;
107}
108
110{
111 // perform initialization
112 LOG(info) << "Start initializing TrackInterpolation";
113 if (mInitDone) {
114 LOG(error) << "Initialization already performed.";
115 return;
116 }
117
118 const auto& elParam = ParameterElectronics::Instance();
119 mTPCTimeBinMUS = elParam.ZbinWidth;
120
121 mFastTransform = std::move(TPCFastTransformHelperO2::instance()->create(0));
122
123 mBz = o2::base::Propagator::Instance()->getNominalBz();
124 mRecoParam.setBfield(mBz);
125 mGeoTRD = o2::trd::Geometry::instance();
127
128 mSourcesConfigured = src;
129 mSourcesConfiguredMap = srcMap;
130 mSingleSourcesConfigured = (mSourcesConfigured == mSourcesConfiguredMap);
131 mTrackTypes.insert({GTrackID::ITSTPC, 0});
132 mTrackTypes.insert({GTrackID::ITSTPCTRD, 1});
133 mTrackTypes.insert({GTrackID::ITSTPCTOF, 2});
134 mTrackTypes.insert({GTrackID::ITSTPCTRDTOF, 3});
135
137 geom->fillMatrixCache(o2::math_utils::bit2Mask(o2::math_utils::TransformType::T2L, o2::math_utils::TransformType::L2G));
138
139 mInitDone = true;
140 LOGP(info, "Done initializing TrackInterpolation. Configured track input: {}. Track input specifically for map: {}",
141 GTrackID::getSourcesNames(mSourcesConfigured), mSingleSourcesConfigured ? "identical" : GTrackID::getSourcesNames(mSourcesConfiguredMap));
142}
143
145{
146 LOGP(debug, "Check if input track {} is accepted", gid.asString());
147 bool hasOuterPoint = gidTable[GTrackID::TRD].isIndexSet() || gidTable[GTrackID::TOF].isIndexSet();
148 if (!hasOuterPoint && !mProcessITSTPConly) {
149 return false; // don't do ITS-only extrapolation through TPC
150 }
151 const auto itsTrk = mRecoCont->getITSTrack(gidTable[GTrackID::ITS]);
152 const auto tpcTrk = mRecoCont->getTPCTrack(gidTable[GTrackID::TPC]);
153
154 if (gidTable[GTrackID::TRD].isIndexSet()) {
155 // TRD specific cuts
156 const auto& trdTrk = mRecoCont->getITSTPCTRDTrack<o2::trd::TrackTRD>(gidTable[GTrackID::ITSTPCTRD]);
157 if (trdTrk.getNtracklets() < mParams->minTRDNTrklts) {
158 return false;
159 }
160 }
161 // reduced chi2 cut is the same for all track types
162 if (itsTrk.getChi2() / itsTrk.getNumberOfClusters() > mParams->maxITSChi2 || tpcTrk.getChi2() / tpcTrk.getNClusterReferences() > mParams->maxTPCChi2) {
163 return false;
164 }
165 if (!hasOuterPoint) {
166 // ITS-TPC track (does not have outer points in TRD or TOF)
167 if (itsTrk.getNumberOfClusters() < mParams->minITSNClsNoOuterPoint || tpcTrk.getNClusterReferences() < mParams->minTPCNClsNoOuterPoint) {
168 return false;
169 }
170 if (itsTrk.getPt() < mParams->minPtNoOuterPoint) {
171 return false;
172 }
173 } else {
174 if (itsTrk.getNumberOfClusters() < mParams->minITSNCls || tpcTrk.getNClusterReferences() < mParams->minTPCNCls) {
175 return false;
176 }
177 }
178
179 auto trc = mRecoCont->getTrackParam(gid);
181 auto prop = o2::base::Propagator::Instance();
182 if (!prop->propagateToDCA(pv, trc, prop->getNominalBz(), 2., o2::base::PropagatorF::MatCorrType::USEMatCorrLUT, &dca)) {
183 return false;
184 }
185 if (dca.getR2() > mParams->maxDCA * mParams->maxDCA) {
186 return false;
187 }
188 return true;
189}
190
192{
193 LOGP(debug, "Trying to find valid source for {} in {}", GTrackID::getSourceName(src), GTrackID::getSourcesNames(mask));
196 return GTrackID::ITSTPCTRD;
197 } else if (mask[GTrackID::ITSTPC]) {
198 return GTrackID::ITSTPC;
199 } else {
200 return GTrackID::NSources;
201 }
202 } else if (src == GTrackID::ITSTPCTRD || src == GTrackID::ITSTPCTOF) {
203 if (mask[GTrackID::ITSTPC]) {
204 return GTrackID::ITSTPC;
205 } else {
206 return GTrackID::NSources;
207 }
208 } else {
209 return GTrackID::NSources;
210 }
211}
212
214{
215 mRecoCont = &inp;
216 uint32_t nTrackSeeds = 0;
217 uint32_t countSeedCandidates[4] = {0};
218 auto pvvec = mRecoCont->getPrimaryVertices();
219 auto trackIndex = mRecoCont->getPrimaryVertexMatchedTracks(); // Global ID's for associated tracks
220 auto vtxRefs = mRecoCont->getPrimaryVertexMatchedTrackRefs(); // references from vertex to these track IDs
221 int nv = vtxRefs.size() - 1;
222 GTrackID::mask_t allowedSources = GTrackID::getSourcesMask("ITS-TPC,ITS-TPC-TRD,ITS-TPC-TOF,ITS-TPC-TRD-TOF") & mSourcesConfigured;
223 constexpr std::array<int, 3> SrcFast = {int(GTrackID::ITSTPCTRD), int(GTrackID::ITSTPCTOF), int(GTrackID::ITSTPCTRDTOF)};
224
225 for (int iv = 0; iv < nv; iv++) {
226 LOGP(debug, "processing PV {} of {}", iv, nv);
227
228 const auto& vtref = vtxRefs[iv];
229 auto pv = pvvec[iv];
230 if (mParams->minTOFTRDPVContributors > 0) { // we want only PVs constrained by fast detectors
231 int nfound = 0;
232 bool usePV = false;
233 for (uint32_t is = 0; is < SrcFast.size() && !usePV; is++) {
234 int src = SrcFast[is], idMin = vtref.getFirstEntryOfSource(src), idMax = idMin + vtref.getEntriesOfSource(src);
235 for (int i = idMin; i < idMax; i++) {
236 if (trackIndex[i].isPVContributor() && (++nfound == mParams->minTOFTRDPVContributors)) {
237 usePV = true;
238 break;
239 }
240 }
241 }
242 if (!usePV) {
243 continue;
244 }
245 }
246
247 for (int is = GTrackID::NSources; is >= 0; is--) {
248 if (!allowedSources[is]) {
249 continue;
250 }
251 LOGP(debug, "Checking source {}", is);
252 int idMin = vtref.getFirstEntryOfSource(is), idMax = idMin + vtref.getEntriesOfSource(is);
253 for (int i = idMin; i < idMax; i++) {
254 auto vid = trackIndex[i];
255 auto vidOrig = vid; // in case only ITS-TPC tracks are configured vid might be overwritten. We need to remember it for the PID
256 if (mParams->ignoreNonPVContrib && !vid.isPVContributor()) {
257 continue;
258 }
259 if (vid.isAmbiguous()) {
260 continue;
261 }
262 auto gidTable = mRecoCont->getSingleDetectorRefs(vid);
263 if (!mSourcesConfigured[is]) {
264 auto src = findValidSource(mSourcesConfigured, static_cast<GTrackID::Source>(vid.getSource()));
266 LOGP(debug, "prepareInputTrackSample: Found valid source {}", GTrackID::getSourceName(src));
267 vid = gidTable[src];
268 gidTable = mRecoCont->getSingleDetectorRefs(vid);
269 } else {
270 break; // no valid source for this vertex track source exists
271 }
272 }
273 ++countSeedCandidates[mTrackTypes[vid.getSource()]];
274 LOGP(debug, "Checking vid {}", vid.asString());
275 if (!isInputTrackAccepted(vid, gidTable, pv)) {
276 continue;
277 }
278 mSeeds.push_back(mRecoCont->getITSTrack(gidTable[GTrackID::ITS]).getParamOut());
279 mSeeds.back().setPID(mRecoCont->getTrackParam(vidOrig).getPID(), true);
280 mGIDs.push_back(vid);
281 mGIDtables.push_back(gidTable);
282 mTrackTimes.push_back(pv.getTimeStamp().getTimeStamp());
283 mTrackIndices[mTrackTypes[vid.getSource()]].push_back(nTrackSeeds++);
284 }
285 }
286 }
287
288 LOGP(info, "Created {} seeds. {} out of {} ITS-TPC-TRD-TOF, {} out of {} ITS-TPC-TRD, {} out of {} ITS-TPC-TOF, {} out of {} ITS-TPC",
289 nTrackSeeds,
290 mTrackIndices[mTrackTypes[GTrackID::ITSTPCTRDTOF]].size(), countSeedCandidates[mTrackTypes[GTrackID::ITSTPCTRDTOF]],
291 mTrackIndices[mTrackTypes[GTrackID::ITSTPCTRD]].size(), countSeedCandidates[mTrackTypes[GTrackID::ITSTPCTRD]],
292 mTrackIndices[mTrackTypes[GTrackID::ITSTPCTOF]].size(), countSeedCandidates[mTrackTypes[GTrackID::ITSTPCTOF]],
293 mTrackIndices[mTrackTypes[GTrackID::ITSTPC]].size(), countSeedCandidates[mTrackTypes[GTrackID::ITSTPC]]);
294}
295
297{
298 std::random_device rd;
299 std::mt19937 g(rd());
300 std::uniform_real_distribution<> distr(0., 1.);
301 float weight = 0;
303}
304
306{
307 // main processing function
308 if (!mInitDone) {
309 LOG(error) << "Initialization not yet done. Aborting...";
310 return;
311 }
312 // set the input containers
313 mTPCTracksClusIdx = mRecoCont->getTPCTracksClusterRefs();
314 mTPCClusterIdxStruct = &mRecoCont->getTPCClusters();
315 {
316 if (!mITSDict) {
317 LOG(error) << "No ITS dictionary available";
318 return;
319 }
320 mITSTrackClusIdx = mRecoCont->getITSTracksClusterRefs();
321 const auto clusITS = mRecoCont->getITSClusters();
322 const auto patterns = mRecoCont->getITSClustersPatterns();
323 auto pattIt = patterns.begin();
324 mITSClustersArray.clear();
325 mITSClustersArray.reserve(clusITS.size());
326 LOGP(info, "We have {} ITS clusters and the number of patterns is {}", clusITS.size(), patterns.size());
327 o2::its::ioutils::convertCompactClusters(clusITS, pattIt, mITSClustersArray, mITSDict);
328 }
329
330 // In case we have more input tracks available than are required per TF
331 // we want to sample them. But we still prefer global ITS-TPC-TRD-TOF tracks
332 // over ITS-TPC-TRD tracks and so on. So we have to shuffle the indices
333 // in blocks.
334 std::random_device rd;
335 std::mt19937 g(rd());
336 std::vector<uint32_t> trackIndices; // here we keep the GIDs for all track types in a single vector to use in loop
337 std::shuffle(mTrackIndices[mTrackTypes[GTrackID::ITSTPCTRDTOF]].begin(), mTrackIndices[mTrackTypes[GTrackID::ITSTPCTRDTOF]].end(), g);
338 std::shuffle(mTrackIndices[mTrackTypes[GTrackID::ITSTPCTRD]].begin(), mTrackIndices[mTrackTypes[GTrackID::ITSTPCTRD]].end(), g);
339 std::shuffle(mTrackIndices[mTrackTypes[GTrackID::ITSTPCTOF]].begin(), mTrackIndices[mTrackTypes[GTrackID::ITSTPCTOF]].end(), g);
340 std::shuffle(mTrackIndices[mTrackTypes[GTrackID::ITSTPC]].begin(), mTrackIndices[mTrackTypes[GTrackID::ITSTPC]].end(), g);
341 trackIndices.insert(trackIndices.end(), mTrackIndices[mTrackTypes[GTrackID::ITSTPCTRDTOF]].begin(), mTrackIndices[mTrackTypes[GTrackID::ITSTPCTRDTOF]].end());
342 trackIndices.insert(trackIndices.end(), mTrackIndices[mTrackTypes[GTrackID::ITSTPCTRD]].begin(), mTrackIndices[mTrackTypes[GTrackID::ITSTPCTRD]].end());
343 trackIndices.insert(trackIndices.end(), mTrackIndices[mTrackTypes[GTrackID::ITSTPCTOF]].begin(), mTrackIndices[mTrackTypes[GTrackID::ITSTPCTOF]].end());
344 trackIndices.insert(trackIndices.end(), mTrackIndices[mTrackTypes[GTrackID::ITSTPC]].begin(), mTrackIndices[mTrackTypes[GTrackID::ITSTPC]].end());
345
346 int nSeeds = mSeeds.size(), lastChecked = 0;
347 mParentID.clear();
348 mParentID.resize(nSeeds, -1);
349
350 int maxOutputTracks = (mMaxTracksPerTF >= 0) ? mMaxTracksPerTF + mAddTracksForMapPerTF : nSeeds;
351 mTrackData.reserve(maxOutputTracks);
352 mClRes.reserve(maxOutputTracks * param::NPadRows);
353 bool maxTracksReached = false;
354 for (int iSeed = 0; iSeed < nSeeds; ++iSeed) {
355 if (mMaxTracksPerTF >= 0 && mTrackDataCompact.size() >= mMaxTracksPerTF + mAddTracksForMapPerTF) {
356 LOG(info) << "Maximum number of tracks per TF reached. Skipping the remaining " << nSeeds - iSeed << " tracks.";
357 break;
358 }
359 int seedIndex = trackIndices[iSeed];
360 if (mParams->enableTrackDownsampling && !isTrackSelected(mSeeds[seedIndex])) {
361 continue;
362 }
363
364 auto addPart = [this, seedIndex](GTrackID::Source src) {
365 this->mGIDs.push_back(this->mGIDtables[seedIndex][src]);
366 this->mGIDtables.push_back(this->mRecoCont->getSingleDetectorRefs(this->mGIDs.back()));
367 this->mTrackTimes.push_back(this->mTrackTimes[seedIndex]);
368 this->mSeeds.push_back(this->mSeeds[seedIndex]);
369 this->mParentID.push_back(seedIndex); // store parent seed id
370 };
371
372 GTrackID::mask_t partsAdded;
373 if (!mSingleSourcesConfigured && !mSourcesConfiguredMap[mGIDs[seedIndex].getSource()]) {
374 auto src = findValidSource(mSourcesConfiguredMap, static_cast<GTrackID::Source>(mGIDs[seedIndex].getSource()));
376 LOGP(debug, "process {}: Found valid source {} for {} | nseeds:{} mSeeds:{} used: {}", iSeed, GTrackID::getSourceName(src), GTrackID::getSourceName(mGIDs[seedIndex].getSource()), nSeeds, mSeeds.size(), mTrackDataCompact.size());
377 addPart(src);
378 }
379 }
380 if (mMaxTracksPerTF >= 0 && mTrackDataCompact.size() >= mMaxTracksPerTF) {
381 if (!maxTracksReached) {
382 LOGP(info, "We already have reached mMaxTracksPerTF={}, but we continue to create seeds until mAddTracksForMapPerTF={} is also reached, iSeed: {} of {} inital seeds", mMaxTracksPerTF, mAddTracksForMapPerTF, iSeed, nSeeds);
383 }
384 maxTracksReached = true;
385 continue;
386 }
387 if (mGIDs[seedIndex].includesDet(DetID::TRD) || mGIDs[seedIndex].includesDet(DetID::TOF)) {
388 interpolateTrack(seedIndex);
389 LOGP(debug, "interpolateTrack {} {}, accepted: {}", iSeed, GTrackID::getSourceName(mGIDs[seedIndex].getSource()), mTrackDataCompact.size());
390 if (mProcessSeeds) {
391 if (mGIDs[seedIndex].includesDet(DetID::TRD) && mGIDs[seedIndex].includesDet(DetID::TOF) && !partsAdded[GTrackID::ITSTPCTRD]) {
392 addPart(GTrackID::ITSTPCTRD);
393 }
394 if (!partsAdded[GTrackID::ITSTPC]) {
395 addPart(GTrackID::ITSTPC);
396 }
397 }
398 } else {
399 extrapolateTrack(seedIndex);
400 LOGP(debug, "extrapolateTrack {} {}, accepted: {}", iSeed, GTrackID::getSourceName(mGIDs[seedIndex].getSource()), mTrackDataCompact.size());
401 }
402 lastChecked = iSeed;
403 }
404 std::vector<int> remSeeds;
405 if (mSeeds.size() > ++lastChecked) {
406 remSeeds.resize(mSeeds.size() - lastChecked);
407 std::iota(remSeeds.begin(), remSeeds.end(), lastChecked);
408 std::shuffle(remSeeds.begin(), remSeeds.end(), g);
409 LOGP(info, "Up to {} tracks out of {} additional seeds will be processed in random order, of which {} are stripped versions, accepted seeds: {}", mAddTracksForMapPerTF, remSeeds.size(), mSeeds.size() - nSeeds, mTrackDataCompact.size());
410 }
411 int extraChecked = 0;
412 for (int iSeed : remSeeds) {
413 if (mAddTracksForMapPerTF > 0 && mTrackDataCompact.size() >= mMaxTracksPerTF + mAddTracksForMapPerTF) {
414 LOGP(info, "Maximum number {} of additional tracks per TF reached. Skipping the remaining {} tracks", mAddTracksForMapPerTF, remSeeds.size() - extraChecked);
415 break;
416 }
417 extraChecked++;
418 if (mGIDs[iSeed].includesDet(DetID::TRD) || mGIDs[iSeed].includesDet(DetID::TOF)) {
419 interpolateTrack(iSeed);
420 LOGP(debug, "extra check {} of {}, seed {} interpolateTrack {}, used: {}", extraChecked, remSeeds.size(), iSeed, GTrackID::getSourceName(mGIDs[iSeed].getSource()), mTrackDataCompact.size());
421 } else {
422 LOGP(debug, "extra check {} of {}, seed {} extrapolateTrack {}, used: {}", extraChecked, remSeeds.size(), iSeed, GTrackID::getSourceName(mGIDs[iSeed].getSource()), mTrackDataCompact.size());
423 extrapolateTrack(iSeed);
424 }
425 }
426 LOG(info) << "Could process " << mTrackData.size() << " tracks successfully. " << mRejectedResiduals << " residuals were rejected. " << mClRes.size() << " residuals were accepted.";
427 mRejectedResiduals = 0;
428}
429
431{
432 LOGP(debug, "Starting track interpolation for GID {}", mGIDs[iSeed].asString());
433 TrackData trackData;
434 std::unique_ptr<TrackDataExtended> trackDataExtended;
435 std::vector<TPCClusterResiduals> clusterResiduals;
436 auto propagator = o2::base::Propagator::Instance();
437 const auto& gidTable = mGIDtables[iSeed];
438 const auto& trkTPC = mRecoCont->getTPCTrack(gidTable[GTrackID::TPC]);
439 const auto& trkITS = mRecoCont->getITSTrack(gidTable[GTrackID::ITS]);
440 if (mDumpTrackPoints) {
441 trackDataExtended = std::make_unique<TrackDataExtended>();
442 (*trackDataExtended).gid = mGIDs[iSeed];
443 (*trackDataExtended).clIdx.setFirstEntry(mClRes.size());
444 (*trackDataExtended).trkITS = trkITS;
445 (*trackDataExtended).trkTPC = trkTPC;
446 auto nCl = trkITS.getNumberOfClusters();
447 auto clEntry = trkITS.getFirstClusterEntry();
448 for (int iCl = nCl - 1; iCl >= 0; iCl--) { // clusters are stored from outer to inner layers
449 const auto& clsITS = mITSClustersArray[mITSTrackClusIdx[clEntry + iCl]];
450 (*trackDataExtended).clsITS.push_back(clsITS);
451 }
452 }
453 trackData.gid = mGIDs[iSeed];
454 trackData.par = mSeeds[iSeed];
455 auto& trkWork = mSeeds[iSeed];
456 o2::track::TrackPar trkInner{trkWork};
457 // reset the cache array (sufficient to set cluster available to zero)
458 for (auto& elem : mCache) {
459 elem.clAvailable = 0;
460 }
461 trackData.clIdx.setFirstEntry(mClRes.size()); // reference the first cluster residual belonging to this track
462 float clusterTimeBinOffset = mTrackTimes[iSeed] / mTPCTimeBinMUS;
463
464 // store the TPC cluster positions in the cache
465 for (int iCl = trkTPC.getNClusterReferences(); iCl--;) {
466 uint8_t sector, row;
467 uint32_t clusterIndexInRow;
468 const auto& clTPC = trkTPC.getCluster(mTPCTracksClusIdx, iCl, *mTPCClusterIdxStruct, sector, row);
469 float clTPCX;
470 std::array<float, 2> clTPCYZ;
471 mFastTransform->TransformIdeal(sector, row, clTPC.getPad(), clTPC.getTime(), clTPCX, clTPCYZ[0], clTPCYZ[1], clusterTimeBinOffset);
472 mCache[row].clSec = sector;
473 mCache[row].clAvailable = 1;
474 mCache[row].clY = clTPCYZ[0];
475 mCache[row].clZ = clTPCYZ[1];
476 mCache[row].clAngle = o2::math_utils::sector2Angle(sector);
477 }
478
479 // extrapolate seed through TPC and store track position at each pad row
480 for (int iRow = 0; iRow < param::NPadRows; ++iRow) {
481 if (!mCache[iRow].clAvailable) {
482 continue;
483 }
484 if (!trkWork.rotate(mCache[iRow].clAngle)) {
485 LOG(debug) << "Failed to rotate track during first extrapolation";
486 return;
487 }
488 if (!propagator->PropagateToXBxByBz(trkWork, param::RowX[iRow], mParams->maxSnp, mParams->maxStep, mMatCorr)) {
489 LOG(debug) << "Failed on first extrapolation";
490 return;
491 }
492 mCache[iRow].y[ExtOut] = trkWork.getY();
493 mCache[iRow].z[ExtOut] = trkWork.getZ();
494 mCache[iRow].sy2[ExtOut] = trkWork.getSigmaY2();
495 mCache[iRow].szy[ExtOut] = trkWork.getSigmaZY();
496 mCache[iRow].sz2[ExtOut] = trkWork.getSigmaZ2();
497 mCache[iRow].snp[ExtOut] = trkWork.getSnp();
498 // printf("Track alpha at row %i: %.2f, Y(%.2f), Z(%.2f)\n", iRow, trkWork.getAlpha(), trkWork.getY(), trkWork.getZ());
499 }
500
501 // start from outermost cluster with outer refit and back propagation
502 if (gidTable[GTrackID::TOF].isIndexSet()) {
503 LOG(debug) << "TOF point available";
504 const auto& clTOF = mRecoCont->getTOFClusters()[gidTable[GTrackID::TOF]];
505 if (mDumpTrackPoints) {
506 (*trackDataExtended).clsTOF = clTOF;
507 (*trackDataExtended).matchTOF = mRecoCont->getTOFMatch(mGIDs[iSeed]);
508 }
509 const int clTOFSec = clTOF.getCount();
510 const float clTOFAlpha = o2::math_utils::sector2Angle(clTOFSec);
511 if (!trkWork.rotate(clTOFAlpha)) {
512 LOG(debug) << "Failed to rotate into TOF cluster sector frame";
513 return;
514 }
515 float clTOFxyz[3] = {clTOF.getX(), clTOF.getY(), clTOF.getZ()};
516 if (!clTOF.isInNominalSector()) {
517 o2::tof::Geo::alignedToNominalSector(clTOFxyz, clTOFSec); // go from the aligned to nominal sector frame
518 }
519 std::array<float, 2> clTOFYZ{clTOFxyz[1], clTOFxyz[2]};
520 std::array<float, 3> clTOFCov{mParams->sigYZ2TOF, 0.f, mParams->sigYZ2TOF}; // assume no correlation between y and z and equal cluster error sigma^2 = (3cm)^2 / 12
521 if (!propagator->PropagateToXBxByBz(trkWork, clTOFxyz[0], mParams->maxSnp, mParams->maxStep, mMatCorr)) {
522 LOG(debug) << "Failed final propagation to TOF radius";
523 return;
524 }
525 // TODO: check if reset of covariance matrix is needed here (or, in case TOF point is not available at outermost TRD layer)
526 if (!trkWork.update(clTOFYZ, clTOFCov)) {
527 LOG(debug) << "Failed to update extrapolated ITS track with TOF cluster";
528 // LOGF(info, "trkWork.y=%f, cl.y=%f, trkWork.z=%f, cl.z=%f", trkWork.getY(), clTOFYZ[0], trkWork.getZ(), clTOFYZ[1]);
529 return;
530 }
531 }
532 if (gidTable[GTrackID::TRD].isIndexSet()) {
533 LOG(debug) << "TRD available";
534 const auto& trkTRD = mRecoCont->getITSTPCTRDTrack<o2::trd::TrackTRD>(gidTable[GTrackID::ITSTPCTRD]);
535 if (mDumpTrackPoints) {
536 (*trackDataExtended).trkTRD = trkTRD;
537 }
538 for (int iLayer = o2::trd::constants::NLAYER - 1; iLayer >= 0; --iLayer) {
539 std::array<float, 2> trkltTRDYZ{};
540 std::array<float, 3> trkltTRDCov{};
541 int res = processTRDLayer(trkTRD, iLayer, trkWork, &trkltTRDYZ, &trkltTRDCov);
542 if (res == -1) { // no TRD tracklet in this layer
543 continue;
544 }
545 if (res < -1) { // failed to reach this layer
546 return;
547 }
548 if (!trkWork.update(trkltTRDYZ, trkltTRDCov)) {
549 LOG(debug) << "Failed to update track at TRD layer " << iLayer;
550 return;
551 }
552 }
553 }
554
555 if (mDumpTrackPoints) {
556 (*trackDataExtended).trkOuter = trkWork;
557 }
558 auto trkOuter = trkWork; // outer param
559
560 // go back through the TPC and store updated track positions
561 bool outerParamStored = false;
562 for (int iRow = param::NPadRows; iRow--;) {
563 if (!mCache[iRow].clAvailable) {
564 continue;
565 }
566 if (mProcessSeeds && !outerParamStored) {
567 // for debug purposes we store the track parameters
568 // of the refitted ITS-(TRD)-(TOF) track at the
569 // outermose TPC cluster if we are processing all seeds
570 // i.e. if we in any case also process the ITS-TPC only
571 // part of the same track
572 trackData.par = trkWork;
573 outerParamStored = true;
574 }
575 if (!trkWork.rotate(mCache[iRow].clAngle)) {
576 LOG(debug) << "Failed to rotate track during back propagation";
577 return;
578 }
579 if (!propagator->PropagateToXBxByBz(trkWork, param::RowX[iRow], mParams->maxSnp, mParams->maxStep, mMatCorr)) {
580 LOG(debug) << "Failed on back propagation";
581 // printf("trkX(%.2f), clX(%.2f), clY(%.2f), clZ(%.2f), alphaTOF(%.2f)\n", trkWork.getX(), param::RowX[iRow], clTOFYZ[0], clTOFYZ[1], clTOFAlpha);
582 return;
583 }
584 mCache[iRow].y[ExtIn] = trkWork.getY();
585 mCache[iRow].z[ExtIn] = trkWork.getZ();
586 mCache[iRow].sy2[ExtIn] = trkWork.getSigmaY2();
587 mCache[iRow].szy[ExtIn] = trkWork.getSigmaZY();
588 mCache[iRow].sz2[ExtIn] = trkWork.getSigmaZ2();
589 mCache[iRow].snp[ExtIn] = trkWork.getSnp();
590 }
591
592 // calculate weighted mean at each pad row (assume for now y and z are uncorrelated) and store residuals to TPC clusters
593 unsigned short deltaRow = 0;
594 for (int iRow = 0; iRow < param::NPadRows; ++iRow) {
595 if (!mCache[iRow].clAvailable) {
596 ++deltaRow;
597 continue;
598 }
599 float wTotY = 1.f / mCache[iRow].sy2[ExtOut] + 1.f / mCache[iRow].sy2[ExtIn];
600 float wTotZ = 1.f / mCache[iRow].sz2[ExtOut] + 1.f / mCache[iRow].sz2[ExtIn];
601 mCache[iRow].y[Int] = (mCache[iRow].y[ExtOut] / mCache[iRow].sy2[ExtOut] + mCache[iRow].y[ExtIn] / mCache[iRow].sy2[ExtIn]) / wTotY;
602 mCache[iRow].z[Int] = (mCache[iRow].z[ExtOut] / mCache[iRow].sz2[ExtOut] + mCache[iRow].z[ExtIn] / mCache[iRow].sz2[ExtIn]) / wTotZ;
603
604 // simple average w/o weighting for angle
605 mCache[iRow].snp[Int] = (mCache[iRow].snp[ExtOut] + mCache[iRow].snp[ExtIn]) / 2.f;
606
607 const auto dY = mCache[iRow].clY - mCache[iRow].y[Int];
608 const auto dZ = mCache[iRow].clZ - mCache[iRow].z[Int];
609 const auto y = mCache[iRow].y[Int];
610 const auto z = mCache[iRow].z[Int];
611 const auto snp = mCache[iRow].snp[Int];
612 const auto sec = mCache[iRow].clSec;
613 clusterResiduals.emplace_back(dY, dZ, y, z, snp, sec, deltaRow);
614
615 deltaRow = 1;
616 }
617 trackData.chi2TRD = gidTable[GTrackID::TRD].isIndexSet() ? mRecoCont->getITSTPCTRDTrack<o2::trd::TrackTRD>(gidTable[GTrackID::ITSTPCTRD]).getChi2() : 0;
618 trackData.chi2TPC = trkTPC.getChi2();
619 trackData.chi2ITS = trkITS.getChi2();
620 trackData.nClsTPC = trkTPC.getNClusterReferences();
621 trackData.nClsITS = trkITS.getNumberOfClusters();
622 trackData.nTrkltsTRD = gidTable[GTrackID::TRD].isIndexSet() ? mRecoCont->getITSTPCTRDTrack<o2::trd::TrackTRD>(gidTable[GTrackID::ITSTPCTRD]).getNtracklets() : 0;
623 if (gidTable[GTrackID::TOF].isIndexSet()) {
624 const auto& tofMatch = mRecoCont->getTOFMatch(mGIDs[iSeed]);
625 trackData.deltaTOF = tofMatch.getSignal() - tofMatch.getFT0Best() - tofMatch.getLTIntegralOut().getTOF(trkTPC.getPID().getID());
626 trackData.clAvailTOF = uint16_t(tofMatch.getFT0BestRes());
627 } else {
628 trackData.clAvailTOF = 0;
629 }
630 trackData.dEdxTPC = trkTPC.getdEdx().dEdxTotTPC;
631
632 TrackParams params; // for refitted track parameters and flagging rejected clusters
633 if (mParams->skipOutlierFiltering || validateTrack(trackData, params, clusterResiduals)) {
634 // track is good
635 int nClValidated = 0;
636 int iRow = 0;
637 for (unsigned int iCl = 0; iCl < clusterResiduals.size(); ++iCl) {
638 iRow += clusterResiduals[iCl].dRow;
639 if (params.flagRej[iCl]) {
640 // skip masked cluster residual
641 continue;
642 }
643 ++nClValidated;
644 const float tgPhi = clusterResiduals[iCl].snp / std::sqrt((1.f - clusterResiduals[iCl].snp) * (1.f + clusterResiduals[iCl].snp));
645 const auto dy = clusterResiduals[iCl].dy;
646 const auto dz = clusterResiduals[iCl].dz;
647 const auto y = clusterResiduals[iCl].y;
648 const auto z = clusterResiduals[iCl].z;
649 const auto sec = clusterResiduals[iCl].sec;
650 if ((std::abs(dy) < param::MaxResid) && (std::abs(dz) < param::MaxResid) && (std::abs(y) < param::MaxY) && (std::abs(z) < param::MaxZ) && (std::abs(tgPhi) < param::MaxTgSlp)) {
651 mClRes.emplace_back(dy, dz, tgPhi, y, z, iRow, sec);
652 } else {
653 ++mRejectedResiduals;
654 }
655 }
656 trackData.clIdx.setEntries(nClValidated);
657
658 bool stopPropagation = !mExtDetResid;
659 if (!stopPropagation) {
660 // do we have TRD residuals to add?
661 trkWork = trkOuter;
662 if (gidTable[GTrackID::TRD].isIndexSet()) {
663 const auto& trkTRD = mRecoCont->getITSTPCTRDTrack<o2::trd::TrackTRD>(gidTable[GTrackID::ITSTPCTRD]);
664 for (int iLayer = 0; iLayer < o2::trd::constants::NLAYER; iLayer++) {
665 std::array<float, 2> trkltTRDYZ{};
666 int res = processTRDLayer(trkTRD, iLayer, trkWork, &trkltTRDYZ, nullptr, &trackData);
667 if (res == -1) { // no traklet on this layer
668 continue;
669 }
670 if (res < -1) { // failed to reach this layer
671 stopPropagation = true;
672 break;
673 }
674
675 float tgPhi = trkWork.getSnp() / std::sqrt((1.f - trkWork.getSnp()) * (1.f + trkWork.getSnp()));
676 auto dy = trkltTRDYZ[0] - trkWork.getY();
677 auto dz = trkltTRDYZ[1] - trkWork.getZ();
678 if ((std::abs(dy) < param::MaxResid) && (std::abs(dz) < param::MaxResid) && (std::abs(trkWork.getY()) < param::MaxY) && (std::abs(trkWork.getZ()) < param::MaxZ) && (std::abs(tgPhi) < param::MaxTgSlp)) {
679 mClRes.emplace_back(dy, dz, tgPhi, trkWork.getY(), trkWork.getZ(), 160 + iLayer, o2::math_utils::angle2Sector(trkWork.getAlpha()), (short)res);
680 trackData.nExtDetResid++;
681 }
682 }
683 }
684
685 // do we have TOF residual to add?
686 while (gidTable[GTrackID::TOF].isIndexSet() && !stopPropagation) {
687 const auto& clTOF = mRecoCont->getTOFClusters()[gidTable[GTrackID::TOF]];
688 float clTOFxyz[3] = {clTOF.getX(), clTOF.getY(), clTOF.getZ()};
689 if (!clTOF.isInNominalSector()) {
690 o2::tof::Geo::alignedToNominalSector(clTOFxyz, clTOF.getCount()); // go from the aligned to nominal sector frame
691 }
692 const float clTOFAlpha = o2::math_utils::sector2Angle(clTOF.getCount());
693 if (trkWork.getAlpha() != clTOFAlpha && !trkWork.rotate(clTOFAlpha)) {
694 LOG(debug) << "Failed to rotate into TOF cluster sector frame";
695 stopPropagation = true;
696 break;
697 }
698 if (!propagator->PropagateToXBxByBz(trkWork, clTOFxyz[0], mParams->maxSnp, mParams->maxStep, mMatCorr)) {
699 LOG(debug) << "Failed final propagation to TOF radius";
700 break;
701 }
702
703 float tgPhi = trkWork.getSnp() / std::sqrt((1.f - trkWork.getSnp()) * (1.f + trkWork.getSnp()));
704 auto dy = clTOFxyz[1] - trkWork.getY();
705 auto dz = clTOFxyz[2] - trkWork.getZ();
706 if ((std::abs(dy) < param::MaxResid) && (std::abs(dz) < param::MaxResid) && (std::abs(trkWork.getY()) < param::MaxY) && (std::abs(trkWork.getZ()) < param::MaxZ) && (std::abs(tgPhi) < param::MaxTgSlp)) {
707 mClRes.emplace_back(dy, dz, tgPhi, trkWork.getY(), trkWork.getZ(), 170, clTOF.getCount(), clTOF.getPadInSector());
708 trackData.nExtDetResid++;
709 }
710 break;
711 }
712
713 // add ITS residuals
714 while (!stopPropagation) {
715 auto& trkWorkITS = trkInner; // this is ITS outer param
716 auto nCl = trkITS.getNumberOfClusters();
717 auto clEntry = trkITS.getFirstClusterEntry();
719 for (int iCl = 0; iCl < nCl; iCl++) { // clusters are stored from outer to inner layers
720 const auto& cls = mITSClustersArray[mITSTrackClusIdx[clEntry + iCl]];
721 int chip = cls.getSensorID();
722 float chipX, chipAlpha;
723 geom->getSensorXAlphaRefPlane(cls.getSensorID(), chipX, chipAlpha);
724 if (!trkWorkITS.rotate(chipAlpha) || !propagator->PropagateToXBxByBz(trkWorkITS, chipX, mParams->maxSnp, mParams->maxStep, mMatCorr)) {
725 LOGP(debug, "Failed final propagation to ITS X={} alpha={}", chipX, chipAlpha);
726 stopPropagation = true;
727 break;
728 }
729 float tgPhi = trkWorkITS.getSnp() / std::sqrt((1.f - trkWorkITS.getSnp()) * (1.f + trkWorkITS.getSnp()));
730 auto dy = cls.getY() - trkWorkITS.getY();
731 auto dz = cls.getZ() - trkWorkITS.getZ();
732 if ((std::abs(dy) < param::MaxResid) && (std::abs(dz) < param::MaxResid) && (std::abs(trkWorkITS.getY()) < param::MaxY) && (std::abs(trkWorkITS.getZ()) < param::MaxZ) && (std::abs(tgPhi) < param::MaxTgSlp)) {
733 mClRes.emplace_back(dy, dz, tgPhi, trkWorkITS.getY(), trkWorkITS.getZ(), 180 + geom->getLayer(cls.getSensorID()), -1, cls.getSensorID());
734 trackData.nExtDetResid++;
735 }
736 }
737 break;
738 }
739 }
740
741 mGIDsSuccess.push_back(mGIDs[iSeed]);
742 mTrackDataCompact.emplace_back(trackData.clIdx.getFirstEntry(), nClValidated, mGIDs[iSeed].getSource(), trackData.nExtDetResid);
743 mTrackData.push_back(std::move(trackData));
744 if (mDumpTrackPoints) {
745 (*trackDataExtended).clIdx.setEntries(nClValidated);
746 (*trackDataExtended).nExtDetResid = trackData.nExtDetResid;
747 mTrackDataExtended.push_back(std::move(*trackDataExtended));
748 }
749 }
750 if (mParams->writeUnfiltered) {
751 TrackData trkDataTmp = trackData;
752 trkDataTmp.clIdx.setFirstEntry(mClResUnfiltered.size());
753 trkDataTmp.clIdx.setEntries(clusterResiduals.size());
754 mTrackDataUnfiltered.push_back(std::move(trkDataTmp));
755 mClResUnfiltered.insert(mClResUnfiltered.end(), clusterResiduals.begin(), clusterResiduals.end());
756 }
757}
758
760 std::array<float, 2>* trkltTRDYZ, std::array<float, 3>* trkltTRDCov, TrackData* trkData)
761{
762 // return chamber ID (0:539) in case of successful processing, -1 if there is no TRD tracklet at given layer, -2 if processing failed
763 int trkltIdx = trkTRD.getTrackletIndex(iLayer);
764 if (trkltIdx < 0) {
765 return -1; // no TRD tracklet in this layer
766 }
767 const auto& trdSP = mRecoCont->getTRDCalibratedTracklets()[trkltIdx];
768 const auto& trdTrklt = mRecoCont->getTRDTracklets()[trkltIdx];
769 auto trkltDet = trdTrklt.getDetector();
770 auto trkltSec = trkltDet / (o2::trd::constants::NLAYER * o2::trd::constants::NSTACK);
771 if (trkltSec != o2::math_utils::angle2Sector(trkWork.getAlpha())) {
772 if (!trkWork.rotate(o2::math_utils::sector2Angle(trkltSec))) {
773 LOG(debug) << "Track could not be rotated in TRD tracklet coordinate system in layer " << iLayer;
774 return -2;
775 }
776 }
777 if (!o2::base::Propagator::Instance()->PropagateToXBxByBz(trkWork, trdSP.getX(), mParams->maxSnp, mParams->maxStep, mMatCorr)) {
778 LOG(debug) << "Failed propagation to TRD layer " << iLayer;
779 return -2;
780 }
781 if (trkltTRDYZ) {
782 const auto* pad = mGeoTRD->getPadPlane(trkltDet);
783 float tilt = tan(TMath::DegToRad() * pad->getTiltingAngle()); // tilt is signed! and returned in degrees
784 float tiltCorrUp = tilt * (trdSP.getZ() - trkWork.getZ());
785 float zPosCorrUp = trdSP.getZ() + mRecoParam.getZCorrCoeffNRC() * trkWork.getTgl(); // maybe Z can be corrected on avarage already by the tracklet transformer?
786 float padLength = pad->getRowSize(trdTrklt.getPadRow());
787 if (!((trkWork.getSigmaZ2() < (padLength * padLength / 12.f)) && (std::fabs(trdSP.getZ() - trkWork.getZ()) < padLength))) {
788 tiltCorrUp = 0.f;
789 }
790 (*trkltTRDYZ)[0] = trdSP.getY() - tiltCorrUp;
791 (*trkltTRDYZ)[1] = zPosCorrUp;
792 if (trkltTRDCov) {
793 mRecoParam.recalcTrkltCov(tilt, trkWork.getSnp(), pad->getRowSize(trdTrklt.getPadRow()), *trkltTRDCov);
794 }
795 }
796 if (trkData) {
797 auto slope = trdSP.getDy();
798 if (std::abs(slope) < param::MaxTRDSlope) {
799 trkData->TRDTrkltSlope[iLayer] = slope * 0x7fff / param::MaxTRDSlope;
800 }
801 }
802 return trkltDet;
803}
804
806{
807 // extrapolate ITS-only track through TPC and store residuals to TPC clusters in the output vectors
808 LOGP(debug, "Starting track extrapolation for GID {}", mGIDs[iSeed].asString());
809 const auto& gidTable = mGIDtables[iSeed];
810 TrackData trackData;
811 std::unique_ptr<TrackDataExtended> trackDataExtended;
812 std::vector<TPCClusterResiduals> clusterResiduals;
813 trackData.clIdx.setFirstEntry(mClRes.size());
814 const auto& trkITS = mRecoCont->getITSTrack(gidTable[GTrackID::ITS]);
815 const auto& trkTPC = mRecoCont->getTPCTrack(gidTable[GTrackID::TPC]);
816 if (mDumpTrackPoints) {
817 trackDataExtended = std::make_unique<TrackDataExtended>();
818 (*trackDataExtended).gid = mGIDs[iSeed];
819 (*trackDataExtended).clIdx.setFirstEntry(mClRes.size());
820 (*trackDataExtended).trkITS = trkITS;
821 (*trackDataExtended).trkTPC = trkTPC;
822 auto nCl = trkITS.getNumberOfClusters();
823 auto clEntry = trkITS.getFirstClusterEntry();
824 for (int iCl = nCl - 1; iCl >= 0; iCl--) { // clusters are stored from outer to inner layers
825 const auto& clsITS = mITSClustersArray[mITSTrackClusIdx[clEntry + iCl]];
826 (*trackDataExtended).clsITS.push_back(clsITS);
827 }
828 }
829 trackData.gid = mGIDs[iSeed];
830 trackData.par = mSeeds[iSeed];
831
832 auto trkWork = mSeeds[iSeed];
833 float clusterTimeBinOffset = mTrackTimes[iSeed] / mTPCTimeBinMUS;
834 auto propagator = o2::base::Propagator::Instance();
835 unsigned short rowPrev = 0; // used to calculate dRow of two consecutive cluster residuals
836 unsigned short nMeasurements = 0;
837 uint8_t clRowPrev = constants::MAXGLOBALPADROW; // used to identify and skip split clusters on the same pad row
838 for (int iCl = trkTPC.getNClusterReferences(); iCl--;) {
839 uint8_t sector, row;
840 uint32_t clusterIndexInRow;
841 const auto& cl = trkTPC.getCluster(mTPCTracksClusIdx, iCl, *mTPCClusterIdxStruct, sector, row);
842 if (clRowPrev == row) {
843 // if there are split clusters we only take the first one on the pad row
844 continue;
845 } else if (clRowPrev < constants::MAXGLOBALPADROW && clRowPrev > row) {
846 // we seem to be looping, abort this track
847 LOGP(debug, "TPC track with pT={} GeV and {} clusters has cluster {} on row {} while the previous cluster was on row {}",
848 mSeeds[iSeed].getPt(), trkTPC.getNClusterReferences(), iCl, row, clRowPrev);
849 return;
850 } else {
851 // this is the first cluster we see on this pad row
852 clRowPrev = row;
853 }
854 float x = 0, y = 0, z = 0;
855 mFastTransform->TransformIdeal(sector, row, cl.getPad(), cl.getTime(), x, y, z, clusterTimeBinOffset);
856 if (!trkWork.rotate(o2::math_utils::sector2Angle(sector))) {
857 return;
858 }
859 if (!propagator->PropagateToXBxByBz(trkWork, x, mParams->maxSnp, mParams->maxStep, mMatCorr)) {
860 return;
861 }
862
863 const auto dY = y - trkWork.getY();
864 const auto dZ = z - trkWork.getZ();
865 const auto ty = trkWork.getY();
866 const auto tz = trkWork.getZ();
867 const auto snp = trkWork.getSnp();
868 const auto sec = sector;
869
870 clusterResiduals.emplace_back(dY, dZ, ty, tz, snp, sec, row - rowPrev);
871
872 rowPrev = row;
873 ++nMeasurements;
874 }
875
876 TrackParams params; // for refitted track parameters and flagging rejected clusters
877 if (clusterResiduals.size() > constants::MAXGLOBALPADROW) {
878 LOGP(warn, "Extrapolated ITS-TPC track and found more reesiduals than possible ({})", clusterResiduals.size());
879 return;
880 }
881
882 trackData.chi2TPC = trkTPC.getChi2();
883 trackData.chi2ITS = trkITS.getChi2();
884 trackData.nClsTPC = trkTPC.getNClusterReferences();
885 trackData.nClsITS = trkITS.getNumberOfClusters();
886 trackData.clIdx.setEntries(nMeasurements);
887 trackData.dEdxTPC = trkTPC.getdEdx().dEdxTotTPC;
888 if (mDumpTrackPoints) {
889 (*trackDataExtended).trkOuter = trkWork;
890 }
891
892 if (mParams->skipOutlierFiltering || validateTrack(trackData, params, clusterResiduals)) {
893 // track is good, store TPC part
894
895 int nClValidated = 0, iRow = 0;
896 unsigned int iCl = 0;
897 for (iCl = 0; iCl < clusterResiduals.size(); ++iCl) {
898 iRow += clusterResiduals[iCl].dRow;
899 if (iRow < param::NPadRows && params.flagRej[iCl]) { // skip masked cluster residual
900 continue;
901 }
902 ++nClValidated;
903 const float tgPhi = clusterResiduals[iCl].snp / std::sqrt((1.f - clusterResiduals[iCl].snp) * (1.f + clusterResiduals[iCl].snp));
904 const auto dy = clusterResiduals[iCl].dy;
905 const auto dz = clusterResiduals[iCl].dz;
906 const auto y = clusterResiduals[iCl].y;
907 const auto z = clusterResiduals[iCl].z;
908 if ((std::abs(dy) < param::MaxResid) && (std::abs(dz) < param::MaxResid) && (std::abs(y) < param::MaxY) && (std::abs(z) < param::MaxZ) && (std::abs(tgPhi) < param::MaxTgSlp)) {
909 mClRes.emplace_back(dy, dz, tgPhi, y, z, iRow, clusterResiduals[iCl].sec);
910 } else {
911 ++mRejectedResiduals;
912 }
913 }
914 trackData.clIdx.setEntries(nClValidated);
915
916 bool stopPropagation = !mExtDetResid;
917 if (!stopPropagation) {
918 // do we have TRD residuals to add?
919 int iSeedFull = mParentID[iSeed] == -1 ? iSeed : mParentID[iSeed];
920 auto gidFull = mGIDs[iSeedFull];
921 const auto& gidTableFull = mGIDtables[iSeedFull];
922 if (gidTableFull[GTrackID::TRD].isIndexSet()) {
923 const auto& trkTRD = mRecoCont->getITSTPCTRDTrack<o2::trd::TrackTRD>(gidTableFull[GTrackID::ITSTPCTRD]);
924 trackData.nTrkltsTRD = trkTRD.getNtracklets();
925 for (int iLayer = 0; iLayer < o2::trd::constants::NLAYER; iLayer++) {
926 std::array<float, 2> trkltTRDYZ{};
927 int res = processTRDLayer(trkTRD, iLayer, trkWork, &trkltTRDYZ, nullptr, &trackData);
928 if (res == -1) { // no traklet on this layer
929 continue;
930 }
931 if (res < -1) { // failed to reach this layer
932 stopPropagation = true;
933 break;
934 }
935
936 float tgPhi = trkWork.getSnp() / std::sqrt((1.f - trkWork.getSnp()) * (1.f + trkWork.getSnp()));
937 auto dy = trkltTRDYZ[0] - trkWork.getY();
938 auto dz = trkltTRDYZ[1] - trkWork.getZ();
939 const auto sec = clusterResiduals[iCl].sec;
940 if ((std::abs(dy) < param::MaxResid) && (std::abs(dz) < param::MaxResid) && (std::abs(trkWork.getY()) < param::MaxY) && (std::abs(trkWork.getZ()) < param::MaxZ) && (std::abs(tgPhi) < param::MaxTgSlp)) {
941 mClRes.emplace_back(dy, dz, tgPhi, trkWork.getY(), trkWork.getZ(), 160 + iLayer, o2::math_utils::angle2Sector(trkWork.getAlpha()), (short)res);
942 trackData.nExtDetResid++;
943 }
944 }
945 }
946
947 // do we have TOF residual to add?
948 trackData.clAvailTOF = 0;
949 while (gidTableFull[GTrackID::TOF].isIndexSet() && !stopPropagation) {
950 const auto& tofMatch = mRecoCont->getTOFMatch(gidFull);
951 trackData.deltaTOF = tofMatch.getSignal() - tofMatch.getFT0Best() - tofMatch.getLTIntegralOut().getTOF(trkTPC.getPID().getID());
952 trackData.clAvailTOF = uint16_t(tofMatch.getFT0BestRes());
953 const auto& clTOF = mRecoCont->getTOFClusters()[gidTableFull[GTrackID::TOF]];
954 const float clTOFAlpha = o2::math_utils::sector2Angle(clTOF.getCount());
955 float clTOFxyz[3] = {clTOF.getX(), clTOF.getY(), clTOF.getZ()};
956 if (!clTOF.isInNominalSector()) {
957 o2::tof::Geo::alignedToNominalSector(clTOFxyz, clTOF.getCount()); // go from the aligned to nominal sector frame
958 }
959 if (trkWork.getAlpha() != clTOFAlpha && !trkWork.rotate(clTOFAlpha)) {
960 LOG(debug) << "Failed to rotate into TOF cluster sector frame";
961 stopPropagation = true;
962 break;
963 }
964 if (!propagator->PropagateToXBxByBz(trkWork, clTOFxyz[0], mParams->maxSnp, mParams->maxStep, mMatCorr)) {
965 LOG(debug) << "Failed final propagation to TOF radius";
966 break;
967 }
968
969 float tgPhi = trkWork.getSnp() / std::sqrt((1.f - trkWork.getSnp()) * (1.f + trkWork.getSnp()));
970 auto dy = clTOFxyz[1] - trkWork.getY();
971 auto dz = clTOFxyz[2] - trkWork.getZ();
972 if ((std::abs(dy) < param::MaxResid) && (std::abs(dz) < param::MaxResid) && (std::abs(trkWork.getY()) < param::MaxY) && (std::abs(trkWork.getZ()) < param::MaxZ) && (std::abs(tgPhi) < param::MaxTgSlp)) {
973 mClRes.emplace_back(dy, dz, tgPhi, trkWork.getY(), trkWork.getZ(), 170, clTOF.getCount(), clTOF.getPadInSector());
974 trackData.nExtDetResid++;
975 }
976 break;
977 }
978
979 // add ITS residuals
980 while (!stopPropagation) {
981 o2::track::TrackPar trkWorkITS{trackData.par}; // this is ITS outer param
982 auto nCl = trkITS.getNumberOfClusters();
983 auto clEntry = trkITS.getFirstClusterEntry();
985 for (int iCl = 0; iCl < nCl; iCl++) { // clusters are stored from outer to inner layers
986 const auto& cls = mITSClustersArray[mITSTrackClusIdx[clEntry + iCl]];
987 int chip = cls.getSensorID();
988 float chipX, chipAlpha;
989 geom->getSensorXAlphaRefPlane(cls.getSensorID(), chipX, chipAlpha);
990 if (!trkWorkITS.rotate(chipAlpha) || !propagator->PropagateToXBxByBz(trkWorkITS, chipX, mParams->maxSnp, mParams->maxStep, mMatCorr)) {
991 LOGP(debug, "Failed final propagation to ITS X={} alpha={}", chipX, chipAlpha);
992 stopPropagation = true;
993 break;
994 }
995 float tgPhi = trkWorkITS.getSnp() / std::sqrt((1.f - trkWorkITS.getSnp()) * (1.f + trkWorkITS.getSnp()));
996 auto dy = cls.getY() - trkWorkITS.getY();
997 auto dz = cls.getZ() - trkWorkITS.getZ();
998 if ((std::abs(dy) < param::MaxResid) && (std::abs(dz) < param::MaxResid) && (std::abs(trkWorkITS.getY()) < param::MaxY) && (std::abs(trkWorkITS.getZ()) < param::MaxZ) && (std::abs(tgPhi) < param::MaxTgSlp)) {
999 mClRes.emplace_back(dy, dz, tgPhi, trkWorkITS.getY(), trkWorkITS.getZ(), 180 + geom->getLayer(cls.getSensorID()), -1, cls.getSensorID());
1000 trackData.nExtDetResid++;
1001 }
1002 }
1003 break;
1004 }
1005 }
1006 mTrackData.push_back(std::move(trackData));
1007 mGIDsSuccess.push_back(mGIDs[iSeed]);
1008 mTrackDataCompact.emplace_back(trackData.clIdx.getFirstEntry(), nClValidated, mGIDs[iSeed].getSource(), trackData.nExtDetResid);
1009 if (mDumpTrackPoints) {
1010 (*trackDataExtended).clIdx.setEntries(nClValidated);
1011 (*trackDataExtended).nExtDetResid = trackData.nExtDetResid;
1012 mTrackDataExtended.push_back(std::move(*trackDataExtended));
1013 }
1014 }
1015 if (mParams->writeUnfiltered) {
1016 TrackData trkDataTmp = trackData;
1017 trkDataTmp.clIdx.setFirstEntry(mClResUnfiltered.size());
1018 trkDataTmp.clIdx.setEntries(clusterResiduals.size());
1019 mTrackDataUnfiltered.push_back(std::move(trkDataTmp));
1020 mClResUnfiltered.insert(mClResUnfiltered.end(), clusterResiduals.begin(), clusterResiduals.end());
1021 }
1022}
1023
1024bool TrackInterpolation::validateTrack(const TrackData& trk, TrackParams& params, const std::vector<TPCClusterResiduals>& clsRes) const
1025{
1026 if (clsRes.size() < mParams->minNCl) {
1027 // no enough clusters for this track to be considered
1028 LOG(debug) << "Skipping track with too few clusters: " << clsRes.size();
1029 return false;
1030 }
1031
1032 bool resHelix = compareToHelix(trk, params, clsRes);
1033 if (!resHelix) {
1034 LOG(debug) << "Skipping track too far from helix approximation";
1035 return false;
1036 }
1037 if (fabsf(mBz) > 0.01 && fabsf(params.qpt) > mParams->maxQ2Pt) {
1038 LOG(debug) << "Skipping track with too high q/pT: " << params.qpt;
1039 return false;
1040 }
1041 if (!outlierFiltering(trk, params, clsRes)) {
1042 return false;
1043 }
1044 return true;
1045}
1046
1047bool TrackInterpolation::compareToHelix(const TrackData& trk, TrackParams& params, const std::vector<TPCClusterResiduals>& clsRes) const
1048{
1049 std::array<float, param::NPadRows> residHelixY;
1050 std::array<float, param::NPadRows> residHelixZ;
1051
1052 std::array<float, param::NPadRows> xLab;
1053 std::array<float, param::NPadRows> yLab;
1054 std::array<float, param::NPadRows> sPath;
1055
1056 float curvature = fabsf(trk.par.getQ2Pt() * mBz * o2::constants::physics::LightSpeedCm2S * 1e-14f);
1057 int secFirst = clsRes[0].sec;
1058 float phiSect = (secFirst + .5f) * o2::constants::math::SectorSpanRad;
1059 float snPhi = sin(phiSect);
1060 float csPhi = cos(phiSect);
1061 sPath[0] = 0.f;
1062
1063 int iRow = 0;
1064 int nCl = clsRes.size();
1065 for (unsigned int iP = 0; iP < nCl; ++iP) {
1066 iRow += clsRes[iP].dRow;
1067 float yTrk = clsRes[iP].y;
1068 // LOGF(info, "iRow(%i), yTrk(%f)", iRow, yTrk);
1069 xLab[iP] = param::RowX[iRow];
1070 if (clsRes[iP].sec != secFirst) {
1071 float phiSectCurrent = (clsRes[iP].sec + .5f) * o2::constants::math::SectorSpanRad;
1072 float cs = cos(phiSectCurrent - phiSect);
1073 float sn = sin(phiSectCurrent - phiSect);
1074 xLab[iP] = param::RowX[iRow] * cs - yTrk * sn;
1075 yLab[iP] = yTrk * cs + param::RowX[iRow] * sn;
1076 } else {
1077 xLab[iP] = param::RowX[iRow];
1078 yLab[iP] = yTrk;
1079 }
1080 // this is needed only later, but we retrieve it already now to save another loop
1081 params.zTrk[iP] = clsRes[iP].z;
1082 params.xTrk[iP] = param::RowX[iRow];
1083 params.dy[iP] = clsRes[iP].dy;
1084 params.dz[iP] = clsRes[iP].dz;
1085 // done retrieving values for later
1086 if (iP > 0) {
1087 float dx = xLab[iP] - xLab[iP - 1];
1088 float dy = yLab[iP] - yLab[iP - 1];
1089 float ds2 = dx * dx + dy * dy;
1090 float ds = sqrt(ds2); // circular path (linear approximation)
1091 // if the curvature of the track or the (approximated) chord length is too large the more exact formula is used:
1092 // chord length = 2r * asin(ds/(2r))
1093 // using the first two terms of the tailer expansion for asin(x) ~ x + x^3 / 6
1094 if (ds * curvature > 0.05) {
1095 ds *= (1.f + ds2 * curvature * curvature / 24.f);
1096 }
1097 sPath[iP] = sPath[iP - 1] + ds;
1098 }
1099 }
1100 if (fabsf(mBz) < 0.01) {
1101 // for B=0 we don't need to try a circular fit...
1102 return true;
1103 }
1104 float xcSec = 0.f;
1105 float ycSec = 0.f;
1106 float r = 0.f;
1107 TrackResiduals::fitCircle(nCl, xLab, yLab, xcSec, ycSec, r, residHelixY);
1108 // LOGF(info, "Done with circle fit. nCl(%i), xcSec(%f), ycSec(%f), r(%f).", nCl, xcSec, ycSec, r);
1109 /*
1110 for (int i=0; i<nCl; ++i) {
1111 LOGF(info, "i(%i), xLab(%f), yLab(%f).", i, xLab[i], yLab[i]);
1112 }
1113 */
1114 // determine curvature
1115 float phiI = TMath::ATan2(yLab[0], xLab[0]);
1116 float phiF = TMath::ATan2(yLab[nCl - 1], xLab[nCl - 1]);
1117 if (phiI < 0) {
1119 }
1120 if (phiF < 0) {
1122 }
1123 float dPhi = phiF - phiI;
1124 float curvSign = -1.f;
1125 if (dPhi > 0) {
1126 if (dPhi < o2::constants::math::PI) {
1127 curvSign = 1.f;
1128 }
1129 } else if (dPhi < -o2::constants::math::PI) {
1130 curvSign = 1.f;
1131 }
1132 params.qpt = std::copysign(1.f / (r * mBz * o2::constants::physics::LightSpeedCm2S * 1e-14f), curvSign);
1133
1134 // calculate circle coordinates in the lab frame
1135 float xc = xcSec * csPhi - ycSec * snPhi;
1136 float yc = xcSec * snPhi + ycSec * csPhi;
1137
1138 std::array<float, 2> pol1Z;
1139 TrackResiduals::fitPoly1(nCl, sPath, params.zTrk, pol1Z);
1140
1141 params.tgl = pol1Z[0];
1142
1143 // max deviations in both directions from helix fit in y and z
1144 float hMinY = 1e9f;
1145 float hMaxY = -1e9f;
1146 float hMinZ = 1e9f;
1147 float hMaxZ = -1e9f;
1148 // extract residuals in Z and fill track slopes in sector frame
1149 int secCurr = secFirst;
1150 iRow = 0;
1151 for (unsigned int iCl = 0; iCl < nCl; ++iCl) {
1152 iRow += clsRes[iCl].dRow;
1153 float resZ = params.zTrk[iCl] - (pol1Z[1] + sPath[iCl] * pol1Z[0]);
1154 residHelixZ[iCl] = resZ;
1155 if (resZ < hMinZ) {
1156 hMinZ = resZ;
1157 }
1158 if (resZ > hMaxZ) {
1159 hMaxZ = resZ;
1160 }
1161 if (residHelixY[iCl] < hMinY) {
1162 hMinY = residHelixY[iCl];
1163 }
1164 if (residHelixY[iCl] > hMaxY) {
1165 hMaxY = residHelixY[iCl];
1166 }
1167 int sec = clsRes[iCl].sec;
1168 if (sec != secCurr) {
1169 secCurr = sec;
1170 phiSect = (.5f + sec) * o2::constants::math::SectorSpanRad;
1171 snPhi = sin(phiSect);
1172 csPhi = cos(phiSect);
1173 xcSec = xc * csPhi + yc * snPhi; // recalculate circle center in the sector frame
1174 }
1175 float cstalp = (param::RowX[iRow] - xcSec) / r;
1176 if (fabsf(cstalp) > 1.f - sFloatEps) {
1177 // track cannot reach this pad row
1178 cstalp = std::copysign(1.f - sFloatEps, cstalp);
1179 }
1180 params.tglArr[iCl] = cstalp / sqrt((1 - cstalp) * (1 + cstalp)); // 1 / tan(acos(cstalp)) = cstalp / sqrt(1 - cstalp^2)
1181
1182 // In B+ the slope of q- should increase with x. Just look on q * B
1183 if (params.qpt * mBz > 0) {
1184 params.tglArr[iCl] *= -1.f;
1185 }
1186 }
1187 // LOGF(info, "CompareToHelix: hMaxY(%f), hMinY(%f), hMaxZ(%f), hMinZ(%f). Max deviation allowed: y(%.2f), z(%.2f)", hMaxY, hMinY, hMaxZ, hMinZ, mParams->maxDevHelixY, mParams->maxDevHelixZ);
1188 // LOGF(info, "New pt/Q (%f), old pt/Q (%f)", 1./params.qpt, 1./trk.qPt);
1189 return fabsf(hMaxY - hMinY) < mParams->maxDevHelixY && fabsf(hMaxZ - hMinZ) < mParams->maxDevHelixZ;
1190}
1191
1192bool TrackInterpolation::outlierFiltering(const TrackData& trk, TrackParams& params, const std::vector<TPCClusterResiduals>& clsRes) const
1193{
1194 if (clsRes.size() < mParams->nMALong) {
1195 LOG(debug) << "Skipping track with too few clusters for long moving average: " << clsRes.size();
1196 return false;
1197 }
1198 float rmsLong = checkResiduals(trk, params, clsRes);
1199 if (static_cast<float>(params.flagRej.count()) / clsRes.size() > mParams->maxRejFrac) {
1200 LOGP(debug, "Skipping track with too many clusters rejected: {} out of {}", params.flagRej.count(), clsRes.size());
1201 return false;
1202 }
1203 if (rmsLong > mParams->maxRMSLong) {
1204 LOG(debug) << "Skipping track with too large RMS: " << rmsLong;
1205 return false;
1206 }
1207 return true;
1208}
1209
1210float TrackInterpolation::checkResiduals(const TrackData& trk, TrackParams& params, const std::vector<TPCClusterResiduals>& clsRes) const
1211{
1212 float rmsLong = 0.f;
1213
1214 int nCl = clsRes.size();
1215 int iClFirst = 0;
1216 int iClLast = nCl - 1;
1217 int secStart = clsRes[0].sec;
1218
1219 // arrays with differences / abs(differences) of points to their neighbourhood, initialized to zero
1220 std::array<float, param::NPadRows> yDiffLL{};
1221 std::array<float, param::NPadRows> zDiffLL{};
1222 std::array<float, param::NPadRows> absDevY{};
1223 std::array<float, param::NPadRows> absDevZ{};
1224
1225 for (unsigned int iCl = 0; iCl < nCl; ++iCl) {
1226 if (iCl < iClLast && clsRes[iCl].sec == secStart) {
1227 continue;
1228 }
1229 // sector changed or last cluster reached
1230 // now run estimators for all points in the same sector
1231 int nClSec = iCl - iClFirst;
1232 if (iCl == iClLast) {
1233 ++nClSec;
1234 }
1235 diffToLocLine(nClSec, iClFirst, params.xTrk, params.dy, yDiffLL);
1236 diffToLocLine(nClSec, iClFirst, params.xTrk, params.dz, zDiffLL);
1237 iClFirst = iCl;
1238 secStart = clsRes[iCl].sec;
1239 }
1240 // store abs deviations
1241 int nAccY = 0;
1242 int nAccZ = 0;
1243 for (int iCl = nCl; iCl--;) {
1244 if (fabsf(yDiffLL[iCl]) > param::sEps) {
1245 absDevY[nAccY++] = fabsf(yDiffLL[iCl]);
1246 }
1247 if (fabsf(zDiffLL[iCl]) > param::sEps) {
1248 absDevZ[nAccZ++] = fabsf(zDiffLL[iCl]);
1249 }
1250 }
1251 if (nAccY < mParams->minNumberOfAcceptedResiduals || nAccZ < mParams->minNumberOfAcceptedResiduals) {
1252 // mask all clusters
1253 LOGP(debug, "Accepted {} clusters for dY {} clusters for dZ, but required at least {} for both", nAccY, nAccZ, mParams->minNumberOfAcceptedResiduals);
1254 params.flagRej.set();
1255 return 0.f;
1256 }
1257 // estimate rms on 90% of the smallest deviations
1258 int nKeepY = static_cast<int>(.9 * nAccY);
1259 int nKeepZ = static_cast<int>(.9 * nAccZ);
1260 std::nth_element(absDevY.begin(), absDevY.begin() + nKeepY, absDevY.begin() + nAccY);
1261 std::nth_element(absDevZ.begin(), absDevZ.begin() + nKeepZ, absDevZ.begin() + nAccZ);
1262 float rmsYkeep = 0.f;
1263 float rmsZkeep = 0.f;
1264 for (int i = nKeepY; i--;) {
1265 rmsYkeep += absDevY[i] * absDevY[i];
1266 }
1267 for (int i = nKeepZ; i--;) {
1268 rmsZkeep += absDevZ[i] * absDevZ[i];
1269 }
1270 rmsYkeep = std::sqrt(rmsYkeep / nKeepY);
1271 rmsZkeep = std::sqrt(rmsZkeep / nKeepZ);
1272 if (rmsYkeep < param::sEps || rmsZkeep < param::sEps) {
1273 LOG(warning) << "Too small RMS: " << rmsYkeep << "(y), " << rmsZkeep << "(z).";
1274 params.flagRej.set();
1275 return 0.f;
1276 }
1277 float rmsYkeepI = 1.f / rmsYkeep;
1278 float rmsZkeepI = 1.f / rmsZkeep;
1279 int nAcc = 0;
1280 std::array<float, param::NPadRows> yAcc;
1281 std::array<float, param::NPadRows> yDiffLong;
1282 for (int iCl = 0; iCl < nCl; ++iCl) {
1283 yDiffLL[iCl] *= rmsYkeepI;
1284 zDiffLL[iCl] *= rmsZkeepI;
1285 if (yDiffLL[iCl] * yDiffLL[iCl] + zDiffLL[iCl] * zDiffLL[iCl] > mParams->maxStdDevMA) {
1286 params.flagRej.set(iCl);
1287 } else {
1288 yAcc[nAcc++] = params.dy[iCl];
1289 }
1290 }
1291 if (nAcc > mParams->nMALong) {
1292 diffToMA(nAcc, yAcc, yDiffLong);
1293 float average = 0.f;
1294 float rms = 0.f;
1295 for (int i = 0; i < nAcc; ++i) {
1296 average += yDiffLong[i];
1297 rms += yDiffLong[i] * yDiffLong[i];
1298 }
1299 average /= nAcc;
1300 rmsLong = rms / nAcc - average * average;
1301 rmsLong = (rmsLong > 0) ? std::sqrt(rmsLong) : 0.f;
1302 }
1303 return rmsLong;
1304}
1305
1306void TrackInterpolation::diffToLocLine(const int np, int idxOffset, const std::array<float, param::NPadRows>& x, const std::array<float, param::NPadRows>& y, std::array<float, param::NPadRows>& diffY) const
1307{
1308 // Calculate the difference between the points and the linear extrapolations from the neighbourhood.
1309 // Nothing more than multiple 1-d fits at once. Instead of building 4 sums (x, x^2, y, xy), 4 * nPoints sums are calculated at once
1310 // compare to TrackResiduals::fitPoly1() method
1311
1312 // adding one entry to the vectors saves an additional if statement when calculating the cumulants
1313 std::vector<float> sumX1vec(np + 1);
1314 std::vector<float> sumX2vec(np + 1);
1315 std::vector<float> sumY1vec(np + 1);
1316 std::vector<float> sumXYvec(np + 1);
1317 auto sumX1 = &(sumX1vec[1]);
1318 auto sumX2 = &(sumX2vec[1]);
1319 auto sumY1 = &(sumY1vec[1]);
1320 auto sumXY = &(sumXYvec[1]);
1321
1322 // accumulate sums for all points
1323 for (int iCl = 0; iCl < np; ++iCl) {
1324 int idx = iCl + idxOffset;
1325 sumX1[iCl] = sumX1[iCl - 1] + x[idx];
1326 sumX2[iCl] = sumX2[iCl - 1] + x[idx] * x[idx];
1327 sumY1[iCl] = sumY1[iCl - 1] + y[idx];
1328 sumXY[iCl] = sumXY[iCl - 1] + x[idx] * y[idx];
1329 }
1330
1331 for (int iCl = 0; iCl < np; ++iCl) {
1332 int iClLeft = iCl - mParams->nMAShort;
1333 int iClRight = iCl + mParams->nMAShort;
1334 if (iClLeft < 0) {
1335 iClLeft = 0;
1336 }
1337 if (iClRight >= np) {
1338 iClRight = np - 1;
1339 }
1340 int nPoints = iClRight - iClLeft;
1341 if (nPoints < mParams->nMAShort) {
1342 continue;
1343 }
1344 float nPointsInv = 1.f / nPoints;
1345 int iClLeftP = iClLeft - 1;
1346 int iClCurrP = iCl - 1;
1347 // extract sum from iClLeft to iClRight from cumulants, excluding iCl from the fit
1348 float sX1 = sumX1[iClRight] - sumX1[iClLeftP] - (sumX1[iCl] - sumX1[iClCurrP]);
1349 float sX2 = sumX2[iClRight] - sumX2[iClLeftP] - (sumX2[iCl] - sumX2[iClCurrP]);
1350 float sY1 = sumY1[iClRight] - sumY1[iClLeftP] - (sumY1[iCl] - sumY1[iClCurrP]);
1351 float sXY = sumXY[iClRight] - sumXY[iClLeftP] - (sumXY[iCl] - sumXY[iClCurrP]);
1352 float det = sX2 - nPointsInv * sX1 * sX1;
1353 if (fabsf(det) < 1e-12f) {
1354 continue;
1355 }
1356 float slope = (sXY - nPointsInv * sX1 * sY1) / det;
1357 float offset = nPointsInv * sY1 - nPointsInv * slope * sX1;
1358 diffY[iCl + idxOffset] = y[iCl + idxOffset] - slope * x[iCl + idxOffset] - offset;
1359 }
1360}
1361
1362void TrackInterpolation::diffToMA(const int np, const std::array<float, param::NPadRows>& y, std::array<float, param::NPadRows>& diffMA) const
1363{
1364 // Calculate
1365 std::vector<float> sumVec(np + 1);
1366 auto sum = &(sumVec[1]);
1367 for (int i = 0; i < np; ++i) {
1368 sum[i] = sum[i - 1] + y[i];
1369 }
1370 for (int i = 0; i < np; ++i) {
1371 diffMA[i] = 0;
1372 int iLeft = i - mParams->nMALong;
1373 int iRight = i + mParams->nMALong;
1374 if (iLeft < 0) {
1375 iLeft = 0;
1376 }
1377 if (iRight >= np) {
1378 iRight = np - 1;
1379 }
1380 int nPoints = iRight - iLeft;
1381 if (nPoints < mParams->nMALong) {
1382 // this cannot happen, since at least mParams->nMALong points are required as neighbours for this function to be called
1383 continue;
1384 }
1385 float movingAverage = (sum[iRight] - sum[iLeft - 1] - (sum[i] - sum[i - 1])) / nPoints;
1386 diffMA[i] = y[i] - movingAverage;
1387 }
1388}
1389
1391{
1392 mTrackData.clear();
1393 mTrackDataCompact.clear();
1394 mTrackDataExtended.clear();
1395 mClRes.clear();
1396 mTrackDataUnfiltered.clear();
1397 mClResUnfiltered.clear();
1398 mGIDsSuccess.clear();
1399 for (auto& vec : mTrackIndices) {
1400 vec.clear();
1401 }
1402 mGIDs.clear();
1403 mGIDtables.clear();
1404 mTrackTimes.clear();
1405 mSeeds.clear();
1406}
1407
1408//______________________________________________
1410{
1411 // Attention! For the refit we are using reference VDrift and TDriftOffest rather than high-rate calibrated, since we want to have fixed reference over the run
1412 if (v.refVDrift != mTPCVDriftRef) {
1413 mTPCVDriftRef = v.refVDrift;
1414 mTPCDriftTimeOffsetRef = v.refTimeOffset;
1415 LOGP(info, "Imposing reference VDrift={}/TDrift={} for TPC residuals extraction", mTPCVDriftRef, mTPCDriftTimeOffsetRef);
1416 o2::tpc::TPCFastTransformHelperO2::instance()->updateCalibration(*mFastTransform, 0, 1.0, mTPCVDriftRef, mTPCDriftTimeOffsetRef);
1417 }
1418}
std::string asString(TDataMember const &dm, char *pointer)
Global TRD definitions and constants.
int32_t i
Definition of the GeometryTGeo class.
Definition of the parameter class for the detector electronics.
uint16_t pos
Definition RawData.h:3
uint16_t slope
Definition RawData.h:1
uint32_t res
Definition RawData.h:0
Definition of the TrackInterpolation class.
Definition of the TrackResiduals class.
calibration data from laser track calibration
std::ostringstream debug
Referenc on track indices contributing to the vertex, with possibility chose tracks from specific sou...
GPUd() value_type estimateLTFast(o2 static GPUd() float estimateLTIncrement(const o2 PropagatorImpl * Instance(bool uninitialized=false)
Definition Propagator.h:147
static BasicCCDBManager & instance()
T * getSpecific(std::string const &path, long timestamp=-1, MD metaData=MD(), std::map< std::string, std::string > *headers=nullptr)
retrieve an object of type T from CCDB as stored under path, timestamp and metaData
static mask_t getSourcesMask(const std::string_view srcList)
static std::string getSourcesNames(mask_t srcm)
Static class with identifiers, bitmasks and names for ALICE detectors.
Definition DetID.h:58
static constexpr ID TRD
Definition DetID.h:65
static constexpr ID TOF
Definition DetID.h:66
float getSensorRefAlpha(int isn) const
static GeometryTGeo * Instance()
float getSensorRefX(int isn) const
static void rotateToSector(Float_t *xyz, Int_t isector)
Definition Geo.cxx:1029
static void getPos(Int_t *det, Float_t *pos)
Definition Geo.cxx:491
static constexpr Int_t NPADSXSECTOR
Definition Geo.h:120
static void getVolumeIndices(Int_t index, Int_t *detId)
Definition Geo.cxx:543
static void alignedToNominalSector(Float_t *xyz, Int_t isector)
Definition Geo.cxx:1049
int updateCalibration(TPCFastTransform &transform, Long_t TimeStamp, float vDriftFactor=1.f, float vDriftRef=0.f, float driftTimeOffset=0.f)
Updates the transformation with the new time stamp.
static TPCFastTransformHelperO2 * instance()
Singleton.
void prepareInputTrackSample(const o2::globaltracking::RecoContainer &inp)
Prepare input track sample (not relying on CreateTracksVariadic functionality)
bool isInputTrackAccepted(const o2::dataformats::GlobalTrackID &gid, const o2::globaltracking::RecoContainer::GlobalIDSet &gidTable, const o2::dataformats::PrimaryVertex &pv) const
Check if input track passes configured cuts.
void diffToMA(const int np, const std::array< float, param::NPadRows > &y, std::array< float, param::NPadRows > &diffMA) const
For a given set of points, calculate their deviation from the moving average (build from the neighbou...
void diffToLocLine(const int np, int idxOffset, const std::array< float, param::NPadRows > &x, const std::array< float, param::NPadRows > &y, std::array< float, param::NPadRows > &diffY) const
For a given set of points, calculate the differences from each point to the fitted lines from all oth...
@ ExtIn
extrapolation inwards of TRD/TOF track
@ Int
interpolation (mean positions of both extrapolations)
@ ExtOut
extrapolation outwards of ITS track
o2::dataformats::GlobalTrackID::Source findValidSource(const o2::dataformats::GlobalTrackID::mask_t mask, const o2::dataformats::GlobalTrackID::Source src) const
For given vertex track source which is not in mSourcesConfigured find the seeding source which is ena...
float checkResiduals(const TrackData &trk, TrackParams &params, const std::vector< TPCClusterResiduals > &clsRes) const
bool compareToHelix(const TrackData &trk, TrackParams &params, const std::vector< TPCClusterResiduals > &clsRes) const
void process()
Main processing function.
void setTPCVDrift(const o2::tpc::VDriftCorrFact &v)
int processTRDLayer(const o2::trd::TrackTRD &trkTRD, int iLayer, o2::track::TrackParCov &trkWork, std::array< float, 2 > *trkltTRDYZ=nullptr, std::array< float, 3 > *trkltTRDCov=nullptr, TrackData *trkData=nullptr)
bool isTrackSelected(const o2::track::TrackParCov &trk) const
Given the defined downsampling factor tsalisThreshold check if track is selected.
void reset()
Reset cache and output vectors.
bool outlierFiltering(const TrackData &trk, TrackParams &params, const std::vector< TPCClusterResiduals > &clsRes) const
void init(o2::dataformats::GlobalTrackID::mask_t src, o2::dataformats::GlobalTrackID::mask_t srcMap)
Initialize everything, set the requested track sources.
bool validateTrack(const TrackData &trk, TrackParams &params, const std::vector< TPCClusterResiduals > &clsRes) const
static bool fitPoly1(int nCl, std::array< float, param::NPadRows > &x, std::array< float, param::NPadRows > &y, std::array< float, 2 > &res)
static void fitCircle(int nCl, std::array< float, param::NPadRows > &x, std::array< float, param::NPadRows > &y, float &xc, float &yc, float &r, std::array< float, param::NPadRows > &residHelixY)
static Geometry * instance()
Definition Geometry.h:33
void recalcTrkltCov(const float tilt, const float snp, const float rowSize, std::array< float, 3 > &cov) const
Recalculate tracklet covariance based on phi angle of related track.
Definition RecoParam.cxx:55
float getZCorrCoeffNRC() const
Get tracklet z correction coefficient for track-eta based corraction.
Definition RecoParam.h:49
void setBfield(float bz)
Load parameterization for given magnetic field.
Definition RecoParam.cxx:23
GLint GLenum GLint x
Definition glcorearb.h:403
GLenum src
Definition glcorearb.h:1767
GLsizeiptr size
Definition glcorearb.h:659
GLuint GLuint end
Definition glcorearb.h:469
const GLdouble * v
Definition glcorearb.h:832
GLuint GLuint GLfloat weight
Definition glcorearb.h:5477
GLenum const GLfloat * params
Definition glcorearb.h:272
GLintptr offset
Definition glcorearb.h:660
GLboolean GLboolean g
Definition glcorearb.h:1233
GLboolean r
Definition glcorearb.h:1233
GLint GLuint mask
Definition glcorearb.h:291
GLdouble GLdouble GLdouble z
Definition glcorearb.h:843
constexpr float TwoPI
constexpr float SectorSpanRad
constexpr float PI
constexpr float LightSpeedCm2S
constexpr double MassPionCharged
void convertCompactClusters(gsl::span< const itsmft::CompClusterExt > clusters, gsl::span< const unsigned char >::iterator &pattIt, std::vector< o2::BaseCluster< float > > &output, const itsmft::TopologyDictionary *dict)
convert compact clusters to 3D spacepoints
Definition IOUtils.cxx:35
int angle2Sector(float phi)
Definition Utils.h:183
float sector2Angle(int sect)
Definition Utils.h:193
constexpr int MAXGLOBALPADROW
Definition Constants.h:34
Global TPC definitions and constants.
Definition SimTraits.h:167
Enum< T >::Iterator begin(Enum< T >)
Definition Defs.h:173
DataT sum(const Vector< DataT, N > &a)
Definition Vector.h:107
TrackParCovF TrackParCov
Definition Track.h:33
constexpr int NLAYER
the number of layers
Definition Constants.h:27
constexpr int NSTACK
the number of stacks per sector
Definition Constants.h:26
GlobalIDSet getSingleDetectorRefs(GTrackID gidx) const
const o2::tpc::TrackTPC & getTPCTrack(GTrackID id) const
gsl::span< const o2::trd::CalibratedTracklet > getTRDCalibratedTracklets() const
const o2::tpc::ClusterNativeAccess & getTPCClusters() const
const o2::track::TrackParCov & getTrackParam(GTrackID gidx) const
const o2::its::TrackITS & getITSTrack(GTrackID gid) const
auto getITSTPCTRDTrack(GTrackID id) const
std::array< GTrackID, GTrackID::NSources > GlobalIDSet
gsl::span< const o2::trd::Tracklet64 > getTRDTracklets() const
const o2::dataformats::MatchInfoTOF & getTOFMatch(GTrackID id) const
static constexpr int L2G
Definition Cartesian.h:54
static constexpr int T2L
Definition Cartesian.h:55
static bool downsampleTsallisCharged(float pt, float factorPt, float sqrts, float &weight, float rnd, float mass=0.13957)
Definition Tsallis.cxx:31
float maxRejFrac
if the fraction of rejected clusters of a track is higher, the full track is invalidated
float maxDevHelixZ
max deviation in Z for clusters wrt helix fit
int nMALong
number of points to be used for moving average (long range)
float minPtNoOuterPoint
minimum pt for ITS-TPC tracks to be considered for extrapolation
int minTRDNTrklts
min number of TRD space points
int nMAShort
number of points to be used for estimation of distance from local line (short range)
float maxDevHelixY
max deviation in Y for clusters wrt helix fit
float maxStdDevMA
max cluster std. deviation (Y^2 + Z^2) wrt moving average to accept
bool skipOutlierFiltering
if set, the outlier filtering will not be applied at all
bool ignoreNonPVContrib
flag if tracks which did not contribute to the PV should be ignored or not
int minITSNClsNoOuterPoint
min number of ITS clusters if no hit in TRD or TOF exists
int minTPCNCls
min number of TPC clusters
float tsalisThreshold
in case the sampling functions returns a value smaller than this the track is discarded (1....
float maxStep
maximum step for propagation
int minTOFTRDPVContributors
min contributors from TRD or TOF (fast detectors) to consider tracks of this PV
bool enableTrackDownsampling
flag if track sampling shall be enabled or not
float maxQ2Pt
max fitted q/pt for a track to be used for calibration
float sigYZ2TOF
for now assume cluster error for TOF equal for all clusters in both Y and Z
int minNumberOfAcceptedResiduals
min number of accepted residuals for
float maxRMSLong
maximum variance of the cluster residuals wrt moving avarage for a track to be considered
int minITSNCls
min number of ITS clusters
float maxSnp
max snp when propagating tracks
bool writeUnfiltered
if set, all residuals and track parameters will be aggregated and dumped additionally without outlier...
int minNCl
min number of clusters in a track to be used for calibration
int minTPCNClsNoOuterPoint
min number of TPC clusters if no hit in TRD or TOF exists
Structure filled for each track with track quality information and a vector with TPCClusterResiduals.
float chi2TRD
chi2 of TRD track
float dEdxTPC
TPC dEdx information.
unsigned short nTrkltsTRD
number of attached TRD tracklets
unsigned short clAvailTOF
whether or not track seed has a matched TOF cluster, if so, gives the resolution of the T0 in ps
short TRDTrkltSlope[6]
TRD tracklet slope 0x7fff / param::MaxTRDSlope.
unsigned short nClsITS
number of attached ITS clusters
o2::dataformats::RangeReference clIdx
index of first cluster residual and total number of TPC cluster residuals of this track
float chi2TPC
chi2 of TPC track
unsigned short nClsTPC
number of attached TPC clusters
o2::dataformats::GlobalTrackID gid
global track ID for seeding track
float deltaTOF
TOFsignal - T0 - texp(PID), if T0 is available.
uint8_t nExtDetResid
number of external detectors (to TPC) residuals stored, on top of clIdx.getEntries
o2::track::TrackPar par
ITS track at inner TPC radius.
float chi2ITS
chi2 of ITS track
Structure for on-the-fly re-calculated track parameters at the validation stage.
unsigned char sec
TPC sector (0..35)
unsigned char row
TPC pad row.
static void init(long timestamp=-1)
short channel
extra channel info (ITS chip ID, TRD chamber, TOF main pad within the sector)
o2::mch::DsIndex ds
std::vector< o2::ctf::BufferType > vec
LOG(info)<< "Compressed in "<< sw.CpuTime()<< " s"
std::uniform_int_distribution< unsigned long long > distr
std::random_device rd
std::vector< int > row