File indexing completed on 2024-10-25 23:57:05
0001 #include "Geometry/CaloGeometry/interface/CaloGenericDetId.h"
0002 #include "Geometry/CaloGeometry/interface/TruncatedPyramid.h"
0003 #include "Geometry/EcalAlgo/interface/EcalBarrelGeometry.h"
0004 #include "Geometry/CaloGeometry/interface/CaloCellGeometry.h"
0005 #include "FWCore/Utilities/interface/Exception.h"
0006 #include "DataFormats/EcalDetId/interface/EBDetId.h"
0007
0008 #include <CLHEP/Geometry/Point3D.h>
0009 #include <CLHEP/Geometry/Plane3D.h>
0010 #include <CLHEP/Geometry/Vector3D.h>
0011
0012 #include <iomanip>
0013 #include <iostream>
0014
0015 typedef CaloCellGeometry::CCGFloat CCGFloat;
0016 typedef CaloCellGeometry::Pt3D Pt3D;
0017 typedef CaloCellGeometry::Pt3DVec Pt3DVec;
0018 typedef HepGeom::Plane3D<CCGFloat> Pl3D;
0019
0020 EcalBarrelGeometry::EcalBarrelGeometry()
0021 : _nnxtalEta(85),
0022 _nnxtalPhi(360),
0023 _PhiBaskets(18),
0024 m_borderMgr(nullptr),
0025 m_borderPtrVec(nullptr),
0026 m_radius(-1.),
0027 m_check(false),
0028 m_cellVec(k_NumberOfCellsForCorners) {
0029 const int neba[] = {25, 45, 65, 85};
0030 _EtaBaskets = std::vector<int>(neba, neba + 4);
0031 }
0032
0033 EcalBarrelGeometry::~EcalBarrelGeometry() {
0034 if (m_borderPtrVec) {
0035 auto ptr = m_borderPtrVec.load(std::memory_order_acquire);
0036 for (auto& v : (*ptr)) {
0037 delete v;
0038 v = nullptr;
0039 }
0040 delete m_borderPtrVec.load();
0041 }
0042 delete m_borderMgr.load();
0043 }
0044
0045 unsigned int EcalBarrelGeometry::alignmentTransformIndexLocal(const DetId& id) {
0046 const CaloGenericDetId gid(id);
0047
0048 assert(gid.isEB());
0049
0050 unsigned int index(EBDetId(id).ism() - 1);
0051
0052 return index;
0053 }
0054
0055 DetId EcalBarrelGeometry::detIdFromLocalAlignmentIndex(unsigned int iLoc) {
0056 return EBDetId(iLoc + 1, 1, EBDetId::SMCRYSTALMODE);
0057 }
0058
0059 unsigned int EcalBarrelGeometry::alignmentTransformIndexGlobal(const DetId& ) {
0060 return (unsigned int)DetId::Ecal - 1;
0061 }
0062
0063 DetId EcalBarrelGeometry::getClosestCell(const GlobalPoint& r) const {
0064
0065 int leverx = 1;
0066 int levery = 1;
0067 CCGFloat pointz = r.z();
0068 int zbin = 1;
0069 if (pointz < 0)
0070 zbin = -1;
0071
0072
0073 CCGFloat pointeta = r.eta();
0074
0075 CCGFloat deta = 999.;
0076 int etabin = 1;
0077
0078 int guessed_eta = (int)(fabs(pointeta) / 0.0174) + 1;
0079 int guessed_eta_begin = guessed_eta - 1;
0080 int guessed_eta_end = guessed_eta + 1;
0081 if (guessed_eta_begin < 1)
0082 guessed_eta_begin = 1;
0083 if (guessed_eta_end > 85)
0084 guessed_eta_end = 85;
0085
0086 for (int bin = guessed_eta_begin; bin <= guessed_eta_end; bin++) {
0087 try {
0088 if (!present(EBDetId(zbin * bin, 1, EBDetId::ETAPHIMODE)))
0089 continue;
0090
0091 CCGFloat eta = getGeometry(EBDetId(zbin * bin, 1, EBDetId::ETAPHIMODE))->etaPos();
0092
0093 if (fabs(pointeta - eta) < deta) {
0094 deta = fabs(pointeta - eta);
0095 etabin = bin;
0096 } else
0097 break;
0098 } catch (cms::Exception& e) {
0099 }
0100 }
0101
0102
0103 constexpr CCGFloat twopi = M_PI + M_PI;
0104
0105 constexpr CCGFloat tilt = twopi / 36.;
0106
0107 CCGFloat pointphi = r.phi() + tilt;
0108
0109
0110 if (pointphi > twopi)
0111 pointphi -= twopi;
0112 if (pointphi < 0)
0113 pointphi += twopi;
0114
0115
0116 int phibin = static_cast<int>(pointphi / (twopi / _nnxtalPhi)) + 1;
0117
0118
0119
0120
0121
0122
0123 try {
0124 EBDetId myCell(zbin * etabin, phibin, EBDetId::ETAPHIMODE);
0125
0126 if (!present(myCell))
0127 return DetId(0);
0128
0129 Pt3D A;
0130 Pt3D B;
0131 Pt3D C;
0132 Pt3D point(r.x(), r.y(), r.z());
0133
0134
0135
0136
0137
0138
0139 CCGFloat history[4]{0.f};
0140
0141
0142
0143 int start = 1;
0144 int counter = 0;
0145
0146 while (leverx == 1 || levery == 1) {
0147 leverx = 0;
0148 levery = 0;
0149 const CaloCellGeometry::CornersVec& corners(getGeometry(myCell)->getCorners());
0150 CCGFloat SS[4];
0151
0152
0153 for (short i = 0; i < 4; ++i) {
0154 A = Pt3D(corners[i % 4].x(), corners[i % 4].y(), corners[i % 4].z());
0155 B = Pt3D(corners[(i + 1) % 4].x(), corners[(i + 1) % 4].y(), corners[(i + 1) % 4].z());
0156 C = Pt3D(corners[4 + (i + 1) % 4].x(), corners[4 + (i + 1) % 4].y(), corners[4 + (i + 1) % 4].z());
0157 Pl3D plane(A, B, C);
0158 plane.normalize();
0159 CCGFloat distance = plane.distance(point);
0160 if (plane.d() > 0.)
0161 distance = -distance;
0162 if (corners[0].z() < 0.)
0163 distance = -distance;
0164 SS[i] = distance;
0165 }
0166
0167
0168
0169
0170
0171 if ((SS[0] > 0. && SS[2] > 0.) || (SS[0] < 0. && SS[2] < 0.)) {
0172 levery = 1;
0173 if (history[0] > 0. && history[2] > 0. && SS[0] < 0 && SS[2] < 0 &&
0174 (fabs(SS[0]) + fabs(SS[2])) > (fabs(history[0]) + fabs(history[2])))
0175 levery = 0;
0176 if (history[0] < 0. && history[2] < 0. && SS[0] > 0 && SS[2] > 0 &&
0177 (fabs(SS[0]) + fabs(SS[2])) > (fabs(history[0]) + fabs(history[2])))
0178 levery = 0;
0179
0180 if (SS[0] > 0.) {
0181 EBDetId nextPoint;
0182 if (myCell.iphi() == EBDetId::MIN_IPHI)
0183 nextPoint = EBDetId(myCell.ieta(), EBDetId::MAX_IPHI);
0184 else
0185 nextPoint = EBDetId(myCell.ieta(), myCell.iphi() - 1);
0186 if (present(nextPoint))
0187 myCell = nextPoint;
0188 else
0189 levery = 0;
0190 } else {
0191 EBDetId nextPoint;
0192 if (myCell.iphi() == EBDetId::MAX_IPHI)
0193 nextPoint = EBDetId(myCell.ieta(), EBDetId::MIN_IPHI);
0194 else
0195 nextPoint = EBDetId(myCell.ieta(), myCell.iphi() + 1);
0196 if (present(nextPoint))
0197 myCell = nextPoint;
0198 else
0199 levery = 0;
0200 }
0201 }
0202
0203 if (((SS[1] > 0. && SS[3] > 0.) || (SS[1] < 0. && SS[3] < 0.)) && start == 1) {
0204 leverx = 1;
0205
0206 if (history[1] > 0. && history[3] > 0. && SS[1] < 0 && SS[3] < 0 &&
0207 (fabs(SS[1]) + fabs(SS[3])) > (fabs(history[1]) + fabs(history[3]))) {
0208 leverx = 0;
0209 start = 0;
0210 }
0211
0212 if (history[1] < 0. && history[3] < 0. && SS[1] > 0 && SS[3] > 0 &&
0213 (fabs(SS[1]) + fabs(SS[3])) > (fabs(history[1]) + fabs(history[3]))) {
0214 leverx = 0;
0215 start = 0;
0216 }
0217
0218 if (SS[1] > 0.) {
0219 EBDetId nextPoint;
0220 if (myCell.ieta() == -1)
0221 nextPoint = EBDetId(1, myCell.iphi());
0222 else {
0223 int nieta = myCell.ieta() + 1;
0224 if (nieta == 86)
0225 nieta = 85;
0226 nextPoint = EBDetId(nieta, myCell.iphi());
0227 }
0228 if (present(nextPoint))
0229 myCell = nextPoint;
0230 else
0231 leverx = 0;
0232 } else {
0233 EBDetId nextPoint;
0234 if (myCell.ieta() == 1)
0235 nextPoint = EBDetId(-1, myCell.iphi());
0236 else {
0237 int nieta = myCell.ieta() - 1;
0238 if (nieta == -86)
0239 nieta = -85;
0240 nextPoint = EBDetId(nieta, myCell.iphi());
0241 }
0242 if (present(nextPoint))
0243 myCell = nextPoint;
0244 else
0245 leverx = 0;
0246 }
0247 }
0248
0249
0250
0251 std::copy(SS, SS + 4, history);
0252
0253 counter++;
0254 if (counter == 10) {
0255 leverx = 0;
0256 levery = 0;
0257 }
0258 }
0259
0260 return DetId(myCell);
0261 } catch (cms::Exception& e) {
0262 return DetId(0);
0263 }
0264 }
0265
0266 CaloSubdetectorGeometry::DetIdSet EcalBarrelGeometry::getCells(const GlobalPoint& r, double dR) const {
0267 constexpr int maxphi(EBDetId::MAX_IPHI);
0268 constexpr int maxeta(EBDetId::MAX_IETA);
0269 constexpr float scale(maxphi / (2 * M_PI));
0270
0271 CaloSubdetectorGeometry::DetIdSet dis;
0272
0273 if (0.000001 < dR) {
0274 if (dR > M_PI / 2.)
0275 {
0276 dis = CaloSubdetectorGeometry::getCells(r, dR);
0277 } else {
0278 const float dR2(dR * dR);
0279 const float reta(r.eta());
0280 const float rz(r.z());
0281 const float rphi(r.phi());
0282 const float lowEta(reta - dR);
0283 const float highEta(reta + dR);
0284
0285 if (highEta > -1.5 && lowEta < 1.5)
0286 {
0287 const int ieta_center(int(reta * scale + ((rz < 0) ? (-1) : (1))));
0288 const float phi(rphi < 0 ? rphi + float(2 * M_PI) : rphi);
0289 const int iphi_center(int(phi * scale + 11.f));
0290
0291 const float fr(dR * scale);
0292 const float frp(1.08f * fr + 1.f);
0293 const float frm(0.92f * fr - 1.f);
0294 const int idr((int)frp);
0295 const int idr2p((int)(frp * frp));
0296 const int idr2m(frm > 0 ? int(frm * frm) : 0);
0297
0298 for (int de(-idr); de <= idr; ++de)
0299 {
0300 int ieta(de + ieta_center);
0301
0302 if (std::abs(ieta) <= maxeta && ieta != 0)
0303 {
0304 const int de2(de * de);
0305 for (int dp(-idr); dp <= idr; ++dp)
0306 {
0307 const int irange2(dp * dp + de2);
0308
0309 if (irange2 <= idr2p)
0310 {
0311 const int iphi((iphi_center + dp + maxphi - 1) % maxphi + 1);
0312
0313 if (iphi != 0) {
0314 const EBDetId id(ieta, iphi);
0315
0316 bool ok(irange2 < idr2m);
0317
0318 if (!ok)
0319 {
0320 const CaloCellGeometry* cell(&m_cellVec[id.denseIndex()]);
0321 const float eta(cell->etaPos());
0322 const float phi(cell->phiPos());
0323 ok = (reco::deltaR2(eta, phi, reta, rphi) < dR2);
0324 }
0325 if (ok)
0326 dis.insert(id);
0327 }
0328 }
0329 }
0330 }
0331 }
0332 }
0333 }
0334 }
0335 return dis;
0336 }
0337
0338 const EcalBarrelGeometry::OrderedListOfEEDetId* EcalBarrelGeometry::getClosestEndcapCells(EBDetId id) const {
0339 OrderedListOfEEDetId* ptr(nullptr);
0340 auto ptrVec = m_borderPtrVec.load(std::memory_order_acquire);
0341 if (!ptrVec) {
0342 if (0 != id.rawId()) {
0343 const int iPhi(id.iphi());
0344 const int iz(id.ieta() > 0 ? 1 : -1);
0345 const EEDetId eeid(EEDetId::idOuterRing(iPhi, iz));
0346 const int iq(eeid.iquadrant());
0347 const int xout(1 == iq || 4 == iq ? 1 : -1);
0348 const int yout(1 == iq || 2 == iq ? 1 : -1);
0349 if (!m_borderMgr.load(std::memory_order_acquire)) {
0350 EZMgrFL<EEDetId>* expect = nullptr;
0351 auto ptrMgr = new EZMgrFL<EEDetId>(720 * 9, 9);
0352 bool exchanged = m_borderMgr.compare_exchange_strong(expect, ptrMgr, std::memory_order_acq_rel);
0353 if (!exchanged)
0354 delete ptrMgr;
0355 }
0356 VecOrdListEEDetIdPtr* expect = nullptr;
0357 auto ptrVec = new VecOrdListEEDetIdPtr();
0358 ptrVec->reserve(720);
0359 for (unsigned int i(0); i != 720; ++i) {
0360 const int kz(360 > i ? -1 : 1);
0361 const EEDetId eeid(EEDetId::idOuterRing(i % 360 + 1, kz));
0362
0363 const int jx(eeid.ix());
0364 const int jy(eeid.iy());
0365
0366 OrderedListOfEEDetId& olist(*new OrderedListOfEEDetId(m_borderMgr.load(std::memory_order_acquire)));
0367 int il(0);
0368
0369 for (unsigned int k(1); k <= 25; ++k) {
0370 const int kx(1 == k || 2 == k || 3 == k || 12 == k || 13 == k
0371 ? 0
0372 : (4 == k || 6 == k || 8 == k || 15 == k || 20 == k
0373 ? 1
0374 : (5 == k || 7 == k || 9 == k || 16 == k || 19 == k
0375 ? -1
0376 : (10 == k || 14 == k || 21 == k || 22 == k || 25 == k ? 2 : -2))));
0377 const int ky(1 == k || 4 == k || 5 == k || 10 == k || 11 == k
0378 ? 0
0379 : (2 == k || 6 == k || 7 == k || 14 == k || 17 == k
0380 ? 1
0381 : (3 == k || 8 == k || 9 == k || 18 == k || 21 == k
0382 ? -1
0383 : (12 == k || 15 == k || 16 == k || 22 == k || 23 == k ? 2 : -2))));
0384
0385 if (8 >= il && EEDetId::validDetId(jx + kx * xout, jy + ky * yout, kz)) {
0386 olist[il++] = EEDetId(jx + kx * xout, jy + ky * yout, kz);
0387 }
0388 }
0389 ptrVec->emplace_back(&olist);
0390 }
0391 bool exchanged = m_borderPtrVec.compare_exchange_strong(expect, ptrVec, std::memory_order_acq_rel);
0392 if (!exchanged)
0393 delete ptrVec;
0394 ptrVec = m_borderPtrVec.load(std::memory_order_acquire);
0395 ptr = (*ptrVec)[iPhi - 1 + (0 > iz ? 0 : 360)];
0396 }
0397 }
0398 return ptr;
0399 }
0400
0401 void EcalBarrelGeometry::localCorners(Pt3DVec& lc, const CCGFloat* pv, unsigned int i, Pt3D& ref) {
0402 const bool negz(EBDetId::kSizeForDenseIndexing / 2 > i);
0403 const bool odd(1 == i % 2);
0404
0405 if (((negz && !odd) || (!negz && odd))) {
0406 TruncatedPyramid::localCornersReflection(lc, pv, ref);
0407 } else {
0408 TruncatedPyramid::localCornersSwap(lc, pv, ref);
0409 }
0410 }
0411
0412 void EcalBarrelGeometry::newCell(
0413 const GlobalPoint& f1, const GlobalPoint& f2, const GlobalPoint& f3, const CCGFloat* parm, const DetId& detId) {
0414 const unsigned int cellIndex(EBDetId(detId).denseIndex());
0415 m_cellVec[cellIndex] = TruncatedPyramid(cornersMgr(), f1, f2, f3, parm);
0416 addValidID(detId);
0417 }
0418
0419 CCGFloat EcalBarrelGeometry::avgRadiusXYFrontFaceCenter() const {
0420 if (!m_check.load(std::memory_order_acquire)) {
0421 CCGFloat sum(0);
0422 for (uint32_t i(0); i != m_cellVec.size(); ++i) {
0423 auto cell(cellGeomPtr(i));
0424 if (nullptr != cell) {
0425 const GlobalPoint& pos(cell->getPosition());
0426 sum += pos.perp();
0427 }
0428 }
0429 m_radius = sum / m_cellVec.size();
0430 m_check.store(true, std::memory_order_release);
0431 }
0432 return m_radius;
0433 }
0434
0435 CaloCellGeometryPtr EcalBarrelGeometry::getGeometryRawPtr(uint32_t index) const {
0436
0437 return CaloCellGeometryPtr(m_cellVec.size() <= index || nullptr == m_cellVec[index].param() ? nullptr
0438 : &m_cellVec[index]);
0439 }
0440
0441 bool EcalBarrelGeometry::present(const DetId& id) const {
0442 if (id.det() == DetId::Ecal && id.subdetId() == EcalBarrel) {
0443 EBDetId ebId(id);
0444 if (EBDetId::validDetId(ebId.ieta(), ebId.iphi()))
0445 return true;
0446 }
0447 return false;
0448 }