diff --git a/Detectors/TPC/calibration/CMakeLists.txt b/Detectors/TPC/calibration/CMakeLists.txt index b0b2704d7ea00..aaec52b71e691 100644 --- a/Detectors/TPC/calibration/CMakeLists.txt +++ b/Detectors/TPC/calibration/CMakeLists.txt @@ -165,6 +165,10 @@ o2_add_test_root_macro(macro/drawCMV.C COMPILE_ONLY PUBLIC_LINK_LIBRARIES O2::TPCCalibration O2::TPCBase LABELS tpc) +o2_add_test_root_macro(macro/calculatedEdx.C + COMPILE_ONLY + PUBLIC_LINK_LIBRARIES O2::TPCCalibration O2::TPCBase + LABELS tpc) o2_add_test(IDCFourierTransform COMPONENT_NAME calibration diff --git a/Detectors/TPC/calibration/include/TPCCalibration/CalculatedEdx.h b/Detectors/TPC/calibration/include/TPCCalibration/CalculatedEdx.h index 4d8c4e89322a8..fc896b198956b 100644 --- a/Detectors/TPC/calibration/include/TPCCalibration/CalculatedEdx.h +++ b/Detectors/TPC/calibration/include/TPCCalibration/CalculatedEdx.h @@ -20,16 +20,40 @@ // o2 includes #include "DataFormatsTPC/TrackTPC.h" #include "DataFormatsTPC/dEdxInfo.h" +#include "DataFormatsTPC/VDriftCorrFact.h" +#include "TPCBase/Mapper.h" #include "GPUO2InterfaceRefit.h" #include "CalibdEdxContainer.h" +#include "CorrectionMapsHelper.h" #include "CommonUtils/TreeStreamRedirector.h" #include "TPCCalibration/CorrectdEdxDistortions.h" #include "TPCFastTransformPOD.h" +#include "GPUCommonRtypes.h" +#include "SimulationDataFormat/MCCompLabel.h" #include +#include +#include +#include +#include +#include + +namespace o2::gpu +{ +class TPCFastTransform; +} namespace o2::tpc { +/// \brief average cluster occupancy of a track, per TPC region +struct AverageOccupancy { + double IROC = 0.; + double OROC1 = 0.; + double OROC2 = 0.; + double OROC3 = 0.; + ClassDefNV(AverageOccupancy, 1); +}; + /// \brief dEdx calculation class /// /// This class is used to calculate dEdx of reconstructed tracks. @@ -46,7 +70,7 @@ namespace o2::tpc /// c.setMembers(tpcTrackClIdxVecInput, clusterIndex, tpcTracks); // set the member variables: TrackTPC, TPCClRefElem, o2::tpc::ClusterNativeAccess /// c.setRefit(); // set the refit pointer to perform refitting of tracks, otherwise setPropagateTrack to true /// start looping over the tracks -/// c.calculatedEdx(track, output, 0.015, 0.60, CorrectionFlags::TopologyPol | CorrectionFlags::dEdxResidual, ClusterFlags::ExcludeEdgeCl) // this will fill the dEdxInfo output for given track +/// c.calculatedEdx(track, output, averageOcc, 0.015, 0.60, CorrectionFlags::TopologyPol | CorrectionFlags::dEdxResidual, ClusterFlags::ExcludeEdgeCl) // this will fill the dEdxInfo output and per-region average track occupancy averageOcc for given track enum class CorrectionFlags : unsigned short { None = 0, @@ -61,11 +85,14 @@ enum class CorrectionFlags : unsigned short { enum class ClusterFlags : unsigned short { None = 0, ExcludeSingleCl = 1 << 0, ///< flag to exclude single clusters in dEdx calculation - ExcludeSplitCl = 1 << 1, ///< flag to exclude split clusters in dEdx calculation - ExcludeEdgeCl = 1 << 2, ///< flag to exclude sector edge clusters in dEdx calculation - ExcludeSubthresholdCl = 1 << 3, ///< flag to exclude subthreshold clusters in dEdx calculation - ExcludeSectorBoundaries = 1 << 4, ///< flag to exclude sector boundary clusters in subthreshold cluster treatment - ExcludeSharedCl = 1 << 5, ///< flag to exclude clusters shared between tracks + ExcludeSplitPadCl = 1 << 1, ///< flag to exclude split pad clusters in dEdx calculation + ExcludeSplitTimeCl = 1 << 2, ///< flag to exclude split time clusters in dEdx calculation + ExcludeSplitCl = 1 << 3, ///< flag to exclude split pad or time clusters in dEdx calculation + ExcludeEdgeCl = 1 << 4, ///< flag to exclude sector edge clusters in dEdx calculation + ExcludeSubthresholdCl = 1 << 5, ///< flag to exclude subthreshold clusters in dEdx calculation + ExcludeSectorBoundaries = 1 << 6, ///< flag to exclude sector boundary clusters in subthreshold cluster treatment + ExcludeSharedCl = 1 << 7, ///< flag to exclude clusters shared between tracks in dEdx calculation + ExcludeSamePadRowCl = 1 << 8, ///< flag to exclude clusters in the same pad row in dEdx calculation }; inline CorrectionFlags operator&(CorrectionFlags a, CorrectionFlags b) { return static_cast(static_cast(a) & static_cast(b)); } @@ -76,6 +103,32 @@ inline ClusterFlags operator&(ClusterFlags a, ClusterFlags b) { return static_ca inline ClusterFlags operator~(ClusterFlags a) { return static_cast(~static_cast(a)); } inline ClusterFlags operator|(ClusterFlags a, ClusterFlags b) { return static_cast(static_cast(a) | static_cast(b)); } +/// \brief bundles the settings of one calculatedEdx() call (everything except the track/output/averageOcc) +/// used by calculatedEdxMultipleSettings() to evaluate several settings for the same track without repeating +/// the track refit/propagation for every setting +struct dEdxSettings { + ClusterFlags clusterMask = ClusterFlags::None; ///< clusters to exclude + CorrectionFlags correctionMask = CorrectionFlags::TopologyPol | CorrectionFlags::dEdxResidual; ///< corrections to apply + unsigned short subthresholdMethod = 0; ///< subthreshold cluster charge filling method + unsigned short stackBoundaryMethod = 0; ///< stack boundary cluster exclusion method + unsigned short sameRowClusterMethod = 2; ///< same (sector,row) cluster group combination method: 0 = do not merge (one sample per cluster), 1 = merge and sum the qTot, 2 = merge and use the largest qTot + float maxSubthresholdChargeTot = 100000.f; ///< upper limit for the per-region minimum qTot used as the virtual charge of a subthreshold cluster (default effectively disables the cap) + float maxSubthresholdChargeMax = 100000.f; ///< upper limit for the per-region minimum qMax used as the virtual charge of a subthreshold cluster (default effectively disables the cap) + float low = 0.015f; ///< lower cluster cut + float high = 0.6f; ///< higher cluster cut + std::string debugRootFile = "dEdxDebug.root"; ///< debug streamer output file used if mDebug is set +}; + +/// \brief per-cluster info needed by the calculatedEdx()/calculatedEdxMultipleSettings() overloads that take clusters directly instead of extracting them from the track via mTPCTrackClIdxVecInput/mClusterIndex +/// isShared cannot be looked up from mTPCRefitterShMap for externally supplied clusters, so it must be supplied by the caller +struct ClInfo { + unsigned char sectorIndex = 0; + unsigned char rowIndex = 0; + bool isShared = false; +}; + +using ClInfoVec = std::vector; + class CalculatedEdx { public: @@ -92,9 +145,21 @@ class CalculatedEdx /// set the refitter void setRefit(const unsigned int nHbfPerTf = 32); + /// Should be called before setRefit(), call setTPCVDrift() after this; if called after setRefit() instead, the + /// existing refitter (which holds a pointer into the correction map buffer being replaced here) is dropped and + /// must be recreated with a new setRefit() call before further use + void setTPCCorrMap(const o2::gpu::TPCFastTransform& corrMap); + + /// Should be called before setRefit(), and after setTPCCorrMap() if that is used too; if called after setRefit() + /// instead, the existing refitter is dropped and must be recreated with a new setRefit() call before further use + void setTPCVDrift(const o2::tpc::VDriftCorrFact& v); + /// \param propagate propagate the tracks to extract the track parameters instead of performing a refit void setPropagateTrack(const bool propagate) { mPropagateTrack = propagate; } + /// \param propagate propagate the tracks to extract the track parameters instead of performing a refit + void setPropagateParams(const bool propagate) { mPropagateParams = propagate; } + /// \param debug use debug streamer and set debug vectors void setDebug(const bool debug) { mDebug = debug; } @@ -104,42 +169,114 @@ class CalculatedEdx /// \param maxMissingCl maximum number of missing clusters for subthreshold check void setMaxMissingCl(int maxMissingCl) { mMaxMissingCl = maxMissingCl; } - /// \param minChargeTotThreshold upper limit for the possible minimum charge tot in subthreshold treatment - void setMinChargeTotThreshold(float minChargeTotThreshold) { mMinChargeTotThreshold = minChargeTotThreshold; } + /// \param n subthreshold clusters are not filled within min(nRows/2, n) rows of the track's outer end, mirroring the online tracker's allowChangeClusters gate, + // set to 0 to fill every 1-row gap regardless of its position on the track + void setSubThreshEdgeRows(int n) { mSubThreshEdgeRows = n; } + + /// \param d max |pad_i - pad_j| within a same (sector,row) cluster group for it to be eligible for merging; a larger spread means the group is looper legs, always kept separate + void setSameRowMaxPadDiff(float d) { mSameRowMaxPadDiff = d; } - /// \param minChargeMaxThreshold upper limit for the possible minimum charge max in subthreshold treatment - void setMinChargeMaxThreshold(float minChargeMaxThreshold) { mMinChargeMaxThreshold = minChargeMaxThreshold; } + /// \param d max |time_i - time_j| (time bins) within a same (sector,row) cluster group for it to be eligible for merging; a larger spread means looper legs, always kept separate + void setSameRowMaxTimeDiff(float d) { mSameRowMaxTimeDiff = d; } - /// set the debug streamer - void setStreamer(const char* debugRootFile) { mStreamer = std::make_unique(debugRootFile, "recreate"); }; + /// set the debug streamer for a given output file; a new streamer is only created the first time a given debugRootFile is seen, + /// so different calculatedEdx() calls using different debugRootFile names each get their own independent debug file + void setStreamer(const char* debugRootFile) + { + auto& streamer = mStreamers[debugRootFile]; + if (!streamer) { + streamer = std::make_unique(debugRootFile, "recreate"); + } + }; /// set the debug streamer of the space-charge dedx correction void setSCStreamer(const char* debugRootFile = "debug_sc_corrections.root") { mSCdEdxCorrection.setStreamer(debugRootFile); } + /// \param lumi set luminosity for space-charge correction map scaling + void setLumi(const float lumi) { mSCdEdxCorrection.setLumi(lumi); } + /// \return returns magnetic field in kG float getFieldNominalGPUBz() { return mFieldNominalGPUBz; } /// \return returns maxMissingCl for subthreshold cluster treatment int getMaxMissingCl() { return mMaxMissingCl; } - /// \return returns the upper limit for the possible minimum charge tot in subthreshold treatment - float getMinChargeTotThreshold() { return mMinChargeTotThreshold; } + /// \return returns the outer-end row exclusion for subthreshold cluster treatment + int getSubThreshEdgeRows() const { return mSubThreshEdgeRows; } + + /// \return returns the max pad spread for a same (sector,row) cluster group to be eligible for merging + float getSameRowMaxPadDiff() const { return mSameRowMaxPadDiff; } + + /// \return returns the max time spread for a same (sector,row) cluster group to be eligible for merging + float getSameRowMaxTimeDiff() const { return mSameRowMaxTimeDiff; } + + /// \return returns the number of rows where refit/propagation failed (row.propagationFailed) since the last resetDebugCounters(); with setRefit(), this only counts rows where the propagation fallback was also unable to recover the row + long getNPropagationFailed() const { return mNPropagationFailed; } + + /// \return returns the number of rows where setRefit()'s RefitTrackAsGPU() could not reach the row and the row's track state instead came from the propagation fallback since the last resetDebugCounters(); always 0 outside setRefit() mode + long getNRefitFallback() const { return mNRefitFallback; } - /// \return returns the upper limit for the possible minimum charge max in subthreshold treatment - float getMinChargeMaxThreshold() { return mMinChargeMaxThreshold; } + /// \return returns the number of rows gathered by gatherRowClusterData() (processed for refit/propagation) since the last resetDebugCounters() + long getNRowsProcessed() const { return mNRowsProcessed; } - /// fill missing clusters with minimum charge (method=0) or minimum charge/2 (method=1) or Landau (method=2) - void fillMissingClusters(int missingClusters[4], float minChargeTot, float minChargeMax, int method, std::array, 5>& chargeTotROC, std::array, 5>& chargeMaxROC); + /// \return returns the number of row gaps filled as subthreshold clusters by calculatedEdxFromRowData() since the last resetDebugCounters() per setting + const std::vector& getNSubThresholdFilledPerSettings() const { return mNSubThresholdFilledPerSettings; } + + /// reset the running counters returned by getNPropagationFailed()/getNRefitFallback()/getNRowsProcessed()/getNSubThresholdFilledPerSettings() + void resetDebugCounters() + { + mNPropagationFailed = 0; + mNRefitFallback = 0; + mNRowsProcessed = 0; + mNSubThresholdFilledPerSettings.clear(); + } + + /// fill missing clusters per region with that region's running minimum charge (method=0) or half of it (method=1), + /// \param missingClusters number of row gaps to fill, per region (IROC, OROC1, OROC2, OROC3) + /// \param minChargeTot per-region running minimum qTot among the accepted clusters of that region + /// \param minChargeMax per-region running minimum qMax among the accepted clusters of that region + void fillMissingClusters(int missingClusters[4], const float minChargeTot[4], const float minChargeMax[4], int method, std::array, 5>& chargeTotROC, std::array, 5>& chargeMaxROC); + + /// \param rowOrder (sector, row) keys in the order they are first encountered while scanning the track's native cluster references (0..nClusterReferences-1), i.e. the track's true physical row-traversal order + /// \param mergeableRows (sector, row) keys of groups with >1 cluster that pass the setSameRowMaxPadDiff()/setSameRowMaxTimeDiff() proximity gate, i.e. are eligible for dEdxSettings::sameRowClusterMethod to merge them; a group missing from here (be it size 1, or size >1 but far apart -- looper legs) is always emitted as one sample per raw cluster, regardless of sameRowClusterMethod + void handleSameRowClusters(o2::tpc::TrackTPC& track, std::vector>& rowOrder, std::map, std::vector>& clustersByRow, std::set>& mergeableRows, std::map>& clusterReferencesByIndex); + + /// same as handleSameRowClusters() above, but groups externally supplied clusters instead of the track's clusters accessed via mTPCTrackClIdxVecInput/mClusterIndex + /// \param rowOrder (sector, row) keys in the order they are first encountered while scanning clusters (0..clusters.size()-1), i.e. the order they were supplied in + /// \param mergeableRows (sector, row) keys of groups with >1 cluster that pass the setSameRowMaxPadDiff()/setSameRowMaxTimeDiff() proximity gate, i.e. are eligible for dEdxSettings::sameRowClusterMethod to merge them; a group missing from here (be it size 1, or size >1 but far apart -- looper legs) is always emitted as one sample per raw cluster, regardless of sameRowClusterMethod + void handleSameRowClusters(const std::vector& clusters, const ClInfoVec& clusterInfos, std::vector>& rowOrder, std::map, std::vector>& clustersByRow, std::set>& mergeableRows); /// get the truncated mean for the input track with the truncation range, charge type, region and corrections /// the cluster charge is normalized by effective length*gain, you can turn off the normalization by setting all corrections to false /// \param track input track /// \param output output dEdxInfo + /// \param averageOcc output average cluster occupancy of the track, per TPC region /// \param low lower cluster cut /// \param high higher cluster cut - /// \param mask to apply different corrections: TopologySimple = simple analytical topology correction, TopologyPol = topology correction from polynomials, GainFull = full gain map from calibration container, + /// \param correctionMask to apply different corrections: TopologySimple = simple analytical topology correction, TopologyPol = topology correction from polynomials, GainFull = full gain map from calibration container, /// GainResidual = residuals gain map from calibration container, dEdxResidual = residual dEdx correction - void calculatedEdx(TrackTPC& track, dEdxInfo& output, float low = 0.015f, float high = 0.6f, CorrectionFlags correctionMask = CorrectionFlags::TopologyPol | CorrectionFlags::dEdxResidual, ClusterFlags clusterMask = ClusterFlags::None, int subthresholdMethod = 0, const char* debugRootFile = "dEdxDebug.root"); + /// \param maxSubthresholdChargeTot upper limit for the per-region minimum qTot used as the virtual charge of a subthreshold cluster + /// \param maxSubthresholdChargeMax upper limit for the per-region minimum qMax used as the virtual charge of a subthreshold cluster + /// \param sameRowClusterMethod same (sector,row) cluster group combination method: 0 = do not merge (one sample per cluster), 1 = merge and sum the qTot, 2 = merge and use the largest qTot + void calculatedEdx(TrackTPC& track, dEdxInfo& output, AverageOccupancy& averageOcc, float low = 0.015f, float high = 0.6f, CorrectionFlags correctionMask = CorrectionFlags::TopologyPol | CorrectionFlags::dEdxResidual, ClusterFlags clusterMask = ClusterFlags::None, int subthresholdMethod = 0, int stackBoundaryMethod = 0, const char* debugRootFile = "dEdxDebug.root", float maxSubthresholdChargeTot = 100000.f, float maxSubthresholdChargeMax = 100000.f, int sameRowClusterMethod = 2); + + /// evaluate several dEdx settings for the same track while performing the track refit/propagation to each cluster row only once + /// \param track input track + /// \param outputs output dEdxInfo, filled with one entry per entry in settingsList, in the same order + /// \param averageOcc output average cluster occupancy of the track, per TPC region; a single value, since occupancy does not depend on the dEdx settings and is therefore the same for every entry in settingsList + /// \param settingsList list of dEdx settings to evaluate for this track + /// \param mcLabel if non-null and mDebug is set, written to the "dEdxDebugTrack" row of every settingsList entry so debug rows can be matched back to the true MC track + void calculatedEdxMultipleSettings(TrackTPC& track, std::vector& outputs, AverageOccupancy& averageOcc, const std::vector& settingsList, const MCCompLabel* mcLabel = nullptr); + + /// same as calculatedEdx() above, but takes the track's clusters and per-cluster info directly instead of extracting them from the track via mTPCTrackClIdxVecInput/mClusterIndex + /// \param clusters clusters of the track, one entry per entry in clusterInfos + /// \param clusterInfos per-cluster (sectorIndex, rowIndex, isShared), one entry per entry in clusters + void calculatedEdx(TrackTPC& track, const std::vector& clusters, const ClInfoVec& clusterInfos, dEdxInfo& output, AverageOccupancy& averageOcc, float low = 0.015f, float high = 0.6f, CorrectionFlags correctionMask = CorrectionFlags::TopologyPol | CorrectionFlags::dEdxResidual, ClusterFlags clusterMask = ClusterFlags::None, int subthresholdMethod = 0, int stackBoundaryMethod = 0, const char* debugRootFile = "dEdxDebug.root", float maxSubthresholdChargeTot = 100000.f, float maxSubthresholdChargeMax = 100000.f, int sameRowClusterMethod = 2); + + /// same as calculatedEdxMultipleSettings() above, but takes the track's clusters and per-cluster info directly instead of extracting them from the track via mTPCTrackClIdxVecInput/mClusterIndex + /// \param clusters clusters of the track, one entry per entry in clusterInfos + /// \param clusterInfos per-cluster (sectorIndex, rowIndex, isShared), one entry per entry in clusters + void calculatedEdxMultipleSettings(TrackTPC& track, const std::vector& clusters, const ClInfoVec& clusterInfos, std::vector& outputs, AverageOccupancy& averageOcc, const std::vector& settingsList, const MCCompLabel* mcLabel = nullptr); /// get the truncated mean for the input charge vector and the truncation range low*nCl1 raw cluster (split hit or looper legs) has one entry per raw cluster here, so + /// dEdxSettings::sameRowClusterMethod == 0 (never merge) can still expose each fragment as its own sample + /// with its own threshold/gain/occupancy, exactly as if it were gathered on its own + struct RowFragment { + o2::tpc::ClusterNative cl; + unsigned char pad; ///< clamped local pad index used for the threshold/gain/dead-channel lookups below + float threshold; + float gain; + float gainResidual; + unsigned int occupancy; + bool isShared; + bool isDeadRegion; + }; + + /// dEdxSettings::sameRowClusterMethod == 2: the fragment with the largest qTot in a mergeable group + static const RowFragment& pickDominantFragment(const std::vector& fragments); + + /// \brief per (sector,row) cluster/track data gathered once per track by gatherRowClusterData(), independent of the dEdx settings reused by calculatedEdxFromRowData() for each entry in a settingsList so the track refit/propagation done in gatherRowClusterData() is not repeated per setting + struct RowClusterData { + std::vector fragments; ///< usually size 1; >1 for a same (sector,row) group. calculatedEdxFromRowData() reduces this to 1..N effective samples per dEdxSettings::sameRowClusterMethod + bool mergeable; ///< true if fragments.size()>1 and they pass the setSameRowMaxPadDiff()/setSameRowMaxTimeDiff() proximity gate; meaningless (never read) when fragments.size()<=1 + RowFragment mergedFragment; ///< dEdxSettings::sameRowClusterMethod==1's synthesized sum-merged sample (charge-weighted pad/time, summed qTot, max qMax) and its own threshold/gain/gainResidual/occupancy/isDeadRegion lookups; computed once here (like the per-fragment quantities above) instead of once per settings entry. Only valid when mergeable && fragments.size()>1 + o2::tpc::TrackTPC trackSnapshot; ///< track state after refit/propagation to this row; identical for every fragment since it only depends on the row's X, not on which cluster + unsigned char sectorIndex; + unsigned char rowIndex; + unsigned int region; + GEMstack stack; + int stackNumber; + StackID stackID; + bool propagationFailed; ///< true if refit/propagation to this row failed, or the resulting track param is NaN, and no fallback recovered it either + bool refitFellBack; ///< true if setRefit() mode's RefitTrackAsGPU() could not reach this row + int missingClusters; ///< number of skipped rows since the previous entry in rowData (i.e. rowIndex - previous rowIndex - 1); same for every settings entry since rowOrder does not depend on the settings + bool sameSectorAsPrevRow; ///< true if this row's sector equals the previous entry in rowData's sector + bool missingClusterGapDeadOrEdge; ///< true if any of the missingClusters skipped row(s) would land on a dead channel or off the padrow edge + std::vector inputClusterIndices; ///< indices of the input clusters merged into this row: positions in the caller-supplied clusters vector (externally-supplied-cluster overload) or the track's cluster-reference list (reference overload). Streamed to "dEdxDebugCl" so external tooling can map a row back to its input cluster(s) + }; + + /// gather, for every (sector, row) of the track's row-traversal order, performing the refit/propagation to each cluster row exactly once + /// \param track input track, mutated in place by refit/propagation + /// \param rowData output per-row data + /// \param averageOcc output average cluster occupancy of the track, per TPC region + void gatherRowClusterData(o2::tpc::TrackTPC& track, std::vector& rowData, AverageOccupancy& averageOcc); + + /// same as gatherRowClusterData() above, but sources clusters from externally supplied clusters/clusterInfos instead of the track's own cluster references + /// \param track input track, mutated in place by refit/propagation + /// \param clusters clusters of the track, one entry per entry in clusterInfos + /// \param clusterInfos per-cluster (sectorIndex, rowIndex, isShared), one entry per entry in clusters + /// \param rowData output per-row data + /// \param averageOcc output average cluster occupancy of the track, per TPC region + void gatherRowClusterData(o2::tpc::TrackTPC& track, const std::vector& clusters, const ClInfoVec& clusterInfos, std::vector& rowData, AverageOccupancy& averageOcc); + + /// per-row processing shared by both gatherRowClusterData() overloads, once the row's raw cluster fragments/sector/row are known regardless of where they came from: + /// refits/propagates the track to this row exactly once (shared by every fragment, since it only depends on the row's X), looks up each fragment's threshold/gain/occupancy/dead-channel status, checks the missing-cluster gap, and appends the resulting group entry to rowData + /// \param track input track, mutated in place by refit/propagation + /// \param fragmentClusters raw native cluster(s) grouped into this (sector,row); usually size 1 + /// \param fragmentIsShared per-fragment isShared, one entry per fragmentClusters + /// \param sectorIndex sector of this row + /// \param rowIndex TPC row index + /// \param mergeable true if fragmentClusters.size()>1 and they pass the pad/time proximity gate (see setSameRowMaxPadDiff()/setSameRowMaxTimeDiff()); ignored when fragmentClusters.size()<=1 + /// \param rowIndexOld rowIndex of the previous entry appended to rowData (255 if this is the first row) + /// \param sectorIndexOld sectorIndex of the previous entry appended to rowData (255 if this is the first row) + /// \param occupancyROC per-region occupancy accumulator, updated in place + /// \param rowData output per-row data; the new row is appended, and rowData.back() (if non-empty) is read as the previous row for the missing-cluster-gap check + /// \param refitAbandoned setRefit() mode only: false as long as RefitTrackAsGPU() keeps succeeding; the first + /// time it fails for this track, set to true and stays true for the rest of the track. On that first + /// failure, the propagation fallback resumes from the track's state just before the failed attempt + /// (i.e. its state after the last successfully refit row), not from the track's pristine pre-loop state. + void gatherRowClusterDataForRow(o2::tpc::TrackTPC& track, const std::vector& fragmentClusters, const std::vector& fragmentIsShared, unsigned char sectorIndex, unsigned char rowIndex, bool mergeable, unsigned char rowIndexOld, unsigned char sectorIndexOld, std::array, 4>& occupancyROC, std::vector& rowData, bool& refitAbandoned); + + /// geometrically propagate track (rotating into the row's sector frame first) to xPosition + /// \return true if any of the three attempts succeeded (track left at xPosition); false if all failed (track left unchanged, at its state on entry) + bool propagateTrackToX(o2::track::TrackParCov& track, float xPosition, unsigned char sectorIndex) const; + + /// re-flatten mTPCCorrMapFull into mTPCCorrMap/mTPCCorrMapBuffer + void rebuildTPCCorrMapPOD(); + + /// compute the dEdx output for one dEdx settings entry from the row data previously gathered by gatherRowClusterData() + /// \param rowData per row data gathered by gatherRowClusterData() for the track being processed + /// \param settings dEdx settings to apply + /// \param settingsIndex index of settings within its settingsList + /// \param trackTime0 track.getTime0() of the track being processed, captured before refit/propagation (unaffected by it) + /// \param trackOrig pristine track (before refit/propagation mutated it), used for the debug "dEdxDebugTrack" row; ignored if mDebug is false + /// \param averageOcc average cluster occupancy of the track as computed by gatherRowClusterData(), only used for the debug "dEdxDebugTrack" row; ignored if mDebug is false + /// \param output output dEdxInfo + /// \param mcLabel if non-null and mDebug is set, written to the "dEdxDebugTrack" row as the "mcLabel" branch + void calculatedEdxFromRowData(const std::vector& rowData, const dEdxSettings& settings, size_t settingsIndex, float trackTime0, const o2::tpc::TrackTPC& trackOrig, const AverageOccupancy& averageOcc, dEdxInfo& output, const MCCompLabel* mcLabel = nullptr); + std::vector* mTracks{nullptr}; ///< vector containing the tpc tracks which will be processed std::vector* mTPCTrackClIdxVecInput{nullptr}; ///< input vector with TPC tracks cluster indicies const o2::tpc::ClusterNativeAccess* mClusterIndex{nullptr}; ///< needed to access clusternative with tpctracks - const o2::gpu::TPCFastTransformPOD* mTPCCorrMap{nullptr}; ///< cluster correction maps helper + const o2::gpu::TPCFastTransformPOD* mTPCCorrMap{nullptr}; ///< cluster correction maps helper; flattened from mTPCCorrMapFull by rebuildTPCCorrMapPOD(); this is the only pointer setRefit()'s GPUO2InterfaceRefit ever sees o2::gpu::aligned_unique_buffer_ptr mTPCCorrMapBuffer; + std::unique_ptr mTPCCorrMapFull; ///< regular (non-POD), exclusively-owned transform mutated by setTPCCorrMap()/setTPCVDrift(); rebuildTPCCorrMapPOD() re-flattens it into mTPCCorrMap/mTPCCorrMapBuffer after every change std::vector mTPCRefitterShMap; ///< externally set TPC clusters sharing map std::vector mTPCRefitterOccMap; ///< externally set TPC clusters occupancy map std::unique_ptr mRefit{nullptr}; ///< TPC refitter used for TPC tracks refit during the reconstruction - int mMaxMissingCl{1}; ///< maximum number of missing clusters for subthreshold check - float mMinChargeTotThreshold{50}; ///< upper limit for minimum charge tot value in subthreshold treatment, i.e for a high dEdx track adding a minimum value of 500 to track as a virtual charge doesn't make sense - float mMinChargeMaxThreshold{50}; ///< upper limit for minimum charge max value in subthreshold treatment, i.e for a high dEdx track adding a minimum value of 500 to track as a virtual charge doesn't make sense - float mFieldNominalGPUBz{5}; ///< magnetic field in kG, used for track propagation - bool mPropagateTrack{false}; ///< propagating the track instead of performing a refit - bool mDebug{false}; ///< use the debug streamer - CalibdEdxContainer mCalibCont; ///< calibration container - std::unique_ptr mStreamer{nullptr}; ///< debug streamer + int mMaxMissingCl{1}; ///< maximum number of missing clusters for subthreshold check + int mSubThreshEdgeRows{30}; ///< no subthreshold fill within min(nRows/2, this) rows of the track's outer end 0 disables + float mSameRowMaxPadDiff{3.f}; ///< max pad spread within a same (sector,row) group for it to be eligible for dEdxSettings::sameRowClusterMethod to merge it as one split hit (else = looper legs, always kept separate) + float mSameRowMaxTimeDiff{4.f}; ///< max time-bin spread within a same (sector,row) group for it to be eligible for merging + float mFieldNominalGPUBz{5}; ///< magnetic field in kG, used for track propagation + bool mPropagateTrack{false}; ///< propagating the track instead of performing a refit (faster than refit) + bool mPropagateParams{false}; ///< propagating the parameters instead of full propagation (faster than track propagation) + bool mDebug{false}; ///< use the debug streamer + CalibdEdxContainer mCalibCont; ///< calibration container + std::unordered_map> mStreamers; ///< debug streamers, keyed by output file name so each debugRootFile gets its own tree + long mDebugTrackIndex{-1}; ///< running index of the track being processed, written to the debug trees so per-cluster rows can be grouped back into tracks + long mNPropagationFailed{0}; ///< number of rows where refit/propagation failed (and the fallback below, if applicable, also failed) since the last resetDebugCounters() + long mNRefitFallback{0}; ///< number of rows where setRefit()'s RefitTrackAsGPU() failed but the propagation fallback recovered the row, since the last resetDebugCounters() + long mNRowsProcessed{0}; ///< number of rows gathered by gatherRowClusterData() since the last resetDebugCounters() + std::vector mNSubThresholdFilledPerSettings; ///< number of row gaps filled as subthreshold clusters, per dEdxSettings list index, since the last resetDebugCounters() CorrectdEdxDistortions mSCdEdxCorrection; ///< for space-charge correction of dE/dx + + std::array, 4> mStackBoundaries = {{{0, 62}, {63, 96}, {97, 126}, {127, 151}}}; // for excluding stack boundaries in dEdx calculation }; } // namespace o2::tpc -#endif +#endif \ No newline at end of file diff --git a/Detectors/TPC/calibration/macro/calculatedEdx.C b/Detectors/TPC/calibration/macro/calculatedEdx.C new file mode 100644 index 0000000000000..315a3b5d61c22 --- /dev/null +++ b/Detectors/TPC/calibration/macro/calculatedEdx.C @@ -0,0 +1,373 @@ +// 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 calculatedEdx.C +/// \brief Example macro showing how to use o2::tpc::CalculatedEdx to calculate TPC dE/dx from TPC tracks and native clusters. +/// Supports real data (CTF- or TF-reconstructed) and MC productions, and optionally restricts the calculation to TPC tracks matched to an ITS track. +/// Accepts a list of dEdx settings (truncation range, correction mask, cluster mask, subthreshold/stack-boundary method). +/// The track refit/propagation is performed only once per track no matter how many settings entries are given. +/// The output tree has one row per track, with one "dEdx" branch per entry in settingsList (see the "dEdx: low=..., high=..." log lines for what each index means). + +#if !defined(__CLING__) || defined(__ROOTCLING__) +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include "TFile.h" +#include "TROOT.h" +#include "TTree.h" +#include "Framework/Logger.h" +#include "CommonUtils/TreeStreamRedirector.h" +#include "CommonDataFormat/TFIDInfo.h" +#include "ReconstructionDataFormats/TrackTPCITS.h" +#include "DataFormatsITS/TrackITS.h" +#include "DataFormatsTPC/TrackTPC.h" +#include "DataFormatsTPC/ClusterNative.h" +#include "DataFormatsTPC/ClusterNativeHelper.h" +#include "SimulationDataFormat/MCCompLabel.h" +#include "TPCCalibration/CalculatedEdx.h" +#endif + +using namespace o2::tpc; +namespace fs = std::filesystem; + +namespace +{ +TFile* openOrNull(const std::string& fileName) +{ + auto f = TFile::Open(fileName.data()); + if (!f || !f->IsOpen() || f->IsZombie()) { + LOGP(error, "Could not open file {}", fileName); + return nullptr; + } + return f; +} + +/// open fileName and retrieve treeName from it; logs why and returns {nullptr, nullptr} if either step fails +std::pair, std::unique_ptr> openTreeOrNull(const std::string& fileName, const char* treeName) +{ + std::unique_ptr file(openOrNull(fileName)); + if (!file) { + return {nullptr, nullptr}; + } + std::unique_ptr tree((TTree*)file->Get(treeName)); + if (!tree) { + LOGP(error, "Could not find tree '{}' in {}", treeName, fileName); + return {nullptr, nullptr}; + } + return {std::move(file), std::move(tree)}; +} + +/// default single-entry settings list +std::vector defaultSettingsList() +{ + dEdxSettings s; + s.low = 0.015f; + s.high = 0.6f; + s.correctionMask = CorrectionFlags::TopologyPol | CorrectionFlags::dEdxResidual; + s.clusterMask = ClusterFlags::ExcludeEdgeCl; + s.subthresholdMethod = 0; + s.stackBoundaryMethod = 0; + s.sameRowClusterMethod = 2; + return {s}; +} +} // namespace + +/// \param dir directory containing tpctracks.root, tpc-native-clusters.root and, if isMatchedToITS, o2trac_its.root/o2match_itstpc.root +/// \param runNumberOrTimeStamp run number or timestamp used to load the calibration objects from CCDB for every timeframe, overridden per timeframe whenever a tfIDFileName file is found in dir +/// \param outFile name of the output file with the calculated dE/dx tree +/// \param localCCDBFolder if non-empty, load calibration objects from local CCDB snapshots in this folder instead of from the CCDB server +/// \param settingsList list of dEdx settings to evaluate for every track (truncation range, correction/cluster mask, +/// subthreshold/stack-boundary/same-row-cluster method) +/// \param useRefit refit the tracks at each cluster row using the GPU refitter (default) +/// \param propagateTrack propagate the full track, including material corrections, instead of refitting; only used if useRefit is false +/// \param propagateParams propagate only the track parameters (fastest option, no material corrections); only used if useRefit and propagateTrack are both false +/// \param debug enable the CalculatedEdx debug streamer, additionally writing one dEdxDebug.root file with per-cluster information (or dEdxDebug_s.root per settings entry j, if settingsList has more than one entry) +/// \param isMC set to true for MC productions to read the true MC track label (TPCTracksMCTruth branch) for each track into an additional "mcLabel" branch +/// \param isMatchedToITS set to true to also read o2trac_its.root/o2match_itstpc.root and restrict the dE/dx calculation to TPC tracks matched to an ITS track +/// \param tfIDFileName name of the optional timeframe ID file; if found in dir, its per-entry time stamp is used for the CCDB instead of runNumberOrTimeStamp +/// \param tpcTracksFileName name of the file with the TPC tracks +/// \param clusterNativeFileName name of the file with the TPC native clusters +/// \param itsTracksFileName name of the file with the ITS tracks; only used if isMatchedToITS is true +/// \param matchFileName name of the file with the TPC-ITS match information; only used if isMatchedToITS is true +/// \param maxEvents if >= 0, process at most this many events, instead of all available ones, default -1 processes all events +/// \param applySCToRefitTransform put the space-charge-corrected cluster map into the refit transform; only used if useRefit is true +/// \param applyVDriftToRefitTransform apply the calibrated drift velocity/t0 to the refit transform; only used if useRefit is true +void calculatedEdx(const std::string dir = ".", + const long runNumberOrTimeStamp = 0, + const std::string outFile = "dEdxCalc.root", + const std::string localCCDBFolder = "", + const std::vector settingsList = defaultSettingsList(), + const bool useRefit = true, + const bool propagateTrack = false, + const bool propagateParams = false, + const bool debug = false, + const bool isMC = false, + const bool isMatchedToITS = false, + const std::string tfIDFileName = "o2_tfidinfo.root", + const std::string tpcTracksFileName = "tpctracks.root", + const std::string clusterNativeFileName = "tpc-native-clusters.root", + const std::string itsTracksFileName = "o2trac_its.root", + const std::string matchFileName = "o2match_itstpc.root", + const long long maxEvents = -1, + const bool applySCToRefitTransform = false, + const bool applyVDriftToRefitTransform = true) +{ + if (settingsList.empty()) { + LOGP(error, "settingsList must not be empty"); + return; + } + for (size_t i = 0; i < settingsList.size(); i++) { + const auto& s = settingsList[i]; + if (s.low < 0.f || s.high > 1.f || s.low >= s.high) { + LOGP(error, "settingsList[{}]: invalid truncation range [{}, {}); expected 0 <= low < high <= 1", i, s.low, s.high); + return; + } + if (s.subthresholdMethod != 0 && s.subthresholdMethod != 1) { + LOGP(error, "settingsList[{}]: invalid subthresholdMethod {}; expected 0 (minimum charge) or 1 (minimum charge / 2)", i, s.subthresholdMethod); + return; + } + if (s.stackBoundaryMethod > 2) { + LOGP(error, "settingsList[{}]: invalid stackBoundaryMethod {}; expected 0 (disabled), 1 (exclude boundary row) or 2 (also exclude the adjacent row)", i, s.stackBoundaryMethod); + return; + } + if (s.sameRowClusterMethod > 2) { + LOGP(error, "settingsList[{}]: invalid sameRowClusterMethod {}; expected 0 (do not merge), 1 (merge, sum qTot) or 2 (merge, highest qTot)", i, s.sameRowClusterMethod); + return; + } + // the output tree stores one "dEdx" branch per settingsList entry; log what each index means here + LOGP(info, "dEdx{}: low={}, high={}, correctionMask={}, clusterMask={}, subthresholdMethod={}, stackBoundaryMethod={}, sameRowClusterMethod={}", i, s.low, s.high, static_cast(s.correctionMask), static_cast(s.clusterMask), s.subthresholdMethod, s.stackBoundaryMethod, s.sameRowClusterMethod); + } + + const std::clock_t c_start = std::clock(); + const auto t_start = std::chrono::high_resolution_clock::now(); + + CalculatedEdx calcdEdx; + calcdEdx.setDebug(debug); + calcdEdx.setPropagateTrack(propagateTrack); + calcdEdx.setPropagateParams(propagateParams); + + // own copy of settingsList, with debugRootFile made unique per settings entry (if there is more than one) so debug streams from different settings never collide + std::vector activeSettingsList = settingsList; + if (debug) { + for (size_t iSettings = 0; iSettings < settingsList.size(); iSettings++) { + activeSettingsList[iSettings].debugRootFile = (settingsList.size() == 1) ? "dEdxDebug.root" : fmt::format("dEdxDebug_s{}.root", iSettings); + } + } + + auto [tpcFile, tpcTree] = openTreeOrNull(fmt::format("{}/{}", dir, tpcTracksFileName), "tpcrec"); + if (!tpcTree) { + return; + } + + std::vector tpcTracks, *tpcTracksPtr = &tpcTracks; + std::vector* tpcTrackClIdxVecInput{nullptr}; + tpcTree->SetBranchAddress("TPCTracks", &tpcTracksPtr); + tpcTree->SetBranchAddress("ClusRefs", &tpcTrackClIdxVecInput); + + std::vector tpcMCTruth, *tpcMCTruthPtr = &tpcMCTruth; + if (isMC) { + if (!tpcTree->GetBranch("TPCTracksMCTruth")) { + LOGP(error, "Branch 'TPCTracksMCTruth' not found in {}/{}, cannot resolve MC truth", dir, tpcTracksFileName); + return; + } + tpcTree->SetBranchAddress("TPCTracksMCTruth", &tpcMCTruthPtr); + } + + std::unique_ptr itsFile; + std::unique_ptr itsTree; + std::unique_ptr matchFile; + std::unique_ptr matchTree; + std::vector itsTracks, *itsTracksPtr = &itsTracks; + std::vector matchTracks, *matchTracksPtr = &matchTracks; + + if (isMatchedToITS) { + std::tie(itsFile, itsTree) = openTreeOrNull(fmt::format("{}/{}", dir, itsTracksFileName), "o2sim"); + std::tie(matchFile, matchTree) = openTreeOrNull(fmt::format("{}/{}", dir, matchFileName), "matchTPCITS"); + if (!itsTree || !matchTree) { + return; + } + itsTree->SetBranchAddress("ITSTrack", &itsTracksPtr); + matchTree->SetBranchAddress("TPCITS", &matchTracksPtr); + } + + std::unique_ptr tfIDFile; + std::unique_ptr tfIDTree; + o2::dataformats::TFIDInfo* tfIDInfo{nullptr}; + Long64_t timeStamp = runNumberOrTimeStamp; + const auto tfIDFullName = fmt::format("{}/{}", dir, tfIDFileName); + if (fs::exists(tfIDFullName)) { + std::tie(tfIDFile, tfIDTree) = openTreeOrNull(tfIDFullName, "tfidTree"); + if (tfIDTree) { + tfIDTree->SetBranchAddress("tfidinfo", &tfIDInfo); + tfIDTree->SetBranchAddress("ts", &timeStamp); + LOGP(info, "Using per-time-frame CCDB time stamps from {}", tfIDFullName); + } + } + + const auto clName = fmt::format("{}/{}", dir, clusterNativeFileName); + if (!fs::exists(clName)) { + LOGP(error, "Cluster file {} does not exist", clName); + return; + } + ClusterNativeHelper::Reader tpcClusterReader{}; + tpcClusterReader.init(clName.data()); + if (tpcClusterReader.getTreeSize() == 0) { + LOGP(error, "Could not read a native cluster tree from {}", clName); + return; + } + + o2::utils::TreeStreamRedirector stream(outFile.data(), "recreate"); + + ClusterNativeAccess clusterIndex{}; + std::unique_ptr clusterBuffer{}; + ClusterNativeHelper::ConstMCLabelContainerViewWithBuffer clusterMCBuffer; + memset(&clusterIndex, 0, sizeof(clusterIndex)); + + long long nEvents = tpcTree->GetEntriesFast(); + if (isMatchedToITS) { + nEvents = std::min(nEvents, std::min(itsTree->GetEntriesFast(), matchTree->GetEntriesFast())); + } + if (tfIDTree && tfIDTree->GetEntriesFast() < nEvents) { + LOGP(error, "tfIDInfo tree has fewer entries ({}) than the data trees ({}); ignoring it and using runNumberOrTimeStamp for all events", + tfIDTree->GetEntriesFast(), nEvents); + tfIDTree.reset(); + tfIDFile.reset(); + } + if (maxEvents >= 0 && maxEvents < nEvents) { + LOGP(info, "Limiting to the first {} of {} available events (maxEvents)", maxEvents, nEvents); + nEvents = maxEvents; + } + + for (long long iEvent = 0; iEvent < nEvents; iEvent++) { + tpcTree->GetEntry(iEvent); + if (isMC && tpcMCTruth.size() != tpcTracks.size()) { + LOGP(error, "TPCTracksMCTruth size ({}) does not match TPCTracks size ({}) for event {}, skipping event", + tpcMCTruth.size(), tpcTracks.size(), iEvent); + continue; + } + tpcClusterReader.read(iEvent); + tpcClusterReader.fillIndex(clusterIndex, clusterBuffer, clusterMCBuffer); + if (isMatchedToITS) { + itsTree->GetEntry(iEvent); + matchTree->GetEntry(iEvent); + } + if (tfIDTree) { + tfIDTree->GetEntry(iEvent); + } + + // setMembers()/loadCalibs/setRefit()... depend on tracks, clusters and timestamp, so they must be redone per event + calcdEdx.setMembers(tpcTrackClIdxVecInput, clusterIndex, &tpcTracks); + if (localCCDBFolder.empty()) { + bool loadSCCorrMap = false; + for (const auto& s : settingsList) { + loadSCCorrMap |= (s.correctionMask & CorrectionFlags::dEdxSC) == CorrectionFlags::dEdxSC; + } + calcdEdx.loadCalibsFromCCDB(timeStamp, isMC, loadSCCorrMap, applySCToRefitTransform && useRefit, applyVDriftToRefitTransform && useRefit); + } else { + calcdEdx.loadCalibsFromLocalCCDBFolder(localCCDBFolder.data()); + } + if (useRefit) { + calcdEdx.setRefit(); + } + + const size_t nSelectable = isMatchedToITS ? matchTracks.size() : tpcTracks.size(); + LOGP(info, "Processing event {} with {} {} and {} settings", iEvent, nSelectable, isMatchedToITS ? "matched tracks" : "tracks", settingsList.size()); + + std::vector tpcOut; + std::vector itsTracksOut; + std::vector matchTracksOut; + std::vector> dEdxOut; // [track][settings] + std::vector averageOccOut; + std::vector mcLabelOut; + + for (size_t i = 0; i < nSelectable; i++) { + size_t tpcIndex = i; + if (isMatchedToITS) { + const auto& itstpc = matchTracks[i]; + if (itstpc.getRefITS().getSource() != o2::dataformats::GlobalTrackID::ITS) { + continue; + } + tpcIndex = itstpc.getRefTPC().getIndex(); + itsTracksOut.emplace_back(itsTracks[itstpc.getRefITS().getIndex()]); + matchTracksOut.emplace_back(itstpc); + } + + TrackTPC track(tpcTracks[tpcIndex]); // local copy: refit/propagation inside calculatedEdxMultipleSettings mutate the track in place + std::vector dEdxVec; + AverageOccupancy averageOcc; + calcdEdx.calculatedEdxMultipleSettings(track, dEdxVec, averageOcc, activeSettingsList, isMC ? &tpcMCTruth[tpcIndex] : nullptr); + + tpcOut.emplace_back(track); + dEdxOut.emplace_back(std::move(dEdxVec)); + averageOccOut.emplace_back(averageOcc); + if (isMC) { + mcLabelOut.emplace_back(tpcMCTruth[tpcIndex]); + } + } + + // per-event summary: refit/propagation failures and how many row gaps were filled as subthreshold clusters per settingsList entry + const long nPropagationFailed = calcdEdx.getNPropagationFailed(); + const long nRowsProcessed = calcdEdx.getNRowsProcessed(); + const auto& nSubThresholdFilledPerSettings = calcdEdx.getNSubThresholdFilledPerSettings(); + calcdEdx.resetDebugCounters(); + std::string subThresholdBreakdown; + for (size_t i = 0; i < nSubThresholdFilledPerSettings.size(); i++) { + subThresholdBreakdown += fmt::format("{}dEdx{}={}", i > 0 ? ", " : "", i, nSubThresholdFilledPerSettings[i]); + } + LOGP(info, "Event {}: refit/propagation failed for {}/{} rows ({:.2f}%); gap-cluster(s) filled as subthreshold per settings entry: {}", + iEvent, nPropagationFailed, nRowsProcessed, nRowsProcessed > 0 ? 100. * nPropagationFailed / nRowsProcessed : 0., subThresholdBreakdown); + + // one row per track, with one "dEdx" branch per entry in settingsList + for (size_t i = 0; i < dEdxOut.size(); i++) { + auto& row = stream << "tree" + << "iEvent=" << iEvent + << "timeStamp=" << timeStamp + << "tpc=" << tpcOut[i] + << "averageOcc=" << averageOccOut[i]; + for (size_t iSettings = 0; iSettings < settingsList.size(); iSettings++) { + row << fmt::format("dEdx{}=", iSettings).c_str() << dEdxOut[i][iSettings]; + } + if (tfIDTree) { + row << "tfIDInfo=" << tfIDInfo; + } + if (isMatchedToITS) { + row << "its=" << itsTracksOut[i] + << "itstpc=" << matchTracksOut[i]; + } + if (isMC) { + const auto& label = mcLabelOut[i]; + row << "mcLabel=" << label; + } + row << "\n"; + } + } + + stream.Close(); + + const std::clock_t c_end = std::clock(); + const auto t_end = std::chrono::high_resolution_clock::now(); + + std::cout << std::fixed << std::setprecision(2) + << "CPU time used: " + << (1000.0 * (c_end - c_start) / CLOCKS_PER_SEC) / 60000.0 << " minutes\n" + << "Wall clock time passed: " + << std::chrono::duration(t_end - t_start).count() / 60000.0 << " minutes\n"; +} diff --git a/Detectors/TPC/calibration/src/CalculatedEdx.cxx b/Detectors/TPC/calibration/src/CalculatedEdx.cxx index 18b2f6e3010c7..c949d149343f4 100644 --- a/Detectors/TPC/calibration/src/CalculatedEdx.cxx +++ b/Detectors/TPC/calibration/src/CalculatedEdx.cxx @@ -28,16 +28,102 @@ #include "GPUO2InterfaceUtils.h" #include "GPUTPCGMMergedTrackHit.h" +#include +#include + using namespace o2::tpc; +namespace +{ +// sentinel for minChargeTotROC/minChargeMaxROC +constexpr float kNoValidCharge = std::numeric_limits::max(); + +// dEdxSettings::sameRowClusterMethod == 1: one synthesized cluster summing a mergeable group's fragments' +// qTot/qMax and charge-weighted pad/time +o2::tpc::ClusterNative buildMergedClusterSum(const std::vector& fragments) +{ + float weightedPadSum = 0.f; + float weightedTimeSum = 0.f; + float totalCharge = 0.f; + uint16_t maxCharge = 0; + + const o2::tpc::ClusterNative& firstCluster = fragments[0]; + o2::tpc::ClusterNative combinedCluster = firstCluster; + + for (const auto& cl : fragments) { + const float clPad = cl.getPad(); + const float clTime = cl.getTime(); + const uint16_t clqTot = cl.getQtot(); + const uint16_t clqMax = cl.qMax; + + weightedPadSum += clPad * clqTot; + weightedTimeSum += clTime * clqTot; + totalCharge += clqTot; + maxCharge = std::max(maxCharge, clqMax); + } + + if (totalCharge > o2::tpc::ClusterNative::maxRegularQtot) { + combinedCluster.setSaturatedQtot(static_cast(totalCharge)); + } else { + combinedCluster.qTotPacked = static_cast(totalCharge); + } + combinedCluster.qMax = maxCharge; + combinedCluster.padPacked = static_cast(weightedPadSum / totalCharge * o2::tpc::ClusterNative::scalePadPacked); + combinedCluster.timeFlagsPacked = (static_cast(weightedTimeSum / totalCharge * o2::tpc::ClusterNative::scaleTimePacked) & 0xFFFFFF) | (firstCluster.timeFlagsPacked & 0xFF000000); + return combinedCluster; +} +} // namespace + +const CalculatedEdx::RowFragment& CalculatedEdx::pickDominantFragment(const std::vector& fragments) +{ + const RowFragment* dominant = &fragments[0]; + for (const auto& frag : fragments) { + if (frag.cl.getQtot() > dominant->cl.getQtot()) { + dominant = &frag; + } + } + return *dominant; +} + CalculatedEdx::CalculatedEdx() { + mTPCCorrMapFull = TPCFastTransformHelperO2::instance()->create(0); + rebuildTPCCorrMapPOD(); +} + +void CalculatedEdx::rebuildTPCCorrMapPOD() +{ + // re-flatten the regular mTPCCorrMapFull into the POD buffer setRefit()'s GPU refitter reads gpu::aligned_unique_buffer_ptr buffer; - gpu::TPCFastTransformPOD::create(buffer, *TPCFastTransformHelperO2::instance()->create(0)); + gpu::TPCFastTransformPOD::create(buffer, *mTPCCorrMapFull); mTPCCorrMapBuffer = std::move(buffer); mTPCCorrMap = mTPCCorrMapBuffer.get(); } +void CalculatedEdx::setTPCCorrMap(const o2::gpu::TPCFastTransform& corrMap) +{ + if (mRefit) { + // the existing refitter holds a pointer into the correction map buffer we are about to replace; drop it + // rather than leave it dangling. Normal per-event usage (loadCalibsFromCCDB()/loadCalibsFromLocalCCDBFolder() + // followed by setRefit()) re-creates it right after this call with the new map, so this is not an error. + LOGP(warning, "CalculatedEdx::setTPCCorrMap() called after setRefit(); invalidating the existing refitter, call setRefit() again before using it."); + mRefit.reset(); + } + mTPCCorrMapFull = TPCFastTransformHelperO2::instance()->create(0, corrMap.getCorrection()); + rebuildTPCCorrMapPOD(); +} + +void CalculatedEdx::setTPCVDrift(const o2::tpc::VDriftCorrFact& v) +{ + if (mRefit) { + // see setTPCCorrMap() above: drop the now-stale refitter instead of leaving it dangling + LOGP(warning, "CalculatedEdx::setTPCVDrift() called after setRefit(); invalidating the existing refitter, call setRefit() again before using it."); + mRefit.reset(); + } + TPCFastTransformHelperO2::instance()->updateCalibration(*mTPCCorrMapFull, 0, v.corrFact, v.refVDrift, v.getTimeOffset()); + rebuildTPCCorrMapPOD(); +} + void CalculatedEdx::setMembers(std::vector* tpcTrackClIdxVecInput, const o2::tpc::ClusterNativeAccess& clIndex, std::vector* vTPCTracksArrayInp) { mTracks = vTPCTracksArrayInp; @@ -47,7 +133,7 @@ void CalculatedEdx::setMembers(std::vector* tpcTrackClIdx void CalculatedEdx::setRefit(const unsigned int nHbfPerTf) { - mTPCRefitterShMap.reserve(mClusterIndex->nClustersTotal); + mTPCRefitterShMap.resize(mClusterIndex->nClustersTotal); auto sizeOcc = o2::gpu::GPUO2InterfaceRefit::fillOccupancyMapGetSize(nHbfPerTf, nullptr); mTPCRefitterOccMap.resize(sizeOcc); std::fill(mTPCRefitterOccMap.begin(), mTPCRefitterOccMap.end(), 0); @@ -55,7 +141,7 @@ void CalculatedEdx::setRefit(const unsigned int nHbfPerTf) mRefit = std::make_unique(mClusterIndex, mTPCCorrMap, mFieldNominalGPUBz, mTPCTrackClIdxVecInput->data(), nHbfPerTf, mTPCRefitterShMap.data(), mTPCRefitterOccMap.data(), mTPCRefitterOccMap.size()); } -void CalculatedEdx::fillMissingClusters(int missingClusters[4], float minChargeTot, float minChargeMax, int method, std::array, 5>& chargeTotROC, std::array, 5>& chargeMaxROC) +void CalculatedEdx::fillMissingClusters(int missingClusters[4], const float minChargeTot[4], const float minChargeMax[4], int method, std::array, 5>& chargeTotROC, std::array, 5>& chargeMaxROC) { if (method != 0 && method != 1) { LOGP(info, "Unrecognized subthreshold cluster treatment. Not adding virtual charges to the track!"); @@ -63,9 +149,15 @@ void CalculatedEdx::fillMissingClusters(int missingClusters[4], float minChargeT } for (int roc = 0; roc < 4; roc++) { + // minChargeTot/MaxROC[roc] is only ever updated from an accepted real cluster in that ROC + // a region this track never accepted a single cluster in leaves it at kNoValidCharge, so skip + // this ROC's fill entirely rather than injecting the sentinel as if it were a measured charge + if (minChargeTot[roc] >= kNoValidCharge || minChargeMax[roc] >= kNoValidCharge) { + continue; + } + const float chargeTot = (method == 1) ? minChargeTot[roc] / 2.f : minChargeTot[roc]; + const float chargeMax = (method == 1) ? minChargeMax[roc] / 2.f : minChargeMax[roc]; for (int i = 0; i < missingClusters[roc]; i++) { - float chargeTot = (method == 1) ? minChargeTot / 2.f : minChargeTot; - float chargeMax = (method == 1) ? minChargeMax / 2.f : minChargeMax; chargeTotROC[roc].emplace_back(chargeTot); chargeTotROC[4].emplace_back(chargeTot); @@ -76,329 +168,710 @@ void CalculatedEdx::fillMissingClusters(int missingClusters[4], float minChargeT } } -void CalculatedEdx::calculatedEdx(o2::tpc::TrackTPC& track, dEdxInfo& output, float low, float high, CorrectionFlags correctionMask, ClusterFlags clusterMask, int subthresholdMethod, const char* debugRootFile) +void CalculatedEdx::handleSameRowClusters(o2::tpc::TrackTPC& track, std::vector>& rowOrder, std::map, std::vector>& clustersByRow, std::set>& mergeableRows, std::map>& clusterReferencesByIndex) { // get number of clusters const int nClusters = track.getNClusterReferences(); - int nClsROC[4] = {0, 0, 0, 0}; - int nClsSubThreshROC[4] = {0, 0, 0, 0}; - - const int nType = 5; - std::array, nType> chargeTotROC; - std::array, nType> chargeMaxROC; - for (int i = 0; i < nType; ++i) { - chargeTotROC[i].reserve(Mapper::PADROWS); - chargeMaxROC[i].reserve(Mapper::PADROWS); - } - - // debug vectors - std::vector excludeClVector; - std::vector regionVector; - std::vector rowIndexVector; - std::vector padVector; - std::vector sectorVector; - std::vector stackVector; - std::vector localXVector; - std::vector localYVector; - std::vector offsPadVector; - - std::vector topologyCorrVector; - std::vector topologyCorrTotVector; - std::vector topologyCorrMaxVector; - std::vector gainVector; - std::vector gainResidualVector; - std::vector residualCorrTotVector; - std::vector residualCorrMaxVector; - std::vector scCorrVector; - - std::vector trackVector; - std::vector clVector; - std::vector occupancyVector; - std::vector isClusterShared; - - if (mDebug) { - excludeClVector.reserve(nClusters); - regionVector.reserve(nClusters); - rowIndexVector.reserve(nClusters); - padVector.reserve(nClusters); - stackVector.reserve(nClusters); - sectorVector.reserve(nClusters); - localXVector.reserve(nClusters); - localYVector.reserve(nClusters); - offsPadVector.reserve(nClusters); - topologyCorrVector.reserve(nClusters); - topologyCorrTotVector.reserve(nClusters); - topologyCorrMaxVector.reserve(nClusters); - gainVector.reserve(nClusters); - gainResidualVector.reserve(nClusters); - residualCorrTotVector.reserve(nClusters); - residualCorrMaxVector.reserve(nClusters); - trackVector.reserve(nClusters); - clVector.reserve(nClusters); - scCorrVector.reserve(nClusters); - occupancyVector.reserve(nClusters); - isClusterShared.reserve(nClusters); - } - - // for missing clusters - unsigned char rowIndexOld = 0; - unsigned char sectorIndexOld = 0; - float minChargeTot = 100000.f; - float minChargeMax = 100000.f; - - // loop over the clusters + // group clusters by (sector, row) for (int iCl = 0; iCl < nClusters; iCl++) { - const o2::tpc::ClusterNative& cl = track.getCluster(*mTPCTrackClIdxVecInput, iCl, *mClusterIndex); unsigned char sectorIndex = 0; unsigned char rowIndex = 0; unsigned int clusterIndexNumb = 0; - // set sectorIndex, rowIndex, clusterIndexNumb track.getClusterReference(*mTPCTrackClIdxVecInput, iCl, sectorIndex, rowIndex, clusterIndexNumb); - // check if the cluster is shared - const unsigned int absoluteIndex = mClusterIndex->clusterOffset[sectorIndex][rowIndex] + clusterIndexNumb; - const bool isShared = mRefit ? (mTPCRefitterShMap[absoluteIndex] & o2::gpu::GPUTPCGMMergedTrackHit::flagShared) : 0; + const auto rowKey = std::make_pair(sectorIndex, rowIndex); + if (clustersByRow.find(rowKey) == clustersByRow.end()) { + rowOrder.emplace_back(rowKey); + } - // get region, pad, stack and stack ID - const int region = Mapper::REGION[rowIndex]; - const unsigned char pad = std::clamp(static_cast(cl.getPad() + 0.5f), static_cast(0), Mapper::PADSPERROW[region][Mapper::getLocalRowFromGlobalRow(rowIndex)] - 1); // the left side of the pad is defined at e.g. 3.5 and the right side at 4.5 - const CRU cru(Sector(sectorIndex), region); - const auto stack = cru.gemStack(); - StackID stackID{sectorIndex, stack}; - // the stack number for debugging - const int stackNumber = static_cast(stack); + // add the cluster index to the corresponding (sector, row) key in clustersByRow + clustersByRow[rowKey].emplace_back(iCl); - // get local coordinates, offset and flags - const float localX = o2::tpc::Mapper::instance().getPadCentre(PadPos(rowIndex, pad)).X(); - const float localY = Mapper::instance().getPadCentre(PadPos(rowIndex, pad)).Y(); - const float offsPad = (cl.getPad() - pad) * o2::tpc::Mapper::instance().getPadRegionInfo(Mapper::REGION[rowIndex]).getPadWidth(); - const auto flagsCl = cl.getFlags(); + // store the reference data in clusterReferencesByIndex + clusterReferencesByIndex[iCl] = std::make_tuple(sectorIndex, rowIndex, clusterIndexNumb); + } - int excludeCl = 0; // works as a bit mask - if (((clusterMask & ClusterFlags::ExcludeSingleCl) == ClusterFlags::ExcludeSingleCl) && ((flagsCl & ClusterNative::flagSingle) == ClusterNative::flagSingle)) { - excludeCl += 0b001; // 1 for single cluster + // flag groups of several clusters that ended up in the same (sector, row) as eligible for + // dEdxSettings::sameRowClusterMethod to merge but ONLY when they are close in pad and time + for (const auto& [rowKey, clusterIndices] : clustersByRow) { + if (clusterIndices.size() <= 1) { + continue; } - if (((clusterMask & ClusterFlags::ExcludeSplitCl) == ClusterFlags::ExcludeSplitCl) && (((flagsCl & ClusterNative::flagSplitPad) == ClusterNative::flagSplitPad) || ((flagsCl & ClusterNative::flagSplitTime) == ClusterNative::flagSplitTime))) { - excludeCl += 0b010; // 2 for split cluster + + // proximity gate: skip the merge if any pair in the group is farther apart than the pad/time window + float minPad = 1e9f, maxPad = -1e9f, minTime = 1e9f, maxTime = -1e9f; + for (int clusterIdx : clusterIndices) { + const o2::tpc::ClusterNative& cl = track.getCluster(*mTPCTrackClIdxVecInput, clusterIdx, *mClusterIndex); + minPad = std::min(minPad, cl.getPad()); + maxPad = std::max(maxPad, cl.getPad()); + minTime = std::min(minTime, cl.getTime()); + maxTime = std::max(maxTime, cl.getTime()); } - if (((clusterMask & ClusterFlags::ExcludeEdgeCl) == ClusterFlags::ExcludeEdgeCl) && ((flagsCl & ClusterNative::flagEdge) == ClusterNative::flagEdge)) { - excludeCl += 0b100; // 4 for edge cluster + if ((maxPad - minPad) > mSameRowMaxPadDiff || (maxTime - minTime) > mSameRowMaxTimeDiff) { + continue; // looper legs / distinct crossings -> keep as separate samples } - if (((clusterMask & ClusterFlags::ExcludeSharedCl) == ClusterFlags::ExcludeSharedCl) && isShared) { - excludeCl += 0b10000; // for shared cluster + + mergeableRows.insert(rowKey); + } +} + +void CalculatedEdx::handleSameRowClusters(const std::vector& clusters, const ClInfoVec& clusterInfos, std::vector>& rowOrder, std::map, std::vector>& clustersByRow, std::set>& mergeableRows) +{ + const int nClusters = static_cast(clusters.size()); + + // group clusters by (sector, row) + for (int iCl = 0; iCl < nClusters; iCl++) { + const auto rowKey = std::make_pair(clusterInfos[iCl].sectorIndex, clusterInfos[iCl].rowIndex); + if (clustersByRow.find(rowKey) == clustersByRow.end()) { + rowOrder.emplace_back(rowKey); } - // get the x position of the track - const float xPosition = Mapper::instance().getPadCentre(PadPos(rowIndex, 0)).X(); + // add the cluster index to the corresponding (sector, row) key in clustersByRow + clustersByRow[rowKey].emplace_back(iCl); + } + + // flag groups of several clusters in the same (sector, row) as eligible for dEdxSettings::sameRowClusterMethod to merge, only when they are close in pad and time + for (const auto& [rowKey, clusterIndices] : clustersByRow) { + if (clusterIndices.size() <= 1) { + continue; + } - bool check = true; - if (!mPropagateTrack) { - if (mRefit == nullptr) { - LOGP(error, "mRefit is a nullptr, call the function setRefit() before looping over the tracks."); + // proximity gate: skip the merge if any pair in the group is farther apart than the pad/time window + float minPad = 1e9f, maxPad = -1e9f, minTime = 1e9f, maxTime = -1e9f; + for (int clusterIdx : clusterIndices) { + const o2::tpc::ClusterNative& cl = clusters[clusterIdx]; + minPad = std::min(minPad, cl.getPad()); + maxPad = std::max(maxPad, cl.getPad()); + minTime = std::min(minTime, cl.getTime()); + maxTime = std::max(maxTime, cl.getTime()); + } + if ((maxPad - minPad) > mSameRowMaxPadDiff || (maxTime - minTime) > mSameRowMaxTimeDiff) { + continue; // looper legs / distinct crossings -> keep as separate samples + } + + mergeableRows.insert(rowKey); + } +} + +void CalculatedEdx::gatherRowClusterData(o2::tpc::TrackTPC& track, std::vector& rowData, AverageOccupancy& averageOcc) +{ + rowData.clear(); + + bool refitAbandoned = false; + + // handle same (sector, row) clusters + std::vector> rowOrder; + std::map, std::vector> clustersByRow; + std::set> mergeableRows; + std::map> clusterReferencesByIndex; + + handleSameRowClusters(track, rowOrder, clustersByRow, mergeableRows, clusterReferencesByIndex); + + rowData.reserve(rowOrder.size()); + + // per-region occupancy, for the average occupancy output + std::array, 4> occupancyROC; + + // for tracking missing clusters + unsigned char rowIndexOld = 255; + unsigned char sectorIndexOld = 255; + + // loop over the clusters in the track's row-traversal order (rowOrder) + for (const auto& rowKey : rowOrder) { + const auto& clusterIndices = clustersByRow.at(rowKey); + const unsigned char rowIndex = rowKey.second; + const unsigned char sectorIndex = rowKey.first; + + std::vector fragmentClusters; + std::vector fragmentIsShared; + fragmentClusters.reserve(clusterIndices.size()); + fragmentIsShared.reserve(clusterIndices.size()); + for (const int clusterIdx : clusterIndices) { + fragmentClusters.emplace_back(track.getCluster(*mTPCTrackClIdxVecInput, clusterIdx, *mClusterIndex)); + const auto& [fragSectorIndex, fragRowIndex, clusterIndexNumb] = clusterReferencesByIndex[clusterIdx]; + const unsigned int absoluteIndex = mClusterIndex->clusterOffset[fragSectorIndex][fragRowIndex] + clusterIndexNumb; + fragmentIsShared.emplace_back(mRefit ? (mTPCRefitterShMap[absoluteIndex] & o2::gpu::GPUTPCGMMergedTrackHit::flagShared) : 0); + } + + const bool mergeable = mergeableRows.count(rowKey) > 0; + gatherRowClusterDataForRow(track, fragmentClusters, fragmentIsShared, sectorIndex, rowIndex, mergeable, rowIndexOld, sectorIndexOld, occupancyROC, rowData, refitAbandoned); + rowIndexOld = rowIndex; + sectorIndexOld = sectorIndex; + } + + // calculate average cl occupancy for the track per TPC region; skip clusters where getOccupancy() had no data (sentinel -1) + double* const averageOccROC[4] = {&averageOcc.IROC, &averageOcc.OROC1, &averageOcc.OROC2, &averageOcc.OROC3}; + for (int roc = 0; roc < 4; roc++) { + unsigned int sumOcc = 0; + size_t nValidOcc = 0; + for (const unsigned int occ : occupancyROC[roc]) { + if (occ != static_cast(-1)) { + sumOcc += occ; + ++nValidOcc; } + } + if (nValidOcc > 0) { + *averageOccROC[roc] = static_cast(sumOcc) / nValidOcc; + } + } +} + +bool CalculatedEdx::propagateTrackToX(o2::track::TrackParCov& track, float xPosition, unsigned char sectorIndex) const +{ + const o2::track::TrackParCov trackBackup = track; + bool check = track.rotate(o2::math_utils::detail::sector2Angle(sectorIndex)); + if (check) { + check = o2::base::Propagator::Instance()->PropagateToXBxByBz(track, xPosition, 0.999f, 0.5f, o2::base::Propagator::MatCorrType::USEMatCorrLUT); + } + if (!check) { + track = trackBackup; + check = track.rotate(o2::math_utils::detail::sector2Angle(sectorIndex)); + if (check) { + check = o2::base::Propagator::Instance()->PropagateToXBxByBz(track, xPosition, 0.999f, 0.5f, o2::base::Propagator::MatCorrType::USEMatCorrNONE); + } + } + if (!check) { + track = trackBackup; + check = track.rotateParam(o2::math_utils::detail::sector2Angle(sectorIndex)); + if (check) { + check = track.propagateParamTo(xPosition, mFieldNominalGPUBz); + } + } + if (!check) { + track = trackBackup; + } + return check; +} + +void CalculatedEdx::gatherRowClusterDataForRow(o2::tpc::TrackTPC& track, const std::vector& fragmentClusters, const std::vector& fragmentIsShared, unsigned char sectorIndex, unsigned char rowIndex, bool mergeable, unsigned char rowIndexOld, unsigned char sectorIndexOld, std::array, 4>& occupancyROC, std::vector& rowData, bool& refitAbandoned) +{ + RowClusterData row; + row.sectorIndex = sectorIndex; + row.rowIndex = rowIndex; + row.mergeable = mergeable; + + // get region and stack + const int region = Mapper::REGION[rowIndex]; + const CRU cru(Sector(sectorIndex), region); + const auto stack = cru.gemStack(); + StackID stackID{sectorIndex, stack}; + const int stackNumber = static_cast(stack); + + row.region = region; + row.stack = stack; + row.stackID = stackID; + row.stackNumber = stackNumber; + + // per-fragment quantities: pad depends on the individual cluster's own pad, so threshold/gain/gainResidual/occupancy/isDeadRegion + // all looked up by pad, are computed once per fragment here + row.fragments.reserve(fragmentClusters.size()); + for (size_t iFrag = 0; iFrag < fragmentClusters.size(); ++iFrag) { + const o2::tpc::ClusterNative& cl = fragmentClusters[iFrag]; + RowFragment frag; + frag.cl = cl; + frag.isShared = fragmentIsShared[iFrag]; + frag.pad = std::clamp(static_cast(cl.getPad() + 0.5f), static_cast(0), Mapper::PADSPERROW[region][Mapper::getLocalRowFromGlobalRow(rowIndex)] - 1); // the left side of the pad is defined at e.g. 3.5 and the right side at 4.5 + frag.occupancy = getOccupancy(cl.getTime()); + frag.threshold = mCalibCont.getZeroSupressionThreshold(sectorIndex, rowIndex, frag.pad); + frag.gain = mCalibCont.getGain(sectorIndex, rowIndex, frag.pad); + frag.gainResidual = mCalibCont.getResidualGain(sectorIndex, rowIndex, frag.pad); + frag.isDeadRegion = mCalibCont.isDead(static_cast(sectorIndex), static_cast(rowIndex), static_cast(frag.pad)); + + row.fragments.emplace_back(std::move(frag)); + } + + // dEdxSettings::sameRowClusterMethod==1's merged sample + if (mergeable && fragmentClusters.size() > 1) { + row.mergedFragment.cl = buildMergedClusterSum(fragmentClusters); + row.mergedFragment.pad = std::clamp(static_cast(row.mergedFragment.cl.getPad() + 0.5f), static_cast(0), Mapper::PADSPERROW[region][Mapper::getLocalRowFromGlobalRow(rowIndex)] - 1); + row.mergedFragment.occupancy = getOccupancy(row.mergedFragment.cl.getTime()); + row.mergedFragment.threshold = mCalibCont.getZeroSupressionThreshold(sectorIndex, rowIndex, row.mergedFragment.pad); + row.mergedFragment.gain = mCalibCont.getGain(sectorIndex, rowIndex, row.mergedFragment.pad); + row.mergedFragment.gainResidual = mCalibCont.getResidualGain(sectorIndex, rowIndex, row.mergedFragment.pad); + row.mergedFragment.isDeadRegion = mCalibCont.isDead(static_cast(sectorIndex), static_cast(rowIndex), static_cast(row.mergedFragment.pad)); + row.mergedFragment.isShared = row.fragments[0].isShared; // preserves the historical "first fragment's isShared" choice for a merged sample + } + + // occupancy sample(s) for this row's average occupancy contribution + std::array, 4>::size_type occStackIdx; + if (stack == GEMstack::IROCgem) { + occStackIdx = 0; + } else if (stack == GEMstack::OROC1gem) { + occStackIdx = 1; + } else if (stack == GEMstack::OROC2gem) { + occStackIdx = 2; + } else { + occStackIdx = 3; + } + if (mergeable && fragmentClusters.size() > 1) { + occupancyROC[occStackIdx].emplace_back(row.mergedFragment.occupancy); + } else { + for (const auto& frag : row.fragments) { + occupancyROC[occStackIdx].emplace_back(frag.occupancy); + } + } + + // get the x position of the track + const float xPosition = Mapper::instance().getPadCentre(PadPos(rowIndex, 0)).X(); + bool check = true; + bool refitFellBack = false; + if (mRefit) { + if (!refitAbandoned) { + // snapshot the track's state as it stood before this row's refit attempt (i.e. after the previous row's + // successful refit) -- RefitTrackAsGPU() writes back into `track` even on failure, so on failure this is + // the most recent known-good state to fall back from, not the track's pristine pre-loop state + const o2::track::TrackParCov trackBeforeRefit = track; + // refit this track mRefit->setTrackReferenceX(xPosition); - check = (mRefit->RefitTrackAsGPU(track, false, true) < 0) ? false : true; + // RefitTrackAsGPU() returns < 0 when it fails; reachedReference is false when the fit succeeded but the final move-to-reference step could not reach xPosition, so both trigger the fallback. + bool reachedReference = true; + check = (mRefit->RefitTrackAsGPU(track, false, true, &reachedReference) < 0) ? false : reachedReference; + if (!check || std::isnan(track.getParam(1))) { + refitAbandoned = true; + static_cast(track) = trackBeforeRefit; + check = propagateTrackToX(track, xPosition, sectorIndex); + refitFellBack = check; + } } else { - // propagate this track to the plane X=xk (cm) in the field "b" (kG) - track.rotate(o2::math_utils::detail::sector2Angle(sectorIndex)); - check = o2::base::Propagator::Instance()->PropagateToXBxByBz(track, xPosition, 0.999f, 2., o2::base::Propagator::MatCorrType::USEMatCorrLUT); + // already abandoned refit for this track (the previous row's RefitTrackAsGPU() failed) + // keep propagating incrementally from `track`'s current state + check = propagateTrackToX(track, xPosition, sectorIndex); + refitFellBack = check; + } + } else if (mPropagateTrack) { + // propagate this track to the plane X=xk (cm) in the field "b" (kG) + check = propagateTrackToX(track, xPosition, sectorIndex); + } else if (mPropagateParams) { + // propagate the params of the track instead of full propagation; same rollback rationale as mPropagateTrack above + const o2::track::TrackParCov trackBackup = track; + check = track.rotateParam(o2::math_utils::detail::sector2Angle(sectorIndex)); + if (check) { + check = track.propagateParamTo(xPosition, mFieldNominalGPUBz); } + if (!check) { + static_cast(track) = trackBackup; + } + } - if (!check || std::isnan(track.getParam(1))) { - excludeCl += 0b1000; // 8 for failure of track propagation or refit + row.propagationFailed = (!check || std::isnan(track.getParam(1))); + row.refitFellBack = refitFellBack && !row.propagationFailed; + ++mNRowsProcessed; + if (row.propagationFailed) { + ++mNPropagationFailed; + } else if (row.refitFellBack) { + ++mNRefitFallback; + } + + // snapshot of the track state after refit/propagation to this row; reused by calculatedEdxFromRowData() for every settings entry + row.trackSnapshot = track; + + // number of rows skipped between this row and the previous entry in rowOrder + row.sameSectorAsPrevRow = (sectorIndexOld == sectorIndex); + row.missingClusters = (rowIndexOld == 255) ? 0 : (std::abs(static_cast(rowIndex) - static_cast(rowIndexOld)) - 1); + + // veto the gap as a subthreshold candidate if any of its missing row(s) would land on a dead channel or off the padrow edge + row.missingClusterGapDeadOrEdge = false; + if (row.missingClusters > 0 && row.missingClusters <= mMaxMissingCl) { + const o2::gpu::GPUTPCGeometry gpuGeom; + const RowClusterData& prevRow = rowData.back(); + // bracket the gap by its lower/upper real row, independent of the rowOrder direction + const int rowLo = std::min(rowIndex, rowIndexOld); + const int rowHi = std::max(rowIndex, rowIndexOld); + const float padLo = (rowLo == static_cast(rowIndexOld)) ? prevRow.fragments[0].cl.getPad() : row.fragments[0].cl.getPad(); + const float padHi = (rowHi == static_cast(rowIndexOld)) ? prevRow.fragments[0].cl.getPad() : row.fragments[0].cl.getPad(); + const float yLo = gpuGeom.LinearPad2Y(sectorIndex, rowLo, padLo); + const float yHi = gpuGeom.LinearPad2Y(sectorIndex, rowHi, padHi); + for (int k = 1; k <= row.missingClusters; ++k) { + const unsigned char missingRow = static_cast(rowLo + k); + const float frac = static_cast(k) / (row.missingClusters + 1); + const float missingPad = gpuGeom.LinearY2Pad(sectorIndex, missingRow, yLo + (yHi - yLo) * frac); + if (missingPad < 0.f || missingPad >= gpuGeom.NPads(missingRow)) { + row.missingClusterGapDeadOrEdge = true; + break; + } + const int missingRegion = Mapper::REGION[missingRow]; + const unsigned char missingPadClamped = std::clamp(static_cast(missingPad + 0.5f), static_cast(0), Mapper::PADSPERROW[missingRegion][Mapper::getLocalRowFromGlobalRow(missingRow)] - 1); + if (mCalibCont.isDead(static_cast(sectorIndex), static_cast(missingRow), static_cast(missingPadClamped))) { + row.missingClusterGapDeadOrEdge = true; + break; + } } + } - if (excludeCl != 0) { - // for debugging - if (mDebug) { - excludeClVector.emplace_back(excludeCl); - regionVector.emplace_back(region); - rowIndexVector.emplace_back(rowIndex); - padVector.emplace_back(pad); - sectorVector.emplace_back(sectorIndex); - stackVector.emplace_back(stackNumber); - localXVector.emplace_back(localX); - localYVector.emplace_back(localY); - offsPadVector.emplace_back(offsPad); - trackVector.emplace_back(track); - clVector.emplace_back(cl); - occupancyVector.emplace_back(getOccupancy(cl)); - isClusterShared.emplace_back(isShared); - - topologyCorrVector.emplace_back(-999.f); - topologyCorrTotVector.emplace_back(-999.f); - topologyCorrMaxVector.emplace_back(-999.f); - gainVector.emplace_back(-999.f); - gainResidualVector.emplace_back(-999.f); - residualCorrTotVector.emplace_back(-999.f); - residualCorrMaxVector.emplace_back(-999.f); - scCorrVector.emplace_back(-999.f); + rowData.emplace_back(std::move(row)); +} + +void CalculatedEdx::gatherRowClusterData(o2::tpc::TrackTPC& track, const std::vector& clusters, const ClInfoVec& clusterInfos, std::vector& rowData, AverageOccupancy& averageOcc) +{ + rowData.clear(); + + bool refitAbandoned = false; + + // handle same (sector, row) clusters + std::vector> rowOrder; + std::map, std::vector> clustersByRow; + std::set> mergeableRows; + + handleSameRowClusters(clusters, clusterInfos, rowOrder, clustersByRow, mergeableRows); + + rowData.reserve(rowOrder.size()); + + // per-region occupancy, for the average occupancy output + std::array, 4> occupancyROC; + + // for tracking missing clusters + unsigned char rowIndexOld = 255; + unsigned char sectorIndexOld = 255; + + // loop over the clusters in the rowOrder; the refit/propagation to a row is done exactly once here regardless of how many raw clusters the row has + for (const auto& rowKey : rowOrder) { + const auto& clusterIndices = clustersByRow.at(rowKey); + const unsigned char rowIndex = rowKey.second; + const unsigned char sectorIndex = rowKey.first; + + std::vector fragmentClusters; + std::vector fragmentIsShared; + fragmentClusters.reserve(clusterIndices.size()); + fragmentIsShared.reserve(clusterIndices.size()); + for (const int clusterIdx : clusterIndices) { + fragmentClusters.emplace_back(clusters[clusterIdx]); + // isShared cannot be looked up from mTPCRefitterShMap for externally supplied clusters, so it is taken from the info directly + fragmentIsShared.emplace_back(clusterInfos[clusterIdx].isShared); + } + + const bool mergeable = mergeableRows.count(rowKey) > 0; + gatherRowClusterDataForRow(track, fragmentClusters, fragmentIsShared, sectorIndex, rowIndex, mergeable, rowIndexOld, sectorIndexOld, occupancyROC, rowData, refitAbandoned); + rowData.back().inputClusterIndices = clusterIndices; // positions in the externally supplied clusters vector grouped into this row + rowIndexOld = rowIndex; + sectorIndexOld = sectorIndex; + } + + // calculate average cl occupancy for the track per TPC region; skip clusters where getOccupancy() had no data (sentinel -1) + double* const averageOccROC[4] = {&averageOcc.IROC, &averageOcc.OROC1, &averageOcc.OROC2, &averageOcc.OROC3}; + for (int roc = 0; roc < 4; roc++) { + unsigned int sumOcc = 0; + size_t nValidOcc = 0; + for (const unsigned int occ : occupancyROC[roc]) { + if (occ != static_cast(-1)) { + sumOcc += occ; + ++nValidOcc; } - // to avoid counting the skipped cluster as a subthreshold cluster - rowIndexOld = rowIndex; - sectorIndexOld = sectorIndex; - continue; } + if (nValidOcc > 0) { + *averageOccROC[roc] = static_cast(sumOcc) / nValidOcc; + } + } +} - // get charge values - float chargeTot = cl.getQtot(); - float chargeMax = cl.qMax; +void CalculatedEdx::calculatedEdxFromRowData(const std::vector& rowData, const dEdxSettings& settings, size_t settingsIndex, float trackTime0, const o2::tpc::TrackTPC& trackOrig, const AverageOccupancy& averageOcc, dEdxInfo& output, const MCCompLabel* mcLabel) +{ + // NHits and NHitsSubthreshold values per region + int nClsROC[4] = {0, 0, 0, 0}; + int nClsSubThreshROC[4] = {0, 0, 0, 0}; - // get threshold - const float threshold = mCalibCont.getZeroSupressionThreshold(sectorIndex, rowIndex, pad); + const unsigned short sameRowClusterMethod = (settings.sameRowClusterMethod <= 2) ? settings.sameRowClusterMethod : 0; + if (settings.sameRowClusterMethod > 2) { + LOGP(warning, "Unrecognized sameRowClusterMethod {} (expected 0, 1, or 2); treating same-row cluster groups as method 0 (do not merge)", settings.sameRowClusterMethod); + } + + // corrected qTot and qMax values per region + const int nType = 5; + std::array, nType> chargeTotROC; + std::array, nType> chargeMaxROC; + for (int i = 0; i < nType; ++i) { + chargeTotROC[i].reserve(Mapper::PADROWS); + chargeMaxROC[i].reserve(Mapper::PADROWS); + } + + // per-region (IROC, OROC1, OROC2, OROC3) running minimum charge among accepted clusters, used as the virtual charge for that region's subthreshold clusters below + float minChargeTotROC[4] = {kNoValidCharge, kNoValidCharge, kNoValidCharge, kNoValidCharge}; + float minChargeMaxROC[4] = {kNoValidCharge, kNoValidCharge, kNoValidCharge, kNoValidCharge}; + + o2::utils::TreeStreamRedirector* debugStreamer = nullptr; + std::vector occupancyVector; + if (mDebug) { + setStreamer(settings.debugRootFile.c_str()); + debugStreamer = mStreamers.at(settings.debugRootFile).get(); + ++mDebugTrackIndex; + occupancyVector.reserve(rowData.size()); + } + + // a gap is not filled as a subthreshold cluster when the row that closes it sits within the outermost min(nRows/2, mSubThreshEdgeRows) rows + const int edgeRowCut = std::min(static_cast(rowData.size()) / 2, mSubThreshEdgeRows); + + for (size_t iRowData = 0; iRowData < rowData.size(); ++iRowData) { + const auto& row = rowData[iRowData]; + + // one effective sample per row for this settings entry, per dEdxSettings::sameRowClusterMethod; pointers into + // row.fragments/row.mergedFragment (both owned by rowData, alive for the whole calculatedEdxFromRowData() call) + std::vector samples; + + const bool doMerge = row.mergeable && row.fragments.size() > 1 && sameRowClusterMethod != 0; + if (!doMerge) { + samples.reserve(row.fragments.size()); + for (const auto& frag : row.fragments) { + samples.push_back(&frag); + } + } else if (sameRowClusterMethod == 2) { + samples.push_back(&pickDominantFragment(row.fragments)); + } else { // sameRowClusterMethod == 1: use the sum-merged sample gatherRowClusterDataForRow() already computed once for this row + samples.push_back(&row.mergedFragment); + } + + // ExcludeSamePadRowCl: true whenever the row-group HAD >1 raw fragment, regardless of whether this settings entry actually merged them + const bool isCombined = row.fragments.size() > 1; // find missing clusters - int missingClusters = rowIndexOld - rowIndex - 1; - if ((missingClusters > 0) && (missingClusters <= mMaxMissingCl)) { - if ((clusterMask & ClusterFlags::ExcludeSectorBoundaries) == ClusterFlags::ExcludeSectorBoundaries) { - if (sectorIndexOld == sectorIndex) { - if (stack == GEMstack::IROCgem) { + const int missingClusters = row.missingClusters; + if ((missingClusters > 0) && (missingClusters <= mMaxMissingCl) && !row.missingClusterGapDeadOrEdge && (static_cast(iRowData) >= edgeRowCut)) { + if ((settings.clusterMask & ClusterFlags::ExcludeSectorBoundaries) == ClusterFlags::ExcludeSectorBoundaries) { + if (row.sameSectorAsPrevRow) { + if (row.stack == GEMstack::IROCgem) { nClsSubThreshROC[0] += missingClusters; nClsROC[0] += missingClusters; - } else if (stack == GEMstack::OROC1gem) { + } else if (row.stack == GEMstack::OROC1gem) { nClsSubThreshROC[1] += missingClusters; nClsROC[1] += missingClusters; - } else if (stack == GEMstack::OROC2gem) { + } else if (row.stack == GEMstack::OROC2gem) { nClsSubThreshROC[2] += missingClusters; nClsROC[2] += missingClusters; - } else if (stack == GEMstack::OROC3gem) { + } else if (row.stack == GEMstack::OROC3gem) { nClsSubThreshROC[3] += missingClusters; nClsROC[3] += missingClusters; } } } else { - if (stack == GEMstack::IROCgem) { + if (row.stack == GEMstack::IROCgem) { nClsSubThreshROC[0] += missingClusters; nClsROC[0] += missingClusters; - } else if (stack == GEMstack::OROC1gem) { + } else if (row.stack == GEMstack::OROC1gem) { nClsSubThreshROC[1] += missingClusters; nClsROC[1] += missingClusters; - } else if (stack == GEMstack::OROC2gem) { + } else if (row.stack == GEMstack::OROC2gem) { nClsSubThreshROC[2] += missingClusters; nClsROC[2] += missingClusters; - } else if (stack == GEMstack::OROC3gem) { + } else if (row.stack == GEMstack::OROC3gem) { nClsSubThreshROC[3] += missingClusters; nClsROC[3] += missingClusters; } } }; - rowIndexOld = rowIndex; - sectorIndexOld = sectorIndex; - // get effective length - float effectiveLength = 1.0f; - float effectiveLengthTot = 1.0f; - float effectiveLengthMax = 1.0f; - if ((correctionMask & CorrectionFlags::TopologySimple) == CorrectionFlags::TopologySimple) { - effectiveLength = getTrackTopologyCorrection(track, region); - chargeTot /= effectiveLength; - chargeMax /= effectiveLength; - }; - if ((correctionMask & CorrectionFlags::TopologyPol) == CorrectionFlags::TopologyPol) { - effectiveLengthTot = getTrackTopologyCorrectionPol(track, cl, region, chargeTot, ChargeType::Tot, threshold); - effectiveLengthMax = getTrackTopologyCorrectionPol(track, cl, region, chargeMax, ChargeType::Max, threshold); - chargeTot /= effectiveLengthTot; - chargeMax /= effectiveLengthMax; - }; + for (const auto& sample : samples) { + if (mDebug) { + occupancyVector.emplace_back(sample->occupancy); + } - // get gain - float gain = 1.0f; - float gainResidual = 1.0f; - if ((correctionMask & CorrectionFlags::GainFull) == CorrectionFlags::GainFull) { - gain = mCalibCont.getGain(sectorIndex, rowIndex, pad); - }; - if ((correctionMask & CorrectionFlags::GainResidual) == CorrectionFlags::GainResidual) { - gainResidual = mCalibCont.getResidualGain(sectorIndex, rowIndex, pad); - }; - chargeTot /= gain * gainResidual; - chargeMax /= gain * gainResidual; - - // get dEdx correction on tgl and sector plane - float corrTot = 1.0f; - float corrMax = 1.0f; - if ((correctionMask & CorrectionFlags::dEdxResidual) == CorrectionFlags::dEdxResidual) { - corrTot = mCalibCont.getResidualCorrection(stackID, ChargeType::Tot, track.getTgl(), track.getSnp()); - corrMax = mCalibCont.getResidualCorrection(stackID, ChargeType::Max, track.getTgl(), track.getSnp()); - if (corrTot > 0) { - chargeTot /= corrTot; + // get cluster values + float chargeTot = sample->cl.getQtot(); + float chargeMax = sample->cl.getQmax(); + + // corrections + float effectiveLength = 1.0f; + float effectiveLengthTot = 1.0f; + float effectiveLengthMax = 1.0f; + float gain = 1.0f; + float gainResidual = 1.0f; + float corrTot = 1.0f; + float corrMax = 1.0f; + float scCorr = 1.0f; + + int excludeCl = 0; // works as a bit mask + const uint8_t flagsCl = sample->cl.getFlags(); + if (((settings.clusterMask & ClusterFlags::ExcludeSingleCl) == ClusterFlags::ExcludeSingleCl) && ((flagsCl & ClusterNative::flagSingle) == ClusterNative::flagSingle)) { + excludeCl += 0b001; // 1 for single cluster + } + if (((settings.clusterMask & ClusterFlags::ExcludeSplitPadCl) == ClusterFlags::ExcludeSplitPadCl) && ((flagsCl & ClusterNative::flagSplitPad) == ClusterNative::flagSplitPad)) { + excludeCl += 0b010; // 2 for split pad cluster + } + if (((settings.clusterMask & ClusterFlags::ExcludeSplitTimeCl) == ClusterFlags::ExcludeSplitTimeCl) && ((flagsCl & ClusterNative::flagSplitTime) == ClusterNative::flagSplitTime)) { + excludeCl += 0b0100; // 4 for split time cluster + } + if (((settings.clusterMask & ClusterFlags::ExcludeSplitCl) == ClusterFlags::ExcludeSplitCl) && (((flagsCl & ClusterNative::flagSplitPad) == ClusterNative::flagSplitPad) || ((flagsCl & ClusterNative::flagSplitTime) == ClusterNative::flagSplitTime))) { + excludeCl += 0b01000; // 8 for split cluster + } + if (((settings.clusterMask & ClusterFlags::ExcludeEdgeCl) == ClusterFlags::ExcludeEdgeCl) && ((flagsCl & ClusterNative::flagEdge) == ClusterNative::flagEdge)) { + excludeCl += 0b010000; // 16 for edge cluster + } + if (((settings.clusterMask & ClusterFlags::ExcludeSharedCl) == ClusterFlags::ExcludeSharedCl) && sample->isShared) { + excludeCl += 0b0100000; // 32 for shared cluster + } + if (((settings.clusterMask & ClusterFlags::ExcludeSamePadRowCl) == ClusterFlags::ExcludeSamePadRowCl) && isCombined) { + excludeCl += 0b01000000; // 64 for combined cluster + } + if ((settings.stackBoundaryMethod == 1 || settings.stackBoundaryMethod == 2) && isInStackBoundaries(row.stackNumber, row.rowIndex, settings.stackBoundaryMethod)) { + excludeCl += 0b010000000; // 128 for stack boundary cluster + } + if (sample->isDeadRegion) { + excludeCl += 0b0100000000; // 256 for dead region + } + if (row.propagationFailed) { + excludeCl += 0b01000000000; // 512 for failure of track propagation or refit + } + + // get effective length + if ((settings.correctionMask & CorrectionFlags::TopologySimple) == CorrectionFlags::TopologySimple) { + effectiveLength = getTrackTopologyCorrection(row.trackSnapshot, row.region); + chargeTot /= effectiveLength; + chargeMax /= effectiveLength; }; - if (corrMax > 0) { - chargeMax /= corrMax; + + const bool gainFullApplied = (settings.correctionMask & CorrectionFlags::GainFull) == CorrectionFlags::GainFull; + float topoChargeTot = chargeTot; + float topoChargeMax = chargeMax; + if (gainFullApplied) { + gain = sample->gain; + chargeTot /= gain; + chargeMax /= gain; + } else { + topoChargeTot *= sample->gain; + topoChargeMax *= sample->gain; + } + + // topology correction + if ((settings.correctionMask & CorrectionFlags::TopologyPol) == CorrectionFlags::TopologyPol) { + effectiveLengthTot = getTrackTopologyCorrectionPol(row.trackSnapshot, sample->cl, row.region, topoChargeTot, ChargeType::Tot, sample->threshold); + effectiveLengthMax = getTrackTopologyCorrectionPol(row.trackSnapshot, sample->cl, row.region, topoChargeMax, ChargeType::Max, sample->threshold); + chargeTot /= effectiveLengthTot; + chargeMax /= effectiveLengthMax; }; - }; - // set the min charge - if (chargeTot < minChargeTot) { - minChargeTot = chargeTot; - }; + // residual dE/dx correction on tgl and sector plane + if ((settings.correctionMask & CorrectionFlags::dEdxResidual) == CorrectionFlags::dEdxResidual) { + corrTot = mCalibCont.getResidualCorrection(row.stackID, ChargeType::Tot, row.trackSnapshot.getTgl(), row.trackSnapshot.getSnp()); + corrMax = mCalibCont.getResidualCorrection(row.stackID, ChargeType::Max, row.trackSnapshot.getTgl(), row.trackSnapshot.getSnp()); + if (corrTot > 0) { + chargeTot /= corrTot; + }; + if (corrMax > 0) { + chargeMax /= corrMax; + }; + }; - if (chargeMax < minChargeMax) { - minChargeMax = chargeMax; - }; + // residual gain map + if ((settings.correctionMask & CorrectionFlags::GainResidual) == CorrectionFlags::GainResidual) { + gainResidual = sample->gainResidual; + chargeTot /= gainResidual; + chargeMax /= gainResidual; + }; - // space-charge dEdx corrections - const float time = cl.getTime() - track.getTime0(); // ToDo: get correct time from ITS-TPC track if possible - float scCorr = 1.0f; - if ((correctionMask & CorrectionFlags::dEdxSC) == CorrectionFlags::dEdxSC) { - scCorr = mSCdEdxCorrection.getCorrection(time, sectorIndex, rowIndex, pad); - if (scCorr > 0) { - chargeTot /= scCorr; + // space-charge dEdx corrections + const float time = sample->cl.getTime() - trackTime0; // ToDo: get correct time from ITS-TPC track if possible + if ((settings.correctionMask & CorrectionFlags::dEdxSC) == CorrectionFlags::dEdxSC) { + scCorr = mSCdEdxCorrection.getCorrection(time, row.sectorIndex, row.rowIndex, sample->pad); + if (scCorr > 0) { + chargeTot /= scCorr; + }; + if (scCorr > 0) { + chargeMax /= scCorr; + }; + } + + // for debugging + if (mDebug) { + const o2::gpu::GPUTPCGeometry gpuGeom; + const float localX = gpuGeom.Row2X(row.rowIndex); + const float localY = gpuGeom.LinearPad2Y(row.sectorIndex, row.rowIndex, sample->cl.getPad()); + const LocalPosition2D l2D{localX, localY}; + const auto g2D = Mapper::LocalToGlobal(l2D, Sector(row.sectorIndex)); + const float globalX = g2D.x(); + const float globalY = g2D.y(); + + // slice to the base parametrization (X, alpha, params, covariance) instead of the full TrackTPC, since only the parametrization changes cluster-to-cluster after refit/propagation + const o2::track::TrackParCov trackParam = row.trackSnapshot; + + (*debugStreamer) << "dEdxDebugCl" + << "trackIndex=" << mDebugTrackIndex + << "trackParam=" << trackParam + << "cl=" << sample->cl + << "chargeTot=" << chargeTot + << "chargeMax=" << chargeMax + << "excludeCl=" << excludeCl + << "region=" << row.region + << "rowIndex=" << row.rowIndex + << "sectorIndex=" << row.sectorIndex + << "stack=" << row.stackNumber + << "localX=" << localX + << "localY=" << localY + << "globalX=" << globalX + << "globalY=" << globalY + << "isShared=" << sample->isShared + << "isCombined=" << isCombined + << "refitFellBack=" << row.refitFellBack + << "topologyCorr=" << effectiveLength + << "topologyCorrTot=" << effectiveLengthTot + << "topologyCorrMax=" << effectiveLengthMax + << "gain=" << gain + << "gainResidual=" << gainResidual + << "residualCorrTot=" << corrTot + << "residualCorrMax=" << corrMax + << "scCorr=" << scCorr + << "occupancy=" << sample->occupancy + << "inputClusterIndices=" << row.inputClusterIndices + << "\n"; }; - if (corrMax > 0) { - chargeMax /= scCorr; + + if (excludeCl != 0) { + continue; + } + + // set the region's min charge, only from clusters actually included in the dEdx calculation, + // so excluded clusters (dead region, edge, failed propagation, ...) don't bias the virtual charge used for subthreshold filling + if (chargeTot < minChargeTotROC[row.stackNumber]) { + minChargeTotROC[row.stackNumber] = chargeTot; }; - } - if (stack == GEMstack::IROCgem) { - chargeTotROC[0].emplace_back(chargeTot); - chargeMaxROC[0].emplace_back(chargeMax); - nClsROC[0]++; - } else if (stack == GEMstack::OROC1gem) { - chargeTotROC[1].emplace_back(chargeTot); - chargeMaxROC[1].emplace_back(chargeMax); - nClsROC[1]++; - } else if (stack == GEMstack::OROC2gem) { - chargeTotROC[2].emplace_back(chargeTot); - chargeMaxROC[2].emplace_back(chargeMax); - nClsROC[2]++; - } else if (stack == GEMstack::OROC3gem) { - chargeTotROC[3].emplace_back(chargeTot); - chargeMaxROC[3].emplace_back(chargeMax); - nClsROC[3]++; - }; + if (chargeMax < minChargeMaxROC[row.stackNumber]) { + minChargeMaxROC[row.stackNumber] = chargeMax; + }; - chargeTotROC[4].emplace_back(chargeTot); - chargeMaxROC[4].emplace_back(chargeMax); - - // for debugging - if (mDebug) { - excludeClVector.emplace_back(0); // cl is successfully processed - regionVector.emplace_back(region); - rowIndexVector.emplace_back(rowIndex); - padVector.emplace_back(pad); - sectorVector.emplace_back(sectorIndex); - stackVector.emplace_back(stackNumber); - localXVector.emplace_back(localX); - localYVector.emplace_back(localY); - offsPadVector.emplace_back(offsPad); - trackVector.emplace_back(track); - clVector.emplace_back(cl); - occupancyVector.emplace_back(getOccupancy(cl)); - isClusterShared.emplace_back(isShared); - - topologyCorrVector.emplace_back(effectiveLength); - topologyCorrTotVector.emplace_back(effectiveLengthTot); - topologyCorrMaxVector.emplace_back(effectiveLengthMax); - gainVector.emplace_back(gain); - gainResidualVector.emplace_back(gainResidual); - residualCorrTotVector.emplace_back(corrTot); - residualCorrMaxVector.emplace_back(corrMax); - scCorrVector.emplace_back(scCorr); - }; + if (row.stack == GEMstack::IROCgem) { + chargeTotROC[0].emplace_back(chargeTot); + chargeMaxROC[0].emplace_back(chargeMax); + nClsROC[0]++; + } else if (row.stack == GEMstack::OROC1gem) { + chargeTotROC[1].emplace_back(chargeTot); + chargeMaxROC[1].emplace_back(chargeMax); + nClsROC[1]++; + } else if (row.stack == GEMstack::OROC2gem) { + chargeTotROC[2].emplace_back(chargeTot); + chargeMaxROC[2].emplace_back(chargeMax); + nClsROC[2]++; + } else if (row.stack == GEMstack::OROC3gem) { + chargeTotROC[3].emplace_back(chargeTot); + chargeMaxROC[3].emplace_back(chargeMax); + nClsROC[3]++; + }; + + chargeTotROC[4].emplace_back(chargeTot); + chargeMaxROC[4].emplace_back(chargeMax); + } + } + + // fill subthreshold clusters if not excluded + if (((settings.clusterMask & ClusterFlags::ExcludeSubthresholdCl) == ClusterFlags::None)) { + float cappedMinChargeTotROC[4], cappedMinChargeMaxROC[4]; + for (int roc = 0; roc < 4; roc++) { + cappedMinChargeTotROC[roc] = (minChargeTotROC[roc] >= kNoValidCharge) ? minChargeTotROC[roc] : std::min(minChargeTotROC[roc], settings.maxSubthresholdChargeTot); + cappedMinChargeMaxROC[roc] = (minChargeMaxROC[roc] >= kNoValidCharge) ? minChargeMaxROC[roc] : std::min(minChargeMaxROC[roc], settings.maxSubthresholdChargeMax); + // a ROC with no valid accepted-cluster charge at all makes fillMissingClusters() below skip it entirely + // (nothing is pushed into chargeTotROC/chargeMaxROC for it) -- so the gaps counted for this ROC earlier + // in the row loop must be un-counted here too, or NHits*/NHitsSubThreshold* and the + // mNSubThresholdFilledPerSettings diagnostic would report fills that never actually happened + if (minChargeTotROC[roc] >= kNoValidCharge || minChargeMaxROC[roc] >= kNoValidCharge) { + nClsROC[roc] -= nClsSubThreshROC[roc]; + nClsSubThreshROC[roc] = 0; + } + } + fillMissingClusters(nClsSubThreshROC, cappedMinChargeTotROC, cappedMinChargeMaxROC, settings.subthresholdMethod, chargeTotROC, chargeMaxROC); + if (mNSubThresholdFilledPerSettings.size() <= settingsIndex) { + mNSubThresholdFilledPerSettings.resize(settingsIndex + 1, 0); + } + mNSubThresholdFilledPerSettings[settingsIndex] += nClsSubThreshROC[0] + nClsSubThreshROC[1] + nClsSubThreshROC[2] + nClsSubThreshROC[3]; } // number of clusters @@ -407,77 +880,152 @@ void CalculatedEdx::calculatedEdx(o2::tpc::TrackTPC& track, dEdxInfo& output, fl output.NHitsSubThresholdOROC2 = nClsROC[2]; output.NHitsSubThresholdOROC3 = nClsROC[3]; - // check if the lost clusters are subthreshold clusters based on the charge thresholds - if (minChargeTot <= mMinChargeTotThreshold && minChargeMax <= mMinChargeMaxThreshold) { - output.NHitsIROC = nClsROC[0] - nClsSubThreshROC[0]; - output.NHitsOROC1 = nClsROC[1] - nClsSubThreshROC[1]; - output.NHitsOROC2 = nClsROC[2] - nClsSubThreshROC[2]; - output.NHitsOROC3 = nClsROC[3] - nClsSubThreshROC[3]; - - // fill subthreshold clusters if not excluded - if (((clusterMask & ClusterFlags::ExcludeSubthresholdCl) == ClusterFlags::None)) { - fillMissingClusters(nClsSubThreshROC, minChargeTot, minChargeMax, subthresholdMethod, chargeTotROC, chargeMaxROC); - } - } else { - output.NHitsIROC = nClsROC[0]; - output.NHitsOROC1 = nClsROC[1]; - output.NHitsOROC2 = nClsROC[2]; - output.NHitsOROC3 = nClsROC[3]; - } + // the gaps found above are always treated as subthreshold clusters (except a ROC with no valid charge to fill + // them with at all, un-counted above so it isn't double-reported as both "hit" and "subthreshold hit") + output.NHitsIROC = nClsROC[0] - nClsSubThreshROC[0]; + output.NHitsOROC1 = nClsROC[1] - nClsSubThreshROC[1]; + output.NHitsOROC2 = nClsROC[2] - nClsSubThreshROC[2]; + output.NHitsOROC3 = nClsROC[3] - nClsSubThreshROC[3]; // copy corrected cluster charges auto chargeTotVector = mDebug ? chargeTotROC[4] : std::vector(); auto chargeMaxVector = mDebug ? chargeMaxROC[4] : std::vector(); // calculate dEdx - output.dEdxTotIROC = getTruncMean(chargeTotROC[0], low, high); - output.dEdxTotOROC1 = getTruncMean(chargeTotROC[1], low, high); - output.dEdxTotOROC2 = getTruncMean(chargeTotROC[2], low, high); - output.dEdxTotOROC3 = getTruncMean(chargeTotROC[3], low, high); - output.dEdxTotTPC = getTruncMean(chargeTotROC[4], low, high); - - output.dEdxMaxIROC = getTruncMean(chargeMaxROC[0], low, high); - output.dEdxMaxOROC1 = getTruncMean(chargeMaxROC[1], low, high); - output.dEdxMaxOROC2 = getTruncMean(chargeMaxROC[2], low, high); - output.dEdxMaxOROC3 = getTruncMean(chargeMaxROC[3], low, high); - output.dEdxMaxTPC = getTruncMean(chargeMaxROC[4], low, high); - - // for debugging + output.dEdxTotIROC = getTruncMean(chargeTotROC[0], settings.low, settings.high); + output.dEdxTotOROC1 = getTruncMean(chargeTotROC[1], settings.low, settings.high); + output.dEdxTotOROC2 = getTruncMean(chargeTotROC[2], settings.low, settings.high); + output.dEdxTotOROC3 = getTruncMean(chargeTotROC[3], settings.low, settings.high); + output.dEdxTotTPC = getTruncMean(chargeTotROC[4], settings.low, settings.high); + + output.dEdxMaxIROC = getTruncMean(chargeMaxROC[0], settings.low, settings.high); + output.dEdxMaxOROC1 = getTruncMean(chargeMaxROC[1], settings.low, settings.high); + output.dEdxMaxOROC2 = getTruncMean(chargeMaxROC[2], settings.low, settings.high); + output.dEdxMaxOROC3 = getTruncMean(chargeMaxROC[3], settings.low, settings.high); + output.dEdxMaxTPC = getTruncMean(chargeMaxROC[4], settings.low, settings.high); + + // for debugging: one row per track, with the track as it was before refit/propagation touched it, per-cluster rows were already written to the "dEdxDebugCl" tree above (each with its own propagated track parameters) and can be grouped back to this row via trackIndex + if (mDebug) { + float minChargeTot = minChargeTotROC[0], minChargeMax = minChargeMaxROC[0]; + for (int roc = 1; roc < 4; roc++) { + minChargeTot = (minChargeTotROC[roc] < minChargeTot) ? minChargeTotROC[roc] : minChargeTot; + minChargeMax = (minChargeMaxROC[roc] < minChargeMax) ? minChargeMaxROC[roc] : minChargeMax; + } + const MCCompLabel label = mcLabel ? *mcLabel : MCCompLabel{}; + (*debugStreamer) << "dEdxDebugTrack" + << "trackIndex=" << mDebugTrackIndex + << "track=" << trackOrig + << "output=" << output + << "averageOcc=" << averageOcc + << "nCl=" << rowData.size() + << "minChargeTot=" << minChargeTot + << "minChargeMax=" << minChargeMax + << "chargeTotVector=" << chargeTotVector + << "chargeMaxVector=" << chargeMaxVector + << "occupancy=" << occupancyVector + << "mcLabel=" << label + << "\n"; + } +} + +void CalculatedEdx::calculatedEdxMultipleSettings(o2::tpc::TrackTPC& track, std::vector& outputs, AverageOccupancy& averageOcc, const std::vector& settingsList, const MCCompLabel* mcLabel) +{ + outputs.clear(); + if (settingsList.empty()) { + return; + } + + o2::tpc::TrackTPC trackOrig; + if (mDebug) { + trackOrig = track; // pristine track, before refit/propagation mutates it cluster-by-cluster below + } + const float trackTime0 = track.getTime0(); // unaffected by refit/propagation, so it is the same for every row and every settings entry + + // gather the per-row cluster/track data once, performing the refit/propagation to each cluster row exactly once; this also fills averageOcc, which does not depend on the dEdx settings and is therefore shared by every settings entry + std::vector rowData; + gatherRowClusterData(track, rowData, averageOcc); + + // evaluate each settings entry against the shared row data + outputs.resize(settingsList.size()); + for (size_t i = 0; i < settingsList.size(); ++i) { + calculatedEdxFromRowData(rowData, settingsList[i], i, trackTime0, trackOrig, averageOcc, outputs[i], mcLabel); + } +} + +void CalculatedEdx::calculatedEdx(o2::tpc::TrackTPC& track, const std::vector& clusters, const ClInfoVec& clusterInfos, dEdxInfo& output, AverageOccupancy& averageOcc, float low, float high, CorrectionFlags correctionMask, ClusterFlags clusterMask, int subthresholdMethod, int stackBoundaryMethod, const char* debugRootFile, float maxSubthresholdChargeTot, float maxSubthresholdChargeMax, int sameRowClusterMethod) +{ + dEdxSettings settings; + settings.low = low; + settings.high = high; + settings.correctionMask = correctionMask; + settings.clusterMask = clusterMask; + settings.subthresholdMethod = subthresholdMethod; + settings.stackBoundaryMethod = stackBoundaryMethod; + settings.debugRootFile = debugRootFile; + settings.maxSubthresholdChargeTot = maxSubthresholdChargeTot; + settings.maxSubthresholdChargeMax = maxSubthresholdChargeMax; + settings.sameRowClusterMethod = sameRowClusterMethod; + + o2::tpc::TrackTPC trackOrig; if (mDebug) { - if (mStreamer == nullptr) { - setStreamer(debugRootFile); - } - - (*mStreamer) << "dEdxDebug" - << "Ncl=" << nClusters - << "excludeClVector=" << excludeClVector - << "regionVector=" << regionVector - << "rowIndexVector=" << rowIndexVector - << "padVector=" << padVector - << "sectorVector=" << sectorVector - << "stackVector=" << stackVector - << "topologyCorrVector=" << topologyCorrVector - << "topologyCorrTotVector=" << topologyCorrTotVector - << "topologyCorrMaxVector=" << topologyCorrMaxVector - << "gainVector=" << gainVector - << "gainResidualVector=" << gainResidualVector - << "residualCorrTotVector=" << residualCorrTotVector - << "residualCorrMaxVector=" << residualCorrMaxVector - << "scCorrVector=" << scCorrVector - << "localXVector=" << localXVector - << "localYVector=" << localYVector - << "offsPadVector=" << offsPadVector - << "trackVector=" << trackVector - << "clVector=" << clVector - << "minChargeTot=" << minChargeTot - << "minChargeMax=" << minChargeMax - << "output=" << output - << "occupancy=" << occupancyVector - << "chargeTotVector=" << chargeTotVector - << "chargeMaxVector=" << chargeMaxVector - << "isClusterShared=" << isClusterShared - << "\n"; + trackOrig = track; // pristine track, before refit/propagation mutates it cluster-by-cluster below } + const float trackTime0 = track.getTime0(); + + std::vector rowData; + gatherRowClusterData(track, clusters, clusterInfos, rowData, averageOcc); + + calculatedEdxFromRowData(rowData, settings, 0, trackTime0, trackOrig, averageOcc, output); +} + +void CalculatedEdx::calculatedEdxMultipleSettings(o2::tpc::TrackTPC& track, const std::vector& clusters, const ClInfoVec& clusterInfos, std::vector& outputs, AverageOccupancy& averageOcc, const std::vector& settingsList, const MCCompLabel* mcLabel) +{ + outputs.clear(); + if (settingsList.empty()) { + return; + } + + o2::tpc::TrackTPC trackOrig; + if (mDebug) { + trackOrig = track; // pristine track, before refit/propagation mutates it cluster-by-cluster below + } + const float trackTime0 = track.getTime0(); // unaffected by refit/propagation, so it is the same for every row and every settings entry + + // gather the per-row cluster/track data once, performing the refit/propagation to each cluster row exactly once; this also fills averageOcc, which does not depend on the dEdx settings and is therefore shared by every settings entry + std::vector rowData; + gatherRowClusterData(track, clusters, clusterInfos, rowData, averageOcc); + + // evaluate each settings entry against the shared row data + outputs.resize(settingsList.size()); + for (size_t i = 0; i < settingsList.size(); ++i) { + calculatedEdxFromRowData(rowData, settingsList[i], i, trackTime0, trackOrig, averageOcc, outputs[i], mcLabel); + } +} + +void CalculatedEdx::calculatedEdx(o2::tpc::TrackTPC& track, dEdxInfo& output, AverageOccupancy& averageOcc, float low, float high, CorrectionFlags correctionMask, ClusterFlags clusterMask, int subthresholdMethod, int stackBoundaryMethod, const char* debugRootFile, float maxSubthresholdChargeTot, float maxSubthresholdChargeMax, int sameRowClusterMethod) +{ + dEdxSettings settings; + settings.low = low; + settings.high = high; + settings.correctionMask = correctionMask; + settings.clusterMask = clusterMask; + settings.subthresholdMethod = subthresholdMethod; + settings.stackBoundaryMethod = stackBoundaryMethod; + settings.debugRootFile = debugRootFile; + settings.maxSubthresholdChargeTot = maxSubthresholdChargeTot; + settings.maxSubthresholdChargeMax = maxSubthresholdChargeMax; + settings.sameRowClusterMethod = sameRowClusterMethod; + + o2::tpc::TrackTPC trackOrig; + if (mDebug) { + trackOrig = track; // pristine track, before refit/propagation mutates it cluster-by-cluster below + } + const float trackTime0 = track.getTime0(); + + std::vector rowData; + gatherRowClusterData(track, rowData, averageOcc); + + calculatedEdxFromRowData(rowData, settings, 0, trackTime0, trackOrig, averageOcc, output); } float CalculatedEdx::getTruncMean(std::vector& charge, float low, float high) const @@ -499,6 +1047,7 @@ float CalculatedEdx::getTruncMean(std::vector& charge, float low, float h if (nCl > 0) { sum /= nCl; } + // if nCl == 0 (charge was empty, or too few entries for low/high to select any index), sum stays 0 return sum; } @@ -518,7 +1067,11 @@ float CalculatedEdx::getTrackTopologyCorrectionPol(const o2::tpc::TrackTPC& trac { const float snp = std::abs(track.getSnp()); const float tgl = track.getTgl(); - const float snp2 = snp * snp; + constexpr float maxSnp2 = 0.99f; + float snp2 = snp * snp; + if (snp2 > maxSnp2) { + snp2 = maxSnp2; + } const float tgl2 = tgl * tgl; const float sec2 = 1.f / (1.f - snp2); const float tanTheta = std::sqrt(tgl2 * sec2); @@ -532,7 +1085,37 @@ float CalculatedEdx::getTrackTopologyCorrectionPol(const o2::tpc::TrackTPC& trac return effectiveLength; } -void CalculatedEdx::loadCalibsFromCCDB(long runNumberOrTimeStamp, const bool isMC) +unsigned int CalculatedEdx::getOccupancy(float clTime) const +{ + // occupancy is only meaningful when the refit method is used, since mTPCRefitterOccMap is only filled by setRefit() + const int nTimeBinsPerOccupBin = 16; + const int iBinOcc = clTime / nTimeBinsPerOccupBin + 2; + if (!mRefit || iBinOcc < 0 || static_cast(iBinOcc) >= mTPCRefitterOccMap.size()) { + return -1; + } + return mTPCRefitterOccMap[iBinOcc]; +} + +bool CalculatedEdx::isInStackBoundaries(int stackNumber, unsigned char rowIndex, int stackBoundaryMethod) +{ + // retrieve boundaries for the given stack + const auto& boundaries = mStackBoundaries[stackNumber]; + // check direct match for method 1 or 2 + for (unsigned char boundary : boundaries) { + if (rowIndex == boundary) { + return true; + } + } + // additional checks for method 2 + if (stackBoundaryMethod == 2) { + if (rowIndex == boundaries[0] + 1 || rowIndex == boundaries[1] - 1) { + return true; + } + } + return false; +} + +void CalculatedEdx::loadCalibsFromCCDB(long runNumberOrTimeStamp, const bool isMC, const bool loadSCCorrMap, const bool loadSCCorrMapForRefit, const bool loadVDriftForRefit) { // setup CCDB manager auto& cm = o2::ccdb::BasicCCDBManager::instance(); @@ -562,7 +1145,7 @@ void CalculatedEdx::loadCalibsFromCCDB(long runNumberOrTimeStamp, const bool isM mCalibCont.setGainMapResidual(gainMapResidual); // set the residual dEdx correction - o2::tpc::CalibdEdxCorrection* residualObj = cm.getForTimeStamp(o2::tpc::CDBTypeMap.at(o2::tpc::CDBType::CalTimeGain), tRun); + o2::tpc::CalibdEdxCorrection* residualObj = isMC ? cm.getForTimeStamp(o2::tpc::CDBTypeMap.at(o2::tpc::CDBType::CalTimeGainMC), tRun) : cm.getForTimeStamp(o2::tpc::CDBTypeMap.at(o2::tpc::CDBType::CalTimeGain), tRun); const auto* residualCorr = static_cast(residualObj); mCalibCont.setResidualCorrection(*residualCorr); @@ -583,14 +1166,44 @@ void CalculatedEdx::loadCalibsFromCCDB(long runNumberOrTimeStamp, const bool isM const o2::base::MatLayerCylSet* matLut = o2::base::MatLayerCylSet::rectifyPtrFromFile(cm.get("GLO/Param/MatLUT")); propagator->setMatLUT(matLut); - // load sc correction maps - auto avgMap = isMC ? cm.getForTimeStamp(o2::tpc::CDBTypeMap.at(o2::tpc::CDBType::CalCorrMapMC), tRun) : cm.getForTimeStamp(o2::tpc::CDBTypeMap.at(o2::tpc::CDBType::CalCorrMap), tRun); - avgMap->rectifyAfterReadingFromFile(); + // load the space-charge correction maps + if (loadSCCorrMap || loadSCCorrMapForRefit) { + auto avgMap = isMC ? cm.getForTimeStamp(o2::tpc::CDBTypeMap.at(o2::tpc::CDBType::CalCorrMapMC), tRun) : cm.getForTimeStamp(o2::tpc::CDBTypeMap.at(o2::tpc::CDBType::CalCorrMap), tRun); + avgMap->rectifyAfterReadingFromFile(); + + if (loadSCCorrMap) { + auto derMap = isMC ? cm.getForTimeStamp(o2::tpc::CDBTypeMap.at(o2::tpc::CDBType::CalCorrDerivMapMC), tRun) : cm.getForTimeStamp(o2::tpc::CDBTypeMap.at(o2::tpc::CDBType::CalCorrDerivMap), tRun); + derMap->rectifyAfterReadingFromFile(); + mSCdEdxCorrection.setCorrectionMaps(avgMap, derMap); + } + + if (loadSCCorrMapForRefit) { + // feed the space-charge-corrected map into the refit transform + setTPCCorrMap(*avgMap); + LOGP(info, "refit transform: using CCDB space-charge correction map {}", + o2::tpc::CDBTypeMap.at(isMC ? o2::tpc::CDBType::CalCorrMapMC : o2::tpc::CDBType::CalCorrMap)); + } + } - auto derMap = isMC ? cm.getForTimeStamp(o2::tpc::CDBTypeMap.at(o2::tpc::CDBType::CalCorrDerivMapMC), tRun) : cm.getForTimeStamp(o2::tpc::CDBTypeMap.at(o2::tpc::CDBType::CalCorrDerivMap), tRun); - derMap->rectifyAfterReadingFromFile(); + // apply the calibrated drift velocity + time offset to the refit transform + if (loadVDriftForRefit) { + const bool prevFatalWhenNull = cm.getFatalWhenNull(); + cm.setFatalWhenNull(false); + if (auto* vd = cm.getForTimeStamp(o2::tpc::CDBTypeMap.at(o2::tpc::CDBType::CalVDriftTgl), tRun)) { + setTPCVDrift(*vd); + LOGP(info, "refit transform vDrift calib: corrFact={:.5f} refVDrift={:.5f} timeOffset={:.4f}us", vd->corrFact, vd->refVDrift, vd->getTimeOffset()); + } else { + LOGP(warning, "no TPC/Calib/VDriftTgl at ts {} -- refit transform stays at nominal vDrift/t0", tRun); + } + cm.setFatalWhenNull(prevFatalWhenNull); + } - mSCdEdxCorrection.setCorrectionMaps(avgMap, derMap); + // set the dead channel map + o2::tpc::DeadChannelMapCreator deadCMCreator; + deadCMCreator.init(); + deadCMCreator.load(tRun); + const o2::tpc::CalDet& deadMap = deadCMCreator.getDeadChannelMap(); + mCalibCont.setDeadChannelMap(deadMap); } void CalculatedEdx::loadCalibsFromLocalCCDBFolder(const char* localCCDBFolder) @@ -602,6 +1215,20 @@ void CalculatedEdx::loadCalibsFromLocalCCDBFolder(const char* localCCDBFolder) setZeroSuppressionThresholdFromFile(localCCDBFolder, "/TPC/Config/FEEPad/snapshot.root", "ccdb_object"); setMagneticFieldFromFile(localCCDBFolder, "/GLO/Config/GRPMagField/snapshot.root", "ccdb_object"); setPropagatorFromFile(localCCDBFolder, "/GLO/Param/MatLUT/snapshot.root", "ccdb_object"); + setVDriftFromFile(localCCDBFolder, "/TPC/Calib/VDriftTgl/snapshot.root", "ccdb_object"); // optional: skipped if absent +} + +void CalculatedEdx::setVDriftFromFile(const char* folder, const char* file, const char* object) +{ + std::unique_ptr vdFile(TFile::Open(fmt::format("{}{}", folder, file).data())); + if (!vdFile || vdFile->IsZombie()) { + LOGP(warning, "no {}{} -- refit transform stays at nominal vDrift/t0", folder, file); + return; + } + if (auto* vd = (o2::tpc::VDriftCorrFact*)vdFile->Get(object)) { + setTPCVDrift(*vd); + LOGP(info, "refit transform vDrift calib from {}: corrFact={:.5f} refVDrift={:.5f} timeOffset={:.4f}us", vdFile->GetName(), vd->corrFact, vd->refVDrift, vd->getTimeOffset()); + } } void CalculatedEdx::setTrackTopologyCorrectionFromFile(const char* folder, const char* file, const char* object) @@ -673,12 +1300,4 @@ void CalculatedEdx::setPropagatorFromFile(const char* folder, const char* file, o2::base::MatLayerCylSet* matLut = o2::base::MatLayerCylSet::rectifyPtrFromFile((o2::base::MatLayerCylSet*)matLutFile->Get(object)); propagator->setMatLUT(matLut); } -} - -unsigned int CalculatedEdx::getOccupancy(const o2::tpc::ClusterNative& cl) const -{ - const int nTimeBinsPerOccupBin = 16; - const int iBinOcc = cl.getTime() / nTimeBinsPerOccupBin + 2; - const unsigned int occupancy = mTPCRefitterOccMap.empty() ? -1 : mTPCRefitterOccMap[iBinOcc]; - return occupancy; -} +} \ No newline at end of file diff --git a/Detectors/TPC/calibration/src/TPCCalibrationLinkDef.h b/Detectors/TPC/calibration/src/TPCCalibrationLinkDef.h index 740d3f9138e57..08998a696cfeb 100644 --- a/Detectors/TPC/calibration/src/TPCCalibrationLinkDef.h +++ b/Detectors/TPC/calibration/src/TPCCalibrationLinkDef.h @@ -115,6 +115,7 @@ #pragma link C++ class o2::tpc::TPCFastSpaceChargeCorrectionHelper + ; #pragma link C++ class o2::tpc::CalculatedEdx + ; +#pragma link C++ struct o2::tpc::AverageOccupancy + ; #pragma link C++ class o2::tpc::TPCScaler + ; #pragma link C++ struct o2::tpc::TPCScalerWeights + ; #pragma link C++ class o2::tpc::TPCMShapeCorrection + ; diff --git a/GPU/GPUTracking/Interface/GPUO2InterfaceRefit.cxx b/GPU/GPUTracking/Interface/GPUO2InterfaceRefit.cxx index cd184b3820533..855ae10fc4228 100644 --- a/GPU/GPUTracking/Interface/GPUO2InterfaceRefit.cxx +++ b/GPU/GPUTracking/Interface/GPUO2InterfaceRefit.cxx @@ -133,10 +133,10 @@ void GPUO2InterfaceRefit::updateCalib(const TPCFastTransformPOD* trans, float bz mRefit->SetFastTransform(trans); } -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::RefitTrackAsGPU(o2::tpc::TrackTPC& trk, bool outward, bool resetCov, bool* reachedReferenceOut) { return mRefit->RefitTrackAsGPU(trk, outward, resetCov, reachedReferenceOut); } +int32_t GPUO2InterfaceRefit::RefitTrackAsTrackParCov(o2::tpc::TrackTPC& trk, bool outward, bool resetCov, bool* reachedReferenceOut) { return mRefit->RefitTrackAsTrackParCov(trk, outward, resetCov, reachedReferenceOut); } +int32_t GPUO2InterfaceRefit::RefitTrackAsGPU(o2::track::TrackParCov& trk, const o2::tpc::TrackTPCClusRef& clusRef, float time0, float* chi2, bool outward, bool resetCov, bool* reachedReferenceOut) { return mRefit->RefitTrackAsGPU(trk, clusRef, time0, chi2, outward, resetCov, reachedReferenceOut); } +int32_t GPUO2InterfaceRefit::RefitTrackAsTrackParCov(o2::track::TrackParCov& trk, const o2::tpc::TrackTPCClusRef& clusRef, float time0, float* chi2, bool outward, bool resetCov, bool* reachedReferenceOut) { return mRefit->RefitTrackAsTrackParCov(trk, clusRef, time0, chi2, outward, resetCov, reachedReferenceOut); } 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..351fee25e865f 100644 --- a/GPU/GPUTracking/Interface/GPUO2InterfaceRefit.h +++ b/GPU/GPUTracking/Interface/GPUO2InterfaceRefit.h @@ -64,10 +64,10 @@ class GPUO2InterfaceRefit GPUO2InterfaceRefit(const o2::tpc::ClusterNativeAccess* cl, const o2::gpu::TPCFastTransformPOD* trans, float bzNominalGPU, const o2::tpc::TPCClRefElem* trackRef, uint32_t nHbfPerTf = 0, const uint8_t* sharedmap = nullptr, const uint32_t* occupancymap = nullptr, int32_t occupancyMapSize = -1, const std::vector* trks = nullptr, o2::base::Propagator* p = nullptr); ~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); + int32_t RefitTrackAsGPU(o2::tpc::TrackTPC& trk, bool outward = false, bool resetCov = false, bool* reachedReferenceOut = nullptr); + int32_t RefitTrackAsTrackParCov(o2::tpc::TrackTPC& trk, bool outward = false, bool resetCov = false, bool* reachedReferenceOut = nullptr); + int32_t RefitTrackAsGPU(o2::track::TrackParCov& trk, const o2::tpc::TrackTPCClusRef& clusRef, float time0, float* chi2 = nullptr, bool outward = false, bool resetCov = false, bool* reachedReferenceOut = nullptr); + int32_t RefitTrackAsTrackParCov(o2::track::TrackParCov& trk, const o2::tpc::TrackTPCClusRef& clusRef, float time0, float* chi2 = nullptr, bool outward = false, bool resetCov = false, bool* reachedReferenceOut = nullptr); void setTrackReferenceX(float v); void setIgnoreErrorsAtTrackEnds(bool v); void updateCalib(const o2::gpu::TPCFastTransformPOD* trans, float bzNominalGPU); diff --git a/GPU/GPUTracking/Merger/GPUTPCGMO2Output.cxx b/GPU/GPUTracking/Merger/GPUTPCGMO2Output.cxx index d7c8eb9c44aab..1ba55d98f5ca1 100644 --- a/GPU/GPUTracking/Merger/GPUTPCGMO2Output.cxx +++ b/GPU/GPUTracking/Merger/GPUTPCGMO2Output.cxx @@ -124,6 +124,12 @@ GPUdii() void GPUTPCGMO2Output::Thread(int32_t nBlocks TrackTPC oTrack; const int32_t i = trackSort[iTmp].x; const auto& track = tracks[i]; + GPUCA_DEBUG_STREAMER_CHECK(if (o2::utils::DebugStreamer::checkStream(o2::utils::StreamFlags::streamdEdx)) { + o2::utils::DebugStreamer::instance()->getStreamer("debug_dedx", "UPDATE") << o2::utils::DebugStreamer::instance()->getUniqueTreeName("tree_indices").data() + << "trackID=" << iTmp + << "iTrk=" << i + << "\n"; + }) auto snpIn = track.GetParam().GetSinPhi(); if (snpIn > SNPThresh) { snpIn = SNPThresh; @@ -185,7 +191,7 @@ GPUdii() void GPUTPCGMO2Output::Thread(int32_t nBlocks uint32_t nOutCl2 = 0; float t1 = 0, t2 = 0; int32_t sector1 = 0, sector2 = 0; - const o2::tpc::ClusterNativeAccess* GPUrestrict() clusters = merger.GetConstantMem()->ioPtrs.clustersNative; + const o2::tpc::ClusterNativeAccess* GPUrestrict() clusters = merger.GetConstantMem() -> ioPtrs.clustersNative; for (uint32_t j = 0; j < track.NClusters(); j++) { if ((trackClusters[track.FirstClusterRef() + j].state & flagsReject)) { continue; @@ -278,7 +284,7 @@ template <> GPUdii() void GPUTPCGMO2Output::Thread(int32_t nBlocks, int32_t nThreads, int32_t iBlock, int32_t iThread, GPUsharedref() GPUSharedMemory& smem, processorType& GPUrestrict() merger) { #ifndef GPUCA_GPUCODE - const o2::tpc::ClusterNativeAccess* GPUrestrict() clusters = merger.GetConstantMem()->ioPtrs.clustersNative; + const o2::tpc::ClusterNativeAccess* GPUrestrict() clusters = merger.GetConstantMem() -> ioPtrs.clustersNative; if (clusters == nullptr || clusters->clustersMCTruth == nullptr) { return; } diff --git a/GPU/GPUTracking/Merger/GPUTPCGMTrackParam.cxx b/GPU/GPUTracking/Merger/GPUTPCGMTrackParam.cxx index 23842d8a1f859..7f0c5f12e6baf 100644 --- a/GPU/GPUTracking/Merger/GPUTPCGMTrackParam.cxx +++ b/GPU/GPUTracking/Merger/GPUTPCGMTrackParam.cxx @@ -52,7 +52,7 @@ GPUd() bool GPUTPCGMTrackParam::Fit(GPUTPCGMMerger* GPUrestrict() merger, int32_ static constexpr float kDeg2Rad = M_PI / 180.f; CADEBUG(static constexpr float kSectAngle = 2 * M_PI / 18.f); - const GPUParam& GPUrestrict() param = merger->Param(); + const GPUParam& GPUrestrict() param = merger -> Param(); GPUdEdx dEdx, dEdxAlt; GPUTPCGMPropagator prop; @@ -126,7 +126,7 @@ GPUd() bool GPUTPCGMTrackParam::Fit(GPUTPCGMMerger* GPUrestrict() merger, int32_ const float clAlpha = param.Alpha(clusters[ihit].sector); float xx, yy, zz; { - const ClusterNative& GPUrestrict() cl = merger->GetConstantMem()->ioPtrs.clustersNative->clustersLinear[clusters[ihit].num]; + const ClusterNative& GPUrestrict() cl = merger -> GetConstantMem()->ioPtrs.clustersNative->clustersLinear[clusters[ihit].num]; merger->GetConstantMem()->calibObjects.fastTransform->Transform(clusters[ihit].sector, clusters[ihit].row, cl.getPad(), cl.getTime(), xx, yy, zz, mTOffset); } // clang-format off @@ -304,11 +304,11 @@ GPUd() bool GPUTPCGMTrackParam::Fit(GPUTPCGMMerger* GPUrestrict() merger, int32_ relTime /= clusterCount; relTime = relTime - CAMath::Round(relTime); if (acc) { - dEdx.fillCluster(qtot, qmax, cluster.row, cluster.sector, mP[2], mP[3], merger->GetConstantMem()->calibObjects, zz, pad, relTime); + dEdx.fillCluster(qtot, qmax, cluster.row, cluster.sector, mP[2], mP[3], merger->GetConstantMem()->calibObjects, zz, pad, relTime, iTrk, clusterState); } if GPUCA_RTC_CONSTEXPR (GPUCA_GET_CONSTEXPR(param.rec.tpc, dEdxClusterRejectionFlagMask) != GPUCA_GET_CONSTEXPR(param.rec.tpc, dEdxClusterRejectionFlagMaskAlt)) { if (accAlt) { - dEdxAlt.fillCluster(qtot, qmax, cluster.row, cluster.sector, mP[2], mP[3], merger->GetConstantMem()->calibObjects, zz, pad, relTime); + dEdxAlt.fillCluster(qtot, qmax, cluster.row, cluster.sector, mP[2], mP[3], merger->GetConstantMem()->calibObjects, zz, pad, relTime, iTrk, clusterState); } } } @@ -364,14 +364,16 @@ GPUd() bool GPUTPCGMTrackParam::Fit(GPUTPCGMMerger* GPUrestrict() merger, int32_ return true; } -GPUdni() void GPUTPCGMTrackParam::MoveToReference(GPUTPCGMPropagator& prop, const GPUParam& param, float& Alpha) +GPUdni() bool GPUTPCGMTrackParam::MoveToReference(GPUTPCGMPropagator& prop, const GPUParam& param, float& Alpha) { static constexpr float kDeg2Rad = M_PI / 180.f; static constexpr float kSectAngle = 2 * M_PI / 18.f; + bool reachedReference = true; if (param.rec.tpc.trackReferenceX <= 500) { GPUTPCGMTrackParam save = *this; float saveAlpha = Alpha; + reachedReference = false; for (int32_t attempt = 0; attempt < 3; attempt++) { float dAngle = CAMath::Round(CAMath::ATan2(mP[0], mX) / kDeg2Rad / 20.f) * kSectAngle; Alpha += dAngle; @@ -380,7 +382,7 @@ GPUdni() void GPUTPCGMTrackParam::MoveToReference(GPUTPCGMPropagator& prop, cons } ConstrainSinPhi(); if (CAMath::Abs(mP[0]) <= mX * CAMath::Tan(kSectAngle / 2.f)) { - return; + return true; } } *this = save; @@ -392,6 +394,7 @@ GPUdni() void GPUTPCGMTrackParam::MoveToReference(GPUTPCGMPropagator& prop, cons ConstrainSinPhi(); Alpha += dAngle; } + return reachedReference; } GPUd() void GPUTPCGMTrackParam::MirrorTo(GPUTPCGMPropagator& GPUrestrict() prop, float toY, float toZ, bool inFlyDirection, const GPUParam& param, uint8_t row, uint8_t clusterState, bool mirrorParameters, int8_t sector) @@ -440,7 +443,7 @@ GPUd() int32_t GPUTPCGMTrackParam::MergeDoubleRowClusters(int32_t& ihit, int32_t xx = yy = zz = 0.f; clusterState = 0; while (true) { - const ClusterNative& GPUrestrict() cl = merger->GetConstantMem()->ioPtrs.clustersNative->clustersLinear[clusters[ihit].num]; + const ClusterNative& GPUrestrict() cl = merger -> GetConstantMem()->ioPtrs.clustersNative->clustersLinear[clusters[ihit].num]; float clamp = cl.getQtot(); float clx, cly, clz; merger->GetConstantMem()->calibObjects.fastTransform->Transform(clusters[ihit].sector, clusters[ihit].row, cl.getPad(), cl.getTime(), clx, cly, clz, mTOffset); @@ -493,7 +496,7 @@ GPUd() float GPUTPCGMTrackParam::AttachClusters(const GPUTPCGMMerger* GPUrestric if (param.rec.tpc.disableRefitAttachment & 1) { return -1e6f; } - const GPUTPCTracker& GPUrestrict() tracker = *(Merger->GetConstantMem()->tpcTrackers + sector); + const GPUTPCTracker& GPUrestrict() tracker = *(Merger -> GetConstantMem()->tpcTrackers + sector); const GPUTPCRow& GPUrestrict() row = tracker.Row(iRow); GPUglobalref() const cahit2* hits = tracker.HitData(row); GPUglobalref() const calink* firsthit = tracker.FirstHitInBin(row); @@ -679,7 +682,7 @@ GPUdi() void GPUTPCGMTrackParam::AttachClustersLooperFollow(const GPUTPCGMMerger bool inFlyDirection = (Merger->MergedTracks()[iTrack].Leg() & 1) ^ up; static constexpr float kSectAngle = 2 * M_PI / 18.f; - const GPUParam& GPUrestrict() param = Merger->Param(); + const GPUParam& GPUrestrict() param = Merger -> Param(); bool right = (mP[2] < 0) ^ up; const int32_t sectorSide = sector >= (int32_t)(GPUTPCGeometry::NSECTORS / 2) ? (GPUTPCGeometry::NSECTORS / 2) : 0; float lrFactor = right ^ !up ? 1.f : -1.f; @@ -803,7 +806,7 @@ GPUd() float GPUTPCGMTrackParam::ShiftZ(const GPUTPCGMMergedTrackHit* clusters, if (N == 0) { N = 1; } - const auto& GPUrestrict() cls = merger->GetConstantMem()->ioPtrs.clustersNative->clustersLinear; + const auto& GPUrestrict() cls = merger -> GetConstantMem()->ioPtrs.clustersNative->clustersLinear; float z0 = cls[clusters[0].num].getTime(), zn = cls[clusters[N - 1].num].getTime(); const auto tmp = zn > z0 ? std::array{zn, z0, GPUTPCGeometry::Row2X(clusters[N - 1].row)} : std::array{z0, zn, GPUTPCGeometry::Row2X(clusters[0].row)}; return ShiftZ(merger, clusters[0].sector, tmp[0], tmp[1], tmp[2]); diff --git a/GPU/GPUTracking/Merger/GPUTPCGMTrackParam.h b/GPU/GPUTracking/Merger/GPUTPCGMTrackParam.h index 51689753f1ca5..6426c1c2ad776 100644 --- a/GPU/GPUTracking/Merger/GPUTPCGMTrackParam.h +++ b/GPU/GPUTracking/Merger/GPUTPCGMTrackParam.h @@ -142,7 +142,7 @@ class GPUTPCGMTrackParam GPUd() bool CheckCov() const; GPUd() bool Fit(GPUTPCGMMerger* merger, int32_t iTrk, GPUTPCGMMergedTrackHit* clusters, int32_t& N, int32_t& NTolerated, float& Alpha, int32_t attempt, float maxSinPhi, GPUTPCGMMergedTrack& track); - GPUd() void MoveToReference(GPUTPCGMPropagator& prop, const GPUParam& param, float& alpha); + GPUd() bool MoveToReference(GPUTPCGMPropagator& prop, const GPUParam& param, float& alpha); ///< returns false if the track could not be propagated to param.rec.tpc.trackReferenceX (position/momentum then rolled back to their pre-call values, but the trailing sector-normalization step may still rotate the track and update alpha); true if it reached it or the reference is disabled (trackReferenceX > 500) GPUd() void MirrorTo(GPUTPCGMPropagator& prop, float toY, float toZ, bool inFlyDirection, const GPUParam& param, uint8_t row, uint8_t clusterState, bool mirrorParameters, int8_t sector); GPUd() int32_t MergeDoubleRowClusters(int32_t& ihit, int32_t wayDirection, GPUTPCGMMergedTrackHit* clusters, const GPUTPCGMMerger* merger, GPUTPCGMPropagator& prop, float& xx, float& yy, float& zz, int32_t maxN, float clAlpha, uint8_t& clusterState, bool rejectChi2); @@ -214,12 +214,12 @@ class GPUTPCGMTrackParam private: GPUd() int32_t initResetT0(); - float mX; // x position - float mTOffset; // Z offset with early transform, T offset otherwise - float mP[5]; // 'active' track parameters: Y, Z, SinPhi, DzDs, q/Pt - float mC[15]; // the covariance matrix for Y,Z,SinPhi,.. - float mChi2; // the chi^2 value - int32_t mNDF; // the Number of Degrees of Freedom + float mX; // x position + float mTOffset; // Z offset with early transform, T offset otherwise + float mP[5]; // 'active' track parameters: Y, Z, SinPhi, DzDs, q/Pt + float mC[15]; // the covariance matrix for Y,Z,SinPhi,.. + float mChi2; // the chi^2 value + int32_t mNDF; // the Number of Degrees of Freedom }; struct GPUTPCGMLoopData { diff --git a/GPU/GPUTracking/Refit/GPUTrackingRefit.cxx b/GPU/GPUTracking/Refit/GPUTrackingRefit.cxx index de0525edcce2a..e10bc97b917dc 100644 --- a/GPU/GPUTracking/Refit/GPUTrackingRefit.cxx +++ b/GPU/GPUTracking/Refit/GPUTrackingRefit.cxx @@ -212,7 +212,7 @@ GPUd() static const float* getPar(const GPUTPCGMTrackParam& trk) { return trk.Ge GPUd() static const float* getPar(const TrackParCov& trk) { return trk.getParams(); } template -GPUd() int32_t GPUTrackingRefit::RefitTrack(T& trkX, bool outward, bool resetCov) +GPUd() int32_t GPUTrackingRefit::RefitTrack(T& trkX, bool outward, bool resetCov, bool* reachedReferenceOut) { CADEBUG(int32_t ii; printf("\nRefitting track\n")); typename internal::refitTrackTypes::propagator prop; @@ -392,20 +392,22 @@ GPUd() int32_t GPUTrackingRefit::RefitTrack(T& trkX, bool outward, bool resetCov resetCov = false; nFitted++; } + bool reachedReference = true; if constexpr (std::is_same_v) { float alpha = prop.GetAlpha(); - trk.MoveToReference(prop, *mPparam, alpha); + reachedReference = trk.MoveToReference(prop, *mPparam, alpha); trk.NormalizeAlpha(alpha); prop.SetAlpha(alpha); } else if constexpr (std::is_same_v) { static constexpr float kDeg2Rad = M_PI / 180.f; static constexpr float kSectAngle = 2 * M_PI / 18.f; if (mPparam->rec.tpc.trackReferenceX <= 500) { - if (prop->PropagateToXBxByBz(trk, mPparam->rec.tpc.trackReferenceX)) { + reachedReference = prop->PropagateToXBxByBz(trk, mPparam->rec.tpc.trackReferenceX); + if (reachedReference) { 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); + reachedReference = trk.rotate(newAlpha) && prop->PropagateToXBxByBz(trk, mPparam->rec.tpc.trackReferenceX); } } } @@ -414,16 +416,19 @@ GPUd() int32_t GPUTrackingRefit::RefitTrack(T& trkX, bool outward, bool resetCov } convertTrack::propagator>(trkX, trk, prop, &TrackParCovChi2); + if (reachedReferenceOut) { + *reachedReferenceOut = reachedReference; + } return nFitted; } #if !defined(GPUCA_GPUCODE) || defined(GPUCA_GPUCODE_DEVICE) // FIXME: DR: WORKAROUND to avoid CUDA bug creating host symbols for device code. -template GPUdni() int32_t GPUTrackingRefit::RefitTrack(GPUTPCGMMergedTrack& trk, bool outward, bool resetCov); -template GPUdni() int32_t GPUTrackingRefit::RefitTrack(GPUTPCGMMergedTrack& trk, bool outward, bool resetCov); -template GPUdni() int32_t GPUTrackingRefit::RefitTrack(TrackTPC& trk, bool outward, bool resetCov); -template GPUdni() int32_t GPUTrackingRefit::RefitTrack(TrackTPC& trk, bool outward, bool resetCov); -template GPUdni() int32_t GPUTrackingRefit::RefitTrack(GPUTrackingRefit::TrackParCovWithArgs& trk, bool outward, bool resetCov); -template GPUdni() int32_t GPUTrackingRefit::RefitTrack(GPUTrackingRefit::TrackParCovWithArgs& trk, bool outward, bool resetCov); +template GPUdni() int32_t GPUTrackingRefit::RefitTrack(GPUTPCGMMergedTrack& trk, bool outward, bool resetCov, bool* reachedReferenceOut); +template GPUdni() int32_t GPUTrackingRefit::RefitTrack(GPUTPCGMMergedTrack& trk, bool outward, bool resetCov, bool* reachedReferenceOut); +template GPUdni() int32_t GPUTrackingRefit::RefitTrack(TrackTPC& trk, bool outward, bool resetCov, bool* reachedReferenceOut); +template GPUdni() int32_t GPUTrackingRefit::RefitTrack(TrackTPC& trk, bool outward, bool resetCov, bool* reachedReferenceOut); +template GPUdni() int32_t GPUTrackingRefit::RefitTrack(GPUTrackingRefit::TrackParCovWithArgs& trk, bool outward, bool resetCov, bool* reachedReferenceOut); +template GPUdni() int32_t GPUTrackingRefit::RefitTrack(GPUTrackingRefit::TrackParCovWithArgs& trk, bool outward, bool resetCov, bool* reachedReferenceOut); #endif #ifndef GPUCA_GPUCODE diff --git a/GPU/GPUTracking/Refit/GPUTrackingRefit.h b/GPU/GPUTracking/Refit/GPUTrackingRefit.h index 70c9fd47d90f6..e2538756fee13 100644 --- a/GPU/GPUTracking/Refit/GPUTrackingRefit.h +++ b/GPU/GPUTracking/Refit/GPUTrackingRefit.h @@ -63,10 +63,10 @@ class GPUTrackingRefit void SetTrackHitReferences(const uint32_t* v) { mPtrackHitReferences = v; } void SetFastTransform(const TPCFastTransformPOD* v) { mPfastTransform = v; } void SetGPUParam(const GPUParam* v) { mPparam = v; } - GPUd() int32_t RefitTrackAsGPU(GPUTPCGMMergedTrack& trk, bool outward = false, bool resetCov = false) { return RefitTrack(trk, outward, resetCov); } - GPUd() int32_t RefitTrackAsTrackParCov(GPUTPCGMMergedTrack& trk, bool outward = false, bool resetCov = false) { return RefitTrack(trk, outward, resetCov); } - GPUd() int32_t RefitTrackAsGPU(o2::tpc::TrackTPC& trk, bool outward = false, bool resetCov = false) { return RefitTrack(trk, outward, resetCov); } - GPUd() int32_t RefitTrackAsTrackParCov(o2::tpc::TrackTPC& trk, bool outward = false, bool resetCov = false) { return RefitTrack(trk, outward, resetCov); } + GPUd() int32_t RefitTrackAsGPU(GPUTPCGMMergedTrack& trk, bool outward = false, bool resetCov = false, bool* reachedReferenceOut = nullptr) { return RefitTrack(trk, outward, resetCov, reachedReferenceOut); } + GPUd() int32_t RefitTrackAsTrackParCov(GPUTPCGMMergedTrack& trk, bool outward = false, bool resetCov = false, bool* reachedReferenceOut = nullptr) { return RefitTrack(trk, outward, resetCov, reachedReferenceOut); } + GPUd() int32_t RefitTrackAsGPU(o2::tpc::TrackTPC& trk, bool outward = false, bool resetCov = false, bool* reachedReferenceOut = nullptr) { return RefitTrack(trk, outward, resetCov, reachedReferenceOut); } + GPUd() int32_t RefitTrackAsTrackParCov(o2::tpc::TrackTPC& trk, bool outward = false, bool resetCov = false, bool* reachedReferenceOut = nullptr) { return RefitTrack(trk, outward, resetCov, reachedReferenceOut); } struct TrackParCovWithArgs { o2::track::TrackParCov& trk; @@ -74,15 +74,15 @@ class GPUTrackingRefit float time0; float* chi2; }; - GPUd() int32_t RefitTrackAsGPU(o2::track::TrackParCov& trk, const o2::tpc::TrackTPCClusRef& clusRef, float time0, float* chi2 = nullptr, bool outward = false, bool resetCov = false) + GPUd() int32_t RefitTrackAsGPU(o2::track::TrackParCov& trk, const o2::tpc::TrackTPCClusRef& clusRef, float time0, float* chi2 = nullptr, bool outward = false, bool resetCov = false, bool* reachedReferenceOut = nullptr) { TrackParCovWithArgs x{trk, clusRef, time0, chi2}; - return RefitTrack(x, outward, resetCov); + return RefitTrack(x, outward, resetCov, reachedReferenceOut); } - 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, bool* reachedReferenceOut = nullptr) { TrackParCovWithArgs x{trk, clusRef, time0, chi2}; - return RefitTrack(x, outward, resetCov); + return RefitTrack(x, outward, resetCov, reachedReferenceOut); } bool mIgnoreErrorsOnTrackEnds = true; // Ignore errors during propagation / update at the beginning / end of tracks for int16_t tracks / tracks with high incl. angle @@ -97,7 +97,7 @@ class GPUTrackingRefit const TPCFastTransformPOD* mPfastTransform = nullptr; // Ptr to TPC fast transform object helper const GPUParam* mPparam = nullptr; // Ptr to GPUParam template - GPUd() int32_t RefitTrack(T& trk, bool outward, bool resetCov); + GPUd() int32_t RefitTrack(T& trk, bool outward, bool resetCov, bool* reachedReferenceOut = nullptr); template GPUd() void convertTrack(T& trk, const S& trkX, U& prop, float* chi2); template diff --git a/GPU/GPUTracking/dEdx/GPUdEdx.h b/GPU/GPUTracking/dEdx/GPUdEdx.h index dad62c1decb53..5b0363addb3f4 100644 --- a/GPU/GPUTracking/dEdx/GPUdEdx.h +++ b/GPU/GPUTracking/dEdx/GPUdEdx.h @@ -33,7 +33,7 @@ class GPUdEdx public: // The driver must call clear(), fill clusters row by row outside-in, then run computedEdx() to get the result GPUd() void clear(); - GPUd() void fillCluster(float qtot, float qmax, int32_t padRow, uint8_t sector, float trackSnp, float trackTgl, const GPUCalibObjectsConst& calib, float z, float pad, float relTime); + GPUd() void fillCluster(float qtot, float qmax, int32_t padRow, uint8_t sector, float trackSnp, float trackTgl, const GPUCalibObjectsConst& calib, float z, float pad, float relTime, int32_t iTrk, uint8_t flags); GPUd() void fillSubThreshold(int32_t padRow); GPUd() void computedEdx(GPUdEdxInfo& output, const GPUParam& param); @@ -73,7 +73,7 @@ GPUdi() void GPUdEdx::checkSubThresh(int32_t roc) mLastROC = roc; } -GPUdnii() void GPUdEdx::fillCluster(float qtot, float qmax, int32_t padRow, uint8_t sector, float trackSnp, float trackTgl, const GPUCalibObjectsConst& calib, float z, float pad, float relTime) +GPUdnii() void GPUdEdx::fillCluster(float qtot, float qmax, int32_t padRow, uint8_t sector, float trackSnp, float trackTgl, const GPUCalibObjectsConst& calib, float z, float pad, float relTime, int32_t iTrk, uint8_t flags) { // container containing all the dE/dx corrections auto calibContainer = calib.dEdxCalibContainer; @@ -84,6 +84,8 @@ GPUdnii() void GPUdEdx::fillCluster(float qtot, float qmax, int32_t padRow, uint if (mCount >= MAX_NCL) { return; } + const float clqTot = qtot; + const float clqMax = qmax; float snp2 = trackSnp * trackSnp; if (snp2 > constants::MAX_SIN_PHI_LOW) { snp2 = constants::MAX_SIN_PHI_LOW; @@ -166,6 +168,10 @@ GPUdnii() void GPUdEdx::fillCluster(float qtot, float qmax, int32_t padRow, uint << "qTotResidualCorr=" << qTotResidualCorr << "residualGainMapGain=" << residualGainMapGain << "fullGainMapGain=" << fullGainMapGain + << "iTrk=" << iTrk + << "flags=" << flags + << "clqTot=" << clqTot + << "clqMax=" << clqMax << "\n"; }) }