Back to home page

Project CMSSW displayed by LXR

 
 

    


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

0001 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0002 #include "Geometry/CaloGeometry/interface/FlatHexagon.h"
0003 #include "Geometry/CaloGeometry/interface/TruncatedPyramid.h"
0004 #include <algorithm>
0005 #include <iostream>
0006 
0007 //#define EDM_ML_DEBUG
0008 
0009 typedef FlatHexagon::CCGFloat CCGFloat;
0010 typedef FlatHexagon::Pt3D Pt3D;
0011 typedef FlatHexagon::Pt3DVec Pt3DVec;
0012 typedef FlatHexagon::Tr3D Tr3D;
0013 
0014 typedef HepGeom::Vector3D<CCGFloat> FVec3D;
0015 typedef HepGeom::Plane3D<CCGFloat> Plane3D;
0016 
0017 typedef HepGeom::Vector3D<double> DVec3D;
0018 typedef HepGeom::Plane3D<double> DPlane3D;
0019 typedef HepGeom::Point3D<double> DPt3D;
0020 
0021 static const float tolmin = 1.e-12;
0022 
0023 //----------------------------------------------------------------------
0024 
0025 FlatHexagon::FlatHexagon()
0026     : CaloCellGeometry(), m_axis(0., 0., 0.), m_corOne(0., 0., 0.), m_local(0., 0., 0.), m_global(0., 0., 0.) {}
0027 
0028 FlatHexagon::FlatHexagon(const FlatHexagon& tr) : CaloCellGeometry(tr) { *this = tr; }
0029 
0030 FlatHexagon& FlatHexagon::operator=(const FlatHexagon& tr) {
0031   CaloCellGeometry::operator=(tr);
0032   if (this != &tr) {
0033     m_axis = tr.m_axis;
0034     m_corOne = tr.m_corOne;
0035     m_local = tr.m_local;
0036     m_global = tr.m_global;
0037     m_tr = tr.m_tr;
0038   }
0039 #ifdef EDM_ML_DEBUG
0040   edm::LogVerbatim("CaloGeometry") << "FlatHexagon(Copy): Local " << m_local << " Global " << m_global << " eta "
0041                                    << etaPos() << " phi " << phiPos() << " Translation " << m_tr.getTranslation()
0042                                    << " and rotation " << m_tr.getRotation();
0043 #endif
0044   return *this;
0045 }
0046 
0047 FlatHexagon::FlatHexagon(
0048     CornersMgr* cMgr, const GlobalPoint& fCtr, const GlobalPoint& bCtr, const GlobalPoint& cor1, const CCGFloat* parV)
0049     : CaloCellGeometry(fCtr, cMgr, parV),
0050       m_axis((bCtr - fCtr).unit()),
0051       m_corOne(cor1.x(), cor1.y(), cor1.z()),
0052       m_local(0., 0., 0.) {
0053   getTransform(m_tr, nullptr);
0054   Pt3D glb = m_tr * m_local;
0055   m_global = GlobalPoint(glb.x(), glb.y(), glb.z());
0056 #ifdef EDM_ML_DEBUG
0057   edm::LogVerbatim("CaloGeometry") << "FlatHexagon: Local " << m_local << " Global " << glb << " eta " << etaPos()
0058                                    << " phi " << phiPos() << " Translation " << m_tr.getTranslation()
0059                                    << " and rotation " << m_tr.getRotation();
0060 #endif
0061 }
0062 
0063 FlatHexagon::FlatHexagon(const CornersVec& corn, const CCGFloat* par)
0064     : CaloCellGeometry(corn, par), m_corOne(corn[0].x(), corn[0].y(), corn[0].z()), m_local(0., 0., 0.) {
0065   getTransform(m_tr, nullptr);
0066   m_axis = makeAxis();
0067   Pt3D glb = m_tr * m_local;
0068   m_global = GlobalPoint(glb.x(), glb.y(), glb.z());
0069 #ifdef EDM_ML_DEBUG
0070   edm::LogVerbatim("CaloGeometry") << "FlatHexagon: Local " << m_local << " Global " << glb << " eta " << etaPos()
0071                                    << " phi " << phiPos() << " Translation " << m_tr.getTranslation()
0072                                    << " and rotation " << m_tr.getRotation();
0073 #endif
0074 }
0075 
0076 FlatHexagon::FlatHexagon(const FlatHexagon& tr, const Pt3D& local) : CaloCellGeometry(tr), m_local(local) {
0077   *this = tr;
0078   Pt3D glb = m_tr * m_local;
0079   m_global = GlobalPoint(glb.x(), glb.y(), glb.z());
0080 #ifdef EDM_ML_DEBUG
0081   edm::LogVerbatim("CaloGeometry") << "FlatHexagon: Local " << m_local << " Global " << glb << " eta " << etaPos()
0082                                    << " phi " << phiPos() << " Translation " << m_tr.getTranslation()
0083                                    << " and rotation " << m_tr.getRotation();
0084 #endif
0085 }
0086 
0087 FlatHexagon::~FlatHexagon() {}
0088 
0089 GlobalPoint FlatHexagon::getPosition(const Pt3D& local) const {
0090   Pt3D glb = m_tr * local;
0091 #ifdef EDM_ML_DEBUG
0092   edm::LogVerbatim("CaloGeometry") << "FlatHexagon::Local " << local.x() << ":" << local.y() << ":" << local.z()
0093                                    << " Global " << glb.x() << ":" << glb.y() << ":" << glb.z() << " TR " << m_tr.xx()
0094                                    << ":" << m_tr.xy() << ":" << m_tr.xz() << ":" << m_tr.yx() << ":" << m_tr.yy()
0095                                    << ":" << m_tr.yz() << ":" << m_tr.zx() << ":" << m_tr.zy() << ":" << m_tr.zz()
0096                                    << ":" << m_tr.dx() << ":" << m_tr.dy() << ":" << m_tr.dz();
0097 #endif
0098   return GlobalPoint(glb.x(), glb.y(), glb.z());
0099 }
0100 
0101 float FlatHexagon::etaSpan() const {
0102   assert(param() != nullptr);
0103   float eta1 =
0104       -std::log(std::tan(0.5 * std::atan((m_global.perp() + param()[k_R]) / std::max(tolmin, std::abs(m_global.z())))));
0105   float eta2 =
0106       -std::log(std::tan(0.5 * std::atan((m_global.perp() - param()[k_R]) / std::max(tolmin, std::abs(m_global.z())))));
0107   float dEta = std::abs(eta1 - eta2);
0108   return dEta;
0109 }
0110 
0111 float FlatHexagon::phiSpan() const {
0112   assert(param() != nullptr);
0113   float dPhi = 2.0 * std::atan(param()[k_r] / std::max(tolmin, m_global.perp()));
0114   return dPhi;
0115 }
0116 
0117 Pt3D FlatHexagon::getLocal(const GlobalPoint& global) const {
0118   Pt3D local = m_tr.inverse() * Pt3D(global.x(), global.y(), global.z());
0119 #ifdef EDM_ML_DEBUG
0120   edm::LogVerbatim("CaloGeometry") << "FlatHexagon::Global " << global.x() << ":" << global.y() << ":" << global.z()
0121                                    << " Local " << local.x() << ":" << local.y() << ":" << local.z();
0122 #endif
0123   return local;
0124 }
0125 
0126 CCGFloat FlatHexagon::getThetaAxis() const { return m_axis.theta(); }
0127 
0128 CCGFloat FlatHexagon::getPhiAxis() const { return m_axis.phi(); }
0129 
0130 const GlobalVector& FlatHexagon::axis() const { return m_axis; }
0131 
0132 void FlatHexagon::vocalCorners(Pt3DVec& vec, const CCGFloat* pv, Pt3D& ref) const { localCorners(vec, pv, ref); }
0133 
0134 void FlatHexagon::createCorners(const std::vector<CCGFloat>& pv, const Tr3D& tr, std::vector<GlobalPoint>& co) {
0135   assert(3 == pv.size());
0136   assert(ncorner_ == co.size());
0137 
0138   Pt3DVec ko(ncorner_, Pt3D(0, 0, 0));
0139 
0140   Pt3D tmp;
0141   Pt3DVec to(ncorner_, Pt3D(0, 0, 0));
0142   localCorners(to, &pv.front(), tmp);
0143 
0144   for (unsigned int i(0); i != ncorner_; ++i) {
0145     ko[i] = tr * to[i];  // apply transformation
0146     const Pt3D& p(ko[i]);
0147     co[i] = GlobalPoint(p.x(), p.y(), p.z());
0148 #ifdef EDM_ML_DEBUG
0149     edm::LogVerbatim("CaloGeometry") << "Corner[" << i << "] = " << co[i];
0150 #endif
0151   }
0152 }
0153 
0154 void FlatHexagon::localCorners(Pt3DVec& lc, const CCGFloat* pv, Pt3D& ref) {
0155   assert(nullptr != pv);
0156   assert(ncorner_ == lc.size());
0157 
0158   const CCGFloat dz(pv[FlatHexagon::k_dZ]);
0159   const CCGFloat r(pv[FlatHexagon::k_r]);
0160   const CCGFloat R(pv[FlatHexagon::k_R]);
0161 
0162   lc[0] = Pt3D(0, -2 * R, -dz);  // (0,-,-)
0163   lc[1] = Pt3D(-r, -R, -dz);     // (-,-,-)
0164   lc[2] = Pt3D(-r, R, -dz);      // (-,+,-)
0165   lc[3] = Pt3D(0, 2 * R, -dz);   // (0,+,-)
0166   lc[4] = Pt3D(r, R, -dz);       // (+,+,-)
0167   lc[5] = Pt3D(r, -R, -dz);      // (+,-,-)
0168   lc[6] = Pt3D(0, -2 * R, dz);   // (0,-,+)
0169   lc[7] = Pt3D(-r, -R, dz);      // (-,-,+)
0170   lc[8] = Pt3D(-r, R, dz);       // (-,+,+)
0171   lc[9] = Pt3D(0, 2 * R, dz);    // (0,+,+)
0172   lc[10] = Pt3D(r, R, dz);       // (+,+,+)
0173   lc[11] = Pt3D(r, -R, dz);      // (+,-,+)
0174 
0175   ref = oneBySix_ * (lc[0] + lc[1] + lc[2] + lc[3] + lc[4] + lc[5]);
0176 #ifdef EDM_ML_DEBUG
0177   edm::LogVerbatim("CaloGeometry") << "Ref " << ref << " Local Corners " << lc[0] << "|" << lc[1] << "|" << lc[2] << "|"
0178                                    << lc[3] << "|" << lc[4] << "|" << lc[5] << "|" << lc[6] << "|" << lc[7] << "|"
0179                                    << lc[8] << "|" << lc[9] << "|" << lc[10] << "|" << lc[11];
0180 #endif
0181 }
0182 
0183 void FlatHexagon::getTransform(Tr3D& tr, Pt3DVec* lptr) const {
0184   const GlobalPoint& p(CaloCellGeometry::getPosition());
0185   const Pt3D gFront(p.x(), p.y(), p.z());
0186   const DPt3D dgFront(p.x(), p.y(), p.z());
0187 
0188   Pt3D lFront;
0189   assert(nullptr != param());
0190   std::vector<Pt3D> lc(ncorner_, Pt3D(0, 0, 0));
0191   localCorners(lc, param(), lFront);
0192 
0193   // figure out if reflction volume or not
0194 
0195   Pt3D lBack((lc[6] + lc[7] + lc[8] + lc[9] + lc[10] + lc[11]) / 6.0);
0196 
0197   const DPt3D dlFront(lFront.x(), lFront.y(), lFront.z());
0198   const DPt3D dlBack(lBack.x(), lBack.y(), lBack.z());
0199   const DPt3D dlOne(lc[0].x(), lc[0].y(), lc[0].z());
0200 
0201   const FVec3D dgAxis(axis().x(), axis().y(), axis().z());
0202 
0203   const DPt3D dmOne(m_corOne.x(), m_corOne.y(), m_corOne.z());
0204 
0205   const DPt3D dgBack(dgFront + (dlBack - dlFront).mag() * dgAxis);
0206   DPt3D dgOne(dgFront + (dlOne - dlFront).mag() * (dmOne - dgFront).unit());
0207 
0208   const double dlangle((dlBack - dlFront).angle(dlOne - dlFront));
0209   const double dgangle((dgBack - dgFront).angle(dgOne - dgFront));
0210   const double dangle(dlangle - dgangle);
0211 
0212   if (1.e-6 < fabs(dangle)) {  //guard against precision problems
0213     const DPlane3D dgPl(dgFront, dgOne, dgBack);
0214     const DPt3D dp2(dgFront + dgPl.normal().unit());
0215 
0216     DPt3D dgOld(dgOne);
0217 
0218     dgOne = (dgFront + HepGeom::Rotate3D(-dangle, dgFront, dp2) * DVec3D(dgOld - dgFront));
0219   }
0220 
0221   tr = Tr3D(dlFront, dlBack, dlOne, dgFront, dgBack, dgOne);
0222 
0223   if (nullptr != lptr)
0224     (*lptr) = lc;
0225 }
0226 
0227 void FlatHexagon::initCorners(CaloCellGeometry::CornersVec& co) {
0228   if (co.uninitialized()) {
0229     CornersVec& corners(co);
0230     Pt3DVec lc;
0231     Tr3D tr;
0232     getTransform(tr, &lc);
0233 
0234     for (unsigned int i(0); i != ncorner_; ++i) {
0235       const Pt3D corn(tr * lc[i]);
0236       corners[i] = GlobalPoint(corn.x(), corn.y(), corn.z());
0237     }
0238   }
0239 }
0240 
0241 GlobalVector FlatHexagon::makeAxis() { return GlobalVector(backCtr() - getPosition()).unit(); }
0242 
0243 GlobalPoint FlatHexagon::backCtr() const {
0244   float dz =
0245       (getCorners()[ncornerBy2_].z() > getCorners()[0].z()) ? param()[FlatHexagon::k_dZ] : -param()[FlatHexagon::k_dZ];
0246   Pt3D local_b(m_local.x(), m_local.y(), m_local.z() + dz);
0247   Pt3D global_b = m_tr * local_b;
0248   GlobalPoint global(global_b.x(), global_b.y(), global_b.z());
0249   return global;
0250 }
0251 
0252 //----------------------------------------------------------------------
0253 //----------------------------------------------------------------------
0254 
0255 std::ostream& operator<<(std::ostream& s, const FlatHexagon& cell) {
0256   s << "Center: " << cell.getPosition() << " eta " << cell.etaPos() << " phi " << cell.phiPos() << std::endl;
0257   s << "Axis: " << cell.getThetaAxis() << " " << cell.getPhiAxis() << std::endl;
0258   const CaloCellGeometry::CornersVec& corners(cell.getCorners());
0259   for (unsigned int i = 0; i != corners.size(); ++i) {
0260     s << "Corner: " << corners[i] << std::endl;
0261   }
0262   return s;
0263 }