diff --git a/HLTrigger/Configuration/python/HLT_75e33/eventsetup/hltLSTGeometry_cfi.py b/HLTrigger/Configuration/python/HLT_75e33/eventsetup/hltLSTGeometry_cfi.py new file mode 100644 index 0000000000000..b6a964749abf9 --- /dev/null +++ b/HLTrigger/Configuration/python/HLT_75e33/eventsetup/hltLSTGeometry_cfi.py @@ -0,0 +1,5 @@ +import FWCore.ParameterSet.Config as cms + +hltLSTGeometry = cms.ESProducer('LSTGeometryESProducer', + ptCut = cms.double(0.8), +) diff --git a/HLTrigger/Configuration/python/HLT_75e33_cff.py b/HLTrigger/Configuration/python/HLT_75e33_cff.py index 0abb71a072bc8..d57129a7c56f9 100644 --- a/HLTrigger/Configuration/python/HLT_75e33_cff.py +++ b/HLTrigger/Configuration/python/HLT_75e33_cff.py @@ -95,6 +95,7 @@ fragment.load("HLTrigger/Configuration/HLT_75e33/eventsetup/hltESPPixelCPEFastParams_cfi") +fragment.load("HLTrigger/Configuration/HLT_75e33/eventsetup/hltLSTGeometry_cfi") fragment.load("HLTrigger/Configuration/HLT_75e33/eventsetup/hltESPModulesDevLST_cfi") fragment.load("HLTrigger/Configuration/HLT_75e33/eventsetup/hltESPTTRHBuilderWithoutRefit_cfi") diff --git a/HLTrigger/Configuration/python/HLT_75e33_timing_cff.py b/HLTrigger/Configuration/python/HLT_75e33_timing_cff.py index 3bc1b1f1ffff5..c98bf9b4715df 100644 --- a/HLTrigger/Configuration/python/HLT_75e33_timing_cff.py +++ b/HLTrigger/Configuration/python/HLT_75e33_timing_cff.py @@ -100,6 +100,7 @@ fragment.load("HLTrigger/Configuration/HLT_75e33/eventsetup/hltESPPixelCPEFastParams_cfi") +fragment.load("HLTrigger/Configuration/HLT_75e33/eventsetup/hltLSTGeometry_cfi") fragment.load("HLTrigger/Configuration/HLT_75e33/eventsetup/hltESPModulesDevLST_cfi") fragment.load("HLTrigger/Configuration/HLT_75e33/eventsetup/hltESPTTRHBuilderWithoutRefit_cfi") diff --git a/HLTrigger/Configuration/python/HLT_NGTScouting_cff.py b/HLTrigger/Configuration/python/HLT_NGTScouting_cff.py index 2702107d0cad7..ec1d2171cd7e2 100644 --- a/HLTrigger/Configuration/python/HLT_NGTScouting_cff.py +++ b/HLTrigger/Configuration/python/HLT_NGTScouting_cff.py @@ -94,6 +94,7 @@ fragment.load("HLTrigger/Configuration/HLT_75e33/eventsetup/hltESPKFTrajectorySmootherForL2Muon_cfi") fragment.load("HLTrigger/Configuration/HLT_75e33/eventsetup/trackdnn_source_cfi") fragment.load("HLTrigger/Configuration/HLT_75e33/eventsetup/hltESPPixelCPEFastParams_cfi") +fragment.load("HLTrigger/Configuration/HLT_75e33/eventsetup/hltLSTGeometry_cfi") fragment.load("HLTrigger/Configuration/HLT_75e33/eventsetup/hltESPModulesDevLST_cfi") fragment.load("HLTrigger/Configuration/HLT_75e33/eventsetup/hltESPTTRHBuilderWithoutRefit_cfi") fragment.load("HLTrigger/Configuration/HLT_75e33/eventsetup/hltESPMkFit_cfi") diff --git a/RecoTracker/IterativeTracking/python/HighPtTripletStep_cff.py b/RecoTracker/IterativeTracking/python/HighPtTripletStep_cff.py index 753e5562b76cc..9df0bbcde8ca8 100644 --- a/RecoTracker/IterativeTracking/python/HighPtTripletStep_cff.py +++ b/RecoTracker/IterativeTracking/python/HighPtTripletStep_cff.py @@ -427,11 +427,12 @@ _HighPtTripletStepTask_LST = HighPtTripletStepTask.copy() from RecoLocalTracker.Phase2TrackerRecHits.Phase2TrackerRecHits_cfi import siPhase2RecHits +from RecoTracker.LSTGeometry.lstGeometryESProducer_cfi import lstGeometryESProducer from RecoTracker.LST.lstSeedTracks_cff import lstInitialStepSeedTracks,lstHighPtTripletStepSeedTracks from RecoTracker.LST.lstInputProducer_cfi import lstInputProducer from RecoTracker.LST.lstProducerTask_cff import * -_HighPtTripletStepTask_LST.add(siPhase2RecHits, lstInitialStepSeedTracks, lstHighPtTripletStepSeedTracks, lstInputProducer, +_HighPtTripletStepTask_LST.add(siPhase2RecHits, lstInitialStepSeedTracks, lstHighPtTripletStepSeedTracks, lstGeometryESProducer, lstInputProducer, lstProducerTask, highPtTripletStepLSTpTracks, highPtTripletStepLSTT4T5Tracks, highPtTripletStepSelectorLSTT4T5) (trackingPhase2PU140 & trackingLST).toReplaceWith(HighPtTripletStepTask, _HighPtTripletStepTask_LST) diff --git a/RecoTracker/LST/plugins/alpaka/LSTModulesDevESProducer.cc b/RecoTracker/LST/plugins/alpaka/LSTModulesDevESProducer.cc index 7152da9ed13c7..e9327eaf2e5f8 100644 --- a/RecoTracker/LST/plugins/alpaka/LSTModulesDevESProducer.cc +++ b/RecoTracker/LST/plugins/alpaka/LSTModulesDevESProducer.cc @@ -1,5 +1,7 @@ // LST includes #include "RecoTracker/LSTCore/interface/alpaka/LST.h" +#include "RecoTracker/LSTGeometry/interface/Common.h" +#include "RecoTracker/LSTGeometry/interface/Geometry.h" #include "FWCore/ParameterSet/interface/ParameterSet.h" @@ -14,22 +16,27 @@ namespace ALPAKA_ACCELERATOR_NAMESPACE { class LSTModulesDevESProducer : public ESProducer { private: - std::string ptCutLabel_; + double ptCut_; + edm::ESGetToken lstGeoToken_; public: LSTModulesDevESProducer(edm::ParameterSet const& iConfig) - : ESProducer(iConfig), ptCutLabel_(iConfig.getParameter("ptCutLabel")) { - setWhatProduced(this, ptCutLabel_); + : ESProducer(iConfig), ptCut_(iConfig.getParameter("ptCut")) { + std::string ptCutStr = lst::floatToStr(ptCut_, 1); + + auto cc = setWhatProduced(this, ptCutStr); + lstGeoToken_ = cc.consumes(edm::ESInputTag("", ptCutStr)); } static void fillDescriptions(edm::ConfigurationDescriptions& descriptions) { edm::ParameterSetDescription desc; - desc.add("ptCutLabel", "0.8"); + desc.add("ptCut", 0.8); descriptions.addWithDefaultLabel(desc); } std::unique_ptr> produce(TrackerRecoGeometryRecord const& iRecord) { - return lst::loadAndFillESHost(ptCutLabel_); + const auto& lstg = iRecord.get(lstGeoToken_); + return lst::fillESDataHost(lstg); } }; diff --git a/RecoTracker/LST/plugins/alpaka/LSTProducer.cc b/RecoTracker/LST/plugins/alpaka/LSTProducer.cc index 496c1ae1182b0..ddb1c0916c8ee 100644 --- a/RecoTracker/LST/plugins/alpaka/LSTProducer.cc +++ b/RecoTracker/LST/plugins/alpaka/LSTProducer.cc @@ -1,6 +1,7 @@ #include #include "RecoTracker/LSTCore/interface/alpaka/LST.h" +#include "RecoTracker/LSTGeometry/interface/Common.h" #include "FWCore/MessageLogger/interface/MessageLogger.h" #include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h" @@ -25,13 +26,14 @@ namespace ALPAKA_ACCELERATOR_NAMESPACE { public: LSTProducer(edm::ParameterSet const& config) : EDProducer(config), - lstInputToken_{consumes(config.getParameter("lstInput"))}, - lstESToken_{esConsumes(edm::ESInputTag("", config.getParameter("ptCutLabel")))}, verbose_(config.getParameter("verbose")), ptCut_(config.getParameter("ptCut")), + ptCutStr_(lst::floatToStr(ptCut_, 1)), clustSizeCut_(static_cast(config.getParameter("clustSizeCut"))), nopLSDupClean_(config.getParameter("nopLSDupClean")), tcpLSTriplets_(config.getParameter("tcpLSTriplets")), + lstInputToken_{consumes(config.getParameter("lstInput"))}, + lstESToken_{esConsumes(edm::ESInputTag("", ptCutStr_))}, lstOutputToken_{produces()} {} void produce(edm::StreamID sid, device::Event& iEvent, const device::EventSetup& iSetup) const override { @@ -42,7 +44,7 @@ namespace ALPAKA_ACCELERATOR_NAMESPACE { lst.run(iEvent.queue(), verbose_, - static_cast(ptCut_), + ptCut_, clustSizeCut_, &lstESDeviceData, &lstInputDC, @@ -60,20 +62,20 @@ namespace ALPAKA_ACCELERATOR_NAMESPACE { desc.add("verbose", false); desc.add("ptCut", 0.8); desc.add("clustSizeCut", 16); - desc.add("ptCutLabel", "0.8"); desc.add("nopLSDupClean", false); desc.add("tcpLSTriplets", false); descriptions.addWithDefaultLabel(desc); } private: - const device::EDGetToken lstInputToken_; - const device::ESGetToken, TrackerRecoGeometryRecord> lstESToken_; const bool verbose_; const double ptCut_; + const std::string ptCutStr_; const uint16_t clustSizeCut_; const bool nopLSDupClean_; const bool tcpLSTriplets_; + const device::EDGetToken lstInputToken_; + const device::ESGetToken, TrackerRecoGeometryRecord> lstESToken_; const device::EDPutToken lstOutputToken_; }; diff --git a/RecoTracker/LST/python/lstProducerTask_cff.py b/RecoTracker/LST/python/lstProducerTask_cff.py index 588b354788635..3de184640b97c 100644 --- a/RecoTracker/LST/python/lstProducerTask_cff.py +++ b/RecoTracker/LST/python/lstProducerTask_cff.py @@ -4,4 +4,8 @@ from RecoTracker.LST.lstModulesDevESProducer_cfi import lstModulesDevESProducer -lstProducerTask = cms.Task(lstModulesDevESProducer, lstProducer) +from RecoTracker.LST.lstInputProducer_cfi import lstInputProducer + +from RecoTracker.LSTGeometry.lstGeometryESProducer_cfi import lstGeometryESProducer + +lstProducerTask = cms.Task(lstGeometryESProducer, lstModulesDevESProducer, lstInputProducer, lstProducer) diff --git a/RecoTracker/LSTCore/BuildFile.xml b/RecoTracker/LSTCore/BuildFile.xml index bdc37903ba99b..2225bd153aa94 100644 --- a/RecoTracker/LSTCore/BuildFile.xml +++ b/RecoTracker/LSTCore/BuildFile.xml @@ -6,6 +6,7 @@ + diff --git a/RecoTracker/LSTCore/interface/Common.h b/RecoTracker/LSTCore/interface/Common.h index b5dece7656b2d..b675caa29b221 100644 --- a/RecoTracker/LSTCore/interface/Common.h +++ b/RecoTracker/LSTCore/interface/Common.h @@ -4,6 +4,8 @@ #include "DataFormats/Common/interface/StdArray.h" #include "HeterogeneousCore/AlpakaInterface/interface/config.h" +#include "RecoTracker/LSTGeometry/interface/Common.h" + #if defined(FP16_Base) #if defined ALPAKA_ACC_GPU_CUDA_ENABLED #include @@ -42,6 +44,9 @@ namespace lst { constexpr uint16_t kTCEmptyLowerModule = 0xFFFF; // Sentinel for empty lowerModule index constexpr unsigned int kTCEmptyHitIdx = 0xFFFFFFFF; // Sentinel for empty hit slots + constexpr float kB = lstgeometry::kB; + constexpr float kC = lstgeometry::kC; + // Half precision wrapper functions. #if defined(FP16_Base) #define __F2H __float2half diff --git a/RecoTracker/LSTCore/interface/EndcapGeometry.h b/RecoTracker/LSTCore/interface/EndcapGeometry.h index b8c44c14fb143..59ce62d24eccf 100644 --- a/RecoTracker/LSTCore/interface/EndcapGeometry.h +++ b/RecoTracker/LSTCore/interface/EndcapGeometry.h @@ -1,6 +1,9 @@ #ifndef RecoTracker_LSTCore_interface_EndcapGeometry_h #define RecoTracker_LSTCore_interface_EndcapGeometry_h +#include "RecoTracker/LSTGeometry/interface/Sensor.h" +#include "RecoTracker/LSTGeometry/interface/Slope.h" + #include #include #include @@ -18,9 +21,11 @@ namespace lst { unsigned int nEndCapMap; EndcapGeometry() = default; - EndcapGeometry(std::string const& filename); + EndcapGeometry(std::string const&); + EndcapGeometry(lstgeometry::Slopes const&, lstgeometry::Sensors const& sensors); void load(std::string const&); + void load(lstgeometry::Slopes const&, lstgeometry::Sensors const& sensors); void fillGeoMapArraysExplicit(); float getdxdy_slope(unsigned int detid) const; }; diff --git a/RecoTracker/LSTCore/interface/LSTESData.h b/RecoTracker/LSTCore/interface/LSTESData.h index 5131d013a4279..2c0a0458b3530 100644 --- a/RecoTracker/LSTCore/interface/LSTESData.h +++ b/RecoTracker/LSTCore/interface/LSTESData.h @@ -5,6 +5,7 @@ #include "RecoTracker/LSTCore/interface/EndcapGeometryDevHostCollection.h" #include "RecoTracker/LSTCore/interface/ModulesHostCollection.h" #include "RecoTracker/LSTCore/interface/PixelMap.h" +#include "RecoTracker/LSTGeometry/interface/Geometry.h" #include "HeterogeneousCore/AlpakaInterface/interface/CopyToDevice.h" @@ -40,7 +41,8 @@ namespace lst { pixelMapping(pixelMappingIn) {} }; - std::unique_ptr> loadAndFillESHost(std::string& ptCutLabel); + std::unique_ptr> loadAndFillESDataHost(std::string& ptCutLabel); + std::unique_ptr> fillESDataHost(lstgeometry::Geometry const& lstg); } // namespace lst diff --git a/RecoTracker/LSTCore/interface/ModuleConnectionMap.h b/RecoTracker/LSTCore/interface/ModuleConnectionMap.h index 63c3496523c0d..ecc04f2ad5d0b 100644 --- a/RecoTracker/LSTCore/interface/ModuleConnectionMap.h +++ b/RecoTracker/LSTCore/interface/ModuleConnectionMap.h @@ -12,10 +12,12 @@ namespace lst { std::map> moduleConnections_; public: - ModuleConnectionMap(); - ModuleConnectionMap(std::string const& filename); + ModuleConnectionMap() = default; + ModuleConnectionMap(std::string const&); + ModuleConnectionMap(std::map> const&); void load(std::string const&); + void load(std::map> const&); void add(std::string const&); void print(); diff --git a/RecoTracker/LSTCore/interface/TiltedGeometry.h b/RecoTracker/LSTCore/interface/TiltedGeometry.h index 7a17106195522..b5ce91de0d08c 100644 --- a/RecoTracker/LSTCore/interface/TiltedGeometry.h +++ b/RecoTracker/LSTCore/interface/TiltedGeometry.h @@ -1,6 +1,8 @@ #ifndef RecoTracker_LSTCore_interface_TiltedGeometry_h #define RecoTracker_LSTCore_interface_TiltedGeometry_h +#include "RecoTracker/LSTGeometry/interface/Slope.h" + #include #include #include @@ -13,9 +15,11 @@ namespace lst { public: TiltedGeometry() = default; - TiltedGeometry(std::string const& filename); + TiltedGeometry(std::string const&); + TiltedGeometry(lstgeometry::Slopes const&); void load(std::string const&); + void load(lstgeometry::Slopes const&); float getDrDz(unsigned int detid) const; float getDxDy(unsigned int detid) const; diff --git a/RecoTracker/LSTCore/interface/alpaka/Common.h b/RecoTracker/LSTCore/interface/alpaka/Common.h index a380fc3888583..83a5619974c0c 100644 --- a/RecoTracker/LSTCore/interface/alpaka/Common.h +++ b/RecoTracker/LSTCore/interface/alpaka/Common.h @@ -32,8 +32,8 @@ namespace ALPAKA_ACCELERATOR_NAMESPACE::lst { 25.007152356, 37.2186993757, 52.3104270826, 68.6658656666, 85.9770373007, 108.301772384}; HOST_DEVICE_CONSTANT float kMiniRminMeanEndcap[5] = { 130.992832231, 154.813883559, 185.352604327, 221.635123002, 265.022076742}; - HOST_DEVICE_CONSTANT float k2Rinv1GeVf = (2.99792458e-3 * 3.8) / 2; - HOST_DEVICE_CONSTANT float kR1GeVf = 1. / (2.99792458e-3 * 3.8); + HOST_DEVICE_CONSTANT float k2Rinv1GeVf = (kC * kB) / 2; + HOST_DEVICE_CONSTANT float kR1GeVf = 1. / (kC * kB); HOST_DEVICE_CONSTANT float kSinAlphaMax = 0.95; HOST_DEVICE_CONSTANT float kDeltaZLum = 15.0; HOST_DEVICE_CONSTANT float kPixelPSZpitch = 0.15; @@ -44,8 +44,6 @@ namespace ALPAKA_ACCELERATOR_NAMESPACE::lst { HOST_DEVICE_CONSTANT float kWidthPS = 0.01; HOST_DEVICE_CONSTANT float kPt_betaMax = 7.0; HOST_DEVICE_CONSTANT int kNTripletThreshold = 1000; - // To be updated with std::numeric_limits::infinity() in the code and data files - HOST_DEVICE_CONSTANT float kVerticalModuleSlope = 123456789.0; HOST_DEVICE_CONSTANT int kLogicalOTLayers = 11; // logical OT layers are 1..11 HOST_DEVICE_CONSTANT float kMiniDeltaTilted[3] = {0.26f, 0.26f, 0.26f}; diff --git a/RecoTracker/LSTCore/src/EndcapGeometry.cc b/RecoTracker/LSTCore/src/EndcapGeometry.cc index 17e72379bb2ec..5d50194a8c04f 100644 --- a/RecoTracker/LSTCore/src/EndcapGeometry.cc +++ b/RecoTracker/LSTCore/src/EndcapGeometry.cc @@ -2,11 +2,16 @@ #include #include +#include #include #include lst::EndcapGeometry::EndcapGeometry(std::string const& filename) { load(filename); } +lst::EndcapGeometry::EndcapGeometry(lstgeometry::Slopes const& slopes, lstgeometry::Sensors const& sensors) { + load(slopes, sensors); +} + void lst::EndcapGeometry::load(std::string const& filename) { dxdy_slope_.clear(); centroid_phis_.clear(); @@ -26,7 +31,9 @@ void lst::EndcapGeometry::load(std::string const& filename) { ifile.read(reinterpret_cast(¢roid_phi), sizeof(centroid_phi)); if (ifile) { - dxdy_slope_[detid] = dxdy_slope; + // TODO: This is needed for now in case old geometry files are used. Remove once deemed unnecessary. + constexpr float kLegacyVerticalSlope = 123456789.0f; + dxdy_slope_[detid] = (dxdy_slope == kLegacyVerticalSlope) ? std::numeric_limits::infinity() : dxdy_slope; centroid_phis_[detid] = centroid_phi; } else { // End of file or read failed @@ -39,6 +46,18 @@ void lst::EndcapGeometry::load(std::string const& filename) { fillGeoMapArraysExplicit(); } +void lst::EndcapGeometry::load(lstgeometry::Slopes const& slopes, lstgeometry::Sensors const& sensors) { + dxdy_slope_.clear(); + centroid_phis_.clear(); + + for (const auto& [detId, slope] : slopes) { + dxdy_slope_[detId] = slope.dxdy; + centroid_phis_[detId] = sensors.at(detId).centerPhi; + } + + fillGeoMapArraysExplicit(); +} + void lst::EndcapGeometry::fillGeoMapArraysExplicit() { nEndCapMap = centroid_phis_.size(); diff --git a/RecoTracker/LSTCore/src/LSTESData.cc b/RecoTracker/LSTCore/src/LSTESData.cc index 93cf76be6312f..2a840418bc276 100644 --- a/RecoTracker/LSTCore/src/LSTESData.cc +++ b/RecoTracker/LSTCore/src/LSTESData.cc @@ -1,3 +1,5 @@ +#include + #include "RecoTracker/LSTCore/interface/LSTESData.h" #include "RecoTracker/LSTCore/interface/EndcapGeometry.h" #include "RecoTracker/LSTCore/interface/ModuleConnectionMap.h" @@ -6,8 +8,6 @@ #include "ModuleMethods.h" -#include - namespace { std::string geometryDataDir() { std::string path_str, path; @@ -79,7 +79,7 @@ namespace { } } // namespace -std::unique_ptr> lst::loadAndFillESHost(std::string& ptCutLabel) { +std::unique_ptr> lst::loadAndFillESDataHost(std::string& ptCutLabel) { uint16_t nModules; uint16_t nLowerModules; unsigned int nPixels; @@ -119,3 +119,83 @@ std::unique_ptr> lst::loadAndFillESHost(s std::move(endcapGeometryDev), pixelMappingPtr); } + +std::unique_ptr> lst::fillESDataHost(lstgeometry::Geometry const& lstg) { + uint16_t nModules; + uint16_t nLowerModules; + unsigned int nPixels; + MapPLStoLayer pLStoLayer; + EndcapGeometry endcapGeometry; + TiltedGeometry tiltedGeometry; + PixelMap pixelMapping; + ModuleConnectionMap moduleConnectionMap; + + endcapGeometry.load(lstg.endcap_slopes, lstg.sensors); + auto endcapGeometryDev = + std::make_shared(cms::alpakatools::host(), endcapGeometry.nEndCapMap); + std::memcpy(endcapGeometryDev->view().geoMapDetId().data(), + endcapGeometry.geoMapDetId_buf.data(), + endcapGeometry.nEndCapMap * sizeof(unsigned int)); + std::memcpy(endcapGeometryDev->view().geoMapPhi().data(), + endcapGeometry.geoMapPhi_buf.data(), + endcapGeometry.nEndCapMap * sizeof(float)); + + tiltedGeometry.load(lstg.barrel_slopes); + + std::map> final_modulemap; + for (auto const& [detId, connections] : lstg.module_map) { + if (connections.empty()) + continue; + final_modulemap[detId] = connections; + } + moduleConnectionMap.load(final_modulemap); + + for (auto& [layersubdetcharge, map] : lstg.pixel_map) { + auto& [layer, subdet, charge] = layersubdetcharge; + + std::map> final_pixelmap; + for (unsigned int isuperbin = 0; isuperbin < map.size(); isuperbin++) { + auto const& vec = map.at(isuperbin); + if (vec.empty()) + continue; + final_pixelmap[isuperbin] = vec; + } + + int charge_index = (charge == 0 ? 0 : (charge > 0 ? 1 : 2)); + int layer_subdet_index = layer - 1 + (subdet == Endcap ? 2 : 0); + pLStoLayer[charge_index][layer_subdet_index] = lst::ModuleConnectionMap(final_pixelmap); + } + + ModuleMetaData mmd; + unsigned int counter = 0; + for (auto const& [detId, sensor] : lstg.sensors) { + mmd.detIdToIndex[detId] = counter; + mmd.module_x[detId] = sensor.centerX; + mmd.module_y[detId] = sensor.centerY; + mmd.module_z[detId] = sensor.centerZ; + mmd.module_type[detId] = static_cast(sensor.moduleType); + counter++; + } + mmd.detIdToIndex[kPixelModuleId] = counter; //pixel module is the last module in the module list + counter++; + nModules = counter; + + auto modulesBuffers = constructModuleCollection(pLStoLayer, + mmd, + nModules, + nLowerModules, + nPixels, + pixelMapping, + endcapGeometry, + tiltedGeometry, + moduleConnectionMap); + + auto pixelMappingPtr = std::make_shared(std::move(pixelMapping)); + return std::make_unique>(nModules, + nLowerModules, + nPixels, + endcapGeometry.nEndCapMap, + std::move(modulesBuffers), + std::move(endcapGeometryDev), + pixelMappingPtr); +} diff --git a/RecoTracker/LSTCore/src/ModuleConnectionMap.cc b/RecoTracker/LSTCore/src/ModuleConnectionMap.cc index 0da0f4cc4ac6f..64fcb0769886f 100644 --- a/RecoTracker/LSTCore/src/ModuleConnectionMap.cc +++ b/RecoTracker/LSTCore/src/ModuleConnectionMap.cc @@ -5,10 +5,12 @@ #include #include -lst::ModuleConnectionMap::ModuleConnectionMap() {} - lst::ModuleConnectionMap::ModuleConnectionMap(std::string const& filename) { load(filename); } +lst::ModuleConnectionMap::ModuleConnectionMap(std::map> const& map) { + load(map); +} + void lst::ModuleConnectionMap::load(std::string const& filename) { moduleConnections_.clear(); @@ -53,6 +55,10 @@ void lst::ModuleConnectionMap::load(std::string const& filename) { } } +void lst::ModuleConnectionMap::load(std::map> const& map) { + moduleConnections_ = map; +} + void lst::ModuleConnectionMap::add(std::string const& filename) { std::ifstream ifile; ifile.open(filename.c_str()); diff --git a/RecoTracker/LSTCore/src/ModuleMethods.h b/RecoTracker/LSTCore/src/ModuleMethods.h index 1d5d75f18ac2b..b237333158204 100644 --- a/RecoTracker/LSTCore/src/ModuleMethods.h +++ b/RecoTracker/LSTCore/src/ModuleMethods.h @@ -16,6 +16,9 @@ #include "HeterogeneousCore/AlpakaInterface/interface/memory.h" namespace lst { + // TODO: At some point it might be good to move some of these things to LSTGeometry + // or replace the functions with values computed from standard geometry/topology methods + struct ModuleMetaData { std::map detIdToIndex; std::map module_x; @@ -138,6 +141,15 @@ namespace lst { uint16_t index = it->second; auto& connectedModules = moduleConnectionMap.getConnectedModuleDetIds(detId); nConnectedModules[index] = connectedModules.size(); + if (nConnectedModules[index] > max_connected_modules) { +#ifdef WARNINGS + printf("Warning: module %u has %u connections, exceeding max_connected_modules=%u. Truncating.\n", + detId, + nConnectedModules[index], + max_connected_modules); +#endif + nConnectedModules[index] = max_connected_modules; + } for (uint16_t i = 0; i < nConnectedModules[index]; i++) { moduleMap[index][i] = mmd.detIdToIndex.at(connectedModules[i]); } @@ -213,24 +225,21 @@ namespace lst { } } - mmd.detIdToIndex[1] = counter; //pixel module is the last module in the module list + mmd.detIdToIndex[kPixelModuleId] = counter; //pixel module is the last module in the module list counter++; nModules = counter; } - inline std::shared_ptr loadModulesFromFile(MapPLStoLayer const& pLStoLayer, - const char* moduleMetaDataFilePath, - uint16_t& nModules, - uint16_t& nLowerModules, - unsigned int& nPixels, - PixelMap& pixelMapping, - const EndcapGeometry& endcapGeometry, - const TiltedGeometry& tiltedGeometry, - const ModuleConnectionMap& moduleConnectionMap) { - ModuleMetaData mmd; - - loadCentroidsFromFile(moduleMetaDataFilePath, mmd, nModules); - + inline std::shared_ptr constructModuleCollection( + MapPLStoLayer const& pLStoLayer, + ModuleMetaData& mmd, + uint16_t& nModules, + uint16_t& nLowerModules, + unsigned int& nPixels, + PixelMap& pixelMapping, + const EndcapGeometry& endcapGeometry, + const TiltedGeometry& tiltedGeometry, + const ModuleConnectionMap& moduleConnectionMap) { // TODO: this whole section could use some refactoring auto [totalSizes, connectedModuleDetIds, @@ -382,7 +391,7 @@ namespace lst { *host_nLowerModules = nLowerModules; // Fill pixel part - pixelMapping.pixelModuleIndex = mmd.detIdToIndex.at(1); + pixelMapping.pixelModuleIndex = mmd.detIdToIndex.at(kPixelModuleId); auto modulesPixel_view = modulesHC->view().modulesPixel(); auto connectedPixels = @@ -402,5 +411,29 @@ namespace lst { return modulesHC; } + + inline std::shared_ptr loadModulesFromFile(MapPLStoLayer const& pLStoLayer, + const char* moduleMetaDataFilePath, + uint16_t& nModules, + uint16_t& nLowerModules, + unsigned int& nPixels, + PixelMap& pixelMapping, + const EndcapGeometry& endcapGeometry, + const TiltedGeometry& tiltedGeometry, + const ModuleConnectionMap& moduleConnectionMap) { + ModuleMetaData mmd; + + loadCentroidsFromFile(moduleMetaDataFilePath, mmd, nModules); + return constructModuleCollection(pLStoLayer, + mmd, + nModules, + nLowerModules, + nPixels, + pixelMapping, + endcapGeometry, + tiltedGeometry, + moduleConnectionMap); + } + } // namespace lst #endif diff --git a/RecoTracker/LSTCore/src/TiltedGeometry.cc b/RecoTracker/LSTCore/src/TiltedGeometry.cc index d65a9a4a5f7b9..6c2f10b32b35d 100644 --- a/RecoTracker/LSTCore/src/TiltedGeometry.cc +++ b/RecoTracker/LSTCore/src/TiltedGeometry.cc @@ -2,11 +2,14 @@ #include #include +#include #include #include lst::TiltedGeometry::TiltedGeometry(std::string const& filename) { load(filename); } +lst::TiltedGeometry::TiltedGeometry(lstgeometry::Slopes const& slopes) { load(slopes); } + void lst::TiltedGeometry::load(std::string const& filename) { drdzs_.clear(); dxdys_.clear(); @@ -26,8 +29,10 @@ void lst::TiltedGeometry::load(std::string const& filename) { ifile.read(reinterpret_cast(&dxdy), sizeof(dxdy)); if (ifile) { - drdzs_[detid] = drdz; - dxdys_[detid] = dxdy; + // TODO: This is needed for now in case old geometry files are used. Remove once deemed unnecessary. + constexpr float kLegacyVerticalSlope = 123456789.0f; + drdzs_[detid] = (drdz == kLegacyVerticalSlope) ? std::numeric_limits::infinity() : drdz; + dxdys_[detid] = (dxdy == kLegacyVerticalSlope) ? std::numeric_limits::infinity() : dxdy; } else { // End of file or read failed if (!ifile.eof()) { @@ -37,6 +42,16 @@ void lst::TiltedGeometry::load(std::string const& filename) { } } +void lst::TiltedGeometry::load(lstgeometry::Slopes const& slopes) { + drdzs_.clear(); + dxdys_.clear(); + + for (const auto& [detId, slope] : slopes) { + drdzs_[detId] = slope.drdz; + dxdys_[detId] = slope.dxdy; + } +} + float lst::TiltedGeometry::getDrDz(unsigned int detid) const { auto res = drdzs_.find(detid); return res == drdzs_.end() ? 0.f : res->second; diff --git a/RecoTracker/LSTCore/src/alpaka/MiniDoublet.h b/RecoTracker/LSTCore/src/alpaka/MiniDoublet.h index 5b2effc827bf0..ba7f1b8d9c596 100644 --- a/RecoTracker/LSTCore/src/alpaka/MiniDoublet.h +++ b/RecoTracker/LSTCore/src/alpaka/MiniDoublet.h @@ -307,7 +307,7 @@ namespace ALPAKA_ACCELERATOR_NAMESPACE::lst { float drprime_x, drprime_y; // drprime * {sin,cos}(atan(slope)) // Algebraic: sin(atan(slope)) = |slope|/sqrt(1+slope^2), cos(atan(slope)) = 1/sqrt(1+slope^2) const float slope = mod.slope; - if (slope != kVerticalModuleSlope && edm::isFinite(slope)) { + if (edm::isFinite(slope)) { const float inv_hypot_slope = 1.f / alpaka::math::sqrt(acc, 1.f + slope * slope); drprime_x = drprime * ((xp > 0.f) - (xp < 0.f)) * alpaka::math::abs(acc, slope) * inv_hypot_slope; drprime_y = drprime * ((yp > 0.f) - (yp < 0.f)) * inv_hypot_slope; @@ -321,7 +321,7 @@ namespace ALPAKA_ACCELERATOR_NAMESPACE::lst { // Compute the new strip hit position (handle slope = infinity and slope = 0 cases) float xn, yn; - if (slope == kVerticalModuleSlope || edm::isNotFinite(slope)) { + if (edm::isNotFinite(slope)) { xn = xa; yn = yo; } else if (slope == 0) { diff --git a/RecoTracker/LSTCore/src/alpaka/Quintuplet.h b/RecoTracker/LSTCore/src/alpaka/Quintuplet.h index 15b0aa1f7f405..70fdaef7d8ec6 100644 --- a/RecoTracker/LSTCore/src/alpaka/Quintuplet.h +++ b/RecoTracker/LSTCore/src/alpaka/Quintuplet.h @@ -667,9 +667,8 @@ namespace ALPAKA_ACCELERATOR_NAMESPACE::lst { // Computing sigmas is a very tricky affair // if the module is tilted or endcap, we need to use the slopes properly! - absArctanSlope = ((slopes[i] != kVerticalModuleSlope && edm::isFinite(slopes[i])) - ? alpaka::math::abs(acc, alpaka::math::atan(acc, slopes[i])) - : kPi / 2.f); + absArctanSlope = + (edm::isFinite(slopes[i]) ? alpaka::math::abs(acc, alpaka::math::atan(acc, slopes[i])) : kPi / 2.f); if (xs[i] > 0 and ys[i] > 0) { angleM = kPi / 2.f - absArctanSlope; @@ -751,9 +750,8 @@ namespace ALPAKA_ACCELERATOR_NAMESPACE::lst { float chiSquared = 0.f; float absArctanSlope, angleM, xPrime, yPrime, sigma2; for (size_t i = 0; i < nPoints; i++) { - absArctanSlope = ((slopes[i] != kVerticalModuleSlope && edm::isFinite(slopes[i])) - ? alpaka::math::abs(acc, alpaka::math::atan(acc, slopes[i])) - : kPi / 2.f); + absArctanSlope = + (edm::isFinite(slopes[i]) ? alpaka::math::abs(acc, alpaka::math::atan(acc, slopes[i])) : kPi / 2.f); if (xs[i] > 0 and ys[i] > 0) { angleM = kPi / 2.f - absArctanSlope; } else if (xs[i] < 0 and ys[i] > 0) { diff --git a/RecoTracker/LSTCore/standalone/bin/lst.cc b/RecoTracker/LSTCore/standalone/bin/lst.cc index ee65f14d378ea..851f24df59cb7 100644 --- a/RecoTracker/LSTCore/standalone/bin/lst.cc +++ b/RecoTracker/LSTCore/standalone/bin/lst.cc @@ -369,7 +369,7 @@ void run_lst() { full_timer.Start(); // Determine which maps to use based on given pt cut for standalone. std::string ptCutString = (ana.ptCut >= 0.8) ? "0.8" : "0.6"; - auto hostESData = lst::loadAndFillESHost(ptCutString); + auto hostESData = lst::loadAndFillESDataHost(ptCutString); auto deviceESData = cms::alpakatools::CopyToDevice>::copyAsync(queues[0], *hostESData.get()); float timeForMapLoading = full_timer.RealTime() * 1000; diff --git a/RecoTracker/LSTCore/standalone/code/core/lst_math.h b/RecoTracker/LSTCore/standalone/code/core/lst_math.h index 3d31efcabb21b..cd4ef9d57f206 100644 --- a/RecoTracker/LSTCore/standalone/code/core/lst_math.h +++ b/RecoTracker/LSTCore/standalone/code/core/lst_math.h @@ -4,7 +4,12 @@ #include #include +#include "RecoTracker/LSTGeometry/interface/Helix.h" + namespace lst_math { + + using Helix = lstgeometry::Helix; + inline float Phi_mpi_pi(float x) { if (std::isnan(x)) { std::cout << "MathUtil::Phi_mpi_pi() function called with NaN" << std::endl; @@ -33,85 +38,6 @@ namespace lst_math { inline float ptEstimateFromRadius(float radius) { return 2.99792458e-3 * 3.8 * radius; }; - class Helix { - public: - std::vector center_; - float radius_; - float phi_; - float lam_; - float charge_; - - Helix(std::vector c, float r, float p, float l, float q) { - center_ = c; - radius_ = r; - phi_ = Phi_mpi_pi(p); - lam_ = l; - charge_ = q; - } - - Helix(float pt, float eta, float phi, float vx, float vy, float vz, float charge) { - // Radius based on pt - radius_ = pt / (2.99792458e-3 * 3.8); - phi_ = phi; - charge_ = charge; - - // reference point vector which for sim track is the vertex point - float ref_vec_x = vx; - float ref_vec_y = vy; - float ref_vec_z = vz; - - // The reference to center vector - float inward_radial_vec_x = charge_ * radius_ * sin(phi_); - float inward_radial_vec_y = charge_ * radius_ * -cos(phi_); - float inward_radial_vec_z = 0; - - // Center point - float center_vec_x = ref_vec_x + inward_radial_vec_x; - float center_vec_y = ref_vec_y + inward_radial_vec_y; - float center_vec_z = ref_vec_z + inward_radial_vec_z; - center_.push_back(center_vec_x); - center_.push_back(center_vec_y); - center_.push_back(center_vec_z); - - // Lambda - lam_ = copysign(M_PI / 2. - 2. * atan(exp(-std::abs(eta))), eta); - } - - const std::vector center() { return center_; } - const float radius() { return radius_; } - const float phi() { return phi_; } - const float lam() { return lam_; } - const float charge() { return charge_; } - - std::tuple get_helix_point(float t) { - float x = center()[0] - charge() * radius() * sin(phi() - (charge()) * t); - float y = center()[1] + charge() * radius() * cos(phi() - (charge()) * t); - float z = center()[2] + radius() * tan(lam()) * t; - float r = sqrt(x * x + y * y); - return std::make_tuple(x, y, z, r); - } - - float infer_t(const std::vector point) { - // Solve for t based on z position - float t = (point[2] - center()[2]) / (radius() * tan(lam())); - return t; - } - - float compare_radius(const std::vector point) { - float t = infer_t(point); - auto [x, y, z, r] = get_helix_point(t); - float point_r = sqrt(point[0] * point[0] + point[1] * point[1]); - return (point_r - r); - } - - float compare_xy(const std::vector point) { - float t = infer_t(point); - auto [x, y, z, r] = get_helix_point(t); - float xy_dist = sqrt(pow(point[0] - x, 2) + pow(point[1] - y, 2)); - return xy_dist; - } - }; - class Hit { public: float x_, y_, z_; diff --git a/RecoTracker/LSTCore/standalone/code/core/trkCore.cc b/RecoTracker/LSTCore/standalone/code/core/trkCore.cc index 7dc9026ea766b..343d78b60b85a 100644 --- a/RecoTracker/LSTCore/standalone/code/core/trkCore.cc +++ b/RecoTracker/LSTCore/standalone/code/core/trkCore.cc @@ -658,20 +658,20 @@ float drfracSimHitConsistentWithHelix(lst_math::Helix& helix, std::vector const& trk_simhit_y, std::vector const& trk_simhit_z) { // Sim hit vector - std::vector point = {trk_simhit_x[isimhitidx], trk_simhit_y[isimhitidx], trk_simhit_z[isimhitidx]}; + std::array point = {trk_simhit_x[isimhitidx], trk_simhit_y[isimhitidx], trk_simhit_z[isimhitidx]}; // Inferring parameter t of helix parametric function via z position - float t = helix.infer_t(point); + float t = helix.inferT(point); // If the best fit is more than pi parameter away then it's a returning hit and should be ignored if (not(t <= M_PI)) return 999; // Expected hit position with given z - auto [x, y, z, r] = helix.get_helix_point(t); + auto [x, y, z, r] = helix.point(t); // ( expected_r - simhit_r ) / expected_r - float drfrac = std::abs(helix.compare_radius(point)) / r; + float drfrac = std::abs(helix.compareRadius(point)) / r; return drfrac; } @@ -711,10 +711,10 @@ float distxySimHitConsistentWithHelix(lst_math::Helix& helix, std::vector const& trk_simhit_y, std::vector const& trk_simhit_z) { // Sim hit vector - std::vector point = {trk_simhit_x[isimhitidx], trk_simhit_y[isimhitidx], trk_simhit_z[isimhitidx]}; + std::array point = {trk_simhit_x[isimhitidx], trk_simhit_y[isimhitidx], trk_simhit_z[isimhitidx]}; // Inferring parameter t of helix parametric function via z position - float t = helix.infer_t(point); + float t = helix.inferT(point); // If the best fit is more than pi parameter away then it's a returning hit and should be ignored if (not(t <= M_PI)) @@ -724,7 +724,7 @@ float distxySimHitConsistentWithHelix(lst_math::Helix& helix, //auto [x, y, z, r] = helix.get_helix_point(t); // ( expected_r - simhit_r ) / expected_r - float distxy = helix.compare_xy(point); + float distxy = helix.compareXY(point); return distxy; } diff --git a/RecoTracker/LSTGeometry/BuildFile.xml b/RecoTracker/LSTGeometry/BuildFile.xml new file mode 100644 index 0000000000000..be1beb06995d9 --- /dev/null +++ b/RecoTracker/LSTGeometry/BuildFile.xml @@ -0,0 +1,8 @@ + + + + + + + + diff --git a/RecoTracker/LSTGeometry/interface/Common.h b/RecoTracker/LSTGeometry/interface/Common.h new file mode 100644 index 0000000000000..05e40da472dbd --- /dev/null +++ b/RecoTracker/LSTGeometry/interface/Common.h @@ -0,0 +1,45 @@ +#ifndef RecoTracker_LSTGeometry_interface_Common_h +#define RecoTracker_LSTGeometry_interface_Common_h + +#include +#include +#include +#include +#include +#include +#include +#include + +namespace lstgeometry { + + constexpr float kB = 3.8; + constexpr float kC = 0.00299792458; + + constexpr unsigned int kBarrelLayers = 6; + constexpr unsigned int kEndcapLayers = 5; + + // For pixel maps + constexpr unsigned int kNEta = 25; + constexpr unsigned int kNPhi = 72; + constexpr unsigned int kNZ = 25; + constexpr std::array kPtBounds = {{2.0, 10'000.0}}; + + // This is defined as a constant in case the legacy value (123456789) needs to be used + constexpr float kDefaultSlope = std::numeric_limits::infinity(); + + float degToRad(float degrees); + float phi_mpi_pi(float phi); + float roundAngle(float angle, float tol = 1e-3); + float roundCoordinate(float coord, float tol = 1e-3); + std::pair getEtaPhi(float x, float y, float z, float refphi = 0); +} // namespace lstgeometry + +namespace lst { + inline std::string floatToStr(float num, unsigned int precision = 1) { + std::ostringstream outSS; + outSS << std::setprecision(precision) << num; + return outSS.str(); + } +} // namespace lst + +#endif diff --git a/RecoTracker/LSTGeometry/interface/Geometry.h b/RecoTracker/LSTGeometry/interface/Geometry.h new file mode 100644 index 0000000000000..786909f5e6a7c --- /dev/null +++ b/RecoTracker/LSTGeometry/interface/Geometry.h @@ -0,0 +1,29 @@ +#ifndef RecoTracker_LSTGeometry_interface_Geometry_h +#define RecoTracker_LSTGeometry_interface_Geometry_h + +#include +#include + +#include "RecoTracker/LSTGeometry/interface/ModuleMap.h" +#include "RecoTracker/LSTGeometry/interface/PixelMap.h" +#include "RecoTracker/LSTGeometry/interface/Sensor.h" +#include "RecoTracker/LSTGeometry/interface/Slope.h" + +namespace lstgeometry { + + struct Geometry { + Sensors sensors; + Slopes barrel_slopes; + Slopes endcap_slopes; + PixelMap pixel_map; + ModuleMap module_map; + + Geometry(Sensors sensors, + std::array const &average_r_barrel, + std::array const &average_z_endcap, + float ptCut); + }; + +} // namespace lstgeometry + +#endif diff --git a/RecoTracker/LSTGeometry/interface/Helix.h b/RecoTracker/LSTGeometry/interface/Helix.h new file mode 100644 index 0000000000000..d3ff36a022297 --- /dev/null +++ b/RecoTracker/LSTGeometry/interface/Helix.h @@ -0,0 +1,143 @@ +#ifndef RecoTracker_LSTGeometry_interface_Helix_h +#define RecoTracker_LSTGeometry_interface_Helix_h + +// This header contains implementations so that the standalone part of LSTCore can use it without having to link some extra library. + +#include +#include +#include + +#include "RecoTracker/LSTGeometry/interface/Common.h" + +namespace lstgeometry { + + struct Helix { + std::array center; + float radius; + float phi; + float lambda; + int charge; + + Helix(std::array center, float radius, float phi, float lambda, int charge) + : center(center), radius(radius), phi(phi_mpi_pi(phi)), lambda(lambda), charge(charge) {} + + // Clarification : phi was derived assuming a negatively charged particle would start + // at the first quadrant. However the way signs are set up in the get_track_point function + // implies the particle actually starts out in the fourth quadrant, and phi is measured from + // the y axis as opposed to x axis in the expression provided in this function. Hence I tucked + // in an extra pi/2 to account for these effects + Helix(float pt, float vx, float vy, float vz, float mx, float my, float mz, int charge) : charge(charge) { + radius = pt / (kC * kB); + + float t = 2. * std::asin(std::sqrt((vx - mx) * (vx - mx) + (vy - my) * (vy - my)) / (2. * radius)); + phi = std::numbers::pi_v / 2. + std::atan((vy - my) / (vx - mx)) + + ((vy - my) / (vx - mx) < 0) * (std::numbers::pi_v)+charge * t / 2. + + (my - vy < 0) * (std::numbers::pi_v / 2.) - (my - vy > 0) * (std::numbers::pi_v / 2.); + + center[0] = vx + charge * radius * std::sin(phi); + center[1] = vy - charge * radius * std::cos(phi); + center[2] = vz; + lambda = std::atan((mz - vz) / (radius * t)); + } + + Helix(float pt, float eta, float phi, float vx, float vy, float vz, float charge) : phi(phi), charge(charge) { + // Radius based on pt + radius = pt / (kC * kB); + + // reference point vector which for sim track is the vertex point + float ref_vec_x = vx; + float ref_vec_y = vy; + float ref_vec_z = vz; + + // The reference to center vector + float inward_radial_vec_x = charge * radius * std::sin(phi); + float inward_radial_vec_y = charge * radius * -std::cos(phi); + float inward_radial_vec_z = 0; + + // Center point + center[0] = ref_vec_x + inward_radial_vec_x; + center[1] = ref_vec_y + inward_radial_vec_y; + center[2] = ref_vec_z + inward_radial_vec_z; + + // Lambda + lambda = std::atan(std::sinh(eta)); + } + + std::array point(float t) const { + float x = center[0] - charge * radius * std::sin(phi - charge * t); + float y = center[1] + charge * radius * std::cos(phi - charge * t); + float z = center[2] + radius * std::tan(lambda) * t; + float r = std::sqrt(x * x + y * y); + return {x, y, z, r}; + } + + float inferT(std::array const& point) const { + return (point[2] - center[2]) / (radius * std::tan(lambda)); + } + + std::tuple pointFromRadius(float target_r) const { + float cx = center[0], cy = center[1]; + float center_r = std::sqrt(cx * cx + cy * cy); + + float sin_val = (cx * cx + cy * cy + radius * radius - target_r * target_r) / (2.f * charge * radius * center_r); + sin_val = std::clamp(sin_val, -1.f, 1.f); + + float offset = std::atan2(-cy, cx); + float alpha = std::asin(sin_val); + + // Two solutions: theta = alpha - offset or pi - alpha - offset + // theta = phi - charge*t => t = charge*(phi - theta) + auto to_t = [this](float theta) { + float t = charge * (phi - theta); + const float two_pi = 2.f * std::numbers::pi_v; + t = std::fmod(t, two_pi); + if (t < 0.f) + t += two_pi; + return t; + }; + + float t1 = to_t(alpha - offset); + float t2 = to_t(std::numbers::pi_v - alpha - offset); + + // Prefer the smallest t in [0, pi], matching the original search range + bool t1_ok = t1 <= std::numbers::pi_v; + bool t2_ok = t2 <= std::numbers::pi_v; + float t = (t1_ok && t2_ok) ? std::min(t1, t2) : t1_ok ? t1 : t2_ok ? t2 : std::min(t1, t2); + + float x = center[0] - charge * radius * std::sin(phi - charge * t); + float y = center[1] + charge * radius * std::cos(phi - charge * t); + float z = center[2] + radius * std::tan(lambda) * t; + float r = std::sqrt(x * x + y * y); + + return std::make_tuple(x, y, z, r); + } + + std::tuple pointFromZ(float target_z) const { + float t = (target_z - center[2]) / (radius * std::tan(lambda)); + + float x = center[0] - charge * radius * std::sin(phi - charge * t); + float y = center[1] + charge * radius * std::cos(phi - charge * t); + float z = center[2] + radius * std::tan(lambda) * t; + float r = std::sqrt(x * x + y * y); + + return std::make_tuple(x, y, z, r); + } + + float compareRadius(std::array const& pt) const { + float t = inferT(pt); + auto [x, y, z, r] = point(t); + float point_r = std::sqrt(pt[0] * pt[0] + pt[1] * pt[1]); + return (point_r - r); + } + + float compareXY(std::array const& pt) const { + float t = inferT(pt); + auto [x, y, z, r] = point(t); + float xy_dist = std::sqrt(std::pow(pt[0] - x, 2) + std::pow(pt[1] - y, 2)); + return xy_dist; + } + }; + +} // namespace lstgeometry + +#endif diff --git a/RecoTracker/LSTGeometry/interface/IO.h b/RecoTracker/LSTGeometry/interface/IO.h new file mode 100644 index 0000000000000..cf110eca35171 --- /dev/null +++ b/RecoTracker/LSTGeometry/interface/IO.h @@ -0,0 +1,21 @@ +#ifndef RecoTracker_LSTGeometry_interface_IO_h +#define RecoTracker_LSTGeometry_interface_IO_h + +#include +#include + +#include "RecoTracker/LSTGeometry/interface/ModuleMap.h" +#include "RecoTracker/LSTGeometry/interface/PixelMap.h" +#include "RecoTracker/LSTGeometry/interface/Sensor.h" +#include "RecoTracker/LSTGeometry/interface/Slope.h" + +namespace lstgeometry { + + void writeSensorCentroids(Sensors const& sensors, std::string const& base_filename, bool binary = true); + void writeSlopes(Slopes const& slopes, Sensors const& sensors, std::string const& base_filename, bool binary = true); + void writeModuleConnections(ModuleMap const& connections, std::string const& base_filename, bool binary = true); + void writePixelMaps(PixelMap const& maps, std::string const& base_filename, bool binary = true); + +} // namespace lstgeometry + +#endif diff --git a/RecoTracker/LSTGeometry/interface/ModuleMap.h b/RecoTracker/LSTGeometry/interface/ModuleMap.h new file mode 100644 index 0000000000000..6e75888403dce --- /dev/null +++ b/RecoTracker/LSTGeometry/interface/ModuleMap.h @@ -0,0 +1,22 @@ +#ifndef RecoTracker_LSTGeometry_interface_ModuleMapMethods_h +#define RecoTracker_LSTGeometry_interface_ModuleMapMethods_h + +#include +#include + +#include "RecoTracker/LSTGeometry/interface/Sensor.h" +#include "RecoTracker/LSTGeometry/interface/SensorBinning.h" + +namespace lstgeometry { + + using ModuleMap = std::unordered_map>; + + ModuleMap buildModuleMap(Sensors const& sensors, + BinnedDetIds const& binned_detids, + std::array const& average_r_barrel, + std::array const& average_z_endcap, + float pt_cut); + +} // namespace lstgeometry + +#endif diff --git a/RecoTracker/LSTGeometry/interface/PixelMap.h b/RecoTracker/LSTGeometry/interface/PixelMap.h new file mode 100644 index 0000000000000..07a96265e5168 --- /dev/null +++ b/RecoTracker/LSTGeometry/interface/PixelMap.h @@ -0,0 +1,22 @@ +#ifndef RecoTracker_LSTGeometry_interface_PixelMap_h +#define RecoTracker_LSTGeometry_interface_PixelMap_h + +#include +#include +#include + +#include "RecoTracker/LSTGeometry/interface/Sensor.h" +#include "RecoTracker/LSTGeometry/interface/SensorBinning.h" + +namespace lstgeometry { + + using LayerSubdetChargeKey = std::tuple; + using LayerSubdetChargeMap = + std::unordered_map>, boost::hash>; + using PixelMap = LayerSubdetChargeMap; + + PixelMap buildPixelMap(Sensors const& sensors, float pt_cut); + +} // namespace lstgeometry + +#endif diff --git a/RecoTracker/LSTGeometry/interface/Sensor.h b/RecoTracker/LSTGeometry/interface/Sensor.h new file mode 100644 index 0000000000000..45100825126f1 --- /dev/null +++ b/RecoTracker/LSTGeometry/interface/Sensor.h @@ -0,0 +1,102 @@ +#ifndef RecoTracker_LSTGeometry_interface_Sensor_h +#define RecoTracker_LSTGeometry_interface_Sensor_h + +// Some parts of this file are guarded by the LST_STANDALONE flag to avoid depending on Eigen for standalone LST + +#include +#include + +#ifndef LST_STANDALONE +#include +#endif + +#include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h" +#include "Geometry/CommonTopologies/interface/GeomDetEnumerators.h" +#include "DataFormats/SiStripDetId/interface/SiStripEnums.h" +#include "DataFormats/GeometrySurface/interface/Surface.h" + +#include "RecoTracker/LSTGeometry/interface/Common.h" + +namespace lstgeometry { + +#ifndef LST_STANDALONE + using MatrixF3x3 = Eigen::Matrix; + using MatrixF4x2 = Eigen::Matrix; + using MatrixF4x3 = Eigen::Matrix; + using MatrixF8x3 = Eigen::Matrix; + using MatrixFNx2 = Eigen::Matrix; + using MatrixFNx3 = Eigen::Matrix; + using RowVectorF2 = Eigen::Matrix; + using ColVectorF3 = Eigen::Matrix; + using RowVectorF3 = Eigen::Matrix; +#endif + + using ModuleType = TrackerGeometry::ModuleType; + using SubDetector = GeomDetEnumerators::SubDetector; + using Location = GeomDetEnumerators::Location; + + enum SubDet { InnerPixel = 0, Barrel = 5, Endcap = 4 }; + enum Side { NegZ = 1, PosZ = 2, Center = 3 }; + + struct SensorExtraData { + // Module-level properties + SubDetector subdet; + Location location; + Side side; + unsigned int moduleId; + unsigned int layer; + unsigned int ring; + bool inverted; + // Sensor-level properties + float centerRho; + float centerPhi; + bool lower; + bool strip; +#ifndef LST_STANDALONE // Extra data is not used in standalone, so it doesn't matter that there is a struct mismatch + MatrixF4x3 corners; +#endif + // Redundant, but convenient to avoid repeated computations + float centerEta; + float centerTheta; + float minR; + float maxR; + float minZ; + float maxZ; + float minPhi; + float maxPhi; + }; + + struct Sensor { + // Info that is always needed + ModuleType moduleType; + float centerPhi; + float centerX; + float centerY; + float centerZ; + // Info that can be dropped after map construction + std::unique_ptr extra; + + Sensor() = default; + Sensor(const Sensor&) = delete; + Sensor& operator=(const Sensor&) = delete; + Sensor(Sensor&&) = default; + Sensor& operator=(Sensor&&) = default; + Sensor(unsigned int detId, + ModuleType moduleType, + SubDetector subdet, + Location location, + Side side, + unsigned int moduleId, + unsigned int layer, + unsigned int ring, + float centerRho, + float centerZ, + float centerPhi, + Surface const& surface); + }; + + using Sensors = std::unordered_map; + +} // namespace lstgeometry + +#endif diff --git a/RecoTracker/LSTGeometry/interface/SensorBinning.h b/RecoTracker/LSTGeometry/interface/SensorBinning.h new file mode 100644 index 0000000000000..b8b02d6688b45 --- /dev/null +++ b/RecoTracker/LSTGeometry/interface/SensorBinning.h @@ -0,0 +1,38 @@ +#ifndef RecoTracker_LSTGeometry_interface_SensorBinning_h +#define RecoTracker_LSTGeometry_interface_SensorBinning_h + +#include +#include +#include +#include +#include +#include + +#include "RecoTracker/LSTGeometry/interface/Common.h" +#include "RecoTracker/LSTGeometry/interface/Sensor.h" + +namespace lstgeometry { + + using LocationLayerThetaBinPhiBinKey = std::tuple; + + using BinnedDetIds = std::unordered_map, + boost::hash>; + + // We split modules into overlapping theta-phi bins so that it's easier to construct module maps + static constexpr unsigned int kNThetaBins = 4; + static constexpr float kThetaBinRad = std::numbers::pi_v / kNThetaBins; + static constexpr unsigned int kNPhiBins = 6; + static constexpr float kPhiBinWidth = 2 * std::numbers::pi_v / kNPhiBins; + + bool isInThetaPhiBin(float theta, float phi, unsigned int theta_bin, unsigned int phi_bin); + std::pair getThetaPhiBins(float theta, float phi); + std::pair getCompatibleEtaRange(Sensor const& sensor, float zmin_bound, float zmax_bound); + std::pair, std::pair> getCompatiblePhiRange(Sensor const& sensor, + float ptmin, + float ptmax); + BinnedDetIds binDetIds(Sensors const& sensors); + +} // namespace lstgeometry + +#endif diff --git a/RecoTracker/LSTGeometry/interface/Slope.h b/RecoTracker/LSTGeometry/interface/Slope.h new file mode 100644 index 0000000000000..d873f853c71d3 --- /dev/null +++ b/RecoTracker/LSTGeometry/interface/Slope.h @@ -0,0 +1,24 @@ +#ifndef RecoTracker_LSTGeometry_interface_Slope_h +#define RecoTracker_LSTGeometry_interface_Slope_h + +#include + +#include "RecoTracker/LSTGeometry/interface/Sensor.h" + +namespace lstgeometry { + + struct Slope { + float drdz; + float dxdy; + + Slope() = default; + Slope(float dx, float dy, float dz); + }; + + using Slopes = std::unordered_map; + + std::tuple computeSlopes(Sensors const& sensors); + +} // namespace lstgeometry + +#endif diff --git a/RecoTracker/LSTGeometry/plugins/BuildFile.xml b/RecoTracker/LSTGeometry/plugins/BuildFile.xml new file mode 100644 index 0000000000000..b09736f9f7259 --- /dev/null +++ b/RecoTracker/LSTGeometry/plugins/BuildFile.xml @@ -0,0 +1,14 @@ + + + + + + + + + + + + + + diff --git a/RecoTracker/LSTGeometry/plugins/LSTGeometryESProducer.cc b/RecoTracker/LSTGeometry/plugins/LSTGeometryESProducer.cc new file mode 100644 index 0000000000000..de12474729a42 --- /dev/null +++ b/RecoTracker/LSTGeometry/plugins/LSTGeometryESProducer.cc @@ -0,0 +1,128 @@ +#include +#include +#include +#include + +#include "FWCore/Framework/interface/ModuleFactory.h" +#include "FWCore/Framework/interface/ESProducer.h" + +#include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h" +#include "Geometry/Records/interface/TrackerTopologyRcd.h" +#include "DataFormats/TrackerCommon/interface/TrackerTopology.h" +#include "RecoTracker/Record/interface/TrackerRecoGeometryRecord.h" +#include "DataFormats/GeometrySurface/interface/RectangularPlaneBounds.h" +#include "Geometry/CommonTopologies/interface/GeomDetEnumerators.h" + +#include "RecoTracker/LSTGeometry/interface/Geometry.h" +#include "RecoTracker/LSTGeometry/interface/Sensor.h" + +class LSTGeometryESProducer : public edm::ESProducer { +public: + LSTGeometryESProducer(const edm::ParameterSet &iConfig); + + static void fillDescriptions(edm::ConfigurationDescriptions &descriptions); + + std::unique_ptr produce(const TrackerRecoGeometryRecord &iRecord); + +private: + double ptCut_; + + edm::ESGetToken geomToken_; + edm::ESGetToken ttopoToken_; + + const TrackerGeometry *trackerGeom_ = nullptr; + const TrackerTopology *trackerTopo_ = nullptr; +}; + +LSTGeometryESProducer::LSTGeometryESProducer(const edm::ParameterSet &iConfig) + : ptCut_(iConfig.getParameter("ptCut")) { + std::string ptCutStr = lst::floatToStr(ptCut_, 1); + + auto cc = setWhatProduced(this, ptCutStr); + geomToken_ = cc.consumes(); + ttopoToken_ = cc.consumes(); +} + +void LSTGeometryESProducer::fillDescriptions(edm::ConfigurationDescriptions &descriptions) { + edm::ParameterSetDescription desc; + desc.add("ptCut", 0.8); + descriptions.addWithDefaultLabel(desc); +} + +std::unique_ptr LSTGeometryESProducer::produce(const TrackerRecoGeometryRecord &iRecord) { + trackerGeom_ = &iRecord.get(geomToken_); + trackerTopo_ = &iRecord.get(ttopoToken_); + + lstgeometry::Sensors sensors; + + std::array avg_r_barrel{}; + std::array avg_z_endcap{}; + std::array avg_r_barrel_counter{}; + std::array avg_z_endcap_counter{}; + + for (auto &det : trackerGeom_->dets()) { + const DetId detId = det->geographicalId(); + const auto moduleType = trackerGeom_->getDetectorType(detId); + + // TODO: Is there a more straightforward way to only loop through these? + if (moduleType != TrackerGeometry::ModuleType::Ph2PSP && moduleType != TrackerGeometry::ModuleType::Ph2PSS && + moduleType != TrackerGeometry::ModuleType::Ph2SS) { + continue; + } + + const unsigned int detid = detId(); + + const auto &surface = det->surface(); + const auto &position = surface.position(); + + const float rho_cm = position.perp(); + const float z_cm = lstgeometry::roundCoordinate(position.z()); + const float phi_rad = lstgeometry::roundAngle(position.phi()); + + const auto subdet = trackerGeom_->geomDetSubDetector(detId.subdetId()); + const auto location = + GeomDetEnumerators::isBarrel(subdet) ? lstgeometry::Location::barrel : lstgeometry::Location::endcap; + const auto side = static_cast( + location == lstgeometry::Location::barrel ? static_cast(trackerTopo_->barrelTiltTypeP2(detId)) + : trackerTopo_->side(detId)); + const unsigned int moduleId = trackerTopo_->module(detId); + const unsigned int layer = trackerTopo_->layer(detId); + const unsigned int ring = trackerTopo_->endcapRingP2(detId); + + if (layer == 0 || (location == lstgeometry::Location::barrel && layer > lstgeometry::kBarrelLayers) || + (location == lstgeometry::Location::endcap && layer > lstgeometry::kEndcapLayers)) { + continue; + } + + if (det->isLeaf()) { + // Leafs are the sensors + sensors[detid] = lstgeometry::Sensor( + detid, moduleType, subdet, location, side, moduleId, layer, ring, rho_cm, z_cm, phi_rad, surface); + + continue; + } + + if (location == lstgeometry::Location::barrel) { + avg_r_barrel[layer - 1] += rho_cm; + avg_r_barrel_counter[layer - 1] += 1; + } else { + avg_z_endcap[layer - 1] += std::fabs(z_cm); + avg_z_endcap_counter[layer - 1] += 1; + } + } + + for (size_t i = 0; i < avg_r_barrel.size(); ++i) { + if (avg_r_barrel_counter[i] > 0) + avg_r_barrel[i] /= avg_r_barrel_counter[i]; + } + for (size_t i = 0; i < avg_z_endcap.size(); ++i) { + if (avg_z_endcap_counter[i] > 0) + avg_z_endcap[i] /= avg_z_endcap_counter[i]; + } + + auto lstGeometry = std::make_unique(std::move(sensors), avg_r_barrel, avg_z_endcap, ptCut_); + + return lstGeometry; +} + +DEFINE_FWK_EVENTSETUP_MODULE(LSTGeometryESProducer); diff --git a/RecoTracker/LSTGeometry/src/Common.cc b/RecoTracker/LSTGeometry/src/Common.cc new file mode 100644 index 0000000000000..5d2bdbf809eff --- /dev/null +++ b/RecoTracker/LSTGeometry/src/Common.cc @@ -0,0 +1,42 @@ +#include "RecoTracker/LSTGeometry/interface/Common.h" + +namespace lstgeometry { + + float degToRad(float degrees) { return degrees * (std::numbers::pi_v / 180); } + + float phi_mpi_pi(float phi) { + while (phi >= std::numbers::pi_v) + phi -= 2 * std::numbers::pi_v; + while (phi < -std::numbers::pi_v) + phi += 2 * std::numbers::pi_v; + return phi; + } + + float roundAngle(float angle, float tol) { + const float pi = std::numbers::pi_v; + if (std::fabs(angle) < tol) { + return 0.0; + } else if (std::fabs(angle - pi / 2) < tol) { + return pi / 2; + } else if (std::fabs(angle + pi / 2) < tol) { + return -pi / 2; + } else if (std::fabs(angle - pi) < tol || std::fabs(angle + pi) < tol) { + return -pi; + } + return angle; + } + + float roundCoordinate(float coord, float tol) { + if (std::fabs(coord) < tol) { + return 0.0; + } + return coord; + } + + std::pair getEtaPhi(float x, float y, float z, float refphi) { + float phi = phi_mpi_pi(std::atan2(y, x) - refphi); + float eta = std::asinh(z / std::sqrt(x * x + y * y)); + return std::make_pair(eta, phi); + } + +} // namespace lstgeometry diff --git a/RecoTracker/LSTGeometry/src/Geometry.cc b/RecoTracker/LSTGeometry/src/Geometry.cc new file mode 100644 index 0000000000000..caea8d31a1e70 --- /dev/null +++ b/RecoTracker/LSTGeometry/src/Geometry.cc @@ -0,0 +1,30 @@ +#include "RecoTracker/LSTGeometry/interface/Geometry.h" +#include "RecoTracker/LSTGeometry/interface/SensorBinning.h" +#include "RecoTracker/LSTGeometry/interface/Slope.h" + +using namespace lstgeometry; + +Geometry::Geometry(Sensors sensors, + std::array const &average_r_barrel, + std::array const &average_z_endcap, + float pt_cut) { + auto slopes = computeSlopes(sensors); + barrel_slopes = std::move(std::get<0>(slopes)); + endcap_slopes = std::move(std::get<1>(slopes)); + + auto binned_detids = binDetIds(sensors); + + pixel_map = buildPixelMap(sensors, pt_cut); + + module_map = buildModuleMap(sensors, binned_detids, average_r_barrel, average_z_endcap, pt_cut); + + // Drop all the extra data that is no longer needed + for (auto &[detId, sensor] : sensors) { + sensor.extra.reset(); + } + + this->sensors = std::move(sensors); +} + +#include "FWCore/Utilities/interface/typelookup.h" +TYPELOOKUP_DATA_REG(lstgeometry::Geometry); diff --git a/RecoTracker/LSTGeometry/src/IO.cc b/RecoTracker/LSTGeometry/src/IO.cc new file mode 100644 index 0000000000000..8d4c9bb9b3878 --- /dev/null +++ b/RecoTracker/LSTGeometry/src/IO.cc @@ -0,0 +1,137 @@ +#include +#include +#include +#include +#include + +#include "RecoTracker/LSTGeometry/interface/IO.h" + +namespace lstgeometry { + + void writeSensorCentroids(Sensors const& sensors, std::string const& base_filename, bool binary) { + std::filesystem::path filepath(base_filename); + std::filesystem::create_directories(filepath.parent_path()); + + std::string filename = base_filename + (binary ? ".bin" : ".txt"); + std::ofstream file(filename, binary ? std::ios::binary : std::ios::out); + + if (binary) { + for (auto const& [detid, sensor] : sensors) { + float x = sensor.centerX; + float y = sensor.centerY; + float z = sensor.centerZ; + unsigned int moduleType = static_cast(sensor.moduleType); + file.write(reinterpret_cast(&detid), sizeof(detid)); + file.write(reinterpret_cast(&x), sizeof(x)); + file.write(reinterpret_cast(&y), sizeof(y)); + file.write(reinterpret_cast(&z), sizeof(z)); + file.write(reinterpret_cast(&moduleType), sizeof(moduleType)); + } + } else { + for (auto const& [detid, sensor] : sensors) { + file << detid << "," << sensor.centerX << "," << sensor.centerY << "," << sensor.centerZ << "," + << static_cast(sensor.moduleType) << std::endl; + } + } + } + + void writeSlopes(Slopes const& slopes, Sensors const& sensors, std::string const& base_filename, bool binary) { + std::filesystem::path filepath(base_filename); + std::filesystem::create_directories(filepath.parent_path()); + + std::string filename = base_filename + (binary ? ".bin" : ".txt"); + std::ofstream file(filename, binary ? std::ios::binary : std::ios::out); + + if (binary) { + for (auto const& [detid, slope] : slopes) { + float drdz_slope = slope.drdz; + float dxdy_slope = slope.dxdy; + float phi = sensors.at(detid).centerPhi; + file.write(reinterpret_cast(&detid), sizeof(detid)); + if (drdz_slope != kDefaultSlope) { + file.write(reinterpret_cast(&drdz_slope), sizeof(drdz_slope)); + file.write(reinterpret_cast(&dxdy_slope), sizeof(dxdy_slope)); + } else { + file.write(reinterpret_cast(&dxdy_slope), sizeof(dxdy_slope)); + file.write(reinterpret_cast(&phi), sizeof(phi)); + } + } + } else { + for (auto const& [detid, slope] : slopes) { + float drdz_slope = slope.drdz; + float dxdy_slope = slope.dxdy; + float phi = sensors.at(detid).centerPhi; + file << detid << ","; + if (drdz_slope != kDefaultSlope) { + file << drdz_slope << "," << dxdy_slope << std::endl; + } else { + file << dxdy_slope << "," << phi << std::endl; + } + } + } + } + + void writeModuleConnections(ModuleMap const& connections, std::string const& base_filename, bool binary) { + std::filesystem::path filepath(base_filename); + std::filesystem::create_directories(filepath.parent_path()); + + std::string filename = base_filename + (binary ? ".bin" : ".txt"); + std::ofstream file(filename, binary ? std::ios::binary : std::ios::out); + + if (binary) { + for (auto const& [detid, list] : connections) { + file.write(reinterpret_cast(&detid), sizeof(detid)); + unsigned int length = list.size(); + file.write(reinterpret_cast(&length), sizeof(length)); + for (unsigned int i : list) { + file.write(reinterpret_cast(&i), sizeof(i)); + } + } + } else { + for (auto const& [detid, list] : connections) { + file << detid << "," << list.size(); + for (unsigned int i : list) { + file << "," << i; + } + file << std::endl; + } + } + } + + void writePixelMaps(PixelMap const& maps, std::string const& base_filename, bool binary) { + std::filesystem::path filepath(base_filename); + std::filesystem::create_directories(filepath.parent_path()); + + for (auto const& [layersubdetcharge, vec_of_vecs] : maps) { + auto const& [layer, subdet, charge] = layersubdetcharge; + + std::string charge_str = charge > 0 ? "_pos" : (charge < 0 ? "_neg" : ""); + std::string filename = + std::format("{}{}_layer{}_subdet{}.{}", base_filename, charge_str, layer, subdet, binary ? "bin" : "txt"); + + std::ofstream file(filename, binary ? std::ios::binary : std::ios::out); + + for (unsigned int isuperbin = 0; isuperbin < vec_of_vecs.size(); isuperbin++) { + auto const& vec = vec_of_vecs.at(isuperbin); + if (vec.empty()) + continue; + + unsigned int length = vec.size(); + if (binary) { + file.write(reinterpret_cast(&isuperbin), sizeof(isuperbin)); + file.write(reinterpret_cast(&length), sizeof(length)); + for (unsigned int i : vec) { + file.write(reinterpret_cast(&i), sizeof(i)); + } + } else { + file << isuperbin << "," << length; + for (unsigned int i : vec) { + file << "," << i; + } + file << std::endl; + } + } + } + } + +} // namespace lstgeometry diff --git a/RecoTracker/LSTGeometry/src/ModuleMap.cc b/RecoTracker/LSTGeometry/src/ModuleMap.cc new file mode 100644 index 0000000000000..01bc9c7547737 --- /dev/null +++ b/RecoTracker/LSTGeometry/src/ModuleMap.cc @@ -0,0 +1,406 @@ +#include +#include +#include +#include +#include +#include + +#include "RecoTracker/LSTGeometry/interface/Helix.h" +#include "RecoTracker/LSTGeometry/interface/ModuleMap.h" + +namespace lstgeometry { + + using Point = boost::geometry::model::d2::point_xy; + using Polygon = boost::geometry::model::polygon; + + Polygon getEtaPhiPolygon(MatrixF4x3 const& corners, float refphi, float zshift = 0) { + MatrixF4x2 mod_boundaries_etaphi; + for (int i = 0; i < 4; ++i) { + auto ref_etaphi = getEtaPhi(corners(i, 1), corners(i, 2), corners(i, 0) + zshift, refphi); + mod_boundaries_etaphi(i, 0) = ref_etaphi.first; + mod_boundaries_etaphi(i, 1) = ref_etaphi.second; + } + + Polygon poly; + // <= because we need to close the polygon with the first point + for (int i = 0; i <= 4; ++i) { + boost::geometry::append(poly, Point(mod_boundaries_etaphi(i % 4, 0), mod_boundaries_etaphi(i % 4, 1))); + } + boost::geometry::correct(poly); + return poly; + } + + bool moduleOverlapsInEtaPhi(MatrixF4x3 const& ref_mod_boundaries, + MatrixF4x3 const& tar_mod_boundaries, + float refphi = 0, + float zshift = 0, + float ref_center_phi = std::numeric_limits::max(), + float tar_center_phi = std::numeric_limits::max()) { + if (ref_center_phi < -std::numbers::pi_v || ref_center_phi > std::numbers::pi_v) { + RowVectorF3 ref_center = ref_mod_boundaries.colwise().sum(); + ref_center /= 4.; + ref_center_phi = std::atan2(ref_center(2), ref_center(1)); + } + if (tar_center_phi < -std::numbers::pi_v || tar_center_phi > std::numbers::pi_v) { + RowVectorF3 tar_center = tar_mod_boundaries.colwise().sum(); + tar_center /= 4.; + tar_center_phi = std::atan2(tar_center(2), tar_center(1)); + } + + if (std::fabs(phi_mpi_pi(ref_center_phi - tar_center_phi)) > std::numbers::pi_v / 2.) + return false; + + MatrixF4x2 ref_mod_boundaries_etaphi; + MatrixF4x2 tar_mod_boundaries_etaphi; + + for (int i = 0; i < 4; ++i) { + auto ref_etaphi = + getEtaPhi(ref_mod_boundaries(i, 1), ref_mod_boundaries(i, 2), ref_mod_boundaries(i, 0) + zshift, refphi); + auto tar_etaphi = + getEtaPhi(tar_mod_boundaries(i, 1), tar_mod_boundaries(i, 2), tar_mod_boundaries(i, 0) + zshift, refphi); + ref_mod_boundaries_etaphi(i, 0) = ref_etaphi.first; + ref_mod_boundaries_etaphi(i, 1) = ref_etaphi.second; + tar_mod_boundaries_etaphi(i, 0) = tar_etaphi.first; + tar_mod_boundaries_etaphi(i, 1) = tar_etaphi.second; + } + + // Quick cut + RowVectorF2 diff = ref_mod_boundaries_etaphi.row(0) - tar_mod_boundaries_etaphi.row(0); + if (std::fabs(diff(0)) > 0.5) + return false; + if (std::fabs(phi_mpi_pi(diff(1))) > 1.) + return false; + + Polygon ref_poly, tar_poly; + + // <= 4 because we need to close the polygon with the first point + for (int i = 0; i <= 4; ++i) { + boost::geometry::append(ref_poly, + Point(ref_mod_boundaries_etaphi(i % 4, 0), ref_mod_boundaries_etaphi(i % 4, 1))); + boost::geometry::append(tar_poly, + Point(tar_mod_boundaries_etaphi(i % 4, 0), tar_mod_boundaries_etaphi(i % 4, 1))); + } + boost::geometry::correct(ref_poly); + boost::geometry::correct(tar_poly); + + return boost::geometry::intersects(ref_poly, tar_poly); + } + + std::vector getStraightLineConnections(unsigned int ref_detid, + Sensors const& sensors, + BinnedDetIds const& binned_detids) { + auto const& ref_sensor = sensors.at(ref_detid); + + float refphi = std::atan2(ref_sensor.centerY, ref_sensor.centerX); + unsigned short ref_layer = ref_sensor.extra->layer; + auto ref_location = ref_sensor.extra->location; + + auto thetaphibins = getThetaPhiBins(ref_sensor.extra->centerTheta, ref_sensor.centerPhi); + + auto const& tar_detids_to_be_considered = + binned_detids.at(std::make_tuple(ref_location, ref_layer + 1, thetaphibins.first, thetaphibins.second)); + + std::vector list_of_detids_etaphi_layer_tar; + for (unsigned int tar_detid : tar_detids_to_be_considered) { + auto const& tar_sensor = sensors.at(tar_detid); + if (moduleOverlapsInEtaPhi(ref_sensor.extra->corners, + tar_sensor.extra->corners, + refphi, + 0, + ref_sensor.centerPhi, + tar_sensor.centerPhi) || + moduleOverlapsInEtaPhi(ref_sensor.extra->corners, + tar_sensor.extra->corners, + refphi, + 10, + ref_sensor.centerPhi, + tar_sensor.centerPhi) || + moduleOverlapsInEtaPhi(ref_sensor.extra->corners, + tar_sensor.extra->corners, + refphi, + -10, + ref_sensor.centerPhi, + tar_sensor.centerPhi)) + list_of_detids_etaphi_layer_tar.push_back(tar_detid); + } + + // Consider barrel to endcap connections if the intersection area is > 0 + // We construct the reference polygon as a vector of polygons because the boost::geometry::difference + // function can return multiple polygons if the difference results in disjoint pieces + if (ref_location == Location::barrel) { + std::vector barrel_endcap_connected_tar_detids; + + for (float zshift : {0, 10, -10}) { + std::vector ref_polygon; + ref_polygon.push_back(getEtaPhiPolygon(ref_sensor.extra->corners, refphi, zshift)); + + // Check whether there is still significant non-zero area + for (unsigned int tar_detid : list_of_detids_etaphi_layer_tar) { + if (ref_polygon.empty()) + break; + Polygon tar_polygon = getEtaPhiPolygon(sensors.at(tar_detid).extra->corners, refphi, zshift); + + std::vector difference; + for (auto const& ref_polygon_piece : ref_polygon) { + std::vector tmp_difference; + boost::geometry::difference(ref_polygon_piece, tar_polygon, tmp_difference); + difference.insert(difference.end(), tmp_difference.begin(), tmp_difference.end()); + } + + ref_polygon = std::move(difference); + } + + float area = 0.; + for (auto const& ref_polygon_piece : ref_polygon) + area += boost::geometry::area(ref_polygon_piece); + + if (area <= 5e-3) + continue; + + auto const& new_tar_detids_to_be_considered = + binned_detids.at(std::make_tuple(Location::endcap, 1, thetaphibins.first, thetaphibins.second)); + + for (unsigned int tar_detid : new_tar_detids_to_be_considered) { + auto const& sensor_target = sensors.at(tar_detid); + float tarphi = std::atan2(sensor_target.centerY, sensor_target.centerX); + + if (std::fabs(phi_mpi_pi(tarphi - refphi)) > std::numbers::pi_v / 2.) + continue; + + Polygon tar_polygon = getEtaPhiPolygon(sensor_target.extra->corners, refphi, zshift); + + bool intersects = false; + for (auto const& ref_polygon_piece : ref_polygon) { + if (boost::geometry::intersects(ref_polygon_piece, tar_polygon)) { + intersects = true; + break; + } + } + + if (intersects) + barrel_endcap_connected_tar_detids.push_back(tar_detid); + } + } + list_of_detids_etaphi_layer_tar.insert(list_of_detids_etaphi_layer_tar.end(), + barrel_endcap_connected_tar_detids.begin(), + barrel_endcap_connected_tar_detids.end()); + } + + return list_of_detids_etaphi_layer_tar; + } + + MatrixF4x3 boundsAfterCurved(unsigned int ref_detid, + Sensors const& sensors, + std::array const& average_r_barrel, + std::array const& average_z_endcap, + float ptCut, + bool doR = true) { + auto const& ref_sensor = sensors.at(ref_detid); + auto const& bounds = ref_sensor.extra->corners; + int charge = 1; + float z_r = ref_sensor.centerZ / + std::sqrt(ref_sensor.centerX * ref_sensor.centerX + ref_sensor.centerY * ref_sensor.centerY); + float refphi = std::atan2(ref_sensor.centerY, ref_sensor.centerX); + unsigned short ref_layer = ref_sensor.extra->layer; + auto ref_location = ref_sensor.extra->location; + MatrixF4x3 next_layer_bound_points; + + for (int i = 0; i < bounds.rows(); i++) { + auto helix_p10 = Helix(ptCut, 0, 0, 10, bounds(i, 1), bounds(i, 2), bounds(i, 0), -charge); + auto helix_m10 = Helix(ptCut, 0, 0, -10, bounds(i, 1), bounds(i, 2), bounds(i, 0), -charge); + auto helix_p10_pos = Helix(ptCut, 0, 0, 10, bounds(i, 1), bounds(i, 2), bounds(i, 0), charge); + auto helix_m10_pos = Helix(ptCut, 0, 0, -10, bounds(i, 1), bounds(i, 2), bounds(i, 0), charge); + float bound_z_r = bounds(i, 0) / std::sqrt(bounds(i, 1) * bounds(i, 1) + bounds(i, 2) * bounds(i, 2)); + float bound_phi = std::atan2(bounds(i, 2), bounds(i, 1)); + float phi_diff = phi_mpi_pi(bound_phi - refphi); + + std::tuple next_point; + if (ref_location == Location::barrel) { + if (doR) { + float tar_layer_radius = average_r_barrel[ref_layer]; + if (bound_z_r < z_r) { + auto const& h = phi_diff > 0 ? helix_p10 : helix_p10_pos; + next_point = h.pointFromRadius(tar_layer_radius); + } else { + auto const& h = phi_diff > 0 ? helix_m10 : helix_m10_pos; + next_point = h.pointFromRadius(tar_layer_radius); + } + } else { + float tar_layer_z = average_z_endcap[0]; + if (bound_z_r < z_r) { + if (phi_diff > 0) { + next_point = helix_p10.pointFromZ(std::copysign(tar_layer_z, helix_p10.lambda)); + } else { + next_point = helix_p10_pos.pointFromZ(std::copysign(tar_layer_z, helix_p10_pos.lambda)); + } + } else { + if (phi_diff > 0) { + next_point = helix_m10.pointFromZ(std::copysign(tar_layer_z, helix_m10.lambda)); + } else { + next_point = helix_m10_pos.pointFromZ(std::copysign(tar_layer_z, helix_m10_pos.lambda)); + } + } + } + } else { + float tar_layer_z = average_z_endcap[ref_layer]; + if (bound_z_r < z_r) { + if (phi_diff > 0) { + next_point = helix_p10.pointFromZ(std::copysign(tar_layer_z, helix_p10.lambda)); + } else { + next_point = helix_p10_pos.pointFromZ(std::copysign(tar_layer_z, helix_p10_pos.lambda)); + } + } else { + if (phi_diff > 0) { + next_point = helix_m10.pointFromZ(std::copysign(tar_layer_z, helix_m10.lambda)); + } else { + next_point = helix_m10_pos.pointFromZ(std::copysign(tar_layer_z, helix_m10_pos.lambda)); + } + } + } + next_layer_bound_points(i, 0) = std::get<2>(next_point); + next_layer_bound_points(i, 1) = std::get<0>(next_point); + next_layer_bound_points(i, 2) = std::get<1>(next_point); + } + + return next_layer_bound_points; + } + + std::vector getCurvedLineConnections(unsigned int ref_detid, + Sensors const& sensors, + BinnedDetIds const& binned_detids, + std::array const& average_r_barrel, + std::array const& average_z_endcap, + float ptCut) { + auto const& ref_sensor = sensors.at(ref_detid); + + float refphi = std::atan2(ref_sensor.centerY, ref_sensor.centerX); + + unsigned short ref_layer = ref_sensor.extra->layer; + auto ref_location = ref_sensor.extra->location; + + auto thetaphibins = getThetaPhiBins(ref_sensor.extra->centerTheta, ref_sensor.centerPhi); + + auto const& tar_detids_to_be_considered = + binned_detids.at({ref_location, ref_layer + 1, thetaphibins.first, thetaphibins.second}); + + auto next_layer_bound_points = boundsAfterCurved(ref_detid, sensors, average_r_barrel, average_z_endcap, ptCut); + + std::vector list_of_detids_etaphi_layer_tar; + for (unsigned int tar_detid : tar_detids_to_be_considered) { + auto const& tar_sensor = sensors.at(tar_detid); + if (moduleOverlapsInEtaPhi(next_layer_bound_points, + tar_sensor.extra->corners, + refphi, + 0, + std::numeric_limits::max(), + tar_sensor.centerPhi)) + list_of_detids_etaphi_layer_tar.push_back(tar_detid); + } + + // Consider barrel to endcap connections if the intersection area is > 0 + // We construct the reference polygon as a vector of polygons because the boost::geometry::difference + // function can return multiple polygons if the difference results in disjoint pieces + if (ref_location == Location::barrel) { + std::vector barrel_endcap_connected_tar_detids; + + int zshift = 0; + + std::vector ref_polygon; + ref_polygon.push_back(getEtaPhiPolygon(next_layer_bound_points, refphi, zshift)); + + // Check whether there is still significant non-zero area + for (unsigned int tar_detid : list_of_detids_etaphi_layer_tar) { + if (ref_polygon.empty()) + break; + Polygon tar_polygon = getEtaPhiPolygon(sensors.at(tar_detid).extra->corners, refphi, zshift); + + std::vector difference; + for (auto const& ref_polygon_piece : ref_polygon) { + std::vector tmp_difference; + boost::geometry::difference(ref_polygon_piece, tar_polygon, tmp_difference); + difference.insert(difference.end(), tmp_difference.begin(), tmp_difference.end()); + } + + ref_polygon = std::move(difference); + } + + float area = 0.; + for (auto const& ref_polygon_piece : ref_polygon) + area += boost::geometry::area(ref_polygon_piece); + + if (area > 5e-3) { + auto const& new_tar_detids_to_be_considered = + binned_detids.at({Location::endcap, 1, thetaphibins.first, thetaphibins.second}); + + for (unsigned int tar_detid : new_tar_detids_to_be_considered) { + auto const& sensor_target = sensors.at(tar_detid); + float tarphi = std::atan2(sensor_target.centerY, sensor_target.centerX); + + if (std::fabs(phi_mpi_pi(tarphi - refphi)) > std::numbers::pi_v / 2.) + continue; + + Polygon tar_polygon = getEtaPhiPolygon(sensor_target.extra->corners, refphi, zshift); + + bool intersects = false; + for (auto const& ref_polygon_piece : ref_polygon) { + if (boost::geometry::intersects(ref_polygon_piece, tar_polygon)) { + intersects = true; + break; + } + } + + if (intersects) + barrel_endcap_connected_tar_detids.push_back(tar_detid); + } + } + + list_of_detids_etaphi_layer_tar.insert(list_of_detids_etaphi_layer_tar.end(), + barrel_endcap_connected_tar_detids.begin(), + barrel_endcap_connected_tar_detids.end()); + } + + return list_of_detids_etaphi_layer_tar; + } + + ModuleMap mergeLineConnections( + std::initializer_list>*> connections_list) { + ModuleMap merged; + + for (auto* connections : connections_list) { + for (auto const& [detid, list] : *connections) { + auto& target = merged[detid]; + target.insert(target.end(), list.begin(), list.end()); + } + } + + for (auto& [detid, list] : merged) { + std::sort(list.begin(), list.end()); + list.erase(std::unique(list.begin(), list.end()), list.end()); + } + + return merged; + } + + ModuleMap buildModuleMap(Sensors const& sensors, + BinnedDetIds const& binned_detids, + std::array const& average_r_barrel, + std::array const& average_z_endcap, + float pt_cut) { + std::unordered_map> straight_line_connections; + std::unordered_map> curved_line_connections; + + for (auto const& [ref_detid, s] : sensors) { + // exclude the outermost modules that do not have connections to other modules + if (!((s.extra->location == Location::barrel && s.extra->lower && s.extra->layer != 6) || + (s.extra->location == Location::endcap && s.extra->lower && s.extra->layer != 5 && + !(s.extra->ring == 15 && s.extra->layer == 1) && !(s.extra->ring == 15 && s.extra->layer == 2) && + !(s.extra->ring == 12 && s.extra->layer == 3) && !(s.extra->ring == 12 && s.extra->layer == 4)))) + continue; + straight_line_connections[ref_detid] = getStraightLineConnections(ref_detid, sensors, binned_detids); + curved_line_connections[ref_detid] = + getCurvedLineConnections(ref_detid, sensors, binned_detids, average_r_barrel, average_z_endcap, pt_cut); + } + return mergeLineConnections({&straight_line_connections, &curved_line_connections}); + } + +} // namespace lstgeometry diff --git a/RecoTracker/LSTGeometry/src/PixelMap.cc b/RecoTracker/LSTGeometry/src/PixelMap.cc new file mode 100644 index 0000000000000..d91e3e3ccef2a --- /dev/null +++ b/RecoTracker/LSTGeometry/src/PixelMap.cc @@ -0,0 +1,137 @@ +#include +#include + +#include "RecoTracker/LSTGeometry/interface/Common.h" +#include "RecoTracker/LSTGeometry/interface/PixelMap.h" + +namespace lstgeometry { + + PixelMap buildPixelMap(Sensors const& sensors, float pt_cut) { + // Charge 0 is the union of charge 1 and charge -1 + PixelMap maps; + + std::size_t nSuperbin = kPtBounds.size() * kNPhi * kNEta * kNZ; + + // Initialize empty lists for the pixel map + for (unsigned int layer : {1, 2}) { + for (unsigned int subdet : {SubDet::Barrel, SubDet::Endcap}) { + for (int charge : {-1, 0, 1}) { + maps.try_emplace({layer, subdet, charge}, nSuperbin); + } + } + } + + // Loop over the detids and for each detid compute which superbins it is connected to + for (auto const& [detId, sensor] : sensors) { + // Phase-2 enum differs from the legacy one used here + unsigned int subdet = sensor.extra->subdet == SubDetector::P2OTB ? SubDet::Barrel : SubDet::Endcap; + auto layer = sensor.extra->layer; + if (layer > 2) + continue; + auto location = sensor.extra->location; + + // Skip if the module is not PS module and is not lower sensor + if (sensor.moduleType == ModuleType::Ph2SS || !sensor.extra->lower) + continue; + + // For this module, now compute which super bins they belong to + // To compute which super bins it belongs to, one needs to provide at least pt and z window to compute compatible eta and phi range + // So we have a loop in pt and Z + for (unsigned int ipt = 0; ipt < kPtBounds.size(); ipt++) { + for (unsigned int iz = 0; iz < kNZ; iz++) { + // The zmin, zmax of consideration + float zmin = -30 + iz * (60. / kNZ); + float zmax = -30 + (iz + 1) * (60. / kNZ); + + zmin -= 0.05; + zmax += 0.05; + + // The ptmin, ptmax of consideration + float pt_lo = ipt == 0 ? pt_cut : kPtBounds[ipt - 1]; + float pt_hi = kPtBounds[ipt]; + + auto [etamin, etamax] = getCompatibleEtaRange(sensor, zmin, zmax); + + etamin -= 0.05; + etamax += 0.05; + + if (layer == 2 && location == Location::endcap) { + if (etamax < 2.3) + continue; + if (etamin < 2.3) + etamin = 2.3; + } + + // Compute the indices of the compatible eta range + unsigned int ietamin = static_cast(std::max((etamin + 2.6f) / (5.2f / kNEta), 0.0f)); + unsigned int ietamax = + static_cast(std::min((etamax + 2.6f) / (5.2f / kNEta), static_cast(kNEta - 1))); + + auto phi_ranges = getCompatiblePhiRange(sensor, pt_lo, pt_hi); + + unsigned int iphimin_pos = static_cast((phi_ranges.first.first + std::numbers::pi_v) / + (2. * std::numbers::pi_v / kNPhi)); + unsigned int iphimax_pos = static_cast((phi_ranges.first.second + std::numbers::pi_v) / + (2. * std::numbers::pi_v / kNPhi)); + unsigned int iphimin_neg = static_cast((phi_ranges.second.first + std::numbers::pi_v) / + (2. * std::numbers::pi_v / kNPhi)); + unsigned int iphimax_neg = static_cast((phi_ranges.second.second + std::numbers::pi_v) / + (2. * std::numbers::pi_v / kNPhi)); + + // <= to cover some inefficiencies + for (unsigned int ieta = ietamin; ieta <= ietamax; ieta++) { + // if the range is crossing the -pi v. pi boundary special care is needed + if (iphimin_pos <= iphimax_pos) { + for (unsigned int iphi = iphimin_pos; iphi < iphimax_pos; iphi++) { + unsigned int isuperbin = (ipt * kNPhi * kNEta * kNZ) + (ieta * kNPhi * kNZ) + (iphi * kNZ) + iz; + maps[{layer, subdet, 1}][isuperbin].push_back(detId); + maps[{layer, subdet, 0}][isuperbin].push_back(detId); + } + } else { + for (unsigned int iphi = 0; iphi < iphimax_pos; iphi++) { + unsigned int isuperbin = (ipt * kNPhi * kNEta * kNZ) + (ieta * kNPhi * kNZ) + (iphi * kNZ) + iz; + maps[{layer, subdet, 1}][isuperbin].push_back(detId); + maps[{layer, subdet, 0}][isuperbin].push_back(detId); + } + for (unsigned int iphi = iphimin_pos; iphi < kNPhi; iphi++) { + unsigned int isuperbin = (ipt * kNPhi * kNEta * kNZ) + (ieta * kNPhi * kNZ) + (iphi * kNZ) + iz; + maps[{layer, subdet, 1}][isuperbin].push_back(detId); + maps[{layer, subdet, 0}][isuperbin].push_back(detId); + } + } + if (iphimin_neg <= iphimax_neg) { + for (unsigned int iphi = iphimin_neg; iphi < iphimax_neg; iphi++) { + unsigned int isuperbin = (ipt * kNPhi * kNEta * kNZ) + (ieta * kNPhi * kNZ) + (iphi * kNZ) + iz; + maps[{layer, subdet, -1}][isuperbin].push_back(detId); + maps[{layer, subdet, 0}][isuperbin].push_back(detId); + } + } else { + for (unsigned int iphi = 0; iphi < iphimax_neg; iphi++) { + unsigned int isuperbin = (ipt * kNPhi * kNEta * kNZ) + (ieta * kNPhi * kNZ) + (iphi * kNZ) + iz; + maps[{layer, subdet, -1}][isuperbin].push_back(detId); + maps[{layer, subdet, 0}][isuperbin].push_back(detId); + } + for (unsigned int iphi = iphimin_neg; iphi < kNPhi; iphi++) { + unsigned int isuperbin = (ipt * kNPhi * kNEta * kNZ) + (ieta * kNPhi * kNZ) + (iphi * kNZ) + iz; + maps[{layer, subdet, -1}][isuperbin].push_back(detId); + maps[{layer, subdet, 0}][isuperbin].push_back(detId); + } + } + } + } + } + } + + for (auto& [key, vec_of_vecs] : maps) { + for (auto& vec : vec_of_vecs) { + if (vec.empty()) + continue; + std::sort(vec.begin(), vec.end()); + vec.erase(std::unique(vec.begin(), vec.end()), vec.end()); + } + } + + return maps; + } + +} // namespace lstgeometry diff --git a/RecoTracker/LSTGeometry/src/Sensor.cc b/RecoTracker/LSTGeometry/src/Sensor.cc new file mode 100644 index 0000000000000..18a2c664b66d0 --- /dev/null +++ b/RecoTracker/LSTGeometry/src/Sensor.cc @@ -0,0 +1,142 @@ +#include + +#include "DataFormats/GeometrySurface/interface/RectangularPlaneBounds.h" + +#include "RecoTracker/LSTGeometry/interface/Sensor.h" + +using namespace lstgeometry; + +// Not sure if there is functionality for this already in CMSSW +bool isInverted(unsigned int moduleId, Location location, Side side, unsigned int layer) { + bool moduleIdIsEven = moduleId % 2 == 0; + if (location == Location::endcap) { + if (side == Side::NegZ) { + return !moduleIdIsEven; + } else if (side == Side::PosZ) { + return moduleIdIsEven; + } + } else if (location == Location::barrel) { + if (side == Side::Center) { + if (layer <= 3) { + return !moduleIdIsEven; + } else if (layer >= 4) { + return moduleIdIsEven; + } + } else if (side == Side::NegZ || side == Side::PosZ) { + if (layer <= 2) { + return !moduleIdIsEven; + } else if (layer == 3) { + return moduleIdIsEven; + } + } + } + return false; +} + +// This differs from TrackerTopology::isLower since it considers if the module is inverted +bool isLower(unsigned int moduleId, Location location, Side side, unsigned int layer, unsigned int detId) { + return isInverted(moduleId, location, side, layer) ? !(detId & 1) : (detId & 1); +} + +bool isStrip(ModuleType moduleType, bool isInverted, bool isLower) { + if (moduleType == ModuleType::Ph2SS) + return true; + if (isInverted) + return isLower; + return !isLower; +} + +Sensor::Sensor(unsigned int detId, + ModuleType moduleType, + SubDetector subdet, + Location location, + Side side, + unsigned int moduleId, + unsigned int layer, + unsigned int ring, + float centerRho, + float centerZ, + float centerPhi, + Surface const &surface) + : moduleType(moduleType), + centerPhi(centerPhi), + centerX(centerRho * std::cos(centerPhi)), + centerY(centerRho * std::sin(centerPhi)), + centerZ(centerZ), + extra(std::make_unique()) { + extra->subdet = subdet; + extra->location = location; + extra->side = side; + extra->moduleId = moduleId; + extra->layer = layer; + extra->ring = ring; + extra->inverted = isInverted(moduleId, location, side, layer); + extra->centerRho = centerRho; + extra->lower = isLower(moduleId, location, side, layer, detId); + extra->strip = isStrip(moduleType, extra->inverted, extra->lower); + extra->centerEta = std::asinh(centerZ / centerRho); + extra->centerTheta = std::numbers::pi_v / 2. - std::atan(centerZ / centerRho); + + const Bounds &bounds = surface.bounds(); + const RectangularPlaneBounds *rectangular_bounds = dynamic_cast(&bounds); + + float wid, len; + if (rectangular_bounds) { + wid = rectangular_bounds->width(); + len = rectangular_bounds->length(); + } else { + throw std::runtime_error("Sensor::Sensor: Surface bounds are not rectangular"); + } + + auto c1 = GloballyPositioned::LocalPoint(-wid / 2, -len / 2, 0); + auto c2 = GloballyPositioned::LocalPoint(-wid / 2, len / 2, 0); + auto c3 = GloballyPositioned::LocalPoint(wid / 2, len / 2, 0); + auto c4 = GloballyPositioned::LocalPoint(wid / 2, -len / 2, 0); + auto c1g = surface.toGlobal(c1); + auto c2g = surface.toGlobal(c2); + auto c3g = surface.toGlobal(c3); + auto c4g = surface.toGlobal(c4); + // store corners with z, x, y ordering + extra->corners << c1g.z(), c1g.x(), c1g.y(), c2g.z(), c2g.x(), c2g.y(), c3g.z(), c3g.x(), c3g.y(), c4g.z(), c4g.x(), + c4g.y(); + + // Precompute min/max R, Z, Phi for convenience + extra->minR = std::numeric_limits::max(); + extra->maxR = std::numeric_limits::lowest(); + extra->minZ = std::numeric_limits::max(); + extra->maxZ = std::numeric_limits::lowest(); + extra->minPhi = std::numeric_limits::max(); + float minPosPhi = std::numeric_limits::max(); + extra->maxPhi = std::numeric_limits::lowest(); + float maxNegPhi = std::numeric_limits::lowest(); + unsigned int nPos = 0; + unsigned int nOverPi2 = 0; + for (int i = 0; i < extra->corners.rows(); i++) { + float x = extra->corners(i, 1); + float y = extra->corners(i, 2); + float r = std::sqrt(x * x + y * y); + float z = extra->corners(i, 0); + float phi = phi_mpi_pi(std::numbers::pi_v + std::atan2(-extra->corners(i, 2), -extra->corners(i, 1))); + + extra->minR = std::min(extra->minR, r); + extra->maxR = std::max(extra->maxR, r); + extra->minZ = std::min(extra->minZ, z); + extra->maxZ = std::max(extra->maxZ, z); + extra->minPhi = std::min(extra->minPhi, phi); + extra->maxPhi = std::max(extra->maxPhi, phi); + + if (phi > 0) { + minPosPhi = std::min(minPosPhi, phi); + nPos++; + } else { + maxNegPhi = std::max(maxNegPhi, phi); + } + if (std::fabs(phi) > std::numbers::pi_v / 2.) { + nOverPi2++; + } + } + if (nOverPi2 == 4 && nPos != 4 && nPos != 0) { + extra->minPhi = minPosPhi; + extra->maxPhi = maxNegPhi; + } +} diff --git a/RecoTracker/LSTGeometry/src/SensorBinning.cc b/RecoTracker/LSTGeometry/src/SensorBinning.cc new file mode 100644 index 0000000000000..8c4d5deb525e1 --- /dev/null +++ b/RecoTracker/LSTGeometry/src/SensorBinning.cc @@ -0,0 +1,101 @@ +#include "RecoTracker/LSTGeometry/interface/SensorBinning.h" + +namespace lstgeometry { + + bool isInThetaPhiBin(float theta, float phi, unsigned int theta_bin, unsigned int phi_bin) { + if (theta_bin == 0) { + if (theta > 3. * kThetaBinRad / 2.) + return false; + } else if (theta_bin == kNThetaBins - 1) { + if (theta < (2 * (kNThetaBins - 1) - 1) * kThetaBinRad / 2.) + return false; + } else if (theta < (2 * theta_bin - 1) * kThetaBinRad / 2. || + theta > (2 * (theta_bin + 1) + 1) * kThetaBinRad / 2.) { + return false; + } + + float pi = std::numbers::pi_v; + if (phi_bin == 0) { + if (phi > -pi + kPhiBinWidth && phi < pi - kPhiBinWidth) + return false; + } else { + if (phi < -pi + (phi_bin - 1) * kPhiBinWidth || phi > -pi + (phi_bin + 1) * kPhiBinWidth) + return false; + } + + return true; + } + + std::pair getThetaPhiBins(float theta, float phi) { + unsigned int theta_bin = std::clamp(static_cast(theta / kThetaBinRad), 0u, kNThetaBins - 1); + + float pi = std::numbers::pi_v; + unsigned int phi_bin = + std::clamp(static_cast((phi + pi + kPhiBinWidth / 2.) / kPhiBinWidth), 0u, kNPhiBins - 1); + if (phi >= pi - kPhiBinWidth / 2) + phi_bin = 0; // The 0 bin wraps around + + return std::make_pair(theta_bin, phi_bin); + } + + std::pair getCompatibleEtaRange(Sensor const& sensor, float zmin_bound, float zmax_bound) { + float minr = sensor.extra->minR; + float maxr = sensor.extra->maxR; + float minz = sensor.extra->minZ; + float maxz = sensor.extra->maxZ; + float mineta = -std::log(std::tan(std::atan2(minz > 0 ? maxr : minr, minz - zmin_bound) / 2.)); + float maxeta = -std::log(std::tan(std::atan2(maxz > 0 ? minr : maxr, maxz - zmax_bound) / 2.)); + + if (maxeta < mineta) + std::swap(maxeta, mineta); + return std::make_pair(mineta, maxeta); + } + + std::pair, std::pair> getCompatiblePhiRange(Sensor const& sensor, + float ptmin, + float ptmax) { + float minr = sensor.extra->minR; + float maxr = sensor.extra->maxR; + float minphi = sensor.extra->minPhi; + float maxphi = sensor.extra->maxPhi; + float A = kC * kB / 2.; + float pos_q_phi_lo_bound = phi_mpi_pi(A * minr / ptmax + minphi); + float pos_q_phi_hi_bound = phi_mpi_pi(A * maxr / ptmin + maxphi); + float neg_q_phi_lo_bound = phi_mpi_pi(-A * maxr / ptmin + minphi); + float neg_q_phi_hi_bound = phi_mpi_pi(-A * minr / ptmax + maxphi); + return std::make_pair(std::make_pair(pos_q_phi_lo_bound, pos_q_phi_hi_bound), + std::make_pair(neg_q_phi_lo_bound, neg_q_phi_hi_bound)); + } + + BinnedDetIds binDetIds(Sensors const& sensors) { + BinnedDetIds binned_detids; + + // Initialize all vectors + for (unsigned int thetabin = 0; thetabin < kNThetaBins; thetabin++) { + for (unsigned int phibin = 0; phibin < kNPhiBins; phibin++) { + for (unsigned int layer = 1; layer <= kBarrelLayers; layer++) { + binned_detids[{Location::barrel, layer, thetabin, phibin}] = {}; + } + for (unsigned int layer = 1; layer <= kEndcapLayers; layer++) { + binned_detids[{Location::endcap, layer, thetabin, phibin}] = {}; + } + } + } + + for (auto const& [detid, sensor] : sensors) { + if (!sensor.extra->lower) + continue; + // We need to loop over all bins instead of using getThetaPhiBins because of the overlapping bins + for (unsigned int thetabin = 0; thetabin < kNThetaBins; thetabin++) { + for (unsigned int phibin = 0; phibin < kNPhiBins; phibin++) { + if (isInThetaPhiBin(sensor.extra->centerTheta, sensor.centerPhi, thetabin, phibin)) { + binned_detids[{sensor.extra->location, sensor.extra->layer, thetabin, phibin}].push_back(detid); + } + } + } + } + + return binned_detids; + } + +} // namespace lstgeometry diff --git a/RecoTracker/LSTGeometry/src/Slope.cc b/RecoTracker/LSTGeometry/src/Slope.cc new file mode 100644 index 0000000000000..e655398adc99d --- /dev/null +++ b/RecoTracker/LSTGeometry/src/Slope.cc @@ -0,0 +1,42 @@ +#include +#include + +#include "RecoTracker/LSTGeometry/interface/Common.h" +#include "RecoTracker/LSTGeometry/interface/Slope.h" + +namespace lstgeometry { + + Slope::Slope(float dx, float dy, float dz) { + float dr = sqrt(dx * dx + dy * dy); + drdz = dz != 0 ? dr / dz : kDefaultSlope; + dxdy = dy != 0 ? -dx / dy : kDefaultSlope; + } + + // Use each sensor's corners to calculate and categorize drdz and dxdy slopes. + std::tuple computeSlopes(Sensors const& sensors) { + Slopes barrel_slopes; + Slopes endcap_slopes; + + for (auto const& [detId, sensor] : sensors) { + float dx = roundCoordinate(sensor.extra->corners(1, 1) - sensor.extra->corners(0, 1)); + float dy = roundCoordinate(sensor.extra->corners(1, 2) - sensor.extra->corners(0, 2)); + float dz = roundCoordinate(sensor.extra->corners(1, 0) - sensor.extra->corners(0, 0)); + + Slope slope(dx, dy, dz); + + auto location = sensor.extra->location; + bool is_tilted = sensor.extra->side != Side::Center; + bool is_strip = sensor.extra->strip; + + if (!is_strip) + continue; + + if (location == Location::barrel and is_tilted) + barrel_slopes[detId] = slope; + else if (location == Location::endcap) + endcap_slopes[detId] = slope; + } + + return std::make_tuple(barrel_slopes, endcap_slopes); + } +} // namespace lstgeometry diff --git a/RecoTracker/LSTGeometry/test/BuildFile.xml b/RecoTracker/LSTGeometry/test/BuildFile.xml new file mode 100644 index 0000000000000..339b55f111d07 --- /dev/null +++ b/RecoTracker/LSTGeometry/test/BuildFile.xml @@ -0,0 +1,7 @@ + + + + + + + diff --git a/RecoTracker/LSTGeometry/test/DumpLSTGeometry.cc b/RecoTracker/LSTGeometry/test/DumpLSTGeometry.cc new file mode 100644 index 0000000000000..72988740c3108 --- /dev/null +++ b/RecoTracker/LSTGeometry/test/DumpLSTGeometry.cc @@ -0,0 +1,49 @@ +#include "FWCore/Framework/interface/one/EDAnalyzer.h" +#include "FWCore/Framework/interface/MakerMacros.h" +#include "FWCore/Framework/interface/ESTransientHandle.h" +#include "FWCore/Framework/interface/ESHandle.h" +#include "FWCore/Framework/interface/EventSetup.h" + +#include "RecoTracker/Record/interface/TrackerRecoGeometryRecord.h" + +#include "RecoTracker/LSTGeometry/interface/Common.h" +#include "RecoTracker/LSTGeometry/interface/Geometry.h" +#include "RecoTracker/LSTGeometry/interface/IO.h" + +class DumpLSTGeometry : public edm::one::EDAnalyzer<> { +public: + explicit DumpLSTGeometry(const edm::ParameterSet& config); + +private: + void analyze(const edm::Event& event, const edm::EventSetup& eventSetup) override; + + double ptCut_; + std::string ptCutStr_; + std::string outputDirectory_; + bool binaryOutput_; + + edm::ESGetToken lstGeoToken_; +}; + +DumpLSTGeometry::DumpLSTGeometry(const edm::ParameterSet& config) + : ptCut_(config.getParameter("ptCut")), + ptCutStr_(lst::floatToStr(ptCut_, 1)), + outputDirectory_(config.getUntrackedParameter("outputDirectory", "data/")), + binaryOutput_(config.getUntrackedParameter("outputAsBinary", true)), + lstGeoToken_{esConsumes(edm::ESInputTag("", ptCutStr_))} {} + +void DumpLSTGeometry::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) { + auto const& lstg = iSetup.getData(lstGeoToken_); + + lstgeometry::writeSensorCentroids(lstg.sensors, outputDirectory_ + "sensor_centroids", binaryOutput_); + lstgeometry::writeSlopes( + lstg.barrel_slopes, lstg.sensors, outputDirectory_ + "tilted_barrel_orientation", binaryOutput_); + lstgeometry::writeSlopes(lstg.endcap_slopes, lstg.sensors, outputDirectory_ + "endcap_orientation", binaryOutput_); + lstgeometry::writePixelMaps(lstg.pixel_map, outputDirectory_ + "pixelmap/pLS_map", binaryOutput_); + lstgeometry::writeModuleConnections( + lstg.module_map, outputDirectory_ + "module_connection_tracing_merged", binaryOutput_); + + edm::LogInfo("DumpLSTGeometry") << "Geometry data was successfully dumped." << std::endl; +} + +DEFINE_FWK_MODULE(DumpLSTGeometry); diff --git a/RecoTracker/LSTGeometry/test/dumpLSTGeometry.py b/RecoTracker/LSTGeometry/test/dumpLSTGeometry.py new file mode 100644 index 0000000000000..2cc7b45c2ae10 --- /dev/null +++ b/RecoTracker/LSTGeometry/test/dumpLSTGeometry.py @@ -0,0 +1,66 @@ +from argparse import ArgumentParser + +import FWCore.ParameterSet.Config as cms +from Configuration.AlCa.GlobalTag import GlobalTag +import Configuration.Geometry.defaultPhase2ConditionsEra_cff as _settings +from RecoTracker.LSTGeometry.lstGeometryESProducer_cfi import lstGeometryESProducer as _lstGeom + +trackingLSTCommon = cms.Modifier() +trackingLST = cms.ModifierChain(trackingLSTCommon) + +################################################################### +# Set default phase-2 settings +################################################################### +_PH2_GLOBAL_TAG, _PH2_ERA = _settings.get_era_and_conditions(_settings.DEFAULT_VERSION) + +# No era in Fireworks/Geom reco dumper +process = cms.Process("DUMP", _PH2_ERA, trackingLST) + +# import of standard configurations +process.load("Configuration.StandardSequences.Services_cff") +process.load("FWCore.MessageService.MessageLogger_cfi") +process.load("Configuration.Geometry.GeometryExtendedRun4DefaultReco_cff") +process.load("Configuration.StandardSequences.MagneticField_cff") +process.load("Configuration.StandardSequences.Reconstruction_cff") +process.load("Configuration.StandardSequences.FrontierConditions_GlobalTag_cff") + +process.GlobalTag = GlobalTag(process.GlobalTag, _PH2_GLOBAL_TAG, "") + +process.MessageLogger.cerr.threshold = "INFO" +process.MessageLogger.cerr.LSTGeometryESProducer = dict(limit=-1) + +process.source = cms.Source("EmptySource") +process.maxEvents.input = 1 + +parser = ArgumentParser() +parser.add_argument( + "--outputDirectory", + default="data/", + help="Output directory for LST geometry files" +) +parser.add_argument( + "--ptCut", + type=float, + default=0.8, + help="pT cut for LST module maps" +) +parser.add_argument( + "--binaryOutput", + action="store_true", + help="Dump LST geometry as binary files" +) +options = parser.parse_args() + +process.dump = cms.EDAnalyzer( + "DumpLSTGeometry", + outputDirectory = cms.untracked.string(options.outputDirectory + "/"), + ptCut = cms.double(options.ptCut), + outputAsBinary = cms.untracked.bool(options.binaryOutput), +) + +process.lstGeometryESProducer = _lstGeom.clone(ptCut = cms.double(options.ptCut)) +process.dTask = cms.Task(process.lstGeometryESProducer) +process.dSeq = cms.Sequence(process.dump,process.dTask) +process.p = cms.Path(process.dSeq) + +print(f"Requesting LST geometry dump into directory: {options.outputDirectory}\n") diff --git a/RecoTracker/LSTGeometry/test/testDumpLSTGeometry.sh b/RecoTracker/LSTGeometry/test/testDumpLSTGeometry.sh new file mode 100755 index 0000000000000..33a0762f58de0 --- /dev/null +++ b/RecoTracker/LSTGeometry/test/testDumpLSTGeometry.sh @@ -0,0 +1,10 @@ +#!/bin/sh + +function die { echo $1: status $2; exit $2; } + +if [ "${SCRAM_TEST_NAME}" != "" ] ; then + mkdir ${SCRAM_TEST_NAME} + cd ${SCRAM_TEST_NAME} +fi + +(cmsRun ${SCRAM_TEST_PATH}/dumpLSTGeometry.py --outputDirectory data --ptCut 0.8 --binaryOutput) || die "failed to run dumpLSTGeometry.py" $?