From d44ccffaef4b55b08f92783d6ed053c05b620be2 Mon Sep 17 00:00:00 2001 From: GiorgioAlbertoLucia Date: Thu, 8 Oct 2026 13:21:04 +0200 Subject: [PATCH 01/12] change to unique function to load map. Add check to ensure that the minimum time resolution is not zero for rescaling --- .../IOTOFSimulation/DPLDigitizerParam.h | 12 +- .../include/IOTOFSimulation/Digitizer.h | 9 +- .../ALICE3/IOTOF/simulation/src/Digitizer.cxx | 106 +++++++++++++----- 3 files changed, 95 insertions(+), 32 deletions(-) diff --git a/Detectors/Upgrades/ALICE3/IOTOF/simulation/include/IOTOFSimulation/DPLDigitizerParam.h b/Detectors/Upgrades/ALICE3/IOTOF/simulation/include/IOTOFSimulation/DPLDigitizerParam.h index 3d2ed2c995c30..b487b3c47c327 100644 --- a/Detectors/Upgrades/ALICE3/IOTOF/simulation/include/IOTOFSimulation/DPLDigitizerParam.h +++ b/Detectors/Upgrades/ALICE3/IOTOF/simulation/include/IOTOFSimulation/DPLDigitizerParam.h @@ -30,8 +30,16 @@ struct DPLDigitizerParam : public o2::conf::ConfigurableParamHelper mChips; //! Chips in the detector, indexed by chip ID std::deque>> mExtraLabelBuffer; //! buffer for multiple mc labels to the same pixel diff --git a/Detectors/Upgrades/ALICE3/IOTOF/simulation/src/Digitizer.cxx b/Detectors/Upgrades/ALICE3/IOTOF/simulation/src/Digitizer.cxx index e4bcde1255b89..ff237079c2156 100644 --- a/Detectors/Upgrades/ALICE3/IOTOF/simulation/src/Digitizer.cxx +++ b/Detectors/Upgrades/ALICE3/IOTOF/simulation/src/Digitizer.cxx @@ -19,6 +19,7 @@ #include "IOTOFSimulation/Digitizer.h" #include "IOTOFSimulation/DPLDigitizerParam.h" #include "DetectorsRaw/HBFUtils.h" +#include "CCDB/BasicCCDBManager.h" #include #include @@ -26,6 +27,7 @@ #include #include +#include #include #include #include @@ -66,15 +68,28 @@ void Digitizer::init() } } - if (!digitizerParams.efficiencyFilePath.empty()) { - loadEfficiencyMap(digitizerParams.efficiencyFilePath); - } + const auto& digitizerParams = o2::iotof::DPLDigitizerParam::Instance(); LOG(info) << "Initializing IOTOF digitizer"; LOG(info) << " Time resolution: " << digitizerParams.timeResolution * 1e3 << " ps"; LOG(info) << " Charge threshold: " << digitizerParams.chargeThreshold << " electrons"; - LOG(info) << " Detection efficiency: " << digitizerParams.efficiency * 100 << " %"; LOG(info) << " Continuous mode: " << (mContinuous ? "ON" : "OFF"); + + loadMap(mEfficiencyMap, digitizerParams.efficiencyMapPath, "hEfficiencyMap"); + if (!mEfficiencyMap) { + LOG(info) << "No efficiency map loaded, using uniform efficiency: " << digitizerParams.efficiency * 100 << " %"; + } + + loadMap(mResolutionMap, digitizerParams.resolutionMapPath, "hResolutionMap"); + if (!mResolutionMap) { + LOG(info) << "No resolution map loaded, using uniform time resolution: " << digitizerParams.timeResolution * 1e3 << " ps"; + } + + loadMap(mTimeOfArrivalMap, digitizerParams.timeOfArrivalMapPath, "hTimeOfArrivalMap"); + if (!mTimeOfArrivalMap) { + LOG(info) << "No time of arrival map loaded"; + } + sSegmentation = o2::iotof::Segmentation::Instance(); } @@ -173,12 +188,11 @@ void Digitizer::processHit(const o2::itsmft::Hit& hit, int evID, int srcID) LOG(debug) << "Hit rejected by efficiency cut at pixel (row,col) = (" << rowIS << ", " << colIS << ")"; continue; } + double smearedTime = smearTime(hitTimeWrtBC, avgHitLocalX[irow][icol], avgHitLocalZ[irow][icol]); const int nElectronsSampled = gRandom->Poisson(electronsPerStep * nEleResp); // Noise can be added here if needed - double smearedTime = smearTime(hitTimeWrtBC); - registerDigits(chip, roFrameAbs, smearedTime, nROF, static_cast(rowIS), static_cast(colIS), nElectronsSampled, label); } @@ -326,12 +340,22 @@ void Digitizer::stepping(const o2::itsmft::Hit& hit, float**& respMatrix, float* } //_______________________________________________________________________ -double Digitizer::smearTime(double time) const +double Digitizer::smearTime(double time, const float x, const float y) const { // Apply Gaussian smearing to simulate detector time resolution const auto& digitizerParams = o2::iotof::DPLDigitizerParam::Instance(); + + float resolutionScaling = 1.; + if (mResolutionMap) { + const float minimumResolution = mResolutionMap->GetMinimum(); + int bin = mResolutionMap->FindBin(x * o2::iotof::Digitizer::cm2um, y * o2::iotof::Digitizer::cm2um); + resolutionScaling = minimumResolution > 0. ? mResolutionMap->GetBinContent(bin) / minimumResolution : 1.; + LOG(debug) << "Time resolution map check: x=" << x * o2::iotof::Digitizer::cm2um << ", y=" << y * o2::iotof::Digitizer::cm2um << ", bin=" << bin << ", resolution=" << minimumResolution; + LOG(debug) << "Time resolution scaling: " << resolutionScaling; + } + if (digitizerParams.timeResolution > 0) { - return time + gRandom->Gaus(0, digitizerParams.timeResolution); + return time + gRandom->Gaus(0, digitizerParams.timeResolution * resolutionScaling); } return time; } @@ -347,31 +371,59 @@ int Digitizer::energyToCharge(float energyLoss) const } //_______________________________________________________________________ -void Digitizer::loadEfficiencyMap(const std::string& filePath) + +void Digitizer::loadMap(TH2D*& map, const std::string& path, const char* mapName) { - // Load the efficiency map from a file - TFile* file = TFile::Open(filePath.c_str()); - if (!file || !file->IsOpen()) { - LOG(error) << "Failed to open efficiency map file: " << filePath; + // Load a map either from CCDB (path prefixed with "ccdb://") or from a ROOT file + // (anything TFile::Open understands: local path, alien://, root://, http://, ...) + static constexpr std::string_view ccdbPrefix = "ccdb://"; + + if (path.empty()) { return; } - auto* rawMap = dynamic_cast(file->Get("hEfficiencyMap")); - if (!rawMap) { - LOG(error) << "Failed to retrieve efficiency map from file: " << filePath; - LOG(error) << "Available keys in the file:"; - TIter next(file->GetListOfKeys()); - TKey* key; - while ((key = dynamic_cast(next()))) { - LOG(error) << " " << key->GetName() << " (" << key->GetClassName() << ")"; + TH2D* rawMap = nullptr; + TFile* file = nullptr; + + if (path.rfind(ccdbPrefix, 0) == 0) { + const std::string ccdbPath = path.substr(ccdbPrefix.size()); + LOG(info) << "Loading " << mapName << " from CCDB: " << ccdbPath; + rawMap = o2::ccdb::BasicCCDBManager::instance().get(ccdbPath); // owned by the CCDB manager + if (!rawMap) { + LOG(error) << "Failed to retrieve " << mapName << " from CCDB path: " << ccdbPath; + return; + } + } else { + LOG(info) << "Loading " << mapName << " from file: " << path; + file = TFile::Open(path.c_str()); + if (!file || !file->IsOpen()) { + LOG(error) << "Failed to open file: " << path; + delete file; + return; + } + rawMap = dynamic_cast(file->Get(mapName)); + if (!rawMap) { + LOG(error) << "Failed to retrieve " << mapName << " from file: " << path; + LOG(error) << "Available keys in the file:"; + TIter next(file->GetListOfKeys()); + TKey* key; + while ((key = dynamic_cast(next()))) { + LOG(error) << " " << key->GetName() << " (" << key->GetClassName() << ")"; + } + file->Close(); + delete file; + return; } - file->Close(); - return; } - mEfficiencyMap = dynamic_cast(rawMap->Clone("mEfficiencyMap")); - mEfficiencyMap->SetDirectory(nullptr); // Detach from file to avoid deletion when file is closed - file->Close(); + map = dynamic_cast(rawMap->Clone()); + map->SetDirectory(nullptr); // Detach from file to avoid deletion when file is closed + LOG(info) << "Loaded " << mapName << " (" << map->GetNbinsX() << " x " << map->GetNbinsY() << " bins)"; + + if (file) { + file->Close(); + delete file; + } } //_______________________________________________________________________ @@ -469,7 +521,7 @@ void Digitizer::registerDigits(Chip& chip, uint32_t roFrame, double time, int nR int tdc = int((time - nbc * o2::constants::lhc::LHCBunchSpacingNS) / digitizerParams.tdcBin); nbc += mEventTime.toLong(); - double absoluteTime = tdc * digitizerParams.tdcBin + nbc * o2::constants::lhc::LHCBunchSpacingNS; + double absoluteTime = tdc * digitizerParams.tdcBin * 1.e-9 + nbc * o2::constants::lhc::LHCBunchSpacingNS; auto key = o2::iotof::Digit::getOrderingKey(nbc, tdc, row, col); o2::iotof::LabeledDigit* existingDigit = chip.findDigit(key); From 676caa8d664a2b279f06963b2ae772bbb08b227b Mon Sep 17 00:00:00 2001 From: GiorgioAlbertoLucia Date: Thu, 8 Oct 2026 13:25:18 +0200 Subject: [PATCH 02/12] add time of arrival map in the smearing process --- .../Upgrades/ALICE3/IOTOF/simulation/src/Digitizer.cxx | 8 +++++++- 1 file changed, 7 insertions(+), 1 deletion(-) diff --git a/Detectors/Upgrades/ALICE3/IOTOF/simulation/src/Digitizer.cxx b/Detectors/Upgrades/ALICE3/IOTOF/simulation/src/Digitizer.cxx index ff237079c2156..5424a01baa63d 100644 --- a/Detectors/Upgrades/ALICE3/IOTOF/simulation/src/Digitizer.cxx +++ b/Detectors/Upgrades/ALICE3/IOTOF/simulation/src/Digitizer.cxx @@ -353,9 +353,15 @@ double Digitizer::smearTime(double time, const float x, const float y) const LOG(debug) << "Time resolution map check: x=" << x * o2::iotof::Digitizer::cm2um << ", y=" << y * o2::iotof::Digitizer::cm2um << ", bin=" << bin << ", resolution=" << minimumResolution; LOG(debug) << "Time resolution scaling: " << resolutionScaling; } + float timeOfArrivalOffset = 0.; + if (mTimeOfArrivalMap) { + int bin = mTimeOfArrivalMap->FindBin(x * o2::iotof::Digitizer::cm2um, y * o2::iotof::Digitizer::cm2um); + timeOfArrivalOffset = mTimeOfArrivalMap->GetBinContent(bin); + LOG(debug) << "Time of arrival map check: x=" << x * o2::iotof::Digitizer::cm2um << ", y=" << y * o2::iotof::Digitizer::cm2um << ", bin=" << bin << ", time offset=" << timeOfArrivalOffset; + } if (digitizerParams.timeResolution > 0) { - return time + gRandom->Gaus(0, digitizerParams.timeResolution * resolutionScaling); + return time + gRandom->Gaus(timeOfArrivalOffset, digitizerParams.timeResolution * resolutionScaling); } return time; } From 67154077116114c6bcc916d5a79e5bd83c02050a Mon Sep 17 00:00:00 2001 From: GiorgioAlbertoLucia Date: Thu, 8 Oct 2026 13:32:24 +0200 Subject: [PATCH 03/12] precompute the resolution map based on input and configured resolution value --- .../include/IOTOFSimulation/Digitizer.h | 3 ++ .../ALICE3/IOTOF/simulation/src/Digitizer.cxx | 28 +++++++++++++++++++ 2 files changed, 31 insertions(+) diff --git a/Detectors/Upgrades/ALICE3/IOTOF/simulation/include/IOTOFSimulation/Digitizer.h b/Detectors/Upgrades/ALICE3/IOTOF/simulation/include/IOTOFSimulation/Digitizer.h index f9d373feb3638..afdcaf303b5f1 100644 --- a/Detectors/Upgrades/ALICE3/IOTOF/simulation/include/IOTOFSimulation/Digitizer.h +++ b/Detectors/Upgrades/ALICE3/IOTOF/simulation/include/IOTOFSimulation/Digitizer.h @@ -95,6 +95,8 @@ class Digitizer : public TObject /// \param mapName Name of the histogram inside the ROOT file (ignored for CCDB) void loadMap(TH2D*& map, const std::string& path, const char* mapName); + void prepareScaledResolutionMap(); + /// Check if the hit passes efficiency cut /// \param x Detector local coordinate x in cm with respect to the center of the sensitive volume. /// \param z Detector local coordinate z in cm with respect to the center of the sensitive volume. @@ -120,6 +122,7 @@ class Digitizer : public TObject const o2::iotof::GeometryTGeo* mGeometry = nullptr; ///< IOTOF geometry TH2D* mEfficiencyMap = nullptr; ///< Efficiency map for the detector TH2D* mResolutionMap = nullptr; ///< Resolution map for the detector + TH2D* mSaledResolutionMap = nullptr; ///< Scaled resolution map for the detector TH2D* mTimeOfArrivalMap = nullptr; ///< Time of arrival map for the detector std::vector mChips; //! Chips in the detector, indexed by chip ID diff --git a/Detectors/Upgrades/ALICE3/IOTOF/simulation/src/Digitizer.cxx b/Detectors/Upgrades/ALICE3/IOTOF/simulation/src/Digitizer.cxx index 5424a01baa63d..64d87f4d7f436 100644 --- a/Detectors/Upgrades/ALICE3/IOTOF/simulation/src/Digitizer.cxx +++ b/Detectors/Upgrades/ALICE3/IOTOF/simulation/src/Digitizer.cxx @@ -81,6 +81,7 @@ void Digitizer::init() } loadMap(mResolutionMap, digitizerParams.resolutionMapPath, "hResolutionMap"); + prepareScaledResolutionMap(); if (!mResolutionMap) { LOG(info) << "No resolution map loaded, using uniform time resolution: " << digitizerParams.timeResolution * 1e3 << " ps"; } @@ -432,6 +433,33 @@ void Digitizer::loadMap(TH2D*& map, const std::string& path, const char* mapName } } +void Digitizer::prepareScaledResolutionMap() +{ + if (!mResolutionMap) { + LOG(warn) << "No resolution map available to prepare scaled resolution map."; + return; + } + + const auto& digitizerParams = o2::iotof::DPLDigitizerParam::Instance(); + const float nominalTimeResolution = digitizerParams.timeResolution; + const float minimumResolution = mResolutionMap->GetMinimum(); + if (minimumResolution <= 0) { + LOG(warn) << "Minimum resolution in the map is non-positive, cannot prepare scaled resolution map."; + return; + } + + mResolutionScalingMap = dynamic_cast(mResolutionMap->Clone("hScaledResolutionMap")); + mResolutionScalingMap->SetDirectory(nullptr); // Detach from file to avoid deletion when file is closed + + for (int binX = 1; binX <= mResolutionScalingMap->GetNbinsX(); ++binX) { + for (int binY = 1; binY <= mResolutionScalingMap->GetNbinsY(); ++binY) { + float originalValue = mResolutionScalingMap->GetBinContent(binX, binY); + float scalingValue = originalValue / minimumResolution; + mResolutionScalingMap->SetBinContent(binX, binY, scalingValue * nominalTimeResolution); + } + } +} + //_______________________________________________________________________ bool Digitizer::isEfficient(const float x, const float z) const { From ac201eb2e398fefafc6c60d8053426e6ba5b55c0 Mon Sep 17 00:00:00 2001 From: GiorgioAlbertoLucia Date: Thu, 8 Oct 2026 13:41:07 +0200 Subject: [PATCH 04/12] add changes required by previous commit --- .../include/IOTOFSimulation/Digitizer.h | 2 +- .../ALICE3/IOTOF/simulation/src/Digitizer.cxx | 28 +++++++++---------- 2 files changed, 14 insertions(+), 16 deletions(-) diff --git a/Detectors/Upgrades/ALICE3/IOTOF/simulation/include/IOTOFSimulation/Digitizer.h b/Detectors/Upgrades/ALICE3/IOTOF/simulation/include/IOTOFSimulation/Digitizer.h index afdcaf303b5f1..2b1d33d8a9fe7 100644 --- a/Detectors/Upgrades/ALICE3/IOTOF/simulation/include/IOTOFSimulation/Digitizer.h +++ b/Detectors/Upgrades/ALICE3/IOTOF/simulation/include/IOTOFSimulation/Digitizer.h @@ -122,7 +122,7 @@ class Digitizer : public TObject const o2::iotof::GeometryTGeo* mGeometry = nullptr; ///< IOTOF geometry TH2D* mEfficiencyMap = nullptr; ///< Efficiency map for the detector TH2D* mResolutionMap = nullptr; ///< Resolution map for the detector - TH2D* mSaledResolutionMap = nullptr; ///< Scaled resolution map for the detector + TH2D* mScaledResolutionMap = nullptr; ///< Scaled resolution map for the detector TH2D* mTimeOfArrivalMap = nullptr; ///< Time of arrival map for the detector std::vector mChips; //! Chips in the detector, indexed by chip ID diff --git a/Detectors/Upgrades/ALICE3/IOTOF/simulation/src/Digitizer.cxx b/Detectors/Upgrades/ALICE3/IOTOF/simulation/src/Digitizer.cxx index 64d87f4d7f436..c76daa6ee2df1 100644 --- a/Detectors/Upgrades/ALICE3/IOTOF/simulation/src/Digitizer.cxx +++ b/Detectors/Upgrades/ALICE3/IOTOF/simulation/src/Digitizer.cxx @@ -346,13 +346,11 @@ double Digitizer::smearTime(double time, const float x, const float y) const // Apply Gaussian smearing to simulate detector time resolution const auto& digitizerParams = o2::iotof::DPLDigitizerParam::Instance(); - float resolutionScaling = 1.; - if (mResolutionMap) { - const float minimumResolution = mResolutionMap->GetMinimum(); - int bin = mResolutionMap->FindBin(x * o2::iotof::Digitizer::cm2um, y * o2::iotof::Digitizer::cm2um); - resolutionScaling = minimumResolution > 0. ? mResolutionMap->GetBinContent(bin) / minimumResolution : 1.; - LOG(debug) << "Time resolution map check: x=" << x * o2::iotof::Digitizer::cm2um << ", y=" << y * o2::iotof::Digitizer::cm2um << ", bin=" << bin << ", resolution=" << minimumResolution; - LOG(debug) << "Time resolution scaling: " << resolutionScaling; + float resolution = digitizerParams.timeResolution; + if (mScaledResolutionMap) { + int bin = mScaledResolutionMap->FindBin(x * o2::iotof::Digitizer::cm2um, y * o2::iotof::Digitizer::cm2um); + resolution = mScaledResolutionMap->GetBinContent(bin); + LOG(debug) << "Time resolution map check: x=" << x * o2::iotof::Digitizer::cm2um << ", y=" << y * o2::iotof::Digitizer::cm2um << ", bin=" << bin << ", resolution=" << resolution; } float timeOfArrivalOffset = 0.; if (mTimeOfArrivalMap) { @@ -362,7 +360,7 @@ double Digitizer::smearTime(double time, const float x, const float y) const } if (digitizerParams.timeResolution > 0) { - return time + gRandom->Gaus(timeOfArrivalOffset, digitizerParams.timeResolution * resolutionScaling); + return time + gRandom->Gaus(timeOfArrivalOffset, resolution); } return time; } @@ -448,14 +446,14 @@ void Digitizer::prepareScaledResolutionMap() return; } - mResolutionScalingMap = dynamic_cast(mResolutionMap->Clone("hScaledResolutionMap")); - mResolutionScalingMap->SetDirectory(nullptr); // Detach from file to avoid deletion when file is closed + mScaledResolutionMap = dynamic_cast(mResolutionMap->Clone("hScaledResolutionMap")); + mScaledResolutionMap->SetDirectory(nullptr); // Detach from file to avoid deletion when file is closed - for (int binX = 1; binX <= mResolutionScalingMap->GetNbinsX(); ++binX) { - for (int binY = 1; binY <= mResolutionScalingMap->GetNbinsY(); ++binY) { - float originalValue = mResolutionScalingMap->GetBinContent(binX, binY); + for (int binX = 1; binX <= mScaledResolutionMap->GetNbinsX(); ++binX) { + for (int binY = 1; binY <= mScaledResolutionMap->GetNbinsY(); ++binY) { + float originalValue = mScaledResolutionMap->GetBinContent(binX, binY); float scalingValue = originalValue / minimumResolution; - mResolutionScalingMap->SetBinContent(binX, binY, scalingValue * nominalTimeResolution); + mScaledResolutionMap->SetBinContent(binX, binY, scalingValue * nominalTimeResolution); } } } @@ -555,7 +553,7 @@ void Digitizer::registerDigits(Chip& chip, uint32_t roFrame, double time, int nR int tdc = int((time - nbc * o2::constants::lhc::LHCBunchSpacingNS) / digitizerParams.tdcBin); nbc += mEventTime.toLong(); - double absoluteTime = tdc * digitizerParams.tdcBin * 1.e-9 + nbc * o2::constants::lhc::LHCBunchSpacingNS; + double absoluteTime = tdc * digitizerParams.tdcBin + nbc * o2::constants::lhc::LHCBunchSpacingNS; auto key = o2::iotof::Digit::getOrderingKey(nbc, tdc, row, col); o2::iotof::LabeledDigit* existingDigit = chip.findDigit(key); From 16cd140e18314af3058c095db2f0a177128f8528 Mon Sep 17 00:00:00 2001 From: GiorgioAlbertoLucia Date: Thu, 8 Oct 2026 17:31:06 +0200 Subject: [PATCH 05/12] add macro to study time resolution --- .../ALICE3/IOTOF/macros/CMakeLists.txt | 11 + .../IOTOF/macros/CheckTimeResolutionIOTOF.C | 559 ++++++++++++++++++ 2 files changed, 570 insertions(+) create mode 100644 Detectors/Upgrades/ALICE3/IOTOF/macros/CheckTimeResolutionIOTOF.C diff --git a/Detectors/Upgrades/ALICE3/IOTOF/macros/CMakeLists.txt b/Detectors/Upgrades/ALICE3/IOTOF/macros/CMakeLists.txt index 8b08fabf6f477..8d6e72b489b18 100644 --- a/Detectors/Upgrades/ALICE3/IOTOF/macros/CMakeLists.txt +++ b/Detectors/Upgrades/ALICE3/IOTOF/macros/CMakeLists.txt @@ -31,3 +31,14 @@ o2_add_test_root_macro(CheckClustersIOTOF.C o2_add_test_root_macro(CheckTopologiesIOTOF.C LABELS iotof COMPILE_ONLY) + +o2_add_test_root_macro(CheckTimeResolutionIOTOF.C + PUBLIC_LINK_LIBRARIES O2::ITSMFTBase + O2::ITSMFTSimulation + O2::IOTOFBase + O2::IOTOFSimulation + O2::MathUtils + O2::SimulationDataFormat + O2::DetectorsBase + O2::Steer + LABELS iotof COMPILE_ONLY) diff --git a/Detectors/Upgrades/ALICE3/IOTOF/macros/CheckTimeResolutionIOTOF.C b/Detectors/Upgrades/ALICE3/IOTOF/macros/CheckTimeResolutionIOTOF.C new file mode 100644 index 0000000000000..703198a727229 --- /dev/null +++ b/Detectors/Upgrades/ALICE3/IOTOF/macros/CheckTimeResolutionIOTOF.C @@ -0,0 +1,559 @@ +// 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 CheckTimeResolutionIOTOF.C +/// \brief Simple macro to check the time resolution of TF3 digits +/// +/// For each digit with a valid MC label, the true time is computed as the +/// MC hit time (relative to the collision) plus the collision time taken from +/// the digitization context. The difference between the digit time and the +/// true time is the time residual, whose width is the time resolution. + +#if !defined(__CLING__) || defined(__ROOTCLING__) +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include + +#include "IOTOFBase/Segmentation.h" +#include "IOTOFBase/IOTOFBaseParam.h" +#include "IOTOFBase/GeometryTGeo.h" +#include "IOTOFSimulation/DPLDigitizerParam.h" +#include "DataFormatsIOTOF/Digit.h" +#include "ITSMFTSimulation/Hit.h" +#include "MathUtils/Utils.h" +#include "SimulationDataFormat/ConstMCTruthContainer.h" +#include "SimulationDataFormat/IOMCTruthContainerView.h" +#include "SimulationDataFormat/MCCompLabel.h" +#include "SimulationDataFormat/MCTrack.h" +#include "SimulationDataFormat/DigitizationContext.h" +#include "CommonDataFormat/InteractionRecord.h" +#include "DetectorsBase/GeometryManager.h" + +#include "DataFormatsITSMFT/ROFRecord.h" + +#endif + +#define ENABLE_UPGRADES + +namespace +{ +constexpr int kNLayers = 2; +constexpr int kNEtaRegions = 2; +const char* kLayerName[kNLayers] = {"ITOF", "OTOF"}; +} // namespace + +/// Fit the time residual distribution with a gaussian in +-2 RMS around the mean and report the resolution +void fitTimeResidual(TH1* h, const char* name, float expectedSigmaPs) +{ + if (!h || h->GetEntries() < 10) { + Warning(name, "Not enough entries to fit the time residual distribution"); + return; + } + const double mean = h->GetMean(), rms = h->GetRMS(); + h->Fit("gaus", "QR", "", mean - 2 * rms, mean + 2 * rms); + auto fit = h->GetFunction("gaus"); + if (!fit) { + return; + } + fit->SetLineColor(kRed); + Info(name, "mean(dt)=%.2f ps, RMS(dt)=%.2f ps, fitted sigma(dt)=%.2f +- %.2f ps, expected sigma(dt)=%.2f ps", + mean, rms, fit->GetParameter(2), fit->GetParError(2), expectedSigmaPs); +} + +/// Convert a TProfile2D filled with the "s" option to a TH2 with the spread of each bin +TH2F* profileToSpread(TProfile2D* prof, const char* name, const char* title) +{ + auto h = new TH2F(name, title, + prof->GetNbinsX(), prof->GetXaxis()->GetXmin(), prof->GetXaxis()->GetXmax(), + prof->GetNbinsY(), prof->GetYaxis()->GetXmin(), prof->GetYaxis()->GetXmax()); + for (int ix = 1; ix <= prof->GetNbinsX(); ++ix) { + for (int iy = 1; iy <= prof->GetNbinsY(); ++iy) { + if (prof->GetBinEntries(prof->GetBin(ix, iy)) < 2) { + continue; + } + h->SetBinContent(ix, iy, prof->GetBinError(ix, iy)); + } + } + return h; +} + +/// Set pad margins leaving room for the axis titles (and for the colour palette and z title of 2D histograms) +void setPadStyle(bool is2D) +{ + gPad->SetTicks(1, 1); + gPad->SetLeftMargin(0.16); + gPad->SetBottomMargin(0.12); + gPad->SetTopMargin(0.08); + gPad->SetRightMargin(is2D ? 0.2 : 0.05); +} + +/// Draw a histogram alone on the canvas and print it as a new page of the given pdf file +void printPage(TCanvas* canv, const char* pdf, TH1* h, const char* opt = "", bool logy = false, bool logz = false) +{ + canv->Clear(); + canv->cd(); + const bool is2D = h->GetDimension() > 1; + setPadStyle(is2D); + gPad->SetLogy(logy); + gPad->SetLogz(logz); + h->GetYaxis()->SetTitleOffset(1.7); + if (is2D) { + h->GetZaxis()->SetTitleOffset(1.6); + } + h->Draw(opt); + canv->Print(pdf, Form("Title:%s", h->GetName())); +} + +/// Human readable name of the selected particle species +TString speciesLabel(int pdg) +{ + switch (std::abs(pdg)) { + case 0: + return "all particles"; + case 11: + return "e^{#pm}"; + case 13: + return "#mu^{#pm}"; + case 211: + return "#pi^{#pm}"; + case 321: + return "K^{#pm}"; + case 2212: + return "p, #bar{p}"; + default: + return Form("|PDG| = %d", std::abs(pdg)); + } +} + +/// Draw the time spectra of ITOF and OTOF split in central and forward pseudorapidity regions, with the total in black. +/// The dashed lines mark the arrival time of a beta = 1 particle at eta = 0. +void drawTimeSpectra(TCanvas* canv, std::array, kNLayers>& hists, const std::array& refTime, + const char* header, float etaCut) +{ + canv->Clear(); + canv->cd(); + setPadStyle(false); + gPad->SetLogy(false); + gPad->SetLogz(false); + + auto hTotal = (TH1F*)hists[0][0]->Clone(Form("%s_total", hists[0][0]->GetName())); + hTotal->Reset(); + for (auto& layerHists : hists) { + for (auto& h : layerHists) { + hTotal->Add(h); + } + } + hTotal->SetLineColor(kBlack); + hTotal->SetFillStyle(0); + hTotal->SetStats(0); + hTotal->GetYaxis()->SetTitleOffset(1.7); + const double peak = hTotal->GetMaximum(); + hTotal->SetMaximum(1.2 * peak); + hTotal->Draw("hist"); + + const int fillColor[kNLayers][kNEtaRegions] = {{kBlue, kRed}, {kBlue + 3, kRed + 3}}; + const int fillStyle[kNLayers] = {3004, 3005}; + auto leg = new TLegend(0.5, 0.55, 0.88, 0.88); + leg->SetBorderSize(0); + leg->SetFillStyle(0); + leg->SetHeader(header); + for (int ie = 0; ie < kNEtaRegions; ++ie) { + for (int il = 0; il < kNLayers; ++il) { + auto h = hists[il][ie]; + h->SetLineColor(fillColor[il][ie]); + h->SetFillColor(fillColor[il][ie]); + h->SetFillStyle(fillStyle[il]); + h->SetStats(0); + h->Draw("hist same"); + leg->AddEntry(h, Form("%s, |#eta| %s %.1f", kLayerName[il], ie == 0 ? "<" : ">", etaCut), "f"); + } + } + hTotal->Draw("hist same"); + leg->Draw(); + + gPad->Update(); + for (int il = 0; il < kNLayers; ++il) { + if (refTime[il] <= 0.f) { + continue; + } + auto line = new TLine(refTime[il], gPad->GetUymin(), refTime[il], 1.03 * peak); + line->SetLineStyle(2); + line->SetLineColor(kGray + 2); + line->Draw(); + auto txt = new TLatex(refTime[il], 1.05 * peak, kLayerName[il]); + txt->SetTextAlign(21); + txt->SetTextSize(0.035); + txt->SetTextColor(kGray + 2); + txt->Draw(); + } +} + +void CheckTimeResolutionIOTOF(std::string digifile = "tf3digits.root", std::string hitfile = "o2sim_HitsTF3.root", std::string kinefile = "o2sim_Kine.root", + std::string inputGeom = "o2sim_geometry.root", std::string collContextFile = "collisioncontext.root", + int pdgSel = 211, float ptMin = 1.f, float ptMax = 10.f, float etaCut = 0.5f, + std::string cfgStr = "IOTOFBase.segmentedInnerTOF=true;IOTOFBase.segmentedOuterTOF=true;IOTOFBase.enableForwardTOF=false;IOTOFBase.enableBackwardTOF=false;") +{ + gStyle->SetPalette(55); + gStyle->SetOptStat(0); // only the fitted histograms get a box, with the fit results + gStyle->SetOptFit(1); + + using namespace o2::base; + using namespace o2::iotof; + + using o2::iotof::Digit; + using o2::itsmft::Hit; + + constexpr float sec2ns = 1e9f; + constexpr float ns2ps = 1e3f; + constexpr float cm2um = 1e4f; + const float speedOfLightCmNs = TMath::C() * 1e-7; // m/s -> cm/ns + + o2::conf::ConfigurableParam::updateFromString(cfgStr); + + const auto& chipInfo = o2::iotof::ChipSpecificsParam::Instance(); + const auto& digiPars = o2::iotof::DPLDigitizerParam::Instance(); + auto seg = o2::iotof::Segmentation::Instance(); + + // Expected resolution: gaussian smearing convoluted with the TDC quantisation (uniform, flooring adds a -tdcBin/2 bias) + const float expectedSigmaPs = std::sqrt(digiPars.timeResolution * digiPars.timeResolution + digiPars.tdcBin * digiPars.tdcBin / 12.f) * ns2ps; + const float dtRangePs = std::max(200.f, 8.f * expectedSigmaPs); + const float tdcBinPs = digiPars.tdcBin * ns2ps; + Info("CheckTimeResolutionIOTOF", "Nominal time resolution %.1f ps, TDC bin %.1f ps -> expected sigma(dt) = %.2f ps", + digiPars.timeResolution * ns2ps, tdcBinPs, expectedSigmaPs); + + TFile* f = TFile::Open("CheckTimeResolution.root", "recreate"); + + TNtuple* nt = new TNtuple("ntt", "digit time ntuple", "id:layer:x:y:z:xLoc:zLoc:dxPix:dzPix:eta:pt:pdg:tTrue:tDig:dt"); + + // Histograms + const float tMaxNs = 13.f; // time window (relative to the collision) for the time spectra + const int nTdcBins = int(tMaxNs / digiPars.tdcBin); + const float halfSizeRow = 0.5f * chipInfo.ActiveMatrixSizeRows(); + const float halfSizeCol = 0.5f * chipInfo.ActiveMatrixSizeCols(); + const float halfPitchRowUm = 0.5f * chipInfo.PitchRow * cm2um; + const float halfPitchColUm = 0.5f * chipInfo.PitchCol * cm2um; + + std::array hTTrue, hTDig, hDt; + std::array hTDigVsTTrue, hDtVsTTrue, hDtVsXLoc, hDtVsZLoc, hDtVsZGlo; + std::array pDtChip, pDtInPixel; + std::array, kNLayers> hTdcSpectrumDig, hTdcSpectrumTrue; + for (int il = 0; il < kNLayers; ++il) { + const char* ln = kLayerName[il]; + hTTrue[il] = new TH1F(Form("h_tTrue_%s", ln), Form("%s: MC true time;t_{true} - t_{collision} (ns);Digits", ln), 260, 0, tMaxNs); + hTDig[il] = new TH1F(Form("h_tDig_%s", ln), Form("%s: digitized time;t_{digit} - t_{collision} (ns);Digits", ln), 260, 0, tMaxNs); + hTDigVsTTrue[il] = new TH2F(Form("h_tDig_vs_tTrue_%s", ln), Form("%s: digitized vs MC true time;t_{true} - t_{collision} (ns);t_{digit} - t_{collision} (ns);Digits", ln), + 260, 0, tMaxNs, 260, 0, tMaxNs); + hDt[il] = new TH1F(Form("h_dt_%s", ln), Form("%s: time residual;t_{digit} - t_{true} (ps);Digits", ln), 400, -dtRangePs, dtRangePs); + hDtVsTTrue[il] = new TH2F(Form("h_dt_vs_tTrue_%s", ln), Form("%s: time residual vs MC true time;t_{true} - t_{collision} (ns);t_{digit} - t_{true} (ps);Digits", ln), + 130, 0, tMaxNs, 200, -dtRangePs, dtRangePs); + hDtVsXLoc[il] = new TH2F(Form("h_dt_vs_xLoc_%s", ln), Form("%s: time residual vs local x in the chip;x_{local} of the pixel (cm);t_{digit} - t_{true} (ps);Digits", ln), + 100, -halfSizeRow, halfSizeRow, 200, -dtRangePs, dtRangePs); + hDtVsZLoc[il] = new TH2F(Form("h_dt_vs_zLoc_%s", ln), Form("%s: time residual vs local z in the chip;z_{local} of the pixel (cm);t_{digit} - t_{true} (ps);Digits", ln), + 100, -halfSizeCol, halfSizeCol, 200, -dtRangePs, dtRangePs); + hDtVsZGlo[il] = new TH2F(Form("h_dt_vs_zGlo_%s", ln), Form("%s: time residual vs global z;z_{global} of the pixel (cm);t_{digit} - t_{true} (ps);Digits", ln), + 200, -400, 400, 200, -dtRangePs, dtRangePs); + pDtChip[il] = new TProfile2D(Form("p_dt_vs_chip_%s", ln), Form("%s: mean time residual vs position in the chip;x_{local} of the pixel (cm);z_{local} of the pixel (cm);#LTt_{digit} - t_{true}#GT (ps)", ln), + 50, -halfSizeRow, halfSizeRow, 50, -halfSizeCol, halfSizeCol, "s"); + pDtInPixel[il] = new TProfile2D(Form("p_dt_vs_inpixel_%s", ln), Form("%s: mean time residual vs position in the pixel;x_{hit} - x_{pixel} (#mum);z_{hit} - z_{pixel} (#mum);#LTt_{digit} - t_{true}#GT (ps)", ln), + 40, -halfPitchRowUm, halfPitchRowUm, 40, -halfPitchColUm, halfPitchColUm, "s"); + for (int ie = 0; ie < kNEtaRegions; ++ie) { + const char* en = ie == 0 ? "central" : "forward"; + hTdcSpectrumDig[il][ie] = new TH1F(Form("h_tdc_digit_%s_%s", ln, en), Form("Digitized time;TDC (%g ps/bin);Entries", tdcBinPs), nTdcBins, 0, nTdcBins); + hTdcSpectrumTrue[il][ie] = new TH1F(Form("h_tdc_true_%s_%s", ln, en), Form("MC true time;t_{true} - t_{collision} (%g ps/bin);Entries", tdcBinPs), nTdcBins, 0, nTdcBins); + } + } + + // Geometry + o2::base::GeometryManager::loadGeometry(inputGeom); + auto* gman = o2::iotof::GeometryTGeo::Instance(); + gman->fillMatrixCache(o2::math_utils::bit2Mask(o2::math_utils::TransformType::L2G)); + + // Collision context: (source, entry) -> collision time + std::map, double> collisionTimeNS; + auto* context = o2::steer::DigitizationContext::loadFromFile(collContextFile); + if (context) { + const auto& records = context->getEventRecords(); + const auto& parts = context->getEventParts(); + for (size_t iColl = 0; iColl < records.size() && iColl < parts.size(); ++iColl) { + for (const auto& part : parts[iColl]) { + collisionTimeNS[{part.sourceID, part.entryID}] = records[iColl].getTimeNS(); + } + } + Info("CheckTimeResolutionIOTOF", "Loaded %zu collisions from %s", records.size(), collContextFile.data()); + } else { + Warning("CheckTimeResolutionIOTOF", "Could not load the collision context from %s: using the ROF BC as collision time (valid only in triggered mode)", collContextFile.data()); + } + + // Hits + TFile* hitFile = TFile::Open(hitfile.data()); + TTree* hitTree = (TTree*)hitFile->Get("o2sim"); + int nevH = hitTree->GetEntries(); // hits are stored as one event per entry + std::vector*> hitArray(nevH, nullptr); + + std::vector> mc2hitVec(nevH); + + // Kinematics (optional): used for the particle species, pT and eta selections + TFile* kineFile = TFile::Open(kinefile.data()); + TTree* kineTree = (kineFile && !kineFile->IsZombie()) ? (TTree*)kineFile->Get("o2sim") : nullptr; + std::vector*> mcTrackArray(nevH, nullptr); + if (!kineTree) { + Warning("CheckTimeResolutionIOTOF", "Could not load the kinematics from %s: no species/pT selection, eta from the digit position", kinefile.data()); + } + + // Digits + TFile* digFile = TFile::Open(digifile.data()); + TTree* digTree = (TTree*)digFile->Get("o2sim"); + + std::vector* digArr{nullptr}; + std::vector* rofRecordsArr{nullptr}; + o2::dataformats::IOMCTruthContainerView* plabelsArr{nullptr}; + + digTree->SetBranchAddress("TF3Digit", &digArr); + digTree->SetBranchAddress("TF3DigitROF", &rofRecordsArr); + digTree->SetBranchAddress("TF3DigitMCTruth", &plabelsArr); + + digTree->GetEntry(0); + + // Load all MC hit (and kinematics) events upfront and build the hit lookup map. + for (int im = 0; im < nevH; ++im) { + hitTree->SetBranchAddress("TF3Hit", &hitArray[im]); + hitTree->GetEntry(im); + auto& mc2hit = mc2hitVec[im]; + for (int ih = hitArray[im]->size(); ih--;) { + const auto& hit = (*hitArray[im])[ih]; + uint64_t key = (uint64_t(hit.GetTrackID()) << 32) + hit.GetDetectorID(); + mc2hit.emplace(key, ih); + } + if (kineTree && im < kineTree->GetEntries()) { + kineTree->SetBranchAddress("MCTrack", &mcTrackArray[im]); + kineTree->GetEntry(im); + } + } + + auto& rofArr = *rofRecordsArr; + + o2::dataformats::ConstMCTruthContainer labels; + plabelsArr->copyandflatten(labels); + + int nNoHit = 0, nNoCollision = 0; + std::array sumRadius{0., 0.}; + std::array nRadius{0, 0}; + + // LOOP on : ROFRecord array + for (unsigned int iROF = 0; iROF < rofArr.size(); ++iROF) { + + const unsigned int rofIndex = rofArr[iROF].getFirstEntry(); + const unsigned int rofNEntries = rofArr[iROF].getNEntries(); + const double rofTimeNS = rofArr[iROF].getBCData().bc2ns(); + + // LOOP on : digits array + for (unsigned int iDigit = rofIndex; iDigit < rofIndex + rofNEntries; iDigit++) { + if (iDigit % 1000 == 0) { + std::cout << "Reading digit " << iDigit << " / " << digArr->size() << std::endl; + } + + const auto& digit = (*digArr)[iDigit]; + Int_t ix = digit.getRow(), iz = digit.getColumn(); + Int_t chipID = digit.getChipIndex(); + Int_t subDetID = gman->getIOTOFLayer(chipID); + if (subDetID < 0 || subDetID >= kNLayers) { + continue; + } + + auto lab = (labels.getLabels(iDigit))[0]; + if (!lab.isValid() || lab.getSourceID() != 0) { // noise or not from the loaded hit file + continue; + } + + // collision time + double tCollNS = rofTimeNS; + if (context) { + auto collEntry = collisionTimeNS.find({lab.getSourceID(), lab.getEventID()}); + if (collEntry == collisionTimeNS.end()) { + nNoCollision++; + continue; + } + tCollNS = collEntry->second; + } + + // get MC info + const int trID = lab.getTrackID(); + std::unordered_map* mc2hit = &mc2hitVec[lab.getEventID()]; + uint64_t key = (uint64_t(trID) << 32) + chipID; + auto hitEntry = mc2hit->find(key); + if (hitEntry == mc2hit->end()) { + LOG(error) << "Failed to find MC hit entry for Tr" << trID << " chipID" << chipID; + nNoHit++; + continue; + } + Hit& hit = (*hitArray[lab.getEventID()])[hitEntry->second]; + + // Local position of the digit (pixel center) and of the hit (mid point between start and end) + Float_t xD = 0.f, zD = 0.f; + seg->detectorToLocal(ix, iz, xD, zD, subDetID); + o2::math_utils::Point3D locD(xD, 0.f, zD); + const auto gloD = gman->getMatrixL2G(chipID)(locD); + + auto xyzLocE = gman->getMatrixL2G(chipID) ^ (hit.GetPos()); + auto xyzLocS = gman->getMatrixL2G(chipID) ^ (hit.GetPosStart()); + const float xH = 0.5f * (xyzLocS.X() + xyzLocE.X()); + const float zH = 0.5f * (xyzLocS.Z() + xyzLocE.Z()); + + // Particle kinematics + const float radius = std::hypot(gloD.X(), gloD.Y()); + float eta = -std::log(std::tan(0.5f * std::atan2(radius, gloD.Z()))); + float pt = -1.f; + int pdg = 0; + const auto* mcTracks = mcTrackArray[lab.getEventID()]; + if (mcTracks && trID >= 0 && trID < (int)mcTracks->size()) { + const auto& mcTrack = (*mcTracks)[trID]; + eta = mcTrack.GetEta(); + pt = mcTrack.GetPt(); + pdg = mcTrack.GetPdgCode(); + } + + // Times (ns) + const double tTrueNS = hit.GetTime() * sec2ns; /// true time relative to the collision + const double tDigNS = digit.getTime() - tCollNS; /// digit time relative to the collision + const float dtPs = (tDigNS - tTrueNS) * ns2ps; /// time residual + + const float dxPixUm = (xH - xD) * cm2um; + const float dzPixUm = (zH - zD) * cm2um; + + float ntVars[] = {float(chipID), float(subDetID), float(gloD.X()), float(gloD.Y()), float(gloD.Z()), xD, zD, dxPixUm, dzPixUm, + eta, pt, float(pdg), float(tTrueNS), float(tDigNS), dtPs}; + nt->Fill(ntVars); + + hTTrue[subDetID]->Fill(tTrueNS); + hTDig[subDetID]->Fill(tDigNS); + hTDigVsTTrue[subDetID]->Fill(tTrueNS, tDigNS); + hDt[subDetID]->Fill(dtPs); + hDtVsTTrue[subDetID]->Fill(tTrueNS, dtPs); + hDtVsXLoc[subDetID]->Fill(xD, dtPs); + hDtVsZLoc[subDetID]->Fill(zD, dtPs); + hDtVsZGlo[subDetID]->Fill(gloD.Z(), dtPs); + pDtChip[subDetID]->Fill(xD, zD, dtPs); + pDtInPixel[subDetID]->Fill(dxPixUm, dzPixUm, dtPs); + sumRadius[subDetID] += radius; + nRadius[subDetID]++; + + // Time spectra for the selected species in the given pT range + const bool passSpecies = !kineTree || pdgSel == 0 || std::abs(pdg) == std::abs(pdgSel); + const bool passPt = !kineTree || (pt > ptMin && pt < ptMax); + if (passSpecies && passPt) { + const int etaRegion = std::abs(eta) < etaCut ? 0 : 1; + hTdcSpectrumDig[subDetID][etaRegion]->Fill(tDigNS / digiPars.tdcBin); + hTdcSpectrumTrue[subDetID][etaRegion]->Fill(tTrueNS / digiPars.tdcBin); + } + + } // end loop on digits array + + } // end loop on ROFRecords + + if (nNoHit || nNoCollision) { + Warning("CheckTimeResolutionIOTOF", "Skipped digits: %d without matching hit, %d without matching collision", nNoHit, nNoCollision); + } + + // Arrival time of a beta = 1 particle at eta = 0, in TDC units + std::array refTimeTdc{0.f, 0.f}; + for (int il = 0; il < kNLayers; ++il) { + if (nRadius[il]) { + refTimeTdc[il] = sumRadius[il] / nRadius[il] / speedOfLightCmNs / digiPars.tdcBin; + } + } + TString header = speciesLabel(kineTree ? pdgSel : 0); + if (kineTree) { + header += Form(", %g < #it{p}_{T} < %g GeV/#it{c}", ptMin, ptMax); + } + + // Histograms created from here on must be written to the output file + f->cd(); + + // One plot per page: each pdf is opened with "[" and closed with "]" + auto canv = new TCanvas("canv", "", 900, 800); + auto openPdf = [&](const char* pdf) { canv->Print(Form("%s[", pdf)); }; + auto closePdf = [&](const char* pdf) { canv->Print(Form("%s]", pdf)); }; + + // Time spectra split in layers and pseudorapidity regions: digitized time and MC true time + const char* pdfSpectra = "tf3digits_time_spectra.pdf"; + openPdf(pdfSpectra); + drawTimeSpectra(canv, hTdcSpectrumDig, refTimeTdc, header, etaCut); + canv->Print(pdfSpectra, "Title:digitized time spectra"); + drawTimeSpectra(canv, hTdcSpectrumTrue, refTimeTdc, header, etaCut); + canv->Print(pdfSpectra, "Title:MC true time spectra"); + closePdf(pdfSpectra); + + // Digitized and MC true time distributions + const char* pdfTime = "tf3digits_time.pdf"; + openPdf(pdfTime); + for (int il = 0; il < kNLayers; ++il) { + printPage(canv, pdfTime, hTTrue[il], "", true); + printPage(canv, pdfTime, hTDig[il], "", true); + printPage(canv, pdfTime, hTDigVsTTrue[il], "colz", false, true); + } + closePdf(pdfTime); + + // Time residual distributions (digit time - true time) + const char* pdfDt = "tf3digits_dt.pdf"; + openPdf(pdfDt); + for (int il = 0; il < kNLayers; ++il) { + fitTimeResidual(hDt[il], kLayerName[il], expectedSigmaPs); + printPage(canv, pdfDt, hDt[il], "", true); + printPage(canv, pdfDt, hDtVsTTrue[il], "colz"); + printPage(canv, pdfDt, hDtVsZGlo[il], "colz"); + } + closePdf(pdfDt); + + // Time residual as a function of the position in the chip + const char* pdfChip = "tf3digits_dt_vs_chip_position.pdf"; + openPdf(pdfChip); + for (int il = 0; il < kNLayers; ++il) { + printPage(canv, pdfChip, hDtVsXLoc[il], "colz"); + printPage(canv, pdfChip, hDtVsZLoc[il], "colz"); + printPage(canv, pdfChip, pDtChip[il], "colz"); + auto hSpread = profileToSpread(pDtChip[il], Form("h_sigmadt_vs_chip_%s", kLayerName[il]), + Form("%s: RMS of the time residual vs position in the chip;x_{local} of the pixel (cm);z_{local} of the pixel (cm);RMS(t_{digit} - t_{true}) (ps)", kLayerName[il])); + printPage(canv, pdfChip, hSpread, "colz"); + } + closePdf(pdfChip); + + // Time residual as a function of the position inside the pixel + const char* pdfInPixel = "tf3digits_dt_inpixel.pdf"; + openPdf(pdfInPixel); + for (int il = 0; il < kNLayers; ++il) { + printPage(canv, pdfInPixel, pDtInPixel[il], "colz"); + auto hSpread = profileToSpread(pDtInPixel[il], Form("h_sigmadt_vs_inpixel_%s", kLayerName[il]), + Form("%s: RMS of the time residual vs position in the pixel;x_{hit} - x_{pixel} (#mum);z_{hit} - z_{pixel} (#mum);RMS(t_{digit} - t_{true}) (ps)", kLayerName[il])); + printPage(canv, pdfInPixel, hSpread, "colz"); + } + closePdf(pdfInPixel); + + f->Write(); + f->Close(); +} From d6271ec518665e403d03d21af1fbb2bcb311fe7f Mon Sep 17 00:00:00 2001 From: GiorgioAlbertoLucia Date: Thu, 8 Oct 2026 17:31:18 +0200 Subject: [PATCH 06/12] fix bug --- Detectors/Upgrades/ALICE3/IOTOF/simulation/src/Digitizer.cxx | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/Detectors/Upgrades/ALICE3/IOTOF/simulation/src/Digitizer.cxx b/Detectors/Upgrades/ALICE3/IOTOF/simulation/src/Digitizer.cxx index c76daa6ee2df1..a6db9ab262955 100644 --- a/Detectors/Upgrades/ALICE3/IOTOF/simulation/src/Digitizer.cxx +++ b/Detectors/Upgrades/ALICE3/IOTOF/simulation/src/Digitizer.cxx @@ -80,13 +80,13 @@ void Digitizer::init() LOG(info) << "No efficiency map loaded, using uniform efficiency: " << digitizerParams.efficiency * 100 << " %"; } - loadMap(mResolutionMap, digitizerParams.resolutionMapPath, "hResolutionMap"); + loadMap(mResolutionMap, digitizerParams.resolutionMapPath, "hSigmaPixel"); prepareScaledResolutionMap(); if (!mResolutionMap) { LOG(info) << "No resolution map loaded, using uniform time resolution: " << digitizerParams.timeResolution * 1e3 << " ps"; } - loadMap(mTimeOfArrivalMap, digitizerParams.timeOfArrivalMapPath, "hTimeOfArrivalMap"); + loadMap(mTimeOfArrivalMap, digitizerParams.timeOfArrivalMapPath, "toa_pixel"); if (!mTimeOfArrivalMap) { LOG(info) << "No time of arrival map loaded"; } From 235a5464822245826d5d991149d086083fe29308 Mon Sep 17 00:00:00 2001 From: GiorgioAlbertoLucia Date: Fri, 9 Oct 2026 10:31:26 +0200 Subject: [PATCH 07/12] fix compilation and execution bug --- Detectors/Upgrades/ALICE3/IOTOF/simulation/src/Digitizer.cxx | 4 +--- Steer/DigitizerWorkflow/src/IOTOFDigitizerSpec.cxx | 5 +++++ 2 files changed, 6 insertions(+), 3 deletions(-) diff --git a/Detectors/Upgrades/ALICE3/IOTOF/simulation/src/Digitizer.cxx b/Detectors/Upgrades/ALICE3/IOTOF/simulation/src/Digitizer.cxx index a6db9ab262955..c0d559ba75add 100644 --- a/Detectors/Upgrades/ALICE3/IOTOF/simulation/src/Digitizer.cxx +++ b/Detectors/Upgrades/ALICE3/IOTOF/simulation/src/Digitizer.cxx @@ -68,8 +68,6 @@ void Digitizer::init() } } - const auto& digitizerParams = o2::iotof::DPLDigitizerParam::Instance(); - LOG(info) << "Initializing IOTOF digitizer"; LOG(info) << " Time resolution: " << digitizerParams.timeResolution * 1e3 << " ps"; LOG(info) << " Charge threshold: " << digitizerParams.chargeThreshold << " electrons"; @@ -189,7 +187,7 @@ void Digitizer::processHit(const o2::itsmft::Hit& hit, int evID, int srcID) LOG(debug) << "Hit rejected by efficiency cut at pixel (row,col) = (" << rowIS << ", " << colIS << ")"; continue; } - double smearedTime = smearTime(hitTimeWrtBC, avgHitLocalX[irow][icol], avgHitLocalZ[irow][icol]); + double smearedTime = smearTime(hitTimeWrtBC, avgHitLocalX[irow][icol] - xPixelCenter, avgHitLocalZ[irow][icol] - zPixelCenter); const int nElectronsSampled = gRandom->Poisson(electronsPerStep * nEleResp); // Noise can be added here if needed diff --git a/Steer/DigitizerWorkflow/src/IOTOFDigitizerSpec.cxx b/Steer/DigitizerWorkflow/src/IOTOFDigitizerSpec.cxx index 008b2bff841c5..1fa43b72684a3 100644 --- a/Steer/DigitizerWorkflow/src/IOTOFDigitizerSpec.cxx +++ b/Steer/DigitizerWorkflow/src/IOTOFDigitizerSpec.cxx @@ -60,6 +60,8 @@ class IOTOFDPLDigitizerTask : o2::base::BaseDPLDigitizer mDigitizer.setGeometry(geom); mDigitizer.init(); + + mROMode = mDigitizer.isContinuous() ? o2::parameters::GRPObject::CONTINUOUS : o2::parameters::GRPObject::PRESENT; } void run(framework::ProcessingContext& pc) @@ -149,6 +151,9 @@ class IOTOFDPLDigitizerTask : o2::base::BaseDPLDigitizer pc.outputs().snapshot(Output{mOrigin, "DIGITSMC2ROF", 0}, dummyMC2ROF); } + LOG(info) << mID.getName() << ": Sending ROMode= " << mROMode << " to GRPUpdater"; + pc.outputs().snapshot(Output{mOrigin, "ROMode", 0}, mROMode); + timer.Stop(); LOG(info) << "Digitization took " << timer.CpuTime() << "s"; From 4897009d32add8df3461b0a01436af277c8fa1b2 Mon Sep 17 00:00:00 2001 From: GiorgioAlbertoLucia Date: Fri, 9 Oct 2026 10:32:01 +0200 Subject: [PATCH 08/12] remove redundant operation --- Detectors/Upgrades/ALICE3/IOTOF/simulation/src/Digitizer.cxx | 2 -- 1 file changed, 2 deletions(-) diff --git a/Detectors/Upgrades/ALICE3/IOTOF/simulation/src/Digitizer.cxx b/Detectors/Upgrades/ALICE3/IOTOF/simulation/src/Digitizer.cxx index c0d559ba75add..e884378374dfc 100644 --- a/Detectors/Upgrades/ALICE3/IOTOF/simulation/src/Digitizer.cxx +++ b/Detectors/Upgrades/ALICE3/IOTOF/simulation/src/Digitizer.cxx @@ -556,8 +556,6 @@ void Digitizer::registerDigits(Chip& chip, uint32_t roFrame, double time, int nR auto key = o2::iotof::Digit::getOrderingKey(nbc, tdc, row, col); o2::iotof::LabeledDigit* existingDigit = chip.findDigit(key); - chip.addDigit(row, col, nElectrons, absoluteTime, nbc, tdc, label); - if (!existingDigit) { // No existing digit, create a new one chip.addDigit(row, col, nElectrons, absoluteTime, nbc, tdc, label); From 1786c8f6344e1476b3da104280704a1b4c9b58aa Mon Sep 17 00:00:00 2001 From: GiorgioAlbertoLucia Date: Fri, 9 Oct 2026 10:52:44 +0200 Subject: [PATCH 09/12] remove unnecessary plot --- .../IOTOF/macros/CheckTimeResolutionIOTOF.C | 33 +++++++------------ 1 file changed, 12 insertions(+), 21 deletions(-) diff --git a/Detectors/Upgrades/ALICE3/IOTOF/macros/CheckTimeResolutionIOTOF.C b/Detectors/Upgrades/ALICE3/IOTOF/macros/CheckTimeResolutionIOTOF.C index 703198a727229..1cd5efb3d1abe 100644 --- a/Detectors/Upgrades/ALICE3/IOTOF/macros/CheckTimeResolutionIOTOF.C +++ b/Detectors/Upgrades/ALICE3/IOTOF/macros/CheckTimeResolutionIOTOF.C @@ -216,7 +216,7 @@ void drawTimeSpectra(TCanvas* canv, std::array, void CheckTimeResolutionIOTOF(std::string digifile = "tf3digits.root", std::string hitfile = "o2sim_HitsTF3.root", std::string kinefile = "o2sim_Kine.root", std::string inputGeom = "o2sim_geometry.root", std::string collContextFile = "collisioncontext.root", - int pdgSel = 211, float ptMin = 1.f, float ptMax = 10.f, float etaCut = 0.5f, + int pdgSel = 211, float ptMin = 1.f, float ptMax = 10.f, float etaCut = 0.5f, float dtMaxPs = 1e4f, std::string cfgStr = "IOTOFBase.segmentedInnerTOF=true;IOTOFBase.segmentedOuterTOF=true;IOTOFBase.enableForwardTOF=false;IOTOFBase.enableBackwardTOF=false;") { gStyle->SetPalette(55); @@ -261,7 +261,7 @@ void CheckTimeResolutionIOTOF(std::string digifile = "tf3digits.root", std::stri std::array hTTrue, hTDig, hDt; std::array hTDigVsTTrue, hDtVsTTrue, hDtVsXLoc, hDtVsZLoc, hDtVsZGlo; - std::array pDtChip, pDtInPixel; + std::array pDtInPixel; std::array, kNLayers> hTdcSpectrumDig, hTdcSpectrumTrue; for (int il = 0; il < kNLayers; ++il) { const char* ln = kLayerName[il]; @@ -278,8 +278,6 @@ void CheckTimeResolutionIOTOF(std::string digifile = "tf3digits.root", std::stri 100, -halfSizeCol, halfSizeCol, 200, -dtRangePs, dtRangePs); hDtVsZGlo[il] = new TH2F(Form("h_dt_vs_zGlo_%s", ln), Form("%s: time residual vs global z;z_{global} of the pixel (cm);t_{digit} - t_{true} (ps);Digits", ln), 200, -400, 400, 200, -dtRangePs, dtRangePs); - pDtChip[il] = new TProfile2D(Form("p_dt_vs_chip_%s", ln), Form("%s: mean time residual vs position in the chip;x_{local} of the pixel (cm);z_{local} of the pixel (cm);#LTt_{digit} - t_{true}#GT (ps)", ln), - 50, -halfSizeRow, halfSizeRow, 50, -halfSizeCol, halfSizeCol, "s"); pDtInPixel[il] = new TProfile2D(Form("p_dt_vs_inpixel_%s", ln), Form("%s: mean time residual vs position in the pixel;x_{hit} - x_{pixel} (#mum);z_{hit} - z_{pixel} (#mum);#LTt_{digit} - t_{true}#GT (ps)", ln), 40, -halfPitchRowUm, halfPitchRowUm, 40, -halfPitchColUm, halfPitchColUm, "s"); for (int ie = 0; ie < kNEtaRegions; ++ie) { @@ -361,7 +359,7 @@ void CheckTimeResolutionIOTOF(std::string digifile = "tf3digits.root", std::stri o2::dataformats::ConstMCTruthContainer labels; plabelsArr->copyandflatten(labels); - int nNoHit = 0, nNoCollision = 0; + int nNoHit = 0, nNoCollision = 0, nOutliers = 0; std::array sumRadius{0., 0.}; std::array nRadius{0, 0}; @@ -450,6 +448,12 @@ void CheckTimeResolutionIOTOF(std::string digifile = "tf3digits.root", std::stri eta, pt, float(pdg), float(tTrueNS), float(tDigNS), dtPs}; nt->Fill(ntVars); + // Reject digits with an unphysical residual (e.g. a corrupted digit time): kept in the ntuple, excluded from the histograms + if (std::abs(dtPs) > dtMaxPs) { + nOutliers++; + continue; + } + hTTrue[subDetID]->Fill(tTrueNS); hTDig[subDetID]->Fill(tDigNS); hTDigVsTTrue[subDetID]->Fill(tTrueNS, tDigNS); @@ -458,7 +462,6 @@ void CheckTimeResolutionIOTOF(std::string digifile = "tf3digits.root", std::stri hDtVsXLoc[subDetID]->Fill(xD, dtPs); hDtVsZLoc[subDetID]->Fill(zD, dtPs); hDtVsZGlo[subDetID]->Fill(gloD.Z(), dtPs); - pDtChip[subDetID]->Fill(xD, zD, dtPs); pDtInPixel[subDetID]->Fill(dxPixUm, dzPixUm, dtPs); sumRadius[subDetID] += radius; nRadius[subDetID]++; @@ -476,8 +479,9 @@ void CheckTimeResolutionIOTOF(std::string digifile = "tf3digits.root", std::stri } // end loop on ROFRecords - if (nNoHit || nNoCollision) { - Warning("CheckTimeResolutionIOTOF", "Skipped digits: %d without matching hit, %d without matching collision", nNoHit, nNoCollision); + if (nNoHit || nNoCollision || nOutliers) { + Warning("CheckTimeResolutionIOTOF", "Skipped digits: %d without matching hit, %d without matching collision, %d with |dt| > %g ps", + nNoHit, nNoCollision, nOutliers, dtMaxPs); } // Arrival time of a beta = 1 particle at eta = 0, in TDC units @@ -530,19 +534,6 @@ void CheckTimeResolutionIOTOF(std::string digifile = "tf3digits.root", std::stri } closePdf(pdfDt); - // Time residual as a function of the position in the chip - const char* pdfChip = "tf3digits_dt_vs_chip_position.pdf"; - openPdf(pdfChip); - for (int il = 0; il < kNLayers; ++il) { - printPage(canv, pdfChip, hDtVsXLoc[il], "colz"); - printPage(canv, pdfChip, hDtVsZLoc[il], "colz"); - printPage(canv, pdfChip, pDtChip[il], "colz"); - auto hSpread = profileToSpread(pDtChip[il], Form("h_sigmadt_vs_chip_%s", kLayerName[il]), - Form("%s: RMS of the time residual vs position in the chip;x_{local} of the pixel (cm);z_{local} of the pixel (cm);RMS(t_{digit} - t_{true}) (ps)", kLayerName[il])); - printPage(canv, pdfChip, hSpread, "colz"); - } - closePdf(pdfChip); - // Time residual as a function of the position inside the pixel const char* pdfInPixel = "tf3digits_dt_inpixel.pdf"; openPdf(pdfInPixel); From 500a5a11a0b3ff33c3db12d5269f2f21d6c9d227 Mon Sep 17 00:00:00 2001 From: GiorgioAlbertoLucia Date: Fri, 9 Oct 2026 11:39:18 +0200 Subject: [PATCH 10/12] fix time of arrival units --- Detectors/Upgrades/ALICE3/IOTOF/simulation/src/Digitizer.cxx | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Detectors/Upgrades/ALICE3/IOTOF/simulation/src/Digitizer.cxx b/Detectors/Upgrades/ALICE3/IOTOF/simulation/src/Digitizer.cxx index e884378374dfc..9a7abec64fb50 100644 --- a/Detectors/Upgrades/ALICE3/IOTOF/simulation/src/Digitizer.cxx +++ b/Detectors/Upgrades/ALICE3/IOTOF/simulation/src/Digitizer.cxx @@ -353,7 +353,7 @@ double Digitizer::smearTime(double time, const float x, const float y) const float timeOfArrivalOffset = 0.; if (mTimeOfArrivalMap) { int bin = mTimeOfArrivalMap->FindBin(x * o2::iotof::Digitizer::cm2um, y * o2::iotof::Digitizer::cm2um); - timeOfArrivalOffset = mTimeOfArrivalMap->GetBinContent(bin); + timeOfArrivalOffset = mTimeOfArrivalMap->GetBinContent(bin) / o2::iotof::Digitizer::ns2ps; // convert to ns LOG(debug) << "Time of arrival map check: x=" << x * o2::iotof::Digitizer::cm2um << ", y=" << y * o2::iotof::Digitizer::cm2um << ", bin=" << bin << ", time offset=" << timeOfArrivalOffset; } From da9e451dc17758c5457ed0a609b37990d2e31337 Mon Sep 17 00:00:00 2001 From: GiorgioAlbertoLucia Date: Fri, 9 Oct 2026 11:41:25 +0200 Subject: [PATCH 11/12] fix syntax --- .../ALICE3/IOTOF/simulation/include/IOTOFSimulation/Digitizer.h | 1 + Detectors/Upgrades/ALICE3/IOTOF/simulation/src/Digitizer.cxx | 2 +- 2 files changed, 2 insertions(+), 1 deletion(-) diff --git a/Detectors/Upgrades/ALICE3/IOTOF/simulation/include/IOTOFSimulation/Digitizer.h b/Detectors/Upgrades/ALICE3/IOTOF/simulation/include/IOTOFSimulation/Digitizer.h index 2b1d33d8a9fe7..2b2a549b37d59 100644 --- a/Detectors/Upgrades/ALICE3/IOTOF/simulation/include/IOTOFSimulation/Digitizer.h +++ b/Detectors/Upgrades/ALICE3/IOTOF/simulation/include/IOTOFSimulation/Digitizer.h @@ -116,6 +116,7 @@ class Digitizer : public TObject return mExtraLabelBuffer[index].get(); } + static constexpr float ps2ns = 1e-3f; ///< picoseconds to nanoseconds conversion static constexpr float sec2ns = 1e9f; ///< seconds to nanoseconds conversion static constexpr float cm2um = 1e4f; ///< centimeters to micrometers conversion diff --git a/Detectors/Upgrades/ALICE3/IOTOF/simulation/src/Digitizer.cxx b/Detectors/Upgrades/ALICE3/IOTOF/simulation/src/Digitizer.cxx index 9a7abec64fb50..65be36ac62cae 100644 --- a/Detectors/Upgrades/ALICE3/IOTOF/simulation/src/Digitizer.cxx +++ b/Detectors/Upgrades/ALICE3/IOTOF/simulation/src/Digitizer.cxx @@ -353,7 +353,7 @@ double Digitizer::smearTime(double time, const float x, const float y) const float timeOfArrivalOffset = 0.; if (mTimeOfArrivalMap) { int bin = mTimeOfArrivalMap->FindBin(x * o2::iotof::Digitizer::cm2um, y * o2::iotof::Digitizer::cm2um); - timeOfArrivalOffset = mTimeOfArrivalMap->GetBinContent(bin) / o2::iotof::Digitizer::ns2ps; // convert to ns + timeOfArrivalOffset = mTimeOfArrivalMap->GetBinContent(bin) * o2::iotof::Digitizer::ps2ns; // convert to ns LOG(debug) << "Time of arrival map check: x=" << x * o2::iotof::Digitizer::cm2um << ", y=" << y * o2::iotof::Digitizer::cm2um << ", bin=" << bin << ", time offset=" << timeOfArrivalOffset; } From fef03486396ab450cdee6dbf3c648915557557d4 Mon Sep 17 00:00:00 2001 From: GiorgioAlbertoLucia Date: Fri, 9 Oct 2026 12:00:14 +0200 Subject: [PATCH 12/12] clang-format --- .../IOTOF/macros/CheckTimeResolutionIOTOF.C | 6 ++--- .../IOTOFSimulation/DPLDigitizerParam.h | 22 +++++++++---------- .../ALICE3/IOTOF/simulation/src/Digitizer.cxx | 2 +- 3 files changed, 15 insertions(+), 15 deletions(-) diff --git a/Detectors/Upgrades/ALICE3/IOTOF/macros/CheckTimeResolutionIOTOF.C b/Detectors/Upgrades/ALICE3/IOTOF/macros/CheckTimeResolutionIOTOF.C index 1cd5efb3d1abe..1e4a309775adb 100644 --- a/Detectors/Upgrades/ALICE3/IOTOF/macros/CheckTimeResolutionIOTOF.C +++ b/Detectors/Upgrades/ALICE3/IOTOF/macros/CheckTimeResolutionIOTOF.C @@ -437,9 +437,9 @@ void CheckTimeResolutionIOTOF(std::string digifile = "tf3digits.root", std::stri } // Times (ns) - const double tTrueNS = hit.GetTime() * sec2ns; /// true time relative to the collision - const double tDigNS = digit.getTime() - tCollNS; /// digit time relative to the collision - const float dtPs = (tDigNS - tTrueNS) * ns2ps; /// time residual + const double tTrueNS = hit.GetTime() * sec2ns; /// true time relative to the collision + const double tDigNS = digit.getTime() - tCollNS; /// digit time relative to the collision + const float dtPs = (tDigNS - tTrueNS) * ns2ps; /// time residual const float dxPixUm = (xH - xD) * cm2um; const float dzPixUm = (zH - zD) * cm2um; diff --git a/Detectors/Upgrades/ALICE3/IOTOF/simulation/include/IOTOFSimulation/DPLDigitizerParam.h b/Detectors/Upgrades/ALICE3/IOTOF/simulation/include/IOTOFSimulation/DPLDigitizerParam.h index b487b3c47c327..48c873d86dc2f 100644 --- a/Detectors/Upgrades/ALICE3/IOTOF/simulation/include/IOTOFSimulation/DPLDigitizerParam.h +++ b/Detectors/Upgrades/ALICE3/IOTOF/simulation/include/IOTOFSimulation/DPLDigitizerParam.h @@ -26,19 +26,19 @@ struct DPLDigitizerParam : public o2::conf::ConfigurableParamHelperFindBin(x * o2::iotof::Digitizer::cm2um, y * o2::iotof::Digitizer::cm2um);