EcalAlgo/src/EcalPreshowerGeometry.cc

0001 #include "Geometry/CaloGeometry/interface/PreshowerStrip.h"
0002 #include "Geometry/CaloGeometry/interface/CaloGenericDetId.h"
0003 #include "Geometry/EcalAlgo/interface/EcalPreshowerGeometry.h"
0004 #include "DataFormats/EcalDetId/interface/ESDetId.h"
0005 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0006 #include "FWCore/Utilities/interface/Exception.h"
0007 #include <iostream>
0008
0009 typedef CaloCellGeometry::CCGFloat CCGFloat;
0010 typedef CaloCellGeometry::Pt3D Pt3D;
0011 typedef CaloCellGeometry::Pt3DVec Pt3DVec;
0012 typedef HepGeom::Plane3D<CCGFloat> Pl3D;
0013
0014 //#define EDM_ML_DEBUG
0015
0016 EcalPreshowerGeometry::EcalPreshowerGeometry()
0017     : m_xWidWaf(6.3),
0018       m_xInterLadGap(0.05),  // additional gap between wafers in adj ladders
0019       m_xIntraLadGap(0.04),  // gap between wafers in same ladder
0020       m_yWidAct(6.1),
0021       m_yCtrOff(0.05),  // gap at center
0022       m_cellVec(k_NumberOfCellsForCorners) {
0023   m_zplane[0] = 0.;
0024   m_zplane[1] = 0.;
0025   m_zplane[2] = 0.;
0026   m_zplane[3] = 0.;
0027 #ifdef EDM_ML_DEBUG
0028   edm::LogVerbatim("EcalGeom") << "EcalPreshowerGeometry::Creating an instance";
0029 #endif
0030 }
0031
0032 EcalPreshowerGeometry::~EcalPreshowerGeometry() {}
0033
0034 unsigned int EcalPreshowerGeometry::alignmentTransformIndexLocal(const DetId& id) {
0035   const CaloGenericDetId gid(id);
0036
0037   assert(gid.isES());
0038
0039   // plane 2 is split into 2 dees along x=0 for both signs of z
0040
0041   // plane 1 at zsign=-1 is split into 2 dees between six=19 and six=20 for siy<=20,
0042   //                                             and six=21 and 22 for siy>=21
0043
0044   // plane 1 at zsign=+1 is split into 2 dees between six=20 and six=21 for siy<=20,
0045   //                                             and six=19 and 20 for siy>=21
0046
0047   // Desired numbering
0048   //                LEFT    RIGHT (as one faces the Dee from the IP)
0049   //  ES-  pl=2     0       1
0050   //       pl=1     2       3    the reversal of pl=2 and pl=1 is intentional here (CM Kuo)
0051   //  ES+  pl=1     4       5
0052   //       pl=2     6       7
0053
0054   const ESDetId esid(id);
0055   const int jx(esid.six() - 1);
0056   const int jy(esid.siy() - 1);
0057   const int jz(esid.zside() + 1);
0058   const int pl(esid.plane() - 1);
0059   const bool second(1 == pl);
0060   const bool top(19 < jy);
0061   const bool negz(0 == jz);
0062   const int lrl(19 > jx ? 0 : 1);
0063   const int lrr(21 > jx ? 0 : 1);
0064
0065   return (second ? jx / 20 + 3 * jz :        // 2nd plane split along middle
0066               (negz && !top ? lrl + 2 :      // 1st plane at neg z and bottom half split at six=19&20
0067                    (negz && top ? lrr + 2 :  // 1st plane at neg z and top half split at six=21&22
0068                         (!negz && !top ? lrr + 4 : lrl + 4))));  // opposite at positive z
0069 }
0070
0071 DetId EcalPreshowerGeometry::detIdFromLocalAlignmentIndex(unsigned int iLoc) {
0072   return ESDetId(1, 10 + 20 * (iLoc % 2), 10, 2 > iLoc || 5 < iLoc ? 2 : 1, 2 * (iLoc / 4) - 1);
0073 }
0074
0075 unsigned int EcalPreshowerGeometry::alignmentTransformIndexGlobal(const DetId& /*id*/) {
0076   return (unsigned int)DetId::Ecal - 1;
0077 }
0078
0079 void EcalPreshowerGeometry::initializeParms() {
0080   unsigned int n1minus(0);
0081   unsigned int n2minus(0);
0082   unsigned int n1plus(0);
0083   unsigned int n2plus(0);
0084   CCGFloat z1minus(0);
0085   CCGFloat z2minus(0);
0086   CCGFloat z1plus(0);
0087   CCGFloat z2plus(0);
0088   const std::vector<DetId>& esDetIds(getValidDetIds());
0089 #ifdef EDM_ML_DEBUG
0090   edm::LogVerbatim("EcalGeom") << "EcalPreshowerGeometry:: Get " << esDetIds.size() << " valid DetIds";
0091 #endif
0092
0093   for (unsigned int i(0); i != esDetIds.size(); ++i) {
0094     const ESDetId esid(esDetIds[i]);
0095     auto cell = getGeometry(esid);
0096     if (nullptr != cell) {
0097       const CCGFloat zz(cell->getPosition().z());
0098       if (1 == esid.plane()) {
0099         if (0 > esid.zside()) {
0100           z1minus += zz;
0101           ++n1minus;
0102         } else {
0103           z1plus += zz;
0104           ++n1plus;
0105         }
0106       }
0107       if (2 == esid.plane()) {
0108         if (0 > esid.zside()) {
0109           z2minus += zz;
0110           ++n2minus;
0111         } else {
0112           z2plus += zz;
0113           ++n2plus;
0114         }
0115       }
0116     }
0117   }
0118   assert(0 != n1minus && 0 != n2minus && 0 != n1plus && 0 != n2plus);
0119   z1minus /= (1. * n1minus);
0120   z2minus /= (1. * n2minus);
0121   z1plus /= (1. * n1plus);
0122   z2plus /= (1. * n2plus);
0123   assert(0 != z1minus && 0 != z2minus && 0 != z1plus && 0 != z2plus);
0124   setzPlanes(z1minus, z2minus, z1plus, z2plus);
0125 }
0126
0127 void EcalPreshowerGeometry::setzPlanes(CCGFloat z1minus, CCGFloat z2minus, CCGFloat z1plus, CCGFloat z2plus) {
0128   assert(0 > z1minus && 0 > z2minus && 0 < z1plus && 0 < z2plus);
0129
0130   m_zplane[0] = z1minus;
0131   m_zplane[1] = z2minus;
0132   m_zplane[2] = z1plus;
0133   m_zplane[3] = z2plus;
0134 }
0135
0136 // Get closest cell, etc...
0137 DetId EcalPreshowerGeometry::getClosestCell(const GlobalPoint& point) const { return getClosestCellInPlane(point, 2); }
0138
0139 DetId EcalPreshowerGeometry::getClosestCellInPlane(const GlobalPoint& point, int plane) const {
0140   const CCGFloat x(point.x());
0141   const CCGFloat y(point.y());
0142   const CCGFloat z(point.z());
0143
0144   if (0 == z || 1 > plane || 2 < plane)
0145     return DetId(0);
0146
0147   const unsigned int iz((0 > z ? 0 : 2) + plane - 1);
0148
0149   const CCGFloat ze(m_zplane[iz]);
0150   const CCGFloat xe(x * ze / z);
0151   const CCGFloat ye(y * ze / z);
0152
0153   const CCGFloat x0(1 == plane ? xe : ye);
0154   const CCGFloat y0(1 == plane ? ye : xe);
0155
0156   static const CCGFloat xWid(m_xWidWaf + m_xIntraLadGap + m_xInterLadGap);
0157
0158   const int row(1 + int(y0 + 20. * m_yWidAct - m_yCtrOff) / m_yWidAct);
0159   const int col(1 + int((x0 + 20. * xWid) / xWid));
0160
0161   CCGFloat closest(1e9);
0162
0163   DetId detId(0);
0164
0165   const int jz(0 > ze ? -1 : 1);
0166
0167 #ifdef EDM_ML_DEBUG
0168   edm::LogVerbatim("EcalGeom") << "** p=" << point << ", (" << xe << ", " << ye << ", " << ze << "), row=" << row
0169                                << ", col=" << col;
0170 #endif
0171   for (int ix(-1); ix != 2; ++ix)  // search within +-1 in row and col
0172   {
0173     for (int iy(-1); iy != 2; ++iy) {
0174       for (int jstrip(ESDetId::ISTRIP_MIN); jstrip <= ESDetId::ISTRIP_MAX; ++jstrip) {
0175         const int jx(1 == plane ? col + ix : row + iy);
0176         const int jy(1 == plane ? row + iy : col + ix);
0177         if (ESDetId::validDetId(jstrip, jx, jy, plane, jz)) {
0178           const ESDetId esId(jstrip, jx, jy, plane, jz);
0179           auto cell = getGeometry(esId);
0180           if (nullptr != cell) {
0181             const GlobalPoint& p(cell->getPosition());
0182             const CCGFloat dist2((p.x() - xe) * (p.x() - xe) + (p.y() - ye) * (p.y() - ye));
0183             if (dist2 < closest && present(esId)) {
0184               closest = dist2;
0185               detId = esId;
0186             }
0187           }
0188         }
0189       }
0190     }
0191   }
0192   return detId;
0193 }
0194
0195 void EcalPreshowerGeometry::localCorners(Pt3DVec& lc, const CCGFloat* pv, unsigned int /*i*/, Pt3D& ref) {
0196   PreshowerStrip::localCorners(lc, pv, ref);
0197 }
0198
0199 void EcalPreshowerGeometry::newCell(const GlobalPoint& f1,
0200                                     const GlobalPoint& /*f2*/,
0201                                     const GlobalPoint& /*f3*/,
0202                                     const CCGFloat* parm,
0203                                     const DetId& detId) {
0204   const unsigned int cellIndex(ESDetId(detId).denseIndex());
0205   m_cellVec[cellIndex] = PreshowerStrip(f1, cornersMgr(), parm);
0206   addValidID(detId);
0207 #ifdef EDM_ML_DEBUG
0208   edm::LogVerbatim("EcalGeom") << "EcalPreshowerGeometry::Add a new DetId " << ESDetId(detId);
0209 #endif
0210 }
0211
0212 CaloCellGeometryPtr EcalPreshowerGeometry::getGeometryRawPtr(uint32_t index) const {
0213   return CaloCellGeometryPtr(m_cellVec.size() <= index || nullptr == m_cellVec[index].param() ? nullptr
0214                                                                                               : &m_cellVec[index]);
0215 }