Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-10-25 23:57:06 

0001 #include "Geometry/CaloGeometry/interface/CaloGenericDetId.h"
0002 #include "Geometry/EcalAlgo/interface/EcalEndcapGeometry.h"
0003 #include "Geometry/CaloGeometry/interface/CaloCellGeometry.h"
0004 #include "Geometry/CaloGeometry/interface/TruncatedPyramid.h"
0005 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0006 #include "FWCore/Utilities/interface/Exception.h"
0007 #include <CLHEP/Geometry/Point3D.h>
0008 #include <CLHEP/Geometry/Plane3D.h>
0009 
0010 typedef CaloCellGeometry::CCGFloat CCGFloat;
0011 typedef CaloCellGeometry::Pt3D Pt3D;
0012 typedef CaloCellGeometry::Pt3DVec Pt3DVec;
0013 typedef HepGeom::Plane3D<CCGFloat> Pl3D;
0014 
0015 //#define EDM_ML_DEBUG
0016 
0017 EcalEndcapGeometry::EcalEndcapGeometry(void)
0018     : _nnmods(316),
0019       _nncrys(25),
0020       zeP(0.),
0021       zeN(0.),
0022       m_wref(0.),
0023       m_del(0.),
0024       m_nref(0),
0025       m_borderMgr(nullptr),
0026       m_borderPtrVec(nullptr),
0027       m_avgZ(-1),
0028       m_check(false),
0029       m_cellVec(k_NumberOfCellsForCorners) {
0030   m_xlo[0] = 999.;
0031   m_xlo[1] = 999.;
0032   m_xhi[0] = -999.;
0033   m_xhi[1] = -999.;
0034   m_ylo[0] = 999.;
0035   m_ylo[1] = 999.;
0036   m_yhi[0] = -999.;
0037   m_yhi[1] = -999.;
0038   m_xoff[0] = 0.;
0039   m_xoff[1] = 0.;
0040   m_yoff[0] = 0.;
0041   m_yoff[1] = 0.;
0042 }
0043 
0044 EcalEndcapGeometry::~EcalEndcapGeometry() {
0045   if (m_borderPtrVec) {
0046     auto ptr = m_borderPtrVec.load(std::memory_order_acquire);
0047     for (auto& v : (*ptr)) {
0048       delete v;
0049       v = nullptr;
0050     }
0051     delete m_borderPtrVec.load();
0052   }
0053   delete m_borderMgr.load();
0054 }
0055 
0056 unsigned int EcalEndcapGeometry::alignmentTransformIndexLocal(const DetId& id) {
0057   const CaloGenericDetId gid(id);
0058 
0059   assert(gid.isEE());
0060   unsigned int index(EEDetId(id).ix() / 51 + (EEDetId(id).zside() < 0 ? 0 : 2));
0061 
0062   return index;
0063 }
0064 
0065 DetId EcalEndcapGeometry::detIdFromLocalAlignmentIndex(unsigned int iLoc) {
0066   return EEDetId(20 + 50 * (iLoc % 2), 50, 2 * (iLoc / 2) - 1);
0067 }
0068 
0069 unsigned int EcalEndcapGeometry::alignmentTransformIndexGlobal(const DetId& /*id*/) {
0070   return (unsigned int)DetId::Ecal - 1;
0071 }
0072 
0073 void EcalEndcapGeometry::initializeParms() {
0074   zeP = 0.;
0075   zeN = 0.;
0076   unsigned nP = 0;
0077   unsigned nN = 0;
0078   m_nref = 0;
0079 
0080   for (uint32_t i(0); i != m_cellVec.size(); ++i) {
0081     auto cell = cellGeomPtr(i);
0082     if (nullptr != cell) {
0083       const CCGFloat z(cell->getPosition().z());
0084       if (z > 0.) {
0085         zeP += z;
0086         ++nP;
0087       } else {
0088         zeN += z;
0089         ++nN;
0090       }
0091       const EEDetId myId(EEDetId::detIdFromDenseIndex(i));
0092       const unsigned int ix(myId.ix());
0093       const unsigned int iy(myId.iy());
0094       if (ix > m_nref)
0095         m_nref = ix;
0096       if (iy > m_nref)
0097         m_nref = iy;
0098     }
0099   }
0100   if (0 < nP)
0101     zeP /= (CCGFloat)nP;
0102   if (0 < nN)
0103     zeN /= (CCGFloat)nN;
0104 
0105   m_xlo[0] = 999;
0106   m_xhi[0] = -999;
0107   m_ylo[0] = 999;
0108   m_yhi[0] = -999;
0109   m_xlo[1] = 999;
0110   m_xhi[1] = -999;
0111   m_ylo[1] = 999;
0112   m_yhi[1] = -999;
0113   for (uint32_t i(0); i != m_cellVec.size(); ++i) {
0114     auto cell = cellGeomPtr(i);
0115     if (nullptr != cell) {
0116       const GlobalPoint& p(cell->getPosition());
0117       const CCGFloat z(p.z());
0118       const CCGFloat zz(0 > z ? zeN : zeP);
0119       const CCGFloat x(p.x() * zz / z);
0120       const CCGFloat y(p.y() * zz / z);
0121 
0122       if (0 > z && x < m_xlo[0])
0123         m_xlo[0] = x;
0124       if (0 < z && x < m_xlo[1])
0125         m_xlo[1] = x;
0126       if (0 > z && y < m_ylo[0])
0127         m_ylo[0] = y;
0128       if (0 < z && y < m_ylo[1])
0129         m_ylo[1] = y;
0130 
0131       if (0 > z && x > m_xhi[0])
0132         m_xhi[0] = x;
0133       if (0 < z && x > m_xhi[1])
0134         m_xhi[1] = x;
0135       if (0 > z && y > m_yhi[0])
0136         m_yhi[0] = y;
0137       if (0 < z && y > m_yhi[1])
0138         m_yhi[1] = y;
0139     }
0140   }
0141 
0142   m_xoff[0] = (m_xhi[0] + m_xlo[0]) / 2.;
0143   m_xoff[1] = (m_xhi[1] + m_xlo[1]) / 2.;
0144   m_yoff[0] = (m_yhi[0] + m_ylo[0]) / 2.;
0145   m_yoff[1] = (m_yhi[1] + m_ylo[1]) / 2.;
0146 
0147   m_del = (m_xhi[0] - m_xlo[0] + m_xhi[1] - m_xlo[1] + m_yhi[0] - m_ylo[0] + m_yhi[1] - m_ylo[1]);
0148 
0149   if (1 != m_nref)
0150     m_wref = m_del / (4. * (m_nref - 1));
0151 
0152   m_xlo[0] -= m_wref / 2;
0153   m_xlo[1] -= m_wref / 2;
0154   m_xhi[0] += m_wref / 2;
0155   m_xhi[1] += m_wref / 2;
0156 
0157   m_ylo[0] -= m_wref / 2;
0158   m_ylo[1] -= m_wref / 2;
0159   m_yhi[0] += m_wref / 2;
0160   m_yhi[1] += m_wref / 2;
0161 
0162   m_del += m_wref;
0163 
0164 #ifdef EDM_ML_DEBUG
0165   edm::LogVerbatim("EcalGeom") << "zeP=" << zeP << ", zeN=" << zeN << ", nP=" << nP << ", nN=" << nN;
0166 
0167   edm::LogVerbatim("EcalGeom") << "xlo[0]=" << m_xlo[0] << ", xlo[1]=" << m_xlo[1] << ", xhi[0]=" << m_xhi[0]
0168                                << ", xhi[1]=" << m_xhi[1] << "\nylo[0]=" << m_ylo[0] << ", ylo[1]=" << m_ylo[1]
0169                                << ", yhi[0]=" << m_yhi[0] << ", yhi[1]=" << m_yhi[1];
0170 
0171   edm::LogVerbatim("EcalGeom") << "xoff[0]=" << m_xoff[0] << ", xoff[1]" << m_xoff[1] << ", yoff[0]=" << m_yoff[0]
0172                                << ", yoff[1]" << m_yoff[1];
0173 
0174   edm::LogVerbatim("EcalGeom") << "nref=" << m_nref << ", m_wref=" << m_wref;
0175 #endif
0176 }
0177 
0178 unsigned int EcalEndcapGeometry::xindex(CCGFloat x, CCGFloat z) const {
0179   const CCGFloat xlo(0 > z ? m_xlo[0] : m_xlo[1]);
0180   const int i(1 + int((x - xlo) / m_wref));
0181 
0182   return (1 > i ? 1 : (m_nref < (unsigned int)i ? m_nref : (unsigned int)i));
0183 }
0184 
0185 unsigned int EcalEndcapGeometry::yindex(CCGFloat y, CCGFloat z) const {
0186   const CCGFloat ylo(0 > z ? m_ylo[0] : m_ylo[1]);
0187   const int i(1 + int((y - ylo) / m_wref));
0188 
0189   return (1 > i ? 1 : (m_nref < (unsigned int)i ? m_nref : (unsigned int)i));
0190 }
0191 
0192 EEDetId EcalEndcapGeometry::gId(float x, float y, float z) const {
0193   const CCGFloat fac(fabs((0 > z ? zeN : zeP) / z));
0194   const unsigned int ix(xindex(x * fac, z));
0195   const unsigned int iy(yindex(y * fac, z));
0196   const unsigned int iz(z > 0 ? 1 : -1);
0197 
0198   if (EEDetId::validDetId(ix, iy, iz)) {
0199     return EEDetId(ix, iy, iz);  // first try is on target
0200   } else                         // try nearby coordinates, spiraling out from center
0201   {
0202     for (unsigned int i(1); i != 6; ++i) {
0203       for (unsigned int k(0); k != 8; ++k) {
0204         const int jx(0 == k || 4 == k || 5 == k ? +i : (1 == k || 5 < k ? -i : 0));
0205         const int jy(2 == k || 4 == k || 6 == k ? +i : (3 == k || 5 == k || 7 == k ? -i : 0));
0206         if (EEDetId::validDetId(ix + jx, iy + jy, iz)) {
0207           return EEDetId(ix + jx, iy + jy, iz);
0208         }
0209       }
0210     }
0211   }
0212   return EEDetId();  // nowhere near any crystal
0213 }
0214 
0215 // Get closest cell, etc...
0216 DetId EcalEndcapGeometry::getClosestCell(const GlobalPoint& r) const {
0217   try {
0218     EEDetId mycellID(gId(r.x(), r.y(), r.z()));  // educated guess
0219 
0220     if (EEDetId::validDetId(mycellID.ix(), mycellID.iy(), mycellID.zside())) {
0221       // now get points in convenient ordering
0222 
0223       Pt3D A;
0224       Pt3D B;
0225       Pt3D C;
0226       Pt3D point(r.x(), r.y(), r.z());
0227       // D.K. : equation of plane : AA*x+BB*y+CC*z+DD=0;
0228       // finding equation for each edge
0229 
0230       // ================================================================
0231       CCGFloat x, y, z;
0232       unsigned offset = 0;
0233       int zsign = 1;
0234       //================================================================
0235 
0236       // compute the distance of the point with respect of the 4 crystal lateral planes
0237 
0238       if (nullptr != getGeometry(mycellID)) {
0239         const GlobalPoint& myPosition = getGeometry(mycellID)->getPosition();
0240 
0241         x = myPosition.x();
0242         y = myPosition.y();
0243         z = myPosition.z();
0244 
0245         offset = 0;
0246         // This will disappear when Andre has applied his fix
0247 
0248         if (z > 0) {
0249           if (x > 0 && y > 0)
0250             offset = 1;
0251           else if (x < 0 && y > 0)
0252             offset = 2;
0253           else if (x > 0 && y < 0)
0254             offset = 0;
0255           else if (x < 0 && y < 0)
0256             offset = 3;
0257           zsign = 1;
0258         } else {
0259           if (x > 0 && y > 0)
0260             offset = 3;
0261           else if (x < 0 && y > 0)
0262             offset = 2;
0263           else if (x > 0 && y < 0)
0264             offset = 0;
0265           else if (x < 0 && y < 0)
0266             offset = 1;
0267           zsign = -1;
0268         }
0269         GlobalPoint corners[8];
0270         for (unsigned ic = 0; ic < 4; ++ic) {
0271           corners[ic] = getGeometry(mycellID)->getCorners()[(unsigned)((zsign * ic + offset) % 4)];
0272           corners[4 + ic] = getGeometry(mycellID)->getCorners()[(unsigned)(4 + (zsign * ic + offset) % 4)];
0273         }
0274 
0275         CCGFloat SS[4];
0276         for (short i = 0; i < 4; ++i) {
0277           A = Pt3D(corners[i % 4].x(), corners[i % 4].y(), corners[i % 4].z());
0278           B = Pt3D(corners[(i + 1) % 4].x(), corners[(i + 1) % 4].y(), corners[(i + 1) % 4].z());
0279           C = Pt3D(corners[4 + (i + 1) % 4].x(), corners[4 + (i + 1) % 4].y(), corners[4 + (i + 1) % 4].z());
0280           Pl3D plane(A, B, C);
0281           plane.normalize();
0282           CCGFloat distance = plane.distance(point);
0283           if (corners[0].z() < 0.)
0284             distance = -distance;
0285           SS[i] = distance;
0286         }
0287 
0288         // Only one move in necessary direction
0289 
0290         const bool yout(0 > SS[0] * SS[2]);
0291         const bool xout(0 > SS[1] * SS[3]);
0292 
0293         if (yout || xout) {
0294           const int ydel(!yout ? 0 : (0 < SS[0] ? -1 : 1));
0295           const int xdel(!xout ? 0 : (0 < SS[1] ? -1 : 1));
0296           const unsigned int ix(mycellID.ix() + xdel);
0297           const unsigned int iy(mycellID.iy() + ydel);
0298           const unsigned int iz(mycellID.zside());
0299           if (EEDetId::validDetId(ix, iy, iz))
0300             mycellID = EEDetId(ix, iy, iz);
0301         }
0302 
0303         return mycellID;
0304       }
0305       return DetId(0);
0306     }
0307   } catch (cms::Exception& e) {
0308     return DetId(0);
0309   }
0310   return DetId(0);
0311 }
0312 
0313 CaloSubdetectorGeometry::DetIdSet EcalEndcapGeometry::getCells(const GlobalPoint& r, double dR) const {
0314   CaloSubdetectorGeometry::DetIdSet dis;  // return object
0315   if (0.000001 < dR) {
0316     if (dR > M_PI / 2.)  // this code assumes small dR
0317     {
0318       dis = CaloSubdetectorGeometry::getCells(r, dR);  // base class version
0319     } else {
0320       const float dR2(dR * dR);
0321       const float reta(r.eta());
0322       const float rphi(r.phi());
0323       const float rx(r.x());
0324       const float ry(r.y());
0325       const float rz(r.z());
0326       const float fac(std::abs(zeP / rz));
0327       const float xx(rx * fac);  // xyz at endcap z
0328       const float yy(ry * fac);
0329       const float zz(rz * fac);
0330 
0331       const float xang(std::atan(xx / zz));
0332       const float lowX(zz > 0 ? zz * std::tan(xang - dR) : zz * std::tan(xang + dR));
0333       const float highX(zz > 0 ? zz * std::tan(xang + dR) : zz * std::tan(xang - dR));
0334       const float yang(std::atan(yy / zz));
0335       const float lowY(zz > 0 ? zz * std::tan(yang - dR) : zz * std::tan(yang + dR));
0336       const float highY(zz > 0 ? zz * std::tan(yang + dR) : zz * std::tan(yang - dR));
0337 
0338       const float refxlo(0 > rz ? m_xlo[0] : m_xlo[1]);
0339       const float refxhi(0 > rz ? m_xhi[0] : m_xhi[1]);
0340       const float refylo(0 > rz ? m_ylo[0] : m_ylo[1]);
0341       const float refyhi(0 > rz ? m_yhi[0] : m_yhi[1]);
0342 
0343       if (lowX < refxhi &&  // proceed if any possible overlap with the endcap
0344           lowY < refyhi && highX > refxlo && highY > refylo) {
0345         const int ix_ctr(xindex(xx, rz));
0346         const int iy_ctr(yindex(yy, rz));
0347         const int iz(rz > 0 ? 1 : -1);
0348 
0349         const int ix_hi(ix_ctr + int((highX - xx) / m_wref) + 2);
0350         const int ix_lo(ix_ctr - int((xx - lowX) / m_wref) - 2);
0351 
0352         const int iy_hi(iy_ctr + int((highY - yy) / m_wref) + 2);
0353         const int iy_lo(iy_ctr - int((yy - lowY) / m_wref) - 2);
0354 
0355         for (int kx(ix_lo); kx <= ix_hi; ++kx) {
0356           if (kx > 0 && kx <= (int)m_nref) {
0357             for (int ky(iy_lo); ky <= iy_hi; ++ky) {
0358               if (ky > 0 && ky <= (int)m_nref) {
0359                 if (EEDetId::validDetId(kx, ky, iz))  // reject invalid ids
0360                 {
0361                   const EEDetId id(kx, ky, iz);
0362                   const CaloCellGeometry* cell(&m_cellVec[id.denseIndex()]);
0363                   const float eta(cell->etaPos());
0364                   const float phi(cell->phiPos());
0365                   if (reco::deltaR2(eta, phi, reta, rphi) < dR2)
0366                     dis.insert(id);
0367                 }
0368               }
0369             }
0370           }
0371         }
0372       }
0373     }
0374   }
0375   return dis;
0376 }
0377 
0378 const EcalEndcapGeometry::OrderedListOfEBDetId* EcalEndcapGeometry::getClosestBarrelCells(EEDetId id) const {
0379   OrderedListOfEBDetId* ptr(nullptr);
0380   auto ptrVec = m_borderPtrVec.load(std::memory_order_acquire);
0381   if (!ptrVec) {
0382     if (0 != id.rawId() && nullptr != getGeometry(id)) {
0383       const float phi(370. + getGeometry(id)->getPosition().phi().degrees());
0384       const int iPhi(1 + int(phi) % 360);
0385       const int iz(id.zside());
0386       if (!m_borderMgr.load(std::memory_order_acquire)) {
0387         EZMgrFL<EBDetId>* expect = nullptr;
0388         auto ptrMgr = new EZMgrFL<EBDetId>(720 * 9, 9);
0389         bool exchanged = m_borderMgr.compare_exchange_strong(expect, ptrMgr, std::memory_order_acq_rel);
0390         if (!exchanged)
0391           delete ptrMgr;
0392       }
0393       VecOrdListEBDetIdPtr* expect = nullptr;
0394       auto ptrVec = new VecOrdListEBDetIdPtr();
0395       ptrVec->reserve(720);
0396       for (unsigned int i(0); i != 720; ++i) {
0397         const int kz(360 > i ? -1 : 1);
0398         const int iEta(kz * 85);
0399         const int iEtam1(kz * 84);
0400         const int iEtam2(kz * 83);
0401         const int jPhi(i % 360 + 1);
0402         OrderedListOfEBDetId& olist(*new OrderedListOfEBDetId(m_borderMgr.load(std::memory_order_acquire)));
0403         olist[0] = EBDetId(iEta, jPhi);
0404         olist[1] = EBDetId(iEta, myPhi(jPhi + 1));
0405         olist[2] = EBDetId(iEta, myPhi(jPhi - 1));
0406         olist[3] = EBDetId(iEtam1, jPhi);
0407         olist[4] = EBDetId(iEtam1, myPhi(jPhi + 1));
0408         olist[5] = EBDetId(iEtam1, myPhi(jPhi - 1));
0409         olist[6] = EBDetId(iEta, myPhi(jPhi + 2));
0410         olist[7] = EBDetId(iEta, myPhi(jPhi - 2));
0411         olist[8] = EBDetId(iEtam2, jPhi);
0412         ptrVec->emplace_back(&olist);
0413       }
0414       bool exchanged = m_borderPtrVec.compare_exchange_strong(expect, ptrVec, std::memory_order_acq_rel);
0415       if (!exchanged)
0416         delete ptrVec;
0417       ptrVec = m_borderPtrVec.load(std::memory_order_acquire);
0418       ptr = (*ptrVec)[(iPhi - 1) + (0 > iz ? 0 : 360)];
0419     }
0420   }
0421   return ptr;
0422 }
0423 
0424 void EcalEndcapGeometry::localCorners(Pt3DVec& lc, const CCGFloat* pv, unsigned int /*i*/, Pt3D& ref) {
0425   TruncatedPyramid::localCorners(lc, pv, ref);
0426 }
0427 
0428 void EcalEndcapGeometry::newCell(
0429     const GlobalPoint& f1, const GlobalPoint& f2, const GlobalPoint& f3, const CCGFloat* parm, const DetId& detId) {
0430   const unsigned int cellIndex(EEDetId(detId).denseIndex());
0431   m_cellVec[cellIndex] = TruncatedPyramid(cornersMgr(), f1, f2, f3, parm);
0432   addValidID(detId);
0433 }
0434 
0435 CCGFloat EcalEndcapGeometry::avgAbsZFrontFaceCenter() const {
0436   if (!m_check.load(std::memory_order_acquire)) {
0437     CCGFloat sum(0);
0438     for (unsigned int i(0); i != m_cellVec.size(); ++i) {
0439       auto cell = cellGeomPtr(i);
0440       if (nullptr != cell) {
0441         sum += fabs(cell->getPosition().z());
0442       }
0443     }
0444     m_avgZ = sum / m_cellVec.size();
0445     m_check.store(true, std::memory_order_release);
0446   }
0447   return m_avgZ;
0448 }
0449 
0450 CaloCellGeometryPtr EcalEndcapGeometry::getGeometryRawPtr(uint32_t index) const {
0451   return CaloCellGeometryPtr(m_cellVec.size() <= index || nullptr == m_cellVec[index].param() ? nullptr
0452                                                                                               : &m_cellVec[index]);
0453 }
0454 
0455 bool EcalEndcapGeometry::present(const DetId& id) const {
0456   if (id.det() == DetId::Ecal && id.subdetId() == EcalEndcap) {
0457     EEDetId eeId(id);
0458     if (EEDetId::validDetId(eeId.ix(), eeId.iy(), eeId.zside()))
0459       return true;
0460   }
0461   return false;
0462 }
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!