CaloGeometry/src/TruncatedPyramid.cc

0001 #include "Geometry/CaloGeometry/interface/TruncatedPyramid.h"
0002 #include <algorithm>
0003 #include <iostream>
0004
0005 typedef TruncatedPyramid::CCGFloat CCGFloat;
0006 typedef TruncatedPyramid::Pt3D Pt3D;
0007 typedef TruncatedPyramid::Pt3DVec Pt3DVec;
0008 typedef TruncatedPyramid::Tr3D Tr3D;
0009
0010 typedef HepGeom::Vector3D<CCGFloat> FVec3D;
0011 typedef HepGeom::Plane3D<CCGFloat> Plane3D;
0012
0013 typedef HepGeom::Vector3D<double> DVec3D;
0014 typedef HepGeom::Plane3D<double> DPlane3D;
0015 typedef HepGeom::Point3D<double> DPt3D;
0016
0017 //----------------------------------------------------------------------
0018
0019 TruncatedPyramid::TruncatedPyramid() : CaloCellGeometry(), m_axis(0., 0., 0.), m_corOne(0., 0., 0.) {}
0020
0021 TruncatedPyramid::TruncatedPyramid(const TruncatedPyramid& tr) : CaloCellGeometry(tr) { *this = tr; }
0022
0023 TruncatedPyramid& TruncatedPyramid::operator=(const TruncatedPyramid& tr) {
0024   CaloCellGeometry::operator=(tr);
0025   if (this != &tr) {
0026     m_axis = tr.m_axis;
0027     m_corOne = tr.m_corOne;
0028   }
0029   return *this;
0030 }
0031
0032 TruncatedPyramid::TruncatedPyramid(
0033     CornersMgr* cMgr, const GlobalPoint& fCtr, const GlobalPoint& bCtr, const GlobalPoint& cor1, const CCGFloat* parV)
0034     : CaloCellGeometry(fCtr, cMgr, parV), m_axis((bCtr - fCtr).unit()), m_corOne(cor1.x(), cor1.y(), cor1.z()) {
0035   initSpan();
0036 }
0037
0038 TruncatedPyramid::TruncatedPyramid(const CornersVec& corn, const CCGFloat* par)
0039     : CaloCellGeometry(corn, par), m_axis(makeAxis()), m_corOne(corn[0].x(), corn[0].y(), corn[0].z()) {
0040   initSpan();
0041 }
0042
0043 TruncatedPyramid::~TruncatedPyramid() {}
0044
0045 GlobalPoint TruncatedPyramid::getPosition(CCGFloat depth) const {
0046   return CaloCellGeometry::getPosition() + depth * m_axis;
0047 }
0048
0049 CCGFloat TruncatedPyramid::getThetaAxis() const { return m_axis.theta(); }
0050
0051 CCGFloat TruncatedPyramid::getPhiAxis() const { return m_axis.phi(); }
0052
0053 const GlobalVector& TruncatedPyramid::axis() const { return m_axis; }
0054
0055 void TruncatedPyramid::vocalCorners(Pt3DVec& vec, const CCGFloat* pv, Pt3D& ref) const { localCorners(vec, pv, ref); }
0056
0057 GlobalVector TruncatedPyramid::makeAxis() { return GlobalVector(backCtr() - CaloCellGeometry::getPosition()).unit(); }
0058
0059 const GlobalPoint TruncatedPyramid::backCtr() const {
0060   return GlobalPoint(0.25 * (getCorners()[4].x() + getCorners()[5].x() + getCorners()[6].x() + getCorners()[7].x()),
0061                      0.25 * (getCorners()[4].y() + getCorners()[5].y() + getCorners()[6].y() + getCorners()[7].y()),
0062                      0.25 * (getCorners()[4].z() + getCorners()[5].z() + getCorners()[6].z() + getCorners()[7].z()));
0063 }
0064
0065 void TruncatedPyramid::getTransform(Tr3D& tr, Pt3DVec* lptr) const {
0066   const GlobalPoint& p(CaloCellGeometry::getPosition());
0067   const Pt3D gFront(p.x(), p.y(), p.z());
0068   const DPt3D dgFront(p.x(), p.y(), p.z());
0069
0070   const double dz(param()[0]);
0071
0072   Pt3D lFront;
0073   assert(nullptr != param());
0074   std::vector<Pt3D> lc(8, Pt3D(0, 0, 0));
0075   if (11.2 > dz) {
0076     localCorners(lc, param(), lFront);
0077   } else {
0078     localCornersSwap(lc, param(), lFront);
0079   }
0080
0081   // figure out if reflction volume or not
0082
0083   Pt3D lBack(0.25 * (lc[4] + lc[5] + lc[6] + lc[7]));
0084
0085   const double disl((lFront - lc[0]).mag());
0086   const double disr((lFront - lc[3]).mag());
0087   const double disg((gFront - m_corOne).mag());
0088
0089   const double dell(fabs(disg - disl));
0090   const double delr(fabs(disg - disr));
0091
0092   if (11.2 < dz && delr < dell)  // reflection volume if true
0093   {
0094     localCornersReflection(lc, param(), lFront);
0095     lBack = 0.25 * (lc[4] + lc[5] + lc[6] + lc[7]);
0096   }
0097
0098   const DPt3D dlFront(lFront.x(), lFront.y(), lFront.z());
0099   const DPt3D dlBack(lBack.x(), lBack.y(), lBack.z());
0100   const DPt3D dlOne(lc[0].x(), lc[0].y(), lc[0].z());
0101
0102   const FVec3D dgAxis(axis().x(), axis().y(), axis().z());
0103
0104   const DPt3D dmOne(m_corOne.x(), m_corOne.y(), m_corOne.z());
0105
0106   const DPt3D dgBack(dgFront + (dlBack - dlFront).mag() * dgAxis);
0107   DPt3D dgOne(dgFront + (dlOne - dlFront).mag() * (dmOne - dgFront).unit());
0108
0109   const double dlangle((dlBack - dlFront).angle(dlOne - dlFront));
0110   const double dgangle((dgBack - dgFront).angle(dgOne - dgFront));
0111   const double dangle(dlangle - dgangle);
0112
0113   if (1.e-6 < fabs(dangle))  //guard against precision problems
0114   {
0115     const DPlane3D dgPl(dgFront, dgOne, dgBack);
0116     const DPt3D dp2(dgFront + dgPl.normal().unit());
0117
0118     DPt3D dgOld(dgOne);
0119
0120     dgOne = (dgFront + HepGeom::Rotate3D(-dangle, dgFront, dp2) * DVec3D(dgOld - dgFront));
0121   }
0122
0123   tr = Tr3D(dlFront, dlBack, dlOne, dgFront, dgBack, dgOne);
0124
0125   if (nullptr != lptr)
0126     (*lptr) = lc;
0127 }
0128
0129 void TruncatedPyramid::initCorners(CaloCellGeometry::CornersVec& corners) {
0130   if (corners.uninitialized()) {
0131     Pt3DVec lc;
0132
0133     Tr3D tr;
0134     getTransform(tr, &lc);
0135
0136     for (unsigned int i(0); i != 8; ++i) {
0137       const Pt3D corn(tr * lc[i]);
0138       corners[i] = GlobalPoint(corn.x(), corn.y(), corn.z());
0139     }
0140   }
0141 }
0142
0143 namespace truncPyr {
0144   Pt3D refl(const Pt3D& p) { return Pt3D(-p.x(), p.y(), p.z()); }
0145 }  // namespace truncPyr
0146
0147 void TruncatedPyramid::localCornersReflection(Pt3DVec& lc, const CCGFloat* pv, Pt3D& ref) {
0148   //   using namespace truncPyr ;
0149   localCorners(lc, pv, ref);
0150   Pt3D tmp;
0151   /*
0152    tmp   = lc[0] ;
0153    lc[0] = refl( lc[2] ) ;
0154    lc[2] = refl( tmp   ) ;
0155    tmp   = lc[1] ;
0156    lc[1] = refl( lc[3] ) ;
0157    lc[3] = refl( tmp   ) ;
0158    tmp   = lc[4] ;
0159    lc[4] = refl( lc[6] ) ;
0160    lc[6] = refl( tmp   ) ;
0161    tmp   = lc[5] ;
0162    lc[5] = refl( lc[7] ) ;
0163    lc[7] = refl( tmp   ) ;
0164 */
0165   lc[0] = truncPyr::refl(lc[0]);
0166   lc[1] = truncPyr::refl(lc[1]);
0167   lc[2] = truncPyr::refl(lc[2]);
0168   lc[3] = truncPyr::refl(lc[3]);
0169   lc[4] = truncPyr::refl(lc[4]);
0170   lc[5] = truncPyr::refl(lc[5]);
0171   lc[6] = truncPyr::refl(lc[6]);
0172   lc[7] = truncPyr::refl(lc[7]);
0173
0174   ref = 0.25 * (lc[0] + lc[1] + lc[2] + lc[3]);
0175 }
0176
0177 void TruncatedPyramid::localCorners(Pt3DVec& lc, const CCGFloat* pv, Pt3D& ref) {
0178   assert(nullptr != pv);
0179   assert(8 == lc.size());
0180
0181   const CCGFloat dz(pv[TruncatedPyramid::k_Dz]);
0182   const CCGFloat th(pv[TruncatedPyramid::k_Theta]);
0183   const CCGFloat ph(pv[TruncatedPyramid::k_Phi]);
0184   const CCGFloat h1(pv[TruncatedPyramid::k_Dy1]);
0185   const CCGFloat b1(pv[TruncatedPyramid::k_Dx1]);
0186   const CCGFloat t1(pv[TruncatedPyramid::k_Dx2]);
0187   const CCGFloat a1(pv[TruncatedPyramid::k_Alp1]);
0188   const CCGFloat h2(pv[TruncatedPyramid::k_Dy2]);
0189   const CCGFloat b2(pv[TruncatedPyramid::k_Dx3]);
0190   const CCGFloat t2(pv[TruncatedPyramid::k_Dx4]);
0191   const CCGFloat a2(pv[TruncatedPyramid::k_Alp2]);
0192
0193   const CCGFloat ta1(tan(a1));  // lower plane
0194   const CCGFloat ta2(tan(a2));  // upper plane
0195
0196   const CCGFloat tth(tan(th));
0197   const CCGFloat tthcp(tth * cos(ph));
0198   const CCGFloat tthsp(tth * sin(ph));
0199
0200   const unsigned int off(h1 < h2 ? 0 : 4);
0201
0202   lc[0 + off] = Pt3D(-dz * tthcp - h1 * ta1 - b1, -dz * tthsp - h1, -dz);  // (-,-,-)
0203   lc[1 + off] = Pt3D(-dz * tthcp + h1 * ta1 - t1, -dz * tthsp + h1, -dz);  // (-,+,-)
0204   lc[2 + off] = Pt3D(-dz * tthcp + h1 * ta1 + t1, -dz * tthsp + h1, -dz);  // (+,+,-)
0205   lc[3 + off] = Pt3D(-dz * tthcp - h1 * ta1 + b1, -dz * tthsp - h1, -dz);  // (+,-,-)
0206   lc[4 - off] = Pt3D(dz * tthcp - h2 * ta2 - b2, dz * tthsp - h2, dz);     // (-,-,+)
0207   lc[5 - off] = Pt3D(dz * tthcp + h2 * ta2 - t2, dz * tthsp + h2, dz);     // (-,+,+)
0208   lc[6 - off] = Pt3D(dz * tthcp + h2 * ta2 + t2, dz * tthsp + h2, dz);     // (+,+,+)
0209   lc[7 - off] = Pt3D(dz * tthcp - h2 * ta2 + b2, dz * tthsp - h2, dz);     // (+,-,+)
0210
0211   ref = 0.25 * (lc[0] + lc[1] + lc[2] + lc[3]);
0212 }
0213
0214 void TruncatedPyramid::localCornersSwap(Pt3DVec& lc, const CCGFloat* pv, Pt3D& ref) {
0215   localCorners(lc, pv, ref);
0216
0217   Pt3D tmp;
0218   tmp = lc[0];
0219   lc[0] = lc[3];
0220   lc[3] = tmp;
0221   tmp = lc[1];
0222   lc[1] = lc[2];
0223   lc[2] = tmp;
0224   tmp = lc[4];
0225   lc[4] = lc[7];
0226   lc[7] = tmp;
0227   tmp = lc[5];
0228   lc[5] = lc[6];
0229   lc[6] = tmp;
0230
0231   ref = 0.25 * (lc[0] + lc[1] + lc[2] + lc[3]);
0232 }
0233
0234 // the following function is static and a helper for the endcap & barrel loader classes
0235 // when initializing from DDD: fills corners vector from trap params plus transform
0236
0237 void TruncatedPyramid::createCorners(const std::vector<CCGFloat>& pv, const Tr3D& tr, std::vector<GlobalPoint>& co) {
0238   assert(11 == pv.size());
0239   assert(8 == co.size());
0240   // to get the ordering right for fast sim, we have to use their convention
0241   // which were based on the old static geometry. Some gymnastics required here.
0242
0243   const CCGFloat dz(pv[0]);
0244   const CCGFloat h1(pv[3]);
0245   const CCGFloat h2(pv[7]);
0246   Pt3DVec ko(8, Pt3D(0, 0, 0));
0247
0248   // if reflection, different things for barrel and endcap
0249   static const FVec3D x(1, 0, 0);
0250   static const FVec3D y(0, 1, 0);
0251   static const FVec3D z(0, 0, 1);
0252   const bool refl(((tr * x).cross(tr * y)).dot(tr * z) < 0);  // has reflection!
0253
0254   Pt3D tmp;
0255   Pt3DVec to(8, Pt3D(0, 0, 0));
0256   localCorners(to, &pv.front(), tmp);
0257
0258   for (unsigned int i(0); i != 8; ++i) {
0259     ko[i] = tr * to[i];  // apply transformation
0260   }
0261
0262   if (refl || h1 > h2) {
0263     if (11.2 < dz)  //barrel
0264     {
0265       if (!refl) {
0266         to[0] = ko[3];
0267         to[1] = ko[2];
0268         to[2] = ko[1];
0269         to[3] = ko[0];
0270         to[4] = ko[7];
0271         to[5] = ko[6];
0272         to[6] = ko[5];
0273         to[7] = ko[4];
0274       } else {
0275         to[0] = ko[0];
0276         to[1] = ko[1];
0277         to[2] = ko[2];
0278         to[3] = ko[3];
0279         to[4] = ko[4];
0280         to[5] = ko[5];
0281         to[6] = ko[6];
0282         to[7] = ko[7];
0283       }
0284     } else  //endcap
0285     {
0286       to[0] = ko[0];
0287       to[1] = ko[3];
0288       to[2] = ko[2];
0289       to[3] = ko[1];
0290       to[4] = ko[4];
0291       to[5] = ko[7];
0292       to[6] = ko[6];
0293       to[7] = ko[5];
0294     }
0295     copy(to.begin(), to.end(), ko.begin());  // faster than ko = to ? maybe.
0296   }
0297   for (unsigned int i(0); i != 8; ++i) {
0298     const Pt3D& p(ko[i]);
0299     co[i] = GlobalPoint(p.x(), p.y(), p.z());
0300   }
0301 }
0302 //----------------------------------------------------------------------
0303 //----------------------------------------------------------------------
0304
0305 std::ostream& operator<<(std::ostream& s, const TruncatedPyramid& cell) {
0306   s << "Center: " << ((const CaloCellGeometry&)cell).getPosition() << std::endl;
0307   const float thetaaxis(cell.getThetaAxis());
0308   const float phiaxis(cell.getPhiAxis());
0309   s << "Axis: " << thetaaxis << " " << phiaxis << std::endl;
0310   const CaloCellGeometry::CornersVec& corners(cell.getCorners());
0311   for (unsigned int i = 0; i != corners.size(); ++i) {
0312     s << "Corner: " << corners[i] << std::endl;
0313   }
0314   return s;
0315 }