Skip to content

Commit

Permalink
Merge pull request #35563 from bsunanda/Phase2-hgx291
Browse files Browse the repository at this point in the history
Phase2-hgx291 Bug fix to address rotated layers in D86
  • Loading branch information
cmsbuild committed Oct 19, 2021
2 parents c41e862 + e085563 commit 2e838c0
Show file tree
Hide file tree
Showing 11 changed files with 445 additions and 70 deletions.
12 changes: 8 additions & 4 deletions Geometry/HGCalCommonData/interface/HGCalDDDConstants.h
Expand Up @@ -30,7 +30,7 @@ class HGCalDDDConstants {
~HGCalDDDConstants();

std::pair<int, int> assignCell(float x, float y, int lay, int subSec, bool reco) const;
std::array<int, 5> assignCellHex(float x, float y, int lay, bool reco) const;
std::array<int, 5> assignCellHex(float x, float y, int lay, bool reco, bool extend = false, bool debug = false) const;
std::array<int, 3> assignCellTrap(float x, float y, float z, int lay, bool reco) const;
std::pair<double, double> cellEtaPhiTrap(int type, int irad) const;
bool cellInLayer(int waferU, int waferV, int cellU, int cellV, int lay, bool reco) const;
Expand Down Expand Up @@ -136,6 +136,7 @@ class HGCalDDDConstants {
int& cellV,
int& celltype,
double& wt,
bool extend = false,
bool debug = false) const;
bool waferHexagon6() const {
return ((mode_ == HGCalGeometryMode::Hexagon) || (mode_ == HGCalGeometryMode::HexagonFull));
Expand Down Expand Up @@ -212,13 +213,15 @@ class HGCalDDDConstants {
const double& cellR,
const std::vector<double>& posX,
const std::vector<double>& posY) const;
void cellHex(double xloc, double yloc, int cellType, int& cellU, int& cellV, bool debug = false) const;
void cellHex(
double xloc, double yloc, int cellType, int& cellU, int& cellV, bool extend = false, bool debug = false) const;
std::pair<int, float> getIndex(int lay, bool reco) const;
int layerFromIndex(int index, bool reco) const;
bool isValidCell(int layindex, int wafer, int cell) const;
bool isValidCell8(int lay, int waferU, int waferV, int cellU, int cellV, int type) const;
int32_t waferIndex(int wafer, int index) const;
bool waferInLayerTest(int wafer, int lay, bool full) const;
std::pair<double, double> waferPositionNoRot(int lay, int waferU, int waferV, bool reco, bool debug = false) const;
std::pair<double, double> waferPosition(int waferU, int waferV, bool reco) const;

HGCalGeomTools geomTools_;
Expand All @@ -229,9 +232,10 @@ class HGCalDDDConstants {
const HGCalParameters* hgpar_;
constexpr static double tan30deg_ = 0.5773502693;
const double sqrt3_;
const HGCalGeometryMode::GeometryMode mode_;
const bool fullAndPart_;
double rmax_, hexside_;
HGCalGeometryMode::GeometryMode mode_;
bool fullAndPart_;
double rmaxT_, hexsideT_;
int32_t tot_wafers_, modHalf_;
std::array<uint32_t, 2> tot_layers_;
Simrecovecs max_modules_layer_;
Expand Down
179 changes: 122 additions & 57 deletions Geometry/HGCalCommonData/src/HGCalDDDConstants.cc

Large diffs are not rendered by default.

1 change: 1 addition & 0 deletions Geometry/HGCalCommonData/src/HGCalParametersFromDD.cc
Expand Up @@ -502,6 +502,7 @@ void HGCalParametersFromDD::getCellPosition(HGCalParameters& php, int type) {
php.cellCoarseIndex_ = cellIndex;
else
php.cellFineIndex_ = cellIndex;

#ifdef EDM_ML_DEBUG
if (type == 1) {
edm::LogVerbatim("HGCalGeom") << "CellPosition for type " << type << " for " << php.cellCoarseX_.size()
Expand Down
1 change: 1 addition & 0 deletions Geometry/HGCalGeometry/interface/HGCalGeometry.h
Expand Up @@ -93,6 +93,7 @@ class HGCalGeometry final : public CaloSubdetectorGeometry {

// Get closest cell, etc...
DetId getClosestCell(const GlobalPoint& r) const override;
DetId getClosestCellHex(const GlobalPoint& r, bool extend) const;

/** \brief Get a list of all cells within a dR of the given cell
Expand Down
61 changes: 55 additions & 6 deletions Geometry/HGCalGeometry/src/HGCalGeometry.cc
Expand Up @@ -219,12 +219,17 @@ GlobalPoint HGCalGeometry::getPosition(const DetId& detid) const {
<< " Global " << glob;
#endif
} else {
xy = m_topology.dddConstants().locateCell(id.iLay, id.iSec1, id.iSec2, id.iCell1, id.iCell2, true, false, true);
const HepGeom::Point3D<float> lcoord(xy.first, xy.second, 0);
glob = m_cellVec[cellIndex].getPosition(lcoord);
#ifdef EDM_ML_DEBUG
edm::LogVerbatim("HGCalGeom") << "getPositionWafer:: index " << cellIndex << " Local " << lcoord.x() << ":"
<< lcoord.y() << " ID " << id.iLay << ":" << id.iSec1 << ":" << id.iSec2 << ":"
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);
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
}
Expand Down Expand Up @@ -497,7 +502,51 @@ DetId HGCalGeometry::getClosestCell(const GlobalPoint& r) const {
id.iType = kxy[2];
} else {
id.iLay = m_topology.dddConstants().getLayer(r.z(), true);
const auto& kxy = m_topology.dddConstants().assignCellHex(local.x(), local.y(), id.iLay, true);
#ifdef EDM_ML_DEBUG
edm::LogVerbatim("HGCalGeom") << "ZZ " << r.z() << " Layer " << id.iLay << " Global " << r << " Local " << local;
#endif
const auto& kxy = m_topology.dddConstants().assignCellHex(local.x(), local.y(), id.iLay, false, true);
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.det << ":" << 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 not valid or out of bounds return a null DetId
return DetId();
}

DetId HGCalGeometry::getClosestCellHex(const GlobalPoint& r, bool extend) const {
unsigned int cellIndex = getClosestCellIndex(r);
if (cellIndex < m_cellVec.size() && m_det != DetId::HGCalHSc) {
HGCalTopology::DecodedDetId id = m_topology.decode(m_validGeomIds[cellIndex]);
if (id.det == 0)
id.det = static_cast<int>(m_topology.detector());
HepGeom::Point3D<float> local;
if (r.z() > 0) {
local = HepGeom::Point3D<float>(r.x(), r.y(), 0);
id.zSide = 1;
} else {
local = HepGeom::Point3D<float>(-r.x(), r.y(), 0);
id.zSide = -1;
}
if (m_topology.waferHexagon8()) {
id.iLay = m_topology.dddConstants().getLayer(r.z(), true);
#ifdef EDM_ML_DEBUG
edm::LogVerbatim("HGCalGeom") << "ZZ " << r.z() << " Layer " << id.iLay << " Global " << r << " Local " << local;
#endif
const auto& kxy = m_topology.dddConstants().assignCellHex(local.x(), local.y(), id.iLay, extend, true);
id.iSec1 = kxy[0];
id.iSec2 = kxy[1];
id.iType = kxy[2];
Expand Down
2 changes: 1 addition & 1 deletion Geometry/HGCalGeometry/src/HGCalGeometryLoader.cc
Expand Up @@ -137,7 +137,7 @@ HGCalGeometry* HGCalGeometryLoader::build(const HGCalTopology& topology) {
int type = topology.dddConstants().getTypeHex(layer, u, v);
DetId detId = (topology.isHFNose() ? (DetId)(HFNoseDetId(zside, type, layer, u, v, 0, 0))
: (DetId)(HGCSiliconDetId(det, zside, type, layer, u, v, 0, 0)));
const auto& w = topology.dddConstants().waferPosition(layer, u, v, true);
const auto& w = topology.dddConstants().waferPosition(layer, u, v, true, true);
double xx = (zside > 0) ? w.first : -w.first;
CLHEP::Hep3Vector h3v(xx, w.second, mytr.h3v.z());
const HepGeom::Transform3D ht3d(mytr.hr, h3v);
Expand Down
87 changes: 87 additions & 0 deletions Geometry/HGCalGeometry/test/HGCalGeometryRotCheck.cc
@@ -0,0 +1,87 @@
#include <cmath>
#include <iostream>
#include <string>
#include <vector>

#include "CommonTools/UtilAlgos/interface/TFileService.h"

#include "FWCore/Framework/interface/Frameworkfwd.h"
#include "FWCore/Framework/interface/one/EDAnalyzer.h"
#include "FWCore/Framework/interface/Event.h"
#include "FWCore/Framework/interface/EventSetup.h"
#include "FWCore/MessageLogger/interface/MessageLogger.h"
#include "FWCore/ParameterSet/interface/ParameterSet.h"
#include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h"

#include "Geometry/Records/interface/IdealGeometryRecord.h"
#include "Geometry/HGCalGeometry/interface/HGCalGeometry.h"
#include "DataFormats/ForwardDetId/interface/HGCSiliconDetId.h"

class HGCalGeometryRotCheck : public edm::one::EDAnalyzer<edm::one::WatchRuns> {
public:
explicit HGCalGeometryRotCheck(const edm::ParameterSet&);
~HGCalGeometryRotCheck() override = default;
static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);

void beginRun(edm::Run const&, edm::EventSetup const&) override;
void analyze(edm::Event const& iEvent, edm::EventSetup const&) override {}
void endRun(edm::Run const&, edm::EventSetup const&) override {}

private:
const std::string nameDetector_;
const edm::ESGetToken<HGCalGeometry, IdealGeometryRecord> geomToken_;
const std::vector<int> layers_;
};

HGCalGeometryRotCheck::HGCalGeometryRotCheck(const edm::ParameterSet& iC)
: nameDetector_(iC.getParameter<std::string>("detectorName")),
geomToken_(esConsumes<HGCalGeometry, IdealGeometryRecord, edm::Transition::BeginRun>(
edm::ESInputTag{"", nameDetector_})),
layers_(iC.getParameter<std::vector<int>>("layers")) {}

void HGCalGeometryRotCheck::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
edm::ParameterSetDescription desc;
std::vector<int> layer = {27, 28, 29, 30, 31, 32, 33};
desc.add<std::string>("detectorName", "HGCalHESiliconSensitive");
desc.add<std::vector<int>>("layers", layer);
descriptions.add("hgcalGeometryRotCheck", desc);
}

void HGCalGeometryRotCheck::beginRun(const edm::Run&, const edm::EventSetup& iSetup) {
const auto& geomR = iSetup.getData(geomToken_);
const HGCalGeometry* geom = &geomR;
int layerF = *(layers_.begin());
int layerL = *(--layers_.end());
int layerOff = geom->topology().dddConstants().getLayerOffset();
edm::LogVerbatim("HGCalGeom") << nameDetector_ << " with layers in the range " << layerF << ":" << layerL
<< " Offset " << layerOff;

auto rrF = geom->topology().dddConstants().rangeRLayer(layerF - layerOff, true);
auto rrE = geom->topology().dddConstants().rangeRLayer(layerL - layerOff, true);
edm::LogVerbatim("HGCalGeom") << " RFront " << rrF.first << ":" << rrF.second << " RBack " << rrE.first << ":"
<< rrE.second;
double r = rrE.first + 5.0;
const int nPhi = 10;
while (r <= rrF.second) {
for (int k = 0; k < nPhi; ++k) {
double phi = 2.0 * k * M_PI / nPhi;
for (auto lay : layers_) {
double zz = geom->topology().dddConstants().waferZ(lay - layerOff, true);
GlobalPoint global1(r * cos(phi), r * sin(phi), zz);
DetId id = geom->getClosestCellHex(global1, true);
HGCSiliconDetId detId = HGCSiliconDetId(id);
GlobalPoint global2 = geom->getPosition(id);
double dx = global1.x() - global2.x();
double dy = global1.y() - global2.y();
double dR = std::sqrt(dx * dx + dy * dy);
edm::LogVerbatim("HGCalGeom") << "Layer: " << lay << " ID " << detId << " I/P " << global1 << " O/P " << global2
<< " dR " << dR;
}
}
r += 100.0;
}
}

//define this as a plug-in
#include "FWCore/Framework/interface/MakerMacros.h"
DEFINE_FWK_MODULE(HGCalGeometryRotCheck);
90 changes: 90 additions & 0 deletions Geometry/HGCalGeometry/test/HGCalGeometryRotTest.cc
@@ -0,0 +1,90 @@
#include <cmath>
#include <iostream>
#include <string>
#include <vector>

#include "CommonTools/UtilAlgos/interface/TFileService.h"

#include "FWCore/Framework/interface/Frameworkfwd.h"
#include "FWCore/Framework/interface/one/EDAnalyzer.h"
#include "FWCore/Framework/interface/Event.h"
#include "FWCore/Framework/interface/EventSetup.h"
#include "FWCore/MessageLogger/interface/MessageLogger.h"
#include "FWCore/ParameterSet/interface/ParameterSet.h"
#include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h"

#include "Geometry/Records/interface/IdealGeometryRecord.h"
#include "Geometry/HGCalGeometry/interface/HGCalGeometry.h"
#include "DataFormats/ForwardDetId/interface/HGCSiliconDetId.h"

class HGCalGeometryRotTest : public edm::one::EDAnalyzer<edm::one::WatchRuns> {
public:
explicit HGCalGeometryRotTest(const edm::ParameterSet&);
~HGCalGeometryRotTest() override = default;
static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);

void beginRun(edm::Run const&, edm::EventSetup const&) override;
void analyze(edm::Event const& iEvent, edm::EventSetup const&) override {}
void endRun(edm::Run const&, edm::EventSetup const&) override {}

private:
const std::string nameDetector_;
const edm::ESGetToken<HGCalGeometry, IdealGeometryRecord> geomToken_;
const std::vector<int> layers_;
const std::vector<int> waferU_, waferV_, cellU_, cellV_, types_;
};

HGCalGeometryRotTest::HGCalGeometryRotTest(const edm::ParameterSet& iC)
: nameDetector_(iC.getParameter<std::string>("detectorName")),
geomToken_(esConsumes<HGCalGeometry, IdealGeometryRecord, edm::Transition::BeginRun>(
edm::ESInputTag{"", nameDetector_})),
layers_(iC.getParameter<std::vector<int>>("layers")),
waferU_(iC.getParameter<std::vector<int>>("waferUs")),
waferV_(iC.getParameter<std::vector<int>>("waferVs")),
cellU_(iC.getParameter<std::vector<int>>("cellUs")),
cellV_(iC.getParameter<std::vector<int>>("cellVs")),
types_(iC.getParameter<std::vector<int>>("types")) {}

void HGCalGeometryRotTest::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
edm::ParameterSetDescription desc;
std::vector<int> layer = {27, 28, 29, 30, 31, 32};
std::vector<int> waferU = {-2, -3, 1, -8, 2, 8};
std::vector<int> waferV = {0, -2, 3, 0, 9, 0};
std::vector<int> cellU = {16, 4, 8, 11, 11, 5};
std::vector<int> cellV = {20, 10, 17, 13, 9, 2};
std::vector<int> type = {0, 0, 0, 2, 2, 2};
desc.add<std::string>("detectorName", "HGCalHESiliconSensitive");
desc.add<std::vector<int>>("layers", layer);
desc.add<std::vector<int>>("waferUs", waferU);
desc.add<std::vector<int>>("waferVs", waferV);
desc.add<std::vector<int>>("cellUs", cellU);
desc.add<std::vector<int>>("cellVs", cellV);
desc.add<std::vector<int>>("types", type);
descriptions.add("hgcalGeometryRotTest", desc);
}

void HGCalGeometryRotTest::beginRun(const edm::Run&, const edm::EventSetup& iSetup) {
const auto& geomR = iSetup.getData(geomToken_);
const HGCalGeometry* geom = &geomR;
DetId::Detector det = (nameDetector_ == "HGCalEESensitive") ? DetId::HGCalEE : DetId::HGCalHSi;
int layerF = *(layers_.begin());
int layerL = *(--layers_.end());
int layerOff = geom->topology().dddConstants().getLayerOffset();
edm::LogVerbatim("HGCalGeom") << nameDetector_ << " with layers in the range " << layerF << ":" << layerL
<< " Offset " << layerOff << " and for " << waferU_.size() << " wafers and cells";

for (unsigned int k = 0; k < waferU_.size(); ++k) {
for (auto lay : layers_) {
HGCSiliconDetId detId(det, 1, types_[k], lay - layerOff, waferU_[k], waferV_[k], cellU_[k], cellV_[k]);
GlobalPoint global = geom->getPosition(DetId(detId));
auto xy = geom->topology().dddConstants().waferPosition(lay - layerOff, waferU_[k], waferV_[k], true);
edm::LogVerbatim("HGCalGeom") << "Layer: " << lay << " U " << waferU_[k] << " V " << waferV_[k] << " Position ("
<< xy.first << ", " << xy.second << ")";
edm::LogVerbatim("HGCalGeom") << detId << " Position " << global;
}
}
}

//define this as a plug-in
#include "FWCore/Framework/interface/MakerMacros.h"
DEFINE_FWK_MODULE(HGCalGeometryRotTest);
@@ -0,0 +1,39 @@
import FWCore.ParameterSet.Config as cms

from Configuration.Eras.Era_Phase2C11I13M9_cff import Phase2C11I13M9
process = cms.Process('PROD',Phase2C11I13M9)
process.load('Configuration.Geometry.GeometryExtended2026D86_cff')
process.load('Configuration.Geometry.GeometryExtended2026D86Reco_cff')

process.load("SimGeneral.HepPDTESSource.pdt_cfi")
process.load('Geometry.HGCalGeometry.hgcalGeometryRotCheck_cfi')
process.load('FWCore.MessageService.MessageLogger_cfi')

if hasattr(process,'MessageLogger'):
process.MessageLogger.HGCalGeom=dict()

process.load("IOMC.RandomEngine.IOMC_cff")
process.RandomNumberGeneratorService.generator.initialSeed = 456789

process.source = cms.Source("EmptySource")

process.generator = cms.EDProducer("FlatRandomEGunProducer",
PGunParameters = cms.PSet(
PartID = cms.vint32(14),
MinEta = cms.double(-3.5),
MaxEta = cms.double(3.5),
MinPhi = cms.double(-3.14159265359),
MaxPhi = cms.double(3.14159265359),
MinE = cms.double(9.99),
MaxE = cms.double(10.01)
),
AddAntiParticle = cms.bool(False),
Verbosity = cms.untracked.int32(0),
firstRun = cms.untracked.uint32(1)
)

process.maxEvents = cms.untracked.PSet(
input = cms.untracked.int32(1)
)

process.p1 = cms.Path(process.generator*process.hgcalGeometryRotCheck)

0 comments on commit 2e838c0

Please sign in to comment.