diff --git a/RecoTracker/LST/plugins/alpaka/LSTModulesDevESProducer.cc b/RecoTracker/LST/plugins/alpaka/LSTModulesDevESProducer.cc index 7152da9ed13c7..47771e7041d20 100644 --- a/RecoTracker/LST/plugins/alpaka/LSTModulesDevESProducer.cc +++ b/RecoTracker/LST/plugins/alpaka/LSTModulesDevESProducer.cc @@ -1,5 +1,8 @@ +#include + // LST includes #include "RecoTracker/LSTCore/interface/alpaka/LST.h" +#include "RecoTracker/LSTGeometry/interface/Geometry.h" #include "FWCore/ParameterSet/interface/ParameterSet.h" @@ -14,22 +17,29 @@ 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::ostringstream ptCutOSS; + ptCutOSS << std::setprecision(1) << ptCut_; + std::string ptCutStr = ptCutOSS.str(); + + auto cc = setWhatProduced(this, "LSTModuleMaps_" + ptCutStr); + lstGeoToken_ = cc.consumes(edm::ESInputTag("", "LSTGeometry_" + 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..bf7a5658196a0 100644 --- a/RecoTracker/LST/plugins/alpaka/LSTProducer.cc +++ b/RecoTracker/LST/plugins/alpaka/LSTProducer.cc @@ -1,5 +1,7 @@ #include +#include + #include "RecoTracker/LSTCore/interface/alpaka/LST.h" #include "FWCore/MessageLogger/interface/MessageLogger.h" @@ -25,13 +27,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_(getPtCutStr(ptCut_)), clustSizeCut_(static_cast(config.getParameter("clustSizeCut"))), nopLSDupClean_(config.getParameter("nopLSDupClean")), tcpLSTriplets_(config.getParameter("tcpLSTriplets")), + lstInputToken_{consumes(config.getParameter("lstInput"))}, + lstESToken_{esConsumes(edm::ESInputTag("", "LSTModuleMaps_" + ptCutStr_))}, lstOutputToken_{produces()} {} void produce(edm::StreamID sid, device::Event& iEvent, const device::EventSetup& iSetup) const override { @@ -42,7 +45,7 @@ namespace ALPAKA_ACCELERATOR_NAMESPACE { lst.run(iEvent.queue(), verbose_, - static_cast(ptCut_), + ptCut_, clustSizeCut_, &lstESDeviceData, &lstInputDC, @@ -60,21 +63,27 @@ 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_; + + static std::string getPtCutStr(double ptCut) { + std::ostringstream ptCutOSS; + ptCutOSS << std::setprecision(1) << ptCut; + return ptCutOSS.str(); + } }; } // namespace ALPAKA_ACCELERATOR_NAMESPACE 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 e9cbb2cccf046..4c062cae70ebf 100644 --- a/RecoTracker/LSTCore/BuildFile.xml +++ b/RecoTracker/LSTCore/BuildFile.xml @@ -5,6 +5,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..fe4af4186c287 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&); void load(std::string const&); + void load(lstgeometry::Slopes const&, lstgeometry::Sensors const&); 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 5ed0546351490..f4dc33268f8fd 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; @@ -43,8 +43,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..cfdebea6bd7d2 100644 --- a/RecoTracker/LSTCore/src/EndcapGeometry.cc +++ b/RecoTracker/LSTCore/src/EndcapGeometry.cc @@ -7,6 +7,10 @@ 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(); @@ -39,6 +43,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..204940babf0e0 100644 --- a/RecoTracker/LSTCore/src/LSTESData.cc +++ b/RecoTracker/LSTCore/src/LSTESData.cc @@ -7,6 +7,8 @@ #include "ModuleMethods.h" #include +#include +#include namespace { std::string geometryDataDir() { @@ -79,7 +81,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 +121,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) { + final_modulemap[detId] = std::vector(connections.begin(), connections.end()); + } + 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& set = map.at(isuperbin); + final_pixelmap[isuperbin] = std::vector(set.begin(), set.end()); + } + + if (charge == 0) { + pLStoLayer[0][layer - 1 + (subdet == Endcap ? 2 : 0)] = lst::ModuleConnectionMap(final_pixelmap); + } else if (charge > 0) { + pLStoLayer[1][layer - 1 + (subdet == Endcap ? 2 : 0)] = lst::ModuleConnectionMap(final_pixelmap); + } else { + pLStoLayer[2][layer - 1 + (subdet == Endcap ? 2 : 0)] = 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..5ef07f3594c65 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; @@ -213,24 +216,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 +382,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 +402,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..049bddfbd34d1 100644 --- a/RecoTracker/LSTCore/src/TiltedGeometry.cc +++ b/RecoTracker/LSTCore/src/TiltedGeometry.cc @@ -7,6 +7,8 @@ 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(); @@ -37,6 +39,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 806869522e062..d6bbe192ab07e 100644 --- a/RecoTracker/LSTCore/src/alpaka/MiniDoublet.h +++ b/RecoTracker/LSTCore/src/alpaka/MiniDoublet.h @@ -304,8 +304,7 @@ namespace ALPAKA_ACCELERATOR_NAMESPACE::lst { drprime = (moduleSeparation / alpaka::math::sin(acc, angleA + angleB)) * alpaka::math::sin(acc, angleA); // Compute arctan of the slope and take care of the slope = infinity case - absArctanSlope = - ((slope != kVerticalModuleSlope && edm::isFinite(slope)) ? fabs(alpaka::math::atan(acc, slope)) : kPi / 2.f); + absArctanSlope = (edm::isFinite(slope) ? fabs(alpaka::math::atan(acc, slope)) : kPi / 2.f); // Depending on which quadrant the pixel hit lies, we define the angleM by shifting them slightly differently if (xp > 0 and yp > 0) { @@ -328,8 +327,7 @@ namespace ALPAKA_ACCELERATOR_NAMESPACE::lst { ya = yp + drprime_y; // Compute the new strip hit position (if the slope value is in special condition take care of the exceptions) - if (slope == kVerticalModuleSlope || - edm::isNotFinite(slope)) // Designated for tilted module when the slope is infinity (module lying along y-axis) + if (edm::isNotFinite(slope)) // Designated for tilted module when the slope is infinity (module lying along y-axis) { xn = xa; // New x point is simply where the anchor is yn = yo; // No shift in y 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 bd1d2825309c1..191a601bae0db 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..d4af3f35d997c --- /dev/null +++ b/RecoTracker/LSTGeometry/interface/Common.h @@ -0,0 +1,58 @@ +#ifndef RecoTracker_LSTGeometry_interface_Common_h +#define RecoTracker_LSTGeometry_interface_Common_h + +#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" + +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 }; + + 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 + +#endif diff --git a/RecoTracker/LSTGeometry/interface/DetectorGeometry.h b/RecoTracker/LSTGeometry/interface/DetectorGeometry.h new file mode 100644 index 0000000000000..2e91641026f85 --- /dev/null +++ b/RecoTracker/LSTGeometry/interface/DetectorGeometry.h @@ -0,0 +1,82 @@ +#ifndef RecoTracker_LSTGeometry_interface_DetectorGeometry_h +#define RecoTracker_LSTGeometry_interface_DetectorGeometry_h + +#include +#include +#include +#include +#include + +#include "RecoTracker/LSTGeometry/interface/Common.h" +#include "RecoTracker/LSTGeometry/interface/Sensor.h" + +namespace lstgeometry { + + using LayerEtaBinPhiBinKey = std::tuple; + + class DetectorGeometry { + private: + std::shared_ptr sensors_; + std::array avg_radii_; + std::array avg_z_; + std::unordered_map, boost::hash> + barrel_lower_det_ids_; + std::unordered_map, boost::hash> + endcap_lower_det_ids_; + + public: + DetectorGeometry(std::shared_ptr sensors, + std::array const& avg_radii, + std::array const& avg_z); +#ifndef LST_STANDALONE + MatrixF4x3 const& getCorners(unsigned int detId) const; +#endif + + std::vector getDetIds(std::function&)> filter = + [](auto const&) { return true; }) const; + + void buildByLayer(Sensors const& sensors); + + std::vector const& getBarrelLayerDetIds(unsigned int layer, + unsigned int etabin, + unsigned int phibin) const; + + std::vector const& getEndcapLayerDetIds(unsigned int layer, + unsigned int etabin, + unsigned int phibin) const; + + float getBarrelLayerAverageRadius(unsigned int layer) const; + + float getEndcapLayerAverageAbsZ(unsigned int layer) const; + + float getMinR(unsigned int detId) const; + + float getMaxR(unsigned int detId) const; + + float getMinZ(unsigned int detId) const; + + float getMaxZ(unsigned int detId) const; + + float getMinPhi(unsigned int detId) const; + + float getMaxPhi(unsigned int detId) const; + + std::pair getCompatibleEtaRange(unsigned int detId, float zmin_bound, float zmax_bound) const; + + std::pair, std::pair> getCompatiblePhiRange(unsigned int detId, + float ptmin, + float ptmax) const; + + // We split modules into overlapping eta-phi bins so that it's easier to construct module maps + // These values are just guesses and can be optimized later + static constexpr unsigned int kNEtaBins = 4; + static constexpr float kEtaBinRad = std::numbers::pi_v / kNEtaBins; + static constexpr unsigned int kNPhiBins = 6; + static constexpr float kPhiBinWidth = 2 * std::numbers::pi_v / kNPhiBins; + + static bool isInEtaPhiBin(float eta, float phi, unsigned int eta_bin, unsigned int phi_bin); + static std::pair getEtaPhiBins(float eta, float phi); + }; +} // namespace lstgeometry + +#endif diff --git a/RecoTracker/LSTGeometry/interface/Geometry.h b/RecoTracker/LSTGeometry/interface/Geometry.h new file mode 100644 index 0000000000000..ba466842618b4 --- /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/Slope.h" +#include "RecoTracker/LSTGeometry/interface/PixelMap.h" +#include "RecoTracker/LSTGeometry/interface/ModuleMap.h" +#include "RecoTracker/LSTGeometry/interface/Sensor.h" + +namespace lstgeometry { + + struct Geometry { + std::shared_ptr sensors; + Slopes barrel_slopes; + Slopes endcap_slopes; + PixelMap pixel_map; + ModuleMap module_map; + + Geometry(std::shared_ptr 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..66da3e933aa89 --- /dev/null +++ b/RecoTracker/LSTGeometry/interface/Helix.h @@ -0,0 +1,131 @@ +#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), 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 * sin(phi); + float inward_radial_vec_y = charge * radius * -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 = copysign(M_PI / 2. - 2. * atan(exp(-std::abs(eta))), 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 { + auto objective_function = [this, target_r](float t) { + float x = this->center[0] - this->charge * this->radius * std::sin(this->phi - this->charge * t); + float y = this->center[1] + this->charge * this->radius * std::cos(this->phi - this->charge * t); + return std::fabs(std::sqrt(x * x + y * y) - target_r); + }; + int bits = std::numeric_limits::digits; + auto result = boost::math::tools::brent_find_minima(objective_function, 0.0f, std::numbers::pi_v, bits); + float t = result.first; + + 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 { + auto objective_function = [this, target_z](float t) { + float z = this->center[2] + this->radius * std::tan(this->lambda) * t; + return std::fabs(z - target_z); + }; + int bits = std::numeric_limits::digits; + auto result = boost::math::tools::brent_find_minima(objective_function, 0.0f, std::numbers::pi_v, bits); + float t = result.first; + + 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..62cae47ce652d --- /dev/null +++ b/RecoTracker/LSTGeometry/interface/IO.h @@ -0,0 +1,160 @@ +#ifndef RecoTracker_LSTGeometry_interface_IO_h +#define RecoTracker_LSTGeometry_interface_IO_h + +#include +#include +#include +#include +#include +#include +#include + +#include "RecoTracker/LSTGeometry/interface/Common.h" +#include "RecoTracker/LSTGeometry/interface/Sensor.h" +#include "RecoTracker/LSTGeometry/interface/PixelMap.h" +#include "RecoTracker/LSTGeometry/interface/Slope.h" + +namespace lstgeometry { + + void writeSensorCentroids(Sensors const& sensors, std::string const& base_filename, bool binary = true) { + 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 = true) { + 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(std::unordered_map> const& connections, + std::string const& base_filename, + bool binary = true) { + 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, set] : connections) { + file.write(reinterpret_cast(&detid), sizeof(detid)); + unsigned int length = set.size(); + file.write(reinterpret_cast(&length), sizeof(length)); + for (unsigned int i : set) { + file.write(reinterpret_cast(&i), sizeof(i)); + } + } + } else { + for (auto const& [detid, set] : connections) { + file << detid << "," << set.size(); + for (unsigned int i : set) { + file << "," << i; + } + file << std::endl; + } + } + } + + void writePixelMaps(PixelMap const& maps, std::string const& base_filename, bool binary = true) { + std::filesystem::path filepath(base_filename); + std::filesystem::create_directories(filepath.parent_path()); + + if (binary) { + for (auto const& [layersubdetcharge, map] : maps) { + auto const& [layer, subdet, charge] = layersubdetcharge; + + std::string charge_str = charge > 0 ? "_pos" : (charge < 0 ? "_neg" : ""); + std::string filename = std::format("{}{}_layer{}_subdet{}.bin", base_filename, charge_str, layer, subdet); + + std::ofstream file(filename, std::ios::binary); + + for (unsigned int isuperbin = 0; isuperbin < map.size(); isuperbin++) { + auto const& set = map.at(isuperbin); + + file.write(reinterpret_cast(&isuperbin), sizeof(isuperbin)); + unsigned int length = set.size(); + file.write(reinterpret_cast(&length), sizeof(length)); + for (unsigned int i : set) { + file.write(reinterpret_cast(&i), sizeof(i)); + } + } + } + } else { + for (auto const& [layersubdetcharge, map] : maps) { + auto const& [layer, subdet, charge] = layersubdetcharge; + + std::string charge_str = charge > 0 ? "_pos" : (charge < 0 ? "_neg" : ""); + std::string filename = std::format("{}{}_layer{}_subdet{}.txt", base_filename, charge_str, layer, subdet); + + std::ofstream file(filename); + + for (unsigned int isuperbin = 0; isuperbin < map.size(); isuperbin++) { + auto const& set = map.at(isuperbin); + + unsigned int length = set.size(); + file << isuperbin << "," << length; + for (unsigned int i : set) { + file << "," << i; + } + file << std::endl; + } + } + } + } + +} // namespace lstgeometry + +#endif diff --git a/RecoTracker/LSTGeometry/interface/Module.h b/RecoTracker/LSTGeometry/interface/Module.h new file mode 100644 index 0000000000000..f28f31afe2a8e --- /dev/null +++ b/RecoTracker/LSTGeometry/interface/Module.h @@ -0,0 +1,24 @@ +#ifndef RecoTracker_LSTGeometry_interface_Module_h +#define RecoTracker_LSTGeometry_interface_Module_h + +#include + +#include "RecoTracker/LSTGeometry/interface/Common.h" + +namespace lstgeometry { + + // A module contains 2 sensors. The common properties of the 2 sensors are stored in the Module struct, and the sensor-specific properties are stored in the Sensor struct. + struct Module { + ModuleType moduleType; + SubDetector subdet; + Location location; + Side side; + unsigned int moduleId; + unsigned int layer; + unsigned int ring; + }; + + using Modules = std::unordered_map; +} // namespace lstgeometry + +#endif diff --git a/RecoTracker/LSTGeometry/interface/ModuleMap.h b/RecoTracker/LSTGeometry/interface/ModuleMap.h new file mode 100644 index 0000000000000..53ba71fc6f100 --- /dev/null +++ b/RecoTracker/LSTGeometry/interface/ModuleMap.h @@ -0,0 +1,18 @@ +#ifndef RecoTracker_LSTGeometry_interface_ModuleMapMethods_h +#define RecoTracker_LSTGeometry_interface_ModuleMapMethods_h + +#include +#include + +#include "RecoTracker/LSTGeometry/interface/DetectorGeometry.h" +#include "RecoTracker/LSTGeometry/interface/Sensor.h" + +namespace lstgeometry { + + using ModuleMap = std::unordered_map>; + + ModuleMap buildModuleMap(Sensors const& sensors, DetectorGeometry const& det_geom, 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..a8ced0c0a8a67 --- /dev/null +++ b/RecoTracker/LSTGeometry/interface/PixelMap.h @@ -0,0 +1,24 @@ +#ifndef RecoTracker_LSTGeometry_interface_PixelMap_h +#define RecoTracker_LSTGeometry_interface_PixelMap_h + +#include +#include +#include +#include + +#include "RecoTracker/LSTGeometry/interface/DetectorGeometry.h" +#include "RecoTracker/LSTGeometry/interface/Sensor.h" + +namespace lstgeometry { + + using LayerSubdetChargeKey = std::tuple; + using LayerSubdetChargeMap = std::unordered_map>, + boost::hash>; + using PixelMap = LayerSubdetChargeMap; + + PixelMap buildPixelMap(Sensors const& sensors, DetectorGeometry const& det_geom, 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..15acf01946b49 --- /dev/null +++ b/RecoTracker/LSTGeometry/interface/Sensor.h @@ -0,0 +1,55 @@ +#ifndef RecoTracker_LSTGeometry_interface_Sensor_h +#define RecoTracker_LSTGeometry_interface_Sensor_h + +#include + +#include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h" +#include "DataFormats/GeometrySurface/interface/Surface.h" + +#include "RecoTracker/LSTGeometry/interface/Common.h" + +namespace lstgeometry { + + struct Sensor { + // Module-level properties + ModuleType moduleType; + SubDetector subdet; + Location location; + Side side; + unsigned int moduleId; + unsigned int layer; + unsigned int ring; + bool inverted; + // Sensor-level properties + float centerRho; + float centerZ; + float centerPhi; + bool lower; + bool strip; +#ifndef LST_STANDALONE + MatrixF4x3 corners; +#endif + // Redundant, but convenient to avoid repeated computations + float centerX; + float centerY; + + 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/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..6eaca7a7e2d77 --- /dev/null +++ b/RecoTracker/LSTGeometry/plugins/BuildFile.xml @@ -0,0 +1,15 @@ + + + + + + + + + + + + + + + diff --git a/RecoTracker/LSTGeometry/plugins/LSTGeometryESProducer.cc b/RecoTracker/LSTGeometry/plugins/LSTGeometryESProducer.cc new file mode 100644 index 0000000000000..d0b3f06564971 --- /dev/null +++ b/RecoTracker/LSTGeometry/plugins/LSTGeometryESProducer.cc @@ -0,0 +1,124 @@ +#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" + +#include +#include +#include +#include +#include + +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::ostringstream ptCutOSS; + ptCutOSS << std::setprecision(1) << ptCut_; + std::string ptCutStr = ptCutOSS.str(); + + auto cc = setWhatProduced(this, "LSTGeometry_" + 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_); + + auto sensors = std::make_shared(); + + 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 (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) { + avg_r_barrel[i] /= avg_r_barrel_counter[i]; + } + for (size_t i = 0; i < avg_z_endcap.size(); ++i) { + avg_z_endcap[i] /= avg_z_endcap_counter[i]; + } + + auto lstGeometry = std::make_unique(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/DetectorGeometry.cc b/RecoTracker/LSTGeometry/src/DetectorGeometry.cc new file mode 100644 index 0000000000000..322418f71ce28 --- /dev/null +++ b/RecoTracker/LSTGeometry/src/DetectorGeometry.cc @@ -0,0 +1,260 @@ +#include "RecoTracker/LSTGeometry/interface/DetectorGeometry.h" + +namespace lstgeometry { + + bool DetectorGeometry::isInEtaPhiBin(float eta, float phi, unsigned int eta_bin, unsigned int phi_bin) { + float theta = 2. * std::atan(std::exp(-eta)); + + if (eta_bin == 0) { + if (theta > 3. * kEtaBinRad / 2.) + return false; + } else if (eta_bin == kNEtaBins - 1) { + if (theta < (2 * (kNEtaBins - 1) - 1) * kEtaBinRad / 2.) + return false; + } else if (theta < (2 * eta_bin - 1) * kEtaBinRad / 2. || theta > (2 * (eta_bin + 1) + 1) * kEtaBinRad / 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 DetectorGeometry::getEtaPhiBins(float eta, float phi) { + float theta = 2. * std::atan(std::exp(-eta)); + unsigned int eta_bin = std::clamp(static_cast(theta / kEtaBinRad), 0u, kNEtaBins - 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(eta_bin, phi_bin); + } + + DetectorGeometry::DetectorGeometry(std::shared_ptr sensors, + std::array const& avg_radii, + std::array const& avg_z) + : sensors_(sensors), avg_radii_(avg_radii), avg_z_(avg_z) {} + + MatrixF4x3 const& DetectorGeometry::getCorners(unsigned int detId) const { return sensors_->at(detId).corners; } + + std::vector DetectorGeometry::getDetIds( + std::function&)> filter) const { + std::vector detIds; + detIds.reserve(sensors_->size()); + for (auto const& entry : *sensors_) { + if (filter(entry)) { + detIds.push_back(entry.first); + } + } + return detIds; + } + + void DetectorGeometry::buildByLayer(Sensors const& sensors) { + // Clear just in case they were already built + barrel_lower_det_ids_.clear(); + endcap_lower_det_ids_.clear(); + + // Initialize all vectors + for (unsigned int etabin = 0; etabin < kNEtaBins; etabin++) { + for (unsigned int phibin = 0; phibin < kNPhiBins; phibin++) { + for (unsigned int layer = 1; layer <= kBarrelLayers; layer++) { + barrel_lower_det_ids_[{layer, etabin, phibin}] = {}; + } + for (unsigned int layer = 1; layer <= kEndcapLayers; layer++) { + endcap_lower_det_ids_[{layer, etabin, phibin}] = {}; + } + } + } + + for (unsigned int layer = 1; layer <= kBarrelLayers; layer++) { + auto detids = getDetIds([&sensors, &layer](auto const& x) { + auto const& s = sensors.at(x.first); + return s.location == Location::barrel && s.layer == layer && s.lower; + }); + for (auto detid : detids) { + auto corners = getCorners(detid); + RowVectorF3 center = corners.colwise().mean(); + center /= 4.; + auto etaphi = getEtaPhi(center(1), center(2), center(0)); + for (unsigned int etabin = 0; etabin < kNEtaBins; etabin++) { + for (unsigned int phibin = 0; phibin < kNPhiBins; phibin++) { + if (isInEtaPhiBin(etaphi.first, etaphi.second, etabin, phibin)) { + barrel_lower_det_ids_[{layer, etabin, phibin}].push_back(detid); + } + } + } + } + } + for (unsigned int layer = 1; layer <= kEndcapLayers; layer++) { + auto detids = getDetIds([&sensors, &layer](auto const& x) { + auto const& s = sensors.at(x.first); + return s.location == Location::endcap && s.layer == layer && s.lower; + }); + for (auto detid : detids) { + auto corners = getCorners(detid); + RowVectorF3 center = corners.colwise().mean(); + center /= 4.; + auto etaphi = getEtaPhi(center(1), center(2), center(0)); + for (unsigned int etabin = 0; etabin < kNEtaBins; etabin++) { + for (unsigned int phibin = 0; phibin < kNPhiBins; phibin++) { + if (isInEtaPhiBin(etaphi.first, etaphi.second, etabin, phibin)) { + endcap_lower_det_ids_[{layer, etabin, phibin}].push_back(detid); + } + } + } + } + } + } + + std::vector const& DetectorGeometry::getBarrelLayerDetIds(unsigned int layer, + unsigned int etabin, + unsigned int phibin) const { + return barrel_lower_det_ids_.at({layer, etabin, phibin}); + } + + std::vector const& DetectorGeometry::getEndcapLayerDetIds(unsigned int layer, + unsigned int etabin, + unsigned int phibin) const { + return endcap_lower_det_ids_.at({layer, etabin, phibin}); + } + + float DetectorGeometry::getBarrelLayerAverageRadius(unsigned int layer) const { return avg_radii_[layer - 1]; } + + float DetectorGeometry::getEndcapLayerAverageAbsZ(unsigned int layer) const { return avg_z_[layer - 1]; } + + float DetectorGeometry::getMinR(unsigned int detId) const { + auto const& corners = getCorners(detId); + float minR = std::numeric_limits::max(); + for (int i = 0; i < corners.rows(); i++) { + float x = corners(i, 1); + float y = corners(i, 2); + minR = std::min(minR, std::sqrt(x * x + y * y)); + } + return minR; + } + + float DetectorGeometry::getMaxR(unsigned int detId) const { + auto const& corners = getCorners(detId); + float maxR = std::numeric_limits::min(); + for (int i = 0; i < corners.rows(); i++) { + float x = corners(i, 1); + float y = corners(i, 2); + maxR = std::max(maxR, std::sqrt(x * x + y * y)); + } + return maxR; + } + + float DetectorGeometry::getMinZ(unsigned int detId) const { + auto const& corners = getCorners(detId); + float minZ = std::numeric_limits::max(); + for (int i = 0; i < corners.rows(); i++) { + float z = corners(i, 0); + minZ = std::min(minZ, z); + } + return minZ; + } + + float DetectorGeometry::getMaxZ(unsigned int detId) const { + auto const& corners = getCorners(detId); + float maxZ = std::numeric_limits::lowest(); + for (int i = 0; i < corners.rows(); i++) { + float z = corners(i, 0); + maxZ = std::max(maxZ, z); + } + return maxZ; + } + + float DetectorGeometry::getMinPhi(unsigned int detId) const { + auto const& corners = getCorners(detId); + float minPhi = std::numeric_limits::max(); + float minPosPhi = std::numeric_limits::max(); + float minNegPhi = std::numeric_limits::max(); + unsigned int nPos = 0; + unsigned int nOverPi2 = 0; + for (int i = 0; i < corners.rows(); i++) { + float phi = phi_mpi_pi(std::numbers::pi_v + std::atan2(-corners(i, 2), -corners(i, 1))); + minPhi = std::min(minPhi, phi); + if (phi > 0) { + minPosPhi = std::min(minPosPhi, phi); + nPos++; + } else { + minNegPhi = std::min(minNegPhi, phi); + } + if (std::fabs(phi) > std::numbers::pi_v / 2.) { + nOverPi2++; + } + } + if (nPos == 4 || nPos == 0) + return minPhi; + if (nOverPi2 == 4) + return minPosPhi; + return minPhi; + } + + float DetectorGeometry::getMaxPhi(unsigned int detId) const { + auto const& corners = getCorners(detId); + float maxPhi = std::numeric_limits::lowest(); + float maxPosPhi = std::numeric_limits::lowest(); + float maxNegPhi = std::numeric_limits::lowest(); + unsigned int nPos = 0; + unsigned int nOverPi2 = 0; + for (int i = 0; i < corners.rows(); i++) { + float phi = phi_mpi_pi(std::numbers::pi_v + std::atan2(-corners(i, 2), -corners(i, 1))); + maxPhi = std::max(maxPhi, phi); + if (phi > 0) { + maxPosPhi = std::max(maxPosPhi, phi); + nPos++; + } else { + maxNegPhi = std::max(maxNegPhi, phi); + } + if (std::fabs(phi) > std::numbers::pi_v / 2.) { + nOverPi2++; + } + } + if (nPos == 4 || nPos == 0) + return maxPhi; + if (nOverPi2 == 4) + return maxNegPhi; + return maxPhi; + } + + std::pair DetectorGeometry::getCompatibleEtaRange(unsigned int detId, + float zmin_bound, + float zmax_bound) const { + float minr = getMinR(detId); + float maxr = getMaxR(detId); + float minz = getMinZ(detId); + float maxz = getMaxZ(detId); + 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> DetectorGeometry::getCompatiblePhiRange( + unsigned int detId, float ptmin, float ptmax) const { + float minr = getMinR(detId); + float maxr = getMaxR(detId); + float minphi = getMinPhi(detId); + float maxphi = getMaxPhi(detId); + 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)); + } +} // namespace lstgeometry diff --git a/RecoTracker/LSTGeometry/src/Geometry.cc b/RecoTracker/LSTGeometry/src/Geometry.cc new file mode 100644 index 0000000000000..631760a079c4f --- /dev/null +++ b/RecoTracker/LSTGeometry/src/Geometry.cc @@ -0,0 +1,25 @@ +#include "RecoTracker/LSTGeometry/interface/Geometry.h" +#include "RecoTracker/LSTGeometry/interface/DetectorGeometry.h" +#include "RecoTracker/LSTGeometry/interface/Slope.h" + +using namespace lstgeometry; + +Geometry::Geometry(std::shared_ptr sensors, + std::array const &average_r_barrel, + std::array const &average_z_endcap, + float pt_cut) + : sensors(sensors) { + auto slopes = computeSlopes(*sensors); + barrel_slopes = std::move(std::get<0>(slopes)); + endcap_slopes = std::move(std::get<1>(slopes)); + + auto det_geom = DetectorGeometry(sensors, average_r_barrel, average_z_endcap); + det_geom.buildByLayer(*sensors); + + pixel_map = buildPixelMap(*sensors, det_geom, pt_cut); + + module_map = buildModuleMap(*sensors, det_geom, pt_cut); +} + +#include "FWCore/Utilities/interface/typelookup.h" +TYPELOOKUP_DATA_REG(Geometry); diff --git a/RecoTracker/LSTGeometry/src/ModuleMap.cc b/RecoTracker/LSTGeometry/src/ModuleMap.cc new file mode 100644 index 0000000000000..b58f8a2875c22 --- /dev/null +++ b/RecoTracker/LSTGeometry/src/ModuleMap.cc @@ -0,0 +1,378 @@ +#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(MatrixFNx3 const& mod_boundaries, float refphi, float zshift = 0) { + int npoints = mod_boundaries.rows(); + MatrixFNx2 mod_boundaries_etaphi(npoints, 2); + for (int i = 0; i < npoints; ++i) { + auto ref_etaphi = getEtaPhi(mod_boundaries(i, 1), mod_boundaries(i, 2), mod_boundaries(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 <= npoints; ++i) { + boost::geometry::append(poly, + Point(mod_boundaries_etaphi(i % npoints, 0), mod_boundaries_etaphi(i % npoints, 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) { + RowVectorF3 ref_center = ref_mod_boundaries.colwise().sum(); + ref_center /= 4.; + RowVectorF3 tar_center = tar_mod_boundaries.colwise().sum(); + tar_center /= 4.; + + float ref_center_phi = std::atan2(ref_center(2), ref_center(1)); + float 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, + DetectorGeometry const& det_geom) { + auto const& sensor = sensors.at(ref_detid); + + float refphi = std::atan2(sensor.centerY, sensor.centerX); + + auto refmodule = sensors.at(ref_detid); + + unsigned short ref_layer = refmodule.layer; + auto ref_location = refmodule.location; + + auto etaphi = getEtaPhi(sensor.centerX, sensor.centerY, sensor.centerZ); + auto etaphibins = DetectorGeometry::getEtaPhiBins(etaphi.first, etaphi.second); + + auto const& tar_detids_to_be_considered = + ref_location == Location::barrel + ? det_geom.getBarrelLayerDetIds(ref_layer + 1, etaphibins.first, etaphibins.second) + : det_geom.getEndcapLayerDetIds(ref_layer + 1, etaphibins.first, etaphibins.second); + + std::vector list_of_detids_etaphi_layer_tar; + for (unsigned int tar_detid : tar_detids_to_be_considered) { + if (moduleOverlapsInEtaPhi(det_geom.getCorners(ref_detid), det_geom.getCorners(tar_detid), refphi, 0) || + moduleOverlapsInEtaPhi(det_geom.getCorners(ref_detid), det_geom.getCorners(tar_detid), refphi, 10) || + moduleOverlapsInEtaPhi(det_geom.getCorners(ref_detid), det_geom.getCorners(tar_detid), refphi, -10)) + 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::unordered_set barrel_endcap_connected_tar_detids; + + for (float zshift : {0, 10, -10}) { + std::vector ref_polygon; + ref_polygon.push_back(getEtaPhiPolygon(det_geom.getCorners(ref_detid), 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.size()) + break; + Polygon tar_polygon = getEtaPhiPolygon(det_geom.getCorners(tar_detid), 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 <= 1e-6) + continue; + + auto const& new_tar_detids_to_be_considered = + det_geom.getEndcapLayerDetIds(1, etaphibins.first, etaphibins.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(det_geom.getCorners(tar_detid), 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.insert(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, DetectorGeometry const& det_geom, float ptCut, bool doR = true) { + auto bounds = det_geom.getCorners(ref_detid); + auto const& sensor = sensors.at(ref_detid); + int charge = 1; + float z_r = sensor.centerZ / std::sqrt(sensor.centerX * sensor.centerX + sensor.centerY * sensor.centerY); + float refphi = std::atan2(sensor.centerY, sensor.centerX); + auto const& refmodule = sensors.at(ref_detid); + unsigned short ref_layer = refmodule.layer; + auto ref_location = refmodule.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 = det_geom.getBarrelLayerAverageRadius(ref_layer + 1); + 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 = det_geom.getEndcapLayerAverageAbsZ(1); + 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 = det_geom.getEndcapLayerAverageAbsZ(ref_layer + 1); + 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, + DetectorGeometry const& det_geom, + float ptCut) { + auto const& sensor = sensors.at(ref_detid); + + float refphi = std::atan2(sensor.centerY, sensor.centerX); + + auto const& refmodule = sensors.at(ref_detid); + + unsigned short ref_layer = refmodule.layer; + auto ref_location = refmodule.location; + + auto etaphi = getEtaPhi(sensor.centerX, sensor.centerY, sensor.centerZ); + auto etaphibins = DetectorGeometry::getEtaPhiBins(etaphi.first, etaphi.second); + + auto const& tar_detids_to_be_considered = + ref_location == Location::barrel + ? det_geom.getBarrelLayerDetIds(ref_layer + 1, etaphibins.first, etaphibins.second) + : det_geom.getEndcapLayerDetIds(ref_layer + 1, etaphibins.first, etaphibins.second); + + auto next_layer_bound_points = boundsAfterCurved(ref_detid, sensors, det_geom, ptCut); + + std::vector list_of_detids_etaphi_layer_tar; + for (unsigned int tar_detid : tar_detids_to_be_considered) { + if (moduleOverlapsInEtaPhi(next_layer_bound_points, det_geom.getCorners(tar_detid), refphi, 0)) + 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::unordered_set 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.size()) + break; + Polygon tar_polygon = getEtaPhiPolygon(det_geom.getCorners(tar_detid), 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 > 1e-6) { + auto const& new_tar_detids_to_be_considered = + det_geom.getEndcapLayerDetIds(1, etaphibins.first, etaphibins.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(det_geom.getCorners(tar_detid), 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.insert(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(list.begin(), list.end()); + } + } + + return merged; + } + + ModuleMap buildModuleMap(Sensors const& sensors, DetectorGeometry const& det_geo, float pt_cut) { + auto detids_etaphi_layer_ref = det_geo.getDetIds([&sensors](auto const& x) { + auto const& s = sensors.at(x.first); + // exclude the outermost modules that do not have connections to other modules + return ((s.location == Location::barrel && s.lower && s.layer != 6) || + (s.location == Location::endcap && s.lower && s.layer != 5 && !(s.ring == 15 && s.layer == 1) && + !(s.ring == 15 && s.layer == 2) && !(s.ring == 12 && s.layer == 3) && !(s.ring == 12 && s.layer == 4))); + }); + + std::unordered_map> straight_line_connections; + std::unordered_map> curved_line_connections; + + for (auto ref_detid : detids_etaphi_layer_ref) { + straight_line_connections[ref_detid] = getStraightLineConnections(ref_detid, sensors, det_geo); + curved_line_connections[ref_detid] = getCurvedLineConnections(ref_detid, sensors, det_geo, 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..20b191703bc4b --- /dev/null +++ b/RecoTracker/LSTGeometry/src/PixelMap.cc @@ -0,0 +1,127 @@ +#include "RecoTracker/LSTGeometry/interface/Common.h" +#include "RecoTracker/LSTGeometry/interface/PixelMap.h" + +namespace lstgeometry { + + PixelMap buildPixelMap(Sensors const& sensors, DetectorGeometry const& det_geom, 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 detId : det_geom.getDetIds()) { + auto const& sensor = sensors.at(detId); + + // Phase-2 enum differs from the legacy one used here + unsigned int subdet = sensor.subdet == SubDetector::P2OTB ? SubDet::Barrel : SubDet::Endcap; + auto layer = sensor.layer; + if (layer > 2) + continue; + auto location = sensor.location; + + // Skip if the module is not PS module and is not lower sensor + if (sensor.moduleType == ModuleType::Ph2SS || !sensor.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] = det_geom.getCompatibleEtaRange(detId, 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 = det_geom.getCompatiblePhiRange(detId, 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].insert(detId); + maps[{layer, subdet, 0}][isuperbin].insert(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].insert(detId); + maps[{layer, subdet, 0}][isuperbin].insert(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].insert(detId); + maps[{layer, subdet, 0}][isuperbin].insert(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].insert(detId); + maps[{layer, subdet, 0}][isuperbin].insert(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].insert(detId); + maps[{layer, subdet, 0}][isuperbin].insert(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].insert(detId); + maps[{layer, subdet, 0}][isuperbin].insert(detId); + } + } + } + } + } + } + + return maps; + } + +} // namespace lstgeometry diff --git a/RecoTracker/LSTGeometry/src/Sensor.cc b/RecoTracker/LSTGeometry/src/Sensor.cc new file mode 100644 index 0000000000000..8ad711b9bf8b0 --- /dev/null +++ b/RecoTracker/LSTGeometry/src/Sensor.cc @@ -0,0 +1,90 @@ +#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), + subdet(subdet), + location(location), + side(side), + moduleId(moduleId), + layer(layer), + ring(ring), + inverted(isInverted(moduleId, location, side, layer)), + centerRho(centerRho), + centerZ(centerZ), + centerPhi(centerPhi), + lower(isLower(moduleId, location, side, layer, detId)), + strip(isStrip(moduleType, inverted, lower)), + centerX(centerRho * std::cos(centerPhi)), + centerY(centerRho * std::sin(centerPhi)) { + const Bounds &bounds = surface.bounds(); + const RectangularPlaneBounds &plane_bounds = dynamic_cast(bounds); + float wid = plane_bounds.width(); + float len = plane_bounds.length(); + 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 + 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(); +} diff --git a/RecoTracker/LSTGeometry/src/Slope.cc b/RecoTracker/LSTGeometry/src/Slope.cc new file mode 100644 index 0000000000000..eab4161cbec38 --- /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.corners(1, 1) - sensor.corners(0, 1)); + float dy = roundCoordinate(sensor.corners(1, 2) - sensor.corners(0, 2)); + float dz = roundCoordinate(sensor.corners(1, 0) - sensor.corners(0, 0)); + + Slope slope(dx, dy, dz); + + auto location = sensor.location; + bool is_tilted = sensor.side != Side::Center; + bool is_strip = sensor.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..ed9690e1c1817 --- /dev/null +++ b/RecoTracker/LSTGeometry/test/DumpLSTGeometry.cc @@ -0,0 +1,54 @@ +#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/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_; + + static std::string getPtCutStr(double ptCut) { + std::ostringstream ptCutOSS; + ptCutOSS << std::setprecision(1) << ptCut; + return ptCutOSS.str(); + } +}; + +DumpLSTGeometry::DumpLSTGeometry(const edm::ParameterSet& config) + : ptCut_(config.getParameter("ptCut")), + ptCutStr_(getPtCutStr(ptCut_)), + outputDirectory_(config.getUntrackedParameter("outputDirectory", "data/")), + binaryOutput_(config.getUntrackedParameter("outputAsBinary", true)), + lstGeoToken_{esConsumes(edm::ESInputTag("", "LSTGeometry_" + 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..7082a094333a3 --- /dev/null +++ b/RecoTracker/LSTGeometry/test/dumpLSTGeometry.py @@ -0,0 +1,68 @@ +import FWCore.ParameterSet.Config as cms +from Configuration.AlCa.GlobalTag import GlobalTag +import Configuration.Geometry.defaultPhase2ConditionsEra_cff as _settings +import FWCore.ParameterSet.VarParsing as VarParsing + +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 + +process.add_(cms.ESProducer("LSTGeometryESProducer")) + +options = VarParsing.VarParsing() +options.register( + "outputDirectory", + "data/", + VarParsing.VarParsing.multiplicity.singleton, + VarParsing.VarParsing.varType.string, + "Output directory for LST geometry files", +) +options.register( + "ptCut", + 0.8, + VarParsing.VarParsing.multiplicity.singleton, + VarParsing.VarParsing.varType.float, + "pT cut for LST module maps", +) +options.register( + "binaryOutput", + True, + VarParsing.VarParsing.multiplicity.singleton, + VarParsing.VarParsing.varType.bool, + "Dump LST geometry as binary files", +) +options.parseArguments() + + +process.dump = cms.EDAnalyzer( + "DumpLSTGeometry", + outputDirectory = cms.untracked.string(options.outputDirectory), + ptCut = cms.double(options.ptCut), + outputAsBinary = cms.untracked.bool(options.binaryOutput), +) + +print(f"Requesting LST geometry dump into directory: {options.outputDirectory}\n") +process.p = cms.Path(process.dump) diff --git a/RecoTracker/LSTGeometry/test/testDumpLSTGeometry.sh b/RecoTracker/LSTGeometry/test/testDumpLSTGeometry.sh new file mode 100755 index 0000000000000..f3a7524015b53 --- /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=true) || die "failed to run dumpLSTGeometry.py" $?