Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2022-06-20 00:21:44

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 #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   // loop over modules
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(layer, ring, iphi, true);
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()) {  // front should always point to the center, so swap front and back
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()) {  // front should always point to the center, so swap front and back
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 }