diff --git a/Detectors/Calibration/include/DetectorsCalibration/IntegratedClusterCalibrator.h b/Detectors/Calibration/include/DetectorsCalibration/IntegratedClusterCalibrator.h index 9720142d391b1..30cf2880c00c5 100644 --- a/Detectors/Calibration/include/DetectorsCalibration/IntegratedClusterCalibrator.h +++ b/Detectors/Calibration/include/DetectorsCalibration/IntegratedClusterCalibrator.h @@ -331,17 +331,18 @@ struct TimeSeriesdEdx { }; struct TimeSeriesITSTPC { - float mVDrift = 0; ///< drift velocity in cm/us - float mPressure = 0; ///< pressure - float mTemperature = 0; ///< temperature - TimeSeries mTSTPC; ///< TPC standalone DCAs - TimeSeries mTSITSTPC; ///< ITS-TPC standalone DCAs - ITSTPC_Matching mITSTPCAll; ///< ITS-TPC matching efficiency for ITS standalone + afterburner - ITSTPC_Matching mITSTPCStandalone; ///< ITS-TPC matching efficiency for ITS standalone - ITSTPC_Matching mITSTPCAfterburner; ///< ITS-TPC matchin efficiency fir ITS afterburner - TimeSeriesdEdx mdEdxQTot; ///< time series for dE/dx qTot monitoring - TimeSeriesdEdx mdEdxQMax; ///< time series for dE/dx qMax monitoring - std::vector mOccupancyMapTPC; ///< cluster occupancy map + float mVDrift = 0; ///< drift velocity in cm/us + float mPressure = 0; ///< pressure + float mTemperature = 0; ///< temperature + TimeSeries mTSTPC; ///< TPC standalone DCAs + TimeSeries mTSITSTPC; ///< ITS-TPC standalone DCAs + ITSTPC_Matching mITSTPCAll; ///< ITS-TPC matching efficiency for ITS standalone + afterburner + ITSTPC_Matching mITSTPCStandalone; ///< ITS-TPC matching efficiency for ITS standalone + ITSTPC_Matching mITSTPCAfterburner; ///< ITS-TPC matchin efficiency fir ITS afterburner + TimeSeriesdEdx mdEdxQTot; ///< time series for dE/dx qTot monitoring + TimeSeriesdEdx mdEdxQMax; ///< time series for dE/dx qMax monitoring + std::vector mOccupancyMapTPC; ///< cluster occupancy map + std::vector> mSecEdgeFlucCorr; ///< applied sector edge fluctuation correction std::vector nPrimVertices; ///< number of primary vertices std::vector nPrimVertices_ITS; ///< number of primary vertices selected with ITS cut 0.2 CDBTypeMap{ {CDBType::CalScaler, "TPC/Calib/Scaler"}, {CDBType::CalScalerWeights, "TPC/Calib/ScalerWeights"}, {CDBType::CalMShape, "TPC/Calib/MShapePotential"}, + {CDBType::CalSecEdgeCorrection, "TPC/Calib/CorrectionMapSecEdgeFluc"}, + {CDBType::CalSecEdgeInfo, "TPC/Calib/SecEdgeFlucInfo"}, // correction maps loader params {CDBType::CorrMapParam, "TPC/Calib/CorrMapParam"}, // distortion maps diff --git a/Detectors/TPC/calibration/CMakeLists.txt b/Detectors/TPC/calibration/CMakeLists.txt index 6aeb497c1cf23..b0b2704d7ea00 100644 --- a/Detectors/TPC/calibration/CMakeLists.txt +++ b/Detectors/TPC/calibration/CMakeLists.txt @@ -61,6 +61,7 @@ o2_add_library(TPCCalibration src/CMVContainer.cxx src/CorrectionMapsLoader.cxx src/CMVHelper.cxx + src/SectorEdgeFluctuations.cxx PUBLIC_LINK_LIBRARIES O2::DataFormatsTPC O2::TPCBaseRecSim O2::TPCReconstruction ROOT::Minuit Microsoft.GSL::GSL @@ -121,7 +122,8 @@ o2_target_root_dictionary(TPCCalibration include/TPCCalibration/PressureTemperatureHelper.h include/TPCCalibration/CMVContainer.h include/TPCCalibration/CorrectionMapsLoader.h - include/TPCCalibration/CMVHelper.h) + include/TPCCalibration/CMVHelper.h + include/TPCCalibration/SectorEdgeFluctuations.h) o2_add_test_root_macro(macro/comparePedestalsAndNoise.C PUBLIC_LINK_LIBRARIES O2::TPCBaseRecSim diff --git a/Detectors/TPC/calibration/include/TPCCalibration/CorrectionMapsLoader.h b/Detectors/TPC/calibration/include/TPCCalibration/CorrectionMapsLoader.h index 32a61225fe82f..9db0f3ce6f5c2 100644 --- a/Detectors/TPC/calibration/include/TPCCalibration/CorrectionMapsLoader.h +++ b/Detectors/TPC/calibration/include/TPCCalibration/CorrectionMapsLoader.h @@ -19,6 +19,7 @@ #include #include "CorrectionMapsHelper.h" #include "CorrectionMapsTypes.h" +#include "TPCCalibration/SectorEdgeFluctuations.h" namespace o2 { @@ -47,14 +48,19 @@ class CorrectionMapsLoader : public o2::gpu::CorrectionMapsHelper void checkMeanScaleConsistency(float meanLumi, float threshold) const; static void requestCCDBInputs(std::vector& inputs, const o2::tpc::CorrectionMapsGloOpts& gloOpts); + void enableSecEdgeFlucCorrection(const bool enable = true) { mApplySecEdgeFlucCorr = enable; } + bool applySecEdgeFlucCorrection() const { return mApplySecEdgeFlucCorr; } + const auto& getSectorEdgeFlucInfo() const { return mSecEdgeFlucInfo; } protected: static void addOption(std::vector& options, o2::framework::ConfigParamSpec&& osp); static void addInput(std::vector& inputs, o2::framework::InputSpec&& isp); - float mInstLumiCTPFactor = 1.0; // multiplicative factor for inst. lumi - int mLumiCTPSource = 0; // 0: main, 1: alternative CTP lumi source - bool mIDC2CTPFallbackActive = false; // flag indicating that fallback from IDC to CTP scaling is active + float mInstLumiCTPFactor = 1.0; // multiplicative factor for inst. lumi + int mLumiCTPSource = 0; // 0: main, 1: alternative CTP lumi source + bool mIDC2CTPFallbackActive = false; // flag indicating that fallback from IDC to CTP scaling is active + o2::tpc::SectorEdgeFluctuations mSecEdgeFlucInfo; // definition of sector edge fluctuation distortion map scaling + bool mApplySecEdgeFlucCorr = true; // flag indicating if sector edge fluctuation correction is enabled }; } // namespace tpc diff --git a/Detectors/TPC/calibration/include/TPCCalibration/SectorEdgeFluctuations.h b/Detectors/TPC/calibration/include/TPCCalibration/SectorEdgeFluctuations.h new file mode 100644 index 0000000000000..e7fa7afc958e5 --- /dev/null +++ b/Detectors/TPC/calibration/include/TPCCalibration/SectorEdgeFluctuations.h @@ -0,0 +1,108 @@ +// 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 SectorEdgeFluctuations.h +/// \brief Class to parse and query time-dependent TPC sector edge fluctuation intervals +/// \author Matthias Kleiner + +#ifndef ALICEO2_TPC_SECTOREDGEFLUCTUATIONS_H +#define ALICEO2_TPC_SECTOREDGEFLUCTUATIONS_H + +#include +#include +#include +#include +#include "Rtypes.h" + +class TTree; + +namespace o2::tpc +{ + +/// One time interval during which a set of TPC sectors has edge fluctuations. +/// Each sector carries an optional scaling factor (default 1.0). +struct SectorEdgeInterval { + Long64_t startTimeMS{0}; ///< interval start, Unix time in ms + Long64_t endTimeMS{0}; ///< interval end, Unix time in ms + std::vector> sectors; ///< {o2SectorId, scalingFactor} + ClassDefNV(SectorEdgeInterval, 1); +}; + +/// Parses and queries TPC sector edge fluctuation intervals from a CSV text file. +/// +/// Expected line format (comma-separated, whitespace around tokens is ignored): +/// runNumber, startMS, endMS, durationMS, label[, SectorID[=scale], ...] +/// +/// The sector list is optional. If omitted (or all tokens fail to parse), the +/// interval is applied to all 36 sectors (A0-A17 and C0-C17) with scale 1.0. +/// +/// Examples: +/// 560352,1732244770094,1732244771094,1000,edge distortions,A3 +/// 560352,1732244771344,1732244776344,5000,edge distortions,A3=1.2,C0,C1=0.6,C2,C3 +/// +/// If the same sector appears in multiple overlapping intervals at a queried +/// timestamp, the scale from the interval with the latest end-time is used. +class SectorEdgeFluctuations +{ + public: + /// Load intervals from file. Clears any previously loaded data. + /// Throws std::runtime_error if the file cannot be opened. + bool loadFromCSVFile(const std::string& filename); + + /// dump this object to a file + /// \param file output file + /// \param name name of the output object + void dumpToFile(const char* file, const char* name = "ccdb_object", const char* brName = "SectorEdgeFluctuation"); + + /// load from input file (which were written using the dumpToFile method) + /// \param inpf input file + /// \param name name of the object in the file + void loadFromFile(const char* inpf, const char* name = "ccdb_object", const int iEntry = 0, const char* brName = "SectorEdgeFluctuation"); + + /// set this object from input tree + void setFromTree(TTree& tree, const int iEntry = 0, const char* brName = "SectorEdgeFluctuation"); + + /// Returns all {o2SectorId, scalingFactor} pairs active for the given run + /// at the given Unix timestamp (milliseconds). Returns empty if none are active + /// or if the run is not known. + std::vector> getSectorsAtTime(int run, Long64_t timestampMS) const; + + /// Convert a sector string such as "A3" or "C14" to the O2 integer sector + /// index (0-35). Returns -1 on parse error. + static int parseSectorId(const std::string& sectorStr); + + /// Total number of intervals across all runs. + size_t size() const + { + size_t n = 0; + for (const auto& [run, v] : mIntervals) { + n += v.size(); + } + return n; + } + + /// number of total runs stored + size_t getNRuns() const { return mIntervals.size(); } + bool empty() const { return mIntervals.empty(); } + + /// get stored data + const auto& getIntervals() const { return mIntervals; } + + private: + /// Per-run intervals, each sorted by startTimeMS. + std::map> mIntervals; + + ClassDefNV(SectorEdgeFluctuations, 1); +}; + +} // namespace o2::tpc + +#endif // ALICEO2_TPC_SECTOREDGEFLUCTUATIONS_H diff --git a/Detectors/TPC/calibration/include/TPCCalibration/TPCFastSpaceChargeCorrectionHelper.h b/Detectors/TPC/calibration/include/TPCCalibration/TPCFastSpaceChargeCorrectionHelper.h index 40c5634b4f1e8..08e0c5aa23202 100644 --- a/Detectors/TPC/calibration/include/TPCCalibration/TPCFastSpaceChargeCorrectionHelper.h +++ b/Detectors/TPC/calibration/include/TPCCalibration/TPCFastSpaceChargeCorrectionHelper.h @@ -41,6 +41,9 @@ using namespace o2::gpu; class TPCFastSpaceChargeCorrectionHelper { + public: + using SectorScales = std::array; + public: /// _____________ Constructors / destructors __________________________ @@ -115,15 +118,32 @@ class TPCFastSpaceChargeCorrectionHelper /// initialise inverse transformation from linear combination of several input corrections void initInverse(std::vector& corrections, const std::vector& scaling, bool prn); - /// merge several corrections + /// weighted add of several corrections /// \param mainCorrection main correction /// \param scale scaling factor for the main correction /// \param additionalCorrections vector of pairs of additional corrections and their scaling factors - /// \param prn printout flag /// \return main correction merged with additional corrections + void addCorrections( + o2::gpu::TPCFastSpaceChargeCorrection& mainCorrection, double scale, + const std::vector>& additionalCorrections); + + /// weighted add of several corrections with sector-dependent scaling factors + /// \param mainCorrection main correction + /// \param scale scaling factor for the main correction + /// \param additionalCorrections vector of pairs of additional corrections and their sector-dependent scaling factors + /// \return main correction merged with additional corrections + void addCorrections( + o2::gpu::TPCFastSpaceChargeCorrection& mainCorrection, SectorScales scale, + const std::vector>& additionalCorrections); + + /// merge of two corrections sector-wise + /// \param destinationCorrection main correction to which the source correction will be added + /// \param sourceCorrection correction to be added to the main correction + /// \param sectors vector of sector indices for which the correction will be added + /// \return main correction merged with the source correction void mergeCorrections( - o2::gpu::TPCFastSpaceChargeCorrection& mainCorrection, float scale, - const std::vector>& additionalCorrections, bool prn); + o2::gpu::TPCFastSpaceChargeCorrection& destinationCorrection, const o2::gpu::TPCFastSpaceChargeCorrection& sourceCorrection, + const std::vector& sectors); /// how far the voxel mean is allowed to be outside of the voxel (1.1 means 10%) void setVoxelMeanValidityRange(double range) diff --git a/Detectors/TPC/calibration/src/CorrectionMapsLoader.cxx b/Detectors/TPC/calibration/src/CorrectionMapsLoader.cxx index c8bdfa0f99350..51435bc04da77 100644 --- a/Detectors/TPC/calibration/src/CorrectionMapsLoader.cxx +++ b/Detectors/TPC/calibration/src/CorrectionMapsLoader.cxx @@ -35,6 +35,12 @@ void CorrectionMapsLoader::extractCCDBInputs(ProcessingContext& pc, float tpcSca if (lumiMode != LumiScaleMode::NoCorrection) { pc.inputs().get("tpcCorrMapRef"); } + + if (mApplySecEdgeFlucCorr) { + pc.inputs().get("tpcCorrMapSecFluc"); + pc.inputs().get("tpSecFlucInfo"); + } + const int maxDumRep = 5; int dumRep = 0; o2::ctp::LumiInfo lumiObj; @@ -93,11 +99,11 @@ void CorrectionMapsLoader::requestCCDBInputs(std::vector& inputs, con { LOGP(info, "Requesting CCDB inputs for TPC correction maps with lumiType={} and lumiMode={}", static_cast(gloOpts.lumiType), static_cast(gloOpts.lumiMode)); if (gloOpts.lumiMode == LumiScaleMode::Linear) { - addInput(inputs, {"tpcCorrMap", "TPC", "CorrMap", 0, Lifetime::Condition, ccdbParamSpec(CDBTypeMap.at(CDBType::CalCorrMap), {}, 1)}); // time-dependent - addInput(inputs, {"tpcCorrMapRef", "TPC", "CorrMapRef", 0, Lifetime::Condition, ccdbParamSpec(CDBTypeMap.at(CDBType::CalCorrMapRef), {}, 0)}); // load once + addInput(inputs, {"tpcCorrMap", o2::header::gDataOriginTPC, "CorrMap", 0, Lifetime::Condition, ccdbParamSpec(CDBTypeMap.at(CDBType::CalCorrMap), {}, 1)}); // time-dependent + addInput(inputs, {"tpcCorrMapRef", o2::header::gDataOriginTPC, "CorrMapRef", 0, Lifetime::Condition, ccdbParamSpec(CDBTypeMap.at(CDBType::CalCorrMapRef), {}, 0)}); // load once } else if (gloOpts.lumiMode == LumiScaleMode::DerivativeMap) { - addInput(inputs, {"tpcCorrMap", "TPC", "CorrMap", 0, Lifetime::Condition, ccdbParamSpec(CDBTypeMap.at(CDBType::CalCorrMap), {}, 1)}); // time-dependent - addInput(inputs, {"tpcCorrMapRef", "TPC", "CorrMapRef", 0, Lifetime::Condition, ccdbParamSpec(CDBTypeMap.at(CDBType::CalCorrDerivMap), {}, 1)}); // time-dependent + addInput(inputs, {"tpcCorrMap", o2::header::gDataOriginTPC, "CorrMap", 0, Lifetime::Condition, ccdbParamSpec(CDBTypeMap.at(CDBType::CalCorrMap), {}, 1)}); // time-dependent + addInput(inputs, {"tpcCorrMapRef", o2::header::gDataOriginTPC, "CorrMapRef", 0, Lifetime::Condition, ccdbParamSpec(CDBTypeMap.at(CDBType::CalCorrDerivMap), {}, 1)}); // time-dependent } else if (gloOpts.lumiMode == LumiScaleMode::DerivativeMapMC) { // for MC corrections addInput(inputs, {"tpcCorrMap", "TPC", "CorrMap", 0, Lifetime::Condition, ccdbParamSpec(CDBTypeMap.at(CDBType::CalCorrMapMC), {}, 1)}); // time-dependent @@ -110,11 +116,17 @@ void CorrectionMapsLoader::requestCCDBInputs(std::vector& inputs, con LOG(fatal) << "Correction mode unknown! Choose either 0 (default) or 1 (derivative map) for flag corrmap-lumi-mode."; } + // load sector edge fluctuation correction only for data + if (gloOpts.enableSecEdgeFlucCorrection) { + addInput(inputs, {"tpcCorrMapSecFluc", o2::header::gDataOriginTPC, "CorrMapSecFluc", 0, Lifetime::Condition, ccdbParamSpec(CDBTypeMap.at(CDBType::CalSecEdgeCorrection), {}, 1)}); // time-dependent + addInput(inputs, {"tpSecFlucInfo", o2::header::gDataOriginTPC, "InfoMapSecFluc", 0, Lifetime::Condition, ccdbParamSpec(CDBTypeMap.at(CDBType::CalSecEdgeInfo), {}, 1)}); // time-dependent + } + if (gloOpts.requestCTPLumi) { addInput(inputs, {"CTPLumi", "CTP", "LUMI", 0, Lifetime::Timeframe}); } - addInput(inputs, {"tpcCorrPar", "TPC", "CorrMapParam", 0, Lifetime::Condition, ccdbParamSpec(CDBTypeMap.at(CDBType::CorrMapParam), {}, 0)}); // load once + addInput(inputs, {"tpcCorrPar", o2::header::gDataOriginTPC, "CorrMapParam", 0, Lifetime::Condition, ccdbParamSpec(CDBTypeMap.at(CDBType::CorrMapParam), {}, 0)}); // load once } //________________________________________________________ @@ -136,7 +148,7 @@ void CorrectionMapsLoader::addOption(std::vector& options, Conf //________________________________________________________ bool CorrectionMapsLoader::accountCCDBInputs(const ConcreteDataMatcher& matcher, void* obj) { - if (matcher == ConcreteDataMatcher("TPC", "CorrMap", 0)) { + if (matcher == ConcreteDataMatcher(o2::header::gDataOriginTPC, "CorrMap", 0)) { setCorrMap((o2::gpu::TPCFastTransform*)obj); mCorrMap->rectifyAfterReadingFromFile(); mCorrMap->setCTP2IDCFallBackThreshold(o2::tpc::CorrMapParam::Instance().CTP2IDCFallBackThreshold); @@ -165,7 +177,7 @@ bool CorrectionMapsLoader::accountCCDBInputs(const ConcreteDataMatcher& matcher, setUpdatedMap(); return true; } - if (matcher == ConcreteDataMatcher("TPC", "CorrMapRef", 0)) { + if (matcher == ConcreteDataMatcher(o2::header::gDataOriginTPC, "CorrMapRef", 0)) { setCorrMapRef((o2::gpu::TPCFastTransform*)obj); mCorrMapRef->rectifyAfterReadingFromFile(); mCorrMapRef->setCTP2IDCFallBackThreshold(o2::tpc::CorrMapParam::Instance().CTP2IDCFallBackThreshold); @@ -194,7 +206,7 @@ bool CorrectionMapsLoader::accountCCDBInputs(const ConcreteDataMatcher& matcher, setUpdatedMapRef(); return true; } - if (matcher == ConcreteDataMatcher("TPC", "CorrMapParam", 0)) { + if (matcher == ConcreteDataMatcher(o2::header::gDataOriginTPC, "CorrMapParam", 0)) { const auto& par = o2::tpc::CorrMapParam::Instance(); mMeanLumiOverride = par.lumiMean; // negative value switches off corrections !!! mMeanLumiRefOverride = par.lumiMeanRef; @@ -225,6 +237,18 @@ bool CorrectionMapsLoader::accountCCDBInputs(const ConcreteDataMatcher& matcher, canUseCorrections() ? "ON" : "OFF", lumiS[scaleType], mMeanLumiOverride, mMeanLumiRefOverride, static_cast(getLumiScaleMode()), mLumiCTPSource, mInstCTPLumiOverride, mInstLumiCTPFactor); } + if (matcher == ConcreteDataMatcher(o2::header::gDataOriginTPC, "CorrMapSecFluc", 0)) { + setCorrMapSecEdgeFluc((o2::gpu::TPCFastTransform*)obj); + mCorrMapSecEdgeFluc->rectifyAfterReadingFromFile(); + setUpdatedMapSecEdgeFluc(); + return true; + } + if (matcher == ConcreteDataMatcher(o2::header::gDataOriginTPC, "InfoMapSecFluc", 0)) { + LOGP(info, "Updating TPC sector edge fluctuation info"); + mSecEdgeFlucInfo.setFromTree(*((TTree*)obj)); + LOGP(info, "Loaded sector edge fluctuation information with {} intervals for {} runs", mSecEdgeFlucInfo.size(), mSecEdgeFlucInfo.getNRuns()); + return true; + } return false; } diff --git a/Detectors/TPC/calibration/src/CorrectionMapsOptions.cxx b/Detectors/TPC/calibration/src/CorrectionMapsOptions.cxx index 5518d680420ca..c28b5a5732560 100644 --- a/Detectors/TPC/calibration/src/CorrectionMapsOptions.cxx +++ b/Detectors/TPC/calibration/src/CorrectionMapsOptions.cxx @@ -33,6 +33,7 @@ CorrectionMapsGloOpts CorrectionMapsOptions::parseGlobalOptions(const o2::framew tpcopt.lumiMode = static_cast(lumiModeVal); tpcopt.enableMShapeCorrection = opts.get("enable-M-shape-correction"); + tpcopt.enableSecEdgeFlucCorrection = opts.get("enable-sec-edge-fluc-correction"); tpcopt.requestCTPLumi = !opts.get("disable-ctp-lumi-request"); tpcopt.checkCTPIDCconsistency = !opts.get("disable-lumi-type-consistency-check"); if (!tpcopt.requestCTPLumi && tpcopt.lumiType == LumiScaleType::CTPLumi) { @@ -49,6 +50,7 @@ void CorrectionMapsOptions::addGlobalOptions(std::vector& optio addOption(options, ConfigParamSpec{"enable-M-shape-correction", o2::framework::VariantType::Bool, false, {"Enable M-shape distortion correction"}}); addOption(options, ConfigParamSpec{"disable-ctp-lumi-request", o2::framework::VariantType::Bool, false, {"do not request CTP lumi (regardless what is used for corrections)"}}); addOption(options, ConfigParamSpec{"disable-lumi-type-consistency-check", o2::framework::VariantType::Bool, false, {"disable check of selected CTP or IDC scaling source being consistent with the map"}}); + addOption(options, ConfigParamSpec{"enable-sec-edge-fluc-correction", o2::framework::VariantType::Bool, false, {"Enable sector edge fluctuation correction"}}); } void CorrectionMapsOptions::addOption(std::vector& options, ConfigParamSpec&& osp) diff --git a/Detectors/TPC/calibration/src/SectorEdgeFluctuations.cxx b/Detectors/TPC/calibration/src/SectorEdgeFluctuations.cxx new file mode 100644 index 0000000000000..26357da044ee1 --- /dev/null +++ b/Detectors/TPC/calibration/src/SectorEdgeFluctuations.cxx @@ -0,0 +1,253 @@ +// 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 SectorEdgeFluctuations.cxx +/// \brief Class to parse and query time-dependent TPC sector edge fluctuation intervals + +#include "TPCCalibration/SectorEdgeFluctuations.h" +#include "DataFormatsTPC/Defs.h" +#include "TTree.h" +#include "TFile.h" + +#include +#include +#include +#include +#include + +#include "Framework/Logger.h" + +namespace o2::tpc +{ + +int SectorEdgeFluctuations::parseSectorId(const std::string& sectorStr) +{ + if (sectorStr.size() < 2) { + return -1; + } + + const char side = std::toupper(static_cast(sectorStr[0])); + if (side != 'A' && side != 'C') { + return -1; + } + + int num = -1; + try { + num = std::stoi(sectorStr.substr(1)); + } catch (...) { + return -1; + } + + if (num < 0 || num > 17) { + return -1; + } + + return (side == 'A') ? num : (num + SECTORSPERSIDE); +} + +bool SectorEdgeFluctuations::loadFromCSVFile(const std::string& filename) +{ + mIntervals.clear(); + + std::ifstream file(filename); + if (!file.is_open()) { + LOGP(error, "SectorEdgeFluctuations: cannot open file: {}", filename); + return false; + } + + std::string line; + int lineNum = 0; + int nSkipped = 0; + int nLoaded = 0; + + while (std::getline(file, line)) { + ++lineNum; + + // skip empty lines and comments + const auto firstNonSpace = line.find_first_not_of(" \t\r\n"); + if (firstNonSpace == std::string::npos || line[firstNonSpace] == '#') { + continue; + } + + // tokenise on comma + std::vector tokens; + { + std::stringstream ss(line); + std::string tok; + while (std::getline(ss, tok, ',')) { + // trim leading/trailing whitespace + const auto s = tok.find_first_not_of(" \t\r\n"); + const auto e = tok.find_last_not_of(" \t\r\n"); + tokens.push_back((s == std::string::npos) ? "" : tok.substr(s, e - s + 1)); + } + } + + // minimum: runNum(0), startMS(1), endMS(2), duration(3), label(4); sectors are optional + if (tokens.size() < 5) { + LOGP(warning, "SectorEdgeFluctuations: skipping malformed line {}", lineNum); + ++nSkipped; + continue; + } + + int run = -1; + try { + run = std::stoi(tokens[0]); + } catch (...) { + LOGP(warning, "SectorEdgeFluctuations: cannot parse run number on line {}", lineNum); + ++nSkipped; + continue; + } + + SectorEdgeInterval interval; + try { + interval.startTimeMS = std::stoll(tokens[1]); + interval.endTimeMS = std::stoll(tokens[2]); + } catch (...) { + LOGP(warning, "SectorEdgeFluctuations: cannot parse timestamps on line {}", lineNum); + ++nSkipped; + continue; + } + + if (interval.endTimeMS < interval.startTimeMS) { + LOGP(warning, "SectorEdgeFluctuations: end < start on line {}, skipping", lineNum); + ++nSkipped; + continue; + } + + // tokens[4] is the human-readable label; sectors start at index 5 + for (size_t i = 5; i < tokens.size(); ++i) { + std::string sectorStr = tokens[i]; + float scale = 1.0f; + + // parse optional "SectorID=scale" suffix + const auto eqPos = sectorStr.find('='); + if (eqPos != std::string::npos) { + try { + scale = std::stof(sectorStr.substr(eqPos + 1)); + } catch (...) { + LOGP(warning, "SectorEdgeFluctuations: cannot parse scale in '{}' on line {}, using 1.0", tokens[i], lineNum); + } + sectorStr = sectorStr.substr(0, eqPos); + } + + const int sectorId = parseSectorId(sectorStr); + if (sectorId < 0) { + LOGP(warning, "SectorEdgeFluctuations: unknown sector '{}' on line {}, skipping token", sectorStr, lineNum); + continue; + } + + // deduplicate: last occurrence in the line wins + auto dup = std::find_if(interval.sectors.begin(), interval.sectors.end(), [sectorId](const std::pair& p) { return p.first == sectorId; }); + if (dup != interval.sectors.end()) { + dup->second = scale; + } else { + interval.sectors.emplace_back(sectorId, scale); + } + } + + if (interval.sectors.empty()) { + // no sector tokens (or all invalid): apply interval to all 36 sectors + const int nSec = SECTORSPERSIDE * SIDES; + for (int s = 0; s < nSec; ++s) { + interval.sectors.emplace_back(s, 1.0f); + } + } + mIntervals[run].push_back(std::move(interval)); + ++nLoaded; + } + + // sort each run's intervals by start time so getSectorsAtTime can break early + for (auto& [run, intervals] : mIntervals) { + std::sort(intervals.begin(), intervals.end(), [](const SectorEdgeInterval& a, const SectorEdgeInterval& b) { + return a.startTimeMS < b.startTimeMS; + }); + } + + LOGP(info, "SectorEdgeFluctuations: loaded {} intervals for {} run(s) from '{}' ({} lines skipped)", nLoaded, mIntervals.size(), filename, nSkipped); + return true; +} + +std::vector> SectorEdgeFluctuations::getSectorsAtTime(int run, Long64_t timestampMS) const +{ + const auto runIt = mIntervals.find(run); + if (runIt == mIntervals.end()) { + return {}; + } + + // Collect all sectors whose interval is active at timestampMS. + // When the same sector appears in multiple overlapping intervals, keep the + // scale from the interval with the latest endTimeMS (most specific). + // sectorBestScale: sectorId -> {scale, endTimeMS} + std::map> sectorBestScale; + + const auto& intervals = runIt->second; + const auto endIt = std::upper_bound(intervals.begin(), intervals.end(), timestampMS, [](Long64_t ts, const SectorEdgeInterval& iv) { return ts < iv.startTimeMS; }); + + for (auto it = intervals.begin(); it != endIt; ++it) { + if (it->endTimeMS < timestampMS) { + continue; + } + for (const auto& [sector, scale] : it->sectors) { + auto sit = sectorBestScale.find(sector); + if (sit == sectorBestScale.end() || it->endTimeMS > sit->second.second) { + sectorBestScale[sector] = {scale, it->endTimeMS}; + } + } + } + + std::vector> result; + result.reserve(sectorBestScale.size()); + for (const auto& [sector, scaleAndEnd] : sectorBestScale) { + result.emplace_back(sector, scaleAndEnd.first); + } + return result; +} + +void SectorEdgeFluctuations::dumpToFile(const char* file, const char* name, const char* brName) +{ + TFile out(file, "RECREATE"); + TTree tree(name, name); + tree.SetAutoSave(0); + tree.Branch(brName, this); + tree.Fill(); + tree.Write(); + out.Close(); +} + +void SectorEdgeFluctuations::loadFromFile(const char* inpf, const char* name, const int iEntry, const char* brName) +{ + TFile inp(inpf, "READ"); + if (inp.IsZombie() || !inp.IsOpen()) { + LOGP(error, "SectorEdgeFluctuations: cannot open file '{}'", inpf); + return; + } + TTree* tree = dynamic_cast(inp.Get(name)); + if (!tree) { + LOGP(error, "SectorEdgeFluctuations: object '{}' not found or not a TTree in '{}'", name, inpf); + return; + } + setFromTree(*tree, iEntry, brName); +} + +void SectorEdgeFluctuations::setFromTree(TTree& tree, const int iEntry, const char* brName) +{ + SectorEdgeFluctuations* msecFlucTmp = this; + tree.SetBranchAddress(brName, &msecFlucTmp); + const int entries = tree.GetEntries(); + if (entries > iEntry) { + tree.GetEntry(iEntry); + } else { + LOGP(error, "SectorEdgeFluctuation not found in input file"); + } + tree.SetBranchAddress(brName, nullptr); +} + +} // namespace o2::tpc diff --git a/Detectors/TPC/calibration/src/TPCCalibrationLinkDef.h b/Detectors/TPC/calibration/src/TPCCalibrationLinkDef.h index 847ae5ad7d788..740d3f9138e57 100644 --- a/Detectors/TPC/calibration/src/TPCCalibrationLinkDef.h +++ b/Detectors/TPC/calibration/src/TPCCalibrationLinkDef.h @@ -128,4 +128,7 @@ #pragma link C++ class o2::tpc::CMVPerTF + ; #pragma link C++ class o2::tpc::CMVPerTFCompressed + ; +#pragma link C++ class o2::tpc::SectorEdgeFluctuations + ; +#pragma link C++ class o2::tpc::SectorEdgeInterval + ; + #endif diff --git a/Detectors/TPC/calibration/src/TPCFastSpaceChargeCorrectionHelper.cxx b/Detectors/TPC/calibration/src/TPCFastSpaceChargeCorrectionHelper.cxx index fa2f8434f80ee..2ec64765f243a 100644 --- a/Detectors/TPC/calibration/src/TPCFastSpaceChargeCorrectionHelper.cxx +++ b/Detectors/TPC/calibration/src/TPCFastSpaceChargeCorrectionHelper.cxx @@ -188,7 +188,7 @@ void TPCFastSpaceChargeCorrectionHelper::fillSpaceChargeCorrectionFromMap(TPCFas } } } // row - }; // thread + }; // thread std::vector threads(mNthreads); @@ -301,7 +301,7 @@ std::unique_ptr TPCFastSpaceChargeCorrectionHelper } } } // row - }; // thread + }; // thread std::vector threads(mNthreads); @@ -994,7 +994,7 @@ void TPCFastSpaceChargeCorrectionHelper::initInverse(std::vector threads(mNthreads); @@ -1013,19 +1013,38 @@ void TPCFastSpaceChargeCorrectionHelper::initInverse(std::vector>& additionalCorrections, bool /*prn*/) +void TPCFastSpaceChargeCorrectionHelper::addCorrections( + o2::gpu::TPCFastSpaceChargeCorrection& mainCorrection, double mainScale, + const std::vector>& additionalCorrections) +{ + /// weighted add of several corrections + SectorScales mainSectorScale; + mainSectorScale.fill(mainScale); + std::vector> additionalSectorScales; + for (const auto& corr : additionalCorrections) { + SectorScales sectorScale; + sectorScale.fill(corr.second); + additionalSectorScales.emplace_back(corr.first, sectorScale); + } + + addCorrections(mainCorrection, mainSectorScale, additionalSectorScales); +} + +void TPCFastSpaceChargeCorrectionHelper::addCorrections( + o2::gpu::TPCFastSpaceChargeCorrection& mainCorrection, SectorScales mainScale, + const std::vector>& additionalCorrections) { - /// merge several corrections + /// weighted add of several corrections TStopwatch watch; - LOG(info) << "fast space charge correction helper: Merge corrections"; + LOG(info) << "fast space charge correction helper: Add corrections"; const auto& geo = mainCorrection.getGeometry(); for (int sector = 0; sector < geo.getNumberOfSectors(); sector++) { + float secMainScale = mainScale[sector]; + auto myThread = [&](int iThread) { for (int row = iThread; row < geo.getNumberOfRows(); row += mNthreads) { auto& rowInfo = mainCorrection.getRowInfo(row); @@ -1040,8 +1059,7 @@ void TPCFastSpaceChargeCorrectionHelper::mergeCorrections( constexpr int nKnotPar3d = nKnotPar1d * 3; { // scale the main correction - - double parscale[4] = {mainScale, mainScale, mainScale, mainScale * mainScale}; + double parscale[4] = {secMainScale, secMainScale, secMainScale, secMainScale * secMainScale}; for (int iknot = 0, ind = 0; iknot < spline.getNumberOfKnots(); iknot++) { for (int ipar = 0; ipar < nKnotPar1d; ++ipar) { for (int idim = 0; idim < 3; idim++, ind++) { @@ -1072,7 +1090,11 @@ void TPCFastSpaceChargeCorrectionHelper::mergeCorrections( for (int icorr = 0; icorr < additionalCorrections.size(); ++icorr) { const auto& corr = *(additionalCorrections[icorr].first); - double scale = additionalCorrections[icorr].second; + double scale = additionalCorrections[icorr].second[sector]; + if (scale == 0.) { + continue; + } + auto& linfo = corr.getRowInfo(row); double scaleU = rowInfo.gridMeasured.getYscale() / linfo.gridMeasured.getYscale(); @@ -1150,7 +1172,93 @@ void TPCFastSpaceChargeCorrectionHelper::mergeCorrections( } } // sector - float duration = watch.RealTime(); + double duration = watch.RealTime(); + LOGP(info, "Merge of corrections tooks: {}s", duration); +} + +void TPCFastSpaceChargeCorrectionHelper::mergeCorrections(o2::gpu::TPCFastSpaceChargeCorrection& destinationCorrection, + const o2::gpu::TPCFastSpaceChargeCorrection& sourceCorrection, + const std::vector& sectors) +{ + /// merge of two corrections sector-wise + TStopwatch watch; + LOG(info) << "fast space charge correction helper: Merge corrections"; + + const auto& geo = destinationCorrection.getGeometry(); + + for (int sector : sectors) { + if (sector < 0 || sector >= geo.getNumberOfSectors()) { + LOGP(fatal, "Invalid sector number {}. Valid range is [0, {})", sector, geo.getNumberOfSectors()); + continue; + } + auto myThread = [&](int iThread) { + for (int row = iThread; row < geo.getNumberOfRows(); row += mNthreads) { + + { // replace the direct correction + const auto& destSpline = destinationCorrection.getSplineForRow(row); + float* destSplineParameters = destinationCorrection.getCorrectionData(sector, row); + const auto& sourceSpline = sourceCorrection.getSplineForRow(row); + const float* sourceSplineParameters = sourceCorrection.getCorrectionData(sector, row); + + // ensure the splines are compatible + if (destSpline.getGridX1().getNumberOfKnots() != sourceSpline.getGridX1().getNumberOfKnots() || + destSpline.getGridX2().getNumberOfKnots() != sourceSpline.getGridX2().getNumberOfKnots()) { + LOGP(error, "Splines for sector {} row {} are not compatible: number of knots in U or V direction do not match", sector, row); + continue; + } + // replace the destination correction with the source correction for this sector and row + memcpy(destSplineParameters, sourceSplineParameters, destSpline.getNumberOfParameters() * sizeof(float)); + } + + { // replace the inverse correction X + const auto& destSpline = destinationCorrection.getSplineInvXforRow(row); + float* destSplineParameters = destinationCorrection.getCorrectionDataInvX(sector, row); + const auto& sourceSpline = sourceCorrection.getSplineInvXforRow(row); + const float* sourceSplineParameters = sourceCorrection.getCorrectionDataInvX(sector, row); + // ensure the splines are compatible + if (destSpline.getGridX1().getNumberOfKnots() != sourceSpline.getGridX1().getNumberOfKnots() || + destSpline.getGridX2().getNumberOfKnots() != sourceSpline.getGridX2().getNumberOfKnots()) { + LOGP(error, "Inverse X splines for sector {} row {} are not compatible: number of knots in U or V direction do not match", sector, row); + continue; + } + memcpy(destSplineParameters, sourceSplineParameters, destSpline.getNumberOfParameters() * sizeof(float)); + } + + { // replace the inverse correction YZ + const auto& destSpline = destinationCorrection.getSplineInvYZforRow(row); + float* destSplineParameters = destinationCorrection.getCorrectionDataInvYZ(sector, row); + const auto& sourceSpline = sourceCorrection.getSplineInvYZforRow(row); + const float* sourceSplineParameters = sourceCorrection.getCorrectionDataInvYZ(sector, row); + // ensure the splines are compatible + if (destSpline.getGridX1().getNumberOfKnots() != sourceSpline.getGridX1().getNumberOfKnots() || + destSpline.getGridX2().getNumberOfKnots() != sourceSpline.getGridX2().getNumberOfKnots()) { + LOGP(error, "Inverse YZ splines for sector {} row {} are not compatible: number of knots in U or V direction do not match", sector, row); + continue; + } + memcpy(destSplineParameters, sourceSplineParameters, destSpline.getNumberOfParameters() * sizeof(float)); + } + + // replace the sector row info + auto& destSecRowInfo = destinationCorrection.getRowInfo(row); + const auto& sourceSecRowInfo = sourceCorrection.getRowInfo(row); + destSecRowInfo = sourceSecRowInfo; + } // row + }; // thread + + std::vector threads(mNthreads); + + // run n threads + for (int i = 0; i < mNthreads; i++) { + threads[i] = std::thread(myThread, i); + } + + // wait for the threads to finish + for (auto& th : threads) { + th.join(); + } + + } // sector + double duration = watch.RealTime(); LOGP(info, "Merge of corrections tooks: {}s", duration); } diff --git a/Detectors/TPC/workflow/include/TPCWorkflow/TPCTimeSeriesSpec.h b/Detectors/TPC/workflow/include/TPCWorkflow/TPCTimeSeriesSpec.h index d7da0b9acb343..3b14b1c60cccc 100644 --- a/Detectors/TPC/workflow/include/TPCWorkflow/TPCTimeSeriesSpec.h +++ b/Detectors/TPC/workflow/include/TPCWorkflow/TPCTimeSeriesSpec.h @@ -23,7 +23,7 @@ namespace tpc static constexpr header::DataDescription getDataDescriptionTimeSeries() { return header::DataDescription{"TIMESERIES"}; } static constexpr header::DataDescription getDataDescriptionTPCTimeSeriesTFId() { return header::DataDescription{"ITPCTSTFID"}; } -o2::framework::DataProcessorSpec getTPCTimeSeriesSpec(const bool disableWriter, const o2::base::Propagator::MatCorrType matType, const bool enableUnbinnedWriter, o2::dataformats::GlobalTrackID::mask_t src); +o2::framework::DataProcessorSpec getTPCTimeSeriesSpec(const bool disableWriter, const o2::base::Propagator::MatCorrType matType, const bool enableUnbinnedWriter, o2::dataformats::GlobalTrackID::mask_t src, const bool enableSecEdgeFluc); } // end namespace tpc } // end namespace o2 diff --git a/Detectors/TPC/workflow/src/TPCScalerSpec.cxx b/Detectors/TPC/workflow/src/TPCScalerSpec.cxx index 1df192dd5ec00..9ca483531fbcd 100644 --- a/Detectors/TPC/workflow/src/TPCScalerSpec.cxx +++ b/Detectors/TPC/workflow/src/TPCScalerSpec.cxx @@ -45,6 +45,7 @@ class TPCScalerSpec : public Task mTPCCorrMapsLoader.setLumiScaleType(sclOpts.lumiType); mTPCCorrMapsLoader.setLumiScaleMode(sclOpts.lumiMode); mTPCCorrMapsLoader.setCheckCTPIDCConsistency(sclOpts.checkCTPIDCconsistency); + mTPCCorrMapsLoader.enableSecEdgeFlucCorrection(sclOpts.enableSecEdgeFlucCorrection); }; void init(framework::InitContext& ic) final @@ -178,14 +179,20 @@ class TPCScalerSpec : public Task pc.outputs().snapshot(Output{header::gDataOriginCTP, "LUMICTP"}, lumiCTP); } - buildMap(pc); + buildMap(pc, timestamp); } - void buildMap(ProcessingContext& pc) + void buildMap(ProcessingContext& pc, int64_t timestamp) { const auto lumiMode = mTPCCorrMapsLoader.getLumiScaleMode(); o2::gpu::TPCFastTransform finalMap; - std::vector> additionalCorrections; + std::vector> additionalCorrections; + + auto uniformScale = [](double s) { + o2::tpc::TPCFastSpaceChargeCorrectionHelper::SectorScales ss; + ss.fill(s); + return ss; + }; if (lumiMode == LumiScaleMode::NoCorrection) { std::unique_ptr dummy(TPCFastTransformHelperO2::instance()->create(0)); @@ -201,26 +208,41 @@ class TPCScalerSpec : public Task // if standard scaling is used: map(lumi) = (mean_map - ref_map) * lumiScale + ref_map if (lumiMode == LumiScaleMode::Linear) { - const std::vector> step0{{&(corrMapRef->getCorrection()), -1.f}}; + const std::vector> step0{{&(corrMapRef->getCorrection()), -1.f}}; // finalMap = (mean_map - finalMap) - TPCFastSpaceChargeCorrectionHelper::instance()->mergeCorrections(finalMap.getCorrection(), 1, step0, true); + TPCFastSpaceChargeCorrectionHelper::instance()->addCorrections(finalMap.getCorrection(), 1., step0); // finalMap = finalMap * lumiScale + ref_map - const std::vector> step1{{&(corrMapRef->getCorrection()), 1.f}}; - TPCFastSpaceChargeCorrectionHelper::instance()->mergeCorrections(finalMap.getCorrection(), lumiScale, step1, true); + const std::vector> step1{{&(corrMapRef->getCorrection()), 1.}}; + TPCFastSpaceChargeCorrectionHelper::instance()->addCorrections(finalMap.getCorrection(), lumiScale, step1); } else if (lumiMode == LumiScaleMode::DerivativeMap || lumiMode == LumiScaleMode::DerivativeMapMC) { - additionalCorrections.emplace_back(&(corrMapRef->getCorrection()), lumiScale); + additionalCorrections.emplace_back(&(corrMapRef->getCorrection()), uniformScale(lumiScale)); } // if mshape map valid if (!mTPCCorrMapsLoader.isCorrMapMShapeDummy()) { LOGP(info, "Adding M-shape correction to the final map with scaling factor {}", mMShapeScalingFac); - additionalCorrections.emplace_back(&(mTPCCorrMapsLoader.getCorrMapMShape()->getCorrection()), 1.f); + additionalCorrections.emplace_back(&(mTPCCorrMapsLoader.getCorrMapMShape()->getCorrection()), uniformScale(1.)); + } + + // --- sector-edge fluctuation correction --- + if (mTPCCorrMapsLoader.applySecEdgeFlucCorrection()) { + LOGP(info, "Checking for sector edge fluctuation"); + const int currRun = pc.services().get().runNumber; + const auto activeSectors = mTPCCorrMapsLoader.getSectorEdgeFlucInfo().getSectorsAtTime(currRun, static_cast(timestamp)); + if (!activeSectors.empty()) { + LOGP(info, "Adding edge-sector correction for {} active sector(s)", activeSectors.size()); + o2::tpc::TPCFastSpaceChargeCorrectionHelper::SectorScales sectorScales{}; + for (const auto& [sector, scale] : activeSectors) { + sectorScales[sector] = static_cast(scale); + } + additionalCorrections.emplace_back(&mTPCCorrMapsLoader.getCorrMapSecEdgeFluc()->getCorrection(), sectorScales); + } } if (!additionalCorrections.empty()) { - TPCFastSpaceChargeCorrectionHelper::instance()->mergeCorrections(finalMap.getCorrection(), 1, additionalCorrections, true); + TPCFastSpaceChargeCorrectionHelper::instance()->addCorrections(finalMap.getCorrection(), uniformScale(1.), additionalCorrections); } } @@ -319,7 +341,6 @@ o2::framework::DataProcessorSpec getTPCScalerSpec(bool enableIDCs, bool enableMS LOGP(info, "Publishing M-shape correction map"); inputs.emplace_back("mshape", o2::header::gDataOriginTPC, "MSHAPEPOTCCDB", 0, Lifetime::Condition, ccdbParamSpec(o2::tpc::CDBTypeMap.at(o2::tpc::CDBType::CalMShape), {}, 1)); // time-dependent } - auto ccdbRequest = std::make_shared(true, // orbitResetTime false, // GRPECS=true for nHBF per TF false, // GRPLHCIF diff --git a/Detectors/TPC/workflow/src/TPCTimeSeriesSpec.cxx b/Detectors/TPC/workflow/src/TPCTimeSeriesSpec.cxx index 0c0ae72056318..7248569f800dd 100644 --- a/Detectors/TPC/workflow/src/TPCTimeSeriesSpec.cxx +++ b/Detectors/TPC/workflow/src/TPCTimeSeriesSpec.cxx @@ -19,6 +19,7 @@ #include "Framework/ConfigParamRegistry.h" #include "Framework/DataProcessorSpec.h" #include "Framework/ControlService.h" +#include "Framework/CCDBParamSpec.h" #include "TPCWorkflow/ProcessingHelpers.h" #include "TPCBase/Mapper.h" #include "DetectorsBase/GRPGeomHelper.h" @@ -45,6 +46,8 @@ #include "DataFormatsTOF/Cluster.h" #include "DataFormatsFT0/RecPoints.h" #include "TPCCalibration/PressureTemperatureHelper.h" +#include "TPCBaseRecSim/CDBTypes.h" +#include "TPCCalibration/SectorEdgeFluctuations.h" using namespace o2::globaltracking; using GTrackID = o2::dataformats::GlobalTrackID; @@ -61,7 +64,7 @@ class TPCTimeSeries : public Task { public: /// \constructor - TPCTimeSeries(std::shared_ptr req, const bool disableWriter, const o2::base::Propagator::MatCorrType matType, const bool enableUnbinnedWriter, const bool tpcOnly, std::shared_ptr dr) : mCCDBRequest(req), mDisableWriter(disableWriter), mMatType(matType), mUnbinnedWriter(enableUnbinnedWriter), mTPCOnly(tpcOnly), mDataRequest(dr) {}; + TPCTimeSeries(std::shared_ptr req, const bool disableWriter, const o2::base::Propagator::MatCorrType matType, const bool enableUnbinnedWriter, const bool tpcOnly, std::shared_ptr dr, const bool enableSecEdgeFluc) : mCCDBRequest(req), mDisableWriter(disableWriter), mMatType(matType), mUnbinnedWriter(enableUnbinnedWriter), mTPCOnly(tpcOnly), mDataRequest(dr), mEnableSecEdgeFluc(enableSecEdgeFluc) {}; void init(framework::InitContext& ic) final { @@ -132,12 +135,18 @@ class TPCTimeSeries : public Task mVDrift = mTPCVDriftHelper.getVDriftObject().getVDrift(); LOGP(info, "Updated reference drift velocity to: {}", mVDrift); } + if (mEnableSecEdgeFluc) { + pc.inputs().get("tpcSecFlucInfo"); + } mBufferDCA.mVDrift = mVDrift; const int nBins = getNBins(); mTimeMS = o2::base::GRPGeomHelper::instance().getOrbitResetTimeMS() + processing_helpers::getFirstTForbit(pc) * o2::constants::lhc::LHCOrbitMUS / 1000; mRun = processing_helpers::getRunNumber(pc); + if (mEnableSecEdgeFluc) { + mBufferDCA.mSecEdgeFlucCorr = mSecEdgeFlucInfo.getSectorsAtTime(mRun, static_cast(mTimeMS)); + } mBufferDCA.mTemperature = mPTHelper.getMeanTemperature(mTimeMS); mBufferDCA.mPressure = mPTHelper.getPressure(mTimeMS); @@ -875,6 +884,11 @@ class TPCTimeSeries : public Task mTPCVDriftHelper.accountCCDBInputs(matcher, obj); mPTHelper.accountCCDBInputs(matcher, obj); o2::base::GRPGeomHelper::instance().finaliseCCDB(matcher, obj); + if (matcher == ConcreteDataMatcher(o2::header::gDataOriginTPC, "InfoMapSecFluc", 0)) { + LOGP(info, "Updating TPC sector edge fluctuation info"); + mSecEdgeFlucInfo.setFromTree(*((TTree*)obj)); + LOGP(info, "Loaded sector edge fluctuation information with {} intervals for {} runs", mSecEdgeFlucInfo.size(), mSecEdgeFlucInfo.getNRuns()); + } } private: @@ -1039,6 +1053,7 @@ class TPCTimeSeries : public Task const bool mUnbinnedWriter{false}; /// write out additional unbinned data const bool mTPCOnly{false}; ///< produce only TPC variables std::shared_ptr mDataRequest; ///< steers the input + bool mEnableSecEdgeFluc{false}; ///< enable write out of sector edge fluctuations int mPhiBins = SECTORSPERSIDE; ///< number of phi bins int mTglBins{3}; ///< number of tgl bins int mQPtBins{20}; ///< number of qPt bins @@ -1112,6 +1127,7 @@ class TPCTimeSeries : public Task int mRun{}; ///< run number int mMaxOccupancyHistBins{912}; ///< maximum number of occupancy bins PressureTemperatureHelper mPTHelper; ///< helper to extract pressure and temperature from CCDB + o2::tpc::SectorEdgeFluctuations mSecEdgeFlucInfo; ///< definition of sector edge fluctuation distortion map scaling /// check if track passes coarse cuts bool acceptTrack(const TrackTPC& track) const { return std::abs(track.getTgl()) < mMaxTgl; } @@ -1820,7 +1836,7 @@ class TPCTimeSeries : public Task } }; -o2::framework::DataProcessorSpec getTPCTimeSeriesSpec(const bool disableWriter, const o2::base::Propagator::MatCorrType matType, const bool enableUnbinnedWriter, GTrackID::mask_t src) +o2::framework::DataProcessorSpec getTPCTimeSeriesSpec(const bool disableWriter, const o2::base::Propagator::MatCorrType matType, const bool enableUnbinnedWriter, GTrackID::mask_t src, const bool enableSecEdgeFluc) { auto dataRequest = std::make_shared(); bool useMC = false; @@ -1849,6 +1865,10 @@ o2::framework::DataProcessorSpec getTPCTimeSeriesSpec(const bool disableWriter, o2::tpc::VDriftHelper::requestCCDBInputs(dataRequest->inputs); PressureTemperatureHelper::requestCCDBInputs(dataRequest->inputs); + if (enableSecEdgeFluc) { + dataRequest->inputs.emplace_back("tpcSecFlucInfo", o2::header::gDataOriginTPC, "InfoMapSecFluc", 0, Lifetime::Condition, ccdbParamSpec(CDBTypeMap.at(CDBType::CalSecEdgeInfo), {}, 1)); + } + std::vector outputs; outputs.emplace_back(o2::header::gDataOriginTPC, getDataDescriptionTimeSeries(), 0, Lifetime::Sporadic); if (!disableWriter) { @@ -1859,7 +1879,7 @@ o2::framework::DataProcessorSpec getTPCTimeSeriesSpec(const bool disableWriter, "tpc-time-series", dataRequest->inputs, outputs, - AlgorithmSpec{adaptFromTask(ccdbRequest, disableWriter, matType, enableUnbinnedWriter, tpcOnly, dataRequest)}, + AlgorithmSpec{adaptFromTask(ccdbRequest, disableWriter, matType, enableUnbinnedWriter, tpcOnly, dataRequest, enableSecEdgeFluc)}, Options{ {"min-momentum", VariantType::Float, 0.2f, {"Minimum momentum of the tracks"}}, {"min-cluster", VariantType::Int, 80, {"Minimum number of clusters of the tracks"}}, diff --git a/Detectors/TPC/workflow/src/tpc-time-series.cxx b/Detectors/TPC/workflow/src/tpc-time-series.cxx index f7bcf00cb27ea..bcba768a7034c 100644 --- a/Detectors/TPC/workflow/src/tpc-time-series.cxx +++ b/Detectors/TPC/workflow/src/tpc-time-series.cxx @@ -29,7 +29,9 @@ void customize(std::vector& workflowOptions) {"disable-root-output", VariantType::Bool, false, {"disable root-files output writers"}}, {"enable-unbinned-root-output", VariantType::Bool, false, {"writing out unbinned track data"}}, {"track-sources", VariantType::String, std::string{o2::dataformats::GlobalTrackID::ALL}, {"comma-separated list of sources to use"}}, - {"material-type", VariantType::Int, 2, {"Type for the material budget during track propagation: 0=None, 1=Geo, 2=LUT"}}}; + {"material-type", VariantType::Int, 2, {"Type for the material budget during track propagation: 0=None, 1=Geo, 2=LUT"}}, + {"enable-sec-edge-fluc-correction", VariantType::Bool, false, {"Enable sector edge fluctuation correction output"}}}; + std::swap(workflowOptions, options); } @@ -43,7 +45,8 @@ WorkflowSpec defineDataProcessing(ConfigContext const& config) const bool enableUnbinnedWriter = config.options().get("enable-unbinned-root-output"); auto src = o2::dataformats::GlobalTrackID::getSourcesMask(config.options().get("track-sources")); auto materialType = static_cast(config.options().get("material-type")); - workflow.emplace_back(o2::tpc::getTPCTimeSeriesSpec(disableWriter, materialType, enableUnbinnedWriter, src)); + const bool enableSecEdgeFluc = config.options().get("enable-sec-edge-fluc-correction"); + workflow.emplace_back(o2::tpc::getTPCTimeSeriesSpec(disableWriter, materialType, enableUnbinnedWriter, src, enableSecEdgeFluc)); if (!disableWriter) { workflow.emplace_back(o2::tpc::getTPCTimeSeriesWriterSpec()); } diff --git a/GPU/TPCFastTransformation/CorrectionMapsHelper.h b/GPU/TPCFastTransformation/CorrectionMapsHelper.h index 095bb837eaacd..9fdaa7f0e576a 100644 --- a/GPU/TPCFastTransformation/CorrectionMapsHelper.h +++ b/GPU/TPCFastTransformation/CorrectionMapsHelper.h @@ -33,10 +33,12 @@ class CorrectionMapsHelper const o2::gpu::TPCFastTransform* getCorrMap() const { return mCorrMap; } const o2::gpu::TPCFastTransform* getCorrMapRef() const { return mCorrMapRef; } + const o2::gpu::TPCFastTransform* getCorrMapSecEdgeFluc() const { return mCorrMapSecEdgeFluc; } const o2::gpu::TPCFastTransform* getCorrMapMShape() const { return mCorrMapMShape.get(); } void setCorrMap(o2::gpu::TPCFastTransform* m) { mCorrMap = m; } void setCorrMapRef(o2::gpu::TPCFastTransform* m) { mCorrMapRef = m; } + void setCorrMapSecEdgeFluc(o2::gpu::TPCFastTransform* m) { mCorrMapSecEdgeFluc = m; } void setCorrMapMShape(std::unique_ptr&& m); void reportScaling(); @@ -98,6 +100,7 @@ class CorrectionMapsHelper void setUpdatedMap() { mUpdatedFlags |= UpdateFlags::MapBit; } void setUpdatedMapRef() { mUpdatedFlags |= UpdateFlags::MapRefBit; } void setUpdatedMapMShape() { mUpdatedFlags |= UpdateFlags::MapMShapeBit; } + void setUpdatedMapSecEdgeFluc() { mUpdatedFlags |= UpdateFlags::MapSecEdgeFlucBit; } void setUpdatedLumi() { mUpdatedFlags |= UpdateFlags::LumiBit; } void acknowledgeUpdate() { mUpdatedFlags = 0; } void setLumiCTPAvailable(bool v) { mLumiCTPAvailable = v; } @@ -131,7 +134,8 @@ class CorrectionMapsHelper enum UpdateFlags { MapBit = 0x1, MapRefBit = 0x2, LumiBit = 0x4, - MapMShapeBit = 0x10 }; + MapMShapeBit = 0x8, + MapSecEdgeFlucBit = 0x10 }; bool mLumiCTPAvailable = false; // is CTP Lumi available // these 2 are global options, must be set by the workflow global options tpc::LumiScaleType mLumiScaleType = tpc::LumiScaleType::Unset; // use CTP Lumi (1) or TPCScaler (2) for the correction scaling, 0 - no scaling @@ -149,8 +153,9 @@ class CorrectionMapsHelper bool mCheckCTPIDCConsistency{true}; // check of selected CTP or IDC scaling source being consistent with the map o2::gpu::TPCFastTransform* mCorrMap{nullptr}; // current transform o2::gpu::TPCFastTransform* mCorrMapRef{nullptr}; // reference transform + o2::gpu::TPCFastTransform* mCorrMapSecEdgeFluc{nullptr}; // sector edge fluctuation correction map std::unique_ptr mCorrMapMShape{nullptr}; // correction map for M-shape distortions on A-side - ClassDefNV(CorrectionMapsHelper, 6); + ClassDefNV(CorrectionMapsHelper, 7); }; } // namespace o2::gpu diff --git a/GPU/TPCFastTransformation/CorrectionMapsTypes.h b/GPU/TPCFastTransformation/CorrectionMapsTypes.h index 092a2927ebe3e..94ea96c5ac7e5 100644 --- a/GPU/TPCFastTransformation/CorrectionMapsTypes.h +++ b/GPU/TPCFastTransformation/CorrectionMapsTypes.h @@ -42,6 +42,7 @@ struct CorrectionMapsGloOpts { bool enableMShapeCorrection = false; bool requestCTPLumi = true; ///< request CTP Lumi regardless of what is used for corrections scaling bool checkCTPIDCconsistency = true; ///< check the selected CTP or IDC scaling source being consistent with mean scaler of the map + bool enableSecEdgeFlucCorrection = true; ///< enable correction of sector edge fluctuations }; } // namespace o2::tpc #endif