From 850a4dc9ed9a7f8db8b5c53ff1372f500f5d306a Mon Sep 17 00:00:00 2001 From: Archil Durglishvili Date: Wed, 30 Oct 2024 15:50:08 +0100 Subject: [PATCH 01/17] HCal new readout implementation --- .../compact/ALLEGRO_o1_v03/ALLEGRO_o1_v03.xml | 4 +- .../ALLEGRO_o1_v03/HCalBarrel_TileCal_v03.xml | 133 ++++++++ .../HCalEndcaps_ThreeParts_TileCal_v03.xml | 145 +++++++++ .../FCCSWHCalPhiThetaHandle_k4geo.h | 131 ++++++++ .../FCCSWHCalPhiTheta_k4geo.h | 193 ++++++++++++ .../src/FCCSWHCalPhiThetaHandle_k4geo.cpp | 4 + .../src/FCCSWHCalPhiTheta_k4geo.cpp | 284 ++++++++++++++++++ .../src/plugins/SegmentationFactories.cpp | 3 + 8 files changed, 895 insertions(+), 2 deletions(-) create mode 100644 FCCee/ALLEGRO/compact/ALLEGRO_o1_v03/HCalBarrel_TileCal_v03.xml create mode 100644 FCCee/ALLEGRO/compact/ALLEGRO_o1_v03/HCalEndcaps_ThreeParts_TileCal_v03.xml create mode 100644 detectorSegmentations/include/detectorSegmentations/FCCSWHCalPhiThetaHandle_k4geo.h create mode 100644 detectorSegmentations/include/detectorSegmentations/FCCSWHCalPhiTheta_k4geo.h create mode 100644 detectorSegmentations/src/FCCSWHCalPhiThetaHandle_k4geo.cpp create mode 100644 detectorSegmentations/src/FCCSWHCalPhiTheta_k4geo.cpp diff --git a/FCCee/ALLEGRO/compact/ALLEGRO_o1_v03/ALLEGRO_o1_v03.xml b/FCCee/ALLEGRO/compact/ALLEGRO_o1_v03/ALLEGRO_o1_v03.xml index fc2ca3b09..0d478c06e 100644 --- a/FCCee/ALLEGRO/compact/ALLEGRO_o1_v03/ALLEGRO_o1_v03.xml +++ b/FCCee/ALLEGRO/compact/ALLEGRO_o1_v03/ALLEGRO_o1_v03.xml @@ -45,9 +45,9 @@ - + - + diff --git a/FCCee/ALLEGRO/compact/ALLEGRO_o1_v03/HCalBarrel_TileCal_v03.xml b/FCCee/ALLEGRO/compact/ALLEGRO_o1_v03/HCalBarrel_TileCal_v03.xml new file mode 100644 index 000000000..cfcc37eac --- /dev/null +++ b/FCCee/ALLEGRO/compact/ALLEGRO_o1_v03/HCalBarrel_TileCal_v03.xml @@ -0,0 +1,133 @@ + + + + The first FCCee HCal layout based on ATLAS HCal, not optimised yet + 1. Nov 2022, J. Faltova: update material and radial segmentation for FCCee + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + system:4,layer:5,row:9,theta:9,phi:10 + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/FCCee/ALLEGRO/compact/ALLEGRO_o1_v03/HCalEndcaps_ThreeParts_TileCal_v03.xml b/FCCee/ALLEGRO/compact/ALLEGRO_o1_v03/HCalEndcaps_ThreeParts_TileCal_v03.xml new file mode 100644 index 000000000..5d64a9b5a --- /dev/null +++ b/FCCee/ALLEGRO/compact/ALLEGRO_o1_v03/HCalEndcaps_ThreeParts_TileCal_v03.xml @@ -0,0 +1,145 @@ + + + + HCal layout based on ATLAS HCal, with realistic longitudinal segmentation and steel support + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + system:4,type:3,layer:6,row:11,theta:11,phi:10 + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/detectorSegmentations/include/detectorSegmentations/FCCSWHCalPhiThetaHandle_k4geo.h b/detectorSegmentations/include/detectorSegmentations/FCCSWHCalPhiThetaHandle_k4geo.h new file mode 100644 index 000000000..e787bbd45 --- /dev/null +++ b/detectorSegmentations/include/detectorSegmentations/FCCSWHCalPhiThetaHandle_k4geo.h @@ -0,0 +1,131 @@ +#ifndef DETECTORSEGMENTATIONS_HCALPHITHETAHANDLE_K4GEO_H +#define DETECTORSEGMENTATIONS_HCALPHITHETAHANDLE_K4GEO_H + +// FCCSW +#include "detectorSegmentations/FCCSWHCalPhiTheta_k4geo.h" + +// DD4hep +#include "DD4hep/Segmentations.h" +#include "DD4hep/detail/SegmentationsInterna.h" + +/// Namespace for the AIDA detector description toolkit +namespace dd4hep { +/// Namespace for base segmentations + + +// Forward declarations +class Segmentation; +template class SegmentationWrapper; + +/// We need some abbreviation to make the code more readable. +typedef Handle> FCCSWGridPhiThetaHandle_k4geo; + +/// Implementation class for the HCal phi-theta segmentation. +/** + * Concrete user handle to serve specific needs of client code + * which requires access to the base functionality not served + * by the super-class Segmentation. + * + * Note: + * We only check the validity of the underlying handle. + * If for whatever reason the implementation object is not valid + * This is not checked. + * In principle this CANNOT happen unless some brain-dead has + * fiddled with the handled object directly..... + * + * Note: + * The handle base corrsponding to this object in for + * conveniance reasons instantiated in DD4hep/src/Segmentations.cpp. + * + */ +class FCCSWHCalPhiTheta_k4geo : public FCCSWGridPhiThetaHandle_k4geo { +public: + /// Defintion of the basic handled object + typedef FCCSWGridPhiThetaHandle_k4geo::Object Object; + +public: + /// Default constructor + FCCSWHCalPhiTheta_k4geo() = default; + /// Copy constructor + FCCSWHCalPhiTheta_k4geo(const FCCSWHCalPhiTheta_k4geo& e) = default; + /// Copy Constructor from segmentation base object + FCCSWHCalPhiTheta_k4geo(const Segmentation& e) : Handle(e) {} + /// Copy constructor from handle + FCCSWHCalPhiTheta_k4geo(const Handle& e) : Handle(e) {} + /// Copy constructor from other polymorph/equivalent handle + template + FCCSWHCalPhiTheta_k4geo(const Handle& e) : Handle(e) {} + /// Assignment operator + FCCSWHCalPhiTheta_k4geo& operator=(const FCCSWHCalPhiTheta_k4geo& seg) = default; + /// Equality operator + bool operator==(const FCCSWHCalPhiTheta_k4geo& seg) const { return m_element == seg.m_element; } + + /// determine the position based on the cell ID + inline Position position(const CellID& id) const { return Position(access()->implementation->position(id)); } + + /// determine the cell ID based on the position + inline dd4hep::CellID cellID(const Position& local, const Position& global, const VolumeID& volID) const { + return access()->implementation->cellID(local, global, volID); + } + + /// access the grid size in theta + inline double gridSizeTheta() const { return access()->implementation->gridSizeTheta(); } + + /// access the grid size in Phi + inline int phiBins() const { return access()->implementation->phiBins(); } + + /// access the coordinate offset in theta + inline double offsetTheta() const { return access()->implementation->offsetTheta(); } + + /// access the coordinate offset in Phi + inline double offsetPhi() const { return access()->implementation->offsetPhi(); } + + /// set the coordinate offset in theta + inline void setOffsetTheta(double offset) const { access()->implementation->setOffsetTheta(offset); } + + /// set the coordinate offset in Phi + inline void setOffsetPhi(double offset) const { access()->implementation->setOffsetPhi(offset); } + + /// set the coordinate offset in z-axis + inline void setOffsetZ(std::vector const&offset) const { access()->implementation->setOffsetZ(offset); } + + /// set the z width + inline void setWidthZ(std::vector const&width) const { access()->implementation->setWidthZ(width); } + + /// set the offset in radius + inline void setOffsetR(std::vector const&offset) const { access()->implementation->setOffsetR(offset); } + + /// set the number of layers with different dR + inline void setNumLayers(std::vector const&num) const { access()->implementation->setNumLayers(num); } + + /// set the dR of each layer + inline void setdRlayer(std::vector const&dRlayer) const { access()->implementation->setdRlayer(dRlayer); } + + /// set the grid size in theta + inline void setGridSizeTheta(double cellSize) const { access()->implementation->setGridSizeTheta(cellSize); } + + /// set the grid size in Phi + inline void setPhiBins(int cellSize) const { access()->implementation->setPhiBins(cellSize); } + + /// access the field name used for theta + inline const std::string& fieldNameTheta() const { return access()->implementation->fieldNameTheta(); } + + /// access the field name used for Phi + inline const std::string& fieldNamePhi() const { return access()->implementation->fieldNamePhi(); } + + /** \brief Returns a std::vector of the cellDimensions of the given cell ID + in natural order of dimensions (dPhi, dTheta) + + Returns a std::vector of the cellDimensions of the given cell ID + \param cellID is ignored as all cells have the same dimension + \return std::vector size 2: + -# size in phi + -# size in theta + */ + inline std::vector cellDimensions(const CellID& /*id*/) const { + return {access()->implementation->gridSizePhi(), access()->implementation->gridSizeTheta()}; + } +}; + +} /* End namespace dd4hep */ +#endif // DETECTORSEGMENTATIONS_HCALPHITHETAHANDLE_K4GEO_H diff --git a/detectorSegmentations/include/detectorSegmentations/FCCSWHCalPhiTheta_k4geo.h b/detectorSegmentations/include/detectorSegmentations/FCCSWHCalPhiTheta_k4geo.h new file mode 100644 index 000000000..439d56ab7 --- /dev/null +++ b/detectorSegmentations/include/detectorSegmentations/FCCSWHCalPhiTheta_k4geo.h @@ -0,0 +1,193 @@ +#ifndef DETECTORSEGMENTATIONS_HCALPHITHETA_K4GEO_H +#define DETECTORSEGMENTATIONS_HCALPHITHETA_K4GEO_H + +// FCCSW +#include "detectorSegmentations/GridTheta_k4geo.h" + +/** FCCSWHCalPhiTheta_k4geo Detector/detectorSegmentations/detectorSegmentations/FCCSWHCalPhiTheta_k4geo.h FCCSWHCalPhiTheta_k4geo.h + * + * Segmentation in theta and phi. + * Based on GridTheta_k4geo, addition of azimuthal angle coordinate. + * + * Rectangular shape cells are defined in r-z plan and all the hits within a defined cell boundaries are assigned the same cellID. + * + */ + +namespace dd4hep { + namespace DDSegmentation { + class FCCSWHCalPhiTheta_k4geo : public GridTheta_k4geo { + public: + /// default constructor using an arbitrary type + FCCSWHCalPhiTheta_k4geo(const std::string& aCellEncoding); + /// Default constructor used by derived classes passing an existing decoder + FCCSWHCalPhiTheta_k4geo(const BitFieldCoder* decoder); + + /// destructor + virtual ~FCCSWHCalPhiTheta_k4geo() = default; + + /** Determine the position of HCal cell based on the cellID. + * @warning This segmentation has no knowledge of radius, so radius = 1 is taken into calculations. + * @param[in] aCellId ID of a cell. + * return Position (radius = 1). + */ + virtual Vector3D position(const CellID& aCellID) const; + + /** Assign a cellID based on the global position. + * @param[in] aLocalPosition (not used). + * @param[in] aGlobalPosition position in the global coordinates. + * @param[in] aVolumeId ID of a volume. + * return Cell ID. + */ + virtual CellID cellID(const Vector3D& aLocalPosition, const Vector3D& aGlobalPosition, + const VolumeID& aVolumeID) const; + + /** Calculate layer radii and edges in z-axis, then define cell edges in each layer using defineCellEdges(). + */ + void defineCellsInRZplan() const; + + /** Define cell edges in z-axis for the given layer. + * Logic: + * 1) Find theta bin centers that fit within the given layer; + * 2) Define a cell edge in z-axis as the middle of each pair of theta bin centers + * @param[in] layer index + */ + void defineCellEdges(const unsigned int layer) const; + + /** Determine the azimuthal angle of HCal cell based on the cellID. + * @param[in] aCellId ID of a cell. + * return Phi. + */ + double phi(const CellID& aCellID) const; + + /** Get the grid size in phi. + * return Grid size in phi. + */ + inline double gridSizePhi() const { return 2 * M_PI / static_cast(m_phiBins); } + + /** Get the number of bins in azimuthal angle. + * return Number of bins in phi. + */ + inline int phiBins() const { return m_phiBins; } + + /** Get the coordinate offset in azimuthal angle. + * return The offset in phi. + */ + inline double offsetPhi() const { return m_offsetPhi; } + + /** Get the vector of theta bins (cells) in a give layer. + */ + inline std::vector thetaBins(const uint layer) const { + if(m_radii.empty()) defineCellsInRZplan(); + if(!m_thetaBins.empty()) return m_thetaBins[layer]; + else return std::vector(); + } + + /** Get the coordinate offset in z-axis. + * return The offset in z. + */ + inline std::vector offsetZ() const { return m_offsetZ; } + + /** Get the z width of the layer. + * return the z width. + */ + inline std::vector widthZ() const { return m_widthZ; } + + /** Get the coordinate offset in radius. + * return the offset in radius. + */ + inline std::vector offsetR() const { return m_offsetR; } + + /** Get the number of layers. + * return the number of layers. + */ + inline std::vector numLayers() const { return m_numLayers; } + + /** Get the dR of layers. + * return the dR. + */ + inline std::vector dRlayer() const { return m_dRlayer; } + + /** Get the field name for azimuthal angle. + * return The field name for phi. + */ + inline const std::string& fieldNamePhi() const { return m_phiID; } + + /** Get the postion og the geometric center of the cell based on the cellID + * @param[in] aCellID + * return the global coordinates of cell center + */ + inline Vector3D centerPosition(const CellID& cID) const; + + + /** Set the number of bins in azimuthal angle. + * @param[in] aNumberBins Number of bins in phi. + */ + inline void setPhiBins(int bins) { m_phiBins = bins; } + + /** Set the coordinate offset in azimuthal angle. + * @param[in] aOffset Offset in phi. + */ + inline void setOffsetPhi(double offset) { m_offsetPhi = offset; } + + /** Set the coordinate offset in z-axis. + * @param[in] offset in z (centre of the layer). + */ + inline void setOffsetZ(std::vector const&offset) { m_offsetZ = offset; } + + /** Set the z width. + * @param[in] width in z. + */ + inline void setWidthZ(std::vector const&width) { m_widthZ = width; } + + /** Set the coordinate offset in radius. + * @param[in] offset in radius. + */ + inline void setOffsetR(std::vector const&offset) { m_offsetR = offset; } + + /** Set the number of layers. + * @param[in] number of layers + */ + inline void setNumLayers(std::vector const&num) { m_numLayers = num; } + + /** Set the dR of layers. + * @param[in] dR of layers. + */ + inline void setdRlayer(std::vector const&dRlayer) { m_dRlayer = dRlayer; } + + + /** Set the field name used for azimuthal angle. + * @param[in] aFieldName Field name for phi. + */ + inline void setFieldNamePhi(const std::string& fieldName) { m_phiID = fieldName; } + + protected: + /// determine the azimuthal angle phi based on the current cell ID + double phi() const; + /// the number of bins in phi + int m_phiBins; + /// the coordinate offset in phi + double m_offsetPhi; + /// the field name used for phi + std::string m_phiID; + /// the z offset of middle of the layer + std::vector m_offsetZ; + /// the z width of the layer + std::vector m_widthZ; + /// the radius of the layer + std::vector m_offsetR; + /// number of layers + std::vector m_numLayers; + /// dR of the layer + std::vector m_dRlayer; + /// radius of each layer + mutable std::vector m_radii; + /// z-min and z-max of each layer + mutable std::vector > m_layerEdges; + /// theta bins (cells) in each layer + mutable std::vector > m_thetaBins; + /// z-min and z-max of each cell (theta bin) in each layer + mutable std::vector > > m_cellEdges; + }; + } +} +#endif /* DETSEGMENTATION_HCALPHITHETA_H */ diff --git a/detectorSegmentations/src/FCCSWHCalPhiThetaHandle_k4geo.cpp b/detectorSegmentations/src/FCCSWHCalPhiThetaHandle_k4geo.cpp new file mode 100644 index 000000000..1e3f8d0e4 --- /dev/null +++ b/detectorSegmentations/src/FCCSWHCalPhiThetaHandle_k4geo.cpp @@ -0,0 +1,4 @@ +#include "detectorSegmentations/FCCSWHCalPhiThetaHandle_k4geo.h" +#include "DD4hep/detail/Handle.inl" + +DD4HEP_INSTANTIATE_HANDLE_UNNAMED(dd4hep::DDSegmentation::FCCSWHCalPhiTheta_k4geo); diff --git a/detectorSegmentations/src/FCCSWHCalPhiTheta_k4geo.cpp b/detectorSegmentations/src/FCCSWHCalPhiTheta_k4geo.cpp new file mode 100644 index 000000000..cc49d9e32 --- /dev/null +++ b/detectorSegmentations/src/FCCSWHCalPhiTheta_k4geo.cpp @@ -0,0 +1,284 @@ +#include "detectorSegmentations/FCCSWHCalPhiTheta_k4geo.h" +#include "DD4hep/Printout.h" + +namespace dd4hep { +namespace DDSegmentation { + +/// default constructor using an encoding string +FCCSWHCalPhiTheta_k4geo::FCCSWHCalPhiTheta_k4geo(const std::string& cellEncoding) : GridTheta_k4geo(cellEncoding) { + // define type and description + _type = "FCCSWHCalPhiTheta_k4geo"; + _description = "Phi-theta segmentation in the global coordinates"; + + // register all necessary parameters (additional to those registered in GridTheta_k4geo) + registerParameter("phi_bins", "Number of bins phi", m_phiBins, 1); + registerParameter("offset_phi", "Angular offset in phi", m_offsetPhi, 0., SegmentationParameter::AngleUnit, true); + registerParameter("offset_z", "Offset in z-axis of the layer center", m_offsetZ, std::vector()); + registerParameter("width_z", "Width in z of the layer", m_widthZ, std::vector()); + registerParameter("offset_r", "Offset in radius of the layer (Rmin)", m_offsetR, std::vector()); + registerParameter("numLayers", "Number of layers", m_numLayers, std::vector()); + registerParameter("dRlayer", "dR of the layer", m_dRlayer, std::vector()); + registerIdentifier("identifier_phi", "Cell ID identifier for phi", m_phiID, "phi"); +} + +FCCSWHCalPhiTheta_k4geo::FCCSWHCalPhiTheta_k4geo(const BitFieldCoder* decoder) : GridTheta_k4geo(decoder) { + // define type and description + _type = "FCCSWHCalPhiTheta_k4geo"; + _description = "Phi-theta segmentation in the global coordinates"; + + // register all necessary parameters (additional to those registered in GridTheta_k4geo) + registerParameter("phi_bins", "Number of bins phi", m_phiBins, 1); + registerParameter("offset_phi", "Angular offset in phi", m_offsetPhi, 0., SegmentationParameter::AngleUnit, true); + registerParameter("offset_z", "Offset in z-axis of the layer center", m_offsetZ, std::vector()); + registerParameter("width_z", "Width in z of the layer", m_widthZ, std::vector()); + registerParameter("offset_r", "Offset in radius of the layer (Rmin)", m_offsetR, std::vector()); + registerParameter("numLayers", "Number of layers", m_numLayers, std::vector()); + registerParameter("dRlayer", "dR of the layer", m_dRlayer, std::vector()); + registerIdentifier("identifier_phi", "Cell ID identifier for phi", m_phiID, "phi"); +} + + +/// determine the global position based on the cell ID +Vector3D FCCSWHCalPhiTheta_k4geo::position(const CellID& cID) const { + uint layer = _decoder->get(cID,"layer"); + double radius = 1.0; + + if(m_radii.empty()) defineCellsInRZplan(); + if(!m_radii.empty()) radius = m_radii[layer]; + + return positionFromRThetaPhi(radius, theta(cID), phi(cID)); +} + + +/// determine the global position based on the cell ID +/// returns the geometric center of the cell +Vector3D FCCSWHCalPhiTheta_k4geo::centerPosition(const CellID& cID) const { + uint layer = _decoder->get(cID,"layer"); + int thetaID = _decoder->get(cID,"theta"); + double zpos = 0.; + double radius = 1.0; + + if(m_radii.empty()) defineCellsInRZplan(); + if(!m_radii.empty()) radius = m_radii[layer]; + if(!m_cellEdges.empty()) zpos = m_cellEdges[layer][thetaID].first + (m_cellEdges[layer][thetaID].second - m_cellEdges[layer][thetaID].first) * 0.5; + + auto pos = positionFromRThetaPhi(radius, theta(cID), phi(cID)); + + // retrun the position with corrected z corrdinate to match to the geometric center + return Vector3D(pos.x(),pos.y(),zpos); +} + + +void FCCSWHCalPhiTheta_k4geo::defineCellsInRZplan() const { + if(m_radii.empty()) + { + // check if all necessary variables are available + if(m_offsetZ.empty() || m_widthZ.empty() || m_offsetR.empty() || m_numLayers.empty() || m_dRlayer.empty()) + { + dd4hep::printout(dd4hep::ERROR, "FCCSWHCalPhiTheta_k4geo","Please check the readout description in the XML file!", + "One of the variables is missing: offset_z | width_z | offset_r | numLayers | dRlayer"); + return; + } + + // calculate the radius for each layer + if(m_offsetZ.size() == 1) // Barrel + { + dd4hep::printout(dd4hep::INFO, "FCCSWHCalPhiTheta_k4geo","Barrel configuration found!"); + uint N_dR = m_numLayers.size(); + double moduleDepth = 0.; + for(uint i_dR = 0; i_dR < N_dR; i_dR++) + { + for(int i_row = 1; i_row <= m_numLayers[i_dR]; i_row++) + { + moduleDepth+=m_dRlayer[i_dR]; + m_radii.push_back(m_offsetR[0] + moduleDepth - m_dRlayer[i_dR]*0.5); + // layer lower and upper edges in z-axis + m_layerEdges.push_back( std::make_pair(m_offsetZ[0] - 0.5*m_widthZ[0], m_offsetZ[0] + 0.5*m_widthZ[0]) ); + } + } + } // Barrel + + if(m_offsetZ.size() > 1) // ThreePartsEndCap + { + dd4hep::printout(dd4hep::INFO, "FCCSWHCalPhiTheta_k4geo","ThreePartsEndCap configuration found!"); + uint N_dR = m_numLayers.size()/3; + double moduleDepth1 = 0.; + double moduleDepth2 = 0.; + double moduleDepth3 = 0.; + // part1 + for(uint i_dR = 0; i_dR < N_dR; i_dR++) + { + for(int i_row = 1; i_row <= m_numLayers[i_dR]; i_row++) + { + moduleDepth1+=m_dRlayer[i_dR]; + m_radii.push_back(m_offsetR[0] + moduleDepth1 - m_dRlayer[i_dR]*0.5); + // layer lower and upper edges in z-axis + m_layerEdges.push_back( std::make_pair(m_offsetZ[0] - 0.5*m_widthZ[0], m_offsetZ[0] + 0.5*m_widthZ[0]) ); + } + } + // part2 + for(uint i_dR = 0; i_dR < N_dR; i_dR++) + { + for(int i_row = 1; i_row <= m_numLayers[i_dR + N_dR]; i_row++) + { + moduleDepth2+=m_dRlayer[i_dR]; + m_radii.push_back(m_offsetR[1] + moduleDepth2 - m_dRlayer[i_dR]*0.5); + // layer lower and upper edges in z-axis + m_layerEdges.push_back( std::make_pair(m_offsetZ[1] - 0.5*m_widthZ[1], m_offsetZ[1] + 0.5*m_widthZ[1]) ); + } + } + // part3 + for(uint i_dR = 0; i_dR < N_dR; i_dR++) + { + for(int i_row = 1; i_row <= m_numLayers[i_dR + 2*N_dR]; i_row++) + { + moduleDepth3+=m_dRlayer[i_dR]; + m_radii.push_back(m_offsetR[2] + moduleDepth3 - m_dRlayer[i_dR]*0.5); + // layer lower and upper edges in z-axis + m_layerEdges.push_back( std::make_pair(m_offsetZ[2] - 0.5*m_widthZ[2], m_offsetZ[2] + 0.5*m_widthZ[2]) ); + } + } + } // ThreePartsEndcap + + // print info of calculated radii and edges + for(uint i_layer = 0; i_layer < m_radii.size(); i_layer++){ + dd4hep::printout(dd4hep::INFO, "FCCSWHCalPhiTheta_k4geo","layer %d radius: %.2f, z range: %.2f - %.2f cm", + i_layer, m_radii[i_layer], m_layerEdges[i_layer].first, m_layerEdges[i_layer].second); + } + + // allocate thetaBins vector for each layer + m_thetaBins.resize(m_radii.size()); + // allocate cellEdges vector for each layer + m_cellEdges.resize(m_radii.size()); + + // determine theta bins and cell edges for each layer + for(uint i_layer = 0; i_layer < m_radii.size(); i_layer++) defineCellEdges(i_layer); + } +} + + +void FCCSWHCalPhiTheta_k4geo::defineCellEdges(const uint layer) const { + if(m_thetaBins[layer].size()==0 && m_radii.size() > 0) + { + //find theta bins that fit within the given layer + int ibin = positionToBin(0.02, m_gridSizeTheta, m_offsetTheta); // <--- start from theta bin outside the HCal theta range + while(m_radii[layer]*std::cos(m_offsetTheta+ibin*m_gridSizeTheta)/std::sin(m_offsetTheta+ibin*m_gridSizeTheta) > m_layerEdges[layer].first) + { + if(m_radii[layer]*std::cos(m_offsetTheta+ibin*m_gridSizeTheta)/std::sin(m_offsetTheta+ibin*m_gridSizeTheta) < m_layerEdges[layer].second) + { + m_thetaBins[layer].push_back(ibin); + } + ibin++; + } + + // find edges of each cell (theta bin) in the given layer + auto prevBin = m_thetaBins[layer][0]; + // set the upper edge of the first cell in the given layer (starting from positive z part) + m_cellEdges[layer][prevBin] = std::make_pair(0.,m_layerEdges[layer].second); + for(auto bin : m_thetaBins[layer]) + { + if(bin!=prevBin) + { + double z1 = m_radii[layer]*std::cos(m_offsetTheta+bin*m_gridSizeTheta)/std::sin(m_offsetTheta+bin*m_gridSizeTheta); + double z2 = m_radii[layer]*std::cos(m_offsetTheta+prevBin*m_gridSizeTheta)/std::sin(m_offsetTheta+prevBin*m_gridSizeTheta); + // set the lower edge of the prevBin cell + m_cellEdges[layer][prevBin].first = z1 + 0.5*(z2 - z1); + // set the upper edge of current bin cell + m_cellEdges[layer][bin] = std::make_pair(0.,m_cellEdges[layer][prevBin].first); + prevBin = bin; + } + } + // set the lower edge of the last cell in the given layer + m_cellEdges[layer][prevBin].first = m_layerEdges[layer].first; + + // for the EndCap, do it again but for negative z part + if(m_offsetZ.size() > 1) + { + while(m_radii[layer]*std::cos(m_offsetTheta+ibin*m_gridSizeTheta)/std::sin(m_offsetTheta+ibin*m_gridSizeTheta) > (-m_layerEdges[layer].second)) + { + if(m_radii[layer]*std::cos(m_offsetTheta+ibin*m_gridSizeTheta)/std::sin(m_offsetTheta+ibin*m_gridSizeTheta) < (-m_layerEdges[layer].first)) + { + m_thetaBins[layer].push_back(ibin); + } + ibin++; + } + + // Create a span view over the theta bins corresponding to the Endcap in negative z part + std::span thetaBins(m_thetaBins[layer].begin() + m_thetaBins[layer].size()/2, m_thetaBins[layer].end()); + prevBin = thetaBins[0]; + + // set the upper edge of the first cell in the given layer at negative z part + m_cellEdges[layer][prevBin] = std::make_pair(0.,-m_layerEdges[layer].first); + for(auto bin : thetaBins) + { + if(bin!=prevBin) + { + double z1 = m_radii[layer]*std::cos(m_offsetTheta+bin*m_gridSizeTheta)/std::sin(m_offsetTheta+bin*m_gridSizeTheta); + double z2 = m_radii[layer]*std::cos(m_offsetTheta+prevBin*m_gridSizeTheta)/std::sin(m_offsetTheta+prevBin*m_gridSizeTheta); + // set the lower edge of the prevBin cell + m_cellEdges[layer][prevBin].first = z1 + 0.5*(z2 - z1); + // set the upper edge of current bin cell + m_cellEdges[layer][bin] = std::make_pair(0.,m_cellEdges[layer][prevBin].first); + prevBin = bin; + } + } + // set the lower edge of the last cell in the given layer + m_cellEdges[layer][prevBin].first = (-m_layerEdges[layer].second); + }// negative-z endcap + + dd4hep::printout(dd4hep::DEBUG, "FCCSWHCalPhiTheta_k4geo","Number of cells in layer %d: %d", layer, m_thetaBins[layer].size()); + for(auto bin : m_thetaBins[layer]) dd4hep::printout(dd4hep::DEBUG, "FCCSWHCalPhiTheta_k4geo","Layer %d cell theta bin: %d, edges: %.2f - %.2f cm", layer, bin, m_cellEdges[layer][bin]); + } +} + +/// create the cell ID based on the position +CellID FCCSWHCalPhiTheta_k4geo::cellID(const Vector3D& /* localPosition */, const Vector3D& globalPosition, + const VolumeID& vID) const { + + CellID cID = vID; + + // The volume ID comes with "row" field information (number of sequences) that would screw up the topo-clustering using cell + // neighbours map produced with RecFCCeeCalorimeter/src/components/CreateFCCeeCaloNeighbours.cpp, + // therefore, lets set it to zero, as it is for the cell IDs in the neighbours map. + _decoder->set(cID, "row", 0); + + // For endcap, the volume ID comes with "type" field information which would screw up the topo-clustering as the "row" field, + // therefore, lets set it to zero, as it is for the cell IDs in the neighbours map. + if(m_offsetZ.size() > 1) _decoder->set(cID, "type", 0); + + double lTheta = thetaFromXYZ(globalPosition); + double lPhi = phiFromXYZ(globalPosition); + uint layer = _decoder->get(vID,"layer"); + + // define cell boundaries in R-z plan + if(m_radii.empty()) defineCellsInRZplan(); + + // check if the cells are defined for the given layer + if(m_thetaBins[layer].empty()) dd4hep::printout(dd4hep::ERROR, "FCCSWHCalPhiTheta_k4geo","No cells are defined for layer %d", layer); + + // find the cell (theta bin) corresponding to the hit and return the cellID + for(auto bin : m_thetaBins[layer]) + { + double posz = globalPosition.z(); + if(posz > m_cellEdges[layer][bin].first && posz < m_cellEdges[layer][bin].second) + { + _decoder->set(cID, m_thetaID, bin); + _decoder->set(cID, m_phiID, positionToBin(lPhi, 2 * M_PI / (double)m_phiBins, m_offsetPhi)); + return cID; + } + } + + dd4hep::printout(dd4hep::WARNING, "FCCSWHCalPhiTheta_k4geo","The hit is outside the defined range of the layer %d", layer); + + _decoder->set(cID, m_thetaID, positionToBin(lTheta, m_gridSizeTheta, m_offsetTheta)); + _decoder->set(cID, m_phiID, positionToBin(lPhi, 2 * M_PI / (double)m_phiBins, m_offsetPhi)); + return cID; +} + +/// determine the azimuthal angle phi based on the cell ID +double FCCSWHCalPhiTheta_k4geo::phi(const CellID& cID) const { + CellID phiValue = _decoder->get(cID, m_phiID); + return binToPosition(phiValue, 2. * M_PI / (double)m_phiBins, m_offsetPhi); +} +} +} diff --git a/detectorSegmentations/src/plugins/SegmentationFactories.cpp b/detectorSegmentations/src/plugins/SegmentationFactories.cpp index d362660de..9b0632c55 100644 --- a/detectorSegmentations/src/plugins/SegmentationFactories.cpp +++ b/detectorSegmentations/src/plugins/SegmentationFactories.cpp @@ -34,3 +34,6 @@ DECLARE_SEGMENTATION(GridDRcalo_k4geo, create_segmentation) + +#include "detectorSegmentations/FCCSWHCalPhiTheta_k4geo.h" +DECLARE_SEGMENTATION(FCCSWHCalPhiTheta_k4geo, create_segmentation) From 3b008db7dca9262f5ada2792d0f148388582d291 Mon Sep 17 00:00:00 2001 From: Archil Durglishvili Date: Thu, 7 Nov 2024 12:01:00 +0100 Subject: [PATCH 02/17] Additional PhiRow segmentation for HCal --- .../ALLEGRO_o1_v03/HCalBarrel_TileCal_v03.xml | 13 + .../HCalEndcaps_ThreeParts_TileCal_v03.xml | 13 + .../FCCSWHCalPhiRowHandle_k4geo.h | 113 +++ .../FCCSWHCalPhiRow_k4geo.h | 240 +++++++ .../src/FCCSWHCalPhiRowHandle_k4geo.cpp | 4 + .../src/FCCSWHCalPhiRow_k4geo.cpp | 674 ++++++++++++++++++ .../src/plugins/SegmentationFactories.cpp | 3 + 7 files changed, 1060 insertions(+) create mode 100644 detectorSegmentations/include/detectorSegmentations/FCCSWHCalPhiRowHandle_k4geo.h create mode 100644 detectorSegmentations/include/detectorSegmentations/FCCSWHCalPhiRow_k4geo.h create mode 100644 detectorSegmentations/src/FCCSWHCalPhiRowHandle_k4geo.cpp create mode 100644 detectorSegmentations/src/FCCSWHCalPhiRow_k4geo.cpp diff --git a/FCCee/ALLEGRO/compact/ALLEGRO_o1_v03/HCalBarrel_TileCal_v03.xml b/FCCee/ALLEGRO/compact/ALLEGRO_o1_v03/HCalBarrel_TileCal_v03.xml index cfcc37eac..a62dab702 100644 --- a/FCCee/ALLEGRO/compact/ALLEGRO_o1_v03/HCalBarrel_TileCal_v03.xml +++ b/FCCee/ALLEGRO/compact/ALLEGRO_o1_v03/HCalBarrel_TileCal_v03.xml @@ -57,6 +57,19 @@ offset_phi="-pi+(pi/BarHCal_n_phi_modules)"/> system:4,layer:5,row:9,theta:9,phi:10 + + + system:4,layer:5,row:9,phi:10 + diff --git a/FCCee/ALLEGRO/compact/ALLEGRO_o1_v03/HCalEndcaps_ThreeParts_TileCal_v03.xml b/FCCee/ALLEGRO/compact/ALLEGRO_o1_v03/HCalEndcaps_ThreeParts_TileCal_v03.xml index 5d64a9b5a..c630263ca 100644 --- a/FCCee/ALLEGRO/compact/ALLEGRO_o1_v03/HCalEndcaps_ThreeParts_TileCal_v03.xml +++ b/FCCee/ALLEGRO/compact/ALLEGRO_o1_v03/HCalEndcaps_ThreeParts_TileCal_v03.xml @@ -56,6 +56,19 @@ offset_phi="-pi+(pi/EndcapHCal_n_phi_modules)"/> system:4,type:3,layer:6,row:11,theta:11,phi:10 + + + system:4,type:3,layer:6,row:-8,phi:10 + returning empty neighbours"); + return cellNeighbours; + } + + // part1 min and max layer index + minLayerIdEndcap[0] = minMaxLayerId[0].first; + maxLayerIdEndcap[0] = minMaxLayerId[0].second; + // part2 min and max layer index + minLayerIdEndcap[1] = minMaxLayerId[1].first; + maxLayerIdEndcap[1] = minMaxLayerId[1].second; + // part3 min and max layer index + minLayerIdEndcap[2] = minMaxLayerId[2].first; + maxLayerIdEndcap[2] = minMaxLayerId[2].second; + + // Part 1 + if(currentLayerId >= minLayerIdEndcap[0] && currentLayerId <= maxLayerIdEndcap[0]) + { + minLayerId = minLayerIdEndcap[0]; + maxLayerId = maxLayerIdEndcap[0]; + EndcapPart1 = true; + } + // Part 2 + if(currentLayerId >= minLayerIdEndcap[1] && currentLayerId <= maxLayerIdEndcap[1]) + { + minLayerId = minLayerIdEndcap[1]; + maxLayerId = maxLayerIdEndcap[1]; + EndcapPart2 = true; + } + // Part 3 + if(currentLayerId >= minLayerIdEndcap[2] && currentLayerId <= maxLayerIdEndcap[2]) + { + minLayerId = minLayerIdEndcap[2]; + maxLayerId = maxLayerIdEndcap[2]; + EndcapPart3 = true; + } + + // correct the min and max CellId for endcap + if(currentCellId < 0) // negative-z part (cell index is negative (-1, -2, ... -N)) + { + // second half of elements in m_cellIndexes[currentLayerId] vector corresponds to the negative-z layer cells + minCellId = m_cellIndexes[currentLayerId].back(); + maxCellId = m_cellIndexes[currentLayerId][m_cellIndexes[currentLayerId].size()/2]; + } + else // positive-z part (cell index is positive (1, 2, ... N)) + { + // first half of elements in m_cellIndexes[currentLayerId] vector corresponds to the positive-z layer cells + minCellId = m_cellIndexes[currentLayerId].front(); + maxCellId = m_cellIndexes[currentLayerId][m_cellIndexes[currentLayerId].size()/2 - 1]; + } + } + else // for Barrel + { + minLayerId = 0; + maxLayerId = m_radii.size()-1; + } + //-------------------------------- + + + //-------------------------------- + // Find neighbours + //-------------------------------- + + // if this is not the first layer then look for neighbours in the previous layer + if(currentLayerId > minLayerId) + { + CellID nID = cID ; + int prevLayerId = currentLayerId - 1; + _decoder->set(nID,"layer",prevLayerId); + + // if the granularity is the same for the previous layer then take the cells with currentCellId, currentCellId - 1, and currentCellId + 1 + if(m_gridSizeRow[prevLayerId] == m_gridSizeRow[currentLayerId]) + { + _decoder->set(nID,m_rowID,currentCellId); + cellNeighbours.push_back(nID); // add the cell from the previous layer of the same phi module + if(currentCellId > minCellId) + { + _decoder->set(nID,m_rowID,currentCellId - 1); + cellNeighbours.push_back(nID); // add the cell from the previous layer of the same phi module + } + if(currentCellId < maxCellId) + { + _decoder->set(nID,m_rowID,currentCellId + 1); + cellNeighbours.push_back(nID); // add the cell from the previous layer of the same phi module + } + } + // if the granularity is different + else + { + // determine the cell index in the previous layer that is below of the current cell + int idx = (currentCellId > 0) ? ((currentCellId-1)/m_gridSizeRow[prevLayerId] + 1) : ((currentCellId+1)/m_gridSizeRow[prevLayerId] - 1); + _decoder->set(nID,m_rowID,idx); + cellNeighbours.push_back(nID); // add the cell from the previous layer of the same phi module + + // + if((m_gridSizeRow[prevLayerId] - abs(currentCellId)%m_gridSizeRow[prevLayerId]) == (m_gridSizeRow[prevLayerId]-1) && currentCellId > minCellId) + { + _decoder->set(nID,m_rowID, (idx > 0) ? (idx - 1) : (idx + 1)); + cellNeighbours.push_back(nID); // add the cell from the previous layer of the same phi module + } + + // + if(abs(currentCellId)%m_gridSizeRow[prevLayerId] == 0 && currentCellId < maxCellId) + { + _decoder->set(nID,m_rowID, (idx > 0) ? (idx + 1) : (idx - 1)); + cellNeighbours.push_back(nID); // add the cell from the previous layer of the same phi module + } + } + } + + // if this is not the last layer then look for neighbours in the next layer + if(currentLayerId < maxLayerId) + { + CellID nID = cID ; + int nextLayerId = currentLayerId + 1; + _decoder->set(nID,"layer",nextLayerId); + + // if the granularity is the same for the next layer then take the cells with currentCellId, currentCellId - 1, and currentCellId + 1 + if(m_gridSizeRow[nextLayerId] == m_gridSizeRow[currentLayerId]) + { + _decoder->set(nID,m_rowID,currentCellId); + cellNeighbours.push_back(nID); // add the cell from the next layer of the same phi module + if(currentCellId > minCellId) + { + _decoder->set(nID,m_rowID,currentCellId - 1); + cellNeighbours.push_back(nID); // add the cell from the next layer of the same phi module + } + if(currentCellId < maxCellId) + { + _decoder->set(nID,m_rowID,currentCellId + 1); + cellNeighbours.push_back(nID); // add the cell from the next layer of the same phi module + } + } + // if the granularity is different + else + { + // determine the cell index in the next layer that is below of the current cell + int idx = (currentCellId > 0) ? ((currentCellId-1)/m_gridSizeRow[nextLayerId] + 1) : ((currentCellId+1)/m_gridSizeRow[nextLayerId] - 1); + _decoder->set(nID,m_rowID,idx); + cellNeighbours.push_back(nID); // add the cell from the next layer of the same phi module + + // + if((m_gridSizeRow[nextLayerId] - abs(currentCellId)%m_gridSizeRow[nextLayerId]) == (m_gridSizeRow[nextLayerId]-1) && currentCellId > minCellId) + { + _decoder->set(nID,m_rowID, (idx > 0) ? (idx - 1) : (idx + 1)); + cellNeighbours.push_back(nID); // add the cell from the next layer of the same phi module + } + + // + if(abs(currentCellId)%m_gridSizeRow[nextLayerId] == 0 && currentCellId < maxCellId) + { + _decoder->set(nID,m_rowID, (idx > 0) ? (idx + 1) : (idx - 1)); + cellNeighbours.push_back(nID); // add the cell from the next layer of the same phi module + } + } + } + + // if this is not the first cell in the given layer then add the previous cell + if(currentCellId > minCellId) + { + CellID nID = cID ; + _decoder->set(nID,m_rowID,currentCellId - 1); + cellNeighbours.push_back(nID); // add the previous cell from current layer of the same phi module + } + // if this is not the last cell in the given layer then add the next cell + if(currentCellId < maxCellId) + { + CellID nID = cID ; + _decoder->set(nID,m_rowID,currentCellId + 1); + cellNeighbours.push_back(nID); // add the next cell from current layer of the same phi module + } + + // if this is the Endcap then look for neighbours in different parts as well + if(m_offsetZ.size() > 1) + { + double currentLayerRmin = m_radii[currentLayerId] - 0.5*m_layerDepth[currentLayerId]; + double currentLayerRmax = m_radii[currentLayerId] + 0.5*m_layerDepth[currentLayerId]; + + // if the cell is in negative-z part, then swap min and max cell indexes + if(currentCellId < 0) + { + minCellId = m_cellIndexes[currentLayerId][m_cellIndexes[currentLayerId].size()/2]; // this should be -1 + maxCellId = m_cellIndexes[currentLayerId].back(); // this should be -N + } + + // if it is the last cell in the part1 + if(EndcapPart1 && currentCellId == maxCellId ) + { + // find the layers in the part2 that share a border with the current layer + for(int part2layerId = minLayerIdEndcap[1]; part2layerId <= maxLayerIdEndcap[1]; part2layerId++) + { + double Rmin = m_radii[part2layerId] - 0.5*m_layerDepth[part2layerId]; + double Rmax = m_radii[part2layerId] + 0.5*m_layerDepth[part2layerId]; + + if( (Rmin >= currentLayerRmin && Rmin <= currentLayerRmax) + || (Rmax >= currentLayerRmin && Rmax <= currentLayerRmax) + || (currentLayerRmin >= Rmin && currentLayerRmin <= Rmax) + || (currentLayerRmax >= Rmin && currentLayerRmax <= Rmax) + ) + { + CellID nID = cID ; + _decoder->set(nID,"layer",part2layerId); + _decoder->set(nID,m_rowID,minCellId); + cellNeighbours.push_back(nID); // add the first cell from part2 layer + } + } + } + + // if it is the first cell in the part2 + if(EndcapPart2 && currentCellId == minCellId) + { + // find the layers in part1 that share a border with the current layer + for(int part1layerId = minLayerIdEndcap[0]; part1layerId <= maxLayerIdEndcap[0]; part1layerId++) + { + double Rmin = m_radii[part1layerId] - 0.5*m_layerDepth[part1layerId]; + double Rmax = m_radii[part1layerId] + 0.5*m_layerDepth[part1layerId]; + + if( (Rmin >= currentLayerRmin && Rmin <= currentLayerRmax) + || (Rmax >= currentLayerRmin && Rmax <= currentLayerRmax) + || (currentLayerRmin >= Rmin && currentLayerRmin <= Rmax) + || (currentLayerRmax >= Rmin && currentLayerRmax <= Rmax) + ) + { + CellID nID = cID ; + _decoder->set(nID,"layer",part1layerId); + _decoder->set(nID,m_rowID,maxCellId); + cellNeighbours.push_back(nID); // add the last cell from the part1 layer + } + } + } + + // if it is the last cell in the part2 + if(EndcapPart2 && currentCellId == maxCellId) + { + // find the layers in the part3 that share a border with the current layer + for(int part3layerId = minLayerIdEndcap[2]; part3layerId <= maxLayerIdEndcap[2]; part3layerId++) + { + double Rmin = m_radii[part3layerId] - 0.5*m_layerDepth[part3layerId]; + double Rmax = m_radii[part3layerId] + 0.5*m_layerDepth[part3layerId]; + + if( (Rmin >= currentLayerRmin && Rmin <= currentLayerRmax) + || (Rmax >= currentLayerRmin && Rmax <= currentLayerRmax) + || (currentLayerRmin >= Rmin && currentLayerRmin <= Rmax) + || (currentLayerRmax >= Rmin && currentLayerRmax <= Rmax) + ) + { + CellID nID = cID ; + _decoder->set(nID,"layer",part3layerId); + _decoder->set(nID,m_rowID,minCellId); + cellNeighbours.push_back(nID); // add the first cell from the part3 layer + } + } + } + + // if it is the first cell in the part3 + if(EndcapPart3 && currentCellId == minCellId) + { + // find the layers in the part2 that share a border with the current layer + for(int part2layerId = minLayerIdEndcap[1]; part2layerId <= maxLayerIdEndcap[1]; part2layerId++) + { + double Rmin = m_radii[part2layerId] - 0.5*m_layerDepth[part2layerId]; + double Rmax = m_radii[part2layerId] + 0.5*m_layerDepth[part2layerId]; + + if( (Rmin >= currentLayerRmin && Rmin <= currentLayerRmax) + || (Rmax >= currentLayerRmin && Rmax <= currentLayerRmax) + || (currentLayerRmin >= Rmin && currentLayerRmin <= Rmax) + || (currentLayerRmax >= Rmin && currentLayerRmax <= Rmax) + ) + { + CellID nID = cID ; + _decoder->set(nID,"layer",part2layerId); + _decoder->set(nID,m_rowID,maxCellId); + cellNeighbours.push_back(nID); // add the last cell from the part2 layer + } + } + } + } + + // Now loop over the neighbours and add the cells from next/previous phi module + for(auto nID : cellNeighbours) + { + CellID newID = nID; + // previous: if the current is 0 then previous is the last bin (id = m_phiBins - 1) else current - 1 + _decoder->set(newID,m_phiID, (_decoder->get(nID,m_phiID) == 0) ? m_phiBins - 1 : _decoder->get(nID,m_phiID) - 1); + cellNeighbours.push_back(newID); + // next: if the current is the last bin (id = m_phiBins - 1) then the next is the first bin (id = 0) else current + 1 + _decoder->set(newID,m_phiID, (_decoder->get(nID,m_phiID) == (m_phiBins - 1)) ? 0 : _decoder->get(nID,m_phiID) + 1); + cellNeighbours.push_back(newID); + } + + // At the end, find neighbours with the same layer/row in next/previous phi module + CellID nID = cID ; + // previous: if the current is 0 then previous is the last bin (id = m_phiBins - 1) else current - 1 + _decoder->set(nID,m_phiID, (_decoder->get(cID,m_phiID) == 0) ? m_phiBins - 1 : _decoder->get(cID,m_phiID) - 1); + cellNeighbours.push_back(nID); + // next: if the current is the last bin (id = m_phiBins - 1) then the next is the first bin (id = 0) else current + 1 + _decoder->set(nID,m_phiID, (_decoder->get(cID,m_phiID) == (m_phiBins - 1)) ? 0 : _decoder->get(cID,m_phiID) + 1); + cellNeighbours.push_back(nID); + + return cellNeighbours; +} + +/// Determine minimum and maximum polar angle of the cell +std::array FCCSWHCalPhiRow_k4geo::cellTheta(const CellID& cID) const { + std::array cTheta = {M_PI,M_PI}; + + // get the cell index + int idx = _decoder->get(cID, m_rowID); + // get the layer index + uint layer = _decoder->get(cID,"layer"); + + if(m_radii.empty()) calculateLayerRadii(); + if(m_cellEdges.empty()) return cTheta; + + double zlow = m_cellEdges[layer][idx].first; + double zhigh = m_cellEdges[layer][idx].second; + + double Rmin = m_radii[layer] - 0.5*m_layerDepth[layer]; + double Rmax = m_radii[layer] + 0.5*m_layerDepth[layer]; + + if( theta(cID) < M_PI/2. ) + { + cTheta[0] = std::atan2(Rmin,zhigh); // theta min + cTheta[1] = std::atan2(Rmax,zlow); // theta max + } + else + { + cTheta[0] = std::atan2(Rmax,zhigh); // theta min + cTheta[1] = std::atan2(Rmin,zlow); // theta max + } + + return cTheta; +} + +} +} diff --git a/detectorSegmentations/src/plugins/SegmentationFactories.cpp b/detectorSegmentations/src/plugins/SegmentationFactories.cpp index 9b0632c55..a3ac9b0e9 100644 --- a/detectorSegmentations/src/plugins/SegmentationFactories.cpp +++ b/detectorSegmentations/src/plugins/SegmentationFactories.cpp @@ -37,3 +37,6 @@ DECLARE_SEGMENTATION(FCCSWEndcapTurbine_k4geo, create_segmentation) + +#include "detectorSegmentations/FCCSWHCalPhiRow_k4geo.h" +DECLARE_SEGMENTATION(FCCSWHCalPhiRow_k4geo, create_segmentation) From 3a79ed0fb621f4e19327d226d8e4c81f7b62a957 Mon Sep 17 00:00:00 2001 From: Archil Durglishvili Date: Thu, 7 Nov 2024 12:05:37 +0100 Subject: [PATCH 03/17] Updates for HCal PhiTheta segmentation: determine neighbours --- .../FCCSWHCalPhiThetaHandle_k4geo.h | 6 +- .../FCCSWHCalPhiTheta_k4geo.h | 22 +- .../src/FCCSWHCalPhiTheta_k4geo.cpp | 453 ++++++++++++++++++ 3 files changed, 476 insertions(+), 5 deletions(-) diff --git a/detectorSegmentations/include/detectorSegmentations/FCCSWHCalPhiThetaHandle_k4geo.h b/detectorSegmentations/include/detectorSegmentations/FCCSWHCalPhiThetaHandle_k4geo.h index e787bbd45..475b39649 100644 --- a/detectorSegmentations/include/detectorSegmentations/FCCSWHCalPhiThetaHandle_k4geo.h +++ b/detectorSegmentations/include/detectorSegmentations/FCCSWHCalPhiThetaHandle_k4geo.h @@ -18,7 +18,7 @@ class Segmentation; template class SegmentationWrapper; /// We need some abbreviation to make the code more readable. -typedef Handle> FCCSWGridPhiThetaHandle_k4geo; +typedef Handle> FCCSWHCalPhiThetaHandle_k4geo; /// Implementation class for the HCal phi-theta segmentation. /** @@ -38,10 +38,10 @@ typedef Handle> FCC * conveniance reasons instantiated in DD4hep/src/Segmentations.cpp. * */ -class FCCSWHCalPhiTheta_k4geo : public FCCSWGridPhiThetaHandle_k4geo { +class FCCSWHCalPhiTheta_k4geo : public FCCSWHCalPhiThetaHandle_k4geo { public: /// Defintion of the basic handled object - typedef FCCSWGridPhiThetaHandle_k4geo::Object Object; + typedef FCCSWHCalPhiThetaHandle_k4geo::Object Object; public: /// Default constructor diff --git a/detectorSegmentations/include/detectorSegmentations/FCCSWHCalPhiTheta_k4geo.h b/detectorSegmentations/include/detectorSegmentations/FCCSWHCalPhiTheta_k4geo.h index 439d56ab7..58338b57a 100644 --- a/detectorSegmentations/include/detectorSegmentations/FCCSWHCalPhiTheta_k4geo.h +++ b/detectorSegmentations/include/detectorSegmentations/FCCSWHCalPhiTheta_k4geo.h @@ -26,9 +26,8 @@ namespace dd4hep { virtual ~FCCSWHCalPhiTheta_k4geo() = default; /** Determine the position of HCal cell based on the cellID. - * @warning This segmentation has no knowledge of radius, so radius = 1 is taken into calculations. * @param[in] aCellId ID of a cell. - * return Position (radius = 1). + * return Position. */ virtual Vector3D position(const CellID& aCellID) const; @@ -41,6 +40,12 @@ namespace dd4hep { virtual CellID cellID(const Vector3D& aLocalPosition, const Vector3D& aGlobalPosition, const VolumeID& aVolumeID) const; + /** Find neighbours of the cell + * @param[in] aCellId ID of a cell. + * return vector of neighbour cellIDs. + */ + std::vector neighbours(const CellID& cID) const; + /** Calculate layer radii and edges in z-axis, then define cell edges in each layer using defineCellEdges(). */ void defineCellsInRZplan() const; @@ -119,6 +124,17 @@ namespace dd4hep { inline Vector3D centerPosition(const CellID& cID) const; + /** Determine the minimum and maximum polar angle of HCal cell based on the cellID. + * @param[in] aCellId ID of a cell. + * return Theta. + */ + std::array cellTheta(const CellID& cID) const; + + /** Get the min and max layer indexes of each HCal part. + * For Endcap, returns the three elements vector, while for Barrel - single element vector. + */ + std::vector > getMinMaxLayerId() const ; + /** Set the number of bins in azimuthal angle. * @param[in] aNumberBins Number of bins in phi. */ @@ -183,6 +199,8 @@ namespace dd4hep { mutable std::vector m_radii; /// z-min and z-max of each layer mutable std::vector > m_layerEdges; + /// dR of each layer + mutable std::vector m_layerDepth; /// theta bins (cells) in each layer mutable std::vector > m_thetaBins; /// z-min and z-max of each cell (theta bin) in each layer diff --git a/detectorSegmentations/src/FCCSWHCalPhiTheta_k4geo.cpp b/detectorSegmentations/src/FCCSWHCalPhiTheta_k4geo.cpp index cc49d9e32..765b9b3fc 100644 --- a/detectorSegmentations/src/FCCSWHCalPhiTheta_k4geo.cpp +++ b/detectorSegmentations/src/FCCSWHCalPhiTheta_k4geo.cpp @@ -94,6 +94,7 @@ void FCCSWHCalPhiTheta_k4geo::defineCellsInRZplan() const { m_radii.push_back(m_offsetR[0] + moduleDepth - m_dRlayer[i_dR]*0.5); // layer lower and upper edges in z-axis m_layerEdges.push_back( std::make_pair(m_offsetZ[0] - 0.5*m_widthZ[0], m_offsetZ[0] + 0.5*m_widthZ[0]) ); + m_layerDepth.push_back(m_dRlayer[i_dR]); } } } // Barrel @@ -114,6 +115,7 @@ void FCCSWHCalPhiTheta_k4geo::defineCellsInRZplan() const { m_radii.push_back(m_offsetR[0] + moduleDepth1 - m_dRlayer[i_dR]*0.5); // layer lower and upper edges in z-axis m_layerEdges.push_back( std::make_pair(m_offsetZ[0] - 0.5*m_widthZ[0], m_offsetZ[0] + 0.5*m_widthZ[0]) ); + m_layerDepth.push_back(m_dRlayer[i_dR]); } } // part2 @@ -125,6 +127,7 @@ void FCCSWHCalPhiTheta_k4geo::defineCellsInRZplan() const { m_radii.push_back(m_offsetR[1] + moduleDepth2 - m_dRlayer[i_dR]*0.5); // layer lower and upper edges in z-axis m_layerEdges.push_back( std::make_pair(m_offsetZ[1] - 0.5*m_widthZ[1], m_offsetZ[1] + 0.5*m_widthZ[1]) ); + m_layerDepth.push_back(m_dRlayer[i_dR]); } } // part3 @@ -136,6 +139,7 @@ void FCCSWHCalPhiTheta_k4geo::defineCellsInRZplan() const { m_radii.push_back(m_offsetR[2] + moduleDepth3 - m_dRlayer[i_dR]*0.5); // layer lower and upper edges in z-axis m_layerEdges.push_back( std::make_pair(m_offsetZ[2] - 0.5*m_widthZ[2], m_offsetZ[2] + 0.5*m_widthZ[2]) ); + m_layerDepth.push_back(m_dRlayer[i_dR]); } } } // ThreePartsEndcap @@ -280,5 +284,454 @@ double FCCSWHCalPhiTheta_k4geo::phi(const CellID& cID) const { CellID phiValue = _decoder->get(cID, m_phiID); return binToPosition(phiValue, 2. * M_PI / (double)m_phiBins, m_offsetPhi); } + + +/// Get the min and max layer indexes for each part of the HCal +std::vector > FCCSWHCalPhiTheta_k4geo::getMinMaxLayerId() const { + /* + * hardcoded numbers would be the following: + * std::vector minLayerIdEndcap = {0, 6, 15}; + * std::vector maxLayerIdEndcap = {5, 14, 36}; + * but lets try to avoid hardcoding: + */ + + std::vector > minMaxLayerId; + + if(m_radii.empty()) defineCellsInRZplan(); + if(m_radii.empty()) return minMaxLayerId; + + + std::vector minLayerIdEndcap = {0, 0, 0}; + std::vector maxLayerIdEndcap = {-1, 0, 0}; + if(m_offsetZ.size() > 1) + { + uint Nl = m_numLayers.size()/3; + + // count the number of layers in the first part of the Endcap + for(uint i=0; i < Nl; i++) maxLayerIdEndcap[0] += m_numLayers[i]; + + minLayerIdEndcap[1] = maxLayerIdEndcap[0]+1; + maxLayerIdEndcap[1] = maxLayerIdEndcap[0]; + // count the number of layers in the second part of the Endcap + for(uint i=0; i < Nl; i++) maxLayerIdEndcap[1] += m_numLayers[i+Nl]; + + minLayerIdEndcap[2] = maxLayerIdEndcap[1]+1; + maxLayerIdEndcap[2] = maxLayerIdEndcap[1]; + // count the number of layers in the third part of the Endcap + for(uint i=0; i < Nl; i++) maxLayerIdEndcap[2] += m_numLayers[i+2*Nl]; + + minMaxLayerId.push_back(std::make_pair(minLayerIdEndcap[0],maxLayerIdEndcap[0])); + minMaxLayerId.push_back(std::make_pair(minLayerIdEndcap[1],maxLayerIdEndcap[1])); + minMaxLayerId.push_back(std::make_pair(minLayerIdEndcap[2],maxLayerIdEndcap[2])); + } + else // for Barrel + { + minMaxLayerId.push_back(std::make_pair(0,m_radii.size()-1)); + } + + return minMaxLayerId; +} + + +/// Calculates the neighbours of the given cell ID and adds them to the list of neighbours +std::vector FCCSWHCalPhiTheta_k4geo::neighbours(const CellID& cID) const { + std::vector cellNeighbours; + + if(m_radii.empty()) defineCellsInRZplan(); + if(m_thetaBins.empty()) return cellNeighbours; + + bool EndcapPart1 = false; + bool EndcapPart2 = false; + bool EndcapPart3 = false; + + int minLayerId = -1; + int maxLayerId = -1; + + int currentLayerId = _decoder->get(cID,"layer"); + int currentCellThetaBin = _decoder->get(cID,m_thetaID); + + int minCellThetaBin = m_thetaBins[currentLayerId].front(); + int maxCellThetaBin = m_thetaBins[currentLayerId].back(); + + //-------------------------------- + // Determine min and max layer Id + //-------------------------------- + // if this is the segmentation of three parts Endcap + std::vector minLayerIdEndcap = {0, 0, 0}; + std::vector maxLayerIdEndcap = {0, 0, 0}; + if(m_offsetZ.size() > 1) + { + std::vector > minMaxLayerId(getMinMaxLayerId()); + if(minMaxLayerId.empty()) + { + dd4hep::printout(dd4hep::ERROR, "FCCSWHCalPhiTheta_k4geo","Can not get ThreePartsEndcap min and max layer indexes! --> returning empty neighbours"); + return cellNeighbours; + } + + // part1 min and max layer index + minLayerIdEndcap[0] = minMaxLayerId[0].first; + maxLayerIdEndcap[0] = minMaxLayerId[0].second; + // part2 min and max layer index + minLayerIdEndcap[1] = minMaxLayerId[1].first; + maxLayerIdEndcap[1] = minMaxLayerId[1].second; + // part3 min and max layer index + minLayerIdEndcap[2] = minMaxLayerId[2].first; + maxLayerIdEndcap[2] = minMaxLayerId[2].second; + + // Part 1 + if(currentLayerId >= minLayerIdEndcap[0] && currentLayerId <= maxLayerIdEndcap[0]) + { + minLayerId = minLayerIdEndcap[0]; + maxLayerId = maxLayerIdEndcap[0]; + EndcapPart1 = true; + } + // Part 2 + if(currentLayerId >= minLayerIdEndcap[1] && currentLayerId <= maxLayerIdEndcap[1]) + { + minLayerId = minLayerIdEndcap[1]; + maxLayerId = maxLayerIdEndcap[1]; + EndcapPart2 = true; + } + // Part 3 + if(currentLayerId >= minLayerIdEndcap[2] && currentLayerId <= maxLayerIdEndcap[2]) + { + minLayerId = minLayerIdEndcap[2]; + maxLayerId = maxLayerIdEndcap[2]; + EndcapPart3 = true; + } + + // correct the min and max theta bin for endcap + if(theta(cID) > M_PI/2) // negative-z part + { + // second half of elements in m_thetaBins[currentLayerId] vector corresponds to the negative-z layer cells + minCellThetaBin = m_thetaBins[currentLayerId][m_thetaBins[currentLayerId].size()/2]; + maxCellThetaBin = m_thetaBins[currentLayerId].back(); + } + else // positive-z part + { + // first half of elements in m_thetaBins[currentLayerId] vector corresponds to the positive-z layer cells + minCellThetaBin = m_thetaBins[currentLayerId].front(); + maxCellThetaBin = m_thetaBins[currentLayerId][m_thetaBins[currentLayerId].size()/2 - 1]; + } + } + else // for Barrel + { + minLayerId = 0; + maxLayerId = m_radii.size()-1; + } + //-------------------------------- + + + //-------------------------------- + // Find neighbours + //-------------------------------- + + //--------------------------------------------- + // this part is same for both Barrel and Endcap + //--------------------------------------------- + // if this is not the first cell in the given layer then add the previous cell + if(currentCellThetaBin > minCellThetaBin) + { + CellID nID = cID ; + _decoder->set(nID,m_thetaID,currentCellThetaBin - 1); + cellNeighbours.push_back(nID); // add the previous cell from current layer of the same phi module + } + // if this is not the last cell in the given layer then add the next cell + if(currentCellThetaBin < maxCellThetaBin) + { + CellID nID = cID ; + _decoder->set(nID,m_thetaID,currentCellThetaBin + 1); + cellNeighbours.push_back(nID); // add the next cell from current layer of the same phi module + } + //---------------------------------------------- + + // deal with the Barrel + if(m_offsetZ.size() == 1) + { + // if this is not the first layer then look for neighbours in the previous layer + if(currentLayerId > minLayerId) + { + CellID nID = cID ; + int prevLayerId = currentLayerId - 1; + _decoder->set(nID,"layer",prevLayerId); + + _decoder->set(nID,m_thetaID,currentCellThetaBin); + cellNeighbours.push_back(nID); // add the cell with the same theta bin from the previous layer of the same phi module + + // if the cID is in the positive-z side and prev layer cell is not in the first theta bin then add the cell from previous theta bin + if(theta(cID) < M_PI/2. && currentCellThetaBin > m_thetaBins[prevLayerId].front() ) + { + _decoder->set(nID,m_thetaID,currentCellThetaBin - 1); + cellNeighbours.push_back(nID); // add the cell from the previous layer of the same phi module + } + // if the cID is in the negative-z side and prev layer cell is not in the last theta bin then add the cell from previous theta bin + if(theta(cID) > M_PI/2. && currentCellThetaBin < m_thetaBins[prevLayerId].back() ) + { + _decoder->set(nID,m_thetaID,currentCellThetaBin + 1); + cellNeighbours.push_back(nID); // add the cell from the previous layer of the same phi module + } + } + + // if this is not the last layer then look for neighbours in the next layer + if(currentLayerId < maxLayerId) + { + double currentCellZmin = m_cellEdges[currentLayerId][currentCellThetaBin].first; + double currentCellZmax = m_cellEdges[currentLayerId][currentCellThetaBin].second; + + CellID nID = cID ; + int nextLayerId = currentLayerId + 1; + _decoder->set(nID,"layer",nextLayerId); + + _decoder->set(nID,m_thetaID,currentCellThetaBin); + cellNeighbours.push_back(nID); // add the cell with the same theta bin from the next layer of the same phi module + + // if the cID is in the positive-z side + if(theta(cID) < M_PI/2.) + { + //add the next layer cell from the next theta bin + _decoder->set(nID,m_thetaID,currentCellThetaBin + 1); + cellNeighbours.push_back(nID); + + //add the next layer cell from the next-to-next theta bin if it overlaps with the current cell in z-coordinate + double zmax = m_cellEdges[nextLayerId][currentCellThetaBin + 2].second; + if(zmax >= currentCellZmin) + { + //add the next layer cell from the next to next theta bin + _decoder->set(nID,m_thetaID,currentCellThetaBin + 2); + cellNeighbours.push_back(nID); + } + } + // if the cID is in the negative-z side + if(theta(cID) > M_PI/2.) + { + //add the next layer cell from the previous theta bin + _decoder->set(nID,m_thetaID,currentCellThetaBin - 1); + cellNeighbours.push_back(nID); + + //add the next layer cell from the next-to-next theta bin if it overlaps with the current cell in z-coordinate + double zmin = m_cellEdges[nextLayerId][currentCellThetaBin - 2].first; + if(zmin <= currentCellZmax) + { + //add the next layer cell from the next to next theta bin + _decoder->set(nID,m_thetaID,currentCellThetaBin - 2); + cellNeighbours.push_back(nID); + } + } + } + } + + // if this is the Endcap then look for neighbours in different parts as well + if(m_offsetZ.size() > 1) + { + double currentCellZmin = m_cellEdges[currentLayerId][currentCellThetaBin].first; + double currentCellZmax = m_cellEdges[currentLayerId][currentCellThetaBin].second; + + // if this is not the first layer then look for neighbours in the previous layer + if(currentLayerId > minLayerId) + { + CellID nID = cID ; + int prevLayerId = currentLayerId - 1; + _decoder->set(nID,"layer",prevLayerId); + // find the ones that share at least part of a border with the current cell + for( auto bin : m_thetaBins[prevLayerId] ) + { + double zmin = m_cellEdges[prevLayerId][bin].first; + double zmax = m_cellEdges[prevLayerId][bin].second; + + if( (zmin >= currentCellZmin && zmin <=currentCellZmax) + || (zmax >= currentCellZmin && zmax <=currentCellZmax) + || (currentCellZmin >= zmin && currentCellZmax <= zmax) + ) + { + _decoder->set(nID,m_thetaID,bin); + cellNeighbours.push_back(nID); // add the cell from the previous layer of the same phi module + } + } + } + // if this is not the last layer then look for neighbours in the next layer + if(currentLayerId < maxLayerId) + { + CellID nID = cID ; + int nextLayerId = currentLayerId + 1; + _decoder->set(nID,"layer",nextLayerId); + // find the ones that share at least part of a border with the current cell + for( auto bin : m_thetaBins[nextLayerId] ) + { + double zmin = m_cellEdges[nextLayerId][bin].first; + double zmax = m_cellEdges[nextLayerId][bin].second; + + if( (zmin >= currentCellZmin && zmin <=currentCellZmax) + || (zmax >= currentCellZmin && zmax <=currentCellZmax) + || (currentCellZmin >= zmin && currentCellZmax <= zmax) + ) + { + _decoder->set(nID,m_thetaID,bin); + cellNeighbours.push_back(nID); // add the cell from the previous layer of the same phi module + } + } + } + + + // + double currentLayerRmin = m_radii[currentLayerId] - 0.5*m_layerDepth[currentLayerId]; + double currentLayerRmax = m_radii[currentLayerId] + 0.5*m_layerDepth[currentLayerId]; + + + // if the cell is in negative-z part, then swap min and max theta bins + if(theta(cID) > M_PI/2.) + { + minCellThetaBin = m_thetaBins[currentLayerId].back(); + maxCellThetaBin = m_thetaBins[currentLayerId][m_thetaBins[currentLayerId].size()/2]; + } + + // if it is the last cell in the part1 + if(EndcapPart1 && currentCellThetaBin == minCellThetaBin ) + { + // find the layers in the part2 that share a border with the current layer + for(int part2layerId = minLayerIdEndcap[1]; part2layerId <= maxLayerIdEndcap[1]; part2layerId++) + { + double Rmin = m_radii[part2layerId] - 0.5*m_layerDepth[part2layerId]; + double Rmax = m_radii[part2layerId] + 0.5*m_layerDepth[part2layerId]; + + if( (Rmin >= currentLayerRmin && Rmin <= currentLayerRmax) + || (Rmax >= currentLayerRmin && Rmax <= currentLayerRmax) + || (currentLayerRmin >= Rmin && currentLayerRmin <= Rmax) + || (currentLayerRmax >= Rmin && currentLayerRmax <= Rmax) + ) + { + CellID nID = cID ; + _decoder->set(nID,"layer",part2layerId); + _decoder->set(nID,m_thetaID,maxCellThetaBin); + cellNeighbours.push_back(nID); // add the last theta bin cell from part2 layer + } + } + } + + // if it is the last theta bin cell in the part2 + if(EndcapPart2 && currentCellThetaBin == maxCellThetaBin) + { + // find the layers in part1 that share a border with the current layer + for(int part1layerId = minLayerIdEndcap[0]; part1layerId <= maxLayerIdEndcap[0]; part1layerId++) + { + double Rmin = m_radii[part1layerId] - 0.5*m_layerDepth[part1layerId]; + double Rmax = m_radii[part1layerId] + 0.5*m_layerDepth[part1layerId]; + + if( (Rmin >= currentLayerRmin && Rmin <= currentLayerRmax) + || (Rmax >= currentLayerRmin && Rmax <= currentLayerRmax) + || (currentLayerRmin >= Rmin && currentLayerRmin <= Rmax) + || (currentLayerRmax >= Rmin && currentLayerRmax <= Rmax) + ) + { + CellID nID = cID ; + _decoder->set(nID,"layer",part1layerId); + _decoder->set(nID,m_thetaID,minCellThetaBin); + cellNeighbours.push_back(nID); // add the first theta bin cell from the part1 layer + } + } + } + + // if it is the first theta bin cell in the part2 + if(EndcapPart2 && currentCellThetaBin == minCellThetaBin) + { + // find the layers in the part3 that share a border with the current layer + for(int part3layerId = minLayerIdEndcap[2]; part3layerId <= maxLayerIdEndcap[2]; part3layerId++) + { + double Rmin = m_radii[part3layerId] - 0.5*m_layerDepth[part3layerId]; + double Rmax = m_radii[part3layerId] + 0.5*m_layerDepth[part3layerId]; + + if( (Rmin >= currentLayerRmin && Rmin <= currentLayerRmax) + || (Rmax >= currentLayerRmin && Rmax <= currentLayerRmax) + || (currentLayerRmin >= Rmin && currentLayerRmin <= Rmax) + || (currentLayerRmax >= Rmin && currentLayerRmax <= Rmax) + ) + { + CellID nID = cID ; + _decoder->set(nID,"layer",part3layerId); + _decoder->set(nID,m_thetaID,maxCellThetaBin); + cellNeighbours.push_back(nID); // add the first cell from the part3 layer + } + } + } + + // if it is the last theta bin cell in the part3 + if(EndcapPart3 && currentCellThetaBin == maxCellThetaBin) + { + // find the layers in the part2 that share a border with the current layer + for(int part2layerId = minLayerIdEndcap[1]; part2layerId <= maxLayerIdEndcap[1]; part2layerId++) + { + double Rmin = m_radii[part2layerId] - 0.5*m_layerDepth[part2layerId]; + double Rmax = m_radii[part2layerId] + 0.5*m_layerDepth[part2layerId]; + + if( (Rmin >= currentLayerRmin && Rmin <= currentLayerRmax) + || (Rmax >= currentLayerRmin && Rmax <= currentLayerRmax) + || (currentLayerRmin >= Rmin && currentLayerRmin <= Rmax) + || (currentLayerRmax >= Rmin && currentLayerRmax <= Rmax) + ) + { + CellID nID = cID ; + _decoder->set(nID,"layer",part2layerId); + _decoder->set(nID,m_thetaID,minCellThetaBin); + cellNeighbours.push_back(nID); // add the first theta bin cell from the part2 layer + } + } + } + } + + // Now loop over the neighbours and add the cells from next/previous phi module + for(auto nID : cellNeighbours) + { + CellID newID = nID; + // previous: if the current is 0 then previous is the last bin (id = m_phiBins - 1) else current - 1 + _decoder->set(newID,m_phiID, (_decoder->get(nID,m_phiID) == 0) ? m_phiBins - 1 : _decoder->get(nID,m_phiID) - 1); + cellNeighbours.push_back(newID); + // next: if the current is the last bin (id = m_phiBins - 1) then the next is the first bin (id = 0) else current + 1 + _decoder->set(newID,m_phiID, (_decoder->get(nID,m_phiID) == (m_phiBins - 1)) ? 0 : _decoder->get(nID,m_phiID) + 1); + cellNeighbours.push_back(newID); + } + + // At the end, find neighbours with the same layer/row in next/previous phi module + CellID nID = cID ; + // previous: if the current is 0 then previous is the last bin (id = m_phiBins - 1) else current - 1 + _decoder->set(nID,m_phiID, (_decoder->get(cID,m_phiID) == 0) ? m_phiBins - 1 : _decoder->get(cID,m_phiID) - 1); + cellNeighbours.push_back(nID); + // next: if the current is the last bin (id = m_phiBins - 1) then the next is the first bin (id = 0) else current + 1 + _decoder->set(nID,m_phiID, (_decoder->get(cID,m_phiID) == (m_phiBins - 1)) ? 0 : _decoder->get(cID,m_phiID) + 1); + cellNeighbours.push_back(nID); + + return cellNeighbours; +} + +/// Determine minimum and maximum polar angle of the cell +std::array FCCSWHCalPhiTheta_k4geo::cellTheta(const CellID& cID) const { + std::array cTheta = {M_PI,M_PI}; + + // get the cell index + int idx = _decoder->get(cID, m_thetaID); + // get the layer index + uint layer = _decoder->get(cID,"layer"); + + if(m_radii.empty()) defineCellsInRZplan(); + if(m_cellEdges.empty()) return cTheta; + + double zlow = m_cellEdges[layer][idx].first; + double zhigh = m_cellEdges[layer][idx].second; + + double Rmin = m_radii[layer] - 0.5*m_layerDepth[layer]; + double Rmax = m_radii[layer] + 0.5*m_layerDepth[layer]; + + if( theta(cID) < M_PI/2. ) + { + cTheta[0] = std::atan2(Rmin,zhigh); // theta min + cTheta[1] = std::atan2(Rmax,zlow); // theta max + } + else + { + cTheta[0] = std::atan2(Rmax,zhigh); // theta min + cTheta[1] = std::atan2(Rmin,zlow); // theta max + } + + return cTheta; +} + } } From 6fac8a8d4f87788895c15b1b6909e0e4f5dede56 Mon Sep 17 00:00:00 2001 From: Archil Durglishvili Date: Sat, 9 Nov 2024 21:01:17 +0100 Subject: [PATCH 04/17] fix in the HCal neighbours --- .../detectorSegmentations/FCCSWHCalPhiTheta_k4geo.h | 13 +++---------- detectorSegmentations/src/FCCSWHCalPhiRow_k4geo.cpp | 6 ++++-- .../src/FCCSWHCalPhiTheta_k4geo.cpp | 9 +++++---- 3 files changed, 12 insertions(+), 16 deletions(-) diff --git a/detectorSegmentations/include/detectorSegmentations/FCCSWHCalPhiTheta_k4geo.h b/detectorSegmentations/include/detectorSegmentations/FCCSWHCalPhiTheta_k4geo.h index 58338b57a..75427523b 100644 --- a/detectorSegmentations/include/detectorSegmentations/FCCSWHCalPhiTheta_k4geo.h +++ b/detectorSegmentations/include/detectorSegmentations/FCCSWHCalPhiTheta_k4geo.h @@ -25,9 +25,9 @@ namespace dd4hep { /// destructor virtual ~FCCSWHCalPhiTheta_k4geo() = default; - /** Determine the position of HCal cell based on the cellID. - * @param[in] aCellId ID of a cell. - * return Position. + /** Get the postion of the geometric center of the cell based on the cellID + * @param[in] aCellID + * return the global coordinates of cell center */ virtual Vector3D position(const CellID& aCellID) const; @@ -117,13 +117,6 @@ namespace dd4hep { */ inline const std::string& fieldNamePhi() const { return m_phiID; } - /** Get the postion og the geometric center of the cell based on the cellID - * @param[in] aCellID - * return the global coordinates of cell center - */ - inline Vector3D centerPosition(const CellID& cID) const; - - /** Determine the minimum and maximum polar angle of HCal cell based on the cellID. * @param[in] aCellId ID of a cell. * return Theta. diff --git a/detectorSegmentations/src/FCCSWHCalPhiRow_k4geo.cpp b/detectorSegmentations/src/FCCSWHCalPhiRow_k4geo.cpp index 405261bdf..511f4adb4 100644 --- a/detectorSegmentations/src/FCCSWHCalPhiRow_k4geo.cpp +++ b/detectorSegmentations/src/FCCSWHCalPhiRow_k4geo.cpp @@ -614,8 +614,10 @@ std::vector FCCSWHCalPhiRow_k4geo::neighbours(const CellID& cID) const } } +/* // Now loop over the neighbours and add the cells from next/previous phi module - for(auto nID : cellNeighbours) + std::vector cellNeighboursCopy(cellNeighbours); + for(auto nID : cellNeighboursCopy) { CellID newID = nID; // previous: if the current is 0 then previous is the last bin (id = m_phiBins - 1) else current - 1 @@ -634,7 +636,7 @@ std::vector FCCSWHCalPhiRow_k4geo::neighbours(const CellID& cID) const // next: if the current is the last bin (id = m_phiBins - 1) then the next is the first bin (id = 0) else current + 1 _decoder->set(nID,m_phiID, (_decoder->get(cID,m_phiID) == (m_phiBins - 1)) ? 0 : _decoder->get(cID,m_phiID) + 1); cellNeighbours.push_back(nID); - +*/ return cellNeighbours; } diff --git a/detectorSegmentations/src/FCCSWHCalPhiTheta_k4geo.cpp b/detectorSegmentations/src/FCCSWHCalPhiTheta_k4geo.cpp index 765b9b3fc..6477d0a5d 100644 --- a/detectorSegmentations/src/FCCSWHCalPhiTheta_k4geo.cpp +++ b/detectorSegmentations/src/FCCSWHCalPhiTheta_k4geo.cpp @@ -38,7 +38,7 @@ FCCSWHCalPhiTheta_k4geo::FCCSWHCalPhiTheta_k4geo(const BitFieldCoder* decoder) : } -/// determine the global position based on the cell ID +/** /// determine the global position based on the cell ID Vector3D FCCSWHCalPhiTheta_k4geo::position(const CellID& cID) const { uint layer = _decoder->get(cID,"layer"); double radius = 1.0; @@ -48,11 +48,11 @@ Vector3D FCCSWHCalPhiTheta_k4geo::position(const CellID& cID) const { return positionFromRThetaPhi(radius, theta(cID), phi(cID)); } - +**/ /// determine the global position based on the cell ID /// returns the geometric center of the cell -Vector3D FCCSWHCalPhiTheta_k4geo::centerPosition(const CellID& cID) const { +Vector3D FCCSWHCalPhiTheta_k4geo::position(const CellID& cID) const { uint layer = _decoder->get(cID,"layer"); int thetaID = _decoder->get(cID,"theta"); double zpos = 0.; @@ -678,7 +678,8 @@ std::vector FCCSWHCalPhiTheta_k4geo::neighbours(const CellID& cID) con } // Now loop over the neighbours and add the cells from next/previous phi module - for(auto nID : cellNeighbours) + std::vector cellNeighboursCopy(cellNeighbours); + for(auto nID : cellNeighboursCopy) { CellID newID = nID; // previous: if the current is 0 then previous is the last bin (id = m_phiBins - 1) else current - 1 From af6c9699a1cde05577a7a329f885b2ff54d44472 Mon Sep 17 00:00:00 2001 From: Archil Durglishvili Date: Sat, 9 Nov 2024 21:57:45 +0100 Subject: [PATCH 05/17] New HCal readouts in ALLEGRO_o1_v04 --- .../compact/ALLEGRO_o1_v04/ALLEGRO_o1_v04.xml | 4 +- .../ALLEGRO_o1_v04/HCalBarrel_TileCal_v03.xml | 146 ++++++++++++++++ .../HCalEndcaps_ThreeParts_TileCal_v03.xml | 158 ++++++++++++++++++ 3 files changed, 306 insertions(+), 2 deletions(-) create mode 100644 FCCee/ALLEGRO/compact/ALLEGRO_o1_v04/HCalBarrel_TileCal_v03.xml create mode 100644 FCCee/ALLEGRO/compact/ALLEGRO_o1_v04/HCalEndcaps_ThreeParts_TileCal_v03.xml diff --git a/FCCee/ALLEGRO/compact/ALLEGRO_o1_v04/ALLEGRO_o1_v04.xml b/FCCee/ALLEGRO/compact/ALLEGRO_o1_v04/ALLEGRO_o1_v04.xml index a0a5dcf1a..4acb0d6f3 100644 --- a/FCCee/ALLEGRO/compact/ALLEGRO_o1_v04/ALLEGRO_o1_v04.xml +++ b/FCCee/ALLEGRO/compact/ALLEGRO_o1_v04/ALLEGRO_o1_v04.xml @@ -45,9 +45,9 @@ - + - + diff --git a/FCCee/ALLEGRO/compact/ALLEGRO_o1_v04/HCalBarrel_TileCal_v03.xml b/FCCee/ALLEGRO/compact/ALLEGRO_o1_v04/HCalBarrel_TileCal_v03.xml new file mode 100644 index 000000000..a62dab702 --- /dev/null +++ b/FCCee/ALLEGRO/compact/ALLEGRO_o1_v04/HCalBarrel_TileCal_v03.xml @@ -0,0 +1,146 @@ + + + + The first FCCee HCal layout based on ATLAS HCal, not optimised yet + 1. Nov 2022, J. Faltova: update material and radial segmentation for FCCee + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + system:4,layer:5,row:9,theta:9,phi:10 + + + + system:4,layer:5,row:9,phi:10 + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/FCCee/ALLEGRO/compact/ALLEGRO_o1_v04/HCalEndcaps_ThreeParts_TileCal_v03.xml b/FCCee/ALLEGRO/compact/ALLEGRO_o1_v04/HCalEndcaps_ThreeParts_TileCal_v03.xml new file mode 100644 index 000000000..c630263ca --- /dev/null +++ b/FCCee/ALLEGRO/compact/ALLEGRO_o1_v04/HCalEndcaps_ThreeParts_TileCal_v03.xml @@ -0,0 +1,158 @@ + + + + HCal layout based on ATLAS HCal, with realistic longitudinal segmentation and steel support + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + system:4,type:3,layer:6,row:11,theta:11,phi:10 + + + + system:4,type:3,layer:6,row:-8,phi:10 + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + From 4cfd76955476ccc0206d27e0d3ee5d68e5cf9cb7 Mon Sep 17 00:00:00 2001 From: Archil-AD Date: Thu, 14 Nov 2024 10:31:04 +0100 Subject: [PATCH 06/17] Update detectorSegmentations/src/FCCSWHCalPhiRow_k4geo.cpp fixing typo Co-authored-by: Andre Sailer --- detectorSegmentations/src/FCCSWHCalPhiRow_k4geo.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/detectorSegmentations/src/FCCSWHCalPhiRow_k4geo.cpp b/detectorSegmentations/src/FCCSWHCalPhiRow_k4geo.cpp index 511f4adb4..a1be32584 100644 --- a/detectorSegmentations/src/FCCSWHCalPhiRow_k4geo.cpp +++ b/detectorSegmentations/src/FCCSWHCalPhiRow_k4geo.cpp @@ -67,7 +67,7 @@ Vector3D FCCSWHCalPhiRow_k4geo::position(const CellID& cID) const { // calculate z-coordinate of the cell center double zpos = minLayerZ + (idx-1) * m_dz_row * m_gridSizeRow[layer] + 0.5 * m_dz_row * m_gridSizeRow[layer]; - // for negative-z Endcap, the index is negetive (starts from -1!) + // for negative-z Endcap, the index is negative (starts from -1!) if(idx < 0) zpos = -minLayerZ + (idx+1) * m_dz_row * m_gridSizeRow[layer] - 0.5 * m_dz_row * m_gridSizeRow[layer]; return Vector3D( radius * std::cos(phi), radius * std::sin(phi), zpos ); From c12027a0140a481cc6ef8f662a1729bcdd12c824 Mon Sep 17 00:00:00 2001 From: Archil-AD Date: Thu, 14 Nov 2024 10:48:34 +0100 Subject: [PATCH 07/17] Update description of grid_size_row parameter --- detectorSegmentations/src/FCCSWHCalPhiRow_k4geo.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/detectorSegmentations/src/FCCSWHCalPhiRow_k4geo.cpp b/detectorSegmentations/src/FCCSWHCalPhiRow_k4geo.cpp index a1be32584..4a4202c11 100644 --- a/detectorSegmentations/src/FCCSWHCalPhiRow_k4geo.cpp +++ b/detectorSegmentations/src/FCCSWHCalPhiRow_k4geo.cpp @@ -16,7 +16,7 @@ FCCSWHCalPhiRow_k4geo::FCCSWHCalPhiRow_k4geo(const std::string& cellEncoding) : // register all necessary parameters registerParameter("phi_bins", "Number of bins phi", m_phiBins, 1); registerParameter("offset_phi", "Angular offset in phi", m_offsetPhi, 0., SegmentationParameter::AngleUnit, true); - registerParameter("grid_size_row", "Cell size in row", m_gridSizeRow, std::vector()); + registerParameter("grid_size_row", "Number of rows combined in a cell", m_gridSizeRow, std::vector()); registerParameter("dz_row", "dz of row", m_dz_row, 0.); registerParameter("offset_z", "Offset in z-axis of the layer center", m_offsetZ, std::vector()); registerParameter("width_z", "Width in z of the layer", m_widthZ, std::vector()); @@ -35,7 +35,7 @@ FCCSWHCalPhiRow_k4geo::FCCSWHCalPhiRow_k4geo(const BitFieldCoder* decoder) : Seg // register all necessary parameters registerParameter("phi_bins", "Number of bins phi", m_phiBins, 1); registerParameter("offset_phi", "Angular offset in phi", m_offsetPhi, 0., SegmentationParameter::AngleUnit, true); - registerParameter("grid_size_row", "Cell size in row", m_gridSizeRow, std::vector()); + registerParameter("grid_size_row", "Number of rows combined in a cell", m_gridSizeRow, std::vector()); registerParameter("dz_row", "dz of row", m_dz_row, 0.); registerParameter("offset_z", "Offset in z-axis of the layer center", m_offsetZ, std::vector()); registerParameter("width_z", "Width in z of the layer", m_widthZ, std::vector()); From 73a3fb666c83e2c039ce30f43f654cf7966f0b27 Mon Sep 17 00:00:00 2001 From: Archil Durglishvili Date: Sun, 17 Nov 2024 14:21:42 +0100 Subject: [PATCH 08/17] Added units in HCal xml files --- .../compact/ALLEGRO_o1_v03/HCalBarrel_TileCal_v03.xml | 3 ++- .../ALLEGRO_o1_v03/HCalEndcaps_ThreeParts_TileCal_v03.xml | 3 ++- .../compact/ALLEGRO_o1_v04/HCalBarrel_TileCal_v03.xml | 5 +++-- .../ALLEGRO_o1_v04/HCalEndcaps_ThreeParts_TileCal_v03.xml | 5 +++-- 4 files changed, 10 insertions(+), 6 deletions(-) diff --git a/FCCee/ALLEGRO/compact/ALLEGRO_o1_v03/HCalBarrel_TileCal_v03.xml b/FCCee/ALLEGRO/compact/ALLEGRO_o1_v03/HCalBarrel_TileCal_v03.xml index a62dab702..f5ce03fe6 100644 --- a/FCCee/ALLEGRO/compact/ALLEGRO_o1_v03/HCalBarrel_TileCal_v03.xml +++ b/FCCee/ALLEGRO/compact/ALLEGRO_o1_v03/HCalBarrel_TileCal_v03.xml @@ -60,7 +60,7 @@ + diff --git a/FCCee/ALLEGRO/compact/ALLEGRO_o1_v03/HCalEndcaps_ThreeParts_TileCal_v03.xml b/FCCee/ALLEGRO/compact/ALLEGRO_o1_v03/HCalEndcaps_ThreeParts_TileCal_v03.xml index c630263ca..42e6904ef 100644 --- a/FCCee/ALLEGRO/compact/ALLEGRO_o1_v03/HCalEndcaps_ThreeParts_TileCal_v03.xml +++ b/FCCee/ALLEGRO/compact/ALLEGRO_o1_v03/HCalEndcaps_ThreeParts_TileCal_v03.xml @@ -59,7 +59,7 @@ + system:4,layer:5,row:9,theta:9,phi:10 - + + diff --git a/FCCee/ALLEGRO/compact/ALLEGRO_o1_v04/HCalEndcaps_ThreeParts_TileCal_v03.xml b/FCCee/ALLEGRO/compact/ALLEGRO_o1_v04/HCalEndcaps_ThreeParts_TileCal_v03.xml index c630263ca..48c59bb81 100644 --- a/FCCee/ALLEGRO/compact/ALLEGRO_o1_v04/HCalEndcaps_ThreeParts_TileCal_v03.xml +++ b/FCCee/ALLEGRO/compact/ALLEGRO_o1_v04/HCalEndcaps_ThreeParts_TileCal_v03.xml @@ -56,10 +56,10 @@ offset_phi="-pi+(pi/EndcapHCal_n_phi_modules)"/> system:4,type:3,layer:6,row:11,theta:11,phi:10 - + + Date: Sun, 17 Nov 2024 14:29:43 +0100 Subject: [PATCH 09/17] Fix in phi position calculation and added phi neigbours for HCal PhiRow segmentation --- detectorSegmentations/src/FCCSWHCalPhiRow_k4geo.cpp | 8 +++----- 1 file changed, 3 insertions(+), 5 deletions(-) diff --git a/detectorSegmentations/src/FCCSWHCalPhiRow_k4geo.cpp b/detectorSegmentations/src/FCCSWHCalPhiRow_k4geo.cpp index 4a4202c11..9cabe64ff 100644 --- a/detectorSegmentations/src/FCCSWHCalPhiRow_k4geo.cpp +++ b/detectorSegmentations/src/FCCSWHCalPhiRow_k4geo.cpp @@ -50,7 +50,6 @@ FCCSWHCalPhiRow_k4geo::FCCSWHCalPhiRow_k4geo(const BitFieldCoder* decoder) : Seg /// determine the global position based on the cell ID Vector3D FCCSWHCalPhiRow_k4geo::position(const CellID& cID) const { uint layer = _decoder->get(cID,"layer"); - double phi = _decoder->get(cID,m_phiID); if(m_radii.empty()) calculateLayerRadii(); if(m_radii.empty() || m_layerEdges.empty()) @@ -70,7 +69,7 @@ Vector3D FCCSWHCalPhiRow_k4geo::position(const CellID& cID) const { // for negative-z Endcap, the index is negative (starts from -1!) if(idx < 0) zpos = -minLayerZ + (idx+1) * m_dz_row * m_gridSizeRow[layer] - 0.5 * m_dz_row * m_gridSizeRow[layer]; - return Vector3D( radius * std::cos(phi), radius * std::sin(phi), zpos ); + return Vector3D( radius * std::cos(phi(cID)), radius * std::sin(phi(cID)), zpos ); } @@ -212,7 +211,7 @@ void FCCSWHCalPhiRow_k4geo::defineCellIndexes(const uint layer) const { double z1 = minLayerZ + (idx-1) * m_dz_row * m_gridSizeRow[layer]; // lower edge double z2 = minLayerZ + (idx) * m_dz_row * m_gridSizeRow[layer]; // upper edge - // for negative-z Endcap, the index is negetive (starts from -1!) + // for negative-z Endcap, the index is negative (starts from -1!) if(idx < 0) { z1 = -minLayerZ + (idx) * m_dz_row * m_gridSizeRow[layer]; // lower edge @@ -614,7 +613,6 @@ std::vector FCCSWHCalPhiRow_k4geo::neighbours(const CellID& cID) const } } -/* // Now loop over the neighbours and add the cells from next/previous phi module std::vector cellNeighboursCopy(cellNeighbours); for(auto nID : cellNeighboursCopy) @@ -636,7 +634,7 @@ std::vector FCCSWHCalPhiRow_k4geo::neighbours(const CellID& cID) const // next: if the current is the last bin (id = m_phiBins - 1) then the next is the first bin (id = 0) else current + 1 _decoder->set(nID,m_phiID, (_decoder->get(cID,m_phiID) == (m_phiBins - 1)) ? 0 : _decoder->get(cID,m_phiID) + 1); cellNeighbours.push_back(nID); -*/ + return cellNeighbours; } From bda1517c9348e753df5d0551ef0270f5417c9a52 Mon Sep 17 00:00:00 2001 From: Archil Durglishvili Date: Sun, 17 Nov 2024 14:34:34 +0100 Subject: [PATCH 10/17] Added an option to add diagonal neighbour cells for HCal PhiTheta segmentation --- .../FCCSWHCalPhiTheta_k4geo.h | 2 +- .../src/FCCSWHCalPhiTheta_k4geo.cpp | 181 ++++++++++++++---- 2 files changed, 145 insertions(+), 38 deletions(-) diff --git a/detectorSegmentations/include/detectorSegmentations/FCCSWHCalPhiTheta_k4geo.h b/detectorSegmentations/include/detectorSegmentations/FCCSWHCalPhiTheta_k4geo.h index 75427523b..607caa2c4 100644 --- a/detectorSegmentations/include/detectorSegmentations/FCCSWHCalPhiTheta_k4geo.h +++ b/detectorSegmentations/include/detectorSegmentations/FCCSWHCalPhiTheta_k4geo.h @@ -44,7 +44,7 @@ namespace dd4hep { * @param[in] aCellId ID of a cell. * return vector of neighbour cellIDs. */ - std::vector neighbours(const CellID& cID) const; + std::vector neighbours(const CellID& cID, bool aDiagonal) const; /** Calculate layer radii and edges in z-axis, then define cell edges in each layer using defineCellEdges(). */ diff --git a/detectorSegmentations/src/FCCSWHCalPhiTheta_k4geo.cpp b/detectorSegmentations/src/FCCSWHCalPhiTheta_k4geo.cpp index 6477d0a5d..04c227ab1 100644 --- a/detectorSegmentations/src/FCCSWHCalPhiTheta_k4geo.cpp +++ b/detectorSegmentations/src/FCCSWHCalPhiTheta_k4geo.cpp @@ -334,7 +334,7 @@ std::vector > FCCSWHCalPhiTheta_k4geo::getMinMaxLayerId() c /// Calculates the neighbours of the given cell ID and adds them to the list of neighbours -std::vector FCCSWHCalPhiTheta_k4geo::neighbours(const CellID& cID) const { +std::vector FCCSWHCalPhiTheta_k4geo::neighbours(const CellID& cID, bool aDiagonal) const { std::vector cellNeighbours; if(m_radii.empty()) defineCellsInRZplan(); @@ -448,6 +448,9 @@ std::vector FCCSWHCalPhiTheta_k4geo::neighbours(const CellID& cID) con // deal with the Barrel if(m_offsetZ.size() == 1) { + double currentCellZmin = m_cellEdges[currentLayerId][currentCellThetaBin].first; + double currentCellZmax = m_cellEdges[currentLayerId][currentCellThetaBin].second; + // if this is not the first layer then look for neighbours in the previous layer if(currentLayerId > minLayerId) { @@ -463,21 +466,40 @@ std::vector FCCSWHCalPhiTheta_k4geo::neighbours(const CellID& cID) con { _decoder->set(nID,m_thetaID,currentCellThetaBin - 1); cellNeighbours.push_back(nID); // add the cell from the previous layer of the same phi module + if(aDiagonal && currentCellThetaBin > (m_thetaBins[prevLayerId].front()+1)) + { + //add the previous layer cell from the prev to prev theta bin if it overlaps with the current cell in z-coordinate + double zmin = m_cellEdges[prevLayerId][currentCellThetaBin - 2].first; + if(zmin <= currentCellZmax) + { + //add the previous layer cell from the prev to prev theta bin + _decoder->set(nID,m_thetaID,currentCellThetaBin - 2); + cellNeighbours.push_back(nID); + } + } } // if the cID is in the negative-z side and prev layer cell is not in the last theta bin then add the cell from previous theta bin if(theta(cID) > M_PI/2. && currentCellThetaBin < m_thetaBins[prevLayerId].back() ) { _decoder->set(nID,m_thetaID,currentCellThetaBin + 1); cellNeighbours.push_back(nID); // add the cell from the previous layer of the same phi module + if(aDiagonal && currentCellThetaBin < (m_thetaBins[prevLayerId].back()-1)) + { + //add the previous layer cell from the next to next theta bin if it overlaps with the current cell in z-coordinate + double zmax = m_cellEdges[prevLayerId][currentCellThetaBin + 2].second; + if(zmax >= currentCellZmin) + { + //add the previous layer cell from the next to next theta bin + _decoder->set(nID,m_thetaID,currentCellThetaBin + 2); + cellNeighbours.push_back(nID); + } + } } } // if this is not the last layer then look for neighbours in the next layer if(currentLayerId < maxLayerId) { - double currentCellZmin = m_cellEdges[currentLayerId][currentCellThetaBin].first; - double currentCellZmax = m_cellEdges[currentLayerId][currentCellThetaBin].second; - CellID nID = cID ; int nextLayerId = currentLayerId + 1; _decoder->set(nID,"layer",nextLayerId); @@ -492,13 +514,16 @@ std::vector FCCSWHCalPhiTheta_k4geo::neighbours(const CellID& cID) con _decoder->set(nID,m_thetaID,currentCellThetaBin + 1); cellNeighbours.push_back(nID); - //add the next layer cell from the next-to-next theta bin if it overlaps with the current cell in z-coordinate - double zmax = m_cellEdges[nextLayerId][currentCellThetaBin + 2].second; - if(zmax >= currentCellZmin) + if(aDiagonal) { - //add the next layer cell from the next to next theta bin - _decoder->set(nID,m_thetaID,currentCellThetaBin + 2); - cellNeighbours.push_back(nID); + //add the next layer cell from the next-to-next theta bin if it overlaps with the current cell in z-coordinate + double zmax = m_cellEdges[nextLayerId][currentCellThetaBin + 2].second; + if(zmax >= currentCellZmin) + { + //add the next layer cell from the next to next theta bin + _decoder->set(nID,m_thetaID,currentCellThetaBin + 2); + cellNeighbours.push_back(nID); + } } } // if the cID is in the negative-z side @@ -508,13 +533,16 @@ std::vector FCCSWHCalPhiTheta_k4geo::neighbours(const CellID& cID) con _decoder->set(nID,m_thetaID,currentCellThetaBin - 1); cellNeighbours.push_back(nID); - //add the next layer cell from the next-to-next theta bin if it overlaps with the current cell in z-coordinate - double zmin = m_cellEdges[nextLayerId][currentCellThetaBin - 2].first; - if(zmin <= currentCellZmax) + if(aDiagonal) { - //add the next layer cell from the next to next theta bin - _decoder->set(nID,m_thetaID,currentCellThetaBin - 2); - cellNeighbours.push_back(nID); + //add the next layer cell from the prev to prev theta bin if it overlaps with the current cell in z-coordinate + double zmin = m_cellEdges[nextLayerId][currentCellThetaBin - 2].first; + if(zmin <= currentCellZmax) + { + //add the next layer cell from the prev to prev theta bin + _decoder->set(nID,m_thetaID,currentCellThetaBin - 2); + cellNeighbours.push_back(nID); + } } } } @@ -538,13 +566,39 @@ std::vector FCCSWHCalPhiTheta_k4geo::neighbours(const CellID& cID) con double zmin = m_cellEdges[prevLayerId][bin].first; double zmax = m_cellEdges[prevLayerId][bin].second; - if( (zmin >= currentCellZmin && zmin <=currentCellZmax) - || (zmax >= currentCellZmin && zmax <=currentCellZmax) - || (currentCellZmin >= zmin && currentCellZmax <= zmax) - ) + // if the cID is in the positive-z side + if(theta(cID) < M_PI/2.) { - _decoder->set(nID,m_thetaID,bin); - cellNeighbours.push_back(nID); // add the cell from the previous layer of the same phi module + if( (zmin >= currentCellZmin && zmin < currentCellZmax) + || (zmax >= currentCellZmin && zmax <= currentCellZmax) + || (currentCellZmin >= zmin && currentCellZmax <= zmax) + ) + { + _decoder->set(nID,m_thetaID,bin); + cellNeighbours.push_back(nID); // add the cell from the previous layer of the same phi module + } + if(aDiagonal && zmin == currentCellZmax) + { + _decoder->set(nID,m_thetaID,bin); + cellNeighbours.push_back(nID); // add the cell from the previous layer of the same phi module + } + } + // if the cID is in the negative-z side + if(theta(cID) > M_PI/2.) + { + if( (zmin >= currentCellZmin && zmin <= currentCellZmax) + || (zmax > currentCellZmin && zmax <= currentCellZmax) + || (currentCellZmin >= zmin && currentCellZmax <= zmax) + ) + { + _decoder->set(nID,m_thetaID,bin); + cellNeighbours.push_back(nID); // add the cell from the previous layer of the same phi module + } + if(aDiagonal && zmax == currentCellZmin) + { + _decoder->set(nID,m_thetaID,bin); + cellNeighbours.push_back(nID); // add the cell from the previous layer of the same phi module + } } } } @@ -559,14 +613,39 @@ std::vector FCCSWHCalPhiTheta_k4geo::neighbours(const CellID& cID) con { double zmin = m_cellEdges[nextLayerId][bin].first; double zmax = m_cellEdges[nextLayerId][bin].second; - - if( (zmin >= currentCellZmin && zmin <=currentCellZmax) - || (zmax >= currentCellZmin && zmax <=currentCellZmax) - || (currentCellZmin >= zmin && currentCellZmax <= zmax) - ) + // if the cID is in the positive-z side + if(theta(cID) < M_PI/2.) + { + if( (zmin >= currentCellZmin && zmin <=currentCellZmax) + || (zmax > currentCellZmin && zmax <=currentCellZmax) + || (currentCellZmin >= zmin && currentCellZmax <= zmax) + ) + { + _decoder->set(nID,m_thetaID,bin); + cellNeighbours.push_back(nID); // add the cell from the next layer of the same phi module + } + if(aDiagonal && zmax == currentCellZmin) + { + _decoder->set(nID,m_thetaID,bin); + cellNeighbours.push_back(nID); // add the cell from the next layer of the same phi module + } + } + // if the cID is in the negative-z side + if(theta(cID) > M_PI/2.) { - _decoder->set(nID,m_thetaID,bin); - cellNeighbours.push_back(nID); // add the cell from the previous layer of the same phi module + if( (zmin >= currentCellZmin && zmin < currentCellZmax) + || (zmax >= currentCellZmin && zmax <=currentCellZmax) + || (currentCellZmin >= zmin && currentCellZmax <= zmax) + ) + { + _decoder->set(nID,m_thetaID,bin); + cellNeighbours.push_back(nID); // add the cell from the next layer of the same phi module + } + if(aDiagonal && zmin == currentCellZmax) + { + _decoder->set(nID,m_thetaID,bin); + cellNeighbours.push_back(nID); // add the cell from the next layer of the same phi module + } } } } @@ -594,8 +673,8 @@ std::vector FCCSWHCalPhiTheta_k4geo::neighbours(const CellID& cID) con double Rmax = m_radii[part2layerId] + 0.5*m_layerDepth[part2layerId]; if( (Rmin >= currentLayerRmin && Rmin <= currentLayerRmax) - || (Rmax >= currentLayerRmin && Rmax <= currentLayerRmax) - || (currentLayerRmin >= Rmin && currentLayerRmin <= Rmax) + || (Rmax > currentLayerRmin && Rmax <= currentLayerRmax) + || (currentLayerRmin >= Rmin && currentLayerRmin < Rmax) || (currentLayerRmax >= Rmin && currentLayerRmax <= Rmax) ) { @@ -604,6 +683,13 @@ std::vector FCCSWHCalPhiTheta_k4geo::neighbours(const CellID& cID) con _decoder->set(nID,m_thetaID,maxCellThetaBin); cellNeighbours.push_back(nID); // add the last theta bin cell from part2 layer } + if(aDiagonal && Rmax == currentLayerRmin) + { + CellID nID = cID ; + _decoder->set(nID,"layer",part2layerId); + _decoder->set(nID,m_thetaID,maxCellThetaBin); + cellNeighbours.push_back(nID); // add the last theta bin cell from part2 layer + } } } @@ -616,10 +702,10 @@ std::vector FCCSWHCalPhiTheta_k4geo::neighbours(const CellID& cID) con double Rmin = m_radii[part1layerId] - 0.5*m_layerDepth[part1layerId]; double Rmax = m_radii[part1layerId] + 0.5*m_layerDepth[part1layerId]; - if( (Rmin >= currentLayerRmin && Rmin <= currentLayerRmax) + if( (Rmin >= currentLayerRmin && Rmin < currentLayerRmax) || (Rmax >= currentLayerRmin && Rmax <= currentLayerRmax) || (currentLayerRmin >= Rmin && currentLayerRmin <= Rmax) - || (currentLayerRmax >= Rmin && currentLayerRmax <= Rmax) + || (currentLayerRmax > Rmin && currentLayerRmax <= Rmax) ) { CellID nID = cID ; @@ -627,6 +713,13 @@ std::vector FCCSWHCalPhiTheta_k4geo::neighbours(const CellID& cID) con _decoder->set(nID,m_thetaID,minCellThetaBin); cellNeighbours.push_back(nID); // add the first theta bin cell from the part1 layer } + if(aDiagonal && Rmin == currentLayerRmax) + { + CellID nID = cID ; + _decoder->set(nID,"layer",part1layerId); + _decoder->set(nID,m_thetaID,minCellThetaBin); + cellNeighbours.push_back(nID); // add the first theta bin cell from the part1 layer + } } } @@ -640,8 +733,8 @@ std::vector FCCSWHCalPhiTheta_k4geo::neighbours(const CellID& cID) con double Rmax = m_radii[part3layerId] + 0.5*m_layerDepth[part3layerId]; if( (Rmin >= currentLayerRmin && Rmin <= currentLayerRmax) - || (Rmax >= currentLayerRmin && Rmax <= currentLayerRmax) - || (currentLayerRmin >= Rmin && currentLayerRmin <= Rmax) + || (Rmax > currentLayerRmin && Rmax <= currentLayerRmax) + || (currentLayerRmin >= Rmin && currentLayerRmin < Rmax) || (currentLayerRmax >= Rmin && currentLayerRmax <= Rmax) ) { @@ -650,6 +743,13 @@ std::vector FCCSWHCalPhiTheta_k4geo::neighbours(const CellID& cID) con _decoder->set(nID,m_thetaID,maxCellThetaBin); cellNeighbours.push_back(nID); // add the first cell from the part3 layer } + if(aDiagonal && Rmax == currentLayerRmin) + { + CellID nID = cID ; + _decoder->set(nID,"layer",part3layerId); + _decoder->set(nID,m_thetaID,maxCellThetaBin); + cellNeighbours.push_back(nID); // add the first cell from the part3 layer + } } } @@ -662,10 +762,10 @@ std::vector FCCSWHCalPhiTheta_k4geo::neighbours(const CellID& cID) con double Rmin = m_radii[part2layerId] - 0.5*m_layerDepth[part2layerId]; double Rmax = m_radii[part2layerId] + 0.5*m_layerDepth[part2layerId]; - if( (Rmin >= currentLayerRmin && Rmin <= currentLayerRmax) + if( (Rmin >= currentLayerRmin && Rmin < currentLayerRmax) || (Rmax >= currentLayerRmin && Rmax <= currentLayerRmax) || (currentLayerRmin >= Rmin && currentLayerRmin <= Rmax) - || (currentLayerRmax >= Rmin && currentLayerRmax <= Rmax) + || (currentLayerRmax > Rmin && currentLayerRmax <= Rmax) ) { CellID nID = cID ; @@ -673,6 +773,13 @@ std::vector FCCSWHCalPhiTheta_k4geo::neighbours(const CellID& cID) con _decoder->set(nID,m_thetaID,minCellThetaBin); cellNeighbours.push_back(nID); // add the first theta bin cell from the part2 layer } + if(aDiagonal && Rmin == currentLayerRmax) + { + CellID nID = cID ; + _decoder->set(nID,"layer",part2layerId); + _decoder->set(nID,m_thetaID,minCellThetaBin); + cellNeighbours.push_back(nID); // add the first theta bin cell from the part2 layer + } } } } From 69c28eb3320c87c2f6302114f4757e7869c739ae Mon Sep 17 00:00:00 2001 From: Archil Durglishvili Date: Tue, 26 Nov 2024 16:07:45 +0100 Subject: [PATCH 11/17] adding a test for ALLEGRO v04 in test/CMakeLists.txt --- test/CMakeLists.txt | 9 +++++++++ 1 file changed, 9 insertions(+) diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index 0b0203867..46576dbf9 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -131,6 +131,15 @@ ADD_TEST( t_${test_name} "${CMAKE_INSTALL_PREFIX}/bin/run_test_${PackageName}.sh SET_TESTS_PROPERTIES( t_${test_name} PROPERTIES FAIL_REGULAR_EXPRESSION " Exception; EXCEPTION;ERROR;Error" ) endif() +#-------------------------------------------------- +# test for ALLEGRO o1 v04 +if(DCH_INFO_H_EXIST) +SET( test_name "test_ALLEGRO_o1_v04" ) +ADD_TEST( t_${test_name} "${CMAKE_INSTALL_PREFIX}/bin/run_test_${PackageName}.sh" + ddsim --compactFile=${CMAKE_CURRENT_SOURCE_DIR}/../FCCee/ALLEGRO/compact/ALLEGRO_o1_v04/ALLEGRO_o1_v04.xml --runType=batch -G -N=1 --outputFile=testALLEGRO_o1_v04.root --gun.direction "1,0,1" ) +SET_TESTS_PROPERTIES( t_${test_name} PROPERTIES FAIL_REGULAR_EXPRESSION " Exception;EXCEPTION;ERROR;Error" ) +endif() + #-------------------------------------------------- # test for ARC o1 v01 SET( test_name "test_ARC_o1_v01_run" ) From 20e98c7b60eab19b756ab9787d59720105b15b05 Mon Sep 17 00:00:00 2001 From: Archil Durglishvili Date: Thu, 28 Nov 2024 11:50:11 +0100 Subject: [PATCH 12/17] creating HCal symlinks in ALLEGRO_o1_v04 from ALLEGRO_o1_v03 --- .../ALLEGRO_o1_v04/HCalBarrel_TileCal_v03.xml | 148 +--------------- .../HCalEndcaps_ThreeParts_TileCal_v03.xml | 160 +----------------- 2 files changed, 2 insertions(+), 306 deletions(-) mode change 100644 => 120000 FCCee/ALLEGRO/compact/ALLEGRO_o1_v04/HCalBarrel_TileCal_v03.xml mode change 100644 => 120000 FCCee/ALLEGRO/compact/ALLEGRO_o1_v04/HCalEndcaps_ThreeParts_TileCal_v03.xml diff --git a/FCCee/ALLEGRO/compact/ALLEGRO_o1_v04/HCalBarrel_TileCal_v03.xml b/FCCee/ALLEGRO/compact/ALLEGRO_o1_v04/HCalBarrel_TileCal_v03.xml deleted file mode 100644 index e108c7bb3..000000000 --- a/FCCee/ALLEGRO/compact/ALLEGRO_o1_v04/HCalBarrel_TileCal_v03.xml +++ /dev/null @@ -1,147 +0,0 @@ - - - - The first FCCee HCal layout based on ATLAS HCal, not optimised yet - 1. Nov 2022, J. Faltova: update material and radial segmentation for FCCee - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - system:4,layer:5,row:9,theta:9,phi:10 - - - - system:4,layer:5,row:9,phi:10 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - diff --git a/FCCee/ALLEGRO/compact/ALLEGRO_o1_v04/HCalBarrel_TileCal_v03.xml b/FCCee/ALLEGRO/compact/ALLEGRO_o1_v04/HCalBarrel_TileCal_v03.xml new file mode 120000 index 000000000..1bde7347d --- /dev/null +++ b/FCCee/ALLEGRO/compact/ALLEGRO_o1_v04/HCalBarrel_TileCal_v03.xml @@ -0,0 +1 @@ +../ALLEGRO_o1_v03/HCalBarrel_TileCal_v03.xml \ No newline at end of file diff --git a/FCCee/ALLEGRO/compact/ALLEGRO_o1_v04/HCalEndcaps_ThreeParts_TileCal_v03.xml b/FCCee/ALLEGRO/compact/ALLEGRO_o1_v04/HCalEndcaps_ThreeParts_TileCal_v03.xml deleted file mode 100644 index 48c59bb81..000000000 --- a/FCCee/ALLEGRO/compact/ALLEGRO_o1_v04/HCalEndcaps_ThreeParts_TileCal_v03.xml +++ /dev/null @@ -1,159 +0,0 @@ - - - - HCal layout based on ATLAS HCal, with realistic longitudinal segmentation and steel support - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - system:4,type:3,layer:6,row:11,theta:11,phi:10 - - - - system:4,type:3,layer:6,row:-8,phi:10 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - diff --git a/FCCee/ALLEGRO/compact/ALLEGRO_o1_v04/HCalEndcaps_ThreeParts_TileCal_v03.xml b/FCCee/ALLEGRO/compact/ALLEGRO_o1_v04/HCalEndcaps_ThreeParts_TileCal_v03.xml new file mode 120000 index 000000000..9fb72eebd --- /dev/null +++ b/FCCee/ALLEGRO/compact/ALLEGRO_o1_v04/HCalEndcaps_ThreeParts_TileCal_v03.xml @@ -0,0 +1 @@ +../ALLEGRO_o1_v03/HCalEndcaps_ThreeParts_TileCal_v03.xml \ No newline at end of file From 8cc31d8ccc4fdb15b1392551cb450fdedee8cffc Mon Sep 17 00:00:00 2001 From: Archil Durglishvili Date: Sat, 30 Nov 2024 16:44:26 +0100 Subject: [PATCH 13/17] HCal Endcap segmentation class is generalized for N-parts endcap --- .../ALLEGRO_o1_v03/HCalBarrel_TileCal_v03.xml | 2 + .../HCalEndcaps_ThreeParts_TileCal_v03.xml | 2 + .../FCCSWHCalPhiRowHandle_k4geo.h | 6 + .../FCCSWHCalPhiRow_k4geo.h | 45 +- .../FCCSWHCalPhiThetaHandle_k4geo.h | 7 + .../FCCSWHCalPhiTheta_k4geo.h | 47 +- .../src/FCCSWHCalPhiRow_k4geo.cpp | 292 +++++------- .../src/FCCSWHCalPhiTheta_k4geo.cpp | 432 ++++++++---------- 8 files changed, 389 insertions(+), 444 deletions(-) diff --git a/FCCee/ALLEGRO/compact/ALLEGRO_o1_v03/HCalBarrel_TileCal_v03.xml b/FCCee/ALLEGRO/compact/ALLEGRO_o1_v03/HCalBarrel_TileCal_v03.xml index f5ce03fe6..1adb2a8fd 100644 --- a/FCCee/ALLEGRO/compact/ALLEGRO_o1_v03/HCalBarrel_TileCal_v03.xml +++ b/FCCee/ALLEGRO/compact/ALLEGRO_o1_v03/HCalBarrel_TileCal_v03.xml @@ -46,6 +46,7 @@ implementation->setOffsetPhi(offset); } + /// set the detector layout + inline void setDetLayout(int detLayout) const { access()->implementation->setDetLayout(detLayout); } + /// set the coordinate offset in z-axis inline void setOffsetZ(std::vector const&offset) const { access()->implementation->setOffsetZ(offset); } @@ -101,6 +104,9 @@ class FCCSWHCalPhiRow_k4geo : public FCCSWHCalPhiRowHandle_k4geo { /// set the grid size in Phi inline void setPhiBins(int cellSize) const { access()->implementation->setPhiBins(cellSize); } + /// access the field name used for layer + inline const std::string& fieldNameLayer() const { return access()->implementation->fieldNameLayer(); } + /// access the field name used for theta inline const std::string& fieldNameRow() const { return access()->implementation->fieldNameRow(); } diff --git a/detectorSegmentations/include/detectorSegmentations/FCCSWHCalPhiRow_k4geo.h b/detectorSegmentations/include/detectorSegmentations/FCCSWHCalPhiRow_k4geo.h index ce9b6f44f..acda782e8 100644 --- a/detectorSegmentations/include/detectorSegmentations/FCCSWHCalPhiRow_k4geo.h +++ b/detectorSegmentations/include/detectorSegmentations/FCCSWHCalPhiRow_k4geo.h @@ -11,6 +11,9 @@ * Segmentation in row and phi. * * Cells are defined in r-z plan by merging the row/sequence of scintillator/absorber. + * Each row consists of 2 * Master_Plate + 1 * Spacer_Plate + 1 * Scintillator. + * Considering a single row as a single cell gives the highest possible granularity. + * More details: https://indico.cern.ch/event/1475808/contributions/6219554/attachments/2966253/5218774/FCC_FullSim_HCal_slides.pdf * */ @@ -42,6 +45,7 @@ namespace dd4hep { const VolumeID& aVolumeID) const; /** Find neighbours of the cell + * Definition of neighbours is explained on slide 7: https://indico.cern.ch/event/1475808/contributions/6219554/attachments/2966253/5218774/FCC_FullSim_HCal_slides.pdf * @param[in] aCellId ID of a cell. * return vector of neighbour cellIDs. */ @@ -49,6 +53,10 @@ namespace dd4hep { /** Calculate layer radii and edges in z-axis. + * Following member variables are calculated: + * m_radii + * m_layerEdges + * m_layerDepth */ void calculateLayerRadii() const; @@ -58,7 +66,9 @@ namespace dd4hep { * In case of a cell with single row/sequence, the index is directly the number of row in the layer. * In case of a cell with several rows/sequences merged, the index is the number of cell in the layer. * For the layers of negative-z Endcap, indexes of cells are negative. - * + * Following member variables are calculated: + * m_cellIndexes + * m_cellEdges * @param[in] layer index */ void defineCellIndexes(const unsigned int layer) const; @@ -79,7 +89,7 @@ namespace dd4hep { */ inline int phiBins() const { return m_phiBins; } - /** Get the coordinate offset in azimuthal angle. + /** Get the coordinate offset in azimuthal angle, which is the center of cell with m_phiID=0. * return The offset in phi. */ inline double offsetPhi() const { return m_offsetPhi; } @@ -114,27 +124,34 @@ namespace dd4hep { */ std::vector > getMinMaxLayerId() const ; - /** Get the coordinate offset in z-axis. + /** Get the coordinate offset in z-axis. + * Offset is the middle position of the Barrel or each section of the Endcap. + * For the Barrel, the vector size is 1, while for the Endcap - number of section. * return The offset in z. */ inline std::vector offsetZ() const { return m_offsetZ; } - /** Get the z width of the layer. - * return the z width. + /** Get the z length of the layer. + * return the z length. */ inline std::vector widthZ() const { return m_widthZ; } /** Get the coordinate offset in radius. + * Offset is the inner radius of the first layer in the Barrel or in each section of the Endcap. + * For the Barrel, the vector size is 1, while for the Endcap - number of sections. * return the offset in radius. */ inline std::vector offsetR() const { return m_offsetR; } - /** Get the number of layers. + /** Get the number of layers for each different thickness retrieved with dRlayer(). + * For the Barrel, the vector size equals to the number of different thicknesses used to form the layers. + * For the Endcap, the vector size equals to the number of sections in the Endcap times the number of different thicknesses used to form the layers. * return the number of layers. */ inline std::vector numLayers() const { return m_numLayers; } - /** Get the dR of layers. + /** Get the dR (thickness) of layers. + * The size of the vector equals to the number of different thicknesses used to form the layers. * return the dR. */ inline std::vector dRlayer() const { return m_dRlayer; } @@ -144,6 +161,11 @@ namespace dd4hep { */ inline const std::string& fieldNamePhi() const { return m_phiID; } + /** Get the field name for layer. + * return The field name for layer. + */ + inline const std::string& fieldNameLayer() const { return m_layerID; } + /** Get the field name for row number. * return The field name for row. */ @@ -164,6 +186,11 @@ namespace dd4hep { */ inline void setGridSizeRow(std::vector const&size) { m_gridSizeRow = size; } + /** Set the detector layout (0 = Barrel; 1 = Endcap). + * @param[in] detLayout + */ + inline void setDetLayout(int detLayout) { m_detLayout = detLayout; } + /** Set the coordinate offset in z-axis. * @param[in] offset in z (centre of the layer). */ @@ -208,12 +235,16 @@ namespace dd4hep { double m_offsetPhi; /// the field name used for phi std::string m_phiID; + /// the field name used for layer + std::string m_layerID; /// the grid size in row for each layer std::vector m_gridSizeRow; /// dz of row double m_dz_row; /// the field name used for row std::string m_rowID; + /// the detector layout (0 = Barrel; 1 = Endcap) + int m_detLayout; /// the z offset of middle of the layer std::vector m_offsetZ; /// the z width of the layer diff --git a/detectorSegmentations/include/detectorSegmentations/FCCSWHCalPhiThetaHandle_k4geo.h b/detectorSegmentations/include/detectorSegmentations/FCCSWHCalPhiThetaHandle_k4geo.h index 475b39649..a9f128147 100644 --- a/detectorSegmentations/include/detectorSegmentations/FCCSWHCalPhiThetaHandle_k4geo.h +++ b/detectorSegmentations/include/detectorSegmentations/FCCSWHCalPhiThetaHandle_k4geo.h @@ -86,6 +86,9 @@ class FCCSWHCalPhiTheta_k4geo : public FCCSWHCalPhiThetaHandle_k4geo { /// set the coordinate offset in Phi inline void setOffsetPhi(double offset) const { access()->implementation->setOffsetPhi(offset); } + /// set the detector layout + inline void setDetLayout(int detLayout) const { access()->implementation->setDetLayout(detLayout); } + /// set the coordinate offset in z-axis inline void setOffsetZ(std::vector const&offset) const { access()->implementation->setOffsetZ(offset); } @@ -113,6 +116,10 @@ class FCCSWHCalPhiTheta_k4geo : public FCCSWHCalPhiThetaHandle_k4geo { /// access the field name used for Phi inline const std::string& fieldNamePhi() const { return access()->implementation->fieldNamePhi(); } + /// access the field name used for layer + inline const std::string& fieldNameLayer() const { return access()->implementation->fieldNameLayer(); } + + /** \brief Returns a std::vector of the cellDimensions of the given cell ID in natural order of dimensions (dPhi, dTheta) diff --git a/detectorSegmentations/include/detectorSegmentations/FCCSWHCalPhiTheta_k4geo.h b/detectorSegmentations/include/detectorSegmentations/FCCSWHCalPhiTheta_k4geo.h index 607caa2c4..e02ad0ede 100644 --- a/detectorSegmentations/include/detectorSegmentations/FCCSWHCalPhiTheta_k4geo.h +++ b/detectorSegmentations/include/detectorSegmentations/FCCSWHCalPhiTheta_k4geo.h @@ -10,6 +10,8 @@ * Based on GridTheta_k4geo, addition of azimuthal angle coordinate. * * Rectangular shape cells are defined in r-z plan and all the hits within a defined cell boundaries are assigned the same cellID. + * Cell borders are defined such that closely follow the theta projective towers. + * More details: https://indico.cern.ch/event/1475808/contributions/6219554/attachments/2966253/5218774/FCC_FullSim_HCal_slides.pdf * */ @@ -40,13 +42,21 @@ namespace dd4hep { virtual CellID cellID(const Vector3D& aLocalPosition, const Vector3D& aGlobalPosition, const VolumeID& aVolumeID) const; - /** Find neighbours of the cell + /** Find neighbours of the cell. + * Definition of neighbours is explained on slide 9: https://indico.cern.ch/event/1475808/contributions/6219554/attachments/2966253/5218774/FCC_FullSim_HCal_slides.pdf * @param[in] aCellId ID of a cell. + * @param[in] aDiagonal if true, will include neighbours from diagonal positions in the next and previous layers. * return vector of neighbour cellIDs. */ std::vector neighbours(const CellID& cID, bool aDiagonal) const; /** Calculate layer radii and edges in z-axis, then define cell edges in each layer using defineCellEdges(). + * Following member variables are calculated: + * m_radii + * m_layerEdges + * m_layerDepth + * m_thetaBins (updated through defineCellEdges()) + * m_cellEdges (updated through defineCellEdges()) */ void defineCellsInRZplan() const; @@ -74,12 +84,12 @@ namespace dd4hep { */ inline int phiBins() const { return m_phiBins; } - /** Get the coordinate offset in azimuthal angle. + /** Get the coordinate offset in azimuthal angle, which is the center of cell with m_phiID=0 * return The offset in phi. */ inline double offsetPhi() const { return m_offsetPhi; } - /** Get the vector of theta bins (cells) in a give layer. + /** Get the vector of theta bins (cells) in a given layer. */ inline std::vector thetaBins(const uint layer) const { if(m_radii.empty()) defineCellsInRZplan(); @@ -87,27 +97,34 @@ namespace dd4hep { else return std::vector(); } - /** Get the coordinate offset in z-axis. + /** Get the coordinate offset in z-axis. + * Offset is the middle position of the Barrel or each section of the Endcap. + * For the Barrel, the vector size is 1, while for the Endcap - number of section. * return The offset in z. */ inline std::vector offsetZ() const { return m_offsetZ; } - /** Get the z width of the layer. - * return the z width. + /** Get the z length of the layer. + * return the z length. */ inline std::vector widthZ() const { return m_widthZ; } /** Get the coordinate offset in radius. + * Offset is the inner radius of the first layer in the Barrel or in each section of the Endcap. + * For the Barrel, the vector size is 1, while for the Endcap - number of sections. * return the offset in radius. */ inline std::vector offsetR() const { return m_offsetR; } - /** Get the number of layers. + /** Get the number of layers for each different thickness retrieved with dRlayer(). + * For the Barrel, the vector size equals to the number of different thicknesses used to form the layers. + * For the Endcap, the vector size equals to the number of sections in the Endcap times the number of different thicknesses used to form the layers. * return the number of layers. */ inline std::vector numLayers() const { return m_numLayers; } - /** Get the dR of layers. + /** Get the dR (thickness) of layers. + * The size of the vector equals to the number of different thicknesses used to form the layers. * return the dR. */ inline std::vector dRlayer() const { return m_dRlayer; } @@ -117,6 +134,11 @@ namespace dd4hep { */ inline const std::string& fieldNamePhi() const { return m_phiID; } + /** Get the field name for layer. + * return The field name for layer. + */ + inline const std::string& fieldNameLayer() const { return m_layerID; } + /** Determine the minimum and maximum polar angle of HCal cell based on the cellID. * @param[in] aCellId ID of a cell. * return Theta. @@ -138,6 +160,11 @@ namespace dd4hep { */ inline void setOffsetPhi(double offset) { m_offsetPhi = offset; } + /** Set the detector layout (0 = Barrel; 1 = Endcap). + * @param[in] detLayout + */ + inline void setDetLayout(int detLayout) { m_detLayout = detLayout; } + /** Set the coordinate offset in z-axis. * @param[in] offset in z (centre of the layer). */ @@ -178,6 +205,10 @@ namespace dd4hep { double m_offsetPhi; /// the field name used for phi std::string m_phiID; + /// the field name used for layer + std::string m_layerID; + /// the detector layout (0 = Barrel; 1 = Endcap) + int m_detLayout; /// the z offset of middle of the layer std::vector m_offsetZ; /// the z width of the layer diff --git a/detectorSegmentations/src/FCCSWHCalPhiRow_k4geo.cpp b/detectorSegmentations/src/FCCSWHCalPhiRow_k4geo.cpp index 9cabe64ff..f17059abb 100644 --- a/detectorSegmentations/src/FCCSWHCalPhiRow_k4geo.cpp +++ b/detectorSegmentations/src/FCCSWHCalPhiRow_k4geo.cpp @@ -18,6 +18,7 @@ FCCSWHCalPhiRow_k4geo::FCCSWHCalPhiRow_k4geo(const std::string& cellEncoding) : registerParameter("offset_phi", "Angular offset in phi", m_offsetPhi, 0., SegmentationParameter::AngleUnit, true); registerParameter("grid_size_row", "Number of rows combined in a cell", m_gridSizeRow, std::vector()); registerParameter("dz_row", "dz of row", m_dz_row, 0.); + registerParameter("detLayout", "The detector layout (0 = Barrel; 1 = Endcap)", m_detLayout, -1); registerParameter("offset_z", "Offset in z-axis of the layer center", m_offsetZ, std::vector()); registerParameter("width_z", "Width in z of the layer", m_widthZ, std::vector()); registerParameter("offset_r", "Offset in radius of the layer (Rmin)", m_offsetR, std::vector()); @@ -25,6 +26,7 @@ FCCSWHCalPhiRow_k4geo::FCCSWHCalPhiRow_k4geo(const std::string& cellEncoding) : registerParameter("dRlayer", "dR of the layer", m_dRlayer, std::vector()); registerIdentifier("identifier_phi", "Cell ID identifier for phi", m_phiID, "phi"); registerIdentifier("identifier_row", "Cell ID identifier for row", m_rowID, "row"); + registerIdentifier("identifier_layer", "Cell ID identifier for layer", m_layerID, "layer"); } FCCSWHCalPhiRow_k4geo::FCCSWHCalPhiRow_k4geo(const BitFieldCoder* decoder) : Segmentation(decoder) { @@ -37,6 +39,7 @@ FCCSWHCalPhiRow_k4geo::FCCSWHCalPhiRow_k4geo(const BitFieldCoder* decoder) : Seg registerParameter("offset_phi", "Angular offset in phi", m_offsetPhi, 0., SegmentationParameter::AngleUnit, true); registerParameter("grid_size_row", "Number of rows combined in a cell", m_gridSizeRow, std::vector()); registerParameter("dz_row", "dz of row", m_dz_row, 0.); + registerParameter("detLayout", "The detector layout (0 = Barrel; 1 = Endcap)", m_detLayout, -1); registerParameter("offset_z", "Offset in z-axis of the layer center", m_offsetZ, std::vector()); registerParameter("width_z", "Width in z of the layer", m_widthZ, std::vector()); registerParameter("offset_r", "Offset in radius of the layer (Rmin)", m_offsetR, std::vector()); @@ -44,12 +47,13 @@ FCCSWHCalPhiRow_k4geo::FCCSWHCalPhiRow_k4geo(const BitFieldCoder* decoder) : Seg registerParameter("dRlayer", "dR of the layer", m_dRlayer, std::vector()); registerIdentifier("identifier_phi", "Cell ID identifier for phi", m_phiID, "phi"); registerIdentifier("identifier_row", "Cell ID identifier for row", m_rowID, "row"); + registerIdentifier("identifier_layer", "Cell ID identifier for layer", m_layerID, "layer"); } /// determine the global position based on the cell ID Vector3D FCCSWHCalPhiRow_k4geo::position(const CellID& cID) const { - uint layer = _decoder->get(cID,"layer"); + uint layer = _decoder->get(cID,m_layerID); if(m_radii.empty()) calculateLayerRadii(); if(m_radii.empty() || m_layerEdges.empty()) @@ -62,7 +66,7 @@ Vector3D FCCSWHCalPhiRow_k4geo::position(const CellID& cID) const { double minLayerZ = m_layerEdges[layer].first; // get index of the cell in the layer (index starts from 1!) - int idx = _decoder->get(cID, "row"); + int idx = _decoder->get(cID, m_rowID); // calculate z-coordinate of the cell center double zpos = minLayerZ + (idx-1) * m_dz_row * m_gridSizeRow[layer] + 0.5 * m_dz_row * m_gridSizeRow[layer]; @@ -77,76 +81,33 @@ void FCCSWHCalPhiRow_k4geo::calculateLayerRadii() const { if(m_radii.empty()) { // check if all necessary variables are available - if(m_offsetZ.empty() || m_widthZ.empty() || m_offsetR.empty() || m_numLayers.empty() || m_dRlayer.empty()) + if(m_detLayout==-1 || m_offsetZ.empty() || m_widthZ.empty() || m_offsetR.empty() || m_numLayers.empty() || m_dRlayer.empty()) { dd4hep::printout(dd4hep::ERROR, "FCCSWHCalPhiRow_k4geo","Please check the readout description in the XML file!", - "One of the variables is missing: offset_z | width_z | offset_r | numLayers | dRlayer"); + "One of the variables is missing: detLayout | offset_z | width_z | offset_r | numLayers | dRlayer"); return; } - // calculate the radius for each layer - if(m_offsetZ.size() == 1) // Barrel - { - dd4hep::printout(dd4hep::INFO, "FCCSWHCalPhiRow_k4geo","Barrel configuration found!"); - uint N_dR = m_numLayers.size(); - double moduleDepth = 0.; - for(uint i_dR = 0; i_dR < N_dR; i_dR++) - { - for(int i_row = 1; i_row <= m_numLayers[i_dR]; i_row++) - { - moduleDepth+=m_dRlayer[i_dR]; - m_radii.push_back(m_offsetR[0] + moduleDepth - m_dRlayer[i_dR]*0.5); - // layer lower and upper edges in z-axis - m_layerEdges.push_back( std::make_pair(m_offsetZ[0] - 0.5*m_widthZ[0], m_offsetZ[0] + 0.5*m_widthZ[0]) ); - m_layerDepth.push_back(m_dRlayer[i_dR]); - } - } - } // Barrel + if(m_detLayout==0) dd4hep::printout(dd4hep::INFO, "FCCSWHCalPhiRow_k4geo","Barrel configuration found!"); + else dd4hep::printout(dd4hep::INFO, "FCCSWHCalPhiRow_k4geo","EndCap configuration found!"); - if(m_offsetZ.size() > 1) // ThreePartsEndCap + // calculate the radius for each layer + uint N_dR = m_numLayers.size()/m_offsetZ.size(); + std::vector moduleDepth(m_offsetZ.size()); + for(uint i_section = 0; i_section < m_offsetZ.size(); i_section++) { - dd4hep::printout(dd4hep::INFO, "FCCSWHCalPhiRow_k4geo","ThreePartsEndCap configuration found!"); - uint N_dR = m_numLayers.size()/3; - double moduleDepth1 = 0.; - double moduleDepth2 = 0.; - double moduleDepth3 = 0.; - // part1 - for(uint i_dR = 0; i_dR < N_dR; i_dR++) - { - for(int i_row = 1; i_row <= m_numLayers[i_dR]; i_row++) - { - moduleDepth1+=m_dRlayer[i_dR]; - m_radii.push_back(m_offsetR[0] + moduleDepth1 - m_dRlayer[i_dR]*0.5); - // layer lower and upper edges in z-axis - m_layerEdges.push_back( std::make_pair(m_offsetZ[0] - 0.5*m_widthZ[0], m_offsetZ[0] + 0.5*m_widthZ[0]) ); - m_layerDepth.push_back(m_dRlayer[i_dR]); - } - } - // part2 for(uint i_dR = 0; i_dR < N_dR; i_dR++) { - for(int i_row = 1; i_row <= m_numLayers[i_dR + N_dR]; i_row++) + for(int i_row = 1; i_row <= m_numLayers[i_dR+i_section*N_dR]; i_row++) { - moduleDepth2+=m_dRlayer[i_dR]; - m_radii.push_back(m_offsetR[1] + moduleDepth2 - m_dRlayer[i_dR]*0.5); + moduleDepth[i_section]+=m_dRlayer[i_dR]; + m_radii.push_back(m_offsetR[i_section] + moduleDepth[i_section] - m_dRlayer[i_dR]*0.5); // layer lower and upper edges in z-axis - m_layerEdges.push_back( std::make_pair(m_offsetZ[1] - 0.5*m_widthZ[1], m_offsetZ[1] + 0.5*m_widthZ[1]) ); + m_layerEdges.push_back( std::make_pair(m_offsetZ[i_section] - 0.5*m_widthZ[i_section], m_offsetZ[i_section] + 0.5*m_widthZ[i_section]) ); m_layerDepth.push_back(m_dRlayer[i_dR]); } } - // part3 - for(uint i_dR = 0; i_dR < N_dR; i_dR++) - { - for(int i_row = 1; i_row <= m_numLayers[i_dR + 2*N_dR]; i_row++) - { - moduleDepth3+=m_dRlayer[i_dR]; - m_radii.push_back(m_offsetR[2] + moduleDepth3 - m_dRlayer[i_dR]*0.5); - // layer lower and upper edges in z-axis - m_layerEdges.push_back( std::make_pair(m_offsetZ[2] - 0.5*m_widthZ[2], m_offsetZ[2] + 0.5*m_widthZ[2]) ); - m_layerDepth.push_back(m_dRlayer[i_dR]); - } - } - } // ThreePartsEndcap + } // print info of calculated radii and edges for(uint i_layer = 0; i_layer < m_radii.size(); i_layer++){ @@ -191,7 +152,7 @@ void FCCSWHCalPhiRow_k4geo::defineCellIndexes(const uint layer) const { } // for the EndCap, do it again but for negative-z part - if(m_offsetZ.size() > 1) + if(m_detLayout == 1) { irow = 0; while( (minLayerZ + (irow+1) * m_dz_row) < (maxLayerZ+0.0001) ) @@ -230,9 +191,9 @@ CellID FCCSWHCalPhiRow_k4geo::cellID(const Vector3D& /* localPosition */, const const VolumeID& vID) const { // get the row number from volumeID (starts from 0!) - int nrow = _decoder->get(vID, "row"); + int nrow = _decoder->get(vID, m_rowID); // get the layer number from volumeID - uint layer = _decoder->get(vID,"layer"); + uint layer = _decoder->get(vID,m_layerID); CellID cID = vID; @@ -240,14 +201,14 @@ CellID FCCSWHCalPhiRow_k4geo::cellID(const Vector3D& /* localPosition */, const int idx = floor(nrow/m_gridSizeRow[layer]) + 1; // if the hit is in the negative-z part of the Endcap then assign negative index - if(m_offsetZ.size() > 1 && globalPosition.z() < 0) idx *= -1; + if(m_detLayout == 1 && globalPosition.z() < 0) idx *= -1; _decoder->set(cID, m_rowID, idx); _decoder->set(cID, m_phiID, positionToBin(dd4hep::DDSegmentation::Util::phiFromXYZ(globalPosition), 2 * M_PI / (double)m_phiBins, m_offsetPhi)); // For endcap, the volume ID comes with "type" field information which would screw up the topo-clustering, // therefore, lets set it to zero, as it is for the cell IDs in the neighbours map. - if(m_offsetZ.size() > 1) _decoder->set(cID, "type", 0); + if(m_detLayout == 1) _decoder->set(cID, "type", 0); return cID; } @@ -266,33 +227,27 @@ std::vector > FCCSWHCalPhiRow_k4geo::getMinMaxLayerId() con if(m_radii.empty()) calculateLayerRadii(); if(m_radii.empty()) return minMaxLayerId; + std::vector minLayerId(m_offsetZ.size(), 0); + std::vector maxLayerId(m_offsetZ.size(), 0); - std::vector minLayerIdEndcap = {0, 0, 0}; - std::vector maxLayerIdEndcap = {-1, 0, 0}; - if(m_offsetZ.size() > 1) - { - uint Nl = m_numLayers.size()/3; + uint Nl = m_numLayers.size()/m_offsetZ.size(); - // count the number of layers in the first part of the Endcap - for(uint i=0; i < Nl; i++) maxLayerIdEndcap[0] += m_numLayers[i]; + for(uint i_section = 0; i_section < m_offsetZ.size(); i_section++) + { + if(i_section > 0) + { + minLayerId[i_section] = maxLayerId[i_section-1] + 1; + maxLayerId[i_section] = maxLayerId[i_section-1]; + } - minLayerIdEndcap[1] = maxLayerIdEndcap[0]+1; - maxLayerIdEndcap[1] = maxLayerIdEndcap[0]; - // count the number of layers in the second part of the Endcap - for(uint i=0; i < Nl; i++) maxLayerIdEndcap[1] += m_numLayers[i+Nl]; + for(uint i=0; i < Nl; i++) + { + maxLayerId[i_section] += m_numLayers[i+i_section*Nl]; + } - minLayerIdEndcap[2] = maxLayerIdEndcap[1]+1; - maxLayerIdEndcap[2] = maxLayerIdEndcap[1]; - // count the number of layers in the third part of the Endcap - for(uint i=0; i < Nl; i++) maxLayerIdEndcap[2] += m_numLayers[i+2*Nl]; + if(i_section==0) maxLayerId[0] -= 1; - minMaxLayerId.push_back(std::make_pair(minLayerIdEndcap[0],maxLayerIdEndcap[0])); - minMaxLayerId.push_back(std::make_pair(minLayerIdEndcap[1],maxLayerIdEndcap[1])); - minMaxLayerId.push_back(std::make_pair(minLayerIdEndcap[2],maxLayerIdEndcap[2])); - } - else // for Barrel - { - minMaxLayerId.push_back(std::make_pair(0,m_radii.size()-1)); + minMaxLayerId.push_back(std::make_pair(minLayerId[i_section], maxLayerId[i_section])); } return minMaxLayerId; @@ -306,14 +261,11 @@ std::vector FCCSWHCalPhiRow_k4geo::neighbours(const CellID& cID) const if(m_radii.empty()) calculateLayerRadii(); if(m_cellIndexes.empty()) return cellNeighbours; - bool EndcapPart1 = false; - bool EndcapPart2 = false; - bool EndcapPart3 = false; - + uint EndcapPart = 0; int minLayerId = -1; int maxLayerId = -1; - int currentLayerId = _decoder->get(cID,"layer"); + int currentLayerId = _decoder->get(cID,m_layerID); int currentCellId = _decoder->get(cID,m_rowID); int minCellId = m_cellIndexes[currentLayerId].front(); @@ -322,54 +274,29 @@ std::vector FCCSWHCalPhiRow_k4geo::neighbours(const CellID& cID) const //-------------------------------- // Determine min and max layer Id //-------------------------------- - // if this is the segmentation of three parts Endcap - /* - * hardcoded numbers would be the following: - * std::vector minLayerIdEndcap = {0, 6, 15}; - * std::vector maxLayerIdEndcap = {5, 14, 36}; - * but lets try to avoid hardcoding: - */ - std::vector minLayerIdEndcap = {0, 0, 0}; - std::vector maxLayerIdEndcap = {0, 0, 0}; - if(m_offsetZ.size() > 1) + std::vector minLayerIdEndcap(m_offsetZ.size(),0); + std::vector maxLayerIdEndcap(m_offsetZ.size(),0); + if(m_detLayout == 1) // Endcap { std::vector > minMaxLayerId(getMinMaxLayerId()); if(minMaxLayerId.empty()) { - dd4hep::printout(dd4hep::ERROR, "FCCSWHCalPhiRow_k4geo","Can not get ThreePartsEndcap min and max layer indexes! --> returning empty neighbours"); + dd4hep::printout(dd4hep::ERROR, "FCCSWHCalPhiRow_k4geo","Can not get Endcap min and max layer indexes! --> returning empty neighbours"); return cellNeighbours; } - // part1 min and max layer index - minLayerIdEndcap[0] = minMaxLayerId[0].first; - maxLayerIdEndcap[0] = minMaxLayerId[0].second; - // part2 min and max layer index - minLayerIdEndcap[1] = minMaxLayerId[1].first; - maxLayerIdEndcap[1] = minMaxLayerId[1].second; - // part3 min and max layer index - minLayerIdEndcap[2] = minMaxLayerId[2].first; - maxLayerIdEndcap[2] = minMaxLayerId[2].second; - - // Part 1 - if(currentLayerId >= minLayerIdEndcap[0] && currentLayerId <= maxLayerIdEndcap[0]) + // determine min and max layer Ids and which section/part it is + for(uint i_section = 0; i_section < minMaxLayerId.size(); i_section++) { - minLayerId = minLayerIdEndcap[0]; - maxLayerId = maxLayerIdEndcap[0]; - EndcapPart1 = true; - } - // Part 2 - if(currentLayerId >= minLayerIdEndcap[1] && currentLayerId <= maxLayerIdEndcap[1]) - { - minLayerId = minLayerIdEndcap[1]; - maxLayerId = maxLayerIdEndcap[1]; - EndcapPart2 = true; - } - // Part 3 - if(currentLayerId >= minLayerIdEndcap[2] && currentLayerId <= maxLayerIdEndcap[2]) - { - minLayerId = minLayerIdEndcap[2]; - maxLayerId = maxLayerIdEndcap[2]; - EndcapPart3 = true; + minLayerIdEndcap[i_section] = minMaxLayerId[i_section].first; + maxLayerIdEndcap[i_section] = minMaxLayerId[i_section].second; + + if(currentLayerId >= minLayerIdEndcap[i_section] && currentLayerId <= maxLayerIdEndcap[i_section]) + { + minLayerId = minLayerIdEndcap[i_section]; + maxLayerId = maxLayerIdEndcap[i_section]; + EndcapPart = i_section; + } } // correct the min and max CellId for endcap @@ -403,7 +330,7 @@ std::vector FCCSWHCalPhiRow_k4geo::neighbours(const CellID& cID) const { CellID nID = cID ; int prevLayerId = currentLayerId - 1; - _decoder->set(nID,"layer",prevLayerId); + _decoder->set(nID,m_layerID,prevLayerId); // if the granularity is the same for the previous layer then take the cells with currentCellId, currentCellId - 1, and currentCellId + 1 if(m_gridSizeRow[prevLayerId] == m_gridSizeRow[currentLayerId]) @@ -450,7 +377,7 @@ std::vector FCCSWHCalPhiRow_k4geo::neighbours(const CellID& cID) const { CellID nID = cID ; int nextLayerId = currentLayerId + 1; - _decoder->set(nID,"layer",nextLayerId); + _decoder->set(nID,m_layerID,nextLayerId); // if the granularity is the same for the next layer then take the cells with currentCellId, currentCellId - 1, and currentCellId + 1 if(m_gridSizeRow[nextLayerId] == m_gridSizeRow[currentLayerId]) @@ -507,8 +434,8 @@ std::vector FCCSWHCalPhiRow_k4geo::neighbours(const CellID& cID) const cellNeighbours.push_back(nID); // add the next cell from current layer of the same phi module } - // if this is the Endcap then look for neighbours in different parts as well - if(m_offsetZ.size() > 1) + // if this is the Endcap (and consists of more than 1 part/section) then look for neighbours in different parts as well + if(m_detLayout == 1 && m_offsetZ.size() > 1) { double currentLayerRmin = m_radii[currentLayerId] - 0.5*m_layerDepth[currentLayerId]; double currentLayerRmax = m_radii[currentLayerId] + 0.5*m_layerDepth[currentLayerId]; @@ -520,8 +447,8 @@ std::vector FCCSWHCalPhiRow_k4geo::neighbours(const CellID& cID) const maxCellId = m_cellIndexes[currentLayerId].back(); // this should be -N } - // if it is the last cell in the part1 - if(EndcapPart1 && currentCellId == maxCellId ) + // if it is the last cell in the first part + if(EndcapPart == 0 && currentCellId == maxCellId ) { // find the layers in the part2 that share a border with the current layer for(int part2layerId = minLayerIdEndcap[1]; part2layerId <= maxLayerIdEndcap[1]; part2layerId++) @@ -536,67 +463,72 @@ std::vector FCCSWHCalPhiRow_k4geo::neighbours(const CellID& cID) const ) { CellID nID = cID ; - _decoder->set(nID,"layer",part2layerId); + _decoder->set(nID,m_layerID,part2layerId); _decoder->set(nID,m_rowID,minCellId); cellNeighbours.push_back(nID); // add the first cell from part2 layer } } } - // if it is the first cell in the part2 - if(EndcapPart2 && currentCellId == minCellId) + // if the Endcap consists of more than 2 parts: + for(uint i_section = 1; i_section < (m_offsetZ.size()-1); i_section++) { - // find the layers in part1 that share a border with the current layer - for(int part1layerId = minLayerIdEndcap[0]; part1layerId <= maxLayerIdEndcap[0]; part1layerId++) - { - double Rmin = m_radii[part1layerId] - 0.5*m_layerDepth[part1layerId]; - double Rmax = m_radii[part1layerId] + 0.5*m_layerDepth[part1layerId]; + if(i_section != EndcapPart) continue; - if( (Rmin >= currentLayerRmin && Rmin <= currentLayerRmax) - || (Rmax >= currentLayerRmin && Rmax <= currentLayerRmax) - || (currentLayerRmin >= Rmin && currentLayerRmin <= Rmax) - || (currentLayerRmax >= Rmin && currentLayerRmax <= Rmax) - ) + // if it is the first cell then look for neighbours in previous part + if(currentCellId == minCellId) + { + // find the layers in previous part that share a border with the current layer + for(int prevPartLayerId = minLayerIdEndcap[i_section-1]; prevPartLayerId <= maxLayerIdEndcap[i_section-1]; prevPartLayerId++) { - CellID nID = cID ; - _decoder->set(nID,"layer",part1layerId); - _decoder->set(nID,m_rowID,maxCellId); - cellNeighbours.push_back(nID); // add the last cell from the part1 layer + double Rmin = m_radii[prevPartLayerId] - 0.5*m_layerDepth[prevPartLayerId]; + double Rmax = m_radii[prevPartLayerId] + 0.5*m_layerDepth[prevPartLayerId]; + + if( (Rmin >= currentLayerRmin && Rmin <= currentLayerRmax) + || (Rmax >= currentLayerRmin && Rmax <= currentLayerRmax) + || (currentLayerRmin >= Rmin && currentLayerRmin <= Rmax) + || (currentLayerRmax >= Rmin && currentLayerRmax <= Rmax) + ) + { + CellID nID = cID ; + _decoder->set(nID,m_layerID,prevPartLayerId); + _decoder->set(nID,m_rowID,maxCellId); + cellNeighbours.push_back(nID); // add the last cell from the previous part layer + } } } - } - - // if it is the last cell in the part2 - if(EndcapPart2 && currentCellId == maxCellId) - { - // find the layers in the part3 that share a border with the current layer - for(int part3layerId = minLayerIdEndcap[2]; part3layerId <= maxLayerIdEndcap[2]; part3layerId++) + // if it is the last cell then look for neighbours in the next part + if(currentCellId == maxCellId) { - double Rmin = m_radii[part3layerId] - 0.5*m_layerDepth[part3layerId]; - double Rmax = m_radii[part3layerId] + 0.5*m_layerDepth[part3layerId]; - - if( (Rmin >= currentLayerRmin && Rmin <= currentLayerRmax) - || (Rmax >= currentLayerRmin && Rmax <= currentLayerRmax) - || (currentLayerRmin >= Rmin && currentLayerRmin <= Rmax) - || (currentLayerRmax >= Rmin && currentLayerRmax <= Rmax) - ) + // find the layers in the next that share a border with the current layer + for(int nextPartLayerId = minLayerIdEndcap[i_section+1]; nextPartLayerId <= maxLayerIdEndcap[i_section+1]; nextPartLayerId++) { - CellID nID = cID ; - _decoder->set(nID,"layer",part3layerId); - _decoder->set(nID,m_rowID,minCellId); - cellNeighbours.push_back(nID); // add the first cell from the part3 layer + double Rmin = m_radii[nextPartLayerId] - 0.5*m_layerDepth[nextPartLayerId]; + double Rmax = m_radii[nextPartLayerId] + 0.5*m_layerDepth[nextPartLayerId]; + + if( (Rmin >= currentLayerRmin && Rmin <= currentLayerRmax) + || (Rmax >= currentLayerRmin && Rmax <= currentLayerRmax) + || (currentLayerRmin >= Rmin && currentLayerRmin <= Rmax) + || (currentLayerRmax >= Rmin && currentLayerRmax <= Rmax) + ) + { + CellID nID = cID ; + _decoder->set(nID,m_layerID,nextPartLayerId); + _decoder->set(nID,m_rowID,minCellId); + cellNeighbours.push_back(nID); // add the first cell from the next part layer + } } } } - // if it is the first cell in the part3 - if(EndcapPart3 && currentCellId == minCellId) + // if it is the first cell in the last part of the Endcap + if(EndcapPart == (m_offsetZ.size() - 1) && currentCellId == minCellId) { - // find the layers in the part2 that share a border with the current layer - for(int part2layerId = minLayerIdEndcap[1]; part2layerId <= maxLayerIdEndcap[1]; part2layerId++) + // find the layers in the previous part that share a border with the current layer + for(int prevPartLayerId = minLayerIdEndcap[m_offsetZ.size()-2]; prevPartLayerId <= maxLayerIdEndcap[m_offsetZ.size()-2]; prevPartLayerId++) { - double Rmin = m_radii[part2layerId] - 0.5*m_layerDepth[part2layerId]; - double Rmax = m_radii[part2layerId] + 0.5*m_layerDepth[part2layerId]; + double Rmin = m_radii[prevPartLayerId] - 0.5*m_layerDepth[prevPartLayerId]; + double Rmax = m_radii[prevPartLayerId] + 0.5*m_layerDepth[prevPartLayerId]; if( (Rmin >= currentLayerRmin && Rmin <= currentLayerRmax) || (Rmax >= currentLayerRmin && Rmax <= currentLayerRmax) @@ -605,7 +537,7 @@ std::vector FCCSWHCalPhiRow_k4geo::neighbours(const CellID& cID) const ) { CellID nID = cID ; - _decoder->set(nID,"layer",part2layerId); + _decoder->set(nID,m_layerID,prevPartLayerId); _decoder->set(nID,m_rowID,maxCellId); cellNeighbours.push_back(nID); // add the last cell from the part2 layer } @@ -645,7 +577,7 @@ std::array FCCSWHCalPhiRow_k4geo::cellTheta(const CellID& cID) const // get the cell index int idx = _decoder->get(cID, m_rowID); // get the layer index - uint layer = _decoder->get(cID,"layer"); + uint layer = _decoder->get(cID,m_layerID); if(m_radii.empty()) calculateLayerRadii(); if(m_cellEdges.empty()) return cTheta; diff --git a/detectorSegmentations/src/FCCSWHCalPhiTheta_k4geo.cpp b/detectorSegmentations/src/FCCSWHCalPhiTheta_k4geo.cpp index 04c227ab1..a6248a3d3 100644 --- a/detectorSegmentations/src/FCCSWHCalPhiTheta_k4geo.cpp +++ b/detectorSegmentations/src/FCCSWHCalPhiTheta_k4geo.cpp @@ -13,12 +13,14 @@ FCCSWHCalPhiTheta_k4geo::FCCSWHCalPhiTheta_k4geo(const std::string& cellEncoding // register all necessary parameters (additional to those registered in GridTheta_k4geo) registerParameter("phi_bins", "Number of bins phi", m_phiBins, 1); registerParameter("offset_phi", "Angular offset in phi", m_offsetPhi, 0., SegmentationParameter::AngleUnit, true); + registerParameter("detLayout", "The detector layout (0 = Barrel; 1 = Endcap)", m_detLayout, -1); registerParameter("offset_z", "Offset in z-axis of the layer center", m_offsetZ, std::vector()); registerParameter("width_z", "Width in z of the layer", m_widthZ, std::vector()); registerParameter("offset_r", "Offset in radius of the layer (Rmin)", m_offsetR, std::vector()); registerParameter("numLayers", "Number of layers", m_numLayers, std::vector()); registerParameter("dRlayer", "dR of the layer", m_dRlayer, std::vector()); registerIdentifier("identifier_phi", "Cell ID identifier for phi", m_phiID, "phi"); + registerIdentifier("identifier_layer", "Cell ID identifier for layer", m_layerID, "layer"); } FCCSWHCalPhiTheta_k4geo::FCCSWHCalPhiTheta_k4geo(const BitFieldCoder* decoder) : GridTheta_k4geo(decoder) { @@ -29,18 +31,20 @@ FCCSWHCalPhiTheta_k4geo::FCCSWHCalPhiTheta_k4geo(const BitFieldCoder* decoder) : // register all necessary parameters (additional to those registered in GridTheta_k4geo) registerParameter("phi_bins", "Number of bins phi", m_phiBins, 1); registerParameter("offset_phi", "Angular offset in phi", m_offsetPhi, 0., SegmentationParameter::AngleUnit, true); + registerParameter("detLayout", "The detector layout (0 = Barrel; 1 = Endcap)", m_detLayout, -1); registerParameter("offset_z", "Offset in z-axis of the layer center", m_offsetZ, std::vector()); registerParameter("width_z", "Width in z of the layer", m_widthZ, std::vector()); registerParameter("offset_r", "Offset in radius of the layer (Rmin)", m_offsetR, std::vector()); registerParameter("numLayers", "Number of layers", m_numLayers, std::vector()); registerParameter("dRlayer", "dR of the layer", m_dRlayer, std::vector()); registerIdentifier("identifier_phi", "Cell ID identifier for phi", m_phiID, "phi"); + registerIdentifier("identifier_layer", "Cell ID identifier for layer", m_layerID, "layer"); } /** /// determine the global position based on the cell ID Vector3D FCCSWHCalPhiTheta_k4geo::position(const CellID& cID) const { - uint layer = _decoder->get(cID,"layer"); + uint layer = _decoder->get(cID,m_layerID); double radius = 1.0; if(m_radii.empty()) defineCellsInRZplan(); @@ -53,8 +57,8 @@ Vector3D FCCSWHCalPhiTheta_k4geo::position(const CellID& cID) const { /// determine the global position based on the cell ID /// returns the geometric center of the cell Vector3D FCCSWHCalPhiTheta_k4geo::position(const CellID& cID) const { - uint layer = _decoder->get(cID,"layer"); - int thetaID = _decoder->get(cID,"theta"); + uint layer = _decoder->get(cID,m_layerID); + int thetaID = _decoder->get(cID,m_thetaID); double zpos = 0.; double radius = 1.0; @@ -73,76 +77,33 @@ void FCCSWHCalPhiTheta_k4geo::defineCellsInRZplan() const { if(m_radii.empty()) { // check if all necessary variables are available - if(m_offsetZ.empty() || m_widthZ.empty() || m_offsetR.empty() || m_numLayers.empty() || m_dRlayer.empty()) + if(m_detLayout==-1 || m_offsetZ.empty() || m_widthZ.empty() || m_offsetR.empty() || m_numLayers.empty() || m_dRlayer.empty()) { - dd4hep::printout(dd4hep::ERROR, "FCCSWHCalPhiTheta_k4geo","Please check the readout description in the XML file!", - "One of the variables is missing: offset_z | width_z | offset_r | numLayers | dRlayer"); + dd4hep::printout(dd4hep::ERROR, "FCCSWHCalPhiRow_k4geo","Please check the readout description in the XML file!", + "One of the variables is missing: detLayout | offset_z | width_z | offset_r | numLayers | dRlayer"); return; } - // calculate the radius for each layer - if(m_offsetZ.size() == 1) // Barrel - { - dd4hep::printout(dd4hep::INFO, "FCCSWHCalPhiTheta_k4geo","Barrel configuration found!"); - uint N_dR = m_numLayers.size(); - double moduleDepth = 0.; - for(uint i_dR = 0; i_dR < N_dR; i_dR++) - { - for(int i_row = 1; i_row <= m_numLayers[i_dR]; i_row++) - { - moduleDepth+=m_dRlayer[i_dR]; - m_radii.push_back(m_offsetR[0] + moduleDepth - m_dRlayer[i_dR]*0.5); - // layer lower and upper edges in z-axis - m_layerEdges.push_back( std::make_pair(m_offsetZ[0] - 0.5*m_widthZ[0], m_offsetZ[0] + 0.5*m_widthZ[0]) ); - m_layerDepth.push_back(m_dRlayer[i_dR]); - } - } - } // Barrel + if(m_detLayout==0) dd4hep::printout(dd4hep::INFO, "FCCSWHCalPhiRow_k4geo","Barrel configuration found!"); + else dd4hep::printout(dd4hep::INFO, "FCCSWHCalPhiRow_k4geo","EndCap configuration found!"); - if(m_offsetZ.size() > 1) // ThreePartsEndCap + // calculate the radius for each layer + uint N_dR = m_numLayers.size()/m_offsetZ.size(); + std::vector moduleDepth(m_offsetZ.size()); + for(uint i_section = 0; i_section < m_offsetZ.size(); i_section++) { - dd4hep::printout(dd4hep::INFO, "FCCSWHCalPhiTheta_k4geo","ThreePartsEndCap configuration found!"); - uint N_dR = m_numLayers.size()/3; - double moduleDepth1 = 0.; - double moduleDepth2 = 0.; - double moduleDepth3 = 0.; - // part1 for(uint i_dR = 0; i_dR < N_dR; i_dR++) { - for(int i_row = 1; i_row <= m_numLayers[i_dR]; i_row++) + for(int i_row = 1; i_row <= m_numLayers[i_dR+i_section*N_dR]; i_row++) { - moduleDepth1+=m_dRlayer[i_dR]; - m_radii.push_back(m_offsetR[0] + moduleDepth1 - m_dRlayer[i_dR]*0.5); + moduleDepth[i_section]+=m_dRlayer[i_dR]; + m_radii.push_back(m_offsetR[i_section] + moduleDepth[i_section] - m_dRlayer[i_dR]*0.5); // layer lower and upper edges in z-axis - m_layerEdges.push_back( std::make_pair(m_offsetZ[0] - 0.5*m_widthZ[0], m_offsetZ[0] + 0.5*m_widthZ[0]) ); + m_layerEdges.push_back( std::make_pair(m_offsetZ[i_section] - 0.5*m_widthZ[i_section], m_offsetZ[i_section] + 0.5*m_widthZ[i_section]) ); m_layerDepth.push_back(m_dRlayer[i_dR]); } } - // part2 - for(uint i_dR = 0; i_dR < N_dR; i_dR++) - { - for(int i_row = 1; i_row <= m_numLayers[i_dR + N_dR]; i_row++) - { - moduleDepth2+=m_dRlayer[i_dR]; - m_radii.push_back(m_offsetR[1] + moduleDepth2 - m_dRlayer[i_dR]*0.5); - // layer lower and upper edges in z-axis - m_layerEdges.push_back( std::make_pair(m_offsetZ[1] - 0.5*m_widthZ[1], m_offsetZ[1] + 0.5*m_widthZ[1]) ); - m_layerDepth.push_back(m_dRlayer[i_dR]); - } - } - // part3 - for(uint i_dR = 0; i_dR < N_dR; i_dR++) - { - for(int i_row = 1; i_row <= m_numLayers[i_dR + 2*N_dR]; i_row++) - { - moduleDepth3+=m_dRlayer[i_dR]; - m_radii.push_back(m_offsetR[2] + moduleDepth3 - m_dRlayer[i_dR]*0.5); - // layer lower and upper edges in z-axis - m_layerEdges.push_back( std::make_pair(m_offsetZ[2] - 0.5*m_widthZ[2], m_offsetZ[2] + 0.5*m_widthZ[2]) ); - m_layerDepth.push_back(m_dRlayer[i_dR]); - } - } - } // ThreePartsEndcap + } // print info of calculated radii and edges for(uint i_layer = 0; i_layer < m_radii.size(); i_layer++){ @@ -196,7 +157,7 @@ void FCCSWHCalPhiTheta_k4geo::defineCellEdges(const uint layer) const { m_cellEdges[layer][prevBin].first = m_layerEdges[layer].first; // for the EndCap, do it again but for negative z part - if(m_offsetZ.size() > 1) + if(m_detLayout == 1) { while(m_radii[layer]*std::cos(m_offsetTheta+ibin*m_gridSizeTheta)/std::sin(m_offsetTheta+ibin*m_gridSizeTheta) > (-m_layerEdges[layer].second)) { @@ -248,11 +209,11 @@ CellID FCCSWHCalPhiTheta_k4geo::cellID(const Vector3D& /* localPosition */, cons // For endcap, the volume ID comes with "type" field information which would screw up the topo-clustering as the "row" field, // therefore, lets set it to zero, as it is for the cell IDs in the neighbours map. - if(m_offsetZ.size() > 1) _decoder->set(cID, "type", 0); + if(m_detLayout == 1) _decoder->set(cID, "type", 0); double lTheta = thetaFromXYZ(globalPosition); double lPhi = phiFromXYZ(globalPosition); - uint layer = _decoder->get(vID,"layer"); + uint layer = _decoder->get(vID,m_layerID); // define cell boundaries in R-z plan if(m_radii.empty()) defineCellsInRZplan(); @@ -288,45 +249,32 @@ double FCCSWHCalPhiTheta_k4geo::phi(const CellID& cID) const { /// Get the min and max layer indexes for each part of the HCal std::vector > FCCSWHCalPhiTheta_k4geo::getMinMaxLayerId() const { - /* - * hardcoded numbers would be the following: - * std::vector minLayerIdEndcap = {0, 6, 15}; - * std::vector maxLayerIdEndcap = {5, 14, 36}; - * but lets try to avoid hardcoding: - */ - std::vector > minMaxLayerId; if(m_radii.empty()) defineCellsInRZplan(); if(m_radii.empty()) return minMaxLayerId; + std::vector minLayerId(m_offsetZ.size(), 0); + std::vector maxLayerId(m_offsetZ.size(), 0); - std::vector minLayerIdEndcap = {0, 0, 0}; - std::vector maxLayerIdEndcap = {-1, 0, 0}; - if(m_offsetZ.size() > 1) - { - uint Nl = m_numLayers.size()/3; + uint Nl = m_numLayers.size()/m_offsetZ.size(); - // count the number of layers in the first part of the Endcap - for(uint i=0; i < Nl; i++) maxLayerIdEndcap[0] += m_numLayers[i]; + for(uint i_section = 0; i_section < m_offsetZ.size(); i_section++) + { + if(i_section > 0) + { + minLayerId[i_section] = maxLayerId[i_section-1] + 1; + maxLayerId[i_section] = maxLayerId[i_section-1]; + } - minLayerIdEndcap[1] = maxLayerIdEndcap[0]+1; - maxLayerIdEndcap[1] = maxLayerIdEndcap[0]; - // count the number of layers in the second part of the Endcap - for(uint i=0; i < Nl; i++) maxLayerIdEndcap[1] += m_numLayers[i+Nl]; + for(uint i=0; i < Nl; i++) + { + maxLayerId[i_section] += m_numLayers[i+i_section*Nl]; + } - minLayerIdEndcap[2] = maxLayerIdEndcap[1]+1; - maxLayerIdEndcap[2] = maxLayerIdEndcap[1]; - // count the number of layers in the third part of the Endcap - for(uint i=0; i < Nl; i++) maxLayerIdEndcap[2] += m_numLayers[i+2*Nl]; + if(i_section == 0) maxLayerId[0] -= 1; - minMaxLayerId.push_back(std::make_pair(minLayerIdEndcap[0],maxLayerIdEndcap[0])); - minMaxLayerId.push_back(std::make_pair(minLayerIdEndcap[1],maxLayerIdEndcap[1])); - minMaxLayerId.push_back(std::make_pair(minLayerIdEndcap[2],maxLayerIdEndcap[2])); - } - else // for Barrel - { - minMaxLayerId.push_back(std::make_pair(0,m_radii.size()-1)); + minMaxLayerId.push_back(std::make_pair(minLayerId[i_section], maxLayerId[i_section])); } return minMaxLayerId; @@ -340,14 +288,11 @@ std::vector FCCSWHCalPhiTheta_k4geo::neighbours(const CellID& cID, boo if(m_radii.empty()) defineCellsInRZplan(); if(m_thetaBins.empty()) return cellNeighbours; - bool EndcapPart1 = false; - bool EndcapPart2 = false; - bool EndcapPart3 = false; - + uint EndcapPart = 0; int minLayerId = -1; int maxLayerId = -1; - int currentLayerId = _decoder->get(cID,"layer"); + int currentLayerId = _decoder->get(cID,m_layerID); int currentCellThetaBin = _decoder->get(cID,m_thetaID); int minCellThetaBin = m_thetaBins[currentLayerId].front(); @@ -356,48 +301,29 @@ std::vector FCCSWHCalPhiTheta_k4geo::neighbours(const CellID& cID, boo //-------------------------------- // Determine min and max layer Id //-------------------------------- - // if this is the segmentation of three parts Endcap - std::vector minLayerIdEndcap = {0, 0, 0}; - std::vector maxLayerIdEndcap = {0, 0, 0}; - if(m_offsetZ.size() > 1) + std::vector minLayerIdEndcap(m_offsetZ.size(),0); + std::vector maxLayerIdEndcap(m_offsetZ.size(),0); + if(m_detLayout == 1) { std::vector > minMaxLayerId(getMinMaxLayerId()); if(minMaxLayerId.empty()) { - dd4hep::printout(dd4hep::ERROR, "FCCSWHCalPhiTheta_k4geo","Can not get ThreePartsEndcap min and max layer indexes! --> returning empty neighbours"); + dd4hep::printout(dd4hep::ERROR, "FCCSWHCalPhiTheta_k4geo","Can not get Endcap min and max layer indexes! --> returning empty neighbours"); return cellNeighbours; } - // part1 min and max layer index - minLayerIdEndcap[0] = minMaxLayerId[0].first; - maxLayerIdEndcap[0] = minMaxLayerId[0].second; - // part2 min and max layer index - minLayerIdEndcap[1] = minMaxLayerId[1].first; - maxLayerIdEndcap[1] = minMaxLayerId[1].second; - // part3 min and max layer index - minLayerIdEndcap[2] = minMaxLayerId[2].first; - maxLayerIdEndcap[2] = minMaxLayerId[2].second; - - // Part 1 - if(currentLayerId >= minLayerIdEndcap[0] && currentLayerId <= maxLayerIdEndcap[0]) + // determine min and max layer Ids and which section/part it is + for(uint i_section = 0; i_section < minMaxLayerId.size(); i_section++) { - minLayerId = minLayerIdEndcap[0]; - maxLayerId = maxLayerIdEndcap[0]; - EndcapPart1 = true; - } - // Part 2 - if(currentLayerId >= minLayerIdEndcap[1] && currentLayerId <= maxLayerIdEndcap[1]) - { - minLayerId = minLayerIdEndcap[1]; - maxLayerId = maxLayerIdEndcap[1]; - EndcapPart2 = true; - } - // Part 3 - if(currentLayerId >= minLayerIdEndcap[2] && currentLayerId <= maxLayerIdEndcap[2]) - { - minLayerId = minLayerIdEndcap[2]; - maxLayerId = maxLayerIdEndcap[2]; - EndcapPart3 = true; + minLayerIdEndcap[i_section] = minMaxLayerId[i_section].first; + maxLayerIdEndcap[i_section] = minMaxLayerId[i_section].second; + + if(currentLayerId >= minLayerIdEndcap[i_section] && currentLayerId <= maxLayerIdEndcap[i_section]) + { + minLayerId = minLayerIdEndcap[i_section]; + maxLayerId = maxLayerIdEndcap[i_section]; + EndcapPart = i_section; + } } // correct the min and max theta bin for endcap @@ -446,7 +372,7 @@ std::vector FCCSWHCalPhiTheta_k4geo::neighbours(const CellID& cID, boo //---------------------------------------------- // deal with the Barrel - if(m_offsetZ.size() == 1) + if(m_detLayout == 0) { double currentCellZmin = m_cellEdges[currentLayerId][currentCellThetaBin].first; double currentCellZmax = m_cellEdges[currentLayerId][currentCellThetaBin].second; @@ -456,7 +382,7 @@ std::vector FCCSWHCalPhiTheta_k4geo::neighbours(const CellID& cID, boo { CellID nID = cID ; int prevLayerId = currentLayerId - 1; - _decoder->set(nID,"layer",prevLayerId); + _decoder->set(nID,m_layerID,prevLayerId); _decoder->set(nID,m_thetaID,currentCellThetaBin); cellNeighbours.push_back(nID); // add the cell with the same theta bin from the previous layer of the same phi module @@ -502,7 +428,7 @@ std::vector FCCSWHCalPhiTheta_k4geo::neighbours(const CellID& cID, boo { CellID nID = cID ; int nextLayerId = currentLayerId + 1; - _decoder->set(nID,"layer",nextLayerId); + _decoder->set(nID,m_layerID,nextLayerId); _decoder->set(nID,m_thetaID,currentCellThetaBin); cellNeighbours.push_back(nID); // add the cell with the same theta bin from the next layer of the same phi module @@ -548,8 +474,8 @@ std::vector FCCSWHCalPhiTheta_k4geo::neighbours(const CellID& cID, boo } } - // if this is the Endcap then look for neighbours in different parts as well - if(m_offsetZ.size() > 1) + // Endcap + if(m_detLayout == 1) { double currentCellZmin = m_cellEdges[currentLayerId][currentCellThetaBin].first; double currentCellZmax = m_cellEdges[currentLayerId][currentCellThetaBin].second; @@ -559,7 +485,7 @@ std::vector FCCSWHCalPhiTheta_k4geo::neighbours(const CellID& cID, boo { CellID nID = cID ; int prevLayerId = currentLayerId - 1; - _decoder->set(nID,"layer",prevLayerId); + _decoder->set(nID,m_layerID,prevLayerId); // find the ones that share at least part of a border with the current cell for( auto bin : m_thetaBins[prevLayerId] ) { @@ -607,7 +533,7 @@ std::vector FCCSWHCalPhiTheta_k4geo::neighbours(const CellID& cID, boo { CellID nID = cID ; int nextLayerId = currentLayerId + 1; - _decoder->set(nID,"layer",nextLayerId); + _decoder->set(nID,m_layerID,nextLayerId); // find the ones that share at least part of a border with the current cell for( auto bin : m_thetaBins[nextLayerId] ) { @@ -650,135 +576,143 @@ std::vector FCCSWHCalPhiTheta_k4geo::neighbours(const CellID& cID, boo } } - - // - double currentLayerRmin = m_radii[currentLayerId] - 0.5*m_layerDepth[currentLayerId]; - double currentLayerRmax = m_radii[currentLayerId] + 0.5*m_layerDepth[currentLayerId]; - - - // if the cell is in negative-z part, then swap min and max theta bins - if(theta(cID) > M_PI/2.) + // if the Endcap consists of more than 1 part/section then look for neighbours in different parts as well + if(m_offsetZ.size() > 1) { - minCellThetaBin = m_thetaBins[currentLayerId].back(); - maxCellThetaBin = m_thetaBins[currentLayerId][m_thetaBins[currentLayerId].size()/2]; - } + // + double currentLayerRmin = m_radii[currentLayerId] - 0.5*m_layerDepth[currentLayerId]; + double currentLayerRmax = m_radii[currentLayerId] + 0.5*m_layerDepth[currentLayerId]; - // if it is the last cell in the part1 - if(EndcapPart1 && currentCellThetaBin == minCellThetaBin ) - { - // find the layers in the part2 that share a border with the current layer - for(int part2layerId = minLayerIdEndcap[1]; part2layerId <= maxLayerIdEndcap[1]; part2layerId++) + // if the cell is in negative-z part, then swap min and max theta bins + if(theta(cID) > M_PI/2.) { - double Rmin = m_radii[part2layerId] - 0.5*m_layerDepth[part2layerId]; - double Rmax = m_radii[part2layerId] + 0.5*m_layerDepth[part2layerId]; - - if( (Rmin >= currentLayerRmin && Rmin <= currentLayerRmax) - || (Rmax > currentLayerRmin && Rmax <= currentLayerRmax) - || (currentLayerRmin >= Rmin && currentLayerRmin < Rmax) - || (currentLayerRmax >= Rmin && currentLayerRmax <= Rmax) - ) - { - CellID nID = cID ; - _decoder->set(nID,"layer",part2layerId); - _decoder->set(nID,m_thetaID,maxCellThetaBin); - cellNeighbours.push_back(nID); // add the last theta bin cell from part2 layer - } - if(aDiagonal && Rmax == currentLayerRmin) - { - CellID nID = cID ; - _decoder->set(nID,"layer",part2layerId); - _decoder->set(nID,m_thetaID,maxCellThetaBin); - cellNeighbours.push_back(nID); // add the last theta bin cell from part2 layer - } + minCellThetaBin = m_thetaBins[currentLayerId].back(); + maxCellThetaBin = m_thetaBins[currentLayerId][m_thetaBins[currentLayerId].size()/2]; } - } - // if it is the last theta bin cell in the part2 - if(EndcapPart2 && currentCellThetaBin == maxCellThetaBin) - { - // find the layers in part1 that share a border with the current layer - for(int part1layerId = minLayerIdEndcap[0]; part1layerId <= maxLayerIdEndcap[0]; part1layerId++) + // if it is the last cell in the part1 + if(EndcapPart == 0 && currentCellThetaBin == minCellThetaBin ) { - double Rmin = m_radii[part1layerId] - 0.5*m_layerDepth[part1layerId]; - double Rmax = m_radii[part1layerId] + 0.5*m_layerDepth[part1layerId]; - - if( (Rmin >= currentLayerRmin && Rmin < currentLayerRmax) - || (Rmax >= currentLayerRmin && Rmax <= currentLayerRmax) - || (currentLayerRmin >= Rmin && currentLayerRmin <= Rmax) - || (currentLayerRmax > Rmin && currentLayerRmax <= Rmax) - ) - { - CellID nID = cID ; - _decoder->set(nID,"layer",part1layerId); - _decoder->set(nID,m_thetaID,minCellThetaBin); - cellNeighbours.push_back(nID); // add the first theta bin cell from the part1 layer - } - if(aDiagonal && Rmin == currentLayerRmax) + // find the layers in the part2 that share a border with the current layer + for(int part2layerId = minLayerIdEndcap[1]; part2layerId <= maxLayerIdEndcap[1]; part2layerId++) { - CellID nID = cID ; - _decoder->set(nID,"layer",part1layerId); - _decoder->set(nID,m_thetaID,minCellThetaBin); - cellNeighbours.push_back(nID); // add the first theta bin cell from the part1 layer + double Rmin = m_radii[part2layerId] - 0.5*m_layerDepth[part2layerId]; + double Rmax = m_radii[part2layerId] + 0.5*m_layerDepth[part2layerId]; + + if( (Rmin >= currentLayerRmin && Rmin <= currentLayerRmax) + || (Rmax > currentLayerRmin && Rmax <= currentLayerRmax) + || (currentLayerRmin >= Rmin && currentLayerRmin < Rmax) + || (currentLayerRmax >= Rmin && currentLayerRmax <= Rmax) + ) + { + CellID nID = cID ; + _decoder->set(nID,m_layerID,part2layerId); + _decoder->set(nID,m_thetaID,maxCellThetaBin); + cellNeighbours.push_back(nID); // add the last theta bin cell from part2 layer + } + if(aDiagonal && Rmax == currentLayerRmin) + { + CellID nID = cID ; + _decoder->set(nID,m_layerID,part2layerId); + _decoder->set(nID,m_thetaID,maxCellThetaBin); + cellNeighbours.push_back(nID); // add the last theta bin cell from part2 layer + } } } - } - // if it is the first theta bin cell in the part2 - if(EndcapPart2 && currentCellThetaBin == minCellThetaBin) - { - // find the layers in the part3 that share a border with the current layer - for(int part3layerId = minLayerIdEndcap[2]; part3layerId <= maxLayerIdEndcap[2]; part3layerId++) + // if the Endcap consists of more than 2 parts: + for(uint i_section = 1; i_section < (m_offsetZ.size()-1); i_section++) { - double Rmin = m_radii[part3layerId] - 0.5*m_layerDepth[part3layerId]; - double Rmax = m_radii[part3layerId] + 0.5*m_layerDepth[part3layerId]; - - if( (Rmin >= currentLayerRmin && Rmin <= currentLayerRmax) - || (Rmax > currentLayerRmin && Rmax <= currentLayerRmax) - || (currentLayerRmin >= Rmin && currentLayerRmin < Rmax) - || (currentLayerRmax >= Rmin && currentLayerRmax <= Rmax) - ) + if(i_section != EndcapPart) continue; + + // if it is the last theta bin cell then look for neighbours in previous part + if(currentCellThetaBin == maxCellThetaBin) { - CellID nID = cID ; - _decoder->set(nID,"layer",part3layerId); - _decoder->set(nID,m_thetaID,maxCellThetaBin); - cellNeighbours.push_back(nID); // add the first cell from the part3 layer + // find the layers in the previous part that share a border with the current layer + for(int prevPartLayerId = minLayerIdEndcap[i_section-1]; prevPartLayerId <= maxLayerIdEndcap[i_section-1]; prevPartLayerId++) + { + double Rmin = m_radii[prevPartLayerId] - 0.5*m_layerDepth[prevPartLayerId]; + double Rmax = m_radii[prevPartLayerId] + 0.5*m_layerDepth[prevPartLayerId]; + + if( (Rmin >= currentLayerRmin && Rmin < currentLayerRmax) + || (Rmax >= currentLayerRmin && Rmax <= currentLayerRmax) + || (currentLayerRmin >= Rmin && currentLayerRmin <= Rmax) + || (currentLayerRmax > Rmin && currentLayerRmax <= Rmax) + ) + { + CellID nID = cID ; + _decoder->set(nID,m_layerID,prevPartLayerId); + _decoder->set(nID,m_thetaID,minCellThetaBin); + cellNeighbours.push_back(nID); // add the first theta bin cell from the part1 layer + } + if(aDiagonal && Rmin == currentLayerRmax) + { + CellID nID = cID ; + _decoder->set(nID,m_layerID,prevPartLayerId); + _decoder->set(nID,m_thetaID,minCellThetaBin); + cellNeighbours.push_back(nID); // add the first theta bin cell from the part1 layer + } + } } - if(aDiagonal && Rmax == currentLayerRmin) + + // if it is the first theta bin cell then look for neighbours in the next part + if(currentCellThetaBin == minCellThetaBin) { - CellID nID = cID ; - _decoder->set(nID,"layer",part3layerId); - _decoder->set(nID,m_thetaID,maxCellThetaBin); - cellNeighbours.push_back(nID); // add the first cell from the part3 layer + // find the layers in the next part that share a border with the current layer + for(int nextPartLayerId = minLayerIdEndcap[i_section+1]; nextPartLayerId <= maxLayerIdEndcap[i_section+1]; nextPartLayerId++) + { + double Rmin = m_radii[nextPartLayerId] - 0.5*m_layerDepth[nextPartLayerId]; + double Rmax = m_radii[nextPartLayerId] + 0.5*m_layerDepth[nextPartLayerId]; + + if( (Rmin >= currentLayerRmin && Rmin <= currentLayerRmax) + || (Rmax > currentLayerRmin && Rmax <= currentLayerRmax) + || (currentLayerRmin >= Rmin && currentLayerRmin < Rmax) + || (currentLayerRmax >= Rmin && currentLayerRmax <= Rmax) + ) + { + CellID nID = cID ; + _decoder->set(nID,m_layerID,nextPartLayerId); + _decoder->set(nID,m_thetaID,maxCellThetaBin); + cellNeighbours.push_back(nID); // add the first cell from the part3 layer + } + if(aDiagonal && Rmax == currentLayerRmin) + { + CellID nID = cID ; + _decoder->set(nID,m_layerID,nextPartLayerId); + _decoder->set(nID,m_thetaID,maxCellThetaBin); + cellNeighbours.push_back(nID); // add the first cell from the part3 layer + } + } } } - } - // if it is the last theta bin cell in the part3 - if(EndcapPart3 && currentCellThetaBin == maxCellThetaBin) - { - // find the layers in the part2 that share a border with the current layer - for(int part2layerId = minLayerIdEndcap[1]; part2layerId <= maxLayerIdEndcap[1]; part2layerId++) + // if it is the last theta bin cell in the last part of the Endcap + if(EndcapPart == (m_offsetZ.size() - 1) && currentCellThetaBin == maxCellThetaBin) { - double Rmin = m_radii[part2layerId] - 0.5*m_layerDepth[part2layerId]; - double Rmax = m_radii[part2layerId] + 0.5*m_layerDepth[part2layerId]; - - if( (Rmin >= currentLayerRmin && Rmin < currentLayerRmax) - || (Rmax >= currentLayerRmin && Rmax <= currentLayerRmax) - || (currentLayerRmin >= Rmin && currentLayerRmin <= Rmax) - || (currentLayerRmax > Rmin && currentLayerRmax <= Rmax) - ) + // find the layers in the previous part that share a border with the current layer + for(int prevPartLayerId = minLayerIdEndcap[m_offsetZ.size()-2]; prevPartLayerId <= maxLayerIdEndcap[m_offsetZ.size()-2]; prevPartLayerId++) { - CellID nID = cID ; - _decoder->set(nID,"layer",part2layerId); - _decoder->set(nID,m_thetaID,minCellThetaBin); - cellNeighbours.push_back(nID); // add the first theta bin cell from the part2 layer - } - if(aDiagonal && Rmin == currentLayerRmax) - { - CellID nID = cID ; - _decoder->set(nID,"layer",part2layerId); - _decoder->set(nID,m_thetaID,minCellThetaBin); - cellNeighbours.push_back(nID); // add the first theta bin cell from the part2 layer + double Rmin = m_radii[prevPartLayerId] - 0.5*m_layerDepth[prevPartLayerId]; + double Rmax = m_radii[prevPartLayerId] + 0.5*m_layerDepth[prevPartLayerId]; + + if( (Rmin >= currentLayerRmin && Rmin < currentLayerRmax) + || (Rmax >= currentLayerRmin && Rmax <= currentLayerRmax) + || (currentLayerRmin >= Rmin && currentLayerRmin <= Rmax) + || (currentLayerRmax > Rmin && currentLayerRmax <= Rmax) + ) + { + CellID nID = cID ; + _decoder->set(nID,m_layerID,prevPartLayerId); + _decoder->set(nID,m_thetaID,minCellThetaBin); + cellNeighbours.push_back(nID); // add the first theta bin cell from the part2 layer + } + if(aDiagonal && Rmin == currentLayerRmax) + { + CellID nID = cID ; + _decoder->set(nID,m_layerID,prevPartLayerId); + _decoder->set(nID,m_thetaID,minCellThetaBin); + cellNeighbours.push_back(nID); // add the first theta bin cell from the part2 layer + } } } } @@ -816,7 +750,7 @@ std::array FCCSWHCalPhiTheta_k4geo::cellTheta(const CellID& cID) cons // get the cell index int idx = _decoder->get(cID, m_thetaID); // get the layer index - uint layer = _decoder->get(cID,"layer"); + uint layer = _decoder->get(cID,m_layerID); if(m_radii.empty()) defineCellsInRZplan(); if(m_cellEdges.empty()) return cTheta; From e752b14f44e0762b51e476a3e9158a7cf0ce849b Mon Sep 17 00:00:00 2001 From: Archil Durglishvili Date: Sat, 30 Nov 2024 21:27:43 +0100 Subject: [PATCH 14/17] Add the method thetaMax() for HCal phi-row segmentation needed for SW clustering --- .../FCCSWHCalPhiRow_k4geo.h | 7 ++++- .../src/FCCSWHCalPhiRow_k4geo.cpp | 27 +++++++++++++++++++ 2 files changed, 33 insertions(+), 1 deletion(-) diff --git a/detectorSegmentations/include/detectorSegmentations/FCCSWHCalPhiRow_k4geo.h b/detectorSegmentations/include/detectorSegmentations/FCCSWHCalPhiRow_k4geo.h index acda782e8..1df1353d9 100644 --- a/detectorSegmentations/include/detectorSegmentations/FCCSWHCalPhiRow_k4geo.h +++ b/detectorSegmentations/include/detectorSegmentations/FCCSWHCalPhiRow_k4geo.h @@ -111,7 +111,7 @@ namespace dd4hep { */ std::array cellTheta(const CellID& cID) const; - /** Get the vector of cell indexes in a give layer. + /** Get the vector of cell indexes in a given layer. */ inline std::vector cellIndexes(const uint layer) const { if(m_radii.empty()) calculateLayerRadii(); @@ -119,6 +119,11 @@ namespace dd4hep { else return std::vector(); } + /** Get the thetaMax needed for SW clustering + * return max theta value of the detector + */ + double thetaMax() const; + /** Get the min and max layer indexes of each HCal part. * For Endcap, returns the three elements vector, while for Barrel - single element vector. */ diff --git a/detectorSegmentations/src/FCCSWHCalPhiRow_k4geo.cpp b/detectorSegmentations/src/FCCSWHCalPhiRow_k4geo.cpp index f17059abb..09d644982 100644 --- a/detectorSegmentations/src/FCCSWHCalPhiRow_k4geo.cpp +++ b/detectorSegmentations/src/FCCSWHCalPhiRow_k4geo.cpp @@ -602,5 +602,32 @@ std::array FCCSWHCalPhiRow_k4geo::cellTheta(const CellID& cID) const return cTheta; } +/// determine maximum theta value of the detector. This is used by SW clustering +double FCCSWHCalPhiRow_k4geo::thetaMax() const { + std::vector > minMaxLayerId(getMinMaxLayerId()); + if(minMaxLayerId.empty()) return 0.; + + // get the first layerId in the Barrel or in the last part of the Endcap + uint layer = minMaxLayerId[minMaxLayerId.size()-1].first; + + if(m_radii.empty()) calculateLayerRadii(); + if(m_cellEdges.empty()) return 0; + + // get the last cell index (which is in the positive-z side) + int idx = abs(m_cellIndexes[layer].back()); + + // get the z-coordinate of the right-hand edge of the last cell + double zhigh = m_cellEdges[layer][idx].second; + + // get the inner radius of the first layer + double Rmin = m_radii[layer] - 0.5*m_layerDepth[layer]; + + // calculate the minimum theta of the last cell in the first layer -> this is the minimum theta of the detector (Barrel or Endcap) + double thetaMin = std::atan2(Rmin,zhigh); // theta min + + return (M_PI - thetaMin); // theta max +} + + } } From 652739174faf08306f2c89d9fd1c0f740c811ff4 Mon Sep 17 00:00:00 2001 From: Archil Durglishvili Date: Mon, 2 Dec 2024 13:04:19 +0100 Subject: [PATCH 15/17] Adding some description in the HCal xml files --- .../ALLEGRO_o1_v03/HCalBarrel_TileCal_v03.xml | 44 ++++++++++++++++-- .../HCalEndcaps_ThreeParts_TileCal_v03.xml | 46 +++++++++++++++++-- .../src/FCCSWHCalPhiTheta_k4geo.cpp | 2 +- 3 files changed, 81 insertions(+), 11 deletions(-) diff --git a/FCCee/ALLEGRO/compact/ALLEGRO_o1_v03/HCalBarrel_TileCal_v03.xml b/FCCee/ALLEGRO/compact/ALLEGRO_o1_v03/HCalBarrel_TileCal_v03.xml index 1adb2a8fd..8b9e9a07e 100644 --- a/FCCee/ALLEGRO/compact/ALLEGRO_o1_v03/HCalBarrel_TileCal_v03.xml +++ b/FCCee/ALLEGRO/compact/ALLEGRO_o1_v03/HCalBarrel_TileCal_v03.xml @@ -27,6 +27,10 @@ + + + + @@ -45,13 +49,28 @@ + system:4,layer:5,row:9,theta:9,phi:10 + diff --git a/FCCee/ALLEGRO/compact/ALLEGRO_o1_v03/HCalEndcaps_ThreeParts_TileCal_v03.xml b/FCCee/ALLEGRO/compact/ALLEGRO_o1_v03/HCalEndcaps_ThreeParts_TileCal_v03.xml index 28a31ba69..32f2978bf 100644 --- a/FCCee/ALLEGRO/compact/ALLEGRO_o1_v03/HCalEndcaps_ThreeParts_TileCal_v03.xml +++ b/FCCee/ALLEGRO/compact/ALLEGRO_o1_v03/HCalEndcaps_ThreeParts_TileCal_v03.xml @@ -25,6 +25,10 @@ + + + + @@ -44,13 +48,29 @@ + system:4,type:3,layer:6,row:11,theta:11,phi:10 + Date: Mon, 2 Dec 2024 15:41:50 +0100 Subject: [PATCH 16/17] Added some sanity checks of the HCal xml --- .../src/FCCSWHCalPhiRow_k4geo.cpp | 42 ++++++++++++++++++- .../src/FCCSWHCalPhiTheta_k4geo.cpp | 34 ++++++++++++++- 2 files changed, 74 insertions(+), 2 deletions(-) diff --git a/detectorSegmentations/src/FCCSWHCalPhiRow_k4geo.cpp b/detectorSegmentations/src/FCCSWHCalPhiRow_k4geo.cpp index 09d644982..3bae5de42 100644 --- a/detectorSegmentations/src/FCCSWHCalPhiRow_k4geo.cpp +++ b/detectorSegmentations/src/FCCSWHCalPhiRow_k4geo.cpp @@ -83,11 +83,51 @@ void FCCSWHCalPhiRow_k4geo::calculateLayerRadii() const { // check if all necessary variables are available if(m_detLayout==-1 || m_offsetZ.empty() || m_widthZ.empty() || m_offsetR.empty() || m_numLayers.empty() || m_dRlayer.empty()) { - dd4hep::printout(dd4hep::ERROR, "FCCSWHCalPhiRow_k4geo","Please check the readout description in the XML file!", + dd4hep::printout(dd4hep::ERROR, "FCCSWHCalPhiRow_k4geo","Please check the readout description in the XML file!\n%s", "One of the variables is missing: detLayout | offset_z | width_z | offset_r | numLayers | dRlayer"); return; } + // some sanity checks of the xml + if( m_offsetZ.size() != m_offsetR.size() ) + { + dd4hep::printout(dd4hep::ERROR, "FCCSWHCalPhiRow_k4geo","Please check the readout description in the XML file!\n%s", + "Number of elements in offsetZ and offsetR must be the same!"); + return; + } + if( m_widthZ.size() != m_offsetR.size() ) + { + dd4hep::printout(dd4hep::ERROR, "FCCSWHCalPhiRow_k4geo","Please check the readout description in the XML file!\n%s", + "Number of elements in widthZ and offsetR must be the same!"); + return; + } + if( m_detLayout == 0 && m_offsetZ.size() != 1) + { + dd4hep::printout(dd4hep::ERROR, "FCCSWHCalPhiRow_k4geo","Please check the readout description in the XML file!\n%s", + "Number of elements in offsetZ/offsetR/widthZ must be 1 for the Barrel!"); + return; + } + if( m_numLayers.size() % m_offsetZ.size() != 0 ) + { + dd4hep::printout(dd4hep::ERROR, "FCCSWHCalPhiRow_k4geo","Please check the readout description in the XML file!\n%s", + "Number of elements in numLayers must be multiple of offsetZ.size()!"); + return; + } + if( m_dRlayer.size() != m_numLayers.size()/m_offsetZ.size() ) + { + dd4hep::printout(dd4hep::ERROR, "FCCSWHCalPhiRow_k4geo","Please check the readout description in the XML file!\n%s", + "Number of elements in dRlayer must be equal to numLayers.size()/offsetZ.size()!"); + return; + } + uint nlayers = 0; + for(auto n : m_numLayers) nlayers+=n; + if( m_gridSizeRow.size() != nlayers ) + { + dd4hep::printout(dd4hep::ERROR, "FCCSWHCalPhiRow_k4geo","Please check the readout description in the XML file!\n%s", + "Number of elements in gridSizeRow must be equal to sum of contents of numLayers!"); + return; + } + if(m_detLayout==0) dd4hep::printout(dd4hep::INFO, "FCCSWHCalPhiRow_k4geo","Barrel configuration found!"); else dd4hep::printout(dd4hep::INFO, "FCCSWHCalPhiRow_k4geo","EndCap configuration found!"); diff --git a/detectorSegmentations/src/FCCSWHCalPhiTheta_k4geo.cpp b/detectorSegmentations/src/FCCSWHCalPhiTheta_k4geo.cpp index 07770e633..08346e4f6 100644 --- a/detectorSegmentations/src/FCCSWHCalPhiTheta_k4geo.cpp +++ b/detectorSegmentations/src/FCCSWHCalPhiTheta_k4geo.cpp @@ -79,11 +79,43 @@ void FCCSWHCalPhiTheta_k4geo::defineCellsInRZplan() const { // check if all necessary variables are available if(m_detLayout==-1 || m_offsetZ.empty() || m_widthZ.empty() || m_offsetR.empty() || m_numLayers.empty() || m_dRlayer.empty()) { - dd4hep::printout(dd4hep::ERROR, "FCCSWHCalPhiRow_k4geo","Please check the readout description in the XML file!", + dd4hep::printout(dd4hep::ERROR, "FCCSWHCalPhiRow_k4geo","Please check the readout description in the XML file!\n%s", "One of the variables is missing: detLayout | offset_z | width_z | offset_r | numLayers | dRlayer"); return; } + // some sanity checks of the xml + if( m_offsetZ.size() != m_offsetR.size() ) + { + dd4hep::printout(dd4hep::ERROR, "FCCSWHCalPhiRow_k4geo","Please check the readout description in the XML file!\n%s", + "Number of elements in offsetZ and offsetR must be the same!"); + return; + } + if( m_widthZ.size() != m_offsetR.size() ) + { + dd4hep::printout(dd4hep::ERROR, "FCCSWHCalPhiRow_k4geo","Please check the readout description in the XML file!\n%s", + "Number of elements in widthZ and offsetR must be the same!"); + return; + } + if( m_detLayout == 0 && m_offsetZ.size() != 1) + { + dd4hep::printout(dd4hep::ERROR, "FCCSWHCalPhiRow_k4geo","Please check the readout description in the XML file!\n%s", + "Number of elements in offsetZ/offsetR/widthZ must be 1 for the Barrel!"); + return; + } + if( m_numLayers.size() % m_offsetZ.size() != 0 ) + { + dd4hep::printout(dd4hep::ERROR, "FCCSWHCalPhiRow_k4geo","Please check the readout description in the XML file!\n%s", + "Number of elements in numLayers must be multiple of offsetZ.size()!"); + return; + } + if( m_dRlayer.size() != m_numLayers.size()/m_offsetZ.size() ) + { + dd4hep::printout(dd4hep::ERROR, "FCCSWHCalPhiRow_k4geo","Please check the readout description in the XML file!\n%s", + "Number of elements in dRlayer must be equal to numLayers.size()/offsetZ.size()!"); + return; + } + if(m_detLayout==0) dd4hep::printout(dd4hep::INFO, "FCCSWHCalPhiRow_k4geo","Barrel configuration found!"); else dd4hep::printout(dd4hep::INFO, "FCCSWHCalPhiRow_k4geo","EndCap configuration found!"); From 40312528252fca0f90812e273bae16e0ced1adb9 Mon Sep 17 00:00:00 2001 From: Archil Durglishvili Date: Mon, 2 Dec 2024 15:57:06 +0100 Subject: [PATCH 17/17] Fixing typo: PhiRow to PhiTheta --- .../src/FCCSWHCalPhiTheta_k4geo.cpp | 16 ++++++++-------- 1 file changed, 8 insertions(+), 8 deletions(-) diff --git a/detectorSegmentations/src/FCCSWHCalPhiTheta_k4geo.cpp b/detectorSegmentations/src/FCCSWHCalPhiTheta_k4geo.cpp index 08346e4f6..e71688951 100644 --- a/detectorSegmentations/src/FCCSWHCalPhiTheta_k4geo.cpp +++ b/detectorSegmentations/src/FCCSWHCalPhiTheta_k4geo.cpp @@ -79,7 +79,7 @@ void FCCSWHCalPhiTheta_k4geo::defineCellsInRZplan() const { // check if all necessary variables are available if(m_detLayout==-1 || m_offsetZ.empty() || m_widthZ.empty() || m_offsetR.empty() || m_numLayers.empty() || m_dRlayer.empty()) { - dd4hep::printout(dd4hep::ERROR, "FCCSWHCalPhiRow_k4geo","Please check the readout description in the XML file!\n%s", + dd4hep::printout(dd4hep::ERROR, "FCCSWHCalPhiTheta_k4geo","Please check the readout description in the XML file!\n%s", "One of the variables is missing: detLayout | offset_z | width_z | offset_r | numLayers | dRlayer"); return; } @@ -87,37 +87,37 @@ void FCCSWHCalPhiTheta_k4geo::defineCellsInRZplan() const { // some sanity checks of the xml if( m_offsetZ.size() != m_offsetR.size() ) { - dd4hep::printout(dd4hep::ERROR, "FCCSWHCalPhiRow_k4geo","Please check the readout description in the XML file!\n%s", + dd4hep::printout(dd4hep::ERROR, "FCCSWHCalPhiTheta_k4geo","Please check the readout description in the XML file!\n%s", "Number of elements in offsetZ and offsetR must be the same!"); return; } if( m_widthZ.size() != m_offsetR.size() ) { - dd4hep::printout(dd4hep::ERROR, "FCCSWHCalPhiRow_k4geo","Please check the readout description in the XML file!\n%s", + dd4hep::printout(dd4hep::ERROR, "FCCSWHCalPhiTheta_k4geo","Please check the readout description in the XML file!\n%s", "Number of elements in widthZ and offsetR must be the same!"); return; } if( m_detLayout == 0 && m_offsetZ.size() != 1) { - dd4hep::printout(dd4hep::ERROR, "FCCSWHCalPhiRow_k4geo","Please check the readout description in the XML file!\n%s", + dd4hep::printout(dd4hep::ERROR, "FCCSWHCalPhiTheta_k4geo","Please check the readout description in the XML file!\n%s", "Number of elements in offsetZ/offsetR/widthZ must be 1 for the Barrel!"); return; } if( m_numLayers.size() % m_offsetZ.size() != 0 ) { - dd4hep::printout(dd4hep::ERROR, "FCCSWHCalPhiRow_k4geo","Please check the readout description in the XML file!\n%s", + dd4hep::printout(dd4hep::ERROR, "FCCSWHCalPhiTheta_k4geo","Please check the readout description in the XML file!\n%s", "Number of elements in numLayers must be multiple of offsetZ.size()!"); return; } if( m_dRlayer.size() != m_numLayers.size()/m_offsetZ.size() ) { - dd4hep::printout(dd4hep::ERROR, "FCCSWHCalPhiRow_k4geo","Please check the readout description in the XML file!\n%s", + dd4hep::printout(dd4hep::ERROR, "FCCSWHCalPhiTheta_k4geo","Please check the readout description in the XML file!\n%s", "Number of elements in dRlayer must be equal to numLayers.size()/offsetZ.size()!"); return; } - if(m_detLayout==0) dd4hep::printout(dd4hep::INFO, "FCCSWHCalPhiRow_k4geo","Barrel configuration found!"); - else dd4hep::printout(dd4hep::INFO, "FCCSWHCalPhiRow_k4geo","EndCap configuration found!"); + if(m_detLayout==0) dd4hep::printout(dd4hep::INFO, "FCCSWHCalPhiTheta_k4geo","Barrel configuration found!"); + else dd4hep::printout(dd4hep::INFO, "FCCSWHCalPhiTheta_k4geo","EndCap configuration found!"); // calculate the radius for each layer uint N_dR = m_numLayers.size()/m_offsetZ.size();