From 66f8031a9b29fa363a43f3f9e20a25579d6dd729 Mon Sep 17 00:00:00 2001 From: tubagundem Date: Thu, 30 Jul 2026 11:52:17 +0200 Subject: [PATCH 1/2] TPC: Improve standalone dEdx calculation class --- Detectors/TPC/calibration/CMakeLists.txt | 4 + .../include/TPCCalibration/CalculatedEdx.h | 100 ++- .../TPC/calibration/macro/calculatedEdx.C | 362 +++++++++++ .../TPC/calibration/src/CalculatedEdx.cxx | 568 ++++++++++-------- .../calibration/src/TPCCalibrationLinkDef.h | 1 + 5 files changed, 766 insertions(+), 269 deletions(-) create mode 100644 Detectors/TPC/calibration/macro/calculatedEdx.C 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..c25850c406abf 100644 --- a/Detectors/TPC/calibration/include/TPCCalibration/CalculatedEdx.h +++ b/Detectors/TPC/calibration/include/TPCCalibration/CalculatedEdx.h @@ -20,16 +20,32 @@ // o2 includes #include "DataFormatsTPC/TrackTPC.h" #include "DataFormatsTPC/dEdxInfo.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 +#include +#include +#include +#include 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 +62,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 +77,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)); } @@ -95,6 +114,9 @@ class CalculatedEdx /// \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; } @@ -110,12 +132,22 @@ class CalculatedEdx /// \param minChargeMaxThreshold upper limit for the possible minimum charge max in subthreshold treatment void setMinChargeMaxThreshold(float minChargeMaxThreshold) { mMinChargeMaxThreshold = minChargeMaxThreshold; } - /// 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; } @@ -128,18 +160,21 @@ class CalculatedEdx /// \return returns the upper limit for the possible minimum charge max in subthreshold treatment float getMinChargeMaxThreshold() { return mMinChargeMaxThreshold; } - /// fill missing clusters with minimum charge (method=0) or minimum charge/2 (method=1) or Landau (method=2) + /// fill missing clusters with minimum charge (method=0) or minimum charge/2 (method=1) void fillMissingClusters(int missingClusters[4], float minChargeTot, float minChargeMax, int method, std::array, 5>& chargeTotROC, std::array, 5>& chargeMaxROC); + void handleSameRowClusters(o2::tpc::TrackTPC& track, std::map, std::vector>& clustersByRow, std::map, o2::tpc::ClusterNative>& combinedClustersByRow, std::map>& clusterReferencesByIndex); + /// 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"); + 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"); /// get the truncated mean for the input charge vector and the truncation range low*nCl* mTracks{nullptr}; ///< vector containing the tpc tracks which will be processed std::vector* mTPCTrackClIdxVecInput{nullptr}; ///< input vector with TPC tracks cluster indicies @@ -233,18 +269,22 @@ class CalculatedEdx 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 + 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 (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 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..2edbea05de5d4 --- /dev/null +++ b/Detectors/TPC/calibration/macro/calculatedEdx.C @@ -0,0 +1,362 @@ +// 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. + +#if !defined(__CLING__) || defined(__ROOTCLING__) +#include +#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 "SimulationDataFormat/MCTrack.h" +#include "Steer/MCKinematicsReader.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)}; +} +} // 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 correctionMask which corrections to apply +/// \param clusterMask which clusters to exclude from the dE/dx calculation +/// \param lowTruncation lower cluster cut of the truncated mean, e.g. 0.015 keeps clusters starting from the 1.5% percentile of sorted charges +/// \param highTruncation higher cluster cut of the truncated mean, e.g. 0.6 keeps clusters up to the 60% percentile of sorted charges +/// \param subthresholdMethod how to fill subthreshold clusters: 0 = minimum charge, 1 = minimum charge / 2 +/// \param stackBoundaryMethod exclude clusters on GEM stack boundary rows from the dE/dx calculation: 0 = disabled, 1 = exclude the boundary row, 2 = also exclude the row next to the boundary +/// \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_t.root file with per-cluster information per worker thread +/// \param nThreads number of worker threads used to process the tracks of one event in parallel +/// \param isMC set to true for MC productions to read the true MC track (TPCTracksMCTruth branch + mcKineFileName) for each track into an additional "mctrack" 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 mcKineFileName base name (without extension) of the MC kinematics file used to resolve MC track labels; only used if isMC is true +void calculatedEdx(const std::string dir = ".", + const long runNumberOrTimeStamp = 0, + const std::string outFile = "dEdxCalc.root", + const std::string localCCDBFolder = "", + const CorrectionFlags correctionMask = CorrectionFlags::TopologyPol | CorrectionFlags::dEdxResidual, + const ClusterFlags clusterMask = ClusterFlags::ExcludeEdgeCl, + const float lowTruncation = 0.015f, + const float highTruncation = 0.6f, + const int subthresholdMethod = 0, + const int stackBoundaryMethod = 0, + 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 size_t nThreads = 8, + 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 std::string mcKineFileName = "o2sim") +{ + ROOT::EnableThreadSafety(); + + if (nThreads == 0) { + LOGP(error, "nThreads must be at least 1"); + return; + } + if (lowTruncation < 0.f || highTruncation > 1.f || lowTruncation >= highTruncation) { + LOGP(error, "Invalid truncation range [{}, {}); expected 0 <= lowTruncation < highTruncation <= 1", lowTruncation, highTruncation); + return; + } + if (subthresholdMethod != 0 && subthresholdMethod != 1) { + LOGP(error, "Invalid subthresholdMethod {}; expected 0 (minimum charge) or 1 (minimum charge / 2)", subthresholdMethod); + return; + } + if (stackBoundaryMethod < 0 || stackBoundaryMethod > 2) { + LOGP(error, "Invalid stackBoundaryMethod {}; expected 0 (disabled), 1 (exclude boundary row) or 2 (also exclude the adjacent row)", stackBoundaryMethod); + return; + } + + const std::clock_t c_start = std::clock(); + const auto t_start = std::chrono::high_resolution_clock::now(); + + // one CalculatedEdx per worker thread + std::vector calcdEdxPerThread(nThreads); + for (auto& calcdEdx : calcdEdxPerThread) { + calcdEdx.setDebug(debug); + calcdEdx.setPropagateTrack(propagateTrack); + calcdEdx.setPropagateParams(propagateParams); + } + + 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; + o2::steer::MCKinematicsReader mcReader; + if (isMC) { + if (!tpcTree->GetBranch("TPCTracksMCTruth")) { + LOGP(error, "Branch 'TPCTracksMCTruth' not found in {}/{}, cannot resolve MC truth", dir, tpcTracksFileName); + return; + } + tpcTree->SetBranchAddress("TPCTracksMCTruth", &tpcMCTruthPtr); + const auto kineFileName = fmt::format("{}/{}", dir, mcKineFileName); + if (!mcReader.initFromKinematics(kineFileName)) { + LOGP(error, "Could not initialize MCKinematicsReader from {}", kineFileName); + return; + } + } + + 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(); + } + + 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, per thread instance + for (auto& calcdEdx : calcdEdxPerThread) { + calcdEdx.setMembers(tpcTrackClIdxVecInput, clusterIndex, &tpcTracks); + if (localCCDBFolder.empty()) { + calcdEdx.loadCalibsFromCCDB(timeStamp, isMC); + } else { + calcdEdx.loadCalibsFromLocalCCDBFolder(localCCDBFolder.data()); + } + if (useRefit) { + calcdEdx.setRefit(); + } + } + + const size_t nSelectable = isMatchedToITS ? matchTracks.size() : tpcTracks.size(); + LOGP(info, "Processing event {} with {} {} using {} threads", iEvent, nSelectable, isMatchedToITS ? "matched tracks" : "tracks", nThreads); + + std::vector> tpcOut(nThreads); + std::vector> itsTracksOut(nThreads); + std::vector> matchTracksOut(nThreads); + std::vector> dEdxOut(nThreads); + std::vector> averageOccOut(nThreads); + std::vector> mcLabelOut(nThreads); + + const size_t chunkSize = (nSelectable + nThreads - 1) / nThreads; + std::vector threads; + for (size_t iThread = 0; iThread < nThreads; iThread++) { + const size_t start = iThread * chunkSize; + const size_t end = std::min(start + chunkSize, nSelectable); + if (start >= end) { + continue; + } + threads.emplace_back([&, iThread, start, end]() { + auto& calcdEdx = calcdEdxPerThread[iThread]; + const std::string debugFile = fmt::format("dEdxDebug_t{}.root", iThread); + for (size_t i = start; i < end; 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[iThread].emplace_back(itsTracks[itstpc.getRefITS().getIndex()]); + matchTracksOut[iThread].emplace_back(itstpc); + } + + TrackTPC track(tpcTracks[tpcIndex]); // local copy: setRefit()/propagation mutate the track in place + dEdxInfo dEdx; + AverageOccupancy averageOcc; + calcdEdx.calculatedEdx(track, dEdx, averageOcc, lowTruncation, highTruncation, correctionMask, clusterMask, subthresholdMethod, stackBoundaryMethod, debugFile.c_str()); + + tpcOut[iThread].emplace_back(track); + dEdxOut[iThread].emplace_back(dEdx); + averageOccOut[iThread].emplace_back(averageOcc); + if (isMC) { + mcLabelOut[iThread].emplace_back(tpcMCTruth[tpcIndex]); + } + } + }); + } + for (auto& th : threads) { + th.join(); + } + + // write out sequentially in the main thread: no locking needed since all worker threads have already joined + for (size_t iThread = 0; iThread < nThreads; iThread++) { + for (size_t i = 0; i < dEdxOut[iThread].size(); i++) { + auto& row = stream << "tree" + << "iEvent=" << iEvent + << "timeStamp=" << timeStamp + << "tpc=" << tpcOut[iThread][i] + << "dEdx=" << dEdxOut[iThread][i] + << "averageOcc=" << averageOccOut[iThread][i]; + if (tfIDTree) { + row << "tfIDInfo=" << tfIDInfo; + } + if (isMatchedToITS) { + row << "its=" << itsTracksOut[iThread][i] + << "itstpc=" << matchTracksOut[iThread][i]; + } + if (isMC) { + const auto& label = mcLabelOut[iThread][i]; + const auto* mcTrack = label.isValid() ? mcReader.getTrack(label) : nullptr; + row << "mcLabel=" << label; + if (mcTrack) { + row << "mctrack=" << *mcTrack; + } + } + 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..79921f2d9000b 100644 --- a/Detectors/TPC/calibration/src/CalculatedEdx.cxx +++ b/Detectors/TPC/calibration/src/CalculatedEdx.cxx @@ -47,7 +47,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); @@ -62,10 +62,10 @@ void CalculatedEdx::fillMissingClusters(int missingClusters[4], float minChargeT return; } + float chargeMax = (method == 1) ? minChargeMax / 2.f : minChargeMax; + float chargeTot = (method == 1) ? minChargeTot / 2.f : minChargeTot; for (int roc = 0; roc < 4; 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,14 +76,81 @@ 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::map, std::vector>& clustersByRow, std::map, o2::tpc::ClusterNative>& combinedClustersByRow, std::map>& clusterReferencesByIndex) { // get number of clusters const int nClusters = track.getNClusterReferences(); + // 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; + + track.getClusterReference(*mTPCTrackClIdxVecInput, iCl, sectorIndex, rowIndex, clusterIndexNumb); + + // add the cluster index to the corresponding (sector, row) key in clustersByRow + clustersByRow[{sectorIndex, rowIndex}].emplace_back(iCl); + + // store the reference data in clusterReferencesByIndex + clusterReferencesByIndex[iCl] = std::make_tuple(sectorIndex, rowIndex, clusterIndexNumb); + } + + // combine clusters in the same (sector, row) and store the result + for (const auto& [rowKey, clusterIndices] : clustersByRow) { + if (clusterIndices.size() > 1) { // only combine if there are multiple clusters in the same row + + // initialize variables for the combined cluster properties + float weightedPadSum = 0; + float weightedTimeSum = 0; + float totalCharge = 0; + uint16_t maxCharge = 0; + + // use the first cluster as a template for other fields + const o2::tpc::ClusterNative& firstCluster = track.getCluster(*mTPCTrackClIdxVecInput, clusterIndices[0], *mClusterIndex); + o2::tpc::ClusterNative combinedCluster = firstCluster; + + // iterate over all the clusters in the current row to combine their properties + for (int clusterIdx : clusterIndices) { + const o2::tpc::ClusterNative& cl = track.getCluster(*mTPCTrackClIdxVecInput, clusterIdx, *mClusterIndex); + + float clPad = cl.getPad(); + float clTime = cl.getTime(); + uint16_t clqTot = cl.getQtot(); + uint16_t clqMax = cl.qMax; + + // calculate weighted sums for pad and time + weightedPadSum += clPad * clqTot; + weightedTimeSum += clTime * clqTot; + totalCharge += clqTot; + maxCharge = std::max(maxCharge, clqMax); + } + + // finalize the combined cluster properties + 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); + + // store the combined cluster in the result map for the (sector, row) + combinedClustersByRow[rowKey] = combinedCluster; + } + } +} + +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) +{ + // NHits and NHitsSubthreshold values per region int nClsROC[4] = {0, 0, 0, 0}; int nClsSubThreshROC[4] = {0, 0, 0, 0}; + // corrected qTot and qMax values per region const int nType = 5; std::array, nType> chargeTotROC; std::array, nType> chargeMaxROC; @@ -92,72 +159,65 @@ void CalculatedEdx::calculatedEdx(o2::tpc::TrackTPC& track, dEdxInfo& output, fl 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; + // occupancy vector for a track (all clusters, for debugging) and per-region (for the average occupancy output) 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); - } + std::array, 4> occupancyROC; - // for missing clusters - unsigned char rowIndexOld = 0; - unsigned char sectorIndexOld = 0; + // for tracking missing clusters + unsigned char rowIndexOld = 255; + unsigned char sectorIndexOld = 255; float minChargeTot = 100000.f; float minChargeMax = 100000.f; - // loop over the clusters - for (int iCl = 0; iCl < nClusters; iCl++) { - - const o2::tpc::ClusterNative& cl = track.getCluster(*mTPCTrackClIdxVecInput, iCl, *mClusterIndex); + // 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; + + // handle same (sector, row) clusters + std::map, std::vector> clustersByRow; + std::map, o2::tpc::ClusterNative> combinedClustersByRow; + std::map> clusterReferencesByIndex; + + handleSameRowClusters(track, clustersByRow, combinedClustersByRow, clusterReferencesByIndex); + + o2::utils::TreeStreamRedirector* debugStreamer = nullptr; + o2::tpc::TrackTPC trackOrig; + if (mDebug) { + setStreamer(debugRootFile); + debugStreamer = mStreamers.at(debugRootFile).get(); + ++mDebugTrackIndex; + trackOrig = track; // pristine track, before refit/propagation mutates it cluster-by-cluster below + } - unsigned char sectorIndex = 0; - unsigned char rowIndex = 0; - unsigned int clusterIndexNumb = 0; + // loop over the clusters + for (const auto& [rowKey, clusterIndices] : clustersByRow) { + const unsigned char rowIndex = rowKey.second; + int clusterIdx = clusterIndices[0]; + const o2::tpc::ClusterNative& clConst = track.getCluster(*mTPCTrackClIdxVecInput, clusterIdx, *mClusterIndex); + const auto& [sectorIndex, rowIndexRef, clusterIndexNumb] = clusterReferencesByIndex[clusterIdx]; + bool isCombined = false; + + o2::tpc::ClusterNative cl = clConst; + + if (clusterIndices.size() > 1) { + cl = combinedClustersByRow[rowKey]; + isCombined = true; + } - // set sectorIndex, rowIndex, clusterIndexNumb - track.getClusterReference(*mTPCTrackClIdxVecInput, iCl, sectorIndex, rowIndex, clusterIndexNumb); + // get cluster values + float chargeTot = cl.getQtot(); + float chargeMax = cl.getQmax(); + const float clPad = cl.getPad(); + const float clTime = cl.getTime(); + const uint8_t flagsCl = cl.getFlags(); + unsigned int occupancy = getOccupancy(clTime); + occupancyVector.emplace_back(occupancy); // check if the cluster is shared const unsigned int absoluteIndex = mClusterIndex->clusterOffset[sectorIndex][rowIndex] + clusterIndexNumb; @@ -165,134 +225,77 @@ void CalculatedEdx::calculatedEdx(o2::tpc::TrackTPC& track, dEdxInfo& output, fl // 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 unsigned char pad = std::clamp(static_cast(clPad + 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); - // 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(); + if (stack == GEMstack::IROCgem) { + occupancyROC[0].emplace_back(occupancy); + } else if (stack == GEMstack::OROC1gem) { + occupancyROC[1].emplace_back(occupancy); + } else if (stack == GEMstack::OROC2gem) { + occupancyROC[2].emplace_back(occupancy); + } else if (stack == GEMstack::OROC3gem) { + occupancyROC[3].emplace_back(occupancy); + } 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 } + if (((clusterMask & ClusterFlags::ExcludeSplitPadCl) == ClusterFlags::ExcludeSplitPadCl) && ((flagsCl & ClusterNative::flagSplitPad) == ClusterNative::flagSplitPad)) { + excludeCl += 0b010; // 2 for split pad cluster + } + if (((clusterMask & ClusterFlags::ExcludeSplitTimeCl) == ClusterFlags::ExcludeSplitTimeCl) && ((flagsCl & ClusterNative::flagSplitTime) == ClusterNative::flagSplitTime)) { + excludeCl += 0b0100; // 4 for split time cluster + } if (((clusterMask & ClusterFlags::ExcludeSplitCl) == ClusterFlags::ExcludeSplitCl) && (((flagsCl & ClusterNative::flagSplitPad) == ClusterNative::flagSplitPad) || ((flagsCl & ClusterNative::flagSplitTime) == ClusterNative::flagSplitTime))) { - excludeCl += 0b010; // 2 for split cluster + excludeCl += 0b01000; // 8 for split cluster } if (((clusterMask & ClusterFlags::ExcludeEdgeCl) == ClusterFlags::ExcludeEdgeCl) && ((flagsCl & ClusterNative::flagEdge) == ClusterNative::flagEdge)) { - excludeCl += 0b100; // 4 for edge cluster + excludeCl += 0b010000; // 16 for edge cluster } if (((clusterMask & ClusterFlags::ExcludeSharedCl) == ClusterFlags::ExcludeSharedCl) && isShared) { - excludeCl += 0b10000; // for shared cluster + excludeCl += 0b0100000; // 32 for shared cluster + } + if (((clusterMask & ClusterFlags::ExcludeSamePadRowCl) == ClusterFlags::ExcludeSamePadRowCl) && isCombined) { + excludeCl += 0b01000000; // 64 for combined cluster + } + if ((stackBoundaryMethod == 1 || stackBoundaryMethod == 2) && isInStackBoundaries(stackNumber, rowIndex, stackBoundaryMethod)) { + excludeCl += 0b010000000; // 128 for stack boundary cluster + } + if (mCalibCont.isDead(static_cast(sectorIndex), static_cast(rowIndex), static_cast(pad))) { + excludeCl += 0b0100000000; // 256 for dead region } // get the x position of the track const float xPosition = Mapper::instance().getPadCentre(PadPos(rowIndex, 0)).X(); - bool check = true; - if (!mPropagateTrack) { - if (mRefit == nullptr) { - LOGP(error, "mRefit is a nullptr, call the function setRefit() before looping over the tracks."); - } + if (mRefit) { + // refit this track mRefit->setTrackReferenceX(xPosition); check = (mRefit->RefitTrackAsGPU(track, false, true) < 0) ? false : true; - } else { + } else if (mPropagateTrack) { // 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); + } else if (mPropagateParams) { + // propagate the params of the track instead of full propagation + track.rotateParam(o2::math_utils::detail::sector2Angle(sectorIndex)); + check = track.propagateParamTo(xPosition, mFieldNominalGPUBz); } if (!check || std::isnan(track.getParam(1))) { - excludeCl += 0b1000; // 8 for failure of track propagation or refit - } - - 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); - } - // to avoid counting the skipped cluster as a subthreshold cluster - rowIndexOld = rowIndex; - sectorIndexOld = sectorIndex; - continue; + excludeCl += 0b01000000000; // 512 for failure of track propagation or refit } - // get charge values - float chargeTot = cl.getQtot(); - float chargeMax = cl.qMax; - // get threshold const float threshold = mCalibCont.getZeroSupressionThreshold(sectorIndex, rowIndex, pad); - // 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) { - nClsSubThreshROC[0] += missingClusters; - nClsROC[0] += missingClusters; - } else if (stack == GEMstack::OROC1gem) { - nClsSubThreshROC[1] += missingClusters; - nClsROC[1] += missingClusters; - } else if (stack == GEMstack::OROC2gem) { - nClsSubThreshROC[2] += missingClusters; - nClsROC[2] += missingClusters; - } else if (stack == GEMstack::OROC3gem) { - nClsSubThreshROC[3] += missingClusters; - nClsROC[3] += missingClusters; - } - } - } else { - if (stack == GEMstack::IROCgem) { - nClsSubThreshROC[0] += missingClusters; - nClsROC[0] += missingClusters; - } else if (stack == GEMstack::OROC1gem) { - nClsSubThreshROC[1] += missingClusters; - nClsROC[1] += missingClusters; - } else if (stack == GEMstack::OROC2gem) { - nClsSubThreshROC[2] += missingClusters; - nClsROC[2] += missingClusters; - } else if (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; @@ -306,8 +309,6 @@ void CalculatedEdx::calculatedEdx(o2::tpc::TrackTPC& track, dEdxInfo& output, fl }; // get gain - float gain = 1.0f; - float gainResidual = 1.0f; if ((correctionMask & CorrectionFlags::GainFull) == CorrectionFlags::GainFull) { gain = mCalibCont.getGain(sectorIndex, rowIndex, pad); }; @@ -318,8 +319,6 @@ void CalculatedEdx::calculatedEdx(o2::tpc::TrackTPC& track, dEdxInfo& output, fl 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()); @@ -331,28 +330,115 @@ void CalculatedEdx::calculatedEdx(o2::tpc::TrackTPC& track, dEdxInfo& output, fl }; }; - // set the min charge - if (chargeTot < minChargeTot) { - minChargeTot = chargeTot; - }; - - if (chargeMax < minChargeMax) { - minChargeMax = chargeMax; - }; - // 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; + const float time = clTime - track.getTime0(); // ToDo: get correct time from ITS-TPC track if possible if ((correctionMask & CorrectionFlags::dEdxSC) == CorrectionFlags::dEdxSC) { scCorr = mSCdEdxCorrection.getCorrection(time, sectorIndex, rowIndex, pad); if (scCorr > 0) { chargeTot /= scCorr; }; - if (corrMax > 0) { + if (scCorr > 0) { chargeMax /= scCorr; }; } + // for debugging + if (mDebug) { + const o2::gpu::GPUTPCGeometry gpuGeom; + const float localX = gpuGeom.Row2X(rowIndex); + const float localY = gpuGeom.LinearPad2Y(sectorIndex, rowIndex, clPad); + const LocalPosition2D l2D{localX, localY}; + const auto g2D = Mapper::LocalToGlobal(l2D, Sector(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 = track; + + // one row per cluster, tagged with the running track index so rows can be grouped back to the track's + // "dEdxDebugTrack" row; trackParam holds the parameters after refit/propagation to this cluster's row, + // so they differ from cluster to cluster (and from the pristine track stored in "dEdxDebugTrack") + (*debugStreamer) << "dEdxDebugCl" + << "trackIndex=" << mDebugTrackIndex + << "trackParam=" << trackParam + << "cl=" << cl + << "excludeCl=" << excludeCl + << "region=" << region + << "rowIndex=" << rowIndex + << "sectorIndex=" << sectorIndex + << "stack=" << stackNumber + << "localX=" << localX + << "localY=" << localY + << "globalX=" << globalX + << "globalY=" << globalY + << "isShared=" << isShared + << "isCombined=" << isCombined + << "topologyCorr=" << effectiveLength + << "topologyCorrTot=" << effectiveLengthTot + << "topologyCorrMax=" << effectiveLengthMax + << "gain=" << gain + << "gainResidual=" << gainResidual + << "residualCorrTot=" << corrTot + << "residualCorrMax=" << corrMax + << "scCorr=" << scCorr + << "occupancy=" << occupancy + << "\n"; + }; + + if (excludeCl != 0) { + // to avoid counting the skipped cluster as a subthreshold cluster + rowIndexOld = rowIndex; + sectorIndexOld = sectorIndex; + continue; + } + + // set the 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 < minChargeTot) { + minChargeTot = chargeTot; + }; + + if (chargeMax < minChargeMax) { + minChargeMax = chargeMax; + }; + + // find missing clusters + int missingClusters = rowIndex - rowIndexOld - 1; + if ((missingClusters > 0) && (missingClusters <= mMaxMissingCl)) { + if ((clusterMask & ClusterFlags::ExcludeSectorBoundaries) == ClusterFlags::ExcludeSectorBoundaries) { + if (sectorIndexOld == sectorIndex) { + if (stack == GEMstack::IROCgem) { + nClsSubThreshROC[0] += missingClusters; + nClsROC[0] += missingClusters; + } else if (stack == GEMstack::OROC1gem) { + nClsSubThreshROC[1] += missingClusters; + nClsROC[1] += missingClusters; + } else if (stack == GEMstack::OROC2gem) { + nClsSubThreshROC[2] += missingClusters; + nClsROC[2] += missingClusters; + } else if (stack == GEMstack::OROC3gem) { + nClsSubThreshROC[3] += missingClusters; + nClsROC[3] += missingClusters; + } + } + } else { + if (stack == GEMstack::IROCgem) { + nClsSubThreshROC[0] += missingClusters; + nClsROC[0] += missingClusters; + } else if (stack == GEMstack::OROC1gem) { + nClsSubThreshROC[1] += missingClusters; + nClsROC[1] += missingClusters; + } else if (stack == GEMstack::OROC2gem) { + nClsSubThreshROC[2] += missingClusters; + nClsROC[2] += missingClusters; + } else if (stack == GEMstack::OROC3gem) { + nClsSubThreshROC[3] += missingClusters; + nClsROC[3] += missingClusters; + } + } + }; + if (stack == GEMstack::IROCgem) { chargeTotROC[0].emplace_back(chargeTot); chargeMaxROC[0].emplace_back(chargeMax); @@ -374,31 +460,8 @@ void CalculatedEdx::calculatedEdx(o2::tpc::TrackTPC& track, dEdxInfo& output, fl 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); - }; + rowIndexOld = rowIndex; + sectorIndexOld = sectorIndex; } // number of clusters @@ -442,41 +505,39 @@ void CalculatedEdx::calculatedEdx(o2::tpc::TrackTPC& track, dEdxInfo& output, fl output.dEdxMaxOROC3 = getTruncMean(chargeMaxROC[3], low, high); output.dEdxMaxTPC = getTruncMean(chargeMaxROC[4], low, high); - // for debugging - if (mDebug) { - if (mStreamer == nullptr) { - setStreamer(debugRootFile); + // calculate average cl occupancy for the track per TPC region; skip clusters where getOccupancy() had no data (sentinel -1), + // otherwise a single such entry would poison the sum via unsigned overflow + 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; } + } - (*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"; + // for debugging: one row per track, with the track as it was before refit/propagation touched it, + // summarizing the accepted clusters; 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) { + (*debugStreamer) << "dEdxDebugTrack" + << "trackIndex=" << mDebugTrackIndex + << "track=" << trackOrig + << "output=" << output + << "averageOcc=" << averageOcc + << "nCl=" << clustersByRow.size() + << "minChargeTot=" << minChargeTot + << "minChargeMax=" << minChargeMax + << "chargeTotVector=" << chargeTotVector + << "chargeMaxVector=" << chargeMaxVector + << "occupancy=" << occupancyVector + << "\n"; } } @@ -532,6 +593,36 @@ float CalculatedEdx::getTrackTopologyCorrectionPol(const o2::tpc::TrackTPC& trac return effectiveLength; } +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) { // setup CCDB manager @@ -591,6 +682,13 @@ void CalculatedEdx::loadCalibsFromCCDB(long runNumberOrTimeStamp, const bool isM derMap->rectifyAfterReadingFromFile(); mSCdEdxCorrection.setCorrectionMaps(avgMap, derMap); + + // set the dead channel map + o2::tpc::DeadChannelMapCreator deadCMCreator; + deadCMCreator.init(); + deadCMCreator.load(tRun); + const o2::tpc::CalDet& deadMap = deadCMCreator.getDeadChannelMap(); + mCalibCont.setDeadChannelMap(deadMap); } void CalculatedEdx::loadCalibsFromLocalCCDBFolder(const char* localCCDBFolder) @@ -673,12 +771,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 + ; From a0af1e3644ecf936b13e73d8a51eb4775c2575f2 Mon Sep 17 00:00:00 2001 From: tubagundem Date: Fri, 31 Jul 2026 12:14:32 +0200 Subject: [PATCH 2/2] Fix cluster traversal order and propagation robustness in CalculatedEdx --- .../include/TPCCalibration/CalculatedEdx.h | 3 +- .../TPC/calibration/src/CalculatedEdx.cxx | 42 ++++++++++++++----- 2 files changed, 34 insertions(+), 11 deletions(-) diff --git a/Detectors/TPC/calibration/include/TPCCalibration/CalculatedEdx.h b/Detectors/TPC/calibration/include/TPCCalibration/CalculatedEdx.h index c25850c406abf..7107f7c65f992 100644 --- a/Detectors/TPC/calibration/include/TPCCalibration/CalculatedEdx.h +++ b/Detectors/TPC/calibration/include/TPCCalibration/CalculatedEdx.h @@ -163,7 +163,8 @@ class CalculatedEdx /// fill missing clusters with minimum charge (method=0) or minimum charge/2 (method=1) void fillMissingClusters(int missingClusters[4], float minChargeTot, float minChargeMax, int method, std::array, 5>& chargeTotROC, std::array, 5>& chargeMaxROC); - void handleSameRowClusters(o2::tpc::TrackTPC& track, std::map, std::vector>& clustersByRow, std::map, o2::tpc::ClusterNative>& combinedClustersByRow, std::map>& clusterReferencesByIndex); + /// \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 + void handleSameRowClusters(o2::tpc::TrackTPC& track, std::vector>& rowOrder, std::map, std::vector>& clustersByRow, std::map, o2::tpc::ClusterNative>& combinedClustersByRow, std::map>& clusterReferencesByIndex); /// 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 diff --git a/Detectors/TPC/calibration/src/CalculatedEdx.cxx b/Detectors/TPC/calibration/src/CalculatedEdx.cxx index 79921f2d9000b..417ac3b9eaa48 100644 --- a/Detectors/TPC/calibration/src/CalculatedEdx.cxx +++ b/Detectors/TPC/calibration/src/CalculatedEdx.cxx @@ -76,7 +76,7 @@ void CalculatedEdx::fillMissingClusters(int missingClusters[4], float minChargeT } } -void CalculatedEdx::handleSameRowClusters(o2::tpc::TrackTPC& track, std::map, std::vector>& clustersByRow, std::map, o2::tpc::ClusterNative>& combinedClustersByRow, std::map>& clusterReferencesByIndex) +void CalculatedEdx::handleSameRowClusters(o2::tpc::TrackTPC& track, std::vector>& rowOrder, std::map, std::vector>& clustersByRow, std::map, o2::tpc::ClusterNative>& combinedClustersByRow, std::map>& clusterReferencesByIndex) { // get number of clusters const int nClusters = track.getNClusterReferences(); @@ -91,8 +91,13 @@ void CalculatedEdx::handleSameRowClusters(o2::tpc::TrackTPC& track, std::map> rowOrder; std::map, std::vector> clustersByRow; std::map, o2::tpc::ClusterNative> combinedClustersByRow; std::map> clusterReferencesByIndex; - handleSameRowClusters(track, clustersByRow, combinedClustersByRow, clusterReferencesByIndex); + handleSameRowClusters(track, rowOrder, clustersByRow, combinedClustersByRow, clusterReferencesByIndex); o2::utils::TreeStreamRedirector* debugStreamer = nullptr; o2::tpc::TrackTPC trackOrig; @@ -195,8 +201,9 @@ void CalculatedEdx::calculatedEdx(o2::tpc::TrackTPC& track, dEdxInfo& output, Av trackOrig = track; // pristine track, before refit/propagation mutates it cluster-by-cluster below } - // loop over the clusters - for (const auto& [rowKey, clusterIndices] : clustersByRow) { + // loop over the clusters in the track's true physical row-traversal order (rowOrder) + for (const auto& rowKey : rowOrder) { + const auto& clusterIndices = clustersByRow.at(rowKey); const unsigned char rowIndex = rowKey.second; int clusterIdx = clusterIndices[0]; const o2::tpc::ClusterNative& clConst = track.getCluster(*mTPCTrackClIdxVecInput, clusterIdx, *mClusterIndex); @@ -280,12 +287,27 @@ void CalculatedEdx::calculatedEdx(o2::tpc::TrackTPC& track, dEdxInfo& output, Av check = (mRefit->RefitTrackAsGPU(track, false, true) < 0) ? false : true; } else if (mPropagateTrack) { // 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); + // snapshot the fit state first and roll it back on failure + // rotate() rejecting the frame change, or PropagateToXBxByBz failing mid-step e.g. its material-LUT lookup driving the state into an unphysical regime for a difficult trajectory + // so that a single bad row doesn't leave the track corrupted for every subsequent row's propagation attempt + const o2::track::TrackParCov trackBackup = track; + 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) { + static_cast(track) = trackBackup; + } } else if (mPropagateParams) { - // propagate the params of the track instead of full propagation - track.rotateParam(o2::math_utils::detail::sector2Angle(sectorIndex)); - check = track.propagateParamTo(xPosition, mFieldNominalGPUBz); + // propagate the params of the track instead of full propagation; same rollback rationale as mPropagateTrack above + const o2::track::TrackParCov trackBackup = track; + check = track.rotateParam(o2::math_utils::detail::sector2Angle(sectorIndex)); + if (check) { + check = track.propagateParamTo(xPosition, mFieldNominalGPUBz); + } + if (!check) { + static_cast(track) = trackBackup; + } } if (!check || std::isnan(track.getParam(1))) {