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& ) {
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
0162
0163
0164
0165
0166
0167
0168
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);
0194 } else
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();
0207 }
0208
0209
0210 DetId EcalEndcapGeometry::getClosestCell(const GlobalPoint& r) const {
0211 try {
0212 EEDetId mycellID(gId(r.x(), r.y(), r.z()));
0213
0214 if (EEDetId::validDetId(mycellID.ix(), mycellID.iy(), mycellID.zside())) {
0215
0216
0217 Pt3D A;
0218 Pt3D B;
0219 Pt3D C;
0220 Pt3D point(r.x(), r.y(), r.z());
0221
0222
0223
0224
0225 CCGFloat x, y, z;
0226 unsigned offset = 0;
0227 int zsign = 1;
0228
0229
0230
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
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
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;
0309 if (0.000001 < dR) {
0310 if (dR > M_PI / 2.)
0311 {
0312 dis = CaloSubdetectorGeometry::getCells(r, dR);
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);
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 &&
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))
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 , 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
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 }