diff --git a/FCCee/ALLEGRO/compact/ALLEGRO_o1_v03/HCalBarrel_TileCal_v03.xml b/FCCee/ALLEGRO/compact/ALLEGRO_o1_v03/HCalBarrel_TileCal_v03.xml index f5ce03fe6..1adb2a8fd 100644 --- a/FCCee/ALLEGRO/compact/ALLEGRO_o1_v03/HCalBarrel_TileCal_v03.xml +++ b/FCCee/ALLEGRO/compact/ALLEGRO_o1_v03/HCalBarrel_TileCal_v03.xml @@ -46,6 +46,7 @@ implementation->setOffsetPhi(offset); } + /// set the detector layout + inline void setDetLayout(int detLayout) const { access()->implementation->setDetLayout(detLayout); } + /// set the coordinate offset in z-axis inline void setOffsetZ(std::vector const&offset) const { access()->implementation->setOffsetZ(offset); } @@ -101,6 +104,9 @@ class FCCSWHCalPhiRow_k4geo : public FCCSWHCalPhiRowHandle_k4geo { /// set the grid size in Phi inline void setPhiBins(int cellSize) const { access()->implementation->setPhiBins(cellSize); } + /// access the field name used for layer + inline const std::string& fieldNameLayer() const { return access()->implementation->fieldNameLayer(); } + /// access the field name used for theta inline const std::string& fieldNameRow() const { return access()->implementation->fieldNameRow(); } diff --git a/detectorSegmentations/include/detectorSegmentations/FCCSWHCalPhiRow_k4geo.h b/detectorSegmentations/include/detectorSegmentations/FCCSWHCalPhiRow_k4geo.h index ce9b6f44f..acda782e8 100644 --- a/detectorSegmentations/include/detectorSegmentations/FCCSWHCalPhiRow_k4geo.h +++ b/detectorSegmentations/include/detectorSegmentations/FCCSWHCalPhiRow_k4geo.h @@ -11,6 +11,9 @@ * Segmentation in row and phi. * * Cells are defined in r-z plan by merging the row/sequence of scintillator/absorber. + * Each row consists of 2 * Master_Plate + 1 * Spacer_Plate + 1 * Scintillator. + * Considering a single row as a single cell gives the highest possible granularity. + * More details: https://indico.cern.ch/event/1475808/contributions/6219554/attachments/2966253/5218774/FCC_FullSim_HCal_slides.pdf * */ @@ -42,6 +45,7 @@ namespace dd4hep { const VolumeID& aVolumeID) const; /** Find neighbours of the cell + * Definition of neighbours is explained on slide 7: https://indico.cern.ch/event/1475808/contributions/6219554/attachments/2966253/5218774/FCC_FullSim_HCal_slides.pdf * @param[in] aCellId ID of a cell. * return vector of neighbour cellIDs. */ @@ -49,6 +53,10 @@ namespace dd4hep { /** Calculate layer radii and edges in z-axis. + * Following member variables are calculated: + * m_radii + * m_layerEdges + * m_layerDepth */ void calculateLayerRadii() const; @@ -58,7 +66,9 @@ namespace dd4hep { * In case of a cell with single row/sequence, the index is directly the number of row in the layer. * In case of a cell with several rows/sequences merged, the index is the number of cell in the layer. * For the layers of negative-z Endcap, indexes of cells are negative. - * + * Following member variables are calculated: + * m_cellIndexes + * m_cellEdges * @param[in] layer index */ void defineCellIndexes(const unsigned int layer) const; @@ -79,7 +89,7 @@ namespace dd4hep { */ inline int phiBins() const { return m_phiBins; } - /** Get the coordinate offset in azimuthal angle. + /** Get the coordinate offset in azimuthal angle, which is the center of cell with m_phiID=0. * return The offset in phi. */ inline double offsetPhi() const { return m_offsetPhi; } @@ -114,27 +124,34 @@ namespace dd4hep { */ std::vector > getMinMaxLayerId() const ; - /** Get the coordinate offset in z-axis. + /** Get the coordinate offset in z-axis. + * Offset is the middle position of the Barrel or each section of the Endcap. + * For the Barrel, the vector size is 1, while for the Endcap - number of section. * return The offset in z. */ inline std::vector offsetZ() const { return m_offsetZ; } - /** Get the z width of the layer. - * return the z width. + /** Get the z length of the layer. + * return the z length. */ inline std::vector widthZ() const { return m_widthZ; } /** Get the coordinate offset in radius. + * Offset is the inner radius of the first layer in the Barrel or in each section of the Endcap. + * For the Barrel, the vector size is 1, while for the Endcap - number of sections. * return the offset in radius. */ inline std::vector offsetR() const { return m_offsetR; } - /** Get the number of layers. + /** Get the number of layers for each different thickness retrieved with dRlayer(). + * For the Barrel, the vector size equals to the number of different thicknesses used to form the layers. + * For the Endcap, the vector size equals to the number of sections in the Endcap times the number of different thicknesses used to form the layers. * return the number of layers. */ inline std::vector numLayers() const { return m_numLayers; } - /** Get the dR of layers. + /** Get the dR (thickness) of layers. + * The size of the vector equals to the number of different thicknesses used to form the layers. * return the dR. */ inline std::vector dRlayer() const { return m_dRlayer; } @@ -144,6 +161,11 @@ namespace dd4hep { */ inline const std::string& fieldNamePhi() const { return m_phiID; } + /** Get the field name for layer. + * return The field name for layer. + */ + inline const std::string& fieldNameLayer() const { return m_layerID; } + /** Get the field name for row number. * return The field name for row. */ @@ -164,6 +186,11 @@ namespace dd4hep { */ inline void setGridSizeRow(std::vector const&size) { m_gridSizeRow = size; } + /** Set the detector layout (0 = Barrel; 1 = Endcap). + * @param[in] detLayout + */ + inline void setDetLayout(int detLayout) { m_detLayout = detLayout; } + /** Set the coordinate offset in z-axis. * @param[in] offset in z (centre of the layer). */ @@ -208,12 +235,16 @@ namespace dd4hep { double m_offsetPhi; /// the field name used for phi std::string m_phiID; + /// the field name used for layer + std::string m_layerID; /// the grid size in row for each layer std::vector m_gridSizeRow; /// dz of row double m_dz_row; /// the field name used for row std::string m_rowID; + /// the detector layout (0 = Barrel; 1 = Endcap) + int m_detLayout; /// the z offset of middle of the layer std::vector m_offsetZ; /// the z width of the layer diff --git a/detectorSegmentations/include/detectorSegmentations/FCCSWHCalPhiThetaHandle_k4geo.h b/detectorSegmentations/include/detectorSegmentations/FCCSWHCalPhiThetaHandle_k4geo.h index 475b39649..a9f128147 100644 --- a/detectorSegmentations/include/detectorSegmentations/FCCSWHCalPhiThetaHandle_k4geo.h +++ b/detectorSegmentations/include/detectorSegmentations/FCCSWHCalPhiThetaHandle_k4geo.h @@ -86,6 +86,9 @@ class FCCSWHCalPhiTheta_k4geo : public FCCSWHCalPhiThetaHandle_k4geo { /// set the coordinate offset in Phi inline void setOffsetPhi(double offset) const { access()->implementation->setOffsetPhi(offset); } + /// set the detector layout + inline void setDetLayout(int detLayout) const { access()->implementation->setDetLayout(detLayout); } + /// set the coordinate offset in z-axis inline void setOffsetZ(std::vector const&offset) const { access()->implementation->setOffsetZ(offset); } @@ -113,6 +116,10 @@ class FCCSWHCalPhiTheta_k4geo : public FCCSWHCalPhiThetaHandle_k4geo { /// access the field name used for Phi inline const std::string& fieldNamePhi() const { return access()->implementation->fieldNamePhi(); } + /// access the field name used for layer + inline const std::string& fieldNameLayer() const { return access()->implementation->fieldNameLayer(); } + + /** \brief Returns a std::vector of the cellDimensions of the given cell ID in natural order of dimensions (dPhi, dTheta) diff --git a/detectorSegmentations/include/detectorSegmentations/FCCSWHCalPhiTheta_k4geo.h b/detectorSegmentations/include/detectorSegmentations/FCCSWHCalPhiTheta_k4geo.h index 607caa2c4..e02ad0ede 100644 --- a/detectorSegmentations/include/detectorSegmentations/FCCSWHCalPhiTheta_k4geo.h +++ b/detectorSegmentations/include/detectorSegmentations/FCCSWHCalPhiTheta_k4geo.h @@ -10,6 +10,8 @@ * Based on GridTheta_k4geo, addition of azimuthal angle coordinate. * * Rectangular shape cells are defined in r-z plan and all the hits within a defined cell boundaries are assigned the same cellID. + * Cell borders are defined such that closely follow the theta projective towers. + * More details: https://indico.cern.ch/event/1475808/contributions/6219554/attachments/2966253/5218774/FCC_FullSim_HCal_slides.pdf * */ @@ -40,13 +42,21 @@ namespace dd4hep { virtual CellID cellID(const Vector3D& aLocalPosition, const Vector3D& aGlobalPosition, const VolumeID& aVolumeID) const; - /** Find neighbours of the cell + /** Find neighbours of the cell. + * Definition of neighbours is explained on slide 9: https://indico.cern.ch/event/1475808/contributions/6219554/attachments/2966253/5218774/FCC_FullSim_HCal_slides.pdf * @param[in] aCellId ID of a cell. + * @param[in] aDiagonal if true, will include neighbours from diagonal positions in the next and previous layers. * return vector of neighbour cellIDs. */ std::vector neighbours(const CellID& cID, bool aDiagonal) const; /** Calculate layer radii and edges in z-axis, then define cell edges in each layer using defineCellEdges(). + * Following member variables are calculated: + * m_radii + * m_layerEdges + * m_layerDepth + * m_thetaBins (updated through defineCellEdges()) + * m_cellEdges (updated through defineCellEdges()) */ void defineCellsInRZplan() const; @@ -74,12 +84,12 @@ namespace dd4hep { */ inline int phiBins() const { return m_phiBins; } - /** Get the coordinate offset in azimuthal angle. + /** Get the coordinate offset in azimuthal angle, which is the center of cell with m_phiID=0 * return The offset in phi. */ inline double offsetPhi() const { return m_offsetPhi; } - /** Get the vector of theta bins (cells) in a give layer. + /** Get the vector of theta bins (cells) in a given layer. */ inline std::vector thetaBins(const uint layer) const { if(m_radii.empty()) defineCellsInRZplan(); @@ -87,27 +97,34 @@ namespace dd4hep { else return std::vector(); } - /** Get the coordinate offset in z-axis. + /** Get the coordinate offset in z-axis. + * Offset is the middle position of the Barrel or each section of the Endcap. + * For the Barrel, the vector size is 1, while for the Endcap - number of section. * return The offset in z. */ inline std::vector offsetZ() const { return m_offsetZ; } - /** Get the z width of the layer. - * return the z width. + /** Get the z length of the layer. + * return the z length. */ inline std::vector widthZ() const { return m_widthZ; } /** Get the coordinate offset in radius. + * Offset is the inner radius of the first layer in the Barrel or in each section of the Endcap. + * For the Barrel, the vector size is 1, while for the Endcap - number of sections. * return the offset in radius. */ inline std::vector offsetR() const { return m_offsetR; } - /** Get the number of layers. + /** Get the number of layers for each different thickness retrieved with dRlayer(). + * For the Barrel, the vector size equals to the number of different thicknesses used to form the layers. + * For the Endcap, the vector size equals to the number of sections in the Endcap times the number of different thicknesses used to form the layers. * return the number of layers. */ inline std::vector numLayers() const { return m_numLayers; } - /** Get the dR of layers. + /** Get the dR (thickness) of layers. + * The size of the vector equals to the number of different thicknesses used to form the layers. * return the dR. */ inline std::vector dRlayer() const { return m_dRlayer; } @@ -117,6 +134,11 @@ namespace dd4hep { */ inline const std::string& fieldNamePhi() const { return m_phiID; } + /** Get the field name for layer. + * return The field name for layer. + */ + inline const std::string& fieldNameLayer() const { return m_layerID; } + /** Determine the minimum and maximum polar angle of HCal cell based on the cellID. * @param[in] aCellId ID of a cell. * return Theta. @@ -138,6 +160,11 @@ namespace dd4hep { */ inline void setOffsetPhi(double offset) { m_offsetPhi = offset; } + /** Set the detector layout (0 = Barrel; 1 = Endcap). + * @param[in] detLayout + */ + inline void setDetLayout(int detLayout) { m_detLayout = detLayout; } + /** Set the coordinate offset in z-axis. * @param[in] offset in z (centre of the layer). */ @@ -178,6 +205,10 @@ namespace dd4hep { double m_offsetPhi; /// the field name used for phi std::string m_phiID; + /// the field name used for layer + std::string m_layerID; + /// the detector layout (0 = Barrel; 1 = Endcap) + int m_detLayout; /// the z offset of middle of the layer std::vector m_offsetZ; /// the z width of the layer diff --git a/detectorSegmentations/src/FCCSWHCalPhiRow_k4geo.cpp b/detectorSegmentations/src/FCCSWHCalPhiRow_k4geo.cpp index 9cabe64ff..f17059abb 100644 --- a/detectorSegmentations/src/FCCSWHCalPhiRow_k4geo.cpp +++ b/detectorSegmentations/src/FCCSWHCalPhiRow_k4geo.cpp @@ -18,6 +18,7 @@ FCCSWHCalPhiRow_k4geo::FCCSWHCalPhiRow_k4geo(const std::string& cellEncoding) : registerParameter("offset_phi", "Angular offset in phi", m_offsetPhi, 0., SegmentationParameter::AngleUnit, true); registerParameter("grid_size_row", "Number of rows combined in a cell", m_gridSizeRow, std::vector()); registerParameter("dz_row", "dz of row", m_dz_row, 0.); + registerParameter("detLayout", "The detector layout (0 = Barrel; 1 = Endcap)", m_detLayout, -1); registerParameter("offset_z", "Offset in z-axis of the layer center", m_offsetZ, std::vector()); registerParameter("width_z", "Width in z of the layer", m_widthZ, std::vector()); registerParameter("offset_r", "Offset in radius of the layer (Rmin)", m_offsetR, std::vector()); @@ -25,6 +26,7 @@ FCCSWHCalPhiRow_k4geo::FCCSWHCalPhiRow_k4geo(const std::string& cellEncoding) : registerParameter("dRlayer", "dR of the layer", m_dRlayer, std::vector()); registerIdentifier("identifier_phi", "Cell ID identifier for phi", m_phiID, "phi"); registerIdentifier("identifier_row", "Cell ID identifier for row", m_rowID, "row"); + registerIdentifier("identifier_layer", "Cell ID identifier for layer", m_layerID, "layer"); } FCCSWHCalPhiRow_k4geo::FCCSWHCalPhiRow_k4geo(const BitFieldCoder* decoder) : Segmentation(decoder) { @@ -37,6 +39,7 @@ FCCSWHCalPhiRow_k4geo::FCCSWHCalPhiRow_k4geo(const BitFieldCoder* decoder) : Seg registerParameter("offset_phi", "Angular offset in phi", m_offsetPhi, 0., SegmentationParameter::AngleUnit, true); registerParameter("grid_size_row", "Number of rows combined in a cell", m_gridSizeRow, std::vector()); registerParameter("dz_row", "dz of row", m_dz_row, 0.); + registerParameter("detLayout", "The detector layout (0 = Barrel; 1 = Endcap)", m_detLayout, -1); registerParameter("offset_z", "Offset in z-axis of the layer center", m_offsetZ, std::vector()); registerParameter("width_z", "Width in z of the layer", m_widthZ, std::vector()); registerParameter("offset_r", "Offset in radius of the layer (Rmin)", m_offsetR, std::vector()); @@ -44,12 +47,13 @@ FCCSWHCalPhiRow_k4geo::FCCSWHCalPhiRow_k4geo(const BitFieldCoder* decoder) : Seg registerParameter("dRlayer", "dR of the layer", m_dRlayer, std::vector()); registerIdentifier("identifier_phi", "Cell ID identifier for phi", m_phiID, "phi"); registerIdentifier("identifier_row", "Cell ID identifier for row", m_rowID, "row"); + registerIdentifier("identifier_layer", "Cell ID identifier for layer", m_layerID, "layer"); } /// determine the global position based on the cell ID Vector3D FCCSWHCalPhiRow_k4geo::position(const CellID& cID) const { - uint layer = _decoder->get(cID,"layer"); + uint layer = _decoder->get(cID,m_layerID); if(m_radii.empty()) calculateLayerRadii(); if(m_radii.empty() || m_layerEdges.empty()) @@ -62,7 +66,7 @@ Vector3D FCCSWHCalPhiRow_k4geo::position(const CellID& cID) const { double minLayerZ = m_layerEdges[layer].first; // get index of the cell in the layer (index starts from 1!) - int idx = _decoder->get(cID, "row"); + int idx = _decoder->get(cID, m_rowID); // calculate z-coordinate of the cell center double zpos = minLayerZ + (idx-1) * m_dz_row * m_gridSizeRow[layer] + 0.5 * m_dz_row * m_gridSizeRow[layer]; @@ -77,76 +81,33 @@ void FCCSWHCalPhiRow_k4geo::calculateLayerRadii() const { if(m_radii.empty()) { // check if all necessary variables are available - if(m_offsetZ.empty() || m_widthZ.empty() || m_offsetR.empty() || m_numLayers.empty() || m_dRlayer.empty()) + if(m_detLayout==-1 || m_offsetZ.empty() || m_widthZ.empty() || m_offsetR.empty() || m_numLayers.empty() || m_dRlayer.empty()) { dd4hep::printout(dd4hep::ERROR, "FCCSWHCalPhiRow_k4geo","Please check the readout description in the XML file!", - "One of the variables is missing: offset_z | width_z | offset_r | numLayers | dRlayer"); + "One of the variables is missing: detLayout | offset_z | width_z | offset_r | numLayers | dRlayer"); return; } - // calculate the radius for each layer - if(m_offsetZ.size() == 1) // Barrel - { - dd4hep::printout(dd4hep::INFO, "FCCSWHCalPhiRow_k4geo","Barrel configuration found!"); - uint N_dR = m_numLayers.size(); - double moduleDepth = 0.; - for(uint i_dR = 0; i_dR < N_dR; i_dR++) - { - for(int i_row = 1; i_row <= m_numLayers[i_dR]; i_row++) - { - moduleDepth+=m_dRlayer[i_dR]; - m_radii.push_back(m_offsetR[0] + moduleDepth - m_dRlayer[i_dR]*0.5); - // layer lower and upper edges in z-axis - m_layerEdges.push_back( std::make_pair(m_offsetZ[0] - 0.5*m_widthZ[0], m_offsetZ[0] + 0.5*m_widthZ[0]) ); - m_layerDepth.push_back(m_dRlayer[i_dR]); - } - } - } // Barrel + if(m_detLayout==0) dd4hep::printout(dd4hep::INFO, "FCCSWHCalPhiRow_k4geo","Barrel configuration found!"); + else dd4hep::printout(dd4hep::INFO, "FCCSWHCalPhiRow_k4geo","EndCap configuration found!"); - if(m_offsetZ.size() > 1) // ThreePartsEndCap + // calculate the radius for each layer + uint N_dR = m_numLayers.size()/m_offsetZ.size(); + std::vector moduleDepth(m_offsetZ.size()); + for(uint i_section = 0; i_section < m_offsetZ.size(); i_section++) { - dd4hep::printout(dd4hep::INFO, "FCCSWHCalPhiRow_k4geo","ThreePartsEndCap configuration found!"); - uint N_dR = m_numLayers.size()/3; - double moduleDepth1 = 0.; - double moduleDepth2 = 0.; - double moduleDepth3 = 0.; - // part1 - for(uint i_dR = 0; i_dR < N_dR; i_dR++) - { - for(int i_row = 1; i_row <= m_numLayers[i_dR]; i_row++) - { - moduleDepth1+=m_dRlayer[i_dR]; - m_radii.push_back(m_offsetR[0] + moduleDepth1 - m_dRlayer[i_dR]*0.5); - // layer lower and upper edges in z-axis - m_layerEdges.push_back( std::make_pair(m_offsetZ[0] - 0.5*m_widthZ[0], m_offsetZ[0] + 0.5*m_widthZ[0]) ); - m_layerDepth.push_back(m_dRlayer[i_dR]); - } - } - // part2 for(uint i_dR = 0; i_dR < N_dR; i_dR++) { - for(int i_row = 1; i_row <= m_numLayers[i_dR + N_dR]; i_row++) + for(int i_row = 1; i_row <= m_numLayers[i_dR+i_section*N_dR]; i_row++) { - moduleDepth2+=m_dRlayer[i_dR]; - m_radii.push_back(m_offsetR[1] + moduleDepth2 - m_dRlayer[i_dR]*0.5); + moduleDepth[i_section]+=m_dRlayer[i_dR]; + m_radii.push_back(m_offsetR[i_section] + moduleDepth[i_section] - m_dRlayer[i_dR]*0.5); // layer lower and upper edges in z-axis - m_layerEdges.push_back( std::make_pair(m_offsetZ[1] - 0.5*m_widthZ[1], m_offsetZ[1] + 0.5*m_widthZ[1]) ); + m_layerEdges.push_back( std::make_pair(m_offsetZ[i_section] - 0.5*m_widthZ[i_section], m_offsetZ[i_section] + 0.5*m_widthZ[i_section]) ); m_layerDepth.push_back(m_dRlayer[i_dR]); } } - // part3 - for(uint i_dR = 0; i_dR < N_dR; i_dR++) - { - for(int i_row = 1; i_row <= m_numLayers[i_dR + 2*N_dR]; i_row++) - { - moduleDepth3+=m_dRlayer[i_dR]; - m_radii.push_back(m_offsetR[2] + moduleDepth3 - m_dRlayer[i_dR]*0.5); - // layer lower and upper edges in z-axis - m_layerEdges.push_back( std::make_pair(m_offsetZ[2] - 0.5*m_widthZ[2], m_offsetZ[2] + 0.5*m_widthZ[2]) ); - m_layerDepth.push_back(m_dRlayer[i_dR]); - } - } - } // ThreePartsEndcap + } // print info of calculated radii and edges for(uint i_layer = 0; i_layer < m_radii.size(); i_layer++){ @@ -191,7 +152,7 @@ void FCCSWHCalPhiRow_k4geo::defineCellIndexes(const uint layer) const { } // for the EndCap, do it again but for negative-z part - if(m_offsetZ.size() > 1) + if(m_detLayout == 1) { irow = 0; while( (minLayerZ + (irow+1) * m_dz_row) < (maxLayerZ+0.0001) ) @@ -230,9 +191,9 @@ CellID FCCSWHCalPhiRow_k4geo::cellID(const Vector3D& /* localPosition */, const const VolumeID& vID) const { // get the row number from volumeID (starts from 0!) - int nrow = _decoder->get(vID, "row"); + int nrow = _decoder->get(vID, m_rowID); // get the layer number from volumeID - uint layer = _decoder->get(vID,"layer"); + uint layer = _decoder->get(vID,m_layerID); CellID cID = vID; @@ -240,14 +201,14 @@ CellID FCCSWHCalPhiRow_k4geo::cellID(const Vector3D& /* localPosition */, const int idx = floor(nrow/m_gridSizeRow[layer]) + 1; // if the hit is in the negative-z part of the Endcap then assign negative index - if(m_offsetZ.size() > 1 && globalPosition.z() < 0) idx *= -1; + if(m_detLayout == 1 && globalPosition.z() < 0) idx *= -1; _decoder->set(cID, m_rowID, idx); _decoder->set(cID, m_phiID, positionToBin(dd4hep::DDSegmentation::Util::phiFromXYZ(globalPosition), 2 * M_PI / (double)m_phiBins, m_offsetPhi)); // For endcap, the volume ID comes with "type" field information which would screw up the topo-clustering, // therefore, lets set it to zero, as it is for the cell IDs in the neighbours map. - if(m_offsetZ.size() > 1) _decoder->set(cID, "type", 0); + if(m_detLayout == 1) _decoder->set(cID, "type", 0); return cID; } @@ -266,33 +227,27 @@ std::vector > FCCSWHCalPhiRow_k4geo::getMinMaxLayerId() con if(m_radii.empty()) calculateLayerRadii(); if(m_radii.empty()) return minMaxLayerId; + std::vector minLayerId(m_offsetZ.size(), 0); + std::vector maxLayerId(m_offsetZ.size(), 0); - std::vector minLayerIdEndcap = {0, 0, 0}; - std::vector maxLayerIdEndcap = {-1, 0, 0}; - if(m_offsetZ.size() > 1) - { - uint Nl = m_numLayers.size()/3; + uint Nl = m_numLayers.size()/m_offsetZ.size(); - // count the number of layers in the first part of the Endcap - for(uint i=0; i < Nl; i++) maxLayerIdEndcap[0] += m_numLayers[i]; + for(uint i_section = 0; i_section < m_offsetZ.size(); i_section++) + { + if(i_section > 0) + { + minLayerId[i_section] = maxLayerId[i_section-1] + 1; + maxLayerId[i_section] = maxLayerId[i_section-1]; + } - minLayerIdEndcap[1] = maxLayerIdEndcap[0]+1; - maxLayerIdEndcap[1] = maxLayerIdEndcap[0]; - // count the number of layers in the second part of the Endcap - for(uint i=0; i < Nl; i++) maxLayerIdEndcap[1] += m_numLayers[i+Nl]; + for(uint i=0; i < Nl; i++) + { + maxLayerId[i_section] += m_numLayers[i+i_section*Nl]; + } - minLayerIdEndcap[2] = maxLayerIdEndcap[1]+1; - maxLayerIdEndcap[2] = maxLayerIdEndcap[1]; - // count the number of layers in the third part of the Endcap - for(uint i=0; i < Nl; i++) maxLayerIdEndcap[2] += m_numLayers[i+2*Nl]; + if(i_section==0) maxLayerId[0] -= 1; - minMaxLayerId.push_back(std::make_pair(minLayerIdEndcap[0],maxLayerIdEndcap[0])); - minMaxLayerId.push_back(std::make_pair(minLayerIdEndcap[1],maxLayerIdEndcap[1])); - minMaxLayerId.push_back(std::make_pair(minLayerIdEndcap[2],maxLayerIdEndcap[2])); - } - else // for Barrel - { - minMaxLayerId.push_back(std::make_pair(0,m_radii.size()-1)); + minMaxLayerId.push_back(std::make_pair(minLayerId[i_section], maxLayerId[i_section])); } return minMaxLayerId; @@ -306,14 +261,11 @@ std::vector FCCSWHCalPhiRow_k4geo::neighbours(const CellID& cID) const if(m_radii.empty()) calculateLayerRadii(); if(m_cellIndexes.empty()) return cellNeighbours; - bool EndcapPart1 = false; - bool EndcapPart2 = false; - bool EndcapPart3 = false; - + uint EndcapPart = 0; int minLayerId = -1; int maxLayerId = -1; - int currentLayerId = _decoder->get(cID,"layer"); + int currentLayerId = _decoder->get(cID,m_layerID); int currentCellId = _decoder->get(cID,m_rowID); int minCellId = m_cellIndexes[currentLayerId].front(); @@ -322,54 +274,29 @@ std::vector FCCSWHCalPhiRow_k4geo::neighbours(const CellID& cID) const //-------------------------------- // Determine min and max layer Id //-------------------------------- - // if this is the segmentation of three parts Endcap - /* - * hardcoded numbers would be the following: - * std::vector minLayerIdEndcap = {0, 6, 15}; - * std::vector maxLayerIdEndcap = {5, 14, 36}; - * but lets try to avoid hardcoding: - */ - std::vector minLayerIdEndcap = {0, 0, 0}; - std::vector maxLayerIdEndcap = {0, 0, 0}; - if(m_offsetZ.size() > 1) + std::vector minLayerIdEndcap(m_offsetZ.size(),0); + std::vector maxLayerIdEndcap(m_offsetZ.size(),0); + if(m_detLayout == 1) // Endcap { std::vector > minMaxLayerId(getMinMaxLayerId()); if(minMaxLayerId.empty()) { - dd4hep::printout(dd4hep::ERROR, "FCCSWHCalPhiRow_k4geo","Can not get ThreePartsEndcap min and max layer indexes! --> returning empty neighbours"); + dd4hep::printout(dd4hep::ERROR, "FCCSWHCalPhiRow_k4geo","Can not get Endcap min and max layer indexes! --> returning empty neighbours"); return cellNeighbours; } - // part1 min and max layer index - minLayerIdEndcap[0] = minMaxLayerId[0].first; - maxLayerIdEndcap[0] = minMaxLayerId[0].second; - // part2 min and max layer index - minLayerIdEndcap[1] = minMaxLayerId[1].first; - maxLayerIdEndcap[1] = minMaxLayerId[1].second; - // part3 min and max layer index - minLayerIdEndcap[2] = minMaxLayerId[2].first; - maxLayerIdEndcap[2] = minMaxLayerId[2].second; - - // Part 1 - if(currentLayerId >= minLayerIdEndcap[0] && currentLayerId <= maxLayerIdEndcap[0]) + // determine min and max layer Ids and which section/part it is + for(uint i_section = 0; i_section < minMaxLayerId.size(); i_section++) { - minLayerId = minLayerIdEndcap[0]; - maxLayerId = maxLayerIdEndcap[0]; - EndcapPart1 = true; - } - // Part 2 - if(currentLayerId >= minLayerIdEndcap[1] && currentLayerId <= maxLayerIdEndcap[1]) - { - minLayerId = minLayerIdEndcap[1]; - maxLayerId = maxLayerIdEndcap[1]; - EndcapPart2 = true; - } - // Part 3 - if(currentLayerId >= minLayerIdEndcap[2] && currentLayerId <= maxLayerIdEndcap[2]) - { - minLayerId = minLayerIdEndcap[2]; - maxLayerId = maxLayerIdEndcap[2]; - EndcapPart3 = true; + minLayerIdEndcap[i_section] = minMaxLayerId[i_section].first; + maxLayerIdEndcap[i_section] = minMaxLayerId[i_section].second; + + if(currentLayerId >= minLayerIdEndcap[i_section] && currentLayerId <= maxLayerIdEndcap[i_section]) + { + minLayerId = minLayerIdEndcap[i_section]; + maxLayerId = maxLayerIdEndcap[i_section]; + EndcapPart = i_section; + } } // correct the min and max CellId for endcap @@ -403,7 +330,7 @@ std::vector FCCSWHCalPhiRow_k4geo::neighbours(const CellID& cID) const { CellID nID = cID ; int prevLayerId = currentLayerId - 1; - _decoder->set(nID,"layer",prevLayerId); + _decoder->set(nID,m_layerID,prevLayerId); // if the granularity is the same for the previous layer then take the cells with currentCellId, currentCellId - 1, and currentCellId + 1 if(m_gridSizeRow[prevLayerId] == m_gridSizeRow[currentLayerId]) @@ -450,7 +377,7 @@ std::vector FCCSWHCalPhiRow_k4geo::neighbours(const CellID& cID) const { CellID nID = cID ; int nextLayerId = currentLayerId + 1; - _decoder->set(nID,"layer",nextLayerId); + _decoder->set(nID,m_layerID,nextLayerId); // if the granularity is the same for the next layer then take the cells with currentCellId, currentCellId - 1, and currentCellId + 1 if(m_gridSizeRow[nextLayerId] == m_gridSizeRow[currentLayerId]) @@ -507,8 +434,8 @@ std::vector FCCSWHCalPhiRow_k4geo::neighbours(const CellID& cID) const cellNeighbours.push_back(nID); // add the next cell from current layer of the same phi module } - // if this is the Endcap then look for neighbours in different parts as well - if(m_offsetZ.size() > 1) + // if this is the Endcap (and consists of more than 1 part/section) then look for neighbours in different parts as well + if(m_detLayout == 1 && m_offsetZ.size() > 1) { double currentLayerRmin = m_radii[currentLayerId] - 0.5*m_layerDepth[currentLayerId]; double currentLayerRmax = m_radii[currentLayerId] + 0.5*m_layerDepth[currentLayerId]; @@ -520,8 +447,8 @@ std::vector FCCSWHCalPhiRow_k4geo::neighbours(const CellID& cID) const maxCellId = m_cellIndexes[currentLayerId].back(); // this should be -N } - // if it is the last cell in the part1 - if(EndcapPart1 && currentCellId == maxCellId ) + // if it is the last cell in the first part + if(EndcapPart == 0 && currentCellId == maxCellId ) { // find the layers in the part2 that share a border with the current layer for(int part2layerId = minLayerIdEndcap[1]; part2layerId <= maxLayerIdEndcap[1]; part2layerId++) @@ -536,67 +463,72 @@ std::vector FCCSWHCalPhiRow_k4geo::neighbours(const CellID& cID) const ) { CellID nID = cID ; - _decoder->set(nID,"layer",part2layerId); + _decoder->set(nID,m_layerID,part2layerId); _decoder->set(nID,m_rowID,minCellId); cellNeighbours.push_back(nID); // add the first cell from part2 layer } } } - // if it is the first cell in the part2 - if(EndcapPart2 && currentCellId == minCellId) + // if the Endcap consists of more than 2 parts: + for(uint i_section = 1; i_section < (m_offsetZ.size()-1); i_section++) { - // find the layers in part1 that share a border with the current layer - for(int part1layerId = minLayerIdEndcap[0]; part1layerId <= maxLayerIdEndcap[0]; part1layerId++) - { - double Rmin = m_radii[part1layerId] - 0.5*m_layerDepth[part1layerId]; - double Rmax = m_radii[part1layerId] + 0.5*m_layerDepth[part1layerId]; + if(i_section != EndcapPart) continue; - if( (Rmin >= currentLayerRmin && Rmin <= currentLayerRmax) - || (Rmax >= currentLayerRmin && Rmax <= currentLayerRmax) - || (currentLayerRmin >= Rmin && currentLayerRmin <= Rmax) - || (currentLayerRmax >= Rmin && currentLayerRmax <= Rmax) - ) + // if it is the first cell then look for neighbours in previous part + if(currentCellId == minCellId) + { + // find the layers in previous part that share a border with the current layer + for(int prevPartLayerId = minLayerIdEndcap[i_section-1]; prevPartLayerId <= maxLayerIdEndcap[i_section-1]; prevPartLayerId++) { - CellID nID = cID ; - _decoder->set(nID,"layer",part1layerId); - _decoder->set(nID,m_rowID,maxCellId); - cellNeighbours.push_back(nID); // add the last cell from the part1 layer + double Rmin = m_radii[prevPartLayerId] - 0.5*m_layerDepth[prevPartLayerId]; + double Rmax = m_radii[prevPartLayerId] + 0.5*m_layerDepth[prevPartLayerId]; + + if( (Rmin >= currentLayerRmin && Rmin <= currentLayerRmax) + || (Rmax >= currentLayerRmin && Rmax <= currentLayerRmax) + || (currentLayerRmin >= Rmin && currentLayerRmin <= Rmax) + || (currentLayerRmax >= Rmin && currentLayerRmax <= Rmax) + ) + { + CellID nID = cID ; + _decoder->set(nID,m_layerID,prevPartLayerId); + _decoder->set(nID,m_rowID,maxCellId); + cellNeighbours.push_back(nID); // add the last cell from the previous part layer + } } } - } - - // if it is the last cell in the part2 - if(EndcapPart2 && currentCellId == maxCellId) - { - // find the layers in the part3 that share a border with the current layer - for(int part3layerId = minLayerIdEndcap[2]; part3layerId <= maxLayerIdEndcap[2]; part3layerId++) + // if it is the last cell then look for neighbours in the next part + if(currentCellId == maxCellId) { - double Rmin = m_radii[part3layerId] - 0.5*m_layerDepth[part3layerId]; - double Rmax = m_radii[part3layerId] + 0.5*m_layerDepth[part3layerId]; - - if( (Rmin >= currentLayerRmin && Rmin <= currentLayerRmax) - || (Rmax >= currentLayerRmin && Rmax <= currentLayerRmax) - || (currentLayerRmin >= Rmin && currentLayerRmin <= Rmax) - || (currentLayerRmax >= Rmin && currentLayerRmax <= Rmax) - ) + // find the layers in the next that share a border with the current layer + for(int nextPartLayerId = minLayerIdEndcap[i_section+1]; nextPartLayerId <= maxLayerIdEndcap[i_section+1]; nextPartLayerId++) { - CellID nID = cID ; - _decoder->set(nID,"layer",part3layerId); - _decoder->set(nID,m_rowID,minCellId); - cellNeighbours.push_back(nID); // add the first cell from the part3 layer + double Rmin = m_radii[nextPartLayerId] - 0.5*m_layerDepth[nextPartLayerId]; + double Rmax = m_radii[nextPartLayerId] + 0.5*m_layerDepth[nextPartLayerId]; + + if( (Rmin >= currentLayerRmin && Rmin <= currentLayerRmax) + || (Rmax >= currentLayerRmin && Rmax <= currentLayerRmax) + || (currentLayerRmin >= Rmin && currentLayerRmin <= Rmax) + || (currentLayerRmax >= Rmin && currentLayerRmax <= Rmax) + ) + { + CellID nID = cID ; + _decoder->set(nID,m_layerID,nextPartLayerId); + _decoder->set(nID,m_rowID,minCellId); + cellNeighbours.push_back(nID); // add the first cell from the next part layer + } } } } - // if it is the first cell in the part3 - if(EndcapPart3 && currentCellId == minCellId) + // if it is the first cell in the last part of the Endcap + if(EndcapPart == (m_offsetZ.size() - 1) && currentCellId == minCellId) { - // find the layers in the part2 that share a border with the current layer - for(int part2layerId = minLayerIdEndcap[1]; part2layerId <= maxLayerIdEndcap[1]; part2layerId++) + // find the layers in the previous part that share a border with the current layer + for(int prevPartLayerId = minLayerIdEndcap[m_offsetZ.size()-2]; prevPartLayerId <= maxLayerIdEndcap[m_offsetZ.size()-2]; prevPartLayerId++) { - double Rmin = m_radii[part2layerId] - 0.5*m_layerDepth[part2layerId]; - double Rmax = m_radii[part2layerId] + 0.5*m_layerDepth[part2layerId]; + double Rmin = m_radii[prevPartLayerId] - 0.5*m_layerDepth[prevPartLayerId]; + double Rmax = m_radii[prevPartLayerId] + 0.5*m_layerDepth[prevPartLayerId]; if( (Rmin >= currentLayerRmin && Rmin <= currentLayerRmax) || (Rmax >= currentLayerRmin && Rmax <= currentLayerRmax) @@ -605,7 +537,7 @@ std::vector FCCSWHCalPhiRow_k4geo::neighbours(const CellID& cID) const ) { CellID nID = cID ; - _decoder->set(nID,"layer",part2layerId); + _decoder->set(nID,m_layerID,prevPartLayerId); _decoder->set(nID,m_rowID,maxCellId); cellNeighbours.push_back(nID); // add the last cell from the part2 layer } @@ -645,7 +577,7 @@ std::array FCCSWHCalPhiRow_k4geo::cellTheta(const CellID& cID) const // get the cell index int idx = _decoder->get(cID, m_rowID); // get the layer index - uint layer = _decoder->get(cID,"layer"); + uint layer = _decoder->get(cID,m_layerID); if(m_radii.empty()) calculateLayerRadii(); if(m_cellEdges.empty()) return cTheta; diff --git a/detectorSegmentations/src/FCCSWHCalPhiTheta_k4geo.cpp b/detectorSegmentations/src/FCCSWHCalPhiTheta_k4geo.cpp index 04c227ab1..a6248a3d3 100644 --- a/detectorSegmentations/src/FCCSWHCalPhiTheta_k4geo.cpp +++ b/detectorSegmentations/src/FCCSWHCalPhiTheta_k4geo.cpp @@ -13,12 +13,14 @@ FCCSWHCalPhiTheta_k4geo::FCCSWHCalPhiTheta_k4geo(const std::string& cellEncoding // register all necessary parameters (additional to those registered in GridTheta_k4geo) registerParameter("phi_bins", "Number of bins phi", m_phiBins, 1); registerParameter("offset_phi", "Angular offset in phi", m_offsetPhi, 0., SegmentationParameter::AngleUnit, true); + registerParameter("detLayout", "The detector layout (0 = Barrel; 1 = Endcap)", m_detLayout, -1); registerParameter("offset_z", "Offset in z-axis of the layer center", m_offsetZ, std::vector()); registerParameter("width_z", "Width in z of the layer", m_widthZ, std::vector()); registerParameter("offset_r", "Offset in radius of the layer (Rmin)", m_offsetR, std::vector()); registerParameter("numLayers", "Number of layers", m_numLayers, std::vector()); registerParameter("dRlayer", "dR of the layer", m_dRlayer, std::vector()); registerIdentifier("identifier_phi", "Cell ID identifier for phi", m_phiID, "phi"); + registerIdentifier("identifier_layer", "Cell ID identifier for layer", m_layerID, "layer"); } FCCSWHCalPhiTheta_k4geo::FCCSWHCalPhiTheta_k4geo(const BitFieldCoder* decoder) : GridTheta_k4geo(decoder) { @@ -29,18 +31,20 @@ FCCSWHCalPhiTheta_k4geo::FCCSWHCalPhiTheta_k4geo(const BitFieldCoder* decoder) : // register all necessary parameters (additional to those registered in GridTheta_k4geo) registerParameter("phi_bins", "Number of bins phi", m_phiBins, 1); registerParameter("offset_phi", "Angular offset in phi", m_offsetPhi, 0., SegmentationParameter::AngleUnit, true); + registerParameter("detLayout", "The detector layout (0 = Barrel; 1 = Endcap)", m_detLayout, -1); registerParameter("offset_z", "Offset in z-axis of the layer center", m_offsetZ, std::vector()); registerParameter("width_z", "Width in z of the layer", m_widthZ, std::vector()); registerParameter("offset_r", "Offset in radius of the layer (Rmin)", m_offsetR, std::vector()); registerParameter("numLayers", "Number of layers", m_numLayers, std::vector()); registerParameter("dRlayer", "dR of the layer", m_dRlayer, std::vector()); registerIdentifier("identifier_phi", "Cell ID identifier for phi", m_phiID, "phi"); + registerIdentifier("identifier_layer", "Cell ID identifier for layer", m_layerID, "layer"); } /** /// determine the global position based on the cell ID Vector3D FCCSWHCalPhiTheta_k4geo::position(const CellID& cID) const { - uint layer = _decoder->get(cID,"layer"); + uint layer = _decoder->get(cID,m_layerID); double radius = 1.0; if(m_radii.empty()) defineCellsInRZplan(); @@ -53,8 +57,8 @@ Vector3D FCCSWHCalPhiTheta_k4geo::position(const CellID& cID) const { /// determine the global position based on the cell ID /// returns the geometric center of the cell Vector3D FCCSWHCalPhiTheta_k4geo::position(const CellID& cID) const { - uint layer = _decoder->get(cID,"layer"); - int thetaID = _decoder->get(cID,"theta"); + uint layer = _decoder->get(cID,m_layerID); + int thetaID = _decoder->get(cID,m_thetaID); double zpos = 0.; double radius = 1.0; @@ -73,76 +77,33 @@ void FCCSWHCalPhiTheta_k4geo::defineCellsInRZplan() const { if(m_radii.empty()) { // check if all necessary variables are available - if(m_offsetZ.empty() || m_widthZ.empty() || m_offsetR.empty() || m_numLayers.empty() || m_dRlayer.empty()) + if(m_detLayout==-1 || m_offsetZ.empty() || m_widthZ.empty() || m_offsetR.empty() || m_numLayers.empty() || m_dRlayer.empty()) { - dd4hep::printout(dd4hep::ERROR, "FCCSWHCalPhiTheta_k4geo","Please check the readout description in the XML file!", - "One of the variables is missing: offset_z | width_z | offset_r | numLayers | dRlayer"); + dd4hep::printout(dd4hep::ERROR, "FCCSWHCalPhiRow_k4geo","Please check the readout description in the XML file!", + "One of the variables is missing: detLayout | offset_z | width_z | offset_r | numLayers | dRlayer"); return; } - // calculate the radius for each layer - if(m_offsetZ.size() == 1) // Barrel - { - dd4hep::printout(dd4hep::INFO, "FCCSWHCalPhiTheta_k4geo","Barrel configuration found!"); - uint N_dR = m_numLayers.size(); - double moduleDepth = 0.; - for(uint i_dR = 0; i_dR < N_dR; i_dR++) - { - for(int i_row = 1; i_row <= m_numLayers[i_dR]; i_row++) - { - moduleDepth+=m_dRlayer[i_dR]; - m_radii.push_back(m_offsetR[0] + moduleDepth - m_dRlayer[i_dR]*0.5); - // layer lower and upper edges in z-axis - m_layerEdges.push_back( std::make_pair(m_offsetZ[0] - 0.5*m_widthZ[0], m_offsetZ[0] + 0.5*m_widthZ[0]) ); - m_layerDepth.push_back(m_dRlayer[i_dR]); - } - } - } // Barrel + if(m_detLayout==0) dd4hep::printout(dd4hep::INFO, "FCCSWHCalPhiRow_k4geo","Barrel configuration found!"); + else dd4hep::printout(dd4hep::INFO, "FCCSWHCalPhiRow_k4geo","EndCap configuration found!"); - if(m_offsetZ.size() > 1) // ThreePartsEndCap + // calculate the radius for each layer + uint N_dR = m_numLayers.size()/m_offsetZ.size(); + std::vector moduleDepth(m_offsetZ.size()); + for(uint i_section = 0; i_section < m_offsetZ.size(); i_section++) { - dd4hep::printout(dd4hep::INFO, "FCCSWHCalPhiTheta_k4geo","ThreePartsEndCap configuration found!"); - uint N_dR = m_numLayers.size()/3; - double moduleDepth1 = 0.; - double moduleDepth2 = 0.; - double moduleDepth3 = 0.; - // part1 for(uint i_dR = 0; i_dR < N_dR; i_dR++) { - for(int i_row = 1; i_row <= m_numLayers[i_dR]; i_row++) + for(int i_row = 1; i_row <= m_numLayers[i_dR+i_section*N_dR]; i_row++) { - moduleDepth1+=m_dRlayer[i_dR]; - m_radii.push_back(m_offsetR[0] + moduleDepth1 - m_dRlayer[i_dR]*0.5); + moduleDepth[i_section]+=m_dRlayer[i_dR]; + m_radii.push_back(m_offsetR[i_section] + moduleDepth[i_section] - m_dRlayer[i_dR]*0.5); // layer lower and upper edges in z-axis - m_layerEdges.push_back( std::make_pair(m_offsetZ[0] - 0.5*m_widthZ[0], m_offsetZ[0] + 0.5*m_widthZ[0]) ); + m_layerEdges.push_back( std::make_pair(m_offsetZ[i_section] - 0.5*m_widthZ[i_section], m_offsetZ[i_section] + 0.5*m_widthZ[i_section]) ); m_layerDepth.push_back(m_dRlayer[i_dR]); } } - // part2 - for(uint i_dR = 0; i_dR < N_dR; i_dR++) - { - for(int i_row = 1; i_row <= m_numLayers[i_dR + N_dR]; i_row++) - { - moduleDepth2+=m_dRlayer[i_dR]; - m_radii.push_back(m_offsetR[1] + moduleDepth2 - m_dRlayer[i_dR]*0.5); - // layer lower and upper edges in z-axis - m_layerEdges.push_back( std::make_pair(m_offsetZ[1] - 0.5*m_widthZ[1], m_offsetZ[1] + 0.5*m_widthZ[1]) ); - m_layerDepth.push_back(m_dRlayer[i_dR]); - } - } - // part3 - for(uint i_dR = 0; i_dR < N_dR; i_dR++) - { - for(int i_row = 1; i_row <= m_numLayers[i_dR + 2*N_dR]; i_row++) - { - moduleDepth3+=m_dRlayer[i_dR]; - m_radii.push_back(m_offsetR[2] + moduleDepth3 - m_dRlayer[i_dR]*0.5); - // layer lower and upper edges in z-axis - m_layerEdges.push_back( std::make_pair(m_offsetZ[2] - 0.5*m_widthZ[2], m_offsetZ[2] + 0.5*m_widthZ[2]) ); - m_layerDepth.push_back(m_dRlayer[i_dR]); - } - } - } // ThreePartsEndcap + } // print info of calculated radii and edges for(uint i_layer = 0; i_layer < m_radii.size(); i_layer++){ @@ -196,7 +157,7 @@ void FCCSWHCalPhiTheta_k4geo::defineCellEdges(const uint layer) const { m_cellEdges[layer][prevBin].first = m_layerEdges[layer].first; // for the EndCap, do it again but for negative z part - if(m_offsetZ.size() > 1) + if(m_detLayout == 1) { while(m_radii[layer]*std::cos(m_offsetTheta+ibin*m_gridSizeTheta)/std::sin(m_offsetTheta+ibin*m_gridSizeTheta) > (-m_layerEdges[layer].second)) { @@ -248,11 +209,11 @@ CellID FCCSWHCalPhiTheta_k4geo::cellID(const Vector3D& /* localPosition */, cons // For endcap, the volume ID comes with "type" field information which would screw up the topo-clustering as the "row" field, // therefore, lets set it to zero, as it is for the cell IDs in the neighbours map. - if(m_offsetZ.size() > 1) _decoder->set(cID, "type", 0); + if(m_detLayout == 1) _decoder->set(cID, "type", 0); double lTheta = thetaFromXYZ(globalPosition); double lPhi = phiFromXYZ(globalPosition); - uint layer = _decoder->get(vID,"layer"); + uint layer = _decoder->get(vID,m_layerID); // define cell boundaries in R-z plan if(m_radii.empty()) defineCellsInRZplan(); @@ -288,45 +249,32 @@ double FCCSWHCalPhiTheta_k4geo::phi(const CellID& cID) const { /// Get the min and max layer indexes for each part of the HCal std::vector > FCCSWHCalPhiTheta_k4geo::getMinMaxLayerId() const { - /* - * hardcoded numbers would be the following: - * std::vector minLayerIdEndcap = {0, 6, 15}; - * std::vector maxLayerIdEndcap = {5, 14, 36}; - * but lets try to avoid hardcoding: - */ - std::vector > minMaxLayerId; if(m_radii.empty()) defineCellsInRZplan(); if(m_radii.empty()) return minMaxLayerId; + std::vector minLayerId(m_offsetZ.size(), 0); + std::vector maxLayerId(m_offsetZ.size(), 0); - std::vector minLayerIdEndcap = {0, 0, 0}; - std::vector maxLayerIdEndcap = {-1, 0, 0}; - if(m_offsetZ.size() > 1) - { - uint Nl = m_numLayers.size()/3; + uint Nl = m_numLayers.size()/m_offsetZ.size(); - // count the number of layers in the first part of the Endcap - for(uint i=0; i < Nl; i++) maxLayerIdEndcap[0] += m_numLayers[i]; + for(uint i_section = 0; i_section < m_offsetZ.size(); i_section++) + { + if(i_section > 0) + { + minLayerId[i_section] = maxLayerId[i_section-1] + 1; + maxLayerId[i_section] = maxLayerId[i_section-1]; + } - minLayerIdEndcap[1] = maxLayerIdEndcap[0]+1; - maxLayerIdEndcap[1] = maxLayerIdEndcap[0]; - // count the number of layers in the second part of the Endcap - for(uint i=0; i < Nl; i++) maxLayerIdEndcap[1] += m_numLayers[i+Nl]; + for(uint i=0; i < Nl; i++) + { + maxLayerId[i_section] += m_numLayers[i+i_section*Nl]; + } - minLayerIdEndcap[2] = maxLayerIdEndcap[1]+1; - maxLayerIdEndcap[2] = maxLayerIdEndcap[1]; - // count the number of layers in the third part of the Endcap - for(uint i=0; i < Nl; i++) maxLayerIdEndcap[2] += m_numLayers[i+2*Nl]; + if(i_section == 0) maxLayerId[0] -= 1; - minMaxLayerId.push_back(std::make_pair(minLayerIdEndcap[0],maxLayerIdEndcap[0])); - minMaxLayerId.push_back(std::make_pair(minLayerIdEndcap[1],maxLayerIdEndcap[1])); - minMaxLayerId.push_back(std::make_pair(minLayerIdEndcap[2],maxLayerIdEndcap[2])); - } - else // for Barrel - { - minMaxLayerId.push_back(std::make_pair(0,m_radii.size()-1)); + minMaxLayerId.push_back(std::make_pair(minLayerId[i_section], maxLayerId[i_section])); } return minMaxLayerId; @@ -340,14 +288,11 @@ std::vector FCCSWHCalPhiTheta_k4geo::neighbours(const CellID& cID, boo if(m_radii.empty()) defineCellsInRZplan(); if(m_thetaBins.empty()) return cellNeighbours; - bool EndcapPart1 = false; - bool EndcapPart2 = false; - bool EndcapPart3 = false; - + uint EndcapPart = 0; int minLayerId = -1; int maxLayerId = -1; - int currentLayerId = _decoder->get(cID,"layer"); + int currentLayerId = _decoder->get(cID,m_layerID); int currentCellThetaBin = _decoder->get(cID,m_thetaID); int minCellThetaBin = m_thetaBins[currentLayerId].front(); @@ -356,48 +301,29 @@ std::vector FCCSWHCalPhiTheta_k4geo::neighbours(const CellID& cID, boo //-------------------------------- // Determine min and max layer Id //-------------------------------- - // if this is the segmentation of three parts Endcap - std::vector minLayerIdEndcap = {0, 0, 0}; - std::vector maxLayerIdEndcap = {0, 0, 0}; - if(m_offsetZ.size() > 1) + std::vector minLayerIdEndcap(m_offsetZ.size(),0); + std::vector maxLayerIdEndcap(m_offsetZ.size(),0); + if(m_detLayout == 1) { std::vector > minMaxLayerId(getMinMaxLayerId()); if(minMaxLayerId.empty()) { - dd4hep::printout(dd4hep::ERROR, "FCCSWHCalPhiTheta_k4geo","Can not get ThreePartsEndcap min and max layer indexes! --> returning empty neighbours"); + dd4hep::printout(dd4hep::ERROR, "FCCSWHCalPhiTheta_k4geo","Can not get Endcap min and max layer indexes! --> returning empty neighbours"); return cellNeighbours; } - // part1 min and max layer index - minLayerIdEndcap[0] = minMaxLayerId[0].first; - maxLayerIdEndcap[0] = minMaxLayerId[0].second; - // part2 min and max layer index - minLayerIdEndcap[1] = minMaxLayerId[1].first; - maxLayerIdEndcap[1] = minMaxLayerId[1].second; - // part3 min and max layer index - minLayerIdEndcap[2] = minMaxLayerId[2].first; - maxLayerIdEndcap[2] = minMaxLayerId[2].second; - - // Part 1 - if(currentLayerId >= minLayerIdEndcap[0] && currentLayerId <= maxLayerIdEndcap[0]) + // determine min and max layer Ids and which section/part it is + for(uint i_section = 0; i_section < minMaxLayerId.size(); i_section++) { - minLayerId = minLayerIdEndcap[0]; - maxLayerId = maxLayerIdEndcap[0]; - EndcapPart1 = true; - } - // Part 2 - if(currentLayerId >= minLayerIdEndcap[1] && currentLayerId <= maxLayerIdEndcap[1]) - { - minLayerId = minLayerIdEndcap[1]; - maxLayerId = maxLayerIdEndcap[1]; - EndcapPart2 = true; - } - // Part 3 - if(currentLayerId >= minLayerIdEndcap[2] && currentLayerId <= maxLayerIdEndcap[2]) - { - minLayerId = minLayerIdEndcap[2]; - maxLayerId = maxLayerIdEndcap[2]; - EndcapPart3 = true; + minLayerIdEndcap[i_section] = minMaxLayerId[i_section].first; + maxLayerIdEndcap[i_section] = minMaxLayerId[i_section].second; + + if(currentLayerId >= minLayerIdEndcap[i_section] && currentLayerId <= maxLayerIdEndcap[i_section]) + { + minLayerId = minLayerIdEndcap[i_section]; + maxLayerId = maxLayerIdEndcap[i_section]; + EndcapPart = i_section; + } } // correct the min and max theta bin for endcap @@ -446,7 +372,7 @@ std::vector FCCSWHCalPhiTheta_k4geo::neighbours(const CellID& cID, boo //---------------------------------------------- // deal with the Barrel - if(m_offsetZ.size() == 1) + if(m_detLayout == 0) { double currentCellZmin = m_cellEdges[currentLayerId][currentCellThetaBin].first; double currentCellZmax = m_cellEdges[currentLayerId][currentCellThetaBin].second; @@ -456,7 +382,7 @@ std::vector FCCSWHCalPhiTheta_k4geo::neighbours(const CellID& cID, boo { CellID nID = cID ; int prevLayerId = currentLayerId - 1; - _decoder->set(nID,"layer",prevLayerId); + _decoder->set(nID,m_layerID,prevLayerId); _decoder->set(nID,m_thetaID,currentCellThetaBin); cellNeighbours.push_back(nID); // add the cell with the same theta bin from the previous layer of the same phi module @@ -502,7 +428,7 @@ std::vector FCCSWHCalPhiTheta_k4geo::neighbours(const CellID& cID, boo { CellID nID = cID ; int nextLayerId = currentLayerId + 1; - _decoder->set(nID,"layer",nextLayerId); + _decoder->set(nID,m_layerID,nextLayerId); _decoder->set(nID,m_thetaID,currentCellThetaBin); cellNeighbours.push_back(nID); // add the cell with the same theta bin from the next layer of the same phi module @@ -548,8 +474,8 @@ std::vector FCCSWHCalPhiTheta_k4geo::neighbours(const CellID& cID, boo } } - // if this is the Endcap then look for neighbours in different parts as well - if(m_offsetZ.size() > 1) + // Endcap + if(m_detLayout == 1) { double currentCellZmin = m_cellEdges[currentLayerId][currentCellThetaBin].first; double currentCellZmax = m_cellEdges[currentLayerId][currentCellThetaBin].second; @@ -559,7 +485,7 @@ std::vector FCCSWHCalPhiTheta_k4geo::neighbours(const CellID& cID, boo { CellID nID = cID ; int prevLayerId = currentLayerId - 1; - _decoder->set(nID,"layer",prevLayerId); + _decoder->set(nID,m_layerID,prevLayerId); // find the ones that share at least part of a border with the current cell for( auto bin : m_thetaBins[prevLayerId] ) { @@ -607,7 +533,7 @@ std::vector FCCSWHCalPhiTheta_k4geo::neighbours(const CellID& cID, boo { CellID nID = cID ; int nextLayerId = currentLayerId + 1; - _decoder->set(nID,"layer",nextLayerId); + _decoder->set(nID,m_layerID,nextLayerId); // find the ones that share at least part of a border with the current cell for( auto bin : m_thetaBins[nextLayerId] ) { @@ -650,135 +576,143 @@ std::vector FCCSWHCalPhiTheta_k4geo::neighbours(const CellID& cID, boo } } - - // - double currentLayerRmin = m_radii[currentLayerId] - 0.5*m_layerDepth[currentLayerId]; - double currentLayerRmax = m_radii[currentLayerId] + 0.5*m_layerDepth[currentLayerId]; - - - // if the cell is in negative-z part, then swap min and max theta bins - if(theta(cID) > M_PI/2.) + // if the Endcap consists of more than 1 part/section then look for neighbours in different parts as well + if(m_offsetZ.size() > 1) { - minCellThetaBin = m_thetaBins[currentLayerId].back(); - maxCellThetaBin = m_thetaBins[currentLayerId][m_thetaBins[currentLayerId].size()/2]; - } + // + double currentLayerRmin = m_radii[currentLayerId] - 0.5*m_layerDepth[currentLayerId]; + double currentLayerRmax = m_radii[currentLayerId] + 0.5*m_layerDepth[currentLayerId]; - // if it is the last cell in the part1 - if(EndcapPart1 && currentCellThetaBin == minCellThetaBin ) - { - // find the layers in the part2 that share a border with the current layer - for(int part2layerId = minLayerIdEndcap[1]; part2layerId <= maxLayerIdEndcap[1]; part2layerId++) + // if the cell is in negative-z part, then swap min and max theta bins + if(theta(cID) > M_PI/2.) { - double Rmin = m_radii[part2layerId] - 0.5*m_layerDepth[part2layerId]; - double Rmax = m_radii[part2layerId] + 0.5*m_layerDepth[part2layerId]; - - if( (Rmin >= currentLayerRmin && Rmin <= currentLayerRmax) - || (Rmax > currentLayerRmin && Rmax <= currentLayerRmax) - || (currentLayerRmin >= Rmin && currentLayerRmin < Rmax) - || (currentLayerRmax >= Rmin && currentLayerRmax <= Rmax) - ) - { - CellID nID = cID ; - _decoder->set(nID,"layer",part2layerId); - _decoder->set(nID,m_thetaID,maxCellThetaBin); - cellNeighbours.push_back(nID); // add the last theta bin cell from part2 layer - } - if(aDiagonal && Rmax == currentLayerRmin) - { - CellID nID = cID ; - _decoder->set(nID,"layer",part2layerId); - _decoder->set(nID,m_thetaID,maxCellThetaBin); - cellNeighbours.push_back(nID); // add the last theta bin cell from part2 layer - } + minCellThetaBin = m_thetaBins[currentLayerId].back(); + maxCellThetaBin = m_thetaBins[currentLayerId][m_thetaBins[currentLayerId].size()/2]; } - } - // if it is the last theta bin cell in the part2 - if(EndcapPart2 && currentCellThetaBin == maxCellThetaBin) - { - // find the layers in part1 that share a border with the current layer - for(int part1layerId = minLayerIdEndcap[0]; part1layerId <= maxLayerIdEndcap[0]; part1layerId++) + // if it is the last cell in the part1 + if(EndcapPart == 0 && currentCellThetaBin == minCellThetaBin ) { - double Rmin = m_radii[part1layerId] - 0.5*m_layerDepth[part1layerId]; - double Rmax = m_radii[part1layerId] + 0.5*m_layerDepth[part1layerId]; - - if( (Rmin >= currentLayerRmin && Rmin < currentLayerRmax) - || (Rmax >= currentLayerRmin && Rmax <= currentLayerRmax) - || (currentLayerRmin >= Rmin && currentLayerRmin <= Rmax) - || (currentLayerRmax > Rmin && currentLayerRmax <= Rmax) - ) - { - CellID nID = cID ; - _decoder->set(nID,"layer",part1layerId); - _decoder->set(nID,m_thetaID,minCellThetaBin); - cellNeighbours.push_back(nID); // add the first theta bin cell from the part1 layer - } - if(aDiagonal && Rmin == currentLayerRmax) + // find the layers in the part2 that share a border with the current layer + for(int part2layerId = minLayerIdEndcap[1]; part2layerId <= maxLayerIdEndcap[1]; part2layerId++) { - CellID nID = cID ; - _decoder->set(nID,"layer",part1layerId); - _decoder->set(nID,m_thetaID,minCellThetaBin); - cellNeighbours.push_back(nID); // add the first theta bin cell from the part1 layer + double Rmin = m_radii[part2layerId] - 0.5*m_layerDepth[part2layerId]; + double Rmax = m_radii[part2layerId] + 0.5*m_layerDepth[part2layerId]; + + if( (Rmin >= currentLayerRmin && Rmin <= currentLayerRmax) + || (Rmax > currentLayerRmin && Rmax <= currentLayerRmax) + || (currentLayerRmin >= Rmin && currentLayerRmin < Rmax) + || (currentLayerRmax >= Rmin && currentLayerRmax <= Rmax) + ) + { + CellID nID = cID ; + _decoder->set(nID,m_layerID,part2layerId); + _decoder->set(nID,m_thetaID,maxCellThetaBin); + cellNeighbours.push_back(nID); // add the last theta bin cell from part2 layer + } + if(aDiagonal && Rmax == currentLayerRmin) + { + CellID nID = cID ; + _decoder->set(nID,m_layerID,part2layerId); + _decoder->set(nID,m_thetaID,maxCellThetaBin); + cellNeighbours.push_back(nID); // add the last theta bin cell from part2 layer + } } } - } - // if it is the first theta bin cell in the part2 - if(EndcapPart2 && currentCellThetaBin == minCellThetaBin) - { - // find the layers in the part3 that share a border with the current layer - for(int part3layerId = minLayerIdEndcap[2]; part3layerId <= maxLayerIdEndcap[2]; part3layerId++) + // if the Endcap consists of more than 2 parts: + for(uint i_section = 1; i_section < (m_offsetZ.size()-1); i_section++) { - double Rmin = m_radii[part3layerId] - 0.5*m_layerDepth[part3layerId]; - double Rmax = m_radii[part3layerId] + 0.5*m_layerDepth[part3layerId]; - - if( (Rmin >= currentLayerRmin && Rmin <= currentLayerRmax) - || (Rmax > currentLayerRmin && Rmax <= currentLayerRmax) - || (currentLayerRmin >= Rmin && currentLayerRmin < Rmax) - || (currentLayerRmax >= Rmin && currentLayerRmax <= Rmax) - ) + if(i_section != EndcapPart) continue; + + // if it is the last theta bin cell then look for neighbours in previous part + if(currentCellThetaBin == maxCellThetaBin) { - CellID nID = cID ; - _decoder->set(nID,"layer",part3layerId); - _decoder->set(nID,m_thetaID,maxCellThetaBin); - cellNeighbours.push_back(nID); // add the first cell from the part3 layer + // find the layers in the previous part that share a border with the current layer + for(int prevPartLayerId = minLayerIdEndcap[i_section-1]; prevPartLayerId <= maxLayerIdEndcap[i_section-1]; prevPartLayerId++) + { + double Rmin = m_radii[prevPartLayerId] - 0.5*m_layerDepth[prevPartLayerId]; + double Rmax = m_radii[prevPartLayerId] + 0.5*m_layerDepth[prevPartLayerId]; + + if( (Rmin >= currentLayerRmin && Rmin < currentLayerRmax) + || (Rmax >= currentLayerRmin && Rmax <= currentLayerRmax) + || (currentLayerRmin >= Rmin && currentLayerRmin <= Rmax) + || (currentLayerRmax > Rmin && currentLayerRmax <= Rmax) + ) + { + CellID nID = cID ; + _decoder->set(nID,m_layerID,prevPartLayerId); + _decoder->set(nID,m_thetaID,minCellThetaBin); + cellNeighbours.push_back(nID); // add the first theta bin cell from the part1 layer + } + if(aDiagonal && Rmin == currentLayerRmax) + { + CellID nID = cID ; + _decoder->set(nID,m_layerID,prevPartLayerId); + _decoder->set(nID,m_thetaID,minCellThetaBin); + cellNeighbours.push_back(nID); // add the first theta bin cell from the part1 layer + } + } } - if(aDiagonal && Rmax == currentLayerRmin) + + // if it is the first theta bin cell then look for neighbours in the next part + if(currentCellThetaBin == minCellThetaBin) { - CellID nID = cID ; - _decoder->set(nID,"layer",part3layerId); - _decoder->set(nID,m_thetaID,maxCellThetaBin); - cellNeighbours.push_back(nID); // add the first cell from the part3 layer + // find the layers in the next part that share a border with the current layer + for(int nextPartLayerId = minLayerIdEndcap[i_section+1]; nextPartLayerId <= maxLayerIdEndcap[i_section+1]; nextPartLayerId++) + { + double Rmin = m_radii[nextPartLayerId] - 0.5*m_layerDepth[nextPartLayerId]; + double Rmax = m_radii[nextPartLayerId] + 0.5*m_layerDepth[nextPartLayerId]; + + if( (Rmin >= currentLayerRmin && Rmin <= currentLayerRmax) + || (Rmax > currentLayerRmin && Rmax <= currentLayerRmax) + || (currentLayerRmin >= Rmin && currentLayerRmin < Rmax) + || (currentLayerRmax >= Rmin && currentLayerRmax <= Rmax) + ) + { + CellID nID = cID ; + _decoder->set(nID,m_layerID,nextPartLayerId); + _decoder->set(nID,m_thetaID,maxCellThetaBin); + cellNeighbours.push_back(nID); // add the first cell from the part3 layer + } + if(aDiagonal && Rmax == currentLayerRmin) + { + CellID nID = cID ; + _decoder->set(nID,m_layerID,nextPartLayerId); + _decoder->set(nID,m_thetaID,maxCellThetaBin); + cellNeighbours.push_back(nID); // add the first cell from the part3 layer + } + } } } - } - // if it is the last theta bin cell in the part3 - if(EndcapPart3 && currentCellThetaBin == maxCellThetaBin) - { - // find the layers in the part2 that share a border with the current layer - for(int part2layerId = minLayerIdEndcap[1]; part2layerId <= maxLayerIdEndcap[1]; part2layerId++) + // if it is the last theta bin cell in the last part of the Endcap + if(EndcapPart == (m_offsetZ.size() - 1) && currentCellThetaBin == maxCellThetaBin) { - double Rmin = m_radii[part2layerId] - 0.5*m_layerDepth[part2layerId]; - double Rmax = m_radii[part2layerId] + 0.5*m_layerDepth[part2layerId]; - - if( (Rmin >= currentLayerRmin && Rmin < currentLayerRmax) - || (Rmax >= currentLayerRmin && Rmax <= currentLayerRmax) - || (currentLayerRmin >= Rmin && currentLayerRmin <= Rmax) - || (currentLayerRmax > Rmin && currentLayerRmax <= Rmax) - ) + // find the layers in the previous part that share a border with the current layer + for(int prevPartLayerId = minLayerIdEndcap[m_offsetZ.size()-2]; prevPartLayerId <= maxLayerIdEndcap[m_offsetZ.size()-2]; prevPartLayerId++) { - CellID nID = cID ; - _decoder->set(nID,"layer",part2layerId); - _decoder->set(nID,m_thetaID,minCellThetaBin); - cellNeighbours.push_back(nID); // add the first theta bin cell from the part2 layer - } - if(aDiagonal && Rmin == currentLayerRmax) - { - CellID nID = cID ; - _decoder->set(nID,"layer",part2layerId); - _decoder->set(nID,m_thetaID,minCellThetaBin); - cellNeighbours.push_back(nID); // add the first theta bin cell from the part2 layer + double Rmin = m_radii[prevPartLayerId] - 0.5*m_layerDepth[prevPartLayerId]; + double Rmax = m_radii[prevPartLayerId] + 0.5*m_layerDepth[prevPartLayerId]; + + if( (Rmin >= currentLayerRmin && Rmin < currentLayerRmax) + || (Rmax >= currentLayerRmin && Rmax <= currentLayerRmax) + || (currentLayerRmin >= Rmin && currentLayerRmin <= Rmax) + || (currentLayerRmax > Rmin && currentLayerRmax <= Rmax) + ) + { + CellID nID = cID ; + _decoder->set(nID,m_layerID,prevPartLayerId); + _decoder->set(nID,m_thetaID,minCellThetaBin); + cellNeighbours.push_back(nID); // add the first theta bin cell from the part2 layer + } + if(aDiagonal && Rmin == currentLayerRmax) + { + CellID nID = cID ; + _decoder->set(nID,m_layerID,prevPartLayerId); + _decoder->set(nID,m_thetaID,minCellThetaBin); + cellNeighbours.push_back(nID); // add the first theta bin cell from the part2 layer + } } } } @@ -816,7 +750,7 @@ std::array FCCSWHCalPhiTheta_k4geo::cellTheta(const CellID& cID) cons // get the cell index int idx = _decoder->get(cID, m_thetaID); // get the layer index - uint layer = _decoder->get(cID,"layer"); + uint layer = _decoder->get(cID,m_layerID); if(m_radii.empty()) defineCellsInRZplan(); if(m_cellEdges.empty()) return cTheta;