Project
Loading...
Searching...
No Matches
MatchTOF.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#include <TTree.h>
12#include <cassert>
13
14#include <fairlogger/Logger.h>
15#include "Field/MagneticField.h"
16#include "Field/MagFieldFast.h"
17#include "TOFBase/Geo.h"
18
20
23#include "MathUtils/Cartesian.h"
24#include "MathUtils/Utils.h"
29
30#include <Math/SMatrix.h>
31#include <Math/SVector.h>
32#include <TFile.h>
33#include <TGeoGlobalMagField.h>
38
40
44
47#include "TOFBase/Utils.h"
48
49#ifdef WITH_OPENMP
50#include <omp.h>
51#endif
52
53using namespace o2::globaltracking;
60
61bool MatchTOF::mHasFillScheme = false;
62bool MatchTOF::mFillScheme[o2::constants::lhc::LHCMaxBunches] = {0};
63
65
66//______________________________________________
67void MatchTOF::run(const o2::globaltracking::RecoContainer& inp, unsigned long firstTForbit)
68{
69 mFirstTForbit = firstTForbit;
70
71 if (!mMatchParams) {
73 mSigmaTimeCut = mMatchParams->nsigmaTimeCut;
74 }
75
77 mRecoCont = &inp;
78 mStartIR = inp.startIR;
79 updateTimeDependentParams();
80
81 mTimerMatchTPC.Reset();
82 mTimerMatchITSTPC.Reset();
83 mTimerTot.Reset();
84
85 mCalibInfoTOF.clear();
86
87 for (int i = 0; i < trkType::SIZEALL; i++) {
88 mMatchedTracks[i].clear();
89 mOutTOFLabels[i].clear();
90 }
91
92 for (int it = 0; it < trkType::SIZE; it++) {
93 for (int sec = o2::constants::math::NSectors - 1; sec > -1; sec--) {
94 mMatchedTracksIndex[sec][it].clear();
95 mTracksWork[sec][it].clear();
96 mTrackGid[sec][it].clear();
97 mLTinfos[sec][it].clear();
98 if (mMCTruthON) {
99 mTracksLblWork[sec][it].clear();
100 }
101 mTracksSectIndexCache[it][sec].clear();
102 mTracksSeed[it][sec].clear();
103 }
104 }
105
106 for (int sec = o2::constants::math::NSectors - 1; sec > -1; sec--) {
107 mSideTPC[sec].clear();
108 mExtraTPCFwdTime[sec].clear();
109 }
110
111 mTimerTot.Start();
112 bool isPrepareTOFClusters = prepareTOFClusters();
113 mTimerTot.Stop();
114 LOGF(info, "Timing prepareTOFCluster: Cpu: %.3e s Real: %.3e s in %d slots", mTimerTot.CpuTime(), mTimerTot.RealTime(), mTimerTot.Counter() - 1);
115
116 if (!isPrepareTOFClusters) { // check cluster before of tracks to see also if MC is required
117 return;
118 }
119
120 mTimerTot.Start();
121 if (!prepareTPCData()) {
122 return;
123 }
124 mTimerTot.Stop();
125 LOGF(info, "Timing prepare TPC tracks: Cpu: %.3e s Real: %.3e s in %d slots", mTimerTot.CpuTime(), mTimerTot.RealTime(), mTimerTot.Counter() - 1);
126
127 mTimerTot.Start();
128 if (!prepareFITData()) {
129 return;
130 }
131 mTimerTot.Stop();
132 LOGF(info, "Timing prepare FIT data: Cpu: %.3e s Real: %.3e s in %d slots", mTimerTot.CpuTime(), mTimerTot.RealTime(), mTimerTot.Counter() - 1);
133
134 mTimerTot.Start();
135 std::array<uint32_t, 18> nMatches = {0};
136 // run matching (this can be multi-thread)
137 mTimerMatchITSTPC.Start();
138
139 for (int sec = o2::constants::math::NSectors - 1; sec > -1; sec--) {
140 mMatchedTracksPairsSec[sec].clear(); // new sector
141 }
142
144
145 if (mIsITSTPCused || mIsTPCTRDused || mIsITSTPCTRDused) {
146#ifdef WITH_OPENMP
147#pragma omp parallel for schedule(dynamic) num_threads(mNlanes)
148#endif
149 for (int sec = o2::constants::math::NSectors - 1; sec > -1; sec--) {
150 doMatching(sec);
151 }
152 }
153 mTimerMatchITSTPC.Stop();
154
155 mTimerMatchTPC.Start();
156 if (mIsTPCused) {
157#ifdef WITH_OPENMP
158#pragma omp parallel for schedule(dynamic) num_threads(mNlanes)
159#endif
160 for (int sec = o2::constants::math::NSectors - 1; sec > -1; sec--) {
161 doMatchingForTPC(sec);
162 }
163 }
164 mTimerMatchTPC.Stop();
165
166 // finalize
167 for (int sec = o2::constants::math::NSectors - 1; sec > -1; sec--) {
168 if (mStoreMatchable) {
169 // if MC check if good or fake matches
170 if (mMCTruthON) {
171 for (auto& matchingPair : mMatchedTracksPairsSec[sec]) {
172 int trkType = (int)matchingPair.getTrackType();
173 int itrk = matchingPair.getIdLocal();
174 const auto& labelsTOF = mTOFClusLabels->getLabels(matchingPair.getTOFClIndex());
175 const auto& labelTrack = mTracksLblWork[sec][trkType][itrk];
176
177 // we have not found the track label among those associated to the TOF cluster --> fake match! We will associate the label of the main channel, but negative
178 bool fake = true;
179 for (auto& lbl : labelsTOF) {
180 if (labelTrack == lbl) { // compares src, evID, trID, ignores fake flag.
181 fake = false;
182 }
183 }
184 if (fake) {
185 matchingPair.setFakeMatch();
186 }
187 }
188 } else {
189 for (auto& matchingPair : mMatchedTracksPairsSec[sec]) {
190 int bct0 = int((matchingPair.getSignal() - matchingPair.getLTIntegralOut().getTOF(0) + 5000) * Geo::BC_TIME_INPS_INV); // bc taken assuming speed of light (el) and 5 ns of margin
191 if (bct0 < 0) { // if negative time (it can happen at the beginng of the TF int was truncated per excess... adjusting)
192 bct0--;
193 }
194 float tof = matchingPair.getSignal() - bct0 * Geo::BC_TIME_INPS;
195 if (std::abs(tof - matchingPair.getLTIntegralOut().getTOF(2)) < 600) {
196 } else if (std::abs(tof - matchingPair.getLTIntegralOut().getTOF(3)) < 600) {
197 } else if (std::abs(tof - matchingPair.getLTIntegralOut().getTOF(4)) < 600) {
198 } else if (std::abs(tof - matchingPair.getLTIntegralOut().getTOF(0)) < 600) {
199 } else if (std::abs(tof - matchingPair.getLTIntegralOut().getTOF(1)) < 600) {
200 } else if (std::abs(tof - matchingPair.getLTIntegralOut().getTOF(5)) < 600) {
201 } else if (std::abs(tof - matchingPair.getLTIntegralOut().getTOF(6)) < 600) {
202 } else if (std::abs(tof - matchingPair.getLTIntegralOut().getTOF(7)) < 600) {
203 } else if (std::abs(tof - matchingPair.getLTIntegralOut().getTOF(8)) < 600) {
204 } else { // no pion, kaon, proton, electron, muon, deuteron, triton, 3He, 4He
205 matchingPair.setFakeMatch();
206 }
207 }
208 }
209 }
210
211 LOG(debug) << "...done. Now check the best matches";
212 nMatches[sec] = mMatchedTracksPairsSec[sec].size();
213 selectBestMatches(sec);
214
215 if (mStoreMatchable) {
216 for (auto& matchingPair : mMatchedTracksPairsSec[sec]) {
217 trkType trkTypeSplitted = trkType::TPC;
218 auto sourceID = matchingPair.getTrackRef().getSource();
220 trkTypeSplitted = trkType::ITSTPC;
221 } else if (sourceID == o2::dataformats::GlobalTrackID::TPCTRD) {
222 trkTypeSplitted = trkType::TPCTRD;
223 } else if (sourceID == o2::dataformats::GlobalTrackID::ITSTPCTRD) {
224 trkTypeSplitted = trkType::ITSTPCTRD;
225 }
226 matchingPair.setTrackType(trkTypeSplitted);
227 }
228 }
229 }
230 std::string nMatchesStr = "Number of pairs matched per sector: ";
231 for (int sec = o2::constants::math::NSectors - 1; sec > -1; sec--) {
232 nMatchesStr += fmt::format("{} : {} ; ", sec, nMatches[sec]);
233 }
234 LOG(info) << nMatchesStr;
235 // re-arrange outputs from constrained/unconstrained to the 4 cases (TPC, ITS-TPC, TPC-TRD, ITS-TPC-TRD) to be implemented as soon as TPC-TRD and ITS-TPC-TRD tracks available
236
237 mIsTPCused = false;
238 mIsITSTPCused = false;
239 mIsTPCTRDused = false;
240 mIsITSTPCTRDused = false;
241
242 mTimerTot.Stop();
243 LOGF(info, "Timing Do Matching: Cpu: %.3e s Real: %.3e s in %d slots", mTimerTot.CpuTime(), mTimerTot.RealTime(), mTimerTot.Counter() - 1);
244 LOGF(info, "Timing Do Matching Constrained: Cpu: %.3e s Real: %.3e s in %d slots", mTimerMatchITSTPC.CpuTime(), mTimerMatchITSTPC.RealTime(), mTimerMatchITSTPC.Counter() - 1);
245 LOGF(info, "Timing Do Matching TPC : Cpu: %.3e s Real: %.3e s in %d slots", mTimerMatchTPC.CpuTime(), mTimerMatchTPC.RealTime(), mTimerMatchTPC.Counter() - 1);
246}
247
248//______________________________________________
250{
251 mTPCVDrift = v.refVDrift * v.corrFact;
252 mTPCVDriftCorrFact = v.corrFact;
253 mTPCVDriftRef = v.refVDrift;
254 mTPCDriftTimeOffset = v.getTimeOffset();
255}
256
257//______________________________________________
259{
260 mTPCCorrMapsHelper = maph;
261}
262
263//______________________________________________
264void MatchTOF::print() const
265{
267
268 LOG(info) << "****** component for the matching of tracks to TOF clusters ******";
269
270 LOG(info) << "MC truth: " << (mMCTruthON ? "on" : "off");
271 LOG(info) << "Time tolerance: " << mTimeTolerance;
272 LOG(info) << "Space tolerance: " << mSpaceTolerance;
273 LOG(info) << "SigmaTimeCut: " << mSigmaTimeCut;
274
275 LOG(info) << "**********************************************************************";
276}
277//______________________________________________
279{
281}
282//_____________________________________________________
283bool MatchTOF::prepareFITData()
284{
285 // If available, read FIT Info
286 if (mIsFIT) {
287 mFITRecPoints = mRecoCont->getFT0RecPoints();
288 // prepareInteractionTimes();
289 }
290 return true;
291}
292//______________________________________________
293int MatchTOF::prepareInteractionTimes()
294{
295 // do nothing. If you think it can be useful have a look at MatchTPCITS
296 return 0;
297}
298//______________________________________________
299void MatchTOF::printGrouping(const std::vector<o2::dataformats::MatchInfoTOFReco>& origin, const std::vector<std::vector<o2::dataformats::MatchInfoTOFReco>>& grouped)
300{
301 printf("Original vector\n");
302 for (const auto& seed : origin) {
303 printf("Pair: track=%d TOFcl=%d\n", seed.getIdLocal(), seed.getTOFClIndex());
304 }
305
306 printf("\nGroups\n");
307 int ngroups = 0;
308 for (const auto& gr : grouped) {
309 ngroups++;
310 printf("Group %d\n", ngroups);
311 for (const auto& seed : gr) {
312 printf("Pair: track=%d TOFcl=%d\n", seed.getIdLocal(), seed.getTOFClIndex());
313 }
314 }
315}
316//______________________________________________
317void MatchTOF::groupingMatch(const std::vector<o2::dataformats::MatchInfoTOFReco>& origin, std::vector<std::vector<o2::dataformats::MatchInfoTOFReco>>& grouped, std::vector<std::vector<int>>& firstEls, std::vector<std::vector<int>>& secondEls)
318{
319 grouped.clear();
320 firstEls.clear();
321 secondEls.clear();
322
323 std::vector<o2::dataformats::MatchInfoTOFReco> dummy;
324 std::vector<int> dummyInt;
325
326 int pos = 0;
327 int posInVector = 0;
328 std::vector<int> alreadyUsed(origin.size(), 0);
329 while (posInVector < origin.size()) { // go ahead if there are elements
330 if (alreadyUsed[posInVector]) {
331 posInVector++;
332 continue;
333 }
334 bool found = true;
335 grouped.push_back(dummy);
336 auto& trkId = firstEls.emplace_back(dummyInt);
337 auto& cluId = secondEls.emplace_back(dummyInt);
338
339 // first element is the seed
340 const auto& seed = origin[posInVector];
341 trkId.push_back(seed.getIdLocal());
342 cluId.push_back(seed.getTOFClIndex());
343 // remove element added to the current group
344 grouped[pos].push_back(seed);
345 alreadyUsed[posInVector] = true;
346 posInVector++;
347
348 while (found) {
349 found = false;
350 for (int i = posInVector; i < origin.size(); i++) {
351 if (alreadyUsed[i]) {
352 continue;
353 }
354 const auto& seed = origin[i];
355 int matchFirst = -1;
356 int matchSecond = -1;
357 for (const int& ind : trkId) {
358 if (seed.getIdLocal() == ind) {
359 matchFirst = ind;
360 break;
361 }
362 }
363 for (const int& ind : cluId) {
364 if (seed.getTOFClIndex() == ind) {
365 matchSecond = ind;
366 break;
367 }
368 }
369
370 if (matchFirst >= 0 || matchSecond >= 0) { // belong to this group
371 if (matchFirst < 0) {
372 trkId.push_back(seed.getIdLocal());
373 }
374 if (matchSecond < 0) {
375 cluId.push_back(seed.getTOFClIndex());
376 }
377 grouped[pos].push_back(seed);
378 alreadyUsed[i] = true;
379 found = true;
380 break;
381 }
382 }
383 }
384 pos++; // move to the next group
385 }
386}
387
388//______________________________________________
389bool MatchTOF::prepareTPCData()
390{
391 mNotPropagatedToTOF[trkType::UNCONS] = 0;
392 mNotPropagatedToTOF[trkType::CONSTR] = 0;
393
394 auto creator = [this](auto& trk, GTrackID gid, float time0, float terr) {
395 const int nclustersMin = 0;
396 if constexpr (isTPCTrack<decltype(trk)>()) {
397 if (trk.getNClusters() < nclustersMin) {
398 return true;
399 }
400
401 if (std::abs(trk.getQ2Pt()) > mMaxInvPt) {
402 return true;
403 }
404 this->addTPCSeed(trk, gid, time0, terr);
405 }
406 if constexpr (isTPCITSTrack<decltype(trk)>()) {
407 if (trk.getParamOut().getX() < o2::constants::geom::XTPCOuterRef - 1.) {
408 return true;
409 }
410 this->addITSTPCSeed(trk, gid, time0, terr);
411 }
412 if constexpr (isTRDTrack<decltype(trk)>()) {
413 this->addTRDSeed(trk, gid, time0, terr);
414 }
415 return true;
416 };
417 mRecoCont->createTracksVariadic(creator);
418
419#ifdef WITH_OPENMP
420#pragma omp parallel for schedule(dynamic) num_threads(mNlanes)
421#endif
422 for (int sec = o2::constants::math::NSectors - 1; sec > -1; sec--) {
423 if (mIsTPCused) {
424 propagateTPCTracks(sec);
425 }
426 if (mIsITSTPCused || mIsTPCTRDused || mIsITSTPCTRDused) {
427 propagateConstrTracks(sec);
428 }
429 }
430
431 // re-assign tracks which change sector (no multi-threading)
432 std::array<float, 3> globalPos;
433 for (int sec = o2::constants::math::NSectors - 1; sec > -1; sec--) {
434 if (mIsTPCused) {
435 for (auto& it : mTracksSeed[trkType::UNCONS][sec]) {
436 auto& pair = mTracksWork[sec][trkType::UNCONS][it];
437 auto& trc = pair.first;
438 trc.getXYZGlo(globalPos);
439 int sector = o2::math_utils::angle2Sector(TMath::ATan2(globalPos[1], globalPos[0]));
440
441 int itnew = mTracksWork[sector][trkType::UNCONS].size();
442
443 mSideTPC[sector].push_back(mSideTPC[sec][it]);
444 mExtraTPCFwdTime[sector].push_back(mExtraTPCFwdTime[sec][it]);
445 mTracksWork[sector][trkType::UNCONS].emplace_back(pair);
446 mTrackGid[sector][trkType::UNCONS].emplace_back(mTrackGid[sec][trkType::UNCONS][it]);
447
448 if (mMCTruthON) {
449 mTracksLblWork[sector][trkType::UNCONS].emplace_back(mTracksLblWork[sec][trkType::UNCONS][it]);
450 }
451 mLTinfos[sector][trkType::UNCONS].emplace_back(mLTinfos[sec][trkType::UNCONS][it]);
452 mVZtpcOnly[sector].push_back(mVZtpcOnly[sec][it]);
453 mTracksSectIndexCache[trkType::UNCONS][sector].push_back(itnew);
454 }
455 }
456
457 for (auto& it : mTracksSeed[trkType::CONSTR][sec]) {
458 auto& pair = mTracksWork[sec][trkType::CONSTR][it];
459 auto& trc = pair.first;
460 trc.getXYZGlo(globalPos);
461 int sector = o2::math_utils::angle2Sector(TMath::ATan2(globalPos[1], globalPos[0]));
462
463 int itnew = mTracksWork[sector][trkType::CONSTR].size();
464
465 mTracksWork[sector][trkType::CONSTR].emplace_back(pair);
466 mTrackGid[sector][trkType::CONSTR].emplace_back(mTrackGid[sec][trkType::CONSTR][it]);
467
468 if (mMCTruthON) {
469 mTracksLblWork[sector][trkType::CONSTR].emplace_back(mTracksLblWork[sec][trkType::CONSTR][it]);
470 }
471 mLTinfos[sector][trkType::CONSTR].emplace_back(mLTinfos[sec][trkType::CONSTR][it]);
472 mTracksSectIndexCache[trkType::CONSTR][sector].push_back(itnew);
473 }
474 }
475
476 for (int it = 0; it < trkType::SIZE; it++) {
477 for (int sec = o2::constants::math::NSectors - 1; sec > -1; sec--) {
478 mMatchedTracksIndex[sec][it].resize(mTracksWork[sec][it].size());
479 std::fill(mMatchedTracksIndex[sec][it].begin(), mMatchedTracksIndex[sec][it].end(), -1); // initializing all to -1
480 }
481 }
482
483 if (mIsTPCused) {
484 LOG(debug) << "Number of UNCONSTRAINED tracks that failed to be propagated to TOF = " << mNotPropagatedToTOF[trkType::UNCONS];
485
486 // sort tracks in each sector according to their time (increasing in time)
487#ifdef WITH_OPENMP
488#pragma omp parallel for schedule(dynamic) num_threads(mNlanes)
489#endif
490 for (int sec = o2::constants::math::NSectors - 1; sec > -1; sec--) {
491 auto& indexCache = mTracksSectIndexCache[trkType::UNCONS][sec];
492 LOG(debug) << "Sorting sector" << sec << " | " << indexCache.size() << " tracks";
493 if (!indexCache.size()) {
494 continue;
495 }
496 std::sort(indexCache.begin(), indexCache.end(), [this, sec](int a, int b) {
497 auto& trcA = mTracksWork[sec][trkType::UNCONS][a].second;
498 auto& trcB = mTracksWork[sec][trkType::UNCONS][b].second;
499 return ((trcA.getTimeStamp() - trcA.getTimeStampError()) - (trcB.getTimeStamp() - trcB.getTimeStampError()) < 0.);
500 });
501 } // loop over tracks of single sector
502 }
503 if (mIsITSTPCused || mIsTPCTRDused || mIsITSTPCTRDused) {
504 LOG(debug) << "Number of CONSTRAINED tracks that failed to be propagated to TOF = " << mNotPropagatedToTOF[trkType::CONSTR];
505
506 // sort tracks in each sector according to their time (increasing in time)
507#ifdef WITH_OPENMP
508#pragma omp parallel for schedule(dynamic) num_threads(mNlanes)
509#endif
510 for (int sec = o2::constants::math::NSectors - 1; sec > -1; sec--) {
511 auto& indexCache = mTracksSectIndexCache[trkType::CONSTR][sec];
512 LOG(debug) << "Sorting sector" << sec << " | " << indexCache.size() << " tracks";
513 if (!indexCache.size()) {
514 continue;
515 }
516 std::sort(indexCache.begin(), indexCache.end(), [this, sec](int a, int b) {
517 auto& trcA = mTracksWork[sec][trkType::CONSTR][a].second;
518 auto& trcB = mTracksWork[sec][trkType::CONSTR][b].second;
519 return ((trcA.getTimeStamp() - mSigmaTimeCut * trcA.getTimeStampError()) - (trcB.getTimeStamp() - mSigmaTimeCut * trcB.getTimeStampError()) < 0.);
520 });
521 } // loop over tracks of single sector
522 }
523
524 if (mRecoCont->inputsTPCclusters) {
525 mTPCClusterIdxStruct = &mRecoCont->inputsTPCclusters->clusterIndex;
526 mTPCTracksArray = mRecoCont->getTPCTracks();
527 mTPCTrackClusIdx = mRecoCont->getTPCTracksClusterRefs();
528 mTPCRefitterShMap = mRecoCont->clusterShMapTPC;
529 mTPCRefitterOccMap = mRecoCont->occupancyMapTPC;
530 }
531
532 return true;
533}
534
535//______________________________________________
536void MatchTOF::propagateTPCTracks(int sec)
537{
538 auto& trkWork = mTracksWork[sec][trkType::UNCONS];
539
540 for (int it = 0; it < trkWork.size(); it++) {
541 o2::track::TrackParCov& trc = trkWork[it].first;
542
543 o2::track::TrackLTIntegral& intLT0 = mLTinfos[sec][trkType::UNCONS][it];
544
545 const auto& trackTune = TrackTuneParams::Instance();
546 if (!trackTune.sourceLevelTPC) { // correct only if TPC track was not corrected at the source level
547 if (trackTune.useTPCOuterCorr) {
548 trc.updateParams(trackTune.tpcParOuter);
549 }
550 if (trackTune.tpcCovOuterType != TrackTuneParams::AddCovType::Disable) { // only TRD-refitted track have cov.matrix already man>
551 trc.updateCov(mCovDiagOuter, trackTune.tpcCovOuterType == TrackTuneParams::AddCovType::WithCorrelations);
552 }
553 }
554 if (!propagateToRefXWithoutCov(trc, mXRef, 10, mBz)) { // we first propagate to 371 cm without considering the covariance matri
555 mNotPropagatedToTOF[trkType::UNCONS]++;
556 continue;
557 }
558
559 if (trc.getX() < o2::constants::geom::XTPCOuterRef - 1.) {
560 if (!propagateToRefXWithoutCov(trc, o2::constants::geom::XTPCOuterRef, 10, mBz) || std::abs(trc.getZ()) > Geo::MAXHZTOF) { // we check that the propagat>
561 mNotPropagatedToTOF[trkType::UNCONS]++;
562 continue;
563 }
564 }
565
566 o2::base::Propagator::Instance()->estimateLTFast(intLT0, trc);
567
568 // the "rough" propagation worked; now we can propagate considering also the cov matrix
569 if (!propagateToRefX(trc, mXRef, 2, intLT0)) { // || std::abs(trc.getZ()) > Geo::MAXHZTOF) { // we check that the propagation with the cov matrix w>
570 mNotPropagatedToTOF[trkType::UNCONS]++;
571 continue;
572 }
573
574 // printf("time0 %f -> %f (diff = %f, err = %f)\n",trackTime0, time0, trackTime0 - time0, terr);
575 // printf("time errors %f,%f -> %f,%f\n",(_tr.getDeltaTBwd() + 5) * mTPCTBinMUS,(_tr.getDeltaTFwd() + 5) * mTPCTBinMUS,trackTime0 - time0 + terr, ti>
576
577 std::array<float, 3> globalPos;
578 trc.getXYZGlo(globalPos);
579 int sector = o2::math_utils::angle2Sector(TMath::ATan2(globalPos[1], globalPos[0]));
580
581 if (sector == sec) {
582 mTracksSectIndexCache[trkType::UNCONS][sector].push_back(it);
583 } else {
584 mTracksSeed[trkType::UNCONS][sec].push_back(it); // to be moved to another sector
585 }
586 }
587}
588//______________________________________________
589void MatchTOF::propagateConstrTracks(int sec)
590{
591 std::array<float, 3> globalPos;
592
593 auto& trkWork = mTracksWork[sec][trkType::CONSTR];
594
595 for (int it = 0; it < trkWork.size(); it++) {
596 o2::track::TrackParCov& trc = trkWork[it].first;
597 o2::track::TrackLTIntegral& intLT0 = mLTinfos[sec][trkType::CONSTR][it];
598
599 // propagate to matching Xref
600 if (!propagateToRefXWithoutCov(trc, mXRef, 2, mBz)) { // we first propagate to 371 cm without considering the covariance matrix
601 mNotPropagatedToTOF[trkType::CONSTR]++;
602 continue;
603 }
604
605 // the "rough" propagation worked; now we can propagate considering also the cov matrix
606 if (!propagateToRefX(trc, mXRef, 2, intLT0) || std::abs(trc.getZ()) > Geo::MAXHZTOF) { // we check that the propagation with the cov matrix worked;>
607 mNotPropagatedToTOF[trkType::CONSTR]++;
608 continue;
609 }
610
611 trc.getXYZGlo(globalPos);
612
613 int sector = o2::math_utils::angle2Sector(TMath::ATan2(globalPos[1], globalPos[0]));
614
615 if (sector == sec) {
616 mTracksSectIndexCache[trkType::CONSTR][sector].push_back(it);
617 } else {
618 mTracksSeed[trkType::CONSTR][sec].push_back(it); // to be moved to another sector
619 }
620 }
621}
622//______________________________________________
623void MatchTOF::addITSTPCSeed(const o2::dataformats::TrackTPCITS& _tr, o2::dataformats::GlobalTrackID srcGID, float time0, float terr)
624{
625 mIsITSTPCused = true;
626
627 auto trc = _tr.getParamOut();
628 o2::track::TrackLTIntegral intLT0 = _tr.getLTIntegralOut();
629
630 timeEst ts(time0, terr);
631 addConstrainedSeed(trc, srcGID, intLT0, ts);
632}
633//______________________________________________
634void MatchTOF::addConstrainedSeed(o2::track::TrackParCov& trc, o2::dataformats::GlobalTrackID srcGID, o2::track::TrackLTIntegral intLT0, timeEst timeMUS)
635{
636 std::array<float, 3> globalPos;
637 trc.getXYZGlo(globalPos);
638
639 int sector = o2::math_utils::angle2Sector(TMath::ATan2(globalPos[1], globalPos[0]));
640
641 // create working copy of track param
642 mTracksWork[sector][trkType::CONSTR].emplace_back(std::make_pair(trc, timeMUS));
643 mTrackGid[sector][trkType::CONSTR].emplace_back(srcGID);
644 mLTinfos[sector][trkType::CONSTR].emplace_back(intLT0);
645
646 if (mMCTruthON) {
647 mTracksLblWork[sector][trkType::CONSTR].emplace_back(mRecoCont->getTPCITSTrackMCLabel(srcGID));
648 }
649
650 //delete trc; // Check: is this needed?
651} //______________________________________________
652void MatchTOF::addTRDSeed(const o2::trd::TrackTRD& _tr, o2::dataformats::GlobalTrackID srcGID, float time0, float terr)
653{
654 if (srcGID.getSource() == o2::dataformats::GlobalTrackID::TPCTRD) {
655 mIsTPCTRDused = true;
656 } else if (srcGID.getSource() == o2::dataformats::GlobalTrackID::ITSTPCTRD) {
657 mIsITSTPCTRDused = true;
658 } else { // shouldn't happen
659 LOG(error) << "MatchTOF::addTRDSee: srcGID.getSource() = " << srcGID.getSource() << " not allowed; expected ones are: " << o2::dataformats::GlobalTrackID::TPCTRD << " and " << o2::dataformats::GlobalTrackID::ITSTPCTRD;
660 }
661
662 auto trc = _tr.getOuterParam();
663
664 o2::track::TrackLTIntegral intLT0 = _tr.getLTIntegralOut();
665
666 // o2::dataformats::TimeStampWithError<float, float>
667 timeEst ts(time0, terr + mExtraTimeToleranceTRD);
668
669 addConstrainedSeed(trc, srcGID, intLT0, ts);
670}
671//______________________________________________
672void MatchTOF::addTPCSeed(const o2::tpc::TrackTPC& _tr, o2::dataformats::GlobalTrackID srcGID, float time0, float terr)
673{
674 mIsTPCused = true;
675
676 std::array<float, 3> globalPos;
677
678 // create working copy of track param
679 timeEst timeInfo;
680 // set
681 float extraErr = 0;
682 if (mIsCosmics) {
683 extraErr = 100;
684 }
685
686 float trackTime0 = _tr.getTime0() * mTPCTBinMUS;
687 timeInfo.setTimeStampError((_tr.getDeltaTBwd() + 5) * mTPCTBinMUS + extraErr);
688 // timeInfo.setTimeStampError(trackTime0 - time0 + terr + extraErr);
689 timeInfo.setTimeStamp(trackTime0);
690
691 auto trc = _tr.getOuterParam();
692 trc.getXYZGlo(globalPos);
693
694 int sector = o2::math_utils::angle2Sector(TMath::ATan2(globalPos[1], globalPos[0]));
695
696 mSideTPC[sector].push_back(_tr.hasASideClustersOnly() ? 1 : (_tr.hasCSideClustersOnly() ? -1 : 0));
697 mExtraTPCFwdTime[sector].push_back((_tr.getDeltaTFwd() + 5) * mTPCTBinMUS + extraErr);
698 // mExtraTPCFwdTime[sector].push_back(time0 + terr - trackTime0 + extraErr);
699 mTracksWork[sector][trkType::UNCONS].emplace_back(std::make_pair(trc, timeInfo));
700 mTrackGid[sector][trkType::UNCONS].emplace_back(srcGID);
701
702 if (mMCTruthON) {
703 mTracksLblWork[sector][trkType::UNCONS].emplace_back(mRecoCont->getTPCTrackMCLabel(srcGID));
704 }
705 o2::track::TrackLTIntegral intLT0; // mTPCTracksWork.back().getLTIntegralOut(); // we get the integrated length from TPC-ITC outward propagation
706 // compute track length up to now
707 mLTinfos[sector][trkType::UNCONS].emplace_back(intLT0);
708 float vz0 = _tr.getZAt(0, mBz);
709 if (std::abs(vz0) > 9000) {
710 vz0 = _tr.getZ() - _tr.getX() * _tr.getTgl();
711 }
712 mVZtpcOnly[sector].push_back(vz0);
713
714 /*
715 const auto& trackTune = TrackTuneParams::Instance();
716 if (!trackTune.sourceLevelTPC) { // correct only if TPC track was not corrected at the source level
717 if (trackTune.useTPCOuterCorr) {
718 trc.updateParams(trackTune.tpcParOuter);
719 }
720 if (trackTune.tpcCovOuterType != TrackTuneParams::AddCovType::Disable) { // only TRD-refitted track have cov.matrix already manipulated
721 trc.updateCov(mCovDiagOuter, trackTune.tpcCovOuterType == TrackTuneParams::AddCovType::WithCorrelations);
722 }
723 }
724 if (!propagateToRefXWithoutCov(trc, mXRef, 10, mBz)) { // we first propagate to 371 cm without considering the covariance matri
725 mNotPropagatedToTOF[trkType::UNCONS]++;
726 return;
727 }
728
729 if (trc.getX() < o2::constants::geom::XTPCOuterRef - 1.) {
730 if (!propagateToRefX(trc, o2::constants::geom::XTPCOuterRef, 10, intLT0) || std::abs(trc.getZ()) > Geo::MAXHZTOF) { // we check that the propagation with the cov matrix worked; CHECK: can it happ
731 mNotPropagatedToTOF[trkType::UNCONS]++;
732 return;
733 }
734 }
735
736 // the "rough" propagation worked; now we can propagate considering also the cov matrix
737 if (!propagateToRefX(trc, mXRef, 2, intLT0)) { // || std::abs(trc.getZ()) > Geo::MAXHZTOF) { // we check that the propagation with the cov matrix worked; CHECK: can it happen that it does not if the prop>
738 mNotPropagatedToTOF[trkType::UNCONS]++;
739 return;
740 }
741
742 // printf("time0 %f -> %f (diff = %f, err = %f)\n",trackTime0, time0, trackTime0 - time0, terr);
743 // printf("time errors %f,%f -> %f,%f\n",(_tr.getDeltaTBwd() + 5) * mTPCTBinMUS,(_tr.getDeltaTFwd() + 5) * mTPCTBinMUS,trackTime0 - time0 + terr, time0 + terr - trackTime0);
744
745 trc.getXYZGlo(globalPos);
746 int sector = o2::math_utils::angle2Sector(TMath::ATan2(globalPos[1], globalPos[0]));
747
748 mTracksSectIndexCache[trkType::UNCONS][sector].push_back(it);
749 */
750 // delete trc; // Check: is this needed?
751}
752//______________________________________________
753bool MatchTOF::prepareTOFClusters()
754{
755 mTOFClustersArrayInp = mRecoCont->getTOFClusters();
756 mTOFClusLabels = mRecoCont->getTOFClustersMCLabels();
757 mMCTruthON = mTOFClusLabels && mTOFClusLabels->getNElements();
758
760
761 // copy the track params, propagate to reference X and build sector tables
762 mTOFClusWork.clear();
763 // mTOFClusWork.reserve(mNumOfClusters); // we cannot do this, we don't have mNumOfClusters yet
764 // if (mMCTruthON) {
765 // mTOFClusLblWork.clear();
766 // mTOFClusLblWork.reserve(mNumOfClusters);
767 // }
768
769 for (int sec = o2::constants::math::NSectors - 1; sec > -1; sec--) {
770 mTOFClusSectIndexCache[sec].clear();
771 //mTOFClusSectIndexCache[sec].reserve(100 + 1.2 * mNumOfClusters / o2::constants::math::NSectors); // we cannot do this, we don't have mNumOfClusters yet
772 }
773
774 mNumOfClusters = 0;
775
776 int nClusterInCurrentChunk = mTOFClustersArrayInp.size();
777 LOG(debug) << "nClusterInCurrentChunk = " << nClusterInCurrentChunk;
778 mNumOfClusters += nClusterInCurrentChunk;
779 mTOFClusWork.reserve(mTOFClusWork.size() + mNumOfClusters);
780 for (int it = 0; it < nClusterInCurrentChunk; it++) {
781 const Cluster& clOrig = mTOFClustersArrayInp[it];
782 // create working copy of track param
783 mTOFClusWork.emplace_back(clOrig);
784 auto& cl = mTOFClusWork.back();
785 // cache work track index
786 mTOFClusSectIndexCache[cl.getSector()].push_back(mTOFClusWork.size() - 1);
787 }
788
789 // sort clusters in each sector according to their time (increasing in time)
790 for (int sec = o2::constants::math::NSectors - 1; sec > -1; sec--) {
791 auto& indexCache = mTOFClusSectIndexCache[sec];
792 LOG(debug) << "Sorting sector" << sec << " | " << indexCache.size() << " TOF clusters";
793 if (!indexCache.size()) {
794 continue;
795 }
796 std::sort(indexCache.begin(), indexCache.end(), [this](int a, int b) {
797 auto& clA = mTOFClusWork[a];
798 auto& clB = mTOFClusWork[b];
799 return (clA.getTime() - clB.getTime()) < 0.;
800 });
801 } // loop over TOF clusters of single sector
802
803 if (mMatchedClustersIndex) {
804 delete[] mMatchedClustersIndex;
805 }
806 mMatchedClustersIndex = new int[mNumOfClusters];
807 std::fill_n(mMatchedClustersIndex, mNumOfClusters, -1); // initializing all to -1
808
809 return true;
810}
811//______________________________________________
812void MatchTOF::doMatching(int sec)
813{
814 trkType type = trkType::CONSTR;
815
817 auto& cacheTOF = mTOFClusSectIndexCache[sec]; // array of cached TOF cluster indices for this sector; reminder: they are ordered in time!
818 auto& cacheTrk = mTracksSectIndexCache[type][sec]; // array of cached tracks indices for this sector; reminder: they are ordered in time!
819 int nTracks = cacheTrk.size(), nTOFCls = cacheTOF.size();
820 LOG(debug) << "Matching sector " << sec << ": number of tracks: " << nTracks << ", number of TOF clusters: " << nTOFCls;
821 if (!nTracks || !nTOFCls) {
822 return;
823 }
824 int itof0 = 0; // starting index in TOF clusters for matching of the track
825 int detId[2][5]; // at maximum one track can fall in 2 strips during the propagation; the second dimention of the array is the TOF det index
826 float deltaPos[2][3]; // at maximum one track can fall in 2 strips during the propagation; the second dimention of the array is the residuals
827 o2::track::TrackLTIntegral trkLTInt[2]; // Here we store the integrated track length and time for the (max 2) matched strips
828 int nStepsInsideSameStrip[2] = {0, 0}; // number of propagation steps in the same strip (since we have maximum 2 strips, it has dimention = 2)
829 float deltaPosTemp[3];
830 std::array<float, 3> pos;
831 std::array<float, 3> posBeforeProp;
832 float posFloat[3];
833
834 // prematching for TPC only tracks (identify BC candidate to correct z for TPC track accordingly to v_drift)
835
836 LOG(debug) << "Trying to match %d tracks" << cacheTrk.size();
837 for (int itrk = 0; itrk < cacheTrk.size(); itrk++) {
838 for (int ii = 0; ii < 2; ii++) {
839 detId[ii][2] = -1; // before trying to match, we need to inizialize the detId corresponding to the strip number to -1; this is the array that we will use to save the det id of the maximum 2 strips matched
840 nStepsInsideSameStrip[ii] = 0;
841 }
842 int nStripsCrossedInPropagation = 0; // how many strips were hit during the propagation
843 auto& trackWork = mTracksWork[sec][type][cacheTrk[itrk]];
844 auto& trefTrk = trackWork.first;
845 float pt = trefTrk.getPt();
846 auto& intLT = mLTinfos[sec][type][cacheTrk[itrk]];
847 float timeShift = intLT.getL() * 33.35641; // integrated time for 0.75 beta particles in ps, to take into account the t.o.f. delay with respect the interaction BC
848 // using beta=0.75 to cover beta range [0.59 , 1.04] also for a 8 m track lenght with a 10 ns track resolution (TRD)
849
850 // Printf("intLT (before doing anything): length = %f, time (Pion) = %f", intLT.getL(), intLT.getTOF(o2::track::PID::Pion));
851 float minTrkTime = (trackWork.second.getTimeStamp() - mSigmaTimeCut * trackWork.second.getTimeStampError()) * 1.E6 + timeShift; // minimum time in ps
852 float maxTrkTime = (trackWork.second.getTimeStamp() + mSigmaTimeCut * trackWork.second.getTimeStampError()) * 1.E6 + timeShift + 100E3; // maximum time in ps + 100 ns for slow tracks (beta->0.2)
853 const float sqrt12inv = 1. / sqrt(12.);
854 float resT = (trackWork.second.getTimeStampError() + 100E-3) * sqrt12inv;
855 int istep = 1; // number of steps
856 float step = 1.0; // step size in cm
857
858 //uncomment for local debug
859 /*
860 //trefTrk.getXYZGlo(posBeforeProp);
861 //float posBeforeProp[3] = {trefTrk.getX(), trefTrk.getY(), trefTrk.getZ()}; // in local ref system
862 //printf("Global coordinates: posBeforeProp[0] = %f, posBeforeProp[1] = %f, posBeforeProp[2] = %f\n", posBeforeProp[0], posBeforeProp[1], posBeforeProp[2]);
863 //Printf("Radius xy = %f", TMath::Sqrt(posBeforeProp[0]*posBeforeProp[0] + posBeforeProp[1]*posBeforeProp[1]));
864 //Printf("Radius xyz = %f", TMath::Sqrt(posBeforeProp[0]*posBeforeProp[0] + posBeforeProp[1]*posBeforeProp[1] + posBeforeProp[2]*posBeforeProp[2]));
865 */
866
867 // initializing
868 for (int ii = 0; ii < 2; ii++) {
869 for (int iii = 0; iii < 5; iii++) {
870 detId[ii][iii] = -1;
871 }
872 }
873
874 int detIdTemp[5] = {-1, -1, -1, -1, -1}; // TOF detector id at the current propagation point
875
876 double reachedPoint = mXRef + istep * step;
877
878 while (propagateToRefX(trefTrk, reachedPoint, step, intLT) && nStripsCrossedInPropagation <= 2 && reachedPoint < Geo::RMAX) {
879 // while (o2::base::Propagator::Instance()->PropagateToXBxByBz(trefTrk, mXRef + istep * step, MAXSNP, step, 1, &intLT) && nStripsCrossedInPropagation <= 2 && mXRef + istep * step < Geo::RMAX) {
880
881 trefTrk.getXYZGlo(pos);
882 for (int ii = 0; ii < 3; ii++) { // we need to change the type...
883 posFloat[ii] = pos[ii];
884 }
885
886 // uncomment below only for local debug; this will produce A LOT of output - one print per propagation step
887 /*
888 Printf("posFloat[0] = %f, posFloat[1] = %f, posFloat[2] = %f", posFloat[0], posFloat[1], posFloat[2]);
889 Printf("radius xy = %f", TMath::Sqrt(posFloat[0]*posFloat[0] + posFloat[1]*posFloat[1]));
890 Printf("radius xyz = %f", TMath::Sqrt(posFloat[0]*posFloat[0] + posFloat[1]*posFloat[1] + posFloat[2]*posFloat[2]));
891 */
892
893 for (int idet = 0; idet < 5; idet++) {
894 detIdTemp[idet] = -1;
895 }
896
897 Geo::getPadDxDyDz(posFloat, detIdTemp, deltaPosTemp, sec);
898
899 reachedPoint += step;
900
901 if (detIdTemp[2] == -1) {
902 continue;
903 }
904
905 // uncomment below only for local debug; this will produce A LOT of output - one print per propagation step
906 //Printf("detIdTemp[0] = %d, detIdTemp[1] = %d, detIdTemp[2] = %d, detIdTemp[3] = %d, detIdTemp[4] = %d", detIdTemp[0], detIdTemp[1], detIdTemp[2], detIdTemp[3], detIdTemp[4]);
907 // if (nStripsCrossedInPropagation == 0) { // print in case you have a useful propagation
908 // LOG(debug) << "*********** We have crossed a strip during propagation!*********";
909 // LOG(debug) << "Global coordinates: pos[0] = " << pos[0] << ", pos[1] = " << pos[1] << ", pos[2] = " << pos[2];
910 // LOG(debug) << "detIdTemp[0] = " << detIdTemp[0] << ", detIdTemp[1] = " << detIdTemp[1] << ", detIdTemp[2] = " << detIdTemp[2] << ", detIdTemp[3] = " << detIdTemp[3] << ", detIdTemp[4] = " << detIdTemp[4];
911 // LOG(debug) << "deltaPosTemp[0] = " << deltaPosTemp[0] << ", deltaPosTemp[1] = " << deltaPosTemp[1] << " deltaPosTemp[2] = " << deltaPosTemp[2];
912 // } else {
913 // LOG(debug) << "*********** We have NOT crossed a strip during propagation!*********";
914 // LOG(debug) << "Global coordinates: pos[0] = " << pos[0] << ", pos[1] = " << pos[1] << ", pos[2] = " << pos[2];
915 // LOG(debug) << "detIdTemp[0] = " << detIdTemp[0] << ", detIdTemp[1] = " << detIdTemp[1] << ", detIdTemp[2] = " << detIdTemp[2] << ", detIdTemp[3] = " << detIdTemp[3] << ", detIdTemp[4] = " << detIdTemp[4];
916 // LOG(debug) << "deltaPosTemp[0] = " << deltaPosTemp[0] << ", deltaPosTemp[1] = " << deltaPosTemp[1] << " deltaPosTemp[2] = " << deltaPosTemp[2];
917 // }
918
919 // check if after the propagation we are in a TOF strip
920 // we ended in a TOF strip
921 // LOG(debug) << "nStripsCrossedInPropagation = " << nStripsCrossedInPropagation << ", detId[nStripsCrossedInPropagation][0] = " << detId[nStripsCrossedInPropagation][0] << ", detIdTemp[0] = " << detIdTemp[0] << ", detId[nStripsCrossedInPropagation][1] = " << detId[nStripsCrossedInPropagation][1] << ", detIdTemp[1] = " << detIdTemp[1] << ", detId[nStripsCrossedInPropagation][2] = " << detId[nStripsCrossedInPropagation][2] << ", detIdTemp[2] = " << detIdTemp[2];
922
923 if (nStripsCrossedInPropagation == 0 || // we are crossing a strip for the first time...
924 (nStripsCrossedInPropagation >= 1 && (detId[nStripsCrossedInPropagation - 1][0] != detIdTemp[0] || detId[nStripsCrossedInPropagation - 1][1] != detIdTemp[1] || detId[nStripsCrossedInPropagation - 1][2] != detIdTemp[2]))) { // ...or we are crossing a new strip
925 if (nStripsCrossedInPropagation == 0) {
926 LOG(debug) << "We cross a strip for the first time";
927 }
928 if (nStripsCrossedInPropagation == 2) {
929 break; // we have already matched 2 strips, we cannot match more
930 }
931 nStripsCrossedInPropagation++;
932 }
933 //Printf("nStepsInsideSameStrip[nStripsCrossedInPropagation-1] = %d", nStepsInsideSameStrip[nStripsCrossedInPropagation - 1]);
934 if (nStepsInsideSameStrip[nStripsCrossedInPropagation - 1] == 0) {
935 // fine propagation inside the strip -> 1 mm step
936 trkLTInt[nStripsCrossedInPropagation - 1] = intLT;
937 // temporary variables since propagation can fail
938 int detIdTemp2[5] = {0, 0, 0, 0, 0};
939 float deltaPosTemp2[3] = {deltaPosTemp[0], deltaPosTemp[1], deltaPosTemp[2]};
940 int nstep = 0;
941 const int maxnstep = 50;
942 float xStart = trefTrk.getX();
943 float xStop = xStart;
944 trefTrk.getXYZGlo(pos);
945 for (int ii = 0; ii < 3; ii++) { // we need to change the type...
946 posFloat[ii] = pos[ii];
947 }
948 while (deltaPosTemp2[1] < -0.05 && detIdTemp2[2] != -1 && nstep < maxnstep) { // continuing propagation if dy is negative and we are still inside the strip volume
949 nstep++;
950 xStop += 0.1;
951 propagateToRefXWithoutCov(trefTrk, xStop, 0.1, mBz, posFloat);
952
953 Geo::getPadDxDyDz(posFloat, detIdTemp2, deltaPosTemp2, sec);
954 if (detIdTemp2[2] != -1) { // if propation was succesful -> update params
955 float dx = deltaPosTemp2[0] - deltaPosTemp[0];
956 float dy = deltaPosTemp2[1] - deltaPosTemp[1];
957 float dz = deltaPosTemp2[2] - deltaPosTemp[2];
958 updateTL(trkLTInt[nStripsCrossedInPropagation - 1], sqrt(dx * dx + dy * dy + dz * dz));
959 detIdTemp[0] = detIdTemp2[0];
960 detIdTemp[1] = detIdTemp2[1];
961 detIdTemp[2] = detIdTemp2[2];
962 detIdTemp[3] = detIdTemp2[3];
963 detIdTemp[4] = detIdTemp2[4];
964 deltaPosTemp[0] = deltaPosTemp2[0];
965 deltaPosTemp[1] = deltaPosTemp2[1];
966 deltaPosTemp[2] = deltaPosTemp2[2];
967 }
968 }
969
970 // adjust accordingly to DeltaY
971 updateTL(trkLTInt[nStripsCrossedInPropagation - 1], -deltaPosTemp[1]);
972
973 detId[nStripsCrossedInPropagation - 1][0] = detIdTemp[0];
974 detId[nStripsCrossedInPropagation - 1][1] = detIdTemp[1];
975 detId[nStripsCrossedInPropagation - 1][2] = detIdTemp[2];
976 detId[nStripsCrossedInPropagation - 1][3] = detIdTemp[3];
977 detId[nStripsCrossedInPropagation - 1][4] = detIdTemp[4];
978 deltaPos[nStripsCrossedInPropagation - 1][0] = deltaPosTemp[0];
979 deltaPos[nStripsCrossedInPropagation - 1][1] = deltaPosTemp[1];
980 deltaPos[nStripsCrossedInPropagation - 1][2] = deltaPosTemp[2];
981 // Printf("intLT (after matching to strip %d): length = %f, time (Pion) = %f", nStripsCrossedInPropagation - 1, trkLTInt[nStripsCrossedInPropagation - 1].getL(), trkLTInt[nStripsCrossedInPropagation - 1].getTOF(o2::track::PID::Pion));
982 nStepsInsideSameStrip[nStripsCrossedInPropagation - 1]++;
983 }
984 /* // obsolete
985 else { // a further propagation step in the same strip -> update info (we sum up on all matching with strip - we will divide for the number of steps a bit below)
986 // N.B. the integrated length and time are taken (at least for now) from the first time we crossed the strip, so here we do nothing with those
987 deltaPos[nStripsCrossedInPropagation - 1][0] += deltaPosTemp[0] + (detIdTemp[4] - detId[nStripsCrossedInPropagation - 1][4]) * Geo::XPAD; // residual in x
988 deltaPos[nStripsCrossedInPropagation - 1][1] += deltaPosTemp[1]; // residual in y
989 deltaPos[nStripsCrossedInPropagation - 1][2] += deltaPosTemp[2] + (detIdTemp[3] - detId[nStripsCrossedInPropagation - 1][3]) * Geo::ZPAD; // residual in z
990 nStepsInsideSameStrip[nStripsCrossedInPropagation - 1]++;
991 }
992 */
993 }
994
995 for (Int_t imatch = 0; imatch < nStripsCrossedInPropagation; imatch++) {
996 // we take as residual the average of the residuals along the propagation in the same strip
997 deltaPos[imatch][0] /= nStepsInsideSameStrip[imatch];
998 deltaPos[imatch][1] /= nStepsInsideSameStrip[imatch];
999 deltaPos[imatch][2] /= nStepsInsideSameStrip[imatch];
1000 }
1001
1002 if (nStripsCrossedInPropagation == 0) {
1003 continue; // the track never hit a TOF strip during the propagation
1004 }
1005 bool foundCluster = false;
1006 for (auto itof = itof0; itof < nTOFCls; itof++) {
1007 // printf("itof = %d\n", itof);
1008 auto& trefTOF = mTOFClusWork[cacheTOF[itof]];
1009 // compare the times of the track and the TOF clusters - remember that they both are ordered in time!
1010 //Printf("trefTOF.getTime() = %f, maxTrkTime = %f, minTrkTime = %f", trefTOF.getTime(), maxTrkTime, minTrkTime);
1011
1012 if (trefTOF.getTime() < minTrkTime) { // this cluster has a time that is too small for the current track, we will get to the next one
1013 //Printf("In trefTOF.getTime() < minTrkTime");
1014 itof0 = itof + 1; // but for the next track that we will check, we will ignore this cluster (the time is anyway too small)
1015 continue;
1016 }
1017 if (trefTOF.getTime() > maxTrkTime) { // no more TOF clusters can be matched to this track
1018 break;
1019 }
1020
1021 int mainChannel = trefTOF.getMainContributingChannel();
1022 int indices[5];
1023 Geo::getVolumeIndices(mainChannel, indices);
1024
1025 // compute fine correction using cluster position instead of pad center
1026 // this because in case of multiple-hit cluster position is averaged on all pads contributing to the cluster (then error position matrix can be used for Chi2 if nedeed)
1027 int ndigits = 1;
1028 float posCorr[3] = {0, 0, 0};
1029
1030 if (trefTOF.isBitSet(Cluster::kLeft)) {
1031 posCorr[0] += Geo::XPAD, ndigits++;
1032 }
1033 if (trefTOF.isBitSet(Cluster::kUpLeft)) {
1034 posCorr[0] += Geo::XPAD, posCorr[2] -= Geo::ZPAD, ndigits++;
1035 }
1036 if (trefTOF.isBitSet(Cluster::kDownLeft)) {
1037 posCorr[0] += Geo::XPAD, posCorr[2] += Geo::ZPAD, ndigits++;
1038 }
1039 if (trefTOF.isBitSet(Cluster::kUp)) {
1040 posCorr[2] -= Geo::ZPAD, ndigits++;
1041 }
1042 if (trefTOF.isBitSet(Cluster::kDown)) {
1043 posCorr[2] += Geo::ZPAD, ndigits++;
1044 }
1045 if (trefTOF.isBitSet(Cluster::kRight)) {
1046 posCorr[0] -= Geo::XPAD, ndigits++;
1047 }
1048 if (trefTOF.isBitSet(Cluster::kUpRight)) {
1049 posCorr[0] -= Geo::XPAD, posCorr[2] -= Geo::ZPAD, ndigits++;
1050 }
1051 if (trefTOF.isBitSet(Cluster::kDownRight)) {
1052 posCorr[0] -= Geo::XPAD, posCorr[2] += Geo::ZPAD, ndigits++;
1053 }
1054
1055 float ndifInv = 1. / ndigits;
1056 if (ndigits > 1) {
1057 posCorr[0] *= ndifInv;
1058 posCorr[1] *= ndifInv;
1059 posCorr[2] *= ndifInv;
1060 }
1061
1062 int trackIdTOF;
1063 int eventIdTOF;
1064 int sourceIdTOF;
1065 for (auto iPropagation = 0; iPropagation < nStripsCrossedInPropagation; iPropagation++) {
1066 int bct0 = int((mTOFClusWork[cacheTOF[itof]].getTime() - trkLTInt[iPropagation].getTOF(0) + 5000) * Geo::BC_TIME_INPS_INV); // bc taken assuming speed of light (el) and 5 ns of margin
1067 if (bct0 < 0) { // if negative time (it can happen at the beginng of the TF int was truncated per excess... adjusting)
1068 bct0--;
1069 }
1070 float tof = mTOFClusWork[cacheTOF[itof]].getTime() - bct0 * Geo::BC_TIME_INPS;
1071 if (tof - trkLTInt[iPropagation].getTOF(6) > 2000) { // reject tracks slower than triton
1072 continue;
1073 }
1074 float cosangle = TMath::Cos(Geo::getAngles(indices[1], indices[2]) * TMath::DegToRad());
1075 const float errXinvMin = 1. / (mMatchParams->maxResX * mMatchParams->maxResX);
1076 const float errZinvMin = 1. / (mMatchParams->maxResZ * mMatchParams->maxResZ);
1077 float errXinv2 = 1. / (trefTrk.getSigmaY2());
1078 float errZinv2 = 1. / (trefTrk.getSigmaZ2() * cosangle); // should be valid only at eta=0
1079
1080 if (errXinv2 < errXinvMin) {
1081 errXinv2 = errXinvMin;
1082 }
1083 if (errZinv2 < errZinvMin) {
1084 errZinv2 = errZinvMin;
1085 }
1086
1087 LOG(debug) << "TOF Cluster [" << itof << ", " << cacheTOF[itof] << "]: indices = " << indices[0] << ", " << indices[1] << ", " << indices[2] << ", " << indices[3] << ", " << indices[4];
1088 LOG(debug) << "Propagated Track [" << itrk << "]: detId[" << iPropagation << "] = " << detId[iPropagation][0] << ", " << detId[iPropagation][1] << ", " << detId[iPropagation][2] << ", " << detId[iPropagation][3] << ", " << detId[iPropagation][4];
1089 float resX = deltaPos[iPropagation][0] - (indices[4] - detId[iPropagation][4]) * Geo::XPAD + posCorr[0]; // readjusting the residuals due to the fact that the propagation fell in a pad that was not exactly the one of the cluster
1090 float resZ = deltaPos[iPropagation][2] - (indices[3] - detId[iPropagation][3]) * Geo::ZPAD + posCorr[2]; // readjusting the residuals due to the fact that the propagation fell in a pad that was not exactly the one of the cluster
1091 float resY = deltaPos[iPropagation][1];
1092 float resXor = resX;
1093 float resZor = resZ;
1094 float res = TMath::Sqrt(resX * resX + resZ * resZ);
1095
1096 if (resX < -1.25) { // distance from closest border
1097 resX += 1.25;
1098 } else if (resX > 1.25) {
1099 resX -= 1.25;
1100 } else {
1101 resX = 1E-3 / (pt + 1E-3); // high-pt should be favoured
1102 }
1103
1104 if (resZ < -1.75) { // distance from closest border
1105 resZ += 1.75;
1106 } else if (resZ > 1.75) {
1107 resZ -= 1.75;
1108 } else {
1109 resZ = 1E-3 / (pt + 1E-3); // high-pt should be favoured
1110 }
1111
1112 LOG(debug) << "resX = " << resX << ", resZ = " << resZ << ", res = " << res;
1113 if (indices[0] != detId[iPropagation][0]) {
1114 continue;
1115 }
1116 if (indices[1] != detId[iPropagation][1]) {
1117 continue;
1118 }
1119 if (indices[2] != detId[iPropagation][2]) {
1120 continue;
1121 }
1122 float chi2 = 0.5 * (resX * resX * errXinv2 + resZ * resZ * errZinv2); // TODO: take into account also the time!
1123
1124 if (res < mSpaceTolerance && chi2 < mMatchParams->maxChi2) { // matching ok!
1125 LOG(debug) << "MATCHING FOUND: We have a match! between track " << mTracksSectIndexCache[type][sec][itrk] << " and TOF cluster " << mTOFClusSectIndexCache[indices[0]][itof];
1126 foundCluster = true;
1127 // set event indexes (to be checked)
1128 int eventIndexTOFCluster = mTOFClusSectIndexCache[indices[0]][itof];
1129 mMatchedTracksPairsSec[sec].emplace_back(cacheTrk[itrk], eventIndexTOFCluster, mTOFClusWork[cacheTOF[itof]].getTime(), chi2, trkLTInt[iPropagation], mTrackGid[sec][type][cacheTrk[itrk]], type, (trefTOF.getTime() - (minTrkTime + maxTrkTime - 100E3) * 0.5) * 1E-6, trefTOF.getZ(), resXor, resZor, resY); // subracting 100 ns to max track which was artificially added
1130 mMatchedTracksPairsSec[sec][mMatchedTracksPairsSec[sec].size() - 1].setPt(pt);
1131 mMatchedTracksPairsSec[sec][mMatchedTracksPairsSec[sec].size() - 1].setResX(sqrt(1. / errXinv2));
1132 mMatchedTracksPairsSec[sec][mMatchedTracksPairsSec[sec].size() - 1].setResZ(sqrt(1. / errZinv2));
1133 mMatchedTracksPairsSec[sec][mMatchedTracksPairsSec[sec].size() - 1].setResT(resT);
1134 mMatchedTracksPairsSec[sec][mMatchedTracksPairsSec[sec].size() - 1].setVz(0.0); // not needed for constrained tracks
1135 mMatchedTracksPairsSec[sec][mMatchedTracksPairsSec[sec].size() - 1].setChannel(mainChannel);
1136 }
1137 }
1138 }
1139 }
1140 return;
1141}
1142//______________________________________________
1143void MatchTOF::doMatchingForTPC(int sec)
1144{
1145 float vdriftInBC = Geo::BC_TIME_INPS * 1E-6 * mTPCVDrift;
1146
1147 int bc_grouping = 40;
1148 int bc_grouping_tolerance = bc_grouping + mTimeTolerance / 25;
1149 int bc_grouping_half = (bc_grouping + 1) / 2;
1150 double BCgranularity = Geo::BC_TIME_INPS * bc_grouping;
1151
1153 auto& cacheTOF = mTOFClusSectIndexCache[sec]; // array of cached TOF cluster indices for this sector; reminder: they are ordered in time!
1154 auto& cacheTrk = mTracksSectIndexCache[trkType::UNCONS][sec]; // array of cached tracks indices for this sector; reminder: they are ordered in time!
1155 int nTracks = cacheTrk.size(), nTOFCls = cacheTOF.size();
1156 LOG(debug) << "Matching sector " << sec << ": number of tracks: " << nTracks << ", number of TOF clusters: " << nTOFCls;
1157 if (!nTracks || !nTOFCls) {
1158 return;
1159 }
1160 int itof0 = 0; // starting index in TOF clusters for matching of the track
1161 float deltaPosTemp[3];
1162 std::array<float, 3> pos;
1163 std::array<float, 3> posBeforeProp;
1164 float posFloat[3];
1165
1166 // prematching for TPC only tracks (identify BC candidate to correct z for TPC track accordingly to v_drift)
1167
1168 std::vector<unsigned long> BCcand;
1169
1170 std::vector<int> nStripsCrossedInPropagation;
1171 std::vector<std::array<std::array<int, 5>, 2>> detId;
1172 std::vector<std::array<o2::track::TrackLTIntegral, 2>> trkLTInt;
1173 std::vector<std::array<std::array<float, 3>, 2>> deltaPos;
1174 std::vector<std::array<int, 2>> nStepsInsideSameStrip;
1175 std::vector<std::array<float, 2>> Zshift;
1176
1177 LOG(debug) << "Trying to match %d tracks" << cacheTrk.size();
1178
1179 for (int itrk = 0; itrk < cacheTrk.size(); itrk++) {
1180 auto& trackWork = mTracksWork[sec][trkType::UNCONS][cacheTrk[itrk]];
1181 auto& trefTrk = trackWork.first;
1182 float pt = trefTrk.getPt();
1183 auto& intLT = mLTinfos[sec][trkType::UNCONS][cacheTrk[itrk]];
1184
1185 float timeShift = intLT.getL() * 33.35641; // integrated time for 0.75 beta particles in ps, to take into account the t.o.f. delay with respect the interaction BC
1186 // using beta=0.75 to cover beta range [0.59 , 1.04] also for a 8 m track lenght with a 10 ns track resolution (TRD)
1187
1188 BCcand.clear();
1189 nStripsCrossedInPropagation.clear();
1190
1191 int side = mSideTPC[sec][cacheTrk[itrk]];
1192 // look at BC candidates for the track
1193 double tpctime = trackWork.second.getTimeStamp(); // in mus
1194 double minTrkTime = (tpctime - trackWork.second.getTimeStampError()) * 1.E6 + timeShift; // minimum time in ps
1195 minTrkTime = int(minTrkTime / BCgranularity) * BCgranularity; // align min to a BC
1196 double maxTrkTime = (tpctime + mExtraTPCFwdTime[sec][cacheTrk[itrk]]) * 1.E6 + timeShift; // maximum time in ps
1197 const float sqrt12inv = 1. / sqrt(12.);
1198 float resT = (maxTrkTime - minTrkTime) * sqrt12inv;
1199
1200 if (mIsCosmics) {
1201 for (double tBC = minTrkTime; tBC < maxTrkTime; tBC += BCgranularity) {
1202 unsigned long ibc = (unsigned long)(tBC * Geo::BC_TIME_INPS_INV);
1203 BCcand.emplace_back(ibc);
1204 nStripsCrossedInPropagation.emplace_back(0);
1205 }
1206 }
1207
1208 int itofMax = nTOFCls;
1209
1210 for (auto itof = itof0; itof < nTOFCls; itof++) {
1211 auto& trefTOF = mTOFClusWork[cacheTOF[itof]];
1212
1213 if (trefTOF.getTime() < minTrkTime) { // this cluster has a time that is too small for the current track, we will get to the next one
1214 itof0 = itof + 1;
1215 continue;
1216 }
1217
1218 if (trefTOF.getTime() > maxTrkTime) { // this cluster has a time that is too large for the current track, close loop
1219 itofMax = itof;
1220 break;
1221 }
1222
1223 if ((trefTOF.getZ() * side < 0) && ((side > 0) != (trackWork.first.getTgl() > 0))) {
1224 continue;
1225 }
1226
1227 unsigned long bc = (unsigned long)(trefTOF.getTime() * Geo::BC_TIME_INPS_INV);
1228
1229 bc = (bc / bc_grouping_half) * bc_grouping_half;
1230
1231 bool isalreadyin = false;
1232
1233 for (int k = 0; k < BCcand.size(); k++) {
1234 if (bc == BCcand[k]) {
1235 isalreadyin = true;
1236 }
1237 }
1238
1239 if (!isalreadyin) {
1240 BCcand.emplace_back(bc);
1241 nStripsCrossedInPropagation.emplace_back(0);
1242 }
1243 }
1244
1245 detId.clear();
1246 detId.reserve(BCcand.size());
1247 trkLTInt.clear();
1248 trkLTInt.reserve(BCcand.size());
1249 Zshift.clear();
1250 Zshift.reserve(BCcand.size());
1251 deltaPos.clear();
1252 deltaPos.reserve(BCcand.size());
1253 nStepsInsideSameStrip.clear();
1254 nStepsInsideSameStrip.reserve(BCcand.size());
1255
1256 // Printf("intLT (before doing anything): length = %f, time (Pion) = %f", intLT.getL(), intLT.getTOF(o2::track::PID::Pion));
1257 int istep = 1; // number of steps
1258 float step = 1.0; // step size in cm
1259
1260 //uncomment for local debug
1261 /*
1262 //trefTrk.getXYZGlo(posBeforeProp);
1263 //float posBeforeProp[3] = {trefTrk.getX(), trefTrk.getY(), trefTrk.getZ()}; // in local ref system
1264 //printf("Global coordinates: posBeforeProp[0] = %f, posBeforeProp[1] = %f, posBeforeProp[2] = %f\n", posBeforeProp[0], posBeforeProp[1], posBeforeProp[2]);
1265 //Printf("Radius xy = %f", TMath::Sqrt(posBeforeProp[0]*posBeforeProp[0] + posBeforeProp[1]*posBeforeProp[1]));
1266 //Printf("Radius xyz = %f", TMath::Sqrt(posBeforeProp[0]*posBeforeProp[0] + posBeforeProp[1]*posBeforeProp[1] + posBeforeProp[2]*posBeforeProp[2]));
1267 */
1268
1269 int detIdTemp[5] = {-1, -1, -1, -1, -1}; // TOF detector id at the current propagation point
1270
1271 double reachedPoint = mXRef + istep * step;
1272
1273 // initializing
1274 for (int ibc = 0; ibc < BCcand.size(); ibc++) {
1275 for (int ii = 0; ii < 2; ii++) {
1276 nStepsInsideSameStrip[ibc][ii] = 0;
1277 for (int iii = 0; iii < 5; iii++) {
1278 detId[ibc][ii][iii] = -1;
1279 }
1280 }
1281 }
1282 while (propagateToRefX(trefTrk, reachedPoint, step, intLT) && reachedPoint < Geo::RMAX) {
1283 // while (o2::base::Propagator::Instance()->PropagateToXBxByBz(trefTrk, mXRef + istep * step, MAXSNP, step, 1, &intLT) && nStripsCrossedInPropagation <= 2 && mXRef + istep * step < Geo::RMAX) {
1284
1285 trefTrk.getXYZGlo(pos);
1286 for (int ii = 0; ii < 3; ii++) { // we need to change the type...
1287 posFloat[ii] = pos[ii];
1288 }
1289
1290 // uncomment below only for local debug; this will produce A LOT of output - one print per propagation step
1291 /*
1292 Printf("posFloat[0] = %f, posFloat[1] = %f, posFloat[2] = %f", posFloat[0], posFloat[1], posFloat[2]);
1293 Printf("radius xy = %f", TMath::Sqrt(posFloat[0]*posFloat[0] + posFloat[1]*posFloat[1]));
1294 Printf("radius xyz = %f", TMath::Sqrt(posFloat[0]*posFloat[0] + posFloat[1]*posFloat[1] + posFloat[2]*posFloat[2]));
1295 */
1296
1297 reachedPoint += step;
1298
1299 // check if you fall in a strip
1300 for (int ibc = 0; ibc < BCcand.size(); ibc++) {
1301 for (int idet = 0; idet < 5; idet++) {
1302 detIdTemp[idet] = -1;
1303 }
1304
1305 if (side > 0) {
1306 posFloat[2] = pos[2] - mTPCVDrift * (trackWork.second.getTimeStamp() - BCcand[ibc] * Geo::BC_TIME_INPS * 1E-6);
1307 } else if (side < 0) {
1308 posFloat[2] = pos[2] + mTPCVDrift * (trackWork.second.getTimeStamp() - BCcand[ibc] * Geo::BC_TIME_INPS * 1E-6);
1309 } else {
1310 posFloat[2] = pos[2];
1311 }
1312
1313 float ZshiftCurrent = posFloat[2] - pos[2];
1314
1315 Geo::getPadDxDyDz(posFloat, detIdTemp, deltaPosTemp, sec);
1316
1317 if (detIdTemp[2] == -1) {
1318 continue;
1319 }
1320
1321 if (nStripsCrossedInPropagation[ibc] == 0 || // we are crossing a strip for the first time...
1322 (nStripsCrossedInPropagation[ibc] >= 1 && (detId[ibc][nStripsCrossedInPropagation[ibc] - 1][0] != detIdTemp[0] || detId[ibc][nStripsCrossedInPropagation[ibc] - 1][1] != detIdTemp[1] || detId[ibc][nStripsCrossedInPropagation[ibc] - 1][2] != detIdTemp[2]))) { // ...or we are crossing a new strip
1323 if (nStripsCrossedInPropagation[ibc] == 0) {
1324 LOG(debug) << "We cross a strip for the first time";
1325 }
1326 if (nStripsCrossedInPropagation[ibc] == 2) {
1327 continue; // we have already matched 2 strips, we cannot match more
1328 }
1329 nStripsCrossedInPropagation[ibc]++;
1330 }
1331
1332 //Printf("nStepsInsideSameStrip[nStripsCrossedInPropagation-1] = %d", nStepsInsideSameStrip[nStripsCrossedInPropagation - 1]);
1333 if (nStepsInsideSameStrip[ibc][nStripsCrossedInPropagation[ibc] - 1] == 0) {
1334 trkLTInt[ibc][nStripsCrossedInPropagation[ibc] - 1] = intLT;
1335 // temporary variables since propagation can fail
1336 int detIdTemp2[5] = {0, 0, 0, 0, 0};
1337 float deltaPosTemp2[3] = {deltaPosTemp[0], deltaPosTemp[1], deltaPosTemp[2]};
1338 int nstep = 0;
1339 const int maxnstep = 50;
1340 float xStart = trefTrk.getX();
1341 float xStop = xStart;
1342 trefTrk.getXYZGlo(pos);
1343 for (int ii = 0; ii < 3; ii++) { // we need to change the type...
1344 posFloat[ii] = pos[ii];
1345 }
1346
1347 while (deltaPosTemp2[1] < -0.05 && detIdTemp2[2] != -1 && nstep < maxnstep) { // continuing propagation if dy is negative and we are still inside the strip volume
1348 nstep++;
1349 xStop += 0.1;
1350 propagateToRefXWithoutCov(trefTrk, xStop, 0.1, mBz, posFloat);
1351
1352 posFloat[2] += ZshiftCurrent;
1353
1354 Geo::getPadDxDyDz(posFloat, detIdTemp2, deltaPosTemp2, sec);
1355 if (detIdTemp2[2] != -1) { // if propation was succesful -> update params
1356 float dx = deltaPosTemp2[0] - deltaPosTemp[0];
1357 float dy = deltaPosTemp2[1] - deltaPosTemp[1];
1358 float dz = deltaPosTemp2[2] - deltaPosTemp[2];
1359 updateTL(trkLTInt[ibc][nStripsCrossedInPropagation[ibc] - 1], sqrt(dx * dx + dy * dy + dz * dz));
1360 detIdTemp[0] = detIdTemp2[0];
1361 detIdTemp[1] = detIdTemp2[1];
1362 detIdTemp[2] = detIdTemp2[2];
1363 detIdTemp[3] = detIdTemp2[3];
1364 detIdTemp[4] = detIdTemp2[4];
1365 deltaPosTemp[0] = deltaPosTemp2[0];
1366 deltaPosTemp[1] = deltaPosTemp2[1];
1367 deltaPosTemp[2] = deltaPosTemp2[2];
1368 }
1369 }
1370
1371 // adjust accordingly to DeltaY
1372 updateTL(trkLTInt[ibc][nStripsCrossedInPropagation[ibc] - 1], -deltaPosTemp[1]);
1373
1374 detId[ibc][nStripsCrossedInPropagation[ibc] - 1][0] = detIdTemp[0];
1375 detId[ibc][nStripsCrossedInPropagation[ibc] - 1][1] = detIdTemp[1];
1376 detId[ibc][nStripsCrossedInPropagation[ibc] - 1][2] = detIdTemp[2];
1377 detId[ibc][nStripsCrossedInPropagation[ibc] - 1][3] = detIdTemp[3];
1378 detId[ibc][nStripsCrossedInPropagation[ibc] - 1][4] = detIdTemp[4];
1379 deltaPos[ibc][nStripsCrossedInPropagation[ibc] - 1][0] = deltaPosTemp[0];
1380 deltaPos[ibc][nStripsCrossedInPropagation[ibc] - 1][1] = deltaPosTemp[1];
1381 deltaPos[ibc][nStripsCrossedInPropagation[ibc] - 1][2] = deltaPosTemp[2];
1382
1383 Zshift[ibc][nStripsCrossedInPropagation[ibc] - 1] = ZshiftCurrent;
1384 // Printf("intLT (after matching to strip %d): length = %f, time (Pion) = %f", nStripsCrossedInPropagation - 1, trkLTInt[nStripsCrossedInPropagation - 1].getL(), trkLTInt[nStripsCrossedInPropagation - 1].getTOF(o2::track::PID::Pion));
1385 nStepsInsideSameStrip[ibc][nStripsCrossedInPropagation[ibc] - 1]++;
1386 }
1387 /* // obsolete
1388 else { // a further propagation step in the same strip -> update info (we sum up on all matching with strip - we will divide for the number of steps a bit below)
1389 // N.B. the integrated length and time are taken (at least for now) from the first time we crossed the strip, so here we do nothing with those
1390 deltaPos[ibc][nStripsCrossedInPropagation[ibc] - 1][0] += deltaPosTemp[0] + (detIdTemp[4] - detId[ibc][nStripsCrossedInPropagation[ibc] - 1][4]) * Geo::XPAD; // residual in x
1391 deltaPos[ibc][nStripsCrossedInPropagation[ibc] - 1][1] += deltaPosTemp[1]; // residual in y
1392 deltaPos[ibc][nStripsCrossedInPropagation[ibc] - 1][2] += deltaPosTemp[2] + (detIdTemp[3] - detId[ibc][nStripsCrossedInPropagation[ibc] - 1][3]) * Geo::ZPAD; // residual in z
1393 nStepsInsideSameStrip[ibc][nStripsCrossedInPropagation[ibc] - 1]++;
1394 }
1395 */
1396 }
1397 }
1398 for (int ibc = 0; ibc < BCcand.size(); ibc++) {
1399 float minTime = (BCcand[ibc] - bc_grouping_tolerance) * Geo::BC_TIME_INPS;
1400 float maxTime = (BCcand[ibc] + bc_grouping_tolerance) * Geo::BC_TIME_INPS;
1401 for (Int_t imatch = 0; imatch < nStripsCrossedInPropagation[ibc]; imatch++) {
1402 // we take as residual the average of the residuals along the propagation in the same strip
1403 deltaPos[ibc][imatch][0] /= nStepsInsideSameStrip[ibc][imatch];
1404 deltaPos[ibc][imatch][1] /= nStepsInsideSameStrip[ibc][imatch];
1405 deltaPos[ibc][imatch][2] /= nStepsInsideSameStrip[ibc][imatch];
1406 }
1407
1408 if (nStripsCrossedInPropagation[ibc] == 0) {
1409 continue; // the track never hit a TOF strip during the propagation
1410 }
1411
1412 bool foundCluster = false;
1413 for (auto itof = itof0; itof < itofMax; itof++) {
1414 // printf("itof = %d\n", itof);
1415 auto& trefTOF = mTOFClusWork[cacheTOF[itof]];
1416 // compare the times of the track and the TOF clusters - remember that they both are ordered in time!
1417
1418 if (trefTOF.getTime() < minTime) { // this cluster has a time that is too small for the current track, we will get to the next one
1419 continue;
1420 }
1421 if (trefTOF.getTime() > maxTime) { // no more TOF clusters can be matched to this track
1422 break;
1423 }
1424
1425 int mainChannel = trefTOF.getMainContributingChannel();
1426 int indices[5];
1427 Geo::getVolumeIndices(mainChannel, indices);
1428
1429 bool isInStrip = false;
1430 for (auto iPropagation = 0; iPropagation < nStripsCrossedInPropagation[ibc]; iPropagation++) {
1431 if (detId[ibc][iPropagation][1] == indices[1] && detId[ibc][iPropagation][2] == indices[2]) {
1432 isInStrip = true;
1433 }
1434 }
1435
1436 if (!isInStrip) {
1437 continue;
1438 }
1439
1440 unsigned long bcClus = trefTOF.getTime() * Geo::BC_TIME_INPS_INV;
1441
1442 // compute fine correction using cluster position instead of pad center
1443 // this because in case of multiple-hit cluster position is averaged on all pads contributing to the cluster (then error position matrix can be used for Chi2 if nedeed)
1444 int ndigits = 1;
1445 float posCorr[3] = {0, 0, 0};
1446
1447 if (trefTOF.isBitSet(Cluster::kLeft)) {
1448 posCorr[0] += Geo::XPAD, ndigits++;
1449 }
1450 if (trefTOF.isBitSet(Cluster::kUpLeft)) {
1451 posCorr[0] += Geo::XPAD, posCorr[2] -= Geo::ZPAD, ndigits++;
1452 }
1453 if (trefTOF.isBitSet(Cluster::kDownLeft)) {
1454 posCorr[0] += Geo::XPAD, posCorr[2] += Geo::ZPAD, ndigits++;
1455 }
1456 if (trefTOF.isBitSet(Cluster::kUp)) {
1457 posCorr[2] -= Geo::ZPAD, ndigits++;
1458 }
1459 if (trefTOF.isBitSet(Cluster::kDown)) {
1460 posCorr[2] += Geo::ZPAD, ndigits++;
1461 }
1462 if (trefTOF.isBitSet(Cluster::kRight)) {
1463 posCorr[0] -= Geo::XPAD, ndigits++;
1464 }
1465 if (trefTOF.isBitSet(Cluster::kUpRight)) {
1466 posCorr[0] -= Geo::XPAD, posCorr[2] -= Geo::ZPAD, ndigits++;
1467 }
1468 if (trefTOF.isBitSet(Cluster::kDownRight)) {
1469 posCorr[0] -= Geo::XPAD, posCorr[2] += Geo::ZPAD, ndigits++;
1470 }
1471
1472 float ndifInv = 1. / ndigits;
1473 if (ndigits > 1) {
1474 posCorr[0] *= ndifInv;
1475 posCorr[1] *= ndifInv;
1476 posCorr[2] *= ndifInv;
1477 }
1478
1479 int trackIdTOF;
1480 int eventIdTOF;
1481 int sourceIdTOF;
1482 for (auto iPropagation = 0; iPropagation < nStripsCrossedInPropagation[ibc]; iPropagation++) {
1483 if (detId[ibc][iPropagation][1] != indices[1] || detId[ibc][iPropagation][2] != indices[2]) {
1484 continue;
1485 }
1486
1487 int bct0 = int((mTOFClusWork[cacheTOF[itof]].getTime() - trkLTInt[ibc][iPropagation].getTOF(0) + 5000) * Geo::BC_TIME_INPS_INV); // bc taken assuming speed of light (el) and 5 ns of margin
1488 if (bct0 < 0) { // if negative time (it can happen at the beginng of the TF int was truncated per excess... adjusting)
1489 bct0--;
1490 }
1491 float tof = mTOFClusWork[cacheTOF[itof]].getTime() - bct0 * Geo::BC_TIME_INPS;
1492 if (tof - trkLTInt[ibc][iPropagation].getTOF(6) > 2000) { // reject tracks slower than triton
1493 continue;
1494 }
1495
1496 if (mMatchParams->applyPIDcutTPConly) { // for TPC only tracks allowing possibility to apply a PID cut
1497 if (std::abs(tof - trkLTInt[ibc][iPropagation].getTOF(2)) < 2000) { // pion hypotesis
1498 } else if (std::abs(tof - trkLTInt[ibc][iPropagation].getTOF(3)) < 2000) { // kaon hypoteis
1499 } else if (std::abs(tof - trkLTInt[ibc][iPropagation].getTOF(4)) < 2000) { // proton hypotesis
1500 } else { // reject matching
1501 continue;
1502 }
1503 }
1504
1505 float cosangle = TMath::Cos(Geo::getAngles(indices[1], indices[2]) * TMath::DegToRad());
1506 const float errXinvMin = 1. / (mMatchParams->maxResX * mMatchParams->maxResX);
1507 const float errZinvMin = 1. / (mMatchParams->maxResZ * mMatchParams->maxResZ);
1508 float errXinv2 = 1. / (trefTrk.getSigmaY2());
1509 float errZinv2 = 1. / (trefTrk.getSigmaZ2() * cosangle); // should be valid only at eta=0
1510
1511 if (errXinv2 < errXinvMin) {
1512 errXinv2 = errXinvMin;
1513 }
1514 if (errZinv2 < errZinvMin) {
1515 errZinv2 = errZinvMin;
1516 }
1517
1518 LOG(debug) << "TOF Cluster [" << itof << ", " << cacheTOF[itof] << "]: indices = " << indices[0] << ", " << indices[1] << ", " << indices[2] << ", " << indices[3] << ", " << indices[4];
1519 LOG(debug) << "Propagated Track [" << itrk << "]: detId[" << iPropagation << "] = " << detId[ibc][iPropagation][0] << ", " << detId[ibc][iPropagation][1] << ", " << detId[ibc][iPropagation][2] << ", " << detId[ibc][iPropagation][3] << ", " << detId[ibc][iPropagation][4];
1520 float resX = deltaPos[ibc][iPropagation][0] - (indices[4] - detId[ibc][iPropagation][4]) * Geo::XPAD + posCorr[0]; // readjusting the residuals due to the fact that the propagation fell in a pad that was not exactly the one of the cluster
1521 float resZ = deltaPos[ibc][iPropagation][2] - (indices[3] - detId[ibc][iPropagation][3]) * Geo::ZPAD + posCorr[2]; // readjusting the residuals due to the fact that the propagation fell in a pad that was not exactly the one of the cluster
1522 float resY = deltaPos[ibc][iPropagation][1];
1523 if (BCcand[ibc] > bcClus) {
1524 resZ += (BCcand[ibc] - bcClus) * vdriftInBC * side; // add bc correction
1525 } else {
1526 resZ -= (bcClus - BCcand[ibc]) * vdriftInBC * side;
1527 }
1528 float resXor = resX;
1529 float resZor = resZ;
1530 float res = TMath::Sqrt(resX * resX + resZ * resZ);
1531
1532 if (resX < -1.25) { // distance from closest border
1533 resX += 1.25;
1534 } else if (resX > 1.25) {
1535 resX -= 1.25;
1536 } else {
1537 resX = 1E-3 / (pt + 1E-3); // high-pt should be favoured
1538 }
1539
1540 if (resZ < -1.75) { // distance from closest border
1541 resZ += 1.75;
1542 } else if (resZ > 1.75) {
1543 resZ -= 1.75;
1544 } else {
1545 resZ = 1E-3 / (pt + 1E-3); // high-pt should be favoured
1546 }
1547
1548 if (indices[0] != detId[ibc][iPropagation][0]) {
1549 continue;
1550 }
1551 if (indices[1] != detId[ibc][iPropagation][1]) {
1552 continue;
1553 }
1554 if (indices[2] != detId[ibc][iPropagation][2]) {
1555 continue;
1556 }
1557
1558 LOG(debug) << "resX = " << resX << ", resZ = " << resZ << ", res = " << res;
1559 float chi2 = mIsCosmics ? resX : 0.5 * (resX * resX * errXinv2 + resZ * resZ * errZinv2); // TODO: take into account also the time!
1560
1561 if (res < mSpaceTolerance && chi2 < mMatchParams->maxChi2) { // matching ok!
1562 LOG(debug) << "MATCHING FOUND: We have a match! between track " << mTracksSectIndexCache[trkType::UNCONS][sec][itrk] << " and TOF cluster " << mTOFClusSectIndexCache[indices[0]][itof];
1563 foundCluster = true;
1564 // set event indexes (to be checked)
1565
1566 int eventIndexTOFCluster = mTOFClusSectIndexCache[indices[0]][itof];
1567 mMatchedTracksPairsSec[sec].emplace_back(cacheTrk[itrk], eventIndexTOFCluster, mTOFClusWork[cacheTOF[itof]].getTime(), chi2, trkLTInt[ibc][iPropagation], mTrackGid[sec][trkType::UNCONS][cacheTrk[itrk]], trkType::UNCONS, trefTOF.getTime() * 1E-6 - tpctime, trefTOF.getZ(), resXor, resZor, resY); // TODO: check if this is correct!
1568 mMatchedTracksPairsSec[sec][mMatchedTracksPairsSec[sec].size() - 1].setPt(pt);
1569 mMatchedTracksPairsSec[sec][mMatchedTracksPairsSec[sec].size() - 1].setResX(sqrt(1. / errXinv2));
1570 mMatchedTracksPairsSec[sec][mMatchedTracksPairsSec[sec].size() - 1].setResZ(sqrt(1. / errZinv2));
1571 mMatchedTracksPairsSec[sec][mMatchedTracksPairsSec[sec].size() - 1].setResT(resT);
1572 mMatchedTracksPairsSec[sec][mMatchedTracksPairsSec[sec].size() - 1].setVz(mVZtpcOnly[sec][itrk] + Zshift[ibc][iPropagation]);
1573 mMatchedTracksPairsSec[sec][mMatchedTracksPairsSec[sec].size() - 1].setChannel(mainChannel);
1574 }
1575 }
1576 }
1577 }
1578 }
1579 return;
1580}
1581//______________________________________________
1582int MatchTOF::findFITIndex(int bc, const gsl::span<const o2::ft0::RecPoints>& FITRecPoints, unsigned long firstOrbit)
1583{
1584 const auto& FT0Params = o2::ft0::InteractionTag::Instance();
1585
1586 if ((!mHasFillScheme) && o2::tof::Utils::hasFillScheme()) {
1587 mHasFillScheme = true;
1588 for (int ibc = 0; ibc < o2::tof::Utils::getNinteractionBC(); ibc++) {
1589 mFillScheme[o2::tof::Utils::getInteractionBC(ibc)] = true;
1590 }
1591 }
1592
1593 if (FITRecPoints.size() == 0) {
1594 return -1;
1595 }
1596
1597 int index = -1;
1598 int distMax = 10;
1599 bool bestQuality = false; // prioritize FT0 BC with good quality (FT0-AC + vertex) to remove umbiguity in Pb-Pb (not a strict cut because inefficient in pp)
1600 const int distThr = 8;
1601
1602 for (unsigned int i = 0; i < FITRecPoints.size(); i++) {
1603 const auto& ft = FITRecPoints[i];
1604 if (!FT0Params.isSelected(ft)) {
1605 continue;
1606 }
1607 const o2::InteractionRecord ir = FITRecPoints[i].getInteractionRecord();
1608 if (mHasFillScheme && !mFillScheme[ir.bc]) {
1609 continue;
1610 }
1611 bool quality = (std::abs(FITRecPoints[i].getCollisionTime(0)) < 1000 && std::abs(FITRecPoints[i].getVertex()) < 1000);
1612 if (bestQuality && !quality) { // if current has no good quality and the one previoulsy selected has -> discard the current one
1613 continue;
1614 }
1615 int bct0 = (ir.orbit - firstOrbit) * o2::constants::lhc::LHCMaxBunches + ir.bc;
1616 int dist = bc - bct0;
1617
1618 bool worseDistance = dist < 0 || dist > distThr || dist > distMax;
1619 if (worseDistance) { // discard if BC is not in the proper range or it is worse than the one already found
1620 continue;
1621 }
1622
1623 bestQuality = quality;
1624 distMax = dist;
1625 index = i;
1626 }
1627
1628 return index;
1629}
1630//______________________________________________
1631void MatchTOF::selectBestMatches(int sec)
1632{
1633 if (mSetHighPurity) {
1634 BestMatchesHP(mMatchedTracksPairsSec[sec], mMatchedTracks, mMatchedTracksIndex[sec], mMatchedClustersIndex, mFITRecPoints, mTOFClusWork, mCalibInfoTOF, mTimestamp, mMCTruthON, mTOFClusLabels, mTracksLblWork[sec], mOutTOFLabels);
1635 return;
1636 }
1637 BestMatches(mMatchedTracksPairsSec[sec], mMatchedTracks, mMatchedTracksIndex[sec], mMatchedClustersIndex, mFITRecPoints, mTOFClusWork, mTracksWork[sec], mCalibInfoTOF, mTimestamp, mMCTruthON, mTOFClusLabels, mTracksLblWork[sec], mOutTOFLabels, mMatchParams->calibMaxChi2);
1638}
1639//______________________________________________
1640void MatchTOF::BestMatches(std::vector<o2::dataformats::MatchInfoTOFReco>& matchedTracksPairs, std::vector<o2::dataformats::MatchInfoTOF>* matchedTracks, std::vector<int>* matchedTracksIndex, int* matchedClustersIndex, const gsl::span<const o2::ft0::RecPoints>& FITRecPoints, const std::vector<Cluster>& TOFClusWork, const std::vector<matchTrack>* TracksWork, std::vector<o2::dataformats::CalibInfoTOF>& CalibInfoTOF, unsigned long Timestamp, bool MCTruthON, const o2::dataformats::MCTruthContainer<o2::MCCompLabel>* TOFClusLabels, const std::vector<o2::MCCompLabel>* TracksLblWork, std::vector<o2::MCCompLabel>* OutTOFLabels, float calibMaxChi2)
1641{
1643
1644 // first, we sort according to the chi2
1645 std::sort(matchedTracksPairs.begin(), matchedTracksPairs.end(), [](const o2::dataformats::MatchInfoTOFReco& a, const o2::dataformats::MatchInfoTOFReco& b) { return (a.getChi2() < b.getChi2()); });
1646 int i = 0;
1647
1648 // then we take discard the pairs if their track or cluster was already matched (since they are ordered in chi2, we will take the best matching)
1649 for (o2::dataformats::MatchInfoTOFReco& matchingPair : matchedTracksPairs) {
1650 int trkType = (int)matchingPair.getTrackType();
1651
1652 int itrk = matchingPair.getIdLocal();
1653
1654 int trkTypeSplitted = trkType;
1655 auto sourceID = matchingPair.getTrackRef().getSource();
1657 trkTypeSplitted = (int)trkType::TPCTRD;
1658 } else if (sourceID == o2::dataformats::GlobalTrackID::ITSTPCTRD) {
1659 trkTypeSplitted = (int)trkType::ITSTPCTRD;
1660 }
1661
1662 if (matchedClustersIndex[matchingPair.getTOFClIndex()] != -1) { // the cluster was already filled
1663 continue;
1664 }
1665 if (matchedTracksIndex[trkType][itrk] != -1) { // the track was already filled
1666 if (trkType) { // not applied for unconstrained tracks (trkType == 0)
1667 // let's check if we can use both matching to improve TOF time
1668 int pairInd = matchedTracksIndex[trkType][itrk];
1669 auto& prevMatching = matchedTracks[trkTypeSplitted][pairInd];
1670
1671 if (pairInd >= matchedTracks[trkTypeSplitted].size()) {
1672 LOG(error) << "Pair index out of range when trying to update TOF time. This should not happen";
1673 continue;
1674 }
1675
1676 int istripOld = TOFClusWork[prevMatching.getTOFClIndex()].getMainContributingChannel() / o2::tof::Geo::NPADS;
1677 int istripNew = TOFClusWork[matchingPair.getTOFClIndex()].getMainContributingChannel() / o2::tof::Geo::NPADS;
1678
1679 if (istripOld != istripNew) { // we accept it only if in a different strip
1680 const o2::track::TrackLTIntegral& intLTnew = matchingPair.getLTIntegralOut();
1681 const o2::track::TrackLTIntegral& intLTold = prevMatching.getLTIntegralOut();
1682
1683 const float cinv = 1.E+10 / TMath::C(); // cm/ps
1684 float deltaT = (intLTnew.getL() - intLTold.getL()) * cinv;
1685
1686 float timeNew = TOFClusWork[matchingPair.getTOFClIndex()].getTime() - deltaT;
1687 float timeOld = TOFClusWork[prevMatching.getTOFClIndex()].getTime();
1688
1689 if (std::abs(timeNew - timeOld) < 200) {
1690 // update time information averaging the two (the second one corrected for the difference in the track length)
1691 prevMatching.setSignal((timeNew + timeOld) * 0.5);
1692 float geanttime = (TOFClusWork[matchingPair.getTOFClIndex()].getTgeant() + TOFClusWork[prevMatching.getTOFClIndex()].getTgeant() - deltaT * 1E-3) * 0.5;
1693 double t0 = (TOFClusWork[matchingPair.getTOFClIndex()].getT0true() + TOFClusWork[prevMatching.getTOFClIndex()].getT0true()) * 0.5;
1694 prevMatching.setTgeant(geanttime);
1695 prevMatching.setT0true(t0);
1696 prevMatching.setChi2(0); // flag such cases with chi2 equal to zero
1697 matchedClustersIndex[matchingPair.getTOFClIndex()] = matchedTracksIndex[trkType][itrk]; // flag also the second cluster as already used
1698 }
1699 }
1700 }
1701
1702 continue;
1703 }
1704 matchedTracksIndex[trkType][itrk] = matchedTracks[trkTypeSplitted].size(); // index of the MatchInfoTOF correspoding to this track
1705 matchedClustersIndex[matchingPair.getTOFClIndex()] = matchedTracksIndex[trkType][itrk]; // index of the track that was matched to this cluster
1706
1707 matchingPair.setTgeant(TOFClusWork[matchingPair.getTOFClIndex()].getTgeant());
1708 matchingPair.setT0true(TOFClusWork[matchingPair.getTOFClIndex()].getT0true());
1709
1710 // let's check if cluster has multiple-hits (noferini)
1711 const auto& tofcl = TOFClusWork[matchingPair.getTOFClIndex()];
1712 if (tofcl.getNumOfContributingChannels() > 1) {
1713 // has an additional hit Up or Down (Z-dir)
1714 matchingPair.setHitPatternUpDown(tofcl.isAdditionalChannelSet(o2::tof::Cluster::kUp) ||
1715 tofcl.isAdditionalChannelSet(o2::tof::Cluster::kUpLeft) ||
1716 tofcl.isAdditionalChannelSet(o2::tof::Cluster::kUpRight) ||
1717 tofcl.isAdditionalChannelSet(o2::tof::Cluster::kDown) ||
1718 tofcl.isAdditionalChannelSet(o2::tof::Cluster::kDownLeft) ||
1719 tofcl.isAdditionalChannelSet(o2::tof::Cluster::kDownRight));
1720 // has an additional hit Left or Right (X-dir)
1721 matchingPair.setHitPatternLeftRight(tofcl.isAdditionalChannelSet(o2::tof::Cluster::kLeft) ||
1722 tofcl.isAdditionalChannelSet(o2::tof::Cluster::kDownLeft) ||
1723 tofcl.isAdditionalChannelSet(o2::tof::Cluster::kUpLeft) ||
1724 tofcl.isAdditionalChannelSet(o2::tof::Cluster::kRight) ||
1725 tofcl.isAdditionalChannelSet(o2::tof::Cluster::kDownRight) ||
1726 tofcl.isAdditionalChannelSet(o2::tof::Cluster::kUpRight));
1727 }
1728
1729 // estimate collision time using FT0 info if available
1730 ULong64_t bclongtofCal = (matchingPair.getSignal() - 10000) * o2::tof::Geo::BC_TIME_INPS_INV;
1731 double t0Best = bclongtofCal * o2::tof::Geo::BC_TIME_INPS; // here just BC
1732 float t0BestRes = 200;
1733 if (FITRecPoints.size() > 0) {
1734 int index = findFITIndex(bclongtofCal, FITRecPoints, mFirstTForbit);
1735 if (index > -1 && FITRecPoints[index].isValidTime(1) && FITRecPoints[index].isValidTime(2)) { // require A and C
1736 t0Best += FITRecPoints[index].getCollisionTime(0);
1737 t0BestRes = 15;
1738 }
1739 }
1740 matchingPair.setFT0Best(t0Best, t0BestRes);
1741 matchedTracks[trkTypeSplitted].push_back(matchingPair); // array of MatchInfoTOF
1742
1743 // get fit info
1744 double t0info = 0;
1745
1746 const o2::track::TrackLTIntegral& intLT = matchingPair.getLTIntegralOut();
1747
1748 if (FITRecPoints.size() > 0 && mIsFIT) {
1749 int index = findFITIndex(TOFClusWork[matchingPair.getTOFClIndex()].getBC() - o2::tof::Geo::LATENCYWINDOW_IN_BC, FITRecPoints, mFirstTForbit);
1750
1751 if (index > -1) {
1752 o2::InteractionRecord ir = FITRecPoints[index].getInteractionRecord();
1753 int bct0 = (ir.orbit - mFirstTForbit) * o2::constants::lhc::LHCMaxBunches + ir.bc;
1754 t0info = bct0 * o2::tof::Geo::BC_TIME_INPS;
1755 }
1756 } else if (!mIsFIT) { // move time to time in orbit to avoid loss of precision when truncating from double to float
1757 int bcStarOrbit = int((TOFClusWork[matchingPair.getTOFClIndex()].getTimeRaw() - intLT.getTOF(o2::track::PID::Pion)) * o2::tof::Geo::BC_TIME_INPS_INV);
1758 bcStarOrbit = (bcStarOrbit / o2::constants::lhc::LHCMaxBunches) * o2::constants::lhc::LHCMaxBunches; // truncation
1759 t0info = bcStarOrbit * o2::tof::Geo::BC_TIME_INPS;
1760 }
1761
1762 // add also calibration infos
1763 int flags = 0;
1764 if (sourceID == o2::dataformats::GlobalTrackID::TPC) {
1766 } else if (sourceID == o2::dataformats::GlobalTrackID::ITSTPC) {
1768 } else if (sourceID == o2::dataformats::GlobalTrackID::TPCTRD) {
1770 } else if (sourceID == o2::dataformats::GlobalTrackID::ITSTPCTRD) {
1772 }
1773
1774 int mask = 0;
1775 float deltat;
1776
1777 if (t0info > 0 && mIsFIT) { // FT0 BC found
1778 deltat = TOFClusWork[matchingPair.getTOFClIndex()].getTimeRaw() - t0info - intLT.getTOF(o2::track::PID::Pion) - o2::tof::Geo::LATENCYWINDOW_IN_BC * o2::tof::Geo::BC_TIME_INPS; // latency subtracted
1779 mask = 1 << 8; // only FT0 BC allowed (-> BC=0 since t0info already subtracted and assumed to be aligned to the true IntBC)
1780 } else if (!mIsFIT) {
1781 deltat = o2::tof::Utils::subtractInteractionBC(TOFClusWork[matchingPair.getTOFClIndex()].getTimeRaw() - t0info - intLT.getTOF(o2::track::PID::Pion), mask, true);
1782 }
1783
1784 const o2::track::TrackParCov& trc = TracksWork[trkType][itrk].first;
1785 float pt = trc.getPt(); // from outer parameters!
1786
1787 if (pt > 1.5) {
1789 }
1790
1791 if (pt < 0.5) {
1793 }
1794
1795 if (mask == 0) {
1797 }
1798
1799 if (TOFClusWork[matchingPair.getTOFClIndex()].getNumOfContributingChannels() != 1) {
1801 }
1802
1803 if (matchingPair.getChi2() < calibMaxChi2 && t0info > 0) { // extra cut in ChiSquare for storing calib info
1804 CalibInfoTOF.emplace_back(TOFClusWork[matchingPair.getTOFClIndex()].getMainContributingChannel(),
1805 Timestamp / 1000 + int(TOFClusWork[matchingPair.getTOFClIndex()].getTimeRaw() * 1E-12), // add time stamp
1806 deltat,
1807 TOFClusWork[matchingPair.getTOFClIndex()].getTot(), mask, flags);
1808 }
1809
1810 if (MCTruthON) {
1811 const auto& labelsTOF = TOFClusLabels->getLabels(matchingPair.getTOFClIndex());
1812 auto& labelTrack = TracksLblWork[trkType][itrk];
1813 // we have not found the track label among those associated to the TOF cluster --> fake match! We will associate the label of the main channel, but negative
1814 bool fake = true;
1815 for (auto& lbl : labelsTOF) {
1816 if (labelTrack == lbl) { // compares src, evID, trID, ignores fake flag.
1817 fake = false;
1818 }
1819 }
1820 OutTOFLabels[trkTypeSplitted].emplace_back(labelTrack).setFakeFlag(fake);
1821 }
1822 i++;
1823 }
1824}
1825//______________________________________________
1826void MatchTOF::BestMatchesHP(std::vector<o2::dataformats::MatchInfoTOFReco>& matchedTracksPairs, std::vector<o2::dataformats::MatchInfoTOF>* matchedTracks, std::vector<int>* matchedTracksIndex, int* matchedClustersIndex, const gsl::span<const o2::ft0::RecPoints>& FITRecPoints, const std::vector<Cluster>& TOFClusWork, std::vector<o2::dataformats::CalibInfoTOF>& CalibInfoTOF, unsigned long Timestamp, bool MCTruthON, const o2::dataformats::MCTruthContainer<o2::MCCompLabel>* TOFClusLabels, const std::vector<o2::MCCompLabel>* TracksLblWork, std::vector<o2::MCCompLabel>* OutTOFLabels)
1827{
1829 float chi2SeparationCut = 2;
1830 float chi2S = 3;
1831
1832 std::vector<o2::dataformats::MatchInfoTOFReco> tmpMatch;
1833
1834 // first, we sort according to the chi2
1835 std::sort(matchedTracksPairs.begin(), matchedTracksPairs.end(), [](const o2::dataformats::MatchInfoTOFReco& a, const o2::dataformats::MatchInfoTOFReco& b) { return (a.getChi2() < b.getChi2()); });
1836 int i = 0;
1837 // then we take discard the pairs if their track or cluster was already matched (since they are ordered in chi2, we will take the best matching)
1838 for (const o2::dataformats::MatchInfoTOFReco& matchingPair : matchedTracksPairs) {
1839 int trkType = (int)matchingPair.getTrackType();
1840
1841 int itrk = matchingPair.getIdLocal();
1842
1843 bool discard = matchingPair.getChi2() > chi2S;
1844
1845 if (matchedTracksIndex[trkType][itrk] != -1) { // the track was already filled, check if this competitor is not too close
1846 auto winnerChi = tmpMatch[matchedTracksIndex[trkType][itrk]].getChi2();
1847 if (winnerChi < 0) { // the winner was already discarded as ambiguous
1848 continue;
1849 }
1850 if (matchingPair.getChi2() - winnerChi < chi2SeparationCut) { // discard previously validated winner and it has too close competitor
1851 tmpMatch[matchedTracksIndex[trkType][itrk]].setChi2(-1);
1852 }
1853 continue;
1854 }
1855
1856 if (matchedClustersIndex[matchingPair.getTOFClIndex()] != -1) { // the cluster was already filled, check if this competitor is not too close
1857 auto winnerChi = tmpMatch[matchedClustersIndex[matchingPair.getTOFClIndex()]].getChi2();
1858 if (winnerChi < 0) { // the winner was already discarded as ambiguous
1859 continue;
1860 }
1861 if (matchingPair.getChi2() - winnerChi < chi2SeparationCut) { // discard previously validated winner and it has too close competitor
1862 tmpMatch[matchedClustersIndex[matchingPair.getTOFClIndex()]].setChi2(-1);
1863 }
1864 continue;
1865 }
1866
1867 if (!discard) {
1868 matchedTracksIndex[trkType][itrk] = tmpMatch.size(); // index of the MatchInfoTOF correspoding to this track
1869 matchedClustersIndex[matchingPair.getTOFClIndex()] = matchedTracksIndex[trkType][itrk]; // index of the track that was matched to this clus
1870 tmpMatch.push_back(matchingPair);
1871 }
1872 }
1873
1874 // now write final matches skipping disabled ones
1875 for (auto& matchingPair : tmpMatch) {
1876 if (matchingPair.getChi2() <= 0) {
1877 continue;
1878 }
1879 int trkType = (int)matchingPair.getTrackType();
1880 int itrk = matchingPair.getIdLocal();
1881
1882 int trkTypeSplitted = trkType;
1883 auto sourceID = matchingPair.getTrackRef().getSource();
1885 trkTypeSplitted = (int)trkType::TPCTRD;
1886 } else if (sourceID == o2::dataformats::GlobalTrackID::ITSTPCTRD) {
1887 trkTypeSplitted = (int)trkType::ITSTPCTRD;
1888 }
1889 matchedTracks[trkTypeSplitted].push_back(matchingPair); // array of MatchInfoTOF
1890
1891 const o2::track::TrackLTIntegral& intLT = matchingPair.getLTIntegralOut();
1892
1893 // get fit info
1894 double t0info = 0;
1895
1896 if (FITRecPoints.size() > 0 && mIsFIT) {
1897 int index = findFITIndex(TOFClusWork[matchingPair.getTOFClIndex()].getBC() - o2::tof::Geo::LATENCYWINDOW_IN_BC, FITRecPoints, mFirstTForbit);
1898
1899 if (index > -1) {
1900 o2::InteractionRecord ir = FITRecPoints[index].getInteractionRecord();
1901 int bct0 = (ir.orbit - mFirstTForbit) * o2::constants::lhc::LHCMaxBunches + ir.bc;
1902 t0info = bct0 * o2::tof::Geo::BC_TIME_INPS;
1903 }
1904 } else { // move time to time in orbit to avoid loss of precision when truncating from double to float
1905 int bcStarOrbit = int((TOFClusWork[matchingPair.getTOFClIndex()].getTimeRaw() - intLT.getTOF(o2::track::PID::Pion)) * o2::tof::Geo::BC_TIME_INPS_INV);
1906 bcStarOrbit = (bcStarOrbit / o2::constants::lhc::LHCMaxBunches) * o2::constants::lhc::LHCMaxBunches; // truncation
1907 t0info = bcStarOrbit * o2::tof::Geo::BC_TIME_INPS;
1908 }
1909
1910 // add also calibration infos
1912 CalibInfoTOF.emplace_back(TOFClusWork[matchingPair.getTOFClIndex()].getMainContributingChannel(),
1913 Timestamp / 1000 + int(TOFClusWork[matchingPair.getTOFClIndex()].getTimeRaw() * 1E-12), // add time stamp
1914 TOFClusWork[matchingPair.getTOFClIndex()].getTimeRaw() - t0info - intLT.getTOF(o2::track::PID::Pion),
1915 TOFClusWork[matchingPair.getTOFClIndex()].getTot(), 0);
1916 }
1917
1918 if (MCTruthON) {
1919 const auto& labelsTOF = TOFClusLabels->getLabels(matchingPair.getTOFClIndex());
1920 auto& labelTrack = TracksLblWork[trkType][itrk];
1921 // we have not found the track label among those associated to the TOF cluster --> fake match! We will associate the label of the main channel, but negative
1922 bool fake = true;
1923 for (auto& lbl : labelsTOF) {
1924 if (labelTrack == lbl) { // compares src, evID, trID, ignores fake flag.
1925 fake = false;
1926 }
1927 }
1928 OutTOFLabels[trkTypeSplitted].emplace_back(labelTrack).setFakeFlag(fake);
1929 }
1930 }
1931}
1932//______________________________________________
1933bool MatchTOF::propagateToRefX(o2::track::TrackParCov& trc, float xRef, float stepInCm, o2::track::TrackLTIntegral& intLT)
1934{
1935 // propagate track to matching reference X
1936 o2::base::Propagator::MatCorrType matCorr = o2::base::Propagator::MatCorrType::USEMatCorrLUT; // material correction method
1937 static const float tanHalfSector = tan(o2::constants::math::SectorSpanRad / 2);
1938 bool refReached = false;
1939 float xStart = trc.getX();
1940 // the first propagation will be from 2m, if the track is not at least at 2m
1941 if (xStart < 50.) {
1942 xStart = 50.;
1943 }
1944 int istep = 1;
1945 bool hasPropagated = o2::base::Propagator::Instance()->PropagateToXBxByBz(trc, xStart + istep * stepInCm, MAXSNP, stepInCm, matCorr, &intLT);
1946 while (hasPropagated) {
1947 if (trc.getX() > xRef) {
1948 refReached = true; // we reached the 371cm reference
1949 }
1950 istep++;
1951 if (std::abs(trc.getY()) > trc.getX() * tanHalfSector) { // we are still in the same sector
1952 // we need to rotate the track to go to the new sector
1953 //Printf("propagateToRefX: changing sector");
1954 auto alphaNew = o2::math_utils::angle2Alpha(trc.getPhiPos());
1955 if (!trc.rotate(alphaNew) != 0) {
1956 // Printf("propagateToRefX: failed to rotate");
1957 break; // failed (this line is taken from MatchTPCITS and the following comment too: RS: check effect on matching tracks to neighbouring sector)
1958 }
1959 }
1960 if (refReached) {
1961 break;
1962 }
1963 hasPropagated = o2::base::Propagator::Instance()->PropagateToXBxByBz(trc, xStart + istep * stepInCm, MAXSNP, stepInCm, matCorr, &intLT);
1964 }
1965
1966 // if (std::abs(trc.getSnp()) > MAXSNP) Printf("propagateToRefX: condition on snp not ok, returning false");
1967 //Printf("propagateToRefX: snp of teh track is %f (--> %f grad)", trc.getSnp(), TMath::ASin(trc.getSnp())*TMath::RadToDeg());
1968 return refReached && std::abs(trc.getSnp()) < 0.95; // Here we need to put MAXSNP
1969}
1970//______________________________________________
1971bool MatchTOF::propagateToRefXWithoutCov(const o2::track::TrackParCov& trc, float xRef, float stepInCm, float bzField)
1972{
1973 // propagate track to matching reference X without using the covariance matrix
1974 // we create the copy of the track in a TrackPar object (no cov matrix)
1975 o2::track::TrackPar trcNoCov(trc);
1976 const float tanHalfSector = tan(o2::constants::math::SectorSpanRad / 2);
1977 bool refReached = false;
1978 float xStart = trcNoCov.getX();
1979 // the first propagation will be from 2m, if the track is not at least at 2m
1980 if (xStart < 50.) {
1981 xStart = 50.;
1982 }
1983 int istep = 1;
1984 bool hasPropagated = trcNoCov.propagateParamTo(xStart + istep * stepInCm, bzField);
1985 while (hasPropagated) {
1986 if (trcNoCov.getX() > xRef) {
1987 refReached = true; // we reached the 371cm reference
1988 }
1989 istep++;
1990 if (std::abs(trcNoCov.getY()) > trcNoCov.getX() * tanHalfSector) { // we are still in the same sector
1991 // we need to rotate the track to go to the new sector
1992 //Printf("propagateToRefX: changing sector");
1993 auto alphaNew = o2::math_utils::angle2Alpha(trcNoCov.getPhiPos());
1994 if (!trcNoCov.rotateParam(alphaNew) != 0) {
1995 // Printf("propagateToRefX: failed to rotate");
1996 break; // failed (this line is taken from MatchTPCITS and the following comment too: RS: check effect on matching tracks to neighbouring sector)
1997 }
1998 }
1999 if (refReached) {
2000 break;
2001 }
2002 hasPropagated = trcNoCov.propagateParamTo(xStart + istep * stepInCm, bzField);
2003 }
2004 // if (std::abs(trc.getSnp()) > MAXSNP) Printf("propagateToRefX: condition on snp not ok, returning false");
2005 //Printf("propagateToRefX: snp of teh track is %f (--> %f grad)", trcNoCov.getSnp(), TMath::ASin(trcNoCov.getSnp())*TMath::RadToDeg());
2006
2007 return refReached && std::abs(trcNoCov.getSnp()) < 0.95 && std::abs(trcNoCov.getZ()) < Geo::MAXHZTOF; // Here we need to put MAXSNP
2008}
2009//______________________________________________
2010void MatchTOF::updateTL(o2::track::TrackLTIntegral& intLT, float deltal)
2011{
2012 for (int i = 0; i < intLT.getNTOFs(); i++) {
2013 float betainv = intLT.getTOF(i) / intLT.getL();
2014 intLT.setTOF(intLT.getTOF(i) + deltal * betainv, i);
2015 }
2016 intLT.setL(intLT.getL() + deltal);
2017}
2018
2019//______________________________________________
2020bool MatchTOF::propagateToRefXWithoutCov(const o2::track::TrackParCov& trc, float xRef, float stepInCm, float bzField, float pos[3])
2021{
2022 // propagate track to matching reference X without using the covariance matrix
2023 // we create the copy of the track in a TrackPar object (no cov matrix)
2024 o2::track::TrackPar trcNoCov(trc);
2025 const float tanHalfSector = tan(o2::constants::math::SectorSpanRad / 2);
2026 bool refReached = false;
2027 float xStart = trcNoCov.getX();
2028 // the first propagation will be from 2m, if the track is not at least at 2m
2029 if (xStart < 50.) {
2030 xStart = 50.;
2031 }
2032 int istep = 1;
2033 bool hasPropagated = trcNoCov.propagateParamTo(xStart + istep * stepInCm, bzField);
2034 while (hasPropagated) {
2035 if (trcNoCov.getX() > xRef) {
2036 refReached = true; // we reached the 371cm reference
2037 }
2038 istep++;
2039 if (fabs(trcNoCov.getY()) > trcNoCov.getX() * tanHalfSector) { // we are still in the same sector
2040 // we need to rotate the track to go to the new sector
2041 // Printf("propagateToRefX: changing sector");
2042 auto alphaNew = o2::math_utils::angle2Alpha(trcNoCov.getPhiPos());
2043 if (!trcNoCov.rotateParam(alphaNew) != 0) {
2044 // Printf("propagateToRefX: failed to rotate");
2045 break; // failed (this line is taken from MatchTPCITS and the following comment too: RS: check effect on matching tracks to neighbouring sector)
2046 }
2047 }
2048 if (refReached) {
2049 break;
2050 }
2051 hasPropagated = trcNoCov.propagateParamTo(xStart + istep * stepInCm, bzField);
2052 }
2053 // if (std::abs(trc.getSnp()) > MAXSNP) Printf("propagateToRefX: condition on snp not ok, returning false");
2054 // Printf("propagateToRefX: snp of teh track is %f (--> %f grad)", trcNoCov.getSnp(), TMath::ASin(trcNoCov.getSnp())*TMath::RadToDeg());
2055
2057 trcNoCov.getXYZGlo(xyz);
2058 pos[0] = xyz[0];
2059 pos[1] = xyz[1];
2060 pos[2] = xyz[2];
2061
2062 return refReached && std::abs(trcNoCov.getSnp()) < 0.95 && TMath::Abs(trcNoCov.getZ()) < Geo::MAXHZTOF; // Here we need to put MAXSNP
2063}
2064
2065//______________________________________________
2066void MatchTOF::setDebugFlag(UInt_t flag, bool on)
2067{
2069 if (on) {
2070 mDBGFlags |= flag;
2071 } else {
2072 mDBGFlags &= ~flag;
2073 }
2074}
2075//______________________________________________
2076void MatchTOF::updateTimeDependentParams()
2077{
2080 mTPCTBinMUS = elParam.ZbinWidth; // TPC bin in microseconds
2081 mTPCTBinMUSInv = 1. / mTPCTBinMUS;
2082 mTPCBin2Z = mTPCTBinMUS * mTPCVDrift;
2083
2084 mBz = o2::base::Propagator::Instance()->getNominalBz();
2085 mMaxInvPt = std::abs(mBz) > 0.1 ? 1. / (std::abs(mBz) * 0.05) : 999.;
2086
2087 const auto& trackTune = TrackTuneParams::Instance();
2088 float scale = mTPCCorrMapsHelper->getInstLumiCTP();
2089 if (scale < 0.f) {
2090 scale = 0.f;
2091 }
2092 mCovDiagInner = trackTune.getCovInnerTotal(scale);
2093 mCovDiagOuter = trackTune.getCovOuterTotal(scale);
2094}
2095
2096//_________________________________________________________
2098{
2099 auto& match = mMatchedTracks[trkType::TPC][matchedID];
2100 const auto& tpcTrOrig = mRecoCont->getTPCTrack(match.getTrackRef());
2101 const auto& tofCl = mTOFClustersArrayInp[match.getTOFClIndex()];
2102 const auto& intLT = match.getLTIntegralOut();
2103
2104 // correct the time of the track
2105 auto timeTOFMUS = (tofCl.getTime() - intLT.getTOF(tpcTrOrig.getPID())) * 1e-6; // tof time in \mus, FIXME: account for time of flight to R TOF
2106 auto timeTOFTB = timeTOFMUS * mTPCTBinMUSInv; // TOF time in TPC timebins, including TPC time offset
2107 auto deltaTBins = timeTOFTB - tpcTrOrig.getTime0(); // time shift in timeBins
2108 float timeErr = 0.010; // assume 10 ns error FIXME
2109 auto dzCorr = deltaTBins * mTPCBin2Z;
2110
2111 if (mTPCClusterIdxStruct) { // refit was requested
2112 trConstr.o2::track::TrackParCov::operator=(tpcTrOrig.getOuterParam()); // seed for inward refit of constrained track, made from the outer param
2113 trConstr.setParamOut(tpcTrOrig); // seed for outward refit of constrained track, made from the inner param
2114 } else {
2115 trConstr.o2::track::TrackParCov::operator=(tpcTrOrig); // inner param, we just correct its position, w/o refit
2116 trConstr.setParamOut(tpcTrOrig.getOuterParam()); // outer param, we just correct its position, w/o refit
2117 }
2118
2119 auto& trConstrOut = trConstr.getParamOut();
2120
2121 auto zTrack = trConstr.getZ();
2122 auto zTrackOut = trConstrOut.getZ();
2123
2124 if (tpcTrOrig.hasASideClustersOnly()) {
2125 zTrack += dzCorr;
2126 zTrackOut += dzCorr;
2127 } else if (tpcTrOrig.hasCSideClustersOnly()) {
2128 zTrack -= dzCorr;
2129 zTrackOut -= dzCorr;
2130 } else {
2131 // TODO : special treatment of tracks crossing the CE
2132 }
2133 trConstr.setZ(zTrack);
2134 trConstrOut.setZ(zTrackOut);
2135 //
2136 trConstr.setTimeMUS(timeTOFMUS, timeErr);
2137 trConstr.setRefMatch(matchedID);
2138 if (mTPCClusterIdxStruct) { // refit was requested
2139 float chi2 = 0;
2140 mTPCRefitter->setTrackReferenceX(o2::constants::geom::XTPCInnerRef);
2141 if (mTPCRefitter->RefitTrackAsTrackParCov(trConstr, tpcTrOrig.getClusterRef(), timeTOFTB, &chi2, false, true) < 0) { // inward refit after resetting cov.mat.
2142 LOGP(debug, "Inward Refit failed {}", trConstr.asString());
2143 return false;
2144 }
2145 const auto& trackTune = TrackTuneParams::Instance(); // if needed, correct TPC track after inward refit
2146 if (!trackTune.useTPCInnerCorr) {
2147 trConstr.updateParams(trackTune.tpcParInner);
2148 }
2149 if (trackTune.tpcCovInnerType != TrackTuneParams::AddCovType::Disable) {
2150 trConstr.updateCov(mCovDiagInner, trackTune.tpcCovInnerType == TrackTuneParams::AddCovType::WithCorrelations);
2151 }
2152
2153 trConstr.setChi2Refit(chi2);
2154 //
2155 mTPCRefitter->setTrackReferenceX(o2::constants::geom::XTPCOuterRef);
2156 if (mTPCRefitter->RefitTrackAsTrackParCov(trConstrOut, tpcTrOrig.getClusterRef(), timeTOFTB, &chi2, true, true) < 0) { // outward refit after resetting cov.mat.
2157 LOGP(debug, "Outward Refit failed {}", trConstrOut.asString());
2158 return false;
2159 }
2160 }
2161
2162 return true;
2163}
2164
2165//_________________________________________________________
2167{
2168 if (mTPCClusterIdxStruct) {
2169 mTPCRefitter = std::make_unique<o2::gpu::GPUO2InterfaceRefit>(mTPCClusterIdxStruct, mTPCCorrMapsHelper, mBz,
2170 mTPCTrackClusIdx.data(), 0, mTPCRefitterShMap.data(),
2171 mTPCRefitterOccMap.data(), mTPCRefitterOccMap.size(), nullptr, o2::base::Propagator::Instance());
2172 }
2173}
General auxilliary methods.
Wrapper container for different reconstructed object types.
particle ids, masses, names class definition
Definition of the GeometryManager class.
uint64_t bc
Definition RawEventData.h:5
int32_t i
Header of the General Run Parameters object.
Some ALICE geometry constants of common interest.
Definition of a container to keep Monte Carlo truth external to simulation objects.
Definition of the fast magnetic field parametrization MagFieldFast.
Definition of the MagF class.
ClassImp(MatchTOF)
Class to perform TOF matching to global tracks.
useful math constants
Definition of the parameter class for the detector electronics.
Definition of the parameter class for the detector gas.
Header to collect physics constants.
uint16_t pos
Definition RawData.h:3
uint32_t side
Definition RawData.h:0
uint32_t res
Definition RawData.h:0
Wrapper container for different reconstructed object types.
class to create TPC fast transformation
Track Length and TOF integral.
Configurable params for tracks ad hoc tuning.
calibration data from laser track calibration
std::ostringstream debug
GPUd() value_type estimateLTFast(o2 static GPUd() float estimateLTIncrement(const o2 PropagatorImpl * Instance(bool uninitialized=false)
Definition Propagator.h:143
A container to hold and manage MC truth information/labels.
gsl::span< TruthElement > getLabels(uint32_t dataindex)
o2::track::TrackParCov & getParamOut()
Definition TrackTPCTOF.h:53
void setTimeMUS(const timeEst &t)
Definition TrackTPCTOF.h:43
void setParamOut(const o2::track::TrackParCov &v)
Definition TrackTPCTOF.h:55
void setTPCVDrift(const o2::tpc::VDriftCorrFact &v)
Definition MatchTOF.cxx:249
void setTPCCorrMaps(o2::gpu::CorrectionMapsHelper *maph)
Definition MatchTOF.cxx:258
void printCandidatesTOF() const
set time tolerance on track-TOF times comparison
Definition MatchTOF.cxx:278
void run(const o2::globaltracking::RecoContainer &inp, unsigned long firstTForbit=0)
< perform matching for provided input
Definition MatchTOF.cxx:67
void setDebugFlag(UInt_t flag, bool on=true)
set the name of output debug file
static void groupingMatch(const std::vector< o2::dataformats::MatchInfoTOFReco > &origin, std::vector< std::vector< o2::dataformats::MatchInfoTOFReco > > &grouped, std::vector< std::vector< int > > &firstEls, std::vector< std::vector< int > > &secondEls)
Definition MatchTOF.cxx:317
static void printGrouping(const std::vector< o2::dataformats::MatchInfoTOFReco > &origin, const std::vector< std::vector< o2::dataformats::MatchInfoTOFReco > > &grouped)
Definition MatchTOF.cxx:299
static int findFITIndex(int bc, const gsl::span< const o2::ft0::RecPoints > &FITRecPoints, unsigned long firstOrbit)
bool makeConstrainedTPCTrack(int matchedID, o2::dataformats::TrackTPCTOF &trConstr)
populate externally provided container by TOF-time-constrained TPC tracks
HMPID cluster implementation.
Definition Cluster.h:27
Cluster class for TOF.
Definition Cluster.h:37
static constexpr Float_t RMAX
Definition Geo.h:136
static constexpr Float_t ZPAD
Definition Geo.h:143
static constexpr Float_t XPAD
Definition Geo.h:141
static constexpr Double_t BC_TIME_INPS
Definition Geo.h:105
static constexpr Double_t BC_TIME_INPS_INV
Definition Geo.h:106
static Float_t getAngles(Int_t iplate, Int_t istrip)
Definition Geo.h:81
static constexpr Int_t NPADS
Definition Geo.h:112
static void Init()
Definition Geo.cxx:43
static constexpr Int_t LATENCYWINDOW_IN_BC
Definition Geo.h:174
static constexpr Float_t MAXHZTOF
Definition Geo.h:129
static void getVolumeIndices(Int_t index, Int_t *detId)
Definition Geo.cxx:543
static void getPadDxDyDz(const Float_t *pos, Int_t *det, Float_t *DeltaPos, int sector=-1)
Definition Geo.cxx:835
static double subtractInteractionBC(double time, int &mask, bool subLatency=false)
Definition Utils.cxx:122
static int getNinteractionBC()
Definition Utils.cxx:117
static bool hasFillScheme()
Definition Utils.cxx:261
static int getInteractionBC(int ibc)
Definition Utils.h:50
static constexpr ID Pion
Definition PID.h:96
bool match(const std::vector< std::string > &queries, const char *pattern)
Definition dcs-ccdb.cxx:229
GLsizeiptr size
Definition glcorearb.h:659
GLuint GLuint end
Definition glcorearb.h:469
const GLdouble * v
Definition glcorearb.h:832
GLuint index
Definition glcorearb.h:781
GLboolean GLboolean GLboolean b
Definition glcorearb.h:1233
GLint GLint GLsizei GLint GLenum GLenum type
Definition glcorearb.h:275
GLsizei GLenum const void * indices
Definition glcorearb.h:400
GLbitfield flags
Definition glcorearb.h:1570
GLuint GLfloat GLfloat GLfloat GLfloat GLfloat GLfloat GLfloat t0
Definition glcorearb.h:5034
GLboolean GLboolean GLboolean GLboolean a
Definition glcorearb.h:1233
GLint GLuint mask
Definition glcorearb.h:291
constexpr float XTPCOuterRef
reference radius to propagate outer TPC track
constexpr float XTPCInnerRef
reference radius at which TPC provides the tracks
constexpr int LHCMaxBunches
constexpr int NSectors
constexpr float SectorSpanRad
float angle2Alpha(float phi)
Definition Utils.h:203
int angle2Sector(float phi)
Definition Utils.h:183
Enum< T >::Iterator begin(Enum< T >)
Definition Defs.h:173
TrackParCovF TrackParCov
Definition Track.h:33
value_T step
Definition TrackUtils.h:42
std::string getTime(uint64_t ts)
uint32_t orbit
LHC orbit.
uint16_t bc
bunch crossing ID of interaction
const o2::tpc::TrackTPC & getTPCTrack(GTrackID id) const
auto getTPCITSTrackMCLabel(GTrackID id) const
void createTracksVariadic(T creator, GTrackID::mask_t srcSel=GTrackID::getSourcesMask("all")) const
auto getTPCTrackMCLabel(GTrackID id) const
gsl::span< const unsigned char > clusterShMapTPC
externally set TPC clusters sharing map
std::unique_ptr< o2::tpc::internal::getWorkflowTPCInput_ret > inputsTPCclusters
gsl::span< const unsigned int > occupancyMapTPC
externally set TPC clusters occupancy map
static int constexpr size
Definition DataHeader.h:211
LOG(info)<< "Compressed in "<< sw.CpuTime()<< " s"
o2::InteractionRecord ir(0, 0)