CaloGeometry/src/FlatTrd.cc

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