Skip to content

Commit

Permalink
Fix a bug for the corner calculations
Browse files Browse the repository at this point in the history
  • Loading branch information
Sunanda committed Oct 26, 2018
1 parent 768b8a5 commit 6c90fb5
Show file tree
Hide file tree
Showing 5 changed files with 128 additions and 64 deletions.
1 change: 1 addition & 0 deletions Geometry/HGCalGeometry/interface/HGCalGeometry.h
View file
Expand Up @@ -79,6 +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;

// avoid sorting set in base class
const std::vector<DetId>& getValidDetIds( DetId::Detector det = DetId::Detector(0), int subdet = 0) const override { return m_validIds; }
Expand Down
159 changes: 108 additions & 51 deletions Geometry/HGCalGeometry/src/HGCalGeometry.cc
View file
Expand Up @@ -269,34 +269,36 @@ GlobalPoint HGCalGeometry::getPosition(const DetId& id) const {
return glob;
}

HGCalGeometry::CornersVec HGCalGeometry::getCorners(const DetId& id) const {
HGCalGeometry::CornersVec HGCalGeometry::getCorners(const DetId& detid) const {

unsigned int ncorner = ((m_det == DetId::HGCalHSc) ? FlatTrd::ncorner_ :
FlatHexagon::ncorner_);
HGCalGeometry::CornersVec co (ncorner, GlobalPoint(0,0,0));
unsigned int cellIndex = indexFor(id);
unsigned int cellIndex = indexFor(detid);
if ((cellIndex < m_cellVec.size() && m_det != DetId::HGCalHSc) ||
(cellIndex < m_cellVec2.size() && m_det == DetId::HGCalHSc)) {
HGCalTopology::DecodedDetId id_ = m_topology.decode(id);
HGCalTopology::DecodedDetId id = m_topology.decode(detid);
std::pair<float,float> xy;
if ((mode_ == HGCalGeometryMode::Hexagon) ||
(mode_ == HGCalGeometryMode::HexagonFull)) {
xy = m_topology.dddConstants().locateCellHex(id_.iCell1,id_.iSec1,true);
xy = m_topology.dddConstants().locateCellHex(id.iCell1,id.iSec1,true);
} else if (mode_ == HGCalGeometryMode::Trapezoid) {
xy = m_topology.dddConstants().locateCellTrap(id_.iLay,id_.iSec1,
id_.iCell1,true);
xy = m_topology.dddConstants().locateCellTrap(id.iLay,id.iSec1,
id.iCell1,true);
} 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);
}
if (m_det == DetId::HGCalHSc) {
float dx = 0.5*m_cellVec2[cellIndex].param()[11];
std::pair<double,double> rr = m_topology.dddConstants().cellSizeTrap(id.iType,id.iSec1);
float dx = 0.5*(rr.second-rr.first);
float dy = 0.5*(rr.second+rr.first)*m_cellVec2[cellIndex].param()[11];
float dz = m_cellVec2[cellIndex].param()[0];
static const int signx[] = {-1,-1,1,1,-1,-1,1,1};
static const int signy[] = {-1,1,1,-1,-1,1,1,-1};
static const int signz[] = {-1,-1,-1,-1,1,1,1,1};
for (unsigned int i = 0; i < ncorner; ++i) {
const HepGeom::Point3D<float> lcoord(xy.first+signx[i]*dx,xy.second+signy[i]*dx,signz[i]*dz);
const HepGeom::Point3D<float> lcoord(signx[i]*dx,signy[i]*dy,signz[i]*dz);
co[i] = m_cellVec2[cellIndex].getPosition(lcoord);
}
} else {
Expand All @@ -315,38 +317,93 @@ HGCalGeometry::CornersVec HGCalGeometry::getCorners(const DetId& id) const {
return co;
}

HGCalGeometry::CornersVec HGCalGeometry::get8Corners(const DetId& id) const {
HGCalGeometry::CornersVec HGCalGeometry::get8Corners(const DetId& detid) const {

unsigned int ncorner = FlatTrd::ncorner_;
HGCalGeometry::CornersVec co (ncorner, GlobalPoint(0,0,0));
unsigned int cellIndex = indexFor(id);
unsigned int cellIndex = indexFor(detid);
if ((cellIndex < m_cellVec.size() && m_det != DetId::HGCalHSc) ||
(cellIndex < m_cellVec2.size() && m_det == DetId::HGCalHSc)) {
HGCalTopology::DecodedDetId id_ = m_topology.decode(id);
HGCalTopology::DecodedDetId id = m_topology.decode(detid);
std::pair<float,float> xy;
if ((mode_ == HGCalGeometryMode::Hexagon) ||
(mode_ == HGCalGeometryMode::HexagonFull)) {
xy = m_topology.dddConstants().locateCellHex(id_.iCell1,id_.iSec1,true);
xy = m_topology.dddConstants().locateCellHex(id.iCell1,id.iSec1,true);
} else if (mode_ == HGCalGeometryMode::Trapezoid) {
xy = m_topology.dddConstants().locateCellTrap(id_.iLay,id_.iSec1,
id_.iCell1,true);
xy = m_topology.dddConstants().locateCellTrap(id.iLay,id.iSec1,
id.iCell1,true);
} 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);
}
float dx = ((m_det == DetId::HGCalHSc) ? 0.5*m_cellVec2[cellIndex].param()[11] :
m_cellVec[cellIndex].param()[1]);
float dz = ((m_det == DetId::HGCalHSc) ? m_cellVec2[cellIndex].param()[0] :
m_cellVec[cellIndex].param()[0]);
static const int signx[] = {-1,-1,1,1,-1,-1,1,1};
static const int signy[] = {-1,1,1,-1,-1,1,1,-1};
static const int signz[] = {-1,-1,-1,-1,1,1,1,1};
for (unsigned int i = 0; i < ncorner; ++i) {
const HepGeom::Point3D<float> lcoord(xy.first+signx[i]*dx,
xy.second+signy[i]*dx,signz[i]*dz);
co[i] = ((m_det == DetId::HGCalHSc) ?
(m_cellVec2[cellIndex].getPosition(lcoord)) :
(m_cellVec[cellIndex].getPosition(lcoord)));
if (m_det == DetId::HGCalHSc) {
std::pair<double,double> rr = m_topology.dddConstants().cellSizeTrap(id.iType,id.iSec1);
float dx = 0.5*(rr.second-rr.first);
float dy = 0.5*(rr.second+rr.first)*m_cellVec2[cellIndex].param()[11];
float dz = m_cellVec2[cellIndex].param()[0];
for (unsigned int i = 0; i < ncorner; ++i) {
const HepGeom::Point3D<float> lcoord(signx[i]*dx,signy[i]*dy,signz[i]*dz);
co[i] = m_cellVec2[cellIndex].getPosition(lcoord);
}
} else {
float dx = m_cellVec[cellIndex].param()[1];
float dz = m_cellVec[cellIndex].param()[0];
for (unsigned int i = 0; i < ncorner; ++i) {
const HepGeom::Point3D<float> lcoord(xy.first+signx[i]*dx,
xy.second+signy[i]*dx,signz[i]*dz);
co[i] = m_cellVec[cellIndex].getPosition(lcoord);
}
}
}
return co;
}

HGCalGeometry::CornersVec HGCalGeometry::getNewCorners(const DetId& detid) const {

unsigned int ncorner = (m_det == DetId::HGCalHSc) ? 5 : 7;
HGCalGeometry::CornersVec co (ncorner, GlobalPoint(0,0,0));
unsigned int cellIndex = indexFor(detid);
if ((cellIndex < m_cellVec.size() && m_det != DetId::HGCalHSc) ||
(cellIndex < m_cellVec2.size() && m_det == DetId::HGCalHSc)) {
HGCalTopology::DecodedDetId id = m_topology.decode(detid);
int signz = (id.zSide >= 0) ? -1 : 1;
std::pair<float,float> xy;
if ((mode_ == HGCalGeometryMode::Hexagon) ||
(mode_ == HGCalGeometryMode::HexagonFull)) {
xy = m_topology.dddConstants().locateCellHex(id.iCell1,id.iSec1,true);
} else if (mode_ == HGCalGeometryMode::Trapezoid) {
xy = m_topology.dddConstants().locateCellTrap(id.iLay,id.iSec1,
id.iCell1,true);
} else {
xy = m_topology.dddConstants().locateCell(id.iLay,id.iSec1,id.iSec2,
id.iCell1,id.iCell2,true,false);
}
if (m_det == DetId::HGCalHSc) {
std::pair<double,double> rr = m_topology.dddConstants().cellSizeTrap(id.iType,id.iSec1);
float dx = 0.5*(rr.second-rr.first);
float dy = 0.5*(rr.second+rr.first)*m_cellVec2[cellIndex].param()[11];
float dz = m_cellVec2[cellIndex].param()[0];
static const int signx[] = {-1,-1,1,1};
static const int signy[] = {-1,1,1,-1};
for (unsigned int i = 0; i < ncorner-1; ++i) {
const HepGeom::Point3D<float> lcoord(signx[i]*dx,signy[i]*dy,signz*dz);
co[i] = m_cellVec2[cellIndex].getPosition(lcoord);
}
co[ncorner-1] = GlobalPoint(0,0,-2*signz*dz);
} else {
float dx = m_cellVec[cellIndex].param()[1];
float dy = m_cellVec[cellIndex].param()[2];
float dz = m_cellVec[cellIndex].param()[0];
static const int signx[] = {0,-1,-1,0,1,1};
static const int signy[] = {-2,-1,1,2,1,-1};
for (unsigned int i = 0; i < ncorner-1; ++i) {
const HepGeom::Point3D<float> lcoord(xy.first+signx[i]*dx,xy.second+signy[i]*dy,signz*dz);
co[i] = m_cellVec[cellIndex].getPosition(lcoord);
}
co[ncorner-1] = GlobalPoint(0,0,-2*signz*dz);
}
}
return co;
Expand All @@ -356,50 +413,50 @@ DetId HGCalGeometry::getClosestCell(const GlobalPoint& r) const {
unsigned int cellIndex = getClosestCellIndex(r);
if ((cellIndex < m_cellVec.size() && m_det != DetId::HGCalHSc) ||
(cellIndex < m_cellVec2.size() && m_det == DetId::HGCalHSc)) {
HGCalTopology::DecodedDetId id_ = m_topology.decode(m_validGeomIds[cellIndex]);
HGCalTopology::DecodedDetId id = m_topology.decode(m_validGeomIds[cellIndex]);
HepGeom::Point3D<float> local;
if (r.z() > 0) {
local = HepGeom::Point3D<float>(r.x(),r.y(),0);
id_.zSide = 1;
id.zSide = 1;
} else {
local = HepGeom::Point3D<float>(-r.x(),r.y(),0);
id_.zSide =-1;
id.zSide =-1;
}
if ((mode_ == HGCalGeometryMode::Hexagon) ||
(mode_ == HGCalGeometryMode::HexagonFull)) {
const auto & kxy =
m_topology.dddConstants().assignCell(local.x(),local.y(),id_.iLay,
id_.iType,true);
id_.iCell1 = kxy.second;
id_.iSec1 = kxy.first;
id_.iType = m_topology.dddConstants().waferTypeT(kxy.first);
if (id_.iType != 1) id_.iType = -1;
m_topology.dddConstants().assignCell(local.x(),local.y(),id.iLay,
id.iType,true);
id.iCell1 = kxy.second;
id.iSec1 = kxy.first;
id.iType = m_topology.dddConstants().waferTypeT(kxy.first);
if (id.iType != 1) id.iType = -1;
} else if (mode_ == HGCalGeometryMode::Trapezoid) {
id_.iLay = m_topology.dddConstants().getLayer(r.z(),true);
id.iLay = m_topology.dddConstants().getLayer(r.z(),true);
const auto & kxy =
m_topology.dddConstants().assignCellTrap(r.x(),r.y(),r.z(),
id_.iLay,true);
id_.iSec1 = kxy[0];
id_.iCell1 = kxy[1];
id_.iType = kxy[2];
id.iLay,true);
id.iSec1 = kxy[0];
id.iCell1 = kxy[1];
id.iType = kxy[2];
} else {
id_.iLay = m_topology.dddConstants().getLayer(r.z(),true);
id.iLay = m_topology.dddConstants().getLayer(r.z(),true);
const auto & kxy =
m_topology.dddConstants().assignCellHex(local.x(),local.y(),id_.iLay,
m_topology.dddConstants().assignCellHex(local.x(),local.y(),id.iLay,
true);
id_.iSec1 = kxy[0]; id_.iSec2 = kxy[1]; id_.iType = kxy[2];
id_.iCell1 = kxy[3]; id_.iCell2 = kxy[4];
id.iSec1 = kxy[0]; id.iSec2 = kxy[1]; id.iType = kxy[2];
id.iCell1 = kxy[3]; id.iCell2 = kxy[4];
}
#ifdef EDM_ML_DEBUG
edm::LogVerbatim("HGCalGeom") << "getClosestCell: local " << local
<< " Id " << id_.zSide << ":" << id_.iLay
<< ":" << id_.iSec1 << ":" << id_.iSec2
<< ":" << id_.iType << ":" << id_.iCell1
<< ":" << id_.iCell2;
<< " Id " << id.zSide << ":" << id.iLay
<< ":" << id.iSec1 << ":" << id.iSec2
<< ":" << id.iType << ":" << id.iCell1
<< ":" << id.iCell2;
#endif

//check if returned cell is valid
if (id_.iCell1>=0) return m_topology.encode(id_);
if (id.iCell1>=0) return m_topology.encode(id);
}

//if not valid or out of bounds return a null DetId
Expand Down
3 changes: 1 addition & 2 deletions Geometry/HGCalGeometry/src/HGCalGeometryLoader.cc
View file
Expand Up @@ -39,8 +39,7 @@ HGCalGeometry* HGCalGeometryLoader::build (const HGCalTopology& topology) {
<< numberExpected << " for sub-detector "
<< topology.subDetector() << " Shapes "
<< numberOfShapes << ":" << parametersPerShape_
<< " mode " << mode << " ifNose "
<< topology.isHFNose();
<< " mode " << mode;
#endif
geom->allocateCorners( numberOfCells ) ;
geom->allocatePar(numberOfShapes, parametersPerShape_);
Expand Down
21 changes: 13 additions & 8 deletions Geometry/HGCalGeometry/test/HGCalGeometryCornerTester.cc
View file
Expand Up @@ -29,22 +29,24 @@ class HGCalGeometryCornerTester : public edm::one::EDAnalyzer<> {

private:
void doTest(const HGCalGeometry* geom);
const std::string name;
const std::string name_;
const int cornerType_;
};

HGCalGeometryCornerTester::HGCalGeometryCornerTester(const edm::ParameterSet& iC) : name(iC.getParameter<std::string>("Detector")) { }

HGCalGeometryCornerTester::HGCalGeometryCornerTester(const edm::ParameterSet& iC) :
name_(iC.getParameter<std::string>("detector")),
cornerType_(iC.getParameter<int>("cornerType")) { }

HGCalGeometryCornerTester::~HGCalGeometryCornerTester() {}

void HGCalGeometryCornerTester::analyze(const edm::Event& ,
const edm::EventSetup& iSetup ) {

edm::ESHandle<HGCalGeometry> geomH;
iSetup.get<IdealGeometryRecord>().get(name,geomH);
iSetup.get<IdealGeometryRecord>().get(name_,geomH);
const HGCalGeometry* geom = (geomH.product());
if (!geomH.isValid()) {
std::cout << "Cannot get valid HGCalGeometry Object for " << name
std::cout << "Cannot get valid HGCalGeometry Object for " << name_
<< std::endl;
} else {
doTest(geom);
Expand All @@ -63,10 +65,13 @@ void HGCalGeometryCornerTester::doTest(const HGCalGeometry* geom) {
else if (id.det() == DetId::HGCalHSc) std::cout <<HGCScintillatorDetId(id);
else std::cout <<HGCSiliconDetId(id);
std::cout << std::endl;
const auto cor = geom->getCorners(id);
std::cout << cor.size() << " Corners:";
const auto cor = ((cornerType_ > 0) ? geom->getCorners(id) :
((cornerType_ < 0) ? geom->get8Corners(id) :
geom->getNewCorners(id)));
GlobalPoint gp = geom->getPosition(id);
std::cout << "Center" << gp << "with " << cor.size() << " Corners: ";
for (unsigned k=0; k<cor.size(); ++k)
std::cout << " [" << k << "] " << cor[k];
std::cout << "[" << k << "]" << cor[k];
std::cout << std::endl;
++kount;
}
Expand Down
8 changes: 5 additions & 3 deletions Geometry/HGCalGeometry/test/python/testHGCalCorner_cfg.py
View file
Expand Up @@ -42,16 +42,18 @@
)

process.prodEE = cms.EDAnalyzer("HGCalGeometryCornerTester",
Detector = cms.string("HGCalEESensitive"),
detector = cms.string("HGCalEESensitive"),
cornerType = cms.int32(1)
)

process.prodHEF = process.prodEE.clone(
Detector = "HGCalHESiliconSensitive",
detector = "HGCalHESiliconSensitive",
)

process.prodHEB = process.prodEE.clone(
Detector = "HGCalHEScintillatorSensitive",
detector = "HGCalHEScintillatorSensitive",
)

#process.p1 = cms.Path(process.generator*process.prodEE*process.prodHEF)
#process.p1 = cms.Path(process.prodHEB)
process.p1 = cms.Path(process.generator*process.prodEE*process.prodHEF*process.prodHEB)

0 comments on commit 6c90fb5

Please sign in to comment.