Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2021-02-14 13:07:58 

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