File indexing completed on 2024-04-06 12:14:17
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
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)
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))
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 }
0146
0147 void TruncatedPyramid::localCornersReflection(Pt3DVec& lc, const CCGFloat* pv, Pt3D& ref) {
0148
0149 localCorners(lc, pv, ref);
0150 Pt3D tmp;
0151
0152
0153
0154
0155
0156
0157
0158
0159
0160
0161
0162
0163
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));
0194 const CCGFloat ta2(tan(a2));
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
0235
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
0241
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
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);
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];
0260 }
0261
0262 if (refl || h1 > h2) {
0263 if (11.2 < dz)
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
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());
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 }