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..8b9e9a07e --- /dev/null +++ b/FCCee/ALLEGRO/compact/ALLEGRO_o1_v03/HCalBarrel_TileCal_v03.xml @@ -0,0 +1,183 @@ + + + + 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_v03/HCalEndcaps_ThreeParts_TileCal_v03.xml b/FCCee/ALLEGRO/compact/ALLEGRO_o1_v03/HCalEndcaps_ThreeParts_TileCal_v03.xml new file mode 100644 index 000000000..32f2978bf --- /dev/null +++ b/FCCee/ALLEGRO/compact/ALLEGRO_o1_v03/HCalEndcaps_ThreeParts_TileCal_v03.xml @@ -0,0 +1,197 @@ + + + + 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/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 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 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 diff --git a/detectorSegmentations/include/detectorSegmentations/FCCSWHCalPhiRowHandle_k4geo.h b/detectorSegmentations/include/detectorSegmentations/FCCSWHCalPhiRowHandle_k4geo.h new file mode 100644 index 000000000..5ca4fe05a --- /dev/null +++ b/detectorSegmentations/include/detectorSegmentations/FCCSWHCalPhiRowHandle_k4geo.h @@ -0,0 +1,119 @@ +#ifndef DETECTORSEGMENTATIONS_HCALPHIROWHANDLE_K4GEO_H +#define DETECTORSEGMENTATIONS_HCALPHIROWHANDLE_K4GEO_H + +// FCCSW +#include "detectorSegmentations/FCCSWHCalPhiRow_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> FCCSWHCalPhiRowHandle_k4geo; + +/// Implementation class for the HCal phi-row 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 FCCSWHCalPhiRow_k4geo : public FCCSWHCalPhiRowHandle_k4geo { +public: + /// Defintion of the basic handled object + typedef FCCSWHCalPhiRowHandle_k4geo::Object Object; + +public: + /// Default constructor + FCCSWHCalPhiRow_k4geo() = default; + /// Copy constructor + FCCSWHCalPhiRow_k4geo(const FCCSWHCalPhiRow_k4geo& e) = default; + /// Copy Constructor from segmentation base object + FCCSWHCalPhiRow_k4geo(const Segmentation& e) : Handle(e) {} + /// Copy constructor from handle + FCCSWHCalPhiRow_k4geo(const Handle& e) : Handle(e) {} + /// Copy constructor from other polymorph/equivalent handle + template + FCCSWHCalPhiRow_k4geo(const Handle& e) : Handle(e) {} + /// Assignment operator + FCCSWHCalPhiRow_k4geo& operator=(const FCCSWHCalPhiRow_k4geo& seg) = default; + /// Equality operator + bool operator==(const FCCSWHCalPhiRow_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 row for each layer + inline std::vector gridSizeRow() const { return access()->implementation->gridSizeRow(); } + + /// access the grid size in Phi + inline int phiBins() const { return access()->implementation->phiBins(); } + + /// access the coordinate offset in Phi + inline double offsetPhi() const { return access()->implementation->offsetPhi(); } + + /// 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); } + + /// 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 setGridSizeRow(std::vector const&cellSize) const { access()->implementation->setGridSizeRow(cellSize); } + + /// 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(); } + + /// access the field name used for Phi + inline const std::string& fieldNamePhi() const { return access()->implementation->fieldNamePhi(); } + +}; + +} /* End namespace dd4hep */ +#endif // DETECTORSEGMENTATIONS_HCALPHITHETAHANDLE_K4GEO_H diff --git a/detectorSegmentations/include/detectorSegmentations/FCCSWHCalPhiRow_k4geo.h b/detectorSegmentations/include/detectorSegmentations/FCCSWHCalPhiRow_k4geo.h new file mode 100644 index 000000000..1df1353d9 --- /dev/null +++ b/detectorSegmentations/include/detectorSegmentations/FCCSWHCalPhiRow_k4geo.h @@ -0,0 +1,276 @@ +#ifndef DETECTORSEGMENTATIONS_HCALPHIROW_K4GEO_H +#define DETECTORSEGMENTATIONS_HCALPHIROW_K4GEO_H + +#include "DDSegmentation/SegmentationUtil.h" +#include "DDSegmentation/Segmentation.h" +#include "TVector3.h" + + +/** FCCSWHCalPhiRow_k4geo Detector/detectorSegmentations/detectorSegmentations/FCCSWHCalPhiRow_k4geo.h FCCSWHCalPhiRow_k4geo.h + * + * 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 + * + */ + +namespace dd4hep { + namespace DDSegmentation { + class FCCSWHCalPhiRow_k4geo : public Segmentation { + public: + /// default constructor using an arbitrary type + FCCSWHCalPhiRow_k4geo(const std::string& aCellEncoding); + /// Default constructor used by derived classes passing an existing decoder + FCCSWHCalPhiRow_k4geo(const BitFieldCoder* decoder); + + /// destructor + virtual ~FCCSWHCalPhiRow_k4geo() = default; + + /** Determine the position of HCal cell based on the cellID. + * @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; + + /** 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. + */ + std::vector neighbours(const CellID& cID) const; + + + /** Calculate layer radii and edges in z-axis. + * Following member variables are calculated: + * m_radii + * m_layerEdges + * m_layerDepth + */ + void calculateLayerRadii() const; + + /** Define cell indexes for the given layer. + * This function fills the m_cellIndexes vector per layer with the cell indexes. + * The cell index is encoded in CellID with "row" field. + * 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; + + /** 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, which is the center of cell with m_phiID=0. + * return The offset in phi. + */ + inline double offsetPhi() const { return m_offsetPhi; } + + /** Get the grid size in row for each layer. + * return Grid size in row. + */ + inline std::vector gridSizeRow() const { return m_gridSizeRow; } + + /** Determine the polar angle of HCal cell center based on the cellID. + * @param[in] aCellId ID of a cell. + * return Theta. + */ + inline double theta(const CellID& aCellID) const { return dd4hep::DDSegmentation::Util::thetaFromXYZ(position(aCellID)); } + + /** 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 vector of cell indexes in a given layer. + */ + inline std::vector cellIndexes(const uint layer) const { + if(m_radii.empty()) calculateLayerRadii(); + if(!m_cellIndexes.empty()) return m_cellIndexes[layer]; + 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. + */ + std::vector > getMinMaxLayerId() const ; + + /** 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 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 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 (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; } + + /** Get the field name for azimuthal angle. + * return The field name for phi. + */ + 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. + */ + inline const std::string& fieldNameRow() const { return m_rowID; } + + /** 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 grid size in theta angle. + * @param[in] aCellSize Cell size in theta. + */ + 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). + */ + 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; } + + /** Set the field name used for row number. + * @param[in] aFieldName Field name for row. + */ + inline void setFieldNameRow(const std::string& fieldName) { m_rowID = fieldName; } + + + protected: + /// 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 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 + 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; + /// dR of each layer + mutable std::vector m_layerDepth; + /// cell indexes in each layer + mutable std::vector > m_cellIndexes; + /// z-min and z-max of each cell in each layer + mutable std::vector > > m_cellEdges; + }; + } +} +#endif /* DETSEGMENTATION_HCALPHITHETA_H */ diff --git a/detectorSegmentations/include/detectorSegmentations/FCCSWHCalPhiThetaHandle_k4geo.h b/detectorSegmentations/include/detectorSegmentations/FCCSWHCalPhiThetaHandle_k4geo.h new file mode 100644 index 000000000..a9f128147 --- /dev/null +++ b/detectorSegmentations/include/detectorSegmentations/FCCSWHCalPhiThetaHandle_k4geo.h @@ -0,0 +1,138 @@ +#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> FCCSWHCalPhiThetaHandle_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 FCCSWHCalPhiThetaHandle_k4geo { +public: + /// Defintion of the basic handled object + typedef FCCSWHCalPhiThetaHandle_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 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); } + + /// 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(); } + + /// 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) + + 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..e02ad0ede --- /dev/null +++ b/detectorSegmentations/include/detectorSegmentations/FCCSWHCalPhiTheta_k4geo.h @@ -0,0 +1,235 @@ +#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. + * 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 + * + */ + +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; + + /** 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; + + /** 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; + + /** 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; + + /** 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, 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 given 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. + * 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 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 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 (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; } + + /** Get the field name for azimuthal angle. + * return The field name for phi. + */ + 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. + */ + 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. + */ + 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 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). + */ + 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 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 + 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; + /// 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 + mutable std::vector > > m_cellEdges; + }; + } +} +#endif /* DETSEGMENTATION_HCALPHITHETA_H */ diff --git a/detectorSegmentations/src/FCCSWHCalPhiRowHandle_k4geo.cpp b/detectorSegmentations/src/FCCSWHCalPhiRowHandle_k4geo.cpp new file mode 100644 index 000000000..4293df1e5 --- /dev/null +++ b/detectorSegmentations/src/FCCSWHCalPhiRowHandle_k4geo.cpp @@ -0,0 +1,4 @@ +#include "detectorSegmentations/FCCSWHCalPhiRowHandle_k4geo.h" +#include "DD4hep/detail/Handle.inl" + +DD4HEP_INSTANTIATE_HANDLE_UNNAMED(dd4hep::DDSegmentation::FCCSWHCalPhiRow_k4geo); diff --git a/detectorSegmentations/src/FCCSWHCalPhiRow_k4geo.cpp b/detectorSegmentations/src/FCCSWHCalPhiRow_k4geo.cpp new file mode 100644 index 000000000..3bae5de42 --- /dev/null +++ b/detectorSegmentations/src/FCCSWHCalPhiRow_k4geo.cpp @@ -0,0 +1,673 @@ +#include "detectorSegmentations/FCCSWHCalPhiRow_k4geo.h" +#include "DD4hep/Printout.h" + +namespace dd4hep { +namespace DDSegmentation { + +using std::runtime_error; + + +/// default constructor using an encoding string +FCCSWHCalPhiRow_k4geo::FCCSWHCalPhiRow_k4geo(const std::string& cellEncoding) : Segmentation(cellEncoding) { + // define type and description + _type = "FCCSWHCalPhiRow_k4geo"; + _description = "Phi-theta segmentation in the global coordinates"; + + // 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", "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()); + 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_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) { + // define type and description + _type = "FCCSWHCalPhiRow_k4geo"; + _description = "Phi-theta segmentation in the global coordinates"; + + // 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", "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()); + 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_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,m_layerID); + + if(m_radii.empty()) calculateLayerRadii(); + if(m_radii.empty() || m_layerEdges.empty()) + { + dd4hep::printout(dd4hep::ERROR, "FCCSWHCalPhiRow_k4geo","Could not calculate layer radii!"); + return Vector3D(0.,0.,0.); + } + + double radius = m_radii[layer]; + double minLayerZ = m_layerEdges[layer].first; + + // get index of the cell in the layer (index starts from 1!) + 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]; + + // 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(cID)), radius * std::sin(phi(cID)), zpos ); +} + + +void FCCSWHCalPhiRow_k4geo::calculateLayerRadii() const { + if(m_radii.empty()) + { + // 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", + "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!"); + + // 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++) + { + for(uint i_dR = 0; i_dR < N_dR; i_dR++) + { + for(int i_row = 1; i_row <= m_numLayers[i_dR+i_section*N_dR]; i_row++) + { + 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[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]); + } + } + } + + // print info of calculated radii and edges + for(uint i_layer = 0; i_layer < m_radii.size(); i_layer++){ + dd4hep::printout(dd4hep::INFO, "FCCSWHCalPhiRow_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 cellIndexes vector for each layer + m_cellIndexes.resize(m_radii.size()); + // allocate cellEdges vector for each layer + m_cellEdges.resize(m_radii.size()); + + // determine row bins for each layer + for(uint i_layer = 0; i_layer < m_radii.size(); i_layer++) defineCellIndexes(i_layer); + } +} + +/* +* This function fills the m_cellIndexes vector per layer with the cell indexes. +* The cell index is encoded in CellID with "row" field. +* 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 (-1, -2, ..., -N). +* +* m_cellIndexes vector can be used to define the neighbours using CreateFCCeeCaloNeighbours tool. +*/ +void FCCSWHCalPhiRow_k4geo::defineCellIndexes(const uint layer) const { + if(m_cellIndexes[layer].size()==0 && m_radii.size() > 0) + { + double minLayerZ = m_layerEdges[layer].first; + double maxLayerZ = m_layerEdges[layer].second; + + // find rows/sequences that fit within the given layer range along z + int irow = 0; + while( (minLayerZ + (irow+1)* m_dz_row) < (maxLayerZ+0.0001) ) + { + // define the cell index + int idx = floor(irow/m_gridSizeRow[layer]) + 1; + // add the index if it is not already there + if(std::find(m_cellIndexes[layer].begin(),m_cellIndexes[layer].end(),idx) == m_cellIndexes[layer].end()) m_cellIndexes[layer].push_back(idx); + irow++; + } + + // for the EndCap, do it again but for negative-z part + if(m_detLayout == 1) + { + irow = 0; + while( (minLayerZ + (irow+1) * m_dz_row) < (maxLayerZ+0.0001) ) + { + // define the cell index with negative sign + int idx = - ( floor(irow/m_gridSizeRow[layer]) + 1 ); + // add the index if it is not already there + if(std::find(m_cellIndexes[layer].begin(),m_cellIndexes[layer].end(),idx) == m_cellIndexes[layer].end()) m_cellIndexes[layer].push_back(idx); + irow++; + } + } + + // find edges of each cell in the given layer along z axis + for(auto idx : m_cellIndexes[layer]) + { + // calculate z-coordinates of the cell edges + 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 negative (starts from -1!) + if(idx < 0) + { + z1 = -minLayerZ + (idx) * m_dz_row * m_gridSizeRow[layer]; // lower edge + z2 = -minLayerZ + (idx+1) * m_dz_row * m_gridSizeRow[layer]; // upper edge + } + + m_cellEdges[layer][idx] = std::make_pair(z1,z2); + } + + dd4hep::printout(dd4hep::DEBUG, "FCCSWHCalPhiRow_k4geo","Number of cells in layer %d: %d", layer, m_cellIndexes[layer].size()); + } +} + +/// create the cell ID based on the position +CellID FCCSWHCalPhiRow_k4geo::cellID(const Vector3D& /* localPosition */, const Vector3D& globalPosition, + const VolumeID& vID) const { + + // get the row number from volumeID (starts from 0!) + int nrow = _decoder->get(vID, m_rowID); + // get the layer number from volumeID + uint layer = _decoder->get(vID,m_layerID); + + CellID cID = vID; + + // get the cell index (start from 1!) + 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_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_detLayout == 1) _decoder->set(cID, "type", 0); + + return cID; +} + +/// determine the azimuthal angle phi based on the cell ID +double FCCSWHCalPhiRow_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 > FCCSWHCalPhiRow_k4geo::getMinMaxLayerId() const { + + std::vector > minMaxLayerId; + + 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); + + uint Nl = m_numLayers.size()/m_offsetZ.size(); + + 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]; + } + + for(uint i=0; i < Nl; i++) + { + maxLayerId[i_section] += m_numLayers[i+i_section*Nl]; + } + + if(i_section==0) maxLayerId[0] -= 1; + + minMaxLayerId.push_back(std::make_pair(minLayerId[i_section], maxLayerId[i_section])); + } + + return minMaxLayerId; +} + + +/// Calculates the neighbours of the given cell ID and adds them to the list of neighbours +std::vector FCCSWHCalPhiRow_k4geo::neighbours(const CellID& cID) const { + std::vector cellNeighbours; + + if(m_radii.empty()) calculateLayerRadii(); + if(m_cellIndexes.empty()) return cellNeighbours; + + uint EndcapPart = 0; + int minLayerId = -1; + int maxLayerId = -1; + + int currentLayerId = _decoder->get(cID,m_layerID); + int currentCellId = _decoder->get(cID,m_rowID); + + int minCellId = m_cellIndexes[currentLayerId].front(); + int maxCellId = m_cellIndexes[currentLayerId].back(); + + //-------------------------------- + // Determine min and max layer Id + //-------------------------------- + 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 Endcap min and max layer indexes! --> returning empty neighbours"); + return cellNeighbours; + } + + // determine min and max layer Ids and which section/part it is + for(uint i_section = 0; i_section < minMaxLayerId.size(); i_section++) + { + 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 + 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,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]) + { + _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,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]) + { + _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 (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]; + + // 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 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++) + { + 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_rowID,minCellId); + cellNeighbours.push_back(nID); // add the first cell from part2 layer + } + } + } + + // if the Endcap consists of more than 2 parts: + for(uint i_section = 1; i_section < (m_offsetZ.size()-1); i_section++) + { + if(i_section != EndcapPart) continue; + + // 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++) + { + 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 then look for neighbours in the next part + if(currentCellId == maxCellId) + { + // 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++) + { + 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 last part of the Endcap + if(EndcapPart == (m_offsetZ.size() - 1) && currentCellId == minCellId) + { + // 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[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 part2 layer + } + } + } + } + + // Now loop over the neighbours and add the cells from next/previous phi module + 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 + _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,m_layerID); + + 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; +} + +/// 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 +} + + +} +} 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..e71688951 --- /dev/null +++ b/detectorSegmentations/src/FCCSWHCalPhiTheta_k4geo.cpp @@ -0,0 +1,811 @@ +#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("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) { + // 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("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,m_layerID); + 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::position(const CellID& cID) const { + uint layer = _decoder->get(cID,m_layerID); + int thetaID = _decoder->get(cID,m_thetaID); + 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)); + + // return the position with corrected z coordinate 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_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!\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, "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, "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, "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, "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, "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, "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(); + std::vector moduleDepth(m_offsetZ.size()); + for(uint i_section = 0; i_section < m_offsetZ.size(); i_section++) + { + for(uint i_dR = 0; i_dR < N_dR; i_dR++) + { + for(int i_row = 1; i_row <= m_numLayers[i_dR+i_section*N_dR]; i_row++) + { + 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[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]); + } + } + } + + // 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_detLayout == 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_detLayout == 1) _decoder->set(cID, "type", 0); + + double lTheta = thetaFromXYZ(globalPosition); + double lPhi = phiFromXYZ(globalPosition); + uint layer = _decoder->get(vID,m_layerID); + + // 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); +} + + +/// Get the min and max layer indexes for each part of the HCal +std::vector > FCCSWHCalPhiTheta_k4geo::getMinMaxLayerId() const { + 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); + + uint Nl = m_numLayers.size()/m_offsetZ.size(); + + 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]; + } + + for(uint i=0; i < Nl; i++) + { + maxLayerId[i_section] += m_numLayers[i+i_section*Nl]; + } + + if(i_section == 0) maxLayerId[0] -= 1; + + minMaxLayerId.push_back(std::make_pair(minLayerId[i_section], maxLayerId[i_section])); + } + + 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, bool aDiagonal) const { + std::vector cellNeighbours; + + if(m_radii.empty()) defineCellsInRZplan(); + if(m_thetaBins.empty()) return cellNeighbours; + + uint EndcapPart = 0; + int minLayerId = -1; + int maxLayerId = -1; + + int currentLayerId = _decoder->get(cID,m_layerID); + 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 + //-------------------------------- + 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 Endcap min and max layer indexes! --> returning empty neighbours"); + return cellNeighbours; + } + + // determine min and max layer Ids and which section/part it is + for(uint i_section = 0; i_section < minMaxLayerId.size(); i_section++) + { + 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 + 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_detLayout == 0) + { + 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,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 + + // 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(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) + { + CellID nID = cID ; + int nextLayerId = currentLayerId + 1; + _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 + + // 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); + + if(aDiagonal) + { + //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); + + if(aDiagonal) + { + //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); + } + } + } + } + } + + // Endcap + if(m_detLayout == 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,m_layerID,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 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 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 + } + } + } + } + // 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,m_layerID,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 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.) + { + 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 + } + } + } + } + + // if the Endcap consists of more than 1 part/section 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 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(EndcapPart == 0 && 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,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 the Endcap consists of more than 2 parts: + for(uint i_section = 1; i_section < (m_offsetZ.size()-1); i_section++) + { + if(i_section != EndcapPart) continue; + + // if it is the last theta bin cell then look for neighbours in previous part + if(currentCellThetaBin == maxCellThetaBin) + { + // 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 it is the first theta bin cell then look for neighbours in the next part + if(currentCellThetaBin == minCellThetaBin) + { + // 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 last part of the Endcap + if(EndcapPart == (m_offsetZ.size() - 1) && currentCellThetaBin == maxCellThetaBin) + { + // 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[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 + } + } + } + } + } + + // Now loop over the neighbours and add the cells from next/previous phi module + 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 + _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,m_layerID); + + 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; +} + +} +} diff --git a/detectorSegmentations/src/plugins/SegmentationFactories.cpp b/detectorSegmentations/src/plugins/SegmentationFactories.cpp index d362660de..a3ac9b0e9 100644 --- a/detectorSegmentations/src/plugins/SegmentationFactories.cpp +++ b/detectorSegmentations/src/plugins/SegmentationFactories.cpp @@ -34,3 +34,9 @@ DECLARE_SEGMENTATION(GridDRcalo_k4geo, create_segmentation) + +#include "detectorSegmentations/FCCSWHCalPhiTheta_k4geo.h" +DECLARE_SEGMENTATION(FCCSWHCalPhiTheta_k4geo, create_segmentation) + +#include "detectorSegmentations/FCCSWHCalPhiRow_k4geo.h" +DECLARE_SEGMENTATION(FCCSWHCalPhiRow_k4geo, create_segmentation) 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" )