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);
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)
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 return trkltDet;
797}
798
800{
801 // extrapolate ITS-only track through TPC and store residuals to TPC clusters in the output vectors
802 LOGP(debug, "Starting track extrapolation for GID {}", mGIDs[iSeed].asString());
803 const auto& gidTable = mGIDtables[iSeed];
804 TrackData trackData;
805 std::unique_ptr<TrackDataExtended> trackDataExtended;
806 std::vector<TPCClusterResiduals> clusterResiduals;
807 trackData.clIdx.setFirstEntry(mClRes.size());
808 const auto& trkITS = mRecoCont->getITSTrack(gidTable[GTrackID::ITS]);
809 const auto& trkTPC = mRecoCont->getTPCTrack(gidTable[GTrackID::TPC]);
810 if (mDumpTrackPoints) {
811 trackDataExtended = std::make_unique<TrackDataExtended>();
812 (*trackDataExtended).gid = mGIDs[iSeed];
813 (*trackDataExtended).clIdx.setFirstEntry(mClRes.size());
814 (*trackDataExtended).trkITS = trkITS;
815 (*trackDataExtended).trkTPC = trkTPC;
816 auto nCl = trkITS.getNumberOfClusters();
817 auto clEntry = trkITS.getFirstClusterEntry();
818 for (int iCl = nCl - 1; iCl >= 0; iCl--) { // clusters are stored from outer to inner layers
819 const auto& clsITS = mITSClustersArray[mITSTrackClusIdx[clEntry + iCl]];
820 (*trackDataExtended).clsITS.push_back(clsITS);
821 }
822 }
823 trackData.gid = mGIDs[iSeed];
824 trackData.par = mSeeds[iSeed];
825
826 auto trkWork = mSeeds[iSeed];
827 float clusterTimeBinOffset = mTrackTimes[iSeed] / mTPCTimeBinMUS;
828 auto propagator = o2::base::Propagator::Instance();
829 unsigned short rowPrev = 0; // used to calculate dRow of two consecutive cluster residuals
830 unsigned short nMeasurements = 0;
831 uint8_t clRowPrev = constants::MAXGLOBALPADROW; // used to identify and skip split clusters on the same pad row
832 for (int iCl = trkTPC.getNClusterReferences(); iCl--;) {
833 uint8_t sector, row;
834 uint32_t clusterIndexInRow;
835 const auto& cl = trkTPC.getCluster(mTPCTracksClusIdx, iCl, *mTPCClusterIdxStruct, sector, row);
836 if (clRowPrev == row) {
837 // if there are split clusters we only take the first one on the pad row
838 continue;
839 } else if (clRowPrev < constants::MAXGLOBALPADROW && clRowPrev > row) {
840 // we seem to be looping, abort this track
841 LOGP(debug, "TPC track with pT={} GeV and {} clusters has cluster {} on row {} while the previous cluster was on row {}",
842 mSeeds[iSeed].getPt(), trkTPC.getNClusterReferences(), iCl, row, clRowPrev);
843 return;
844 } else {
845 // this is the first cluster we see on this pad row
846 clRowPrev = row;
847 }
848 float x = 0, y = 0, z = 0;
849 mFastTransform->TransformIdeal(sector, row, cl.getPad(), cl.getTime(), x, y, z, clusterTimeBinOffset);
850 if (!trkWork.rotate(o2::math_utils::sector2Angle(sector))) {
851 return;
852 }
853 if (!propagator->PropagateToXBxByBz(trkWork, x, mParams->maxSnp, mParams->maxStep, mMatCorr)) {
854 return;
855 }
856
857 const auto dY = y - trkWork.getY();
858 const auto dZ = z - trkWork.getZ();
859 const auto ty = trkWork.getY();
860 const auto tz = trkWork.getZ();
861 const auto snp = trkWork.getSnp();
862 const auto sec = sector;
863
864 clusterResiduals.emplace_back(dY, dZ, ty, tz, snp, sec, row - rowPrev);
865
866 rowPrev = row;
867 ++nMeasurements;
868 }
869
870 TrackParams params; // for refitted track parameters and flagging rejected clusters
871 if (clusterResiduals.size() > constants::MAXGLOBALPADROW) {
872 LOGP(warn, "Extrapolated ITS-TPC track and found more reesiduals than possible ({})", clusterResiduals.size());
873 return;
874 }
875
876 trackData.chi2TPC = trkTPC.getChi2();
877 trackData.chi2ITS = trkITS.getChi2();
878 trackData.nClsTPC = trkTPC.getNClusterReferences();
879 trackData.nClsITS = trkITS.getNumberOfClusters();
880 trackData.clIdx.setEntries(nMeasurements);
881 trackData.dEdxTPC = trkTPC.getdEdx().dEdxTotTPC;
882 if (mDumpTrackPoints) {
883 (*trackDataExtended).trkOuter = trkWork;
884 }
885
886 if (mParams->skipOutlierFiltering || validateTrack(trackData, params, clusterResiduals)) {
887 // track is good, store TPC part
888
889 int nClValidated = 0, iRow = 0;
890 unsigned int iCl = 0;
891 for (iCl = 0; iCl < clusterResiduals.size(); ++iCl) {
892 iRow += clusterResiduals[iCl].dRow;
893 if (iRow < param::NPadRows && params.flagRej[iCl]) { // skip masked cluster residual
894 continue;
895 }
896 ++nClValidated;
897 const float tgPhi = clusterResiduals[iCl].snp / std::sqrt((1.f - clusterResiduals[iCl].snp) * (1.f + clusterResiduals[iCl].snp));
898 const auto dy = clusterResiduals[iCl].dy;
899 const auto dz = clusterResiduals[iCl].dz;
900 const auto y = clusterResiduals[iCl].y;
901 const auto z = clusterResiduals[iCl].z;
902 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)) {
903 mClRes.emplace_back(dy, dz, tgPhi, y, z, iRow, clusterResiduals[iCl].sec);
904 } else {
905 ++mRejectedResiduals;
906 }
907 }
908 trackData.clIdx.setEntries(nClValidated);
909
910 bool stopPropagation = !mExtDetResid;
911 if (!stopPropagation) {
912 // do we have TRD residuals to add?
913 int iSeedFull = mParentID[iSeed] == -1 ? iSeed : mParentID[iSeed];
914 auto gidFull = mGIDs[iSeedFull];
915 const auto& gidTableFull = mGIDtables[iSeedFull];
916 if (gidTableFull[GTrackID::TRD].isIndexSet()) {
917 const auto& trkTRD = mRecoCont->getITSTPCTRDTrack<o2::trd::TrackTRD>(gidTableFull[GTrackID::ITSTPCTRD]);
918 for (int iLayer = 0; iLayer < o2::trd::constants::NLAYER; iLayer++) {
919 std::array<float, 2> trkltTRDYZ{};
920 int res = processTRDLayer(trkTRD, iLayer, trkWork, &trkltTRDYZ);
921 if (res == -1) { // no traklet on this layer
922 continue;
923 }
924 if (res < -1) { // failed to reach this layer
925 stopPropagation = true;
926 break;
927 }
928
929 float tgPhi = trkWork.getSnp() / std::sqrt((1.f - trkWork.getSnp()) * (1.f + trkWork.getSnp()));
930 auto dy = trkltTRDYZ[0] - trkWork.getY();
931 auto dz = trkltTRDYZ[1] - trkWork.getZ();
932 const auto sec = clusterResiduals[iCl].sec;
933 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)) {
934 mClRes.emplace_back(dy, dz, tgPhi, trkWork.getY(), trkWork.getZ(), 160 + iLayer, o2::math_utils::angle2Sector(trkWork.getAlpha()), (short)res);
935 trackData.nTrkltsTRD++;
936 trackData.nExtDetResid++;
937 }
938 }
939 }
940
941 // do we have TOF residual to add?
942 trackData.clAvailTOF = 0;
943 while (gidTableFull[GTrackID::TOF].isIndexSet() && !stopPropagation) {
944 const auto& tofMatch = mRecoCont->getTOFMatch(gidFull);
945 trackData.deltaTOF = tofMatch.getSignal() - tofMatch.getFT0Best() - tofMatch.getLTIntegralOut().getTOF(trkTPC.getPID().getID());
946 trackData.clAvailTOF = uint16_t(tofMatch.getFT0BestRes());
947 const auto& clTOF = mRecoCont->getTOFClusters()[gidTableFull[GTrackID::TOF]];
948 const float clTOFAlpha = o2::math_utils::sector2Angle(clTOF.getCount());
949 float clTOFxyz[3] = {clTOF.getX(), clTOF.getY(), clTOF.getZ()};
950 if (!clTOF.isInNominalSector()) {
951 o2::tof::Geo::alignedToNominalSector(clTOFxyz, clTOF.getCount()); // go from the aligned to nominal sector frame
952 }
953 if (trkWork.getAlpha() != clTOFAlpha && !trkWork.rotate(clTOFAlpha)) {
954 LOG(debug) << "Failed to rotate into TOF cluster sector frame";
955 stopPropagation = true;
956 break;
957 }
958 if (!propagator->PropagateToXBxByBz(trkWork, clTOFxyz[0], mParams->maxSnp, mParams->maxStep, mMatCorr)) {
959 LOG(debug) << "Failed final propagation to TOF radius";
960 break;
961 }
962
963 float tgPhi = trkWork.getSnp() / std::sqrt((1.f - trkWork.getSnp()) * (1.f + trkWork.getSnp()));
964 auto dy = clTOFxyz[1] - trkWork.getY();
965 auto dz = clTOFxyz[2] - trkWork.getZ();
966 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)) {
967 mClRes.emplace_back(dy, dz, tgPhi, trkWork.getY(), trkWork.getZ(), 170, clTOF.getCount(), clTOF.getPadInSector());
968 trackData.nExtDetResid++;
969 }
970 break;
971 }
972
973 // add ITS residuals
974 while (!stopPropagation) {
975 o2::track::TrackPar trkWorkITS{trackData.par}; // this is ITS outer param
976 auto nCl = trkITS.getNumberOfClusters();
977 auto clEntry = trkITS.getFirstClusterEntry();
979 for (int iCl = 0; iCl < nCl; iCl++) { // clusters are stored from outer to inner layers
980 const auto& cls = mITSClustersArray[mITSTrackClusIdx[clEntry + iCl]];
981 int chip = cls.getSensorID();
982 float chipX, chipAlpha;
983 geom->getSensorXAlphaRefPlane(cls.getSensorID(), chipX, chipAlpha);
984 if (!trkWorkITS.rotate(chipAlpha) || !propagator->PropagateToXBxByBz(trkWorkITS, chipX, mParams->maxSnp, mParams->maxStep, mMatCorr)) {
985 LOGP(debug, "Failed final propagation to ITS X={} alpha={}", chipX, chipAlpha);
986 stopPropagation = true;
987 break;
988 }
989 float tgPhi = trkWorkITS.getSnp() / std::sqrt((1.f - trkWorkITS.getSnp()) * (1.f + trkWorkITS.getSnp()));
990 auto dy = cls.getY() - trkWorkITS.getY();
991 auto dz = cls.getZ() - trkWorkITS.getZ();
992 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)) {
993 mClRes.emplace_back(dy, dz, tgPhi, trkWorkITS.getY(), trkWorkITS.getZ(), 180 + geom->getLayer(cls.getSensorID()), -1, cls.getSensorID());
994 trackData.nExtDetResid++;
995 }
996 }
997 break;
998 }
999 }
1000 mTrackData.push_back(std::move(trackData));
1001 mGIDsSuccess.push_back(mGIDs[iSeed]);
1002 mTrackDataCompact.emplace_back(trackData.clIdx.getFirstEntry(), nClValidated, mGIDs[iSeed].getSource(), trackData.nExtDetResid);
1003 if (mDumpTrackPoints) {
1004 (*trackDataExtended).clIdx.setEntries(nClValidated);
1005 (*trackDataExtended).nExtDetResid = trackData.nExtDetResid;
1006 mTrackDataExtended.push_back(std::move(*trackDataExtended));
1007 }
1008 }
1009 if (mParams->writeUnfiltered) {
1010 TrackData trkDataTmp = trackData;
1011 trkDataTmp.clIdx.setFirstEntry(mClResUnfiltered.size());
1012 trkDataTmp.clIdx.setEntries(clusterResiduals.size());
1013 mTrackDataUnfiltered.push_back(std::move(trkDataTmp));
1014 mClResUnfiltered.insert(mClResUnfiltered.end(), clusterResiduals.begin(), clusterResiduals.end());
1015 }
1016}
1017
1018bool TrackInterpolation::validateTrack(const TrackData& trk, TrackParams& params, const std::vector<TPCClusterResiduals>& clsRes) const
1019{
1020 if (clsRes.size() < mParams->minNCl) {
1021 // no enough clusters for this track to be considered
1022 LOG(debug) << "Skipping track with too few clusters: " << clsRes.size();
1023 return false;
1024 }
1025
1026 bool resHelix = compareToHelix(trk, params, clsRes);
1027 if (!resHelix) {
1028 LOG(debug) << "Skipping track too far from helix approximation";
1029 return false;
1030 }
1031 if (fabsf(mBz) > 0.01 && fabsf(params.qpt) > mParams->maxQ2Pt) {
1032 LOG(debug) << "Skipping track with too high q/pT: " << params.qpt;
1033 return false;
1034 }
1035 if (!outlierFiltering(trk, params, clsRes)) {
1036 return false;
1037 }
1038 return true;
1039}
1040
1041bool TrackInterpolation::compareToHelix(const TrackData& trk, TrackParams& params, const std::vector<TPCClusterResiduals>& clsRes) const
1042{
1043 std::array<float, param::NPadRows> residHelixY;
1044 std::array<float, param::NPadRows> residHelixZ;
1045
1046 std::array<float, param::NPadRows> xLab;
1047 std::array<float, param::NPadRows> yLab;
1048 std::array<float, param::NPadRows> sPath;
1049
1050 float curvature = fabsf(trk.par.getQ2Pt() * mBz * o2::constants::physics::LightSpeedCm2S * 1e-14f);
1051 int secFirst = clsRes[0].sec;
1052 float phiSect = (secFirst + .5f) * o2::constants::math::SectorSpanRad;
1053 float snPhi = sin(phiSect);
1054 float csPhi = cos(phiSect);
1055 sPath[0] = 0.f;
1056
1057 int iRow = 0;
1058 int nCl = clsRes.size();
1059 for (unsigned int iP = 0; iP < nCl; ++iP) {
1060 iRow += clsRes[iP].dRow;
1061 float yTrk = clsRes[iP].y;
1062 // LOGF(info, "iRow(%i), yTrk(%f)", iRow, yTrk);
1063 xLab[iP] = param::RowX[iRow];
1064 if (clsRes[iP].sec != secFirst) {
1065 float phiSectCurrent = (clsRes[iP].sec + .5f) * o2::constants::math::SectorSpanRad;
1066 float cs = cos(phiSectCurrent - phiSect);
1067 float sn = sin(phiSectCurrent - phiSect);
1068 xLab[iP] = param::RowX[iRow] * cs - yTrk * sn;
1069 yLab[iP] = yTrk * cs + param::RowX[iRow] * sn;
1070 } else {
1071 xLab[iP] = param::RowX[iRow];
1072 yLab[iP] = yTrk;
1073 }
1074 // this is needed only later, but we retrieve it already now to save another loop
1075 params.zTrk[iP] = clsRes[iP].z;
1076 params.xTrk[iP] = param::RowX[iRow];
1077 params.dy[iP] = clsRes[iP].dy;
1078 params.dz[iP] = clsRes[iP].dz;
1079 // done retrieving values for later
1080 if (iP > 0) {
1081 float dx = xLab[iP] - xLab[iP - 1];
1082 float dy = yLab[iP] - yLab[iP - 1];
1083 float ds2 = dx * dx + dy * dy;
1084 float ds = sqrt(ds2); // circular path (linear approximation)
1085 // if the curvature of the track or the (approximated) chord length is too large the more exact formula is used:
1086 // chord length = 2r * asin(ds/(2r))
1087 // using the first two terms of the tailer expansion for asin(x) ~ x + x^3 / 6
1088 if (ds * curvature > 0.05) {
1089 ds *= (1.f + ds2 * curvature * curvature / 24.f);
1090 }
1091 sPath[iP] = sPath[iP - 1] + ds;
1092 }
1093 }
1094 if (fabsf(mBz) < 0.01) {
1095 // for B=0 we don't need to try a circular fit...
1096 return true;
1097 }
1098 float xcSec = 0.f;
1099 float ycSec = 0.f;
1100 float r = 0.f;
1101 TrackResiduals::fitCircle(nCl, xLab, yLab, xcSec, ycSec, r, residHelixY);
1102 // LOGF(info, "Done with circle fit. nCl(%i), xcSec(%f), ycSec(%f), r(%f).", nCl, xcSec, ycSec, r);
1103 /*
1104 for (int i=0; i<nCl; ++i) {
1105 LOGF(info, "i(%i), xLab(%f), yLab(%f).", i, xLab[i], yLab[i]);
1106 }
1107 */
1108 // determine curvature
1109 float phiI = TMath::ATan2(yLab[0], xLab[0]);
1110 float phiF = TMath::ATan2(yLab[nCl - 1], xLab[nCl - 1]);
1111 if (phiI < 0) {
1113 }
1114 if (phiF < 0) {
1116 }
1117 float dPhi = phiF - phiI;
1118 float curvSign = -1.f;
1119 if (dPhi > 0) {
1120 if (dPhi < o2::constants::math::PI) {
1121 curvSign = 1.f;
1122 }
1123 } else if (dPhi < -o2::constants::math::PI) {
1124 curvSign = 1.f;
1125 }
1126 params.qpt = std::copysign(1.f / (r * mBz * o2::constants::physics::LightSpeedCm2S * 1e-14f), curvSign);
1127
1128 // calculate circle coordinates in the lab frame
1129 float xc = xcSec * csPhi - ycSec * snPhi;
1130 float yc = xcSec * snPhi + ycSec * csPhi;
1131
1132 std::array<float, 2> pol1Z;
1133 TrackResiduals::fitPoly1(nCl, sPath, params.zTrk, pol1Z);
1134
1135 params.tgl = pol1Z[0];
1136
1137 // max deviations in both directions from helix fit in y and z
1138 float hMinY = 1e9f;
1139 float hMaxY = -1e9f;
1140 float hMinZ = 1e9f;
1141 float hMaxZ = -1e9f;
1142 // extract residuals in Z and fill track slopes in sector frame
1143 int secCurr = secFirst;
1144 iRow = 0;
1145 for (unsigned int iCl = 0; iCl < nCl; ++iCl) {
1146 iRow += clsRes[iCl].dRow;
1147 float resZ = params.zTrk[iCl] - (pol1Z[1] + sPath[iCl] * pol1Z[0]);
1148 residHelixZ[iCl] = resZ;
1149 if (resZ < hMinZ) {
1150 hMinZ = resZ;
1151 }
1152 if (resZ > hMaxZ) {
1153 hMaxZ = resZ;
1154 }
1155 if (residHelixY[iCl] < hMinY) {
1156 hMinY = residHelixY[iCl];
1157 }
1158 if (residHelixY[iCl] > hMaxY) {
1159 hMaxY = residHelixY[iCl];
1160 }
1161 int sec = clsRes[iCl].sec;
1162 if (sec != secCurr) {
1163 secCurr = sec;
1164 phiSect = (.5f + sec) * o2::constants::math::SectorSpanRad;
1165 snPhi = sin(phiSect);
1166 csPhi = cos(phiSect);
1167 xcSec = xc * csPhi + yc * snPhi; // recalculate circle center in the sector frame
1168 }
1169 float cstalp = (param::RowX[iRow] - xcSec) / r;
1170 if (fabsf(cstalp) > 1.f - sFloatEps) {
1171 // track cannot reach this pad row
1172 cstalp = std::copysign(1.f - sFloatEps, cstalp);
1173 }
1174 params.tglArr[iCl] = cstalp / sqrt((1 - cstalp) * (1 + cstalp)); // 1 / tan(acos(cstalp)) = cstalp / sqrt(1 - cstalp^2)
1175
1176 // In B+ the slope of q- should increase with x. Just look on q * B
1177 if (params.qpt * mBz > 0) {
1178 params.tglArr[iCl] *= -1.f;
1179 }
1180 }
1181 // 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);
1182 // LOGF(info, "New pt/Q (%f), old pt/Q (%f)", 1./params.qpt, 1./trk.qPt);
1183 return fabsf(hMaxY - hMinY) < mParams->maxDevHelixY && fabsf(hMaxZ - hMinZ) < mParams->maxDevHelixZ;
1184}
1185
1186bool TrackInterpolation::outlierFiltering(const TrackData& trk, TrackParams& params, const std::vector<TPCClusterResiduals>& clsRes) const
1187{
1188 if (clsRes.size() < mParams->nMALong) {
1189 LOG(debug) << "Skipping track with too few clusters for long moving average: " << clsRes.size();
1190 return false;
1191 }
1192 float rmsLong = checkResiduals(trk, params, clsRes);
1193 if (static_cast<float>(params.flagRej.count()) / clsRes.size() > mParams->maxRejFrac) {
1194 LOGP(debug, "Skipping track with too many clusters rejected: {} out of {}", params.flagRej.count(), clsRes.size());
1195 return false;
1196 }
1197 if (rmsLong > mParams->maxRMSLong) {
1198 LOG(debug) << "Skipping track with too large RMS: " << rmsLong;
1199 return false;
1200 }
1201 return true;
1202}
1203
1204float TrackInterpolation::checkResiduals(const TrackData& trk, TrackParams& params, const std::vector<TPCClusterResiduals>& clsRes) const
1205{
1206 float rmsLong = 0.f;
1207
1208 int nCl = clsRes.size();
1209 int iClFirst = 0;
1210 int iClLast = nCl - 1;
1211 int secStart = clsRes[0].sec;
1212
1213 // arrays with differences / abs(differences) of points to their neighbourhood, initialized to zero
1214 std::array<float, param::NPadRows> yDiffLL{};
1215 std::array<float, param::NPadRows> zDiffLL{};
1216 std::array<float, param::NPadRows> absDevY{};
1217 std::array<float, param::NPadRows> absDevZ{};
1218
1219 for (unsigned int iCl = 0; iCl < nCl; ++iCl) {
1220 if (iCl < iClLast && clsRes[iCl].sec == secStart) {
1221 continue;
1222 }
1223 // sector changed or last cluster reached
1224 // now run estimators for all points in the same sector
1225 int nClSec = iCl - iClFirst;
1226 if (iCl == iClLast) {
1227 ++nClSec;
1228 }
1229 diffToLocLine(nClSec, iClFirst, params.xTrk, params.dy, yDiffLL);
1230 diffToLocLine(nClSec, iClFirst, params.xTrk, params.dz, zDiffLL);
1231 iClFirst = iCl;
1232 secStart = clsRes[iCl].sec;
1233 }
1234 // store abs deviations
1235 int nAccY = 0;
1236 int nAccZ = 0;
1237 for (int iCl = nCl; iCl--;) {
1238 if (fabsf(yDiffLL[iCl]) > param::sEps) {
1239 absDevY[nAccY++] = fabsf(yDiffLL[iCl]);
1240 }
1241 if (fabsf(zDiffLL[iCl]) > param::sEps) {
1242 absDevZ[nAccZ++] = fabsf(zDiffLL[iCl]);
1243 }
1244 }
1245 if (nAccY < mParams->minNumberOfAcceptedResiduals || nAccZ < mParams->minNumberOfAcceptedResiduals) {
1246 // mask all clusters
1247 LOGP(debug, "Accepted {} clusters for dY {} clusters for dZ, but required at least {} for both", nAccY, nAccZ, mParams->minNumberOfAcceptedResiduals);
1248 params.flagRej.set();
1249 return 0.f;
1250 }
1251 // estimate rms on 90% of the smallest deviations
1252 int nKeepY = static_cast<int>(.9 * nAccY);
1253 int nKeepZ = static_cast<int>(.9 * nAccZ);
1254 std::nth_element(absDevY.begin(), absDevY.begin() + nKeepY, absDevY.begin() + nAccY);
1255 std::nth_element(absDevZ.begin(), absDevZ.begin() + nKeepZ, absDevZ.begin() + nAccZ);
1256 float rmsYkeep = 0.f;
1257 float rmsZkeep = 0.f;
1258 for (int i = nKeepY; i--;) {
1259 rmsYkeep += absDevY[i] * absDevY[i];
1260 }
1261 for (int i = nKeepZ; i--;) {
1262 rmsZkeep += absDevZ[i] * absDevZ[i];
1263 }
1264 rmsYkeep = std::sqrt(rmsYkeep / nKeepY);
1265 rmsZkeep = std::sqrt(rmsZkeep / nKeepZ);
1266 if (rmsYkeep < param::sEps || rmsZkeep < param::sEps) {
1267 LOG(warning) << "Too small RMS: " << rmsYkeep << "(y), " << rmsZkeep << "(z).";
1268 params.flagRej.set();
1269 return 0.f;
1270 }
1271 float rmsYkeepI = 1.f / rmsYkeep;
1272 float rmsZkeepI = 1.f / rmsZkeep;
1273 int nAcc = 0;
1274 std::array<float, param::NPadRows> yAcc;
1275 std::array<float, param::NPadRows> yDiffLong;
1276 for (int iCl = 0; iCl < nCl; ++iCl) {
1277 yDiffLL[iCl] *= rmsYkeepI;
1278 zDiffLL[iCl] *= rmsZkeepI;
1279 if (yDiffLL[iCl] * yDiffLL[iCl] + zDiffLL[iCl] * zDiffLL[iCl] > mParams->maxStdDevMA) {
1280 params.flagRej.set(iCl);
1281 } else {
1282 yAcc[nAcc++] = params.dy[iCl];
1283 }
1284 }
1285 if (nAcc > mParams->nMALong) {
1286 diffToMA(nAcc, yAcc, yDiffLong);
1287 float average = 0.f;
1288 float rms = 0.f;
1289 for (int i = 0; i < nAcc; ++i) {
1290 average += yDiffLong[i];
1291 rms += yDiffLong[i] * yDiffLong[i];
1292 }
1293 average /= nAcc;
1294 rmsLong = rms / nAcc - average * average;
1295 rmsLong = (rmsLong > 0) ? std::sqrt(rmsLong) : 0.f;
1296 }
1297 return rmsLong;
1298}
1299
1300void 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
1301{
1302 // Calculate the difference between the points and the linear extrapolations from the neighbourhood.
1303 // 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
1304 // compare to TrackResiduals::fitPoly1() method
1305
1306 // adding one entry to the vectors saves an additional if statement when calculating the cumulants
1307 std::vector<float> sumX1vec(np + 1);
1308 std::vector<float> sumX2vec(np + 1);
1309 std::vector<float> sumY1vec(np + 1);
1310 std::vector<float> sumXYvec(np + 1);
1311 auto sumX1 = &(sumX1vec[1]);
1312 auto sumX2 = &(sumX2vec[1]);
1313 auto sumY1 = &(sumY1vec[1]);
1314 auto sumXY = &(sumXYvec[1]);
1315
1316 // accumulate sums for all points
1317 for (int iCl = 0; iCl < np; ++iCl) {
1318 int idx = iCl + idxOffset;
1319 sumX1[iCl] = sumX1[iCl - 1] + x[idx];
1320 sumX2[iCl] = sumX2[iCl - 1] + x[idx] * x[idx];
1321 sumY1[iCl] = sumY1[iCl - 1] + y[idx];
1322 sumXY[iCl] = sumXY[iCl - 1] + x[idx] * y[idx];
1323 }
1324
1325 for (int iCl = 0; iCl < np; ++iCl) {
1326 int iClLeft = iCl - mParams->nMAShort;
1327 int iClRight = iCl + mParams->nMAShort;
1328 if (iClLeft < 0) {
1329 iClLeft = 0;
1330 }
1331 if (iClRight >= np) {
1332 iClRight = np - 1;
1333 }
1334 int nPoints = iClRight - iClLeft;
1335 if (nPoints < mParams->nMAShort) {
1336 continue;
1337 }
1338 float nPointsInv = 1.f / nPoints;
1339 int iClLeftP = iClLeft - 1;
1340 int iClCurrP = iCl - 1;
1341 // extract sum from iClLeft to iClRight from cumulants, excluding iCl from the fit
1342 float sX1 = sumX1[iClRight] - sumX1[iClLeftP] - (sumX1[iCl] - sumX1[iClCurrP]);
1343 float sX2 = sumX2[iClRight] - sumX2[iClLeftP] - (sumX2[iCl] - sumX2[iClCurrP]);
1344 float sY1 = sumY1[iClRight] - sumY1[iClLeftP] - (sumY1[iCl] - sumY1[iClCurrP]);
1345 float sXY = sumXY[iClRight] - sumXY[iClLeftP] - (sumXY[iCl] - sumXY[iClCurrP]);
1346 float det = sX2 - nPointsInv * sX1 * sX1;
1347 if (fabsf(det) < 1e-12f) {
1348 continue;
1349 }
1350 float slope = (sXY - nPointsInv * sX1 * sY1) / det;
1351 float offset = nPointsInv * sY1 - nPointsInv * slope * sX1;
1352 diffY[iCl + idxOffset] = y[iCl + idxOffset] - slope * x[iCl + idxOffset] - offset;
1353 }
1354}
1355
1356void TrackInterpolation::diffToMA(const int np, const std::array<float, param::NPadRows>& y, std::array<float, param::NPadRows>& diffMA) const
1357{
1358 // Calculate
1359 std::vector<float> sumVec(np + 1);
1360 auto sum = &(sumVec[1]);
1361 for (int i = 0; i < np; ++i) {
1362 sum[i] = sum[i - 1] + y[i];
1363 }
1364 for (int i = 0; i < np; ++i) {
1365 diffMA[i] = 0;
1366 int iLeft = i - mParams->nMALong;
1367 int iRight = i + mParams->nMALong;
1368 if (iLeft < 0) {
1369 iLeft = 0;
1370 }
1371 if (iRight >= np) {
1372 iRight = np - 1;
1373 }
1374 int nPoints = iRight - iLeft;
1375 if (nPoints < mParams->nMALong) {
1376 // this cannot happen, since at least mParams->nMALong points are required as neighbours for this function to be called
1377 continue;
1378 }
1379 float movingAverage = (sum[iRight] - sum[iLeft - 1] - (sum[i] - sum[i - 1])) / nPoints;
1380 diffMA[i] = y[i] - movingAverage;
1381 }
1382}
1383
1385{
1386 mTrackData.clear();
1387 mTrackDataCompact.clear();
1388 mTrackDataExtended.clear();
1389 mClRes.clear();
1390 mTrackDataUnfiltered.clear();
1391 mClResUnfiltered.clear();
1392 mGIDsSuccess.clear();
1393 for (auto& vec : mTrackIndices) {
1394 vec.clear();
1395 }
1396 mGIDs.clear();
1397 mGIDtables.clear();
1398 mTrackTimes.clear();
1399 mSeeds.clear();
1400}
1401
1402//______________________________________________
1404{
1405 // 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
1406 if (v.refVDrift != mTPCVDriftRef) {
1407 mTPCVDriftRef = v.refVDrift;
1408 mTPCDriftTimeOffsetRef = v.refTimeOffset;
1409 LOGP(info, "Imposing reference VDrift={}/TDrift={} for TPC residuals extraction", mTPCVDriftRef, mTPCDriftTimeOffsetRef);
1410 o2::tpc::TPCFastTransformHelperO2::instance()->updateCalibration(*mFastTransform, 0, 1.0, mTPCVDriftRef, mTPCDriftTimeOffsetRef);
1411 }
1412}
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:143
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.
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)
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)
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
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