Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2023-03-17 11:24:15

0001 // to make hits in EB/EE/HC
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 //#define EDM_ML_DEBUG
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 }