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