diff --git a/Geometry/HGCalCommonData/interface/HGCalDDDConstants.h b/Geometry/HGCalCommonData/interface/HGCalDDDConstants.h index 7c9ffb888ff92..2194c01af7a7a 100644 --- a/Geometry/HGCalCommonData/interface/HGCalDDDConstants.h +++ b/Geometry/HGCalCommonData/interface/HGCalDDDConstants.h @@ -113,12 +113,21 @@ class HGCalDDDConstants { std::pair rowColumnWafer(const int wafer) const; int sectors() const { return hgpar_->nSectors_; } std::pair simToReco(int cell, int layer, int mod, bool half) const; + int tileCount(int layer, int ring) const; bool tileExist(int zside, int layer, int ring, int phi) const { int indx = HGCalTileIndex::tileIndex(layer, ring, 0); auto itr = hgpar_->tileInfoMap_.find(indx); bool ok = (itr == hgpar_->tileInfoMap_.end()) ? false : HGCalTileIndex::tileExist(itr->second.hex, zside, phi); return ok; } + std::pair tileRings(int layer) const { + if (mode_ == HGCalGeometryMode::TrapezoidFile) { + int ll = layer - hgpar_->firstLayer_; + if (ll >= 0 && ll < static_cast(hgpar_->tileRingRange_.size())) + return hgpar_->tileRingRange_[ll]; + } + return std::make_pair(0, 0); + } int tileSiPM(int sipm) const { return ((sipm > 0) ? HGCalTypes::SiPMSmall : HGCalTypes::SiPMLarge); } bool tileTrapezoid() const { return ((mode_ == HGCalGeometryMode::Trapezoid) || (mode_ == HGCalGeometryMode::TrapezoidFile) || @@ -212,16 +221,17 @@ class HGCalDDDConstants { int wafers() const; int wafers(int layer, int type) const; int waferToCopy(int wafer) const { - return ((wafer >= 0) && (wafer < (int)(hgpar_->waferCopy_.size()))) ? hgpar_->waferCopy_[wafer] - : (int)(hgpar_->waferCopy_.size()); + return ((wafer >= 0) && (wafer < static_cast(hgpar_->waferCopy_.size()))) + ? hgpar_->waferCopy_[wafer] + : static_cast(hgpar_->waferCopy_.size()); } // wafer transverse thickness classification (2 = coarse, 1 = fine) int waferTypeT(int wafer) const { - return ((wafer >= 0) && (wafer < (int)(hgpar_->waferTypeT_.size()))) ? hgpar_->waferTypeT_[wafer] : 0; + return ((wafer >= 0) && (wafer < static_cast(hgpar_->waferTypeT_.size()))) ? hgpar_->waferTypeT_[wafer] : 0; } // wafer longitudinal thickness classification (1 = 100um, 2 = 200um, 3=300um) int waferTypeL(int wafer) const { - return ((wafer >= 0) && (wafer < (int)(hgpar_->waferTypeL_.size()))) ? hgpar_->waferTypeL_[wafer] : 0; + return ((wafer >= 0) && (wafer < static_cast(hgpar_->waferTypeL_.size()))) ? hgpar_->waferTypeL_[wafer] : 0; } int waferType(DetId const& id, bool fromFile = false) const; int waferType(int layer, int waferU, int waferV, bool fromFile = false) const; diff --git a/Geometry/HGCalCommonData/interface/HGCalGeomParameters.h b/Geometry/HGCalCommonData/interface/HGCalGeomParameters.h index 154407eeadad9..717d850530932 100644 --- a/Geometry/HGCalCommonData/interface/HGCalGeomParameters.h +++ b/Geometry/HGCalCommonData/interface/HGCalGeomParameters.h @@ -99,6 +99,10 @@ class HGCalGeomParameters { : half(h), wafer(w), xyz(std::move(p)) {} }; + constexpr static int siliconFileEE = 2; + constexpr static int siliconFileHE = 3; + constexpr static int scintillatorFile = 4; + private: void loadGeometryHexagon(const std::map& layers, std::vector& trforms, @@ -143,9 +147,6 @@ class HGCalGeomParameters { void resetZero(std::vector&); constexpr static double tan30deg_ = 0.5773502693; - constexpr static int siliconFileEE = 2; - constexpr static int siliconFileHE = 3; - constexpr static int scintillatorFile = 4; HGCalGeomTools geomTools_; const double sqrt3_; double waferSize_; diff --git a/Geometry/HGCalCommonData/src/HGCalDDDConstants.cc b/Geometry/HGCalCommonData/src/HGCalDDDConstants.cc index 10286d481690c..aefdcf260eee3 100644 --- a/Geometry/HGCalCommonData/src/HGCalDDDConstants.cc +++ b/Geometry/HGCalCommonData/src/HGCalDDDConstants.cc @@ -7,6 +7,7 @@ #include "DataFormats/ForwardDetId/interface/HGCalDetId.h" #include "FWCore/MessageLogger/interface/MessageLogger.h" #include "FWCore/Utilities/interface/Exception.h" +#include "Geometry/HGCalCommonData/interface/HGCalGeomParameters.h" #include "Geometry/HGCalCommonData/interface/HGCalGeomTools.h" #include "Geometry/HGCalCommonData/interface/HGCalGeometryMode.h" #include "Geometry/HGCalCommonData/interface/HGCalTypes.h" @@ -15,6 +16,7 @@ #include "Geometry/HGCalCommonData/interface/HGCalWaferType.h" #include +#include #include #include #include @@ -588,7 +590,7 @@ bool HGCalDDDConstants::isValidTrap(int layer, int irad, int iphi) const { const auto& indx = getIndex(layer, true); if (indx.first < 0) return false; - return ((irad >= hgpar_->iradMinBH_[indx.first]) && (irad <= hgpar_->iradMaxBH_[indx.first]) && (iphi > 0) && + return ((irad >= hgpar_->iradMinBH_[indx.first]) && (irad <= (hgpar_->iradMaxBH_[indx.first] + 1)) && (iphi > 0) && (iphi <= hgpar_->scintCells(layer))); } @@ -949,11 +951,14 @@ double HGCalDDDConstants::mouseBite(bool reco) const { } int HGCalDDDConstants::numberCells(bool reco) const { - int cells(0); - unsigned int nlayer = (reco) ? hgpar_->depth_.size() : hgpar_->layer_.size(); - for (unsigned k = 0; k < nlayer; ++k) { - std::vector ncells = numberCells(((reco) ? hgpar_->depth_[k] : hgpar_->layer_[k]), reco); - cells = std::accumulate(ncells.begin(), ncells.end(), cells); + int cells = + (tileTrapezoid() && (hgpar_->waferMaskMode_ == HGCalGeomParameters::scintillatorFile)) ? tileCount(0, -1) : 0; + if (cells == 0) { + unsigned int nlayer = (reco) ? hgpar_->depth_.size() : hgpar_->layer_.size(); + for (unsigned k = 0; k < nlayer; ++k) { + std::vector ncells = numberCells(((reco) ? hgpar_->depth_[k] : hgpar_->layer_[k]), reco); + cells = std::accumulate(ncells.begin(), ncells.end(), cells); + } } return cells; } @@ -1111,6 +1116,32 @@ std::pair HGCalDDDConstants::simToReco(int cell, int lay, int mod, boo } } +int HGCalDDDConstants::tileCount(int layer, int ring) const { + int laymin(layer), laymax(layer), ringmin(ring), ringmax(ring), kount(0); + if (layer == 0) { + laymin = hgpar_->firstLayer_; + laymax = lastLayer(true); + } + for (int lay = laymin; lay <= laymax; ++lay) { + if (ring < 0) { + int ll = lay - hgpar_->firstLayer_; + ringmin = hgpar_->tileRingRange_[ll].first; + ringmax = hgpar_->tileRingRange_[ll].second; + } + for (int rin = ringmin; rin <= ringmax; ++rin) { + int indx = HGCalTileIndex::tileIndex(lay, rin + 1, 0); + auto itr = hgpar_->tileInfoMap_.find(indx); + if (itr != hgpar_->tileInfoMap_.end()) { + for (int k = 0; k < 4; ++k) { + std::bitset<24> b(itr->second.hex[k]); + kount += b.count(); + } + } + } + } + return (3 * kount); +} + int HGCalDDDConstants::waferFromCopy(int copy) const { const int ncopies = hgpar_->waferCopy_.size(); int wafer(ncopies); diff --git a/Geometry/HGCalCommonData/src/HGCalGeomParameters.cc b/Geometry/HGCalCommonData/src/HGCalGeomParameters.cc index 75019519e3495..08353a393909c 100644 --- a/Geometry/HGCalCommonData/src/HGCalGeomParameters.cc +++ b/Geometry/HGCalCommonData/src/HGCalGeomParameters.cc @@ -2109,7 +2109,7 @@ void HGCalGeomParameters::loadCellTrapezoid(HGCalParameters& php) { // Minimum and maximum radius index for each layer for (unsigned int k = 0; k < php.zLayerHex_.size(); ++k) { php.iradMinBH_.emplace_back(1 + php.tileRingRange_[k].first); - php.iradMaxBH_.emplace_back(php.tileRingRange_[k].second); + php.iradMaxBH_.emplace_back(1 + php.tileRingRange_[k].second); #ifdef EDM_ML_DEBUG int kk = php.scintType(php.firstLayer_ + (int)(k)); edm::LogVerbatim("HGCalGeom") << "New Layer " << k << " Type " << kk << " Low edge " << php.iradMinBH_.back() diff --git a/Geometry/HGCalGeometry/src/HGCalGeometryLoader.cc b/Geometry/HGCalGeometry/src/HGCalGeometryLoader.cc index 2357f1fd8f888..16719ac17ac4c 100644 --- a/Geometry/HGCalGeometry/src/HGCalGeometryLoader.cc +++ b/Geometry/HGCalGeometry/src/HGCalGeometryLoader.cc @@ -182,7 +182,7 @@ HGCalGeometry* HGCalGeometryLoader::build(const HGCalTopology& topology) { geom->sortDetIds(); if (counter != numberExpected) { - if (test) { + if (topology.tileTrapezoid()) { edm::LogVerbatim("HGCalGeom") << "Inconsistent # of cells: expected " << numberExpected << ":" << numberOfCells << " , inited " << counter; } else { diff --git a/Geometry/HGCalGeometry/test/HGCalValidScintTest.cc b/Geometry/HGCalGeometry/test/HGCalValidScintTest.cc new file mode 100644 index 0000000000000..4ee801f2243a8 --- /dev/null +++ b/Geometry/HGCalGeometry/test/HGCalValidScintTest.cc @@ -0,0 +1,108 @@ +#include +#include +#include +#include +#include + +#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/Framework/interface/MakerMacros.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/ForwardSubdetector.h" +#include "DataFormats/ForwardDetId/interface/HGCalDetId.h" +#include "DataFormats/ForwardDetId/interface/HGCScintillatorDetId.h" +#include "DataFormats/ForwardDetId/interface/HGCSiliconDetId.h" + +class HGCalValidScintTest : public edm::one::EDAnalyzer<> { +public: + explicit HGCalValidScintTest(const edm::ParameterSet&); + ~HGCalValidScintTest() override = default; + + static void fillDescriptions(edm::ConfigurationDescriptions& descriptions); + void analyze(edm::Event const& iEvent, edm::EventSetup const&) override; + +private: + const edm::ESGetToken geomToken_; + struct layerInfo { + int ringMin, ringMax; + double rMin, rMax; + layerInfo(int minR = 100, double rMn = 0, int maxR = 0, double rMx = 0) + : ringMin(minR), ringMax(maxR), rMin(rMn), rMax(rMx){}; + }; +}; + +HGCalValidScintTest::HGCalValidScintTest(const edm::ParameterSet& iC) + : geomToken_(esConsumes(edm::ESInputTag{"", "HGCalHEScintillatorSensitive"})) {} + +void HGCalValidScintTest::fillDescriptions(edm::ConfigurationDescriptions& descriptions) { + edm::ParameterSetDescription desc; + descriptions.add("hgcalValidScintTest", desc); +} + +void HGCalValidScintTest::analyze(const edm::Event&, const edm::EventSetup& iSetup) { + const auto& geom = &iSetup.getData(geomToken_); + DetId::Detector det = DetId::HGCalHSc; + edm::LogVerbatim("HGCalGeom") << "Perform test for HGCalHEScintillatorSensitive Detector " << det << " Mode " + << geom->topology().dddConstants().geomMode(); + + int firstLayer = geom->topology().dddConstants().firstLayer(); + int lastLayer = geom->topology().dddConstants().lastLayer(true); + const std::vector& ids = geom->getValidDetIds(); + edm::LogVerbatim("HGCalGeom") << "doTest: " << ids.size() << " valid ids for " << geom->cellElement(); + + int zside(1); + std::map layerMap; + for (int layer = firstLayer; layer <= lastLayer; ++layer) { + std::vector > done; + for (auto const& id : ids) { + HGCScintillatorDetId hid(id); + if ((hid.zside() != zside) || (hid.layer() != layer)) + continue; + std::pair ring = std::make_pair(hid.ring(), 0); + if (std::find(done.begin(), done.end(), ring) != done.end()) + continue; + done.emplace_back(ring); + edm::LogVerbatim("HGCalGeom") << "Corners for " << hid; + + const auto cor = geom->getNewCorners(id); + std::ostringstream st1; + for (unsigned int k = 0; k < cor.size(); ++k) + st1 << " [" << k << "] " << std::setprecision(4) << cor[k]; + edm::LogVerbatim("HGCalGeom") << st1.str(); + + double r = cor[0].perp(); + auto itr = layerMap.find(layer); + if (itr == layerMap.end()) { + layerInfo info; + info.ringMin = info.ringMax = ring.first; + info.rMin = info.rMax = r; + layerMap[layer] = info; + } else { + layerInfo info = itr->second; + if (info.ringMin > ring.first) { + info.ringMin = ring.first; + info.rMin = r; + } else if (info.ringMax < ring.first) { + info.ringMax = ring.first; + info.rMax = r; + } + layerMap[layer] = info; + } + } + } + edm::LogVerbatim("HGCalGeom") << "\n\nSummary of " << layerMap.size() << " Scintillator layers"; + for (auto itr = layerMap.begin(); itr != layerMap.end(); ++itr) + edm::LogVerbatim("HGCalGeom") << "Layer " << itr->first << " lowest Ring " << (itr->second).ringMin << ":" + << (itr->second).rMin << " largest Ring " << (itr->second).ringMax << ":" + << (itr->second).rMax; +} + +//define this as a plug-in +DEFINE_FWK_MODULE(HGCalValidScintTest); diff --git a/Geometry/HGCalGeometry/test/python/testHGCalValidScintTest_cfg.py b/Geometry/HGCalGeometry/test/python/testHGCalValidScintTest_cfg.py new file mode 100644 index 0000000000000..50b91187c4a1e --- /dev/null +++ b/Geometry/HGCalGeometry/test/python/testHGCalValidScintTest_cfg.py @@ -0,0 +1,39 @@ +import FWCore.ParameterSet.Config as cms + +from Configuration.Eras.Era_Phase2C11_cff import Phase2C11 +process = cms.Process('PROD',Phase2C11) + +process.load("SimGeneral.HepPDTESSource.pdt_cfi") +process.load("Configuration.Geometry.GeometryExtended2026D86Reco_cff") +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.load("Geometry.HGCalGeometry.hgcalValidScintTest_cfi") + +process.p1 = cms.Path(process.generator*process.hgcalValidScintTest)