From 66f8031a9b29fa363a43f3f9e20a25579d6dd729 Mon Sep 17 00:00:00 2001 From: tubagundem Date: Thu, 30 Jul 2026 11:52:17 +0200 Subject: [PATCH 1/5] 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/5] 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))) { From bc2f13426d8d56f926ae2fd085e7e7002df953fa Mon Sep 17 00:00:00 2001 From: tubagundem Date: Tue, 4 Aug 2026 14:26:59 +0200 Subject: [PATCH 3/5] Add loadSCCorrMap option to loadCalibsFromCCDB, fetch the MC time gain calibration when isMC is set, remove MC kinematics reader from dEdx calculation macro --- .../include/TPCCalibration/CalculatedEdx.h | 5 +++-- .../TPC/calibration/macro/calculatedEdx.C | 21 ++++--------------- .../TPC/calibration/src/CalculatedEdx.cxx | 18 +++++++++------- 3 files changed, 17 insertions(+), 27 deletions(-) diff --git a/Detectors/TPC/calibration/include/TPCCalibration/CalculatedEdx.h b/Detectors/TPC/calibration/include/TPCCalibration/CalculatedEdx.h index 7107f7c65f992..c11d30b75eb17 100644 --- a/Detectors/TPC/calibration/include/TPCCalibration/CalculatedEdx.h +++ b/Detectors/TPC/calibration/include/TPCCalibration/CalculatedEdx.h @@ -211,8 +211,9 @@ class CalculatedEdx /// load calibration objects from CCDB /// \param runNumberOrTimeStamp run number or time stamp - /// \param isMC set if dEdx space-charge corrections will be loaded for MC or real data - void loadCalibsFromCCDB(long runNumberOrTimeStamp, const bool isMC = false); + /// \param isMC set if dEdx residual and space-charge corrections will be loaded for MC or real data + /// \param loadSCCorrMap set to false to skip loading the space-charge correction maps + void loadCalibsFromCCDB(long runNumberOrTimeStamp, const bool isMC = false, const bool loadSCCorrMap = true); /// load calibration objects from local CCDB folder /// \param localCCDBFolder local CCDB folder diff --git a/Detectors/TPC/calibration/macro/calculatedEdx.C b/Detectors/TPC/calibration/macro/calculatedEdx.C index 2edbea05de5d4..8dc41d93b6705 100644 --- a/Detectors/TPC/calibration/macro/calculatedEdx.C +++ b/Detectors/TPC/calibration/macro/calculatedEdx.C @@ -40,8 +40,6 @@ #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 @@ -91,14 +89,13 @@ std::pair, std::unique_ptr> openTreeOrNull(const s /// \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 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 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", @@ -120,8 +117,7 @@ void calculatedEdx(const std::string dir = ".", 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") + const std::string matchFileName = "o2match_itstpc.root") { ROOT::EnableThreadSafety(); @@ -164,18 +160,12 @@ void calculatedEdx(const std::string dir = ".", 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; @@ -260,7 +250,8 @@ void calculatedEdx(const std::string dir = ".", for (auto& calcdEdx : calcdEdxPerThread) { calcdEdx.setMembers(tpcTrackClIdxVecInput, clusterIndex, &tpcTracks); if (localCCDBFolder.empty()) { - calcdEdx.loadCalibsFromCCDB(timeStamp, isMC); + const bool loadSCCorrMap = (correctionMask & CorrectionFlags::dEdxSC) == CorrectionFlags::dEdxSC; + calcdEdx.loadCalibsFromCCDB(timeStamp, isMC, loadSCCorrMap); } else { calcdEdx.loadCalibsFromLocalCCDBFolder(localCCDBFolder.data()); } @@ -338,11 +329,7 @@ void calculatedEdx(const std::string dir = ".", } 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"; } diff --git a/Detectors/TPC/calibration/src/CalculatedEdx.cxx b/Detectors/TPC/calibration/src/CalculatedEdx.cxx index 417ac3b9eaa48..81dd01ff117e7 100644 --- a/Detectors/TPC/calibration/src/CalculatedEdx.cxx +++ b/Detectors/TPC/calibration/src/CalculatedEdx.cxx @@ -645,7 +645,7 @@ bool CalculatedEdx::isInStackBoundaries(int stackNumber, unsigned char rowIndex, return false; } -void CalculatedEdx::loadCalibsFromCCDB(long runNumberOrTimeStamp, const bool isMC) +void CalculatedEdx::loadCalibsFromCCDB(long runNumberOrTimeStamp, const bool isMC, const bool loadSCCorrMap) { // setup CCDB manager auto& cm = o2::ccdb::BasicCCDBManager::instance(); @@ -675,7 +675,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); @@ -696,14 +696,16 @@ 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 sc correction maps; skip if not needed + if (loadSCCorrMap) { + 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(); + 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); + } // set the dead channel map o2::tpc::DeadChannelMapCreator deadCMCreator; From b66cc623af15faee3beeb773ea6a1b6f2530ce9f Mon Sep 17 00:00:00 2001 From: tubagundem Date: Fri, 7 Aug 2026 15:14:14 +0200 Subject: [PATCH 4/5] Added a function to calculate multiple dEdx settings for the same track without repeating track propagation and introduced minor improvements --- .../include/TPCCalibration/CalculatedEdx.h | 107 ++- .../TPC/calibration/macro/calculatedEdx.C | 124 +++- .../TPC/calibration/src/CalculatedEdx.cxx | 656 ++++++++++++++++-- 3 files changed, 791 insertions(+), 96 deletions(-) diff --git a/Detectors/TPC/calibration/include/TPCCalibration/CalculatedEdx.h b/Detectors/TPC/calibration/include/TPCCalibration/CalculatedEdx.h index c11d30b75eb17..69bb82d3b87f5 100644 --- a/Detectors/TPC/calibration/include/TPCCalibration/CalculatedEdx.h +++ b/Detectors/TPC/calibration/include/TPCCalibration/CalculatedEdx.h @@ -95,6 +95,21 @@ 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 { + float low = 0.015f; ///< lower cluster cut + float high = 0.6f; ///< higher cluster cut + CorrectionFlags correctionMask = CorrectionFlags::TopologyPol | CorrectionFlags::dEdxResidual; ///< corrections to apply + ClusterFlags clusterMask = ClusterFlags::None; ///< clusters to exclude + int subthresholdMethod = 0; ///< subthreshold cluster charge filling method + int stackBoundaryMethod = 0; ///< stack boundary cluster exclusion method + std::string debugRootFile = "dEdxDebug.root"; ///< debug streamer output file used if mDebug is set + 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) +}; + class CalculatedEdx { public: @@ -126,12 +141,6 @@ 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 minChargeMaxThreshold upper limit for the possible minimum charge max in subthreshold treatment - void setMinChargeMaxThreshold(float minChargeMaxThreshold) { mMinChargeMaxThreshold = minChargeMaxThreshold; } - /// 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) @@ -154,14 +163,28 @@ class CalculatedEdx /// \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 number of rows where refit/propagation failed (row.propagationFailed) since the last resetDebugCounters() + long getNPropagationFailed() const { return mNPropagationFailed; } - /// \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; } + + /// \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()/getNRowsProcessed()/getNSubThresholdFilledPerSettings() + void resetDebugCounters() + { + mNPropagationFailed = 0; + mNRowsProcessed = 0; + mNSubThresholdFilledPerSettings.clear(); + } - /// 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); + /// 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 void handleSameRowClusters(o2::tpc::TrackTPC& track, std::vector>& rowOrder, std::map, std::vector>& clustersByRow, std::map, o2::tpc::ClusterNative>& combinedClustersByRow, std::map>& clusterReferencesByIndex); @@ -175,7 +198,16 @@ class CalculatedEdx /// \param high higher cluster cut /// \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, 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"); + /// \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 + 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); + + /// 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 + void calculatedEdxMultipleSettings(TrackTPC& track, std::vector& outputs, AverageOccupancy& averageOcc, const std::vector& settingsList); /// get the truncated mean for the input charge vector and the truncation range low*nCl& rowData, AverageOccupancy& averageOcc); + + /// 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 + void calculatedEdxFromRowData(const std::vector& rowData, const dEdxSettings& settings, size_t settingsIndex, float trackTime0, const o2::tpc::TrackTPC& trackOrig, const AverageOccupancy& averageOcc, dEdxInfo& output); + 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 @@ -272,8 +348,6 @@ class CalculatedEdx 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 (faster than refit) bool mPropagateParams{false}; ///< propagating the parameters instead of full propagation (faster than track propagation) @@ -281,6 +355,9 @@ class CalculatedEdx 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 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 diff --git a/Detectors/TPC/calibration/macro/calculatedEdx.C b/Detectors/TPC/calibration/macro/calculatedEdx.C index 8dc41d93b6705..214b1fa3ff224 100644 --- a/Detectors/TPC/calibration/macro/calculatedEdx.C +++ b/Detectors/TPC/calibration/macro/calculatedEdx.C @@ -12,6 +12,9 @@ /// \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 @@ -72,22 +75,30 @@ std::pair, std::unique_ptr> openTreeOrNull(const s } 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; + 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 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 settingsList list of dEdx settings to evaluate for every track /// \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 debug enable the CalculatedEdx debug streamer, additionally writing one dEdxDebug_t.root file per worker thread with per-cluster information (or dEdxDebug_t_s.root per worker thread per settings entry j, if settingsList has more than one entry) /// \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 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 @@ -96,16 +107,12 @@ std::pair, std::unique_ptr> openTreeOrNull(const s /// \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 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 std::vector settingsList = defaultSettingsList(), const bool useRefit = true, const bool propagateTrack = false, const bool propagateParams = false, @@ -117,7 +124,8 @@ void calculatedEdx(const std::string dir = ".", 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 matchFileName = "o2match_itstpc.root", + const long long maxEvents = -1) { ROOT::EnableThreadSafety(); @@ -125,17 +133,26 @@ void calculatedEdx(const std::string dir = ".", 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); + if (settingsList.empty()) { + LOGP(error, "settingsList must not be empty"); 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; + 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 < 0 || 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; + } + // 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={}", i, s.low, s.high, static_cast(s.correctionMask), static_cast(s.clusterMask), s.subthresholdMethod, s.stackBoundaryMethod); } const std::clock_t c_start = std::clock(); @@ -149,6 +166,16 @@ void calculatedEdx(const std::string dir = ".", calcdEdx.setPropagateParams(propagateParams); } + // one copy of settingsList per thread, with debugRootFile made unique per thread (and per settings entry, if there is more than one) so debug streams from different threads/settings never collide + std::vector> settingsPerThread(nThreads, settingsList); + if (debug) { + for (size_t iThread = 0; iThread < nThreads; iThread++) { + for (size_t iSettings = 0; iSettings < settingsList.size(); iSettings++) { + settingsPerThread[iThread][iSettings].debugRootFile = (settingsList.size() == 1) ? fmt::format("dEdxDebug_t{}.root", iThread) : fmt::format("dEdxDebug_t{}_s{}.root", iThread, iSettings); + } + } + } + auto [tpcFile, tpcTree] = openTreeOrNull(fmt::format("{}/{}", dir, tpcTracksFileName), "tpcrec"); if (!tpcTree) { return; @@ -228,6 +255,10 @@ void calculatedEdx(const std::string dir = ".", 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); @@ -250,7 +281,10 @@ void calculatedEdx(const std::string dir = ".", for (auto& calcdEdx : calcdEdxPerThread) { calcdEdx.setMembers(tpcTrackClIdxVecInput, clusterIndex, &tpcTracks); if (localCCDBFolder.empty()) { - const bool loadSCCorrMap = (correctionMask & CorrectionFlags::dEdxSC) == CorrectionFlags::dEdxSC; + bool loadSCCorrMap = false; + for (const auto& s : settingsList) { + loadSCCorrMap |= (s.correctionMask & CorrectionFlags::dEdxSC) == CorrectionFlags::dEdxSC; + } calcdEdx.loadCalibsFromCCDB(timeStamp, isMC, loadSCCorrMap); } else { calcdEdx.loadCalibsFromLocalCCDBFolder(localCCDBFolder.data()); @@ -261,13 +295,13 @@ void calculatedEdx(const std::string dir = ".", } const size_t nSelectable = isMatchedToITS ? matchTracks.size() : tpcTracks.size(); - LOGP(info, "Processing event {} with {} {} using {} threads", iEvent, nSelectable, isMatchedToITS ? "matched tracks" : "tracks", nThreads); + LOGP(info, "Processing event {} with {} {} using {} threads and {} settings", iEvent, nSelectable, isMatchedToITS ? "matched tracks" : "tracks", nThreads, settingsList.size()); std::vector> tpcOut(nThreads); std::vector> itsTracksOut(nThreads); std::vector> matchTracksOut(nThreads); - std::vector> dEdxOut(nThreads); - std::vector> averageOccOut(nThreads); + std::vector>> dEdxOut(nThreads); // [thread][track][settings] + std::vector> averageOccOut(nThreads); // [thread][track] std::vector> mcLabelOut(nThreads); const size_t chunkSize = (nSelectable + nThreads - 1) / nThreads; @@ -280,7 +314,7 @@ void calculatedEdx(const std::string dir = ".", } threads.emplace_back([&, iThread, start, end]() { auto& calcdEdx = calcdEdxPerThread[iThread]; - const std::string debugFile = fmt::format("dEdxDebug_t{}.root", iThread); + const auto& threadSettingsList = settingsPerThread[iThread]; for (size_t i = start; i < end; i++) { size_t tpcIndex = i; if (isMatchedToITS) { @@ -293,13 +327,13 @@ void calculatedEdx(const std::string dir = ".", matchTracksOut[iThread].emplace_back(itstpc); } - TrackTPC track(tpcTracks[tpcIndex]); // local copy: setRefit()/propagation mutate the track in place - dEdxInfo dEdx; + TrackTPC track(tpcTracks[tpcIndex]); // local copy: refit/propagation inside calculatedEdxMultipleSettings mutate the track in place + std::vector dEdxVec; AverageOccupancy averageOcc; - calcdEdx.calculatedEdx(track, dEdx, averageOcc, lowTruncation, highTruncation, correctionMask, clusterMask, subthresholdMethod, stackBoundaryMethod, debugFile.c_str()); + calcdEdx.calculatedEdxMultipleSettings(track, dEdxVec, averageOcc, threadSettingsList); tpcOut[iThread].emplace_back(track); - dEdxOut[iThread].emplace_back(dEdx); + dEdxOut[iThread].emplace_back(std::move(dEdxVec)); averageOccOut[iThread].emplace_back(averageOcc); if (isMC) { mcLabelOut[iThread].emplace_back(tpcMCTruth[tpcIndex]); @@ -311,15 +345,37 @@ void calculatedEdx(const std::string dir = ".", th.join(); } - // write out sequentially in the main thread: no locking needed since all worker threads have already joined + // per-event summary: refit/propagation failures and how many row gaps were filled as subthreshold clusters per settingsList entry + long nPropagationFailed = 0, nRowsProcessed = 0; + std::vector nSubThresholdFilledPerSettings(settingsList.size(), 0); + for (auto& calcdEdx : calcdEdxPerThread) { + nPropagationFailed += calcdEdx.getNPropagationFailed(); + nRowsProcessed += calcdEdx.getNRowsProcessed(); + const auto& threadSubThresholdFilled = calcdEdx.getNSubThresholdFilledPerSettings(); + for (size_t i = 0; i < threadSubThresholdFilled.size() && i < nSubThresholdFilledPerSettings.size(); i++) { + nSubThresholdFilledPerSettings[i] += threadSubThresholdFilled[i]; + } + 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); + + // write out sequentially in the main thread: no locking needed since all worker threads have already joined; + // one row per track, with one "dEdx" branch per entry in settingsList 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]; + for (size_t iSettings = 0; iSettings < settingsList.size(); iSettings++) { + row << fmt::format("dEdx{}=", iSettings).c_str() << dEdxOut[iThread][i][iSettings]; + } if (tfIDTree) { row << "tfIDInfo=" << tfIDInfo; } diff --git a/Detectors/TPC/calibration/src/CalculatedEdx.cxx b/Detectors/TPC/calibration/src/CalculatedEdx.cxx index 81dd01ff117e7..875d96ce7057a 100644 --- a/Detectors/TPC/calibration/src/CalculatedEdx.cxx +++ b/Detectors/TPC/calibration/src/CalculatedEdx.cxx @@ -55,16 +55,16 @@ 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!"); return; } - float chargeMax = (method == 1) ? minChargeMax / 2.f : minChargeMax; - float chargeTot = (method == 1) ? minChargeTot / 2.f : minChargeTot; for (int roc = 0; roc < 4; roc++) { + 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++) { chargeTotROC[roc].emplace_back(chargeTot); @@ -149,7 +149,525 @@ void CalculatedEdx::handleSameRowClusters(o2::tpc::TrackTPC& track, std::vector< } } -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) +void CalculatedEdx::gatherRowClusterData(o2::tpc::TrackTPC& track, std::vector& rowData, AverageOccupancy& averageOcc) +{ + rowData.clear(); + + // handle same (sector, row) clusters + std::vector> rowOrder; + std::map, std::vector> clustersByRow; + std::map, o2::tpc::ClusterNative> combinedClustersByRow; + std::map> clusterReferencesByIndex; + + handleSameRowClusters(track, rowOrder, clustersByRow, combinedClustersByRow, 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; + 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; + } + + RowClusterData row; + row.cl = cl; + row.clPad = cl.getPad(); + row.clTime = cl.getTime(); + row.chargeTot = cl.getQtot(); + row.chargeMax = cl.getQmax(); + row.sectorIndex = sectorIndex; + row.rowIndex = rowIndex; + row.isCombined = isCombined; + + const unsigned int occupancy = getOccupancy(row.clTime); + row.occupancy = occupancy; + + // check if the cluster is shared + const unsigned int absoluteIndex = mClusterIndex->clusterOffset[sectorIndex][rowIndex] + clusterIndexNumb; + row.isShared = mRefit ? (mTPCRefitterShMap[absoluteIndex] & o2::gpu::GPUTPCGMMergedTrackHit::flagShared) : 0; + + // get region, pad, stack and stack ID + const int region = Mapper::REGION[rowIndex]; + const unsigned char pad = std::clamp(static_cast(row.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}; + const int stackNumber = static_cast(stack); + + row.region = region; + row.pad = pad; + row.stack = stack; + row.stackID = stackID; + row.stackNumber = stackNumber; + + 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); + } + + row.isDeadRegion = mCalibCont.isDead(static_cast(sectorIndex), static_cast(rowIndex), static_cast(pad)); + + // get the x position of the track + const float xPosition = Mapper::instance().getPadCentre(PadPos(rowIndex, 0)).X(); + bool check = true; + if (mRefit) { + // refit this track + mRefit->setTrackReferenceX(xPosition); + 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); snapshot the fit state first and roll it back before each fallback attempt below. A failure here must not leave the track frozen at this row's stale X + // so retry the same target without material corrections and if that also fails, fall back to the simple analytic parameter-only propagation used by mPropagateParams only if all three fail, give up and roll back, marking the row as propagationFailed below + 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; + 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) { + static_cast(track) = trackBackup; + check = track.rotateParam(o2::math_utils::detail::sector2Angle(sectorIndex)); + if (check) { + check = track.propagateParamTo(xPosition, mFieldNominalGPUBz); + } + } + if (!check) { + static_cast(track) = trackBackup; + } + } 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))); + ++mNRowsProcessed; + if (row.propagationFailed) { + ++mNPropagationFailed; + } + + // snapshot of the track state after refit/propagation to this row; reused by calculatedEdxFromRowData() for every settings entry + row.trackSnapshot = track; + + // get threshold and gain + row.threshold = mCalibCont.getZeroSupressionThreshold(sectorIndex, rowIndex, pad); + row.gain = mCalibCont.getGain(sectorIndex, rowIndex, pad); + row.gainResidual = mCalibCont.getResidualGain(sectorIndex, rowIndex, pad); + + // number of rows skipped since the previous row in rowOrder + row.missingClusters = rowIndex - rowIndexOld - 1; + row.sameSectorAsPrevRow = (sectorIndexOld == sectorIndex); + + // 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 && row.sameSectorAsPrevRow) { + const o2::gpu::GPUTPCGeometry gpuGeom; + const RowClusterData& prevRow = rowData.back(); + const float yPrev = gpuGeom.LinearPad2Y(sectorIndex, prevRow.rowIndex, prevRow.clPad); + const float yCur = gpuGeom.LinearPad2Y(sectorIndex, rowIndex, row.clPad); + for (int k = 1; k <= row.missingClusters; ++k) { + const unsigned char missingRow = prevRow.rowIndex + k; + const float frac = static_cast(k) / (row.missingClusters + 1); + const float missingPad = gpuGeom.LinearY2Pad(sectorIndex, missingRow, yPrev + (yCur - yPrev) * 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; + } + } + } + + rowIndexOld = rowIndex; + sectorIndexOld = sectorIndex; + + rowData.emplace_back(std::move(row)); + } + + // 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) +{ + // 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; + 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] = {100000.f, 100000.f, 100000.f, 100000.f}; + float minChargeMaxROC[4] = {100000.f, 100000.f, 100000.f, 100000.f}; + + 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()); + } + + for (const auto& row : rowData) { + if (mDebug) { + occupancyVector.emplace_back(row.occupancy); + } + + // get cluster values + float chargeTot = row.chargeTot; + float chargeMax = row.chargeMax; + + // 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 = row.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) && row.isShared) { + excludeCl += 0b0100000; // 32 for shared cluster + } + if (((settings.clusterMask & ClusterFlags::ExcludeSamePadRowCl) == ClusterFlags::ExcludeSamePadRowCl) && row.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 (row.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 ((settings.correctionMask & CorrectionFlags::TopologyPol) == CorrectionFlags::TopologyPol) { + effectiveLengthTot = getTrackTopologyCorrectionPol(row.trackSnapshot, row.cl, row.region, chargeTot, ChargeType::Tot, row.threshold); + effectiveLengthMax = getTrackTopologyCorrectionPol(row.trackSnapshot, row.cl, row.region, chargeMax, ChargeType::Max, row.threshold); + chargeTot /= effectiveLengthTot; + chargeMax /= effectiveLengthMax; + }; + + // get gain + if ((settings.correctionMask & CorrectionFlags::GainFull) == CorrectionFlags::GainFull) { + gain = row.gain; + }; + if ((settings.correctionMask & CorrectionFlags::GainResidual) == CorrectionFlags::GainResidual) { + gainResidual = row.gainResidual; + }; + chargeTot /= gain * gainResidual; + chargeMax /= gain * gainResidual; + + // get dEdx 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; + }; + }; + + // space-charge dEdx corrections + const float time = row.clTime - 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, row.pad); + if (scCorr > 0) { + chargeTot /= scCorr; + }; + if (scCorr > 0) { + chargeMax /= scCorr; + }; + } + + // for debugging + if (mDebug) { + const o2::gpu::GPUTPCGeometry gpuGeom; + const float localX = gpuGeom.Row2X(row.rowIndex); + const float localY = gpuGeom.LinearPad2Y(row.sectorIndex, row.rowIndex, row.clPad); + 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=" << row.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=" << row.isShared + << "isCombined=" << row.isCombined + << "topologyCorr=" << effectiveLength + << "topologyCorrTot=" << effectiveLengthTot + << "topologyCorrMax=" << effectiveLengthMax + << "gain=" << gain + << "gainResidual=" << gainResidual + << "residualCorrTot=" << corrTot + << "residualCorrMax=" << corrMax + << "scCorr=" << scCorr + << "occupancy=" << row.occupancy + << "\n"; + }; + + // find missing clusters + const int missingClusters = row.missingClusters; + if ((missingClusters > 0) && (missingClusters <= mMaxMissingCl) && !row.missingClusterGapDeadOrEdge) { + if ((settings.clusterMask & ClusterFlags::ExcludeSectorBoundaries) == ClusterFlags::ExcludeSectorBoundaries) { + if (row.sameSectorAsPrevRow) { + if (row.stack == GEMstack::IROCgem) { + nClsSubThreshROC[0] += missingClusters; + nClsROC[0] += missingClusters; + } else if (row.stack == GEMstack::OROC1gem) { + nClsSubThreshROC[1] += missingClusters; + nClsROC[1] += missingClusters; + } else if (row.stack == GEMstack::OROC2gem) { + nClsSubThreshROC[2] += missingClusters; + nClsROC[2] += missingClusters; + } else if (row.stack == GEMstack::OROC3gem) { + nClsSubThreshROC[3] += missingClusters; + nClsROC[3] += missingClusters; + } + } + } else { + if (row.stack == GEMstack::IROCgem) { + nClsSubThreshROC[0] += missingClusters; + nClsROC[0] += missingClusters; + } else if (row.stack == GEMstack::OROC1gem) { + nClsSubThreshROC[1] += missingClusters; + nClsROC[1] += missingClusters; + } else if (row.stack == GEMstack::OROC2gem) { + nClsSubThreshROC[2] += missingClusters; + nClsROC[2] += missingClusters; + } else if (row.stack == GEMstack::OROC3gem) { + nClsSubThreshROC[3] += missingClusters; + nClsROC[3] += missingClusters; + } + } + }; + + 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 (chargeMax < minChargeMaxROC[row.stackNumber]) { + minChargeMaxROC[row.stackNumber] = chargeMax; + }; + + 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); + } + + // number of clusters + output.NHitsSubThresholdIROC = nClsROC[0]; + output.NHitsSubThresholdOROC1 = nClsROC[1]; + output.NHitsSubThresholdOROC2 = nClsROC[2]; + output.NHitsSubThresholdOROC3 = nClsROC[3]; + + // the gaps found above are always treated as subthreshold clusters + 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 (((settings.clusterMask & ClusterFlags::ExcludeSubthresholdCl) == ClusterFlags::None)) { + float cappedMinChargeTotROC[4], cappedMinChargeMaxROC[4]; + for (int roc = 0; roc < 4; roc++) { + cappedMinChargeTotROC[roc] = std::min(minChargeTotROC[roc], settings.maxSubthresholdChargeTot); + cappedMinChargeMaxROC[roc] = std::min(minChargeMaxROC[roc], settings.maxSubthresholdChargeMax); + } + 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]; + } + + // 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], settings.low, settings.high); + output.dEdxTotOROC1 = getTruncMean(chargeTotROC[1], settings.low, settings.high); + output.dEdxTotOROC2 = getTruncMean(chargeTotROC[2], settings.low, settings.high); + output.dEdxTotOROC3 = getTruncMean(chargeTotROC[3], settings.low, settings.high); + output.dEdxTotTPC = getTruncMean(chargeTotROC[4], settings.low, settings.high); + + output.dEdxMaxIROC = getTruncMean(chargeMaxROC[0], settings.low, settings.high); + output.dEdxMaxOROC1 = getTruncMean(chargeMaxROC[1], settings.low, settings.high); + output.dEdxMaxOROC2 = getTruncMean(chargeMaxROC[2], settings.low, settings.high); + output.dEdxMaxOROC3 = getTruncMean(chargeMaxROC[3], settings.low, settings.high); + output.dEdxMaxTPC = getTruncMean(chargeMaxROC[4], settings.low, settings.high); + + // for debugging: one row per track, with the track as it was before refit/propagation touched it, per-cluster rows were already written to the "dEdxDebugCl" tree above (each with its own propagated track parameters) and can be grouped back to this row via trackIndex + if (mDebug) { + float minChargeTot = minChargeTotROC[0], minChargeMax = minChargeMaxROC[0]; + for (int roc = 1; roc < 4; roc++) { + minChargeTot = (minChargeTotROC[roc] < minChargeTot) ? minChargeTotROC[roc] : minChargeTot; + minChargeMax = (minChargeMaxROC[roc] < minChargeMax) ? minChargeMaxROC[roc] : minChargeMax; + } + (*debugStreamer) << "dEdxDebugTrack" + << "trackIndex=" << mDebugTrackIndex + << "track=" << trackOrig + << "output=" << output + << "averageOcc=" << averageOcc + << "nCl=" << rowData.size() + << "minChargeTot=" << minChargeTot + << "minChargeMax=" << minChargeMax + << "chargeTotVector=" << chargeTotVector + << "chargeMaxVector=" << chargeMaxVector + << "occupancy=" << occupancyVector + << "\n"; + } +} + +void CalculatedEdx::calculatedEdxMultipleSettings(o2::tpc::TrackTPC& track, std::vector& outputs, AverageOccupancy& averageOcc, const std::vector& settingsList) +{ + 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]); + } +} + +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) { // NHits and NHitsSubthreshold values per region int nClsROC[4] = {0, 0, 0, 0}; @@ -171,8 +689,10 @@ void CalculatedEdx::calculatedEdx(o2::tpc::TrackTPC& track, dEdxInfo& output, Av // for tracking missing clusters unsigned char rowIndexOld = 255; unsigned char sectorIndexOld = 255; - float minChargeTot = 100000.f; - float minChargeMax = 100000.f; + float clPadOld = 0.f; // previous row's cluster pad, used by the dead/edge gap veto below + // per-region running minimum charge + float minChargeTotROC[4] = {100000.f, 100000.f, 100000.f, 100000.f}; + float minChargeMaxROC[4] = {100000.f, 100000.f, 100000.f, 100000.f}; // corrections float effectiveLength = 1.0f; @@ -286,15 +806,27 @@ void CalculatedEdx::calculatedEdx(o2::tpc::TrackTPC& track, dEdxInfo& output, Av mRefit->setTrackReferenceX(xPosition); 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) - // 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 + // propagate this track to the plane X=xk (cm) in the field "b" (kG); snapshot the fit state first and roll it back before each fallback attempt below. A failure here must not leave the track frozen at this row's stale X + // so retry the same target without material corrections and if that also fails, fall back to the simple analytic parameter-only propagation used by mPropagateParams only if all three fail, give up and roll back, marking the row as propagationFailed below 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; + 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) { + static_cast(track) = trackBackup; + check = track.rotateParam(o2::math_utils::detail::sector2Angle(sectorIndex)); + if (check) { + check = track.propagateParamTo(xPosition, mFieldNominalGPUBz); + } + } if (!check) { static_cast(track) = trackBackup; } @@ -374,8 +906,7 @@ void CalculatedEdx::calculatedEdx(o2::tpc::TrackTPC& track, dEdxInfo& output, Av 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 + // slice to the base parametrization 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 @@ -385,6 +916,8 @@ void CalculatedEdx::calculatedEdx(o2::tpc::TrackTPC& track, dEdxInfo& output, Av << "trackIndex=" << mDebugTrackIndex << "trackParam=" << trackParam << "cl=" << cl + << "chargeTot=" << chargeTot + << "chargeMax=" << chargeMax << "excludeCl=" << excludeCl << "region=" << region << "rowIndex=" << rowIndex @@ -408,26 +941,33 @@ void CalculatedEdx::calculatedEdx(o2::tpc::TrackTPC& track, dEdxInfo& output, Av << "\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; - }; + // find missing clusters - deliberately evaluated before (independent of) this row's own excludeCl status + int missingClusters = rowIndex - rowIndexOld - 1; - if (chargeMax < minChargeMax) { - minChargeMax = chargeMax; - }; + // 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 + bool missingClusterGapDeadOrEdge = false; + if (missingClusters > 0 && missingClusters <= mMaxMissingCl && sectorIndexOld == sectorIndex) { + const o2::gpu::GPUTPCGeometry gpuGeom; + const float yPrev = gpuGeom.LinearPad2Y(sectorIndex, rowIndexOld, clPadOld); + const float yCur = gpuGeom.LinearPad2Y(sectorIndex, rowIndex, clPad); + for (int k = 1; k <= missingClusters; ++k) { + const unsigned char missingRow = rowIndexOld + k; + const float frac = static_cast(k) / (missingClusters + 1); + const float missingPad = gpuGeom.LinearY2Pad(sectorIndex, missingRow, yPrev + (yCur - yPrev) * frac); + if (missingPad < 0.f || missingPad >= gpuGeom.NPads(missingRow)) { + 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))) { + missingClusterGapDeadOrEdge = true; + break; + } + } + } - // find missing clusters - int missingClusters = rowIndex - rowIndexOld - 1; - if ((missingClusters > 0) && (missingClusters <= mMaxMissingCl)) { + if ((missingClusters > 0) && (missingClusters <= mMaxMissingCl) && !missingClusterGapDeadOrEdge) { if ((clusterMask & ClusterFlags::ExcludeSectorBoundaries) == ClusterFlags::ExcludeSectorBoundaries) { if (sectorIndexOld == sectorIndex) { if (stack == GEMstack::IROCgem) { @@ -461,6 +1001,24 @@ void CalculatedEdx::calculatedEdx(o2::tpc::TrackTPC& track, dEdxInfo& output, Av } }; + if (excludeCl != 0) { + // to avoid counting the skipped cluster itself as a real or subthreshold hit + rowIndexOld = rowIndex; + sectorIndexOld = sectorIndex; + clPadOld = clPad; + 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[stackNumber]) { + minChargeTotROC[stackNumber] = chargeTot; + }; + + if (chargeMax < minChargeMaxROC[stackNumber]) { + minChargeMaxROC[stackNumber] = chargeMax; + }; + if (stack == GEMstack::IROCgem) { chargeTotROC[0].emplace_back(chargeTot); chargeMaxROC[0].emplace_back(chargeMax); @@ -484,6 +1042,7 @@ void CalculatedEdx::calculatedEdx(o2::tpc::TrackTPC& track, dEdxInfo& output, Av rowIndexOld = rowIndex; sectorIndexOld = sectorIndex; + clPadOld = clPad; } // number of clusters @@ -492,22 +1051,19 @@ void CalculatedEdx::calculatedEdx(o2::tpc::TrackTPC& track, dEdxInfo& output, Av 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); + 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 (((clusterMask & ClusterFlags::ExcludeSubthresholdCl) == ClusterFlags::None)) { + float cappedMinChargeTotROC[4], cappedMinChargeMaxROC[4]; + for (int roc = 0; roc < 4; roc++) { + cappedMinChargeTotROC[roc] = std::min(minChargeTotROC[roc], maxSubthresholdChargeTot); + cappedMinChargeMaxROC[roc] = std::min(minChargeMaxROC[roc], maxSubthresholdChargeMax); } - } else { - output.NHitsIROC = nClsROC[0]; - output.NHitsOROC1 = nClsROC[1]; - output.NHitsOROC2 = nClsROC[2]; - output.NHitsOROC3 = nClsROC[3]; + fillMissingClusters(nClsSubThreshROC, cappedMinChargeTotROC, cappedMinChargeMaxROC, subthresholdMethod, chargeTotROC, chargeMaxROC); } // copy corrected cluster charges @@ -545,9 +1101,15 @@ void CalculatedEdx::calculatedEdx(o2::tpc::TrackTPC& track, dEdxInfo& output, Av } // 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 + // 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) { + // minChargeTot/Max are now tracked per region (see minChargeTotROC/minChargeMaxROC above); report the + // smallest of the four here so the "dEdxDebugTrack" tree keeps its existing single-scalar branches + 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; + } (*debugStreamer) << "dEdxDebugTrack" << "trackIndex=" << mDebugTrackIndex << "track=" << trackOrig From 7f2e56a28541556f16ec7bb6d09cab52d79865bc Mon Sep 17 00:00:00 2001 From: tubagundem Date: Fri, 7 Aug 2026 15:18:09 +0200 Subject: [PATCH 5/5] Fix formatting --- .../TPC/calibration/include/TPCCalibration/CalculatedEdx.h | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/Detectors/TPC/calibration/include/TPCCalibration/CalculatedEdx.h b/Detectors/TPC/calibration/include/TPCCalibration/CalculatedEdx.h index 69bb82d3b87f5..a1ea651dbbbe3 100644 --- a/Detectors/TPC/calibration/include/TPCCalibration/CalculatedEdx.h +++ b/Detectors/TPC/calibration/include/TPCCalibration/CalculatedEdx.h @@ -355,9 +355,9 @@ class CalculatedEdx 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 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() + long mNPropagationFailed{0}; ///< number of rows where refit/propagation failed 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