diff --git a/Detectors/Upgrades/ALICE3/TRK/base/CMakeLists.txt b/Detectors/Upgrades/ALICE3/TRK/base/CMakeLists.txt index a237a2d12211d..96ebf4ead4b7b 100644 --- a/Detectors/Upgrades/ALICE3/TRK/base/CMakeLists.txt +++ b/Detectors/Upgrades/ALICE3/TRK/base/CMakeLists.txt @@ -12,8 +12,10 @@ o2_add_library(TRKBase SOURCES src/GeometryTGeo.cxx src/TRKBaseParam.cxx + src/SegmentationChip.cxx PUBLIC_LINK_LIBRARIES O2::DetectorsBase) o2_target_root_dictionary(TRKBase HEADERS include/TRKBase/GeometryTGeo.h - include/TRKBase/TRKBaseParam.h) \ No newline at end of file + include/TRKBase/TRKBaseParam.h + include/TRKBase/SegmentationChip.h) \ No newline at end of file diff --git a/Detectors/Upgrades/ALICE3/TRK/base/include/TRKBase/GeometryTGeo.h b/Detectors/Upgrades/ALICE3/TRK/base/include/TRKBase/GeometryTGeo.h index cfd991728d09b..a1e4b9321130f 100644 --- a/Detectors/Upgrades/ALICE3/TRK/base/include/TRKBase/GeometryTGeo.h +++ b/Detectors/Upgrades/ALICE3/TRK/base/include/TRKBase/GeometryTGeo.h @@ -73,7 +73,7 @@ class GeometryTGeo : public o2::detectors::DetMatrixCache void setOwner(bool v) { mOwner = v; } void Print(Option_t* opt = "") const; - void PrintChipID(int index, int subDetID, int petalcase, int disk, int lay, int stave, int halfstave, int indexRetrieved) const; + void PrintChipID(int index, int subDetID, int petalcase, int disk, int lay, int stave, int halfstave) const; int getLayer(int index) const; int getStave(int index) const; diff --git a/Detectors/Upgrades/ALICE3/TRK/base/include/TRKBase/SegmentationChip.h b/Detectors/Upgrades/ALICE3/TRK/base/include/TRKBase/SegmentationChip.h new file mode 100644 index 0000000000000..100af5be1b4d0 --- /dev/null +++ b/Detectors/Upgrades/ALICE3/TRK/base/include/TRKBase/SegmentationChip.h @@ -0,0 +1,282 @@ +// 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 SegmentationChip.h +/// \brief Definition of the SegmentationChipclass + +#ifndef ALICEO2_TRK_SEGMENTATIONCHIP_H_ +#define ALICEO2_TRK_SEGMENTATIONCHIP_H_ + +#include +#include + +#include "MathUtils/Cartesian.h" +#include "TRKBase/Specs.h" + +namespace o2::trk +{ + +/// Segmentation and response for TRK chips in ALICE3 upgrade +/// This is a work-in-progress code derived from the ITS2 and ITS3 segmentations. +class SegmentationChip +{ + // This class defines the segmenation of the TRK chips in the ALICE3 upgrade. + // The "global coordinate system" refers to the hit position in cm in the global coordinate system centered in 0,0,0 + // The "local coordinate system" refers to the hit position in cm in the coordinate system of the sensor, which + // is centered in 0,0,0 in the case of curved layers, and in the middle of the chip in the case of flat layers + // The "detector coordinate system" refers to the hit position in row,col inside the sensor + // This class provides the transformations from the local and detector coordinate systems + // The conversion between global and local coordinate systems is operated by the transformation matrices + // For the curved VD layers there exist three coordinate systems and one is transient. + // 1. The global (curved) coordinate system. The chip's center of coordinate system is + // defined at the the mid-point of the detector. + // 2. The local (flat) coordinate system. This is the tube segment projected onto a flat + // surface. In the projection we implicitly assume that the inner and outer + // stretch does not depend on the radius. + // 3. The detector coordinate system. Defined by the row and column segmentation + // defined at the upper edge in the flat coord. + // For the flat ML and OT layers, there exist two coordinate systems: + // 1. The global (flat) coordinate system. The chip's center of coordinate system is + // defined at the the mid-point of the detector. + // 2. The detector coordinate system. Defined by the row and column segmentation + // TODO: add segmentation for VD disks + + public: + constexpr SegmentationChip() = default; + ~SegmentationChip() = default; + constexpr SegmentationChip(const SegmentationChip&) = default; + constexpr SegmentationChip(SegmentationChip&&) = delete; + constexpr SegmentationChip& operator=(const SegmentationChip&) = default; + constexpr SegmentationChip& operator=(SegmentationChip&&) = delete; + + static constexpr float PitchColVD{constants::VD::petal::layer::pitchZ}; + static constexpr float PitchRowVD{constants::VD::petal::layer::pitchX}; + + static constexpr float PitchColMLOT{constants::moduleMLOT::chip::pitchZ}; + static constexpr float PitchRowMLOT{constants::moduleMLOT::chip::pitchX}; + + static constexpr float SensorLayerThicknessVD = {constants::VD::petal::layer::totalThickness}; // physical thickness of sensitive part = 30 um + static constexpr float SensorLayerThicknessML = {constants::moduleMLOT::chip::totalThickness}; // physical thickness of sensitive part = 100 um + static constexpr float SensorLayerThicknessOT = {constants::moduleMLOT::chip::totalThickness}; // physical thickness of sensitive part = 100 um + + static constexpr float SiliconThicknessVD = constants::VD::silicon::thickness; // effective thickness of sensitive part + static constexpr float SiliconThicknessMLOT = constants::moduleMLOT::silicon::thickness; // effective thickness of sensitive part + + static constexpr std::array radiiVD = constants::VD::petal::layer::radii; + + /// Transformation from Geant detector centered local coordinates (cm) to + /// Pixel cell numbers iRow and iCol. + /// Returns kTRUE if point x,z is inside sensitive volume, kFALSE otherwise. + /// A value of -1 for iRow or iCol indicates that this point is outside of the + /// detector segmentation as defined. + /// \param float x Detector local coordinate x in cm with respect to + /// the center of the sensitive volume. + /// \param float z Detector local coordinate z in cm with respect to + /// the center of the sensitive volulme. + /// \param int iRow Detector x cell coordinate. Has the range 0 <= iRow < mNumberOfRows + /// \param int iCol Detector z cell coordinate. Has the range 0 <= iCol < mNumberOfColumns + /// \param int subDetID Sub-detector ID (0 for VD, 1 for ML/OT) + /// \param int layer Layer number (0 to 2 for VD, 0 to 7 for ML/OT) + /// \param int disk Disk number (0 to 5 for VD) + static bool localToDetector(float xRow, float zCol, int& iRow, int& iCol, int subDetID, int layer, int disk) noexcept + { + if (!isValidGlob(xRow, zCol, subDetID, layer)) { + LOGP(debug, "Local coordinates not valid: row = {} cm, col = {} cm", xRow, zCol); + return false; + } + localToDetectorUnchecked(xRow, zCol, iRow, iCol, subDetID, layer, disk); + + LOG(debug) << "Result from localToDetectorUnchecked: xRow " << xRow << " -> iRow " << iRow << ", zCol " << zCol << " -> iCol " << iCol << " on subDetID, layer, disk: " << subDetID << " " << layer << " " << disk; + + if (!isValidDet(iRow, iCol, subDetID, layer)) { + iRow = iCol = -1; + LOGP(debug, "Detector coordinates not valid: iRow = {}, iCol = {}", iRow, iCol); + return false; + } + return true; + }; + /// same but w/o check for row/column range + static void localToDetectorUnchecked(float xRow, float zCol, int& iRow, int& iCol, int subDetID, int layer, int disk) noexcept + { + // convert to row/col w/o over/underflow check + float pitchRow(0), pitchCol(0); + float maxWidth(0), maxLength(0); + + if (subDetID == 0) { + pitchRow = PitchRowVD; + pitchCol = PitchColVD; + maxWidth = constants::VD::petal::layer::width[layer]; + maxLength = constants::VD::petal::layer::length; + // TODO: change this to use the layer and disk + } else if (subDetID == 1 && layer <= 3) { // ML + pitchRow = PitchRowMLOT; + pitchCol = PitchColMLOT; + maxWidth = constants::ML::width; + maxLength = constants::ML::length; + } else if (subDetID == 1 && layer >= 4) { // OT + pitchRow = PitchRowMLOT; + pitchCol = PitchColMLOT; + maxWidth = constants::OT::width; + maxLength = constants::OT::length; + } + // convert to row/col + iRow = static_cast(std::floor((maxWidth / 2 - xRow) / pitchRow)); + iCol = static_cast(std::floor((zCol + maxLength / 2) / pitchCol)); + }; + + // Check local coordinates (cm) validity. + static constexpr bool isValidGlob(float x, float z, int subDetID, int layer) noexcept + { + float maxWidth(0), maxLength(0); + if (subDetID == 0) { + maxWidth = constants::VD::petal::layer::width[layer]; + maxLength = constants::VD::petal::layer::length; + // TODO: change this to use the layer and disk + } else if (subDetID == 1 && layer <= 3) { // ML + maxWidth = constants::ML::width; + maxLength = constants::ML::length; + } else if (subDetID == 1 && layer >= 4) { // OT + maxWidth = constants::OT::width; + maxLength = constants::OT::length; + } + return (-maxWidth / 2 < x && x < maxWidth / 2 && -maxLength / 2 < z && z < maxLength / 2); + } + + // Check detector coordinates validity. + static constexpr bool isValidDet(float row, float col, int subDetID, int layer) noexcept + { + // Check if the row and column are within the valid range + int nRows(0), nCols(0); + if (subDetID == 0) { + nRows = constants::VD::petal::layer::nRows[layer]; + nCols = constants::VD::petal::layer::nCols; + // TODO: change this to use the layer and disk + } else if (subDetID == 1 && layer <= 3) { // ML + nRows = constants::ML::nRows; + nCols = constants::ML::nCols; + } else if (subDetID == 1 && layer >= 4) { // OT + nRows = constants::OT::nRows; + nCols = constants::OT::nCols; + } + return (row >= 0 && row < static_cast(nRows) && col >= 0 && col < static_cast(nCols)); + } + + /// Transformation from Detector cell coordinates to Geant detector centered + /// local coordinates (cm) + /// \param int iRow Detector x cell coordinate. + /// \param int iCol Detector z cell coordinate. + /// \param float x Detector local coordinate x in cm with respect to the + /// center of the sensitive volume. + /// \param float z Detector local coordinate z in cm with respect to the + /// center of the sensitive volume. + /// If iRow and or iCol is outside of the segmentation range a value of -0.5*Dx() + /// or -0.5*Dz() is returned. + /// \param int subDetID Sub-detector ID (0 for VD, 1 for ML/OT) + /// \param int layer Layer number (0 to 2 for VD, 0 to 7 for ML/OT) + /// \param int disk Disk number (0 to 5 for VD) + static constexpr bool detectorToLocal(int iRow, int iCol, float& xRow, float& zCol, int subDetID, int layer, int disk) noexcept + { + if (!isValidDet(iRow, iCol, subDetID, layer)) { + LOGP(debug, "Detector coordinates not valid: iRow = {}, iCol = {}", iRow, iCol); + return false; + } + detectorToLocalUnchecked(iRow, iCol, xRow, zCol, subDetID, layer, disk); + LOG(debug) << "Result from detectorToLocalUnchecked: iRow " << iRow << " -> xRow " << xRow << ", iCol " << iCol << " -> zCol " << zCol << " on subDetID, layer, disk: " << subDetID << " " << layer << " " << disk; + + if (!isValidGlob(xRow, zCol, subDetID, layer)) { + LOGP(debug, "Local coordinates not valid: row = {} cm, col = {} cm", xRow, zCol); + return false; + } + return true; + }; + + // Same as detectorToLocal w.o. checks. + // We position ourself in the middle of the pixel. + static void detectorToLocalUnchecked(int row, int col, float& xRow, float& zCol, int subDetID, int layer, int disk) noexcept + { + /// xRow = half chip width - iRow(center) * pitch + /// zCol = iCol * pitch - half chip lenght + if (subDetID == 0) { + xRow = 0.5 * (constants::VD::petal::layer::width[layer] - PitchRowVD) - (row * PitchRowVD); + zCol = col * PitchColVD + 0.5 * (PitchColVD - constants::VD::petal::layer::length); + } else if (subDetID == 1 && layer <= 3) { // ML + xRow = 0.5 * (constants::ML::width - PitchRowMLOT) - (row * PitchRowMLOT); + zCol = col * PitchRowMLOT + 0.5 * (PitchRowMLOT - constants::ML::length); + } else if (subDetID == 1 && layer >= 4) { // OT + xRow = 0.5 * (constants::OT::width - PitchRowMLOT) - (row * PitchRowMLOT); + zCol = col * PitchColMLOT + 0.5 * (PitchColMLOT - constants::OT::length); + } + } + + /// Transformation from the curved surface to a flat surface. + /// Additionally a shift in the flat coordinates must be applied because + /// the center of the TGeoShap when projected will be higher than the + /// physical thickness of the chip. Thus we shift the projected center + /// down by this difference to align the coordinate systems. + /// \param layer VD layer number + /// \param xCurved Detector local curved coordinate x in cm with respect to + /// the center of the sensitive volume. + /// \param yCurved Detector local curved coordinate y in cm with respect to + /// the center of the sensitive volume. + /// \return math_utils::Vector2D: x and y represent the detector local flat coordinates x and y + // in cm with respect to the center of the sensitive volume. + static math_utils::Vector2D curvedToFlat(const int layer, const float xCurved, const float yCurved) noexcept + { + // Align the flat surface with the curved survace of the original chip (and account for metal stack, TODO) + float dist = std::hypot(xCurved, yCurved); + float phi = std::atan2(yCurved, xCurved); + + // the y position is in the silicon volume however we need the chip volume (silicon+metalstack) + // this is accounted by a y shift + float xFlat = constants::VD::petal::layer::radii[layer] * phi; /// this is equal to the circumference segment covered between y=0 and the phi angle + float yFlat = constants::VD::petal::layer::radii[layer] - dist; + return math_utils::Vector2D(xFlat, yFlat); + } + + /// Transformation from the flat surface to a curved surface + /// It works only if the detector is not rototraslated. + /// \param layer VD layer number + /// \param xFlat Detector local flat coordinate x in cm with respect to + /// the center of the sensitive volume. + /// \param yFlat Detector local flat coordinate y in cm with respect to + /// the center of the sensitive volume. + /// \return math_utils::Vector2D: x and y represent the detector local curved coordinates x and y + // in cm with respect to the center of the sensitive volume. + static constexpr math_utils::Vector2D flatToCurved(int layer, float xFlat, float yFlat) noexcept + { + // Revert the curvedToFlat transformation + float dist = constants::VD::petal::layer::radii[layer] - yFlat; + float phi = xFlat / constants::VD::petal::layer::radii[layer]; + // the y position is in the chip volume however we need the silicon volume + // this is accounted by a -y shift + float xCurved = dist * std::cos(phi); + float yCurved = dist * std::sin(phi); + return math_utils::Vector2D(xCurved, yCurved); + } + + /// Print segmentation info + static const void Print() noexcept + { + LOG(info) << "Number of rows:\nVD L0: " << constants::VD::petal::layer::nRows[0] + << "\nVD L1: " << constants::VD::petal::layer::nRows[1] + << "\nVD L2: " << constants::VD::petal::layer::nRows[2] + << "\nML stave: " << constants::ML::nRows + << "\nOT stave: " << constants::OT::nRows; + + LOG(info) << "Number of cols:\nVD: " << constants::VD::petal::layer::nCols + << "\nML stave: " << constants::ML::nCols + << "\nOT stave: " << constants::OT::nCols; + } +}; + +} // namespace o2::trk + +#endif diff --git a/Detectors/Upgrades/ALICE3/TRK/base/include/TRKBase/Specs.h b/Detectors/Upgrades/ALICE3/TRK/base/include/TRKBase/Specs.h new file mode 100644 index 0000000000000..373e9d972656b --- /dev/null +++ b/Detectors/Upgrades/ALICE3/TRK/base/include/TRKBase/Specs.h @@ -0,0 +1,129 @@ +// 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 Specs.h +/// \brief specs of the ALICE3 TRK + +#ifndef O2_ALICE_TRK_SPECS +#define O2_ALICE_TRK_SPECS + +#include +#include +// This is a temporary version with the specs for the ALICE3 TRK +// This files defines the design specifications of the chips for VD, ML, OT. +// Each TGeoShape has the following properties +// length: dimension in z-axis +// width: dimension in xy-axes +// color: for visulisation +namespace o2::trk::constants +{ +// Default unit of TGeo = cm +constexpr double cm{1}; +constexpr double mu{1e-4}; +constexpr double mm{1e-1}; + +namespace VD // TODO: add a primitive segmentation with more granularity wrt 1/4 layer = 1 chip +{ +namespace silicon +{ +constexpr double thickness{30 * mu}; // thickness of the silicon (should be 10 um epitaxial layer + 20 um substrate)? +} // namespace silicon +namespace metalstack +{ +constexpr double thickness{0 * mu}; // thickness of the copper metal stack - for the moment it is not implemented +} // namespace metalstack +namespace petal +{ +constexpr int nLayers{3}; // number of layers in each VD petal +constexpr int nDisks{6}; // number of disks in each VD petal +namespace layer +{ +constexpr double pitchX{10 * mu}; // pitch of the row +constexpr double pitchZ{10 * mu}; // pitch of the column +constexpr double totalThickness{silicon::thickness + metalstack::thickness}; // total thickness of the chip +constexpr std::array radii{0.5 * cm, 1.2 * cm, 2.5 * cm}; // radius of layer in cm +constexpr std::array width{radii[0] * 2 * M_PI / 4, radii[1] * 2 * M_PI / 4, radii[2] * 2 * M_PI / 4}; // width of the quarter of layer in cm +constexpr double length{50 * cm}; // length of the layer +constexpr int nCols{static_cast(length / pitchZ)}; // number of columns in the chip +constexpr std::array nRows{static_cast(width[0] / pitchX), static_cast(width[1] / pitchX), static_cast(width[2] / pitchX)}; // number of rows in the chip. For the moment is different for each layer since a siner segmentation in repetitive units is stil to be implemented + +} // namespace layer +namespace disk +{ //// TODO: to be filled +constexpr double radiusIn{0.5 * cm}; +constexpr double radiusOut{2.5 * cm}; +} // namespace disk +} // namespace petal +} // namespace VD + +namespace moduleMLOT /// same for ML and OT for the moment +{ /// TODO: account for different modules in case of changes +namespace silicon +{ +constexpr double thickness{100 * mu}; // thickness of the silicon (should be 10 um epitaxial layer + 90 um substrate)? +} // namespace silicon +namespace metalstack +{ +constexpr double thickness{0 * mu}; // thickness of the copper metal stack - for the moment it is not implemented +} // namespace metalstack +namespace chip +{ +constexpr double width{25 * mm}; // width of the chip +constexpr double length{32 * mm}; // length of the chip +constexpr double pitchX{50 * mu}; // pitch of the row +constexpr double pitchZ{50 * mu}; // pitch of the column +constexpr int nRows{static_cast(width / pitchX)}; // number of columns in the chip +constexpr int nCols{static_cast(length / pitchZ)}; // number of rows in the chip +constexpr double totalThickness{silicon::thickness + metalstack::thickness}; // total thickness of the chip +/// Set to 0 for the moment, to be adjusted with the actual design of the chip if needed +static constexpr float PassiveEdgeReadOut = 0.f; // width of the readout edge (Passive bottom) +static constexpr float PassiveEdgeTop = 0.f; // Passive area on top +static constexpr float PassiveEdgeSide = 0.f; // width of Passive area on left/right of the sensor +} // namespace chip +namespace gaps +{ +constexpr double interChips{0.2 * mm}; // gap between the chips +constexpr double outerEdgeLongSide{1 * mm}; // gap between the chips and the outer edges (long side) +constexpr double outerEdgeShortSide{0.1 * mm}; // gap between the chips and the outer edges (short side) +} // namespace gaps +constexpr double width{chip::width * 2 + gaps::interChips + 2 * gaps::outerEdgeLongSide}; // width of the module +constexpr double length{chip::length * 4 + 3 * gaps::interChips + 2 * gaps::outerEdgeShortSide}; // length of the module +constexpr int nRows{static_cast(width / chip::pitchX)}; // number of columns in the module +constexpr int nCols{static_cast(length / chip::pitchZ)}; // number of rows in the module +} // namespace moduleMLOT + +namespace ML +{ +constexpr double width{constants::moduleMLOT::width * 1}; // width of the stave +constexpr double length{constants::moduleMLOT::length * 10}; // length of the stave +constexpr int nRows{static_cast(width / constants::moduleMLOT::chip::pitchX)}; // number of rows in the stave +constexpr int nCols{static_cast(length / constants::moduleMLOT::chip::pitchZ)}; // number of columns in the stave +} // namespace ML + +namespace OT +{ //// TODO: add shorter lenght of the stave of L4 +constexpr double width{moduleMLOT::width * 2}; // width of the stave +constexpr double length{moduleMLOT::length * 20}; // length of the stave +constexpr int nRows{static_cast(width / moduleMLOT::chip::pitchX)}; // number of rows in the stave +constexpr int nCols{static_cast(length / moduleMLOT::chip::pitchZ)}; // number of columns in the stave +} // namespace OT + +namespace apts /// parameters for the APTS response +{ +constexpr double pitchX{15.0 * mu}; +constexpr double pitchZ{15.0 * mu}; +constexpr double responseYShift{15.5 * mu}; +constexpr double thickness{45 * mu}; +} // namespace apts + +} // namespace o2::trk::constants + +#endif diff --git a/Detectors/Upgrades/ALICE3/TRK/base/src/GeometryTGeo.cxx b/Detectors/Upgrades/ALICE3/TRK/base/src/GeometryTGeo.cxx index 69bae0fad9cee..20088179f4dcc 100644 --- a/Detectors/Upgrades/ALICE3/TRK/base/src/GeometryTGeo.cxx +++ b/Detectors/Upgrades/ALICE3/TRK/base/src/GeometryTGeo.cxx @@ -11,9 +11,9 @@ #include #include -// #include "TRKBase/SegmentationChip.h" +#include "TRKBase/SegmentationChip.h" -// using Segmentation = o2::trk::SegmentationChip; +using Segmentation = o2::trk::SegmentationChip; namespace o2 { @@ -263,58 +263,35 @@ bool GeometryTGeo::getChipID(int index, int& subDetID, int& petalcase, int& disk TString GeometryTGeo::getMatrixPath(int index) const { - int subDetID, petalcase, disk, lay, stave, halfstave; //// TODO: add chips in a second step - getChipID(index, subDetID, petalcase, disk, lay, stave, halfstave); + int subDetID, petalcase, disk, layer, stave, halfstave; //// TODO: add chips in a second step + getChipID(index, subDetID, petalcase, disk, layer, stave, halfstave); - int indexRetrieved = getChipIndex(subDetID, petalcase, disk, lay, stave, halfstave); + // PrintChipID(index, subDetID, petalcase, disk, layer, stave, halfstave); - PrintChipID(index, subDetID, petalcase, disk, lay, stave, halfstave, indexRetrieved); + // TString path = "/cave_1/barrel_1/TRKV_2/TRKLayer0_1/TRKStave0_1/TRKChip0_1/TRKSensor0_1/"; /// dummy path, to be used for tests + TString path = Form("/cave_1/barrel_1/%s_2/", GeometryTGeo::getTRKVolPattern()); - // TString path = Form("/cave_1/barrel_1/%s_2/", GeometryTGeo::getTRKVolPattern()); - TString path = "/cave_1/barrel_1/TRKV_2/TRKLayer0_1/TRKStave0_1/TRKChip0_1/TRKSensor0_1/"; /// dummy path, to be replaced - - // if (wrID >= 0) { - // path += Form("%s%d_1/", getITSWrapVolPattern(), wrID); - // } - - // if (isVD) { - // path += Form("%s%d_1/", getTRKPetalPattern(), index); - - // } else { - // path += Form("%s%d_1/", getTRKLayerPattern(), index); - // } - - // if (!mIsLayerITS3[lay]) { - // path += - // Form("%s%d_1/", getITSLayerPattern(), lay); - // if (mNumberOfHalfBarrels > 0) { - // path += Form("%s%d_%d/", getITSHalfBarrelPattern(), lay, hba); - // } - // path += - // Form("%s%d_%d/", getITSStavePattern(), lay, stav); - - // if (mNumberOfHalfStaves[lay] > 0) { - // path += Form("%s%d_%d/", getITSHalfStavePattern(), lay, sstav); - // } - // if (mNumberOfModules[lay] > 0) { - // path += Form("%s%d_%d/", getITSModulePattern(), lay, mod); - // } - // path += Form("%s%d_%d/%s%d_1", getITSChipPattern(), lay, chipInMod, getITSSensorPattern(), lay); - // } else { - // // hba = carbonform - // // stav = 0 - // // sstav = segment - // // mod = rsu - // // chipInMod = tile - // // sensor = pixelarray - // path += Form("%s_0/", getITS3LayerPattern(lay)); - // path += Form("%s_%d/", getITS3CarbonFormPattern(lay), hba); - // path += Form("%s_0/", getITS3ChipPattern(lay)); - // path += Form("%s_%d/", getITS3SegmentPattern(lay), sstav); - // path += Form("%s_%d/", getITS3RSUPattern(lay), mod); - // path += Form("%s_%d/", getITS3TilePattern(lay), chipInMod); - // path += Form("%s_0", getITS3PixelArrayPattern(lay)); - // } + if (subDetID == 0) { // VD + if (disk >= 0) { + path += Form("%s%d_%s%d_1/", getTRKPetalPattern(), petalcase, getTRKPetalDiskPattern(), disk); // PETALCASEx_DISKy_1 + path += Form("%s%d_%s%d_%s%d_1/", getTRKPetalPattern(), petalcase, getTRKPetalDiskPattern(), disk, getTRKChipPattern(), disk); // PETALCASEx_DISKy_TRKChipy_1 + path += Form("%s%d_%s%d_%s%d_1/", getTRKPetalPattern(), petalcase, getTRKPetalDiskPattern(), disk, getTRKSensorPattern(), disk); // PETALCASEx_DISKy_TRKSensory_1 + } else if (layer >= 0) { + path += Form("%s%d_%s%d_1/", getTRKPetalPattern(), petalcase, getTRKPetalLayerPattern(), layer); // PETALCASEx_LAYERy_1 + path += Form("%s%d_%s%d_%s%d_1/", getTRKPetalPattern(), petalcase, getTRKPetalLayerPattern(), layer, getTRKStavePattern(), layer); // PETALCASEx_LAYERy_TRKStavey_1 + path += Form("%s%d_%s%d_%s%d_1/", getTRKPetalPattern(), petalcase, getTRKPetalLayerPattern(), layer, getTRKChipPattern(), layer); // PETALCASEx_LAYERy_TRKChipy_1 + path += Form("%s%d_%s%d_%s%d_1/", getTRKPetalPattern(), petalcase, getTRKPetalLayerPattern(), layer, getTRKSensorPattern(), layer); // PETALCASEx_LAYERy_TRKSensory_1 + } + } else if (subDetID == 1) { // MLOT + path += Form("%s%d_1/", getTRKLayerPattern(), layer); // TRKLayerx_1 + path += Form("%s%d_%d/", getTRKStavePattern(), layer, stave); // TRKStavex_y + if (mNumberOfHalfStaves[layer] == 2) { // staggered geometry + path += Form("%s%d_%d/", getTRKChipPattern(), layer, halfstave); // TRKChipx_0/1 + } else if (mNumberOfHalfStaves[layer] == 1) { // turbo geometry + path += Form("%s%d_1/", getTRKChipPattern(), layer); // TRKChipx_1 + } + path += Form("%s%d_1/", getTRKSensorPattern(), layer); // TRKSensorx_1 + } return path; } @@ -325,40 +302,40 @@ TGeoHMatrix* GeometryTGeo::extractMatrixSensor(int index) const // Note, the if the effective sensitive layer thickness is smaller than the // total physical sensor tickness, this matrix is biased and connot be used // directly for transformation from sensor frame to global one. - // // Therefore we need to add a shift + auto path = getMatrixPath(index); static TGeoHMatrix matTmp; - // gGeoManager->PushPath(); // Preserve the modeler state. + gGeoManager->PushPath(); // Preserve the modeler state. - // if (!gGeoManager->cd(path.Data())) { - // gGeoManager->PopPath(); - // LOG(error) << "Error in cd-ing to " << path.Data(); - // return nullptr; - // } // end if !gGeoManager + if (!gGeoManager->cd(path.Data())) { + gGeoManager->PopPath(); + LOG(error) << "Error in cd-ing to " << path.Data(); + return nullptr; + } // end if !gGeoManager matTmp = *gGeoManager->GetCurrentMatrix(); // matrix may change after cd // RSS - // printf("%d/%d/%d %s\n", lay, stav, detInSta, path.Data()); // matTmp.Print(); // Restore the modeler state. gGeoManager->PopPath(); static int chipInGlo{0}; + /// TODO: // account for the difference between physical sensitive layer (where charge collection is simulated) and effective sensor thicknesses - // in the ITS3 case this accounted by specialized functions - // double delta = Segmentation::SensorLayerThickness; - // static TGeoTranslation tra(0., 0.5 * delta, 0.); - // #ifdef ENABLE_UPGRADES // only apply for non ITS3 OB layers - // if (!mIsLayerITS3[getLayer(index)]) { - // matTmp *= tra; - // } - // #else + // in the VD case this will be accounted by specialized functions during the clusterization (following what it is done for ITS3) + // this can be done once the right sensor thickness is in place in the geometry + // double delta = 0.; + // if (getSubDetID(index) == 1){ /// ML/OT + // delta = Segmentation::SensorLayerThicknessVD - Segmentation::SiliconTickness; + // static TGeoTranslation tra(0., 0.5 * delta, 0.); // matTmp *= tra; - // #endif + // } + // std::cout<<"-----"< + +using namespace o2::trk; + +// void SegmentationChip::print() +// { +// printf("++++++++++ VD ++++++++++"); +// printf("Pixel size: %.2f (along %d rows) %.2f (along %d columns) microns\n", PitchRowVD * 1e4, -999, PitchColVD * 1e4, NColsVD); +// printf("++++++++++ ML ++++++++++"); +// printf("Pixel size: %.2f (along %d rows) %.2f (along %d columns) microns\n", PitchRowML * 1e4, -999, PitchColML * 1e4, NColsML); + +// } \ No newline at end of file diff --git a/Detectors/Upgrades/ALICE3/TRK/simulation/CMakeLists.txt b/Detectors/Upgrades/ALICE3/TRK/simulation/CMakeLists.txt index a1cb0279efef8..ab817a3fdaa0d 100644 --- a/Detectors/Upgrades/ALICE3/TRK/simulation/CMakeLists.txt +++ b/Detectors/Upgrades/ALICE3/TRK/simulation/CMakeLists.txt @@ -11,7 +11,10 @@ o2_add_library(TRKSimulation SOURCES src/TRKLayer.cxx + src/ChipDigitsContainer.cxx + src/ChipSimResponse.cxx src/Detector.cxx + src/DigiParams.cxx src/Digitizer.cxx src/TRKServices.cxx src/DPLDigitizerParam.cxx @@ -21,10 +24,14 @@ o2_add_library(TRKSimulation PUBLIC_LINK_LIBRARIES O2::TRKBase O2::FT3Simulation O2::ITSMFTSimulation + O2::DetectorsRaw O2::SimulationDataFormat) o2_target_root_dictionary(TRKSimulation - HEADERS include/TRKSimulation/Digitizer.h + HEADERS include/TRKSimulation/ChipDigitsContainer.h + include/TRKSimulation/ChipSimResponse.h + include/TRKSimulation/DigiParams.h + include/TRKSimulation/Digitizer.h include/TRKSimulation/Detector.h include/TRKSimulation/TRKLayer.h include/TRKSimulation/TRKServices.h diff --git a/Detectors/Upgrades/ALICE3/TRK/simulation/include/TRKSimulation/ChipDigitsContainer.h b/Detectors/Upgrades/ALICE3/TRK/simulation/include/TRKSimulation/ChipDigitsContainer.h new file mode 100644 index 0000000000000..658fb823bb596 --- /dev/null +++ b/Detectors/Upgrades/ALICE3/TRK/simulation/include/TRKSimulation/ChipDigitsContainer.h @@ -0,0 +1,37 @@ +// 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. + +#ifndef ALICEO2_TRK_CHIPDIGITSCONTAINER_ +#define ALICEO2_TRK_CHIPDIGITSCONTAINER_ + +#include "ITSMFTBase/SegmentationAlpide.h" +#include "ITSMFTSimulation/ChipDigitsContainer.h" +#include "TRKBase/SegmentationChip.h" +#include "TRKBase/Specs.h" +#include "TRKSimulation/DigiParams.h" +#include + +namespace o2::trk +{ + +class ChipDigitsContainer : public o2::itsmft::ChipDigitsContainer +{ + public: + explicit ChipDigitsContainer(UShort_t idx = 0); + + using Segmentation = SegmentationChip; + + ClassDefNV(ChipDigitsContainer, 1); +}; + +} // namespace o2::trk + +#endif // ALICEO2_TRK_CHIPDIGITSCONTAINER_ diff --git a/Detectors/Upgrades/ALICE3/TRK/simulation/include/TRKSimulation/ChipSimResponse.h b/Detectors/Upgrades/ALICE3/TRK/simulation/include/TRKSimulation/ChipSimResponse.h new file mode 100644 index 0000000000000..29147997f66bf --- /dev/null +++ b/Detectors/Upgrades/ALICE3/TRK/simulation/include/TRKSimulation/ChipSimResponse.h @@ -0,0 +1,37 @@ +// 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. + +#ifndef ALICEO2_TRKSIMULATION_CHIPSIMRESPONSE_H +#define ALICEO2_TRKSIMULATION_CHIPSIMRESPONSE_H + +#include "ITSMFTSimulation/AlpideSimResponse.h" + +namespace o2 +{ +namespace trk +{ + +class ChipSimResponse : public o2::itsmft::AlpideSimResponse +{ + public: + ChipSimResponse() = default; + ChipSimResponse(const ChipSimResponse& other) = default; + ChipSimResponse(const o2::itsmft::AlpideSimResponse* base) : o2::itsmft::AlpideSimResponse(*base) {} + + void initData(int tableNumber, std::string dataPath, const bool quiet = true); + + ClassDef(ChipSimResponse, 1); +}; + +} // namespace trk +} // namespace o2 + +#endif // ALICEO2_TRKSIMULATION_CHIPSIMRESPONSE_H diff --git a/Detectors/Upgrades/ALICE3/TRK/simulation/include/TRKSimulation/DPLDigitizerParam.h b/Detectors/Upgrades/ALICE3/TRK/simulation/include/TRKSimulation/DPLDigitizerParam.h index 59b3551ecbd32..584ffaa3aff75 100644 --- a/Detectors/Upgrades/ALICE3/TRK/simulation/include/TRKSimulation/DPLDigitizerParam.h +++ b/Detectors/Upgrades/ALICE3/TRK/simulation/include/TRKSimulation/DPLDigitizerParam.h @@ -39,7 +39,7 @@ struct DPLDigitizerParam : public o2::conf::ConfigurableParamHelper +#include +#include "TRKBase/TRKBaseParam.h" +#include "TRKBase/GeometryTGeo.h" + +//////////////////////////////////////////////////////////// +// // +// Simulation params for the TRK digitizer // +// // +// This is a provisionary implementation, until proper // +// microscopic simulation and its configuration will // +// be implemented // +// // +//////////////////////////////////////////////////////////// + +namespace o2 +{ +namespace trk +{ + +class ChipSimResponse; + +class DigiParams +{ + + using SignalShape = o2::itsmft::AlpideSignalTrapezoid; + + public: + DigiParams(); + ~DigiParams() = default; + + 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; } + + void setROFrameLength(float ns); + float getROFrameLength() const { return mROFrameLength; } + float getROFrameLengthInv() const { return mROFrameLengthInv; } + + void setStrobeDelay(float ns) { mStrobeDelay = ns; } + float getStrobeDelay() const { return mStrobeDelay; } + + void setStrobeLength(float ns) { mStrobeLength = ns; } + float getStrobeLength() const { return mStrobeLength; } + + void setTimeOffset(double sec) { mTimeOffset = sec; } + double getTimeOffset() const { return mTimeOffset; } + + void setROFrameBiasInBC(int n) { mROFrameBiasInBC = n; } + int getROFrameBiasInBC() const { return mROFrameBiasInBC; } + + void setChargeThreshold(int v, float frac2Account = 0.1); + void setNSimSteps(int v); + void setEnergyToNElectrons(float v) { mEnergyToNElectrons = v; } + + void setVbb(float v) { mVbb = v; } + void setIBVbb(float v) { mIBVbb = v; } + void setOBVbb(float v) { mOBVbb = v; } + + int getChargeThreshold() const { return mChargeThreshold; } + int getMinChargeToAccount() const { return mMinChargeToAccount; } + int getNSimSteps() const { return mNSimSteps; } + float getNSimStepsInv() const { return mNSimStepsInv; } + float getEnergyToNElectrons() const { return mEnergyToNElectrons; } + + float getVbb() const { return mVbb; } + float getIBVbb() const { return mIBVbb; } + float getOBVbb() const { return mOBVbb; } + + bool isTimeOffsetSet() const { return mTimeOffset > -infTime; } + + const o2::trk::ChipSimResponse* getAlpSimResponse() const { return mAlpSimResponse; } + void setAlpSimResponse(const o2::trk::ChipSimResponse* par) { mAlpSimResponse = par; } + + const SignalShape& getSignalShape() const { return mSignalShape; } + SignalShape& getSignalShape() { return (SignalShape&)mSignalShape; } + + virtual void print() const; + + private: + static constexpr double infTime = 1e99; + bool mIsContinuous = false; ///< flag for continuous simulation + float mNoisePerPixel = 1.e-8; ///< ALPIDE 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 = 150; ///< charge threshold in Nelectrons + int mMinChargeToAccount = 15; ///< minimum charge contribution to account + int mNSimSteps = 18; ///< number of steps in response simulation + float mNSimStepsInv = 0; ///< its inverse + + float mEnergyToNElectrons = 1. / 3.6e-9; // conversion of eloss to Nelectrons + + float mVbb = 0.0; ///< back bias absolute value for MFT (in Volt) + 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) + + o2::itsmft::AlpideSignalTrapezoid mSignalShape; ///< signal timeshape parameterization + + const o2::trk::ChipSimResponse* mAlpSimResponse = nullptr; //!< pointer on external response + + // auxiliary precalculated parameters + float mROFrameLengthInv = 0; ///< inverse length of RO frame in ns + + // ClassDef(DigiParams, 2); +}; +} // namespace trk +} // namespace o2 + +#endif \ No newline at end of file diff --git a/Detectors/Upgrades/ALICE3/TRK/simulation/include/TRKSimulation/Digitizer.h b/Detectors/Upgrades/ALICE3/TRK/simulation/include/TRKSimulation/Digitizer.h index 6863c5392cae3..573217fe9b076 100644 --- a/Detectors/Upgrades/ALICE3/TRK/simulation/include/TRKSimulation/Digitizer.h +++ b/Detectors/Upgrades/ALICE3/TRK/simulation/include/TRKSimulation/Digitizer.h @@ -21,12 +21,12 @@ #include "Rtypes.h" // for Digitizer::Class #include "TObject.h" // for TObject -#include "ITSMFTSimulation/ChipDigitsContainer.h" -// #include "ITSMFTSimulation/AlpideSimResponse.h" -#include "ITSMFTSimulation/DigiParams.h" +#include "TRKSimulation/ChipSimResponse.h" +#include "TRKSimulation/ChipDigitsContainer.h" + +#include "TRKSimulation/DigiParams.h" #include "ITSMFTSimulation/Hit.h" #include "TRKBase/GeometryTGeo.h" -// #include "ITS3Base/SegmentationSuperAlpide.h" #include "DataFormatsITSMFT/Digit.h" #include "DataFormatsITSMFT/ROFRecord.h" #include "CommonDataFormat/InteractionRecord.h" @@ -37,7 +37,7 @@ namespace o2::trk { -class Digitizer : public TObject +class Digitizer { using ExtraDig = std::vector; ///< container for extra contributions to PreDigits @@ -46,11 +46,13 @@ class Digitizer : public TObject void setMCLabels(o2::dataformats::MCTruthContainer* mclb) { mMCLabels = mclb; } void setROFRecords(std::vector* rec) { mROFRecords = rec; } - o2::itsmft::DigiParams& getParams() { return (o2::itsmft::DigiParams&)mParams; } - const o2::itsmft::DigiParams& getParams() const { return mParams; } + o2::trk::DigiParams& getParams() { return (o2::trk::DigiParams&)mParams; } + const o2::trk::DigiParams& getParams() const { return mParams; } void init(); + 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); @@ -64,10 +66,10 @@ class Digitizer : public TObject bool isContinuous() const { return mParams.isContinuous(); } void fillOutputContainer(uint32_t maxFrame = 0xffffffff); - void setDigiParams(const o2::itsmft::DigiParams& par) { mParams = par; } - const o2::itsmft::DigiParams& getDigitParams() const { return mParams; } + void setDigiParams(const o2::trk::DigiParams& par) { mParams = par; } + const o2::trk::DigiParams& getDigitParams() const { return mParams; } - // provide the common itsmft::GeometryTGeo to access matrices and segmentation + // provide the common trk::GeometryTGeo to access matrices and segmentation void setGeometry(const o2::trk::GeometryTGeo* gm) { mGeometry = gm; } uint32_t getEventROFrameMin() const { return mEventROFrameMin; } @@ -82,7 +84,7 @@ class Digitizer : public TObject private: void processHit(const o2::itsmft::Hit& hit, uint32_t& maxFr, int evID, int srcID); - void registerDigits(o2::itsmft::ChipDigitsContainer& chip, uint32_t roFrame, float tInROF, int nROF, + void registerDigits(o2::trk::ChipDigitsContainer& chip, uint32_t roFrame, float tInROF, int nROF, uint16_t row, uint16_t col, int nEle, o2::MCCompLabel& lbl); ExtraDig* getExtraDigBuffer(uint32_t roFrame) @@ -97,9 +99,41 @@ class Digitizer : public TObject return mExtraBuff[ind].get(); } + /// Get the number of columns according to the subdetector + /// \param subDetID 0 for VD, 1 for ML/OT + /// \param layer 0 to 2 for VD, 0 to 7 for ML/OT + /// \return Number of columns (for the moment, in the entire layer(VD) or stave (ML/OT) + int getNCols(int subDetID, int layer) + { + if (subDetID == 0) { // VD + return constants::VD::petal::layer::nCols; + } else if (subDetID == 1 && layer <= 3) { // ML + return constants::ML::nCols; + } else if (subDetID == 1 && layer >= 4) { // OT + return constants::OT::nCols; + } + return 0; + } + + /// Get the number of rows according to the subdetector + /// \param subDetID 0 for VD, 1 for ML/OT + /// \param layer 0 to 2 for VD, 0 to 7 for ML/OT + /// \return Number of rows (for the moment, in the entire layer(VD) or stave (ML/OT) + int getNRows(int subDetID, int layer) + { + if (subDetID == 0) { // VD + return constants::VD::petal::layer::nRows[layer]; + } else if (subDetID == 1 && layer <= 3) { // ML + return constants::ML::nRows; + } else if (subDetID == 1 && layer >= 4) { // OT + return constants::OT::nRows; + } + return 0; + } + static constexpr float sec2ns = 1e9; - o2::itsmft::DigiParams mParams; ///< digitization parameters + o2::trk::DigiParams mParams; ///< digitization parameters o2::InteractionTimeRecord mEventTime; ///< global event time and interaction record o2::InteractionRecord mIRFirstSampledTF; ///< IR of the 1st sampled IR, noise-only ROFs will be inserted till this IR only double mCollisionTimeWrtROF{}; @@ -110,19 +144,35 @@ class Digitizer : public TObject uint32_t mEventROFrameMin = 0xffffffff; ///< lowest RO frame for processed events (w/o automatic noise ROFs) uint32_t mEventROFrameMax = 0; ///< highest RO frame forfor processed events (w/o automatic noise ROFs) - o2::itsmft::AlpideSimResponse* mAlpSimResp = nullptr; // simulated response + int mNumberOfChips = 0; + + o2::trk::ChipSimResponse* mChipSimResp = nullptr; // simulated response + o2::trk::ChipSimResponse* mChipSimRespVD = nullptr; // simulated response for VD chips + o2::trk::ChipSimResponse* mChipSimRespMLOT = nullptr; // simulated response for ML/OT chips + + // std::string mResponseFile = "$(O2_ROOT)/share/Detectors/ITSMFT/data/AlpideResponseData/AlpideResponseData.root"; + std::string mResponseFile = "$(O2_ROOT)/share/Detectors/Upgrades/ITS3/data/ITS3ChipResponseData/APTSResponseData.root"; /// using temporarly the APTS response + + bool mSimRespOrientation{false}; // wether the orientation in the response function is flipped + float mSimRespVDShift{0.f}; // adjusting the Y-shift in the APTS response function to match sensor local coord. + float mSimRespVDScaleX{1.f}; // scale x-local coordinate to response function x-coordinate + float mSimRespVDScaleZ{1.f}; // scale z-local coordinate to response function z-coordinate + float mSimRespMLOTShift{0.f}; // adjusting the Y-shift in the APTS response function to match sensor local coord. + float mSimRespMLOTScaleX{1.f}; // scale x-local coordinate to response function x-coordinate + float mSimRespMLOTScaleZ{1.f}; // scale z-local coordinate to response function z-coordinate + float mSimRespVDScaleDepth{1.f}; // scale depth-local coordinate to response function depth-coordinate + float mSimRespMLOTScaleDepth{1.f}; // scale depth-local coordinate to response function depth-coordinate const o2::trk::GeometryTGeo* mGeometry = nullptr; ///< TRK geometry - std::vector mChips; ///< Array of chips digits containers - std::deque> mExtraBuff; ///< burrer (per roFrame) for extra digits + std::vector mChips; ///< Array of chips digits containers + std::deque> mExtraBuff; ///< buffer (per roFrame) for extra digits std::vector* mDigits = nullptr; //! output digits std::vector* mROFRecords = nullptr; //! output ROF records o2::dataformats::MCTruthContainer* mMCLabels = nullptr; //! output labels const o2::itsmft::NoiseMap* mDeadChanMap = nullptr; - - ClassDef(Digitizer, 1); + const o2::itsmft::NoiseMap* mNoiseMap = nullptr; }; -} // namespace o2::trk \ No newline at end of file +} // namespace o2::trk diff --git a/Detectors/Upgrades/ALICE3/TRK/simulation/src/ChipDigitsContainer.cxx b/Detectors/Upgrades/ALICE3/TRK/simulation/src/ChipDigitsContainer.cxx new file mode 100644 index 0000000000000..9ed4a4bedf5c5 --- /dev/null +++ b/Detectors/Upgrades/ALICE3/TRK/simulation/src/ChipDigitsContainer.cxx @@ -0,0 +1,17 @@ +// 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. + +#include "TRKSimulation/ChipDigitsContainer.h" + +using namespace o2::trk; + +ChipDigitsContainer::ChipDigitsContainer(UShort_t idx) + : o2::itsmft::ChipDigitsContainer(idx) {} diff --git a/Detectors/Upgrades/ALICE3/TRK/simulation/src/ChipSimResponse.cxx b/Detectors/Upgrades/ALICE3/TRK/simulation/src/ChipSimResponse.cxx new file mode 100644 index 0000000000000..70c4f131b9724 --- /dev/null +++ b/Detectors/Upgrades/ALICE3/TRK/simulation/src/ChipSimResponse.cxx @@ -0,0 +1,21 @@ +// 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. + +#include "TRKSimulation/ChipSimResponse.h" +#include +#include + +using namespace o2::trk; + +void ChipSimResponse::initData(int tableNumber, std::string dataPath, const bool quiet) +{ + AlpideSimResponse::initData(tableNumber, dataPath, quiet); +} \ No newline at end of file diff --git a/Detectors/Upgrades/ALICE3/TRK/simulation/src/DigiParams.cxx b/Detectors/Upgrades/ALICE3/TRK/simulation/src/DigiParams.cxx new file mode 100644 index 0000000000000..df6f46ac0ecb0 --- /dev/null +++ b/Detectors/Upgrades/ALICE3/TRK/simulation/src/DigiParams.cxx @@ -0,0 +1,72 @@ +// 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 DigiParams.cxx +/// \brief Implementation of the TRK digitization steering params. Based on the ITS2 code. + +#include // for LOG +#include "TRKSimulation/DigiParams.h" +#include + +using namespace o2::trk; + +DigiParams::DigiParams() +{ + // make sure the defaults are consistent + setNSimSteps(mNSimSteps); +} + +void DigiParams::setROFrameLength(float lNS) +{ + // set ROFrame length in nanosecongs + mROFrameLength = lNS; + assert(mROFrameLength > 1.); + mROFrameLengthInv = 1. / mROFrameLength; +} + +void DigiParams::setNSimSteps(int v) +{ + // set number of sampling steps in silicon + mNSimSteps = v > 0 ? v : 1; + mNSimStepsInv = 1.f / mNSimSteps; +} + +void DigiParams::setChargeThreshold(int v, float frac2Account) +{ + // set charge threshold for digits creation and its fraction to account + // contribution from single hit + mChargeThreshold = v; + mMinChargeToAccount = v * frac2Account; + if (mMinChargeToAccount < 0 || mMinChargeToAccount > mChargeThreshold) { + mMinChargeToAccount = mChargeThreshold; + } + LOG(info) << "Set charge threshold to " << mChargeThreshold + << ", single hit will be accounted from " << mMinChargeToAccount + << " electrons"; +} + +//______________________________________________ +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); + printf("ELoss to N electrons factor : %e\n", mEnergyToNElectrons); + printf("Noise level per pixel : %e\n", mNoisePerPixel); + printf("Charge time-response:\n"); + mSignalShape.print(); +} diff --git a/Detectors/Upgrades/ALICE3/TRK/simulation/src/Digitizer.cxx b/Detectors/Upgrades/ALICE3/TRK/simulation/src/Digitizer.cxx index 21e6e629ec418..cc89f0eff1a54 100644 --- a/Detectors/Upgrades/ALICE3/TRK/simulation/src/Digitizer.cxx +++ b/Detectors/Upgrades/ALICE3/TRK/simulation/src/Digitizer.cxx @@ -12,456 +12,563 @@ /// \file Digitizer.cxx #include "DataFormatsITSMFT/Digit.h" -// #include "ITSMFTBase/SegmentationAlpide.h" +#include "TRKBase/SegmentationChip.h" #include "TRKSimulation/DPLDigitizerParam.h" +#include "TRKSimulation/TRKLayer.h" #include "TRKSimulation/Digitizer.h" -// #include "MathUtils/Cartesian.h" -// #include "SimulationDataFormat/MCTruthContainer.h" -// #include "DetectorsRaw/HBFUtils.h" +#include "DetectorsRaw/HBFUtils.h" -// #include +#include // #include -// #include -// #include +#include +#include +#include #include // for LOG using o2::itsmft::Digit; using o2::itsmft::Hit; -// using Segmentation = o2::itsmft::SegmentationAlpide; +using Segmentation = o2::trk::SegmentationChip; using namespace o2::trk; +using namespace o2::itsmft; // using namespace o2::base; - //_______________________________________________________________________ void Digitizer::init() { - // mNumberOfChips = mGeometry->getNumberOfChips(); - // mChips.resize(mNumberOfChips); - // for (int i = mNumberOfChips; i--;) { - // mChips[i].setChipIndex(i); - // if (mNoiseMap) { - // mChips[i].setNoiseMap(mNoiseMap); - // } - // if (mDeadChanMap) { - // mChips[i].disable(mDeadChanMap->isFullChipMasked(i)); - // mChips[i].setDeadChanMap(mDeadChanMap); - // } - // } - // initializing for both collection tables - /*for (int i = 0; i < 2; i++) { - mAlpSimResp[i].initData(i); - }*/ - + LOG(info) << "Initializing digitizer"; + mNumberOfChips = mGeometry->getNumberOfChips(); + mChips.resize(mNumberOfChips); /// temporary, to not make it crash + for (int i = mNumberOfChips; i--;) { + mChips[i].setChipIndex(i); + if (mNoiseMap) { + mChips[i].setNoiseMap(mNoiseMap); + } + if (mDeadChanMap) { + mChips[i].disable(mDeadChanMap->isFullChipMasked(i)); + mChips[i].setDeadChanMap(mDeadChanMap); + } + } // importing the charge collection tables // (initialized while building O2) - // auto file = TFile::Open(mResponseFile.data()); - // if (!file) { - // LOG(fatal) << "Cannot open response file " << mResponseFile; - // } - /*std::string response = "response"; - for (int i=0; i<2; i++) { - response.append(std::to_string(i)); - mAlpSimResp[i] = *(o2::itsmft::AlpideSimResponse*)file->Get(response.data()); - }*/ - // mAlpSimResp[0] = *(o2::itsmft::AlpideSimResponse*)file->Get("response0"); - // mAlpSimResp[1] = *(o2::itsmft::AlpideSimResponse*)file->Get("response1"); + auto file = TFile::Open(mResponseFile.data()); + if (!file) { + LOG(fatal) << "Cannot open response file " << mResponseFile; + } + + // setting the correct response function (for the moment, for both VD and MLOT the APTS response function is udes) + mChipSimResp = (o2::trk::ChipSimResponse*)file->Get("response1"); + mChipSimRespVD = mChipSimResp; /// for the moment considering the same response + mChipSimRespMLOT = mChipSimResp; /// for the moment considering the same response + + /// setting scale factors to adapt to the APTS response function (adjusting pitch and Y shift) + // TODO: adjust Y shift when the geometry is improved + LOG(debug) << " Depth max: " << mChipSimRespVD->getDepthMax(); + LOG(debug) << " Depth min: " << mChipSimRespVD->getDepthMin(); + + float thicknessVD = 0.0095; // cm --- hardcoded based on geometry currently present + float thicknessMLOT = 0.1; // cm --- hardcoded based on geometry currently present + + mSimRespVDScaleX = o2::trk::constants::apts::pitchX / o2::trk::SegmentationChip::PitchRowVD; + mSimRespVDScaleZ = o2::trk::constants::apts::pitchZ / o2::trk::SegmentationChip::PitchColVD; + mSimRespVDScaleDepth = o2::trk::constants::apts::thickness / (thicknessVD); /// introducing this scaling factor because the silicon thickness for the moment is 1 mm -> rescale to 45 um which is the depth of the APTS response + // mSimRespVDShift = mChipSimRespVD->getDepthMax() - thicknessVD * mSimRespVDScaleDepth / 2.f; // the shift should be done considering the rescaling done to adapt to the wrong silicon thickness. TODO: remove the scaling factor for the depth when the silicon thickness match the simulated response + mSimRespVDShift = mChipSimRespVD->getDepthMax(); // the curved, rescaled, sensors have a width from 0 to -45. Must add 10 um (= max depth) to match the APTS response. + mSimRespMLOTScaleX = o2::trk::constants::apts::pitchX / o2::trk::SegmentationChip::PitchRowMLOT; + mSimRespMLOTScaleZ = o2::trk::constants::apts::pitchZ / o2::trk::SegmentationChip::PitchColMLOT; + mSimRespMLOTScaleDepth = o2::trk::constants::apts::thickness / (thicknessMLOT); /// introducing this scaling factor because the silicon thickness for the moment is 1 mm -> rescale to 45 um which is the depth of the APTS response + mSimRespMLOTShift = mChipSimRespMLOT->getDepthMax() - thicknessMLOT * mSimRespMLOTScaleDepth / 2.f; // the shift should be done considering the rescaling done to adapt to the wrong silicon thickness. TODO: remove the scaling factor for the depth when the silicon thickness match the simulated response + mSimRespOrientation = false; // importing the parameters from DPLDigitizerParam.h auto& dOptTRK = DPLDigitizerParam::Instance(); - LOGP(info, "TRK Digitizer is initalised."); + LOGP(info, "TRK Digitizer is initialised."); + mParams.print(); + LOGP(info, "VD shift = {} ; ML/OT shift = {} = {} - {}", mSimRespVDShift, mSimRespMLOTShift, mChipSimRespMLOT->getDepthMax(), thicknessMLOT * mSimRespMLOTScaleDepth / 2.f); + LOGP(info, "VD pixel scale on x = {} ; z = {}", mSimRespVDScaleX, mSimRespVDScaleZ); + LOGP(info, "ML/OT pixel scale on x = {} ; z = {}", mSimRespMLOTScaleX, mSimRespMLOTScaleZ); + LOGP(info, "Response orientation: {}", mSimRespOrientation ? "flipped" : "normal"); + + mIRFirstSampledTF = o2::raw::HBFUtils::Instance().getFirstSampledTFIR(); } -// auto Digitizer::getChipResponse(int chipID) -// { -// if (mNumberOfChips < 10000) { // in MFT -// return mAlpSimRespMFT; -// } +o2::trk::ChipSimResponse* Digitizer::getChipResponse(int chipID) +{ + if (mGeometry->getSubDetID(chipID) == 0) { /// VD + return mChipSimRespVD; + } -// if (chipID < 432) { // in ITS Inner Barrel -// return mAlpSimRespIB; -// } else { // in ITS Outter Barrel -// return mAlpSimRespOB; -// } -// } + else if (mGeometry->getSubDetID(chipID) == 1) { /// ML/OT + return mChipSimRespMLOT; + } + return nullptr; +}; //_______________________________________________________________________ void Digitizer::process(const std::vector* hits, int evID, int srcID) { // digitize single event, the time must have been set beforehand - // LOG(info) << "Digitizing " << mGeometry->getName() << " hits of entry " << evID << " from source " - // << srcID << " at time " << mEventTime << " ROFrame= " << mNewROFrame << ")" - // << " cont.mode: " << isContinuous() - // << " Min/Max ROFrames " << mROFrameMin << "/" << mROFrameMax; + LOG(info) << " Digitizing " << mGeometry->getName() << " (ID: " << mGeometry->getDetID() + << ") hits of entry " << evID << " from source " << srcID + << " at time " << mEventTime << " ROFrame= " << mNewROFrame << ")" + << " cont.mode: " << isContinuous() + << " Min/Max ROFrames " << mROFrameMin << "/" << mROFrameMax; + + std::cout << "Printing segmentation info: " << std::endl; + SegmentationChip::Print(); // // is there something to flush ? - // if (mNewROFrame > mROFrameMin) { - // fillOutputContainer(mNewROFrame - 1); // flush out all frame preceding the new one - // } - - // int nHits = hits->size(); - // std::vector hitIdx(nHits); - // std::iota(std::begin(hitIdx), std::end(hitIdx), 0); - // // sort hits to improve memory access - // std::sort(hitIdx.begin(), hitIdx.end(), - // [hits](auto lhs, auto rhs) { - // return (*hits)[lhs].GetDetectorID() < (*hits)[rhs].GetDetectorID(); - // }); - // 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); - // } + if (mNewROFrame > mROFrameMin) { + fillOutputContainer(mNewROFrame - 1); // flush out all frames preceding the new one + } + + int nHits = hits->size(); + std::vector hitIdx(nHits); + std::iota(std::begin(hitIdx), std::end(hitIdx), 0); + // sort hits to improve memory access + std::sort(hitIdx.begin(), hitIdx.end(), + [hits](auto lhs, auto rhs) { + 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); + } } //_______________________________________________________________________ void Digitizer::setEventTime(const o2::InteractionTimeRecord& irt) { - // // 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--; - // } - - // // we might get interactions to digitize from before - // // the first sampled IR - // if (nbc < 0) { - // mNewROFrame = 0; - // // this event is before the first RO - // mIsBeforeFirstRO = true; - // } else { - // mNewROFrame = nbc / mParams.getROFrameLengthInBC(); - // mIsBeforeFirstRO = false; - // } - // LOG(info) << " NewROFrame " << mNewROFrame << " nbc " << nbc; - - // // in continuous mode depends on starts of periodic readout frame - // mCollisionTimeWrtROF += (nbc % mParams.getROFrameLengthInBC()) * o2::constants::lhc::LHCBunchSpacingNS; - // } else { - // mNewROFrame = 0; - // } - - // 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) { - // mROFrameMax = mNewROFrame - 1; // all frames up to this are finished - // } + LOG(info) << "Setting event time "; + // 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--; + } + + mNewROFrame = nbc / mParams.getROFrameLengthInBC(); + + LOG(info) << " NewROFrame " << mNewROFrame << " = " << nbc << "/" << mParams.getROFrameLengthInBC() << " (nbc/mParams.getROFrameLengthInBC()"; + + // in continuous mode depends on starts of periodic readout frame + mCollisionTimeWrtROF += (nbc % mParams.getROFrameLengthInBC()) * o2::constants::lhc::LHCBunchSpacingNS; + } else { + mNewROFrame = 0; + } + + 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) { + mROFrameMax = mNewROFrame - 1; // all frames up to this are finished + } } //_______________________________________________________________________ void Digitizer::fillOutputContainer(uint32_t frameLast) { // // fill output with digits from min.cached up to requested frame, generating the noise beforehand - // if (frameLast > mROFrameMax) { - // frameLast = mROFrameMax; - // } + if (frameLast > mROFrameMax) { + frameLast = mROFrameMax; + } // // make sure all buffers for extra digits are created up to the maxFrame - // getExtraDigBuffer(mROFrameMax); - - // LOG(info) << "Filling " << mGeometry->getName() << " digits output for RO frames " << mROFrameMin << ":" - // << frameLast; - - // o2::itsmft::ROFRecord rcROF; - - // // we have to write chips in RO increasing order, therefore have to loop over the frames here - // for (; mROFrameMin <= frameLast; mROFrameMin++) { - // rcROF.setROFrame(mROFrameMin); - // rcROF.setFirstEntry(mDigits->size()); // start of current ROF in digits - - // auto& extra = *(mExtraBuff.front().get()); - // for (auto& chip : mChips) { - // if (chip.isDisabled()) { - // continue; - // } - // chip.addNoise(mROFrameMin, mROFrameMin, &mParams); - // auto& buffer = chip.getPreDigits(); - // if (buffer.empty()) { - // continue; - // } - // auto itBeg = buffer.begin(); - // auto iter = itBeg; - // ULong64_t maxKey = chip.getOrderingKey(mROFrameMin + 1, 0, 0) - 1; // fetch digits with key below that - // for (; iter != buffer.end(); ++iter) { - // if (iter->first > maxKey) { - // break; // is the digit ROFrame from the key > the max requested frame - // } - // auto& preDig = iter->second; // preDigit - // if (preDig.charge >= mParams.getChargeThreshold()) { - // int digID = mDigits->size(); - // mDigits->emplace_back(chip.getChipIndex(), preDig.row, preDig.col, preDig.charge); - // mMCLabels->addElement(digID, preDig.labelRef.label); - // auto& nextRef = preDig.labelRef; // extra contributors are in extra array - // while (nextRef.next >= 0) { - // nextRef = extra[nextRef.next]; - // mMCLabels->addElement(digID, nextRef.label); - // } - // } - // } - // buffer.erase(itBeg, iter); - // } - // // 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? - // } - // if (mROFRecords) { - // mROFRecords->push_back(rcROF); - // } - // extra.clear(); // clear container for extra digits of the mROFrameMin ROFrame - // // and move it as a new slot in the end - // mExtraBuff.emplace_back(mExtraBuff.front().release()); - // mExtraBuff.pop_front(); - // } + getExtraDigBuffer(mROFrameMax); + LOG(info) << "Filling " << mGeometry->getName() << " digits output for RO frames " << mROFrameMin << ":" + << frameLast; + + o2::itsmft::ROFRecord rcROF; /// using temporarly itsmft::ROFRecord + + // we have to write chips in RO increasing order, therefore have to loop over the frames here + for (; mROFrameMin <= frameLast; mROFrameMin++) { + rcROF.setROFrame(mROFrameMin); + rcROF.setFirstEntry(mDigits->size()); // start of current ROF in digits + + auto& extra = *(mExtraBuff.front().get()); + for (auto& chip : mChips) { + if (chip.isDisabled()) { + continue; + } + // chip.addNoise(mROFrameMin, mROFrameMin, &mParams); /// TODO: add noise + auto& buffer = chip.getPreDigits(); + if (buffer.empty()) { + continue; + } + auto itBeg = buffer.begin(); + auto iter = itBeg; + ULong64_t maxKey = chip.getOrderingKey(mROFrameMin + 1, 0, 0) - 1; // fetch digits with key below that + for (; iter != buffer.end(); ++iter) { + if (iter->first > maxKey) { + break; // is the digit ROFrame from the key > the max requested frame + } + auto& preDig = iter->second; // preDigit + if (preDig.charge >= mParams.getChargeThreshold()) { + int digID = mDigits->size(); + mDigits->emplace_back(chip.getChipIndex(), preDig.row, preDig.col, preDig.charge); + LOG(debug) << "Adding digit ID: " << digID << " with chipID: " << chip.getChipIndex() << ", row: " << preDig.row << ", col: " << preDig.col << ", charge: " << preDig.charge; + mMCLabels->addElement(digID, preDig.labelRef.label); + auto& nextRef = preDig.labelRef; // extra contributors are in extra array + while (nextRef.next >= 0) { + nextRef = extra[nextRef.next]; + mMCLabels->addElement(digID, nextRef.label); + } + } + } + buffer.erase(itBeg, iter); + } + // 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? + } + if (mROFRecords) { + mROFRecords->push_back(rcROF); + } + extra.clear(); // clear container for extra digits of the mROFrameMin ROFrame + // and move it as a new slot in the end + mExtraBuff.emplace_back(mExtraBuff.front().release()); + mExtraBuff.pop_front(); + } } //_______________________________________________________________________ void Digitizer::processHit(const o2::itsmft::Hit& hit, uint32_t& maxFr, int evID, int srcID) { - // // convert single hit to digits - // int chipID = hit.GetDetectorID(); - // auto& chip = mChips[chipID]; - // if (chip.isDisabled()) { - // LOG(debug) << "skip disabled chip " << chipID; - // return; - // } - // float timeInROF = hit.GetTime() * sec2ns; - // if (timeInROF > 20e3) { - // const int maxWarn = 10; - // static int warnNo = 0; - // if (warnNo < maxWarn) { - // LOG(warning) << "Ignoring hit with time_in_event = " << timeInROF << " ns" - // << ((++warnNo < maxWarn) ? "" : " (suppressing further warnings)"); - // } - // return; - // } - // if (isContinuous()) { - // timeInROF += mCollisionTimeWrtROF; - // } - // if (mIsBeforeFirstRO && timeInROF < 0) { - // // disregard this hit because it comes from an event before readout starts and it does not effect this RO - // return; - // } - - // // calculate RO Frame for this hit - // if (timeInROF < 0) { - // timeInROF = 0.; - // } - // float tTot = mParams.getSignalShape().getMaxDuration(); - // // frame of the hit signal start wrt event ROFrame - // int roFrameRel = int(timeInROF * mParams.getROFrameLengthInv()); - // // 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; - // int nFrames = roFrameRelMax + 1 - roFrameRel; - // uint32_t roFrameMax = mNewROFrame + roFrameRelMax; - // if (roFrameMax > maxFr) { - // maxFr = roFrameMax; // if signal extends beyond current maxFrame, increase the latter - // } - - // // here we start stepping in the depth of the sensor to generate charge diffusion - // float nStepsInv = mParams.getNSimStepsInv(); - // int nSteps = mParams.getNSimSteps(); - // const auto& matrix = mGeometry->getMatrixL2G(hit.GetDetectorID()); - // math_utils::Vector3D xyzLocS(matrix ^ (hit.GetPosStart())); // start position in sensor frame - // math_utils::Vector3D xyzLocE(matrix ^ (hit.GetPos())); // end position in sensor frame - - // math_utils::Vector3D step(xyzLocE); - // step -= xyzLocS; - // step *= nStepsInv; // position increment at each step - // // the electrons will injected in the middle of each step - // math_utils::Vector3D stepH(step * 0.5); - // xyzLocS += stepH; - // xyzLocE -= stepH; - - // int rowS = -1, colS = -1, rowE = -1, colE = -1, nSkip = 0; - // // get entrance pixel row and col - // while (!Segmentation::localToDetector(xyzLocS.X(), xyzLocS.Z(), rowS, colS)) { // guard-ring ? - // if (++nSkip >= nSteps) { - // return; // did not enter to sensitive matrix - // } - // xyzLocS += step; - // } - // // get exit pixel row and col - // while (!Segmentation::localToDetector(xyzLocE.X(), xyzLocE.Z(), rowE, colE)) { // guard-ring ? - // if (++nSkip >= nSteps) { - // return; // did not enter to sensitive matrix - // } - // xyzLocE -= step; - // } - // // estimate the limiting min/max row and col where the non-0 response is possible - // if (rowS > rowE) { - // std::swap(rowS, rowE); - // } - // if (colS > colE) { - // std::swap(colS, colE); - // } - // rowS -= AlpideRespSimMat::NPix / 2; - // rowE += AlpideRespSimMat::NPix / 2; - // if (rowS < 0) { - // rowS = 0; - // } - // if (rowE >= Segmentation::NRows) { - // rowE = Segmentation::NRows - 1; - // } - // colS -= AlpideRespSimMat::NPix / 2; - // colE += AlpideRespSimMat::NPix / 2; - // if (colS < 0) { - // colS = 0; - // } - // if (colE >= Segmentation::NCols) { - // colE = Segmentation::NCols - 1; - // } - // int rowSpan = rowE - rowS + 1, colSpan = colE - colS + 1; // size of plaquet where some response is expected - - // float respMatrix[rowSpan][colSpan]; // response accumulated here - // std::fill(&respMatrix[0][0], &respMatrix[0][0] + rowSpan * colSpan, 0.f); - - // float nElectrons = hit.GetEnergyLoss() * mParams.getEnergyToNElectrons(); // total number of deposited electrons - // nElectrons *= nStepsInv; // N electrons injected per step - // if (nSkip) { - // nSteps -= nSkip; - // } - // // - // int rowPrev = -1, colPrev = -1, row, col; - // float cRowPix = 0.f, cColPix = 0.f; // local coordinated of the current pixel center - - // const o2::itsmft::AlpideSimResponse* resp = getChipResponse(chipID); - - // // take into account that the AlpideSimResponse depth defintion has different min/max boundaries - // // although the max should coincide with the surface of the epitaxial layer, which in the chip - // // local coordinates has Y = +SensorLayerThickness/2 - - // xyzLocS.SetY(xyzLocS.Y() + resp->getDepthMax() - Segmentation::SensorLayerThickness / 2.); - - // // collect charge in every pixel which might be affected by the hit - // for (int iStep = nSteps; iStep--;) { - // // Get the pixel ID - // Segmentation::localToDetector(xyzLocS.X(), xyzLocS.Z(), row, col); - // if (row != rowPrev || col != colPrev) { // update pixel and coordinates of its center - // if (!Segmentation::detectorToLocal(row, col, cRowPix, cColPix)) { - // continue; // should not happen - // } - // rowPrev = row; - // colPrev = col; - // } - // bool flipCol, flipRow; - // // note that response needs coordinates along column row (locX) (locZ) then depth (locY) - // auto rspmat = resp->getResponse(xyzLocS.X() - cRowPix, xyzLocS.Z() - cColPix, xyzLocS.Y(), flipRow, flipCol); - - // xyzLocS += step; - // if (!rspmat) { - // continue; - // } - - // for (int irow = AlpideRespSimMat::NPix; irow--;) { - // int rowDest = row + irow - AlpideRespSimMat::NPix / 2 - rowS; // destination row in the respMatrix - // if (rowDest < 0 || rowDest >= rowSpan) { - // continue; - // } - // for (int icol = AlpideRespSimMat::NPix; icol--;) { - // int colDest = col + icol - AlpideRespSimMat::NPix / 2 - colS; // destination column in the respMatrix - // if (colDest < 0 || colDest >= colSpan) { - // continue; - // } - // respMatrix[rowDest][colDest] += rspmat->getValue(irow, icol, flipRow, flipCol); - // } - // } - // } - - // // fire the pixels assuming Poisson(n_response_electrons) - // o2::MCCompLabel lbl(hit.GetTrackID(), evID, srcID, false); - // auto roFrameAbs = mNewROFrame + roFrameRel; - // for (int irow = rowSpan; irow--;) { - // uint16_t rowIS = irow + rowS; - // for (int icol = colSpan; icol--;) { - // float nEleResp = respMatrix[irow][icol]; - // if (!nEleResp) { - // continue; - // } - // int nEle = gRandom->Poisson(nElectrons * nEleResp); // total charge in given pixel - // // ignore charge which have no chance to fire the pixel - // if (nEle < mParams.getMinChargeToAccount()) { - // continue; - // } - // uint16_t colIS = icol + colS; - // if (mNoiseMap && mNoiseMap->isNoisy(chipID, rowIS, colIS)) { - // continue; - // } - // if (mDeadChanMap && mDeadChanMap->isNoisy(chipID, rowIS, colIS)) { - // continue; - // } - // // - // registerDigits(chip, roFrameAbs, timeInROF, nFrames, rowIS, colIS, nEle, lbl); - // } - // } + 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); + + int layer = mGeometry->getLayer(chipID); + int disk = mGeometry->getDisk(chipID); + + if (disk != -1) { + LOG(debug) << "Skipping disk " << disk; + return; // skipping hits on disks for the moment + } + + LOG(debug) << "Processing hit for chip " << chipID; + auto& chip = mChips[chipID]; + if (chip.isDisabled()) { + LOG(debug) << "Skipping disabled chip " << chipID; + return; + } + float timeInROF = hit.GetTime() * sec2ns; + LOG(debug) << "timeInROF: " << timeInROF; + if (timeInROF > 20e3) { + const int maxWarn = 10; + static int warnNo = 0; + if (warnNo < maxWarn) { + LOG(warning) << "Ignoring hit with time_in_event = " << timeInROF << " ns" + << ((++warnNo < maxWarn) ? "" : " (suppressing further warnings)"); + } + return; + } + if (isContinuous()) { + timeInROF += mCollisionTimeWrtROF; + } + if (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; + return; + } + + // calculate RO Frame for this hit + if (timeInROF < 0) { + timeInROF = 0.; + } + float tTot = mParams.getSignalShape().getMaxDuration(); + // frame of the hit signal start wrt event ROFrame + int roFrameRel = int(timeInROF * mParams.getROFrameLengthInv()); + // 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; + int nFrames = roFrameRelMax + 1 - roFrameRel; + uint32_t roFrameMax = mNewROFrame + roFrameRelMax; + if (roFrameMax > maxFr) { + maxFr = roFrameMax; // if signal extends beyond current maxFrame, increase the latter + } + + // here we start stepping in the depth of the sensor to generate charge diffusion + float nStepsInv = mParams.getNSimStepsInv(); + int nSteps = mParams.getNSimSteps(); + + const auto& matrix = mGeometry->getMatrixL2G(hit.GetDetectorID()); + // matrix.print(); + + /// transorm from the global detector coordinates to the local detector coordinates + math_utils::Vector3D xyzLocS(matrix ^ (hit.GetPosStart())); // start position in sensor frame + math_utils::Vector3D xyzLocE(matrix ^ (hit.GetPos())); // end position in sensor frame + + if (subDetID == 0) { // VD - need to take into account for the curved layers. TODO: consider the disks + // transform the point on the curved surface to a flat one + math_utils::Vector2D xyFlatS = Segmentation::curvedToFlat(layer, xyzLocS.x(), xyzLocS.y()); + math_utils::Vector2D xyFlatE = Segmentation::curvedToFlat(layer, xyzLocE.x(), xyzLocE.y()); + LOG(debug) << "Called curved to flat: " << xyzLocS.x() << " -> " << xyFlatS.x() << ", " << xyzLocS.y() << " -> " << xyFlatS.y(); + // update the local coordinates with the flattened ones + xyzLocS.SetXYZ(xyFlatS.x(), xyFlatS.y(), xyzLocS.Z()); + xyzLocE.SetXYZ(xyFlatE.x(), xyFlatE.y(), xyzLocE.Z()); + } + + // std::cout<<"Printing example of point in 0.35 0.35 0 in global frame: "< examplehitGlob(0.35, 0.35, 0); + // math_utils::Vector3D exampleLoc(matrix ^ (examplehitGlob)); // start position in sensor frame + // std::cout<< "Example hit in local frame: " << exampleLoc << std::endl; + // std::cout<<"Going back to glob coordinates: " << (matrix * exampleLoc) << std::endl; + + //// adapting the depth (Y) of the chip to the APTS response maximum depth + LOG(debug) << "local original: startPos = " << xyzLocS << ", endPos = " << xyzLocE << std::endl; + if (subDetID == 0) { + xyzLocS.SetY(xyzLocS.Y() * mSimRespVDScaleDepth); + xyzLocE.SetY(xyzLocE.Y() * mSimRespVDScaleDepth); + } else { + xyzLocS.SetY(xyzLocS.Y() * mSimRespMLOTScaleDepth); + xyzLocE.SetY(xyzLocE.Y() * mSimRespMLOTScaleDepth); + } + LOG(debug) << "rescaled Y: startPos = " << xyzLocS << ", endPos = " << xyzLocE << std::endl; + + math_utils::Vector3D step(xyzLocE); + step -= xyzLocS; + step *= nStepsInv; // position increment at each step + // the electrons will injected in the middle of each step + // starting from the middle of the first step + math_utils::Vector3D stepH(step * 0.5); + xyzLocS += stepH; + xyzLocE -= stepH; + + LOG(debug) << "Step into the sensitive volume: " << step << ". Number of steps: " << nSteps; + int rowS = -1, colS = -1, rowE = -1, colE = -1, nSkip = 0; + + /// here it is the control whether the hit is in the sensitive matrix based on the segmentation + // get entrance pixel row and col + while (!Segmentation::localToDetector(xyzLocS.X(), xyzLocS.Z(), rowS, colS, subDetID, layer, disk)) { // guard-ring ? + if (++nSkip >= nSteps) { + LOG(debug) << "Did not enter to sensitive matrix, " << nSkip << " >= " << nSteps; + return; // did not enter to sensitive matrix + } + xyzLocS += step; + } + + // get exit pixel row and col + while (!Segmentation::localToDetector(xyzLocE.X(), xyzLocE.Z(), rowE, colE, subDetID, layer, disk)) { /// for the moment chipID = bigger element + if (++nSkip >= nSteps) { + LOG(debug) << "Did not enter to sensitive matrix, " << nSkip << " >= " << nSteps; + return; // did not enter to sensitive matrix + } + xyzLocE -= step; + } + + int nCols = getNCols(subDetID, layer); + int nRows = getNRows(subDetID, layer); + + // estimate the limiting min/max row and col where the non-0 response is possible + if (rowS > rowE) { + std::swap(rowS, rowE); + } + if (colS > colE) { + std::swap(colS, colE); + } + rowS -= AlpideRespSimMat::NPix / 2; + rowE += AlpideRespSimMat::NPix / 2; + if (rowS < 0) { + rowS = 0; + } + if (rowE >= nRows) { + rowE = nRows - 1; + } + colS -= AlpideRespSimMat::NPix / 2; + colE += AlpideRespSimMat::NPix / 2; + if (colS < 0) { + colS = 0; + } + if (colE >= nCols) { + colE = nCols - 1; + } + int rowSpan = rowE - rowS + 1, colSpan = colE - colS + 1; // size of plaquet where some response is expected + + float respMatrix[rowSpan][colSpan]; // response accumulated here + std::fill(&respMatrix[0][0], &respMatrix[0][0] + rowSpan * colSpan, 0.f); + + float nElectrons = hit.GetEnergyLoss() * mParams.getEnergyToNElectrons(); // total number of deposited electrons + nElectrons *= nStepsInv; // N electrons injected per step + if (nSkip) { + nSteps -= nSkip; + } + + int rowPrev = -1, colPrev = -1, row, col; + float cRowPix = 0.f, cColPix = 0.f; // local coordinate of the current pixel center + + const o2::trk::ChipSimResponse* resp = getChipResponse(chipID); + // std::cout << "Printing chip response:" << std::endl; + // resp->print(); + + // take into account that the ChipSimResponse depth defintion has different min/max boundaries + // although the max should coincide with the surface of the epitaxial layer, which in the chip + // local coordinates has Y = +SensorLayerThickness/2 + xyzLocS.SetY(xyzLocS.Y() + ((subDetID == 0) ? mSimRespVDShift : mSimRespMLOTShift)); + + // collect charge in every pixel which might be affected by the hit + for (int iStep = nSteps; iStep--;) { + // Get the pixel ID + Segmentation::localToDetector(xyzLocS.X(), xyzLocS.Z(), row, col, subDetID, layer, disk); + if (row != rowPrev || col != colPrev) { // update pixel and coordinates of its center + if (!Segmentation::detectorToLocal(row, col, cRowPix, cColPix, subDetID, layer, disk)) { + continue; // should not happen + } + rowPrev = row; + colPrev = col; + } + bool flipCol = false, flipRow = false; + // note that response needs coordinates along column row (locX) (locZ) then depth (locY) + float rowMax{}, colMax{}; + const AlpideRespSimMat* rspmat{nullptr}; + if (subDetID == 0) { // VD + rowMax = 0.5f * Segmentation::PitchRowVD * mSimRespVDScaleX; + colMax = 0.5f * Segmentation::PitchColVD * mSimRespVDScaleZ; + rspmat = resp->getResponse(mSimRespVDScaleX * (xyzLocS.X() - cRowPix), mSimRespVDScaleZ * (xyzLocS.Z() - cColPix), xyzLocS.Y(), flipRow, flipCol, rowMax, colMax); + } else { // ML/OT + rowMax = 0.5f * Segmentation::PitchRowMLOT * mSimRespMLOTScaleX; + colMax = 0.5f * Segmentation::PitchColMLOT * mSimRespMLOTScaleZ; + rspmat = resp->getResponse(mSimRespMLOTScaleX * (xyzLocS.X() - cRowPix), mSimRespMLOTScaleZ * (xyzLocS.Z() - cColPix), xyzLocS.Y(), flipRow, flipCol, rowMax, colMax); + } + + float tempPitchX = 0, tempPitchZ = 0; + if (subDetID == 0) { + tempPitchX = Segmentation::PitchRowVD; + tempPitchZ = Segmentation::PitchColVD; + } else { + tempPitchX = Segmentation::PitchRowMLOT; + tempPitchZ = Segmentation::PitchColMLOT; + } + LOG(debug) << "X and Z inside pixel at start = " << (xyzLocS.X() - cRowPix) << " , " << (xyzLocS.Z() - cColPix) << ", rescaled: " << mSimRespMLOTScaleX * (xyzLocS.X() - cRowPix) << " , " << mSimRespMLOTScaleZ * (xyzLocS.Z() - cColPix); + LOG(debug) << "Hit inside pitch? X: " << ((xyzLocS.X() - cRowPix) < tempPitchX) << " Z: " << ((xyzLocS.Z() - cColPix) < tempPitchZ); + + xyzLocS += step; + + if (rspmat == nullptr) { + LOG(debug) << "Error in rspmat for step " << iStep << " / " << nSteps; + continue; + } + LOG(debug) << "rspmat valid! for step " << iStep << " / " << nSteps << ", (row,col) = (" << row << "," << col << ")"; + // rspmat->print(); // print the response matrix for debugging + + for (int irow = AlpideRespSimMat::NPix; irow--;) { + int rowDest = row + irow - AlpideRespSimMat::NPix / 2 - rowS; // destination row in the respMatrix + if (rowDest < 0 || rowDest >= rowSpan) { + continue; + } + for (int icol = AlpideRespSimMat::NPix; icol--;) { + int colDest = col + icol - AlpideRespSimMat::NPix / 2 - colS; // destination column in the respMatrix + if (colDest < 0 || colDest >= colSpan) { + continue; + } + respMatrix[rowDest][colDest] += rspmat->getValue(irow, icol, mSimRespOrientation ? !flipRow : flipRow, flipCol); + } + } + } + + // fire the pixels assuming Poisson(n_response_electrons) + o2::MCCompLabel lbl(hit.GetTrackID(), evID, srcID, false); + auto roFrameAbs = mNewROFrame + roFrameRel; + LOG(debug) << "Spanning through rows and columns; rowspan = " << rowSpan << " colspan = " << colSpan << " = " << colE << " - " << colS << " +1 " << std::endl; + for (int irow = rowSpan; irow--;) { // irow ranging from 4 to 0 + uint16_t rowIS = irow + rowS; // row distant irow from the row of the hit start + for (int icol = colSpan; icol--;) { // icol ranging from 4 to 0 + float nEleResp = respMatrix[irow][icol]; // value of the probability of the response in this pixel + if (nEleResp <= 1.e-36) { + continue; + } + LOG(debug) << "nEleResp: value " << nEleResp << " for pixel " << irow << " " << icol << std::endl; + int nEle = gRandom->Poisson(nElectrons * nEleResp); // total charge in given pixel = number of electrons generated in the hit multiplied by the probability of being detected in their position + LOG(debug) << "Charge detected in the pixel: " << nEle << " for pixel " << irow << " " << icol << std::endl; + // ignore charge which have no chance to fire the pixel + if (nEle < mParams.getMinChargeToAccount()) { /// TODO: substitute with the threshold? + LOG(debug) << "Ignoring pixel with nEle = " << nEle << " < min charge to account " + << mParams.getMinChargeToAccount() << " for pixel " << irow << " " << icol; + continue; + } + + uint16_t colIS = icol + colS; // col distant icol from the col of the hit start + if (mNoiseMap && mNoiseMap->isNoisy(chipID, rowIS, colIS)) { + continue; + } + if (mDeadChanMap && mDeadChanMap->isNoisy(chipID, rowIS, colIS)) { + continue; + } + + registerDigits(chip, roFrameAbs, timeInROF, nFrames, rowIS, colIS, nEle, lbl); + } + } } //________________________________________________________________________________ -void Digitizer::registerDigits(o2::itsmft::ChipDigitsContainer& chip, uint32_t roFrame, float tInROF, int nROF, +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) { // 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 - - // float tStrobe = mParams.getStrobeDelay() - tInROF; // strobe start wrt signal start - // for (int i = 0; i < nROF; i++) { - // uint32_t roFr = roFrame + i; - // int nEleROF = mParams.getSignalShape().getCollectedCharge(nEle, tStrobe, tStrobe + mParams.getStrobeLength()); - // tStrobe += mParams.getROFrameLength(); // for the next ROF - - // // discard too small contributions, they have no chance to produce a digit - // if (nEleROF < mParams.getMinChargeToAccount()) { - // continue; - // } - // if (roFr > mEventROFrameMax) { - // mEventROFrameMax = roFr; - // } - // if (roFr < mEventROFrameMin) { - // mEventROFrameMin = roFr; - // } - // auto key = chip.getOrderingKey(roFr, row, col); - // PreDigit* pd = chip.findDigit(key); - // if (!pd) { - // chip.addDigit(key, roFr, row, col, nEleROF, lbl); - // } else { // there is already a digit at this slot, account as PreDigitExtra contribution - // pd->charge += nEleROF; - // if (pd->labelRef.label == lbl) { // don't store the same label twice - // continue; - // } - // ExtraDig* extra = getExtraDigBuffer(roFr); - // int& nxt = pd->labelRef.next; - // bool skip = false; - // while (nxt >= 0) { - // if ((*extra)[nxt].label == lbl) { // don't store the same label twice - // skip = true; - // break; - // } - // nxt = (*extra)[nxt].next; - // } - // if (skip) { - // continue; - // } - // // new predigit will be added in the end of the chain - // nxt = extra->size(); - // extra->emplace_back(lbl); - // } - // } + 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++) { + uint32_t roFr = roFrame + i; + int nEleROF = mParams.getSignalShape().getCollectedCharge(nEle, tStrobe, tStrobe + mParams.getStrobeLength()); + tStrobe += mParams.getROFrameLength(); // for the next ROF + + // discard too small contributions, they have no chance to produce a digit + if (nEleROF < mParams.getMinChargeToAccount()) { /// use threshold instead? + continue; + } + if (roFr > mEventROFrameMax) { + mEventROFrameMax = roFr; + } + if (roFr < mEventROFrameMin) { + mEventROFrameMin = roFr; + } + auto key = chip.getOrderingKey(roFr, row, col); + o2::itsmft::PreDigit* pd = chip.findDigit(key); + if (!pd) { + chip.addDigit(key, roFr, row, col, nEleROF, lbl); + LOG(debug) << "Added digit " << key << " " << roFr << " " << row << " " << col << " " << nEleROF; + } else { // there is already a digit at this slot, account as PreDigitExtra contribution + pd->charge += nEleROF; + if (pd->labelRef.label == lbl) { // don't store the same label twice + continue; + } + ExtraDig* extra = getExtraDigBuffer(roFr); + int& nxt = pd->labelRef.next; + bool skip = false; + while (nxt >= 0) { + if ((*extra)[nxt].label == lbl) { // don't store the same label twice + skip = true; + break; + } + nxt = (*extra)[nxt].next; + } + if (skip) { + continue; + } + // new predigit will be added in the end of the chain + nxt = extra->size(); + extra->emplace_back(lbl); + } + } } diff --git a/Steer/DigitizerWorkflow/src/TRKDigitizerSpec.cxx b/Steer/DigitizerWorkflow/src/TRKDigitizerSpec.cxx index cb375936744d5..0ed276237bd86 100644 --- a/Steer/DigitizerWorkflow/src/TRKDigitizerSpec.cxx +++ b/Steer/DigitizerWorkflow/src/TRKDigitizerSpec.cxx @@ -93,7 +93,7 @@ class TRKDPLDigitizerTask : BaseDPLDigitizer timer.Start(); LOG(info) << " CALLING TRK DIGITIZATION "; - // mDigitizer.setDigits(&mDigits); + mDigitizer.setDigits(&mDigits); mDigitizer.setROFRecords(&mROFRecords); mDigitizer.setMCLabels(&mLabels); @@ -104,8 +104,10 @@ class TRKDPLDigitizerTask : BaseDPLDigitizer // 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)); @@ -139,7 +141,7 @@ class TRKDPLDigitizerTask : BaseDPLDigitizer mLabels.clear(); mDigits.clear(); mROFRecords.clear(); - }; // and accumulate lambda + }; // end accumulate lambda auto& eventParts = context->getEventParts(withQED); // loop over all composite collisions given from context (aka loop over all the interaction records) @@ -172,6 +174,7 @@ class TRKDPLDigitizerTask : BaseDPLDigitizer 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)