18 #include <Math/Transform3D.h>
19 #include <Math/EulerAngles.h>
27 : m_topology(topology_),
28 m_validGeomIds(topology_.totalGeomModules()),
30 m_subdet(topology_.subDetector()),
31 twoBysqrt3_(2.0 /
std::
sqrt(3.0)) {
94 for (
int cell = 0; cell <
cells; ++cell) {
109 if (typm.first >= 0) {
112 idc = static_cast<DetId>(hid);
124 unsigned int cellAll(0), cellSelect(0);
126 for (
int u = 0; u < 2 *
cells; ++u) {
127 for (
int v = 0;
v < 2 *
cells; ++
v) {
146 edm::LogVerbatim(
"HGCalGeom") <<
"HGCalGeometry keeps " << cellSelect <<
" out of " << cellAll <<
" for wafer "
147 <<
id.iSec1 <<
":" <<
id.iSec2 <<
" in "
148 <<
" layer " <<
id.iLay;
153 edm::LogVerbatim(
"HGCalGeom") <<
"HGCalGeometry::newCell-> [" << cellIndex <<
"]"
154 <<
" front:" <<
f1.x() <<
'/' <<
f1.y() <<
'/' <<
f1.z() <<
" back:" <<
f2.x() <<
'/'
155 <<
f2.y() <<
'/' <<
f2.z() <<
" eta|phi " <<
m_cellVec2[cellIndex].etaPos() <<
":"
158 edm::LogVerbatim(
"HGCalGeom") <<
"HGCalGeometry::newCell-> [" << cellIndex <<
"]"
159 <<
" front:" <<
f1.x() <<
'/' <<
f1.y() <<
'/' <<
f1.z() <<
" back:" <<
f2.x() <<
'/'
160 <<
f2.y() <<
'/' <<
f2.z() <<
" eta|phi " <<
m_cellVec[cellIndex].etaPos() <<
":"
182 if (detId ==
DetId())
191 if (detId ==
DetId())
199 unsigned int cellIndex =
indexFor(detid);
204 std::pair<float, float>
xy;
207 const HepGeom::Point3D<float> lcoord(
xy.first,
xy.second, 0);
208 glob =
m_cellVec[cellIndex].getPosition(lcoord);
210 edm::LogVerbatim(
"HGCalGeom") <<
"getPosition:: index " << cellIndex <<
" Local " << lcoord.x() <<
":"
211 << lcoord.y() <<
" ID " <<
id.iCell1 <<
":" <<
id.iSec1 <<
" Global " << glob;
214 const HepGeom::Point3D<float> lcoord(0, 0, 0);
215 glob =
m_cellVec2[cellIndex].getPosition(lcoord);
217 edm::LogVerbatim(
"HGCalGeom") <<
"getPositionTrap:: index " << cellIndex <<
" Local " << lcoord.x() <<
":"
218 << lcoord.y() <<
" ID " <<
id.iLay <<
":" <<
id.iSec1 <<
":" <<
id.iCell1
219 <<
" Global " << glob;
223 const HepGeom::Point3D<float> lcoord(
xy.first,
xy.second, 0);
224 glob =
m_cellVec[cellIndex].getPosition(lcoord);
226 edm::LogVerbatim(
"HGCalGeom") <<
"getPositionWafer:: index " << cellIndex <<
" Local " << lcoord.x() <<
":"
227 << lcoord.y() <<
" ID " <<
id.iLay <<
":" <<
id.iSec1 <<
":" <<
id.iSec2 <<
":"
228 <<
id.iCell1 <<
":" <<
id.iCell2 <<
" Global " << glob;
236 unsigned int cellIndex =
indexFor(detid);
240 const HepGeom::Point3D<float> lcoord(0, 0, 0);
242 glob =
m_cellVec2[cellIndex].getPosition(lcoord);
244 glob =
m_cellVec[cellIndex].getPosition(lcoord);
248 << cellIndex <<
" Global " << glob;
257 if (corners.size() > 1) {
258 int n = corners.size() - 1;
260 for (
int i = 0;
i <
n; ++
i) {
261 area += ((corners[
j].x() + corners[
i].x()) * (corners[
i].
y() - corners[
j].y()));
271 unsigned int cellIndex =
indexFor(detid);
281 static const int signr[] = {1, 1, -1, -1, 1, 1, -1, -1};
282 static const int signf[] = {-1, 1, 1, -1, -1, 1, 1, -1};
283 static const int signz[] = {-1, -1, -1, -1, 1, 1, 1, 1};
284 for (
unsigned int i = 0;
i < ncorner; ++
i) {
286 (
r + signr[
i] *
dr) *
sin(fi + signf[
i] * dfi),
287 (
v.z() + signz[
i] *
dz));
290 std::pair<float, float>
xy;
296 static const int signx[] = {0, -1, -1, 0, 1, 1, 0, -1, -1, 0, 1, 1};
297 static const int signy[] = {-2, -1, 1, 2, 1, -1, -2, -1, 1, 2, 1, -1};
298 static const int signz[] = {-1, -1, -1, -1, -1, -1, 1, 1, 1, 1, 1, 1};
299 for (
unsigned int i = 0;
i < ncorner; ++
i) {
300 const HepGeom::Point3D<float> lcoord(
xy.first + signx[
i] *
dx,
xy.second + signy[
i] *
dy, signz[
i] *
dz);
308 static const int signx[] = {1, -1, -2, -1, 1, 2, 1, -1, -2, -1, 1, 2};
309 static const int signy[] = {1, 1, 0, -1, -1, 0, 1, 1, 0, -1, -1, 0};
310 static const int signz[] = {-1, -1, -1, -1, -1, -1, 1, 1, 1, 1, 1, 1};
311 for (
unsigned int i = 0;
i < ncorner; ++
i) {
312 const HepGeom::Point3D<float> lcoord(
xy.first + signx[
i] *
dx,
xy.second + signy[
i] *
dy, signz[
i] *
dz);
323 unsigned int cellIndex =
indexFor(detid);
333 static const int signr[] = {1, 1, -1, -1, 1, 1, -1, -1};
334 static const int signf[] = {-1, 1, 1, -1, -1, 1, 1, -1};
335 static const int signz[] = {-1, -1, -1, -1, 1, 1, 1, 1};
336 for (
unsigned int i = 0;
i < ncorner; ++
i) {
338 (
r + signr[
i] *
dr) *
sin(fi + signf[
i] * dfi),
339 (
v.z() + signz[
i] *
dz));
342 std::pair<float, float>
xy;
351 static const int signx[] = {-1, -1, 1, 1, -1, -1, 1, 1};
352 static const int signy[] = {-1, 1, 1, -1, -1, 1, 1, -1};
353 static const int signz[] = {-1, -1, -1, -1, 1, 1, 1, 1};
355 for (
unsigned int i = 0;
i < ncorner; ++
i) {
356 const HepGeom::Point3D<float> lcoord(
xy.first + signx[
i] *
dx,
xy.second + signy[
i] *
dx, signz[
i] *
dz);
366 unsigned int cellIndex =
indexFor(detid);
376 static const int signr[] = {1, 1, -1, -1};
377 static const int signf[] = {-1, 1, 1, -1};
378 for (
unsigned int i = 0;
i < ncorner - 1; ++
i) {
380 (
r + signr[
i] *
dr) *
cos(fi + signf[
i] * dfi), (
r + signr[
i] *
dr) *
sin(fi + signf[
i] * dfi), (
v.z() +
dz));
384 std::pair<float, float>
xy;
393 static const int signx[] = {1, -1, -2, -1, 1, 2};
394 static const int signy[] = {1, 1, 0, -1, -1, 0};
395 for (
unsigned int i = 0;
i < ncorner - 1; ++
i) {
396 const HepGeom::Point3D<float> lcoord(
xy.first + signx[
i] *
dx,
xy.second + signy[
i] *
dy,
dz);
407 int lay = ((momentum.
z() *
id.zSide > 0) ? (
id.iLay + 1) : (
id.iLay - 1));
409 edm::LogVerbatim(
"HGCalGeom") <<
"neighborz1:: ID " <<
id.iLay <<
":" <<
id.iSec1 <<
":" <<
id.iSec2 <<
":"
410 <<
id.iCell1 <<
":" <<
id.iCell2 <<
" New Layer " << lay <<
" Range "
415 (momentum.
z() != 0.0)) {
418 double grad = (
z -
v.z()) / momentum.
z();
422 if (
r >= rlimit.first &&
r <= rlimit.second)
425 edm::LogVerbatim(
"HGCalGeom") <<
"neighborz1:: Position " <<
v <<
" New Z " <<
z <<
":" << grad <<
" new position "
426 <<
p <<
" r-limit " << rlimit.first <<
":" << rlimit.second;
438 int lay = ((momentum.
z() *
id.zSide > 0) ? (
id.iLay + 1) : (
id.iLay - 1));
440 edm::LogVerbatim(
"HGCalGeom") <<
"neighborz2:: ID " <<
id.iLay <<
":" <<
id.iSec1 <<
":" <<
id.iSec2 <<
":"
441 <<
id.iCell1 <<
":" <<
id.iCell2 <<
" New Layer " << lay <<
" Range "
446 (momentum.
z() != 0.0)) {
458 if (
r >= rlimit.first &&
r <= rlimit.second)
463 << tsos.
isValid() <<
" new position " <<
p <<
" r limits " << rlimit.first <<
":"
477 HepGeom::Point3D<float>
local;
479 local = HepGeom::Point3D<float>(
r.x(),
r.y(), 0);
482 local = HepGeom::Point3D<float>(-
r.x(),
r.y(), 0);
487 id.iCell1 = kxy.second;
488 id.iSec1 = kxy.first;
508 edm::LogVerbatim(
"HGCalGeom") <<
"getClosestCell: local " <<
local <<
" Id " <<
id.det <<
":" <<
id.zSide <<
":"
509 <<
id.iLay <<
":" <<
id.iSec1 <<
":" <<
id.iSec2 <<
":" <<
id.iType <<
":"
510 <<
id.iCell1 <<
":" <<
id.iCell2;
531 return "HGCalHEFront";
533 return "HGCalHEBack";
540 if (detId !=
DetId()) {
545 <<
" index " << cellIndex;
559 return (
nullptr == cell->
param() ? nullptr : cell);
564 return (
nullptr == cell->
param() ? nullptr : cell);
572 static const auto do_not_delete = [](
const void*) {};
574 auto cell = std::shared_ptr<const CaloCellGeometry>(&
m_cellVec2[
index], do_not_delete);
575 if (
nullptr == cell->param())
579 auto cell = std::shared_ptr<const CaloCellGeometry>(&
m_cellVec[
index], do_not_delete);
580 if (
nullptr == cell->param())
594 cell->setPosition(
pos);
598 if (
nullptr == cell->param())
603 cell->setPosition(
pos);
607 if (
nullptr == cell->param())
614 edm::LogError(
"HGCalGeom") <<
"HGCalGeometry::addValidID is not implemented";
623 float phip =
r.phi();
625 float dzmin(9999), dphimin(9999), dphi10(0.175);
626 unsigned int cellIndex = vec.size();
627 for (
unsigned int k = 0;
k < vec.size(); ++
k) {
628 float dphi = phip - vec[
k].phiPos();
631 while (dphi <= -
M_PI)
635 if (
dz < (dzmin + 0.001)) {
637 if (
std::abs(dphi) < (dphimin + 0.01)) {
641 if (cellIndex >= vec.size())
648 edm::LogVerbatim(
"HGCalGeom") <<
"getClosestCellIndex::Input " << zp <<
":" << phip <<
" Index " << cellIndex;
649 if (cellIndex < vec.size())
650 edm::LogVerbatim(
"HGCalGeom") <<
" Cell z " << vec[cellIndex].getPosition().z() <<
":" << dzmin <<
" phi "
651 << vec[cellIndex].phiPos() <<
":" << dphimin;
659 bool operator()(
const DetId&
a,
const DetId&
b) {
return (
a.rawId() <
b.rawId()); }
678 iVector.reserve(numberOfCells);
680 dinsVector.reserve(numberOfCells);
694 dimVector.insert(dimVector.end(),
params.begin(),
params.end());
711 dimVector.insert(dimVector.end(),
params.begin(),
params.end());
721 dimVector.insert(dimVector.end(),
params.begin(),
params.end());
727 for (
unsigned int i(0);
i < numberOfCells; ++
i) {
740 iVector.emplace_back(
layer);
744 if (
nullptr != ptr) {
745 ptr->getTransform(tr, (
Pt3DVec*)
nullptr);
749 tr = HepGeom::Translate3D(
gp.x(),
gp.y(),
gp.z());
752 const CLHEP::Hep3Vector
tt(tr.getTranslation());
753 trVector.emplace_back(
tt.x());
754 trVector.emplace_back(
tt.y());
755 trVector.emplace_back(
tt.z());
757 const CLHEP::HepRotation
rr(tr.getRotation());
758 const ROOT::Math::Transform3D rtr(
759 rr.xx(),
rr.xy(),
rr.xz(),
tt.x(),
rr.yx(),
rr.yy(),
rr.yz(),
tt.y(),
rr.zx(),
rr.zy(),
rr.zz(),
tt.z());
762 trVector.emplace_back(ea.Phi());
763 trVector.emplace_back(ea.Theta());
764 trVector.emplace_back(ea.Psi());
773 geomId = static_cast<DetId>(
HGCalDetId(detId).geometryCell());
777 geomId = static_cast<DetId>(
HFNoseDetId(detId).geometryCell());