Build tubes geometry

source/tdis/tracking/ActsGeometryService.cc

@@ -368,7 +368,7 @@ void tdis::tracking::ActsGeometryService::Init() {
 
 
     /// Return the telescope detector
-    gGeometry = tdis::tracking::buildDetector(
+    gGeometry = tdis::tracking::buildOriginalDetector(
             nominalContext,             // gctx is the detector element dependent geometry context
             m_detector_elements,              // detectorStore is the store for the detector element
             m_plane_positions,          // positions are the positions of different layers in the longitudinal direction

source/tdis/tracking/BuildMtpcDetector.cpp

@@ -22,7 +22,9 @@
 #include "Acts/Geometry/PlaneLayer.hpp"
 #include "Acts/Geometry/TrackingGeometry.hpp"
 #include "Acts/Geometry/TrackingVolume.hpp"
+#include <Acts/Geometry/CylinderLayer.hpp>
 #include "Acts/Material/HomogeneousSurfaceMaterial.hpp"
+#include "Acts/Material/HomogeneousVolumeMaterial.hpp"
 #include "Acts/Material/Material.hpp"
 #include "Acts/Material/MaterialSlab.hpp"
 #include "Acts/Surfaces/RadialBounds.hpp"
@@ -32,8 +34,12 @@
 #include "Acts/Utilities/Logger.hpp"
 #include "BuildMtpcDetector.hpp"
 
-std::unique_ptr<const Acts::TrackingGeometry>
-tdis::tracking::buildDetector(
+
+#include <memory>
+#include <vector>
+/*
+static std::unique_ptr<const Acts::TrackingGeometry>
+tdis::tracking::buildOriginalDetector(
     const typename tdis::tracking::MtpcDetectorElement::ContextType &gctx,
     std::vector<std::shared_ptr<tdis::tracking::MtpcDetectorElement> > &detectorStore,
     const std::vector<double> &positions,
@@ -128,3 +134,127 @@ tdis::tracking::buildDetector(
     // Build and return tracking geometry
     return std::make_unique<Acts::TrackingGeometry>(trackVolume);
 }
+*/
+
+/** This geometry uses cylindrical surfaces with:
+ *    - R (radius) of each cylinder - corresponding to each of mTPC ring center
+ *    - The length of each cylinder in z direction equal to whole detector length
+ *  This way the finding a point on a plane for ACTS Kalman filtering is easy in terms of Z coordinate
+ *  Then each cylinder correspond to MtpcDetectorElement in detectorStore
+ */
+std::unique_ptr<const Acts::TrackingGeometry>
+tdis::tracking::buildCylindricalDetector(
+    const typename MtpcDetectorElement::ContextType& gctx,
+    std::unordered_map<uint32_t, std::shared_ptr<MtpcDetectorElement>>& surfaceStore)
+{
+    using namespace Acts;
+    using namespace Acts::UnitLiterals;
+
+    // Define Argon gas material properties at STP
+    double radiationLength = 19.55_m;       // Radiation length in mm (19.55 m)
+    double interactionLength = 70.0_m;      // Interaction length in mm (70 m)
+    double atomicMass = 39.948;             // Atomic mass of Argon
+    double atomicNumber = 18;               // Atomic number of Argon
+    double massDensity = 1.66e-6_g / 1_mm3; // Mass density in g/mm³
+
+    // Create Argon gas material
+    Material argonGas = Material::fromMassDensity(
+        radiationLength, interactionLength, atomicMass, atomicNumber, massDensity);
+
+    // Define constants
+    double innerRadius = 50_mm;     // 5 cm in mm
+    double outerRadius = 150_mm;    // 15 cm in mm
+    double cylinderLength = 550_mm; // 55 cm in mm
+    double halfLength = cylinderLength / 2.0; // Half-length of the cylinder
+    int numRings = 21;              // Number of concentric rings
+    double radialStep = (outerRadius - innerRadius) / numRings;
+
+    // Create vector to hold cylinder layers
+    std::vector<std::shared_ptr<const CylinderLayer>> cylinderLayers;
+
+    for (int i = 0; i < numRings; ++i) {
+        // Calculate the radius of the current ring (centered within its radial step)
+        double radius = innerRadius + (i + 0.5) * radialStep;
+
+        // Define the cylinder bounds
+        auto cylinderBounds = std::make_shared<const CylinderBounds>(radius, halfLength);
+
+        // Create the detector element
+        uint32_t id = static_cast<uint32_t>(i); // Ring index, 0 for innermost
+
+        auto elementTransform = std::make_shared<const Transform3>(Transform3::Identity());
+
+        // Create the detector element
+        auto detElem = std::make_shared<MtpcDetectorElement>(
+            id,
+            elementTransform,
+            cylinderBounds,
+            radialStep // Thickness of the layer
+            // No surface material assigned
+        );
+
+        // Store the detector element
+        surfaceStore[id] = detElem;
+
+        // Get the surface from the detector element
+        auto cylinderSurface = detElem->surface().getSharedPtr();
+
+        // Create a vector of surfaces
+        std::vector<std::shared_ptr<const Surface>> surfaces;
+        surfaces.push_back(cylinderSurface);
+
+        // Create the SurfaceArray
+        auto surfaceArray = std::make_unique<SurfaceArray>(
+            Transform3(Vector3(0, 0, 0)), // Center at origin
+            std::move(surfaces));
+
+        // Now create the CylinderLayer using the create method
+        auto cylinderLayer = CylinderLayer::create(
+            Transform3::Identity(),        // Transform
+            cylinderBounds,                // Cylinder bounds
+            std::move(surfaceArray),       // Surface array
+            radialStep,                    // Layer thickness
+            nullptr,                       // No approach descriptor
+            LayerType::sensitive);         // Layer type
+
+        // Add the layer to the vector
+        cylinderLayers.push_back(std::move(cylinderLayer));
+    }
+
+    // Create a LayerArray from the cylinder layers
+    Acts::LayerArrayCreator::Config layerArrayCreatorConfig;
+    LayerArrayCreator layerArrayCreator(layerArrayCreatorConfig);
+    auto layerArray = layerArrayCreator.layerArray(
+        GeometryContext(),
+        cylinderLayers,
+        innerRadius,
+        outerRadius,
+        BinningType::arbitrary,
+        BinningValue::binR);
+
+    // Define the cylinder volume bounds (outermost dimensions)
+    auto volumeBounds = std::make_shared<CylinderVolumeBounds>(
+        innerRadius, outerRadius, halfLength);
+
+    // Create the material for the volume
+    auto volumeMaterial = std::make_shared<HomogeneousVolumeMaterial>(argonGas);
+
+    // Create the tracking volume with the layers
+    auto trackingVolume = std::make_shared<TrackingVolume>(
+        Transform3::Identity(),      // No transformation (centered at origin)
+        volumeBounds,                // Volume bounds
+        volumeMaterial,              // Volume material
+        std::move(layerArray),       // Layer array
+        nullptr,                     // No contained volumes
+        MutableTrackingVolumeVector{},
+        "TPCVolume");
+
+    // Create the TrackingGeometry with the tracking volume as the world volume
+    auto trackingGeometry = std::make_unique<TrackingGeometry>(trackingVolume);
+
+    return trackingGeometry;
+}
+
+
+
+

source/tdis/tracking/BuildMtpcDetector.hpp

@@ -23,29 +23,15 @@ class TrackingGeometry;
 
 namespace tdis::tracking {
 
-/// Global method to build the telescope tracking geometry
-///
-/// @param gctx is the detector element dependent geometry context
-/// @param detectorStore is the store for the detector element
-/// @param positions are the positions of different layers in the longitudinal
-///                  direction
-/// @param stereoAngles are the stereo angles of different layers, which are
-///                     rotation angles around the longitudinal (normal)
-///                     direction
-/// @param offsets is the offset (u, v) of the layers in the transverse plane
-/// @param bounds is the surface bound values, i.e. minR and maxR
-/// @param thickness is the material thickness of each layer
-
-/// @param binValue indicates which axis the detector surface normals are
-/// parallel to
-std::unique_ptr<const Acts::TrackingGeometry> buildDetector(
-    const typename MtpcDetectorElement::ContextType& gctx,
-    std::vector<std::shared_ptr<MtpcDetectorElement>>& detectorStore,
-    const std::vector<double>& positions,
-    const std::vector<double>& stereoAngles,
-    const std::array<double, 2>& offsets,
-    const std::array<double, 2>& bounds,
-    double thickness,
-    Acts::BinningValue binValue = Acts::BinningValue::binZ);
+
+/** This geometry uses cylindrical surfaces with:
+ *    - R (radius) of each cylinder - corresponding to each of mTPC ring center
+ *    - The length of each cylinder in z direction equal to whole detector length
+ *  This way the finding a point on a plane for ACTS Kalman filtering is easy in terms of Z coordinate
+ *  Then each cylinder correspond to MtpcDetectorElement in detectorStore
+ */
+std::unique_ptr<const Acts::TrackingGeometry> buildCylindricalDetector(
+        const typename MtpcDetectorElement::ContextType& gctx,
+        std::unordered_map<uint32_t, std::shared_ptr<MtpcDetectorElement>>& surfaceStore);
 
 }  // namespace ActsExamples::Telescope

source/tdis/tracking/MtpcDetectorElement.cpp

@@ -1,27 +1,73 @@
-// This file is part of the Acts project.
-//
-// Copyright (C) 2020-2021 CERN for the benefit of the Acts project
-//
-// This Source Code Form is subject to the terms of the Mozilla Public
-// License, v. 2.0. If a copy of the MPL was not distributed with this
-// file, You can obtain one at http://mozilla.org/MPL/2.0/.
-
-#include "Acts/Surfaces/DiscSurface.hpp"
-#include "Acts/Surfaces/PlaneSurface.hpp"
 #include "MtpcDetectorElement.hpp"
 
+#include "Acts/Geometry/GeometryContext.hpp"
+#include "Acts/Surfaces/CylinderBounds.hpp"
+#include "Acts/Surfaces/Surface.hpp"
 
-tdis::tracking::MtpcDetectorElement::MtpcDetectorElement(
-    std::shared_ptr<const Acts::Transform3> transform,
-    std::shared_ptr<const Acts::DiscBounds> dBounds,
-    double thickness,
-    int gem_plane_id,
-    std::shared_ptr<const Acts::ISurfaceMaterial> material)
-    : Acts::DetectorElementBase(),
-      m_elementTransform(std::move(transform)),
-      m_elementSurface(Acts::Surface::makeShared<Acts::DiscSurface>(dBounds, *this)),
-      m_elementThickness(thickness),
-      m_gem_plane_id(gem_plane_id),
-      m_elementDiscBounds(std::move(dBounds)) {
-  m_elementSurface->assignSurfaceMaterial(std::move(material));
-}
+#include <Acts/Surfaces/CylinderSurface.hpp>
+
+#include "Acts/Definitions/Algebra.hpp"
+
+namespace tdis::tracking {
+    // Implementation
+
+    inline MtpcDetectorElement::MtpcDetectorElement(
+        uint32_t planeId,
+        std::shared_ptr<const Acts::Transform3> transform,
+        std::shared_ptr<const Acts::CylinderBounds> cBounds,
+        double thickness,
+        std::shared_ptr<const Acts::ISurfaceMaterial> material)
+        : m_elementTransform(std::move(transform)),
+          m_elementThickness(thickness),
+          m_elementCylinderBounds(std::move(cBounds)),
+          m_id(planeId)
+    {
+        // Create the cylinder surface
+        m_elementSurface = Acts::Surface::makeShared<Acts::CylinderSurface>(
+            this, // Detector element pointer
+            *m_elementTransform,
+            *m_elementCylinderBounds);
+
+        if (material) {
+            m_elementSurface->assignSurfaceMaterial(material);
+        }
+    }
+
+    inline const Acts::Surface& MtpcDetectorElement::surface() const { return *m_elementSurface; }
+
+    inline Acts::Surface& MtpcDetectorElement::surface() { return *m_elementSurface; }
+
+    inline double MtpcDetectorElement::thickness() const { return m_elementThickness; }
+
+    inline const Acts::Transform3& MtpcDetectorElement::transform(const Acts::GeometryContext& gctx) const
+    {
+        // Check if a different transform than the nominal exists
+        if (!m_alignedTransforms.empty()) {
+            // Cast into the right context object
+            auto alignContext = gctx.get<ContextType>();
+            return (*m_alignedTransforms[alignContext.iov].get());
+        }
+        // Return the standard transform if not found
+        return nominalTransform(gctx);
+    }
+
+    inline const Acts::Transform3& MtpcDetectorElement::nominalTransform(const Acts::GeometryContext& /*gctx*/) const
+    {
+        return *m_elementTransform;
+    }
+
+    inline void MtpcDetectorElement::addAlignedTransform(std::unique_ptr<Acts::Transform3> alignedTransform, unsigned int iov)
+    {
+        // Ensure the vector is large enough
+        auto size = m_alignedTransforms.size();
+        if (iov >= size) {
+            m_alignedTransforms.resize(iov + 1);
+        }
+        m_alignedTransforms[iov] = std::move(alignedTransform);
+    }
+
+    inline const std::vector<std::unique_ptr<Acts::Transform3>>&
+    MtpcDetectorElement::alignedTransforms() const {
+        return m_alignedTransforms;
+    }
+}