HGCalGeometry/src/HGCalTBGeometryLoader.cc

0001 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0002 #include "Geometry/HGCalGeometry/interface/HGCalTBGeometryLoader.h"
0003 #include "Geometry/HGCalGeometry/interface/HGCalTBGeometry.h"
0004 #include "Geometry/CaloGeometry/interface/CaloCellGeometry.h"
0005 #include "Geometry/CaloGeometry/interface/FlatTrd.h"
0006 #include "Geometry/HGCalTBCommonData/interface/HGCalTBDDDConstants.h"
0007 #include "DataFormats/ForwardDetId/interface/HGCalDetId.h"
0008 #include "DataFormats/ForwardDetId/interface/ForwardSubdetector.h"
0009
0010 //#define EDM_ML_DEBUG
0011
0012 typedef CaloCellGeometry::CCGFloat CCGFloat;
0013 typedef std::vector<float> ParmVec;
0014
0015 HGCalTBGeometryLoader::HGCalTBGeometryLoader() : twoBysqrt3_(2.0 / std::sqrt(3.0)) {}
0016
0017 HGCalTBGeometry* HGCalTBGeometryLoader::build(const HGCalTBTopology& topology) {
0018   // allocate geometry
0019   HGCalTBGeometry* geom = new HGCalTBGeometry(topology);
0020   unsigned int numberOfCells = topology.totalGeomModules();  // both sides
0021   unsigned int numberExpected = topology.allGeomModules();
0022   parametersPerShape_ = static_cast<int>(HGCalTBGeometry::k_NumberOfParametersPerHex);
0023   uint32_t numberOfShapes = HGCalTBGeometry::k_NumberOfShapes;
0024 #ifdef EDM_ML_DEBUG
0025   edm::LogVerbatim("HGCalGeom") << "Number of Cells " << numberOfCells << ":" << numberExpected << " for sub-detector "
0026                                 << topology.subDetector() << " Shapes " << numberOfShapes << ":" << parametersPerShape_;
0027 #endif
0028   geom->allocateCorners(numberOfCells);
0029   geom->allocatePar(numberOfShapes, parametersPerShape_);
0030
0031   // loop over modules
0032   ParmVec params(parametersPerShape_, 0);
0033   unsigned int counter(0);
0034 #ifdef EDM_ML_DEBUG
0035   edm::LogVerbatim("HGCalGeom") << "HGCalTBGeometryLoader with # of "
0036                                 << "transformation matrices " << topology.dddConstants().getTrFormN() << " and "
0037                                 << topology.dddConstants().volumes() << ":" << topology.dddConstants().sectors()
0038                                 << " volumes";
0039 #endif
0040   for (unsigned itr = 0; itr < topology.dddConstants().getTrFormN(); ++itr) {
0041     HGCalTBParameters::hgtrform mytr = topology.dddConstants().getTrForm(itr);
0042     int zside = mytr.zp;
0043     int layer = mytr.lay;
0044 #ifdef EDM_ML_DEBUG
0045     unsigned int kount(0);
0046     edm::LogVerbatim("HGCalGeom") << "HGCalTBGeometryLoader:: Z:Layer " << zside << ":" << layer << " z "
0047                                   << mytr.h3v.z();
0048 #endif
0049     ForwardSubdetector subdet = topology.subDetector();
0050     for (int wafer = 0; wafer < topology.dddConstants().sectors(); ++wafer) {
0051       std::string code[2] = {"False", "True"};
0052       if (topology.dddConstants().waferInLayer(wafer, layer, true)) {
0053         int type = topology.dddConstants().waferTypeT(wafer);
0054         if (type != 1)
0055           type = 0;
0056         DetId detId = static_cast<DetId>(HGCalDetId(subdet, zside, layer, type, wafer, 0));
0057         const auto& w = topology.dddConstants().waferPosition(wafer, true);
0058         double xx = (zside > 0) ? w.first : -w.first;
0059         CLHEP::Hep3Vector h3v(xx, w.second, mytr.h3v.z());
0060         const HepGeom::Transform3D ht3d(mytr.hr, h3v);
0061 #ifdef EDM_ML_DEBUG
0062         edm::LogVerbatim("HGCalGeom") << "HGCalTBGeometryLoader:: Wafer:Type " << wafer << ":" << type << " DetId "
0063                                       << HGCalDetId(detId) << std::hex << " " << detId.rawId() << std::dec << " transf "
0064                                       << ht3d.getTranslation() << " and " << ht3d.getRotation();
0065 #endif
0066         HGCalTBParameters::hgtrap vol = topology.dddConstants().getModule(wafer, true, true);
0067         params[FlatHexagon::k_dZ] = vol.dz;
0068         params[FlatHexagon::k_r] = topology.dddConstants().cellSizeHex(type);
0069         params[FlatHexagon::k_R] = twoBysqrt3_ * params[FlatHexagon::k_r];
0070
0071         buildGeom(params, ht3d, detId, geom);
0072         counter++;
0073 #ifdef EDM_ML_DEBUG
0074         ++kount;
0075 #endif
0076       }
0077     }
0078 #ifdef EDM_ML_DEBUG
0079     edm::LogVerbatim("HGCalGeom") << kount << " modules found in Layer " << layer << " Z " << zside;
0080 #endif
0081   }
0082
0083   geom->sortDetIds();
0084
0085   if (counter != numberExpected) {
0086     edm::LogError("HGCalGeom") << "Inconsistent # of cells: expected " << numberExpected << ":" << numberOfCells
0087                                << " , inited " << counter;
0088   }
0089
0090   return geom;
0091 }
0092
0093 void HGCalTBGeometryLoader::buildGeom(const ParmVec& params,
0094                                       const HepGeom::Transform3D& ht3d,
0095                                       const DetId& detId,
0096                                       HGCalTBGeometry* geom) {
0097 #ifdef EDM_ML_DEBUG
0098   for (int i = 0; i < parametersPerShape_; ++i)
0099     edm::LogVerbatim("HGCalGeom") << "Parameter[" << i << "] : " << params[i];
0100 #endif
0101   std::vector<GlobalPoint> corners(FlatHexagon::ncorner_);
0102
0103   FlatHexagon::createCorners(params, ht3d, corners);
0104
0105   const CCGFloat* parmPtr(CaloCellGeometry::getParmPtr(params, geom->parMgr(), geom->parVecVec()));
0106
0107   GlobalPoint front(
0108       FlatHexagon::oneBySix_ *
0109           (corners[0].x() + corners[1].x() + corners[2].x() + corners[3].x() + corners[4].x() + corners[5].x()),
0110       FlatHexagon::oneBySix_ *
0111           (corners[0].y() + corners[1].y() + corners[2].y() + corners[3].y() + corners[4].y() + corners[5].y()),
0112       FlatHexagon::oneBySix_ *
0113           (corners[0].z() + corners[1].z() + corners[2].z() + corners[3].z() + corners[4].z() + corners[5].z()));
0114
0115   GlobalPoint back(
0116       FlatHexagon::oneBySix_ *
0117           (corners[6].x() + corners[7].x() + corners[8].x() + corners[9].x() + corners[10].x() + corners[11].x()),
0118       FlatHexagon::oneBySix_ *
0119           (corners[6].y() + corners[7].y() + corners[8].y() + corners[9].y() + corners[10].y() + corners[11].y()),
0120       FlatHexagon::oneBySix_ *
0121           (corners[6].z() + corners[7].z() + corners[8].z() + corners[9].z() + corners[10].z() + corners[11].z()));
0122
0123   if (front.mag2() > back.mag2()) {  // front should always point to the center, so swap front and back
0124     std::swap(front, back);
0125     std::swap_ranges(
0126         corners.begin(), corners.begin() + FlatHexagon::ncornerBy2_, corners.begin() + FlatHexagon::ncornerBy2_);
0127   }
0128   geom->newCell(front, back, corners[0], parmPtr, detId);
0129 }