TelescopeDetector/src/BuildTelescopeDetectorGen3.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 "Acts/Definitions/Algebra.hpp"
0010 #include "Acts/Definitions/Units.hpp"
0011 #include "Acts/Geometry/Blueprint.hpp"
0012 #include "Acts/Geometry/BlueprintNode.hpp"
0013 #include "Acts/Geometry/ContainerBlueprintNode.hpp"
0014 #include "Acts/Geometry/CuboidVolumeBounds.hpp"
0015 #include "Acts/Geometry/CylinderVolumeBounds.hpp"
0016 #include "Acts/Geometry/DiscLayer.hpp"
0017 #include "Acts/Geometry/GeometryContext.hpp"
0018 #include "Acts/Geometry/ILayerArrayCreator.hpp"
0019 #include "Acts/Geometry/ITrackingVolumeHelper.hpp"
0020 #include "Acts/Geometry/LayerArrayCreator.hpp"
0021 #include "Acts/Geometry/LayerBlueprintNode.hpp"
0022 #include "Acts/Geometry/MaterialDesignatorBlueprintNode.hpp"
0023 #include "Acts/Geometry/PadBlueprintNode.hpp"
0024 #include "Acts/Geometry/PlaneLayer.hpp"
0025 #include "Acts/Geometry/StaticBlueprintNode.hpp"
0026 #include "Acts/Geometry/TrackingGeometry.hpp"
0027 #include "Acts/Geometry/TrackingVolume.hpp"
0028 #include "Acts/Geometry/TrapezoidVolumeBounds.hpp"
0029 #include "Acts/Geometry/VolumeAttachmentStrategy.hpp"
0030 #include "Acts/Material/HomogeneousSurfaceMaterial.hpp"
0031 #include "Acts/Material/Material.hpp"
0032 #include "Acts/Material/MaterialSlab.hpp"
0033 #include "Acts/Material/ProtoSurfaceMaterial.hpp"
0034 #include "Acts/Surfaces/RadialBounds.hpp"
0035 #include "Acts/Surfaces/RectangleBounds.hpp"
0036 #include "Acts/Surfaces/Surface.hpp"
0037 #include "Acts/Surfaces/SurfaceArray.hpp"
0038 #include "Acts/Surfaces/SurfacePlacementBase.hpp"
0039 #include "Acts/Utilities/AxisDefinitions.hpp"
0040 #include "Acts/Utilities/AxisSpec.hpp"
0041 #include "Acts/Utilities/Diagnostics.hpp"
0042 #include "Acts/Utilities/Logger.hpp"
0043 #include "Acts/Utilities/ProtoAxis.hpp"
0044 #include "Acts/Utilities/TransformHelpers.hpp"
0045 #include "Acts/Visualization/GeometryView3D.hpp"
0046 #include "Acts/Visualization/ObjVisualization3D.hpp"
0047 #include "ActsExamples/TelescopeDetector/BuildTelescopeDetector.hpp"
0048 #include "ActsExamples/TelescopeDetector/TelescopeDetectorElement.hpp"
0049
0050 #include <algorithm>
0051 #include <cstddef>
0052 #include <filesystem>
0053 #include <fstream>
0054 #include <iostream>
0055 #include <memory>
0056 #include <stdexcept>
0057 #include <utility>
0058 #include <vector>
0059
0060 using namespace Acts;
0061 using namespace Acts::UnitLiterals;
0062 using Acts::Blueprint;
0063 using Acts::BlueprintNode;
0064 using Acts::BlueprintOptions;
0065 using Acts::LayerBlueprintNode;
0066 using Acts::MaterialDesignatorBlueprintNode;
0067 using Acts::PadBlueprintNode;
0068 using Acts::StaticBlueprintNode;
0069
0070 std::unique_ptr<const Acts::TrackingGeometry>
0071 ActsExamples::buildTelescopeDetectorGen3(
0072     const Acts::GeometryContext& gctx,
0073     std::vector<std::shared_ptr<const Acts::SurfacePlacementBase>>&
0074         detectorStore,
0075     const std::vector<double>& positions,
0076     const std::vector<double>& stereoAngles,
0077     const std::array<double, 2>& offsets, const std::array<double, 2>& bounds,
0078     double thickness, TelescopeSurfaceType surfaceType,
0079     Acts::AxisDirection rotDirection, const std::array<double, 2>& envelope_x,
0080     const std::array<double, 2>& envelope_y,
0081     const std::array<double, 2>& envelope_z, const Acts::Logger& logger) {
0082   using namespace Acts::UnitLiterals;
0083   ACTS_INFO("Building telescope detector Gen3...");
0084
0085   // The rectangle bounds for plane surface
0086   const auto pBounds =
0087       std::make_shared<Acts::RectangleBounds>(bounds[0], bounds[1]);
0088
0089   // Material of the surfaces
0090   Acts::Material silicon = Acts::Material::fromMassDensity(
0091       9.370_cm, 46.52_cm, 28.0855, 14, 2.329_g / 1_cm3);
0092   Acts::MaterialSlab matProp(silicon, thickness);
0093   const auto surfaceMaterial =
0094       std::make_shared<Acts::HomogeneousSurfaceMaterial>(matProp);
0095
0096   // Construct the rotation
0097   // This assumes the direction is AxisX, AxisY or AxisZ. No reset is necessary
0098   // in case of AxisZ
0099   Acts::RotationMatrix3 rotation = Acts::RotationMatrix3::Identity();
0100   if (rotDirection == Acts::AxisDirection::AxisX) {
0101     Acts::AngleAxis3 rot{90._degree, Acts::Vector3::UnitY()};
0102     rotation = rot.toRotationMatrix();
0103   } else if (rotDirection == Acts::AxisDirection::AxisY) {
0104     Acts::AngleAxis3 rot{90._degree, Acts::Vector3::UnitX()};
0105     rotation = rot.toRotationMatrix();
0106   }
0107
0108   Blueprint::Config cfg;
0109   cfg.envelope[AxisDirection::AxisX] = {envelope_x[0] * 1_mm,
0110                                         envelope_x[1] * 1_mm};
0111   cfg.envelope[AxisDirection::AxisY] = {envelope_y[0] * 1_mm,
0112                                         envelope_y[1] * 1_mm};
0113   cfg.envelope[AxisDirection::AxisZ] = {envelope_z[0] * 1_mm,
0114                                         envelope_z[1] * 1_mm};
0115
0116   Blueprint root{cfg};
0117
0118   auto& cubcontainer = root.addCuboidContainer("CuboidContainer", rotDirection);
0119
0120   std::size_t nLayers = positions.size();
0121   for (unsigned int i = 0; i < nLayers; i++) {
0122     // The translation without rotation yet
0123     Acts::Vector3 transVec;
0124     if (rotDirection == Acts::AxisDirection::AxisX) {
0125       transVec = Acts::Vector3(positions[i], offsets[1], -offsets[0]);
0126     } else if (rotDirection == Acts::AxisDirection::AxisY) {
0127       transVec = Acts::Vector3(offsets[0], positions[i], -offsets[1]);
0128     } else {  // AxisZ
0129       transVec = Acts::Vector3(offsets[0], offsets[1], positions[i]);
0130     }
0131     Acts::Translation3 trans(transVec);
0132
0133     // The entire transformation (the coordinate system, whose center is defined
0134     // by trans, will be rotated as well)
0135     Acts::Transform3 trafo(trans * Acts::makeTransform3(rotation));
0136
0137     // rotate around local z axis by stereo angle
0138     auto stereo = stereoAngles[i];
0139     trafo *= Acts::AngleAxis3(stereo, Acts::Vector3::UnitZ());
0140
0141     // Create the detector element
0142     std::shared_ptr<TelescopeDetectorElement> detElement = nullptr;
0143
0144     const auto id =
0145         static_cast<TelescopeDetectorElement::Identifier>(detectorStore.size());
0146
0147     if (surfaceType == TelescopeSurfaceType::Plane) {
0148       detElement = std::make_shared<TelescopeDetectorElement>(
0149           id, std::make_shared<Acts::Transform3>(trafo), pBounds, 1._um,
0150           surfaceMaterial);
0151     } else {
0152       throw std::invalid_argument(
0153           "Only TelescopeSurfaceType::Plane is currently supported, not "
0154           "TelescopeSurfaceType::Disc");
0155     }
0156     detectorStore.push_back(detElement);
0157
0158     std::shared_ptr<CuboidVolumeBounds> layerBounds;
0159
0160     // Get the surface part 1
0161     auto surface = detElement->surface().getSharedPtr();
0162
0163     // Local 2D corner point of the rectangle bounds (4 corners)
0164     std::vector<Acts::Vector2> localCorners = pBounds->vertices();
0165
0166     // Convert into global 3D Coordinates
0167     double maxAbsInPlane1 = 0.0;
0168     double maxAbsInPlane2 = 0.0;
0169
0170     for (const auto& localCorner : localCorners) {
0171       Acts::Vector3 globalCorner =
0172           surface->localToGlobal(gctx, localCorner, Acts::Vector3::Zero());
0173
0174       // Relative to Layer centrum (trans), since layerVol exists unrotated at
0175       // trans
0176       Acts::Vector3 relCorner = globalCorner - trans.translation();
0177
0178       if (rotDirection == Acts::AxisDirection::AxisX) {
0179         maxAbsInPlane1 = std::max(maxAbsInPlane1, std::abs(relCorner.y()));
0180         maxAbsInPlane2 = std::max(maxAbsInPlane2, std::abs(relCorner.z()));
0181       } else if (rotDirection == Acts::AxisDirection::AxisY) {
0182         maxAbsInPlane1 = std::max(maxAbsInPlane1, std::abs(relCorner.x()));
0183         maxAbsInPlane2 = std::max(maxAbsInPlane2, std::abs(relCorner.z()));
0184       } else {
0185         maxAbsInPlane1 = std::max(maxAbsInPlane1, std::abs(relCorner.x()));
0186         maxAbsInPlane2 = std::max(maxAbsInPlane2, std::abs(relCorner.y()));
0187       }
0188     }
0189
0190     constexpr double margin = 5_mm;
0191
0192     if (rotDirection == Acts::AxisDirection::AxisX) {
0193       layerBounds = std::make_shared<CuboidVolumeBounds>(
0194           thickness / 2 + 1_mm, maxAbsInPlane1 + margin,
0195           maxAbsInPlane2 + margin);
0196     } else if (rotDirection == Acts::AxisDirection::AxisY) {
0197       layerBounds = std::make_shared<CuboidVolumeBounds>(
0198           maxAbsInPlane1 + margin, thickness / 2 + 1_mm,
0199           maxAbsInPlane2 + margin);
0200     } else {
0201       layerBounds = std::make_shared<CuboidVolumeBounds>(
0202           maxAbsInPlane1 + margin, maxAbsInPlane2 + margin,
0203           thickness / 2 + 1_mm);
0204     }
0205
0206     auto layerVol = std::make_unique<TrackingVolume>(
0207         Acts::Transform3{trans}, layerBounds, std::format("parent_{:}", i));
0208
0209     // Get the surface part 2
0210     layerVol->addSurface(surface);
0211
0212     auto layerNode = std::make_shared<StaticBlueprintNode>(std::move(layerVol));
0213
0214     cubcontainer.addChild(std::move(layerNode));
0215   }
0216
0217   auto trackingGeometry = root.construct({}, gctx, logger);
0218
0219   return trackingGeometry;
0220 }