HcalCommonData/src/HcalNumberingFromDDD.cc

0001 ///////////////////////////////////////////////////////////////////////////////
0002 // File: HcalNumberingFromDDD.cc
0003 // Description: Usage of DDD to get to numbering scheme for hadron calorimeter
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 //#define EDM_ML_DEBUG
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   //First eta index
0047   if (det == 5) {  // Forward HCal
0048     hsubdet = static_cast<int>(HcalForward);
0049     hR = sqrt(hx * hx + hy * hy);
0050     etaR = (heta >= 0. ? hR : -hR);
0051   } else {  // Barrel or Endcap
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 }