File indexing completed on 2024-04-06 12:14:17
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
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];
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);
0163 lc[1] = Pt3D(-r, -R, -dz);
0164 lc[2] = Pt3D(-r, R, -dz);
0165 lc[3] = Pt3D(0, 2 * R, -dz);
0166 lc[4] = Pt3D(r, R, -dz);
0167 lc[5] = Pt3D(r, -R, -dz);
0168 lc[6] = Pt3D(0, -2 * R, dz);
0169 lc[7] = Pt3D(-r, -R, dz);
0170 lc[8] = Pt3D(-r, R, dz);
0171 lc[9] = Pt3D(0, 2 * R, dz);
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
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)) {
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 }
