HcalNumberingFromDDD.cc

Go to the documentation of this file.

// File: HcalNumberingFromDDD.cc
// Description: Usage of DDD to get to numbering scheme for hadron calorimeter

#include "CondFormats/GeometryObjects/interface/HcalParameters.h"
#include "Geometry/HcalCommonData/interface/HcalNumberingFromDDD.h"
#include "DataFormats/Math/interface/GeantUnits.h"

#include "FWCore/MessageLogger/interface/MessageLogger.h"

#include <iostream>

//#define EDM_ML_DEBUG
using namespace geant_units;
using namespace geant_units::operators;

HcalNumberingFromDDD::HcalNumberingFromDDD(const HcalDDDSimConstants *hcons) :
  hcalConstants(hcons) {
#ifdef EDM_ML_DEBUG
  edm::LogInfo("HCalGeom") << "Creating HcalNumberingFromDDD\n";
#endif
}

HcalNumberingFromDDD::~HcalNumberingFromDDD() {
#ifdef EDM_ML_DEBUG
  edm::LogInfo("HCalGeom") << "Deleting HcalNumberingFromDDD\n";
#endif
}

HcalNumberingFromDDD::HcalID HcalNumberingFromDDD::unitID(int det,
                              const math::XYZVectorD& point,
                              int depth,
                              int lay) const {


  double hx  = point.x();
  double hy  = point.y();
  double hz  = point.z();
  double hR  = sqrt(hx*hx+hy*hy+hz*hz);
  double htheta = (hR == 0. ? 0. : acos(std::max(std::min(hz/hR,1.0),-1.0)));
  double hsintheta = sin(htheta);
  double hphi = (hR*hsintheta == 0. ? 0. :atan2(hy,hx));
  double heta = (fabs(hsintheta) == 1.? 0. : -log(fabs(tan(htheta/2.))) );

  int    hsubdet=0;
  double etaR;

  //First eta index
  if (det == 5) { // Forward HCal
    hsubdet = static_cast<int>(HcalForward);
    hR      = sqrt(hx*hx+hy*hy);
    etaR    = (heta >= 0. ? hR : -hR);
  } else { // Barrel or Endcap
    etaR    = heta;
    if (det == 3) {
      hsubdet = static_cast<int>(HcalBarrel);
      etaR    = hcalConstants->getEtaHO(heta,hx,hy,hz);
    } else {
      hsubdet = static_cast<int>(HcalEndcap);
    }
  }

#ifdef EDM_ML_DEBUG
  edm::LogVerbatim("HCalGeom") << "HcalNumberingFromDDD: point = " << point
                   << " det " << det << ":" << hsubdet << " eta/R "
                   << etaR << " phi " << hphi << " depth " << depth
                   << " layer " << lay;
#endif
  return unitID(hsubdet,etaR,hphi,depth,lay);
}

HcalNumberingFromDDD::HcalID HcalNumberingFromDDD::unitID(double eta,double fi,
                              int depth, 
                              int lay) const {
  std::pair<int,double> detEta = hcalConstants->getDetEta(eta, depth);
  return unitID(detEta.first,detEta.second,fi,depth,lay);
}


HcalNumberingFromDDD::HcalID HcalNumberingFromDDD::unitID(int det,
                              double etaR,
                              double phi,
                              int depth,
                              int lay) const {

  double hetaR = fabs(etaR);
  int    ieta  = hcalConstants->getEta(det, lay, hetaR);
  std::pair<double,double> ficons = hcalConstants->getPhiCons(det, ieta);

  int    nphi  = int((2._pi+0.1*ficons.second)/ficons.second);
  int    zside = etaR>0 ? 1: 0;
  double hphi  = phi+ficons.first;
  if (hphi < 0)    hphi += (2._pi);
  int    iphi  = int(hphi/ficons.second) + 1;
  if (iphi > nphi) iphi = 1;

#ifdef EDM_ML_DEBUG
  edm::LogVerbatim("HCalGeom") << "HcalNumberingFromDDD: etaR = " << etaR 
                   << " : " << zside << "/" << ieta << " phi "
                   << hphi << " : " << iphi << std::endl;
#endif
  return unitID(det,zside,depth,ieta,iphi,lay);
}

HcalNumberingFromDDD::HcalID HcalNumberingFromDDD::unitID(int det, int zside,
                              int depth, int etaR,
                              int phi, 
                              int lay) const {


  if (det == static_cast<int>(HcalBarrel) && lay > 17)
    det = static_cast<int>(HcalOuter);

  int units     = hcalConstants->unitPhi(det, etaR);
  int iphi_skip = hcalConstants->phiNumber(phi, units);

  std::pair<int,int> etaDepth = hcalConstants->getEtaDepth(det, etaR, iphi_skip, zside, depth, lay);

#ifdef EDM_ML_DEBUG
  edm::LogVerbatim("HCalGeom") << "HcalNumberingFromDDD: phi units= " << units
                   << "  iphi_skip= " << iphi_skip; 
#endif
  HcalNumberingFromDDD::HcalID tmp(det,zside,etaDepth.second,etaDepth.first,phi,iphi_skip,lay);

#ifdef EDM_ML_DEBUG
  edm::LogVerbatim("HCalGeom") << "HcalNumberingFromDDD: det = " << det << " " 
                   << tmp.subdet << " zside = " << tmp.zside 
                   << " depth = " << tmp.depth << " eta/R = " 
                   << tmp.etaR << " phi = "   << tmp.phi 
                   << " layer = " << tmp.lay;
#endif
  return tmp;
}