CudaMuonSegmentFitting/src/CudaMuonSpacePoint.cpp

0001 // This file is part of the ACTS project.
0002 //
0003 // Copyright (C) 2016 CERN for the benefit of the ACTS project
0004 //
0005 // This Source Code Form is subject to the terms of the Mozilla Public
0006 // License, v. 2.0. If a copy of the MPL was not distributed with this
0007 // file, You can obtain one at https://mozilla.org/MPL/2.0/.
0008
0009 #include "ActsExamples/EventData/CudaMuonSpacePoint.hpp"
0010
0011 #include <stdexcept>
0012 #include <utility>
0013
0014 #include <cuda_runtime.h>
0015
0016 #include "CudaUtilities.hpp"
0017
0018 namespace {
0019
0020 using namespace ActsExamples;
0021
0022 template <typename T>
0023 void resizeColumn(std::vector<T>& column, std::size_t size) {
0024   column.resize(size);
0025 }
0026
0027 void resizeHostData(ActsExamples::CudaMuonSpacePointHostData& host,
0028                     std::size_t size) {
0029   resizeColumn(host.geometryId, size);
0030   resizeColumn(host.muonId, size);
0031
0032   resizeColumn(host.localPositionX, size);
0033   resizeColumn(host.localPositionY, size);
0034   resizeColumn(host.localPositionZ, size);
0035
0036   resizeColumn(host.sensorDirectionX, size);
0037   resizeColumn(host.sensorDirectionY, size);
0038   resizeColumn(host.sensorDirectionZ, size);
0039
0040   resizeColumn(host.toNextSensorX, size);
0041   resizeColumn(host.toNextSensorY, size);
0042   resizeColumn(host.toNextSensorZ, size);
0043
0044   resizeColumn(host.planeNormalX, size);
0045   resizeColumn(host.planeNormalY, size);
0046   resizeColumn(host.planeNormalZ, size);
0047
0048   resizeColumn(host.covariance0, size);
0049   resizeColumn(host.covariance1, size);
0050   resizeColumn(host.covariance2, size);
0051
0052   resizeColumn(host.driftRadius, size);
0053   resizeColumn(host.time, size);
0054 }
0055
0056 void allocateDeviceData(ActsExamples::CudaMuonSpacePointArrays& device,
0057                         std::size_t spacePointCount, std::size_t bucketCount) {
0058   allocateDeviceColumn(device.geometryId, spacePointCount);
0059   allocateDeviceColumn(device.muonId, spacePointCount);
0060
0061   allocateDeviceColumn(device.localPositionX, spacePointCount);
0062   allocateDeviceColumn(device.localPositionY, spacePointCount);
0063   allocateDeviceColumn(device.localPositionZ, spacePointCount);
0064
0065   allocateDeviceColumn(device.sensorDirectionX, spacePointCount);
0066   allocateDeviceColumn(device.sensorDirectionY, spacePointCount);
0067   allocateDeviceColumn(device.sensorDirectionZ, spacePointCount);
0068
0069   allocateDeviceColumn(device.toNextSensorX, spacePointCount);
0070   allocateDeviceColumn(device.toNextSensorY, spacePointCount);
0071   allocateDeviceColumn(device.toNextSensorZ, spacePointCount);
0072
0073   allocateDeviceColumn(device.planeNormalX, spacePointCount);
0074   allocateDeviceColumn(device.planeNormalY, spacePointCount);
0075   allocateDeviceColumn(device.planeNormalZ, spacePointCount);
0076
0077   allocateDeviceColumn(device.covariance0, spacePointCount);
0078   allocateDeviceColumn(device.covariance1, spacePointCount);
0079   allocateDeviceColumn(device.covariance2, spacePointCount);
0080
0081   allocateDeviceColumn(device.driftRadius, spacePointCount);
0082   allocateDeviceColumn(device.time, spacePointCount);
0083
0084   allocateDeviceColumn(device.bucketStart, bucketCount);
0085   allocateDeviceColumn(device.bucketEnd, bucketCount);
0086 }
0087
0088 void freeDeviceData(ActsExamples::CudaMuonSpacePointArrays& device) noexcept {
0089   freeDeviceColumn(device.geometryId);
0090   freeDeviceColumn(device.muonId);
0091
0092   freeDeviceColumn(device.localPositionX);
0093   freeDeviceColumn(device.localPositionY);
0094   freeDeviceColumn(device.localPositionZ);
0095
0096   freeDeviceColumn(device.sensorDirectionX);
0097   freeDeviceColumn(device.sensorDirectionY);
0098   freeDeviceColumn(device.sensorDirectionZ);
0099
0100   freeDeviceColumn(device.toNextSensorX);
0101   freeDeviceColumn(device.toNextSensorY);
0102   freeDeviceColumn(device.toNextSensorZ);
0103
0104   freeDeviceColumn(device.planeNormalX);
0105   freeDeviceColumn(device.planeNormalY);
0106   freeDeviceColumn(device.planeNormalZ);
0107
0108   freeDeviceColumn(device.covariance0);
0109   freeDeviceColumn(device.covariance1);
0110   freeDeviceColumn(device.covariance2);
0111
0112   freeDeviceColumn(device.driftRadius);
0113   freeDeviceColumn(device.time);
0114
0115   freeDeviceColumn(device.bucketStart);
0116   freeDeviceColumn(device.bucketEnd);
0117 }
0118
0119 void copyHostToDevice(ActsExamples::CudaMuonSpacePointArrays& device,
0120                       const ActsExamples::CudaMuonSpacePointHostData& host) {
0121   copyColumnToDevice(device.geometryId, host.geometryId);
0122   copyColumnToDevice(device.muonId, host.muonId);
0123
0124   copyColumnToDevice(device.localPositionX, host.localPositionX);
0125   copyColumnToDevice(device.localPositionY, host.localPositionY);
0126   copyColumnToDevice(device.localPositionZ, host.localPositionZ);
0127
0128   copyColumnToDevice(device.sensorDirectionX, host.sensorDirectionX);
0129   copyColumnToDevice(device.sensorDirectionY, host.sensorDirectionY);
0130   copyColumnToDevice(device.sensorDirectionZ, host.sensorDirectionZ);
0131
0132   copyColumnToDevice(device.toNextSensorX, host.toNextSensorX);
0133   copyColumnToDevice(device.toNextSensorY, host.toNextSensorY);
0134   copyColumnToDevice(device.toNextSensorZ, host.toNextSensorZ);
0135
0136   copyColumnToDevice(device.planeNormalX, host.planeNormalX);
0137   copyColumnToDevice(device.planeNormalY, host.planeNormalY);
0138   copyColumnToDevice(device.planeNormalZ, host.planeNormalZ);
0139
0140   copyColumnToDevice(device.covariance0, host.covariance0);
0141   copyColumnToDevice(device.covariance1, host.covariance1);
0142   copyColumnToDevice(device.covariance2, host.covariance2);
0143
0144   copyColumnToDevice(device.driftRadius, host.driftRadius);
0145   copyColumnToDevice(device.time, host.time);
0146
0147   copyColumnToDevice(device.bucketStart, host.bucketStart);
0148   copyColumnToDevice(device.bucketEnd, host.bucketEnd);
0149 }
0150
0151 void copyDeviceToHost(ActsExamples::CudaMuonSpacePointHostData& host,
0152                       const ActsExamples::CudaMuonSpacePointArrays& device) {
0153   copyColumnToHost(host.geometryId, device.geometryId);
0154   copyColumnToHost(host.muonId, device.muonId);
0155
0156   copyColumnToHost(host.localPositionX, device.localPositionX);
0157   copyColumnToHost(host.localPositionY, device.localPositionY);
0158   copyColumnToHost(host.localPositionZ, device.localPositionZ);
0159
0160   copyColumnToHost(host.sensorDirectionX, device.sensorDirectionX);
0161   copyColumnToHost(host.sensorDirectionY, device.sensorDirectionY);
0162   copyColumnToHost(host.sensorDirectionZ, device.sensorDirectionZ);
0163
0164   copyColumnToHost(host.toNextSensorX, device.toNextSensorX);
0165   copyColumnToHost(host.toNextSensorY, device.toNextSensorY);
0166   copyColumnToHost(host.toNextSensorZ, device.toNextSensorZ);
0167
0168   copyColumnToHost(host.planeNormalX, device.planeNormalX);
0169   copyColumnToHost(host.planeNormalY, device.planeNormalY);
0170   copyColumnToHost(host.planeNormalZ, device.planeNormalZ);
0171
0172   copyColumnToHost(host.covariance0, device.covariance0);
0173   copyColumnToHost(host.covariance1, device.covariance1);
0174   copyColumnToHost(host.covariance2, device.covariance2);
0175
0176   copyColumnToHost(host.driftRadius, device.driftRadius);
0177   copyColumnToHost(host.time, device.time);
0178
0179   copyColumnToHost(host.bucketStart, device.bucketStart);
0180   copyColumnToHost(host.bucketEnd, device.bucketEnd);
0181 }
0182
0183 std::size_t countSpacePoints(
0184     const ActsExamples::MuonSpacePointContainer& spacePoints) {
0185   std::size_t totalSize = 0;
0186
0187   for (const auto& bucket : spacePoints) {
0188     totalSize += bucket.size();
0189   }
0190
0191   return totalSize;
0192 }
0193
0194 }  // namespace
0195
0196 namespace ActsExamples {
0197
0198 CudaMuonSpacePointProxy::CudaMuonSpacePointProxy(
0199     const CudaMuonSpacePointContainer& container, std::size_t index) noexcept
0200     : m_container{&container}, m_index{index} {}
0201
0202 const MuonSpacePoint::MuonId& CudaMuonSpacePointProxy::id() const {
0203   m_idCache = MuonSpacePoint::MuonId{m_container->m_host.muonId[m_index]};
0204   return m_idCache;
0205 }
0206
0207 const Acts::GeometryIdentifier& CudaMuonSpacePointProxy::geometryId() const {
0208   m_geometryIdCache =
0209       Acts::GeometryIdentifier{m_container->m_host.geometryId[m_index]};
0210   return m_geometryIdCache;
0211 }
0212
0213 const Acts::Vector3& CudaMuonSpacePointProxy::localPosition() const {
0214   m_localPositionCache = Acts::Vector3{
0215       m_container->m_host.localPositionX[m_index],
0216       m_container->m_host.localPositionY[m_index],
0217       m_container->m_host.localPositionZ[m_index],
0218   };
0219
0220   return m_localPositionCache;
0221 }
0222
0223 const Acts::Vector3& CudaMuonSpacePointProxy::sensorDirection() const {
0224   m_sensorDirectionCache = Acts::Vector3{
0225       m_container->m_host.sensorDirectionX[m_index],
0226       m_container->m_host.sensorDirectionY[m_index],
0227       m_container->m_host.sensorDirectionZ[m_index],
0228   };
0229
0230   return m_sensorDirectionCache;
0231 }
0232
0233 const Acts::Vector3& CudaMuonSpacePointProxy::planeNormal() const {
0234   m_planeNormalCache = Acts::Vector3{
0235       m_container->m_host.planeNormalX[m_index],
0236       m_container->m_host.planeNormalY[m_index],
0237       m_container->m_host.planeNormalZ[m_index],
0238   };
0239
0240   return m_planeNormalCache;
0241 }
0242
0243 const Acts::Vector3& CudaMuonSpacePointProxy::toNextSensor() const {
0244   m_toNextSensorCache = Acts::Vector3{
0245       m_container->m_host.toNextSensorX[m_index],
0246       m_container->m_host.toNextSensorY[m_index],
0247       m_container->m_host.toNextSensorZ[m_index],
0248   };
0249
0250   return m_toNextSensorCache;
0251 }
0252
0253 const std::array<double, 3>& CudaMuonSpacePointProxy::covariance() const {
0254   m_covarianceCache = {
0255       m_container->m_host.covariance0[m_index],
0256       m_container->m_host.covariance1[m_index],
0257       m_container->m_host.covariance2[m_index],
0258   };
0259
0260   return m_covarianceCache;
0261 }
0262
0263 double CudaMuonSpacePointProxy::driftRadius() const {
0264   return m_container->m_host.driftRadius[m_index];
0265 }
0266
0267 double CudaMuonSpacePointProxy::time() const {
0268   return m_container->m_host.time[m_index];
0269 }
0270
0271 bool CudaMuonSpacePointProxy::isStraw() const {
0272   return id().technology() == MuonSpacePoint::MuonId::TechField::Mdt;
0273 }
0274
0275 bool CudaMuonSpacePointProxy::hasTime() const {
0276   return id().measuresTime();
0277 }
0278
0279 bool CudaMuonSpacePointProxy::measuresLoc0() const {
0280   return id().measuresPhi();
0281 }
0282
0283 bool CudaMuonSpacePointProxy::measuresLoc1() const {
0284   return id().measuresEta();
0285 }
0286
0287 CudaMuonSpacePointPtr::CudaMuonSpacePointPtr(
0288     const CudaMuonSpacePointContainer& container, std::size_t index) noexcept
0289     : m_proxy{container, index}, m_valid{true} {}
0290
0291 CudaMuonSpacePointContainer::CudaMuonSpacePointContainer(size_type size)
0292     : m_size{size} {
0293   resizeHostData(m_host, m_size);
0294 }
0295
0296 CudaMuonSpacePointContainer::CudaMuonSpacePointContainer(
0297     const MuonSpacePointContainer& spacePoints)
0298     : CudaMuonSpacePointContainer{countSpacePoints(spacePoints)} {
0299   m_host.bucketStart.reserve(spacePoints.size());
0300   m_host.bucketEnd.reserve(spacePoints.size());
0301
0302   size_type index = 0;
0303
0304   for (const MuonSpacePointBucket& bucket : spacePoints) {
0305     const size_type bucketBegin = index;
0306
0307     for (const MuonSpacePoint& spacePoint : bucket) {
0308       setGeometryId(index, spacePoint.geometryId().value());
0309       setId(index, spacePoint.id().toInt());
0310
0311       defineCoordinates(index, spacePoint.localPosition(),
0312                         spacePoint.sensorDirection(),
0313                         spacePoint.toNextSensor());
0314
0315       setRadius(index, spacePoint.driftRadius());
0316       setTime(index, spacePoint.time());
0317
0318       const auto& covariance = spacePoint.covariance();
0319       setCovariance(index, covariance[0], covariance[1], covariance[2]);
0320
0321       ++index;
0322     }
0323
0324     // Empty buckets are preserved as ranges for which begin == end.
0325     addBucket(bucketBegin, index);
0326   }
0327 }
0328
0329 CudaMuonSpacePointContainer::CudaMuonSpacePointContainer(
0330     CudaMuonSpacePointContainer&& other) noexcept
0331     : m_size{std::exchange(other.m_size, 0)},
0332       m_host{std::move(other.m_host)},
0333       m_device{std::exchange(other.m_device, {})},
0334       m_onDevice{std::exchange(other.m_onDevice, false)} {}
0335
0336 CudaMuonSpacePointContainer& CudaMuonSpacePointContainer::operator=(
0337     CudaMuonSpacePointContainer&& other) noexcept {
0338   if (this != &other) {
0339     clearDevice();
0340
0341     m_size = std::exchange(other.m_size, 0);
0342     m_host = std::move(other.m_host);
0343     m_device = std::exchange(other.m_device, {});
0344     m_onDevice = std::exchange(other.m_onDevice, false);
0345   }
0346
0347   return *this;
0348 }
0349
0350 CudaMuonSpacePointContainer::~CudaMuonSpacePointContainer() noexcept {
0351   clearDevice();
0352 }
0353
0354 CudaMuonSpacePointContainer::size_type CudaMuonSpacePointContainer::bucketStart(
0355     size_type bucket) const {
0356   checkBucket(bucket);
0357   return m_host.bucketStart[bucket];
0358 }
0359
0360 CudaMuonSpacePointContainer::size_type CudaMuonSpacePointContainer::bucketEnd(
0361     size_type bucket) const {
0362   checkBucket(bucket);
0363   return m_host.bucketEnd[bucket];
0364 }
0365
0366 void CudaMuonSpacePointContainer::addBucket(size_type start, size_type end) {
0367   if (start > end || end > m_size) {
0368     throw std::out_of_range(
0369         std::format("CudaMuonSpacePointContainer invalid bucket range "
0370                     "[{:};{:}]. Allowed [0;{:})",
0371                     start, end, m_size));
0372   }
0373
0374   m_host.bucketStart.push_back(static_cast<std::uint32_t>(start));
0375   m_host.bucketEnd.push_back(static_cast<std::uint32_t>(end));
0376 }
0377
0378 void CudaMuonSpacePointContainer::setGeometryId(
0379     size_type index, Acts::GeometryIdentifier::Value geometryId) {
0380   checkIndex(index);
0381   m_host.geometryId[index] = geometryId;
0382 }
0383
0384 void CudaMuonSpacePointContainer::setId(size_type index, std::uint32_t muonId) {
0385   checkIndex(index);
0386   m_host.muonId[index] = muonId;
0387 }
0388
0389 void CudaMuonSpacePointContainer::defineCoordinates(
0390     size_type index, const Acts::Vector3& position,
0391     const Acts::Vector3& sensorDirection, const Acts::Vector3& toNextSensor) {
0392   checkIndex(index);
0393
0394   m_host.localPositionX[index] = position.x();
0395   m_host.localPositionY[index] = position.y();
0396   m_host.localPositionZ[index] = position.z();
0397
0398   m_host.sensorDirectionX[index] = sensorDirection.x();
0399   m_host.sensorDirectionY[index] = sensorDirection.y();
0400   m_host.sensorDirectionZ[index] = sensorDirection.z();
0401
0402   m_host.toNextSensorX[index] = toNextSensor.x();
0403   m_host.toNextSensorY[index] = toNextSensor.y();
0404   m_host.toNextSensorZ[index] = toNextSensor.z();
0405
0406   // CompositeSpacePoint expects the measurement plane normal
0407   const Acts::Vector3 planeNormal =
0408       sensorDirection.cross(toNextSensor).normalized();
0409
0410   m_host.planeNormalX[index] = planeNormal.x();
0411   m_host.planeNormalY[index] = planeNormal.y();
0412   m_host.planeNormalZ[index] = planeNormal.z();
0413 }
0414
0415 void CudaMuonSpacePointContainer::setRadius(size_type index, double radius) {
0416   checkIndex(index);
0417   m_host.driftRadius[index] = radius;
0418 }
0419
0420 void CudaMuonSpacePointContainer::setTime(size_type index, double time) {
0421   checkIndex(index);
0422   m_host.time[index] = time;
0423 }
0424
0425 void CudaMuonSpacePointContainer::setCovariance(size_type index, double cov0,
0426                                                 double cov1, double cov2) {
0427   checkIndex(index);
0428
0429   m_host.covariance0[index] = cov0;
0430   m_host.covariance1[index] = cov1;
0431   m_host.covariance2[index] = cov2;
0432 }
0433
0434 void CudaMuonSpacePointContainer::moveToDevice() {
0435   clearDevice();
0436
0437   allocateDeviceData(m_device, m_size, bucketCount());
0438   copyHostToDevice(m_device, m_host);
0439
0440   m_onDevice = true;
0441 }
0442
0443 void CudaMuonSpacePointContainer::moveToHost() {
0444   if (!m_onDevice) {
0445     return;
0446   }
0447
0448   copyDeviceToHost(m_host, m_device);
0449 }
0450
0451 void CudaMuonSpacePointContainer::clearDevice() noexcept {
0452   freeDeviceData(m_device);
0453   m_device = {};
0454   m_onDevice = false;
0455 }
0456
0457 CudaMuonSpacePointContainer::value_type CudaMuonSpacePointContainer::at(
0458     size_type index) const {
0459   checkIndex(index);
0460   return value_type{*this, index};
0461 }
0462
0463 inline void CudaMuonSpacePointContainer::checkIndex(size_type index) const {
0464   if (index >= m_size) {
0465     throw std::out_of_range(std::format(
0466         "CudaMuonSpacePointContainer index {:}out of range. Max allowed: {:}",
0467         index, m_size - 1));
0468   }
0469 }
0470
0471 void CudaMuonSpacePointContainer::checkBucket(size_type bucket) const {
0472   if (bucket >= bucketCount()) {
0473     throw std::out_of_range(std::format(
0474         "CudaMuonSpacePointContainer bucket {:} out of range. Max allowed: {:}",
0475         bucket, bucketCount() - 1));
0476   }
0477 }
0478
0479 // The raw MuonId stores layer - 1 in bits 17..20. Here `layer` is already
0480 // zero-based because CUDA layer masks use layer 0 -> bit 0.
0481 void CudaMuonSpacePointContainer::setLogicalLayer(size_type index,
0482                                                   std::uint32_t layer) {
0483   checkIndex(index);
0484
0485   static constexpr std::uint32_t fourBit = 0xFu;
0486   static constexpr std::uint32_t layerShift = 17u;
0487
0488   if (layer > fourBit) {
0489     throw std::out_of_range(
0490         "CudaMuonSpacePointContainer logical layer out of range");
0491   }
0492
0493   std::uint32_t rawId = m_host.muonId[index];
0494
0495   // Clear old layer bits 17..20.
0496   rawId &= ~(fourBit << layerShift);
0497
0498   // Insert new layer bits 17..20.
0499   rawId |= ((layer & fourBit) << layerShift);
0500
0501   m_host.muonId[index] = rawId;
0502 }
0503
0504 }  // namespace ActsExamples