diff --git a/DataFormats/Detectors/GlobalTracking/CMakeLists.txt b/DataFormats/Detectors/GlobalTracking/CMakeLists.txt index b219de73f5b47..6e7c7025adb18 100644 --- a/DataFormats/Detectors/GlobalTracking/CMakeLists.txt +++ b/DataFormats/Detectors/GlobalTracking/CMakeLists.txt @@ -45,4 +45,5 @@ o2_target_root_dictionary( DataFormatsGlobalTracking HEADERS include/DataFormatsGlobalTracking/FilteredRecoTF.h include/DataFormatsGlobalTracking/TrackTuneParams.h + include/DataFormatsGlobalTracking/TrackCosmicsExtended.h ) diff --git a/DataFormats/Detectors/GlobalTracking/include/DataFormatsGlobalTracking/TrackCosmicsExtended.h b/DataFormats/Detectors/GlobalTracking/include/DataFormatsGlobalTracking/TrackCosmicsExtended.h new file mode 100644 index 0000000000000..a6fbbcc969052 --- /dev/null +++ b/DataFormats/Detectors/GlobalTracking/include/DataFormatsGlobalTracking/TrackCosmicsExtended.h @@ -0,0 +1,127 @@ +// Copyright 2019-2026 CERN and copyright holders of ALICE O2. +// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. +// All rights not expressly granted are reserved. +// +// This software is distributed under the terms of the GNU General Public +// License v3 (GPL Version 3), copied verbatim in the file "COPYING". +// +// In applying this license CERN does not waive the privileges and immunities +// granted to it by virtue of its status as an Intergovernmental Organization +// or submit itself to any jurisdiction. + +/// \file TrackCosmicsExtended.h +/// \brief Matched cosmic track with the raw clusters of its legs and of the road around them, for offline refits + +#ifndef ALICEO2_TRACK_COSMICS_EXTENDED_H +#define ALICEO2_TRACK_COSMICS_EXTENDED_H + +#include +#include +#include +#include "ReconstructionDataFormats/TrackCosmics.h" +#include "DataFormatsTPC/TrackTPC.h" +#include "DataFormatsTPC/ClusterNative.h" +#include "SimulationDataFormat/MCCompLabel.h" + +namespace o2::dataformats +{ + +/// raw TPC cluster with its address; transformed coordinates are not stored (re-transform offline with the calibration of the TF) +struct CosmicTPCCluster { + enum Flags : uint8_t { + Attached = 0x1, ///< attached to the TPC track of this leg + Corridor = 0x2, ///< found in the road around the leg + Used = 0x4, ///< attached to some TPC track (for corridor clusters: another track, e.g. a split piece of the leg) + AbsTime = 0x8 ///< found on the other TPC side with the absolute time of the cosmic instead of the time0 of the leg: TrackCosmicsExtended::timeTOFMUS + ///< if >= 0, else the time of TrackCosmicsExtended::cosmic + }; + o2::tpc::ClusterNative cl{}; ///< raw cluster: time, pad, widths, charges, flags + uint8_t sector = 0; + uint8_t row = 0; + uint8_t flags = 0; + bool isAttached() const { return flags & Attached; } + bool isCorridor() const { return flags & Corridor; } + bool isUsed() const { return flags & Used; } + bool isAbsTime() const { return flags & AbsTime; } + ClassDefNV(CosmicTPCCluster, 1); +}; + +/// origin of the ITS / TOF / TRD hits of a cosmic +enum CosmicHitFlags : uint8_t { + HitMatched = 0x1, ///< part of the leg's matched global track + HitRoad = 0x2, ///< found by the road search around the cosmic + HitTOFFlight = 0x4 ///< TOF road hit of the top/bottom pair whose time difference matches the muon's flight (fixes the cosmic's time) +}; + +/// ITS cluster of a leg: the raw compact cluster (chip, anchor pixel, pattern ID) and its ROF; no coordinates (the dictionary and the +/// geometry are applied offline) +struct CosmicITSCluster { + uint16_t chipID = 0; + uint16_t row = 0; + uint16_t col = 0; + uint16_t pattID = 0; + int32_t pattEntry = -1; ///< start of the pattern bytes in TrackCosmicsExtended::itsPatterns (pattern not in the dictionary or group pattern), -1: none + int32_t rofBC = 0; ///< start of the cluster's ROF in BCs since the start of the TF + uint8_t leg = 0; ///< 0 bottom, 1 top + uint8_t flags = 0; ///< CosmicHitFlags + ClassDefNV(CosmicITSCluster, 1); +}; + +/// TOF cluster of a leg: matched to its track or found on the road +struct CosmicTOFCluster { + double timeRaw = 0.; ///< raw TOF time [ps] (the calibration is applied offline) + float tot = 0.f; ///< time over threshold + int32_t channel = -1; + uint8_t leg = 0; + uint8_t flags = 0; ///< CosmicHitFlags + ClassDefNV(CosmicTOFCluster, 1); +}; + +/// TRD tracklet of a leg: attached to its track or found on the road +struct CosmicTRDTracklet { + uint64_t word = 0; ///< raw Tracklet64 word + int32_t trigBC = 0; ///< BC of its trigger since the start of the TF + uint8_t layer = 0; + uint8_t leg = 0; + uint8_t flags = 0; ///< CosmicHitFlags + ClassDefNV(CosmicTRDTracklet, 1); +}; + +/// matched cosmic with everything needed for an offline refit +struct TrackCosmicsExtended { + o2::dataformats::TrackCosmics cosmic{}; ///< matcher output (time in mus; the leg references are only valid within the TF) + o2::tpc::TrackTPC tpcBottom{}; ///< TPC part of the bottom leg (default if none); z refers to its time0; its cluster + ///< references are only valid within the TF: the clusters are in clTPCBottom + o2::tpc::TrackTPC tpcTop{}; ///< TPC part of the top leg + std::vector clTPCBottom; ///< attached + corridor TPC clusters of the bottom leg + std::vector clTPCTop; ///< same for the top leg + std::vector clITS; ///< ITS clusters of the legs' matched tracks and on the road (per half and layer the closest, in the inner barrel up to 5, best first) + std::vector itsPatterns; ///< pattern bytes (row span, column span, bitmap) of the ITS clusters that need them + std::vector clTOF; ///< TOF clusters of the legs' matched tracks and on the road + std::vector trdTracklets; ///< TRD tracklets of the legs' matched tracks and on the road + float timeTOFMUS = -1.f; ///< time of the cosmic [mus since the TF start, negative before it] from its HitTOFFlight pair (valid if hasTOFTime()) + float scoreTOFPair = -1.f; ///< score of the HitTOFFlight pair: road and flight-time residuals squared in units of the cuts (< 0: none) + float scoreTOFReversed = -1.f; ///< same for the best pair in the impossible order (bottom hit first) = an accidental coincidence: + ///< QA of the flag's background; comparable to scoreTOFPair = ambiguous TOF time (< 0: none) + int duplicateOf = -1; ///< entry (in the TF) of the cosmic this one duplicates: the same muon matched twice (e.g. a leg split + ///< into two TPC tracks): >= 30 % of the smaller one's TPC clusters shared with a better one (TOF time, + ///< more attached clusters) (-1: none) + o2::track::TrackParCov polished{}; ///< refit at the TOF time (TPC-only legs, muon mass, energy loss along the flight), at the closest + ///< approach to the beam line; valid if chi2MatchPolished >= 0 + float chi2MatchPolished = -1.f; ///< chi2 of the polished bottom and top halves at the closest approach (< 0: not polished) + o2::MCCompLabel label{}; ///< MC label of the cosmic (MC only) + + bool hasTOFTime() const { return scoreTOFPair >= 0.f; } ///< a HitTOFFlight pair gave the cosmic its time timeTOFMUS + ClassDefNV(TrackCosmicsExtended, 1); +}; + +/// per-TF quantities of the TPC transformation used in the reconstruction +struct CosmicsTFInfo { + float vDrift = 0.f; ///< drift velocity of the transformation [cm/time bin] + float t0 = 0.f; ///< time offset of the transformation [time bins] + ClassDefNV(CosmicsTFInfo, 1); +}; + +} // namespace o2::dataformats + +#endif diff --git a/DataFormats/Detectors/GlobalTracking/src/DataFormatsGlobalTrackingLinkDef.h b/DataFormats/Detectors/GlobalTracking/src/DataFormatsGlobalTrackingLinkDef.h index d3519bbaad607..06d54368484d0 100644 --- a/DataFormats/Detectors/GlobalTracking/src/DataFormatsGlobalTrackingLinkDef.h +++ b/DataFormats/Detectors/GlobalTracking/src/DataFormatsGlobalTrackingLinkDef.h @@ -22,4 +22,16 @@ #pragma link C++ class o2::globaltracking::TrackTuneParams + ; #pragma link C++ class o2::conf::ConfigurableParamHelper < o2::globaltracking::TrackTuneParams> + ; +#pragma link C++ class o2::dataformats::CosmicTPCCluster + ; +#pragma link C++ class std::vector < o2::dataformats::CosmicTPCCluster> + ; +#pragma link C++ class o2::dataformats::CosmicITSCluster + ; +#pragma link C++ class std::vector < o2::dataformats::CosmicITSCluster> + ; +#pragma link C++ class o2::dataformats::CosmicTOFCluster + ; +#pragma link C++ class std::vector < o2::dataformats::CosmicTOFCluster> + ; +#pragma link C++ class o2::dataformats::CosmicTRDTracklet + ; +#pragma link C++ class std::vector < o2::dataformats::CosmicTRDTracklet> + ; +#pragma link C++ class o2::dataformats::TrackCosmicsExtended + ; +#pragma link C++ class std::vector < o2::dataformats::TrackCosmicsExtended> + ; +#pragma link C++ class o2::dataformats::CosmicsTFInfo + ; + #endif diff --git a/DataFormats/Reconstruction/include/ReconstructionDataFormats/TrackCosmics.h b/DataFormats/Reconstruction/include/ReconstructionDataFormats/TrackCosmics.h index 79a34dc585876..427995b3319fe 100644 --- a/DataFormats/Reconstruction/include/ReconstructionDataFormats/TrackCosmics.h +++ b/DataFormats/Reconstruction/include/ReconstructionDataFormats/TrackCosmics.h @@ -30,6 +30,12 @@ class TrackCosmics : public o2::track::TrackParCov using timeEst = o2::dataformats::TimeStampWithError; public: + enum Flags : uint8_t { + NearBeamBottom = 0x1, ///< bottom leg failed the standard seed DCA cuts (close to the beam line, absolutely or within its errors) and + ///< passed only the looser minSeedDCAxyTOF ones (TOF-confirmed pairs only) + NearBeamTop = 0x2 ///< the same for the top leg + }; + TrackCosmics() = default; ~TrackCosmics() = default; TrackCosmics(const TrackCosmics& src) = default; @@ -59,6 +65,10 @@ class TrackCosmics : public o2::track::TrackParCov int getNClusters() const { return mNClusters; } void setNClusters(int n) { mNClusters = n; } + uint8_t getFlags() const { return mFlags; } + void setFlags(uint8_t flags) { mFlags = flags; } + bool isNearBeam() const { return mFlags & (NearBeamBottom | NearBeamTop); } + o2::track::TrackParCov& getParamOut() { return mParamOut; } const o2::track::TrackParCov& getParamOut() const { return mParamOut; } @@ -72,8 +82,9 @@ class TrackCosmics : public o2::track::TrackParCov int mNClusters = 0; ///< total number of fitted clusters timeEst mTimeMUS; ///< time estimate in ns o2::track::TrackParCov mParamOut; ///< refitted outer parameter + uint8_t mFlags = 0; ///< Flags: NearBeamBottom, NearBeamTop - ClassDefNV(TrackCosmics, 1); + ClassDefNV(TrackCosmics, 2); }; } // namespace dataformats } // namespace o2 diff --git a/Detectors/GlobalTracking/include/GlobalTracking/MatchCosmics.h b/Detectors/GlobalTracking/include/GlobalTracking/MatchCosmics.h index 28cf1ab846412..82621288afccd 100644 --- a/Detectors/GlobalTracking/include/GlobalTracking/MatchCosmics.h +++ b/Detectors/GlobalTracking/include/GlobalTracking/MatchCosmics.h @@ -17,6 +17,8 @@ #define ALICEO2_MATCH_COSMICS #include +#include +#include #include #include "ReconstructionDataFormats/TrackCosmics.h" #include "ReconstructionDataFormats/GlobalTrackID.h" @@ -39,7 +41,8 @@ class VDriftCorrFact; namespace gpu { class TPCFastTransformPOD; -} +class GPUO2InterfaceRefit; +} // namespace gpu namespace globaltracking { @@ -67,7 +70,9 @@ class MatchCosmics RejTime, RejProp, RejChi2, - RejOther + RejOther, + RejSameHalf, + RejNoTOF ///< a leg kept by minSeedDCAxyTOF only, and no TOF flight pair }; using InfoAccessor = o2d::AbstractRefAccessor; // there is no unique structure, so the default return type is dummy (int) @@ -78,6 +83,21 @@ class MatchCosmics int id1 = MinusOne; ///< id of 2nd parnter float chi2 = -1.f; ///< matching chi2 int next = MinusOne; ///< index of eventual next record + float tCommon = 0.f; ///< common time [mus] fixed by z continuity of TPC-only legs on opposite TPC sides + float tCommonErr = -1.f; ///< its 1 sigma error [mus]; < 0: not fixed, the refit uses the centre of the time-bracket overlap and the cosmic's time error is the overlap's half-width + float tofScore = -1.f; ///< score of the top / bottom TOF hit pair matching the muon's flight (tofFlightSelection; < 0: none); a pair with one wins against pairs without + float tCommonNoTOF = 0.f; ///< tCommon before a TOF flight pair replaced it: the refit falls back to it if the refit at the TOF time fails + float tCommonErrNoTOF = -1.f; ///< tCommonErr before a TOF flight pair replaced it + }; + + struct TOFCandidate { ///< TOF cluster along the outward continuation of a TPC-only seed + int index = -1; ///< index of the TOF cluster + double timeNS = 0.; ///< its time since the start of the TF [ns] + float dy = 0.f; ///< cluster - predicted y in the frame of the cluster's sector [cm] + float dz = 0.f; ///< cluster - predicted z, the leg's z taken at its own reference time tRef [cm] + float gx = 0.f; ///< global position of the cluster [cm] + float gy = 0.f; + float gz = 0.f; }; struct TrackSeed : public o2::track::TrackParCov { @@ -86,6 +106,10 @@ class MatchCosmics int matchID = MinusOne; ///< entry (none if MinusOne) of its match in the vector of matches short vtIDMin = -1; ///< id of the 1st compatible vertex short vtIDMax = -1; ///< id of the last compatible vertex + float tRef = 0.f; ///< time [mus] the z of the parameters refers to (TPC-only: the TrackTPC time0; others: bracket centre) + int8_t tpcSide = 0; ///< TPC-only seed with clusters on one side: +1 A, -1 C (z = z(t) - side*vD*(t-tRef)); 0: z absolute + std::array xyzRef{}; ///< global position of the reference point before the propagation to the DCA (same-half veto) + bool nearBeam = false; ///< kept by the looser minSeedDCAxyTOF cuts only: usable in TOF-confirmed pairs only }; void setTPCCorrMaps(const o2::gpu::TPCFastTransformPOD* maph); void setTPCVDrift(const o2::tpc::VDriftCorrFact& v); @@ -94,6 +118,8 @@ class MatchCosmics void process(const o2::globaltracking::RecoContainer& data); void setUseMC(bool mc) { mUseMC = mc; } void setUsePVInfo(bool v) { mUsePVInfo = v; } + void setSeedSources(GTrackID::mask_t src) { mSeedSources = src; } ///< track sources used as legs (other loaded ones resolve PV contributors) + void setPVVetoSources(GTrackID::mask_t src) { mPVVetoSources = src; } ///< sources whose PV contributors veto the legs they contain void init(); void end(); @@ -129,7 +155,11 @@ class MatchCosmics private: void updateTimeDependentParams(); RejFlag checkPair(int i, int j); - void registerMatch(int i, int j, float chi2); + bool refitSeedAtTime(const TrackSeed& seed, float timeMUS, TrackSeed& out); + void registerMatch(int i, int j, float chi2, float tCommon = 0.f, float tCommonErr = -1.f, float tofScore = -1.f, float tCommonNoTOF = 0.f, float tCommonErrNoTOF = -1.f); + void prepareTOFClusters(const o2::globaltracking::RecoContainer& data); + const std::vector& getTOFCandidates(int iseed); + float findTOFFlightPair(int i, int j, float tMinMUS, float tMaxMUS, float& tofTimeMUS); void suppressMatch(int partner0, int partner1); void createSeeds(const o2::globaltracking::RecoContainer& data); bool validateMatch(int partner0); @@ -152,9 +182,24 @@ class MatchCosmics bool mFieldON = true; bool mUsePVInfo = false; bool mUseMC = true; + GTrackID::mask_t mSeedSources{GTrackID::MASK_ALL}; ///< track sources used as legs + GTrackID::mask_t mPVVetoSources{GTrackID::MASK_NONE}; ///< sources whose PV contributors veto the legs they contain float mITSROFrameLengthMUS = 0.; float mQ2PtCutoff = 1e9; + float mQ2PtCutoffOppositeSides = 1e9; const MatchCosmicsParams* mMatchParams = nullptr; + const o2::globaltracking::RecoContainer* mRecoData = nullptr; ///< inputs of the TF being processed + o2::gpu::GPUO2InterfaceRefit* mTPCRefitter = nullptr; ///< TPC refitter of the TF being processed (owned by process()) + size_t mNRefitsCommonTime = 0; ///< seeds refitted at the common time of a same-side pair in this TF + std::vector mTOFClusterOrder; ///< TOF clusters of the TF sorted in time (tofFlightSelection) + std::vector mTOFClusterTimeMUS; ///< their times since the start of the TF [mus], same order + std::vector> mSeedTOFCandidates; ///< TOF candidates per seed, filled on first use + std::vector mSeedTOFDone; ///< the TOF candidates of the seed are filled + size_t mNTOFConfirmed = 0; ///< accepted pairs with a TOF flight pair in this TF + size_t mNTOFFallbacks = 0; ///< TOF-confirmed winners refitted at their time without TOF in this TF + size_t mNSeedsNearBeam = 0; ///< TPC-only seeds kept for TOF-confirmed pairs only (minSeedDCAxyTOF) in this TF + size_t mNNearBeamConfirmed = 0; ///< accepted pairs with such a seed, confirmed by a TOF flight pair, in this TF + size_t mNSeedsPVContributors = 0; ///< seeds rejected as part of a primary-vertex contributor in this TF std::vector mCosmicTracks; std::vector mCosmicTracksLbl; diff --git a/Detectors/GlobalTracking/include/GlobalTracking/MatchCosmicsParams.h b/Detectors/GlobalTracking/include/GlobalTracking/MatchCosmicsParams.h index 4b34135d83693..6a418f216de85 100644 --- a/Detectors/GlobalTracking/include/GlobalTracking/MatchCosmicsParams.h +++ b/Detectors/GlobalTracking/include/GlobalTracking/MatchCosmicsParams.h @@ -18,6 +18,7 @@ #include "CommonUtils/ConfigurableParamHelper.h" #include "DetectorsBase/Propagator.h" #include "ReconstructionDataFormats/GlobalTrackID.h" +#include namespace o2 { @@ -29,10 +30,24 @@ struct MatchCosmicsParams : public o2::conf::ConfigurableParamHelper= this [cm] (0: no cut; rejects collision tracks in physics data) + float minSeedDCAxyNSigma = 0.f; // use only tracks with |DCA_xy| >= this * sigma(DCA_xy) (0: no cut; poorly measured collision tracks) + float minSeedDCAxyTOF = -1.f; // with tofFlightSelection: TPC-only legs failing the two cuts above but with |DCA_xy| >= this [cm] are used in TOF-confirmed pairs only (< 0: off; e.g. cosmics crossing the ITS inner barrel) + float minSeedDCAxyNSigmaTOF = 3.f; // the same for these legs in units of sigma(DCA_xy) + bool constrainTPCOnlyZ = false; // TPC-only legs: test z at a common time (same side, or a leg with known time), else require the time implied by z continuity in both brackets + bool vetoSameHalf = false; // reject pairs whose two legs lie on the same side of the closest approach (two pieces of one leg) + bool refitSameSideAtCommonTime = false; // TPC-only legs on the same side: compare them refitted at the centre of their brackets' overlap instead of at their own time0s + bool tofFlightSelection = false; // needs TOF clusters: accepted pairs of TPC-only legs pointing to a top / bottom TOF hit pair with the muon's flight time win the selection, refit at that time (if that fails, at their time without TOF) + float tofRoad = 5.f; // half-width [cm] in y and z of the road at the TOF around the outward continuation of a TPC-only leg + float tofFlightTolerance = 2.f; // max. deviation [ns] of the top / bottom TOF time difference from the flight time along the helix + float tofTimeError = 0.1f; // error [mus] of the TOF time of a confirmed cosmic, for its refit and time window (covers TPC vs TOF offsets) float nSigmaTError = 4.f; // number of sigmas on track time error for matching (except for TPC which provides an interval) float tpcExtraZError2 = 1.f; // extra error^2 on the TPC-only track Z coordinate float fiducialRIP = 1.0f; // consider track having |Y@x=0|< this as passing DCA cut (if requested) @@ -44,6 +59,12 @@ struct MatchCosmicsParams : public o2::conf::ConfigurableParamHelperpropagateToDCABxByBz(v, trc, mMatchParams->maxStep, mMatchParams->matCorr)) { + // a cosmic muon flies downward: along a leg whose outward direction points up (top leg) the inward propagation follows the flight, + // so the energy is lost; along the bottom leg it goes back in the flight, so the energy is gained + std::array momentum{}; + trc.getPxPyPzGlo(momentum); + const int eLossSign = momentum[1] > 0.f ? ELossLoss : ELossGain; + if (!prop->propagateToDCABxByBz(v, trc, mMatchParams->maxStep, mMatchParams->matCorr, nullptr, nullptr, eLossSign)) { trc.matchID = Reject; // reject track continue; } - if (mMatchParams->dcaCutChi2[trc.origID.getSource()] > 0.f && mUsePVInfo && (std::abs(trc.getY()) < mMatchParams->fiducialRIP && std::abs(trc.getZ()) < mMatchParams->fiducialZIP)) { + if (std::abs(trc.getY()) < mMatchParams->minSeedDCAxy || trc.getY() * trc.getY() < mMatchParams->minSeedDCAxyNSigma * mMatchParams->minSeedDCAxyNSigma * trc.getSigmaY2()) { + // passes close to the beam line, absolutely or within its errors: indistinguishable from collision tracks, unless a TOF flight pair + // confirms the cosmic (TPC-only legs above the looser minSeedDCAxyTOF cuts, used in TOF-confirmed pairs only) + const bool nearBeam = mMatchParams->tofFlightSelection && mMatchParams->minSeedDCAxyTOF >= 0.f && trc.origID.getSource() == GTrackID::TPC && + std::abs(trc.getY()) >= mMatchParams->minSeedDCAxyTOF && + trc.getY() * trc.getY() >= mMatchParams->minSeedDCAxyNSigmaTOF * mMatchParams->minSeedDCAxyNSigmaTOF * trc.getSigmaY2(); + if (!nearBeam) { + trc.matchID = Reject; + continue; + } + trc.nearBeam = true; + mNSeedsNearBeam++; + } + if (mMatchParams->dcaCutChi2[trc.origID.getSource()] > 0.f && mUsePVInfo && trc.vtIDMin >= 0 && (std::abs(trc.getY()) < mMatchParams->fiducialRIP && std::abs(trc.getZ()) < mMatchParams->fiducialZIP)) { // do the propagation only if we are in the fiducial IP range. - for (int iv = trc.vtIDMin; iv < trc.vtIDMax; iv++) { + for (int iv = trc.vtIDMin; iv <= trc.vtIDMax; iv++) { // vtIDMax is the last compatible vertex (inclusive); vtIDMin < 0: no compatible vertex const auto& pv = data.getPrimaryVertex(iv); o2::track::TrackParCov trcatPV(trc); o2::dataformats::DCA dca; @@ -92,6 +121,23 @@ void MatchCosmics::process(const o2::globaltracking::RecoContainer& data) } } } + // TPC refitter of this TF: same-side TPC-only legs compared at a common time (checkPair) and the refit of the winners + std::unique_ptr tpcRefitter; + if (data.inputsTPCclusters) { + tpcRefitter = std::make_unique(&data.inputsTPCclusters->clusterIndex, mTPCCorrMaps, mBz, data.getTPCTracksClusterRefs().data(), 0, + data.clusterShMapTPC.data(), data.occupancyMapTPC.data(), data.occupancyMapTPC.size(), nullptr, + o2::base::Propagator::Instance()); + } + mTPCRefitter = tpcRefitter.get(); + mRecoData = &data; + mNRefitsCommonTime = 0; + mNTOFConfirmed = 0; + mNTOFFallbacks = 0; + mNNearBeamConfirmed = 0; + if (mMatchParams->tofFlightSelection) { + prepareTOFClusters(data); + } + // sort in time bracket lower edge, putting rejected tracks in the end std::vector sortID(ntr); std::iota(sortID.begin(), sortID.end(), 0); @@ -112,9 +158,26 @@ void MatchCosmics::process(const o2::globaltracking::RecoContainer& data) } } } + if (mMatchParams->tofFlightSelection) { + LOGP(info, "{} accepted pairs confirmed by a TOF flight pair", mNTOFConfirmed); + } + if (mNSeedsPVContributors) { + LOGP(info, "{} seeds rejected as (part of) a primary-vertex contributor", mNSeedsPVContributors); + } + if (mNSeedsNearBeam) { + LOGP(info, "{} seeds near the beam line used in TOF-confirmed pairs only, {} such pairs confirmed", mNSeedsNearBeam, mNNearBeamConfirmed); + } selectWinners(); refitWinners(data); + if (mNRefitsCommonTime) { + LOGP(info, "{} seeds refitted at the common time of same-side pairs", mNRefitsCommonTime); + } + if (mNTOFFallbacks) { + LOGP(info, "{} TOF-confirmed winners failed the refit at the TOF time and were refitted at their time without TOF", mNTOFFallbacks); + } + mTPCRefitter = nullptr; + mRecoData = nullptr; mTFCount++; } @@ -125,16 +188,7 @@ void MatchCosmics::refitWinners(const o2::globaltracking::RecoContainer& data) LOG(info) << "Refitting " << mWinners.size() << " winner matches"; int count = 0; auto tpcTBinMUSInv = 1. / mTPCTBinMUS; - const auto& tpcClusRefs = data.getTPCTracksClusterRefs(); - const auto& tpcClusShMap = data.clusterShMapTPC; - const auto& tpcClusOccMap = data.occupancyMapTPC; - std::unique_ptr tpcRefitter; - if (data.inputsTPCclusters) { - tpcRefitter = std::make_unique(&data.inputsTPCclusters->clusterIndex, - mTPCCorrMaps, mBz, - tpcClusRefs.data(), 0, tpcClusShMap.data(), - tpcClusOccMap.data(), tpcClusOccMap.size(), nullptr, o2::base::Propagator::Instance()); - } + auto* tpcRefitter = mTPCRefitter; // created in process() const auto& itsClusters = prepareITSClusters(data); // RS FIXME: this is probably a temporary solution, since ITS tracking over boundaries will likely change the TrackITS format @@ -149,7 +203,7 @@ void MatchCosmics::refitWinners(const o2::globaltracking::RecoContainer& data) } } - auto refitITSTrack = [this, &data, &itsTracksROF, &itsClusters](o2::track::TrackParCov& trFit, GTrackID gidx, float& chi2, bool inward = false) { + auto refitITSTrack = [this, &data, &itsTracksROF, &itsClusters](o2::track::TrackParCov& trFit, GTrackID gidx, float& chi2, bool inward, int eLossSign) { const auto& itsTrOrig = data.getITSTrack(gidx); int nclRefit = 0, ncl = itsTrOrig.getNumberOfClusters(), rof = itsTracksROF[gidx.getIndex()]; const auto& itsTrackClusRefs = data.getITSTracksClusterRefs(); @@ -165,7 +219,7 @@ void MatchCosmics::refitWinners(const o2::globaltracking::RecoContainer& data) for (int icl = from; icl != to; icl += step) { // ITS clusters are referred in layer decreasing order const auto& clus = itsClusters[itsTrackClusRefs[clEntry + icl]]; float alpha = geomITS->getSensorRefAlpha(clus.getSensorID()), x = clus.getX(); - if (!trFit.rotate(alpha) || !propagator->propagateToX(trFit, x, propagator->getNominalBz(), this->mMatchParams->maxSnp, this->mMatchParams->maxStep, this->mMatchParams->matCorr)) { + if (!trFit.rotate(alpha) || !propagator->propagateToX(trFit, x, propagator->getNominalBz(), this->mMatchParams->maxSnp, this->mMatchParams->maxStep, this->mMatchParams->matCorr, nullptr, eLossSign)) { break; } chi2 += trFit.getPredictedChi2(clus); @@ -177,12 +231,15 @@ void MatchCosmics::refitWinners(const o2::globaltracking::RecoContainer& data) return nclRefit == ncl ? ncl : -1; }; - for (auto winRID : mWinners) { + // refit the legs of winner winRID at the time t0 [mus] (error dt) and add the cosmic; false: a refit, propagation or cut failed + auto refitWinner = [&](int winRID, float t0, float dt) { const auto& rec = mRecords[winRID]; int poolEntryID[2] = {rec.id0, rec.id1}; - const o2::track::TrackParCov outerLegs[2] = {data.getTrackParamOut(mSeeds[rec.id0].origID), data.getTrackParamOut(mSeeds[rec.id1].origID)}; + o2::track::TrackParCov outerLegs[2] = {data.getTrackParamOut(mSeeds[rec.id0].origID), data.getTrackParamOut(mSeeds[rec.id1].origID)}; + for (auto& leg : outerLegs) { + leg.setPID(o2::track::PID::Muon, true); // as the seeds + } auto tOverlap = mSeeds[rec.id0].tBracket.getOverlap(mSeeds[rec.id1].tBracket); - float t0 = tOverlap.mean(), dt = tOverlap.delta() * 0.5; auto pnt0 = outerLegs[0].getXYZGlo(), pnt1 = outerLegs[1].getXYZGlo(); int btm = 0, top = 1; // we fit topward from bottom @@ -209,10 +266,10 @@ void MatchCosmics::refitWinners(const o2::globaltracking::RecoContainer& data) if (!mFieldON) { trCosm.setQ2Pt(-o2::track::kMostProbablePt); } - int retVal = tpcRefitter->RefitTrackAsTrackParCov(trCosm, tpcTrOrig.getClusterRef(), t0 * tpcTBinMUSInv, &chi2, false, false); // inward refit, reset + int retVal = tpcRefitter->RefitTrackAsTrackParCov(trCosm, tpcTrOrig.getClusterRef(), t0 * tpcTBinMUSInv, &chi2, false, false, ELossGain); // inward refit, reset if (retVal < 0) { // refit failed LOG(debug) << "Inward refit of btm TPC track failed."; - continue; + return false; } nclTot += retVal; LOG(debug) << "chi2 after btm TPC refit with " << retVal << " clusters : " << chi2 << " orig.chi2 was " << tpcTrOrig.getChi2(); @@ -235,9 +292,9 @@ void MatchCosmics::refitWinners(const o2::globaltracking::RecoContainer& data) } trCosm.invert(); if (!trCosm.rotate(mSeeds[poolEntryID[top]].getAlpha()) || - !o2::base::Propagator::Instance()->PropagateToXBxByBz(trCosm, mSeeds[poolEntryID[top]].getX(), mMatchParams->maxSnp, mMatchParams->maxStep, mMatchParams->matCorr)) { + !o2::base::Propagator::Instance()->PropagateToXBxByBz(trCosm, mSeeds[poolEntryID[top]].getX(), mMatchParams->maxSnp, mMatchParams->maxStep, mMatchParams->matCorr, nullptr, ELossGain)) { LOG(debug) << "Rotation/propagation of btm-track to top-track frame failed."; - continue; + return false; } // save bottom parameter at merging point auto trCosmBtm = trCosm; @@ -248,9 +305,9 @@ void MatchCosmics::refitWinners(const o2::globaltracking::RecoContainer& data) // is there ITS sub-track? if (gidxListTop[GTrackID::ITS].isIndexSet()) { - auto nclfit = refitITSTrack(trCosm, gidxListTop[GTrackID::ITS], chi2, false); + auto nclfit = refitITSTrack(trCosm, gidxListTop[GTrackID::ITS], chi2, false, ELossGain); if (nclfit < 0) { - continue; + return false; } LOG(debug) << "chi2 after top ITS refit with " << nclfit << " clusters : " << chi2 << " orig.chi2 was " << data.getITSTrack(gidxListTop[GTrackID::ITS]).getChi2(); nclTot += nclfit; @@ -261,16 +318,16 @@ void MatchCosmics::refitWinners(const o2::globaltracking::RecoContainer& data) if (trCosm.getX() * trCosm.getX() + trCosm.getY() * trCosm.getY() <= o2::constants::geom::XTPCInnerRef * o2::constants::geom::XTPCInnerRef) { float xtogo = 0; if (!trCosm.getXatLabR(o2::constants::geom::XTPCInnerRef, xtogo, mBz, o2::track::DirOutward) || - !o2::base::Propagator::Instance()->PropagateToXBxByBz(trCosm, xtogo, mMatchParams->maxSnp, mMatchParams->maxStep, mMatchParams->matCorr)) { + !o2::base::Propagator::Instance()->PropagateToXBxByBz(trCosm, xtogo, mMatchParams->maxSnp, mMatchParams->maxStep, mMatchParams->matCorr, nullptr, ELossGain)) { LOG(debug) << "Propagation to inner TPC boundary X=" << xtogo << " failed"; - continue; + return false; } } const auto& tpcTrOrig = data.getTPCTrack(gidxListTop[GTrackID::TPC]); - int retVal = tpcRefitter->RefitTrackAsTrackParCov(trCosm, tpcTrOrig.getClusterRef(), t0 * tpcTBinMUSInv, &chi2, true, false); // outward refit, no reset + int retVal = tpcRefitter->RefitTrackAsTrackParCov(trCosm, tpcTrOrig.getClusterRef(), t0 * tpcTBinMUSInv, &chi2, true, false, ELossGain); // outward refit, no reset if (retVal < 0) { // refit failed LOG(debug) << "Outward refit of top TPC track failed."; - continue; + return false; } // outward refit in TPC LOG(debug) << "chi2 after top TPC refit with " << retVal << " clusters : " << chi2 << " orig.chi2 was " << tpcTrOrig.getChi2(); nclTot += retVal; @@ -281,40 +338,74 @@ void MatchCosmics::refitWinners(const o2::globaltracking::RecoContainer& data) auto trCosmTop = outerLegs[top]; if (gidxListTop[GTrackID::TPC].isIndexSet()) { // inward refit in TPC const auto& tpcTrOrig = data.getTPCTrack(gidxListTop[GTrackID::TPC]); - int retVal = tpcRefitter->RefitTrackAsTrackParCov(trCosmTop, tpcTrOrig.getClusterRef(), t0 * tpcTBinMUSInv, &chi2Dummy, false, true); // inward refit, reset + int retVal = tpcRefitter->RefitTrackAsTrackParCov(trCosmTop, tpcTrOrig.getClusterRef(), t0 * tpcTBinMUSInv, &chi2Dummy, false, true, ELossLoss); // inward refit, reset if (retVal < 0) { // refit failed - LOG(debug) << "Outward refit of top TPC track failed."; - continue; + LOG(debug) << "Inward refit of top TPC track failed."; + return false; } // inward refit in TPC } // is there ITS sub-track ? if (gidxListTop[GTrackID::ITS].isIndexSet()) { - auto nclfit = refitITSTrack(trCosmTop, gidxListTop[GTrackID::ITS], chi2Dummy, true); + auto nclfit = refitITSTrack(trCosmTop, gidxListTop[GTrackID::ITS], chi2Dummy, true, ELossLoss); if (nclfit < 0) { - continue; + return false; } nclTot += nclfit; } // ITS refit // propagate to bottom param if (!trCosmTop.rotate(trCosmBtm.getAlpha()) || - !o2::base::Propagator::Instance()->PropagateToXBxByBz(trCosmTop, trCosmBtm.getX(), mMatchParams->maxSnp, mMatchParams->maxStep, mMatchParams->matCorr)) { + !o2::base::Propagator::Instance()->PropagateToXBxByBz(trCosmTop, trCosmBtm.getX(), mMatchParams->maxSnp, mMatchParams->maxStep, mMatchParams->matCorr, nullptr, ELossLoss)) { LOG(debug) << "Rotation/propagation of top-track to bottom-track frame failed."; - continue; + return false; } // calculate weighted average of 2 legs and chi2 o2::track::TrackParCov::MatrixDSym5 cov5; float chi2Match = trCosmBtm.getPredictedChi2(trCosmTop, cov5); + if (mMatchParams->maxChi2Match >= 0.f && chi2Match > mMatchParams->maxChi2Match) { + LOG(debug) << "Top/Bottom refitted legs disagree, chi2Match " << chi2Match; + return false; + } if (!trCosmBtm.update(trCosmTop, cov5)) { LOG(debug) << "Top/Bottom update failed"; - continue; + return false; + } + // TPC-only legs on opposite sides: the legs' pT is required in checkPair, the refitted cosmic's here + if (mSeeds[rec.id0].tpcSide * mSeeds[rec.id1].tpcSide < 0 && std::abs(trCosmBtm.getQ2Pt()) > mQ2PtCutoffOppositeSides) { + LOG(debug) << "Cosmic with legs on opposite TPC sides below minPtOppositeSides"; + return false; } // create final track - mCosmicTracks.emplace_back(mSeeds[poolEntryID[btm]].origID, mSeeds[poolEntryID[top]].origID, trCosmBtm, trCosmTop, chi2, chi2Match, nclTot, t0, dt); + auto& cosmic = mCosmicTracks.emplace_back(mSeeds[poolEntryID[btm]].origID, mSeeds[poolEntryID[top]].origID, trCosmBtm, trCosmTop, chi2, chi2Match, nclTot, t0, dt); + uint8_t flags = 0; + if (mSeeds[poolEntryID[btm]].nearBeam) { + flags |= o2d::TrackCosmics::NearBeamBottom; + } + if (mSeeds[poolEntryID[top]].nearBeam) { + flags |= o2d::TrackCosmics::NearBeamTop; + } + cosmic.setFlags(flags); if (mUseMC) { o2::MCCompLabel lbl[2] = {data.getTrackMCLabel(mSeeds[poolEntryID[btm]].origID), data.getTrackMCLabel(mSeeds[poolEntryID[top]].origID)}; auto& tlb = mCosmicTracksLbl.emplace_back((nclBtm > nclTot - nclBtm ? lbl[0] : lbl[1])); tlb.setFakeFlag(lbl[0] != lbl[1]); } + return true; + }; + for (auto winRID : mWinners) { + const auto& rec = mRecords[winRID]; + // refit at the common time if one is fixed (z continuity of TPC-only legs on opposite sides, TOF flight pair), else at the centre of + // the overlap of the legs' time brackets + auto refitAt = [&](float tCommon, float tCommonErr) { + if (tCommonErr >= 0.f) { + return refitWinner(winRID, tCommon, tCommonErr); + } + auto tOverlap = mSeeds[rec.id0].tBracket.getOverlap(mSeeds[rec.id1].tBracket); + return refitWinner(winRID, tOverlap.mean(), tOverlap.delta() * 0.5f); + }; + // a TOF-confirmed winner whose refit at the TOF time fails is refitted at the time it has without its TOF flight pair + if (!refitAt(rec.tCommon, rec.tCommonErr) && rec.tofScore >= 0.f && refitAt(rec.tCommonNoTOF, rec.tCommonErrNoTOF)) { + mNTOFFallbacks++; + } } LOG(info) << "Validated " << mCosmicTracks.size() << " top-bottom tracks in TF# " << mTFCount; } @@ -424,10 +515,22 @@ MatchCosmics::RejFlag MatchCosmics::checkPair(int i, int j) return (rej = RejTime); // since the brackets are sorted in tmin, all following tbj will also exceed tbi } float chi2 = 1.e9f; + float tCommon = 0.f; // time fixed by z continuity of TPC-only legs on opposite sides, used by the refit + float tCommonErr = -1.f; // its error (< 0: not fixed) + float tofScore = -1.f; // score of the TOF flight pair of an accepted pair (tofFlightSelection; < 0: none) + TrackSeed seed0Common; // same-side TPC-only legs refitted at a common time (refitSameSideAtCommonTime) + TrackSeed seed1Common; + bool commonTime = false; // check // 1) crude check on tgl and q/pt (if B!=0). Note: back-to-back tracks will have mutually params (see TrackPar::invertParam) while (1) { + // TPC-only legs on opposite sides: their z continuity defines the time, so z does not reject random pairs of collision tracks; require + // the pT of a cosmic for both legs already here, so that such a pair cannot win against the true partner of one of its legs + if (seed0.tpcSide * seed1.tpcSide < 0 && std::max(std::abs(seed0.getQ2Pt()), std::abs(seed1.getQ2Pt())) > mQ2PtCutoffOppositeSides) { + rej = RejQ2Pt; + break; + } auto dTgl = seed0.getTgl() + seed1.getTgl(); if (dTgl * dTgl > (mMatchParams->systSigma2[o2::track::kTgl] + seed0.getSigmaTgl2() + seed1.getSigmaTgl2()) * mMatchParams->crudeNSigma2Cut[o2::track::kTgl]) { rej = RejTgl; @@ -440,33 +543,101 @@ MatchCosmics::RejFlag MatchCosmics::checkPair(int i, int j) break; } } - o2::track::TrackParCov seed1Inv = seed1; + if (mMatchParams->vetoSameHalf) { + // a cosmic has its two legs on opposite sides of its closest approach to the beam line: project the legs' reference points (before + // the propagation to the DCA) on the transverse direction of seed0 at its DCA; pieces of one leg are on the same side. Transverse + // only: the z of a TPC-only track refers to its own time0, so z differences between the legs are meaningless. Skipped for + // reference points close to the DCA (e.g. ITS-containing tracks), where the sign is undefined. + std::array pca{}; + seed0.getXYZGlo(pca); + const float phi = seed0.getAlpha() + std::asin(seed0.getSnp()); + const float dir[2] = {std::cos(phi), std::sin(phi)}; + float proj0 = 0.f; + float proj1 = 0.f; + for (int k = 0; k < 2; k++) { + proj0 += (seed0.xyzRef[k] - pca[k]) * dir[k]; + proj1 += (seed1.xyzRef[k] - pca[k]) * dir[k]; + } + constexpr float MinDist = 20.f; // cm + if (std::abs(proj0) > MinDist && std::abs(proj1) > MinDist && proj0 * proj1 > 0.f) { + rej = RejSameHalf; + break; + } + } + // TPC-only legs on the same side: the tracker transformed each leg's clusters with its own time0, a guess (the two legs of a cosmic + // get guesses tens of mus apart), so the legs were distortion-corrected at different z and disagree although they are one track; + // compare them refitted at one common time, the centre of their brackets' overlap (the time the refit of the winners uses) + if (mMatchParams->refitSameSideAtCommonTime && seed0.tpcSide != 0 && seed0.tpcSide == seed1.tpcSide) { // tpcSide != 0: TPC-only one-side legs + const float tPair = seed0.tBracket.getOverlap(seed1.tBracket).mean(); + commonTime = refitSeedAtTime(seed0, tPair, seed0Common) && refitSeedAtTime(seed1, tPair, seed1Common); + } + const TrackSeed& leg0 = commonTime ? seed0Common : seed0; + const TrackSeed& leg1 = commonTime ? seed1Common : seed1; + o2::track::TrackParCov seed1Inv = leg1; seed1Inv.invert(); for (int i = 0; i < o2::track::kNParams; i++) { // add systematic error seed1Inv.updateCov(mMatchParams->systSigma2[i], o2::track::DiagMap[i]); } - if (!seed1Inv.rotate(seed0.getAlpha()) || - !o2::base::Propagator::Instance()->PropagateToXBxByBz(seed1Inv, seed0.getX(), mMatchParams->maxSnp, mMatchParams->maxStep, mMatchParams->matCorr)) { + if (!seed1Inv.rotate(leg0.getAlpha()) || + !o2::base::Propagator::Instance()->PropagateToXBxByBz(seed1Inv, leg0.getX(), mMatchParams->maxSnp, mMatchParams->maxStep, mMatchParams->matCorr)) { rej = RejProp; break; } - auto dSnp = seed0.getSnp() - seed1Inv.getSnp(); - if (dSnp * dSnp > (seed0.getSigmaSnp2() + seed1Inv.getSigmaSnp2()) * mMatchParams->crudeNSigma2Cut[o2::track::kSnp]) { + auto dSnp = leg0.getSnp() - seed1Inv.getSnp(); + if (dSnp * dSnp > (leg0.getSigmaSnp2() + seed1Inv.getSigmaSnp2()) * mMatchParams->crudeNSigma2Cut[o2::track::kSnp]) { rej = RejSnp; break; } - auto dY = seed0.getY() - seed1Inv.getY(); - if (dY * dY > (seed0.getSigmaY2() + seed1Inv.getSigmaY2()) * mMatchParams->crudeNSigma2Cut[o2::track::kY]) { + auto dY = leg0.getY() - seed1Inv.getY(); + if (dY * dY > (leg0.getSigmaY2() + seed1Inv.getSigmaY2()) * mMatchParams->crudeNSigma2Cut[o2::track::kY]) { rej = RejY; break; } - bool ignoreZ = seed0.origID.getSource() == o2d::GlobalTrackID::TPC || seed1.origID.getSource() == o2d::GlobalTrackID::TPC; - // RSTODO this is simplification: one should constraint the TPC-only track Z by the time of other candidate (at least their difference of both tracks are TPC only). - // If the shift is large, eventually the tracks need to be refitted. - if (!ignoreZ) { // cut on Z makes no sense for TPC only tracks - auto dZ = seed0.getZ() - seed1Inv.getZ(); - if (dZ * dZ > (seed0.getSigmaZ2() + seed1Inv.getSigmaZ2()) * mMatchParams->crudeNSigma2Cut[o2::track::kZ]) { + bool ignoreZ = leg0.origID.getSource() == o2d::GlobalTrackID::TPC || leg1.origID.getSource() == o2d::GlobalTrackID::TPC; + if (ignoreZ && mMatchParams->constrainTPCOnlyZ) { + // a TPC-only track with clusters on one side has z relative to its time0: z(t) = z + side * vD * (t - tRef); a CE-crossing or + // non-TPC-only track has an absolute z (side 0). Bring both legs to a common time where possible and test z; for legs on opposite + // sides, z continuity fixes the common time, which must lie in both time brackets. + const int side0 = leg0.tpcSide; + const int side1 = leg1.tpcSide; + const float sigZ2 = (leg0.getSigmaZ2() + seed1Inv.getSigmaZ2()) * mMatchParams->crudeNSigma2Cut[o2::track::kZ]; + if (side0 == 0 || side1 == 0 || side0 == side1) { + float dZ = leg0.getZ() - seed1Inv.getZ(); + float dZTimeTol = 0.f; // the time of a non-TPC absolute leg is only known within its bracket (tRef is the bracket centre) + if (side0 != 0 && side1 != 0) { // same side: the z offset is fixed by the difference of the reference times + dZ -= side0 * mTPCVDrift * (leg0.tRef - leg1.tRef); + } else if (side1 != 0) { // leg0 absolute: move leg1 to the time of leg0 + dZ -= side1 * mTPCVDrift * (leg0.tRef - leg1.tRef); + if (leg0.origID.getSource() != o2d::GlobalTrackID::TPC) { + dZTimeTol = 0.5f * mTPCVDrift * leg0.tBracket.delta(); + } + } else if (side0 != 0) { // leg1 absolute: move leg0 to the time of leg1 + dZ += side0 * mTPCVDrift * (leg1.tRef - leg0.tRef); + if (leg1.origID.getSource() != o2d::GlobalTrackID::TPC) { + dZTimeTol = 0.5f * mTPCVDrift * leg1.tBracket.delta(); + } + } + const float dZTol = std::sqrt(sigZ2) + dZTimeTol; + if (dZ * dZ > dZTol * dZTol) { + rej = RejZ; + break; + } + } else { // opposite sides + const float t = 0.5f * (side0 * (seed1Inv.getZ() - leg0.getZ()) / mTPCVDrift + leg0.tRef + leg1.tRef); + const float tTol = std::sqrt(sigZ2) / (2.f * mTPCVDrift); + if (t < std::max(leg0.tBracket.getMin(), leg1.tBracket.getMin()) - tTol || t > std::min(leg0.tBracket.getMax(), leg1.tBracket.getMax()) + tTol) { + rej = RejZ; + break; + } + tCommon = t; + tCommonErr = std::sqrt(leg0.getSigmaZ2() + seed1Inv.getSigmaZ2()) / (2.f * mTPCVDrift); + } + } + // the z of a TPC-only leg refers to its own time0: it is tested above at a common time with constrainTPCOnlyZ, otherwise ignored + if (!ignoreZ) { // both legs have an absolute z + auto dZ = leg0.getZ() - seed1Inv.getZ(); + if (dZ * dZ > (leg0.getSigmaZ2() + seed1Inv.getSigmaZ2()) * mMatchParams->crudeNSigma2Cut[o2::track::kZ]) { rej = RejZ; break; } @@ -478,26 +649,51 @@ MatchCosmics::RejFlag MatchCosmics::checkPair(int i, int j) seed1Inv.setCov(0., o2::track::CovarMap[o2::track::kZ][o2::track::kQ2Pt]); } // calculate chi2 (expensive) - chi2 = seed0.getPredictedChi2(seed1Inv); + chi2 = leg0.getPredictedChi2(seed1Inv); if (chi2 > mMatchParams->crudeChi2Cut) { rej = RejChi2; break; } rej = Accept; - registerMatch(i, j, chi2); - registerMatch(j, i, chi2); // the reverse reference can be also done in a separate loop + const float tCommonNoTOF = tCommon; // the pair's time without a TOF flight pair: the fallback of the refit at the TOF time + const float tCommonErrNoTOF = tCommonErr; + if (mMatchParams->tofFlightSelection) { // a top / bottom TOF hit pair with the muon's flight time confirms the pair and gives its time + const bool timeFixed = tCommonErr >= 0.f; + const auto overlap = seed0.tBracket.getOverlap(seed1.tBracket); + const float tMin = timeFixed ? tCommon - mMatchParams->nSigmaTError * tCommonErr : overlap.getMin(); + const float tMax = timeFixed ? tCommon + mMatchParams->nSigmaTError * tCommonErr : overlap.getMax(); + float tofTimeMUS = 0.f; + tofScore = findTOFFlightPair(i, j, tMin, tMax, tofTimeMUS); + if (tofScore >= 0.f) { + tCommon = tofTimeMUS; + tCommonErr = mMatchParams->tofTimeError; + mNTOFConfirmed++; + } + } + if (seed0.nearBeam || seed1.nearBeam) { // a leg near the beam line (minSeedDCAxyTOF): TOF-confirmed pairs only + if (tofScore < 0.f) { + rej = RejNoTOF; + break; + } + mNNearBeamConfirmed++; + } + registerMatch(i, j, chi2, tCommon, tCommonErr, tofScore, tCommonNoTOF, tCommonErrNoTOF); + registerMatch(j, i, chi2, tCommon, tCommonErr, tofScore, tCommonNoTOF, tCommonErrNoTOF); // the reverse reference can be also done in a separate loop LOG(debug) << "Chi2 = " << chi2 << " NMatches " << mRecords.size(); break; } #ifdef _ALLOW_DEBUG_TREES_ if (mDBGOut && ((rej == Accept && isDebugFlag(MatchTreeAccOnly)) || isDebugFlag(MatchTreeAll))) { - auto seed1I = seed1; + const TrackSeed& dbgLeg0 = commonTime ? seed0Common : seed0; // the legs compared (refitted at the common time if they were) + auto seed1I = commonTime ? seed1Common : seed1; seed1I.invert(); - if (seed1I.rotate(seed0.getAlpha()) && o2::base::Propagator::Instance()->PropagateToXBxByBz(seed1I, seed0.getX(), mMatchParams->maxSnp, mMatchParams->maxStep, mMatchParams->matCorr)) { + if (seed1I.rotate(dbgLeg0.getAlpha()) && o2::base::Propagator::Instance()->PropagateToXBxByBz(seed1I, dbgLeg0.getX(), mMatchParams->maxSnp, mMatchParams->maxStep, mMatchParams->matCorr)) { int rejI = int(rej); + int commonTimeI = commonTime; (*mDBGOut) << "match" - << "tf=" << mTFCount << "seed0=" << seed0 << "seed1=" << seed1I << "chi2Match=" << chi2 << "rej=" << rejI << "\n"; + << "tf=" << mTFCount << "seed0=" << dbgLeg0 << "seed1=" << seed1I << "chi2Match=" << chi2 << "rej=" << rejI << "commonTime=" << commonTimeI + << "side0=" << int(seed0.tpcSide) << "side1=" << int(seed1.tpcSide) << "tCommon=" << tCommon << "tCommonErr=" << tCommonErr << "tofScore=" << tofScore << "\n"; } } #endif @@ -506,15 +702,193 @@ MatchCosmics::RejFlag MatchCosmics::checkPair(int i, int j) } //________________________________________________________ -void MatchCosmics::registerMatch(int i, int j, float chi2) +bool MatchCosmics::refitSeedAtTime(const TrackSeed& seed, float timeMUS, TrackSeed& out) +{ + // TPC-only seed refitted with its clusters transformed at the time timeMUS, then treated as the seeds in createSeeds and process() (muon, + // extra z error, propagation to the DCA to the beam line); energy loss along the muon's flight in the refit and the propagation, as in + // refitWinners. Its z then refers to timeMUS + if (!mTPCRefitter || !mRecoData) { + return false; + } + const auto& tpcTrack = mRecoData->getTPCTrack(seed.origID); + o2::track::TrackParCov trk = tpcTrack.getParamOut(); + trk.setPID(o2::track::PID::Muon, true); + trk.resetCovariance(); + // the muon flies downward: inward along a leg pointing up (top leg) is along its flight (loss), along the bottom leg against it (gain) + std::array momentum{}; + trk.getPxPyPzGlo(momentum); + const int eLossSign = momentum[1] > 0.f ? ELossLoss : ELossGain; + if (mTPCRefitter->RefitTrackAsTrackParCov(trk, tpcTrack.getClusterRef(), timeMUS / mTPCTBinMUS, nullptr, false, true, eLossSign) < 0) { + return false; + } + trk.setCov(mMatchParams->tpcExtraZError2 + trk.getSigmaZ2(), o2::track::kSigZ2); + const o2::dataformats::VertexBase v; + if (!o2::base::Propagator::Instance()->propagateToDCABxByBz(v, trk, mMatchParams->maxStep, mMatchParams->matCorr, nullptr, nullptr, eLossSign)) { + return false; + } + out = seed; + static_cast(out) = trk; + out.tRef = timeMUS; + mNRefitsCommonTime++; + return true; +} + +//________________________________________________________ +void MatchCosmics::prepareTOFClusters(const o2::globaltracking::RecoContainer& data) +{ + // TOF clusters of the TF sorted in time; the TOF candidates of a seed are searched on its first use + const auto clusters = data.getTOFClusters(); + mTOFClusterOrder.resize(clusters.size()); + std::iota(mTOFClusterOrder.begin(), mTOFClusterOrder.end(), 0); + std::sort(mTOFClusterOrder.begin(), mTOFClusterOrder.end(), [&clusters](int a, int b) { return clusters[a].getTime() < clusters[b].getTime(); }); + mTOFClusterTimeMUS.resize(clusters.size()); + for (size_t k = 0; k < clusters.size(); k++) { + mTOFClusterTimeMUS[k] = clusters[mTOFClusterOrder[k]].getTime() * 1e-6; // [ps] since the start of the TF + } + mSeedTOFCandidates.clear(); + mSeedTOFCandidates.resize(mSeeds.size()); + mSeedTOFDone.assign(mSeeds.size(), false); +} + +//________________________________________________________ +const std::vector& MatchCosmics::getTOFCandidates(int iseed) +{ + // TOF clusters along the outward continuation of a TPC-only seed (from its outer parameters, in the sector it points to and its two + // neighbours) within its time bracket: |dy| < tofRoad, and dz within tofRoad of the drift of a one-side leg's z over the bracket + auto& candidates = mSeedTOFCandidates[iseed]; + if (mSeedTOFDone[iseed]) { + return candidates; + } + mSeedTOFDone[iseed] = true; + const auto& seed = mSeeds[iseed]; + if (seed.origID.getSource() != GTrackID::TPC) { + return candidates; + } + constexpr float MaxFlightMUS = 0.1f; // flight time of the muon between the TPC and the TOF, slow tails + constexpr float RadiusTOF = 380.f; // [cm], to find the sector the leg points to + const auto& tpcTrack = mRecoData->getTPCTrack(seed.origID); + const o2::track::TrackPar& parOut = tpcTrack.getParamOut(); + std::array xyz{}; + std::array dir{}; + parOut.getXYZGlo(xyz); + parOut.getPxPyPzGlo(dir); + // straight line from the outer parameters to the TOF radius (the sector only) + const float a = dir[0] * dir[0] + dir[1] * dir[1]; + const float b = xyz[0] * dir[0] + xyz[1] * dir[1]; + const float c = xyz[0] * xyz[0] + xyz[1] * xyz[1] - RadiusTOF * RadiusTOF; + const float disc = b * b - a * c; + if (a <= 0.f || disc < 0.f) { + return candidates; + } + const float step = (-b + std::sqrt(disc)) / a; + const int sectorCentre = o2::math_utils::angle2Sector(std::atan2(xyz[1] + step * dir[1], xyz[0] + step * dir[0])); + constexpr int NSectors = 18; + o2::track::TrackPar parSector[3]; + bool okSector[3] = {false, false, false}; + for (int k = 0; k < 3; k++) { + parSector[k] = parOut; + okSector[k] = parSector[k].rotateParam(o2::math_utils::sector2Angle((sectorCentre + k - 1 + NSectors) % NSectors)); + } + const float road = mMatchParams->tofRoad; + // range of z(t) - z(tRef) = side * vD * (t - tRef) over the bracket + const float dzDrift0 = seed.tpcSide * mTPCVDrift * (seed.tBracket.getMin() - seed.tRef); + const float dzDrift1 = seed.tpcSide * mTPCVDrift * (seed.tBracket.getMax() - seed.tRef); + const float dzMin = std::min(dzDrift0, dzDrift1) - road; + const float dzMax = std::max(dzDrift0, dzDrift1) + road; + const auto clusters = mRecoData->getTOFClusters(); + auto first = std::lower_bound(mTOFClusterTimeMUS.begin(), mTOFClusterTimeMUS.end(), seed.tBracket.getMin() - MaxFlightMUS); + for (auto it = first; it != mTOFClusterTimeMUS.end() && *it <= seed.tBracket.getMax() + MaxFlightMUS; ++it) { + const int index = mTOFClusterOrder[it - mTOFClusterTimeMUS.begin()]; + const auto& cl = clusters[index]; + const int k = (cl.getSector() - sectorCentre + NSectors + 1) % NSectors; // 0, 1, 2 for the sectors before, at and after the centre + if (k > 2 || !okSector[k]) { + continue; + } + float y = 0.f; + float z = 0.f; + if (!parSector[k].getYZAt(cl.getX(), mBz, y, z)) { + continue; + } + const float dy = cl.getY() - y; + const float dz = cl.getZ() - z; + if (std::abs(dy) > road || dz < dzMin || dz > dzMax) { + continue; + } + const float alpha = o2::math_utils::sector2Angle(cl.getSector()); + const float sinAlpha = std::sin(alpha); + const float cosAlpha = std::cos(alpha); + candidates.push_back(TOFCandidate{index, cl.getTime() * 1e-3, dy, dz, cl.getX() * cosAlpha - cl.getY() * sinAlpha, cl.getX() * sinAlpha + cl.getY() * cosAlpha, cl.getZ()}); + } + return candidates; +} + +//________________________________________________________ +float MatchCosmics::findTOFFlightPair(int i, int j, float tMinMUS, float tMaxMUS, float& tofTimeMUS) +{ + // best pair of TOF candidates of seeds i and j whose time difference matches the muon's flight between them along the helix (the higher + // hit first), with its mean time within [tMinMUS, tMaxMUS] and the z of both legs in the road at that time. Returns its score, the + // squared residuals in units of their cuts (< 0: no pair); tofTimeMUS is the mean time of the two hits + constexpr float MaxFlightMUS = 0.1f; + constexpr float CmPerNS = 29.9792458f; + const auto& candidates0 = getTOFCandidates(i); + const auto& candidates1 = getTOFCandidates(j); + if (candidates0.empty() || candidates1.empty()) { + return -1.f; + } + const auto& seed0 = mSeeds[i]; + const auto& seed1 = mSeeds[j]; + const float curvature = 0.5f * (std::abs(mRecoData->getTPCTrack(seed0.origID).getCurvature(mBz)) + std::abs(mRecoData->getTPCTrack(seed1.origID).getCurvature(mBz))); + const float road = mMatchParams->tofRoad; + const float tolerance = mMatchParams->tofFlightTolerance; + float bestScore = -1.f; + for (const auto& c0 : candidates0) { + for (const auto& c1 : candidates1) { + if (c0.index == c1.index) { + continue; + } + const auto& top = c0.gy > c1.gy ? c0 : c1; + const auto& bottom = c0.gy > c1.gy ? c1 : c0; + // flight path along the helix: arc in the transverse plane from the chord, then the dip + const float chordXY = std::hypot(top.gx - bottom.gx, top.gy - bottom.gy); + const float halfAngleSin = 0.5f * curvature * chordXY; + const float arcXY = halfAngleSin > 1e-4f && halfAngleSin < 1.f ? 2.f * std::asin(halfAngleSin) / curvature : chordXY; + const float length = std::hypot(arcXY, top.gz - bottom.gz); + const float flightDev = float(top.timeNS - bottom.timeNS) + length / CmPerNS; // the muon crosses the top TOF first + if (std::abs(flightDev) > tolerance) { + continue; + } + const float pairTimeMUS = float(0.5e-3 * (c0.timeNS + c1.timeNS)); + if (pairTimeMUS < tMinMUS - MaxFlightMUS || pairTimeMUS > tMaxMUS + MaxFlightMUS) { + continue; + } + const float dz0 = c0.dz - seed0.tpcSide * mTPCVDrift * (pairTimeMUS - seed0.tRef); + const float dz1 = c1.dz - seed1.tpcSide * mTPCVDrift * (pairTimeMUS - seed1.tRef); + if (std::abs(dz0) > road || std::abs(dz1) > road) { + continue; + } + const float score = (c0.dy * c0.dy + c1.dy * c1.dy + dz0 * dz0 + dz1 * dz1) / (road * road) + flightDev * flightDev / (tolerance * tolerance); + if (bestScore < 0.f || score < bestScore) { + bestScore = score; + tofTimeMUS = pairTimeMUS; + } + } + } + return bestScore; +} + +//________________________________________________________ +void MatchCosmics::registerMatch(int i, int j, float chi2, float tCommon, float tCommonErr, float tofScore, float tCommonNoTOF, float tCommonErrNoTOF) { - /// register track index j as a match for track index i + /// register track index j as a match for track index i; the matches of i are ordered in chi2, those confirmed by a TOF flight pair + /// (tofScore >= 0) in front of all others int newRef = mRecords.size(); - auto& matchRec = mRecords.emplace_back(MatchRecord{i, j, chi2, MinusOne}); + auto& matchRec = mRecords.emplace_back(MatchRecord{i, j, chi2, MinusOne, tCommon, tCommonErr, tofScore, tCommonNoTOF, tCommonErrNoTOF}); + const bool confirmed = tofScore >= 0.f; auto* best = &mSeeds[i].matchID; while (*best > MinusOne) { auto& oldMatchRec = mRecords[*best]; - if (oldMatchRec.chi2 > chi2) { // insert new match in front of the old one + const bool oldConfirmed = oldMatchRec.tofScore >= 0.f; + if ((confirmed && !oldConfirmed) || (confirmed == oldConfirmed && oldMatchRec.chi2 > chi2)) { // insert new match in front of the old one matchRec.next = *best; // new record will refer to the one it is superseding *best = newRef; // the reference on the superseded record should now refer to new one break; @@ -540,7 +914,7 @@ void MatchCosmics::createSeeds(const o2::globaltracking::RecoContainer& data) return true; } if constexpr (isTPCTrack()) { - if (!this->mMatchParams->allowTPCOnly) { + if (!this->mMatchParams->allowTPCOnly || _tr.getNClusters() < this->mMatchParams->minSeedNClTPC) { return true; } // unconstrained TPC track, with t0 = TrackTPC.getTime0+0.5*(DeltaFwd-DeltaBwd) and terr = 0.5*(DeltaFwd+DeltaBwd) in TimeBins @@ -554,9 +928,14 @@ void MatchCosmics::createSeeds(const o2::globaltracking::RecoContainer& data) } terr += this->mMatchParams->timeToleranceMUS; trackEntry[_origID] = mSeeds.size(); - mSeeds.emplace_back(TrackSeed{_tr, {t0 - terr, t0 + terr}, _origID, MinusOne}); + auto& seed = mSeeds.emplace_back(TrackSeed{_tr, {t0 - terr, t0 + terr}, _origID, MinusOne}); + seed.setPID(o2::track::PID::Muon, true); // muon mass and charge for the material corrections, whatever the leg's dE/dx PID + seed.getXYZGlo(seed.xyzRef); + seed.tRef = t0; if constexpr (isTPCTrack()) { - mSeeds.back().setCov(this->mMatchParams->tpcExtraZError2 + _tr.getCov()[o2::track::kSigZ2], o2::track::kSigZ2); + seed.setCov(this->mMatchParams->tpcExtraZError2 + _tr.getCov()[o2::track::kSigZ2], o2::track::kSigZ2); + seed.tRef = _tr.getTime0() * this->mTPCTBinMUS; // the z of a TPC-only track refers to its time0 + seed.tpcSide = _tr.hasASideClustersOnly() ? 1 : (_tr.hasCSideClustersOnly() ? -1 : 0); } return true; } else { @@ -564,8 +943,9 @@ void MatchCosmics::createSeeds(const o2::globaltracking::RecoContainer& data) } }; - data.createTracksVariadic(creator); + data.createTracksVariadic(creator, mSeedSources); // other loaded sources only resolve the primary-vertex contributors below + mNSeedsPVContributors = 0; if (mUsePVInfo) { // if needed, veto with the primary vertex info auto trackIndex = data.getPrimaryVertexMatchedTracks(); // Global ID's for associated tracks auto vtxRefs = data.getPrimaryVertexMatchedTrackRefs(); // references from vertex to these track IDs @@ -578,11 +958,28 @@ void MatchCosmics::createSeeds(const o2::globaltracking::RecoContainer& data) auto tvid = trackIndex[it]; auto entry = trackEntry.find(tvid); if (entry == trackEntry.end()) { + // a contributor of a veto source (not itself a leg): the legs it contains (e.g. its TPC track) come from a collision + if (mMatchParams->discardPVContributors && tvid.isPVContributor() && mPVVetoSources[tvid.getSource()] && data.isTrackSourceLoaded(tvid.getSource())) { + const auto refs = data.getSingleDetectorRefs(tvid); + for (int src = 0; src < GTrackID::NSources; src++) { + if (!mSeedSources[src] || !refs[src].isIndexSet()) { + continue; + } + auto partEntry = trackEntry.find(refs[src]); + if (partEntry != trackEntry.end() && mSeeds[partEntry->second].matchID != Reject) { + mSeeds[partEntry->second].matchID = Reject; + mNSeedsPVContributors++; + } + } + } continue; } auto& seed = mSeeds[entry->second]; - if (seed.matchID == Reject || (mMatchParams->discardPVContributors && tvid.isPVContributor())) { + if (seed.matchID != Reject && mMatchParams->discardPVContributors && tvid.isPVContributor()) { // the leg itself is a contributor seed.matchID = Reject; + mNSeedsPVContributors++; + } + if (seed.matchID == Reject) { continue; } if (seed.vtIDMin < 0) { @@ -603,10 +1000,13 @@ void MatchCosmics::updateTimeDependentParams() mBz = o2::base::Propagator::Instance()->getNominalBz(); mFieldON = std::abs(mBz) > 0.01; mQ2PtCutoff = 1.f / std::max(0.05f, mMatchParams->minSeedPt); + mQ2PtCutoffOppositeSides = mMatchParams->minPtOppositeSides > 0.f ? 1.f / mMatchParams->minPtOppositeSides : 1e9; if (mFieldON) { mQ2PtCutoff *= 5.00668 / std::abs(mBz); + mQ2PtCutoffOppositeSides *= 5.00668 / std::abs(mBz); } else { mQ2PtCutoff = 1e9; + mQ2PtCutoffOppositeSides = 1e9; } } diff --git a/Detectors/GlobalTracking/src/MatchCosmicsParams.cxx b/Detectors/GlobalTracking/src/MatchCosmicsParams.cxx index f14ae04897c68..7886558f7fab0 100644 --- a/Detectors/GlobalTracking/src/MatchCosmicsParams.cxx +++ b/Detectors/GlobalTracking/src/MatchCosmicsParams.cxx @@ -14,4 +14,39 @@ /// \author ruben.shahoyan@cern.ch #include "GlobalTracking/MatchCosmicsParams.h" +#include "Framework/Logger.h" +#include + O2ParamImpl(o2::globaltracking::MatchCosmicsParams); + +std::string o2::globaltracking::getMatchCosmicsPreset(const std::string& name) +{ + // physics-v1, tuned on PbPb 2025 (LHC25an 567939, 19 kHz) for fakes, on the cosmics run 562658 and on a cosmic MC with PbPb-like + // distortions for the efficiency (TPC-only legs): + // - seed cuts against collision tracks: pT, transverse DCA to the beam line in cm and in sigma, TPC clusters + // - realistic systematic errors for the leg comparison (true pairs had y / snp / q/pt pulls 2-4x too wide), so that the crude pair + // chi2 cut is meaningful; per-parameter windows open except tgl (3 sigma: the main discriminant against random pairs) + // - loose cut on the chi2 of the refitted legs; z test of TPC-only legs at a common time and same-half veto switched on + // - TPC-only legs on opposite sides (time from z continuity, ~94 % of the PbPb candidates were random pairs): both legs and the + // refitted cosmic above 2 GeV (offline: PbPb candidates / 10.5, cosmic MC efficiency 81.8 -> 80.6 %) + // - TOF flight pair as selection criterion (needs TOF clusters, added to the inputs): pairs confirmed by a top / bottom TOF hit pair + // with the muon's flight time win the selection and are refitted at the TOF time (568041: fewer cosmics, more of them TOF-tagged); + // the common-time refit of same-side legs (refitSameSideAtCommonTime) is not used: in PbPb it added mostly collision-track pairs + static const std::map presets{ + {"physics-v1", + "cosmicsMatch.minSeedPt=1;cosmicsMatch.minSeedDCAxy=3;cosmicsMatch.minSeedDCAxyNSigma=10;cosmicsMatch.minSeedNClTPC=30;" + "cosmicsMatch.crudeChi2Cut=50;cosmicsMatch.systSigma2[0]=0.25;cosmicsMatch.systSigma2[2]=4e-4;cosmicsMatch.systSigma2[4]=2.5e-3;" + "cosmicsMatch.crudeNSigma2Cut[0]=144;cosmicsMatch.crudeNSigma2Cut[2]=144;cosmicsMatch.crudeNSigma2Cut[3]=9;cosmicsMatch.crudeNSigma2Cut[4]=144;" + "cosmicsMatch.maxChi2Match=1000;cosmicsMatch.constrainTPCOnlyZ=true;cosmicsMatch.vetoSameHalf=true;cosmicsMatch.minPtOppositeSides=2;" + "cosmicsMatch.tofFlightSelection=true"}}; + auto it = presets.find(name); + if (it == presets.end()) { + std::string known; + for (const auto& [k, v] : presets) { + known += " " + k; + } + LOGP(fatal, "Unknown cosmics matching preset {}, known:{}", name, known); + return {}; + } + return it->second; +} diff --git a/Detectors/GlobalTrackingWorkflow/CMakeLists.txt b/Detectors/GlobalTrackingWorkflow/CMakeLists.txt index 6f29ab1930f95..69f5ff9379944 100644 --- a/Detectors/GlobalTrackingWorkflow/CMakeLists.txt +++ b/Detectors/GlobalTrackingWorkflow/CMakeLists.txt @@ -26,6 +26,7 @@ o2_add_library(GlobalTrackingWorkflow src/StrangenessTrackingWriterSpec.cxx src/CosmicsMatchingSpec.cxx src/TrackCosmicsWriterSpec.cxx + src/CosmicsClusterCollectorSpec.cxx src/TOFMatcherSpec.cxx src/TOFMatchChecker.cxx src/HMPMatcherSpec.cxx diff --git a/Detectors/GlobalTrackingWorkflow/include/GlobalTrackingWorkflow/CosmicsClusterCollectorSpec.h b/Detectors/GlobalTrackingWorkflow/include/GlobalTrackingWorkflow/CosmicsClusterCollectorSpec.h new file mode 100644 index 0000000000000..2db4d81c2293e --- /dev/null +++ b/Detectors/GlobalTrackingWorkflow/include/GlobalTrackingWorkflow/CosmicsClusterCollectorSpec.h @@ -0,0 +1,31 @@ +// Copyright 2019-2026 CERN and copyright holders of ALICE O2. +// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. +// All rights not expressly granted are reserved. +// +// This software is distributed under the terms of the GNU General Public +// License v3 (GPL Version 3), copied verbatim in the file "COPYING". +// +// In applying this license CERN does not waive the privileges and immunities +// granted to it by virtue of its status as an Intergovernmental Organization +// or submit itself to any jurisdiction. + +/// @file CosmicsClusterCollectorSpec.h +/// @brief Collects the raw clusters of matched cosmics (attached + road around the legs) for offline refits + +#ifndef O2_COSMICS_CLUSTER_COLLECTOR_SPEC_H +#define O2_COSMICS_CLUSTER_COLLECTOR_SPEC_H + +#include "Framework/DataProcessorSpec.h" +#include "ReconstructionDataFormats/GlobalTrackID.h" +#include "DetectorsCommonDataFormats/DetID.h" + +namespace o2::globaltracking +{ + +/// create a processor spec collecting the clusters of the cosmics found by the cosmics-matcher +/// roadDets: detectors (ITS, TOF, TRD) whose hits along the cosmic are collected besides the TPC road +framework::DataProcessorSpec getCosmicsClusterCollectorSpec(o2::dataformats::GlobalTrackID::mask_t src, bool useMC, bool itsStag, o2::detectors::DetID::mask_t roadDets); + +} // namespace o2::globaltracking + +#endif diff --git a/Detectors/GlobalTrackingWorkflow/include/GlobalTrackingWorkflow/CosmicsMatchingSpec.h b/Detectors/GlobalTrackingWorkflow/include/GlobalTrackingWorkflow/CosmicsMatchingSpec.h index 1b1d9c494dc6c..a04c2e1f81000 100644 --- a/Detectors/GlobalTrackingWorkflow/include/GlobalTrackingWorkflow/CosmicsMatchingSpec.h +++ b/Detectors/GlobalTrackingWorkflow/include/GlobalTrackingWorkflow/CosmicsMatchingSpec.h @@ -24,7 +24,14 @@ namespace globaltracking { /// create a processor spec -framework::DataProcessorSpec getCosmicsMatchingSpec(o2::dataformats::GlobalTrackID::mask_t src, bool usePV, bool useMC, bool itsStag); +/// ITS tracks needed by legs with an ITS part (their refit and ITS clusters) if ITS tracks are no legs themselves +o2::dataformats::GlobalTrackID::mask_t getLegITSSources(o2::dataformats::GlobalTrackID::mask_t src); + +/// srcPVContributors plus the parent matches needed to resolve their single-detector parts (RecoContainer::getSingleDetectorRefs) +o2::dataformats::GlobalTrackID::mask_t addPVContributorParents(o2::dataformats::GlobalTrackID::mask_t srcPVContributors); + +framework::DataProcessorSpec getCosmicsMatchingSpec(o2::dataformats::GlobalTrackID::mask_t src, bool usePV, bool useMC, bool itsStag, bool useTOFClusters = false, + o2::dataformats::GlobalTrackID::mask_t srcPVContributors = {}); } // namespace globaltracking } // namespace o2 diff --git a/Detectors/GlobalTrackingWorkflow/include/GlobalTrackingWorkflow/TrackCosmicsWriterSpec.h b/Detectors/GlobalTrackingWorkflow/include/GlobalTrackingWorkflow/TrackCosmicsWriterSpec.h index 9cae3afa7b1d5..fd149b7c3c401 100644 --- a/Detectors/GlobalTrackingWorkflow/include/GlobalTrackingWorkflow/TrackCosmicsWriterSpec.h +++ b/Detectors/GlobalTrackingWorkflow/include/GlobalTrackingWorkflow/TrackCosmicsWriterSpec.h @@ -27,6 +27,9 @@ namespace globaltracking /// write cosmics tracks to a root file framework::DataProcessorSpec getTrackCosmicsWriterSpec(bool useMC); +/// write cosmics with their raw clusters (cosmics-cluster-collector output) to a root file +framework::DataProcessorSpec getCosmicsFullWriterSpec(); + } // namespace globaltracking } // namespace o2 diff --git a/Detectors/GlobalTrackingWorkflow/src/CosmicsClusterCollectorSpec.cxx b/Detectors/GlobalTrackingWorkflow/src/CosmicsClusterCollectorSpec.cxx new file mode 100644 index 0000000000000..c91318243ec0d --- /dev/null +++ b/Detectors/GlobalTrackingWorkflow/src/CosmicsClusterCollectorSpec.cxx @@ -0,0 +1,1658 @@ +// Copyright 2019-2026 CERN and copyright holders of ALICE O2. +// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. +// All rights not expressly granted are reserved. +// +// This software is distributed under the terms of the GNU General Public +// License v3 (GPL Version 3), copied verbatim in the file "COPYING". +// +// In applying this license CERN does not waive the privileges and immunities +// granted to it by virtue of its status as an Intergovernmental Organization +// or submit itself to any jurisdiction. + +/// @file CosmicsClusterCollectorSpec.cxx +/// @brief Collects the raw clusters of matched cosmics (attached + road around the legs) for offline refits +/// +/// The async reconstruction stores neither TPC tracks nor TPC clusters, so the leg references of TrackCosmics cannot be resolved +/// offline. For every cosmic this device stores the TPC tracks of both legs, the raw TPC clusters attached to them, all raw TPC clusters +/// in a road around each leg (gaps, split pieces, delta electrons; clusters attached to other tracks are flagged), the ITS / TOF / TRD +/// hits of the legs' matched tracks and along roads in these detectors, and per TF the quantities of the TPC transformation. +/// +/// Road: each leg's helix is propagated through all pad rows of all sectors on its TPC side, in the frame of its own time0 (the frame +/// in which its clusters are consistent). Points on the other side of the leg's closest approach to the beam line belong to the other +/// leg and are skipped. The predicted real (y, z) is mapped to nominal coordinates with the inverse correction, and clusters with +/// nominal coordinates within the corridor width of the predicted point (distance perpendicular to the track) are taken. The other TPC +/// side is searched only when the time of the cosmic is known (z-continuity time of legs on opposite sides with a small error, or a TOF +/// time); a CE-crossing leg's own frame already covers both sides. Parts of a leg nearly parallel to the pad rows (|snp| >= 0.8 in the +/// sector frame, e.g. around a closest approach to the beam line inside the TPC, where the track runs along one pad row) cannot be +/// reached row by row: there the helix is followed in steps of path length and the rows within the road width of each point are searched. +/// +/// TOF tag: all TOF clusters close to the outward extrapolation of the legs within the cosmic's time window are candidates; the top / +/// bottom pair whose time difference matches the muon's flight along the helix (HitTOFFlight) gives the cosmic its TOF time. The same +/// search in the impossible order (bottom hit first) only finds accidental pairs and is kept as QA of the flag's background. +/// +/// Roads in the other detectors (--road-detectors): TRD tracklets close to the outward extrapolation of the legs and ITS clusters close +/// to the trajectory near the beam line, within the time window of the cosmic (the TOF time if there is one, else the matcher's time, +/// which is precise for TPC-only legs on opposite sides; otherwise the time window of the legs, with a correspondingly loose z cut). +/// They are flagged as found on the road; hits of the legs' matched global tracks are flagged as matched. In the ITS the road keeps per +/// half of the cosmic and layer the cluster closest to the trajectory, in the inner barrel (dense with collision clusters near the beam +/// line) the 5 closest, best first. Cosmics sharing >= 30 % of their TPC clusters with a better one (a leg split into two TPC tracks) are +/// flagged as duplicates, nothing is removed. +/// +/// Polish: a cosmic with a TOF time and TPC-only legs is refitted at that time (as in the matcher: muon mass, energy loss along the +/// flight); one-side legs then have their real z, hence the right material, which the matcher's TPC time cannot always give. +/// +/// Besides the matcher's and the polished track, only raw detector data are stored (TPC ClusterNative, ITS compact clusters + patterns, +/// raw TOF time, TRD tracklet words); calibrations, the cluster dictionary and the geometry are applied offline. +/// +/// --debug-tree writes cosmics_collector_debug.root for test runs: tree "cosmics" with one entry per cosmic, its matching and timing +/// quantities and vectors per detector: cl* every stored TPC cluster transformed (local x, y, z in the frame used by the road, zCos in the +/// common frame of the cosmic's time, global gx, gy), road* the predicted track points of each leg (same frames; roadFrame 0 / 1 own +/// / cosmic's time frame row by row, 2 / 3 the same along the low-angle walk), its* the +/// local (with the ITS road also global) coordinates of the ITS clusters, tof* the raw and calibrated TOF time and the cluster position, +/// trd* the road tracklets and their residuals. E.g. Draw("clGy:clGx", "scoreTOF >= 0") (TOF-tagged; tTOF can be negative). + +#include +#include +#include +#include +#include +#include +#include +#include +#include "TStopwatch.h" +#include "Framework/Task.h" +#include "Framework/ConfigParamRegistry.h" +#include "Framework/DataProcessorSpec.h" +#include "Framework/DeviceSpec.h" +#include "GlobalTrackingWorkflow/CosmicsClusterCollectorSpec.h" +#include "GlobalTrackingWorkflow/CosmicsMatchingSpec.h" +#include "DataFormatsGlobalTracking/RecoContainer.h" +#include "DataFormatsGlobalTracking/TrackCosmicsExtended.h" +#include "ReconstructionDataFormats/TrackCosmics.h" +#include "DataFormatsTPC/TrackTPC.h" +#include "DataFormatsTPC/ClusterNative.h" +#include "DataFormatsITS/TrackITS.h" +#include "DataFormatsITSMFT/CompCluster.h" +#include "DataFormatsITSMFT/ROFRecord.h" +#include "DataFormatsITSMFT/TopologyDictionary.h" +#include "DataFormatsITSMFT/ClusterPattern.h" +#include "ITSMFTBase/SegmentationAlpide.h" +#include "DataFormatsTOF/Cluster.h" +#include "DataFormatsTRD/TrackTRD.h" +#include "DataFormatsTRD/Tracklet64.h" +#include "DataFormatsTRD/TriggerRecord.h" +#include "DataFormatsTRD/CalibratedTracklet.h" +#include "DataFormatsTRD/Constants.h" +#include "DataFormatsITSMFT/DPLAlpideParam.h" +#include "ITSBase/GeometryTGeo.h" +#include "DetectorsCommonDataFormats/DetID.h" +#include "CommonConstants/LHCConstants.h" +#include "CommonDataFormat/TFIDInfo.h" +#include "CommonDataFormat/InteractionRecord.h" +#include "DetectorsBase/GRPGeomHelper.h" +#include "DetectorsBase/Propagator.h" +#include "GPUO2InterfaceRefit.h" +#include "GlobalTracking/MatchCosmicsParams.h" +#include "DataFormatsTPC/WorkflowHelper.h" +#include "ReconstructionDataFormats/Vertex.h" +#include "DetectorsBase/TFIDInfoHelper.h" +#include "TPCBase/ParameterElectronics.h" +#include "MathUtils/Utils.h" +#include "MathUtils/Primitive2D.h" +#include "TPCFastTransformPOD.h" +#include "CommonUtils/TreeStreamRedirector.h" + +using namespace o2::framework; +using GTrackID = o2::dataformats::GlobalTrackID; +using TPCGeo = o2::gpu::TPCFastTransformGeoPOD; +using DetID = o2::detectors::DetID; + +namespace o2::globaltracking +{ + +namespace +{ +constexpr float TanSector = 0.17632698f; // tan(10 deg): half opening of a sector + +/// z outside the drift volume of one TPC side by more than margin +bool outsideDriftVolume(bool sideA, float z, float margin) +{ + const float zLength = TPCGeo::getTPCzLength(); + return sideA ? (z < -margin || z > zLength + margin) : (z > margin || z < -zLength - margin); +} + +/// TPC side of a track with clusters on one side only (+1 A, -1 C: its z moves with the time assumed for its clusters), 0 otherwise +int tpcSide(const o2::tpc::TrackTPC& trk) +{ + return trk.hasASideClustersOnly() ? 1 : (trk.hasCSideClustersOnly() ? -1 : 0); +} +} // namespace + +class CosmicsClusterCollectorSpec : public Task +{ + public: + CosmicsClusterCollectorSpec(std::shared_ptr dr, std::shared_ptr gr, bool useMC, DetID::mask_t roadDets) : mDataRequest(dr), mGGCCDBRequest(gr), mRoadDets(roadDets), mUseMC(useMC) {} + ~CosmicsClusterCollectorSpec() override = default; + void init(InitContext& ic) final; + void run(ProcessingContext& pc) final; + void endOfStream(EndOfStreamContext& ec) final; + void finaliseCCDB(ConcreteDataMatcher& matcher, void* obj) final; + + private: + /// which side of the closest approach to the beam line (transverse) a point is on, relative to the leg's own clusters + struct LegBranch { + bool isLine = false; ///< straight line (no field / very high pT) instead of a circle + float centerX = 0.f; ///< circle centre + float centerY = 0.f; + float pcaX = 0.f; ///< point of closest approach to the beam line + float pcaY = 0.f; + float dirX = 0.f; ///< direction of the straight line + float dirY = 0.f; + int sign = 0; ///< side of the leg's clusters; 0: the leg spans both sides, accept everything + int side(float x, float y) const + { + const float orientation = isLine ? (x - pcaX) * dirX + (y - pcaY) * dirY : (pcaX - centerX) * (y - centerY) - (pcaY - centerY) * (x - centerX); + return orientation > 0.f ? 1 : -1; + } + bool accept(float x, float y) const { return sign == 0 || side(x, y) == sign; } + void init(const o2::track::TrackPar& inner, const o2::track::TrackPar& outer, float bz); + }; + + /// time of the cosmic in TPC time bins + struct CosmicTime { + bool known = false; ///< known well enough to search the other TPC side of one-side legs + float tb = 0.f; ///< time [TB] + float errTB = 0.f; ///< its error [TB] + }; + + void updateTimeDependentParams(ProcessingContext& pc); + void buildUsedMap(const RecoContainer& data, std::vector>& time0Windows); + void addTPCAttached(const RecoContainer& data, const o2::tpc::TrackTPC& trk, std::vector& out, std::unordered_set& taken) const; + /// time frame of a TPC road: the leg's z is shifted by dz and its clusters are transformed with vertexTime + struct RoadFrame { + float vertexTime; ///< vertex time used for the transformation [TB] + float dz; ///< shift of the leg's z into this frame + float zTolerance; ///< extra z tolerance from the error of the vertex time + int sectorMin; + int sectorMax; + uint8_t flag; + }; + void addTPCCorridor(const RecoContainer& data, const o2::tpc::TrackTPC& trk, const CosmicTime& cosmicTime, std::vector& out, std::unordered_set& taken) const; + void walkLowAngleRoad(const o2::tpc::ClusterNativeAccess& clusters, const o2::track::TrackPar& start, bool innerPart, float rMid, const LegBranch& branch, const RoadFrame& frame, + std::vector& out, std::unordered_set& taken) const; + void searchRow(const o2::tpc::ClusterNativeAccess& clusters, int sector, int row, float y, float z, float snp, float tgl, float vertexTime, float zTolerance, uint8_t flag, + float dxRow, std::vector& out, std::unordered_set& taken) const; + void addDebugRoadPoint(int sector, float x, float y, float z, float snp, float tgl, int row, float vertexTime, int frameCode) const; + struct ITSPattRequest { + int cosmic; ///< entry of the cosmic in the output + int cluster; ///< entry of the cluster in its clITS + int index; ///< index of the cluster in the TF + }; + void addITS(const RecoContainer& data, GTrackID gid, uint8_t leg, std::vector& out, int icosm, std::vector& requests, std::unordered_set& matched) const; + void fillITSPatterns(const RecoContainer& data, std::vector& requests, std::vector& cosmics) const; + void addTOF(const RecoContainer& data, GTrackID gid, uint8_t leg, std::vector& out, std::unordered_set& matched) const; + void addTRD(const RecoContainer& data, GTrackID gid, uint8_t leg, std::vector& out, std::unordered_set& matched) const; + std::pair timeWindowMUS(const CosmicTime& cosmicTime) const; + bool predictOutward(const o2::tpc::TrackTPC& leg, const CosmicTime& cosmicTime, int sector, float x, float& y, float& z) const; + void roadTOF(const RecoContainer& data, const o2::tpc::TrackTPC* const* legs, const CosmicTime& cosmicTime, int icosm, std::vector& out, const std::unordered_set& matched, float& timeTOFMUS, float& scoreTOFPair, float& scoreTOFReversed) const; + void polish(const o2::tpc::TrackTPC* const* legs, float timeMUS, o2::gpu::GPUO2InterfaceRefit& refitter, o2::dataformats::TrackCosmicsExtended& out) const; + void roadTRD(const RecoContainer& data, const o2::tpc::TrackTPC* const* legs, const CosmicTime& cosmicTime, int icosm, std::vector& out, const std::unordered_set& matched) const; + void roadITS(const RecoContainer& data, const o2::dataformats::TrackCosmics& cosm, const CosmicTime& cosmicTime, int legsSide, int icosm, std::vector& out, std::vector& requests, + const std::unordered_set& matched) const; + void cacheITSChipCentres(); + void writeDebug(const o2::dataformats::TrackCosmicsExtended& cosm, int icosm) const; + void flagDuplicates(std::vector& cosmics) const; + void writeDebugTOF(const o2::tof::Cluster& c, int icosm, int leg, uint8_t flags) const; + + std::shared_ptr mDataRequest; + std::shared_ptr mGGCCDBRequest; + const o2::gpu::TPCFastTransformPOD* mCorrMap = nullptr; + const o2::itsmft::TopologyDictionary* mITSDict = nullptr; + std::vector mUsed; ///< TPC clusters attached to any TPC track in this TF + o2::InteractionRecord mTFStart{}; ///< first BC of the TF + float mCorridor = 1.f; ///< road half-width [cm] + float mMaxAbsTimeErr = 0.5f; ///< max. time error of the cosmic [mus] to search the other TPC side of a one-side leg + size_t mMaxCosmicsPerTF = 100; ///< cosmics processed per TF at most (protection against fake-dominated settings) + DetID::mask_t mRoadDets{}; ///< detectors searched along the road besides the TPC + float mRoadTOF = 5.f; ///< road half-width at the TOF [cm] + float mTOFFlightTol = 2.f; ///< max. deviation of the top/bottom TOF time difference from the flight time [ns] + float mTOFTimeErr = 0.1f; ///< error of the TOF time of a cosmic for the later roads [mus] (covers TPC vs TOF offsets) + float mRoadTRD = 5.f; ///< road half-width at the TRD [cm] (in z plus the pad length) + float mRoadITS = 1.5f; ///< road half-width in the ITS [cm] + std::vector> mITSChipCentres; ///< global positions of the ITS chip centres (aligned geometry) + float mTPCTBinMUS = 0.2f; ///< TPC time bin [mus] + float mBz = 0.f; + bool mUseMC = false; + mutable bool mStaggeredWarned = false; + mutable bool mNoDictWarned = false; + mutable bool mTRDMismatchWarned = false; + std::unique_ptr mDebugOut; ///< debug tree "cosmics", one entry per cosmic (--debug-tree) + int mDbgTF = 0; ///< context of the debug output: TF counter + int mDbgCosmic = 0; ///< entry of the cosmic in the TF + int mDbgLeg = 0; ///< leg (0 bottom, 1 top) + float mDbgTauC = 0.f; ///< time of the cosmic [TB]: common frame (zCos) of the debug output + /// debug quantities found while processing a cosmic (road points, TOF and TRD hits), written with its entry + struct DebugCosmic { + std::vector roadLeg; + std::vector roadFrame; + std::vector roadSector; + std::vector roadRow; + std::vector roadX; + std::vector roadY; + std::vector roadZ; + std::vector roadZCos; + std::vector roadGx; + std::vector roadGy; + std::vector roadSnp; + std::vector roadTgl; + std::vector tofLeg; + std::vector tofChannel; + std::vector tofFlags; + std::vector tofTimeRaw; + std::vector tofTime; + std::vector tofTot; + std::vector tofX; + std::vector tofY; + std::vector tofZ; + std::vector tofGx; + std::vector tofGy; + std::vector trdLeg; + std::vector trdLayer; + std::vector trdX; + std::vector trdY; + std::vector trdZ; + std::vector trdDy; + std::vector trdDz; + std::vector trdTrigMUS; + }; + mutable std::vector mDbgCosmics; ///< per cosmic of the TF + size_t mNCosmics = 0; ///< cosmics processed + size_t mNClAttached = 0; ///< attached TPC clusters stored + size_t mNClCorridor = 0; ///< road TPC clusters stored + TStopwatch mTimer; +}; + +void CosmicsClusterCollectorSpec::init(InitContext& ic) +{ + mTimer.Stop(); + mTimer.Reset(); + o2::base::GRPGeomHelper::instance().setRequest(mGGCCDBRequest); + mCorridor = ic.options().get("corridor-width"); + mMaxAbsTimeErr = ic.options().get("max-abs-time-err"); + mMaxCosmicsPerTF = std::max(0, ic.options().get("max-cosmics-per-tf")); + mRoadTOF = ic.options().get("tof-road-width"); + mTOFFlightTol = ic.options().get("tof-flight-tolerance"); + mTOFTimeErr = ic.options().get("tof-time-error"); + mRoadTRD = ic.options().get("trd-road-width"); + mRoadITS = ic.options().get("its-road-width"); + if (ic.options().get("debug-tree")) { + const auto timesliceId = ic.services().get().inputTimesliceId; + const std::string name = timesliceId == 0 ? "cosmics_collector_debug.root" : fmt::format("cosmics_collector_debug_{}.root", timesliceId); + mDebugOut = std::make_unique(name.c_str(), "recreate"); + } +} + +void CosmicsClusterCollectorSpec::run(ProcessingContext& pc) +{ + mTimer.Start(false); + RecoContainer recoData; + recoData.collectData(pc, *mDataRequest.get()); + updateTimeDependentParams(pc); + + o2::dataformats::TFIDInfo tfID; + o2::base::TFIDInfoHelper::fillTFIDInfo(pc, tfID); + mTFStart = {0, tfID.firstTForbit}; + o2::dataformats::CosmicsTFInfo tfInfo; + tfInfo.vDrift = mCorrMap->getVDrift(); + tfInfo.t0 = mCorrMap->getT0(); + + std::vector cosmicsOut; + const auto cosmics = recoData.getCosmicTracks(); + const size_t nCosmics = std::min(cosmics.size(), mMaxCosmicsPerTF); + if (nCosmics < cosmics.size()) { + LOGP(warning, "{} cosmics in TF {}, only the first {} are written (max-cosmics-per-tf)", cosmics.size(), tfID.tfCounter, nCosmics); + } + if (nCosmics) { + // only TPC tracks whose time0 is within two drift times of a leg's time0 can share clusters with the roads + const float maxDistTB = 2.2f * TPCGeo::getTPCzLength() / mCorrMap->getVDrift(); + std::vector> time0Windows; + for (size_t ic = 0; ic < nCosmics; ic++) { + for (const auto leg : {cosmics[ic].getRefBottom(), cosmics[ic].getRefTop()}) { + const auto refs = recoData.getSingleDetectorRefs(leg); + if (refs[GTrackID::TPC].isIndexSet()) { + const float time0 = recoData.getTPCTrack(refs[GTrackID::TPC]).getTime0(); + time0Windows.emplace_back(time0 - maxDistTB, time0 + maxDistTB); + } + } + } + buildUsedMap(recoData, time0Windows); + } + std::vector pattRequests; // ITS clusters whose patterns are not in the dictionary + std::unique_ptr tpcRefitter; // polish of the cosmics with a TOF time + if (nCosmics && mRoadDets[DetID::TOF] && recoData.inputsTPCclusters) { + const auto& clusterShMap = recoData.clusterShMapTPC; + const auto& occupancyMap = recoData.occupancyMapTPC; + tpcRefitter = std::make_unique(&recoData.inputsTPCclusters->clusterIndex, mCorrMap, mBz, recoData.getTPCTracksClusterRefs().data(), 0, + clusterShMap.data(), occupancyMap.data(), occupancyMap.size(), nullptr, o2::base::Propagator::Instance()); + } + mDbgTF = tfID.tfCounter; + if (mDebugOut) { + mDbgCosmics.assign(nCosmics, DebugCosmic{}); + } + for (size_t ic = 0; ic < nCosmics; ic++) { + const auto& cosm = cosmics[ic]; + auto& out = cosmicsOut.emplace_back(); + out.cosmic = cosm; + if (mUseMC) { + out.label = recoData.getCosmicTrackMCLabel(ic); + } + const GTrackID legs[2] = {cosm.getRefBottom(), cosm.getRefTop()}; + const o2::tpc::TrackTPC* tpcLegs[2] = {nullptr, nullptr}; + std::vector* tpcCl[2] = {&out.clTPCBottom, &out.clTPCTop}; + std::unordered_set matchedITS; // hits of the legs' matched tracks, not searched again on the road + std::unordered_set matchedTOF; + std::unordered_set matchedTRD; + for (uint8_t leg = 0; leg < 2; leg++) { + auto refs = recoData.getSingleDetectorRefs(legs[leg]); + if (refs[GTrackID::TPC].isIndexSet()) { + tpcLegs[leg] = &recoData.getTPCTrack(refs[GTrackID::TPC]); + (leg == 0 ? out.tpcBottom : out.tpcTop) = *tpcLegs[leg]; + } + if (refs[GTrackID::ITS].isIndexSet()) { + addITS(recoData, refs[GTrackID::ITS], leg, out.clITS, ic, pattRequests, matchedITS); + } + if (refs[GTrackID::TOF].isIndexSet()) { + addTOF(recoData, refs[GTrackID::TOF], leg, out.clTOF, matchedTOF); + if (mDebugOut) { + writeDebugTOF(recoData.getTOFClusters()[refs[GTrackID::TOF].getIndex()], ic, leg, o2::dataformats::HitMatched); + } + } + if (refs[GTrackID::TRD].isIndexSet()) { + addTRD(recoData, refs[GTrackID::TRD], leg, out.trdTracklets, matchedTRD); + } + } + CosmicTime cosmicTime; + cosmicTime.tb = cosm.getTimeMUS().getTimeStamp() / mTPCTBinMUS; + cosmicTime.errTB = cosm.getTimeMUS().getTimeStampError() / mTPCTBinMUS; + cosmicTime.known = cosm.getTimeMUS().getTimeStampError() < mMaxAbsTimeErr; + mDbgCosmic = ic; + // TOF road first: a top/bottom hit pair matching the muon's flight gives the cosmic's time to ~ns (also for one-side legs on the same + // side, whose brackets leave it open by tens of mus); the TPC corridor of the other side and the TRD / ITS roads then use that time + if (mRoadDets[DetID::TOF]) { + roadTOF(recoData, tpcLegs, cosmicTime, ic, out.clTOF, matchedTOF, out.timeTOFMUS, out.scoreTOFPair, out.scoreTOFReversed); + } + if (out.hasTOFTime() && tpcRefitter && legs[0].getSource() == GTrackID::TPC && legs[1].getSource() == GTrackID::TPC) { + polish(tpcLegs, out.timeTOFMUS, *tpcRefitter, out); + } + CosmicTime roadTime = cosmicTime; + if (out.hasTOFTime()) { + roadTime.tb = out.timeTOFMUS / mTPCTBinMUS; + roadTime.errTB = mTOFTimeErr / mTPCTBinMUS; + roadTime.known = true; + } + mDbgTauC = roadTime.tb; + std::unordered_set taken; + for (int leg = 0; leg < 2; leg++) { // attached clusters of both legs first, so that a road never takes the other leg's clusters + if (tpcLegs[leg]) { + addTPCAttached(recoData, *tpcLegs[leg], *tpcCl[leg], taken); + mNClAttached += tpcCl[leg]->size(); + } + } + for (int leg = 0; leg < 2; leg++) { + if (tpcLegs[leg]) { + const size_t nAttached = tpcCl[leg]->size(); + mDbgLeg = leg; + addTPCCorridor(recoData, *tpcLegs[leg], roadTime, *tpcCl[leg], taken); + mNClCorridor += tpcCl[leg]->size() - nAttached; + } + } + if (mRoadDets[DetID::TRD]) { + roadTRD(recoData, tpcLegs, roadTime, ic, out.trdTracklets, matchedTRD); + } + if (mRoadDets[DetID::ITS]) { + // the refitted cosmic's z moves with the time only if both legs are TPC-only on one side (an ITS part fixes it absolutely) + int legsSide = 0; + if (tpcLegs[0] && tpcLegs[1] && legs[0].getSource() == GTrackID::TPC && legs[1].getSource() == GTrackID::TPC && tpcSide(*tpcLegs[0]) == tpcSide(*tpcLegs[1])) { + legsSide = tpcSide(*tpcLegs[0]); + } + roadITS(recoData, cosm, roadTime, legsSide, ic, out.clITS, pattRequests, matchedITS); + } + } + if (!pattRequests.empty()) { + fillITSPatterns(recoData, pattRequests, cosmicsOut); + } + flagDuplicates(cosmicsOut); + if (mDebugOut) { + for (size_t ic = 0; ic < cosmicsOut.size(); ic++) { + writeDebug(cosmicsOut[ic], ic); + } + } + mNCosmics += cosmicsOut.size(); + LOGP(info, "Collected clusters for {} cosmics in TF {}", cosmicsOut.size(), tfID.tfCounter); + pc.outputs().snapshot(Output{"GLO", "COSMFULL", 0}, cosmicsOut); + pc.outputs().snapshot(Output{"GLO", "COSMFULLTF", 0}, tfInfo); + pc.outputs().snapshot(Output{"GLO", "COSMFULLTFID", 0}, tfID); + mTimer.Stop(); +} + +void CosmicsClusterCollectorSpec::updateTimeDependentParams(ProcessingContext& pc) +{ + o2::base::GRPGeomHelper::instance().checkUpdates(pc); + mCorrMap = &o2::gpu::TPCFastTransformPOD::get(pc.inputs().get("corrMap")); + mTPCTBinMUS = o2::tpc::ParameterElectronics::Instance().ZbinWidth; + mBz = o2::base::Propagator::Instance()->getNominalBz(); + if (mRoadDets[DetID::ITS] && mITSChipCentres.empty()) { + cacheITSChipCentres(); + } +} + +void CosmicsClusterCollectorSpec::cacheITSChipCentres() +{ + auto geom = o2::its::GeometryTGeo::Instance(); + geom->fillMatrixCache(o2::math_utils::bit2Mask(o2::math_utils::TransformType::L2G)); + mITSChipCentres.resize(geom->getNumberOfChips()); + for (int chip = 0; chip < geom->getNumberOfChips(); chip++) { + mITSChipCentres[chip] = geom->getMatrixL2G(chip) * o2::math_utils::Point3D(0.f, 0.f, 0.f); + } +} + +void CosmicsClusterCollectorSpec::buildUsedMap(const RecoContainer& data, std::vector>& time0Windows) +{ + // merge the windows, then flag the clusters of the TPC tracks whose time0 falls into one of them + std::sort(time0Windows.begin(), time0Windows.end()); + std::vector> merged; + for (const auto& window : time0Windows) { + if (!merged.empty() && window.first <= merged.back().second) { + merged.back().second = std::max(merged.back().second, window.second); + } else { + merged.push_back(window); + } + } + const auto& clusters = data.getTPCClusters(); + const auto tracks = data.getTPCTracks(); + const auto refs = data.getTPCTracksClusterRefs(); + mUsed.assign(clusters.nClustersTotal, false); + for (const auto& trk : tracks) { + const float time0 = trk.getTime0(); + const auto window = std::upper_bound(merged.begin(), merged.end(), time0, [](float t, const std::pair& w) { return t < w.first; }); + if (window == merged.begin() || time0 > std::prev(window)->second) { + continue; + } + for (int j = 0; j < trk.getNClusterReferences(); j++) { + uint8_t sector = 0; + uint8_t row = 0; + uint32_t clIdx = 0; + trk.getClusterReference(refs, j, sector, row, clIdx); + mUsed[clusters.clusterOffset[sector][row] + clIdx] = true; + } + } +} + +void CosmicsClusterCollectorSpec::addTPCAttached(const RecoContainer& data, const o2::tpc::TrackTPC& trk, std::vector& out, std::unordered_set& taken) const +{ + const auto& clusters = data.getTPCClusters(); + const auto refs = data.getTPCTracksClusterRefs(); + for (int j = 0; j < trk.getNClusterReferences(); j++) { + uint8_t sector = 0; + uint8_t row = 0; + uint32_t clIdx = 0; + trk.getClusterReference(refs, j, sector, row, clIdx); + if (!taken.insert(clusters.clusterOffset[sector][row] + clIdx).second) { + continue; + } + auto& cl = out.emplace_back(); + cl.cl = clusters.clusters[sector][row][clIdx]; + cl.sector = sector; + cl.row = row; + cl.flags = o2::dataformats::CosmicTPCCluster::Attached | o2::dataformats::CosmicTPCCluster::Used; + } +} + +void CosmicsClusterCollectorSpec::LegBranch::init(const o2::track::TrackPar& inner, const o2::track::TrackPar& outer, float bz) +{ + if (std::abs(inner.getCurvature(bz)) < 1e-5f) { // radius > 1 km: straight line + isLine = true; + const auto point = inner.getXYZGlo(); + const float phi = inner.getPhi(); + dirX = std::cos(phi); + dirY = std::sin(phi); + const float proj = point.X() * dirX + point.Y() * dirY; + pcaX = point.X() - proj * dirX; + pcaY = point.Y() - proj * dirY; + } else { + o2::math_utils::CircleXYf_t circle; + float sinAlpha = 0.f; + float cosAlpha = 0.f; + inner.getCircleParams(bz, circle, sinAlpha, cosAlpha); + centerX = circle.xC; + centerY = circle.yC; + const float centerDist = std::sqrt(centerX * centerX + centerY * centerY); + if (centerDist < 1e-3f) { // circle around the beam line: no closest approach + sign = 0; + return; + } + pcaX = centerX * (1.f - circle.rC / centerDist); + pcaY = centerY * (1.f - circle.rC / centerDist); + } + // side of the leg from its end farther from the closest approach; a leg with ends on both sides (both > 10 cm away) spans the + // closest approach and is accepted on both sides + const auto pointIn = inner.getXYZGlo(); + const auto pointOut = outer.getXYZGlo(); + const float distIn2 = (pointIn.X() - pcaX) * (pointIn.X() - pcaX) + (pointIn.Y() - pcaY) * (pointIn.Y() - pcaY); + const float distOut2 = (pointOut.X() - pcaX) * (pointOut.X() - pcaX) + (pointOut.Y() - pcaY) * (pointOut.Y() - pcaY); + const int sideIn = side(pointIn.X(), pointIn.Y()); + const int sideOut = side(pointOut.X(), pointOut.Y()); + constexpr float MinDist2 = 10.f * 10.f; + if (sideIn == sideOut) { + sign = sideIn; + } else if (std::min(distIn2, distOut2) < MinDist2) { + sign = distIn2 > distOut2 ? sideIn : sideOut; + } else { + sign = 0; + } +} + +void CosmicsClusterCollectorSpec::addTPCCorridor(const RecoContainer& data, const o2::tpc::TrackTPC& trk, const CosmicTime& cosmicTime, std::vector& out, std::unordered_set& taken) const +{ + constexpr int NSectorsA = TPCGeo::getNumberOfSectorsA(); + const float zLength = TPCGeo::getTPCzLength(); + const float vDrift = mCorrMap->getVDrift(); + const auto& clusters = data.getTPCClusters(); + const float time0Leg = trk.getTime0(); + const int legSide = tpcSide(trk); + + RoadFrame frames[2]; + int nFrames = 0; + if (legSide == 0) { // CE-crossing leg: its time0 is absolute + frames[nFrames++] = {time0Leg, 0.f, 0.f, 0, 2 * NSectorsA, 0}; + } else { + frames[nFrames++] = {time0Leg, 0.f, 0.f, legSide > 0 ? 0 : NSectorsA, legSide > 0 ? NSectorsA : 2 * NSectorsA, 0}; + if (cosmicTime.known) { // z of a one-side leg moves by side * vD * (t - time0) when its clusters are transformed with the time t + frames[nFrames++] = {cosmicTime.tb, legSide * (cosmicTime.tb - time0Leg) * vDrift, cosmicTime.errTB * vDrift, + legSide > 0 ? NSectorsA : 0, legSide > 0 ? 2 * NSectorsA : NSectorsA, o2::dataformats::CosmicTPCCluster::AbsTime}; + } + } + + const o2::track::TrackPar refPar[2] = {trk, trk.getParamOut()}; // rows below rMid: inner parameter, above: outer + LegBranch branch; + branch.init(refPar[0], refPar[1], mBz); + const auto pointIn = refPar[0].getXYZGlo(); + const auto pointOut = refPar[1].getXYZGlo(); + const float rIn = std::hypot(pointIn.X(), pointIn.Y()); + const float rOut = std::hypot(pointOut.X(), pointOut.Y()); + const float rMid = 0.5f * (rIn + rOut); + + for (int iFrame = 0; iFrame < nFrames; iFrame++) { + const auto& frame = frames[iFrame]; + for (int sector = frame.sectorMin; sector < frame.sectorMax; sector++) { + const float alpha = o2::math_utils::sector2Angle(sector % NSectorsA); + const float sinAlpha = std::sin(alpha); + const float cosAlpha = std::cos(alpha); + const bool sideA = sector < NSectorsA; + for (int iPar = 0; iPar < 2; iPar++) { + auto par = refPar[iPar]; + par.setZ(par.getZ() + frame.dz); + if (!par.rotateParam(alpha)) { // the leg points away from this sector frame: same helix, opposite direction + par.invertParam(); + if (!par.rotateParam(alpha)) { + continue; + } + } + for (int row = 0; row < o2::tpc::constants::MAXGLOBALPADROW; row++) { + const float xRow = TPCGeo::getRowInfoX(row); + if ((xRow < rMid) != (iPar == 0)) { + continue; + } + auto parRow = par; + if (!parRow.propagateParamTo(xRow, mBz)) { + continue; + } + // coarse acceptance at the nominal x of the row before the more expensive move to its real x + constexpr float CoarseMargin = 10.f; // [cm] covers the shift of y and z between the nominal and the real x + if (std::abs(parRow.getY()) > xRow * TanSector + mCorridor + CoarseMargin || outsideDriftVolume(sideA, parRow.getZ(), mCorridor + frame.zTolerance + CoarseMargin)) { + continue; + } + // the corrected clusters of this row lie at the real x of the row, x + dx(y, z), not at its nominal x + float xReal = xRow; + const float zInside = sideA ? std::clamp(parRow.getZ(), 0.f, zLength) : std::clamp(parRow.getZ(), -zLength, 0.f); + mCorrMap->InverseTransformYZtoX(sector, row, parRow.getY(), zInside, xReal); + if (!parRow.propagateParamTo(xReal, mBz)) { + continue; + } + const float y = parRow.getY(); + const float z = parRow.getZ(); + if (std::abs(y) > xRow * TanSector + mCorridor) { + continue; + } + if (outsideDriftVolume(sideA, z, mCorridor + frame.zTolerance)) { + continue; + } + if (!branch.accept(xRow * cosAlpha - y * sinAlpha, xRow * sinAlpha + y * cosAlpha)) { + continue; + } + if (mDebugOut) { + addDebugRoadPoint(sector, xReal, y, z, parRow.getSnp(), parRow.getTgl(), row, frame.vertexTime, frame.flag != 0); + } + searchRow(clusters, sector, row, y, z, parRow.getSnp(), parRow.getTgl(), frame.vertexTime, frame.zTolerance, frame.flag, 0.f, out, taken); + } + } + } + // parts of the leg nearly parallel to the pad rows, which the propagation row by row cannot reach (e.g. a closest approach to the + // beam line inside the TPC: in the sector of the closest approach the track runs along one pad row); a leg whose inner end is not + // inside its outer end (all its clusters along one pad row) is walked from the inner end alone + const bool splitAtRMid = rOut > rIn; + for (int iPar = 0; iPar < (splitAtRMid ? 2 : 1); iPar++) { + auto par = refPar[iPar]; + par.setZ(par.getZ() + frame.dz); + walkLowAngleRoad(clusters, par, iPar == 0, splitAtRMid ? rMid : std::numeric_limits::max(), branch, frame, out, taken); + } + } +} + +void CosmicsClusterCollectorSpec::walkLowAngleRoad(const o2::tpc::ClusterNativeAccess& clusters, const o2::track::TrackPar& start, bool innerPart, float rMid, const LegBranch& branch, + const RoadFrame& frame, std::vector& out, std::unordered_set& taken) const +{ + constexpr int NSectorsA = TPCGeo::getNumberOfSectorsA(); + constexpr float CosSector = 0.98480775f; // cos(10 deg) + constexpr float Step = 0.5f; // [cm] path length per step, below the road width along the pad row + constexpr float MaxPath = 600.f; // [cm] protection; the walk normally ends at the TPC boundary or the closest approach + constexpr float MaxEntryPath = 20.f; // [cm] a start outside the leg's part (up to 10 cm past the closest approach, see LegBranch) is stepped into + constexpr float MinSnp = 0.8f; // below, the row crossings are well defined and found row by row + constexpr float XMargin = 10.f; // [cm] covers the shift between the nominal and the real x of a row (CoarseMargin of the row loop) + constexpr int NRows = o2::tpc::constants::MAXGLOBALPADROW; + const float zLength = TPCGeo::getTPCzLength(); + const float rLow = TPCGeo::getRowInfoX(0) - mCorridor; + const float rHigh = TPCGeo::getRowInfoX(NRows - 1) / CosSector + mCorridor; + std::array sinSector{}; + std::array cosSector{}; + for (int sector = 0; sector < NSectorsA; sector++) { + const float alpha = o2::math_utils::sector2Angle(sector); + sinSector[sector] = std::sin(alpha); + cosSector[sector] = std::cos(alpha); + } + + bool startInside = false; // the start point lies in the leg's part + for (int direction = 0; direction < 2; direction++) { + auto par = start; + if (direction == 1) { + par.invertParam(); + } + bool entered = direction == 1 && startInside; // the walk has reached the leg's part + // the start point is processed in the first direction only + for (float path = direction == 0 ? 0.f : Step; path < MaxPath; path += Step) { + // step along the helix in the frame of its direction (snp = 0), where the propagation is always defined + if (path > 0.f && (!par.rotateParam(par.getPhi()) || !par.propagateParamTo(par.getX() + Step, mBz))) { + break; + } + const auto pos = par.getXYZGlo(); + const float r = std::hypot(pos.X(), pos.Y()); + if (r > rHigh || std::abs(pos.Z()) > zLength + mCorridor + frame.zTolerance) { + break; + } + // outside the leg's part: past its closest approach to the beam line, beyond rMid or inside the inner radius; the start of the walk + // can lie there (a leg's end a few cm past the closest approach, an inner end inside the nominal inner radius) + if (r < rLow || (r < rMid) != innerPart || !branch.accept(pos.X(), pos.Y())) { + if (entered || path > MaxEntryPath) { + break; + } + continue; + } + entered = true; + startInside = startInside || path == 0.f; + const float phiDir = par.getPhi(); + const int sectorPos = o2::math_utils::angle2Sector(std::atan2(pos.Y(), pos.X())); + for (int dSector = -1; dSector <= 1; dSector++) { // the road can reach into the neighbouring sectors + const int sectorInSide = (sectorPos + dSector + NSectorsA) % NSectorsA; + for (int sector = sectorInSide; sector < 2 * NSectorsA; sector += NSectorsA) { + if (sector < frame.sectorMin || sector >= frame.sectorMax) { + continue; + } + const bool sideA = sector < NSectorsA; + if (outsideDriftVolume(sideA, pos.Z(), mCorridor + frame.zTolerance)) { + continue; + } + const float alpha = o2::math_utils::sector2Angle(sectorInSide); + const float sinAlpha = sinSector[sectorInSide]; + const float cosAlpha = cosSector[sectorInSide]; + const float x = pos.X() * cosAlpha + pos.Y() * sinAlpha; + const float y = -pos.X() * sinAlpha + pos.Y() * cosAlpha; + float snp = std::sin(phiDir - alpha); + float tgl = par.getTgl(); + if (std::cos(phiDir - alpha) < 0.f) { // the same line with its direction along +x of the sector frame, as in the row search + snp = -snp; + tgl = -tgl; + } + if (std::abs(snp) < MinSnp || std::abs(y) > x * TanSector + mCorridor) { + continue; + } + const float zInside = sideA ? std::clamp(pos.Z(), 0.f, zLength) : std::clamp(pos.Z(), -zLength, 0.f); + // rows whose real x lies within the road width of the point: the first by bisection on the nominal x, which grows with the row + int rowFirst = 0; + int rowEnd = NRows; + while (rowFirst < rowEnd) { + const int rowMid = (rowFirst + rowEnd) / 2; + if (TPCGeo::getRowInfoX(rowMid) < x - mCorridor - XMargin) { + rowFirst = rowMid + 1; + } else { + rowEnd = rowMid; + } + } + for (int row = rowFirst; row < NRows && TPCGeo::getRowInfoX(row) <= x + mCorridor + XMargin; row++) { + float xReal = TPCGeo::getRowInfoX(row); + mCorrMap->InverseTransformYZtoX(sector, row, y, zInside, xReal); + const float dxRow = xReal - x; + if (std::abs(dxRow) > mCorridor) { + continue; + } + if (mDebugOut) { + addDebugRoadPoint(sector, x, y, pos.Z(), snp, tgl, row, frame.vertexTime, 2 + (frame.flag != 0)); + } + searchRow(clusters, sector, row, y, pos.Z(), snp, tgl, frame.vertexTime, frame.zTolerance, frame.flag, dxRow, out, taken); + } + } + } + } + } +} + +void CosmicsClusterCollectorSpec::searchRow(const o2::tpc::ClusterNativeAccess& clusters, int sector, int row, float y, float z, float snp, float tgl, float vertexTime, float zTolerance, uint8_t flag, + float dxRow, std::vector& out, std::unordered_set& taken) const +{ + // (y, z) is the predicted point, dxRow the real x of the row minus the x of that point (0 when the point lies on the row) + const float zLength = TPCGeo::getTPCzLength(); + const float vDrift = mCorrMap->getVDrift(); + const float t0 = mCorrMap->getT0(); + // nominal (measured) coordinates of the predicted real point; the correction is evaluated inside the drift volume + const float zInside = sector < TPCGeo::getNumberOfSectorsA() ? std::clamp(z, 0.f, zLength) : std::clamp(z, -zLength, 0.f); + float yNominal = 0.f; + float zNominal = 0.f; + mCorrMap->InverseTransformYZtoNominalYZ(sector, row, y, zInside, yNominal, zNominal); + zNominal += z - zInside; + float padPred = 0.f; + float driftLengthPred = 0.f; + TPCGeo::convLocalToPadDriftLength(sector, row, yNominal, zNominal, padPred, driftLengthPred); + const float timePred = driftLengthPred / vDrift + t0 + vertexTime; + // the road is a cylinder around the track; its section with the pad-row plane is an ellipse with half-axes W / cos(phi) in y and + // W * sqrt(cos^2(phi) + tgl^2) / cos(phi) in z + const float cosPhi = std::sqrt((1.f - snp) * (1.f + snp)); + const float cosPhiWindow = std::max(cosPhi, 0.1f); // limits the window for tracks nearly parallel to the pad row + const float norm = 1.f / std::sqrt(1.f + tgl * tgl); + const float dirX = cosPhi * norm; // track direction + const float dirY = snp * norm; + const float dirZ = tgl * norm; + const float corridor2 = mCorridor * mCorridor; + const float windowPad = mCorridor / (cosPhiWindow * TPCGeo::getRowInfoPadWidth(row)); + const float windowTime = (mCorridor * std::sqrt(cosPhiWindow * cosPhiWindow + tgl * tgl) / cosPhiWindow + zTolerance) / vDrift; + const auto* rowClusters = clusters.clusters[sector][row]; + const uint32_t rowOffset = clusters.clusterOffset[sector][row]; + for (uint32_t k = 0; k < clusters.nClusters[sector][row]; k++) { + const auto& c = rowClusters[k]; + const float time = c.getTime(); + const float pad = c.getPad(); + if (std::abs(time - timePred) > windowTime || std::abs(pad - padPred) > windowPad) { + continue; + } + float yCluster = 0.f; + float zCluster = 0.f; + TPCGeo::convPadDriftLengthToLocal(sector, row, pad, (time - t0 - vertexTime) * vDrift, yCluster, zCluster); + const float dy = yCluster - yNominal; + float dz = zCluster - zNominal; + if (zTolerance > 0.f) { // the vertex time of this frame is known within zTolerance / vD: allow that shift along z + dz = std::abs(dz) > zTolerance ? dz - std::copysign(zTolerance, dz) : 0.f; + } + const float proj = dxRow * dirX + dy * dirY + dz * dirZ; + if (dxRow * dxRow + dy * dy + dz * dz - proj * proj > corridor2) { // distance perpendicular to the track + continue; + } + if (!taken.insert(rowOffset + k).second) { + continue; + } + auto& cl = out.emplace_back(); + cl.cl = c; + cl.sector = sector; + cl.row = row; + cl.flags = o2::dataformats::CosmicTPCCluster::Corridor | flag | (mUsed[rowOffset + k] ? o2::dataformats::CosmicTPCCluster::Used : 0); + } +} + +void CosmicsClusterCollectorSpec::addDebugRoadPoint(int sector, float x, float y, float z, float snp, float tgl, int row, float vertexTime, int frameCode) const +{ + // predicted point, z in the frame of the road and in the common frame of the cosmic's time + const float alpha = o2::math_utils::sector2Angle(sector % TPCGeo::getNumberOfSectorsA()); + const float sinAlpha = std::sin(alpha); + const float cosAlpha = std::cos(alpha); + const bool sideA = sector < TPCGeo::getNumberOfSectorsA(); + auto& dbg = mDbgCosmics[mDbgCosmic]; + dbg.roadLeg.push_back(mDbgLeg); + dbg.roadFrame.push_back(frameCode); + dbg.roadSector.push_back(sector); + dbg.roadRow.push_back(row); + dbg.roadX.push_back(x); + dbg.roadY.push_back(y); + dbg.roadZ.push_back(z); + dbg.roadZCos.push_back(z + (sideA ? 1.f : -1.f) * (mDbgTauC - vertexTime) * mCorrMap->getVDrift()); + dbg.roadGx.push_back(x * cosAlpha - y * sinAlpha); + dbg.roadGy.push_back(x * sinAlpha + y * cosAlpha); + dbg.roadSnp.push_back(snp); + dbg.roadTgl.push_back(tgl); +} + +void CosmicsClusterCollectorSpec::addITS(const RecoContainer& data, GTrackID gid, uint8_t leg, std::vector& out, int icosm, std::vector& requests, std::unordered_set& matched) const +{ + if (gid.getSource() != GTrackID::ITS) { // ITS-AB tracklets are not stored + return; + } + if (data.getITSPerLayer()) { + if (!mStaggeredWarned) { + LOGP(warning, "Staggered ITS clusters are not supported, ITS clusters of cosmics are not stored"); + mStaggeredWarned = true; + } + return; + } + const auto& trk = data.getITSTrack(gid); + const auto refs = data.getITSTracksClusterRefs(); + const auto clusters = data.getITSClusters(); + const auto rofs = data.getITSClustersROFRecords(); + for (int i = 0; i < trk.getNumberOfClusters(); i++) { + const int idx = refs[trk.getFirstClusterEntry() + i]; + const auto& c = clusters[idx]; + auto& cl = out.emplace_back(); + cl.chipID = c.getSensorID(); + cl.row = c.getRow(); + cl.col = c.getCol(); + cl.pattID = c.getPatternID(); + cl.leg = leg; + cl.flags = o2::dataformats::HitMatched; + matched.insert(idx); + if (c.getPatternID() == o2::itsmft::CompCluster::InvalidPatternID || (mITSDict && mITSDict->isGroup(c.getPatternID()))) { + requests.push_back({icosm, int(out.size()) - 1, idx}); // the pattern is in the TF's pattern stream, copied by fillITSPatterns + } + auto rof = std::upper_bound(rofs.begin(), rofs.end(), idx, [](int v, const o2::itsmft::ROFRecord& r) { return v < r.getFirstEntry(); }); + if (rof != rofs.begin()) { + cl.rofBC = std::prev(rof)->getBCData().differenceInBC(mTFStart); + } + } +} + +void CosmicsClusterCollectorSpec::addTOF(const RecoContainer& data, GTrackID gid, uint8_t leg, std::vector& out, std::unordered_set& matched) const +{ + const auto& c = data.getTOFClusters()[gid.getIndex()]; + auto& cl = out.emplace_back(); + cl.timeRaw = c.getTimeRaw(); + cl.tot = c.getTot(); + cl.channel = c.getMainContributingChannel(); + cl.leg = leg; + cl.flags = o2::dataformats::HitMatched; + matched.insert(gid.getIndex()); +} + +std::pair CosmicsClusterCollectorSpec::timeWindowMUS(const CosmicTime& cosmicTime) const +{ + // a time fixed by z continuity is given with a 1 sigma error, otherwise the error is the half-width of the legs' time-bracket overlap + const float timeMUS = cosmicTime.tb * mTPCTBinMUS; + const float errMUS = cosmicTime.errTB * mTPCTBinMUS; + const float halfWidth = cosmicTime.known ? 5.f * errMUS + 0.2f : errMUS; + return {timeMUS - halfWidth, timeMUS + halfWidth}; +} + +bool CosmicsClusterCollectorSpec::predictOutward(const o2::tpc::TrackTPC& leg, const CosmicTime& cosmicTime, int sector, float x, float& y, float& z) const +{ + // outward continuation of a leg in the frame of a sector; z in the frame of the cosmic's time (a one-side TPC-only leg has z relative to + // its time0) + o2::track::TrackPar par = leg.getParamOut(); + if (!par.rotateParam(o2::math_utils::sector2Angle(sector % TPCGeo::getNumberOfSectorsA())) || !par.propagateParamTo(x, mBz)) { + return false; + } + y = par.getY(); + z = par.getZ() + tpcSide(leg) * (cosmicTime.tb - leg.getTime0()) * mCorrMap->getVDrift(); + return true; +} + +void CosmicsClusterCollectorSpec::roadTOF(const RecoContainer& data, const o2::tpc::TrackTPC* const* legs, const CosmicTime& cosmicTime, int icosm, std::vector& out, const std::unordered_set& matched, float& timeTOFMUS, float& scoreTOFPair, float& scoreTOFReversed) const +{ + // TOF clusters along the legs' outward continuations (the road in PbPb also contains hits of collision tracks). Preferred: the top/bottom + // pair whose time difference matches the muon's flight between them; its mean time fixes the cosmic's time, also when the legs' brackets + // leave it open by tens of mus (one-side legs on the same side), and the z of one-side legs is shifted to it. Otherwise per leg the + // cluster closest to its continuation. The same search in the impossible order (bottom hit first) only finds accidental pairs: its best + // score is kept as QA of the flag's background. + constexpr float MaxFlightMUS = 0.1f; // flight time of the muon between the TPC and the TOF, slow tails + constexpr float CmPerNS = 29.9792458f; + const auto window = timeWindowMUS(cosmicTime); + const float zTimeTol = 0.5f * (window.second - window.first) / mTPCTBinMUS * mCorrMap->getVDrift(); // z uncertainty of one-side legs + const float vDriftPerMUS = mCorrMap->getVDrift() / mTPCTBinMUS; + const float cosmicTimeMUS = cosmicTime.tb * mTPCTBinMUS; + int legSide[2] = {0, 0}; + float curvature = 0.f; // |1/R| of the muon [1/cm] for the flight path between the TOF hits + int nLegs = 0; + for (int leg = 0; leg < 2; leg++) { + if (legs[leg]) { + legSide[leg] = tpcSide(*legs[leg]); + curvature += std::abs(legs[leg]->getCurvature(mBz)); + nLegs++; + } + } + curvature = nLegs ? curvature / nLegs : 0.f; + struct Candidate { + int index; + bool matched; // hit of the leg's matched global track (already stored) + double timeNS; // since the start of the TF + float dy; + float dzAtCosmicTime; + float gx; + float gy; + float gz; + }; + std::vector candidates[2]; + const auto clusters = data.getTOFClusters(); + int best[2] = {-1, -1}; + float bestScore[2] = {1.f, 1.f}; + for (int i = 0; i < (int)clusters.size(); i++) { + const auto& c = clusters[i]; + const float timeMUS = c.getTime() * 1e-6f; // [ps] since the start of the TF + if (timeMUS < window.first - MaxFlightMUS || timeMUS > window.second + MaxFlightMUS) { + continue; + } + const bool isMatched = matched.count(i) > 0; + const float alpha = o2::math_utils::sector2Angle(c.getSector()); + const float sinAlpha = std::sin(alpha); + const float cosAlpha = std::cos(alpha); + for (int leg = 0; leg < 2; leg++) { + float y = 0.f; + float z = 0.f; + if (!legs[leg] || !predictOutward(*legs[leg], cosmicTime, c.getSector(), c.getX(), y, z)) { + continue; + } + const float normY = (c.getY() - y) / mRoadTOF; + const float normZ = (c.getZ() - z) / (mRoadTOF + zTimeTol); + const float score = std::max(normY * normY, normZ * normZ); + if (score >= 1.f) { + continue; + } + if (!isMatched && score < bestScore[leg]) { + bestScore[leg] = score; + best[leg] = i; + } + candidates[leg].push_back(Candidate{i, isMatched, c.getTime() * 1e-3, c.getY() - y, c.getZ() - z, c.getX() * cosAlpha - c.getY() * sinAlpha, + c.getX() * sinAlpha + c.getY() * cosAlpha, c.getZ()}); + } + } + float bestPairScore = -1.f; + float bestReversedScore = -1.f; + int bestPair[2] = {-1, -1}; + for (const auto& c0 : candidates[0]) { + for (const auto& c1 : candidates[1]) { + if (c0.index == c1.index) { + continue; + } + // top / bottom by the hits' height: a muon goes down (for near-horizontal cosmics both orders are possible, but rare) + const auto& top = c0.gy > c1.gy ? c0 : c1; + const auto& bottom = c0.gy > c1.gy ? c1 : c0; + // flight path along the helix: arc in the transverse plane from the chord, then the dip + const float chordXY = std::hypot(top.gx - bottom.gx, top.gy - bottom.gy); + const float halfAngleSin = 0.5f * curvature * chordXY; + const float arcXY = halfAngleSin > 1e-4f && halfAngleSin < 1.f ? 2.f * std::asin(halfAngleSin) / curvature : chordXY; + const float length = std::hypot(arcXY, top.gz - bottom.gz); + const float timeDiff = float(top.timeNS - bottom.timeNS); + const float flightDev = timeDiff + length / CmPerNS; // the muon crosses the top TOF first + const float reversedDev = timeDiff - length / CmPerNS; // bottom hit first: impossible for a muon, only accidental pairs + const bool reversed = std::abs(flightDev) > mTOFFlightTol; // the two cannot both pass (L/c ~ 25 ns >> tolerance) + if (reversed && std::abs(reversedDev) > mTOFFlightTol) { + continue; + } + const double pairTimeNS = 0.5 * (c0.timeNS + c1.timeNS); + const float shiftMUS = float(pairTimeNS * 1e-3) - cosmicTimeMUS; + const float dz0 = c0.dzAtCosmicTime - legSide[0] * shiftMUS * vDriftPerMUS; + const float dz1 = c1.dzAtCosmicTime - legSide[1] * shiftMUS * vDriftPerMUS; + if (std::abs(dz0) > mRoadTOF || std::abs(dz1) > mRoadTOF) { + continue; + } + const float timeDev = reversed ? reversedDev : flightDev; + const float score = (c0.dy * c0.dy + c1.dy * c1.dy + dz0 * dz0 + dz1 * dz1) / (mRoadTOF * mRoadTOF) + timeDev * timeDev / (mTOFFlightTol * mTOFFlightTol); + if (reversed) { + if (bestReversedScore < 0.f || score < bestReversedScore) { + bestReversedScore = score; + } + continue; + } + if (bestPairScore < 0.f || score < bestPairScore) { + bestPairScore = score; + bestPair[0] = c0.index; + bestPair[1] = c1.index; + timeTOFMUS = float(pairTimeNS * 1e-3); + } + } + } + scoreTOFPair = bestPairScore; + scoreTOFReversed = bestReversedScore; + uint8_t flags = o2::dataformats::HitRoad; + if (bestPairScore >= 0.f) { + best[0] = bestPair[0]; + best[1] = bestPair[1]; + flags |= o2::dataformats::HitTOFFlight; + } else if (best[0] >= 0 && best[0] == best[1]) { // one TOF cluster in the roads of both legs: keep it for the closer one + best[bestScore[0] <= bestScore[1] ? 1 : 0] = -1; + } + for (int leg = 0; leg < 2; leg++) { + if (best[leg] < 0) { + continue; + } + const auto& c = clusters[best[leg]]; + if (matched.count(best[leg])) { // the leg's matched hit, already stored: flag it as part of the flight pair + for (auto& cl : out) { + if (cl.leg == leg && cl.channel == c.getMainContributingChannel() && cl.timeRaw == c.getTimeRaw()) { + cl.flags |= o2::dataformats::HitTOFFlight; + } + } + if (mDebugOut) { + auto& dbg = mDbgCosmics[icosm]; + for (size_t j = 0; j < dbg.tofFlags.size(); j++) { + if (dbg.tofLeg[j] == leg && dbg.tofChannel[j] == c.getMainContributingChannel() && dbg.tofTimeRaw[j] == c.getTimeRaw()) { + dbg.tofFlags[j] |= o2::dataformats::HitTOFFlight; + } + } + } + continue; + } + auto& cl = out.emplace_back(); + cl.timeRaw = c.getTimeRaw(); + cl.tot = c.getTot(); + cl.channel = c.getMainContributingChannel(); + cl.leg = leg; + cl.flags = flags; + if (mDebugOut) { + writeDebugTOF(c, icosm, leg, flags); + } + } +} + +void CosmicsClusterCollectorSpec::polish(const o2::tpc::TrackTPC* const* legs, float timeMUS, o2::gpu::GPUO2InterfaceRefit& refitter, o2::dataformats::TrackCosmicsExtended& out) const +{ + // refit of a cosmic with TPC-only legs at its TOF time, as MatchCosmics::refitWinners: the bottom leg inward, then to the closest approach + // to the beam line; the top leg inward, then to the same point; the two halves combined. Muon mass, energy loss along the muon's flight + // (top to bottom): going inward along the bottom leg and up to the closest approach is against the flight (gain), the top leg inward + // and down to the closest approach is with it (loss). With a precise time the one-side legs get their real z, hence also the right + // material, which the matcher's time (the legs' bracket overlap for one-side legs on the same side) cannot give. + constexpr int ELossGain = 1; + constexpr int ELossLoss = -1; + const auto& matchParams = o2::globaltracking::MatchCosmicsParams::Instance(); // propagation settings as in the matcher + const auto prop = o2::base::Propagator::Instance(); + const float timeTB = timeMUS / mTPCTBinMUS; + o2::track::TrackParCov bottom = legs[0]->getParamOut(); + bottom.setPID(o2::track::PID::Muon, true); + bottom.resetCovariance(); + if (std::abs(mBz) <= 0.01f) { + // B = 0: both legs carry the same conventional q/pt; as the matcher, flip the bottom one so that it is right after the inversion + bottom.setQ2Pt(-o2::track::kMostProbablePt); + } + if (refitter.RefitTrackAsTrackParCov(bottom, legs[0]->getClusterRef(), timeTB, nullptr, false, false, ELossGain) < 0) { + return; + } + bottom.invert(); + const o2::dataformats::VertexBase origin; + if (!prop->propagateToDCABxByBz(origin, bottom, matchParams.maxStep, matchParams.matCorr, nullptr, nullptr, ELossGain)) { + return; + } + o2::track::TrackParCov top = legs[1]->getParamOut(); + top.setPID(o2::track::PID::Muon, true); + if (refitter.RefitTrackAsTrackParCov(top, legs[1]->getClusterRef(), timeTB, nullptr, false, true, ELossLoss) < 0) { + return; + } + if (!top.rotate(bottom.getAlpha()) || !prop->PropagateToXBxByBz(top, bottom.getX(), matchParams.maxSnp, matchParams.maxStep, matchParams.matCorr, nullptr, ELossLoss)) { + return; + } + o2::track::TrackParCov::MatrixDSym5 cov5; + const float chi2Match = bottom.getPredictedChi2(top, cov5); + if (!bottom.update(top, cov5)) { + return; + } + out.polished = bottom; + out.chi2MatchPolished = chi2Match; +} + +void CosmicsClusterCollectorSpec::roadTRD(const RecoContainer& data, const o2::tpc::TrackTPC* const* legs, const CosmicTime& cosmicTime, int icosm, std::vector& out, const std::unordered_set& matched) const +{ + // per leg and layer the TRD tracklet closest to the leg's outward continuation, from triggers whose readout window can contain the cosmic + constexpr float ReadoutWindowMUS = 3.f; // a cosmic leaves tracklets only if it passes within the readout window after a trigger + constexpr float PadLength = 10.f; // [cm] longest TRD pads: the tracklet z is the pad-row centre + constexpr int NLayers = o2::trd::constants::NLAYER; + const auto window = timeWindowMUS(cosmicTime); + const float zTimeTol = 0.5f * (window.second - window.first) / mTPCTBinMUS * mCorrMap->getVDrift(); + const auto tracklets = data.getTRDTracklets(); + const auto calibrated = data.getTRDCalibratedTracklets(); + if (calibrated.size() != tracklets.size()) { + if (!mTRDMismatchWarned) { + LOGP(warning, "{} calibrated vs {} raw TRD tracklets: TRD road of cosmics skipped", calibrated.size(), tracklets.size()); + mTRDMismatchWarned = true; + } + return; + } + struct Candidate { + int tracklet = -1; + int trigBC = 0; + float score = 1.f; + }; + Candidate best[2][NLayers]; + for (const auto& trig : data.getTRDTriggerRecords()) { + const int trigBC = trig.getBCData().differenceInBC(mTFStart); + const float trigMUS = trigBC * o2::constants::lhc::LHCBunchSpacingMUS; + if (trigMUS < window.first - ReadoutWindowMUS || trigMUS > window.second) { + continue; + } + for (int it = trig.getFirstTracklet(); it < trig.getFirstTracklet() + trig.getNumberOfTracklets(); it++) { + if (matched.count(it)) { + continue; + } + const int detector = tracklets[it].getDetector(); + const int sector = detector / o2::trd::constants::NCHAMBERPERSEC; + const int layer = detector % NLayers; + const auto& point = calibrated[it]; + for (int leg = 0; leg < 2; leg++) { + float y = 0.f; + float z = 0.f; + if (!legs[leg] || !predictOutward(*legs[leg], cosmicTime, sector, point.getX(), y, z)) { + continue; + } + const float normY = (point.getY() - y) / mRoadTRD; + const float normZ = (point.getZ() - z) / (mRoadTRD + PadLength + zTimeTol); + const float score = std::max(normY * normY, normZ * normZ); + if (score < best[leg][layer].score) { + best[leg][layer] = {it, trigBC, score}; + } + } + } + } + for (int leg = 0; leg < 2; leg++) { + for (int layer = 0; layer < NLayers; layer++) { + const auto& cand = best[leg][layer]; + if (cand.tracklet < 0) { + continue; + } + auto& tr = out.emplace_back(); + tr.word = tracklets[cand.tracklet].getTrackletWord(); + tr.trigBC = cand.trigBC; + tr.layer = layer; + tr.leg = leg; + tr.flags = o2::dataformats::HitRoad; + if (mDebugOut) { + const auto& point = calibrated[cand.tracklet]; + const int sector = tracklets[cand.tracklet].getDetector() / o2::trd::constants::NCHAMBERPERSEC; + float y = 0.f; + float z = 0.f; + predictOutward(*legs[leg], cosmicTime, sector, point.getX(), y, z); // succeeded in the search + auto& dbg = mDbgCosmics[icosm]; + dbg.trdLeg.push_back(leg); + dbg.trdLayer.push_back(layer); + dbg.trdX.push_back(point.getX()); + dbg.trdY.push_back(point.getY()); + dbg.trdZ.push_back(point.getZ()); + dbg.trdDy.push_back(point.getY() - y); + dbg.trdDz.push_back(point.getZ() - z); + dbg.trdTrigMUS.push_back(cand.trigBC * float(o2::constants::lhc::LHCBunchSpacingMUS)); + } + } + } +} + +void CosmicsClusterCollectorSpec::roadITS(const RecoContainer& data, const o2::dataformats::TrackCosmics& cosm, const CosmicTime& cosmicTime, int legsSide, int icosm, std::vector& out, + std::vector& requests, const std::unordered_set& matched) const +{ + if (data.getITSPerLayer() || mITSChipCentres.empty()) { + return; + } + // trajectory near the beam line from the combined cosmic, sampled every cm along its direction (closest approach at local x = 0) + constexpr float MaxRadius = 45.f; // [cm] outer ITS barrel + margin + constexpr float ChipHalfDiag = 1.7f; // [cm] half diagonal of an ITS chip + o2::track::TrackPar par = cosm; + if (!par.rotateParam(par.getPhi())) { + return; + } + std::vector> points; + for (float x = -MaxRadius - 5.f; x <= MaxRadius + 5.f; x += 1.f) { + bool ok = false; + const auto point = par.getXYZGloAt(x, mBz, ok); + if (ok) { + points.push_back(point); + } + } + if (points.size() < 2) { + return; + } + // closest segment in the transverse plane: distance, z of the trajectory there + auto closest = [&points](float x, float y, float& dist2, float& zTraj) { + dist2 = 1e10f; + for (size_t i = 0; i + 1 < points.size(); i++) { + const float segX = points[i + 1].X() - points[i].X(); + const float segY = points[i + 1].Y() - points[i].Y(); + const float len2 = segX * segX + segY * segY; + const float frac = len2 > 0.f ? std::clamp(((x - points[i].X()) * segX + (y - points[i].Y()) * segY) / len2, 0.f, 1.f) : 0.f; + const float dx = x - (points[i].X() + frac * segX); + const float dy = y - (points[i].Y() + frac * segY); + if (dx * dx + dy * dy < dist2) { + dist2 = dx * dx + dy * dy; + zTraj = points[i].Z() + frac * (points[i + 1].Z() - points[i].Z()); + } + } + }; + float minR2 = 1e10f; + float pcaY = 0.f; + for (const auto& point : points) { + const float r2 = point.X() * point.X() + point.Y() * point.Y(); + if (r2 < minR2) { + minR2 = r2; + pcaY = point.Y(); + } + } + if (minR2 > MaxRadius * MaxRadius) { // the cosmic does not cross the ITS + return; + } + std::vector candidateChip(mITSChipCentres.size(), false); + bool anyChip = false; + const float chipCut2 = (mRoadITS + ChipHalfDiag) * (mRoadITS + ChipHalfDiag); + for (size_t chip = 0; chip < mITSChipCentres.size(); chip++) { + float dist2 = 0.f; + float zTraj = 0.f; + closest(mITSChipCentres[chip].X(), mITSChipCentres[chip].Y(), dist2, zTraj); + if (dist2 < chipCut2) { + candidateChip[chip] = true; + anyChip = true; + } + } + if (!anyChip) { + return; + } + const auto window = timeWindowMUS(cosmicTime); + const float zTimeTol = 0.5f * (window.second - window.first) / mTPCTBinMUS * mCorrMap->getVDrift(); // z of the cosmic is in the frame of its time + // the refitted cosmic has the z of its TPC time; with legs on one TPC side and a road time from the TOF its z moves by side * vD * dt + const float zShift = legsSide * (cosmicTime.tb - cosm.getTimeMUS().getTimeStamp() / mTPCTBinMUS) * mCorrMap->getVDrift(); + const float rofLengthMUS = o2::itsmft::DPLAlpideParam::Instance().roFrameLengthInBC * o2::constants::lhc::LHCBunchSpacingMUS; + // per half of the cosmic and layer the ITS clusters closest to the trajectory, best first: one in the outer barrel, several in the inner + // barrel, where the road near the beam line is dense with collision clusters and the TPC prediction (~mm) does not single out the hit + constexpr int NLayers = 7; + constexpr int NLayersIB = 3; + constexpr int MaxHitsIB = 5; + struct Candidate { + int index = -1; + int rofBC = 0; + float score = 1.f; + }; + Candidate best[2][NLayers][MaxHitsIB]; + const auto clusters = data.getITSClusters(); + auto geom = o2::its::GeometryTGeo::Instance(); + for (const auto& rof : data.getITSClustersROFRecords()) { + const int rofBC = rof.getBCData().differenceInBC(mTFStart); + const float rofMUS = rofBC * o2::constants::lhc::LHCBunchSpacingMUS; + if (rofMUS + rofLengthMUS < window.first || rofMUS > window.second) { + continue; + } + for (int idx = rof.getFirstEntry(); idx < rof.getFirstEntry() + rof.getNEntries(); idx++) { + const auto& c = clusters[idx]; + if (!candidateChip[c.getSensorID()] || matched.count(idx)) { + continue; + } + o2::math_utils::Point3D local; + if (mITSDict && c.getPatternID() != o2::itsmft::CompCluster::InvalidPatternID && !mITSDict->isGroup(c.getPatternID())) { + local = mITSDict->getClusterCoordinates(c); + } else { // anchor pixel: good enough for a cm road + o2::itsmft::SegmentationAlpide::detectorToLocalUnchecked(c.getRow(), c.getCol(), local); + } + const auto global = geom->getMatrixL2G(c.getSensorID()) * local; + float dist2 = 0.f; + float zTraj = 0.f; + closest(global.X(), global.Y(), dist2, zTraj); + const float normZ = (global.Z() - zTraj - zShift) / (mRoadITS + zTimeTol); + const float score = std::max(dist2 / (mRoadITS * mRoadITS), normZ * normZ); + const int half = global.Y() < pcaY ? 0 : 1; // bottom / top half of the cosmic + const int layer = geom->getLayer(c.getSensorID()); + if (layer < 0 || layer >= NLayers) { + continue; + } + auto* candidates = best[half][layer]; + int slot = (layer < NLayersIB ? MaxHitsIB : 1) - 1; + if (!(score < candidates[slot].score)) { // also rejects a NaN score + continue; + } + for (; slot > 0 && score < candidates[slot - 1].score; slot--) { + candidates[slot] = candidates[slot - 1]; + } + candidates[slot] = {idx, rofBC, score}; + } + } + for (int half = 0; half < 2; half++) { + for (int layer = 0; layer < NLayers; layer++) { + for (const auto& candidate : best[half][layer]) { + if (candidate.index < 0) { + break; + } + const auto& c = clusters[candidate.index]; + auto& cl = out.emplace_back(); + cl.chipID = c.getSensorID(); + cl.row = c.getRow(); + cl.col = c.getCol(); + cl.pattID = c.getPatternID(); + cl.rofBC = candidate.rofBC; + cl.leg = half; + cl.flags = o2::dataformats::HitRoad; + if (c.getPatternID() == o2::itsmft::CompCluster::InvalidPatternID || (mITSDict && mITSDict->isGroup(c.getPatternID()))) { + requests.push_back({icosm, int(out.size()) - 1, candidate.index}); + } + } + } + } +} + +void CosmicsClusterCollectorSpec::fillITSPatterns(const RecoContainer& data, std::vector& requests, std::vector& cosmics) const +{ + // the TF's pattern stream holds, in cluster order, the patterns of the clusters with an invalid or a group pattern ID + if (!mITSDict) { + if (!mNoDictWarned) { + LOGP(warning, "No ITS cluster dictionary: patterns of ITS clusters of cosmics are not stored"); + mNoDictWarned = true; + } + return; + } + std::sort(requests.begin(), requests.end(), [](const ITSPattRequest& a, const ITSPattRequest& b) { return a.index < b.index; }); + const auto clusters = data.getITSClusters(); + auto pattIt = data.getITSClustersPatterns().begin(); + size_t ir = 0; + for (int k = 0; k < (int)clusters.size() && ir < requests.size(); k++) { + const auto pattID = clusters[k].getPatternID(); + if (pattID != o2::itsmft::CompCluster::InvalidPatternID && !mITSDict->isGroup(pattID)) { + continue; + } + const auto start = pattIt; + o2::itsmft::ClusterPattern::skipPattern(pattIt); + for (; ir < requests.size() && requests[ir].index == k; ir++) { + auto& cosm = cosmics[requests[ir].cosmic]; + cosm.clITS[requests[ir].cluster].pattEntry = cosm.itsPatterns.size(); + cosm.itsPatterns.insert(cosm.itsPatterns.end(), start, pattIt); + } + } +} + +void CosmicsClusterCollectorSpec::addTRD(const RecoContainer& data, GTrackID gid, uint8_t leg, std::vector& out, std::unordered_set& matched) const +{ + const auto& trk = data.getTrack(gid); // TPC-TRD or ITS-TPC-TRD track + const auto tracklets = data.getTRDTracklets(); + const auto trigs = data.getTRDTriggerRecords(); + for (int layer = 0; layer < 6; layer++) { + const int it = trk.getTrackletIndex(layer); + if (it < 0) { + continue; + } + auto& tr = out.emplace_back(); + tr.word = tracklets[it].getTrackletWord(); + tr.layer = layer; + tr.leg = leg; + tr.flags = o2::dataformats::HitMatched; + matched.insert(it); + auto trig = std::upper_bound(trigs.begin(), trigs.end(), it, [](int v, const o2::trd::TriggerRecord& t) { return v < t.getFirstTracklet(); }); + if (trig != trigs.begin()) { + tr.trigBC = std::prev(trig)->getBCData().differenceInBC(mTFStart); + } + } +} + +void CosmicsClusterCollectorSpec::flagDuplicates(std::vector& cosmics) const +{ + // the same muon can be matched twice, e.g. when a leg is split into two TPC tracks: the roads then collect largely the same clusters + // (PbPb 567939: 10 of 47 cosmic pairs in a TF share 56-99 % of the clusters of the smaller one, all others none). The best one (TOF time, + // then more attached clusters) stays, the others point to it. + constexpr float MinSharedFraction = 0.3f; + const size_t n = cosmics.size(); + if (n < 2) { + return; + } + std::vector> keys(n); // sorted cluster addresses: sector, row and the packed time and pad of the raw cluster + std::vector nAttached(n, 0); + for (size_t ic = 0; ic < n; ic++) { + for (const auto* legClusters : {&cosmics[ic].clTPCBottom, &cosmics[ic].clTPCTop}) { + for (const auto& c : *legClusters) { + keys[ic].push_back((uint64_t(c.sector) << 56) | (uint64_t(c.row) << 48) | (uint64_t(c.cl.timeFlagsPacked) << 16) | c.cl.padPacked); + nAttached[ic] += c.isAttached(); + } + } + std::sort(keys[ic].begin(), keys[ic].end()); + } + std::vector order(n); + std::iota(order.begin(), order.end(), 0); + std::stable_sort(order.begin(), order.end(), [&](int a, int b) { + const bool timedA = cosmics[a].hasTOFTime(); + const bool timedB = cosmics[b].hasTOFTime(); + return timedA != timedB ? timedA : nAttached[a] > nAttached[b]; + }); + std::vector kept; + for (int ic : order) { + for (int ik : kept) { + std::vector shared; + std::set_intersection(keys[ic].begin(), keys[ic].end(), keys[ik].begin(), keys[ik].end(), std::back_inserter(shared)); + const size_t nMin = std::min(keys[ic].size(), keys[ik].size()); + if (nMin > 0 && shared.size() >= MinSharedFraction * nMin) { + cosmics[ic].duplicateOf = ik; + break; + } + } + if (cosmics[ic].duplicateOf < 0) { + kept.push_back(ic); + } + } +} + +void CosmicsClusterCollectorSpec::writeDebug(const o2::dataformats::TrackCosmicsExtended& cosm, int icosm) const +{ + constexpr int NSectorsA = TPCGeo::getNumberOfSectorsA(); + // the time the other-side corridor used (and the common frame zCos): the TOF time if the TOF road found one, else the TPC time + const float timeCosmic = (cosm.hasTOFTime() ? cosm.timeTOFMUS : cosm.cosmic.getTimeMUS().getTimeStamp()) / mTPCTBinMUS; + const bool absTimeKnown = cosm.cosmic.getTimeMUS().getTimeStampError() < mMaxAbsTimeErr; + int nAttached[2] = {0, 0}; + int nRoad[2] = {0, 0}; + int side[2] = {-2, -2}; + std::vector clLeg; + std::vector clFlags; + std::vector clSector; + std::vector clRow; + std::vector clPad; + std::vector clTime; + std::vector clQMax; + std::vector clQTot; + std::vector clX; + std::vector clY; + std::vector clZ; + std::vector clZCos; + std::vector clGx; + std::vector clGy; + const o2::tpc::TrackTPC* legs[2] = {&cosm.tpcBottom, &cosm.tpcTop}; + const std::vector* legClusters[2] = {&cosm.clTPCBottom, &cosm.clTPCTop}; + for (int leg = 0; leg < 2; leg++) { + if (legs[leg]->getNClusters() == 0) { // no TPC part + continue; + } + side[leg] = tpcSide(*legs[leg]); + for (const auto& c : *legClusters[leg]) { + // transform with the vertex time used by the road: the leg's time0, or the cosmic's time for clusters found on the other side + const float vertexTime = c.isAbsTime() ? timeCosmic : legs[leg]->getTime0(); + float x = 0.f; + float y = 0.f; + float z = 0.f; + mCorrMap->Transform(c.sector, c.row, c.cl.getPad(), c.cl.getTime(), x, y, z, vertexTime); + float xCos = 0.f; + float yCos = 0.f; + float zCos = 0.f; + mCorrMap->Transform(c.sector, c.row, c.cl.getPad(), c.cl.getTime(), xCos, yCos, zCos, timeCosmic); + const float alpha = o2::math_utils::sector2Angle(c.sector % NSectorsA); + const float sinAlpha = std::sin(alpha); + const float cosAlpha = std::cos(alpha); + clLeg.push_back(leg); + clFlags.push_back(c.flags); + clSector.push_back(c.sector); + clRow.push_back(c.row); + clPad.push_back(c.cl.getPad()); + clTime.push_back(c.cl.getTime()); + clQMax.push_back(c.cl.getQmax()); + clQTot.push_back(c.cl.getQtot()); + clX.push_back(x); + clY.push_back(y); + clZ.push_back(z); + clZCos.push_back(zCos); + clGx.push_back(x * cosAlpha - y * sinAlpha); + clGy.push_back(x * sinAlpha + y * cosAlpha); + (c.isAttached() ? nAttached : nRoad)[leg]++; + } + } + std::vector itsLeg; + std::vector itsFlags; + std::vector itsChip; + std::vector itsRow; + std::vector itsCol; + std::vector itsPattID; + std::vector itsHasPatt; + std::vector itsRofBC; + std::vector itsGx; + std::vector itsGy; + std::vector itsGz; + std::vector itsXLoc; + std::vector itsZLoc; + for (const auto& c : cosm.clITS) { // local coordinates on the chip from the dictionary or the stored pattern + const o2::itsmft::CompClusterExt compCluster(c.row, c.col, c.pattID, c.chipID); + o2::math_utils::Point3D local; + if (c.pattEntry >= 0) { + o2::itsmft::ClusterPattern pattern; + auto pattIt = cosm.itsPatterns.begin() + c.pattEntry; + pattern.acquirePattern(pattIt); + local = o2::itsmft::TopologyDictionary::getClusterCoordinates(compCluster, pattern, c.pattID != o2::itsmft::CompCluster::InvalidPatternID); + } else if (mITSDict && c.pattID != o2::itsmft::CompCluster::InvalidPatternID) { + local = mITSDict->getClusterCoordinates(compCluster); + } else { // pattern not available: pixel centre + o2::itsmft::SegmentationAlpide::detectorToLocalUnchecked(c.row, c.col, local); + } + o2::math_utils::Point3D global(0.f, 0.f, 0.f); + if (!mITSChipCentres.empty()) { // geometry loaded for the ITS road + global = o2::its::GeometryTGeo::Instance()->getMatrixL2G(c.chipID) * local; + } + itsLeg.push_back(c.leg); + itsFlags.push_back(c.flags); + itsChip.push_back(c.chipID); + itsRow.push_back(c.row); + itsCol.push_back(c.col); + itsPattID.push_back(c.pattID); + itsHasPatt.push_back(c.pattEntry >= 0); + itsRofBC.push_back(c.rofBC); + itsGx.push_back(global.X()); + itsGy.push_back(global.Y()); + itsGz.push_back(global.Z()); + itsXLoc.push_back(local.X()); + itsZLoc.push_back(local.Z()); + } + const auto& dbg = mDbgCosmics[icosm]; + const auto& time = cosm.cosmic.getTimeMUS(); + (*mDebugOut) << "cosmics" + << "tf=" << mDbgTF << "cosm=" << icosm << "t=" << time.getTimeStamp() << "tErr=" << time.getTimeStampError() << "absTime=" << int(absTimeKnown) + << "chi2Match=" << cosm.cosmic.getChi2Match() << "chi2Refit=" << cosm.cosmic.getChi2Refit() << "q2pt=" << cosm.cosmic.getQ2Pt() + << "tgl=" << cosm.cosmic.getTgl() << "side0=" << side[0] << "side1=" << side[1] << "nAtt0=" << nAttached[0] << "nAtt1=" << nAttached[1] + << "nRoad0=" << nRoad[0] << "nRoad1=" << nRoad[1] << "nITS=" << int(cosm.clITS.size()) << "nTOF=" << int(cosm.clTOF.size()) + << "nTRD=" << int(cosm.trdTracklets.size()) << "tTOF=" << cosm.timeTOFMUS << "scoreTOF=" << cosm.scoreTOFPair << "scoreTOFRev=" << cosm.scoreTOFReversed << "dupOf=" << cosm.duplicateOf << "polChi2Match=" << cosm.chi2MatchPolished << "polQ2Pt=" << cosm.polished.getQ2Pt() << "polZ=" << cosm.polished.getZ() + << "mcEvent=" << (cosm.label.isSet() ? cosm.label.getEventID() : -1) << "mcTrack=" << (cosm.label.isSet() ? cosm.label.getTrackID() : -1) + << "clLeg=" << clLeg << "clFlags=" << clFlags << "clSector=" << clSector << "clRow=" << clRow << "clPad=" << clPad << "clTime=" << clTime + << "clQMax=" << clQMax << "clQTot=" << clQTot << "clX=" << clX << "clY=" << clY << "clZ=" << clZ << "clZCos=" << clZCos << "clGx=" << clGx + << "clGy=" << clGy << "roadLeg=" << dbg.roadLeg << "roadFrame=" << dbg.roadFrame << "roadSector=" << dbg.roadSector << "roadRow=" << dbg.roadRow + << "roadX=" << dbg.roadX << "roadY=" << dbg.roadY << "roadZ=" << dbg.roadZ << "roadZCos=" << dbg.roadZCos << "roadGx=" << dbg.roadGx + << "roadGy=" << dbg.roadGy << "roadSnp=" << dbg.roadSnp << "roadTgl=" << dbg.roadTgl << "itsLeg=" << itsLeg << "itsFlags=" << itsFlags + << "itsChip=" << itsChip << "itsRow=" << itsRow << "itsCol=" << itsCol << "itsPattID=" << itsPattID << "itsHasPatt=" << itsHasPatt + << "itsRofBC=" << itsRofBC << "itsGx=" << itsGx << "itsGy=" << itsGy << "itsGz=" << itsGz << "itsXLoc=" << itsXLoc << "itsZLoc=" << itsZLoc + << "tofLeg=" << dbg.tofLeg << "tofChannel=" << dbg.tofChannel << "tofFlags=" << dbg.tofFlags << "tofTimeRaw=" << dbg.tofTimeRaw + << "tofTime=" << dbg.tofTime << "tofTot=" << dbg.tofTot << "tofX=" << dbg.tofX << "tofY=" << dbg.tofY << "tofZ=" << dbg.tofZ << "tofGx=" << dbg.tofGx + << "tofGy=" << dbg.tofGy << "trdLeg=" << dbg.trdLeg << "trdLayer=" << dbg.trdLayer << "trdX=" << dbg.trdX << "trdY=" << dbg.trdY << "trdZ=" << dbg.trdZ + << "trdDy=" << dbg.trdDy << "trdDz=" << dbg.trdDz << "trdTrigMUS=" << dbg.trdTrigMUS << "\n"; +} + +void CosmicsClusterCollectorSpec::writeDebugTOF(const o2::tof::Cluster& c, int icosm, int leg, uint8_t flags) const +{ + // the TOF cluster position is in the frame of its sector + const float alpha = o2::math_utils::sector2Angle(c.getSector()); + const float sinAlpha = std::sin(alpha); + const float cosAlpha = std::cos(alpha); + auto& dbg = mDbgCosmics[icosm]; + dbg.tofLeg.push_back(leg); + dbg.tofChannel.push_back(c.getMainContributingChannel()); + dbg.tofFlags.push_back(flags); + dbg.tofTimeRaw.push_back(c.getTimeRaw()); + dbg.tofTime.push_back(c.getTime()); + dbg.tofTot.push_back(c.getTot()); + dbg.tofX.push_back(c.getX()); + dbg.tofY.push_back(c.getY()); + dbg.tofZ.push_back(c.getZ()); + dbg.tofGx.push_back(c.getX() * cosAlpha - c.getY() * sinAlpha); + dbg.tofGy.push_back(c.getX() * sinAlpha + c.getY() * cosAlpha); +} + +void CosmicsClusterCollectorSpec::finaliseCCDB(ConcreteDataMatcher& matcher, void* obj) +{ + if (o2::base::GRPGeomHelper::instance().finaliseCCDB(matcher, obj)) { + return; + } + if (matcher == ConcreteDataMatcher("ITS", "CLUSDICT", 0)) { + mITSDict = (const o2::itsmft::TopologyDictionary*)obj; + return; + } +} + +void CosmicsClusterCollectorSpec::endOfStream(EndOfStreamContext& ec) +{ + mDebugOut.reset(); + LOGP(info, "Cosmics cluster collector: {} cosmics, {} attached and {} road TPC clusters; Cpu: {:.3e} Real: {:.3e} s in {} slots", + mNCosmics, mNClAttached, mNClCorridor, mTimer.CpuTime(), mTimer.RealTime(), mTimer.Counter() - 1); +} + +DataProcessorSpec getCosmicsClusterCollectorSpec(GTrackID::mask_t src, bool useMC, bool itsStag, DetID::mask_t roadDets) +{ + if (itsStag && roadDets[DetID::ITS]) { + LOGP(warning, "Staggered ITS clusters are not supported: ITS road of cosmics disabled"); + roadDets &= ~DetID::getMask(DetID::ITS); + } + auto dataRequest = std::make_shared(); + dataRequest->setITSPerLayer(itsStag); + dataRequest->requestTracks(src, false); + dataRequest->requestTracks(getLegITSSources(src), false); // the ITS clusters of legs with an ITS part are read through their ITS track + dataRequest->requestClusters(src, false); + if (roadDets[DetID::ITS]) { + dataRequest->requestITSClusters(false); + } + if (roadDets[DetID::TOF]) { + dataRequest->requestTOFClusters(false); + } + if (roadDets[DetID::TRD]) { + dataRequest->requestTRDTracklets(false); + } + dataRequest->requestCoscmicTracks(useMC); + auto ggRequest = std::make_shared(false, // orbitResetTime + false, // GRPECS + false, // GRPLHCIF + true, // GRPMagField + roadDets[DetID::TOF], // askMatLUT: polish refit of TOF-timed cosmics + roadDets[DetID::ITS] ? o2::base::GRPGeomRequest::Aligned : o2::base::GRPGeomRequest::None, // ITS road: chip positions + dataRequest->inputs, + true); + dataRequest->inputs.emplace_back("corrMap", o2::header::gDataOriginTPC, "TPCCORRMAP", 0, Lifetime::Timeframe); + + std::vector outputs; + outputs.emplace_back("GLO", "COSMFULL", 0, Lifetime::Timeframe); + outputs.emplace_back("GLO", "COSMFULLTF", 0, Lifetime::Timeframe); + outputs.emplace_back("GLO", "COSMFULLTFID", 0, Lifetime::Timeframe); + + return DataProcessorSpec{ + "cosmics-cluster-collector", + dataRequest->inputs, + outputs, + AlgorithmSpec{adaptFromTask(dataRequest, ggRequest, useMC, roadDets)}, + Options{ + {"corridor-width", VariantType::Float, 1.f, {"half-width of the road around each leg [cm]"}}, + {"max-abs-time-err", VariantType::Float, 0.5f, {"max. time error of a cosmic [mus] to search the other TPC side of its one-side legs"}}, + {"max-cosmics-per-tf", VariantType::Int, 100, {"write at most this many cosmics per TF (with their clusters), the others are dropped"}}, + {"tof-road-width", VariantType::Float, 5.f, {"half-width of the road at the TOF [cm]"}}, + {"tof-flight-tolerance", VariantType::Float, 2.f, {"max. deviation of the top/bottom TOF time difference from the muon's flight time [ns]"}}, + {"tof-time-error", VariantType::Float, 0.1f, {"error of a cosmic's TOF time for the TPC other-side corridor and the TRD / ITS roads [mus]"}}, + {"trd-road-width", VariantType::Float, 5.f, {"half-width of the road at the TRD [cm] (in z plus the pad length)"}}, + {"its-road-width", VariantType::Float, 1.5f, {"half-width of the road in the ITS [cm]"}}, + {"debug-tree", VariantType::Bool, false, {"write cosmics_collector_debug.root with transformed clusters and road points (test runs)"}}}}; +} + +} // namespace o2::globaltracking diff --git a/Detectors/GlobalTrackingWorkflow/src/CosmicsMatchingSpec.cxx b/Detectors/GlobalTrackingWorkflow/src/CosmicsMatchingSpec.cxx index 20941d8333e27..6f3cf0529e218 100644 --- a/Detectors/GlobalTrackingWorkflow/src/CosmicsMatchingSpec.cxx +++ b/Detectors/GlobalTrackingWorkflow/src/CosmicsMatchingSpec.cxx @@ -66,7 +66,7 @@ namespace globaltracking class CosmicsMatchingSpec : public Task { public: - CosmicsMatchingSpec(std::shared_ptr dr, std::shared_ptr gr, bool usePV, bool useMC) : mDataRequest(dr), mGGCCDBRequest(gr), mUsePVInfo(usePV), mUseMC(useMC) {} + CosmicsMatchingSpec(std::shared_ptr dr, std::shared_ptr gr, GTrackID::mask_t src, GTrackID::mask_t srcPVVeto, bool usePV, bool useMC) : mDataRequest(dr), mGGCCDBRequest(gr), mSeedSources(src), mPVVetoSources(srcPVVeto), mUsePVInfo(usePV), mUseMC(useMC) {} ~CosmicsMatchingSpec() override = default; void init(InitContext& ic) final; void run(ProcessingContext& pc) final; @@ -81,6 +81,8 @@ class CosmicsMatchingSpec : public Task o2::tpc::VDriftHelper mTPCVDriftHelper{}; const o2::gpu::TPCFastTransformPOD* mCorrMap{nullptr}; o2::globaltracking::MatchCosmics mMatching; // matching engine + GTrackID::mask_t mSeedSources; // track sources used as legs + GTrackID::mask_t mPVVetoSources; // sources whose PV contributors veto the legs they contain bool mUseMC = true; bool mUsePVInfo = false; TStopwatch mTimer; @@ -94,6 +96,8 @@ void CosmicsMatchingSpec::init(InitContext& ic) mMatching.setDebugFlag(ic.options().get("debug-tree-flags")); mMatching.setUseMC(mUseMC); mMatching.setUsePVInfo(mUsePVInfo); + mMatching.setSeedSources(mSeedSources); + mMatching.setPVVetoSources(mPVVetoSources); // } @@ -183,7 +187,27 @@ void CosmicsMatchingSpec::endOfStream(EndOfStreamContext& ec) mTimer.CpuTime(), mTimer.RealTime(), mTimer.Counter() - 1); } -DataProcessorSpec getCosmicsMatchingSpec(GTrackID::mask_t src, bool usePV, bool useMC, bool itsStag) +GTrackID::mask_t getLegITSSources(GTrackID::mask_t src) +{ + return GTrackID::includesDet(o2::detectors::DetID::ITS, src) ? GTrackID::getSourceMask(GTrackID::ITS) & ~src : GTrackID::mask_t{}; +} + +GTrackID::mask_t addPVContributorParents(GTrackID::mask_t srcPVContributors) +{ + auto src = srcPVContributors; + if (src[GTrackID::ITSTPCTRDTOF]) { + src |= GTrackID::getSourceMask(GTrackID::ITSTPCTRD); + } + if (src[GTrackID::TPCTRDTOF]) { + src |= GTrackID::getSourceMask(GTrackID::TPCTRD); + } + if (src[GTrackID::ITSTPCTRD] || src[GTrackID::ITSTPCTOF]) { + src |= GTrackID::getSourceMask(GTrackID::ITSTPC); + } + return src; +} + +DataProcessorSpec getCosmicsMatchingSpec(GTrackID::mask_t src, bool usePV, bool useMC, bool itsStag, bool useTOFClusters, GTrackID::mask_t srcPVContributors) { std::vector outputs; Options opts{ @@ -194,9 +218,15 @@ DataProcessorSpec getCosmicsMatchingSpec(GTrackID::mask_t src, bool usePV, bool dataRequest->setITSPerLayer(itsStag); dataRequest->requestTracks(src, useMC); + dataRequest->requestTracks(getLegITSSources(src), false); // the refit of legs with an ITS part reads their ITS track dataRequest->requestClusters(src, false); // no MC labels for clusters needed for refit only + if (useTOFClusters) { + dataRequest->requestTOFClusters(false); // TOF flight pairs of the candidate pairs (MatchCosmicsParams::tofFlightSelection) + } if (usePV) { dataRequest->requestPrimaryVertices(useMC); + // global tracks whose primary-vertex contributors veto their TPC track as a leg (MatchCosmicsParams::discardPVContributors), not legs + dataRequest->requestTracks(addPVContributorParents(srcPVContributors) & ~src, useMC); // same MC flag: some requests are shared with the legs } outputs.emplace_back("GLO", "COSMICTRC", 0, Lifetime::Timeframe); @@ -221,7 +251,7 @@ DataProcessorSpec getCosmicsMatchingSpec(GTrackID::mask_t src, bool usePV, bool "cosmics-matcher", dataRequest->inputs, outputs, - AlgorithmSpec{adaptFromTask(dataRequest, ggRequest, usePV, useMC)}, + AlgorithmSpec{adaptFromTask(dataRequest, ggRequest, src, usePV ? srcPVContributors : GTrackID::mask_t{}, usePV, useMC)}, opts}; } diff --git a/Detectors/GlobalTrackingWorkflow/src/TrackCosmicsWriterSpec.cxx b/Detectors/GlobalTrackingWorkflow/src/TrackCosmicsWriterSpec.cxx index 800978f7a4db3..bb7bbc5c371aa 100644 --- a/Detectors/GlobalTrackingWorkflow/src/TrackCosmicsWriterSpec.cxx +++ b/Detectors/GlobalTrackingWorkflow/src/TrackCosmicsWriterSpec.cxx @@ -16,6 +16,8 @@ #include "DPLUtils/MakeRootTreeWriterSpec.h" #include "ReconstructionDataFormats/TrackCosmics.h" #include "SimulationDataFormat/MCCompLabel.h" +#include "DataFormatsGlobalTracking/TrackCosmicsExtended.h" +#include "CommonDataFormat/TFIDInfo.h" using namespace o2::framework; @@ -51,5 +53,20 @@ DataProcessorSpec getTrackCosmicsWriterSpec(bool useMC) ""})(); } +DataProcessorSpec getCosmicsFullWriterSpec() +{ + auto logger = [](std::vector const& cosmics) { + LOG(info) << "Writing " << cosmics.size() << " cosmics with clusters"; + }; + return MakeRootTreeWriterSpec("cosmics-full-writer", + "o2_cosmics_full.root", + "cosmicsFull", + -1, // do not limit number of events to store + 100, // periodically autosave + BranchDefinition>{InputSpec{"cosmics", "GLO", "COSMFULL", 0}, "cosmics", 1, logger}, + BranchDefinition{InputSpec{"tfinfo", "GLO", "COSMFULLTF", 0}, "tfInfo", 1}, + BranchDefinition{InputSpec{"tfid", "GLO", "COSMFULLTFID", 0}, "tfID", 1})(); +} + } // namespace globaltracking } // namespace o2 diff --git a/Detectors/GlobalTrackingWorkflow/src/cosmics-match-workflow.cxx b/Detectors/GlobalTrackingWorkflow/src/cosmics-match-workflow.cxx index 67e2fd6cdc4cd..a14cad38fb5d2 100644 --- a/Detectors/GlobalTrackingWorkflow/src/cosmics-match-workflow.cxx +++ b/Detectors/GlobalTrackingWorkflow/src/cosmics-match-workflow.cxx @@ -26,7 +26,9 @@ #include "DetectorsCommonDataFormats/DetID.h" #include "GlobalTrackingWorkflowReaders/TrackTPCITSReaderSpec.h" #include "GlobalTrackingWorkflow/CosmicsMatchingSpec.h" +#include "GlobalTracking/MatchCosmicsParams.h" #include "GlobalTrackingWorkflow/TrackCosmicsWriterSpec.h" +#include "GlobalTrackingWorkflow/CosmicsClusterCollectorSpec.h" #include "Algorithm/RangeTokenizer.h" #include "DetectorsRaw/HBFUtilsInitializer.h" #include "Framework/CallbacksPolicy.h" @@ -52,6 +54,10 @@ void customize(std::vector& workflowOptions) {"disable-root-input", o2::framework::VariantType::Bool, false, {"disable root-files input reader"}}, {"disable-root-output", o2::framework::VariantType::Bool, false, {"disable root-files output writer"}}, {"use-pv-info", o2::framework::VariantType::Bool, false, {"request primary vertex for relevant cuts in the collision/cosmics interleaved data"}}, + {"pv-contributor-sources", VariantType::String, "", {"global track sources loaded only to reject the TPC track of a primary-vertex contributor as a leg (cosmicsMatch.discardPVContributors, implies --use-pv-info); sources without TPC are ignored, so the vertexing sources can be passed as they are"}}, + {"enable-cluster-output", o2::framework::VariantType::Bool, false, {"collect the raw clusters of the cosmics (legs + road around them) and write them with the cosmics"}}, + {"road-detectors", VariantType::String, "ITS,TOF,TRD", {"with --enable-cluster-output: detectors whose hits along the cosmic are collected besides the TPC road"}}, + {"cosmics-preset", VariantType::String, "", {"named set of cosmicsMatch settings applied before --configKeyValues (which can override single keys): physics-v1 = cosmics in collision data"}}, {"track-sources", VariantType::String, std::string{GID::ALL}, {"comma-separated list of sources to use"}}, {"configKeyValues", VariantType::String, "", {"Semicolon separated key=value strings ..."}}}; o2::itsmft::DPLAlpideParamInitializer::addITSConfigOption(options); @@ -70,6 +76,9 @@ void customize(std::vector& policies) policies.push_back(o2::tpc::TPCSectorCompletionPolicy("cosmics-matcher", o2::tpc::TPCSectorCompletionPolicy::Config::RequireAll, InputSpec{"cluster", o2::framework::ConcreteDataTypeMatcher{"TPC", "CLUSTERNATIVE"}})()); + policies.push_back(o2::tpc::TPCSectorCompletionPolicy("cosmics-cluster-collector", + o2::tpc::TPCSectorCompletionPolicy::Config::RequireAll, + InputSpec{"cluster", o2::framework::ConcreteDataTypeMatcher{"TPC", "CLUSTERNATIVE"}})()); policies.push_back(CompletionPolicyHelpers::consumeWhenAllOrdered(".*cosm.*[W,w]riter.*")); } @@ -82,7 +91,11 @@ WorkflowSpec defineDataProcessing(ConfigContext const& configcontext) WorkflowSpec specs; GID::mask_t alowedSources = GID::getSourcesMask("ITS,TPC,ITS-TPC,TPC-TRD,TPC-TOF,TPC-TRD-TOF,ITS-TPC-TOF,ITS-TPC-TRD-TOF"); - // Update the (declared) parameters if changed from the command line + // Update the (declared) parameters if changed from the command line: first an eventual preset, then the explicit key=values + auto preset = configcontext.options().get("cosmics-preset"); + if (!preset.empty()) { + o2::conf::ConfigurableParam::updateFromString(o2::globaltracking::getMatchCosmicsPreset(preset)); + } o2::conf::ConfigurableParam::updateFromString(configcontext.options().get("configKeyValues")); // write the configuration used for the workflow o2::conf::ConfigurableParam::writeINI("o2match-cosmics-workflow_configuration.ini"); @@ -105,21 +118,51 @@ WorkflowSpec defineDataProcessing(ConfigContext const& configcontext) src = src | GID::getSourcesMask("CTP"); } - GID::mask_t srcCl = src; + // legs with an ITS part need the ITS tracks and clusters for their refit, also if standalone ITS tracks are no legs (not added to src) + const GID::mask_t srcITS = o2::globaltracking::getLegITSSources(src); + GID::mask_t srcCl = src | srcITS; + const bool useTOFClusters = o2::globaltracking::MatchCosmicsParams::Instance().tofFlightSelection; + if (useTOFClusters) { + srcCl |= GID::getSourceMask(GID::TOF); + } GID::mask_t dummy; if (!configcontext.options().get("disable-root-input")) { specs.emplace_back(o2::tpc::getTPCScalerSpec(sclOpt)); } bool usePV = configcontext.options().get("use-pv-info"); - specs.emplace_back(o2::globaltracking::getCosmicsMatchingSpec(src, usePV, useMC, doStag)); + GID::mask_t srcPV = GID::getSourcesMask("ITS-TPC,TPC-TRD,TPC-TOF,ITS-TPC-TRD,TPC-TRD-TOF,ITS-TPC-TOF,ITS-TPC-TRD-TOF") & + GID::getSourcesMask(configcontext.options().get("pv-contributor-sources")); + if (srcPV.any() && !o2::globaltracking::MatchCosmicsParams::Instance().discardPVContributors) { + LOG(warning) << "--pv-contributor-sources ignored: cosmicsMatch.discardPVContributors is off"; + srcPV.reset(); + } + usePV |= srcPV.any(); + specs.emplace_back(o2::globaltracking::getCosmicsMatchingSpec(src, usePV, useMC, doStag, useTOFClusters, srcPV)); + const auto srcPVLoaded = o2::globaltracking::addPVContributorParents(srcPV); // with the parents needed to resolve the contributors + bool clusterOutput = configcontext.options().get("enable-cluster-output"); + if (clusterOutput) { + if (!src[GID::TPC]) { + LOG(fatal) << "--enable-cluster-output needs TPC tracks in --track-sources"; + } + const auto roadDets = DetID::getMask(configcontext.options().get("road-detectors")) & DetID::getMask("ITS,TOF,TRD"); + for (auto det : {DetID::ITS, DetID::TOF, DetID::TRD}) { + if (roadDets[det]) { + srcCl |= GID::getSourceMask(det == DetID::ITS ? GID::ITS : (det == DetID::TOF ? GID::TOF : GID::TRD)); + } + } + specs.emplace_back(o2::globaltracking::getCosmicsClusterCollectorSpec(src, useMC, doStag, roadDets)); + } - o2::globaltracking::InputHelper::addInputSpecs(configcontext, specs, src, src, src, useMC, dummy); // clusters MC is not needed + o2::globaltracking::InputHelper::addInputSpecs(configcontext, specs, srcCl, src | srcPVLoaded | srcITS, src | srcPVLoaded | srcITS, useMC, dummy); // clusters MC is not needed if (usePV) { o2::globaltracking::InputHelper::addInputSpecsPVertex(configcontext, specs, useMC); // P-vertex is always needed } if (!disableRootOut) { specs.emplace_back(o2::globaltracking::getTrackCosmicsWriterSpec(useMC)); + if (clusterOutput) { + specs.emplace_back(o2::globaltracking::getCosmicsFullWriterSpec()); + } } // configure dpl timer to inject correct firstTForbit: start from the 1st orbit of TF containing 1st sampled orbit diff --git a/GPU/GPUTracking/Interface/GPUO2InterfaceRefit.cxx b/GPU/GPUTracking/Interface/GPUO2InterfaceRefit.cxx index cd184b3820533..aa9a036c4fb41 100644 --- a/GPU/GPUTracking/Interface/GPUO2InterfaceRefit.cxx +++ b/GPU/GPUTracking/Interface/GPUO2InterfaceRefit.cxx @@ -136,7 +136,7 @@ void GPUO2InterfaceRefit::updateCalib(const TPCFastTransformPOD* trans, float bz int32_t GPUO2InterfaceRefit::RefitTrackAsGPU(o2::tpc::TrackTPC& trk, bool outward, bool resetCov) { return mRefit->RefitTrackAsGPU(trk, outward, resetCov); } int32_t GPUO2InterfaceRefit::RefitTrackAsTrackParCov(o2::tpc::TrackTPC& trk, bool outward, bool resetCov) { return mRefit->RefitTrackAsTrackParCov(trk, outward, resetCov); } int32_t GPUO2InterfaceRefit::RefitTrackAsGPU(o2::track::TrackParCov& trk, const o2::tpc::TrackTPCClusRef& clusRef, float time0, float* chi2, bool outward, bool resetCov) { return mRefit->RefitTrackAsGPU(trk, clusRef, time0, chi2, outward, resetCov); } -int32_t GPUO2InterfaceRefit::RefitTrackAsTrackParCov(o2::track::TrackParCov& trk, const o2::tpc::TrackTPCClusRef& clusRef, float time0, float* chi2, bool outward, bool resetCov) { return mRefit->RefitTrackAsTrackParCov(trk, clusRef, time0, chi2, outward, resetCov); } +int32_t GPUO2InterfaceRefit::RefitTrackAsTrackParCov(o2::track::TrackParCov& trk, const o2::tpc::TrackTPCClusRef& clusRef, float time0, float* chi2, bool outward, bool resetCov, int32_t eLossSign) { return mRefit->RefitTrackAsTrackParCov(trk, clusRef, time0, chi2, outward, resetCov, eLossSign); } void GPUO2InterfaceRefit::setIgnoreErrorsAtTrackEnds(bool v) { mRefit->mIgnoreErrorsOnTrackEnds = v; } void GPUO2InterfaceRefit::setTrackReferenceX(float v) { mParam->rec.tpc.trackReferenceX = v; } diff --git a/GPU/GPUTracking/Interface/GPUO2InterfaceRefit.h b/GPU/GPUTracking/Interface/GPUO2InterfaceRefit.h index f85a376b9185a..052f728821cac 100644 --- a/GPU/GPUTracking/Interface/GPUO2InterfaceRefit.h +++ b/GPU/GPUTracking/Interface/GPUO2InterfaceRefit.h @@ -67,7 +67,8 @@ class GPUO2InterfaceRefit int32_t RefitTrackAsGPU(o2::tpc::TrackTPC& trk, bool outward = false, bool resetCov = false); int32_t RefitTrackAsTrackParCov(o2::tpc::TrackTPC& trk, bool outward = false, bool resetCov = false); int32_t RefitTrackAsGPU(o2::track::TrackParCov& trk, const o2::tpc::TrackTPCClusRef& clusRef, float time0, float* chi2 = nullptr, bool outward = false, bool resetCov = false); - int32_t RefitTrackAsTrackParCov(o2::track::TrackParCov& trk, const o2::tpc::TrackTPCClusRef& clusRef, float time0, float* chi2 = nullptr, bool outward = false, bool resetCov = false); + // eLossSign: energy-loss sign of the propagations, 0 along the track parameters (loss forward, gain backward), +1 gain, -1 loss + int32_t RefitTrackAsTrackParCov(o2::track::TrackParCov& trk, const o2::tpc::TrackTPCClusRef& clusRef, float time0, float* chi2 = nullptr, bool outward = false, bool resetCov = false, int32_t eLossSign = 0); void setTrackReferenceX(float v); void setIgnoreErrorsAtTrackEnds(bool v); void updateCalib(const o2::gpu::TPCFastTransformPOD* trans, float bzNominalGPU); diff --git a/GPU/GPUTracking/Refit/GPUTrackingRefit.cxx b/GPU/GPUTracking/Refit/GPUTrackingRefit.cxx index 9e86e4627a2fb..c8d6b00d916cd 100644 --- a/GPU/GPUTracking/Refit/GPUTrackingRefit.cxx +++ b/GPU/GPUTracking/Refit/GPUTrackingRefit.cxx @@ -221,6 +221,7 @@ GPUd() int32_t GPUTrackingRefit::RefitTrack(T& trkX, bool outward, bool resetCov convertTrack::propagator>(trk, trkX, prop, &TrackParCovChi2); int32_t begin = 0, count; float tOffset; + [[maybe_unused]] int32_t eLossSign = 0; // energy-loss sign of the TrackParCov propagations (0: from the direction, outward = loss) if constexpr (std::is_same_v) { count = trkX.NClusters(); tOffset = trkX.GetParam().GetTOffset(); @@ -230,6 +231,7 @@ GPUd() int32_t GPUTrackingRefit::RefitTrack(T& trkX, bool outward, bool resetCov } else if constexpr (std::is_same_v) { count = trkX.clusRef.getEntries(); tOffset = trkX.time0; + eLossSign = trkX.eLossSign; } else { static_assert("Invalid template"); } @@ -357,7 +359,7 @@ GPUd() int32_t GPUTrackingRefit::RefitTrack(T& trkX, bool outward, bool resetCov IgnoreErrors(trk.getSnp()); return -1; } - if (!prop->PropagateToXBxByBz(trk, x, constants::MAX_SIN_PHI_LOW)) { + if (!prop->PropagateToXBxByBz(trk, x, constants::MAX_SIN_PHI_LOW, Propagator::MAX_STEP, Propagator::MatCorrType::USEMatCorrLUT, nullptr, eLossSign)) { IgnoreErrors(trk.getSnp()); return -2; } @@ -401,11 +403,13 @@ GPUd() int32_t GPUTrackingRefit::RefitTrack(T& trkX, bool outward, bool resetCov constexpr float kDeg2Rad = M_PI / 180.f; constexpr float kSectAngle = 2 * M_PI / 18.f; if (mPparam->rec.tpc.trackReferenceX <= 500) { - if (prop->PropagateToXBxByBz(trk, mPparam->rec.tpc.trackReferenceX)) { + // a forced energy-loss sign holds along the refit direction: the way to the reference X can run against it + const int32_t eLossSignRef = (mPparam->rec.tpc.trackReferenceX > trk.getX()) == outward ? eLossSign : -eLossSign; + if (prop->PropagateToXBxByBz(trk, mPparam->rec.tpc.trackReferenceX, Propagator::MAX_SIN_PHI, Propagator::MAX_STEP, Propagator::MatCorrType::USEMatCorrLUT, nullptr, eLossSignRef)) { if (CAMath::Abs(trk.getY()) > trk.getX() * CAMath::Tan(kSectAngle / 2.f)) { float newAlpha = trk.getAlpha() + CAMath::Round(CAMath::ATan2(trk.getY(), trk.getX()) / kDeg2Rad / 20.f) * kSectAngle; GPUTPCGMTrackParam::NormalizeAlpha(newAlpha); - trk.rotate(newAlpha) && prop->PropagateToXBxByBz(trk, mPparam->rec.tpc.trackReferenceX); + trk.rotate(newAlpha) && prop->PropagateToXBxByBz(trk, mPparam->rec.tpc.trackReferenceX, Propagator::MAX_SIN_PHI, Propagator::MAX_STEP, Propagator::MatCorrType::USEMatCorrLUT, nullptr, eLossSignRef); } } } diff --git a/GPU/GPUTracking/Refit/GPUTrackingRefit.h b/GPU/GPUTracking/Refit/GPUTrackingRefit.h index 70c9fd47d90f6..71f8f11bee3eb 100644 --- a/GPU/GPUTracking/Refit/GPUTrackingRefit.h +++ b/GPU/GPUTracking/Refit/GPUTrackingRefit.h @@ -73,15 +73,16 @@ class GPUTrackingRefit const o2::tpc::TrackTPCClusRef& clusRef; float time0; float* chi2; + int32_t eLossSign = 0; // energy-loss sign of the TrackParCov propagations: 0 along the track parameters (loss forward), +1 gain, -1 loss }; GPUd() int32_t RefitTrackAsGPU(o2::track::TrackParCov& trk, const o2::tpc::TrackTPCClusRef& clusRef, float time0, float* chi2 = nullptr, bool outward = false, bool resetCov = false) { TrackParCovWithArgs x{trk, clusRef, time0, chi2}; return RefitTrack(x, outward, resetCov); } - GPUd() int32_t RefitTrackAsTrackParCov(o2::track::TrackParCov& trk, const o2::tpc::TrackTPCClusRef& clusRef, float time0, float* chi2 = nullptr, bool outward = false, bool resetCov = false) + GPUd() int32_t RefitTrackAsTrackParCov(o2::track::TrackParCov& trk, const o2::tpc::TrackTPCClusRef& clusRef, float time0, float* chi2 = nullptr, bool outward = false, bool resetCov = false, int32_t eLossSign = 0) { - TrackParCovWithArgs x{trk, clusRef, time0, chi2}; + TrackParCovWithArgs x{trk, clusRef, time0, chi2, eLossSign}; return RefitTrack(x, outward, resetCov); } diff --git a/prodtests/full-system-test/calib-workflow.sh b/prodtests/full-system-test/calib-workflow.sh index a14ff3b620d45..a8b67e6939857 100644 --- a/prodtests/full-system-test/calib-workflow.sh +++ b/prodtests/full-system-test/calib-workflow.sh @@ -79,6 +79,23 @@ if [[ $CALIB_ASYNC_EXTRACTTIMESERIES == 1 ]] ; then CONFIG_TPCTIMESERIES+=" --mult-max ${TPCTIMESERIES_MULT_MAX}" add_W o2-tpc-time-series-workflow "$DISABLE_ROOT_INPUT ${CONFIG_TPCTIMESERIES}" fi +if [[ $CALIB_ASYNC_EXTRACTCOSMICS == 1 ]] ; then + # cosmic muons in collision data: TPC-only legs matched with the preset's selection, raw clusters of each cosmic -> o2_cosmics_full.root + : ${COSMICS_PRESET:=physics-v1} + COSMICS_ROAD_DETECTORS= + for det in ITS TOF TRD; do + has_detector_reco $det && COSMICS_ROAD_DETECTORS+="${COSMICS_ROAD_DETECTORS:+,}$det" + done + COSMICS_CONFIG= + has_detector_reco TOF || COSMICS_CONFIG="cosmicsMatch.tofFlightSelection=false" # the TOF flight selection of the preset needs TOF clusters + # the TPC part of a primary-vertex contributor comes from a collision: veto it as a leg, with the vertexing sources (those without TPC are ignored) + COSMICS_OPT= + : ${COSMICS_PV_SOURCES:=${VERTEXING_SOURCES:-}} + if [[ ${COSMICS_PV_VETO:-1} == 1 ]] && [[ $BEAMTYPE != "cosmic" ]] && has_detector_matching PRIMVTX && [[ -n ${VERTEXING_SOURCES:-} ]] && [[ -n $COSMICS_PV_SOURCES ]]; then + COSMICS_OPT+=" --pv-contributor-sources $COSMICS_PV_SOURCES" + fi + add_W o2-cosmics-match-workflow "$DISABLE_ROOT_INPUT $DISABLE_MC --track-sources TPC --cosmics-preset ${COSMICS_PRESET} --enable-cluster-output --road-detectors ${COSMICS_ROAD_DETECTORS:-none}$COSMICS_OPT" "$COSMICS_CONFIG" +fi # output-proxy for aggregator if workflow_has_parameter CALIB_PROXIES; then