EcalAlgo/src/EcalBarrelGeometry.cc

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& /*id*/) {
0060   return (unsigned int)DetId::Ecal - 1;
0061 }
0062 // Get closest cell, etc...
0063 DetId EcalBarrelGeometry::getClosestCell(const GlobalPoint& r) const {
0064   // z is the easy one
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   // Now find the closest eta
0073   CCGFloat pointeta = r.eta();
0074   //  double eta;
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   // Now the closest phi. always same number of phi bins(!?)
0103   constexpr CCGFloat twopi = M_PI + M_PI;
0104   // 10 degree tilt
0105   constexpr CCGFloat tilt = twopi / 36.;
0106
0107   CCGFloat pointphi = r.phi() + tilt;
0108
0109   // put phi in correct range (0->2pi)
0110   if (pointphi > twopi)
0111     pointphi -= twopi;
0112   if (pointphi < 0)
0113     pointphi += twopi;
0114
0115   //calculate phi bin, distinguish + and - eta
0116   int phibin = static_cast<int>(pointphi / (twopi / _nnxtalPhi)) + 1;
0117   //   if(point.z()<0.0)
0118   //     {
0119   //       phibin = nxtalPhi/2 - 1 - phibin;
0120   //       if(phibin<0)
0121   //         phibin += nxtalPhi;
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     // D.K. : equation of plane : AA*x+BB*y+CC*z+DD=0;
0135     // finding equation for each edge
0136
0137     // Since the point can lie between crystals, it is necessary to keep track of the movements
0138     // to avoid infinite loops
0139     CCGFloat history[4]{0.f};
0140
0141     //
0142     // stop movement in eta direction when closest cell was found (point between crystals)
0143     int start = 1;
0144     int counter = 0;
0145     // Moving until find closest crystal in eta and phi directions (leverx and levery)
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       // compute the distance of the point with respect of the 4 crystal lateral planes
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       // SS's - normals
0168       // check position of the point with respect to opposite side of crystal
0169       // if SS's have opposite sign, the  point lies inside that crystal
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       // Update the history. If the point lies between crystals, the closest one
0250       // is returned
0251       std::copy(SS, SS + 4, history);
0252
0253       counter++;
0254       if (counter == 10) {
0255         leverx = 0;
0256         levery = 0;
0257       }
0258     }
0259     // D.K. if point lies netween cells, take a closest cell.
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));  // angle to index
0270
0271   CaloSubdetectorGeometry::DetIdSet dis;  // this is the return object
0272
0273   if (0.000001 < dR) {
0274     if (dR > M_PI / 2.)  // this version needs "small" dR
0275     {
0276       dis = CaloSubdetectorGeometry::getCells(r, dR);  // base class version
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)  // in barrel
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));  // phi=-9.4deg is iphi=1
0290
0291         const float fr(dR * scale);         // # crystal widths in dR
0292         const float frp(1.08f * fr + 1.f);  // conservatively above fr
0293         const float frm(0.92f * fr - 1.f);  // conservatively below fr
0294         const int idr((int)frp);            // integerize
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)  // over eta limits
0299         {
0300           int ieta(de + ieta_center);
0301
0302           if (std::abs(ieta) <= maxeta && ieta != 0)  // eta is in EB
0303           {
0304             const int de2(de * de);
0305             for (int dp(-idr); dp <= idr; ++dp)  // over phi limits
0306             {
0307               const int irange2(dp * dp + de2);
0308
0309               if (irange2 <= idr2p)  // cut off corners that must be too far away
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);  // no more calculation necessary if inside this radius
0317
0318                   if (!ok)  // if not ok, then we have to test this cell for being inside cone
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   // Modify the RawPtr class
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 }