diff --git a/Detectors/Upgrades/ALICE3/IOTOF/base/include/IOTOFBase/GeometryTGeo.h b/Detectors/Upgrades/ALICE3/IOTOF/base/include/IOTOFBase/GeometryTGeo.h index d466cc2987f23..616c6409f4577 100644 --- a/Detectors/Upgrades/ALICE3/IOTOF/base/include/IOTOFBase/GeometryTGeo.h +++ b/Detectors/Upgrades/ALICE3/IOTOF/base/include/IOTOFBase/GeometryTGeo.h @@ -39,7 +39,7 @@ class GeometryTGeo : public o2::detectors::DetMatrixCache static const char* getIOTOFVolPattern() { return sIOTOFVolumeName.c_str(); } // Inner TOF - const int getITOFNumberOfChips() { return mNumberOfChipsIOTOF[0]; } + const int getITOFNumberOfChips() const { return mNumberOfChipsIOTOF[0]; } static const char* getITOFLayerPattern() { return sITOFLayerName.c_str(); } static const char* getITOFStavePattern() { return sITOFStaveName.c_str(); } static const char* getITOFModulePattern() { return sITOFModuleName.c_str(); } @@ -47,7 +47,7 @@ class GeometryTGeo : public o2::detectors::DetMatrixCache static const char* getITOFSensorPattern() { return sITOFSensorName.c_str(); } // Outer TOF - const int getOTOFNumberOfChips() { return mNumberOfChipsIOTOF[1]; } + const int getOTOFNumberOfChips() const { return mNumberOfChipsIOTOF[1]; } static const char* getOTOFLayerPattern() { return sOTOFLayerName.c_str(); } static const char* getOTOFStavePattern() { return sOTOFStaveName.c_str(); } static const char* getOTOFModulePattern() { return sOTOFModuleName.c_str(); } @@ -55,13 +55,13 @@ class GeometryTGeo : public o2::detectors::DetMatrixCache static const char* getOTOFSensorPattern() { return sOTOFSensorName.c_str(); } // Forward TOF - const int getFTOFNumberOfChips() { return mNumberOfChipsFTOF; } + const int getFTOFNumberOfChips() const { return mNumberOfChipsFTOF; } static const char* getFTOFLayerPattern() { return sFTOFLayerName.c_str(); } static const char* getFTOFChipPattern() { return sFTOFChipName.c_str(); } static const char* getFTOFSensorPattern() { return sFTOFSensorName.c_str(); } // Backward TOF - const int getBTOFNumberOfChips() { return mNumberOfChipsBTOF; } + const int getBTOFNumberOfChips() const { return mNumberOfChipsBTOF; } static const char* getBTOFLayerPattern() { return sBTOFLayerName.c_str(); } static const char* getBTOFChipPattern() { return sBTOFChipName.c_str(); } static const char* getBTOFSensorPattern() { return sBTOFSensorName.c_str(); } diff --git a/Detectors/Upgrades/ALICE3/IOTOF/simulation/include/IOTOFSimulation/DPLDigitizerParam.h b/Detectors/Upgrades/ALICE3/IOTOF/simulation/include/IOTOFSimulation/DPLDigitizerParam.h index 66014e1c7f06e..679ace26da68c 100644 --- a/Detectors/Upgrades/ALICE3/IOTOF/simulation/include/IOTOFSimulation/DPLDigitizerParam.h +++ b/Detectors/Upgrades/ALICE3/IOTOF/simulation/include/IOTOFSimulation/DPLDigitizerParam.h @@ -31,6 +31,8 @@ struct DPLDigitizerParam : public o2::conf::ConfigurableParamHelper @@ -30,7 +31,6 @@ namespace o2::iotof { o2::iotof::Segmentation* Digitizer::sSegmentation = nullptr; - //_______________________________________________________________________ void Digitizer::init() { @@ -100,7 +100,7 @@ void Digitizer::processHit(const o2::itsmft::Hit& hit, int evID, int srcID) } // Get detector element ID - int chipID = hit.GetDetectorID(); + const int chipID = hit.GetDetectorID(); auto& chip = mChips[chipID]; if (chip.isDisabled()) { LOG(debug) << "Hit rejected because chip " << chipID << " is disabled"; @@ -110,6 +110,8 @@ void Digitizer::processHit(const o2::itsmft::Hit& hit, int evID, int srcID) // Convert energy loss to charge (number of electrons) float energyLoss = hit.GetEnergyLoss(); // in GeV int charge = energyToCharge(energyLoss); + const auto& digitizerParams = o2::iotof::DPLDigitizerParam::Instance(); + int electronsPerStep = static_cast(charge / digitizerParams.nSimSteps); // Apply charge threshold if (charge < mChargeThreshold) { @@ -128,25 +130,150 @@ void Digitizer::processHit(const o2::itsmft::Hit& hit, int evID, int srcID) LOG(debug) << "Invalid detector ID: " << chipID << ", geometry size: " << mGeometry->getSize(); return; // invalid detector ID } + + // Create the digit with time information + o2::MCCompLabel label(hit.GetTrackID(), evID, srcID, false); + const int roFrameAbs = 0; // For now, we can set this to 0 or calculate based on time if needed + const int nROF = 1; // For now, we can assume the signal is contained in one ROF, this can be extended to multiple ROFs based on the time + + if (digitizerParams.performStepping) { + float** respMatrix = nullptr; + int rowStart = 0, colStart = 0, rowSpan = 0, colSpan = 0; + stepping(hit, respMatrix, rowStart, colStart, rowSpan, colSpan); + + for (int irow = rowSpan; irow--;) { + uint16_t rowIS = irow + rowStart; + for (int icol = colSpan; icol--;) { + uint16_t colIS = icol + colStart; + float nEleResp = respMatrix[irow][icol]; + if (!nEleResp) { + continue; + } + const int nElectronsSampled = gRandom->Poisson(electronsPerStep * nEleResp); + // Noise can be added here if needed + + registerDigits(chip, roFrameAbs, smearedTime, nROF, + static_cast(rowIS), static_cast(colIS), nElectronsSampled, label); + } + } + + for (int irow = 0; irow < rowSpan; ++irow) { + delete[] respMatrix[irow]; + } + delete[] respMatrix; + } else { + int row = 0, col = 0; + responseInTheMiddle(hit, row, col); + const int nElectronsSampled = gRandom->Poisson(charge); + registerDigits(chip, roFrameAbs, smearedTime, nROF, static_cast(row), static_cast(col), nElectronsSampled, label); + } +} + +void Digitizer::stepping(const o2::itsmft::Hit& hit, float**& respMatrix, int& rowStart, int& colStart, int& rowSpan, int& colSpan) +{ + const auto& matrix = mGeometry->getMatrixL2G(hit.GetDetectorID()); + const int chipID = hit.GetDetectorID(); + const int subdetectorID = mGeometry->getIOTOFLayer(chipID); + + auto xyzPositionStart(matrix ^ (hit.GetPosStart())); // start position in sensor frame + auto xyzPositionEnd(matrix ^ (hit.GetPos())); // end position in sensor frame + + const auto& digitizerParams = o2::iotof::DPLDigitizerParam::Instance(); + const auto stepVector = (xyzPositionEnd - xyzPositionStart) / digitizerParams.nSimSteps; + xyzPositionStart = xyzPositionStart + stepVector * 0.5f; // center the start position in the middle of the step + xyzPositionEnd = xyzPositionEnd - stepVector * 0.5f; // center the end position in the middle of the step + + rowStart = -1; + colStart = -1; + int rowEnd = -1, colEnd = -1, nSkip = 0, nSteps = digitizerParams.nSimSteps; + while (!sSegmentation->localToDetector(xyzPositionStart.X(), xyzPositionStart.Z(), rowStart, colStart, mGeometry->getIOTOFLayer(chipID))) { + if (++nSkip > digitizerParams.nSimSteps) { // additional check to add: should we exclude something? + LOG(debug) << "Hit position out of bounds for detector ID " << chipID; + return; // hit is outside the active area + } + xyzPositionStart += stepVector; + } + + while (!sSegmentation->localToDetector(xyzPositionEnd.X(), xyzPositionEnd.Z(), rowEnd, colEnd, mGeometry->getIOTOFLayer(chipID))) { + if (++nSkip > digitizerParams.nSimSteps) { // additional check to add: should we exclude something? + LOG(debug) << "Hit position out of bounds for detector ID " << chipID; + return; // hit is outside the active area + } + xyzPositionEnd += stepVector; + } + + if (rowStart > rowEnd) { + std::swap(rowStart, rowEnd); + } + if (colStart > colEnd) { + std::swap(colStart, colEnd); + } + + // Expand the range to take into account the effects of charge sharing + rowStart -= digitizerParams.responseMatrixSize / 2; + rowEnd += digitizerParams.responseMatrixSize / 2; + rowStart = std::max(rowStart, 0); + colStart = std::max(colStart, 0); + + rowEnd = std::min(rowEnd, (subdetectorID == 0 ? sSegmentation->mITofSpecsConfig.NRows : sSegmentation->mOTofSpecsConfig.NRows) - 1); + colEnd = std::min(colEnd, (subdetectorID == 0 ? sSegmentation->mITofSpecsConfig.NCols : sSegmentation->mOTofSpecsConfig.NCols) - 1); + rowSpan = rowEnd - rowStart + 1; + colSpan = colEnd - colStart + 1; + + respMatrix = new float*[rowSpan]; + for (int i = 0; i < rowSpan; ++i) { + respMatrix[i] = new float[colSpan](); + } + + int rowPrev = -1, colPrev = -1, row = 0, col = 0; + if (!respMatrix || rowSpan <= 0 || colSpan <= 0) { + return; + } + if (nSkip) { + nSteps -= nSkip; + } + + auto& currentPosLocal = xyzPositionStart; + for (int iStep = nSteps; iStep--;) { + sSegmentation->localToDetector(currentPosLocal.X(), currentPosLocal.Z(), row, col, subdetectorID); + if (row != rowPrev || col != colPrev) { + rowPrev = row; + colPrev = col; + } + + currentPosLocal += stepVector; // Move to the next step position + + for (int irow = digitizerParams.responseMatrixSize; irow--;) { + int rowDest = row + irow - (digitizerParams.responseMatrixSize / 2) - rowStart; // destination row in the respMatrix + if (rowDest < 0 || rowDest >= rowSpan) { + continue; + } + for (int icol = digitizerParams.responseMatrixSize; icol--;) { + int colDest = col + icol - (digitizerParams.responseMatrixSize / 2) - colStart; // destination column in the respMatrix + if (colDest < 0 || colDest >= colSpan) { + continue; + } + respMatrix[rowDest][colDest] += 1.; + } + } + } +} + +void Digitizer::responseInTheMiddle(const o2::itsmft::Hit& hit, int& row, int& col) +{ const auto& matrix = mGeometry->getMatrixL2G(hit.GetDetectorID()); + const int chipID = hit.GetDetectorID(); math_utils::Vector3D xyzPositionStart(matrix ^ (hit.GetPosStart())); // start position in sensor frame - // math_utils::Vector3D xyzPositionEnd(matrix ^ (hit.GetPos())); // end position in sensor frame + math_utils::Vector3D xyzPositionEnd(matrix ^ (hit.GetPos())); // end position in sensor frame - int row = 0; // Will be determined from start hit position - int col = 0; // Will be determined from start hit position + // Calculate the middle position of the hit + math_utils::Vector3D xyzPositionMiddle = (xyzPositionStart + xyzPositionEnd) * 0.5f; - if (!sSegmentation->localToDetector(xyzPositionStart.X(), xyzPositionStart.Z(), row, col, mGeometry->getIOTOFLayer(chipID))) { + if (!sSegmentation->localToDetector(xyzPositionMiddle.X(), xyzPositionMiddle.Z(), row, col, mGeometry->getIOTOFLayer(chipID))) { LOG(debug) << "Hit position out of bounds for detector ID " << chipID; return; // hit is outside the active area } - - // Create the digit with time information - o2::MCCompLabel label(hit.GetTrackID(), evID, srcID, false); - const int roFrameAbs = 0; // For now, we can set this to 0 or calculate based on time if needed - const int nROF = 1; // For now, we can assume the signal is contained in one ROF, this can be extended to multiple ROFs based on the time - - registerDigits(chip, roFrameAbs, smearedTime, nROF, static_cast(row), static_cast(col), charge, label); } //_______________________________________________________________________ @@ -200,10 +327,9 @@ void Digitizer::fillOutputContainer() auto& chipDigits = chip.getDigits(); for (const auto& [key, digit] : chipDigits) { - /// Charge threshold not implemented yet - /// if (digit.getCharge() < mChargeThreshold) { - /// continue; // skip digits below threshold - /// } + if (digit.getCharge() < mChargeThreshold) { + continue; // skip digits below threshold + } int digitID = mDigits->size(); mDigits->emplace_back(digit.getChipIndex(), digit.getRow(), digit.getColumn(), digit.getCharge(), digit.getTime());