diff --git a/DataFormats/Detectors/GlobalTracking/include/DataFormatsGlobalTracking/RecoContainer.h b/DataFormats/Detectors/GlobalTracking/include/DataFormatsGlobalTracking/RecoContainer.h index db072ef3a421e..5d91c63258e57 100644 --- a/DataFormats/Detectors/GlobalTracking/include/DataFormatsGlobalTracking/RecoContainer.h +++ b/DataFormats/Detectors/GlobalTracking/include/DataFormatsGlobalTracking/RecoContainer.h @@ -27,6 +27,7 @@ #include "SimulationDataFormat/MCTruthContainer.h" #include "SimulationDataFormat/ConstMCTruthContainer.h" #include "DataFormatsCTP/LumiInfo.h" +#include "DataFormatsITSMFT/ClusterID.h" #include #include @@ -188,8 +189,8 @@ namespace globaltracking { // max number of layers for which the ITS/MFT clusters, ROF records and patterns can be provided separately -constexpr int MaxITSLayers = 7; -constexpr int MaxMFTLayers = 10; +constexpr int MaxITSLayers = o2::itsmft::MaxITSClusLayers; +constexpr int MaxMFTLayers = o2::itsmft::MaxMFTClusLayers; // helper class to request DPL input data from the processor specs definition struct DataRequest { diff --git a/DataFormats/Detectors/ITSMFT/ITS/include/DataFormatsITS/TrackITS.h b/DataFormats/Detectors/ITSMFT/ITS/include/DataFormatsITS/TrackITS.h index 9b63509cc9424..2f026fce7a588 100644 --- a/DataFormats/Detectors/ITSMFT/ITS/include/DataFormatsITS/TrackITS.h +++ b/DataFormats/Detectors/ITSMFT/ITS/include/DataFormatsITS/TrackITS.h @@ -221,13 +221,6 @@ class TrackITSExt : public TrackITS GPUhdDefault() TrackITSExt(const TrackITSExt& t) = default; - void setClusterIndex(int l, int i) - { - int ncl = getNumberOfClusters(); - mIndex[ncl++] = (l << 28) + i; - getClusterRefs().setEntries(ncl); - } - GPUhdi() int getClusterIndex(int lr) const { return mIndex[lr]; } GPUh() int getFirstLayerClusterIndex() const diff --git a/DataFormats/Detectors/ITSMFT/common/include/DataFormatsITSMFT/ClusterID.h b/DataFormats/Detectors/ITSMFT/common/include/DataFormatsITSMFT/ClusterID.h new file mode 100644 index 0000000000000..a8d1da83533ce --- /dev/null +++ b/DataFormats/Detectors/ITSMFT/common/include/DataFormatsITSMFT/ClusterID.h @@ -0,0 +1,43 @@ +// Copyright 2019-2020 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 ClusterID.h +/// \brief Composition/decomposition of the ITS/MFT cluster IDs referring to per-layer cluster arrays +/// \author ruben.shahoyan@cern.ch + +#ifndef ALICEO2_ITSMFT_CLUSTERID_H +#define ALICEO2_ITSMFT_CLUSTERID_H + +namespace o2::itsmft +{ + +// max number of layers for which the ITS/MFT clusters, ROF records and patterns can be provided separately +constexpr int MaxITSClusLayers = 7; +constexpr int MaxMFTClusLayers = 10; +constexpr int MaxClusLayers = MaxITSClusLayers > MaxMFTClusLayers ? MaxITSClusLayers : MaxMFTClusLayers; + +///< With the per-layer (staggered readout) ITS/MFT clusters input the clusters are referred to by the +///< composed ID (layer << ClusLayerShift) + index_in_layer. With a single (monolithic) clusters input +///< all clusters sit in the layer slot 0, hence the composed ID coincides with the flat cluster index +///< and the same decoding works for both cases. +///< Note: the bit 31 is excluded from the layer field, since the negative values of the composed ID +///< are reserved for the "no cluster" flags. +constexpr int ClusLayerShift = 27; +constexpr int ClusIndexMask = (0x1 << ClusLayerShift) - 1; +static_assert((1 << (31 - ClusLayerShift)) >= MaxClusLayers, "ClusLayerShift leaves no room for the layer ID"); + +constexpr int composeClusID(int lr, int idx) { return (lr << ClusLayerShift) + idx; } +constexpr int clusID2Layer(int id) { return id >> ClusLayerShift; } +constexpr int clusID2Index(int id) { return id & ClusIndexMask; } + +} // namespace o2::itsmft + +#endif diff --git a/Detectors/GlobalTracking/include/GlobalTracking/MatchTPCITS.h b/Detectors/GlobalTracking/include/GlobalTracking/MatchTPCITS.h index e736f0c9c8a42..88d176d2f29bb 100644 --- a/Detectors/GlobalTracking/include/GlobalTracking/MatchTPCITS.h +++ b/Detectors/GlobalTracking/include/GlobalTracking/MatchTPCITS.h @@ -34,6 +34,7 @@ #include "ReconstructionDataFormats/GlobalTrackID.h" #include "MathUtils/Primitive2D.h" #include "CommonDataFormat/EvIndex.h" +#include #include "CommonDataFormat/InteractionRecord.h" #include "CommonDataFormat/RangeReference.h" #include "CommonDataFormat/BunchFilling.h" @@ -42,6 +43,7 @@ #include "CommonUtils/TreeStreamRedirector.h" #include "DataFormatsITSMFT/Cluster.h" #include "DataFormatsITSMFT/ROFRecord.h" +#include "DataFormatsITSMFT/DPLAlpideParam.h" #include "DataFormatsITS/TrackITS.h" #include "DataFormatsFT0/RecPoints.h" #include "FT0Reconstruction/InteractionTag.h" @@ -53,6 +55,7 @@ #include "GlobalTracking/MatchTPCITSParams.h" #include "DataFormatsITSMFT/TopologyDictionary.h" #include "DataFormatsITSMFT/TrkClusRef.h" +#include "DataFormatsITSMFT/ClusterID.h" #include "ITSMFTReconstruction/ChipMappingITS.h" #include "TPCFastTransformPOD.h" #if !defined(__CINT__) && !defined(__MAKECINT__) && !defined(__ROOTCLING__) && !defined(__CLING__) @@ -107,6 +110,9 @@ constexpr int MinusOne = -1; constexpr int MinusTen = -10; constexpr int Validated = -2; +///< per-layer status of ITS clusters (e.g. for the AfterBurner) +using ITSClusStatus = std::array, o2::its::RecoGeomHelper::getNLayers()>; + ///< flags to tell the status of TPC-ITS tracks comparison enum TrackRejFlag : int { Accept = 0, @@ -264,25 +270,25 @@ struct TPCABSeed { { return lowestLayer < o2::its::RecoGeomHelper::getNLayers() ? firstInLr[lowestLayer] : -1; } - bool checkLinkHasUsedClusters(int linkID, const std::vector& clStatus) const + bool checkLinkHasUsedClusters(int linkID, const ITSClusStatus& clStatus) const { // check if some clusters used by the link or its parents are forbidden (already used by validatet track) while (linkID > MinusOne) { const auto& link = getLink(linkID); - if (link.clID > MinusOne && clStatus[link.clID] != MinusOne) { + if (link.clID > MinusOne && clStatus[o2::itsmft::clusID2Layer(link.clID)][o2::itsmft::clusID2Index(link.clID)] != MinusOne) { return true; } linkID = link.parentID; } return false; } - void flagLinkUsedClusters(int linkID, std::vector& clStatus) const + void flagLinkUsedClusters(int linkID, ITSClusStatus& clStatus) const { // check if some clusters used by the link or its parents are forbidden (already used by validated track) while (linkID > MinusOne) { const auto& link = getLink(linkID); if (link.clID > MinusOne) { - clStatus[link.clID] = MinusTen; + clStatus[o2::itsmft::clusID2Layer(link.clID)][o2::itsmft::clusID2Index(link.clID)] = MinusTen; } linkID = link.parentID; } @@ -293,40 +299,59 @@ struct TPCABSeed { struct InteractionCandidate : public o2::InteractionRecord { o2::math_utils::Bracketf_t tBracket; // interaction time - int rofITS; // corresponding ITS ROF entry (in the ROFRecord vectors) + int rofITS; // corresponding ITS clock-layer cluster ROF entry (in the ROFRecord vectors) uint32_t flag; // origin, etc. o2::dataformats::RangeReference seedsRef; // references to AB seeds - InteractionCandidate(const o2::InteractionRecord& ir, float t, float dt, int rof, uint32_t f = 0) : o2::InteractionRecord(ir), tBracket(t - dt, t + dt), rofITS(rof), flag(f) {} -}; - -struct ITSChipClustersRefs { - ///< contaner for sorted cluster indices for certain time window (usually ROF) and reference on the start and N clusters - ///< for every chip - using ClusRange = o2::dataformats::RangeReference; - std::vector clusterID; // indices of sorted clusters - -#ifndef ENABLE_UPGRADES - std::array chipRefs; // offset and number of clusters in each chip - ITSChipClustersRefs(int nclIni = 50000) - { - clusterID.reserve(nclIni); - } -#else - std::vector chipRefs; // offset and number of clusters in each chip - ITSChipClustersRefs(int nchips = o2::its::RecoGeomHelper::getNChips(), int nclIni = 50000) + std::array rofLr; // per ITS layer: cluster ROF entry compatible with the candidate time, -1: none + uint8_t rofNextLr = 0; // bit lr set: the candidate time is compatible also with the ROF rofLr[lr]+1 of the layer + InteractionCandidate(const o2::InteractionRecord& ir, float t, float dt, int rof, uint32_t f = 0) : o2::InteractionRecord(ir), tBracket(t - dt, t + dt), rofITS(rof), flag(f) { - clusterID.reserve(nclIni); - chipRefs.resize(nchips, ClusRange()); + rofLr.fill(MinusOne); } -#endif +}; + +struct ABClusterInfo { + ///< compact info on an ITS cluster usable by the AfterBurner + float y = 0.f, z = 0.f; ///< Y, Z of the cluster in the tracking frame of its sensor + int id = MinusOne; ///< composed cluster ID, see o2::itsmft::composeClusID + int chip = -1; ///< global chip (sensor) ID +}; +struct ABLayerClusters { + ///< clusters of one ITS layer usable by the AfterBurner (not attached to ITS tracks), stored in per-ROF + ///< blocks sorted in (chip, Z). Only the blocks (ROFs) referenced by some interaction candidate are built. + using ClusRange = o2::dataformats::RangeReference; + std::vector rofTimes; ///< time brackets of all cluster ROFs of the layer + std::vector rofRefs; ///< block of every ROF in the clus vector, firstEntry < 0 if the block was not built + std::vector clus; ///< clusters of the built blocks + std::vector needROF; ///< scratch buffer: blocks to build for the current TF void clear() { - clusterID.clear(); - std::memset(chipRefs.data(), 0, chipRefs.size() * sizeof(ClusRange)); // reset chip->cluster references + rofTimes.clear(); + rofRefs.clear(); + clus.clear(); + needROF.clear(); } - size_t sizeInternal() const { return sizeof(int) * clusterID.size(); } - size_t capInternal() const { return sizeof(int) * clusterID.capacity(); } + size_t sizeInternal() const { return sizeof(ABClusterInfo) * clus.size() + sizeof(o2::math_utils::Bracketf_t) * rofTimes.size() + sizeof(ClusRange) * rofRefs.size() + needROF.size(); } + size_t capInternal() const { return sizeof(ABClusterInfo) * clus.capacity() + sizeof(o2::math_utils::Bracketf_t) * rofTimes.capacity() + sizeof(ClusRange) * rofRefs.capacity() + needROF.capacity(); } +}; + +struct ABLayerView { + ///< thread-local view of the ABLayerClusters block(s) covering the time of the interaction candidate being processed + using ClusRange = o2::dataformats::RangeReference; + int rofKey = -2; ///< 1st ROF of the loaded block(s): -1: no compatible ROF, -2: nothing loaded yet + int nROFsKey = 0; ///< number of consecutive ROFs loaded (1 or 2) + int chipOffs = 0; ///< global ID of the 1st chip of the layer + int nData = 0; ///< number of clusters in the view + const ABClusterInfo* data = nullptr; ///< clusters of the loaded block(s), sorted in (chip, Z) + std::vector chipRefs; ///< cluster range in data for every chip of the layer + std::vector touchedChips; ///< chips with non-empty chipRefs, for fast reset + std::vector merged; ///< buffer to merge 2 blocks when the candidate time is compatible with 2 ROFs +}; + +struct ABThreadClusterViews { + ///< per-thread views of the AB clusters of all layers + std::array layers; }; class MatchTPCITS @@ -340,6 +365,7 @@ class MatchTPCITS using Params = o2::globaltracking::MatchTPCITSParams; using MatCorrType = o2::base::Propagator::MatCorrType; using VDTriplet = o2::dataformats::Triplet; + using AlpParamITS = o2::itsmft::DPLAlpideParam; MatchTPCITS(); // std::unique_ptr to forward declared type needs constructor / destructor in .cxx ~MatchTPCITS(); @@ -348,6 +374,7 @@ class MatchTPCITS static constexpr int MaxLadderCand = 2 * MaxUpDnLadders + 1; // max ladders to check for matching clusters static constexpr int MaxSeedsPerLayer = 50; // TODO static constexpr int NITSLayers = o2::its::RecoGeomHelper::getNLayers(); + static_assert(NITSLayers == AlpParamITS::getNLayers(), "ITS layers count mismatch between geometry helper and DPLAlpideParam"); ///< perform matching for provided input #if !defined(__CINT__) && !defined(__MAKECINT__) && !defined(__ROOTCLING__) && !defined(__CLING__) void run(const o2::globaltracking::RecoContainer& inp, @@ -394,7 +421,13 @@ class MatchTPCITS void setNHBPerTF(int n) { mNHBPerTF = n; } ///< ITS readout mode - void setITSTriggered(bool v) { mITSTriggered = v; } + void setITSTriggered(bool v) + { + mITSTriggered = v; + if (mAlpParams) { // the ROF length depends on the readout mode + setAlpideParam(mAlpParams); + } + } bool isITSTriggered() const { return mITSTriggered; } void setUseFT0(bool v) { mUseFT0 = v; } @@ -403,14 +436,29 @@ class MatchTPCITS void setUseBCFilling(bool v) { mUseBCFilling = v; } bool getUseBCFilling() const { return mUseBCFilling; } - ///< set ITS ROFrame duration in microseconds - void setITSROFrameLengthMUS(float fums); - ///< set ITS ROFrame duration in BC (continuous mode only) - void setITSROFrameLengthInBC(int nbc); - - void setITSTimeBiasInBC(int n); - int getITSTimeBiasInBC() const { return mITSTimeBiasInBC; } - float getITSTimeBiasMUS() const { return mITSTimeBiasMUS; } + ///< set ITS Alpide parameters: all per-layer ITS ROF lengths, biases and the clock layer + ///< are derived from them, no other ITS timing setter is needed. Can be called in any order + ///< wrt setITSTriggered(). + void setAlpideParam(const AlpParamITS* p); + const AlpParamITS* getAlpideParam() const { return mAlpParams; } + + ///< layer whose ROF defines the granularity of the ITS tracks ROFRecords (0 if all layers share the same ROF) + int getITSClockLayer() const { return mITSClockLayer; } + + ///< per-layer ITS ROF length and bias, always defined for all NITSLayers + int getITSROFrameLengthInBC(int lr) const { return mITSROFrameLengthInBC[lr]; } + float getITSROFrameLengthMUS(int lr) const { return mITSROFrameLengthMUS[lr]; } + float getITSROFrameLengthMUSInv(int lr) const { return mITSROFrameLengthMUSInv[lr]; } + float getITSTimeResMUS(int lr) const { return mITSTimeResMUS[lr]; } + int getITSTimeBiasInBC(int lr) const { return mITSTimeBiasInBC[lr]; } + float getITSTimeBiasMUS(int lr) const { return mITSTimeBiasMUS[lr]; } + + ///< time bracket (wrt TF start) of the ROF starting at nBC of the given layer + BracketF getITSROFTimeBracket(long nBC, int lr) const + { + float tMin = (nBC + mITSTimeBiasInBC[lr]) * o2::constants::lhc::LHCBunchSpacingMUS; + return {tMin, tMin + mITSROFrameLengthMUS[lr]}; + } // ==================== >> DPL-driven input >> ======================= void setITSDictionary(const o2::itsmft::TopologyDictionary* d) { mITSDict = d; } @@ -493,9 +541,8 @@ class MatchTPCITS int prepareTPCTracksAfterBurner(); int addTPCSeed(const o2::track::TrackParCov& _tr, float t0, float terr, o2::dataformats::GlobalTrackID srcGID, int tpcID); - int preselectChipClusters(std::vector& clVecOut, const ClusRange& clRange, const ITSChipClustersRefs& itsChipClRefs, + int preselectChipClusters(std::vector& clVecOut, const ClusRange& clRange, const ABLayerView& clView, float trackY, float trackZ, float tolerY, float tolerZ) const; - void fillClustersForAfterBurner(int rofStart, int nROFs, ITSChipClustersRefs& itsChipClRefs); void flagUsedITSClusters(const o2::its::TrackITS& track); void doMatching(int sec); @@ -533,7 +580,7 @@ class MatchTPCITS ///< convert time to ITS ROFrame units in case of continuous ITS readout int time2ITSROFrameCont(float t) const { - int rof = (t - mITSTimeBiasMUS) * mITSROFrameLengthMUSInv; + int rof = (t - mITSTimeBiasMUS[mITSClockLayer]) * mITSROFrameLengthMUSInv[mITSClockLayer]; if (rof < 0) { rof = 0; } @@ -544,7 +591,7 @@ class MatchTPCITS ///< convert time to ITS ROFrame units in case of triggered ITS readout int time2ITSROFrameTrig(float t, int start) const { - t -= mITSTimeBiasMUS; + t -= mITSTimeBiasMUS[mITSClockLayer]; while (start < int(mITSROFTimes.size())) { if (mITSROFTimes[start].getMax() > t) { return start; @@ -575,10 +622,17 @@ class MatchTPCITS return delta > toler ? rejFlag : (delta < -toler ? -rejFlag : Accept); } + const ITSCluster& getITSCluster(int composedID) const + { + return mITSClustersArray[o2::itsmft::clusID2Layer(composedID)][o2::itsmft::clusID2Index(composedID)]; + } + // ========================= AFTERBURNER ========================= int prepareABSeeds(); - void processABSeed(int sid, const ITSChipClustersRefs& itsChipClRefs, uint8_t tID); - int followABSeed(const o2::track::TrackParCov& seed, const ITSChipClustersRefs& itsChipClRefs, int seedID, int lrID, TPCABSeed& ABSeed); + void prepareABClusters(); + void updateABLayerView(ABLayerView& view, int lr, int rof, int nROFs) const; + void processABSeed(int sid, const ABThreadClusterViews& itsClViews, uint8_t tID); + int followABSeed(const o2::track::TrackParCov& seed, const ABLayerView& clView, int seedID, int lrID, TPCABSeed& ABSeed); int registerABTrackLink(TPCABSeed& ABSeed, const o2::track::TrackParCov& trc, int clID, int parentID, int lr, int laddID, float chi2Cl); bool isBetter(float chi2A, float chi2B) { return chi2A < chi2B; } // RS FIMXE TODO void accountForOverlapsAB(int lrSeed); @@ -600,6 +654,7 @@ class MatchTPCITS ///========== Parameters to be set externally, e.g. from CCDB ==================== const Params* mParams = nullptr; const o2::ft0::InteractionTag* mFT0Params = nullptr; + const AlpParamITS* mAlpParams = nullptr; ///< ITS Alpide parameters, set externally MatCorrType mUseMatCorrFlag = MatCorrType::USEMatCorrTGeo; bool mUseBCFilling = false; ///< use BC filling for candidates validation @@ -618,12 +673,15 @@ class MatchTPCITS ///< assigned time0 and its track Z position (converted from mTPCTimeEdgeZSafeMargin) float mTPCTimeEdgeTSafeMargin = 0.f; float mTPCExtConstrainedNSigmaInv = 0.f; // inverse for NSigmas for TPC time-interval from external constraint time sigma - int mITSROFrameLengthInBC = 0; ///< ITS RO frame in BC (for ITS cont. mode only) - float mITSROFrameLengthMUS = -1.; ///< ITS RO frame in \mus - float mITSTimeResMUS = -1.; ///< nominal ITS time resolution derived from ROF - float mITSROFrameLengthMUSInv = -1.; ///< ITS RO frame in \mus inverse - int mITSTimeBiasInBC = 0; ///< ITS RO frame shift in BCs, i.e. t_i = (I_ROF*mITSROFrameLengthInBC + mITSTimeBiasInBC)*BCLength_MUS - float mITSTimeBiasMUS = 0.; ///< ITS RO frame shift in \mus, i.e. t_i = (I_ROF*mITSROFrameLengthInBC)*BCLength_MUS + mITSTimeBiasMUS + ///< ITS ROF timings, always filled for all NITSLayers, also when all layers share the same ROF. + ///< t_i(lr) = (I_ROF*mITSROFrameLengthInBC[lr] + mITSTimeBiasInBC[lr])*BCLength_MUS + int mITSClockLayer = 0; ///< layer defining the ITS tracks ROFRecords granularity + std::array mITSROFrameLengthInBC{}; ///< ITS RO frame in BC per layer (cont. mode only) + std::array mITSROFrameLengthMUS{}; ///< ITS RO frame in \mus per layer + std::array mITSROFrameLengthMUSInv{}; ///< inverse ITS RO frame in \mus per layer + std::array mITSTimeResMUS{}; ///< nominal ITS time resolution derived from per-layer ROF + std::array mITSTimeBiasInBC{}; ///< ITS RO frame shift in BC per layer + std::array mITSTimeBiasMUS{}; ///< ITS RO frame shift in \mus per layer float mTPCVDrift = -1.; ///< TPC drift speed in cm/microseconds float mTPCVDriftInv = -1.; ///< inverse TPC nominal drift speed in cm/microseconds float mTPCDriftTimeOffset = 0; ///< drift time offset in mus @@ -662,15 +720,14 @@ class MatchTPCITS gsl::span mITSTrackROFRec; ///< input ITS tracks ROFRecord span gsl::span mITSTracksArray; ///< input ITS tracks span gsl::span mITSTrackClusIdx; ///< input ITS track cluster indices span - std::vector mITSClustersArray; ///< ITS clusters created in loadInput - std::vector mITSClusterSizes; ///< ITS cluster sizes created in loadInput - - gsl::span mITSClusterROFRec; ///< input ITS clusters ROFRecord span + std::array, NITSLayers> mITSClustersArray{}; ///< ITS clusters created in loadInput + std::array, NITSLayers> mITSClusterSizes{}; ///< ITS cluster sizes created in loadInput + std::array, NITSLayers> mITSClusterROFRec; ///< input ITS clusters ROFRecord span + int mNITSClusters = 0; gsl::span mFITInfo; ///< optional input FIT info span gsl::span mTPCRefitterShMap; ///< externally set TPC clusters sharing map gsl::span mTPCRefitterOccMap; ///< externally set TPC clusters occupancy map - const o2::itsmft::TopologyDictionary* mITSDict{nullptr}; // cluster patterns dictionary #ifdef ENABLE_UPGRADES const o2::its3::TopologyDictionary* mIT3Dict{nullptr}; // cluster patterns dictionary @@ -678,7 +735,7 @@ class MatchTPCITS const o2::tpc::ClusterNativeAccess* mTPCClusterIdxStruct = nullptr; ///< struct holding the TPC cluster indices - const o2::dataformats::MCTruthContainer* mITSClsLabels = nullptr; ///< input ITS Cluster MC labels + std::array*, NITSLayers> mITSClsLabels{}; ///< input ITS Cluster MC labels gsl::span mITSTrkLabels; ///< input ITS Track MC labels gsl::span mTPCTrkLabels; ///< input TPC Track MC labels /// <<<----- @@ -711,7 +768,13 @@ class MatchTPCITS ///< indices of selected track entries in mTPCWork (for tracks selected by AfterBurner) std::vector mTPCABIndexCache; std::vector mABWinnersIDs; - std::vector mABClusterLinkIndex; ///< index of 1st ABClusterLink for every cluster used by AfterBurner, -1: unused, -10: used by external ITS tracks + ///< per storage-slot status of ITS clusters wrt AfterBurner: -1: free, -10: used by an external ITS track or a validated AB track. + ///< The slot is the layer for per-layer clusters input, 0 for the monolithic one; only the slots used by the AfterBurner are booked + ITSClusStatus mABClusterStatus; + std::array mABLayerClusters; ///< per (layer, ROF) blocks of AB-usable clusters; filled for the AB layers and the ROF times also for the clock layer + std::array mABChipsBounds{}; ///< the layer lr owns the global chip IDs [mABChipsBounds[lr], mABChipsBounds[lr+1]) + float mABROFMarginMUS = 0.f; ///< effective margin for candidate time to ITS ROF matching: abROFMarginMUS clamped to below half of the shortest AB layer ROF + float mITSMaxROFOverhangMUS = 0.f; ///< max excess of ITS track time brackets over the end of their clock-layer ROF in the current TF LinksPoolMT mABLinksPool; ///< per sector indices of TPC track entry in mTPCWork diff --git a/Detectors/GlobalTracking/include/GlobalTracking/MatchTPCITSParams.h b/Detectors/GlobalTracking/include/GlobalTracking/MatchTPCITSParams.h index 3ec189deff54b..a45b8e135a20a 100644 --- a/Detectors/GlobalTracking/include/GlobalTracking/MatchTPCITSParams.h +++ b/Detectors/GlobalTracking/include/GlobalTracking/MatchTPCITSParams.h @@ -85,6 +85,8 @@ struct MatchTPCITSParams : public o2::conf::ConfigurableParamHelperrunAfterBurner) { // only used in AfterBurner mRGHelper.init(mParams->lowestLayerAB); // prepare helper for TPC track / ITS clusters matching } + { + const auto* geomITS = o2::its::GeometryTGeo::Instance(); + for (int lr = 0; lr < NITSLayers; lr++) { + mABChipsBounds[lr] = geomITS->getFirstChipIndex(lr); + } + mABChipsBounds[NITSLayers] = geomITS->getLastChipIndex(NITSLayers - 1) + 1; + } clear(); @@ -262,30 +276,53 @@ void MatchTPCITS::init() //______________________________________________ void MatchTPCITS::updateTimeDependentParams() { - ///< update parameters depending on time (once per TF) - auto& elParam = o2::tpc::ParameterElectronics::Instance(); - auto& detParam = o2::tpc::ParameterDetector::Instance(); - mTPCTBinMUS = elParam.ZbinWidth; - mTPCTBinNS = mTPCTBinMUS * 1e3; - mTPCZMax = detParam.TPClength; - mTPCTBinMUSInv = 1. / mTPCTBinMUS; - assert(mITSROFrameLengthMUS > 0.0f); + ///< update parameters depending on time (once per TF or beginning of run) + mBz = o2::base::Propagator::Instance()->getNominalBz(); + mFieldON = std::abs(mBz) > 0.01; + + static bool initOnceDone = false; + if (!initOnceDone) { + initOnceDone = true; + if (mParams->runAfterBurner) { + // margin for matching interaction candidate times to per-layer ITS ROFs in the AfterBurner: must stay + // below half of the shortest AB layer ROF length to guarantee at most 2 compatible ROFs per layer + float margin = std::max(0.f, mParams->abROFMarginMUS), marginMax = mITSROFrameLengthMUS[mParams->lowestLayerAB]; + for (int lr = mParams->lowestLayerAB + 1; lr < NITSLayers; lr++) { + if (mITSROFrameLengthMUS[lr] < marginMax) { + marginMax = mITSROFrameLengthMUS[lr]; + } + } + marginMax *= 0.45f; + if (margin > marginMax) { + if (mABROFMarginMUS != marginMax) { + LOGP(warn, "abROFMarginMUS={} exceeds 45% of the shortest AfterBurner layer ROF length, clamping to {}", margin, marginMax); + } + margin = marginMax; + } + mABROFMarginMUS = margin; + } + auto& elParam = o2::tpc::ParameterElectronics::Instance(); + auto& detParam = o2::tpc::ParameterDetector::Instance(); + mTPCTBinMUS = elParam.ZbinWidth; + mTPCTBinNS = mTPCTBinMUS * 1e3; + mTPCZMax = detParam.TPClength; + mTPCTBinMUSInv = 1. / mTPCTBinMUS; + assert(mITSROFrameLengthMUS[mITSClockLayer] > 0.0f); + + o2::math_utils::Point3D p0(90., 1., 1), p1(90., 100., 100.); + auto matbd = o2::base::Propagator::Instance()->getMatBudget(mParams->matCorr, p0, p1); + mTPCmeanX0Inv = matbd.meanX2X0 / matbd.length; + } mTPCBin2Z = mTPCTBinMUS * mTPCVDrift; mZ2TPCBin = 1. / mTPCBin2Z; mTPCVDriftInv = 1. / mTPCVDrift; mNTPCBinsFullDrift = mTPCZMax * mZ2TPCBin; mTPCTimeEdgeTSafeMargin = z2TPCBin(mParams->safeMarginTPCTimeEdge); mTPCExtConstrainedNSigmaInv = 1.f / mParams->tpcExtConstrainedNSigma; - mBz = o2::base::Propagator::Instance()->getNominalBz(); - mFieldON = std::abs(mBz) > 0.01; mMinTPCTrackPtInv = (mFieldON && mParams->minTPCTrackR > 0) ? 1. / std::abs(mParams->minTPCTrackR * mBz * o2::constants::math::B2C) : 999.; mMinITSTrackPtInv = (mFieldON && mParams->minITSTrackR > 0) ? 1. / std::abs(mParams->minITSTrackR * mBz * o2::constants::math::B2C) : 999.; - o2::math_utils::Point3D p0(90., 1., 1), p1(90., 100., 100.); - auto matbd = o2::base::Propagator::Instance()->getMatBudget(mParams->matCorr, p0, p1); - mTPCmeanX0Inv = matbd.meanX2X0 / matbd.length; - const auto& trackTune = TrackTuneParams::Instance(); float scale = mLumiCTP; if (scale < 0.f) { @@ -645,57 +682,83 @@ bool MatchTPCITS::prepareITSData() const auto& inp = *mRecoCont; // ITS clusters - mITSClusterROFRec = inp.getITSClustersROFRecords(); - const auto clusITS = inp.getITSClusters(); - if (mITSClusterROFRec.empty() || clusITS.empty()) { - LOG(info) << "No ITS clusters"; - return false; - } - const auto patterns = inp.getITSClustersPatterns(); - auto pattIt = patterns.begin(); - mITSClustersArray.reserve(clusITS.size()); + int nClLr = inp.getITSPerLayer() ? NITSLayers : 1; + mNITSClusters = 0; + for (int lr = 0; lr < nClLr; lr++) { + mITSClusterROFRec[lr] = inp.getITSClustersROFRecords(lr); + const auto clusITS = inp.getITSClusters(lr); + mNITSClusters += clusITS.size(); + const auto patterns = inp.getITSClustersPatterns(lr); + auto pattIt = patterns.begin(); + mITSClustersArray[lr].reserve(clusITS.size()); #ifdef ENABLE_UPGRADES - bool withITS3 = o2::GlobalParams::Instance().withITS3; - if (withITS3) { - o2::its3::ioutils::convertCompactClusters(clusITS, pattIt, mITSClustersArray, mIT3Dict); - } else { - o2::its::ioutils::convertCompactClusters(clusITS, pattIt, mITSClustersArray, mITSDict); - } + bool withITS3 = o2::GlobalParams::Instance().withITS3; + if (withITS3) { + o2::its3::ioutils::convertCompactClusters(clusITS, pattIt, mITSClustersArray[lr], mIT3Dict); + } else { + o2::its::ioutils::convertCompactClusters(clusITS, pattIt, mITSClustersArray[lr], mITSDict); + } #else - o2::its::ioutils::convertCompactClusters(clusITS, pattIt, mITSClustersArray, mITSDict); + o2::its::ioutils::convertCompactClusters(clusITS, pattIt, mITSClustersArray[lr], mITSDict); #endif - - // ITS clusters sizes - mITSClusterSizes.reserve(clusITS.size()); - auto pattIt2 = patterns.begin(); - for (auto& clus : clusITS) { - auto pattID = clus.getPatternID(); - unsigned int npix; + // ITS clusters sizes + mITSClusterSizes[lr].reserve(clusITS.size()); + auto pattIt2 = patterns.begin(); + for (auto& clus : clusITS) { + auto pattID = clus.getPatternID(); + unsigned int npix; #ifdef ENABLE_UPGRADES - auto ib = o2::its3::constants::detID::isDetITS3(clus.getChipID()); - if ((pattID == o2::itsmft::CompCluster::InvalidPatternID) || ((withITS3) ? mIT3Dict->isGroup(pattID, ib) : mITSDict->isGroup(pattID))) { // braces guarantee evaluation order + auto ib = o2::its3::constants::detID::isDetITS3(clus.getChipID()); + if ((pattID == o2::itsmft::CompCluster::InvalidPatternID) || ((withITS3) ? mIT3Dict->isGroup(pattID, ib) : mITSDict->isGroup(pattID))) { // braces guarantee evaluation order #else - if (pattID == o2::itsmft::CompCluster::InvalidPatternID || mITSDict->isGroup(pattID)) { + if (pattID == o2::itsmft::CompCluster::InvalidPatternID || mITSDict->isGroup(pattID)) { #endif - o2::itsmft::ClusterPattern patt; - patt.acquirePattern(pattIt2); - npix = patt.getNPixels(); - } else { -#ifdef ENABLE_UPGRADES - if (withITS3) { - npix = mIT3Dict->getNpixels(pattID, ib); + o2::itsmft::ClusterPattern patt; + patt.acquirePattern(pattIt2); + npix = patt.getNPixels(); } else { - npix = mITSDict->getNpixels(pattID); - } +#ifdef ENABLE_UPGRADES + if (withITS3) { + npix = mIT3Dict->getNpixels(pattID, ib); + } else { + npix = mITSDict->getNpixels(pattID); + } #else - npix = mITSDict->getNpixels(pattID); + npix = mITSDict->getNpixels(pattID); #endif + } + mITSClusterSizes[lr].push_back(std::clamp(npix, 0u, 255u)); + } + if (mMCTruthON) { + mITSClsLabels[lr] = inp.getITSClustersMCLabels(lr); + } + } + if (!mNITSClusters) { + LOGP(warn, "No ITS clusters"); + return false; + } + // time brackets of the cluster ROFs: the clock-layer ones relate interaction candidates to ITS ROFs, + // the AfterBurner layers ones select the clusters compatible with the candidate time + for (int lr = 0; lr < NITSLayers; lr++) { + auto& abClus = mABLayerClusters[lr]; + abClus.clear(); + if (lr != mITSClockLayer && !(mParams->runAfterBurner && lr >= mParams->lowestLayerAB)) { + continue; + } + const auto& rofs = mITSClusterROFRec[inp.getITSPerLayer() ? lr : 0]; + abClus.rofTimes.reserve(rofs.size()); + for (const auto& rofRec : rofs) { + abClus.rofTimes.push_back(getITSROFTimeBracket(rofRec.getBCData().differenceInBC(mStartIR), lr)); } - mITSClusterSizes.push_back(std::clamp(npix, 0u, 255u)); } - if (mMCTruthON) { - mITSClsLabels = inp.getITSClustersMCLabels(); + // the AfterBurner uses only layers from mParams->lowestLayerAB, with single ITS clusters input all of them are in the slot 0 + if (inp.getITSPerLayer()) { + for (int lr = mParams->lowestLayerAB; lr < NITSLayers; lr++) { + mABClusterStatus[lr].resize(mITSClustersArray[lr].size(), MinusOne); + } + } else { + mABClusterStatus[0].resize(mITSClustersArray[0].size(), MinusOne); } // ITS tracks @@ -708,10 +771,6 @@ bool MatchTPCITS::prepareITSData() int nROFs = mITSTrackROFRec.size(); mITSWork.reserve(mITSTracksArray.size()); - // total N ITS clusters in TF - const auto& lastClROF = mITSClusterROFRec.back(); - int nITSClus = lastClROF.getFirstEntry() + lastClROF.getNEntries(); - mABClusterLinkIndex.resize(nITSClus, MinusOne); for (int sec = o2::constants::math::NSectors; sec--;) { mITSTimeStart[sec].resize(nROFs, -1); // start of ITS work tracks in every sector } @@ -719,6 +778,8 @@ bool MatchTPCITS::prepareITSData() long maxBCs = nHBF * long(o2::constants::lhc::LHCMaxBunches); o2::track::TrackLTIntegral trackLTInt; trackLTInt.setTimeNotNeeded(); + mITSMaxROFOverhangMUS = 0.f; + const float trackTimeMarginBC = std::max(0.f, mParams->itsTimeStampMarginBC); for (int irof = 0; irof < nROFs; irof++) { const auto& rofRec = mITSTrackROFRec[irof]; @@ -730,17 +791,17 @@ bool MatchTPCITS::prepareITSData() } break; } - float tMin = nBC * o2::constants::lhc::LHCBunchSpacingMUS + mITSTimeBiasMUS; - float tMax = (nBC + mITSROFrameLengthInBC) * o2::constants::lhc::LHCBunchSpacingMUS + mITSTimeBiasMUS; + auto tBracket = getITSROFTimeBracket(nBC, mITSClockLayer); + float tMin = tBracket.getMin(), tMax = tBracket.getMax(); if (!mITSTriggered) { - size_t irofCont = nBC / mITSROFrameLengthInBC; + size_t irofCont = nBC / mITSROFrameLengthInBC[mITSClockLayer]; if (mITSTrackROFContMapping.size() <= irofCont) { // there might be gaps in the non-empty rofs, this will map continuous ROFs index to non empty ones mITSTrackROFContMapping.resize((1 + irofCont / 128) * 128, 0); } mITSTrackROFContMapping[irofCont] = irof; } - mITSROFTimes.emplace_back(tMin, tMax); // ITS ROF min/max time + mITSROFTimes.emplace_back(tMin, tMax); // nominal ITS ROF min/max time, to be extended to the envelope of the per-track time brackets for (int sec = o2::constants::math::NSectors; sec--;) { // start of sector's tracks for this ROF mITSTimeStart[sec][irof] = mITSSectIndexCache[sec].size(); // The sorting does not affect this @@ -758,9 +819,31 @@ bool MatchTPCITS::prepareITSData() if (std::abs(trcOrig.getQ2Pt()) > mMinITSTrackPtInv) { continue; } + // per-track time bracket from the tracker time stamp (bias-corrected BC since the TF start), + // widened by the itsTimeStampMarginBC safety margin; the nominal ROF bracket is used as a + // fallback when the time stamp is invalid (legacy input) + float tMinTrc = tMin, tMaxTrc = tMax; + const auto& tstamp = trcOrig.getTimeStamp(); + if (tstamp.getTimeStampError() > 0.f) { + // the tracker guarantees the raw lower edge of the time stamp to be within the assigned clock-layer ROF + assert(tstamp.getTimeStamp() - tstamp.getTimeStampError() >= nBC + mITSTimeBiasInBC[mITSClockLayer] - 0.5f); + float errBC = tstamp.getTimeStampError() + trackTimeMarginBC; + tMinTrc = (tstamp.getTimeStamp() - errBC) * o2::constants::lhc::LHCBunchSpacingMUS; + tMaxTrc = (tstamp.getTimeStamp() + errBC) * o2::constants::lhc::LHCBunchSpacingMUS; + auto& rofEnv = mITSROFTimes.back(); // extend the ROF envelope used by the TPC-side and triggered-mode entry caches + if (tMinTrc < rofEnv.getMin()) { + rofEnv.setMin(tMinTrc); + } + if (tMaxTrc > rofEnv.getMax()) { + rofEnv.setMax(tMaxTrc); + if (tMaxTrc - tMax > mITSMaxROFOverhangMUS) { + mITSMaxROFOverhangMUS = tMaxTrc - tMax; // max excess of the track brackets over their ROF end + } + } + } int nWorkTracks = mITSWork.size(); // working copy of outer track param - auto& trc = mITSWork.emplace_back(TrackLocITS{trcOrig.getParamOut(), {tMin, tMax}, it, irof, MinusOne}); + auto& trc = mITSWork.emplace_back(TrackLocITS{trcOrig.getParamOut(), {tMinTrc, tMaxTrc}, it, irof, MinusOne}); if (!trc.rotate(o2::math_utils::angle2Alpha(trc.getPhiPos()))) { mITSWork.pop_back(); // discard failed track continue; @@ -810,8 +893,10 @@ bool MatchTPCITS::prepareITSData() } } - // sort tracks in each sector according to their min time, then tgl - // RSTODO: sorting in tgl will be dangerous once the tracks with different time uncertaincies will be added + // Sort tracks in each sector according to their bracket min time (tgl serves only as a deterministic tie-break). + // Since the raw lower edge of every track time stamp is guaranteed to be within its clock-layer ROF and the + // safety margin shifts all tracks alike, the sorting cannot mix tracks of different ROFs, hence the + // mITSTimeStart entries assigned at the filling stage above remain valid. for (int sec = o2::constants::math::NSectors; sec--;) { auto& indexCache = mITSSectIndexCache[sec]; if (mParams->verbosity > 0) { @@ -834,7 +919,7 @@ bool MatchTPCITS::prepareITSData() mMatchRecordsITS.reserve(mITSWork.size() * mParams->maxMatchCandidates); mTimer[SWPrepITS].Stop(); - return nITSClus > 0; + return mNITSClusters > 0; } //_____________________________________________________ @@ -887,7 +972,9 @@ void MatchTPCITS::doMatching(int sec) // estimate ITS 1st ROframe bin this track may match to: TPC track are sorted according to their // timeMax, hence the timeMax - MaxmNTPCBinsFullDrift are non-decreasing auto tmn = trefTPC.tBracket.getMax() - maxTDriftSafe; - itsROBin = mITSTriggered ? time2ITSROFrameTrig(tmn, itsROBin) : time2ITSROFrameCont(tmn); + // in continuous mode the lookup time is decreased by the max excess of the ITS track brackets over their + // ROF end, since the mITSTimeStart binning is in ROF units while the brackets may extend beyond the ROF + itsROBin = mITSTriggered ? time2ITSROFrameTrig(tmn, itsROBin) : time2ITSROFrameCont(tmn - mITSMaxROFOverhangMUS); if (itsROBin >= int(timeStartITS.size())) { // time of TPC track exceeds the max time of ITS in the cache break; @@ -942,27 +1029,6 @@ void MatchTPCITS::doMatching(int sec) fillTPCITSmatchTree(cacheITS[iits], cacheTPC[itpc], rejFlag, chi2, timeCorr); } #endif - /* - // RS: this might be dangerous for ITS tracks with different time coverages. - if (rejFlag == RejectOnTgl) { - // ITS tracks in each ROFrame are ordered in Tgl, hence if this check failed on Tgl check - // (i.e. tgl_its>tgl_tpc+tolerance), then all other ITS tracks in this ROFrame will also have tgl too large. - // Jump on the 1st ITS track of the next ROFrame - int rof = trefITS.roFrame; - bool stop = false; - do { - if (++rof >= int(timeStartITS.size())) { - stop = true; - break; // no more ITS ROFrames in cache - } - iits = timeStartITS[rof] - 1; // next track to be checked -1 - } while (iits <= timeStartITS[trefITS.roFrame]); // skip empty bins - if (stop) { - break; - } - continue; - } - */ if (rejFlag != Accept) { continue; } @@ -1569,7 +1635,7 @@ void MatchTPCITS::fillCalibDebug(int ifit, int iTPC, const o2::dataformats::Trac (*mDBGOut) << "refit" << "multTPC=" << mltTPC << "multITSTr=" << mITSTrackROFRec[tITS.roFrame].getNEntries() - << "multITSCl=" << mITSClusterROFRec[tITS.roFrame].getNEntries() + << "multITSCl=" << mNITSClusters << "tf=" << mTFCount << "\n"; } #endif @@ -1599,8 +1665,9 @@ bool MatchTPCITS::refitTrackTPCITS(int slot, int iTPC, int& iITS, pmr::vector mITSTimeResMUS && tTPC.constraint != TrackLocTPC::Constrained) { - timeErr = mITSTimeResMUS; // chose smallest error + float itsTimeRes = tITS.tBracket.delta() / std::sqrt(12.f); // ITS time resolution from the track time-stamp bracket (uniform distribution) + if (timeErr > itsTimeRes && tTPC.constraint != TrackLocTPC::Constrained) { + timeErr = itsTimeRes; // chose smallest error deltaT = tTPC.constraint == TrackLocTPC::ASide ? tITS.tBracket.mean() - tTPC.time0 : tTPC.time0 - tITS.tBracket.mean(); } timeErr += mParams->globalTimeExtraErrorMUS; @@ -1638,7 +1705,7 @@ bool MatchTPCITS::refitTrackTPCITS(int slot, int iTPC, int& iITS, pmr::vectorgetSensorRefAlpha(clus.getSensorID()), x = clus.getX(); if (!trfit.rotate(alpha) || // note: here we also calculate the L,T integral (in the inward direction, but this is irrelevant) @@ -1763,7 +1830,7 @@ bool MatchTPCITS::refitABTrack(int iITSAB, const TPCABSeed& seed, pmr::vectorgetSensorRefAlpha(clus.getSensorID()), x = clus.getX(); if (!tracOut.rotate(alpha, refLin, propagator->getNominalBz()) || // note: here we also calculate the L,T integral @@ -1921,10 +1988,15 @@ int MatchTPCITS::prepareABSeeds() //______________________________________________ int MatchTPCITS::prepareInteractionTimes() { - // guess interaction times from various sources and relate with ITS rofs + // Guess interaction times from various sources and relate them to the ITS cluster ROFs of the clock layer. + // For every candidate define also the per-layer cluster ROFs compatible with its time within the + // mABROFMarginMUS margin: they define the clusters the AfterBurner will check for this candidate. const float ft0Uncertainty = 0.5e-3; - int nITSROFs = mITSROFTimes.size(); - if (mFITInfo.size()) { + const auto& clockROFTimes = mABLayerClusters[mITSClockLayer].rofTimes; + int nClockROFs = clockROFTimes.size(); + int lowestAB = mParams->runAfterBurner ? mParams->lowestLayerAB : NITSLayers; + std::array ptrLr{}; // monotonic pointers on the per-layer ROFs (candidates are time-ordered) + if (mFITInfo.size() && nClockROFs) { int rof = 0; for (const auto& ft : mFITInfo) { if (!mFT0Params->isSelected(ft)) { @@ -1934,22 +2006,33 @@ int MatchTPCITS::prepareInteractionTimes() if (fitTime < 0) { // should not happen continue; } + while (rof < nClockROFs && clockROFTimes[rof].getMax() + mABROFMarginMUS < fitTime) { + rof++; + } + if (rof >= nClockROFs) { // no candidates beyond the last cluster ROF of the clock layer + break; + } if (size_t(fitTime) >= mInteractionMUSLUT.size()) { mInteractionMUSLUT.resize(size_t(fitTime) + 1, -1); } if (mInteractionMUSLUT[fitTime] < 0) { mInteractionMUSLUT[fitTime] = mInteractions.size(); } - for (; rof < nITSROFs; rof++) { - if (mITSROFTimes[rof] < fitTime) { - continue; + auto& intCand = mInteractions.emplace_back(ft.getInteractionRecord(), fitTime, ft0Uncertainty, rof, o2::detectors::DetID::FT0); + for (int lr = lowestAB; lr < NITSLayers; lr++) { // relate the candidate time to compatible ROFs of the AB layers + const auto& rofTimes = mABLayerClusters[lr].rofTimes; + int nROFsLr = (int)rofTimes.size(); + int& ptr = ptrLr[lr]; + while (ptr < nROFsLr && rofTimes[ptr].getMax() + mABROFMarginMUS < fitTime) { + ptr++; + } + if (ptr < nROFsLr && fitTime >= rofTimes[ptr].getMin() - mABROFMarginMUS) { + intCand.rofLr[lr] = ptr; + if (ptr + 1 < nROFsLr && fitTime >= rofTimes[ptr + 1].getMin() - mABROFMarginMUS) { + intCand.rofNextLr |= 0x1 << lr; // compatible also with the next ROF of the layer + } } - break; - } - if (rof >= nITSROFs) { - break; } - mInteractions.emplace_back(ft.getInteractionRecord(), fitTime, ft0Uncertainty, rof, o2::detectors::DetID::FT0); } } int ent = 0; @@ -1963,6 +2046,164 @@ int MatchTPCITS::prepareInteractionTimes() return mInteractions.size(); } +//______________________________________________ +void MatchTPCITS::prepareABClusters() +{ + // Build per (layer, ROF) blocks of the clusters usable by the AfterBurner (i.e. not attached to ITS + // tracks), sorted in (chip, Z). Only the blocks referenced by interaction candidates with seeds are built. + const bool perLayer = mRecoCont->getITSPerLayer(); + const int lowestAB = mParams->lowestLayerAB; + for (int lr = lowestAB; lr < NITSLayers; lr++) { + auto& abClus = mABLayerClusters[lr]; + abClus.needROF.assign(abClus.rofTimes.size(), 0); + abClus.rofRefs.assign(abClus.rofTimes.size(), ClusRange(-1, 0)); + } + for (const auto& intCand : mInteractions) { // mark the blocks to build + if (!intCand.seedsRef.getEntries()) { + continue; + } + for (int lr = lowestAB; lr < NITSLayers; lr++) { + int rof = intCand.rofLr[lr]; + if (rof < 0) { + continue; + } + auto& need = mABLayerClusters[lr].needROF; + need[rof] = 1; + if ((intCand.rofNextLr >> lr) & 0x1) { + need[rof + 1] = 1; + } + } + } + struct ABBlock { + int lr = 0, rof = 0, nCl = 0, offs = 0; + }; + std::vector blocks; + for (int lr = lowestAB; lr < NITSLayers; lr++) { + const auto& need = mABLayerClusters[lr].needROF; + for (int rof = 0; rof < (int)need.size(); rof++) { + if (need[rof]) { + blocks.push_back({lr, rof}); + } + } + } + if (blocks.empty()) { + return; + } + int nBlocks = blocks.size(); + // count AB-usable clusters of every block +#ifdef WITH_OPENMP +#pragma omp parallel for schedule(dynamic) num_threads(mNThreads) +#endif + for (int ib = 0; ib < nBlocks; ib++) { + auto& blk = blocks[ib]; + int slot = perLayer ? blk.lr : 0; + const auto& rofRec = mITSClusterROFRec[slot][blk.rof]; + const auto& status = mABClusterStatus[slot]; + int last = rofRec.getFirstEntry() + rofRec.getNEntries(), nCl = 0; + if (perLayer) { + for (int icl = rofRec.getFirstEntry(); icl < last; icl++) { + if (status[icl] != MinusTen) { + nCl++; + } + } + } else { // select only the clusters of the block layer + const auto& clusArr = mITSClustersArray[0]; + int chipMin = mABChipsBounds[blk.lr], chipMax = mABChipsBounds[blk.lr + 1]; + for (int icl = rofRec.getFirstEntry(); icl < last; icl++) { + int chip = clusArr[icl].getSensorID(); + if (chip >= chipMin && chip < chipMax && status[icl] != MinusTen) { + nCl++; + } + } + } + blk.nCl = nCl; + } + // assign per-layer storage offsets + std::array layerSizes{}; + for (auto& blk : blocks) { + blk.offs = layerSizes[blk.lr]; + layerSizes[blk.lr] += blk.nCl; + } + for (int lr = lowestAB; lr < NITSLayers; lr++) { + mABLayerClusters[lr].clus.resize(layerSizes[lr]); + } + // fill and sort the blocks +#ifdef WITH_OPENMP +#pragma omp parallel for schedule(dynamic) num_threads(mNThreads) +#endif + for (int ib = 0; ib < nBlocks; ib++) { + const auto& blk = blocks[ib]; + auto& abClus = mABLayerClusters[blk.lr]; + int slot = perLayer ? blk.lr : 0; + const auto& rofRec = mITSClusterROFRec[slot][blk.rof]; + const auto& status = mABClusterStatus[slot]; + const auto& clusArr = mITSClustersArray[slot]; + int last = rofRec.getFirstEntry() + rofRec.getNEntries(); + int chipMin = mABChipsBounds[blk.lr], chipMax = mABChipsBounds[blk.lr + 1]; + ABClusterInfo* dst = abClus.clus.data() + blk.offs; + int nCl = 0; + for (int icl = rofRec.getFirstEntry(); icl < last; icl++) { + const auto& cls = clusArr[icl]; + int chip = cls.getSensorID(); + if ((perLayer || (chip >= chipMin && chip < chipMax)) && status[icl] != MinusTen) { + assert(chip >= chipMin && chip < chipMax); // clusters of a per-layer slot must belong to its layer + dst[nCl++] = {cls.getY(), cls.getZ(), o2::itsmft::composeClusID(slot, icl), chip}; + } + } + assert(nCl == blk.nCl); + std::sort(dst, dst + nCl, [](const ABClusterInfo& a, const ABClusterInfo& b) { return a.chip < b.chip || (a.chip == b.chip && a.z < b.z); }); + abClus.rofRefs[blk.rof].set(blk.offs, nCl); + } +} + +//______________________________________________ +void MatchTPCITS::updateABLayerView(ABLayerView& view, int lr, int rof, int nROFs) const +{ + // Point the thread-local view to the AB clusters block(s) of the layer covering the requested ROF(s), + // refreshing the chip->clusters references. Does nothing if the requested block(s) are already loaded. + if (view.rofKey == rof && view.nROFsKey == nROFs) { + return; + } + for (int chip : view.touchedChips) { // reset previously registered references + view.chipRefs[chip].setEntries(0); + } + view.touchedChips.clear(); + view.rofKey = rof; + view.nROFsKey = nROFs; + view.data = nullptr; + view.nData = 0; + if (rof < 0) { // no ROF of this layer is compatible with the candidate time + return; + } + const auto& abClus = mABLayerClusters[lr]; + const auto& ref0 = abClus.rofRefs[rof]; + assert(ref0.getFirstEntry() >= 0); // the block must have been built in prepareABClusters + if (nROFs == 1) { + view.data = abClus.clus.data() + ref0.getFirstEntry(); + view.nData = ref0.getEntries(); + } else { // merge the 2 blocks, each sorted in (chip, Z) + const auto& ref1 = abClus.rofRefs[rof + 1]; + assert(ref1.getFirstEntry() >= 0); + view.merged.resize(ref0.getEntries() + ref1.getEntries()); + std::merge(abClus.clus.begin() + ref0.getFirstEntry(), abClus.clus.begin() + ref0.getEntriesBound(), + abClus.clus.begin() + ref1.getFirstEntry(), abClus.clus.begin() + ref1.getEntriesBound(), + view.merged.begin(), [](const ABClusterInfo& a, const ABClusterInfo& b) { return a.chip < b.chip || (a.chip == b.chip && a.z < b.z); }); + view.data = view.merged.data(); + view.nData = (int)view.merged.size(); + } + int icl = 0; + while (icl < view.nData) { // register the cluster ranges of the chips + int chip = view.data[icl].chip, jcl = icl + 1; + while (jcl < view.nData && view.data[jcl].chip == chip) { + jcl++; + } + int chipLoc = chip - view.chipOffs; + view.chipRefs[chipLoc].set(icl, jcl - icl); + view.touchedChips.push_back(chipLoc); + icl = jcl; + } +} + //______________________________________________ bool MatchTPCITS::runAfterBurner(pmr::vector& matchedTracks, pmr::vector& matchLabels, pmr::vector& ABTrackletLabels, pmr::vector& ABTrackletClusterIDs, pmr::vector& ABTrackletRefs, pmr::vector>& calib) @@ -1977,33 +2218,77 @@ bool MatchTPCITS::runAfterBurner(pmr::vector& matc return false; } mTimer[SWABMatch].Start(false); - - std::vector itsChipClRefsBuff(mNThreads); -#ifdef ENABLE_UPGRADES - // with upgrades the datatype changed, hence we need to initialize - // each element individually - std::generate(itsChipClRefsBuff.begin(), itsChipClRefsBuff.end(), []() { - return ITSChipClustersRefs(o2::its::GeometryTGeo::Instance()->getNumberOfChips()); - }); -#endif - -#ifdef WITH_OPENMP -#pragma omp parallel for schedule(dynamic) num_threads(mNThreads) -#endif + prepareABClusters(); // build the blocks of AB-usable clusters for the (layer, ROF)s referenced by the candidates + + // Group together consecutive candidates with seeds sharing the same per-layer ROFs: within a group the + // thread-local cluster views stay valid, so they are (re)loaded at most once per group and layer + int lowestAB = mParams->lowestLayerAB; + std::vector seededIC; + std::vector> icGroups; // ranges (1st, last entry) in seededIC + seededIC.reserve(nIntCand); for (int ic = 0; ic < nIntCand; ic++) { const auto& intCand = mInteractions[ic]; - LOGP(debug, "cand T {} Entries: {} : {} : {} | ITS ROF: {}", intCand.tBracket.mean(), intCand.seedsRef.getEntries(), intCand.seedsRef.getFirstEntry(), intCand.seedsRef.getEntriesBound(), intCand.rofITS); if (!intCand.seedsRef.getEntries()) { continue; } + bool anyROF = false; + for (int lr = lowestAB; lr < NITSLayers; lr++) { + if (intCand.rofLr[lr] >= 0) { + anyROF = true; + break; + } + } + if (!anyROF) { // no AB layer has clusters at the candidate time: discard the seeds of this candidate + for (int is = intCand.seedsRef.getFirstEntry(); is < intCand.seedsRef.getEntriesBound(); is++) { + mTPCABSeeds[is].disable(); + } + continue; + } + bool sameROFs = false; + if (!seededIC.empty()) { + const auto& prevCand = mInteractions[seededIC.back()]; + sameROFs = prevCand.rofNextLr == intCand.rofNextLr; + for (int lr = lowestAB; sameROFs && lr < NITSLayers; lr++) { + sameROFs = prevCand.rofLr[lr] == intCand.rofLr[lr]; + } + } + if (sameROFs) { + icGroups.back().second = (int)seededIC.size(); + } else { + icGroups.emplace_back((int)seededIC.size(), (int)seededIC.size()); + } + seededIC.push_back(ic); + } + + std::vector itsClViewsBuff(mNThreads); + for (auto& views : itsClViewsBuff) { // book the chip->clusters references of the AB layers + for (int lr = lowestAB; lr < NITSLayers; lr++) { + auto& view = views.layers[lr]; + view.chipOffs = mABChipsBounds[lr]; + view.chipRefs.resize(mABChipsBounds[lr + 1] - mABChipsBounds[lr], ClusRange(0, 0)); + } + } + +#ifdef WITH_OPENMP +#pragma omp parallel for schedule(dynamic) num_threads(mNThreads) +#endif + for (int ig = 0; ig < (int)icGroups.size(); ig++) { #ifdef WITH_OPENMP uint8_t tid = (uint8_t)omp_get_thread_num(); #else uint8_t tid = 0; #endif - fillClustersForAfterBurner(intCand.rofITS, 1, itsChipClRefsBuff[tid]); // RS FIXME account for possibility of filling 2 ROFs - for (int is = intCand.seedsRef.getFirstEntry(); is < intCand.seedsRef.getEntriesBound(); is++) { // loop over all seeds of this interaction candidate - processABSeed(is, itsChipClRefsBuff[tid], tid); + auto& views = itsClViewsBuff[tid]; + const auto& groupCand = mInteractions[seededIC[icGroups[ig].first]]; // all candidates of the group share the same ROFs + for (int lr = lowestAB; lr < NITSLayers; lr++) { + updateABLayerView(views.layers[lr], lr, groupCand.rofLr[lr], (groupCand.rofNextLr >> lr) & 0x1 ? 2 : 1); + } + for (int j = icGroups[ig].first; j <= icGroups[ig].second; j++) { + const auto& intCand = mInteractions[seededIC[j]]; + LOGP(debug, "cand T {} Entries: {} : {} : {} | ITS clock ROF: {}", intCand.tBracket.mean(), intCand.seedsRef.getEntries(), intCand.seedsRef.getFirstEntry(), intCand.seedsRef.getEntriesBound(), intCand.rofITS); + for (int is = intCand.seedsRef.getFirstEntry(); is < intCand.seedsRef.getEntriesBound(); is++) { // loop over all seeds of this interaction candidate + processABSeed(is, views, tid); + } } } mTimer[SWABMatch].Stop(); @@ -2039,13 +2324,13 @@ bool MatchTPCITS::runAfterBurner(pmr::vector& matc continue; } auto bestID = ABSeed.getBestLinkID(); - if (ABSeed.checkLinkHasUsedClusters(bestID, mABClusterLinkIndex)) { + if (ABSeed.checkLinkHasUsedClusters(bestID, mABClusterStatus)) { ABSeed.setNeedAlternative(); // flag for later processing // RSTMP LOG(info) << "Iter: " << iter << " seed has used clusters " << i << "[" << candAB[i].seedID << "/" << candAB[i].chi2 << "]" << " last lr: " << int(ABSeed.lowestLayer) << " Ncont: " << int(link.nContLayers);; continue; } ABSeed.validate(bestID); - ABSeed.flagLinkUsedClusters(bestID, mABClusterLinkIndex); + ABSeed.flagLinkUsedClusters(bestID, mABClusterStatus); mABWinnersIDs.push_back(tTPC.matchID = candAB[i].seedID); mNABRefsClus += ABSeed.getNLayers(); nwin++; @@ -2079,8 +2364,8 @@ void MatchTPCITS::refitABWinners(pmr::vector& matc } std::map labelOccurence; - auto accountClusterLabel = [&labelOccurence, itsClLabs = mITSClsLabels](int clID) { - auto labels = itsClLabs->getLabels(clID); + auto accountClusterLabel = [&labelOccurence, this](int clID) { + auto labels = mITSClsLabels[o2::itsmft::clusID2Layer(clID)]->getLabels(o2::itsmft::clusID2Index(clID)); for (auto lab : labels) { // check all labels of the cluster if (lab.isSet()) { labelOccurence[lab]++; @@ -2099,7 +2384,7 @@ void MatchTPCITS::refitABWinners(pmr::vector& matc ABTrackletClusterIDs.push_back(winL.clID); ncl++; clref.pattern |= 0x1 << winL.layerID; - clref.setClusterSize(winL.layerID, mITSClusterSizes[winL.clID]); + clref.setClusterSize(winL.layerID, mITSClusterSizes[o2::itsmft::clusID2Layer(winL.clID)][o2::itsmft::clusID2Index(winL.clID)]); if (mMCTruthON) { accountClusterLabel(winL.clID); } @@ -2139,13 +2424,13 @@ void MatchTPCITS::refitABWinners(pmr::vector& matc } //______________________________________________ -void MatchTPCITS::processABSeed(int sid, const ITSChipClustersRefs& itsChipClRefs, uint8_t tID) +void MatchTPCITS::processABSeed(int sid, const ABThreadClusterViews& itsClViews, uint8_t tID) { // prepare matching hypothesis tree for given seed auto& ABSeed = mTPCABSeeds[sid]; ABSeed.threadID = tID; ABSeed.linksEntry = mABLinksPool.threadPool[tID].size(); - followABSeed(ABSeed.track, itsChipClRefs, MinusTen, NITSLayers - 1, ABSeed); // check matches on outermost layer + followABSeed(ABSeed.track, itsClViews.layers[NITSLayers - 1], MinusTen, NITSLayers - 1, ABSeed); // check matches on outermost layer for (int ilr = NITSLayers - 1; ilr > mParams->lowestLayerAB; ilr--) { int nextLinkID = ABSeed.firstInLr[ilr]; if (nextLinkID < 0) { @@ -2158,7 +2443,7 @@ void MatchTPCITS::processABSeed(int sid, const ITSChipClustersRefs& itsChipClRef continue; } int next2nextLinkID = seedLink.nextOnLr; // fetch now since the seedLink may change due to the relocation - followABSeed(seedLink, itsChipClRefs, nextLinkID, ilr - 1, ABSeed); // check matches on the next layer + followABSeed(seedLink, itsClViews.layers[ilr - 1], nextLinkID, ilr - 1, ABSeed); // check matches on the next layer nextLinkID = next2nextLinkID; } } @@ -2185,9 +2470,11 @@ void MatchTPCITS::processABSeed(int sid, const ITSChipClustersRefs& itsChipClRef } //______________________________________________ -int MatchTPCITS::followABSeed(const o2::track::TrackParCov& seed, const ITSChipClustersRefs& itsChipClRefs, int seedID, int lrID, TPCABSeed& ABSeed) +int MatchTPCITS::followABSeed(const o2::track::TrackParCov& seed, const ABLayerView& clView, int seedID, int lrID, TPCABSeed& ABSeed) { - + if (!clView.nData) { // no AB-usable clusters on this layer at the interaction candidate time + return 0; + } auto propagator = o2::base::Propagator::Instance(); float xTgt; const auto& lr = mRGHelper.layers[lrID]; @@ -2248,7 +2535,7 @@ int MatchTPCITS::followABSeed(const o2::track::TrackParCov& seed, const ITSChipC if (lad.chips[chipID].zRange.isOutside(zCross, mParams->nABSigmaZ * errZ)) { continue; } - const auto& clRange = itsChipClRefs.chipRefs[lad.chips[chipID].id]; + const auto& clRange = clView.chipRefs[lad.chips[chipID].id - clView.chipOffs]; if (!clRange.getEntries()) { LOG(debug) << "No clusters in chip range"; continue; @@ -2257,7 +2544,7 @@ int MatchTPCITS::followABSeed(const o2::track::TrackParCov& seed, const ITSChipC float errYcalp = errY * (csa * chipC.csAlp + sna * chipC.snAlp); // sigY_rotate(from alpha0 to alpha1) = sigY * cos(alpha1 - alpha0); float tolerZ = errZ * mParams->nABSigmaZ, tolerY = errYcalp * mParams->nABSigmaY; float yTrack = -xCross * chipC.snAlp + yCross * chipC.csAlp; // track-chip crossing Y in chip frame - if (!preselectChipClusters(chipSelClusters, clRange, itsChipClRefs, yTrack, zCross, tolerY, tolerZ)) { // select candidate clusters for this chip + if (!preselectChipClusters(chipSelClusters, clRange, clView, yTrack, zCross, tolerY, tolerZ)) { // select candidate clusters for this chip LOG(debug) << "No compatible clusters found"; continue; } @@ -2270,7 +2557,7 @@ int MatchTPCITS::followABSeed(const o2::track::TrackParCov& seed, const ITSChipC } for (auto clID : chipSelClusters) { - const auto& cls = mITSClustersArray[clID]; + const auto& cls = getITSCluster(clID); auto chi2 = trcLC.getPredictedChi2(cls); if (chi2 > mParams->cutABTrack2ClChi2) { continue; @@ -2470,79 +2757,44 @@ float MatchTPCITS::correctTPCTrack(o2::track::TrackParCov& trc, const TrackLocTP } //______________________________________________ -void MatchTPCITS::fillClustersForAfterBurner(int rofStart, int nROFs, ITSChipClustersRefs& itsChipClRefs) -{ - // Prepare unused clusters of given ROFs range for matching in the afterburner - // Note: normally only 1 ROF needs to be filled (nROFs==1 ) unless we want - // to account for interaction on the boundary of 2 rofs, which then may contribute to both ROFs. - int first = mITSClusterROFRec[rofStart].getFirstEntry(), last = first; - for (int ir = nROFs; ir--;) { - last += mITSClusterROFRec[rofStart + ir].getNEntries(); - } - itsChipClRefs.clear(); - auto& idxSort = itsChipClRefs.clusterID; - for (int icl = first; icl < last; icl++) { - if (mABClusterLinkIndex[icl] != MinusTen) { // clusters with MinusOne are used in main matching - idxSort.push_back(icl); - } - } - // sort in chip, Z - const auto& clusArr = mITSClustersArray; - std::sort(idxSort.begin(), idxSort.end(), [&clusArr](int i, int j) { - const auto &clI = clusArr[i], &clJ = clusArr[j]; - if (clI.getSensorID() < clJ.getSensorID()) { - return true; - } - if (clI.getSensorID() == clJ.getSensorID()) { - return clI.getZ() < clJ.getZ(); - } - return false; - }); - - int ncl = idxSort.size(); - int lastSens = -1, nClInSens = 0; - ClusRange* chipClRefs = nullptr; - for (int icl = 0; icl < ncl; icl++) { - const auto& clus = mITSClustersArray[idxSort[icl]]; - int sens = clus.getSensorID(); - if (sens != lastSens) { - if (chipClRefs) { // finalize chip reference - chipClRefs->setEntries(nClInSens); - nClInSens = 0; - } - chipClRefs = &itsChipClRefs.chipRefs[(lastSens = sens)]; - chipClRefs->setFirstEntry(icl); - } - nClInSens++; - } - if (chipClRefs) { - chipClRefs->setEntries(nClInSens); // finalize last chip reference - } -} - -//______________________________________________ -void MatchTPCITS::setITSTimeBiasInBC(int n) +void MatchTPCITS::setAlpideParam(const AlpParamITS* p) { - mITSTimeBiasInBC = n; - mITSTimeBiasMUS = mITSTimeBiasInBC * o2::constants::lhc::LHCBunchSpacingNS * 1e-3; -} - -//______________________________________________ -void MatchTPCITS::setITSROFrameLengthMUS(float fums) -{ - mITSROFrameLengthMUS = fums; - mITSTimeResMUS = mITSROFrameLengthMUS / std::sqrt(12.f); - mITSROFrameLengthMUSInv = 1. / mITSROFrameLengthMUS; - mITSROFrameLengthInBC = std::max(1, int(mITSROFrameLengthMUS / (o2::constants::lhc::LHCBunchSpacingNS * 1e-3))); -} - -//______________________________________________ -void MatchTPCITS::setITSROFrameLengthInBC(int nbc) -{ - mITSROFrameLengthInBC = nbc; - mITSROFrameLengthMUS = nbc * o2::constants::lhc::LHCBunchSpacingNS * 1e-3; - mITSTimeResMUS = mITSROFrameLengthMUS / std::sqrt(12.f); - mITSROFrameLengthMUSInv = 1. / mITSROFrameLengthMUS; + ///< derive all ITS ROF timings from the Alpide parameters. The per-layer arrays are always + ///< filled for all NITSLayers, also when every layer shares the same ROF length and bias. + if (!p) { + LOG(fatal) << "ITS Alpide parameters pointer is null"; + } + mAlpParams = p; + constexpr float BCLenMUS = o2::constants::lhc::LHCBunchSpacingMUS; + unsigned int maxNROFsPerOrbit = 0; + mITSClockLayer = 0; + for (int lr = 0; lr < NITSLayers; lr++) { + if (mITSTriggered) { // in the triggered mode all layers are read out together + mITSROFrameLengthMUS[lr] = p->roFrameLengthTrig * 1e-3; + mITSROFrameLengthInBC[lr] = std::max(1, int(mITSROFrameLengthMUS[lr] / BCLenMUS)); + } else { + mITSROFrameLengthInBC[lr] = p->getROFLengthInBC(lr); + mITSROFrameLengthMUS[lr] = mITSROFrameLengthInBC[lr] * BCLenMUS; + } + mITSROFrameLengthMUSInv[lr] = 1.f / mITSROFrameLengthMUS[lr]; + mITSTimeResMUS[lr] = mITSROFrameLengthMUS[lr] / std::sqrt(12.f); + mITSTimeBiasInBC[lr] = p->getROFBiasInBC(lr); + mITSTimeBiasMUS[lr] = mITSTimeBiasInBC[lr] * BCLenMUS; + // the ITS tracks ROFRecords are defined by the layer with the largest number of ROFs per orbit, + // see o2::its::ITSTrackingInterface (o2::its::ROFOverlapView::getClock()). If all layers have the + // same timing, this is the layer 0. + unsigned int nROFsPerOrbit = o2::constants::lhc::LHCMaxBunches / mITSROFrameLengthInBC[lr]; + if (nROFsPerOrbit > maxNROFsPerOrbit) { + maxNROFsPerOrbit = nROFsPerOrbit; + mITSClockLayer = lr; + } + } + std::string rep; + for (int lr = 0; lr < NITSLayers; lr++) { + rep += fmt::format(" L{}:{}/{}", lr, mITSROFrameLengthInBC[lr], mITSTimeBiasInBC[lr]); + } + LOGP(info, "ITS {} readout, per-layer ROFLength/Bias in BC:{} | clock layer {}", + mITSTriggered ? "triggered" : "continuous", rep, mITSClockLayer); } //___________________________________________________________________ @@ -2655,35 +2907,37 @@ void MatchTPCITS::flagUsedITSClusters(const o2::its::TrackITS& track) // flag clusters used by this track int clEntry = track.getFirstClusterEntry(); for (int icl = track.getNumberOfClusters(); icl--;) { - mABClusterLinkIndex[mITSTrackClusIdx[clEntry++]] = MinusTen; + const int clID = mITSTrackClusIdx[clEntry++]; // composed ID: (layer << o2::itsmft::ClusLayerShift) + index_in_layer + auto& clStatus = mABClusterStatus[o2::itsmft::clusID2Layer(clID)]; + if (!clStatus.empty()) { // layers not used by the AfterBurner are not booked + clStatus[o2::itsmft::clusID2Index(clID)] = MinusTen; + } } } //__________________________________________________________ -int MatchTPCITS::preselectChipClusters(std::vector& clVecOut, const ClusRange& clRange, const ITSChipClustersRefs& itsChipClRefs, +int MatchTPCITS::preselectChipClusters(std::vector& clVecOut, const ClusRange& clRange, const ABLayerView& clView, float trackY, float trackZ, float tolerY, float tolerZ) const { clVecOut.clear(); int icID = clRange.getFirstEntry(); for (int icl = clRange.getEntries(); icl--;) { // note: clusters within a chip are sorted in Z - int clID = itsChipClRefs.clusterID[icID++]; // so, we go in clusterID increasing direction - const auto& cls = mITSClustersArray[clID]; - float dz = cls.getZ() - trackZ; - LOG(debug) << "cl" << icl << '/' << clID << " " - << " dZ: " << dz << " [" << tolerZ << "| dY: " << trackY - cls.getY() << " [" << tolerY << "]"; + const auto& cls = clView.data[icID++]; + float dz = cls.z - trackZ; + LOG(debug) << "cl" << icl << '/' << cls.id << " " + << " dZ: " << dz << " [" << tolerZ << "| dY: " << trackY - cls.y << " [" << tolerY << "]"; if (dz > tolerZ) { - float clsZ = cls.getZ(); - LOG(debug) << "Skip the rest since " << trackZ << " < " << clsZ << "\n"; + LOG(debug) << "Skip the rest since " << trackZ << " < " << cls.z << "\n"; break; } else if (dz < -tolerZ) { - LOG(debug) << "Skip cluster dz=" << dz << " Ztr=" << trackZ << " zCl=" << cls.getZ(); + LOG(debug) << "Skip cluster dz=" << dz << " Ztr=" << trackZ << " zCl=" << cls.z; continue; } - if (fabs(trackY - cls.getY()) > tolerY) { - LOG(debug) << "Skip cluster dy= " << trackY - cls.getY() << " Ytr=" << trackY << " yCl=" << cls.getY(); + if (fabs(trackY - cls.y) > tolerY) { + LOG(debug) << "Skip cluster dy= " << trackY - cls.y << " Ytr=" << trackY << " yCl=" << cls.y; continue; } - clVecOut.push_back(clID); + clVecOut.push_back(cls.id); } return clVecOut.size(); } @@ -2735,12 +2989,24 @@ void MatchTPCITS::reportSizes(pmr::vector& matched LOGP(info, "Size SHM, calib : size {:9} cap {:9}", siz, cap); } { - siz = mITSClustersArray.size() * sizeof(ITSCluster); - cap = mITSClustersArray.capacity() * sizeof(ITSCluster); + siz = cap = 0; + for (const auto& clusV : mITSClustersArray) { + siz += clusV.size() * sizeof(ITSCluster); + cap += clusV.capacity() * sizeof(ITSCluster); + } sizTot += siz; capTot += cap; LOGP(info, "Size RSS, mITSClustersArray : size {:9} cap {:9}", siz, cap); // + siz = cap = 0; + for (const auto& abClus : mABLayerClusters) { + siz += abClus.sizeInternal(); + cap += abClus.capInternal(); + } + sizTot += siz; + capTot += cap; + LOGP(info, "Size RSS, mABLayerClusters : size {:9} cap {:9}", siz, cap); + // siz = mMatchRecordsTPC.size() * sizeof(MatchRecord); cap = mMatchRecordsTPC.capacity() * sizeof(MatchRecord); sizTot += siz; @@ -2803,11 +3069,14 @@ void MatchTPCITS::reportSizes(pmr::vector& matched capTot += cap; LOGP(info, "Size RSS, mABWinnersIDs : size {:9} cap {:9}", siz, cap); // - siz = mABClusterLinkIndex.size() * sizeof(int); - cap = mABClusterLinkIndex.capacity() * sizeof(int); + siz = cap = 0; + for (const auto& clStatus : mABClusterStatus) { + siz += clStatus.size() * sizeof(int); + cap += clStatus.capacity() * sizeof(int); + } sizTot += siz; capTot += cap; - LOGP(info, "Size RSS, mABClusterLinkIndex : size {:9} cap {:9}", siz, cap); + LOGP(info, "Size RSS, mABClusterStatus : size {:9} cap {:9}", siz, cap); // for (int is = 0; is < o2::constants::math::NSectors; is++) { siz += mTPCSectIndexCache[is].size() * sizeof(int); @@ -2951,7 +3220,7 @@ void MatchTPCITS::fillTPCITSmatchTree(int itsID, int tpcID, int rejFlag, float c << "rejFlag=" << rejFlag << "multTPC=" << mltTPC << "multITSTr=" << mITSTrackROFRec[trackITS.roFrame].getNEntries() - << "multITSCl=" << mITSClusterROFRec[trackITS.roFrame].getNEntries() + << "multITSCl=" << mITSClusterROFRec[mRecoCont->getITSPerLayer() ? mITSClockLayer : 0][trackITS.roFrame].getNEntries() << "\n"; mTimer[SWDBG].Stop(); @@ -2986,7 +3255,7 @@ void MatchTPCITS::dumpWinnerMatches() (*mDBGOut) << "matchWin" << "multTPC=" << mltTPC << "multITSTr=" << mITSTrackROFRec[tITS.roFrame].getNEntries() - << "multITSCl=" << mITSClusterROFRec[tITS.roFrame].getNEntries() + << "multITSCl=" << mITSClusterROFRec[mRecoCont->getITSPerLayer() ? mITSClockLayer : 0][tITS.roFrame].getNEntries() << "\n"; } mTimer[SWDBG].Stop(); diff --git a/Detectors/GlobalTrackingWorkflow/include/GlobalTrackingWorkflow/TPCITSMatchingSpec.h b/Detectors/GlobalTrackingWorkflow/include/GlobalTrackingWorkflow/TPCITSMatchingSpec.h index 56240fd2c8f98..775e7b8c676f2 100644 --- a/Detectors/GlobalTrackingWorkflow/include/GlobalTrackingWorkflow/TPCITSMatchingSpec.h +++ b/Detectors/GlobalTrackingWorkflow/include/GlobalTrackingWorkflow/TPCITSMatchingSpec.h @@ -23,7 +23,7 @@ namespace o2 namespace globaltracking { /// create a processor spec -framework::DataProcessorSpec getTPCITSMatchingSpec(o2::dataformats::GlobalTrackID::mask_t src, bool useFT0, bool calib, bool skipTPCOnly, bool useGeom, bool useMC, bool requestCTPLumi); +framework::DataProcessorSpec getTPCITSMatchingSpec(o2::dataformats::GlobalTrackID::mask_t src, bool useFT0, bool calib, bool skipTPCOnly, bool useGeom, bool useMC, bool requestCTPLumi, bool itsStag); } // namespace globaltracking } // namespace o2 diff --git a/Detectors/GlobalTrackingWorkflow/src/TPCITSMatchingSpec.cxx b/Detectors/GlobalTrackingWorkflow/src/TPCITSMatchingSpec.cxx index 7f63b61e02be0..90a02fd6fc493 100644 --- a/Detectors/GlobalTrackingWorkflow/src/TPCITSMatchingSpec.cxx +++ b/Detectors/GlobalTrackingWorkflow/src/TPCITSMatchingSpec.cxx @@ -210,15 +210,10 @@ void TPCITSMatchingDPL::updateTimeDependentParams(ProcessingContext& pc) pc.inputs().get("MatchParam"); // Note: ITS/CLUSDICT and ITS/ALPIDEPARAM are requested/loaded by the recocontainer - const auto& alpParams = o2::itsmft::DPLAlpideParam::Instance(); - if (mMatching.isITSTriggered()) { - mMatching.setITSROFrameLengthMUS(alpParams.roFrameLengthTrig / 1.e3); // ITS ROFrame duration in \mus - } else { - mMatching.setITSROFrameLengthInBC(alpParams.roFrameLengthInBC); // ITS ROFrame duration in \mus - } - mMatching.setITSTimeBiasInBC(alpParams.roFrameBiasInBC); - mMatching.setSkipTPCOnly(mSkipTPCOnly); + // all per-layer ITS ROF lengths and biases are derived in MatchTPCITS::setAlpideParam() mMatching.setITSTriggered(!o2::base::GRPGeomHelper::instance().getGRPECS()->isDetContinuousReadOut(o2::detectors::DetID::ITS)); + mMatching.setAlpideParam(&o2::itsmft::DPLAlpideParam::Instance()); + mMatching.setSkipTPCOnly(mSkipTPCOnly); mMatching.setNHBPerTF(o2::base::GRPGeomHelper::instance().getGRPECS()->getNHBFPerTF()); mMatching.setMCTruthOn(mUseMC); mMatching.setUseFT0(mUseFT0); @@ -250,10 +245,11 @@ void TPCITSMatchingDPL::updateTimeDependentParams(ProcessingContext& pc) } } -DataProcessorSpec getTPCITSMatchingSpec(GTrackID::mask_t src, bool useFT0, bool calib, bool skipTPCOnly, bool useGeom, bool useMC, bool requestCTPLumi) +DataProcessorSpec getTPCITSMatchingSpec(GTrackID::mask_t src, bool useFT0, bool calib, bool skipTPCOnly, bool useGeom, bool useMC, bool requestCTPLumi, bool itsStag) { std::vector outputs; auto dataRequest = std::make_shared(); + dataRequest->setITSPerLayer(itsStag); if ((src & GTrackID::getSourcesMask("TPC-TRD,TPC-TOF,TPC-TRD-TOF")).any()) { // preliminary stage of extended workflow ? dataRequest->setMatchingInputStrict(); } diff --git a/Detectors/GlobalTrackingWorkflow/src/tpcits-match-workflow.cxx b/Detectors/GlobalTrackingWorkflow/src/tpcits-match-workflow.cxx index 79ca13430ccd9..05c3f15d6df51 100644 --- a/Detectors/GlobalTrackingWorkflow/src/tpcits-match-workflow.cxx +++ b/Detectors/GlobalTrackingWorkflow/src/tpcits-match-workflow.cxx @@ -96,7 +96,7 @@ WorkflowSpec defineDataProcessing(o2::framework::ConfigContext const& configcont if (!configcontext.options().get("disable-root-input")) { specs.emplace_back(o2::tpc::getTPCScalerSpec(sclOpt)); } - specs.emplace_back(o2::globaltracking::getTPCITSMatchingSpec(srcL, useFT0, calib, !GID::includesSource(GID::TPC, src), useGeom, useMC, sclOpt.requestCTPLumi)); + specs.emplace_back(o2::globaltracking::getTPCITSMatchingSpec(srcL, useFT0, calib, !GID::includesSource(GID::TPC, src), useGeom, useMC, sclOpt.requestCTPLumi, doStag)); if (!configcontext.options().get("disable-root-output")) { specs.emplace_back(o2::globaltracking::getTrackWriterTPCITSSpec(useMC)); diff --git a/Detectors/ITSMFT/ITS/tracking/src/TrackingInterface.cxx b/Detectors/ITSMFT/ITS/tracking/src/TrackingInterface.cxx index a0e8d708cffa2..732a68208732a 100644 --- a/Detectors/ITSMFT/ITS/tracking/src/TrackingInterface.cxx +++ b/Detectors/ITSMFT/ITS/tracking/src/TrackingInterface.cxx @@ -26,6 +26,7 @@ #include "ITSMFTTracking/ITSTrackingConfigParam.h" #include "ITStracking/TrackingInterface.h" +#include "DataFormatsITSMFT/ClusterID.h" #include "DataFormatsITSMFT/ROFRecord.h" #include "DataFormatsITSMFT/PhysTrigger.h" #include "DataFormatsTRD/TriggerRecord.h" @@ -333,7 +334,9 @@ void ITSTrackingInterface::run(framework::ProcessingContext& pc) auto clid = trc.getClusterIndex(ic); if (clid >= 0) { trc.setClusterSize(ic, mTimeFrame->getClusterSize((mDoStaggering) ? ic : 0, clid)); - allClusIdx.push_back(clid); + // with the per-layer clusters input the index is local to the layer, hence the layer must be + // encoded into the stored reference; with the monolithic input the composed ID is just the index + allClusIdx.push_back(o2::itsmft::composeClusID((mDoStaggering) ? ic : 0, clid)); nclf++; } }