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..811684ce48094 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 of tracks where setRefit()'s RefitTrackAsGPU() failed and the row's track state instead came from propagating the original track 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 refit of this track failed and the row's track state came from propagating the original track + 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 refitFailed setRefit() mode only: true if the single per-track refit in refitTrack() failed, in which case the rows are propagated from the track's original 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 refitFailed); + + /// setRefit() mode only: refit the whole track once, outward, so that it ends at the outermost cluster; the track is then propagated row by row in gatherRowClusterDataForRow() + /// \param track input track, mutated in place by the refit; left unchanged if the refit fails + /// \return true if the refit failed + bool refitTrack(o2::tpc::TrackTPC& track); + + /// 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 of tracks where setRefit()'s RefitTrackAsGPU() failed but propagating the original track 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..0f70a7ef68a42 --- /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(); // copy: resetDebugCounters() below clears the member + 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..46ef2b010a1a7 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,702 @@ 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; + } + + // 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(); + + const bool refitFailed = refitTrack(track); + + // 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, refitFailed); + rowIndexOld = rowIndex; + sectorIndexOld = sectorIndex; + } - bool check = true; - if (!mPropagateTrack) { - if (mRefit == nullptr) { - LOGP(error, "mRefit is a nullptr, call the function setRefit() before looping over the tracks."); + // 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; } - mRefit->setTrackReferenceX(xPosition); - check = (mRefit->RefitTrackAsGPU(track, false, true) < 0) ? false : true; - } 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); } + if (nValidOcc > 0) { + *averageOccROC[roc] = static_cast(sumOcc) / nValidOcc; + } + } +} + +bool CalculatedEdx::refitTrack(o2::tpc::TrackTPC& track) +{ + if (!mRefit) { + return false; + } + // RefitTrackAsGPU() refits the whole track from its clusters (not just up to the current row) and may write into the track even if it fails, so keep the original state to fall back from + const o2::track::TrackParCov trackBeforeRefit = track; + // outward refit: the track ends at the outermost cluster, i.e. at the start of the cluster reference list that is traversed row by row + const bool failed = (mRefit->RefitTrackAsGPU(track, true, true) < 0) || std::isnan(track.getParam(1)); + if (failed) { + static_cast(track) = trackBeforeRefit; + } + return failed; +} - if (!check || std::isnan(track.getParam(1))) { - excludeCl += 0b1000; // 8 for failure of track propagation or refit +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; +} - 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); +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 refitFailed) +{ + 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; + if (mRefit || mPropagateTrack) { + // mRefit: the track was refit once in gatherRowClusterData() and is only propagated from row to row here, if that refit failed, the propagation starts from the track's original state + // 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; + } + } + + row.propagationFailed = (!check || std::isnan(track.getParam(1))); + row.refitFellBack = mRefit && refitFailed && !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; } - // to avoid counting the skipped cluster as a subthreshold cluster - rowIndexOld = rowIndex; - sectorIndexOld = sectorIndex; - continue; } + } + + 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(); + + const bool refitFailed = refitTrack(track); - // get charge values - float chargeTot = cl.getQtot(); - float chargeMax = cl.qMax; + // handle same (sector, row) clusters + std::vector> rowOrder; + std::map, std::vector> clustersByRow; + std::set> mergeableRows; - // get threshold - const float threshold = mCalibCont.getZeroSupressionThreshold(sectorIndex, rowIndex, pad); + 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, refitFailed); + 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; + } + } + if (nValidOcc > 0) { + *averageOccROC[roc] = static_cast(sumOcc) / nValidOcc; + } + } +} + +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}; + + 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 = 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; + }; + } - // 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; + // 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 +872,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) { - 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"; + 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) { + 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 +1039,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 +1059,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 +1077,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 +1137,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 +1158,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(); - 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(); + 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); + } - 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)); + } + } + + // 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); + } + + // 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 +1207,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 +1292,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/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..faf44492d4b21 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 GPUCA_DEBUG_STREAMER_CHECK(, 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 GPUCA_DEBUG_STREAMER_CHECK(, iTrk, clusterState)); } } } @@ -440,7 +440,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 +493,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 +679,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 +803,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..c9d8f2017af52 100644 --- a/GPU/GPUTracking/Merger/GPUTPCGMTrackParam.h +++ b/GPU/GPUTracking/Merger/GPUTPCGMTrackParam.h @@ -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/dEdx/GPUdEdx.h b/GPU/GPUTracking/dEdx/GPUdEdx.h index dad62c1decb53..a53ad186ed95f 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 GPUCA_DEBUG_STREAMER_CHECK(, 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 GPUCA_DEBUG_STREAMER_CHECK(, int32_t iTrk, uint8_t flags)) { // container containing all the dE/dx corrections auto calibContainer = calib.dEdxCalibContainer; @@ -84,6 +84,7 @@ GPUdnii() void GPUdEdx::fillCluster(float qtot, float qmax, int32_t padRow, uint if (mCount >= MAX_NCL) { return; } + GPUCA_DEBUG_STREAMER_CHECK(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 +167,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"; }) }