From c739f2de2915e4c1ccea116783d7805d838d8974 Mon Sep 17 00:00:00 2001 From: Maximiliano Puccio Date: Fri, 24 Apr 2026 08:25:46 +0200 Subject: [PATCH 1/6] ALICE3: Add method to get the TRK layer id --- .../ALICE3/TRK/base/include/TRKBase/GeometryTGeo.h | 3 ++- Detectors/Upgrades/ALICE3/TRK/base/src/GeometryTGeo.cxx | 9 +++++++++ 2 files changed, 11 insertions(+), 1 deletion(-) diff --git a/Detectors/Upgrades/ALICE3/TRK/base/include/TRKBase/GeometryTGeo.h b/Detectors/Upgrades/ALICE3/TRK/base/include/TRKBase/GeometryTGeo.h index 576dbf434f757..bfc0b779eb1ca 100644 --- a/Detectors/Upgrades/ALICE3/TRK/base/include/TRKBase/GeometryTGeo.h +++ b/Detectors/Upgrades/ALICE3/TRK/base/include/TRKBase/GeometryTGeo.h @@ -89,7 +89,8 @@ class GeometryTGeo : public o2::detectors::DetMatrixCache int getSubDetID(int index) const; int getPetalCase(int index) const; int getDisk(int index) const; - int getLayer(int index) const; + int getLayer(int index) const; ///< local layer index within the sub-detector (0-based per VD/MLOT) + int getLayerTRK(int index) const; ///< global layer index across the full TRK (VD layers 0..nVD-1, MLOT layers nVD..nTotal-1) int getStave(int index) const; int getHalfStave(int index) const; int getModule(int index) const; diff --git a/Detectors/Upgrades/ALICE3/TRK/base/src/GeometryTGeo.cxx b/Detectors/Upgrades/ALICE3/TRK/base/src/GeometryTGeo.cxx index 10c1c63615d35..fe4abda499239 100644 --- a/Detectors/Upgrades/ALICE3/TRK/base/src/GeometryTGeo.cxx +++ b/Detectors/Upgrades/ALICE3/TRK/base/src/GeometryTGeo.cxx @@ -203,6 +203,15 @@ int GeometryTGeo::getLayer(int index) const return -1; /// -1 if not found } //__________________________________________________________________________ +int GeometryTGeo::getLayerTRK(int index) const +{ + if (getDisk(index) != -1) { + return -1; /// disks do not have a global layer index + } + int subDetID = getSubDetID(index); + return subDetID * o2::trk::constants::VD::petal::nLayers + getLayer(index); // MLOT: offset by number of VD layers +} +//__________________________________________________________________________ int GeometryTGeo::getStave(int index) const { int subDetID = getSubDetID(index); From 6c57378206d096ac0c40f0284d26e27a82a5fd41 Mon Sep 17 00:00:00 2001 From: Maximiliano Puccio Date: Thu, 23 Apr 2026 22:56:46 +0200 Subject: [PATCH 2/6] ALICE3: introduce staggering logic in the digitizer --- .../TRK/base/include/TRKBase/AlmiraParam.h | 29 +- .../ALICE3/TRK/base/include/TRKBase/Specs.h | 2 + .../include/TRKSimulation/DigiParams.h | 42 +-- .../include/TRKSimulation/Digitizer.h | 24 +- .../ALICE3/TRK/simulation/src/DigiParams.cxx | 12 +- .../ALICE3/TRK/simulation/src/Digitizer.cxx | 93 +++--- .../include/TRKWorkflow/DigitReaderSpec.h | 16 +- .../TRK/workflow/src/DigitReaderSpec.cxx | 90 ++++-- .../TRK/workflow/src/DigitWriterSpec.cxx | 105 +++--- .../src/TRKDigitizerSpec.cxx | 302 +++++++++--------- 10 files changed, 391 insertions(+), 324 deletions(-) diff --git a/Detectors/Upgrades/ALICE3/TRK/base/include/TRKBase/AlmiraParam.h b/Detectors/Upgrades/ALICE3/TRK/base/include/TRKBase/AlmiraParam.h index 2048666e21c00..987b459c04764 100644 --- a/Detectors/Upgrades/ALICE3/TRK/base/include/TRKBase/AlmiraParam.h +++ b/Detectors/Upgrades/ALICE3/TRK/base/include/TRKBase/AlmiraParam.h @@ -12,21 +12,40 @@ #ifndef O2_TRK_ALMIRAPARAM_H #define O2_TRK_ALMIRAPARAM_H +#include + #include "CommonConstants/LHCConstants.h" #include "CommonUtils/ConfigurableParam.h" #include "CommonUtils/ConfigurableParamHelper.h" +#include "TRKBase/Specs.h" namespace o2 { namespace trk { -constexpr float DEFAlmiraStrobeDelay = 0.f; ///< default strobe delay in ns wrt ROF start, to be tuned with the real chip response struct AlmiraParam : public o2::conf::ConfigurableParamHelper { - int roFrameLengthInBC = o2::constants::lhc::LHCMaxBunches / 198; ///< ROF length in BC for continuous mode - float strobeDelay = DEFAlmiraStrobeDelay; ///< strobe start in ns wrt ROF start - float strobeLengthCont = -1.; ///< if < 0, full ROF length minus delay - int roFrameBiasInBC = 0; ///< ROF start bias in BC wrt orbit start + static constexpr size_t kNLayers = constants::VD::petal::nLayers + constants::ML::nLayers + constants::OT::nLayers; + static constexpr size_t getNLayers() { return kNLayers; } + + int roFrameLengthInBCPerLayer[kNLayers] = {0}; ///< ROF length in BC per layer + float strobeDelayPerLayer[kNLayers] = {0}; ///< strobe delay in ns per layer + float strobeLengthContPerLayer[kNLayers] = {0}; ///< strobe length in ns per layer + int roFrameBiasInBCPerLayer[kNLayers] = {0}; ///< ROF start bias in BC per layer + int roFrameDelayInBCPerLayer[kNLayers] = {0}; ///< extra ROF delay in BC per layer + + int getROFLengthInBC(int layer) const { + if (roFrameLengthInBCPerLayer[layer] > 0) { + return roFrameLengthInBCPerLayer[layer]; + } else { + return o2::constants::lhc::LHCMaxBunches / 198; + } + } + float getStrobeDelay(int layer) const { return strobeDelayPerLayer[layer]; } + float getStrobeLengthCont(int layer) const { return strobeLengthContPerLayer[layer]; } + int getROFBiasInBC(int layer) const { return roFrameBiasInBCPerLayer[layer]; } + int getROFDelayInBC(int layer) const { return roFrameDelayInBCPerLayer[layer]; } + O2ParamDef(AlmiraParam, "TRKAlmiraParam"); }; diff --git a/Detectors/Upgrades/ALICE3/TRK/base/include/TRKBase/Specs.h b/Detectors/Upgrades/ALICE3/TRK/base/include/TRKBase/Specs.h index 91d6f5669ef33..aad9d045ac6f4 100644 --- a/Detectors/Upgrades/ALICE3/TRK/base/include/TRKBase/Specs.h +++ b/Detectors/Upgrades/ALICE3/TRK/base/include/TRKBase/Specs.h @@ -100,6 +100,7 @@ constexpr int nCols{static_cast(length / chip::pitchZ)}; namespace ML { +constexpr int nLayers{3}; // number of layers in the ML constexpr double width{constants::moduleMLOT::width * 1}; // width of the stave // constexpr double length{constants::moduleMLOT::length * 10}; // length of the stave constexpr double length{124 * cm}; // length of the stave, hardcoded to fit the implemented geometry @@ -117,6 +118,7 @@ constexpr double length{258 * cm}; // len constexpr int nRows{static_cast(width / moduleMLOT::chip::pitchX)}; // number of rows in the halfstave constexpr int nCols{static_cast(length / moduleMLOT::chip::pitchZ)}; // number of columns in the halfstave } // namespace halfstave +constexpr int nLayers{5}; // number of layers in the OT constexpr double width{halfstave::width * 2}; // width of the stave constexpr double length{halfstave::length}; // length of the stave constexpr int nRows{static_cast(width / moduleMLOT::chip::pitchX)}; // number of rows in the stave diff --git a/Detectors/Upgrades/ALICE3/TRK/simulation/include/TRKSimulation/DigiParams.h b/Detectors/Upgrades/ALICE3/TRK/simulation/include/TRKSimulation/DigiParams.h index 3bb58f21dd33b..0253f8bf2bb5d 100644 --- a/Detectors/Upgrades/ALICE3/TRK/simulation/include/TRKSimulation/DigiParams.h +++ b/Detectors/Upgrades/ALICE3/TRK/simulation/include/TRKSimulation/DigiParams.h @@ -15,9 +15,12 @@ #ifndef ALICEO2_TRK_DIGIPARAMS_H #define ALICEO2_TRK_DIGIPARAMS_H +#include + #include #include "ITSMFTSimulation/AlpideSignalTrapezoid.h" #include "ITSMFTSimulation/AlpideSimResponse.h" +#include "TRKBase/AlmiraParam.h" #include "TRKBase/TRKBaseParam.h" #include "TRKBase/GeometryTGeo.h" @@ -50,27 +53,24 @@ class DigiParams void setNoisePerPixel(float v) { mNoisePerPixel = v; } float getNoisePerPixel() const { return mNoisePerPixel; } - void setContinuous(bool v) { mIsContinuous = v; } - bool isContinuous() const { return mIsContinuous; } - - int getROFrameLengthInBC() const { return mROFrameLengthInBC; } - void setROFrameLengthInBC(int n) { mROFrameLengthInBC = n; } + int getROFrameLengthInBC(int layer) const { return mROFrameLayerLengthInBC[layer]; } + void setROFrameLengthInBC(int n, int layer) { mROFrameLayerLengthInBC[layer] = n; } - void setROFrameLength(float ns); - float getROFrameLength() const { return mROFrameLength; } - float getROFrameLengthInv() const { return mROFrameLengthInv; } + void setROFrameLength(float ns, int layer); + float getROFrameLength(int layer) const { return mROFrameLayerLength[layer]; } + float getROFrameLengthInv(int layer) const { return mROFrameLayerLengthInv[layer]; } - void setStrobeDelay(float ns) { mStrobeDelay = ns; } - float getStrobeDelay() const { return mStrobeDelay; } + void setStrobeDelay(float ns, int layer) { mStrobeLayerDelay[layer] = ns; } + float getStrobeDelay(int layer) const { return mStrobeLayerDelay[layer]; } - void setStrobeLength(float ns) { mStrobeLength = ns; } - float getStrobeLength() const { return mStrobeLength; } + void setStrobeLength(float ns, int layer) { mStrobeLayerLength[layer] = ns; } + float getStrobeLength(int layer) const { return mStrobeLayerLength[layer]; } void setTimeOffset(double sec) { mTimeOffset = sec; } double getTimeOffset() const { return mTimeOffset; } - void setROFrameBiasInBC(int n) { mROFrameBiasInBC = n; } - int getROFrameBiasInBC() const { return mROFrameBiasInBC; } + void setROFrameBiasInBC(int n, int layer) { mROFrameLayerBiasInBC[layer] = n; } + int getROFrameBiasInBC(int layer) const { return mROFrameLayerBiasInBC[layer]; } void setChargeThreshold(int v, float frac2Account = 0.1); void setNSimSteps(int v); @@ -102,14 +102,8 @@ class DigiParams private: static constexpr double infTime = 1e99; - bool mIsContinuous = false; ///< flag for continuous simulation float mNoisePerPixel = 1.e-7; ///< Noise per chip - int mROFrameLengthInBC = 0; ///< ROF length in BC for continuos mode - float mROFrameLength = 0; ///< length of RO frame in ns - float mStrobeDelay = 0.; ///< strobe start (in ns) wrt ROF start - float mStrobeLength = 0; ///< length of the strobe in ns (sig. over threshold checked in this window only) double mTimeOffset = -2 * infTime; ///< time offset (in seconds!) to calculate ROFrame from hit time - int mROFrameBiasInBC = 0; ///< misalignment of the ROF start in BC int mChargeThreshold = 75; ///< charge threshold in Nelectrons int mMinChargeToAccount = 7; ///< minimum charge contribution to account int mNSimSteps = 475; ///< number of steps in response simulation @@ -121,12 +115,18 @@ class DigiParams float mIBVbb = 0.0; ///< back bias absolute value for ITS Inner Barrel (in Volt) float mOBVbb = 0.0; ///< back bias absolute value for ITS Outter Barrel (in Volt) + std::array mROFrameLayerLengthInBC; ///< staggering ROF length in BC for continuous mode per layer + std::array mROFrameLayerBiasInBC; ///< staggering ROF bias in BC for continuous mode per layer + std::array mROFrameLayerLength; ///< staggering ROF length in ns for continuous mode per layer + std::array mStrobeLayerLength; ///< staggering strobe length in ns per layer + std::array mStrobeLayerDelay; ///< staggering strobe delay in ns per layer + o2::itsmft::AlpideSignalTrapezoid mSignalShape; ///< signal timeshape parameterization std::unique_ptr mResponse; //!< pointer on external response // auxiliary precalculated parameters - float mROFrameLengthInv = 0; ///< inverse length of RO frame in ns + std::array mROFrameLayerLengthInv; ///< inverse length of RO frame in ns per layer // ClassDef(DigiParams, 2); }; diff --git a/Detectors/Upgrades/ALICE3/TRK/simulation/include/TRKSimulation/Digitizer.h b/Detectors/Upgrades/ALICE3/TRK/simulation/include/TRKSimulation/Digitizer.h index 362de63fb8cb6..5910fc98134aa 100644 --- a/Detectors/Upgrades/ALICE3/TRK/simulation/include/TRKSimulation/Digitizer.h +++ b/Detectors/Upgrades/ALICE3/TRK/simulation/include/TRKSimulation/Digitizer.h @@ -55,18 +55,20 @@ class Digitizer const o2::trk::ChipSimResponse* getChipResponse(int chipID); /// Steer conversion of hits to digits - void process(const std::vector* hits, int evID, int srcID); - void setEventTime(const o2::InteractionTimeRecord& irt); - double getEndTimeOfROFMax() const + void process(const std::vector* hits, int evID, int srcID, int layer); + void setEventTime(const o2::InteractionTimeRecord& irt, int layer); + + void fillOutputContainer(uint32_t maxFrame, int layer); + + void resetROFrameBounds() { - ///< return the time corresponding to end of the last reserved ROFrame : mROFrameMax - return mParams.getROFrameLength() * (mROFrameMax + 1) + mParams.getTimeOffset(); + mROFrameMin = 0; + mROFrameMax = 0; + mNewROFrame = 0; + mIsBeforeFirstRO = false; + mExtraBuff.clear(); } - void setContinuous(bool v) { mParams.setContinuous(v); } - bool isContinuous() const { return mParams.isContinuous(); } - void fillOutputContainer(uint32_t maxFrame = 0xffffffff); - const o2::trk::DigiParams& getDigitParams() const { return mParams; } // provide the common trk::GeometryTGeo to access matrices and segmentation @@ -83,9 +85,9 @@ class Digitizer void setDeadChannelsMap(const o2::itsmft::NoiseMap* mp) { mDeadChanMap = mp; } private: - void processHit(const o2::trk::Hit& hit, uint32_t& maxFr, int evID, int srcID); + void processHit(const o2::trk::Hit& hit, uint32_t& maxFr, int evID, int srcID, int rofLayer); void registerDigits(o2::trk::ChipDigitsContainer& chip, uint32_t roFrame, float tInROF, int nROF, - uint16_t row, uint16_t col, int nEle, o2::MCCompLabel& lbl); + uint16_t row, uint16_t col, int nEle, o2::MCCompLabel& lbl, int layer); ExtraDig* getExtraDigBuffer(uint32_t roFrame) { diff --git a/Detectors/Upgrades/ALICE3/TRK/simulation/src/DigiParams.cxx b/Detectors/Upgrades/ALICE3/TRK/simulation/src/DigiParams.cxx index d5d47b3658b04..3558a6a87ce71 100644 --- a/Detectors/Upgrades/ALICE3/TRK/simulation/src/DigiParams.cxx +++ b/Detectors/Upgrades/ALICE3/TRK/simulation/src/DigiParams.cxx @@ -25,12 +25,12 @@ DigiParams::DigiParams() setNSimSteps(mNSimSteps); } -void DigiParams::setROFrameLength(float lNS) +void DigiParams::setROFrameLength(float lNS, int layer) { // set ROFrame length in nanosecongs - mROFrameLength = lNS; - assert(mROFrameLength > 1.); - mROFrameLengthInv = 1. / mROFrameLength; + mROFrameLayerLength[layer] = lNS; + assert(mROFrameLayerLength[layer] > 1.); + mROFrameLayerLengthInv[layer] = 1. / mROFrameLayerLength[layer]; } void DigiParams::setNSimSteps(int v) @@ -59,10 +59,6 @@ void DigiParams::print() const { // print settings printf("TRK digitization params:\n"); - printf("Continuous readout : %s\n", mIsContinuous ? "ON" : "OFF"); - printf("Readout Frame Length(ns) : %f\n", mROFrameLength); - printf("Strobe delay (ns) : %f\n", mStrobeDelay); - printf("Strobe length (ns) : %f\n", mStrobeLength); printf("Threshold (N electrons) : %d\n", mChargeThreshold); printf("Min N electrons to account : %d\n", mMinChargeToAccount); printf("Number of charge sharing steps : %d\n", mNSimSteps); diff --git a/Detectors/Upgrades/ALICE3/TRK/simulation/src/Digitizer.cxx b/Detectors/Upgrades/ALICE3/TRK/simulation/src/Digitizer.cxx index 31b9a25b7e5f8..890c272fefbc2 100644 --- a/Detectors/Upgrades/ALICE3/TRK/simulation/src/Digitizer.cxx +++ b/Detectors/Upgrades/ALICE3/TRK/simulation/src/Digitizer.cxx @@ -23,6 +23,7 @@ #include #include #include +#include #include // for LOG using o2::itsmft::Digit; @@ -113,14 +114,13 @@ const o2::trk::ChipSimResponse* Digitizer::getChipResponse(int chipID) }; //_______________________________________________________________________ -void Digitizer::process(const std::vector* hits, int evID, int srcID) +void Digitizer::process(const std::vector* hits, int evID, int srcID, int layer) { // digitize single event, the time must have been set beforehand LOG(info) << " Digitizing " << mGeometry->getName() << " (ID: " << mGeometry->getDetID() << ") hits of event " << evID << " from source " << srcID << " at time " << mEventTime.getTimeNS() << " ROFrame = " << mNewROFrame - << " cont.mode: " << isContinuous() << " Min/Max ROFrames " << mROFrameMin << "/" << mROFrameMax; std::cout << "Printing segmentation info: " << std::endl; @@ -128,7 +128,7 @@ void Digitizer::process(const std::vector* hits, int evID, int srcID) // // is there something to flush ? if (mNewROFrame > mROFrameMin) { - fillOutputContainer(mNewROFrame - 1); // flush out all frames preceding the new one + fillOutputContainer(mNewROFrame - 1, layer); // flush out all frames preceding the new one } int nHits = hits->size(); @@ -140,66 +140,55 @@ void Digitizer::process(const std::vector* hits, int evID, int srcID) return (*hits)[lhs].GetDetectorID() < (*hits)[rhs].GetDetectorID(); }); LOG(info) << "Processing " << nHits << " hits"; - for (int i : hitIdx) { - processHit((*hits)[i], mROFrameMax, evID, srcID); - } - - // in the triggered mode store digits after every MC event - // TODO: in the real triggered mode this will not be needed, this is actually for the - // single event processing only - if (!mParams.isContinuous()) { - fillOutputContainer(mROFrameMax); + for (int i : hitIdx | std::views::filter([&](int idx) { + if (layer < 0) { + return true; + } + return mGeometry->getLayerTRK((*hits)[idx].GetDetectorID()) == layer; + })) { + processHit((*hits)[i], mROFrameMax, evID, srcID, layer); } } //_______________________________________________________________________ -void Digitizer::setEventTime(const o2::InteractionTimeRecord& irt) +void Digitizer::setEventTime(const o2::InteractionTimeRecord& irt, int layer) { LOG(info) << "Setting event time to " << irt.getTimeNS() << " ns after orbit 0 bc 0"; // assign event time in ns mEventTime = irt; - if (!mParams.isContinuous()) { - mROFrameMin = 0; // in triggered mode reset the frame counters - mROFrameMax = 0; - } // RO frame corresponding to provided time mCollisionTimeWrtROF = mEventTime.timeInBCNS; // in triggered mode the ROF starts at BC (is there a delay?) - if (mParams.isContinuous()) { - auto nbc = mEventTime.differenceInBC(mIRFirstSampledTF); - - if (mCollisionTimeWrtROF < 0 && nbc > 0) { - nbc--; - } - - if (nbc < 0) { - mNewROFrame = 0; - mIsBeforeFirstRO = true; - } else { - mNewROFrame = nbc / mParams.getROFrameLengthInBC(); - mIsBeforeFirstRO = false; - } + auto nbc = mEventTime.differenceInBC(mIRFirstSampledTF); - LOG(debug) << " NewROFrame " << mNewROFrame << " = " << nbc << "/" << mParams.getROFrameLengthInBC() << " (nbc/mParams.getROFrameLengthInBC()"; + if (mCollisionTimeWrtROF < 0 && nbc > 0) { + nbc--; + } - // in continuous mode depends on starts of periodic readout frame - mCollisionTimeWrtROF += (nbc % mParams.getROFrameLengthInBC()) * o2::constants::lhc::LHCBunchSpacingNS; - } else { + if (nbc < 0) { mNewROFrame = 0; + mIsBeforeFirstRO = true; + } else { + mNewROFrame = nbc / mParams.getROFrameLengthInBC(layer); mIsBeforeFirstRO = false; } + LOG(debug) << " NewROFrame " << mNewROFrame << " = " << nbc << "/" << mParams.getROFrameLengthInBC(layer) << " (nbc/mParams.getROFrameLengthInBC()"; + + // in continuous mode depends on starts of periodic readout frame + mCollisionTimeWrtROF += (nbc % mParams.getROFrameLengthInBC(layer)) * o2::constants::lhc::LHCBunchSpacingNS; + if (mNewROFrame < mROFrameMin) { LOG(error) << "New ROFrame " << mNewROFrame << " (" << irt << ") precedes currently cashed " << mROFrameMin; throw std::runtime_error("deduced ROFrame precedes already processed one"); } - if (mParams.isContinuous() && mROFrameMax < mNewROFrame) { + if (mROFrameMax < mNewROFrame) { mROFrameMax = mNewROFrame - 1; // all frames up to this are finished } } //_______________________________________________________________________ -void Digitizer::fillOutputContainer(uint32_t frameLast) +void Digitizer::fillOutputContainer(uint32_t frameLast, int layer) { // // fill output with digits from min.cached up to requested frame, generating the noise beforehand if (frameLast > mROFrameMax) { @@ -219,7 +208,7 @@ void Digitizer::fillOutputContainer(uint32_t frameLast) auto& extra = *(mExtraBuff.front().get()); for (auto& chip : mChips) { - if (chip.isDisabled()) { + if (chip.isDisabled() || (layer >= 0 && mGeometry->getLayerTRK(chip.getChipIndex()) != layer)) { continue; } chip.addNoise(mROFrameMin, mROFrameMin, &mParams, mGeometry->getSubDetID(chip.getChipIndex()), mGeometry->getLayer(chip.getChipIndex())); /// TODO: add noise @@ -251,11 +240,7 @@ void Digitizer::fillOutputContainer(uint32_t frameLast) } // finalize ROF record rcROF.setNEntries(mDigits->size() - rcROF.getFirstEntry()); // number of digits - if (isContinuous()) { - rcROF.getBCData().setFromLong(mIRFirstSampledTF.toLong() + mROFrameMin * mParams.getROFrameLengthInBC()); - } else { - rcROF.getBCData() = mEventTime; // RSTODO do we need to add trigger delay? - } + rcROF.getBCData().setFromLong(mIRFirstSampledTF.toLong() + mROFrameMin * mParams.getROFrameLengthInBC(layer)); if (mROFRecords) { mROFRecords->push_back(rcROF); } @@ -267,7 +252,7 @@ void Digitizer::fillOutputContainer(uint32_t frameLast) } //_______________________________________________________________________ -void Digitizer::processHit(const o2::trk::Hit& hit, uint32_t& maxFr, int evID, int srcID) +void Digitizer::processHit(const o2::trk::Hit& hit, uint32_t& maxFr, int evID, int srcID, int rofLayer) { int chipID = hit.GetDetectorID(); //// the chip ID at the moment is not referred to the chip but to a wider detector element (e.g. quarter of layer or disk in VD, stave in ML, half stave in OT) int subDetID = mGeometry->getSubDetID(chipID); @@ -297,9 +282,7 @@ void Digitizer::processHit(const o2::trk::Hit& hit, uint32_t& maxFr, int evID, i } return; } - if (isContinuous()) { - timeInROF += mCollisionTimeWrtROF; - } + timeInROF += mCollisionTimeWrtROF; if (mIsBeforeFirstRO && timeInROF < 0) { // disregard this hit because it comes from an event byefore readout starts and it does not effect this RO LOG(debug) << "Ignoring hit with timeInROF = " << timeInROF; @@ -312,9 +295,9 @@ void Digitizer::processHit(const o2::trk::Hit& hit, uint32_t& maxFr, int evID, i } float tTot = mParams.getSignalShape().getMaxDuration(); // frame of the hit signal start wrt event ROFrame - int roFrameRel = int(timeInROF * mParams.getROFrameLengthInv()); + int roFrameRel = int(timeInROF * mParams.getROFrameLengthInv(rofLayer)); // frame of the hit signal end wrt event ROFrame: in the triggered mode we read just 1 frame - uint32_t roFrameRelMax = mParams.isContinuous() ? (timeInROF + tTot) * mParams.getROFrameLengthInv() : roFrameRel; + uint32_t roFrameRelMax = (timeInROF + tTot) * mParams.getROFrameLengthInv(rofLayer); int nFrames = roFrameRelMax + 1 - roFrameRel; uint32_t roFrameMax = mNewROFrame + roFrameRelMax; if (roFrameMax > maxFr) { @@ -509,25 +492,25 @@ void Digitizer::processHit(const o2::trk::Hit& hit, uint32_t& maxFr, int evID, i if (mDeadChanMap && mDeadChanMap->isNoisy(chipID, rowIS, colIS)) { continue; } - registerDigits(chip, roFrameAbs, timeInROF, nFrames, rowIS, colIS, nEle, lbl); + registerDigits(chip, roFrameAbs, timeInROF, nFrames, rowIS, colIS, nEle, lbl, rofLayer); } } } //________________________________________________________________________________ void Digitizer::registerDigits(o2::trk::ChipDigitsContainer& chip, uint32_t roFrame, float tInROF, int nROF, - uint16_t row, uint16_t col, int nEle, o2::MCCompLabel& lbl) + uint16_t row, uint16_t col, int nEle, o2::MCCompLabel& lbl, int layer) { // Register digits for given pixel, accounting for the possible signal contribution to // multiple ROFrame. The signal starts at time tInROF wrt the start of provided roFrame // In every ROFrame we check the collected signal during strobe LOG(debug) << "Registering digits for chip " << chip.getChipIndex() << " at ROFrame " << roFrame << " row " << row << " col " << col << " nEle " << nEle << " label " << lbl; - float tStrobe = mParams.getStrobeDelay() - tInROF; // strobe start wrt signal start - for (int i = 0; i < nROF; i++) { // loop on all the ROFs occupied by the same signal to calculate the charge accumulated in that ROF + float tStrobe = mParams.getStrobeDelay(layer) - tInROF; // strobe start wrt signal start + for (int i = 0; i < nROF; i++) { // loop on all the ROFs occupied by the same signal to calculate the charge accumulated in that ROF uint32_t roFr = roFrame + i; - int nEleROF = mParams.getSignalShape().getCollectedCharge(nEle, tStrobe, tStrobe + mParams.getStrobeLength()); - tStrobe += mParams.getROFrameLength(); // for the next ROF + int nEleROF = mParams.getSignalShape().getCollectedCharge(nEle, tStrobe, tStrobe + mParams.getStrobeLength(layer)); + tStrobe += mParams.getROFrameLength(layer); // for the next ROF // discard too small contributions, they have no chance to produce a digit if (nEleROF < mParams.getMinChargeToAccount()) { /// use threshold instead? diff --git a/Detectors/Upgrades/ALICE3/TRK/workflow/include/TRKWorkflow/DigitReaderSpec.h b/Detectors/Upgrades/ALICE3/TRK/workflow/include/TRKWorkflow/DigitReaderSpec.h index 2a0acd792f4a9..4acea7bd74914 100644 --- a/Detectors/Upgrades/ALICE3/TRK/workflow/include/TRKWorkflow/DigitReaderSpec.h +++ b/Detectors/Upgrades/ALICE3/TRK/workflow/include/TRKWorkflow/DigitReaderSpec.h @@ -12,16 +12,20 @@ #ifndef O2_TRK_DIGITREADER #define O2_TRK_DIGITREADER +#include + #include "TFile.h" #include "TTree.h" #include "DataFormatsITSMFT/Digit.h" #include "DataFormatsITSMFT/GBTCalibData.h" #include "DataFormatsITSMFT/ROFRecord.h" +#include "SimulationDataFormat/IOMCTruthContainerView.h" #include "Framework/DataProcessorSpec.h" #include "Framework/Task.h" #include "Headers/DataHeader.h" #include "DataFormatsITSMFT/ROFRecord.h" #include "DetectorsCommonDataFormats/DetID.h" +#include "TRKBase/AlmiraParam.h" using namespace o2::framework; @@ -41,11 +45,16 @@ class DigitReader : public Task protected: void connectTree(const std::string& filename); + template + void setBranchAddress(const std::string& base, Ptr& addr, int layer = -1); + std::string getBranchName(const std::string& base, int index) const; + + static constexpr int mLayers = o2::trk::AlmiraParam::kNLayers; - std::vector mDigits, *mDigitsPtr = &mDigits; + std::vector*> mDigits{nullptr}; std::vector mCalib, *mCalibPtr = &mCalib; - std::vector mDigROFRec, *mDigROFRecPtr = &mDigROFRec; - std::vector mDigMC2ROFs, *mDigMC2ROFsPtr = &mDigMC2ROFs; + std::vector*> mDigROFRec{nullptr}; + std::vector mPLabels{nullptr}; o2::header::DataOrigin mOrigin = o2::header::gDataOriginInvalid; @@ -64,7 +73,6 @@ class DigitReader : public Task std::string mCalibBranchName = "Calib"; std::string mDigtMCTruthBranchName = "DigitMCTruth"; - std::string mDigtMC2ROFBranchName = "DigitMC2ROF"; }; class TRKDigitReader : public DigitReader diff --git a/Detectors/Upgrades/ALICE3/TRK/workflow/src/DigitReaderSpec.cxx b/Detectors/Upgrades/ALICE3/TRK/workflow/src/DigitReaderSpec.cxx index 09bb1f12a48e4..ec2b6d4d66192 100644 --- a/Detectors/Upgrades/ALICE3/TRK/workflow/src/DigitReaderSpec.cxx +++ b/Detectors/Upgrades/ALICE3/TRK/workflow/src/DigitReaderSpec.cxx @@ -36,12 +36,15 @@ DigitReader::DigitReader(o2::detectors::DetID id, bool useMC, bool useCalib) mDetNameLC = mDetName = id.getName(); mDigTreeName = "o2sim"; + mDigits.resize(mLayers, nullptr); + mDigROFRec.resize(mLayers, nullptr); + mPLabels.resize(mLayers, nullptr); + mDigitBranchName = mDetName + mDigitBranchName; mDigROFBranchName = mDetName + mDigROFBranchName; mCalibBranchName = mDetName + mCalibBranchName; mDigtMCTruthBranchName = mDetName + mDigtMCTruthBranchName; - mDigtMC2ROFBranchName = mDetName + mDigtMC2ROFBranchName; mUseMC = useMC; mUseCalib = useCalib; @@ -58,30 +61,27 @@ void DigitReader::run(ProcessingContext& pc) { auto ent = mTree->GetReadEntry() + 1; assert(ent < mTree->GetEntries()); // this should not happen + mTree->GetEntry(ent); + + for (int iLayer = 0; iLayer < mLayers; ++iLayer) { + LOG(info) << mDetName << "DigitReader on layer " << iLayer << " pushes " << mDigROFRec[iLayer]->size() << " ROFRecords, " + << mDigits[iLayer]->size() << " digits at entry " << ent; - o2::dataformats::IOMCTruthContainerView* plabels = nullptr; - if (mUseMC) { - mTree->SetBranchAddress(mDigtMCTruthBranchName.c_str(), &plabels); + pc.outputs().snapshot(Output{mOrigin, "DIGITSROF", static_cast(iLayer)}, *mDigROFRec[iLayer]); + pc.outputs().snapshot(Output{mOrigin, "DIGITS", static_cast(iLayer)}, *mDigits[iLayer]); + + if (mUseMC) { + auto& sharedlabels = pc.outputs().make>(Output{mOrigin, "DIGITSMCTR", static_cast(iLayer)}); + mPLabels[iLayer]->copyandflatten(sharedlabels); + delete mPLabels[iLayer]; + mPLabels[iLayer] = nullptr; + } } - mTree->GetEntry(ent); - LOG(info) << mDetName << "DigitReader pushes " << mDigROFRec.size() << " ROFRecords, " - << mDigits.size() << " digits at entry " << ent; - // This is a very ugly way of providing DataDescription, which anyway does not need to contain detector name. - // To be fixed once the names-definition class is ready - pc.outputs().snapshot(Output{mOrigin, "DIGITSROF", 0}, mDigROFRec); - pc.outputs().snapshot(Output{mOrigin, "DIGITS", 0}, mDigits); if (mUseCalib) { pc.outputs().snapshot(Output{mOrigin, "GBTCALIB", 0}, mCalib); } - if (mUseMC) { - auto& sharedlabels = pc.outputs().make>(Output{mOrigin, "DIGITSMCTR", 0}); - plabels->copyandflatten(sharedlabels); - delete plabels; - pc.outputs().snapshot(Output{mOrigin, "DIGITSMC2ROF", 0}, mDigMC2ROFs); - } - if (mTree->GetReadEntry() + 1 >= mTree->GetEntries()) { pc.services().get().endOfStream(); pc.services().get().readyToQuit(QuitRequest::Me); @@ -96,35 +96,59 @@ void DigitReader::connectTree(const std::string& filename) mTree.reset((TTree*)mFile->Get(mDigTreeName.c_str())); assert(mTree); - mTree->SetBranchAddress(mDigROFBranchName.c_str(), &mDigROFRecPtr); - mTree->SetBranchAddress(mDigitBranchName.c_str(), &mDigitsPtr); + for (int iLayer = 0; iLayer < mLayers; ++iLayer) { + setBranchAddress(mDigROFBranchName, mDigROFRec[iLayer], iLayer); + setBranchAddress(mDigitBranchName, mDigits[iLayer], iLayer); + if (mUseMC) { + const auto mctruthBranch = getBranchName(mDigtMCTruthBranchName, iLayer); + if (!mTree->GetBranch(mctruthBranch.c_str())) { + throw std::runtime_error("MC data requested but missing branch(es) at layer " + std::to_string(iLayer) + + ": " + mctruthBranch); + } + setBranchAddress(mDigtMCTruthBranchName, mPLabels[iLayer], iLayer); + } + } + if (mUseCalib) { if (!mTree->GetBranch(mCalibBranchName.c_str())) { throw std::runtime_error("GBT calibration data requested but not found in the tree"); } - mTree->SetBranchAddress(mCalibBranchName.c_str(), &mCalibPtr); - } - if (mUseMC) { - if (!mTree->GetBranch(mDigtMC2ROFBranchName.c_str()) || !mTree->GetBranch(mDigtMCTruthBranchName.c_str())) { - throw std::runtime_error("MC data requested but not found in the tree"); - } - mTree->SetBranchAddress(mDigtMC2ROFBranchName.c_str(), &mDigMC2ROFsPtr); + setBranchAddress(mCalibBranchName, mCalibPtr); } LOG(info) << "Loaded tree from " << filename << " with " << mTree->GetEntries() << " entries"; } +std::string DigitReader::getBranchName(const std::string& base, int index) const +{ + if (index >= 0) { + return base + "_" + std::to_string(index); + } + return base; +} + +template +void DigitReader::setBranchAddress(const std::string& base, Ptr& addr, int layer) +{ + const auto name = getBranchName(base, layer); + if (Int_t ret = mTree->SetBranchAddress(name.c_str(), &addr); ret != 0) { + LOGP(fatal, "failed to set branch address for {} ret={}", name, ret); + } +} + DataProcessorSpec getTRKDigitReaderSpec(bool useMC, bool useCalib, std::string defname) { + static constexpr int nLayers = o2::trk::AlmiraParam::kNLayers; std::vector outputSpec; - outputSpec.emplace_back("TRK", "DIGITS", 0, Lifetime::Timeframe); - outputSpec.emplace_back("TRK", "DIGITSROF", 0, Lifetime::Timeframe); + for (int iLayer = 0; iLayer < nLayers; ++iLayer) { + outputSpec.emplace_back("TRK", "DIGITS", iLayer, Lifetime::Timeframe); + outputSpec.emplace_back("TRK", "DIGITSROF", iLayer, Lifetime::Timeframe); + if (useMC) { + outputSpec.emplace_back("TRK", "DIGITSMCTR", iLayer, Lifetime::Timeframe); + } + } if (useCalib) { outputSpec.emplace_back("TRK", "GBTCALIB", 0, Lifetime::Timeframe); } - if (useMC) { - outputSpec.emplace_back("TRK", "DIGITSMCTR", 0, Lifetime::Timeframe); - outputSpec.emplace_back("TRK", "DIGITSMC2ROF", 0, Lifetime::Timeframe); - } return DataProcessorSpec{ "trk-digit-reader", diff --git a/Detectors/Upgrades/ALICE3/TRK/workflow/src/DigitWriterSpec.cxx b/Detectors/Upgrades/ALICE3/TRK/workflow/src/DigitWriterSpec.cxx index 2a743551adddb..591b084aee3ba 100644 --- a/Detectors/Upgrades/ALICE3/TRK/workflow/src/DigitWriterSpec.cxx +++ b/Detectors/Upgrades/ALICE3/TRK/workflow/src/DigitWriterSpec.cxx @@ -9,9 +9,12 @@ // granted to it by virtue of its status as an Intergovernmental Organization // or submit itself to any jurisdiction. -/// @brief Processor spec for a ROOT file writer for ITSMFT digits +/// @brief Processor spec for a ROOT file writer for TRK digits (per-layer) #include "TRKWorkflow/DigitWriterSpec.h" +#include "Framework/ConcreteDataMatcher.h" +#include "Framework/DataRef.h" +#include "TRKBase/AlmiraParam.h" #include "DPLUtils/MakeRootTreeWriterSpec.h" #include "DataFormatsITSMFT/Digit.h" #include "DataFormatsITSMFT/GBTCalibData.h" @@ -24,6 +27,7 @@ #include #include #include +#include using namespace o2::framework; using SubSpecificationType = o2::framework::DataAllocator::SubSpecificationType; @@ -37,16 +41,22 @@ template using BranchDefinition = MakeRootTreeWriterSpec::BranchDefinition; using MCCont = o2::dataformats::ConstMCTruthContainer; -/// create the processor spec -/// describing a processor receiving digits for ITS/MFT and writing them to file -DataProcessorSpec getDigitWriterSpec(bool mctruth, bool dec, bool calib, o2::header::DataOrigin detOrig, o2::detectors::DetID detId) +DataProcessorSpec getTRKDigitWriterSpec(bool mctruth, bool dec, bool calib) { - std::string detStr = o2::detectors::DetID::getName(detId); - std::string detStrL = dec ? "o2_" : ""; // for decoded digits prepend by o2 - detStrL += detStr; - std::transform(detStrL.begin(), detStrL.end(), detStrL.begin(), ::tolower); - auto logger = [](std::vector const& inDigits) { - LOG(info) << "RECEIVED DIGITS SIZE " << inDigits.size(); + static constexpr o2::header::DataOrigin Origin = o2::header::gDataOriginTRK; + const int mLayers = o2::trk::AlmiraParam::kNLayers; + std::string detStr = "TRK"; + std::string detStrL = dec ? "o2_trk" : "trk"; + + auto digitSizes = std::make_shared>(mLayers, 0); + auto digitSizeGetter = [digitSizes](std::vector const& inDigits, DataRef const& ref) { + auto const* dh = DataRefUtils::getHeader(ref); + (*digitSizes)[dh->subSpecification] = inDigits.size(); + }; + auto rofSizes = std::make_shared>(mLayers, 0); + auto rofSizeGetter = [rofSizes](std::vector const& inROFs, DataRef const& ref) { + auto const* dh = DataRefUtils::getHeader(ref); + (*rofSizes)[dh->subSpecification] = inROFs.size(); }; // the callback to be set as hook for custom action when the writer is closed @@ -61,16 +71,18 @@ DataProcessorSpec getDigitWriterSpec(bool mctruth, bool dec, bool calib, o2::hea nent = n; } outputtree->SetEntries(nent); - outputtree->Write("", TObject::kOverwrite); + outputfile->Write("", TObject::kOverwrite); outputfile->Close(); }; // handler for labels - // This is necessary since we can't store the original label buffer in a ROOT entry -- as is -- if it exceeds a certain size. - // We therefore convert it to a special split class. - auto fillLabels = [](TBranch& branch, std::vector const& labelbuffer, DataRef const& /*ref*/) { + auto fillLabels = [detStr, digitSizes, rofSizes](TBranch& branch, std::vector const& labelbuffer, DataRef const& ref) { o2::dataformats::ConstMCTruthContainerView labels(labelbuffer); - LOG(info) << "WRITING " << labels.getNElements() << " LABELS "; + auto const* dh = DataRefUtils::getHeader(ref); + auto layer = static_cast(dh->subSpecification); + LOG(info) << detStr << ": WRITING " << labels.getNElements() << " LABELS" + << std::format(" FOR LAYER {}", layer) << " WITH " << (*digitSizes)[layer] + << " DIGITS IN " << (*rofSizes)[layer] << " ROFS"; o2::dataformats::IOMCTruthContainerView outputcontainer; auto ptr = &outputcontainer; @@ -80,30 +92,49 @@ DataProcessorSpec getDigitWriterSpec(bool mctruth, bool dec, bool calib, o2::hea br->ResetAddress(); }; - return MakeRootTreeWriterSpec((detStr + "DigitWriter" + (dec ? "_dec" : "")).c_str(), + auto getIndex = [](DataRef const& ref) -> size_t { + auto const* dh = DataRefUtils::getHeader(ref); + return static_cast(dh->subSpecification); + }; + auto getName = [](std::string base, size_t index) -> std::string { + return base + "_" + std::to_string(index); + }; + + std::vector vecInpSpecDig, vecInpSpecROF, vecInpSpecLbl; + vecInpSpecDig.reserve(mLayers); + vecInpSpecROF.reserve(mLayers); + vecInpSpecLbl.reserve(mLayers); + for (int iLayer = 0; iLayer < mLayers; iLayer++) { + vecInpSpecDig.emplace_back(getName(detStr + "digits", iLayer), Origin, "DIGITS", iLayer); + vecInpSpecROF.emplace_back(getName(detStr + "digitsROF", iLayer), Origin, "DIGITSROF", iLayer); + vecInpSpecLbl.emplace_back(getName(detStr + "_digitsMCTR", iLayer), Origin, "DIGITSMCTR", iLayer); + } + + return MakeRootTreeWriterSpec(("TRKDigitWriter" + std::string(dec ? "_dec" : "")).c_str(), (detStrL + "digits.root").c_str(), - MakeRootTreeWriterSpec::TreeAttributes{"o2sim", "Digits tree"}, + MakeRootTreeWriterSpec::TreeAttributes{.name = "o2sim", .title = detStr + " Digits tree"}, MakeRootTreeWriterSpec::CustomClose(finishWriting), - // in case of labels we first read them as std::vector and process them correctly in the fillLabels hook - BranchDefinition>{InputSpec{"digitsMCTR", detOrig, "DIGITSMCTR", 0}, - (detStr + "DigitMCTruth").c_str(), - (mctruth ? 1 : 0), fillLabels}, - BranchDefinition>{InputSpec{"digitsMC2ROF", detOrig, "DIGITSMC2ROF", 0}, - (detStr + "DigitMC2ROF").c_str(), - (mctruth ? 1 : 0)}, - BranchDefinition>{InputSpec{"digits", detOrig, "DIGITS", 0}, - (detStr + "Digit").c_str(), - logger}, - BranchDefinition>{InputSpec{"calib", detOrig, "GBTCALIB", 0}, - (detStr + "Calib").c_str(), - (calib ? 1 : 0)}, - BranchDefinition>{InputSpec{"digitsROF", detOrig, "DIGITSROF", 0}, - (detStr + "DigitROF").c_str()})(); -} - -DataProcessorSpec getTRKDigitWriterSpec(bool mctruth, bool dec, bool calib) -{ - return getDigitWriterSpec(mctruth, dec, calib, o2::header::gDataOriginTRK, o2::detectors::DetID::TRK); + BranchDefinition>{vecInpSpecDig, + detStr + "Digit", "digit-branch", + mLayers, + digitSizeGetter, + getIndex, + getName}, + BranchDefinition>{vecInpSpecROF, + detStr + "DigitROF", "digit-rof-branch", + mLayers, + rofSizeGetter, + getIndex, + getName}, + BranchDefinition>{vecInpSpecLbl, + detStr + "DigitMCTruth", "digit-mctruth-branch", + (mctruth ? mLayers : 0), + fillLabels, + getIndex, + getName}, + BranchDefinition>{InputSpec{detStr + "calib", ConcreteDataTypeMatcher{Origin, "GBTCALIB"}}, + detStr + "Calib", "digit-calib-branch", + (calib ? 1 : 0)})(); } } // end namespace trk diff --git a/Steer/DigitizerWorkflow/src/TRKDigitizerSpec.cxx b/Steer/DigitizerWorkflow/src/TRKDigitizerSpec.cxx index 8957ebed223b2..06d922cc1a117 100644 --- a/Steer/DigitizerWorkflow/src/TRKDigitizerSpec.cxx +++ b/Steer/DigitizerWorkflow/src/TRKDigitizerSpec.cxx @@ -30,6 +30,7 @@ #include "TRKSimulation/DPLDigitizerParam.h" #include "TRKBase/AlmiraParam.h" #include "TRKBase/GeometryTGeo.h" +#include "TRKBase/Specs.h" #include "TRKBase/TRKBaseParam.h" #include @@ -47,11 +48,13 @@ namespace std::vector makeOutChannels(o2::header::DataOrigin detOrig, bool mctruth) { std::vector outputs; - outputs.emplace_back(detOrig, "DIGITS", 0, Lifetime::Timeframe); - outputs.emplace_back(detOrig, "DIGITSROF", 0, Lifetime::Timeframe); - if (mctruth) { - outputs.emplace_back(detOrig, "DIGITSMC2ROF", 0, Lifetime::Timeframe); - outputs.emplace_back(detOrig, "DIGITSMCTR", 0, Lifetime::Timeframe); + for (uint32_t iLayer = 0; iLayer < o2::trk::AlmiraParam::getNLayers(); ++iLayer) { + outputs.emplace_back(detOrig, "DIGITS", iLayer, Lifetime::Timeframe); + outputs.emplace_back(detOrig, "DIGITSROF", iLayer, Lifetime::Timeframe); + if (mctruth) { + outputs.emplace_back(detOrig, "DIGITSMC2ROF", iLayer, Lifetime::Timeframe); + outputs.emplace_back(detOrig, "DIGITSMCTR", iLayer, Lifetime::Timeframe); + } } outputs.emplace_back(detOrig, "ROMode", 0, Lifetime::Timeframe); return outputs; @@ -99,163 +102,151 @@ class TRKDPLDigitizerTask : BaseDPLDigitizer timer.Start(); LOG(info) << " CALLING TRK DIGITIZATION "; - mDigitizer.setDigits(&mDigits); - mDigitizer.setROFRecords(&mROFRecords); - mDigitizer.setMCLabels(&mLabels); + auto& eventParts = context->getEventParts(withQED); + uint64_t nDigits{0}; + for (uint32_t iLayer = 0; iLayer < static_cast(mLayers); ++iLayer) { + mDigits[iLayer].clear(); + mROFRecords[iLayer].clear(); + mROFRecordsAccum[iLayer].clear(); + if (mWithMCTruth) { + mLabels[iLayer].clear(); + mLabelsAccum[iLayer].clear(); + mMC2ROFRecordsAccum[iLayer].clear(); + } + + mDigitizer.setDigits(&mDigits[iLayer]); + mDigitizer.setROFRecords(&mROFRecords[iLayer]); + mDigitizer.setMCLabels(&mLabels[iLayer]); + mDigitizer.resetROFrameBounds(); - // digits are directly put into DPL owned resource - auto& digitsAccum = pc.outputs().make>(Output{mOrigin, "DIGITS", 0}); + // digits are directly put into DPL owned resource + auto& digitsAccum = pc.outputs().make>(Output{mOrigin, "DIGITS", iLayer}); - const int roFrameLengthInBC = mDigitizer.getParams().getROFrameLengthInBC(); - const int nROFsPerOrbit = o2::constants::lhc::LHCMaxBunches / roFrameLengthInBC; - const int nROFsTF = nROFsPerOrbit * raw::HBFUtils::Instance().getNOrbitsPerTF(); - mROFRecordsAccum.reserve(nROFsTF); + const int roFrameLengthInBC = mDigitizer.getParams().getROFrameLengthInBC(iLayer); + const int nROFsPerOrbit = o2::constants::lhc::LHCMaxBunches / roFrameLengthInBC; + const int nROFsTF = nROFsPerOrbit * raw::HBFUtils::Instance().getNOrbitsPerTF(); + mROFRecordsAccum[iLayer].reserve(nROFsTF); - auto accumulate = [this, &digitsAccum]() { - // accumulate result of single event processing, called after processing every event supplied - // AND after the final flushing via digitizer::fillOutputContainer - if (mDigits.empty()) { - LOG(debug) << "No digits to accumulate"; - return; // no digits were flushed, nothing to accumulate - } - LOG(debug) << "Accumulating " << mDigits.size() << " digits "; - auto ndigAcc = digitsAccum.size(); - std::copy(mDigits.begin(), mDigits.end(), std::back_inserter(digitsAccum)); - - // fix ROFrecords references on ROF entries - auto nROFRecsOld = mROFRecordsAccum.size(); - - for (int i = 0; i < mROFRecords.size(); i++) { - auto& rof = mROFRecords[i]; - rof.setFirstEntry(ndigAcc + rof.getFirstEntry()); - rof.print(); - - if (mFixMC2ROF < mMC2ROFRecordsAccum.size()) { // fix ROFRecord entry in MC2ROF records - for (int m2rid = mFixMC2ROF; m2rid < mMC2ROFRecordsAccum.size(); m2rid++) { - // need to register the ROFRecors entry for MC event starting from this entry - auto& mc2rof = mMC2ROFRecordsAccum[m2rid]; - if (rof.getROFrame() == mc2rof.minROF) { - mFixMC2ROF++; - mc2rof.rofRecordID = nROFRecsOld + i; - mc2rof.print(); - } - } + auto accumulate = [this, &digitsAccum, &iLayer]() { + // accumulate result of single event processing on one layer, called after each collision + // and after the final flushing via digitizer::fillOutputContainer + if (mDigits[iLayer].empty()) { + return; } - } + auto ndigAcc = digitsAccum.size(); + std::copy(mDigits[iLayer].begin(), mDigits[iLayer].end(), std::back_inserter(digitsAccum)); - std::copy(mROFRecords.begin(), mROFRecords.end(), std::back_inserter(mROFRecordsAccum)); - if (mWithMCTruth) { - mLabelsAccum.mergeAtBack(mLabels); - } - LOG(info) << "Added " << mDigits.size() << " digits "; - // clean containers from already accumulated stuff - mLabels.clear(); - mDigits.clear(); - mROFRecords.clear(); - }; // end accumulate lambda + for (auto& rof : mROFRecords[iLayer]) { + rof.setFirstEntry(ndigAcc + rof.getFirstEntry()); + } - auto& eventParts = context->getEventParts(withQED); - // loop over all composite collisions given from context (aka loop over all the interaction records) - const int bcShift = mDigitizer.getParams().getROFrameBiasInBC(); - // loop over all composite collisions given from context (aka loop over all the interaction records) - for (size_t collID = 0; collID < timesview.size(); ++collID) { - auto irt = timesview[collID]; - if (irt.toLong() < bcShift) { // due to the ROF misalignment the collision would go to negative ROF ID, discard - continue; - } - irt -= bcShift; // account for the ROF start shift - - mDigitizer.setEventTime(irt); - mDigitizer.resetEventROFrames(); // to estimate min/max ROF for this collID - // for each collision, loop over the constituents event and source IDs - // (background signal merging is basically taking place here) - for (auto& part : eventParts[collID]) { - - // get the hits for this event and this source - mHits.clear(); - context->retrieveHits(mSimChains, o2::detectors::SimTraits::DETECTORBRANCHNAMES[mID][0].c_str(), part.sourceID, part.entryID, &mHits); - - if (!mHits.empty()) { - LOG(debug) << "For collision " << collID << " eventID " << part.entryID - << " found " << mHits.size() << " hits "; - mDigitizer.process(&mHits, part.entryID, part.sourceID); // call actual digitization procedure + std::copy(mROFRecords[iLayer].begin(), mROFRecords[iLayer].end(), std::back_inserter(mROFRecordsAccum[iLayer])); + if (mWithMCTruth) { + mLabelsAccum[iLayer].mergeAtBack(mLabels[iLayer]); } + LOG(info) << "Added " << mDigits[iLayer].size() << " digits on layer " << iLayer; + mLabels[iLayer].clear(); + mDigits[iLayer].clear(); + mROFRecords[iLayer].clear(); + }; + + const int bcShift = mDigitizer.getParams().getROFrameBiasInBC(iLayer); + for (size_t collID = 0; collID < timesview.size(); ++collID) { + auto irt = timesview[collID]; + if (irt.toLong() < bcShift) { + continue; + } + irt -= bcShift; + + mDigitizer.setEventTime(irt, iLayer); + mDigitizer.resetEventROFrames(); + for (auto& part : eventParts[collID]) { + mHits.clear(); + context->retrieveHits(mSimChains, o2::detectors::SimTraits::DETECTORBRANCHNAMES[mID][0].c_str(), part.sourceID, part.entryID, &mHits); + + if (!mHits.empty()) { + LOG(debug) << "For collision " << collID << " eventID " << part.entryID + << " found " << mHits.size() << " hits on layer " << iLayer; + mDigitizer.process(&mHits, part.entryID, part.sourceID, iLayer); + } + } + if (mWithMCTruth) { + mMC2ROFRecordsAccum[iLayer].emplace_back(collID, -1, mDigitizer.getEventROFrameMin(), mDigitizer.getEventROFrameMax()); + } + accumulate(); } - mMC2ROFRecordsAccum.emplace_back(collID, -1, mDigitizer.getEventROFrameMin(), mDigitizer.getEventROFrameMax()); + mDigitizer.fillOutputContainer(0xffffffff, iLayer); accumulate(); - } - mDigitizer.fillOutputContainer(); - LOG(debug) << "mDigits size after fill: " << mDigits.size(); - accumulate(); - - // here we have all digits and labels and we can send them to consumer (aka snapshot it onto output) - std::vector expDigitRofVec(nROFsTF); - for (int iROF = 0; iROF < nROFsTF; ++iROF) { - auto& rof = expDigitRofVec[iROF]; - const int orb = iROF * roFrameLengthInBC / o2::constants::lhc::LHCMaxBunches + mFirstOrbitTF; - const int bc = iROF * roFrameLengthInBC % o2::constants::lhc::LHCMaxBunches; - rof.setBCData(o2::InteractionRecord(bc, orb)); - rof.setROFrame(iROF); - rof.setNEntries(0); - rof.setFirstEntry(-1); - } - - for (const auto& rof : mROFRecordsAccum) { - const auto& ir = rof.getBCData(); - const auto irToFirst = ir - firstIR; - const auto irROF = irToFirst.toLong() / roFrameLengthInBC; - if (irROF < 0 || irROF >= nROFsTF) { - continue; + nDigits += digitsAccum.size(); + + std::vector expDigitRofVec(nROFsTF); + for (int iROF = 0; iROF < nROFsTF; ++iROF) { + auto& rof = expDigitRofVec[iROF]; + const int orb = iROF * roFrameLengthInBC / o2::constants::lhc::LHCMaxBunches + mFirstOrbitTF; + const int bc = iROF * roFrameLengthInBC % o2::constants::lhc::LHCMaxBunches; + rof.setBCData(o2::InteractionRecord(bc, orb)); + rof.setROFrame(iROF); + rof.setNEntries(0); + rof.setFirstEntry(-1); } - auto& expROF = expDigitRofVec[irROF]; - expROF.setFirstEntry(rof.getFirstEntry()); - expROF.setNEntries(rof.getNEntries()); - if (expROF.getBCData() != rof.getBCData()) { - LOGP(fatal, "detected mismatch between expected {} and received {}", expROF.asString(), rof.asString()); + + for (const auto& rof : mROFRecordsAccum[iLayer]) { + const auto& ir = rof.getBCData(); + const auto irToFirst = ir - firstIR; + const auto irROF = irToFirst.toLong() / roFrameLengthInBC; + if (irROF < 0 || irROF >= nROFsTF) { + continue; + } + auto& expROF = expDigitRofVec[irROF]; + expROF.setFirstEntry(rof.getFirstEntry()); + expROF.setNEntries(rof.getNEntries()); + if (expROF.getBCData() != rof.getBCData()) { + LOGP(fatal, "detected mismatch between expected {} and received {}", expROF.asString(), rof.asString()); + } } - } - int prevFirst = 0; - for (auto& rof : expDigitRofVec) { - if (rof.getFirstEntry() < 0) { - rof.setFirstEntry(prevFirst); + int prevFirst = 0; + for (auto& rof : expDigitRofVec) { + if (rof.getFirstEntry() < 0) { + rof.setFirstEntry(prevFirst); + } + prevFirst = rof.getFirstEntry(); } - prevFirst = rof.getFirstEntry(); - } - pc.outputs().snapshot(Output{mOrigin, "DIGITSROF", 0}, expDigitRofVec); - if (mWithMCTruth) { - std::vector clippedMC2ROFRecords; - clippedMC2ROFRecords.reserve(mMC2ROFRecordsAccum.size()); - for (auto mc2rof : mMC2ROFRecordsAccum) { - if (mc2rof.rofRecordID < 0 || mc2rof.minROF >= static_cast(nROFsTF)) { - mc2rof.rofRecordID = -1; - mc2rof.minROF = 0; - mc2rof.maxROF = 0; - } else { - mc2rof.maxROF = std::min(mc2rof.maxROF, nROFsTF - 1); - if (mc2rof.minROF > mc2rof.maxROF) { + pc.outputs().snapshot(Output{mOrigin, "DIGITSROF", iLayer}, expDigitRofVec); + if (mWithMCTruth) { + std::vector clippedMC2ROFRecords; + clippedMC2ROFRecords.reserve(mMC2ROFRecordsAccum[iLayer].size()); + for (auto mc2rof : mMC2ROFRecordsAccum[iLayer]) { + if (mc2rof.minROF >= static_cast(nROFsTF) || mc2rof.minROF > mc2rof.maxROF) { mc2rof.rofRecordID = -1; mc2rof.minROF = 0; mc2rof.maxROF = 0; } else { - mc2rof.rofRecordID = mc2rof.minROF; + mc2rof.maxROF = std::min(mc2rof.maxROF, nROFsTF - 1); + if (mc2rof.minROF > mc2rof.maxROF) { + mc2rof.rofRecordID = -1; + mc2rof.minROF = 0; + mc2rof.maxROF = 0; + } else { + mc2rof.rofRecordID = mc2rof.minROF; + } } + clippedMC2ROFRecords.push_back(mc2rof); } - clippedMC2ROFRecords.push_back(mc2rof); + pc.outputs().snapshot(Output{mOrigin, "DIGITSMC2ROF", iLayer}, clippedMC2ROFRecords); + auto& sharedlabels = pc.outputs().make>(Output{mOrigin, "DIGITSMCTR", iLayer}); + mLabelsAccum[iLayer].flatten_to(sharedlabels); + mLabels[iLayer].clear_andfreememory(); + mLabelsAccum[iLayer].clear_andfreememory(); } - pc.outputs().snapshot(Output{mOrigin, "DIGITSMC2ROF", 0}, clippedMC2ROFRecords); - auto& sharedlabels = pc.outputs().make>(Output{mOrigin, "DIGITSMCTR", 0}); - mLabelsAccum.flatten_to(sharedlabels); - // free space of existing label containers - mLabels.clear_andfreememory(); - mLabelsAccum.clear_andfreememory(); } 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"; + LOG(info) << "Produced " << nDigits << " digits"; // we should be only called once; tell DPL that this process is ready to exit pc.services().get().readyToQuit(QuitRequest::Me); @@ -288,13 +279,25 @@ class TRKDPLDigitizerTask : BaseDPLDigitizer const auto& dopt = o2::trk::DPLDigitizerParam::Instance(); // pc.inputs().get("TRK_almiraparam"); const auto& aopt = o2::trk::AlmiraParam::Instance(); - auto frameNS = aopt.roFrameLengthInBC * o2::constants::lhc::LHCBunchSpacingNS; - digipar.setContinuous(true); - digipar.setROFrameBiasInBC(aopt.roFrameBiasInBC); - digipar.setROFrameLengthInBC(aopt.roFrameLengthInBC); - digipar.setROFrameLength(frameNS); // RO frame in ns - digipar.setStrobeDelay(aopt.strobeDelay); - digipar.setStrobeLength(aopt.strobeLengthCont > 0 ? aopt.strobeLengthCont : frameNS - aopt.strobeDelay); + mLayers = constants::VD::petal::nLayers + geom->getNumberOfLayersMLOT(); + mDigits.resize(mLayers); + mROFRecords.resize(mLayers); + mROFRecordsAccum.resize(mLayers); + mLabels.resize(mLayers); + mLabelsAccum.resize(mLayers); + mMC2ROFRecordsAccum.resize(mLayers); + + for (int iLayer = 0; iLayer < mLayers; ++iLayer) { + const auto roFrameLengthInBC = aopt.getROFLengthInBC(iLayer); + const auto frameNS = roFrameLengthInBC * o2::constants::lhc::LHCBunchSpacingNS; + digipar.setROFrameLengthInBC(roFrameLengthInBC, iLayer); + // ROF delay is treated as an additional bias from the digitizer point of view. + digipar.setROFrameBiasInBC(aopt.getROFBiasInBC(iLayer) + aopt.getROFDelayInBC(iLayer), iLayer); + digipar.setStrobeDelay(aopt.getStrobeDelay(iLayer), iLayer); + const auto strobeLengthCont = aopt.getStrobeLengthCont(iLayer); + digipar.setStrobeLength(strobeLengthCont > 0 ? strobeLengthCont : frameNS - aopt.getStrobeDelay(iLayer), iLayer); + digipar.setROFrameLength(frameNS, iLayer); + } // parameters of signal time response: flat-top duration, max rise time and q @ which rise time is 0 digipar.getSignalShape().setParameters(dopt.strobeFlatTop, dopt.strobeMaxRiseTime, dopt.strobeQRiseTime0); digipar.setChargeThreshold(dopt.chargeThreshold); // charge threshold in electrons @@ -352,17 +355,16 @@ class TRKDPLDigitizerTask : BaseDPLDigitizer const o2::detectors::DetID mID{o2::detectors::DetID::TRK}; const o2::header::DataOrigin mOrigin{o2::header::gDataOriginTRK}; o2::trk::Digitizer mDigitizer{}; - std::vector mDigits{}; - std::vector mROFRecords{}; - std::vector mROFRecordsAccum{}; + int mLayers{0}; + std::vector> mDigits{}; + std::vector> mROFRecords{}; + std::vector> mROFRecordsAccum{}; std::vector mHits{}; std::vector* mHitsP{&mHits}; - o2::dataformats::MCTruthContainer mLabels{}; - o2::dataformats::MCTruthContainer mLabelsAccum{}; - std::vector mMC2ROFRecordsAccum{}; + std::vector> mLabels{}; + std::vector> mLabelsAccum{}; + std::vector> mMC2ROFRecordsAccum{}; std::vector mSimChains{}; - - int mFixMC2ROF = 0; // 1st entry in mc2rofRecordsAccum to be fixed for ROFRecordID o2::parameters::GRPObject::ROMode mROMode = o2::parameters::GRPObject::PRESENT; // readout mode }; From c5b2f4f9ed6fe258f8d32f5a40dc6c8c7e982ea2 Mon Sep 17 00:00:00 2001 From: Maximiliano Puccio Date: Thu, 23 Apr 2026 23:11:45 +0200 Subject: [PATCH 3/6] ALICE3: Adapt test macros to new staggering structure --- .../ALICE3/TRK/macros/test/CheckBandwidth.C | 205 +++++---- .../ALICE3/TRK/macros/test/CheckDigitsTRK.C | 429 ++++++++++++++++++ 2 files changed, 544 insertions(+), 90 deletions(-) create mode 100644 Detectors/Upgrades/ALICE3/TRK/macros/test/CheckDigitsTRK.C diff --git a/Detectors/Upgrades/ALICE3/TRK/macros/test/CheckBandwidth.C b/Detectors/Upgrades/ALICE3/TRK/macros/test/CheckBandwidth.C index 06d24361c7721..cda07447d58f5 100644 --- a/Detectors/Upgrades/ALICE3/TRK/macros/test/CheckBandwidth.C +++ b/Detectors/Upgrades/ALICE3/TRK/macros/test/CheckBandwidth.C @@ -41,13 +41,9 @@ namespace { constexpr double DigitBits = 24.; -constexpr double BunchCrossingNS = 25.; -constexpr int ReadoutCycleBC = 18; -constexpr int ReadoutCycleSimBC = 18; -constexpr double ReadoutCycleSeconds = ReadoutCycleBC * BunchCrossingNS * 1.e-9; } // namespace -void CheckBandwidth(std::string digifile = "trkdigits.root", std::string inputGeom = "o2sim_geometry.root", std::string collContextFile = "collisioncontext.root") +void CheckBandwidth(std::string digifile = "trkdigits.root", std::string inputGeom = "sgn_geometry.root", std::string collContextFile = "collisioncontext.root") { gStyle->SetPalette(55); gStyle->SetOptStat(0); @@ -73,15 +69,14 @@ void CheckBandwidth(std::string digifile = "trkdigits.root", std::string inputGe latex.DrawLatex(0.04, 0.06, Form("avg non-empty: %.3f collisions/ROF", nonEmptyAverageValue)); }; - auto drawCollisionInfoBox = [](double averageValue) { - const double effectiveIRRateHz = ReadoutCycleSeconds > 0. ? averageValue / ReadoutCycleSeconds : 0.; + auto drawCollisionInfoBox = [](double effectiveIRRateHz, double rofLengthBC) { TPaveText infoBox(0.55, 0.79, 0.88, 0.9, "NDC"); infoBox.SetFillColor(0); infoBox.SetBorderSize(1); infoBox.SetTextAlign(12); infoBox.SetTextSize(0.028); infoBox.AddText(Form("effective IR: %.3f MHz", effectiveIRRateHz * 1.e-6)); - infoBox.AddText(Form("ROF length: %d BC", ReadoutCycleBC)); + infoBox.AddText(Form("ROF length: %d BC", rofLengthBC)); infoBox.DrawClone(); }; @@ -168,13 +163,43 @@ void CheckBandwidth(std::string digifile = "trkdigits.root", std::string inputGe } } + // --- Digits --- + + TFile* digFile = TFile::Open(digifile.data()); + TTree* digTree = (TTree*)digFile->Get("o2sim"); + const int nDigitTreeEntries = digTree->GetEntries(); + + std::vector*> digArr(nTotalLayers, nullptr); + std::vector*> rofRecords(nTotalLayers, nullptr); + for (int nDigitsLayer{0}; nDigitsLayer < nTotalLayers; ++nDigitsLayer) { + if (!digTree->GetBranch(Form("TRKDigit_%i", nDigitsLayer))) { + break; + } + digTree->SetBranchAddress(Form("TRKDigit_%i", nDigitsLayer), &digArr[nDigitsLayer]); + digTree->SetBranchAddress(Form("TRKDigitROF_%i", nDigitsLayer), &rofRecords[nDigitsLayer]); + } + + digTree->GetEntry(0); + if (nDigitTreeEntries > 1) { + LOG(warning) << "Digit tree has " << nDigitTreeEntries << " entries, but this macro processes entry 0 only."; + } + + std::vector rofLengthBC(nTotalLayers, 0u); + for (int iLayer = 0; iLayer < nTotalLayers; ++iLayer) { + if (rofRecords[iLayer]->size() < 2) { + LOG(fatal) << "ROF record tree for layer " << iLayer << " has " << rofRecords[iLayer]->size() + << " entries, but at least 2 are expected (one per ROF + one empty at the end). Check input files."; + } + rofLengthBC[iLayer] = (*rofRecords[iLayer])[1].getBCData().bc - (*rofRecords[iLayer])[0].getBCData().bc; + } + + // --- Collision context --- TFile* ccFile = TFile::Open(collContextFile.data()); auto* digiContext = (o2::steer::DigitizationContext*)ccFile->Get("DigitizationContext"); const o2::InteractionRecord firstSampledIR{0, digiContext->getFirstOrbitForSampling()}; - std::vector collisionsPerROF; - + std::vector> collisionsPerROF(nTotalLayers); for (const auto& record : digiContext->getEventRecords()) { auto nbc = record.differenceInBC(firstSampledIR); if (record.getTimeOffsetWrtBC() < 0. && nbc > 0) { @@ -183,63 +208,48 @@ void CheckBandwidth(std::string digifile = "trkdigits.root", std::string inputGe if (nbc < 0) { continue; } - const size_t rofID = nbc / ReadoutCycleSimBC; - if (rofID >= collisionsPerROF.size()) { - collisionsPerROF.resize(rofID + 1, 0u); + for (int iLayer = 0; iLayer < nTotalLayers; ++iLayer) { + if (rofLengthBC[iLayer] == 0) { + LOG(fatal) << "ROF length in BC for layer " << iLayer << " is zero. Check input files."; + } + const size_t rofID = nbc / rofLengthBC[iLayer]; + if (rofID >= collisionsPerROF[iLayer].size()) { + collisionsPerROF[iLayer].resize(rofID + 1, 0u); + } + ++collisionsPerROF[iLayer][rofID]; } - ++collisionsPerROF[rofID]; } - // --- Digits --- - - TFile* digFile = TFile::Open(digifile.data()); - TTree* digTree = (TTree*)digFile->Get("o2sim"); - const int nDigitTreeEntries = digTree->GetEntries(); - - std::vector* digArr = nullptr; - std::vector* rofRecords = nullptr; - digTree->SetBranchAddress("TRKDigit", &digArr); - digTree->SetBranchAddress("TRKDigitROF", &rofRecords); - - digTree->GetEntry(0); - if (nDigitTreeEntries > 1) { - LOG(warning) << "Digit tree has " << nDigitTreeEntries << " entries, but this macro processes entry 0 only."; - } - - const int nROFRec = (int)rofRecords->size(); - if (nROFRec != (int)collisionsPerROF.size()) { + const int nROFRec = (int)rofRecords[0]->size(); + if (nROFRec != (int)collisionsPerROF[0].size()) { LOG(fatal) << "Mismatch between number of ROF records in digit tree (" << nROFRec - << ") and number of ROFs computed from collisioncontext.root (" << collisionsPerROF.size() + << ") and number of ROFs computed from collisioncontext.root (" << collisionsPerROF[0].size() << "). Check input files."; } // --- Accumulate per-chip digit counts across all ROFs --- - const double rofNorm = nROFRec > 0 ? 1. / nROFRec : 0.; - const double bitsToGbps = ReadoutCycleSeconds > 0. ? DigitBits / ReadoutCycleSeconds / 1.e9 : 0.; std::vector digitsPerChip(nChips, 0ull); std::vector maxDigitsPerROFPerChip(nChips, 0u); std::vector digitsInCurrentROFPerChip(nChips, 0u); - for (unsigned int iROF = 0; iROF < rofRecords->size(); ++iROF) { + for (unsigned int iROF = 0; iROF < (unsigned int)nROFRec; ++iROF) { std::vector touchedChips; - const unsigned int rofStart = (*rofRecords)[iROF].getFirstEntry(); - const unsigned int rofEnd = rofStart + (*rofRecords)[iROF].getNEntries(); - - for (unsigned int iDigit = rofStart; iDigit < rofEnd; ++iDigit) { - if (iDigit % 1000 == 0) { - std::cout << "Reading digit " << iDigit << " / " << digArr->size() << "\r" << std::flush; - } - const int chipID = (*digArr)[iDigit].getChipIndex(); - if (chipGeom[chipID].disk != -1) { - continue; - } - if (digitsInCurrentROFPerChip[chipID] == 0) { - touchedChips.push_back(chipID); + for (int iLayer = 0; iLayer < nTotalLayers; ++iLayer) { + const unsigned int rofStart = (*rofRecords[iLayer])[iROF].getFirstEntry(); + const unsigned int rofEnd = rofStart + (*rofRecords[iLayer])[iROF].getNEntries(); + for (unsigned int iDigit = rofStart; iDigit < rofEnd; ++iDigit) { + const int chipID = (*digArr[iLayer])[iDigit].getChipIndex(); + if (chipGeom[chipID].disk != -1) { + continue; + } + if (digitsInCurrentROFPerChip[chipID] == 0) { + touchedChips.push_back(chipID); + } + ++digitsPerChip[chipID]; + ++digitsInCurrentROFPerChip[chipID]; } - ++digitsPerChip[chipID]; - ++digitsInCurrentROFPerChip[chipID]; } for (const int chipID : touchedChips) { @@ -273,19 +283,21 @@ void CheckBandwidth(std::string digifile = "trkdigits.root", std::string inputGe // digitsInCurrentROFPerChip is all zeros after the first scan — reuse it here. { std::vector touchedChips; - for (unsigned int iROF = 0; iROF < rofRecords->size(); ++iROF) { + for (unsigned int iROF = 0; iROF < (unsigned int)nROFRec; ++iROF) { touchedChips.clear(); - const unsigned int rofStart = (*rofRecords)[iROF].getFirstEntry(); - const unsigned int rofEnd = rofStart + (*rofRecords)[iROF].getNEntries(); - for (unsigned int iDigit = rofStart; iDigit < rofEnd; ++iDigit) { - const int chipID = (*digArr)[iDigit].getChipIndex(); - if (chipGeom[chipID].disk != -1) { - continue; + for (int iLayer = 0; iLayer < nTotalLayers; ++iLayer) { + const unsigned int rofStart = (*rofRecords[iLayer])[iROF].getFirstEntry(); + const unsigned int rofEnd = rofStart + (*rofRecords[iLayer])[iROF].getNEntries(); + for (unsigned int iDigit = rofStart; iDigit < rofEnd; ++iDigit) { + const int chipID = (*digArr[iLayer])[iDigit].getChipIndex(); + if (chipGeom[chipID].disk != -1) { + continue; + } + if (digitsInCurrentROFPerChip[chipID] == 0) { + touchedChips.push_back(chipID); + } + ++digitsInCurrentROFPerChip[chipID]; } - if (digitsInCurrentROFPerChip[chipID] == 0) { - touchedChips.push_back(chipID); - } - ++digitsInCurrentROFPerChip[chipID]; } for (const int chipID : touchedChips) { const int l = chipGeom[chipID].globalLayer; @@ -321,7 +333,7 @@ void CheckBandwidth(std::string digifile = "trkdigits.root", std::string inputGe continue; } const int l = g.globalLayer; - const double avgDigits = digitsPerChip[chipID] * rofNorm; + const double avgDigits = digitsPerChip[chipID] / collisionsPerROF[l].size(); const double maxDigits = (double)maxDigitsPerROFPerChip[chipID]; layerStats[l].avgDigitsPerROF += avgDigits; layerStats[l].avgMaxDigitsPerROF += maxDigits; @@ -334,37 +346,50 @@ void CheckBandwidth(std::string digifile = "trkdigits.root", std::string inputGe layerStats[l].avgDigitsPerROF *= norm; layerStats[l].avgMaxDigitsPerROF *= norm; } - layerStats[l].avgBandwidthGbps = layerStats[l].avgDigitsPerROF * bitsToGbps; - layerStats[l].peakBandwidthGbps = layerStats[l].peakAvgDigitsPerROF * bitsToGbps; + layerStats[l].avgBandwidthGbps = layerStats[l].avgDigitsPerROF * DigitBits / rofLengthBC[l] / o2::constants::lhc::LHCBunchSpacingNS * 1.e9; + layerStats[l].peakBandwidthGbps = layerStats[l].peakAvgDigitsPerROF * DigitBits / rofLengthBC[l] / o2::constants::lhc::LHCBunchSpacingNS * 1.e9; } // --- Collision plots --- if (nROFRec > 0) { - auto* hCollisionsPerROF = new TH1D("h_collisions_per_rof", "Collisions per ROF;ROF id;N collisions", - nROFRec, -0.5, nROFRec - 0.5); - double totalCollisionsPerROF = 0.; - double peakCollisionsPerROF = 0.; - int nNonEmptyROFs = 0; - - for (int rofID = 0; rofID < nROFRec; ++rofID) { - const double nColl = collisionsPerROF[rofID]; - hCollisionsPerROF->SetBinContent(rofID + 1, nColl); - totalCollisionsPerROF += nColl; - peakCollisionsPerROF = std::max(peakCollisionsPerROF, nColl); - if (nColl > 0.) { - ++nNonEmptyROFs; + std::vector totalCollisionsPerROF(nTotalLayers, 0.); + std::vector peakCollisionsPerROF(nTotalLayers, 0.); + std::vector nNonEmptyROFs(nTotalLayers, 0); + std::vector hCollisionsPerROFPerLayer(nTotalLayers, nullptr); + + for (int iLayer = 0; iLayer < nTotalLayers; ++iLayer) { + hCollisionsPerROFPerLayer[iLayer] = new TH1D(Form("h_collisions_per_rof_layer%d", iLayer), + Form("Layer %d;ROF id;N collisions", iLayer), + nROFRec, -0.5, nROFRec - 0.5); + for (int rofID = 0; rofID < nROFRec; ++rofID) { + const double nColl = collisionsPerROF[iLayer][rofID]; + hCollisionsPerROFPerLayer[iLayer]->SetBinContent(rofID + 1, nColl); + totalCollisionsPerROF[iLayer] += nColl; + peakCollisionsPerROF[iLayer] = std::max(peakCollisionsPerROF[iLayer], nColl); + if (nColl > 0.) { + ++nNonEmptyROFs[iLayer]; + } } } - const double avgCollisionsPerROF = totalCollisionsPerROF / nROFRec; - auto* canvCollisionsPerROF = new TCanvas("canvCollisionsPerROF", "Collisions per ROF", 1050, 1050); - canvCollisionsPerROF->SetTopMargin(0.08); - hCollisionsPerROF->Draw("hist"); - drawCollisionSummary(avgCollisionsPerROF, - nNonEmptyROFs > 0 ? totalCollisionsPerROF / nNonEmptyROFs : 0., - peakCollisionsPerROF); - drawCollisionInfoBox(avgCollisionsPerROF); + const int nCols = std::max(1, (int)std::ceil(std::sqrt((double)nTotalLayers))); + const int nRows = (nTotalLayers + nCols - 1) / nCols; + auto* canvCollisionsPerROF = new TCanvas("canvCollisionsPerROF", "Collisions per ROF", 350 * nCols, 300 * nRows); + canvCollisionsPerROF->Divide(nCols, nRows); + for (int iLayer = 0; iLayer < nTotalLayers; ++iLayer) { + canvCollisionsPerROF->cd(iLayer + 1); + gPad->SetTopMargin(0.10); + gPad->SetBottomMargin(0.14); + gPad->SetLeftMargin(0.14); + hCollisionsPerROFPerLayer[iLayer]->Draw("hist"); + const double avgCollisionsPerROF = totalCollisionsPerROF[iLayer] / collisionsPerROF[iLayer].size(); + drawCollisionSummary(avgCollisionsPerROF, + nNonEmptyROFs[iLayer] > 0 ? totalCollisionsPerROF[iLayer] / nNonEmptyROFs[iLayer] : 0., + peakCollisionsPerROF[iLayer]); + const double effectiveIRRateHz = avgCollisionsPerROF / rofLengthBC[iLayer] / o2::constants::lhc::LHCBunchSpacingNS * 1.e9; + drawCollisionInfoBox(effectiveIRRateHz, rofLengthBC[iLayer]); + } appendCanvasToPdf(canvCollisionsPerROF); } @@ -405,9 +430,9 @@ void CheckBandwidth(std::string digifile = "trkdigits.root", std::string inputGe if (g.petal < 0 || g.petal >= nVDPetals) { continue; } - const double avgDigits = digitsPerChip[chipID] * rofNorm; + const double avgDigits = double(digitsPerChip[chipID]) / collisionsPerROF[g.globalLayer].size(); const double maxDigits = (double)maxDigitsPerROFPerChip[chipID]; - const double bandwidth = avgDigits * bitsToGbps; + const double bandwidth = avgDigits * DigitBits / rofLengthBC[g.globalLayer] / o2::constants::lhc::LHCBunchSpacingNS * 1.e9; hVDDigitsPerROF->SetBinContent(g.petal + 1, g.localLayer + 1, avgDigits); hVDMaxDigitsPerROF->SetBinContent(g.petal + 1, g.localLayer + 1, maxDigits); @@ -483,12 +508,12 @@ void CheckBandwidth(std::string digifile = "trkdigits.root", std::string inputGe continue; } const double staveBinX = g.stave + (g.halfStave + 0.5) / nHalfStaves - 0.5; - const double avgDigits = digitsPerChip[chipID] * rofNorm; + const double avgDigits = double(digitsPerChip[chipID]) / collisionsPerROF[g.globalLayer].size(); const double maxDigits = (double)maxDigitsPerROFPerChip[chipID]; hDigitsPerROF->Fill(staveBinX, sensorID, avgDigits); hMaxDigitsPerROF->Fill(staveBinX, sensorID, maxDigits); - hBandwidth->Fill(staveBinX, sensorID, avgDigits * bitsToGbps); + hBandwidth->Fill(staveBinX, sensorID, avgDigits * DigitBits / rofLengthBC[g.globalLayer] / o2::constants::lhc::LHCBunchSpacingNS * 1.e9); } const auto& ls = layerStats[outputLayer]; diff --git a/Detectors/Upgrades/ALICE3/TRK/macros/test/CheckDigitsTRK.C b/Detectors/Upgrades/ALICE3/TRK/macros/test/CheckDigitsTRK.C new file mode 100644 index 0000000000000..f29a38f554200 --- /dev/null +++ b/Detectors/Upgrades/ALICE3/TRK/macros/test/CheckDigitsTRK.C @@ -0,0 +1,429 @@ +// 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 CheckDigits.C +/// \brief Simple macro to check TRK digits + +#if !defined(__CLING__) || defined(__ROOTCLING__) +#include +#include +#include +#include +#include +#include +#include +#include + +#include "TRKBase/SegmentationChip.h" +#include "TRKBase/GeometryTGeo.h" +#include "DataFormatsITSMFT/Digit.h" +#include "TRKSimulation/Hit.h" +#include "MathUtils/Utils.h" +#include "SimulationDataFormat/ConstMCTruthContainer.h" +#include "SimulationDataFormat/IOMCTruthContainerView.h" +#include "SimulationDataFormat/MCCompLabel.h" +#include "DetectorsBase/GeometryManager.h" +#include "ITSMFTSimulation/AlpideSimResponse.h" +#include "CCDB/BasicCCDBManager.h" + +#include "DataFormatsITSMFT/ROFRecord.h" + +#endif + +#define ENABLE_UPGRADES + +void addTLines(float pitch) +{ + // Add grid lines at multiples of pitch on the current pad + if (!gPad) + return; + + gPad->Update(); + + Double_t xmin = gPad->GetUxmin(); + Double_t xmax = gPad->GetUxmax(); + Double_t ymin = gPad->GetUymin(); + Double_t ymax = gPad->GetUymax(); + + // Calculate the first vertical line position (multiple of pitch) + int nLinesX = 0; + for (float x = xmin; x <= xmax && nLinesX < 1000; x += pitch, nLinesX++) { + TLine* line = new TLine(x, ymin, x, ymax); + line->SetLineStyle(2); + line->SetLineColor(kGray); + line->Draw("same"); + } + + // Calculate the first horizontal line position (multiple of pitch) + int nLinesY = 0; + for (float y = ymin; y <= ymax && nLinesY < 1000; y += pitch, nLinesY++) { + TLine* line = new TLine(xmin, y, xmax, y); + line->SetLineStyle(2); + line->SetLineColor(kGray); + line->Draw("same"); + } + + gPad->Modified(); + gPad->Update(); +} + +void CheckDigits(std::string digifile = "trkdigits.root", std::string hitfile = "o2sim_HitsTRK.root", std::string inputGeom = "o2sim_geometry.root") +{ + gStyle->SetPalette(55); + + using namespace o2::base; + using namespace o2::trk; + + using o2::itsmft::Digit; + using o2::trk::Hit; + + using o2::trk::SegmentationChip; + + TFile* f = TFile::Open("CheckDigits.root", "recreate"); + + TNtuple* nt = new TNtuple("ntd", "digit ntuple", "id:x:y:z:rowD:colD:rowH:colH:xlH:zlH:xlcH:zlcH:dx:dz"); + TNtuple* nt2 = new TNtuple("ntd2", "digit ntuple", "id:z:dxH:dzH"); /// maximum number of elements in a tuple = 15: doing a new tuple to store more variables + + // Geometry + o2::base::GeometryManager::loadGeometry(inputGeom); + auto* gman = o2::trk::GeometryTGeo::Instance(); + gman->fillMatrixCache(o2::math_utils::bit2Mask(o2::math_utils::TransformType::L2G)); + + const int nVDLayers = gman->extractNumberOfLayersVD(); + const int nMLOTLayers = gman->getNumberOfLayersMLOT(); + const int nTotalLayers = nVDLayers + nMLOTLayers; + + SegmentationChip seg; + // seg.Print(); + + // MLOT response plane: y = halfThickness - depthMax. + float depthMax = (float)o2::trk::constants::apts::thickness; // fallback (no CCDB) + auto& ccdbMgr = o2::ccdb::BasicCCDBManager::instance(); + ccdbMgr.setURL("http://alice-ccdb.cern.ch"); + if (auto* alpResp = ccdbMgr.get("IT3/Calib/APTSResponse")) { + depthMax = alpResp->getDepthMax(); + } + const float yPlaneMLOT = o2::trk::SegmentationChip::SiliconThicknessMLOT / 2.f - depthMax; + const float yPlaneVD = -o2::trk::SegmentationChip::SiliconThicknessVD; // VD reference plane in local flat y + // 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); + + // Digits — per-layer branches + TFile* digFile = TFile::Open(digifile.data()); + TTree* digTree = (TTree*)digFile->Get("o2sim"); + + int nDigitLayers = 0; + std::vector*> digArr(nTotalLayers, nullptr); + std::vector*> rofRecordsArr(nTotalLayers, nullptr); + std::vector plabelsArr(nTotalLayers, nullptr); + + for (int iLayer = 0; iLayer < nTotalLayers; ++iLayer) { + if (!digTree->GetBranch(Form("TRKDigit_%i", iLayer))) { + break; + } + digTree->SetBranchAddress(Form("TRKDigit_%i", iLayer), &digArr[iLayer]); + digTree->SetBranchAddress(Form("TRKDigitROF_%i", iLayer), &rofRecordsArr[iLayer]); + digTree->SetBranchAddress(Form("TRKDigitMCTruth_%i", iLayer), &plabelsArr[iLayer]); + ++nDigitLayers; + } + + digTree->GetEntry(0); + + // Load all MC hit events upfront and build the hit lookup map. + for (int im = 0; im < nevH; ++im) { + hitTree->SetBranchAddress("TRKHit", &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); + } + } + + // LOOP over layers, then ROFRecords within each layer + for (int iLayer = 0; iLayer < nDigitLayers; ++iLayer) { + auto& rofArr = *rofRecordsArr[iLayer]; + const int nROFRec = (int)rofArr.size(); + + o2::dataformats::ConstMCTruthContainer labels; + plabelsArr[iLayer]->copyandflatten(labels); + + // 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(); + + // LOOP on : digits array + for (unsigned int iDigit = rofIndex; iDigit < rofIndex + rofNEntries; iDigit++) { + if (iDigit % 1000 == 0) + std::cout << "Layer " << iLayer << ": reading digit " << iDigit << " / " << digArr[iLayer]->size() << std::endl; + + Int_t ix = (*digArr[iLayer])[iDigit].getRow(), iz = (*digArr[iLayer])[iDigit].getColumn(); + Int_t iDetID = (*digArr[iLayer])[iDigit].getChipIndex(); + Int_t layer = gman->getLayer(iDetID); + Int_t disk = gman->getDisk(iDetID); + Int_t subDetID = gman->getSubDetID(iDetID); + Int_t petalCase = gman->getPetalCase(iDetID); + Int_t stave = gman->getStave(iDetID); + Int_t halfstave = gman->getHalfStave(iDetID); + + Float_t x = 0.f, y = 0.f, z = 0.f; + Float_t x_flat = 0.f, z_flat = 0.f; + + if (disk != -1) { + continue; // skip disks for the moment + } + + if (subDetID != 0) { + seg.detectorToLocal(ix, iz, x, z, subDetID, layer, disk); + } else if (subDetID == 0) { + seg.detectorToLocal(ix, iz, x_flat, z_flat, subDetID, layer, disk); + o2::math_utils::Vector2D xyCurved = seg.flatToCurved(layer, x_flat, 0.); + x = xyCurved.X(); + y = xyCurved.Y(); + z = z_flat; + } + + o2::math_utils::Point3D locD(x, y, z); // local Digit curved + o2::math_utils::Point3D locDF(-1, -1, -1); // local Digit flat + + Int_t chipID = (*digArr[iLayer])[iDigit].getChipIndex(); + auto lab = (labels.getLabels(iDigit))[0]; + + int trID = lab.getTrackID(); + + if (!lab.isValid()) { // not a noise + continue; + } + + const auto gloD = gman->getMatrixL2G(chipID)(locD); // convert to global + + std::unordered_map* mc2hit = &mc2hitVec[lab.getEventID()]; + + // get MC info + 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; + continue; + } + + ////// HITS + Hit& hit = (*hitArray[lab.getEventID()])[hitEntry->second]; + + auto xyzLocE = gman->getMatrixL2G(chipID) ^ (hit.GetPos()); // inverse conversion from global to local + auto xyzLocS = gman->getMatrixL2G(chipID) ^ (hit.GetPosStart()); + + // Hit local reference: Both VD and MLOT use response-plane interpolation (in flat local frame). + // For VD, transform curved → flat first, then interpolate. + o2::math_utils::Vector3D locH; /// Hit reference (at response plane) + o2::math_utils::Vector3D locHS; /// Hit, start pos + locHS.SetCoordinates(xyzLocS.X(), xyzLocS.Y(), xyzLocS.Z()); + o2::math_utils::Vector3D locHE; /// Hit, end pos + locHE.SetCoordinates(xyzLocE.X(), xyzLocE.Y(), xyzLocE.Z()); + o2::math_utils::Vector3D locHF; + + if (subDetID == 0) { + // VD: Interpolate to VD reference plane in flat frame; apply same r to X and Z + auto flatSta = seg.curvedToFlat(layer, locHS.X(), locHS.Y()); + auto flatEnd = seg.curvedToFlat(layer, locHE.X(), locHE.Y()); + float x0 = flatSta.X(), y0 = flatSta.Y(), z0 = locHS.Z(); + float dltx = flatEnd.X() - x0, dlty = flatEnd.Y() - y0, dltz = locHE.Z() - z0; + float r = (std::abs(dlty) > 1e-9f) ? (yPlaneVD - y0) / dlty : 0.5f; + locH.SetCoordinates(x0 + r * dltx, yPlaneVD, z0 + r * dltz); + } else { + // MLOT: Interpolate to response plane + float x0 = locHS.X(), y0 = locHS.Y(), z0 = locHS.Z(); + float dltx = locHE.X() - x0, dlty = locHE.Y() - y0, dltz = locHE.Z() - z0; + float r = (std::abs(dlty) > 1e-9f) ? (yPlaneMLOT - y0) / dlty : 0.5f; + locH.SetCoordinates(x0 + r * dltx, yPlaneMLOT, z0 + r * dltz); + } + + int row = 0, col = 0; + float xlc = 0., zlc = 0.; + + if (subDetID == 0) { + Float_t x_flat = 0.f, y_flat = 0.f; + // locH is already in flat frame from interpolation above; convert digit to flat for comparison + o2::math_utils::Vector2D xyFlatD = seg.curvedToFlat(layer, locD.X(), locD.Y()); + locDF.SetCoordinates(xyFlatD.X(), xyFlatD.Y(), locD.Z()); + locHF.SetCoordinates(locH.X(), locH.Y(), locH.Z()); // locH already in flat frame + seg.localToDetector(locHF.X(), locHF.Z(), row, col, subDetID, layer, disk); + } else { + seg.localToDetector(locH.X(), locH.Z(), row, col, subDetID, layer, disk); + } + + seg.detectorToLocal(row, col, xlc, zlc, subDetID, layer, disk); + + if (subDetID == 0) { + nt->Fill(chipID, /// detector ID + gloD.X(), gloD.Y(), gloD.Z(), /// global position retrieved from the digit: digit (row, col) ->local position -> global potision + ix, iz, /// row and column of the digit + row, col, /// row and col retrieved from the hit: hit global position -> hit local position -> detector position (row, col) + locH.X(), locH.Z(), /// x and z of the hit in the local reference frame: hit global position -> hit local position + xlc, zlc, /// x and z of the hit in the local frame: hit global position -> hit local position -> detector position (row, col) -> local position + locHF.X() - locDF.X(), locHF.Z() - locDF.Z()); /// difference in x and z between the hit and the digit in the local frame + nt2->Fill(chipID, gloD.Z(), locHS.X() - locHE.X(), locHS.Z() - locHE.Z()); /// differences between local hit start and hit end positions + } else { + nt->Fill(chipID, /// detector ID + gloD.X(), gloD.Y(), gloD.Z(), /// global position retrieved from the digit: digit (row, col) ->local position -> global potision + ix, iz, /// row and column of the digit + row, col, /// row and col retrieved from the hit: hit global position -> hit local position -> detector position (row, col) + locH.X(), locH.Z(), /// x and z of the hit in the local reference frame: hit global position -> hit local position + xlc, zlc, /// x and z of the hit in the local frame: hit global position -> hit local position -> detector position (row, col) -> local position + locH.X() - locD.X(), locH.Z() - locD.Z()); /// difference in x and z between the hit and the digit in the local frame + nt2->Fill(chipID, gloD.Z(), locHS.X() - locHE.X(), locHS.Z() - locHE.Z()); /// differences between local hit start and hit end positions + } + + } // end loop on digits array + + } // end loop on ROFRecords + + } // end loop on layers + + // digit maps in the xy and yz planes + auto canvXY = new TCanvas("canvXY", "", 1600, 2400); + canvXY->Divide(2, 3); + canvXY->cd(1); + nt->Draw("y:x >>h_y_vs_x_VD(1000, -3, 3, 1000, -3, 3)", "id < 12 ", "colz"); + canvXY->cd(2); + nt->Draw("y:z>>h_y_vs_z_VD(1000, -26, 26, 1000, -3, 3)", "id < 12 ", "colz"); + canvXY->cd(3); + nt->Draw("y:x>>h_y_vs_x_ML(1000, -25, 25, 1000, -25, 25)", "id >= 12 && id < 5132 ", "colz"); + canvXY->cd(4); + nt->Draw("y:z>>h_y_vs_z_ML(1000, -70, 70, 1000, -25, 25)", "id >= 12 && id < 5132 ", "colz"); + canvXY->cd(5); + nt->Draw("y:x>>h_y_vs_x_OT(1000, -85, 85, 1000, -85, 85)", "id >= 5132 ", "colz"); + canvXY->cd(6); + nt->Draw("y:z>>h_y_vs_z_OT(1000, -85, 85, 1000, -130, 130)", "id >= 5132 ", "colz"); + canvXY->SaveAs("trkdigits_y_vs_x_vs_z.pdf"); + + // z distributions + auto canvZ = new TCanvas("canvZ", "", 800, 2400); + canvZ->Divide(1, 3); + canvZ->cd(1); + nt->Draw("z>>h_z_VD(500, -26, 26)", "id < 12 "); + canvZ->cd(2); + nt->Draw("z>>h_z_ML(500, -70, 70)", "id >= 12 && id < 5132 "); + canvZ->cd(3); + nt->Draw("z>>h_z_OT(500, -85, 85)", "id >= 5132 "); + canvZ->SaveAs("trkdigits_z.pdf"); + + // dz distributions (difference between local position of digits and hits in x and z) + auto canvdZ = new TCanvas("canvdZ", "", 800, 2400); + canvdZ->Divide(1, 3); + canvdZ->cd(1); + nt->Draw("dz>>h_dz_VD(500, -0.05, 0.05)", "id < 12 "); + canvdZ->cd(2); + nt->Draw("dz>>h_dz_ML(500, -0.05, 0.05)", "id >= 12 && id < 5132 "); + canvdZ->cd(3); + nt->Draw("dz>>h_dz_OT(500, -0.05, 0.05)", "id >= 5132 "); + canvdZ->SaveAs("trkdigits_dz.pdf"); + canvdZ->SaveAs("trkdigits_dz.root"); + + // distributions of differences between local positions of digits and hits in x and z + auto canvdXdZ = new TCanvas("canvdXdZ", "", 1600, 2400); + canvdXdZ->Divide(2, 3); + canvdXdZ->cd(1); + nt->Draw("dx:dz>>h_dx_vs_dz_VD(500, -0.005, 0.005, 500, -0.005, 0.005)", "id < 12", "colz"); + addTLines(o2::trk::SegmentationChip::PitchRowVD); + auto h = (TH2F*)gPad->GetPrimitive("h_dx_vs_dz_VD"); + LOG(info) << "dx, dz"; + Info("VD", "RMS(dx)=%.1f mu", h->GetRMS(2) * 1e4); + Info("VD", "RMS(dz)=%.1f mu", h->GetRMS(1) * 1e4); + canvdXdZ->cd(2); + nt->Draw("dx:dz>>h_dx_vs_dz_VD_z(500, -0.005, 0.005, 500, -0.005, 0.005)", "id < 12 && abs(z)<0.5", "colz"); + addTLines(o2::trk::SegmentationChip::PitchRowVD); + h = (TH2F*)gPad->GetPrimitive("h_dx_vs_dz_VD_z"); + Info("VD |z|<1", "RMS(dx)=%.1f mu", h->GetRMS(2) * 1e4); + Info("VD |z|<1", "RMS(dz)=%.1f mu", h->GetRMS(1) * 1e4); + canvdXdZ->cd(3); + nt->Draw("dx:dz>>h_dx_vs_dz_ML(600, -0.03, 0.03, 600, -0.03, 0.03)", "id >= 12 && id < 5132", "colz"); + addTLines(o2::trk::SegmentationChip::PitchRowMLOT); + h = (TH2F*)gPad->GetPrimitive("h_dx_vs_dz_ML"); + Info("ML", "RMS(dx)=%.1f mu", h->GetRMS(2) * 1e4); + Info("ML", "RMS(dz)=%.1f mu", h->GetRMS(1) * 1e4); + canvdXdZ->cd(4); + nt->Draw("dx:dz>>h_dx_vs_dz_ML_z(600, -0.03, 0.03, 600, -0.03, 0.03)", "id >= 12 && id < 5132 && abs(z)<2", "colz"); + addTLines(o2::trk::SegmentationChip::PitchRowMLOT); + h = (TH2F*)gPad->GetPrimitive("h_dx_vs_dz_ML_z"); + Info("ML |z|<2", "RMS(dx)=%.1f mu", h->GetRMS(2) * 1e4); + Info("ML |z|<2", "RMS(dz)=%.1f mu", h->GetRMS(1) * 1e4); + canvdXdZ->SaveAs("trkdigits_dx_vs_dz.pdf"); + canvdXdZ->cd(5); + nt->Draw("dx:dz>>h_dx_vs_dz_OT(600, -0.03, 0.03, 600, -0.03, 0.03)", "id >= 5132", "colz"); + addTLines(o2::trk::SegmentationChip::PitchRowMLOT); + h = (TH2F*)gPad->GetPrimitive("h_dx_vs_dz_OT"); + Info("OT", "RMS(dx)=%.1f mu", h->GetRMS(2) * 1e4); + Info("OT", "RMS(dz)=%.1f mu", h->GetRMS(1) * 1e4); + canvdXdZ->cd(6); + nt->Draw("dx:dz>>h_dx_vs_dz_OT_z(600, -0.03, 0.03, 600, -0.03, 0.03)", "id >= 5132 && abs(z)<2", "colz"); + h = (TH2F*)gPad->GetPrimitive("h_dx_vs_dz_OT_z"); + addTLines(o2::trk::SegmentationChip::PitchRowMLOT); + Info("OT |z|<2", "RMS(dx)=%.1f mu", h->GetRMS(2) * 1e4); + Info("OT |z|<2", "RMS(dz)=%.1f mu", h->GetRMS(1) * 1e4); + canvdXdZ->SaveAs("trkdigits_dx_vs_dz.pdf"); + canvdXdZ->SaveAs("trkdigits_dx_vs_dz.root"); + + // distribution of differences between hit start and hit end in local coordinates + auto canvdXdZHit = new TCanvas("canvdXdZHit", "", 1600, 2400); + canvdXdZHit->Divide(2, 3); + canvdXdZHit->cd(1); + nt2->Draw("dxH:dzH>>h_dxH_vs_dzH_VD(300, -0.03, 0.03, 300, -0.03, 0.03)", "id < 12", "colz"); + addTLines(o2::trk::SegmentationChip::PitchRowVD); + LOG(info) << "dxH, dzH"; + h = (TH2F*)gPad->GetPrimitive("h_dxH_vs_dzH_VD"); + Info("VD", "RMS(dxH)=%.1f mu", h->GetRMS(2) * 1e4); + Info("VD", "RMS(dzH)=%.1f mu", h->GetRMS(1) * 1e4); + canvdXdZHit->cd(2); + nt2->Draw("dxH:dzH>>h_dxH_vs_dzH_VD_z(300, -0.03, 0.03, 300, -0.03, 0.03)", "id < 12 && abs(z)<2", "colz"); + addTLines(o2::trk::SegmentationChip::PitchRowVD); + h = (TH2F*)gPad->GetPrimitive("h_dxH_vs_dzH_VD_z"); + Info("VD |z|<2", "RMS(dxH)=%.1f mu", h->GetRMS(2) * 1e4); + Info("VD |z|<2", "RMS(dzH)=%.1f mu", h->GetRMS(1) * 1e4); + canvdXdZHit->cd(3); + nt2->Draw("dxH:dzH>>h_dxH_vs_dzH_ML(300, -0.03, 0.03, 300, -0.03, 0.03)", "id >= 12 && id < 5132", "colz"); + addTLines(o2::trk::SegmentationChip::PitchRowMLOT); + h = (TH2F*)gPad->GetPrimitive("h_dxH_vs_dzH_ML"); + Info("ML", "RMS(dxH)=%.1f mu", h->GetRMS(2) * 1e4); + Info("ML", "RMS(dzH)=%.1f mu", h->GetRMS(1) * 1e4); + canvdXdZHit->cd(4); + nt2->Draw("dxH:dzH>>h_dxH_vs_dzH_ML_z(300, -0.03, 0.03, 300, -0.03, 0.03)", "id >= 12 && id < 5132 && abs(z)<2", "colz"); + addTLines(o2::trk::SegmentationChip::PitchRowMLOT); + h = (TH2F*)gPad->GetPrimitive("h_dxH_vs_dzH_ML_z"); + Info("ML |z|<2", "RMS(dxH)=%.1f mu", h->GetRMS(2) * 1e4); + Info("ML |z|<2", "RMS(dzH)=%.1f mu", h->GetRMS(1) * 1e4); + canvdXdZHit->SaveAs("trkdigits_dxH_vs_dzH.pdf"); + canvdXdZHit->cd(5); + nt2->Draw("dxH:dzH>>h_dxH_vs_dzH_OT(300, -0.03, 0.03, 300, -0.03, 0.03)", "id >= 5132", "colz"); + addTLines(o2::trk::SegmentationChip::PitchRowMLOT); + h = (TH2F*)gPad->GetPrimitive("h_dxH_vs_dzH_OT"); + Info("OT", "RMS(dxH)=%.1f mu", h->GetRMS(2) * 1e4); + Info("OT", "RMS(dzH)=%.1f mu", h->GetRMS(1) * 1e4); + canvdXdZHit->cd(6); + nt2->Draw("dxH:dzH>>h_dxH_vs_dzH_OT_z(300, -0.03, 0.03, 300, -0.03, 0.03)", "id >= 5132 && abs(z)<2", "colz"); + addTLines(o2::trk::SegmentationChip::PitchRowMLOT); + h = (TH2F*)gPad->GetPrimitive("h_dxH_vs_dzH_OT_z"); + Info("OT |z|<2", "RMS(dxH)=%.1f mu", h->GetRMS(2) * 1e4); + Info("OT |z|<2", "RMS(dzH)=%.1f mu", h->GetRMS(1) * 1e4); + canvdXdZHit->SaveAs("trkdigits_dxH_vs_dzH.pdf"); + + f->Write(); + f->Close(); +} From af16e0e16d19c0458906f0c75b7f541d2423005d Mon Sep 17 00:00:00 2001 From: Maximiliano Puccio Date: Fri, 24 Apr 2026 08:27:11 +0200 Subject: [PATCH 4/6] ALICE3: adapt cluster finding to separated streams per layer Fix CheckCluster macro --- .../TRK/base/include/TRKBase/AlmiraParam.h | 8 +- .../TRK/base/include/TRKBase/GeometryTGeo.h | 4 +- .../ALICE3/TRK/base/include/TRKBase/Specs.h | 4 +- .../ALICE3/TRK/macros/test/CheckBandwidth.C | 2 - .../ALICE3/TRK/macros/test/CheckClusters.C | 380 +++++++------- .../ALICE3/TRK/macros/test/CheckDigits.C | 467 ------------------ .../ALICE3/TRK/macros/test/CheckDigitsTRK.C | 250 +++++----- .../include/TRKReconstruction/Clusterer.h | 5 +- .../include/TRKReconstruction/ClustererACTS.h | 4 +- .../TRK/reconstruction/src/Clusterer.cxx | 11 +- .../TRK/reconstruction/src/ClustererACTS.cxx | 4 +- .../include/TRKSimulation/DigiParams.h | 2 +- .../include/TRKWorkflow/ClustererSpec.h | 2 + .../include/TRKWorkflow/DigitReaderSpec.h | 14 +- .../TRK/workflow/src/ClusterWriterSpec.cxx | 86 +++- .../ALICE3/TRK/workflow/src/ClustererSpec.cxx | 126 ++--- 16 files changed, 473 insertions(+), 896 deletions(-) delete mode 100644 Detectors/Upgrades/ALICE3/TRK/macros/test/CheckDigits.C diff --git a/Detectors/Upgrades/ALICE3/TRK/base/include/TRKBase/AlmiraParam.h b/Detectors/Upgrades/ALICE3/TRK/base/include/TRKBase/AlmiraParam.h index 987b459c04764..9929a14c4e39c 100644 --- a/Detectors/Upgrades/ALICE3/TRK/base/include/TRKBase/AlmiraParam.h +++ b/Detectors/Upgrades/ALICE3/TRK/base/include/TRKBase/AlmiraParam.h @@ -31,10 +31,11 @@ struct AlmiraParam : public o2::conf::ConfigurableParamHelper { int roFrameLengthInBCPerLayer[kNLayers] = {0}; ///< ROF length in BC per layer float strobeDelayPerLayer[kNLayers] = {0}; ///< strobe delay in ns per layer float strobeLengthContPerLayer[kNLayers] = {0}; ///< strobe length in ns per layer - int roFrameBiasInBCPerLayer[kNLayers] = {0}; ///< ROF start bias in BC per layer - int roFrameDelayInBCPerLayer[kNLayers] = {0}; ///< extra ROF delay in BC per layer + int roFrameBiasInBCPerLayer[kNLayers] = {0}; ///< ROF start bias in BC per layer + int roFrameDelayInBCPerLayer[kNLayers] = {0}; ///< extra ROF delay in BC per layer - int getROFLengthInBC(int layer) const { + int getROFLengthInBC(int layer) const + { if (roFrameLengthInBCPerLayer[layer] > 0) { return roFrameLengthInBCPerLayer[layer]; } else { @@ -46,7 +47,6 @@ struct AlmiraParam : public o2::conf::ConfigurableParamHelper { int getROFBiasInBC(int layer) const { return roFrameBiasInBCPerLayer[layer]; } int getROFDelayInBC(int layer) const { return roFrameDelayInBCPerLayer[layer]; } - O2ParamDef(AlmiraParam, "TRKAlmiraParam"); }; diff --git a/Detectors/Upgrades/ALICE3/TRK/base/include/TRKBase/GeometryTGeo.h b/Detectors/Upgrades/ALICE3/TRK/base/include/TRKBase/GeometryTGeo.h index bfc0b779eb1ca..53ad7662cbfcd 100644 --- a/Detectors/Upgrades/ALICE3/TRK/base/include/TRKBase/GeometryTGeo.h +++ b/Detectors/Upgrades/ALICE3/TRK/base/include/TRKBase/GeometryTGeo.h @@ -89,8 +89,8 @@ class GeometryTGeo : public o2::detectors::DetMatrixCache int getSubDetID(int index) const; int getPetalCase(int index) const; int getDisk(int index) const; - int getLayer(int index) const; ///< local layer index within the sub-detector (0-based per VD/MLOT) - int getLayerTRK(int index) const; ///< global layer index across the full TRK (VD layers 0..nVD-1, MLOT layers nVD..nTotal-1) + int getLayer(int index) const; ///< local layer index within the sub-detector (0-based per VD/MLOT) + int getLayerTRK(int index) const; ///< global layer index across the full TRK (VD layers 0..nVD-1, MLOT layers nVD..nTotal-1) int getStave(int index) const; int getHalfStave(int index) const; int getModule(int index) const; diff --git a/Detectors/Upgrades/ALICE3/TRK/base/include/TRKBase/Specs.h b/Detectors/Upgrades/ALICE3/TRK/base/include/TRKBase/Specs.h index aad9d045ac6f4..6be2e0bd0f827 100644 --- a/Detectors/Upgrades/ALICE3/TRK/base/include/TRKBase/Specs.h +++ b/Detectors/Upgrades/ALICE3/TRK/base/include/TRKBase/Specs.h @@ -100,7 +100,7 @@ constexpr int nCols{static_cast(length / chip::pitchZ)}; namespace ML { -constexpr int nLayers{3}; // number of layers in the ML +constexpr int nLayers{3}; // number of layers in the ML constexpr double width{constants::moduleMLOT::width * 1}; // width of the stave // constexpr double length{constants::moduleMLOT::length * 10}; // length of the stave constexpr double length{124 * cm}; // length of the stave, hardcoded to fit the implemented geometry @@ -118,7 +118,7 @@ constexpr double length{258 * cm}; // len constexpr int nRows{static_cast(width / moduleMLOT::chip::pitchX)}; // number of rows in the halfstave constexpr int nCols{static_cast(length / moduleMLOT::chip::pitchZ)}; // number of columns in the halfstave } // namespace halfstave -constexpr int nLayers{5}; // number of layers in the OT +constexpr int nLayers{5}; // number of layers in the OT constexpr double width{halfstave::width * 2}; // width of the stave constexpr double length{halfstave::length}; // length of the stave constexpr int nRows{static_cast(width / moduleMLOT::chip::pitchX)}; // number of rows in the stave diff --git a/Detectors/Upgrades/ALICE3/TRK/macros/test/CheckBandwidth.C b/Detectors/Upgrades/ALICE3/TRK/macros/test/CheckBandwidth.C index cda07447d58f5..c071a06516d30 100644 --- a/Detectors/Upgrades/ALICE3/TRK/macros/test/CheckBandwidth.C +++ b/Detectors/Upgrades/ALICE3/TRK/macros/test/CheckBandwidth.C @@ -193,7 +193,6 @@ void CheckBandwidth(std::string digifile = "trkdigits.root", std::string inputGe rofLengthBC[iLayer] = (*rofRecords[iLayer])[1].getBCData().bc - (*rofRecords[iLayer])[0].getBCData().bc; } - // --- Collision context --- TFile* ccFile = TFile::Open(collContextFile.data()); @@ -229,7 +228,6 @@ void CheckBandwidth(std::string digifile = "trkdigits.root", std::string inputGe // --- Accumulate per-chip digit counts across all ROFs --- - std::vector digitsPerChip(nChips, 0ull); std::vector maxDigitsPerROFPerChip(nChips, 0u); std::vector digitsInCurrentROFPerChip(nChips, 0u); diff --git a/Detectors/Upgrades/ALICE3/TRK/macros/test/CheckClusters.C b/Detectors/Upgrades/ALICE3/TRK/macros/test/CheckClusters.C index 327577102d86e..28dc61aed9c8b 100644 --- a/Detectors/Upgrades/ALICE3/TRK/macros/test/CheckClusters.C +++ b/Detectors/Upgrades/ALICE3/TRK/macros/test/CheckClusters.C @@ -26,6 +26,7 @@ #include "DataFormatsTRK/Cluster.h" #include "DataFormatsTRK/ROFRecord.h" +#include "TRKBase/AlmiraParam.h" #include "TRKBase/GeometryTGeo.h" #include "TRKBase/SegmentationChip.h" #include "TRKSimulation/Hit.h" @@ -142,47 +143,63 @@ void CheckClusters(const std::string& clusfile = "o2clus_trk.root", return; } - std::vector* clusArr = nullptr; - std::vector* rofRecVecP = nullptr; - std::vector* patternsPtr = nullptr; - clusTree->SetBranchAddress("TRKClusterComp", &clusArr); - clusTree->SetBranchAddress("TRKClustersROF", &rofRecVecP); - if (clusTree->GetBranch("TRKClusterPatt") != nullptr) { - clusTree->SetBranchAddress("TRKClusterPatt", &patternsPtr); + // Read per-layer cluster branches and accumulate + static constexpr int nLayers = o2::trk::AlmiraParam::kNLayers; + std::vector*> clusArrPerLayer(nLayers, nullptr); + std::vector*> rofRecVecPerLayer(nLayers, nullptr); + std::vector*> patternsPerLayer(nLayers, nullptr); + std::vector*> clusLabArrPerLayer(nLayers, nullptr); + + bool hasMC = true; + for (int iLayer = 0; iLayer < nLayers; iLayer++) { + std::string brClus = std::string("TRKClusterComp_") + std::to_string(iLayer); + std::string brROF = std::string("TRKClustersROF_") + std::to_string(iLayer); + std::string brPatt = std::string("TRKClusterPatt_") + std::to_string(iLayer); + std::string brMCTruth = std::string("TRKClusterMCTruth_") + std::to_string(iLayer); + + if (clusTree->GetBranch(brClus.c_str()) == nullptr) { + LOGP(warning, "Branch {} not found, skipping layer {}", brClus, iLayer); + continue; + } + clusTree->SetBranchAddress(brClus.c_str(), &clusArrPerLayer[iLayer]); + clusTree->SetBranchAddress(brROF.c_str(), &rofRecVecPerLayer[iLayer]); + if (clusTree->GetBranch(brPatt.c_str()) != nullptr) { + clusTree->SetBranchAddress(brPatt.c_str(), &patternsPerLayer[iLayer]); + } + if (clusTree->GetBranch(brMCTruth.c_str()) != nullptr) { + clusTree->SetBranchAddress(brMCTruth.c_str(), &clusLabArrPerLayer[iLayer]); + } else { + hasMC = false; + } } - o2::dataformats::MCTruthContainer* clusLabArr = nullptr; - std::vector mc2rofVec, *mc2rofVecP = &mc2rofVec; - bool hasMC = (clusTree->GetBranch("TRKClusterMCTruth") != nullptr); - if (hasMC) { - clusTree->SetBranchAddress("TRKClusterMCTruth", &clusLabArr); - clusTree->SetBranchAddress("TRKClustersMC2ROF", &mc2rofVecP); + // Read entry and accumulate all layers + clusTree->GetEntry(0); + // Print total clusters per layer + for (int iLayer = 0; iLayer < nLayers; iLayer++) { + LOGP(info, "Layer {}: {} clusters", iLayer, clusArrPerLayer[iLayer]->size()); } - clusTree->GetEntry(0); - const unsigned int nROFRec = rofRecVecP ? (unsigned int)rofRecVecP->size() : 0u; + // Accumulate max ROF count across all layers + unsigned int nROFRec = 0; + for (int iLayer = 0; iLayer < nLayers; iLayer++) { + nROFRec = std::max(nROFRec, (unsigned int)rofRecVecPerLayer[iLayer]->size()); + } LOGP(info, "Number of ROF records: {}", nROFRec); - auto pattIt = patternsPtr ? patternsPtr->cbegin() : std::vector::const_iterator{}; - // ── Build per-ROF MC event range ─────────────────────────────────────────── - std::vector mcEvMin(nROFRec, (int)hitTree->GetEntries()); - std::vector mcEvMax(nROFRec, -1); + // ── Load all MC hit events upfront (TRK has no MC2ROF mapping) ────────────── if (hasMC) { - for (int imc = (int)mc2rofVec.size(); imc--;) { - const auto& mc2rof = mc2rofVec[imc]; - if (mc2rof.rofRecordID < 0) { - continue; - } - for (unsigned int irfd = mc2rof.maxROF - mc2rof.minROF + 1; irfd--;) { - unsigned int irof = mc2rof.rofRecordID + irfd; - if (irof >= nROFRec) { - continue; - } - if (mcEvMin[irof] > imc) { - mcEvMin[irof] = imc; - } - if (mcEvMax[irof] < imc) { - mcEvMax[irof] = imc; + LOGP(info, "Pre-loading {} MC events", hitTree->GetEntries()); + for (int im = 0; im < (int)hitTree->GetEntries(); im++) { + if (hitVecPool[im] == nullptr) { + hitTree->SetBranchAddress("TRKHit", &hitVecPool[im]); + hitTree->GetEntry(im); + auto& mc2hit = mc2hitVec[im]; + const auto* hv = hitVecPool[im]; + for (int ih = (int)hv->size(); ih--;) { + const auto& hit = (*hv)[ih]; + uint64_t key = (uint64_t(hit.GetTrackID()) << 32) + hit.GetDetectorID(); + mc2hit.emplace(key, ih); } } } @@ -203,170 +220,165 @@ void CheckClusters(const std::string& clusfile = "o2clus_trk.root", // ── Main loop ────────────────────────────────────────────────────────────── for (unsigned int irof = 0; irof < nROFRec; irof++) { - const auto& rofRec = (*rofRecVecP)[irof]; - - // Cache MC hit events for this ROF - if (hasMC) { - for (int im = mcEvMin[irof]; im <= mcEvMax[irof]; im++) { - if (hitVecPool[im] == nullptr) { - hitTree->SetBranchAddress("TRKHit", &hitVecPool[im]); - hitTree->GetEntry(im); - auto& mc2hit = mc2hitVec[im]; - const auto* hv = hitVecPool[im]; - for (int ih = (int)hv->size(); ih--;) { - const auto& hit = (*hv)[ih]; - uint64_t key = (uint64_t(hit.GetTrackID()) << 32) + hit.GetDetectorID(); - mc2hit.emplace(key, ih); - } - } + // Process each layer + for (int iLayer = 0; iLayer < nLayers; iLayer++) { + if (rofRecVecPerLayer[iLayer]->empty() || irof >= rofRecVecPerLayer[iLayer]->size()) { + continue; } - } - - for (int icl = 0; icl < rofRec.getNEntries(); icl++) { - const int clEntry = rofRec.getFirstEntry() + icl; - const auto& cluster = (*clusArr)[clEntry]; - nTot++; - - // ── Parse pattern → center-of-gravity within bounding box ────────── - // The cluster stores the bounding-box top-left pixel (row, col). - // The pattern stream encodes [rowSpan, colSpan, bitmap...] for each cluster. - // We accumulate pixel row/col offsets to obtain a sub-pixel CoG correction. - float cogDr{0.f}, cogDc{0.f}; // mean offsets from bbox origin (pixels) - if (patternsPtr) { - const uint8_t rowSpan = *pattIt++; - const uint8_t colSpan = *pattIt++; - const int nBytes = (rowSpan * colSpan + 7) / 8; - int nPix{0}, pixIdx{0}; - for (int ib = 0; ib < nBytes; ib++) { - const uint8_t byte = *pattIt++; - for (int bit = 7; bit >= 0 && pixIdx < rowSpan * colSpan; bit--, pixIdx++) { - if (byte & (1 << bit)) { - cogDr += pixIdx / colSpan; - cogDc += pixIdx % colSpan; - nPix++; + const auto& rofRec = (*rofRecVecPerLayer[iLayer])[irof]; + const auto& clusArr = *clusArrPerLayer[iLayer]; + const auto& patternsPtr = (patternsPerLayer[iLayer] == nullptr) ? nullptr : patternsPerLayer[iLayer]; + const auto& clusLabArr = clusLabArrPerLayer[iLayer]; + + // Create per-layer pattern iterator + auto pattIt = patternsPtr ? patternsPtr->cbegin() : std::vector::const_iterator{}; + + for (int icl = 0; icl < rofRec.getNEntries(); icl++) { + const int clEntry = rofRec.getFirstEntry() + icl; + const auto& cluster = clusArr[clEntry]; + nTot++; + + // ── Parse pattern → center-of-gravity within bounding box ────────── + // The cluster stores the bounding-box top-left pixel (row, col). + // The pattern stream encodes [rowSpan, colSpan, bitmap...] for each cluster. + // We accumulate pixel row/col offsets to obtain a sub-pixel CoG correction. + float cogDr{0.f}, cogDc{0.f}; // mean offsets from bbox origin (pixels) + if (patternsPtr) { + const uint8_t rowSpan = *pattIt++; + const uint8_t colSpan = *pattIt++; + const int nBytes = (rowSpan * colSpan + 7) / 8; + int nPix{0}, pixIdx{0}; + for (int ib = 0; ib < nBytes; ib++) { + const uint8_t byte = *pattIt++; + for (int bit = 7; bit >= 0 && pixIdx < rowSpan * colSpan; bit--, pixIdx++) { + if (byte & (1 << bit)) { + cogDr += pixIdx / colSpan; + cogDc += pixIdx % colSpan; + nPix++; + } } } + if (nPix > 1) { + cogDr /= nPix; + cogDc /= nPix; + } } - if (nPix > 1) { - cogDr /= nPix; - cogDc /= nPix; + + // ── Cluster local → global (CoG position) ───────────────────────────── + // Get local coords of the bounding-box corner pixel, then apply the + // fractional CoG displacement using the pixel pitch. + // Formula from detectorToLocalUnchecked: + // VD : xRow = 0.5*(width[lay]-pitchRow) - row*pitchRow → row↑ xRow↓ + // zCol = col*pitchCol + 0.5*(pitchCol-length) → col↑ zCol↑ + // MLOT: same structure with MLOT pitches + float clLocX{0.f}, clLocZ{0.f}; + o2::trk::SegmentationChip::detectorToLocalUnchecked( + cluster.row, cluster.col, clLocX, clLocZ, + cluster.subDetID, cluster.layer, cluster.disk); + const float pitchRow = (cluster.subDetID == 0) + ? o2::trk::SegmentationChip::PitchRowVD + : o2::trk::SegmentationChip::PitchRowMLOT; + const float pitchCol = (cluster.subDetID == 0) + ? o2::trk::SegmentationChip::PitchColVD + : o2::trk::SegmentationChip::PitchColMLOT; + clLocX -= cogDr * pitchRow; // increasing row → decreasing xRow + clLocZ += cogDc * pitchCol; // increasing col → increasing zCol + const float yResponse = (cluster.subDetID == 0) ? yPlaneVD : yPlaneMLOT; + // For VD the L2G matrix is built in the *curved* local frame (quasi-Cartesian, + // origin at the beam axis). Convert flat (clLocX, 0) → curved (xC, yC) first. + // For MLOT (flat sensors) the local frame is already Cartesian: pass directly. + // clLocX is already in the flat frame from detectorToLocalUnchecked + CoG and + // does NOT need any further transformation for the residual comparison. + o2::math_utils::Point3D locC; + if (cluster.subDetID == 0) { + auto cv = o2::trk::SegmentationChip::flatToCurved(cluster.layer, clLocX, 0.f); + locC = {cv.X(), cv.Y(), clLocZ}; + } else { + locC = {clLocX, yResponse, clLocZ}; + } + auto gloC = gman->getMatrixL2G(cluster.chipID)(locC); + + if (!hasMC || clusLabArr == nullptr) { + // No MC info: just fill geometry columns, leave residuals as 0 + std::array data = { + -1.f, -1.f, + 0.f, 0.f, 0.f, 0.f, 0.f, + (float)gloC.X(), (float)gloC.Y(), (float)gloC.Z(), + clLocX, clLocZ, + (float)rofRec.getROFrame(), (float)cluster.size, (float)cluster.chipID, + (float)cluster.layer, (float)cluster.disk, (float)cluster.subDetID, + (float)cluster.row, (float)cluster.col, -1.f}; + nt.Fill(data.data()); + continue; } - } - // ── Cluster local → global (CoG position) ───────────────────────────── - // Get local coords of the bounding-box corner pixel, then apply the - // fractional CoG displacement using the pixel pitch. - // Formula from detectorToLocalUnchecked: - // VD : xRow = 0.5*(width[lay]-pitchRow) - row*pitchRow → row↑ xRow↓ - // zCol = col*pitchCol + 0.5*(pitchCol-length) → col↑ zCol↑ - // MLOT: same structure with MLOT pitches - float clLocX{0.f}, clLocZ{0.f}; - o2::trk::SegmentationChip::detectorToLocalUnchecked( - cluster.row, cluster.col, clLocX, clLocZ, - cluster.subDetID, cluster.layer, cluster.disk); - const float pitchRow = (cluster.subDetID == 0) - ? o2::trk::SegmentationChip::PitchRowVD - : o2::trk::SegmentationChip::PitchRowMLOT; - const float pitchCol = (cluster.subDetID == 0) - ? o2::trk::SegmentationChip::PitchColVD - : o2::trk::SegmentationChip::PitchColMLOT; - clLocX -= cogDr * pitchRow; // increasing row → decreasing xRow - clLocZ += cogDc * pitchCol; // increasing col → increasing zCol - const float yResponse = (cluster.subDetID == 0) ? yPlaneVD : yPlaneMLOT; - // For VD the L2G matrix is built in the *curved* local frame (quasi-Cartesian, - // origin at the beam axis). Convert flat (clLocX, 0) → curved (xC, yC) first. - // For MLOT (flat sensors) the local frame is already Cartesian: pass directly. - // clLocX is already in the flat frame from detectorToLocalUnchecked + CoG and - // does NOT need any further transformation for the residual comparison. - o2::math_utils::Point3D locC; - if (cluster.subDetID == 0) { - auto cv = o2::trk::SegmentationChip::flatToCurved(cluster.layer, clLocX, 0.f); - locC = {cv.X(), cv.Y(), clLocZ}; - } else { - locC = {clLocX, yResponse, clLocZ}; - } - auto gloC = gman->getMatrixL2G(cluster.chipID)(locC); + // ── MC label ─────────────────────────────────────────────────────── + const auto& labels = clusLabArr->getLabels(clEntry); + if (labels.empty() || !labels[0].isValid()) { + nInvalidLabel++; + continue; + } + const auto& lab = labels[0]; + const int trID = lab.getTrackID(); + const int evID = lab.getEventID(); + + // ── Find matching MC hit ──────────────────────────────────────────── + const auto& mc2hit = mc2hitVec[evID]; + uint64_t key = (uint64_t(trID) << 32) + cluster.chipID; + auto hitEntry = mc2hit.find(key); + if (hitEntry == mc2hit.end()) { + nNoMCHit++; + continue; + } + const auto& hit = (*hitVecPool[evID])[hitEntry->second]; + const float pt = TMath::Hypot(hit.GetPx(), hit.GetPy()); + + // ── Hit global midpoint ──────────────────────────────────────────── + const auto& gloHend = hit.GetPos(); + const auto& gloHsta = hit.GetPosStart(); + o2::math_utils::Point3D gloHmid( + 0.5f * (gloHend.X() + gloHsta.X()), + 0.5f * (gloHend.Y() + gloHsta.Y()), + 0.5f * (gloHend.Z() + gloHsta.Z())); + + // ── Hit global → local ───────────────────────────── + o2::math_utils::Point3D locHsta = gman->getMatrixL2G(cluster.chipID) ^ (gloHsta); // inverse L2G + o2::math_utils::Point3D locHend = gman->getMatrixL2G(cluster.chipID) ^ (gloHend); // inverse L2G + + // ── Propagate hit segment to the sensor response surface ─────────────── + // Rather than the geometric midpoint, find where the track segment crosses + // the response plane (y = responseYShift in the flat local frame). + // For VD (curved): convert both endpoints to flat frame first. + // For ML/OT (flat): use local coordinates directly. + float hitLocX{0.f}, hitLocZ{0.f}; + if (cluster.subDetID == 0) { // VD – curved sensor + auto flatSta = o2::trk::SegmentationChip::curvedToFlat(cluster.layer, locHsta.X(), locHsta.Y()); + auto flatEnd = o2::trk::SegmentationChip::curvedToFlat(cluster.layer, locHend.X(), locHend.Y()); + float x0 = flatSta.X(), y0 = flatSta.Y(), z0 = locHsta.Z(); + float dltx = flatEnd.X() - x0, dlty = flatEnd.Y() - y0, dltz = locHend.Z() - z0; + float r = (std::abs(dlty) > 1e-9f) ? (yPlaneVD - y0) / dlty : 0.5f; + hitLocX = x0 + r * dltx; + hitLocZ = z0 + r * dltz; + } else { // ML/OT – flat sensor + float x0 = locHsta.X(), y0 = locHsta.Y(), z0 = locHsta.Z(); + float dltx = locHend.X() - x0, dlty = locHend.Y() - y0, dltz = locHend.Z() - z0; + float r = (std::abs(dlty) > 1e-9f) ? (yPlaneMLOT - y0) / dlty : 0.5f; + hitLocX = x0 + r * dltx; + hitLocZ = z0 + r * dltz; + } - if (!hasMC || clusLabArr == nullptr) { - // No MC info: just fill geometry columns, leave residuals as 0 + nValid++; std::array data = { - -1.f, -1.f, - 0.f, 0.f, 0.f, 0.f, 0.f, + (float)evID, (float)trID, + hitLocX, hitLocZ, + (float)gloHmid.X(), (float)gloHmid.Y(), (float)gloHmid.Z(), (float)gloC.X(), (float)gloC.Y(), (float)gloC.Z(), clLocX, clLocZ, (float)rofRec.getROFrame(), (float)cluster.size, (float)cluster.chipID, (float)cluster.layer, (float)cluster.disk, (float)cluster.subDetID, - (float)cluster.row, (float)cluster.col, -1.f}; + (float)cluster.row, (float)cluster.col, pt}; nt.Fill(data.data()); - continue; - } - - // ── MC label ─────────────────────────────────────────────────────── - const auto& labels = clusLabArr->getLabels(clEntry); - if (labels.empty() || !labels[0].isValid()) { - nInvalidLabel++; - continue; } - const auto& lab = labels[0]; - const int trID = lab.getTrackID(); - const int evID = lab.getEventID(); - - // ── Find matching MC hit ──────────────────────────────────────────── - const auto& mc2hit = mc2hitVec[evID]; - uint64_t key = (uint64_t(trID) << 32) + cluster.chipID; - auto hitEntry = mc2hit.find(key); - if (hitEntry == mc2hit.end()) { - nNoMCHit++; - continue; - } - const auto& hit = (*hitVecPool[evID])[hitEntry->second]; - const float pt = TMath::Hypot(hit.GetPx(), hit.GetPy()); - - // ── Hit global midpoint ──────────────────────────────────────────── - const auto& gloHend = hit.GetPos(); - const auto& gloHsta = hit.GetPosStart(); - o2::math_utils::Point3D gloHmid( - 0.5f * (gloHend.X() + gloHsta.X()), - 0.5f * (gloHend.Y() + gloHsta.Y()), - 0.5f * (gloHend.Z() + gloHsta.Z())); - - // ── Hit global → local ───────────────────────────── - o2::math_utils::Point3D locHsta = gman->getMatrixL2G(cluster.chipID) ^ (gloHsta); // inverse L2G - o2::math_utils::Point3D locHend = gman->getMatrixL2G(cluster.chipID) ^ (gloHend); // inverse L2G - - // ── Propagate hit segment to the sensor response surface ─────────────── - // Rather than the geometric midpoint, find where the track segment crosses - // the response plane (y = responseYShift in the flat local frame). - // For VD (curved): convert both endpoints to flat frame first. - // For ML/OT (flat): use local coordinates directly. - float hitLocX{0.f}, hitLocZ{0.f}; - if (cluster.subDetID == 0) { // VD – curved sensor - auto flatSta = o2::trk::SegmentationChip::curvedToFlat(cluster.layer, locHsta.X(), locHsta.Y()); - auto flatEnd = o2::trk::SegmentationChip::curvedToFlat(cluster.layer, locHend.X(), locHend.Y()); - float x0 = flatSta.X(), y0 = flatSta.Y(), z0 = locHsta.Z(); - float dltx = flatEnd.X() - x0, dlty = flatEnd.Y() - y0, dltz = locHend.Z() - z0; - float r = (std::abs(dlty) > 1e-9f) ? (yPlaneVD - y0) / dlty : 0.5f; - hitLocX = x0 + r * dltx; - hitLocZ = z0 + r * dltz; - } else { // ML/OT – flat sensor - float x0 = locHsta.X(), y0 = locHsta.Y(), z0 = locHsta.Z(); - float dltx = locHend.X() - x0, dlty = locHend.Y() - y0, dltz = locHend.Z() - z0; - float r = (std::abs(dlty) > 1e-9f) ? (yPlaneMLOT - y0) / dlty : 0.5f; - hitLocX = x0 + r * dltx; - hitLocZ = z0 + r * dltz; - } - - nValid++; - std::array data = { - (float)evID, (float)trID, - hitLocX, hitLocZ, - (float)gloHmid.X(), (float)gloHmid.Y(), (float)gloHmid.Z(), - (float)gloC.X(), (float)gloC.Y(), (float)gloC.Z(), - clLocX, clLocZ, - (float)rofRec.getROFrame(), (float)cluster.size, (float)cluster.chipID, - (float)cluster.layer, (float)cluster.disk, (float)cluster.subDetID, - (float)cluster.row, (float)cluster.col, pt}; - nt.Fill(data.data()); } } diff --git a/Detectors/Upgrades/ALICE3/TRK/macros/test/CheckDigits.C b/Detectors/Upgrades/ALICE3/TRK/macros/test/CheckDigits.C deleted file mode 100644 index ec1adf500f562..0000000000000 --- a/Detectors/Upgrades/ALICE3/TRK/macros/test/CheckDigits.C +++ /dev/null @@ -1,467 +0,0 @@ -// 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 CheckDigits.C -/// \brief Simple macro to check TRK digits - -#if !defined(__CLING__) || defined(__ROOTCLING__) -#include -#include -#include -#include -#include -#include -#include -#include - -#include "TRKBase/SegmentationChip.h" -#include "TRKBase/GeometryTGeo.h" -#include "DataFormatsITSMFT/Digit.h" -#include "TRKSimulation/Hit.h" -#include "MathUtils/Utils.h" -#include "SimulationDataFormat/ConstMCTruthContainer.h" -#include "SimulationDataFormat/IOMCTruthContainerView.h" -#include "SimulationDataFormat/MCCompLabel.h" -#include "DetectorsBase/GeometryManager.h" -#include "ITSMFTSimulation/AlpideSimResponse.h" -#include "CCDB/BasicCCDBManager.h" - -#include "DataFormatsITSMFT/ROFRecord.h" - -#endif - -#define ENABLE_UPGRADES - -void addTLines(float pitch) -{ - // Add grid lines at multiples of pitch on the current pad - if (!gPad) - return; - - gPad->Update(); - - Double_t xmin = gPad->GetUxmin(); - Double_t xmax = gPad->GetUxmax(); - Double_t ymin = gPad->GetUymin(); - Double_t ymax = gPad->GetUymax(); - - // Calculate the first vertical line position (multiple of pitch) - int nLinesX = 0; - for (float x = xmin; x <= xmax && nLinesX < 1000; x += pitch, nLinesX++) { - TLine* line = new TLine(x, ymin, x, ymax); - line->SetLineStyle(2); - line->SetLineColor(kGray); - line->Draw("same"); - } - - // Calculate the first horizontal line position (multiple of pitch) - int nLinesY = 0; - for (float y = ymin; y <= ymax && nLinesY < 1000; y += pitch, nLinesY++) { - TLine* line = new TLine(xmin, y, xmax, y); - line->SetLineStyle(2); - line->SetLineColor(kGray); - line->Draw("same"); - } - - gPad->Modified(); - gPad->Update(); -} - -void CheckDigits(std::string digifile = "trkdigits.root", std::string hitfile = "o2sim_HitsTRK.root", std::string inputGeom = "o2sim_geometry.root", std::string paramfile = "o2sim_par.root") -{ - gStyle->SetPalette(55); - - using namespace o2::base; - using namespace o2::trk; - - using o2::itsmft::Digit; - using o2::trk::Hit; - - using o2::trk::SegmentationChip; - - TFile* f = TFile::Open("CheckDigits.root", "recreate"); - - TNtuple* nt = new TNtuple("ntd", "digit ntuple", "id:x:y:z:rowD:colD:rowH:colH:xlH:zlH:xlcH:zlcH:dx:dz"); - TNtuple* nt2 = new TNtuple("ntd2", "digit ntuple", "id:z:dxH:dzH"); /// maximum number of elements in a tuple = 15: doing a new tuple to store more variables - - // Geometry - o2::base::GeometryManager::loadGeometry(inputGeom); - auto* gman = o2::trk::GeometryTGeo::Instance(); - gman->fillMatrixCache(o2::math_utils::bit2Mask(o2::math_utils::TransformType::L2G)); - - SegmentationChip seg; - // seg.Print(); - - // MLOT response plane: y = halfThickness - depthMax. - float depthMax = (float)o2::trk::constants::apts::thickness; // fallback (no CCDB) - auto& ccdbMgr = o2::ccdb::BasicCCDBManager::instance(); - ccdbMgr.setURL("http://alice-ccdb.cern.ch"); - if (auto* alpResp = ccdbMgr.get("IT3/Calib/APTSResponse")) { - depthMax = alpResp->getDepthMax(); - } - const float yPlaneMLOT = o2::trk::SegmentationChip::SiliconThicknessMLOT / 2.f - depthMax; - const float yPlaneVD = -o2::trk::SegmentationChip::SiliconThicknessVD; // VD reference plane in local flat y - // 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); - - // Digits - TFile* digFile = TFile::Open(digifile.data()); - TTree* digTree = (TTree*)digFile->Get("o2sim"); - - std::vector* digArr = nullptr; - digTree->SetBranchAddress("TRKDigit", &digArr); - - o2::dataformats::IOMCTruthContainerView* plabels = nullptr; - digTree->SetBranchAddress("TRKDigitMCTruth", &plabels); - - // Get Read Out Frame arrays - std::vector* ROFRecordArrray = nullptr; - digTree->SetBranchAddress("TRKDigitROF", &ROFRecordArrray); - std::vector& ROFRecordArrrayRef = *ROFRecordArrray; - - std::vector* MC2ROFRecordArrray = nullptr; - digTree->SetBranchAddress("TRKDigitMC2ROF", &MC2ROFRecordArrray); - std::vector& MC2ROFRecordArrrayRef = *MC2ROFRecordArrray; - - digTree->GetEntry(0); - - int nROFRec = (int)ROFRecordArrrayRef.size(); - std::vector mcEvMin(nROFRec, hitTree->GetEntries()); - std::vector mcEvMax(nROFRec, -1); - o2::dataformats::ConstMCTruthContainer labels; - plabels->copyandflatten(labels); - delete plabels; - - // >> build min and max MC events used by each ROF - for (int imc = MC2ROFRecordArrrayRef.size(); imc--;) { - const auto& mc2rof = MC2ROFRecordArrrayRef[imc]; - // printf("MCRecord: "); - // mc2rof.print(); - - if (mc2rof.rofRecordID < 0) { - continue; // this MC event did not contribute to any ROF - } - - for (int irfd = mc2rof.maxROF - mc2rof.minROF + 1; irfd--;) { - - int irof = mc2rof.rofRecordID + irfd; - - if (irof >= nROFRec) { - LOG(error) << "ROF=" << irof << " from MC2ROF record is >= N ROFs=" << nROFRec; - } - if (mcEvMin[irof] > imc) { - mcEvMin[irof] = imc; - } - if (mcEvMax[irof] < imc) { - mcEvMax[irof] = imc; - } - } - } // << build min and max MC events used by each ROF - - unsigned int rofIndex = 0; - unsigned int rofNEntries = 0; - - // LOOP on : ROFRecord array - for (unsigned int iROF = 0; iROF < ROFRecordArrrayRef.size(); iROF++) { - - rofIndex = ROFRecordArrrayRef[iROF].getFirstEntry(); - rofNEntries = ROFRecordArrrayRef[iROF].getNEntries(); - - // >> read and map MC events contributing to this ROF - for (int im = mcEvMin[iROF]; im <= mcEvMax[iROF]; im++) { - - if (!hitArray[im]) { - - hitTree->SetBranchAddress("TRKHit", &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); - } - } - } - - // LOOP on : digits array - for (unsigned int iDigit = rofIndex; iDigit < rofIndex + rofNEntries; iDigit++) { - // if (iDigit % 10000 != 0) /// looking only at a small sample - // continue; - - if (iDigit % 1000 == 0) - std::cout << "Reading digit " << iDigit << " / " << digArr->size() << std::endl; - - Int_t ix = (*digArr)[iDigit].getRow(), iz = (*digArr)[iDigit].getColumn(); - Int_t iDetID = (*digArr)[iDigit].getChipIndex(); - Int_t layer = gman->getLayer(iDetID); - Int_t disk = gman->getDisk(iDetID); - Int_t subDetID = gman->getSubDetID(iDetID); - Int_t petalCase = gman->getPetalCase(iDetID); - Int_t stave = gman->getStave(iDetID); - Int_t halfstave = gman->getHalfStave(iDetID); - - Float_t x = 0.f, y = 0.f, z = 0.f; - Float_t x_flat = 0.f, z_flat = 0.f; - - if (disk != -1) { - continue; // skip disks for the moment - } - - if (subDetID != 0) { - seg.detectorToLocal(ix, iz, x, z, subDetID, layer, disk); - } else if (subDetID == 0) { - seg.detectorToLocal(ix, iz, x_flat, z_flat, subDetID, layer, disk); - o2::math_utils::Vector2D xyCurved = seg.flatToCurved(layer, x_flat, 0.); - x = xyCurved.X(); - y = xyCurved.Y(); - z = z_flat; - } - - o2::math_utils::Point3D locD(x, y, z); // local Digit curved - o2::math_utils::Point3D locDF(-1, -1, -1); // local Digit flat - - Int_t chipID = (*digArr)[iDigit].getChipIndex(); - auto lab = (labels.getLabels(iDigit))[0]; - - int trID = lab.getTrackID(); - - if (!lab.isValid()) { // not a noise - continue; - } - - const auto gloD = gman->getMatrixL2G(chipID)(locD); // convert to global - - std::unordered_map* mc2hit = &mc2hitVec[lab.getEventID()]; - - // get MC info - 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; - continue; - } - - ////// HITS - Hit& hit = (*hitArray[lab.getEventID()])[hitEntry->second]; - - auto xyzLocE = gman->getMatrixL2G(chipID) ^ (hit.GetPos()); // inverse conversion from global to local - auto xyzLocS = gman->getMatrixL2G(chipID) ^ (hit.GetPosStart()); - - // Hit local reference: Both VD and MLOT use response-plane interpolation (in flat local frame). - // For VD, transform curved → flat first, then interpolate. - o2::math_utils::Vector3D locH; /// Hit reference (at response plane) - o2::math_utils::Vector3D locHS; /// Hit, start pos - locHS.SetCoordinates(xyzLocS.X(), xyzLocS.Y(), xyzLocS.Z()); - o2::math_utils::Vector3D locHE; /// Hit, end pos - locHE.SetCoordinates(xyzLocE.X(), xyzLocE.Y(), xyzLocE.Z()); - o2::math_utils::Vector3D locHF; - - if (subDetID == 0) { - // VD: Interpolate to VD reference plane in flat frame; apply same r to X and Z - auto flatSta = seg.curvedToFlat(layer, locHS.X(), locHS.Y()); - auto flatEnd = seg.curvedToFlat(layer, locHE.X(), locHE.Y()); - float x0 = flatSta.X(), y0 = flatSta.Y(), z0 = locHS.Z(); - float dltx = flatEnd.X() - x0, dlty = flatEnd.Y() - y0, dltz = locHE.Z() - z0; - float r = (std::abs(dlty) > 1e-9f) ? (yPlaneVD - y0) / dlty : 0.5f; - locH.SetCoordinates(x0 + r * dltx, yPlaneVD, z0 + r * dltz); - } else { - // MLOT: Interpolate to response plane - float x0 = locHS.X(), y0 = locHS.Y(), z0 = locHS.Z(); - float dltx = locHE.X() - x0, dlty = locHE.Y() - y0, dltz = locHE.Z() - z0; - float r = (std::abs(dlty) > 1e-9f) ? (yPlaneMLOT - y0) / dlty : 0.5f; - locH.SetCoordinates(x0 + r * dltx, yPlaneMLOT, z0 + r * dltz); - } - - int row = 0, col = 0; - float xlc = 0., zlc = 0.; - - if (subDetID == 0) { - Float_t x_flat = 0.f, y_flat = 0.f; - // locH is already in flat frame from interpolation above; convert digit to flat for comparison - o2::math_utils::Vector2D xyFlatD = seg.curvedToFlat(layer, locD.X(), locD.Y()); - locDF.SetCoordinates(xyFlatD.X(), xyFlatD.Y(), locD.Z()); - locHF.SetCoordinates(locH.X(), locH.Y(), locH.Z()); // locH already in flat frame - seg.localToDetector(locHF.X(), locHF.Z(), row, col, subDetID, layer, disk); - } - - else { - seg.localToDetector(locH.X(), locH.Z(), row, col, subDetID, layer, disk); - } - - seg.detectorToLocal(row, col, xlc, zlc, subDetID, layer, disk); - - if (subDetID == 0) { - nt->Fill(chipID, /// detector ID - gloD.X(), gloD.Y(), gloD.Z(), /// global position retrieved from the digit: digit (row, col) ->local position -> global potision - ix, iz, /// row and column of the digit - row, col, /// row and col retrieved from the hit: hit global position -> hit local position -> detector position (row, col) - locH.X(), locH.Z(), /// x and z of the hit in the local reference frame: hit global position -> hit local position - xlc, zlc, /// x and z of the hit in the local frame: hit global position -> hit local position -> detector position (row, col) -> local position - locHF.X() - locDF.X(), locHF.Z() - locDF.Z()); /// difference in x and z between the hit and the digit in the local frame - - nt2->Fill(chipID, gloD.Z(), locHS.X() - locHE.X(), locHS.Z() - locHE.Z()); /// differences between local hit start and hit end positions - } else { - - nt->Fill(chipID, /// detector ID - gloD.X(), gloD.Y(), gloD.Z(), /// global position retrieved from the digit: digit (row, col) ->local position -> global potision - ix, iz, /// row and column of the digit - row, col, /// row and col retrieved from the hit: hit global position -> hit local position -> detector position (row, col) - locH.X(), locH.Z(), /// x and z of the hit in the local reference frame: hit global position -> hit local position - xlc, zlc, /// x and z of the hit in the local frame: hit global position -> hit local position -> detector position (row, col) -> local position - locH.X() - locD.X(), locH.Z() - locD.Z()); /// difference in x and z between the hit and the digit in the local frame - // locHS.X() - locHE.X(), locHS.Z() - locHE.Z()); /// difference in x and z between the hit and the digit in the local frame - nt2->Fill(chipID, gloD.Z(), locHS.X() - locHE.X(), locHS.Z() - locHE.Z()); /// differences between local hit start and hit end positions - } - - } // end loop on digits array - - } // end loop on ROFRecords array - - // digit maps in the xy and yz planes - auto canvXY = new TCanvas("canvXY", "", 1600, 2400); - canvXY->Divide(2, 3); - canvXY->cd(1); - nt->Draw("y:x >>h_y_vs_x_VD(1000, -3, 3, 1000, -3, 3)", "id < 12 ", "colz"); - canvXY->cd(2); - nt->Draw("y:z>>h_y_vs_z_VD(1000, -26, 26, 1000, -3, 3)", "id < 12 ", "colz"); - canvXY->cd(3); - nt->Draw("y:x>>h_y_vs_x_ML(1000, -25, 25, 1000, -25, 25)", "id >= 12 && id < 5132 ", "colz"); - canvXY->cd(4); - nt->Draw("y:z>>h_y_vs_z_ML(1000, -70, 70, 1000, -25, 25)", "id >= 12 && id < 5132 ", "colz"); - canvXY->cd(5); - nt->Draw("y:x>>h_y_vs_x_OT(1000, -85, 85, 1000, -85, 85)", "id >= 5132 ", "colz"); - canvXY->cd(6); - nt->Draw("y:z>>h_y_vs_z_OT(1000, -85, 85, 1000, -130, 130)", "id >= 5132 ", "colz"); - canvXY->SaveAs("trkdigits_y_vs_x_vs_z.pdf"); - - // z distributions - auto canvZ = new TCanvas("canvZ", "", 800, 2400); - canvZ->Divide(1, 3); - canvZ->cd(1); - nt->Draw("z>>h_z_VD(500, -26, 26)", "id < 12 "); - canvZ->cd(2); - nt->Draw("z>>h_z_ML(500, -70, 70)", "id >= 12 && id < 5132 "); - canvZ->cd(3); - nt->Draw("z>>h_z_OT(500, -85, 85)", "id >= 5132 "); - canvZ->SaveAs("trkdigits_z.pdf"); - - // dz distributions (difference between local position of digits and hits in x and z) - auto canvdZ = new TCanvas("canvdZ", "", 800, 2400); - canvdZ->Divide(1, 3); - canvdZ->cd(1); - nt->Draw("dz>>h_dz_VD(500, -0.05, 0.05)", "id < 12 "); - canvdZ->cd(2); - nt->Draw("dz>>h_dz_ML(500, -0.05, 0.05)", "id >= 12 && id < 5132 "); - canvdZ->cd(3); - nt->Draw("dz>>h_dz_OT(500, -0.05, 0.05)", "id >= 5132 "); - canvdZ->SaveAs("trkdigits_dz.pdf"); - canvdZ->SaveAs("trkdigits_dz.root"); - - // distributions of differences between local positions of digits and hits in x and z - auto canvdXdZ = new TCanvas("canvdXdZ", "", 1600, 2400); - canvdXdZ->Divide(2, 3); - canvdXdZ->cd(1); - nt->Draw("dx:dz>>h_dx_vs_dz_VD(500, -0.005, 0.005, 500, -0.005, 0.005)", "id < 12", "colz"); - addTLines(o2::trk::SegmentationChip::PitchRowVD); - auto h = (TH2F*)gPad->GetPrimitive("h_dx_vs_dz_VD"); - LOG(info) << "dx, dz"; - Info("VD", "RMS(dx)=%.1f mu", h->GetRMS(2) * 1e4); - Info("VD", "RMS(dz)=%.1f mu", h->GetRMS(1) * 1e4); - canvdXdZ->cd(2); - nt->Draw("dx:dz>>h_dx_vs_dz_VD_z(500, -0.005, 0.005, 500, -0.005, 0.005)", "id < 12 && abs(z)<0.5", "colz"); - addTLines(o2::trk::SegmentationChip::PitchRowVD); - h = (TH2F*)gPad->GetPrimitive("h_dx_vs_dz_VD_z"); - Info("VD |z|<1", "RMS(dx)=%.1f mu", h->GetRMS(2) * 1e4); - Info("VD |z|<1", "RMS(dz)=%.1f mu", h->GetRMS(1) * 1e4); - canvdXdZ->cd(3); - nt->Draw("dx:dz>>h_dx_vs_dz_ML(600, -0.03, 0.03, 600, -0.03, 0.03)", "id >= 12 && id < 5132", "colz"); - addTLines(o2::trk::SegmentationChip::PitchRowMLOT); - h = (TH2F*)gPad->GetPrimitive("h_dx_vs_dz_ML"); - Info("ML", "RMS(dx)=%.1f mu", h->GetRMS(2) * 1e4); - Info("ML", "RMS(dz)=%.1f mu", h->GetRMS(1) * 1e4); - canvdXdZ->cd(4); - nt->Draw("dx:dz>>h_dx_vs_dz_ML_z(600, -0.03, 0.03, 600, -0.03, 0.03)", "id >= 12 && id < 5132 && abs(z)<2", "colz"); - addTLines(o2::trk::SegmentationChip::PitchRowMLOT); - h = (TH2F*)gPad->GetPrimitive("h_dx_vs_dz_ML_z"); - Info("ML |z|<2", "RMS(dx)=%.1f mu", h->GetRMS(2) * 1e4); - Info("ML |z|<2", "RMS(dz)=%.1f mu", h->GetRMS(1) * 1e4); - canvdXdZ->SaveAs("trkdigits_dx_vs_dz.pdf"); - canvdXdZ->cd(5); - nt->Draw("dx:dz>>h_dx_vs_dz_OT(600, -0.03, 0.03, 600, -0.03, 0.03)", "id >= 5132", "colz"); - addTLines(o2::trk::SegmentationChip::PitchRowMLOT); - h = (TH2F*)gPad->GetPrimitive("h_dx_vs_dz_OT"); - Info("OT", "RMS(dx)=%.1f mu", h->GetRMS(2) * 1e4); - Info("OT", "RMS(dz)=%.1f mu", h->GetRMS(1) * 1e4); - canvdXdZ->cd(6); - nt->Draw("dx:dz>>h_dx_vs_dz_OT_z(600, -0.03, 0.03, 600, -0.03, 0.03)", "id >= 5132 && abs(z)<2", "colz"); - h = (TH2F*)gPad->GetPrimitive("h_dx_vs_dz_OT_z"); - addTLines(o2::trk::SegmentationChip::PitchRowMLOT); - Info("OT |z|<2", "RMS(dx)=%.1f mu", h->GetRMS(2) * 1e4); - Info("OT |z|<2", "RMS(dz)=%.1f mu", h->GetRMS(1) * 1e4); - canvdXdZ->SaveAs("trkdigits_dx_vs_dz.pdf"); - canvdXdZ->SaveAs("trkdigits_dx_vs_dz.root"); - - // distribution of differences between hit start and hit end in local coordinates - auto canvdXdZHit = new TCanvas("canvdXdZHit", "", 1600, 2400); - canvdXdZHit->Divide(2, 3); - canvdXdZHit->cd(1); - nt2->Draw("dxH:dzH>>h_dxH_vs_dzH_VD(300, -0.03, 0.03, 300, -0.03, 0.03)", "id < 12", "colz"); - addTLines(o2::trk::SegmentationChip::PitchRowVD); - LOG(info) << "dxH, dzH"; - h = (TH2F*)gPad->GetPrimitive("h_dxH_vs_dzH_VD"); - Info("VD", "RMS(dxH)=%.1f mu", h->GetRMS(2) * 1e4); - Info("VD", "RMS(dzH)=%.1f mu", h->GetRMS(1) * 1e4); - canvdXdZHit->cd(2); - nt2->Draw("dxH:dzH>>h_dxH_vs_dzH_VD_z(300, -0.03, 0.03, 300, -0.03, 0.03)", "id < 12 && abs(z)<2", "colz"); - addTLines(o2::trk::SegmentationChip::PitchRowVD); - h = (TH2F*)gPad->GetPrimitive("h_dxH_vs_dzH_VD_z"); - Info("VD |z|<2", "RMS(dxH)=%.1f mu", h->GetRMS(2) * 1e4); - Info("VD |z|<2", "RMS(dzH)=%.1f mu", h->GetRMS(1) * 1e4); - canvdXdZHit->cd(3); - nt2->Draw("dxH:dzH>>h_dxH_vs_dzH_ML(300, -0.03, 0.03, 300, -0.03, 0.03)", "id >= 12 && id < 5132", "colz"); - addTLines(o2::trk::SegmentationChip::PitchRowMLOT); - h = (TH2F*)gPad->GetPrimitive("h_dxH_vs_dzH_ML"); - Info("ML", "RMS(dxH)=%.1f mu", h->GetRMS(2) * 1e4); - Info("ML", "RMS(dzH)=%.1f mu", h->GetRMS(1) * 1e4); - canvdXdZHit->cd(4); - nt2->Draw("dxH:dzH>>h_dxH_vs_dzH_ML_z(300, -0.03, 0.03, 300, -0.03, 0.03)", "id >= 12 && id < 5132 && abs(z)<2", "colz"); - addTLines(o2::trk::SegmentationChip::PitchRowMLOT); - h = (TH2F*)gPad->GetPrimitive("h_dxH_vs_dzH_ML_z"); - Info("ML |z|<2", "RMS(dxH)=%.1f mu", h->GetRMS(2) * 1e4); - Info("ML |z|<2", "RMS(dzH)=%.1f mu", h->GetRMS(1) * 1e4); - canvdXdZHit->SaveAs("trkdigits_dxH_vs_dzH.pdf"); - canvdXdZHit->cd(5); - nt2->Draw("dxH:dzH>>h_dxH_vs_dzH_OT(300, -0.03, 0.03, 300, -0.03, 0.03)", "id >= 5132", "colz"); - addTLines(o2::trk::SegmentationChip::PitchRowMLOT); - h = (TH2F*)gPad->GetPrimitive("h_dxH_vs_dzH_OT"); - Info("OT", "RMS(dxH)=%.1f mu", h->GetRMS(2) * 1e4); - Info("OT", "RMS(dzH)=%.1f mu", h->GetRMS(1) * 1e4); - canvdXdZHit->cd(6); - nt2->Draw("dxH:dzH>>h_dxH_vs_dzH_OT_z(300, -0.03, 0.03, 300, -0.03, 0.03)", "id >= 5132 && abs(z)<2", "colz"); - addTLines(o2::trk::SegmentationChip::PitchRowMLOT); - h = (TH2F*)gPad->GetPrimitive("h_dxH_vs_dzH_OT_z"); - Info("OT |z|<2", "RMS(dxH)=%.1f mu", h->GetRMS(2) * 1e4); - Info("OT |z|<2", "RMS(dzH)=%.1f mu", h->GetRMS(1) * 1e4); - canvdXdZHit->SaveAs("trkdigits_dxH_vs_dzH.pdf"); - - f->Write(); - f->Close(); -} diff --git a/Detectors/Upgrades/ALICE3/TRK/macros/test/CheckDigitsTRK.C b/Detectors/Upgrades/ALICE3/TRK/macros/test/CheckDigitsTRK.C index f29a38f554200..400457fc98585 100644 --- a/Detectors/Upgrades/ALICE3/TRK/macros/test/CheckDigitsTRK.C +++ b/Detectors/Upgrades/ALICE3/TRK/macros/test/CheckDigitsTRK.C @@ -298,131 +298,131 @@ void CheckDigits(std::string digifile = "trkdigits.root", std::string hitfile = } // end loop on layers // digit maps in the xy and yz planes - auto canvXY = new TCanvas("canvXY", "", 1600, 2400); - canvXY->Divide(2, 3); - canvXY->cd(1); - nt->Draw("y:x >>h_y_vs_x_VD(1000, -3, 3, 1000, -3, 3)", "id < 12 ", "colz"); - canvXY->cd(2); - nt->Draw("y:z>>h_y_vs_z_VD(1000, -26, 26, 1000, -3, 3)", "id < 12 ", "colz"); - canvXY->cd(3); - nt->Draw("y:x>>h_y_vs_x_ML(1000, -25, 25, 1000, -25, 25)", "id >= 12 && id < 5132 ", "colz"); - canvXY->cd(4); - nt->Draw("y:z>>h_y_vs_z_ML(1000, -70, 70, 1000, -25, 25)", "id >= 12 && id < 5132 ", "colz"); - canvXY->cd(5); - nt->Draw("y:x>>h_y_vs_x_OT(1000, -85, 85, 1000, -85, 85)", "id >= 5132 ", "colz"); - canvXY->cd(6); - nt->Draw("y:z>>h_y_vs_z_OT(1000, -85, 85, 1000, -130, 130)", "id >= 5132 ", "colz"); - canvXY->SaveAs("trkdigits_y_vs_x_vs_z.pdf"); - - // z distributions - auto canvZ = new TCanvas("canvZ", "", 800, 2400); - canvZ->Divide(1, 3); - canvZ->cd(1); - nt->Draw("z>>h_z_VD(500, -26, 26)", "id < 12 "); - canvZ->cd(2); - nt->Draw("z>>h_z_ML(500, -70, 70)", "id >= 12 && id < 5132 "); - canvZ->cd(3); - nt->Draw("z>>h_z_OT(500, -85, 85)", "id >= 5132 "); - canvZ->SaveAs("trkdigits_z.pdf"); - - // dz distributions (difference between local position of digits and hits in x and z) - auto canvdZ = new TCanvas("canvdZ", "", 800, 2400); - canvdZ->Divide(1, 3); - canvdZ->cd(1); - nt->Draw("dz>>h_dz_VD(500, -0.05, 0.05)", "id < 12 "); - canvdZ->cd(2); - nt->Draw("dz>>h_dz_ML(500, -0.05, 0.05)", "id >= 12 && id < 5132 "); - canvdZ->cd(3); - nt->Draw("dz>>h_dz_OT(500, -0.05, 0.05)", "id >= 5132 "); - canvdZ->SaveAs("trkdigits_dz.pdf"); - canvdZ->SaveAs("trkdigits_dz.root"); - - // distributions of differences between local positions of digits and hits in x and z - auto canvdXdZ = new TCanvas("canvdXdZ", "", 1600, 2400); - canvdXdZ->Divide(2, 3); - canvdXdZ->cd(1); - nt->Draw("dx:dz>>h_dx_vs_dz_VD(500, -0.005, 0.005, 500, -0.005, 0.005)", "id < 12", "colz"); - addTLines(o2::trk::SegmentationChip::PitchRowVD); - auto h = (TH2F*)gPad->GetPrimitive("h_dx_vs_dz_VD"); - LOG(info) << "dx, dz"; - Info("VD", "RMS(dx)=%.1f mu", h->GetRMS(2) * 1e4); - Info("VD", "RMS(dz)=%.1f mu", h->GetRMS(1) * 1e4); - canvdXdZ->cd(2); - nt->Draw("dx:dz>>h_dx_vs_dz_VD_z(500, -0.005, 0.005, 500, -0.005, 0.005)", "id < 12 && abs(z)<0.5", "colz"); - addTLines(o2::trk::SegmentationChip::PitchRowVD); - h = (TH2F*)gPad->GetPrimitive("h_dx_vs_dz_VD_z"); - Info("VD |z|<1", "RMS(dx)=%.1f mu", h->GetRMS(2) * 1e4); - Info("VD |z|<1", "RMS(dz)=%.1f mu", h->GetRMS(1) * 1e4); - canvdXdZ->cd(3); - nt->Draw("dx:dz>>h_dx_vs_dz_ML(600, -0.03, 0.03, 600, -0.03, 0.03)", "id >= 12 && id < 5132", "colz"); - addTLines(o2::trk::SegmentationChip::PitchRowMLOT); - h = (TH2F*)gPad->GetPrimitive("h_dx_vs_dz_ML"); - Info("ML", "RMS(dx)=%.1f mu", h->GetRMS(2) * 1e4); - Info("ML", "RMS(dz)=%.1f mu", h->GetRMS(1) * 1e4); - canvdXdZ->cd(4); - nt->Draw("dx:dz>>h_dx_vs_dz_ML_z(600, -0.03, 0.03, 600, -0.03, 0.03)", "id >= 12 && id < 5132 && abs(z)<2", "colz"); - addTLines(o2::trk::SegmentationChip::PitchRowMLOT); - h = (TH2F*)gPad->GetPrimitive("h_dx_vs_dz_ML_z"); - Info("ML |z|<2", "RMS(dx)=%.1f mu", h->GetRMS(2) * 1e4); - Info("ML |z|<2", "RMS(dz)=%.1f mu", h->GetRMS(1) * 1e4); - canvdXdZ->SaveAs("trkdigits_dx_vs_dz.pdf"); - canvdXdZ->cd(5); - nt->Draw("dx:dz>>h_dx_vs_dz_OT(600, -0.03, 0.03, 600, -0.03, 0.03)", "id >= 5132", "colz"); - addTLines(o2::trk::SegmentationChip::PitchRowMLOT); - h = (TH2F*)gPad->GetPrimitive("h_dx_vs_dz_OT"); - Info("OT", "RMS(dx)=%.1f mu", h->GetRMS(2) * 1e4); - Info("OT", "RMS(dz)=%.1f mu", h->GetRMS(1) * 1e4); - canvdXdZ->cd(6); - nt->Draw("dx:dz>>h_dx_vs_dz_OT_z(600, -0.03, 0.03, 600, -0.03, 0.03)", "id >= 5132 && abs(z)<2", "colz"); - h = (TH2F*)gPad->GetPrimitive("h_dx_vs_dz_OT_z"); - addTLines(o2::trk::SegmentationChip::PitchRowMLOT); - Info("OT |z|<2", "RMS(dx)=%.1f mu", h->GetRMS(2) * 1e4); - Info("OT |z|<2", "RMS(dz)=%.1f mu", h->GetRMS(1) * 1e4); - canvdXdZ->SaveAs("trkdigits_dx_vs_dz.pdf"); - canvdXdZ->SaveAs("trkdigits_dx_vs_dz.root"); - - // distribution of differences between hit start and hit end in local coordinates - auto canvdXdZHit = new TCanvas("canvdXdZHit", "", 1600, 2400); - canvdXdZHit->Divide(2, 3); - canvdXdZHit->cd(1); - nt2->Draw("dxH:dzH>>h_dxH_vs_dzH_VD(300, -0.03, 0.03, 300, -0.03, 0.03)", "id < 12", "colz"); - addTLines(o2::trk::SegmentationChip::PitchRowVD); - LOG(info) << "dxH, dzH"; - h = (TH2F*)gPad->GetPrimitive("h_dxH_vs_dzH_VD"); - Info("VD", "RMS(dxH)=%.1f mu", h->GetRMS(2) * 1e4); - Info("VD", "RMS(dzH)=%.1f mu", h->GetRMS(1) * 1e4); - canvdXdZHit->cd(2); - nt2->Draw("dxH:dzH>>h_dxH_vs_dzH_VD_z(300, -0.03, 0.03, 300, -0.03, 0.03)", "id < 12 && abs(z)<2", "colz"); - addTLines(o2::trk::SegmentationChip::PitchRowVD); - h = (TH2F*)gPad->GetPrimitive("h_dxH_vs_dzH_VD_z"); - Info("VD |z|<2", "RMS(dxH)=%.1f mu", h->GetRMS(2) * 1e4); - Info("VD |z|<2", "RMS(dzH)=%.1f mu", h->GetRMS(1) * 1e4); - canvdXdZHit->cd(3); - nt2->Draw("dxH:dzH>>h_dxH_vs_dzH_ML(300, -0.03, 0.03, 300, -0.03, 0.03)", "id >= 12 && id < 5132", "colz"); - addTLines(o2::trk::SegmentationChip::PitchRowMLOT); - h = (TH2F*)gPad->GetPrimitive("h_dxH_vs_dzH_ML"); - Info("ML", "RMS(dxH)=%.1f mu", h->GetRMS(2) * 1e4); - Info("ML", "RMS(dzH)=%.1f mu", h->GetRMS(1) * 1e4); - canvdXdZHit->cd(4); - nt2->Draw("dxH:dzH>>h_dxH_vs_dzH_ML_z(300, -0.03, 0.03, 300, -0.03, 0.03)", "id >= 12 && id < 5132 && abs(z)<2", "colz"); - addTLines(o2::trk::SegmentationChip::PitchRowMLOT); - h = (TH2F*)gPad->GetPrimitive("h_dxH_vs_dzH_ML_z"); - Info("ML |z|<2", "RMS(dxH)=%.1f mu", h->GetRMS(2) * 1e4); - Info("ML |z|<2", "RMS(dzH)=%.1f mu", h->GetRMS(1) * 1e4); - canvdXdZHit->SaveAs("trkdigits_dxH_vs_dzH.pdf"); - canvdXdZHit->cd(5); - nt2->Draw("dxH:dzH>>h_dxH_vs_dzH_OT(300, -0.03, 0.03, 300, -0.03, 0.03)", "id >= 5132", "colz"); - addTLines(o2::trk::SegmentationChip::PitchRowMLOT); - h = (TH2F*)gPad->GetPrimitive("h_dxH_vs_dzH_OT"); - Info("OT", "RMS(dxH)=%.1f mu", h->GetRMS(2) * 1e4); - Info("OT", "RMS(dzH)=%.1f mu", h->GetRMS(1) * 1e4); - canvdXdZHit->cd(6); - nt2->Draw("dxH:dzH>>h_dxH_vs_dzH_OT_z(300, -0.03, 0.03, 300, -0.03, 0.03)", "id >= 5132 && abs(z)<2", "colz"); - addTLines(o2::trk::SegmentationChip::PitchRowMLOT); - h = (TH2F*)gPad->GetPrimitive("h_dxH_vs_dzH_OT_z"); - Info("OT |z|<2", "RMS(dxH)=%.1f mu", h->GetRMS(2) * 1e4); - Info("OT |z|<2", "RMS(dzH)=%.1f mu", h->GetRMS(1) * 1e4); - canvdXdZHit->SaveAs("trkdigits_dxH_vs_dzH.pdf"); + auto canvXY = new TCanvas("canvXY", "", 1600, 2400); + canvXY->Divide(2, 3); + canvXY->cd(1); + nt->Draw("y:x >>h_y_vs_x_VD(1000, -3, 3, 1000, -3, 3)", "id < 12 ", "colz"); + canvXY->cd(2); + nt->Draw("y:z>>h_y_vs_z_VD(1000, -26, 26, 1000, -3, 3)", "id < 12 ", "colz"); + canvXY->cd(3); + nt->Draw("y:x>>h_y_vs_x_ML(1000, -25, 25, 1000, -25, 25)", "id >= 12 && id < 5132 ", "colz"); + canvXY->cd(4); + nt->Draw("y:z>>h_y_vs_z_ML(1000, -70, 70, 1000, -25, 25)", "id >= 12 && id < 5132 ", "colz"); + canvXY->cd(5); + nt->Draw("y:x>>h_y_vs_x_OT(1000, -85, 85, 1000, -85, 85)", "id >= 5132 ", "colz"); + canvXY->cd(6); + nt->Draw("y:z>>h_y_vs_z_OT(1000, -85, 85, 1000, -130, 130)", "id >= 5132 ", "colz"); + canvXY->SaveAs("trkdigits_y_vs_x_vs_z.pdf"); + + // z distributions + auto canvZ = new TCanvas("canvZ", "", 800, 2400); + canvZ->Divide(1, 3); + canvZ->cd(1); + nt->Draw("z>>h_z_VD(500, -26, 26)", "id < 12 "); + canvZ->cd(2); + nt->Draw("z>>h_z_ML(500, -70, 70)", "id >= 12 && id < 5132 "); + canvZ->cd(3); + nt->Draw("z>>h_z_OT(500, -85, 85)", "id >= 5132 "); + canvZ->SaveAs("trkdigits_z.pdf"); + + // dz distributions (difference between local position of digits and hits in x and z) + auto canvdZ = new TCanvas("canvdZ", "", 800, 2400); + canvdZ->Divide(1, 3); + canvdZ->cd(1); + nt->Draw("dz>>h_dz_VD(500, -0.05, 0.05)", "id < 12 "); + canvdZ->cd(2); + nt->Draw("dz>>h_dz_ML(500, -0.05, 0.05)", "id >= 12 && id < 5132 "); + canvdZ->cd(3); + nt->Draw("dz>>h_dz_OT(500, -0.05, 0.05)", "id >= 5132 "); + canvdZ->SaveAs("trkdigits_dz.pdf"); + canvdZ->SaveAs("trkdigits_dz.root"); + + // distributions of differences between local positions of digits and hits in x and z + auto canvdXdZ = new TCanvas("canvdXdZ", "", 1600, 2400); + canvdXdZ->Divide(2, 3); + canvdXdZ->cd(1); + nt->Draw("dx:dz>>h_dx_vs_dz_VD(500, -0.005, 0.005, 500, -0.005, 0.005)", "id < 12", "colz"); + addTLines(o2::trk::SegmentationChip::PitchRowVD); + auto h = (TH2F*)gPad->GetPrimitive("h_dx_vs_dz_VD"); + LOG(info) << "dx, dz"; + Info("VD", "RMS(dx)=%.1f mu", h->GetRMS(2) * 1e4); + Info("VD", "RMS(dz)=%.1f mu", h->GetRMS(1) * 1e4); + canvdXdZ->cd(2); + nt->Draw("dx:dz>>h_dx_vs_dz_VD_z(500, -0.005, 0.005, 500, -0.005, 0.005)", "id < 12 && abs(z)<0.5", "colz"); + addTLines(o2::trk::SegmentationChip::PitchRowVD); + h = (TH2F*)gPad->GetPrimitive("h_dx_vs_dz_VD_z"); + Info("VD |z|<1", "RMS(dx)=%.1f mu", h->GetRMS(2) * 1e4); + Info("VD |z|<1", "RMS(dz)=%.1f mu", h->GetRMS(1) * 1e4); + canvdXdZ->cd(3); + nt->Draw("dx:dz>>h_dx_vs_dz_ML(600, -0.03, 0.03, 600, -0.03, 0.03)", "id >= 12 && id < 5132", "colz"); + addTLines(o2::trk::SegmentationChip::PitchRowMLOT); + h = (TH2F*)gPad->GetPrimitive("h_dx_vs_dz_ML"); + Info("ML", "RMS(dx)=%.1f mu", h->GetRMS(2) * 1e4); + Info("ML", "RMS(dz)=%.1f mu", h->GetRMS(1) * 1e4); + canvdXdZ->cd(4); + nt->Draw("dx:dz>>h_dx_vs_dz_ML_z(600, -0.03, 0.03, 600, -0.03, 0.03)", "id >= 12 && id < 5132 && abs(z)<2", "colz"); + addTLines(o2::trk::SegmentationChip::PitchRowMLOT); + h = (TH2F*)gPad->GetPrimitive("h_dx_vs_dz_ML_z"); + Info("ML |z|<2", "RMS(dx)=%.1f mu", h->GetRMS(2) * 1e4); + Info("ML |z|<2", "RMS(dz)=%.1f mu", h->GetRMS(1) * 1e4); + canvdXdZ->SaveAs("trkdigits_dx_vs_dz.pdf"); + canvdXdZ->cd(5); + nt->Draw("dx:dz>>h_dx_vs_dz_OT(600, -0.03, 0.03, 600, -0.03, 0.03)", "id >= 5132", "colz"); + addTLines(o2::trk::SegmentationChip::PitchRowMLOT); + h = (TH2F*)gPad->GetPrimitive("h_dx_vs_dz_OT"); + Info("OT", "RMS(dx)=%.1f mu", h->GetRMS(2) * 1e4); + Info("OT", "RMS(dz)=%.1f mu", h->GetRMS(1) * 1e4); + canvdXdZ->cd(6); + nt->Draw("dx:dz>>h_dx_vs_dz_OT_z(600, -0.03, 0.03, 600, -0.03, 0.03)", "id >= 5132 && abs(z)<2", "colz"); + h = (TH2F*)gPad->GetPrimitive("h_dx_vs_dz_OT_z"); + addTLines(o2::trk::SegmentationChip::PitchRowMLOT); + Info("OT |z|<2", "RMS(dx)=%.1f mu", h->GetRMS(2) * 1e4); + Info("OT |z|<2", "RMS(dz)=%.1f mu", h->GetRMS(1) * 1e4); + canvdXdZ->SaveAs("trkdigits_dx_vs_dz.pdf"); + canvdXdZ->SaveAs("trkdigits_dx_vs_dz.root"); + + // distribution of differences between hit start and hit end in local coordinates + auto canvdXdZHit = new TCanvas("canvdXdZHit", "", 1600, 2400); + canvdXdZHit->Divide(2, 3); + canvdXdZHit->cd(1); + nt2->Draw("dxH:dzH>>h_dxH_vs_dzH_VD(300, -0.03, 0.03, 300, -0.03, 0.03)", "id < 12", "colz"); + addTLines(o2::trk::SegmentationChip::PitchRowVD); + LOG(info) << "dxH, dzH"; + h = (TH2F*)gPad->GetPrimitive("h_dxH_vs_dzH_VD"); + Info("VD", "RMS(dxH)=%.1f mu", h->GetRMS(2) * 1e4); + Info("VD", "RMS(dzH)=%.1f mu", h->GetRMS(1) * 1e4); + canvdXdZHit->cd(2); + nt2->Draw("dxH:dzH>>h_dxH_vs_dzH_VD_z(300, -0.03, 0.03, 300, -0.03, 0.03)", "id < 12 && abs(z)<2", "colz"); + addTLines(o2::trk::SegmentationChip::PitchRowVD); + h = (TH2F*)gPad->GetPrimitive("h_dxH_vs_dzH_VD_z"); + Info("VD |z|<2", "RMS(dxH)=%.1f mu", h->GetRMS(2) * 1e4); + Info("VD |z|<2", "RMS(dzH)=%.1f mu", h->GetRMS(1) * 1e4); + canvdXdZHit->cd(3); + nt2->Draw("dxH:dzH>>h_dxH_vs_dzH_ML(300, -0.03, 0.03, 300, -0.03, 0.03)", "id >= 12 && id < 5132", "colz"); + addTLines(o2::trk::SegmentationChip::PitchRowMLOT); + h = (TH2F*)gPad->GetPrimitive("h_dxH_vs_dzH_ML"); + Info("ML", "RMS(dxH)=%.1f mu", h->GetRMS(2) * 1e4); + Info("ML", "RMS(dzH)=%.1f mu", h->GetRMS(1) * 1e4); + canvdXdZHit->cd(4); + nt2->Draw("dxH:dzH>>h_dxH_vs_dzH_ML_z(300, -0.03, 0.03, 300, -0.03, 0.03)", "id >= 12 && id < 5132 && abs(z)<2", "colz"); + addTLines(o2::trk::SegmentationChip::PitchRowMLOT); + h = (TH2F*)gPad->GetPrimitive("h_dxH_vs_dzH_ML_z"); + Info("ML |z|<2", "RMS(dxH)=%.1f mu", h->GetRMS(2) * 1e4); + Info("ML |z|<2", "RMS(dzH)=%.1f mu", h->GetRMS(1) * 1e4); + canvdXdZHit->SaveAs("trkdigits_dxH_vs_dzH.pdf"); + canvdXdZHit->cd(5); + nt2->Draw("dxH:dzH>>h_dxH_vs_dzH_OT(300, -0.03, 0.03, 300, -0.03, 0.03)", "id >= 5132", "colz"); + addTLines(o2::trk::SegmentationChip::PitchRowMLOT); + h = (TH2F*)gPad->GetPrimitive("h_dxH_vs_dzH_OT"); + Info("OT", "RMS(dxH)=%.1f mu", h->GetRMS(2) * 1e4); + Info("OT", "RMS(dzH)=%.1f mu", h->GetRMS(1) * 1e4); + canvdXdZHit->cd(6); + nt2->Draw("dxH:dzH>>h_dxH_vs_dzH_OT_z(300, -0.03, 0.03, 300, -0.03, 0.03)", "id >= 5132 && abs(z)<2", "colz"); + addTLines(o2::trk::SegmentationChip::PitchRowMLOT); + h = (TH2F*)gPad->GetPrimitive("h_dxH_vs_dzH_OT_z"); + Info("OT |z|<2", "RMS(dxH)=%.1f mu", h->GetRMS(2) * 1e4); + Info("OT |z|<2", "RMS(dzH)=%.1f mu", h->GetRMS(1) * 1e4); + canvdXdZHit->SaveAs("trkdigits_dxH_vs_dzH.pdf"); f->Write(); f->Close(); diff --git a/Detectors/Upgrades/ALICE3/TRK/reconstruction/include/TRKReconstruction/Clusterer.h b/Detectors/Upgrades/ALICE3/TRK/reconstruction/include/TRKReconstruction/Clusterer.h index 70518b2ace593..bcd95155f533f 100644 --- a/Detectors/Upgrades/ALICE3/TRK/reconstruction/include/TRKReconstruction/Clusterer.h +++ b/Detectors/Upgrades/ALICE3/TRK/reconstruction/include/TRKReconstruction/Clusterer.h @@ -48,7 +48,6 @@ class Clusterer using Digit = o2::itsmft::Digit; using DigROFRecord = o2::itsmft::ROFRecord; - using DigMC2ROFRecord = o2::itsmft::MC2ROFRecord; using ClusterTruth = o2::dataformats::MCTruthContainer; using ConstDigitTruth = o2::dataformats::ConstMCTruthContainerView; using Label = o2::MCCompLabel; @@ -167,9 +166,7 @@ class Clusterer std::vector& patterns, std::vector& clusterROFs, const ConstDigitTruth* digitLabels = nullptr, - ClusterTruth* clusterLabels = nullptr, - gsl::span digMC2ROFs = {}, - std::vector* clusterMC2ROFs = nullptr); + ClusterTruth* clusterLabels = nullptr); protected: int mNHugeClus = 0; diff --git a/Detectors/Upgrades/ALICE3/TRK/reconstruction/include/TRKReconstruction/ClustererACTS.h b/Detectors/Upgrades/ALICE3/TRK/reconstruction/include/TRKReconstruction/ClustererACTS.h index 37a148aa78afb..5d68193e5e375 100644 --- a/Detectors/Upgrades/ALICE3/TRK/reconstruction/include/TRKReconstruction/ClustererACTS.h +++ b/Detectors/Upgrades/ALICE3/TRK/reconstruction/include/TRKReconstruction/ClustererACTS.h @@ -35,9 +35,7 @@ class ClustererACTS : public Clusterer std::vector& patterns, std::vector& clusterROFs, const ConstDigitTruth* digitLabels = nullptr, - ClusterTruth* clusterLabels = nullptr, - gsl::span digMC2ROFs = {}, - std::vector* clusterMC2ROFs = nullptr) override; + ClusterTruth* clusterLabels = nullptr) override; private: }; diff --git a/Detectors/Upgrades/ALICE3/TRK/reconstruction/src/Clusterer.cxx b/Detectors/Upgrades/ALICE3/TRK/reconstruction/src/Clusterer.cxx index bdaa76319c1f2..e0d689e4db5ed 100644 --- a/Detectors/Upgrades/ALICE3/TRK/reconstruction/src/Clusterer.cxx +++ b/Detectors/Upgrades/ALICE3/TRK/reconstruction/src/Clusterer.cxx @@ -28,9 +28,7 @@ void Clusterer::process(gsl::span digits, std::vector& patterns, std::vector& clusterROFs, const ConstDigitTruth* digitLabels, - ClusterTruth* clusterLabels, - gsl::span digMC2ROFs, - std::vector* clusterMC2ROFs) + ClusterTruth* clusterLabels) { if (!mThread) { mThread = std::make_unique(this); @@ -81,13 +79,6 @@ void Clusterer::process(gsl::span digits, clusterROFs.emplace_back(inROF.getBCData(), inROF.getROFrame(), outFirst, static_cast(clusters.size()) - outFirst); } - - if (clusterMC2ROFs && !digMC2ROFs.empty()) { - clusterMC2ROFs->reserve(clusterMC2ROFs->size() + digMC2ROFs.size()); - for (const auto& in : digMC2ROFs) { - clusterMC2ROFs->emplace_back(in.eventRecordID, in.rofRecordID, in.minROF, in.maxROF); - } - } } //__________________________________________________ diff --git a/Detectors/Upgrades/ALICE3/TRK/reconstruction/src/ClustererACTS.cxx b/Detectors/Upgrades/ALICE3/TRK/reconstruction/src/ClustererACTS.cxx index 2dbf56ae610e3..b764fcdd1cd79 100644 --- a/Detectors/Upgrades/ALICE3/TRK/reconstruction/src/ClustererACTS.cxx +++ b/Detectors/Upgrades/ALICE3/TRK/reconstruction/src/ClustererACTS.cxx @@ -162,9 +162,7 @@ void ClustererACTS::process(gsl::span digits, std::vector& patterns, std::vector& clusterROFs, const ConstDigitTruth* digitLabels, - ClusterTruth* clusterLabels, - gsl::span digMC2ROFs, - std::vector* clusterMC2ROFs) + ClusterTruth* clusterLabels) { if (!mThread) { mThread = std::make_unique(this); diff --git a/Detectors/Upgrades/ALICE3/TRK/simulation/include/TRKSimulation/DigiParams.h b/Detectors/Upgrades/ALICE3/TRK/simulation/include/TRKSimulation/DigiParams.h index 0253f8bf2bb5d..d7d1ea28bfcf7 100644 --- a/Detectors/Upgrades/ALICE3/TRK/simulation/include/TRKSimulation/DigiParams.h +++ b/Detectors/Upgrades/ALICE3/TRK/simulation/include/TRKSimulation/DigiParams.h @@ -126,7 +126,7 @@ class DigiParams std::unique_ptr mResponse; //!< pointer on external response // auxiliary precalculated parameters - std::array mROFrameLayerLengthInv; ///< inverse length of RO frame in ns per layer + std::array mROFrameLayerLengthInv; ///< inverse length of RO frame in ns per layer // ClassDef(DigiParams, 2); }; diff --git a/Detectors/Upgrades/ALICE3/TRK/workflow/include/TRKWorkflow/ClustererSpec.h b/Detectors/Upgrades/ALICE3/TRK/workflow/include/TRKWorkflow/ClustererSpec.h index 18cc6d245025a..9d072e85d574a 100644 --- a/Detectors/Upgrades/ALICE3/TRK/workflow/include/TRKWorkflow/ClustererSpec.h +++ b/Detectors/Upgrades/ALICE3/TRK/workflow/include/TRKWorkflow/ClustererSpec.h @@ -14,6 +14,7 @@ #include "Framework/DataProcessorSpec.h" #include "Framework/Task.h" +#include "TRKBase/AlmiraParam.h" #include "TRKReconstruction/Clusterer.h" #ifdef O2_WITH_ACTS #include "TRKReconstruction/ClustererACTS.h" @@ -30,6 +31,7 @@ class ClustererDPL : public o2::framework::Task void run(o2::framework::ProcessingContext& pc) final; private: + static constexpr int mLayers = o2::trk::AlmiraParam::kNLayers; bool mUseMC = true; int mNThreads = 1; o2::trk::Clusterer mClusterer; diff --git a/Detectors/Upgrades/ALICE3/TRK/workflow/include/TRKWorkflow/DigitReaderSpec.h b/Detectors/Upgrades/ALICE3/TRK/workflow/include/TRKWorkflow/DigitReaderSpec.h index 4acea7bd74914..92b64e0815cfb 100644 --- a/Detectors/Upgrades/ALICE3/TRK/workflow/include/TRKWorkflow/DigitReaderSpec.h +++ b/Detectors/Upgrades/ALICE3/TRK/workflow/include/TRKWorkflow/DigitReaderSpec.h @@ -45,16 +45,16 @@ class DigitReader : public Task protected: void connectTree(const std::string& filename); - template - void setBranchAddress(const std::string& base, Ptr& addr, int layer = -1); - std::string getBranchName(const std::string& base, int index) const; + template + void setBranchAddress(const std::string& base, Ptr& addr, int layer = -1); + std::string getBranchName(const std::string& base, int index) const; - static constexpr int mLayers = o2::trk::AlmiraParam::kNLayers; + static constexpr int mLayers = o2::trk::AlmiraParam::kNLayers; - std::vector*> mDigits{nullptr}; + std::vector*> mDigits{nullptr}; std::vector mCalib, *mCalibPtr = &mCalib; - std::vector*> mDigROFRec{nullptr}; - std::vector mPLabels{nullptr}; + std::vector*> mDigROFRec{nullptr}; + std::vector mPLabels{nullptr}; o2::header::DataOrigin mOrigin = o2::header::gDataOriginInvalid; diff --git a/Detectors/Upgrades/ALICE3/TRK/workflow/src/ClusterWriterSpec.cxx b/Detectors/Upgrades/ALICE3/TRK/workflow/src/ClusterWriterSpec.cxx index bc3a75c646198..863915bac0572 100644 --- a/Detectors/Upgrades/ALICE3/TRK/workflow/src/ClusterWriterSpec.cxx +++ b/Detectors/Upgrades/ALICE3/TRK/workflow/src/ClusterWriterSpec.cxx @@ -11,9 +11,16 @@ /// @file ClusterWriterSpec.cxx +#include +#include +#include #include +#include #include "TRKWorkflow/ClusterWriterSpec.h" +#include "Framework/ConcreteDataMatcher.h" +#include "Framework/DataRef.h" +#include "TRKBase/AlmiraParam.h" #include "DPLUtils/MakeRootTreeWriterSpec.h" #include "DataFormatsTRK/Cluster.h" #include "DataFormatsTRK/ROFRecord.h" @@ -35,31 +42,68 @@ using ROFRecLblType = std::vector; DataProcessorSpec getClusterWriterSpec(bool useMC) { - auto clustersSize = std::make_shared(0); - auto clustersSizeGetter = [clustersSize](ClustersType const& clusters) { - *clustersSize = clusters.size(); + static constexpr o2::header::DataOrigin Origin{o2::header::gDataOriginTRK}; + static constexpr int nLayers = o2::trk::AlmiraParam::kNLayers; + const auto detName = Origin.as(); + + auto compClusterSizes = std::make_shared>(nLayers, 0); + auto compClustersSizeGetter = [compClusterSizes](ClustersType const& compClusters, DataRef const& ref) { + auto const* dh = DataRefUtils::getHeader(ref); + (*compClusterSizes)[dh->subSpecification] = compClusters.size(); + }; + auto logger = [detName, compClusterSizes](ROFrameType const& rofs, DataRef const& ref) { + auto const* dh = DataRefUtils::getHeader(ref); + const auto i = dh->subSpecification; + LOG(info) << detName << "ClusterWriter on layer " << i + << " pulled " << (*compClusterSizes)[i] << " clusters, in " << rofs.size() << " RO frames"; }; - auto logger = [clustersSize](ROFrameType const& rofs) { - LOG(info) << "TRKClusterWriter pulled " << *clustersSize << " clusters, in " << rofs.size() << " RO frames"; + auto getIndex = [](DataRef const& ref) -> size_t { + auto const* dh = DataRefUtils::getHeader(ref); + return static_cast(dh->subSpecification); }; + auto getName = [](std::string base, size_t index) -> std::string { + return base + "_" + std::to_string(index); + }; + auto detNameLC = detName; + std::transform(detNameLC.begin(), detNameLC.end(), detNameLC.begin(), [](unsigned char c) { return std::tolower(c); }); + + std::vector vecInpSpecClus, vecInpSpecPatt, vecInpSpecROF, vecInpSpecLbl; + vecInpSpecClus.reserve(nLayers); + vecInpSpecPatt.reserve(nLayers); + vecInpSpecROF.reserve(nLayers); + vecInpSpecLbl.reserve(nLayers); + for (int iLayer = 0; iLayer < nLayers; iLayer++) { + vecInpSpecClus.emplace_back(getName("compclus", iLayer), Origin, "COMPCLUSTERS", iLayer); + vecInpSpecPatt.emplace_back(getName("patterns", iLayer), Origin, "PATTERNS", iLayer); + vecInpSpecROF.emplace_back(getName("ROframes", iLayer), Origin, "CLUSTERSROF", iLayer); + vecInpSpecLbl.emplace_back(getName("labels", iLayer), Origin, "CLUSTERSMCTR", iLayer); + } - return MakeRootTreeWriterSpec("trk-cluster-writer", + return MakeRootTreeWriterSpec(std::format("{}-cluster-writer", detNameLC).c_str(), "o2clus_trk.root", - MakeRootTreeWriterSpec::TreeAttributes{"o2sim", "Tree with TRK clusters"}, - BranchDefinition{InputSpec{"compclus", "TRK", "COMPCLUSTERS", 0}, - "TRKClusterComp", - clustersSizeGetter}, - BranchDefinition{InputSpec{"patterns", "TRK", "PATTERNS", 0}, - "TRKClusterPatt"}, - BranchDefinition{InputSpec{"ROframes", "TRK", "CLUSTERSROF", 0}, - "TRKClustersROF", - logger}, - BranchDefinition{InputSpec{"labels", "TRK", "CLUSTERSMCTR", 0}, - "TRKClusterMCTruth", - (useMC ? 1 : 0)}, - BranchDefinition{InputSpec{"MC2ROframes", "TRK", "CLUSTERSMC2ROF", 0}, - "TRKClustersMC2ROF", - (useMC ? 1 : 0)})(); + MakeRootTreeWriterSpec::TreeAttributes{.name = "o2sim", .title = "Tree with TRK clusters"}, + BranchDefinition{vecInpSpecClus, + "TRKClusterComp", "compact-cluster-branch", + nLayers, + compClustersSizeGetter, + getIndex, + getName}, + BranchDefinition{vecInpSpecPatt, + "TRKClusterPatt", "cluster-pattern-branch", + nLayers, + getIndex, + getName}, + BranchDefinition{vecInpSpecROF, + "TRKClustersROF", "cluster-rof-branch", + nLayers, + logger, + getIndex, + getName}, + BranchDefinition{vecInpSpecLbl, + "TRKClusterMCTruth", "cluster-label-branch", + (useMC ? nLayers : 0), + getIndex, + getName})(); } } // namespace o2::trk diff --git a/Detectors/Upgrades/ALICE3/TRK/workflow/src/ClustererSpec.cxx b/Detectors/Upgrades/ALICE3/TRK/workflow/src/ClustererSpec.cxx index 5d9ac463b3f54..f91262e021a55 100644 --- a/Detectors/Upgrades/ALICE3/TRK/workflow/src/ClustererSpec.cxx +++ b/Detectors/Upgrades/ALICE3/TRK/workflow/src/ClustererSpec.cxx @@ -17,6 +17,8 @@ #include "Framework/Logger.h" #include "SimulationDataFormat/ConstMCTruthContainer.h" +#include + namespace o2::trk { @@ -30,82 +32,84 @@ void ClustererDPL::init(o2::framework::InitContext& ic) void ClustererDPL::run(o2::framework::ProcessingContext& pc) { - auto digits = pc.inputs().get>("digits"); - auto rofs = pc.inputs().get>("ROframes"); + o2::base::GeometryManager::loadGeometry("sgn_geometry.root", false, true); - gsl::span mc2rofs; - gsl::span labelbuffer; - if (mUseMC) { - labelbuffer = pc.inputs().get>("labels"); - mc2rofs = pc.inputs().get>("MC2ROframes"); - } - o2::dataformats::ConstMCTruthContainerView labels(labelbuffer); + uint64_t totalClusters = 0; + for (int iLayer = 0; iLayer < mLayers; ++iLayer) { + auto digits = pc.inputs().get>(std::format("digits_{}", iLayer)); + auto rofs = pc.inputs().get>(std::format("ROframes_{}", iLayer)); - std::vector clusters; - std::vector patterns; - std::vector clusterROFs; - std::unique_ptr> clusterLabels; - std::vector clusterMC2ROFs; - if (mUseMC) { - clusterLabels = std::make_unique>(); - } - o2::base::GeometryManager::loadGeometry("o2sim_geometry.root", false, true); + gsl::span labelbuffer; + if (mUseMC) { + labelbuffer = pc.inputs().get>(std::format("labels_{}", iLayer)); + } + o2::dataformats::ConstMCTruthContainerView labels(labelbuffer); + + std::vector clusters; + std::vector patterns; + std::vector clusterROFs; + std::unique_ptr> clusterLabels; + if (mUseMC) { + clusterLabels = std::make_unique>(); + } #ifdef O2_WITH_ACTS - if (mUseACTS) { - LOG(info) << "Running TRKClusterer with ACTS"; - mClustererACTS.process(digits, - rofs, - clusters, - patterns, - clusterROFs, - mUseMC ? &labels : nullptr, - clusterLabels.get(), - mc2rofs, - mUseMC ? &clusterMC2ROFs : nullptr); - } else + if (mUseACTS) { + LOG(info) << "Running TRKClusterer with ACTS on layer " << iLayer; + mClustererACTS.process(digits, + rofs, + clusters, + patterns, + clusterROFs, + mUseMC ? &labels : nullptr, + clusterLabels.get()); + } else #endif - { - LOG(info) << "Running TRKClusterer"; - mClusterer.process(digits, - rofs, - clusters, - patterns, - clusterROFs, - mUseMC ? &labels : nullptr, - clusterLabels.get(), - mc2rofs, - mUseMC ? &clusterMC2ROFs : nullptr); - } - - pc.outputs().snapshot(o2::framework::Output{"TRK", "COMPCLUSTERS", 0}, clusters); - pc.outputs().snapshot(o2::framework::Output{"TRK", "PATTERNS", 0}, patterns); - pc.outputs().snapshot(o2::framework::Output{"TRK", "CLUSTERSROF", 0}, clusterROFs); + { + LOG(info) << "Running TRKClusterer on layer " << iLayer; + mClusterer.process(digits, + rofs, + clusters, + patterns, + clusterROFs, + mUseMC ? &labels : nullptr, + clusterLabels.get()); + } - if (mUseMC) { - pc.outputs().snapshot(o2::framework::Output{"TRK", "CLUSTERSMCTR", 0}, *clusterLabels); - pc.outputs().snapshot(o2::framework::Output{"TRK", "CLUSTERSMC2ROF", 0}, clusterMC2ROFs); + const auto subspec = static_cast(iLayer); + pc.outputs().snapshot(o2::framework::Output{"TRK", "COMPCLUSTERS", subspec}, clusters); + pc.outputs().snapshot(o2::framework::Output{"TRK", "PATTERNS", subspec}, patterns); + pc.outputs().snapshot(o2::framework::Output{"TRK", "CLUSTERSROF", subspec}, clusterROFs); + if (mUseMC) { + pc.outputs().snapshot(o2::framework::Output{"TRK", "CLUSTERSMCTR", subspec}, *clusterLabels); + } + totalClusters += clusters.size(); + LOGP(info, "TRKClusterer layer {} pushed {} clusters in {} ROFs", iLayer, clusters.size(), clusterROFs.size()); } - LOGP(info, "TRKClusterer pushed {} clusters in {} ROFs", clusters.size(), clusterROFs.size()); + LOGP(info, "TRKClusterer produced {} clusters", totalClusters); } o2::framework::DataProcessorSpec getClustererSpec(bool useMC) { + static constexpr int nLayers = o2::trk::AlmiraParam::kNLayers; std::vector inputs; - inputs.emplace_back("digits", "TRK", "DIGITS", 0, o2::framework::Lifetime::Timeframe); - inputs.emplace_back("ROframes", "TRK", "DIGITSROF", 0, o2::framework::Lifetime::Timeframe); + for (int iLayer = 0; iLayer < nLayers; ++iLayer) { + inputs.emplace_back(std::format("digits_{}", iLayer), "TRK", "DIGITS", iLayer, o2::framework::Lifetime::Timeframe); + inputs.emplace_back(std::format("ROframes_{}", iLayer), "TRK", "DIGITSROF", iLayer, o2::framework::Lifetime::Timeframe); + if (useMC) { + inputs.emplace_back(std::format("labels_{}", iLayer), "TRK", "DIGITSMCTR", iLayer, o2::framework::Lifetime::Timeframe); + } + } std::vector outputs; - outputs.emplace_back("TRK", "COMPCLUSTERS", 0, o2::framework::Lifetime::Timeframe); - outputs.emplace_back("TRK", "PATTERNS", 0, o2::framework::Lifetime::Timeframe); - outputs.emplace_back("TRK", "CLUSTERSROF", 0, o2::framework::Lifetime::Timeframe); - - if (useMC) { - inputs.emplace_back("labels", "TRK", "DIGITSMCTR", 0, o2::framework::Lifetime::Timeframe); - inputs.emplace_back("MC2ROframes", "TRK", "DIGITSMC2ROF", 0, o2::framework::Lifetime::Timeframe); - outputs.emplace_back("TRK", "CLUSTERSMCTR", 0, o2::framework::Lifetime::Timeframe); - outputs.emplace_back("TRK", "CLUSTERSMC2ROF", 0, o2::framework::Lifetime::Timeframe); + for (int iLayer = 0; iLayer < nLayers; ++iLayer) { + outputs.emplace_back("TRK", "COMPCLUSTERS", iLayer, o2::framework::Lifetime::Timeframe); + outputs.emplace_back("TRK", "PATTERNS", iLayer, o2::framework::Lifetime::Timeframe); + outputs.emplace_back("TRK", "CLUSTERSROF", iLayer, o2::framework::Lifetime::Timeframe); + if (useMC) { + outputs.emplace_back("TRK", "CLUSTERSMCTR", iLayer, o2::framework::Lifetime::Timeframe); + } } return o2::framework::DataProcessorSpec{ From 135bdf737a12442a7f9d6915939f4124292667ea Mon Sep 17 00:00:00 2001 From: Maximiliano Puccio Date: Mon, 27 Apr 2026 22:41:53 +0200 Subject: [PATCH 5/6] Change macro name in CMakeLists.txt --- Detectors/Upgrades/ALICE3/TRK/macros/test/CMakeLists.txt | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Detectors/Upgrades/ALICE3/TRK/macros/test/CMakeLists.txt b/Detectors/Upgrades/ALICE3/TRK/macros/test/CMakeLists.txt index 33d1b4a5afdc6..cdae7c9c379fd 100644 --- a/Detectors/Upgrades/ALICE3/TRK/macros/test/CMakeLists.txt +++ b/Detectors/Upgrades/ALICE3/TRK/macros/test/CMakeLists.txt @@ -21,7 +21,7 @@ o2_add_test_root_macro(CheckBandwidth.C O2::Steer LABELS trk COMPILE_ONLY) -o2_add_test_root_macro(CheckDigits.C +o2_add_test_root_macro(CheckDigitsTRK.C PUBLIC_LINK_LIBRARIES O2::ITSMFTBase O2::ITSMFTSimulation O2::TRKBase From 0a64b97ff77001b1bcf6ec3a004abd7850e0ffdd Mon Sep 17 00:00:00 2001 From: Maximiliano Puccio Date: Wed, 29 Apr 2026 09:27:54 +0200 Subject: [PATCH 6/6] Fix Specs.h and use that in geometry building --- .../ALICE3/TRK/base/include/TRKBase/Specs.h | 4 ++-- .../Upgrades/ALICE3/TRK/base/src/GeometryTGeo.cxx | 7 ++++--- .../Upgrades/ALICE3/TRK/simulation/src/Detector.cxx | 13 +++++++------ 3 files changed, 13 insertions(+), 11 deletions(-) diff --git a/Detectors/Upgrades/ALICE3/TRK/base/include/TRKBase/Specs.h b/Detectors/Upgrades/ALICE3/TRK/base/include/TRKBase/Specs.h index 6be2e0bd0f827..0ed7ca6a8a8d4 100644 --- a/Detectors/Upgrades/ALICE3/TRK/base/include/TRKBase/Specs.h +++ b/Detectors/Upgrades/ALICE3/TRK/base/include/TRKBase/Specs.h @@ -100,7 +100,7 @@ constexpr int nCols{static_cast(length / chip::pitchZ)}; namespace ML { -constexpr int nLayers{3}; // number of layers in the ML +constexpr int nLayers{5}; // number of layers in the ML constexpr double width{constants::moduleMLOT::width * 1}; // width of the stave // constexpr double length{constants::moduleMLOT::length * 10}; // length of the stave constexpr double length{124 * cm}; // length of the stave, hardcoded to fit the implemented geometry @@ -118,7 +118,7 @@ constexpr double length{258 * cm}; // len constexpr int nRows{static_cast(width / moduleMLOT::chip::pitchX)}; // number of rows in the halfstave constexpr int nCols{static_cast(length / moduleMLOT::chip::pitchZ)}; // number of columns in the halfstave } // namespace halfstave -constexpr int nLayers{5}; // number of layers in the OT +constexpr int nLayers{3}; // number of layers in the OT constexpr double width{halfstave::width * 2}; // width of the stave constexpr double length{halfstave::length}; // length of the stave constexpr int nRows{static_cast(width / moduleMLOT::chip::pitchX)}; // number of rows in the stave diff --git a/Detectors/Upgrades/ALICE3/TRK/base/src/GeometryTGeo.cxx b/Detectors/Upgrades/ALICE3/TRK/base/src/GeometryTGeo.cxx index fe4abda499239..ddfc844cc964d 100644 --- a/Detectors/Upgrades/ALICE3/TRK/base/src/GeometryTGeo.cxx +++ b/Detectors/Upgrades/ALICE3/TRK/base/src/GeometryTGeo.cxx @@ -12,6 +12,7 @@ #include #include #include "TRKBase/SegmentationChip.h" +#include "TRKBase/Specs.h" #include #include @@ -64,7 +65,7 @@ GeometryTGeo::GeometryTGeo(bool build, int loadTrans) : DetMatrixCache(detectors void GeometryTGeo::Build(int loadTrans) { ///// current geometry organization: - ///// total elements = x staves (*2 half staves if staggered geometry) * 8 layers ML+OT + 4 petal cases * (3 layers + 6 disks) + ///// total elements = x staves (*2 half staves if staggered geometry) * ML+OT layers + 4 petal cases * (3 layers + 6 disks) ///// indexing from 0 to 35: VD petals -> layers -> disks ///// indexing from 36 to y: MLOT staves @@ -1130,7 +1131,7 @@ void GeometryTGeo::Print(Option_t*) const std::cout << "Detector ID: " << sInstance.get()->getDetID() << std::endl; LOGF(info, "Summary of GeometryTGeo: %s", getName()); - LOGF(info, "Number of layers ML + OL: %d", mNumberOfLayersMLOT); + LOGF(info, "Number of layers ML + OT: %d", mNumberOfLayersMLOT); LOGF(info, "Number of active parts VD: %d", mNumberOfActivePartsVD); LOGF(info, "Number of layers VD: %d", mNumberOfLayersVD); LOGF(info, "Number of petals VD: %d", mNumberOfPetalsVD); @@ -1142,7 +1143,7 @@ void GeometryTGeo::Print(Option_t*) const LOGF(info, "Number of staves and half staves per layer MLOT: "); for (int i = 0; i < mNumberOfLayersMLOT; i++) { std::string mlot = ""; - mlot = (i < 4) ? "ML" : "OT"; + mlot = (i < constants::ML::nLayers) ? "ML" : "OT"; LOGF(info, "Layer: %d, %s, %d staves, %d half staves per stave", i, mlot.c_str(), mNumberOfStaves[i], mNumberOfHalfStaves[i]); } LOGF(info, "Number of modules per stave (half stave) in each ML(OT) layer: "); diff --git a/Detectors/Upgrades/ALICE3/TRK/simulation/src/Detector.cxx b/Detectors/Upgrades/ALICE3/TRK/simulation/src/Detector.cxx index 66ace4746d399..3fa51afe3ba2b 100644 --- a/Detectors/Upgrades/ALICE3/TRK/simulation/src/Detector.cxx +++ b/Detectors/Upgrades/ALICE3/TRK/simulation/src/Detector.cxx @@ -13,6 +13,7 @@ #include "DetectorsBase/Stack.h" +#include "TRKBase/Specs.h" #include "TRKBase/TRKBaseParam.h" #include "TRKSimulation/Hit.h" #include "TRKSimulation/VDGeometryBuilder.h" @@ -99,7 +100,7 @@ void Detector::configMLOT() case kCylindrical: { const std::vector length{128.35f, 128.35f, 128.35f, 128.35f, 128.35f, 256.7f, 256.7f, 256.7f}; LOGP(warning, "Loading cylindrical configuration for ALICE3 TRK"); - for (int i{0}; i < 8; ++i) { + for (int i{0}; i < constants::ML::nLayers + constants::OT::nLayers; ++i) { std::string name = GeometryTGeo::getTRKLayerPattern() + std::to_string(i); mLayers.push_back(std::make_unique(i, name, rInn[i], length[i], thick, MatBudgetParamMode::Thickness)); } @@ -115,9 +116,9 @@ void Detector::configMLOT() const std::vector stagOffsets{0.f, 0.f, 0.f, 1.17f, 0.89f}; LOGP(warning, "Loading segmented configuration for ALICE3 TRK"); - for (int i{0}; i < 8; ++i) { + for (int i{0}; i < constants::ML::nLayers + constants::OT::nLayers; ++i) { std::string name = GeometryTGeo::getTRKLayerPattern() + std::to_string(i); - if (i < 5) { + if (i < constants::ML::nLayers) { mLayers.push_back(std::make_unique(i, name, rInn[i], stagOffsets[i], tiltAngles[i], nStaves[i], nMods[i], thick, MatBudgetParamMode::Thickness)); } else { mLayers.push_back(std::make_unique(i, name, rInn[i], tiltAngles[i], nStaves[i], nMods[i], thick, MatBudgetParamMode::Thickness)); @@ -202,8 +203,8 @@ void Detector::configFromFile(std::string fileName) // Default mode is Thickness MatBudgetParamMode mode = MatBudgetParamMode::Thickness; - if (layerCount < 5) { - // ML layers (0 to 4) require stagOffset (index 5) + if (layerCount < constants::ML::nLayers) { + // ML layers require stagOffset (index 5) if (tmpBuff.size() < 6) { LOGP(fatal, "Invalid configuration for ML layer {}: stagOffset is missing.", layerCount); } @@ -215,7 +216,7 @@ void Detector::configFromFile(std::string fileName) mLayers.push_back(std::make_unique(layerCount, name, rInn, stagOffset, tiltAngle, nStaves, nMods, thick, mode)); } else { - // OT layers (5+) do NOT have stagOffset. The optional mode is at index 5. + // OT layers do NOT have stagOffset. The optional mode is at index 5. if (tmpBuff.size() >= 6) { mode = static_cast(static_cast(tmpBuff[5])); }