Skip to content

Commit

Permalink
Backport the bugfix to HGCal geometry from the PR's 37330 and 37460
Browse files Browse the repository at this point in the history
  • Loading branch information
Sunanda committed Apr 5, 2022
1 parent 106580e commit 250803c
Show file tree
Hide file tree
Showing 4 changed files with 70 additions and 54 deletions.
11 changes: 9 additions & 2 deletions Geometry/HGCalCommonData/interface/HGCalDDDConstants.h
Expand Up @@ -84,8 +84,15 @@ class HGCalDDDConstants {
std::pair<float, float> localToGlobal8(
int lay, int waferU, int waferV, double localX, double localY, bool reco, bool debug) const;
std::pair<float, float> locateCell(int cell, int lay, int type, bool reco) const;
std::pair<float, float> locateCell(
int lay, int waferU, int waferV, int cellU, int cellV, bool reco, bool all, bool debug = false) const;
std::pair<float, float> locateCell(int lay,
int waferU,
int waferV,
int cellU,
int cellV,
bool reco,
bool all,
bool norot = false,
bool debug = false) const;
std::pair<float, float> locateCell(const HGCSiliconDetId&, bool debug = false) const;
std::pair<float, float> locateCell(const HGCScintillatorDetId&, bool debug = false) const;
std::pair<float, float> locateCellHex(int cell, int wafer, bool reco) const;
Expand Down
24 changes: 9 additions & 15 deletions Geometry/HGCalCommonData/src/HGCalDDDConstants.cc
Expand Up @@ -10,6 +10,7 @@
#include "Geometry/HGCalCommonData/interface/HGCalGeomParameters.h"
#include "Geometry/HGCalCommonData/interface/HGCalGeomTools.h"
#include "Geometry/HGCalCommonData/interface/HGCalGeometryMode.h"
#include "Geometry/HGCalCommonData/interface/HGCalCell.h"
#include "Geometry/HGCalCommonData/interface/HGCalTypes.h"
#include "Geometry/HGCalCommonData/interface/HGCalWaferIndex.h"
#include "Geometry/HGCalCommonData/interface/HGCalWaferMask.h"
Expand Down Expand Up @@ -213,7 +214,7 @@ bool HGCalDDDConstants::cellInLayer(int waferU, int waferV, int cellU, int cellV
const auto& indx = getIndex(lay, true);
if (indx.first >= 0) {
if (waferHexagon8() || waferHexagon6()) {
const auto& xy = ((waferHexagon8()) ? locateCell(lay, waferU, waferV, cellU, cellV, reco, true, true)
const auto& xy = ((waferHexagon8()) ? locateCell(lay, waferU, waferV, cellU, cellV, reco, true, false, false)
: locateCell(cellU, lay, waferU, reco));
double rpos = sqrt(xy.first * xy.first + xy.second * xy.second);
return ((rpos >= hgpar_->rMinLayHex_[indx.first]) && (rpos <= hgpar_->rMaxLayHex_[indx.first]));
Expand Down Expand Up @@ -621,7 +622,7 @@ std::pair<float, float> HGCalDDDConstants::localToGlobal8(
bool rotx =
((!hgpar_->layerType_.empty()) && (hgpar_->layerType_[lay - hgpar_->firstLayer_] == HGCalTypes::WaferCenterR));
if (debug)
edm::LogVerbatim("HGCalGeom") << "LocalToGlobal " << lay << ":" << (lay - hgpar_->firstLayer_) << ":" << rotx
edm::LogVerbatim("HGCalGeom") << "LocalToGlobal8 " << lay << ":" << (lay - hgpar_->firstLayer_) << ":" << rotx
<< " Local (" << x << ":" << y << ") Reco " << reco;
if (!reco) {
x *= HGCalParameters::k_ScaleToDDD;
Expand All @@ -631,7 +632,7 @@ std::pair<float, float> HGCalDDDConstants::localToGlobal8(
x += xy.first;
y += xy.second;
if (debug)
edm::LogVerbatim("HGCalGeom") << "With wafer " << x << ":" << y << " by addong " << xy.first << ":" << xy.second;
edm::LogVerbatim("HGCalGeom") << "With wafer " << x << ":" << y << " by adding " << xy.first << ":" << xy.second;
return (rotx ? getXY(lay, x, y, false) : std::make_pair(x, y));
}

Expand Down Expand Up @@ -667,44 +668,39 @@ std::pair<float, float> HGCalDDDConstants::locateCell(int cell, int lay, int typ
}

std::pair<float, float> HGCalDDDConstants::locateCell(
int lay, int waferU, int waferV, int cellU, int cellV, bool reco, bool all, bool debug) const {
int lay, int waferU, int waferV, int cellU, int cellV, bool reco, bool all, bool norot, bool debug) const {
double x(0), y(0);
int indx = HGCalWaferIndex::waferIndex(lay, waferU, waferV);
auto itr = hgpar_->typesInLayers_.find(indx);
int type = ((itr == hgpar_->typesInLayers_.end()) ? 2 : hgpar_->waferTypeL_[itr->second]);
bool rotx =
((!hgpar_->layerType_.empty()) && (hgpar_->layerType_[lay - hgpar_->firstLayer_] == HGCalTypes::WaferCenterR));
#ifdef EDM_ML_DEBUG
bool rotx = (norot) ? false
: ((!hgpar_->layerType_.empty()) &&
(hgpar_->layerType_[lay - hgpar_->firstLayer_] == HGCalTypes::WaferCenterR));
if (debug) {
edm::LogVerbatim("HGCalGeom") << "LocateCell " << lay << ":" << (lay - hgpar_->firstLayer_) << ":" << rotx << ":"
<< waferU << ":" << waferV << ":" << indx << ":"
<< (itr == hgpar_->typesInLayers_.end()) << ":" << type << " Flags " << reco << ":"
<< all;
}
#endif
int kndx = cellV * 100 + cellU;
if (type == 0) {
auto ktr = hgpar_->cellFineIndex_.find(kndx);
if (ktr != hgpar_->cellFineIndex_.end()) {
x = hgpar_->cellFineX_[ktr->second];
y = hgpar_->cellFineY_[ktr->second];
}
#ifdef EDM_ML_DEBUG
if (debug)
edm::LogVerbatim("HGCalGeom") << "Fine " << cellU << ":" << cellV << ":" << kndx << ":" << x << ":" << y << ":"
<< (ktr != hgpar_->cellFineIndex_.end());
#endif
} else {
auto ktr = hgpar_->cellCoarseIndex_.find(kndx);
if (ktr != hgpar_->cellCoarseIndex_.end()) {
x = hgpar_->cellCoarseX_[ktr->second];
y = hgpar_->cellCoarseY_[ktr->second];
}
#ifdef EDM_ML_DEBUG
if (debug)
edm::LogVerbatim("HGCalGeom") << "Coarse " << cellU << ":" << cellV << ":" << kndx << ":" << x << ":" << y << ":"
<< (ktr != hgpar_->cellCoarseIndex_.end());
#endif
}
if (!reco) {
x *= HGCalParameters::k_ScaleToDDD;
Expand All @@ -714,10 +710,8 @@ std::pair<float, float> HGCalDDDConstants::locateCell(
const auto& xy = waferPositionNoRot(lay, waferU, waferV, reco, debug);
x += xy.first;
y += xy.second;
#ifdef EDM_ML_DEBUG
if (debug)
edm::LogVerbatim("HGCalGeom") << "With wafer " << x << ":" << y << " by addong " << xy.first << ":" << xy.second;
#endif
edm::LogVerbatim("HGCalGeom") << "With wafer " << x << ":" << y << " by adding " << xy.first << ":" << xy.second;
}
return (rotx ? getXY(lay, x, y, false) : std::make_pair(x, y));
}
Expand Down
6 changes: 3 additions & 3 deletions Geometry/HGCalGeometry/interface/HGCalGeometry.h
Expand Up @@ -50,7 +50,7 @@ class HGCalGeometry final : public CaloSubdetectorGeometry {

HGCalGeometry(const HGCalTopology& topology);

~HGCalGeometry() override;
~HGCalGeometry() override = default;

void localCorners(Pt3DVec& lc, const CCGFloat* pv, unsigned int i, Pt3D& ref);

Expand All @@ -70,7 +70,7 @@ class HGCalGeometry final : public CaloSubdetectorGeometry {
CaloSubdetectorGeometry::DimVec& dimVector,
CaloSubdetectorGeometry::IVec& dinsVector) const override;

GlobalPoint getPosition(const DetId& id) const;
GlobalPoint getPosition(const DetId& id, bool debug = false) const;
GlobalPoint getWaferPosition(const DetId& id) const;

/// Returns area of a cell
Expand All @@ -79,7 +79,7 @@ class HGCalGeometry final : public CaloSubdetectorGeometry {
/// Returns the corner points of this cell's volume.
CornersVec getCorners(const DetId& id) const;
CornersVec get8Corners(const DetId& id) const;
CornersVec getNewCorners(const DetId& id) const;
CornersVec getNewCorners(const DetId& id, bool debug = false) const;

// Get neighbor in z along a direction
DetId neighborZ(const DetId& idin, const GlobalVector& p) const;
Expand Down
83 changes: 49 additions & 34 deletions Geometry/HGCalGeometry/src/HGCalGeometry.cc
Expand Up @@ -23,6 +23,8 @@ typedef std::vector<float> ParmVec;

//#define EDM_ML_DEBUG

const bool debugLocate = false;

HGCalGeometry::HGCalGeometry(const HGCalTopology& topology_)
: m_topology(topology_),
m_validGeomIds(topology_.totalGeomModules()),
Expand All @@ -40,8 +42,6 @@ HGCalGeometry::HGCalGeometry(const HGCalTopology& topology_)
#endif
}

HGCalGeometry::~HGCalGeometry() {}

void HGCalGeometry::fillNamedParams(DDFilteredView fv) {}

void HGCalGeometry::initializeParms() {}
Expand Down Expand Up @@ -183,7 +183,7 @@ std::shared_ptr<const CaloCellGeometry> HGCalGeometry::getGeometry(const DetId&
return nullptr; // nothing to get
DetId geomId = getGeometryDetId(detId);
const uint32_t cellIndex(m_topology.detId2denseGeomId(geomId));
const GlobalPoint pos = (detId != geomId) ? getPosition(detId) : GlobalPoint();
const GlobalPoint pos = (detId != geomId) ? getPosition(detId, false) : GlobalPoint();
return cellGeomPtr(cellIndex, pos);
}

Expand All @@ -195,7 +195,7 @@ bool HGCalGeometry::present(const DetId& detId) const {
return (nullptr != getGeometryRawPtr(index));
}

GlobalPoint HGCalGeometry::getPosition(const DetId& detid) const {
GlobalPoint HGCalGeometry::getPosition(const DetId& detid, bool debug) const {
unsigned int cellIndex = indexFor(detid);
GlobalPoint glob;
unsigned int maxSize = (m_topology.tileTrapezoid() ? m_cellVec2.size() : m_cellVec.size());
Expand All @@ -206,32 +206,30 @@ GlobalPoint HGCalGeometry::getPosition(const DetId& detid) const {
xy = m_topology.dddConstants().locateCellHex(id.iCell1, id.iSec1, true);
const HepGeom::Point3D<float> lcoord(xy.first, xy.second, 0);
glob = m_cellVec[cellIndex].getPosition(lcoord);
#ifdef EDM_ML_DEBUG
edm::LogVerbatim("HGCalGeom") << "getPosition:: index " << cellIndex << " Local " << lcoord.x() << ":"
<< lcoord.y() << " ID " << id.iCell1 << ":" << id.iSec1 << " Global " << glob;
#endif
if (debug)
edm::LogVerbatim("HGCalGeom") << "getPosition:: index " << cellIndex << " Local " << lcoord.x() << ":"
<< lcoord.y() << " ID " << id.iCell1 << ":" << id.iSec1 << " Global " << glob;
} else if (m_topology.tileTrapezoid()) {
const HepGeom::Point3D<float> lcoord(0, 0, 0);
glob = m_cellVec2[cellIndex].getPosition(lcoord);
#ifdef EDM_ML_DEBUG
edm::LogVerbatim("HGCalGeom") << "getPositionTrap:: index " << cellIndex << " Local " << lcoord.x() << ":"
<< lcoord.y() << " ID " << id.iLay << ":" << id.iSec1 << ":" << id.iCell1
<< " Global " << glob;
#endif
if (debug)
edm::LogVerbatim("HGCalGeom") << "getPositionTrap:: index " << cellIndex << " Local " << lcoord.x() << ":"
<< lcoord.y() << " ID " << id.iLay << ":" << id.iSec1 << ":" << id.iCell1
<< " Global " << glob;
} else {
#ifdef EDM_ML_DEBUG
edm::LogVerbatim("HGCalGeom") << "getPosition for " << HGCSiliconDetId(detid) << " Layer " << id.iLay << " Wafer "
<< id.iSec1 << ":" << id.iSec2 << " Cell " << id.iCell1 << ":" << id.iCell2;
#endif
xy = m_topology.dddConstants().locateCell(id.iLay, id.iSec1, id.iSec2, id.iCell1, id.iCell2, true, true, true);
if (debug)
edm::LogVerbatim("HGCalGeom") << "getPosition for " << HGCSiliconDetId(detid) << " Layer " << id.iLay
<< " Wafer " << id.iSec1 << ":" << id.iSec2 << " Cell " << id.iCell1 << ":"
<< id.iCell2;
xy = m_topology.dddConstants().locateCell(
id.iLay, id.iSec1, id.iSec2, id.iCell1, id.iCell2, true, true, false, debug);
double xx = id.zSide * xy.first;
double zz = id.zSide * m_topology.dddConstants().waferZ(id.iLay, true);
glob = GlobalPoint(xx, xy.second, zz);
#ifdef EDM_ML_DEBUG
edm::LogVerbatim("HGCalGeom") << "getPositionWafer:: index " << cellIndex << " Local " << xy.first << ":"
<< xy.second << " ID " << id.iLay << ":" << id.iSec1 << ":" << id.iSec2 << ":"
<< id.iCell1 << ":" << id.iCell2 << " Global " << glob;
#endif
if (debug)
edm::LogVerbatim("HGCalGeom") << "getPositionWafer:: index " << cellIndex << " Local " << xy.first << ":"
<< xy.second << " ID " << id.iLay << ":" << id.iSec1 << ":" << id.iSec2 << ":"
<< id.iCell1 << ":" << id.iCell2 << " Global " << glob;
}
}
return glob;
Expand Down Expand Up @@ -276,7 +274,7 @@ HGCalGeometry::CornersVec HGCalGeometry::getCorners(const DetId& detid) const {
unsigned int cellIndex = indexFor(detid);
HGCalTopology::DecodedDetId id = m_topology.decode(detid);
if (cellIndex < m_cellVec2.size() && m_det == DetId::HGCalHSc) {
GlobalPoint v = getPosition(detid);
GlobalPoint v = getPosition(detid, false);
int type = std::min(id.iType, 1);
std::pair<double, double> rr = m_topology.dddConstants().cellSizeTrap(type, id.iSec1);
float dr = k_half * (rr.second - rr.first);
Expand Down Expand Up @@ -307,7 +305,8 @@ HGCalGeometry::CornersVec HGCalGeometry::getCorners(const DetId& detid) const {
co[i] = m_cellVec[cellIndex].getPosition(lcoord);
}
} else {
xy = m_topology.dddConstants().locateCell(id.iLay, id.iSec1, id.iSec2, id.iCell1, id.iCell2, true, false);
xy = m_topology.dddConstants().locateCell(
id.iLay, id.iSec1, id.iSec2, id.iCell1, id.iCell2, true, false, true, debugLocate);
float zz = m_topology.dddConstants().waferZ(id.iLay, true);
float dx = k_fac2 * m_cellVec[cellIndex].param()[FlatHexagon::k_r];
float dy = k_fac1 * m_cellVec[cellIndex].param()[FlatHexagon::k_R];
Expand All @@ -332,7 +331,7 @@ HGCalGeometry::CornersVec HGCalGeometry::get8Corners(const DetId& detid) const {
unsigned int cellIndex = indexFor(detid);
HGCalTopology::DecodedDetId id = m_topology.decode(detid);
if (cellIndex < m_cellVec2.size() && m_det == DetId::HGCalHSc) {
GlobalPoint v = getPosition(detid);
GlobalPoint v = getPosition(detid, false);
int type = std::min(id.iType, 1);
std::pair<double, double> rr = m_topology.dddConstants().cellSizeTrap(type, id.iSec1);
float dr = k_half * (rr.second - rr.first);
Expand Down Expand Up @@ -363,7 +362,8 @@ HGCalGeometry::CornersVec HGCalGeometry::get8Corners(const DetId& detid) const {
co[i] = m_cellVec[cellIndex].getPosition(lcoord);
}
} else {
xy = m_topology.dddConstants().locateCell(id.iLay, id.iSec1, id.iSec2, id.iCell1, id.iCell2, true, false);
xy = m_topology.dddConstants().locateCell(
id.iLay, id.iSec1, id.iSec2, id.iCell1, id.iCell2, true, false, true, debugLocate);
dx = k_fac2 * m_cellVec[cellIndex].param()[FlatHexagon::k_r];
float dy = k_fac1 * m_cellVec[cellIndex].param()[FlatHexagon::k_R];
float dz = -id.zSide * m_cellVec[cellIndex].param()[FlatHexagon::k_dZ];
Expand All @@ -379,13 +379,16 @@ HGCalGeometry::CornersVec HGCalGeometry::get8Corners(const DetId& detid) const {
return co;
}

HGCalGeometry::CornersVec HGCalGeometry::getNewCorners(const DetId& detid) const {
HGCalGeometry::CornersVec HGCalGeometry::getNewCorners(const DetId& detid, bool debug) const {
unsigned int ncorner = (m_det == DetId::HGCalHSc) ? 5 : 7;
HGCalGeometry::CornersVec co(ncorner, GlobalPoint(0, 0, 0));
unsigned int cellIndex = indexFor(detid);
HGCalTopology::DecodedDetId id = m_topology.decode(detid);
if (debug)
edm::LogVerbatim("HGCalGeom") << "NewCorners for Layer " << id.iLay << " Wafer " << id.iSec1 << ":" << id.iSec2
<< " Cell " << id.iCell1 << ":" << id.iCell2;
if (cellIndex < m_cellVec2.size() && m_det == DetId::HGCalHSc) {
GlobalPoint v = getPosition(detid);
GlobalPoint v = getPosition(detid, false);
int type = std::min(id.iType, 1);
std::pair<double, double> rr = m_topology.dddConstants().cellSizeTrap(type, id.iSec1);
float dr = k_half * (rr.second - rr.first);
Expand All @@ -407,18 +410,30 @@ HGCalGeometry::CornersVec HGCalGeometry::getNewCorners(const DetId& detid) const
float dz = -id.zSide * m_cellVec[cellIndex].param()[FlatHexagon::k_dZ];
static const int signx[] = {1, -1, -2, -1, 1, 2};
static const int signy[] = {1, 1, 0, -1, -1, 0};
#ifdef EDM_ML_DEBUG
if (debug)
edm::LogVerbatim("HGCalGeom") << "kfac " << k_fac1 << ":" << k_fac2 << " dx:dy:dz " << dx << ":" << dy << ":"
<< dz;
#endif
if (m_topology.waferHexagon6()) {
xy = m_topology.dddConstants().locateCellHex(id.iCell1, id.iSec1, true);
for (unsigned int i = 0; i < ncorner - 1; ++i) {
const HepGeom::Point3D<float> lcoord(xy.first + signx[i] * dx, xy.second + signy[i] * dy, dz);
co[i] = m_cellVec[cellIndex].getPosition(lcoord);
}
} else {
xy = m_topology.dddConstants().locateCell(id.iLay, id.iSec1, id.iSec2, id.iCell1, id.iCell2, true, false);
xy = m_topology.dddConstants().locateCell(
id.iLay, id.iSec1, id.iSec2, id.iCell1, id.iCell2, true, false, true, debug);
float zz = m_topology.dddConstants().waferZ(id.iLay, true);
for (unsigned int i = 0; i < ncorner; ++i) {
auto xyglob = m_topology.dddConstants().localToGlobal8(
id.iLay, id.iSec1, id.iSec2, (xy.first + signx[i] * dx), (xy.second + signy[i] * dy), true, false);
double xloc = xy.first + signx[i] * dx;
double yloc = xy.second + signy[i] * dy;
#ifdef EDM_ML_DEBUG
if (debug)
edm::LogVerbatim("HGCalGeom") << "Corner " << i << " x " << xy.first << ":" << xloc << " y " << xy.second
<< ":" << yloc << " z " << zz << ":" << id.zSide * (zz + dz);
#endif
auto xyglob = m_topology.dddConstants().localToGlobal8(id.iLay, id.iSec1, id.iSec2, xloc, yloc, true, debug);
double xx = id.zSide * xyglob.first;
co[i] = GlobalPoint(xx, xyglob.second, id.zSide * (zz + dz));
}
Expand All @@ -440,7 +455,7 @@ DetId HGCalGeometry::neighborZ(const DetId& idin, const GlobalVector& momentum)
#endif
if ((lay >= m_topology.dddConstants().firstLayer()) && (lay <= m_topology.dddConstants().lastLayer(true)) &&
(momentum.z() != 0.0)) {
GlobalPoint v = getPosition(idin);
GlobalPoint v = getPosition(idin, false);
double z = id.zSide * m_topology.dddConstants().waferZ(lay, true);
double grad = (z - v.z()) / momentum.z();
GlobalPoint p(v.x() + grad * momentum.x(), v.y() + grad * momentum.y(), z);
Expand Down Expand Up @@ -471,7 +486,7 @@ DetId HGCalGeometry::neighborZ(const DetId& idin,
#endif
if ((lay >= m_topology.dddConstants().firstLayer()) && (lay <= m_topology.dddConstants().lastLayer(true)) &&
(momentum.z() != 0.0)) {
GlobalPoint v = getPosition(idin);
GlobalPoint v = getPosition(idin, false);
double z = id.zSide * m_topology.dddConstants().waferZ(lay, true);
FreeTrajectoryState fts(v, momentum, charge, bField);
Plane::PlanePointer nPlane = Plane::build(Plane::PositionType(0, 0, z), Plane::RotationType());
Expand Down

0 comments on commit 250803c

Please sign in to comment.