Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Phase2-hgx167 Fix a bug for the corner calculations #25027

Merged
merged 1 commit into from Nov 1, 2018
Merged

Phase2-hgx167 Fix a bug for the corner calculations #25027

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Jump to
Jump to file
Failed to load files.
Diff view
Diff view
1 change: 1 addition & 0 deletions Geometry/HGCalGeometry/interface/HGCalGeometry.h
View file
Open in desktop
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
Open in desktop
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];
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Please don't use magic numbers like 5,7, and 11. Assign constants and use them instead.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

These numbers are standard in various types of calorimeter cells in CaloGeometry. I shall define the names in CaloGeometry package later on and transmit them to use cases in derived geometries.

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I created an issue for this task: #25078

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
Open in desktop
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
Open in desktop
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
Open in desktop
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)