File indexing completed on 2024-04-06 12:14:47
0001
0002
0003
0004
0005
0006 #include "CondFormats/GeometryObjects/interface/HcalParameters.h"
0007 #include "Geometry/HcalCommonData/interface/HcalNumberingFromDDD.h"
0008 #include "DataFormats/Math/interface/GeantUnits.h"
0009
0010 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0011
0012 #include <iostream>
0013
0014
0015 using namespace geant_units;
0016 using namespace geant_units::operators;
0017
0018 HcalNumberingFromDDD::HcalNumberingFromDDD(const HcalDDDSimConstants* hcons) : hcalConstants(hcons) {
0019 #ifdef EDM_ML_DEBUG
0020 edm::LogVerbatim("HCalGeom") << "Creating HcalNumberingFromDDD";
0021 #endif
0022 }
0023
0024 HcalNumberingFromDDD::~HcalNumberingFromDDD() {
0025 #ifdef EDM_ML_DEBUG
0026 edm::LogVerbatim("HCalGeom") << "Deleting HcalNumberingFromDDD";
0027 #endif
0028 }
0029
0030 HcalNumberingFromDDD::HcalID HcalNumberingFromDDD::unitID(int det,
0031 const math::XYZVectorD& point,
0032 int depth,
0033 int lay) const {
0034 double hx = point.x();
0035 double hy = point.y();
0036 double hz = point.z();
0037 double hR = sqrt(hx * hx + hy * hy + hz * hz);
0038 double htheta = (hR == 0. ? 0. : acos(std::max(std::min(hz / hR, 1.0), -1.0)));
0039 double hsintheta = sin(htheta);
0040 double hphi = (hR * hsintheta == 0. ? 0. : atan2(hy, hx));
0041 double heta = (fabs(hsintheta) == 1. ? 0. : -log(fabs(tan(htheta / 2.))));
0042
0043 int hsubdet = 0;
0044 double etaR;
0045
0046
0047 if (det == 5) {
0048 hsubdet = static_cast<int>(HcalForward);
0049 hR = sqrt(hx * hx + hy * hy);
0050 etaR = (heta >= 0. ? hR : -hR);
0051 } else {
0052 etaR = heta;
0053 if (det == 3) {
0054 hsubdet = static_cast<int>(HcalBarrel);
0055 etaR = hcalConstants->getEtaHO(heta, hx, hy, hz);
0056 } else {
0057 hsubdet = static_cast<int>(HcalEndcap);
0058 }
0059 }
0060
0061 #ifdef EDM_ML_DEBUG
0062 edm::LogVerbatim("HCalGeom") << "HcalNumberingFromDDD: point = " << point << " det " << det << ":" << hsubdet
0063 << " eta/R " << etaR << " phi " << hphi << " depth " << depth << " layer " << lay;
0064 #endif
0065 return unitID(hsubdet, etaR, hphi, depth, lay);
0066 }
0067
0068 HcalNumberingFromDDD::HcalID HcalNumberingFromDDD::unitID(double eta, double fi, int depth, int lay) const {
0069 std::pair<int, double> detEta = hcalConstants->getDetEta(eta, depth);
0070 return unitID(detEta.first, detEta.second, fi, depth, lay);
0071 }
0072
0073 HcalNumberingFromDDD::HcalID HcalNumberingFromDDD::unitID(int det, double etaR, double phi, int depth, int lay) const {
0074 double hetaR = fabs(etaR);
0075 int ieta = hcalConstants->getEta(det, lay, hetaR);
0076 std::pair<double, double> ficons = hcalConstants->getPhiCons(det, ieta);
0077
0078 int nphi = int((2._pi + 0.1 * ficons.second) / ficons.second);
0079 int zside = etaR > 0 ? 1 : 0;
0080 double hphi = phi + ficons.first;
0081 if (hphi < 0)
0082 hphi += (2._pi);
0083 int iphi = int(hphi / ficons.second) + 1;
0084 if (iphi > nphi)
0085 iphi = 1;
0086
0087 #ifdef EDM_ML_DEBUG
0088 edm::LogVerbatim("HCalGeom") << "HcalNumberingFromDDD: etaR = " << etaR << " : " << zside << "/" << ieta << " phi "
0089 << hphi << " : " << iphi;
0090 #endif
0091 return unitID(det, zside, depth, ieta, iphi, lay);
0092 }
0093
0094 HcalNumberingFromDDD::HcalID HcalNumberingFromDDD::unitID(
0095 int det, int zside, int depth, int etaR, int phi, int lay) const {
0096 if (det == static_cast<int>(HcalBarrel) && lay > 17)
0097 det = static_cast<int>(HcalOuter);
0098
0099 int units = hcalConstants->unitPhi(det, etaR);
0100 int iphi_skip = hcalConstants->phiNumber(phi, units);
0101
0102 if ((lay == 1) && (etaR == 16))
0103 etaR = 15;
0104
0105 std::pair<int, int> etaDepth = hcalConstants->getEtaDepth(det, etaR, iphi_skip, zside, depth, lay);
0106
0107 #ifdef EDM_ML_DEBUG
0108 edm::LogVerbatim("HCalGeom") << "HcalNumberingFromDDD: phi units= " << units << " iphi_skip= " << iphi_skip;
0109 #endif
0110 HcalNumberingFromDDD::HcalID tmp(det, zside, etaDepth.second, etaDepth.first, phi, iphi_skip, lay);
0111
0112 #ifdef EDM_ML_DEBUG
0113 edm::LogVerbatim("HCalGeom") << "HcalNumberingFromDDD: det = " << det << " " << tmp.subdet << " zside = " << tmp.zside
0114 << " depth = " << tmp.depth << " eta/R = " << tmp.etaR << " phi = " << tmp.phi
0115 << " layer = " << tmp.lay;
0116 #endif
0117 return tmp;
0118 }