File indexing completed on 2025-04-26 02:01:04
0001 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0002 #include "Geometry/HGCalCommonData/interface/HGCalWaferIndex.h"
0003 #include "Geometry/HGCalGeometry/interface/HGCalGeometryLoader.h"
0004 #include "Geometry/HGCalGeometry/interface/HGCalGeometry.h"
0005 #include "Geometry/CaloGeometry/interface/CaloCellGeometry.h"
0006 #include "Geometry/CaloGeometry/interface/FlatTrd.h"
0007 #include "Geometry/HGCalCommonData/interface/HGCalDDDConstants.h"
0008 #include "DataFormats/ForwardDetId/interface/HFNoseDetId.h"
0009 #include "DataFormats/ForwardDetId/interface/HGCalDetId.h"
0010 #include "DataFormats/ForwardDetId/interface/HGCSiliconDetId.h"
0011 #include "DataFormats/ForwardDetId/interface/HGCScintillatorDetId.h"
0012 #include "DataFormats/ForwardDetId/interface/ForwardSubdetector.h"
0013
0014
0015
0016 typedef CaloCellGeometry::CCGFloat CCGFloat;
0017 typedef std::vector<float> ParmVec;
0018
0019 HGCalGeometryLoader::HGCalGeometryLoader() : twoBysqrt3_(2.0 / std::sqrt(3.0)) {}
0020
0021 HGCalGeometry* HGCalGeometryLoader::build(const HGCalTopology& topology) {
0022
0023 HGCalGeometry* geom = new HGCalGeometry(topology);
0024 unsigned int numberOfCells = topology.totalGeomModules();
0025 unsigned int numberExpected = topology.allGeomModules();
0026 parametersPerShape_ = (topology.tileTrapezoid() ? (int)HGCalGeometry::k_NumberOfParametersPerTrd
0027 : (int)HGCalGeometry::k_NumberOfParametersPerHex);
0028 uint32_t numberOfShapes =
0029 (topology.tileTrapezoid() ? HGCalGeometry::k_NumberOfShapesTrd : HGCalGeometry::k_NumberOfShapes);
0030 HGCalGeometryMode::GeometryMode mode = topology.geomMode();
0031 bool test = ((mode == HGCalGeometryMode::TrapezoidModule) || (mode == HGCalGeometryMode::TrapezoidCassette) ||
0032 (mode == HGCalGeometryMode::TrapezoidFineCell));
0033 #ifdef EDM_ML_DEBUG
0034 edm::LogVerbatim("HGCalGeom") << "Number of Cells " << numberOfCells << ":" << numberExpected << " for sub-detector "
0035 << topology.subDetector() << " Shapes " << numberOfShapes << ":" << parametersPerShape_
0036 << " mode " << mode;
0037 #endif
0038 geom->allocateCorners(numberOfCells);
0039 geom->allocatePar(numberOfShapes, parametersPerShape_);
0040
0041
0042 ParmVec params(parametersPerShape_, 0);
0043 unsigned int counter(0);
0044 #ifdef EDM_ML_DEBUG
0045 edm::LogVerbatim("HGCalGeom") << "HGCalGeometryLoader with # of "
0046 << "transformation matrices " << topology.dddConstants().getTrFormN() << " and "
0047 << topology.dddConstants().volumes() << ":" << topology.dddConstants().sectors()
0048 << " volumes";
0049 #endif
0050 for (unsigned itr = 0; itr < topology.dddConstants().getTrFormN(); ++itr) {
0051 HGCalParameters::hgtrform mytr = topology.dddConstants().getTrForm(itr);
0052 int zside = mytr.zp;
0053 int layer = mytr.lay;
0054 #ifdef EDM_ML_DEBUG
0055 unsigned int kount(0);
0056 edm::LogVerbatim("HGCalGeom") << "HGCalGeometryLoader:: Z:Layer " << zside << ":" << layer << " z " << mytr.h3v.z();
0057 #endif
0058 if (topology.waferHexagon6()) {
0059 ForwardSubdetector subdet = topology.subDetector();
0060 for (int wafer = 0; wafer < topology.dddConstants().sectors(); ++wafer) {
0061 std::string code[2] = {"False", "True"};
0062 if (topology.dddConstants().waferInLayer(wafer, layer, true)) {
0063 int type = topology.dddConstants().waferTypeT(wafer);
0064 if (type != 1)
0065 type = 0;
0066 DetId detId = static_cast<DetId>(HGCalDetId(subdet, zside, layer, type, wafer, 0));
0067 const auto& w = topology.dddConstants().waferPosition(wafer, true);
0068 double xx = (zside > 0) ? w.first : -w.first;
0069 CLHEP::Hep3Vector h3v(xx, w.second, mytr.h3v.z());
0070 const HepGeom::Transform3D ht3d(mytr.hr, h3v);
0071 #ifdef EDM_ML_DEBUG
0072 edm::LogVerbatim("HGCalGeom") << "HGCalGeometryLoader:: Wafer:Type " << wafer << ":" << type << " DetId "
0073 << HGCalDetId(detId) << std::hex << " " << detId.rawId() << std::dec
0074 << " transf " << ht3d.getTranslation() << " and " << ht3d.getRotation();
0075 #endif
0076 HGCalParameters::hgtrap vol = topology.dddConstants().getModule(wafer, true, true);
0077 params[FlatHexagon::k_dZ] = vol.dz;
0078 params[FlatHexagon::k_r] = topology.dddConstants().cellSizeHex(type);
0079 params[FlatHexagon::k_R] = twoBysqrt3_ * params[FlatHexagon::k_r];
0080
0081 buildGeom(params, ht3d, detId, geom, 0);
0082 counter++;
0083 #ifdef EDM_ML_DEBUG
0084 ++kount;
0085 #endif
0086 }
0087 }
0088 } else if (topology.tileTrapezoid()) {
0089 int indx = topology.dddConstants().layerIndex(layer, true);
0090 int ring = topology.dddConstants().getParameter()->iradMinBH_[indx];
0091 int nphi = topology.dddConstants().getParameter()->scintCells(layer);
0092 int type = topology.dddConstants().getParameter()->scintType(layer);
0093 #ifdef EDM_ML_DEBUG
0094 edm::LogVerbatim("HGCalGeom") << "Layer " << layer << ":" << indx << " Ring " << ring << ":"
0095 << topology.dddConstants().getParameter()->iradMaxBH_[indx] << " Phi " << nphi
0096 << " Type " << type;
0097 #endif
0098 for (int md = topology.dddConstants().getParameter()->firstModule_[indx];
0099 md <= topology.dddConstants().getParameter()->lastModule_[indx];
0100 ++md) {
0101 for (int iphi = 1; iphi <= nphi; ++iphi) {
0102 HGCScintillatorDetId id(type, layer, zside * ring, iphi);
0103 std::pair<int, int> typm = topology.dddConstants().tileType(layer, ring, 0);
0104 if (typm.first >= 0) {
0105 id.setType(typm.first);
0106 id.setSiPM(typm.second);
0107 }
0108 bool ok = test ? topology.dddConstants().tileExist(zside, layer, ring, iphi) : true;
0109 #ifdef EDM_ML_DEBUG
0110 edm::LogVerbatim("HGCalGeom") << "HGCalGeometryLoader::layer:rad:phi:type:sipm " << layer << ":"
0111 << ring * zside << ":" << iphi << ":" << type << ":" << typm.first << ":"
0112 << typm.second << " Test " << test << ":" << ok << " ID " << id;
0113 #endif
0114 if (ok) {
0115 DetId detId = static_cast<DetId>(id);
0116 const auto& w = topology.dddConstants().locateCellTrap(zside, layer, ring, iphi, true, false);
0117 double xx = (zside > 0) ? w.first : -w.first;
0118 CLHEP::Hep3Vector h3v(xx, w.second, mytr.h3v.z());
0119 const HepGeom::Transform3D ht3d(mytr.hr, h3v);
0120 #ifdef EDM_ML_DEBUG
0121 edm::LogVerbatim("HGCalGeom")
0122 << "HGCalGeometryLoader::rad:phi:type " << ring * zside << ":" << iphi << ":" << type << " DetId "
0123 << HGCScintillatorDetId(detId) << " " << std::hex << detId.rawId() << std::dec << " transf "
0124 << ht3d.getTranslation() << " R " << ht3d.getTranslation().perp() << " and " << ht3d.getRotation();
0125 #endif
0126 HGCalParameters::hgtrap vol = topology.dddConstants().getModule(md, false, true);
0127 params[FlatTrd::k_dZ] = vol.dz;
0128 params[FlatTrd::k_Theta] = params[FlatTrd::k_Phi] = 0;
0129 params[FlatTrd::k_dY1] = params[FlatTrd::k_dY2] = vol.h;
0130 params[FlatTrd::k_dX1] = params[FlatTrd::k_dX3] = vol.bl;
0131 params[FlatTrd::k_dX2] = params[FlatTrd::k_dX4] = vol.tl;
0132 params[FlatTrd::k_Alp1] = params[FlatTrd::k_Alp2] = 0;
0133 params[FlatTrd::k_Cell] = topology.dddConstants().cellSizeHex(type);
0134
0135 buildGeom(params, ht3d, detId, geom, 1);
0136 counter++;
0137 #ifdef EDM_ML_DEBUG
0138 ++kount;
0139 #endif
0140 }
0141 }
0142 ++ring;
0143 }
0144 } else {
0145 DetId::Detector det = topology.detector();
0146 for (int wafer = 0; wafer < topology.dddConstants().sectors(); ++wafer) {
0147 if (topology.dddConstants().waferInLayer(wafer, layer, true)) {
0148 int copy = topology.dddConstants().getParameter()->waferCopy_[wafer];
0149 int u = HGCalWaferIndex::waferU(copy);
0150 int v = HGCalWaferIndex::waferV(copy);
0151 int type = topology.dddConstants().getTypeHex(layer, u, v);
0152 DetId detId =
0153 (topology.isHFNose() ? static_cast<DetId>(HFNoseDetId(zside, type, layer, u, v, 0, 0))
0154 : static_cast<DetId>(HGCSiliconDetId(det, zside, type, layer, u, v, 0, 0)));
0155 const auto& w = topology.dddConstants().waferPosition(layer, u, v, true, true);
0156 double xx = (zside > 0) ? w.first : -w.first;
0157 CLHEP::Hep3Vector h3v(xx, w.second, mytr.h3v.z());
0158 const HepGeom::Transform3D ht3d(mytr.hr, h3v);
0159 #ifdef EDM_ML_DEBUG
0160 if (topology.isHFNose())
0161 edm::LogVerbatim("HGCalGeom") << "HGCalGeometryLoader::Wafer:Type " << wafer << ":" << type << " DetId "
0162 << HFNoseDetId(detId) << std::hex << " " << detId.rawId() << std::dec
0163 << " trans " << ht3d.getTranslation() << " and " << ht3d.getRotation();
0164 else
0165 edm::LogVerbatim("HGCalGeom") << "HGCalGeometryLoader::Wafer:Type " << wafer << ":" << type << " DetId "
0166 << HGCSiliconDetId(detId) << std::hex << " " << detId.rawId() << std::dec
0167 << " trans " << ht3d.getTranslation() << " and " << ht3d.getRotation();
0168 #endif
0169 HGCalParameters::hgtrap vol = topology.dddConstants().getModule(type, false, true);
0170 params[FlatHexagon::k_dZ] = vol.dz;
0171 params[FlatHexagon::k_r] = topology.dddConstants().cellSizeHex(type);
0172 params[FlatHexagon::k_R] = twoBysqrt3_ * params[FlatHexagon::k_r];
0173
0174 buildGeom(params, ht3d, detId, geom, 0);
0175 counter++;
0176 #ifdef EDM_ML_DEBUG
0177 ++kount;
0178 #endif
0179 }
0180 }
0181 }
0182 #ifdef EDM_ML_DEBUG
0183 edm::LogVerbatim("HGCalGeom") << kount << " modules found in Layer " << layer << " Z " << zside;
0184 #endif
0185 }
0186
0187 geom->sortDetIds();
0188
0189 if (counter != numberExpected) {
0190 if (topology.tileTrapezoid()) {
0191 edm::LogVerbatim("HGCalGeom") << "Inconsistent # of cells: expected " << numberExpected << ":" << numberOfCells
0192 << " , inited " << counter;
0193 } else {
0194 edm::LogError("HGCalGeom") << "Inconsistent # of cells: expected " << numberExpected << ":" << numberOfCells
0195 << " , inited " << counter;
0196 assert(counter == numberExpected);
0197 }
0198 }
0199
0200 return geom;
0201 }
0202
0203 void HGCalGeometryLoader::buildGeom(
0204 const ParmVec& params, const HepGeom::Transform3D& ht3d, const DetId& detId, HGCalGeometry* geom, int mode) {
0205 #ifdef EDM_ML_DEBUG
0206 for (int i = 0; i < parametersPerShape_; ++i)
0207 edm::LogVerbatim("HGCalGeom") << "Parameter[" << i << "] : " << params[i];
0208 #endif
0209 if (mode == 1) {
0210 std::vector<GlobalPoint> corners(FlatTrd::ncorner_);
0211
0212 FlatTrd::createCorners(params, ht3d, corners);
0213
0214 const CCGFloat* parmPtr(CaloCellGeometry::getParmPtr(params, geom->parMgr(), geom->parVecVec()));
0215
0216 GlobalPoint front(0.25 * (corners[0].x() + corners[1].x() + corners[2].x() + corners[3].x()),
0217 0.25 * (corners[0].y() + corners[1].y() + corners[2].y() + corners[3].y()),
0218 0.25 * (corners[0].z() + corners[1].z() + corners[2].z() + corners[3].z()));
0219
0220 GlobalPoint back(0.25 * (corners[4].x() + corners[5].x() + corners[6].x() + corners[7].x()),
0221 0.25 * (corners[4].y() + corners[5].y() + corners[6].y() + corners[7].y()),
0222 0.25 * (corners[4].z() + corners[5].z() + corners[6].z() + corners[7].z()));
0223
0224 if (front.mag2() > back.mag2()) {
0225 std::swap(front, back);
0226 std::swap_ranges(corners.begin(), corners.begin() + FlatTrd::ncornerBy2_, corners.begin() + FlatTrd::ncornerBy2_);
0227 }
0228 geom->newCell(front, back, corners[0], parmPtr, detId);
0229 } else {
0230 std::vector<GlobalPoint> corners(FlatHexagon::ncorner_);
0231
0232 FlatHexagon::createCorners(params, ht3d, corners);
0233
0234 const CCGFloat* parmPtr(CaloCellGeometry::getParmPtr(params, geom->parMgr(), geom->parVecVec()));
0235
0236 GlobalPoint front(
0237 FlatHexagon::oneBySix_ *
0238 (corners[0].x() + corners[1].x() + corners[2].x() + corners[3].x() + corners[4].x() + corners[5].x()),
0239 FlatHexagon::oneBySix_ *
0240 (corners[0].y() + corners[1].y() + corners[2].y() + corners[3].y() + corners[4].y() + corners[5].y()),
0241 FlatHexagon::oneBySix_ *
0242 (corners[0].z() + corners[1].z() + corners[2].z() + corners[3].z() + corners[4].z() + corners[5].z()));
0243
0244 GlobalPoint back(
0245 FlatHexagon::oneBySix_ *
0246 (corners[6].x() + corners[7].x() + corners[8].x() + corners[9].x() + corners[10].x() + corners[11].x()),
0247 FlatHexagon::oneBySix_ *
0248 (corners[6].y() + corners[7].y() + corners[8].y() + corners[9].y() + corners[10].y() + corners[11].y()),
0249 FlatHexagon::oneBySix_ *
0250 (corners[6].z() + corners[7].z() + corners[8].z() + corners[9].z() + corners[10].z() + corners[11].z()));
0251
0252 if (front.mag2() > back.mag2()) {
0253 std::swap(front, back);
0254 std::swap_ranges(
0255 corners.begin(), corners.begin() + FlatHexagon::ncornerBy2_, corners.begin() + FlatHexagon::ncornerBy2_);
0256 }
0257 geom->newCell(front, back, corners[0], parmPtr, detId);
0258 }
0259 }