Back to home page

Project CMSSW displayed by LXR

 
 

    


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 //#define EDM_ML_DEBUG
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   // allocate geometry
0023   HGCalGeometry* geom = new HGCalGeometry(topology);
0024   unsigned int numberOfCells = topology.totalGeomModules();  // both sides
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   // loop over modules
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()) {  // front should always point to the center, so swap front and back
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()) {  // front should always point to the center, so swap front and back
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 }