File indexing completed on 2024-04-06 12:29:53
0001
0002 #include "SimG4CMS/Calo/interface/HcalNumberingFromPS.h"
0003 #include "DataFormats/Math/interface/GeantUnits.h"
0004 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0005
0006 #include <cmath>
0007 #include <iostream>
0008 #include <iomanip>
0009
0010
0011 using namespace geant_units::operators;
0012
0013 HcalNumberingFromPS::HcalNumberingFromPS(const edm::ParameterSet& conf) {
0014 etaTable_ = conf.getParameter<std::vector<double> >("EtaTable");
0015 phibin_ = conf.getParameter<std::vector<double> >("PhiBin");
0016 phioff_ = conf.getParameter<std::vector<double> >("PhiOffset");
0017 etaMin_ = conf.getParameter<std::vector<int> >("EtaMin");
0018 etaMax_ = conf.getParameter<std::vector<int> >("EtaMax");
0019 etaHBHE_ = conf.getParameter<int>("EtaHBHE");
0020 depthHBHE_ = conf.getParameter<std::vector<int> >("DepthHBHE");
0021 depth29Mx_ = conf.getParameter<int>("Depth29Max");
0022 rMinHO_ = conf.getParameter<double>("RMinHO");
0023 zHO_ = conf.getParameter<std::vector<double> >("ZHO");
0024 const double deg = M_PI / 180.0;
0025 for (auto& phi : phibin_)
0026 phi *= deg;
0027 for (auto& phi : phioff_)
0028 phi *= deg;
0029 #ifdef EDM_ML_DEBUG
0030 edm::LogVerbatim("HcalSim") << "HcalNumberingFromPS:: EtaTable with " << etaTable_.size() << " elements";
0031 for (unsigned k = 0; k < etaTable_.size(); ++k)
0032 edm::LogVerbatim("HcalSim") << "EtaTable[" << k << "] = " << etaTable_[k];
0033 edm::LogVerbatim("HcalSim") << "HcalNumberingFromPS:: PhiBin with " << phibin_.size() << " elements";
0034 for (unsigned k = 0; k < phibin_.size(); ++k)
0035 edm::LogVerbatim("HcalSim") << "PhiBin[" << k << "] = " << phibin_[k];
0036 edm::LogVerbatim("HcalSim") << "HcalNumberingFromPS:: PhiOff with " << phioff_.size() << " elements";
0037 for (unsigned k = 0; k < phioff_.size(); ++k)
0038 edm::LogVerbatim("HcalSim") << "PhiOff[" << k << "] = " << phioff_[k];
0039 edm::LogVerbatim("HcalSim") << "HcalNumberingFromPS:: EtaMin/EtaMax with " << etaMin_.size() << ":" << etaMax_.size()
0040 << " elements";
0041 for (unsigned k = 0; k < etaMin_.size(); ++k)
0042 edm::LogVerbatim("HcalSim") << "EtaMin[" << k << "] = " << etaMin_[k] << " EtaMax[" << k << "] = " << etaMax_[k];
0043 edm::LogVerbatim("HcalSim") << "HcalNumberingFromPS:: EtaHBHE " << etaHBHE_ << " DepthHBHE " << depthHBHE_[0] << ":"
0044 << depthHBHE_[1] << " RMinHO " << rMinHO_ << " zHO with " << zHO_.size() << " elements";
0045 for (unsigned k = 0; k < zHO_.size(); ++k)
0046 edm::LogVerbatim("HcalSim") << "ZHO[" << k << "] = " << zHO_[k];
0047 #endif
0048
0049 segmentation_.resize(nEtas_);
0050 for (int ring = 0; ring < nEtas_; ++ring) {
0051 char name[10];
0052 snprintf(name, 10, "Eta%d", ring + 1);
0053 if (ring > 0) {
0054 segmentation_[ring] = conf.getUntrackedParameter<std::vector<int> >(name, segmentation_[ring - 1]);
0055 } else {
0056 segmentation_[ring] = conf.getUntrackedParameter<std::vector<int> >(name);
0057 }
0058 #ifdef EDM_ML_DEBUG
0059 edm::LogVerbatim("HcalSim") << "HcalNumberingFromPS:: Ring " << ring + 1 << " with " << segmentation_[ring].size()
0060 << " layers";
0061 for (unsigned int k = 0; k < segmentation_[ring].size(); ++k)
0062 edm::LogVerbatim("HcalSim") << "Layer[" << k << "] = " << segmentation_[ring][k];
0063 #endif
0064 }
0065 }
0066
0067 HcalNumberingFromDDD::HcalID HcalNumberingFromPS::unitID(int det,
0068 int layer,
0069 int depth,
0070 const math::XYZVectorD& pos) const {
0071 int subdet = ((det == 3) ? static_cast<int>(HcalBarrel) : static_cast<int>(HcalEndcap));
0072 std::pair<int, int> deteta = getEta(subdet, pos);
0073 std::pair<int, int> iphi = getPhi(deteta.first, deteta.second, pos.Phi());
0074 int newDepth(depth);
0075 int zside = ((pos.z() > 0) ? 1 : 0);
0076 if (deteta.first == static_cast<int>(HcalBarrel)) {
0077 newDepth = segmentation_[deteta.second - 1][layer - 1];
0078 if ((deteta.second == etaHBHE_) && (newDepth > depthHBHE_[0]))
0079 newDepth = depthHBHE_[0];
0080 } else if (deteta.first == static_cast<int>(HcalEndcap)) {
0081 newDepth = segmentation_[deteta.second - 1][layer - 1];
0082 if ((deteta.second == etaHBHE_) && (newDepth < depthHBHE_[1]))
0083 newDepth = depthHBHE_[1];
0084 if ((deteta.second == etaMax_[1]) && (newDepth > depth29Mx_))
0085 --(deteta.second);
0086 } else if (deteta.first == static_cast<int>(HcalOuter)) {
0087 newDepth = 4;
0088 }
0089 #ifdef EDM_ML_DEBUG
0090 edm::LogVerbatim("HcalSim") << "HcalNumberingFromPS:: det " << det << ":" << subdet << ":" << deteta.first
0091 << "\t Eta " << pos.Eta() << ":" << deteta.second << "\t Phi " << pos.Phi() << ":"
0092 << iphi.first << ":" << iphi.second << "\t Layer|Depth " << layer << ":" << depth << ":"
0093 << newDepth;
0094 #endif
0095 return HcalNumberingFromDDD::HcalID(deteta.first, zside, newDepth, deteta.second, iphi.first, iphi.second, layer);
0096 }
0097
0098 std::pair<int, int> HcalNumberingFromPS::getEta(const int& det, const math::XYZVectorD& pos) const {
0099 int ieta(1);
0100 int subdet(det);
0101 double eta = std::abs(pos.Eta());
0102 if (pos.Rho() > rMinHO_) {
0103 subdet = static_cast<int>(HcalOuter);
0104 double z = std::abs(pos.z());
0105 if (z > zHO_[3]) {
0106 if (eta <= etaTable_[10])
0107 eta = etaTable_[10] + 0.001;
0108 } else if (z > zHO_[1]) {
0109 if (eta <= etaTable_[4])
0110 eta = etaTable_[4] + 0.001;
0111 }
0112 }
0113 for (unsigned int i = 1; i < etaTable_.size(); i++) {
0114 if (eta < etaTable_[i]) {
0115 ieta = i;
0116 break;
0117 }
0118 }
0119
0120 if ((subdet == static_cast<int>(HcalBarrel)) || (subdet == static_cast<int>(HcalOuter))) {
0121 if (ieta > etaMax_[0])
0122 ieta = etaMax_[0];
0123 } else if (det == static_cast<int>(HcalEndcap)) {
0124 if (ieta <= etaMin_[1])
0125 ieta = etaMin_[1];
0126 }
0127 return std::make_pair(subdet, ieta);
0128 }
0129
0130 std::pair<int, int> HcalNumberingFromPS::getPhi(const int& det, const int& ieta, const double& phi) const {
0131 double fioff = ((det == static_cast<int>(HcalEndcap)) ? phioff_[1] : phioff_[0]);
0132 double fibin = phibin_[ieta - 1];
0133 int nphi = int((2._pi + 0.1 * fibin) / fibin);
0134 double hphi = phi + fioff;
0135 if (hphi < 0)
0136 hphi += (2._pi);
0137 int iphi = int(hphi / fibin) + 1;
0138 if (iphi > nphi)
0139 iphi = 1;
0140 const double fiveDegInRad = 5._deg;
0141 int units = int(fibin / fiveDegInRad + 0.5);
0142 if (units < 1)
0143 units = 1;
0144 int iphi_skip = iphi;
0145 if (units == 2)
0146 iphi_skip = (iphi - 1) * 2 + 1;
0147 else if (units == 4)
0148 iphi_skip = (iphi - 1) * 4 - 1;
0149 if (iphi_skip < 0)
0150 iphi_skip += 72;
0151 return std::make_pair(iphi, iphi_skip);
0152 }
