From 8ed5bbfacf9092ecace8710a9e29ad6ca31f1317 Mon Sep 17 00:00:00 2001 From: Matthias Kleiner Date: Thu, 7 Dec 2023 17:19:25 +0100 Subject: [PATCH] TPC: Adding possibility to correct for V-shape distortions The time dependent distortions that we observe in the DCAs on the A-side can be similarly corrected as the space-charge distortions. For this a new model is added which we can linearly scale up or down according to current V-shape observed in DCA. Further changes: - TPC IDC scaler adding helper function to split the scalers, when writing to file - fixing timestamp in TPC time series Adding v-shape correction object to several places Adding option to recalculate the inverse correction during reconstruction Pass needed options to tpc-scalers workflow from dpl-workflow.sh connect TPCScalerSpec to dependent standalone workflows Adding option to load correction map from local macro Adding updated v-shape scaler with z-pos Updating v-shape calibration procedure temporary fix compiling Renaming V-shape tp M-shape --- .../src/barrel-alignment-workflow.cxx | 5 + .../src/cosmics-match-workflow.cxx | 4 + .../src/secondary-vertexing-workflow.cxx | 5 +- .../src/tof-matcher-workflow.cxx | 5 +- .../src/tpcits-match-workflow.cxx | 4 + .../study/CMakeLists.txt | 1 + .../study/src/tpc-track-study-workflow.cxx | 4 + Detectors/TPC/base/include/TPCBase/CDBTypes.h | 2 + Detectors/TPC/calibration/CMakeLists.txt | 4 +- .../TPCCalibration/CorrectionMapsLoader.h | 16 +- .../TPCFastSpaceChargeCorrectionHelper.h | 9 +- .../include/TPCCalibration/TPCMShapeScaler.h | 99 +++++++++++ .../include/TPCCalibration/TPCScaler.h | 22 ++- .../calibration/src/CorrectionMapsLoader.cxx | 88 +++++++++- .../calibration/src/TPCCalibrationLinkDef.h | 2 + .../TPCFastSpaceChargeCorrectionHelper.cxx | 28 ++- .../TPC/calibration/src/TPCMShapeScaler.cxx | 161 ++++++++++++++++++ Detectors/TPC/calibration/src/TPCScaler.cxx | 67 ++++++++ .../include/TPCWorkflow/TPCScalerSpec.h | 2 +- Detectors/TPC/workflow/src/RecoWorkflow.cxx | 8 + Detectors/TPC/workflow/src/TPCScalerSpec.cxx | 68 ++++++-- .../TPC/workflow/src/TPCTimeSeriesSpec.cxx | 7 +- .../workflow/src/tpc-calib-gainmap-tracks.cxx | 7 +- Detectors/TPC/workflow/src/tpc-scaler.cxx | 10 +- .../workflow/src/trd-tracking-workflow.cxx | 4 + GPU/GPUTracking/DataTypes/GPUDataTypes.h | 1 + GPU/GPUTracking/Global/GPUChainTracking.cxx | 12 ++ GPU/GPUTracking/Global/GPUChainTracking.h | 2 + GPU/GPUTracking/Global/GPUChainTrackingIO.cxx | 10 ++ .../CorrectionMapsHelper.cxx | 25 ++- .../CorrectionMapsHelper.h | 56 ++++-- GPU/TPCFastTransformation/TPCFastTransform.h | 77 ++++++--- .../include/GPUWorkflow/GPUWorkflowSpec.h | 2 + .../include/GPUWorkflow/O2GPUDPLDisplay.h | 1 + GPU/Workflow/src/GPUWorkflowSpec.cxx | 3 +- GPU/Workflow/src/GPUWorkflowTPC.cxx | 12 ++ GPU/Workflow/src/O2GPUDPLDisplay.cxx | 3 + GPU/Workflow/src/gpu-reco-workflow.cxx | 1 + prodtests/full-system-test/dpl-workflow.sh | 7 +- 39 files changed, 766 insertions(+), 78 deletions(-) create mode 100644 Detectors/TPC/calibration/include/TPCCalibration/TPCMShapeScaler.h create mode 100644 Detectors/TPC/calibration/src/TPCMShapeScaler.cxx diff --git a/Detectors/Align/Workflow/src/barrel-alignment-workflow.cxx b/Detectors/Align/Workflow/src/barrel-alignment-workflow.cxx index 467e0af9d6de5..fa3842eb53a49 100644 --- a/Detectors/Align/Workflow/src/barrel-alignment-workflow.cxx +++ b/Detectors/Align/Workflow/src/barrel-alignment-workflow.cxx @@ -19,6 +19,7 @@ #include "TPCReaderWorkflow/TrackReaderSpec.h" #include "TPCReaderWorkflow/ClusterReaderSpec.h" #include "TPCWorkflow/ClusterSharingMapSpec.h" +#include "TPCWorkflow/TPCScalerSpec.h" #include "TPCCalibration/CorrectionMapsLoader.h" #include "TOFWorkflowIO/ClusterReaderSpec.h" #include "TOFWorkflowIO/TOFMatchedReaderSpec.h" @@ -147,6 +148,10 @@ WorkflowSpec defineDataProcessing(ConfigContext const& configcontext) o2::conf::ConfigurableParam::writeINI("o2_barrel_alignment_configuration.ini"); } + if (sclOpt.needTPCScalersWorkflow() && !configcontext.options().get("disable-root-input")) { + specs.emplace_back(o2::tpc::getTPCScalerSpec(sclOpt.lumiType == 2, sclOpt.enableMShapeCorrection)); + } + specs.emplace_back(o2::align::getBarrelAlignmentSpec(srcMP, src, dets, skipDetClusters, enableCosmic, postprocess, useMC, sclOpt)); // RS FIXME: check which clusters are really needed if (!postprocess) { diff --git a/Detectors/GlobalTrackingWorkflow/src/cosmics-match-workflow.cxx b/Detectors/GlobalTrackingWorkflow/src/cosmics-match-workflow.cxx index 5adadcad3f0eb..212930fb99975 100644 --- a/Detectors/GlobalTrackingWorkflow/src/cosmics-match-workflow.cxx +++ b/Detectors/GlobalTrackingWorkflow/src/cosmics-match-workflow.cxx @@ -17,6 +17,7 @@ #include "ITSMFTWorkflow/ClusterReaderSpec.h" #include "TPCReaderWorkflow/TrackReaderSpec.h" #include "TPCReaderWorkflow/ClusterReaderSpec.h" +#include "TPCWorkflow/TPCScalerSpec.h" #include "TPCWorkflow/ClusterSharingMapSpec.h" #include "TOFWorkflowIO/ClusterReaderSpec.h" #include "TOFWorkflowIO/TOFMatchedReaderSpec.h" @@ -101,6 +102,9 @@ WorkflowSpec defineDataProcessing(ConfigContext const& configcontext) } GID::mask_t srcCl = src; GID::mask_t dummy; + if (sclOpt.needTPCScalersWorkflow() && !configcontext.options().get("disable-root-input")) { + specs.emplace_back(o2::tpc::getTPCScalerSpec(sclOpt.lumiType == 2, sclOpt.enableMShapeCorrection)); + } specs.emplace_back(o2::globaltracking::getCosmicsMatchingSpec(src, useMC, sclOpt)); o2::globaltracking::InputHelper::addInputSpecs(configcontext, specs, src, src, src, useMC, dummy); // clusters MC is not needed diff --git a/Detectors/GlobalTrackingWorkflow/src/secondary-vertexing-workflow.cxx b/Detectors/GlobalTrackingWorkflow/src/secondary-vertexing-workflow.cxx index 2fe0321142642..13f2050577e50 100644 --- a/Detectors/GlobalTrackingWorkflow/src/secondary-vertexing-workflow.cxx +++ b/Detectors/GlobalTrackingWorkflow/src/secondary-vertexing-workflow.cxx @@ -17,6 +17,7 @@ #include "GlobalTrackingWorkflowHelpers/InputHelper.h" #include "ITSWorkflow/TrackReaderSpec.h" #include "TPCReaderWorkflow/TrackReaderSpec.h" +#include "TPCWorkflow/TPCScalerSpec.h" #include "TOFWorkflowIO/TOFMatchedReaderSpec.h" #include "TOFWorkflowIO/ClusterReaderSpec.h" #include "ReconstructionDataFormats/GlobalTrackID.h" @@ -96,7 +97,9 @@ WorkflowSpec defineDataProcessing(ConfigContext const& configcontext) src = src | GID::getSourcesMask("CTP"); } WorkflowSpec specs; - + if (sclOpt.needTPCScalersWorkflow() && !configcontext.options().get("disable-root-input")) { + specs.emplace_back(o2::tpc::getTPCScalerSpec(sclOpt.lumiType == 2, sclOpt.enableMShapeCorrection)); + } specs.emplace_back(o2::vertexing::getSecondaryVertexingSpec(src, enableCasc, enable3body, enableStrTr, enableCCDBParams, useMC, sclOpt)); // only TOF clusters are needed if TOF is involved, no clusters MC needed diff --git a/Detectors/GlobalTrackingWorkflow/src/tof-matcher-workflow.cxx b/Detectors/GlobalTrackingWorkflow/src/tof-matcher-workflow.cxx index c63e35deda2ed..f80d21e777815 100644 --- a/Detectors/GlobalTrackingWorkflow/src/tof-matcher-workflow.cxx +++ b/Detectors/GlobalTrackingWorkflow/src/tof-matcher-workflow.cxx @@ -33,6 +33,7 @@ #include "TSystem.h" #include "DetectorsBase/DPLWorkflowUtils.h" #include "TPCCalibration/CorrectionMapsLoader.h" +#include "TPCWorkflow/TPCScalerSpec.h" using namespace o2::framework; using DetID = o2::detectors::DetID; @@ -167,7 +168,9 @@ WorkflowSpec defineDataProcessing(ConfigContext const& configcontext) specs.push_back(s); } } - + if (sclOpt.needTPCScalersWorkflow() && !configcontext.options().get("disable-root-input")) { + specs.emplace_back(o2::tpc::getTPCScalerSpec(sclOpt.lumiType == 2, sclOpt.enableMShapeCorrection)); + } specs.emplace_back(o2::globaltracking::getTOFMatcherSpec(src, useMC, useFIT, refitTPCTOF, strict, extratolerancetrd, writeMatchable, sclOpt, nLanes)); // doTPCrefit not yet supported (need to load TPC clusters?) if (!disableRootOut) { diff --git a/Detectors/GlobalTrackingWorkflow/src/tpcits-match-workflow.cxx b/Detectors/GlobalTrackingWorkflow/src/tpcits-match-workflow.cxx index 4548e6988fac3..364bfbd220aa2 100644 --- a/Detectors/GlobalTrackingWorkflow/src/tpcits-match-workflow.cxx +++ b/Detectors/GlobalTrackingWorkflow/src/tpcits-match-workflow.cxx @@ -10,6 +10,7 @@ // or submit itself to any jurisdiction. #include "TPCWorkflow/ClusterSharingMapSpec.h" +#include "TPCWorkflow/TPCScalerSpec.h" #include "GlobalTrackingWorkflow/TPCITSMatchingSpec.h" #include "GlobalTrackingWorkflow/TrackWriterTPCITSSpec.h" #include "GlobalTrackingWorkflowHelpers/InputHelper.h" @@ -87,6 +88,9 @@ WorkflowSpec defineDataProcessing(o2::framework::ConfigContext const& configcont } o2::framework::WorkflowSpec specs; + if (sclOpt.needTPCScalersWorkflow() && !configcontext.options().get("disable-root-input")) { + specs.emplace_back(o2::tpc::getTPCScalerSpec(sclOpt.lumiType == 2, sclOpt.enableMShapeCorrection)); + } specs.emplace_back(o2::globaltracking::getTPCITSMatchingSpec(srcL, useFT0, calib, !GID::includesSource(GID::TPC, src), useMC, sclOpt)); if (!configcontext.options().get("disable-root-output")) { diff --git a/Detectors/GlobalTrackingWorkflow/study/CMakeLists.txt b/Detectors/GlobalTrackingWorkflow/study/CMakeLists.txt index 21142374ba25f..6e36c5ba6e6eb 100644 --- a/Detectors/GlobalTrackingWorkflow/study/CMakeLists.txt +++ b/Detectors/GlobalTrackingWorkflow/study/CMakeLists.txt @@ -22,6 +22,7 @@ o2_add_library(GlobalTrackingStudy O2::GlobalTrackingWorkflowHelpers O2::DataFormatsGlobalTracking O2::DetectorsVertexing + O2::TPCWorkflow O2::SimulationDataFormat) o2_add_executable(study-workflow diff --git a/Detectors/GlobalTrackingWorkflow/study/src/tpc-track-study-workflow.cxx b/Detectors/GlobalTrackingWorkflow/study/src/tpc-track-study-workflow.cxx index 10ee68d39e5e1..32418df76476b 100644 --- a/Detectors/GlobalTrackingWorkflow/study/src/tpc-track-study-workflow.cxx +++ b/Detectors/GlobalTrackingWorkflow/study/src/tpc-track-study-workflow.cxx @@ -21,6 +21,7 @@ #include "GlobalTrackingWorkflowHelpers/InputHelper.h" #include "DetectorsRaw/HBFUtilsInitializer.h" #include "TPCCalibration/CorrectionMapsLoader.h" +#include "TPCWorkflow/TPCScalerSpec.h" using namespace o2::framework; using GID = o2::dataformats::GlobalTrackID; @@ -70,6 +71,9 @@ WorkflowSpec defineDataProcessing(ConfigContext const& configcontext) } o2::globaltracking::InputHelper::addInputSpecs(configcontext, specs, srcCls, srcTrc, srcTrc, useMC); o2::globaltracking::InputHelper::addInputSpecsPVertex(configcontext, specs, useMC); // P-vertex is always needed + if (sclOpt.needTPCScalersWorkflow() && !configcontext.options().get("disable-root-input")) { + specs.emplace_back(o2::tpc::getTPCScalerSpec(sclOpt.lumiType == 2, sclOpt.enableMShapeCorrection)); + } specs.emplace_back(o2::trackstudy::getTPCTrackStudySpec(srcTrc, srcCls, useMC, sclOpt)); // configure dpl timer to inject correct firstTForbit: start from the 1st orbit of TF containing 1st sampled orbit diff --git a/Detectors/TPC/base/include/TPCBase/CDBTypes.h b/Detectors/TPC/base/include/TPCBase/CDBTypes.h index 48a79b36d96b8..4c4d163676f52 100644 --- a/Detectors/TPC/base/include/TPCBase/CDBTypes.h +++ b/Detectors/TPC/base/include/TPCBase/CDBTypes.h @@ -79,6 +79,7 @@ enum class CDBType { CalTimeSeries, ///< integrated DCAs for longer time interval CalScaler, ///< Scaler from IDCs or combined estimator CalScalerWeights, ///< Weights for scalers + CalMShape, ///< calibration object for M-shape distortions /// CorrMapParam, ///< parameters for CorrectionMapsLoader configuration /// @@ -144,6 +145,7 @@ const std::unordered_map CDBTypeMap{ {CDBType::CalTimeSeries, "TPC/Calib/TimeSeries"}, {CDBType::CalScaler, "TPC/Calib/Scaler"}, {CDBType::CalScalerWeights, "TPC/Calib/ScalerWeights"}, + {CDBType::CalMShape, "TPC/Calib/MShapeScaler"}, // 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 20e72b28d62cc..e2d19c57573f1 100644 --- a/Detectors/TPC/calibration/CMakeLists.txt +++ b/Detectors/TPC/calibration/CMakeLists.txt @@ -54,6 +54,7 @@ o2_add_library(TPCCalibration src/CalculatedEdx.cxx src/TPCScaler.cxx src/CorrMapParam.cxx + src/TPCMShapeScaler.cxx PUBLIC_LINK_LIBRARIES O2::DataFormatsTPC O2::TPCBase O2::TPCReconstruction ROOT::Minuit Microsoft.GSL::GSL @@ -107,7 +108,8 @@ o2_target_root_dictionary(TPCCalibration include/TPCCalibration/TPCFastSpaceChargeCorrectionHelper.h include/TPCCalibration/CalculatedEdx.h include/TPCCalibration/TPCScaler.h - include/TPCCalibration/CorrMapParam.h) + include/TPCCalibration/CorrMapParam.h + include/TPCCalibration/TPCMShapeScaler.h) o2_add_test_root_macro(macro/comparePedestalsAndNoise.C PUBLIC_LINK_LIBRARIES O2::TPCBase diff --git a/Detectors/TPC/calibration/include/TPCCalibration/CorrectionMapsLoader.h b/Detectors/TPC/calibration/include/TPCCalibration/CorrectionMapsLoader.h index d72f69eeca422..6e06dfad09b34 100644 --- a/Detectors/TPC/calibration/include/TPCCalibration/CorrectionMapsLoader.h +++ b/Detectors/TPC/calibration/include/TPCCalibration/CorrectionMapsLoader.h @@ -38,8 +38,14 @@ namespace tpc { struct CorrectionMapsLoaderGloOpts { - int lumiType = 0; - int lumiMode = 0; + int lumiType = 0; ///< 0: no scaling, 1: CTP, 2: IDC + int lumiMode = 0; ///< 0: classical scaling, 1: Using of the derivative map, 2: Using of the derivative map for MC + bool enableMShapeCorrection = false; + + bool needTPCScalersWorkflow() const + { + return lumiType == 2 || enableMShapeCorrection; + } }; class CorrectionMapsLoader : public o2::gpu::CorrectionMapsHelper @@ -55,6 +61,7 @@ class CorrectionMapsLoader : public o2::gpu::CorrectionMapsHelper void updateVDrift(float vdriftCorr, float vdrifRef, float driftTimeOffset = 0); void init(o2::framework::InitContext& ic); void copySettings(const CorrectionMapsLoader& src); + void updateInverse(); /// recalculate inverse correction static void requestCCDBInputs(std::vector& inputs, std::vector& options, const CorrectionMapsLoaderGloOpts& gloOpts); static void addGlobalOptions(std::vector& options); @@ -65,8 +72,9 @@ class CorrectionMapsLoader : public o2::gpu::CorrectionMapsHelper static void addOption(std::vector& options, o2::framework::ConfigParamSpec&& osp); static void addInput(std::vector& inputs, o2::framework::InputSpec&& isp); - float mInstLumiFactor = 1.0; // multiplicative factor for inst. lumi - int mCTPLumiSource = 0; // 0: main, 1: alternative CTP lumi source + float mInstLumiFactor = 1.0; // multiplicative factor for inst. lumi + int mCTPLumiSource = 0; // 0: main, 1: alternative CTP lumi source + int mNthreadsInv = 1; // number of threads used for calculating the inverse correction #endif }; diff --git a/Detectors/TPC/calibration/include/TPCCalibration/TPCFastSpaceChargeCorrectionHelper.h b/Detectors/TPC/calibration/include/TPCCalibration/TPCFastSpaceChargeCorrectionHelper.h index 2f4291a54a086..6aefd7ffa6b6a 100644 --- a/Detectors/TPC/calibration/include/TPCCalibration/TPCFastSpaceChargeCorrectionHelper.h +++ b/Detectors/TPC/calibration/include/TPCCalibration/TPCFastSpaceChargeCorrectionHelper.h @@ -97,6 +97,12 @@ class TPCFastSpaceChargeCorrectionHelper void testGeometry(const TPCFastTransformGeo& geo) const; + /// initialise inverse transformation + void initInverse(o2::gpu::TPCFastSpaceChargeCorrection& correction, bool prn); + + /// initialise inverse transformation from linear combination of several input corrections + void initInverse(std::vector& corrections, const std::vector scaling, bool prn); + private: /// geometry initialization void initGeometry(); @@ -107,9 +113,6 @@ class TPCFastSpaceChargeCorrectionHelper /// initialise max drift length void initMaxDriftLength(o2::gpu::TPCFastSpaceChargeCorrection& correction, bool prn); - /// initialise inverse transformation - void initInverse(o2::gpu::TPCFastSpaceChargeCorrection& correction, bool prn); - static TPCFastSpaceChargeCorrectionHelper* sInstance; ///< singleton instance bool mIsInitialized = 0; ///< initialization flag int mNthreads = 1; ///< n of threads to use diff --git a/Detectors/TPC/calibration/include/TPCCalibration/TPCMShapeScaler.h b/Detectors/TPC/calibration/include/TPCCalibration/TPCMShapeScaler.h new file mode 100644 index 0000000000000..ad28197218726 --- /dev/null +++ b/Detectors/TPC/calibration/include/TPCCalibration/TPCMShapeScaler.h @@ -0,0 +1,99 @@ +// 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 TPCMShapeScaler.h +/// \author Matthias Kleiner + +#ifndef ALICEO2_TPC_TPCMSHAPESCALER +#define ALICEO2_TPC_TPCMSHAPESCALER + +#include "DataFormatsTPC/Defs.h" +#include + +class TTree; + +namespace o2::tpc +{ + +/* +Class for storing the scalers (obtained from DCAs) which are used to correct for the M-shape distortions +*/ + +struct TPCMShape { + double getEndTime(float sampling) const { return mStartTimeMS + mScalerA.size() * sampling; } + double mStartTimeMS; ///< start time of the M-shape + std::vector mScalerA; ///< DCA or whatever + std::vector mZPosition; ///< z-position of the delta potential + bool operator<(const TPCMShape& r) const { return mStartTimeMS < r.mStartTimeMS; } + ClassDefNV(TPCMShape, 1); +}; + +class TPCMShapeScaler +{ + public: + /// default constructor + TPCMShapeScaler() = default; + + /// \return returns number of stored M-shapes + int getNValues() const { return mMShape.size(); } + + /// setting the scalers + void setScaler(const std::vector& time, const std::vector& scaler, const std::vector& zpos); + + /// \return returns run number for which this object is valid + void setRun(int run) { mRun = run; } + + /// dump this object to a file + /// \param file output file + void dumpToFile(const char* file, const char* name); + + /// load parameters from input file (which were written using the writeToFile method) + /// \param inpf input file + void loadFromFile(const char* inpf, const char* name); + + /// set this object from input tree + void setFromTree(TTree& tpcScalerTree); + + /// set sampling time of the stored values + float setSamplingTimeMS(float t) { return mSamplingTimeMS = t; } + + /// adding a M-shape + void addMShape(const TPCMShape& vshape) { mMShape.emplace_back(); } + + /// \return returns run number for which this object is valid + int getRun() const { return mRun; } + + /// \return returns M-shape scaling value for given timestamp + /// \param timestamp timestamp for which the scaler is queried + std::pair getScaler(const double timestamp) const; + + /// \return returns sampling time of the stored values + float getSamplingTimeMS() const { return mSamplingTimeMS; } + + /// \return returns all stored M-shapes + const auto& getMShapes() const { return mMShape; } + + private: + int mRun{}; ///< run for which this object is valid + float mSamplingTimeMS = 1; ///< sampling time of the M-shape scalers + std::vector mMShape; ///< vector containing the M-shapes + + /// if distance to neighbouring data is larger than sampling time return 0 for the scaling + bool checkDeltaTime(double deltaTime) const { return deltaTime < 1.5 * mSamplingTimeMS; } + + /// return linear interpolated value + float getVal(float x, float x0, float x1, float y0, float y1) const { return (y0 * (x1 - x) + y1 * (x - x0)) / (x1 - x0); } + + ClassDefNV(TPCMShapeScaler, 1); +}; + +} // namespace o2::tpc +#endif diff --git a/Detectors/TPC/calibration/include/TPCCalibration/TPCScaler.h b/Detectors/TPC/calibration/include/TPCCalibration/TPCScaler.h index 532c263244394..dbe163ee4a115 100644 --- a/Detectors/TPC/calibration/include/TPCCalibration/TPCScaler.h +++ b/Detectors/TPC/calibration/include/TPCCalibration/TPCScaler.h @@ -50,8 +50,8 @@ class TPCScaler /// \return returns number of stored TPC scaler values int getNValues(o2::tpc::Side side) const { return (side == o2::tpc::Side::A ? mScalerA.size() : mScalerC.size()); } - /// set the parameters for the coefficients of the polynomial - /// \param params parameter for the coefficients + /// set TPC scalers + /// \param values scaling values void setScaler(const std::vector& values, const o2::tpc::Side side) { (side == o2::tpc::Side::A ? (mScalerA = values) : (mScalerC = values)); } /// \return returns ion drift time in ms @@ -73,6 +73,19 @@ class TPCScaler /// \param file output file void dumpToFile(const char* file, const char* name); + /// dump this object to a file for a given time range + /// \param file output file + /// \param startTime start time stamp in ms of the dumped object + /// \param endTime end time stamp in ms of the dumped object (if endTimeMS=-1 the end of the stored data is used) + void dumpToFile(const char* file, const char* name, double startTimeMS, double endTimeMS); + + /// dump this object to several file each containing the scalers for n minutes + /// \param file output file + /// \param minutesPerObject minutes each output file is covering + /// \param marginMS additional margin to the beginning of each object (should be >= getIonDriftTimeMS()) + /// \param marginCCDBMinutes margin for CCDB timestamp + void dumpToFileSlices(const char* file, const char* name, int minutesPerObject, float marginMS = 500, float marginCCDBMinutes = 1); + /// load parameters from input file (which were written using the writeToFile method) /// \param inpf input file void loadFromFile(const char* inpf, const char* name); @@ -114,6 +127,11 @@ class TPCScaler /// \return returns duration in ms for which the scalers are defined float getDurationMS(o2::tpc::Side side) const { return mIntegrationTimeMS * getNValues(side); } + /// clamp scaler values + /// \param minThreshold minimum accepted scaler value + /// \param maxThreshold maximum accepted scaler value + void clampScalers(float minThreshold, float maxThreshold, Side side); + /// setting the weights for the scalers void setScalerWeights(const TPCScalerWeights& weights) { mScalerWeights = weights; } diff --git a/Detectors/TPC/calibration/src/CorrectionMapsLoader.cxx b/Detectors/TPC/calibration/src/CorrectionMapsLoader.cxx index beb0fd72af601..71cdde756ae1a 100644 --- a/Detectors/TPC/calibration/src/CorrectionMapsLoader.cxx +++ b/Detectors/TPC/calibration/src/CorrectionMapsLoader.cxx @@ -23,6 +23,9 @@ #include "Framework/DeviceSpec.h" #include "Framework/ConfigParamRegistry.h" #include "DataFormatsCTP/LumiInfo.h" +#include "TPCCalibration/TPCFastSpaceChargeCorrectionHelper.h" +#include "TPCSpaceCharge/SpaceCharge.h" +#include "CommonDataFormat/Pair.h" using namespace o2::tpc; using namespace o2::framework; @@ -36,6 +39,9 @@ void CorrectionMapsLoader::updateVDrift(float vdriftCorr, float vdrifRef, float if (mCorrMapRef) { o2::tpc::TPCFastTransformHelperO2::instance()->updateCalibration(*mCorrMapRef, 0, vdriftCorr, vdrifRef, driftTimeOffset); } + if (mCorrMapMShape) { + o2::tpc::TPCFastTransformHelperO2::instance()->updateCalibration(*mCorrMapMShape, 0, vdriftCorr, vdrifRef, driftTimeOffset); + } } //________________________________________________________ @@ -57,10 +63,53 @@ void CorrectionMapsLoader::extractCCDBInputs(ProcessingContext& pc) } lumiObj = lumiPrev; } - setInstLumi(mInstLumiFactor * (mCTPLumiSource == 0 ? lumiObj.getLumi() : lumiObj.getLumiAlt())); + setInstLumi(mInstLumiFactor * (mCTPLumiSource == 0 ? lumiObj.getLumi() : lumiObj.getLumiAlt()), !getUseMShapeCorrection()); } else if (getLumiScaleType() == 2 && mInstLumiOverride <= 0.) { float tpcScaler = pc.inputs().get("tpcscaler"); - setInstLumi(mInstLumiFactor * tpcScaler); + setInstLumi(mInstLumiFactor * tpcScaler, !getUseMShapeCorrection()); + } + if (getUseMShapeCorrection()) { + auto mshapescaler = pc.inputs().get>("mshape"); + float potential = mshapescaler.first; + float deltaPotPos = mshapescaler.second; + setMShapeScaler(potential); + if (potential != 0) { + LOGP(info, "Recalculating M-shape correction map for z position of delta potential {}", deltaPotPos); + // update analytical map! + const int nthreads = 5; + using namespace o2::tpc; + o2::tpc::SpaceCharge::setNThreads(nthreads); + const Side side = Side::A; + int field = 5; + o2::tpc::SpaceCharge sc = o2::tpc::SpaceCharge(field, 257, 257, 45); + + float deltaPot = 100; + sc.setIFCChargeUpFallingPot(deltaPot, deltaPotPos, 0, 250, 0, side); // 1/x + sc.setIFCChargeUpRisingPot(deltaPot, deltaPotPos, 0, 0, 0, side); + + sc.poissonSolver(side); + sc.calcEField(side); + const auto numEFields = sc.getElectricFieldsInterpolator(side); + sc.calcGlobalCorrections(numEFields); + + std::function getCorrections = [&sc = sc](const int roc, double x, double y, double z, double& dx, double& dy, double& dz) { + Side side = roc < 18 ? Side::A : Side::C; + sc.getCorrections(x, y, z, side, dx, dy, dz); + }; + + TPCFastSpaceChargeCorrectionHelper::instance()->setNthreads(nthreads); + const int nKnotsY = 4; + const int nKnotsZ = 4; + std::unique_ptr spCorrection = TPCFastSpaceChargeCorrectionHelper::instance()->createFromGlobalCorrection(getCorrections, nKnotsY, nKnotsZ); + std::unique_ptr fastTransform(TPCFastTransformHelperO2::instance()->create(0, *spCorrection)); + LOGP(info, "Setting M-shape correction from calculated model"); + setCorrMapMShape(std::move(fastTransform)); + setUpdatedMapMShape(); + } + } + // update inverse in case it is requested + if (!mScaleInverse) { + updateInverse(); } } @@ -87,8 +136,11 @@ void CorrectionMapsLoader::requestCCDBInputs(std::vector& inputs, std addInput(inputs, {"tpcscaler", o2::header::gDataOriginTPC, "TPCSCALER", 0, Lifetime::Timeframe}); } - addInput(inputs, {"tpcCorrPar", "TPC", "CorrMapParam", 0, Lifetime::Condition, ccdbParamSpec(CDBTypeMap.at(CDBType::CorrMapParam), {}, 0)}); // load once + if (gloOpts.enableMShapeCorrection) { + addInput(inputs, {"mshape", o2::header::gDataOriginTPC, "TPCMShapeScaler", 0, Lifetime::Timeframe}); + } + addInput(inputs, {"tpcCorrPar", "TPC", "CorrMapParam", 0, Lifetime::Condition, ccdbParamSpec(CDBTypeMap.at(CDBType::CorrMapParam), {}, 0)}); // load once addOptions(options); } @@ -97,6 +149,8 @@ void CorrectionMapsLoader::addOptions(std::vector& options) { // these are options which should be added at the level of device using TPC corrections // At the moment - nothing, all options are moved to configurable param CorrMapParam + addOption(options, ConfigParamSpec{"recalculate-inverse-correction", o2::framework::VariantType::Bool, false, {"Recalculate the inverse correction in case linear scaling of corrections is used"}}); + addOption(options, ConfigParamSpec{"nthreads-inverse-correction", o2::framework::VariantType::Int, -1, {"Number of threads used for calculating the inverse correction (-1=all threads)"}}); } //________________________________________________________ @@ -105,6 +159,7 @@ void CorrectionMapsLoader::addGlobalOptions(std::vector& option // these are options which should be added at the workflow level, since they modify the inputs of the devices addOption(options, ConfigParamSpec{"lumi-type", o2::framework::VariantType::Int, 0, {"1 = require CTP lumi for TPC correction scaling, 2 = require TPC scalers for TPC correction scaling"}}); addOption(options, ConfigParamSpec{"corrmap-lumi-mode", o2::framework::VariantType::Int, 0, {"scaling mode: (default) 0 = static + scale * full; 1 = full + scale * derivative; 2 = full + scale * derivative (for MC)"}}); + addOption(options, ConfigParamSpec{"enable-M-shape-correction", o2::framework::VariantType::Bool, false, {"Enable M-shape distortion correction"}}); } //________________________________________________________ @@ -113,6 +168,7 @@ CorrectionMapsLoaderGloOpts CorrectionMapsLoader::parseGlobalOptions(const o2::f CorrectionMapsLoaderGloOpts tpcopt; tpcopt.lumiType = opts.get("lumi-type"); tpcopt.lumiMode = opts.get("corrmap-lumi-mode"); + tpcopt.enableMShapeCorrection = opts.get("enable-M-shape-correction"); return tpcopt; } @@ -202,8 +258,16 @@ void CorrectionMapsLoader::init(o2::framework::InitContext& ic) LOGP(fatal, "Lumi scaling source TPCScaler is not compatible with TPC correction lumi scaler type {}", getLumiScaleType()); } break; + setLumiScaleType(2); + } else if (route.matcher == InputSpec{"mshape", o2::header::gDataOriginTPC, "TPCMShapeScaler", 0, Lifetime::Timeframe}) { + enableMShapeCorrection(true); } } + mScaleInverse = !(ic.options().get("recalculate-inverse-correction")); + mNthreadsInv = (ic.options().get("nthreads-inverse-correction")); + if (!mScaleInverse) { + LOGP(info, "Recalculating the inverse correction for every TF"); + } } //________________________________________________________ @@ -217,8 +281,26 @@ void CorrectionMapsLoader::copySettings(const CorrectionMapsLoader& src) setMeanLumiRefOverride(src.getMeanLumiRefOverride()); setInstLumiOverride(src.getInstLumiOverride()); setLumiScaleMode(src.getLumiScaleMode()); + enableMShapeCorrection(src.getUseMShapeCorrection()); + setMShapeScaler(src.getMShapeScaler(), false); mInstLumiFactor = src.mInstLumiFactor; mCTPLumiSource = src.mCTPLumiSource; + mLumiScaleMode = src.mLumiScaleMode; + mScaleInverse = src.getScaleInverse(); + mNthreadsInv = src.mNthreadsInv; +} + +void CorrectionMapsLoader::updateInverse() +{ + (mNthreadsInv < 0) ? TPCFastSpaceChargeCorrectionHelper::instance()->setNthreadsToMaximum() : TPCFastSpaceChargeCorrectionHelper::instance()->setNthreads(mNthreadsInv); + if (mLumiScaleMode == 1 || mLumiScaleMode == 2) { + setUpdatedMap(); + std::vector scaling{1, mLumiScale, mMShapeScaler}; + std::vector corr{&(mCorrMap->getCorrection()), &(mCorrMapRef->getCorrection()), &(mCorrMapMShape->getCorrection())}; + TPCFastSpaceChargeCorrectionHelper::instance()->initInverse(corr, scaling, false); + } else { + LOGP(info, "Reinitializing inverse correction with lumi scale mode {} not supported for now", mLumiScaleMode); + } } #endif // #ifndef GPUCA_GPUCODE_DEVICE diff --git a/Detectors/TPC/calibration/src/TPCCalibrationLinkDef.h b/Detectors/TPC/calibration/src/TPCCalibrationLinkDef.h index 9b0fd7ce71e0e..13fd86267b841 100644 --- a/Detectors/TPC/calibration/src/TPCCalibrationLinkDef.h +++ b/Detectors/TPC/calibration/src/TPCCalibrationLinkDef.h @@ -114,4 +114,6 @@ #pragma link C++ class o2::tpc::CalculatedEdx + ; #pragma link C++ class o2::tpc::TPCScaler + ; #pragma link C++ struct o2::tpc::TPCScalerWeights + ; +#pragma link C++ class o2::tpc::TPCMShapeScaler + ; +#pragma link C++ struct o2::tpc::TPCMShape + ; #endif diff --git a/Detectors/TPC/calibration/src/TPCFastSpaceChargeCorrectionHelper.cxx b/Detectors/TPC/calibration/src/TPCFastSpaceChargeCorrectionHelper.cxx index 0d20928eef451..980f1bda2aad1 100644 --- a/Detectors/TPC/calibration/src/TPCFastSpaceChargeCorrectionHelper.cxx +++ b/Detectors/TPC/calibration/src/TPCFastSpaceChargeCorrectionHelper.cxx @@ -28,6 +28,7 @@ #include "Riostream.h" #include #include +#include "TStopwatch.h" using namespace o2::gpu; @@ -653,10 +654,22 @@ void TPCFastSpaceChargeCorrectionHelper::initMaxDriftLength(o2::gpu::TPCFastSpac void TPCFastSpaceChargeCorrectionHelper::initInverse(o2::gpu::TPCFastSpaceChargeCorrection& correction, bool prn) { - /// initialise inverse transformation + std::vector corr{&correction}; + initInverse(corr, std::vector{1}, prn); +} +void TPCFastSpaceChargeCorrectionHelper::initInverse(std::vector& corrections, const std::vector scaling, bool prn) +{ + /// initialise inverse transformation + TStopwatch watch; LOG(info) << "fast space charge correction helper: init inverse"; + if (corrections.size() != scaling.size()) { + LOGP(error, "Input corrections and scaling values have different size"); + return; + } + + auto& correction = *(corrections.front()); initMaxDriftLength(correction, prn); double tpcR2min = mGeo.getRowInfo(0).x - 1.; @@ -709,6 +722,17 @@ void TPCFastSpaceChargeCorrectionHelper::initInverse(o2::gpu::TPCFastSpaceCharge for (double v = v0; v < v1 + stepV; v += stepV) { float dx, du, dv; correction.getCorrection(slice, row, u, v, dx, du, dv); + dx *= scaling[0]; + du *= scaling[0]; + dv *= scaling[0]; + // add remaining corrections + for (int i = 1; i < corrections.size(); ++i) { + float dxTmp, duTmp, dvTmp; + corrections[i]->getCorrection(slice, row, u, v, dxTmp, duTmp, dvTmp); + dx += dxTmp * scaling[i]; + du += duTmp * scaling[i]; + dv += dvTmp * scaling[i]; + } double cx = x + dx; double cu = u + du; double cv = v + dv; @@ -789,6 +813,8 @@ void TPCFastSpaceChargeCorrectionHelper::initInverse(o2::gpu::TPCFastSpaceCharge } } // slice + float duration = watch.RealTime(); + LOGP(info, "Inverse took: {}s", duration); } } // namespace tpc diff --git a/Detectors/TPC/calibration/src/TPCMShapeScaler.cxx b/Detectors/TPC/calibration/src/TPCMShapeScaler.cxx new file mode 100644 index 0000000000000..0943834eb0c50 --- /dev/null +++ b/Detectors/TPC/calibration/src/TPCMShapeScaler.cxx @@ -0,0 +1,161 @@ +// 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 TPCMShapeScaler.cxx +/// \brief Definition of TPCMShapeScaler class +/// +/// \author Matthias Kleiner + +#include "TPCCalibration/TPCMShapeScaler.h" +#include +#include +#include "Framework/Logger.h" +#include "CommonConstants/LHCConstants.h" + +using namespace o2::tpc; + +void TPCMShapeScaler::dumpToFile(const char* file, const char* name) +{ + TFile out(file, "RECREATE"); + TTree tree(name, name); + tree.SetAutoSave(0); + tree.Branch("TPCMShapeScaler", this); + tree.Fill(); + tree.SetAlias("relTime", "(Iteration$ % mMShape.@mScalerA.size()) * mSamplingTimeMS"); + tree.SetAlias("time", "mMShape.mStartTimeMS + (Iteration$ % mMShape.@mScalerA.size()) * mSamplingTimeMS"); + out.WriteObject(&tree, name); +} + +void TPCMShapeScaler::loadFromFile(const char* inpf, const char* name) +{ + TFile out(inpf, "READ"); + TTree* tree = (TTree*)out.Get(name); + setFromTree(*tree); +} + +void TPCMShapeScaler::setFromTree(TTree& tpcScalerTree) +{ + TPCMShapeScaler* scalerTmp = this; + tpcScalerTree.SetBranchAddress("TPCMShapeScaler", &scalerTmp); + const int entries = tpcScalerTree.GetEntries(); + if (entries > 0) { + tpcScalerTree.GetEntry(0); + } else { + LOGP(error, "TPCMShapeScaler not found in input file"); + } + tpcScalerTree.SetBranchAddress("TPCMShapeScaler", nullptr); +} + +std::pair TPCMShapeScaler::getScaler(const double timestamp) const +{ + auto idx = std::distance(mMShape.begin(), std::upper_bound(mMShape.begin(), mMShape.end(), TPCMShape{timestamp})); + // LOGP(info, "idx: {} timestamp: {} first: {}", idx, timestamp, mMShape.front().mStartTimeMS); + + if (idx == 0) { + return {0.f, 0.f}; + } else { + const int idxV = idx - 1; + const auto& vs = mMShape[idxV]; + const double timeEndMShape = vs.getEndTime(mSamplingTimeMS); + if (timestamp > timeEndMShape) { + return {0.f, 0.f}; + } + + // perform linear interpolation + float relTime = timestamp - vs.mStartTimeMS; + + // get value to the left and right + float idxData = relTime / mSamplingTimeMS; + + // LOGP(info, "idx: {} timestamp: {} vs.mStartTimeMS: {} idxData: {}", idx, timestamp, vs.mStartTimeMS, idxData); + + // if the time is very close to the sampled value -> just return the sampled value + const float epsilon = 0.001; + const int nearestIdx = static_cast(relTime / mSamplingTimeMS + 0.5); + if (std::abs(idxData - nearestIdx) < epsilon) { + // LOGP(info, "return nearest value for relTime {} mSamplingTimeMS {} idxData: {} nearestIdx: {}", relTime, mSamplingTimeMS, idxData, nearestIdx); + return {vs.mScalerA[nearestIdx], vs.mZPosition[nearestIdx]}; + } + + const int idxLow = static_cast(idxData); + const int idxUp = idxLow + 1; + float scaler = getVal(idxData, idxLow, idxUp, vs.mScalerA[idxLow], vs.mScalerA[idxUp]); + float zPos = getVal(idxData, idxLow, idxUp, vs.mZPosition[idxLow], vs.mZPosition[idxUp]); + // LOGP(info, "Using idxData: {} idxLow: {} idxUp: {} Got: scaler {} zPos: {} vs.mScalerA[idxLow]: {} vs.mScalerA[idxUp]: {} vs.mZPosition[idxLow]: {} vs.mZPosition[idxUp]: {} vs.mScalerA.size() {}", idxData, idxLow, idxUp, scaler, zPos, vs.mScalerA[idxLow], vs.mScalerA[idxUp], vs.mZPosition[idxLow], vs.mZPosition[idxUp], vs.mScalerA.size()); + return {scaler, zPos}; + } + + // o2::utils::DebugStreamer::instance()->getStreamer("debug_fasttransform", "UPDATE") << o2::utils::DebugStreamer::instance()->getUniqueTreeName("tree_Transform").data() + // // corrections in x, u, v + // << "dxOrig=" << dxOrig + // << "\n"; + + // GPUCA_DEBUG_STREAMER_CHECK(if (o2::utils::DebugStreamer::checkStream(o2::utils::StreamFlags::streamFastTransform)) { + // + // } + + // mMShape[idx] = getEndTime(float sampling) + // if (idx == mScalerA.size()) { + // // check end range + // const double deltaTime = std::abs(mScalerATime.back() - timestamp); + // if (checkDeltaTime(deltaTime)) { + // return {mScalerA.back(), mZPosition.back()}; + // } + // } else if (idx == 0) { + // const double deltaTime = std::abs(mScalerATime.front() - timestamp); + // if (checkDeltaTime(deltaTime)) { + // return {mScalerA.front(), mZPosition.front()}; + // } + // } else { + // // check if upper or lower value is closer + // const double deltaTimeLow = std::abs(mScalerATime[idx - 1] - timestamp); + // const double deltaTimeUp = std::abs(mScalerATime[idx] - timestamp); + + // // return closes value + // if (checkDeltaTime(deltaTimeLow) && checkDeltaTime(deltaTimeUp)) { + // // do linear interpolation + // float scaler = getVal(deltaTimeLow, 0, deltaTimeLow + deltaTimeUp, mScalerA[idx - 1], mScalerA[idx]); + // float zPos = getVal(deltaTimeLow, 0, deltaTimeLow + deltaTimeUp, mZPosition[idx - 1], mZPosition[idx]); + // return {scaler, zPos}; + // } else if (checkDeltaTime(deltaTimeLow)) { + // return {mScalerA[idx - 1], mZPosition[idx - 1]}; + // } else if (checkDeltaTime(deltaTimeUp)) { + // return {mScalerA[idx], mZPosition[idx]}; + // } + // } + + return {0.f, 0.f}; +} + +void TPCMShapeScaler::setScaler(const std::vector& time, const std::vector& scaler, const std::vector& zpos) +{ + TPCMShape mshape; + mSamplingTimeMS = time[1] - time[0]; + + int idxStart = 0; + for (int i = 1; i < time.size(); ++i) { + const double deltaT = time[i] - time[i - 1]; + if (deltaT > 1.1 * mSamplingTimeMS) { + // next v-shape + mshape.mZPosition = std::vector(zpos.begin() + idxStart, zpos.begin() + i); + mshape.mScalerA = std::vector(scaler.begin() + idxStart, scaler.begin() + i); + mshape.mStartTimeMS = time[idxStart]; + mMShape.emplace_back(mshape); + idxStart = i; + } else if (i == time.size() - 1) { + // last v-shape + mshape.mZPosition = std::vector(zpos.begin() + idxStart, zpos.end()); + mshape.mScalerA = std::vector(scaler.begin() + idxStart, scaler.end()); + mshape.mStartTimeMS = time[idxStart]; + mMShape.emplace_back(mshape); + } + } +} \ No newline at end of file diff --git a/Detectors/TPC/calibration/src/TPCScaler.cxx b/Detectors/TPC/calibration/src/TPCScaler.cxx index f10555892bb33..d387f1d0b74af 100644 --- a/Detectors/TPC/calibration/src/TPCScaler.cxx +++ b/Detectors/TPC/calibration/src/TPCScaler.cxx @@ -18,6 +18,7 @@ #include #include #include "Framework/Logger.h" +#include "CommonConstants/LHCConstants.h" using namespace o2::tpc; @@ -31,6 +32,66 @@ void TPCScaler::dumpToFile(const char* file, const char* name) out.WriteObject(&tree, name); } +void TPCScaler::dumpToFile(const char* file, const char* name, double startTimeMS, double endTimeMS) +{ + TPCScaler scaler; + scaler.mIonDriftTimeMS = mIonDriftTimeMS; + scaler.mRun = mRun; + scaler.mFirstTFOrbit = (startTimeMS <= 0) ? mFirstTFOrbit : (mFirstTFOrbit + (startTimeMS - mTimeStampMS) / (o2::constants::lhc::LHCOrbitMUS * 0.001)); + scaler.mTimeStampMS = (startTimeMS <= 0) ? mTimeStampMS : startTimeMS; + scaler.mIntegrationTimeMS = mIntegrationTimeMS; + + const int dataIdx = (startTimeMS <= 0) ? 0 : getDataIdx(startTimeMS); + const int idxDataStart = (startTimeMS <= 0) ? 0 : dataIdx; + const int idxDataEndA = ((endTimeMS < 0) || (dataIdx > getNValues(Side::A))) ? getNValues(Side::A) : dataIdx; + const int idxDataEndC = ((endTimeMS < 0) || (dataIdx > getNValues(Side::C))) ? getNValues(Side::C) : dataIdx; + scaler.mScalerA = std::vector(mScalerA.begin() + idxDataStart, mScalerA.begin() + idxDataEndA); + scaler.mScalerC = std::vector(mScalerC.begin() + idxDataStart, mScalerC.begin() + idxDataEndC); + scaler.dumpToFile(file, name); +} + +void TPCScaler::dumpToFileSlices(const char* file, const char* name, int minutesPerObject, float marginMS, float marginCCDBMinutes) +{ + if (getNValues(o2::tpc::Side::A) != getNValues(o2::tpc::Side::C)) { + LOGP(error, "Number of points stored for A-side and C-side is different"); + return; + } + const long marginCCDBMS = marginCCDBMinutes * 60 * 1000; + const float msPerObjectTmp = minutesPerObject * 60 * 1000; + const int valuesPerSlice = static_cast(msPerObjectTmp / mIntegrationTimeMS); + const int marginPerSlice = static_cast(marginMS / mIntegrationTimeMS); + int nSlices = getNValues(Side::A) / valuesPerSlice; // number of output objects + + LOGP(info, "Producing {} objects with a CCDB margin of {} ms and {} margin per slice", nSlices, marginCCDBMS, marginPerSlice); + if (nSlices == 0) { + nSlices = 1; + } + + for (int i = 0; i < nSlices; ++i) { + const int idxDataStart = (i == 0) ? 0 : (i * valuesPerSlice - marginPerSlice); + int idxDataEnd = (i == nSlices - 1) ? getNValues(Side::A) : ((i + 1) * valuesPerSlice + marginPerSlice); + if (idxDataEnd > getNValues(Side::A)) { + idxDataEnd = getNValues(Side::A); + } + + TPCScaler scaler; + scaler.mIonDriftTimeMS = mIonDriftTimeMS; + scaler.mRun = mRun; + scaler.mTimeStampMS = (i == 0) ? mTimeStampMS : (mTimeStampMS + idxDataStart * static_cast(mIntegrationTimeMS)); + scaler.mFirstTFOrbit = mFirstTFOrbit + (scaler.mTimeStampMS - mTimeStampMS) / (o2::constants::lhc::LHCOrbitMUS * 0.001); + scaler.mIntegrationTimeMS = mIntegrationTimeMS; + scaler.mScalerA = std::vector(mScalerA.begin() + idxDataStart, mScalerA.begin() + idxDataEnd); + scaler.mScalerC = std::vector(mScalerC.begin() + idxDataStart, mScalerC.begin() + idxDataEnd); + + const long timePerSliceMS = valuesPerSlice * mIntegrationTimeMS; + const long tsCCDBStart = mTimeStampMS + i * timePerSliceMS; + const long tsCCDBStartMargin = (i == 0) ? (tsCCDBStart - marginCCDBMS) : tsCCDBStart; + const long tsCCDBEnd = (i == nSlices - 1) ? (getEndTimeStampMS(o2::tpc::Side::A) + marginCCDBMS) : (tsCCDBStart + timePerSliceMS); + const std::string fileOut = fmt::format("{}_{}_{}_{}.root", file, i, tsCCDBStartMargin, tsCCDBEnd - 1); + scaler.dumpToFile(fileOut.data(), name); + } +} + void TPCScaler::loadFromFile(const char* inpf, const char* name) { TFile out(inpf, "READ"); @@ -109,3 +170,9 @@ float TPCScalerWeights::getWeight(float deltaTime) const return y; } } + +void TPCScaler::clampScalers(float minThreshold, float maxThreshold, Side side) +{ + auto& scaler = (side == o2::tpc::Side::A) ? mScalerA : mScalerC; + std::transform(std::begin(scaler), std::end(scaler), std::begin(scaler), [minThreshold, maxThreshold](auto val) { return std::clamp(val, minThreshold, maxThreshold); }); +} diff --git a/Detectors/TPC/workflow/include/TPCWorkflow/TPCScalerSpec.h b/Detectors/TPC/workflow/include/TPCWorkflow/TPCScalerSpec.h index 3cea2e0175772..b85a882870ecb 100644 --- a/Detectors/TPC/workflow/include/TPCWorkflow/TPCScalerSpec.h +++ b/Detectors/TPC/workflow/include/TPCWorkflow/TPCScalerSpec.h @@ -19,7 +19,7 @@ namespace o2 namespace tpc { -o2::framework::DataProcessorSpec getTPCScalerSpec(bool enableWeights); +o2::framework::DataProcessorSpec getTPCScalerSpec(bool enableIDCs, bool enableMShape); } // end namespace tpc } // end namespace o2 diff --git a/Detectors/TPC/workflow/src/RecoWorkflow.cxx b/Detectors/TPC/workflow/src/RecoWorkflow.cxx index fd77845ff844b..417e9a68c5b29 100644 --- a/Detectors/TPC/workflow/src/RecoWorkflow.cxx +++ b/Detectors/TPC/workflow/src/RecoWorkflow.cxx @@ -22,6 +22,7 @@ #include "TPCWorkflow/TPCTriggerWriterSpec.h" #include "TPCReaderWorkflow/PublisherSpec.h" #include "TPCWorkflow/ClustererSpec.h" +#include "TPCWorkflow/TPCScalerSpec.h" #include "TPCWorkflow/ClusterDecoderRawSpec.h" #include "GPUWorkflow/GPUWorkflowSpec.h" #include "TPCWorkflow/EntropyEncoderSpec.h" @@ -196,6 +197,9 @@ framework::WorkflowSpec getWorkflow(CompletionPolicyData* policyData, std::vecto laneConfiguration, &hook}, propagateMC)); + if (sclOpts.needTPCScalersWorkflow()) { // for standalone tpc-reco workflow + specs.emplace_back(o2::tpc::getTPCScalerSpec(sclOpts.lumiType == 2, sclOpts.enableMShapeCorrection)); + } } else if (inputType == InputType::ClustersHardware) { specs.emplace_back(o2::tpc::getPublisherSpec(PublisherConf{ "tpc-clusterhardware-reader", @@ -214,6 +218,9 @@ framework::WorkflowSpec getWorkflow(CompletionPolicyData* policyData, std::vecto if (!getenv("DPL_DISABLE_TPC_TRIGGER_READER") || atoi(getenv("DPL_DISABLE_TPC_TRIGGER_READER")) != 1) { specs.emplace_back(o2::tpc::getTPCTriggerReaderSpec()); } + if (sclOpts.needTPCScalersWorkflow()) { // for standalone tpc-reco workflow + specs.emplace_back(o2::tpc::getTPCScalerSpec(sclOpts.lumiType == 2, sclOpts.enableMShapeCorrection)); + } } else if (inputType == InputType::CompClusters) { // TODO: need to check if we want to store the MC labels alongside with compressed clusters // for the moment reading of labels is disabled (last parameter is false) @@ -446,6 +453,7 @@ framework::WorkflowSpec getWorkflow(CompletionPolicyData* policyData, std::vecto cfg.runTPCTracking = true; cfg.lumiScaleType = sclOpts.lumiType; cfg.lumiScaleMode = sclOpts.lumiMode; + cfg.enableMShape = sclOpts.enableMShapeCorrection; cfg.decompressTPC = decompressTPC; cfg.decompressTPCFromROOT = decompressTPC && inputType == InputType::CompClusters; cfg.caClusterer = caClusterer; diff --git a/Detectors/TPC/workflow/src/TPCScalerSpec.cxx b/Detectors/TPC/workflow/src/TPCScalerSpec.cxx index 752d5bfd7d637..9578ffdbf8d0c 100644 --- a/Detectors/TPC/workflow/src/TPCScalerSpec.cxx +++ b/Detectors/TPC/workflow/src/TPCScalerSpec.cxx @@ -22,6 +22,8 @@ #include "TPCBase/CDBInterface.h" #include "DetectorsBase/GRPGeomHelper.h" #include "TPCCalibration/TPCScaler.h" +#include "TPCCalibration/TPCMShapeScaler.h" +#include "CommonDataFormat/Pair.h" #include "TTree.h" using namespace o2::framework; @@ -34,19 +36,21 @@ namespace tpc class TPCScalerSpec : public Task { public: - TPCScalerSpec(std::shared_ptr req, bool enableWeights) : mCCDBRequest(req), mEnableWeights(enableWeights){}; + TPCScalerSpec(std::shared_ptr req, bool enableIDCs, bool enableMShape) : mCCDBRequest(req), mEnableIDCs(enableIDCs), mEnableMShape(enableMShape){}; void init(framework::InitContext& ic) final { o2::base::GRPGeomHelper::instance().setRequest(mCCDBRequest); mIonDriftTimeMS = ic.options().get("ion-drift-time"); mMaxTimeWeightsMS = ic.options().get("max-time-for-weights"); + mVShapeScalingFac = ic.options().get("v-shape-scaling-factor"); + mEnableWeights = !(ic.options().get("disableWeights")); } void run(ProcessingContext& pc) final { o2::base::GRPGeomHelper::instance().checkUpdates(pc); - if (pc.inputs().isValid("tpcscaler")) { + if (mEnableIDCs && pc.inputs().isValid("tpcscaler")) { pc.inputs().get("tpcscaler"); } @@ -55,6 +59,9 @@ class TPCScalerSpec : public Task pc.inputs().get("tpcscalerw"); } } + if (mEnableMShape && pc.inputs().isValid("vshape")) { + pc.inputs().get("vshape"); + } if (pc.services().get().runNumber != mTPCScaler.getRun()) { LOGP(error, "Run number {} of processed data and run number {} of loaded TPC scaler doesnt match!", pc.services().get().runNumber, mTPCScaler.getRun()); @@ -63,11 +70,27 @@ class TPCScalerSpec : public Task const auto orbitResetTimeMS = o2::base::GRPGeomHelper::instance().getOrbitResetTimeMS(); const auto firstTFOrbit = pc.services().get().firstTForbit; const double timestamp = orbitResetTimeMS + firstTFOrbit * o2::constants::lhc::LHCOrbitMUS * 0.001; - float scalerA = mTPCScaler.getMeanScaler(timestamp, o2::tpc::Side::A); - float scalerC = mTPCScaler.getMeanScaler(timestamp, o2::tpc::Side::C); - float meanScaler = (scalerA + scalerC) / 2; - LOGP(info, "Publishing TPC scaler: {}", meanScaler); - pc.outputs().snapshot(Output{header::gDataOriginTPC, "TPCSCALER"}, meanScaler); + + if (mEnableIDCs) { + float scalerA = mTPCScaler.getMeanScaler(timestamp, o2::tpc::Side::A); + float scalerC = mTPCScaler.getMeanScaler(timestamp, o2::tpc::Side::C); + float meanScaler = (scalerA + scalerC) / 2; + LOGP(info, "Publishing TPC scaler: {} for timestamp: {}, firstTFOrbit: {}", meanScaler, timestamp, firstTFOrbit); + pc.outputs().snapshot(Output{header::gDataOriginTPC, "TPCSCALER"}, meanScaler); + } + + if (mEnableMShape) { + // double tsMiddleTF = timestamp + o2::base::GRPGeomHelper::instance().getNHBFPerTF() / 2 * o2::constants::lhc::LHCOrbitMUS * 0.001; + // // float vshapescaler = mVShapeScalingFac * mVShapeTPCScaler.getScaler(tsMiddleTF); + // // LOGP(info, "Publishing TPC V-shape scaler: {} for timestamp {}, firstTFOrbit: {} with VShapeScalingFac: {}", vshapescaler, long(timestamp), firstTFOrbit, mVShapeScalingFac); + // // pc.outputs().snapshot(Output{header::gDataOriginTPC, "TPCMShapeScaler"}, vshapescaler); + double tsMiddleTF = timestamp + o2::base::GRPGeomHelper::instance().getNHBFPerTF() / 2 * o2::constants::lhc::LHCOrbitMUS * 0.001; + const auto scalers = mVShapeTPCScaler.getScaler(tsMiddleTF); + float vshapescaler = mVShapeScalingFac * scalers.first; + float zPos = scalers.first; + LOGP(info, "Publishing TPC V-shape scaler: {} for timestamp {}, firstTFOrbit: {} with VShapeScalingFac: {}", vshapescaler, long(timestamp), firstTFOrbit, mVShapeScalingFac); + pc.outputs().snapshot(Output{header::gDataOriginTPC, "TPCMShapeScaler"}, dataformats::Pair{vshapescaler, zPos}); + } } void finaliseCCDB(o2::framework::ConcreteDataMatcher& matcher, void* obj) final @@ -98,16 +121,23 @@ class TPCScalerSpec : public Task if (mIonDriftTimeMS == -1) { overWriteIntegrationTime(); } + } else if (matcher == ConcreteDataMatcher(o2::header::gDataOriginTPC, "VSHAPESCALERCCDB", 0)) { + LOGP(info, "Updating V-shape TPC scaler"); + mVShapeTPCScaler.setFromTree(*((TTree*)obj)); } } private: std::shared_ptr mCCDBRequest; ///< info for CCDB request - const bool mEnableWeights{false}; ///< use weights for TPC scalers + const bool mEnableIDCs{true}; ///< enable IDCs + const bool mEnableMShape{false}; ///< enable v shape scalers + bool mEnableWeights{false}; ///< use weights for TPC scalers TPCScalerWeights mScalerWeights{}; ///< scaler weights float mIonDriftTimeMS{-1}; ///< ion drift time float mMaxTimeWeightsMS{500}; ///< maximum integration time when weights are used TPCScaler mTPCScaler; ///< tpc scaler + float mVShapeScalingFac{0}; ///< scale v-shape scalers with this value + TPCMShapeScaler mVShapeTPCScaler; ///< TPC V-shape scalers void overWriteIntegrationTime() { @@ -123,16 +153,21 @@ class TPCScalerSpec : public Task } }; -o2::framework::DataProcessorSpec getTPCScalerSpec(bool enableWeights) +o2::framework::DataProcessorSpec getTPCScalerSpec(bool enableIDCs, bool enableMShape) { std::vector inputs; - inputs.emplace_back("tpcscaler", o2::header::gDataOriginTPC, "TPCSCALERCCDB", 0, Lifetime::Condition, ccdbParamSpec(o2::tpc::CDBTypeMap.at(o2::tpc::CDBType::CalScaler), {}, 1)); // time-dependent - if (enableWeights) { + if (enableIDCs) { + LOGP(info, "Publishing IDC scalers for space-charge distortion fluctuation correction"); + inputs.emplace_back("tpcscaler", o2::header::gDataOriginTPC, "TPCSCALERCCDB", 0, Lifetime::Condition, ccdbParamSpec(o2::tpc::CDBTypeMap.at(o2::tpc::CDBType::CalScaler), {}, 1)); // time-dependent inputs.emplace_back("tpcscalerw", o2::header::gDataOriginTPC, "TPCSCALERWCCDB", 0, Lifetime::Condition, ccdbParamSpec(o2::tpc::CDBTypeMap.at(o2::tpc::CDBType::CalScalerWeights), {}, 0)); // non time-dependent } + if (enableMShape) { + LOGP(info, "Publishing V-shape scalers for V-shape distortion fluctuation correction"); + inputs.emplace_back("vshape", o2::header::gDataOriginTPC, "VSHAPESCALERCCDB", 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 + enableMShape, // GRPECS=true for nHBF per TF false, // GRPLHCIF false, // GRPMagField false, // askMatLUT @@ -141,15 +176,20 @@ o2::framework::DataProcessorSpec getTPCScalerSpec(bool enableWeights) std::vector outputs; outputs.emplace_back(o2::header::gDataOriginTPC, "TPCSCALER", 0, Lifetime::Timeframe); + if (enableMShape) { + outputs.emplace_back(o2::header::gDataOriginTPC, "TPCMShapeScaler", 0, Lifetime::Timeframe); + } return DataProcessorSpec{ "tpc-scaler", inputs, outputs, - AlgorithmSpec{adaptFromTask(ccdbRequest, enableWeights)}, + AlgorithmSpec{adaptFromTask(ccdbRequest, enableIDCs, enableMShape)}, Options{ {"ion-drift-time", VariantType::Float, -1.f, {"Overwrite ion drift time if a value >0 is provided"}}, - {"max-time-for-weights", VariantType::Float, 500.f, {"Maximum possible integration time in ms when weights are used"}}}}; + {"max-time-for-weights", VariantType::Float, 500.f, {"Maximum possible integration time in ms when weights are used"}}, + {"v-shape-scaling-factor", VariantType::Float, 1.f, {"Scale V-shape scaler with this value"}}, + {"disableWeights", VariantType::Bool, false, {"Disable weights for TPC scalers"}}}}; } } // namespace tpc diff --git a/Detectors/TPC/workflow/src/TPCTimeSeriesSpec.cxx b/Detectors/TPC/workflow/src/TPCTimeSeriesSpec.cxx index 4dfa67a7b826b..205234770bd28 100644 --- a/Detectors/TPC/workflow/src/TPCTimeSeriesSpec.cxx +++ b/Detectors/TPC/workflow/src/TPCTimeSeriesSpec.cxx @@ -1288,14 +1288,15 @@ class TPCTimeSeries : public Task void sendOutput(ProcessingContext& pc) { + const long timeMS = o2::base::GRPGeomHelper::instance().getOrbitResetTimeMS() + processing_helpers::getFirstTForbit(pc) * o2::constants::lhc::LHCOrbitMUS / 1000; + mBufferDCA.mTSTPC.setStartTime(timeMS); + mBufferDCA.mTSITSTPC.setStartTime(timeMS); + pc.outputs().snapshot(Output{header::gDataOriginTPC, getDataDescriptionTimeSeries()}, mBufferDCA); // in case of ROOT output also store the TFinfo in the TTree if (!mDisableWriter) { o2::dataformats::TFIDInfo tfinfo; o2::base::TFIDInfoHelper::fillTFIDInfo(pc, tfinfo); - const long timeMS = o2::base::GRPGeomHelper::instance().getOrbitResetTimeMS() + processing_helpers::getFirstTForbit(pc) * o2::constants::lhc::LHCOrbitMUS / 1000; - mBufferDCA.mTSTPC.setStartTime(timeMS); - mBufferDCA.mTSITSTPC.setStartTime(timeMS); pc.outputs().snapshot(Output{header::gDataOriginTPC, getDataDescriptionTPCTimeSeriesTFId()}, tfinfo); } } diff --git a/Detectors/TPC/workflow/src/tpc-calib-gainmap-tracks.cxx b/Detectors/TPC/workflow/src/tpc-calib-gainmap-tracks.cxx index 90dd77391661d..5475995437113 100644 --- a/Detectors/TPC/workflow/src/tpc-calib-gainmap-tracks.cxx +++ b/Detectors/TPC/workflow/src/tpc-calib-gainmap-tracks.cxx @@ -21,6 +21,7 @@ #include "TPCWorkflow/TPCCalibPadGainTracksSpec.h" #include "TPCReaderWorkflow/TPCSectorCompletionPolicy.h" #include "TPCCalibration/CorrectionMapsLoader.h" +#include "TPCWorkflow/TPCScalerSpec.h" using namespace o2::framework; @@ -63,6 +64,10 @@ WorkflowSpec defineDataProcessing(ConfigContext const& config) const std::string polynomialsFile = config.options().get("polynomialsFile"); const auto disablePolynomialsCCDB = config.options().get("disablePolynomialsCCDB"); const auto sclOpt = o2::tpc::CorrectionMapsLoader::parseGlobalOptions(config.options()); - WorkflowSpec workflow{getTPCCalibPadGainTracksSpec(publishAfterTFs, debug, useLastExtractedMapAsReference, polynomialsFile, disablePolynomialsCCDB, sclOpt)}; + WorkflowSpec workflow; + if (sclOpt.needTPCScalersWorkflow()) { + workflow.emplace_back(o2::tpc::getTPCScalerSpec(sclOpt.lumiType == 2, sclOpt.enableMShapeCorrection)); + } + workflow.emplace_back(o2::tpc::getTPCCalibPadGainTracksSpec(publishAfterTFs, debug, useLastExtractedMapAsReference, polynomialsFile, disablePolynomialsCCDB, sclOpt)); return workflow; } diff --git a/Detectors/TPC/workflow/src/tpc-scaler.cxx b/Detectors/TPC/workflow/src/tpc-scaler.cxx index e2c2c8ed79c43..cf0b0e21c1736 100644 --- a/Detectors/TPC/workflow/src/tpc-scaler.cxx +++ b/Detectors/TPC/workflow/src/tpc-scaler.cxx @@ -23,7 +23,10 @@ void customize(std::vector& workflowOptions) // option allowing to set parameters std::vector options{ ConfigParamSpec{"configKeyValues", VariantType::String, "", {"Semicolon separated key=value strings"}}, - {"disableWeights", VariantType::Bool, false, {"Disable weights for TPC scalers"}}}; + {"enable-M-shape-correction", VariantType::Bool, false, {"Enable M-shape distortion correction"}}, + {"disable-IDC-scalers", VariantType::Bool, false, {"Disable TPC scalers for space-charge distortion fluctuation correction"}}, + {"configKeyValues", VariantType::String, "", {"Semicolon separated key=value strings"}}}; + std::swap(workflowOptions, options); } @@ -33,7 +36,8 @@ WorkflowSpec defineDataProcessing(ConfigContext const& config) { WorkflowSpec workflow; o2::conf::ConfigurableParam::updateFromString(config.options().get("configKeyValues")); - const bool enableWeights = !(config.options().get("disableWeights")); - workflow.emplace_back(o2::tpc::getTPCScalerSpec(enableWeights)); + const auto enableMShape = config.options().get("enable-M-shape-correction"); + const auto enableIDCs = !config.options().get("disable-IDC-scalers"); + workflow.emplace_back(o2::tpc::getTPCScalerSpec(enableIDCs, enableMShape)); return workflow; } diff --git a/Detectors/TRD/workflow/src/trd-tracking-workflow.cxx b/Detectors/TRD/workflow/src/trd-tracking-workflow.cxx index c7a489b0fa6e1..04f882d01851a 100644 --- a/Detectors/TRD/workflow/src/trd-tracking-workflow.cxx +++ b/Detectors/TRD/workflow/src/trd-tracking-workflow.cxx @@ -25,6 +25,7 @@ #include "TRDWorkflow/TRDPulseHeightSpec.h" #include "GlobalTrackingWorkflowHelpers/InputHelper.h" #include "TPCCalibration/CorrectionMapsLoader.h" +#include "TPCWorkflow/TPCScalerSpec.h" using namespace o2::framework; using GTrackID = o2::dataformats::GlobalTrackID; @@ -112,6 +113,9 @@ WorkflowSpec defineDataProcessing(ConfigContext const& configcontext) // processing devices o2::framework::WorkflowSpec specs; + if (sclOpt.needTPCScalersWorkflow() && !configcontext.options().get("disable-root-input")) { + specs.emplace_back(o2::tpc::getTPCScalerSpec(sclOpt.lumiType == 2, sclOpt.enableMShapeCorrection)); + } specs.emplace_back(o2::trd::getTRDGlobalTrackingSpec(useMC, srcTRD, trigRecFilterActive, strict, pid, policy, sclOpt)); if (vdexb || gain) { specs.emplace_back(o2::trd::getTRDTrackBasedCalibSpec(srcTRD, vdexb, gain)); diff --git a/GPU/GPUTracking/DataTypes/GPUDataTypes.h b/GPU/GPUTracking/DataTypes/GPUDataTypes.h index c00796e3fd7df..7eeab55b5c242 100644 --- a/GPU/GPUTracking/DataTypes/GPUDataTypes.h +++ b/GPU/GPUTracking/DataTypes/GPUDataTypes.h @@ -216,6 +216,7 @@ template