Geant4/src/Geant4Converters.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 "ActsPlugins/Geant4/Geant4Converters.hpp"
0010
0011 #include "Acts/Definitions/Algebra.hpp"
0012 #include "Acts/Geometry/CylinderVolumeBounds.hpp"
0013 #include "Acts/Material/HomogeneousSurfaceMaterial.hpp"
0014 #include "Acts/Material/Material.hpp"
0015 #include "Acts/Material/MaterialSlab.hpp"
0016 #include "Acts/Surfaces/CylinderBounds.hpp"
0017 #include "Acts/Surfaces/CylinderSurface.hpp"
0018 #include "Acts/Surfaces/DiscSurface.hpp"
0019 #include "Acts/Surfaces/LineBounds.hpp"
0020 #include "Acts/Surfaces/PlaneSurface.hpp"
0021 #include "Acts/Surfaces/RadialBounds.hpp"
0022 #include "Acts/Surfaces/RectangleBounds.hpp"
0023 #include "Acts/Surfaces/StrawSurface.hpp"
0024 #include "Acts/Surfaces/TrapezoidBounds.hpp"
0025
0026 #include <algorithm>
0027 #include <cmath>
0028 #include <cstdlib>
0029 #include <iterator>
0030 #include <numbers>
0031 #include <stdexcept>
0032 #include <utility>
0033 #include <vector>
0034
0035 #include "G4Box.hh"
0036 #include "G4LogicalVolume.hh"
0037 #include "G4Material.hh"
0038 #include "G4RotationMatrix.hh"
0039 #include "G4ThreeVector.hh"
0040 #include "G4Transform3D.hh"
0041 #include "G4Trap.hh"
0042 #include "G4Trd.hh"
0043 #include "G4Tubs.hh"
0044 #include "G4VPhysicalVolume.hh"
0045 #include "G4VSolid.hh"
0046
0047 using namespace Acts;
0048
0049 Transform3 ActsPlugins::Geant4AlgebraConverter::transform(
0050     const G4ThreeVector& g4Trans) {
0051   Transform3 gTransform = Transform3::Identity();
0052   Vector3 scaledTrans =
0053       Vector3(scale * g4Trans[0], scale * g4Trans[1], scale * g4Trans[2]);
0054   gTransform.pretranslate(scaledTrans);
0055   return gTransform;
0056 }
0057
0058 Transform3 ActsPlugins::Geant4AlgebraConverter::transform(
0059     const G4RotationMatrix& g4Rot, const G4ThreeVector& g4Trans) {
0060   // Create the translation
0061   Vector3 translation(scale * g4Trans[0], scale * g4Trans[1],
0062                       scale * g4Trans[2]);
0063   // And the rotation to it
0064   RotationMatrix3 rotation;
0065   rotation << g4Rot.xx(), g4Rot.xy(), g4Rot.xz(), g4Rot.yx(), g4Rot.yy(),
0066       g4Rot.yz(), g4Rot.zx(), g4Rot.zy(), g4Rot.zz();
0067   Transform3 transform = Transform3::Identity();
0068   transform.rotate(rotation);
0069   transform.pretranslate(translation);
0070   return transform;
0071 }
0072
0073 Transform3 ActsPlugins::Geant4AlgebraConverter::transform(
0074     const G4Transform3D& g4Trf) {
0075   auto g4Rot = g4Trf.getRotation();
0076   auto g4Trans = g4Trf.getTranslation();
0077   return transform(g4Rot, g4Trans);
0078 }
0079
0080 Transform3 ActsPlugins::Geant4AlgebraConverter::transform(
0081     const G4VPhysicalVolume& g4PhysVol) {
0082   // Get Rotation and translation
0083   auto g4Translation = g4PhysVol.GetTranslation();
0084   auto g4Rotation = g4PhysVol.GetRotation();
0085
0086   G4Transform3D g4Transform =
0087       (g4Rotation == nullptr)
0088           ? G4Transform3D(CLHEP::HepRotation(), g4Translation)
0089           : G4Transform3D(*g4Rotation, g4Translation);
0090
0091   return transform(g4Transform);
0092 }
0093
0094 std::tuple<std::shared_ptr<CylinderBounds>, double>
0095 ActsPlugins::Geant4ShapeConverter::cylinderBounds(const G4Tubs& g4Tubs) {
0096   using B = CylinderBounds;
0097
0098   std::array<double, B::eSize> tArray = {};
0099   tArray[B::eR] =
0100       static_cast<double>(g4Tubs.GetInnerRadius() + g4Tubs.GetOuterRadius()) *
0101       0.5;
0102   tArray[B::eHalfLengthZ] = static_cast<double>(g4Tubs.GetZHalfLength());
0103   tArray[B::eHalfPhiSector] =
0104       0.5 * static_cast<double>(g4Tubs.GetDeltaPhiAngle());
0105   // Geant fiddles around with user given values, i.e. it would not allow [-PI,
0106   // +PI) as a full segment (has to be [0, 2PI)])
0107   if (std::abs(tArray[B::eHalfPhiSector] - std::numbers::pi) <
0108       std::numeric_limits<double>::epsilon()) {
0109     tArray[B::eAveragePhi] = 0.;
0110   } else {
0111     tArray[B::eAveragePhi] = static_cast<double>(g4Tubs.GetStartPhiAngle()) +
0112                              tArray[B::eHalfPhiSector];
0113   }
0114   double thickness = g4Tubs.GetOuterRadius() - g4Tubs.GetInnerRadius();
0115   auto cBounds = std::make_shared<CylinderBounds>(tArray);
0116   return {std::move(cBounds), thickness};
0117 }
0118
0119 std::tuple<std::shared_ptr<RadialBounds>, double>
0120 ActsPlugins::Geant4ShapeConverter::radialBounds(const G4Tubs& g4Tubs) {
0121   using B = RadialBounds;
0122
0123   std::array<double, B::eSize> tArray = {};
0124   tArray[B::eMinR] = static_cast<double>(g4Tubs.GetInnerRadius());
0125   tArray[B::eMaxR] = static_cast<double>(g4Tubs.GetOuterRadius());
0126   tArray[B::eHalfPhiSector] =
0127       0.5 * static_cast<double>(g4Tubs.GetDeltaPhiAngle());
0128   // Geant fiddles around with user given values, i.e. it would not allow [-PI,
0129   // +PI) as a full segment (has to be [0, 2PI)])
0130   if (std::abs(tArray[B::eHalfPhiSector] - std::numbers::pi) <
0131       std::numeric_limits<double>::epsilon()) {
0132     tArray[B::eAveragePhi] = 0.;
0133   } else {
0134     tArray[B::eAveragePhi] = static_cast<double>(g4Tubs.GetStartPhiAngle()) +
0135                              tArray[B::eHalfPhiSector];
0136   }
0137   double thickness = g4Tubs.GetZHalfLength() * 2;
0138   auto rBounds = std::make_shared<RadialBounds>(tArray);
0139   return {std::move(rBounds), thickness};
0140 }
0141
0142 std::shared_ptr<LineBounds> ActsPlugins::Geant4ShapeConverter::lineBounds(
0143     const G4Tubs& g4Tubs) {
0144   auto r = static_cast<double>(g4Tubs.GetOuterRadius());
0145   auto hlZ = static_cast<double>(g4Tubs.GetZHalfLength());
0146   return std::make_shared<LineBounds>(r, hlZ);
0147 }
0148
0149 std::tuple<std::shared_ptr<RectangleBounds>, std::array<int, 2u>, double>
0150 ActsPlugins::Geant4ShapeConverter::rectangleBounds(const G4Box& g4Box) {
0151   std::vector<double> hG4XYZ = {static_cast<double>(g4Box.GetXHalfLength()),
0152                                 static_cast<double>(g4Box.GetYHalfLength()),
0153                                 static_cast<double>(g4Box.GetZHalfLength())};
0154
0155   auto minAt = std::min_element(hG4XYZ.begin(), hG4XYZ.end());
0156   std::size_t minPos = std::distance(hG4XYZ.begin(), minAt);
0157   double thickness = 2. * hG4XYZ[minPos];
0158
0159   std::array<int, 2u> rAxes = {};
0160   switch (minPos) {
0161     case 0: {
0162       rAxes = {1, 2};
0163     } break;
0164     case 1: {
0165       if (keepAxisOrder) {
0166         rAxes = {0, -2};  // flip for right-handed
0167       } else {
0168         rAxes = {2, 0};  // cyclic positive
0169       }
0170     } break;
0171     case 2: {
0172       rAxes = {0, 1};
0173     } break;
0174     default:  // do nothing
0175       break;
0176   }
0177   auto rBounds = std::make_shared<RectangleBounds>(hG4XYZ[std::abs(rAxes[0u])],
0178                                                    hG4XYZ[std::abs(rAxes[1u])]);
0179   return {std::move(rBounds), rAxes, thickness};
0180 }
0181
0182 std::tuple<std::shared_ptr<TrapezoidBounds>, std::array<int, 2u>, double>
0183 ActsPlugins::Geant4ShapeConverter::trapezoidBounds(const G4Trd& g4Trd) {
0184   // primary parameters
0185   double hlX0 = static_cast<double>(g4Trd.GetXHalfLength1());
0186   double hlX1 = static_cast<double>(g4Trd.GetXHalfLength2());
0187   double hlY0 = static_cast<double>(g4Trd.GetYHalfLength1());
0188   double hlY1 = static_cast<double>(g4Trd.GetYHalfLength2());
0189   double hlZ = static_cast<double>(g4Trd.GetZHalfLength());
0190
0191   std::vector<double> dXYZ = {(hlX0 + hlX1) * 0.5, (hlY0 + hlY1) * 0.5, hlZ};
0192
0193   auto minAt = std::min_element(dXYZ.begin(), dXYZ.end());
0194   std::size_t minPos = std::distance(dXYZ.begin(), minAt);
0195   double thickness = 2. * dXYZ[minPos];
0196
0197   double halfLengthXminY = 0.;
0198   double halfLengthXmaxY = 0.;
0199   double halfLengthY = 0.;
0200
0201   std::array<int, 2u> rAxes = {};
0202   switch (minPos) {
0203     case 0: {
0204       halfLengthXminY = hlY0;
0205       halfLengthXmaxY = hlY1;
0206       halfLengthY = hlZ;
0207       rAxes = {1, 2};
0208     } break;
0209     case 1: {
0210       halfLengthXminY = hlX0;
0211       halfLengthXmaxY = hlX1;
0212       halfLengthY = hlZ;
0213       rAxes = {0, -2};
0214     } break;
0215     case 2: {
0216       if (std::abs(hlY0 - hlY1) < std::abs(hlX0 - hlX1)) {
0217         halfLengthXminY = hlX0;
0218         halfLengthXmaxY = hlX1;
0219         halfLengthY = (hlY0 + hlY1) * 0.5;
0220         rAxes = {0, 1};
0221       } else {
0222         halfLengthXminY = hlY0;
0223         halfLengthXmaxY = hlY1;
0224         halfLengthY = (hlX0 + hlX1) * 0.5;
0225         rAxes = {-1, 0};
0226       }
0227     } break;
0228   }
0229
0230   auto tBounds = std::make_shared<TrapezoidBounds>(
0231       halfLengthXminY, halfLengthXmaxY, halfLengthY);
0232   return {std::move(tBounds), rAxes, thickness};
0233 }
0234
0235 std::tuple<std::shared_ptr<TrapezoidBounds>, std::array<int, 2u>, double>
0236 ActsPlugins::Geant4ShapeConverter::trapezoidBounds(const G4Trap& g4Trap) {
0237   // primary parameters
0238   auto y1 = static_cast<double>(g4Trap.GetYHalfLength1());
0239   auto y2 = static_cast<double>(g4Trap.GetYHalfLength2());
0240   auto x1 = static_cast<double>(g4Trap.GetXHalfLength1());
0241   auto x2 = static_cast<double>(g4Trap.GetXHalfLength2());
0242   auto x3 = static_cast<double>(g4Trap.GetXHalfLength3());
0243   auto x4 = static_cast<double>(g4Trap.GetXHalfLength4());
0244   auto phi = static_cast<double>(g4Trap.GetPhi());
0245   auto theta = static_cast<double>(g4Trap.GetTheta());
0246   auto z = static_cast<double>(g4Trap.GetZHalfLength());
0247
0248   double hlX0 = (x1 + x2) * 0.5;
0249   double hlX1 = 2 * z * tan(theta) * cos(phi) + (x3 + x4) * 0.5;
0250   double hlY0 = y1;
0251   double hlY1 = y2 + 2 * z * tan(theta) * sin(phi);
0252   double hlZ = z;
0253
0254   std::vector<double> dXYZ = {(hlX0 + hlX1) * 0.5, (hlY0 + hlY1) * 0.5, hlZ};
0255
0256   auto minAt = std::ranges::min_element(dXYZ);
0257   std::size_t minPos = std::distance(dXYZ.begin(), minAt);
0258   double thickness = 2. * dXYZ[minPos];
0259
0260   double halfLengthXminY = 0.;
0261   double halfLengthXmaxY = 0.;
0262   double halfLengthY = 0.;
0263
0264   std::array<int, 2u> rAxes = {};
0265   switch (minPos) {
0266     case 0: {
0267       halfLengthXminY = std::min(hlY0, hlY1);
0268       halfLengthXmaxY = std::max(hlY0, hlY1);
0269       halfLengthY = hlZ;
0270       rAxes = {1, 2};
0271     } break;
0272     case 1: {
0273       halfLengthXminY = std::min(hlX0, hlX1);
0274       halfLengthXmaxY = std::max(hlX0, hlX1);
0275       halfLengthY = hlZ;
0276       rAxes = {0, -2};
0277     } break;
0278     case 2: {
0279       if (std::abs(hlY0 - hlY1) < std::abs(hlX0 - hlX1)) {
0280         halfLengthXminY = std::min(hlX0, hlX1);
0281         halfLengthXmaxY = std::max(hlX0, hlX1);
0282         halfLengthY = (hlY0 + hlY1) * 0.5;
0283         rAxes = {0, 1};
0284       } else {
0285         halfLengthXminY = std::min(hlY0, hlY1);
0286         halfLengthXmaxY = std::max(hlY0, hlY1);
0287         halfLengthY = (hlX0 + hlX1) * 0.5;
0288         rAxes = {-1, 0};
0289       }
0290     } break;
0291     default: {
0292       throw std::runtime_error("Geant4Converters: could not convert G4Trap.");
0293     }
0294   }
0295
0296   auto tBounds = std::make_shared<TrapezoidBounds>(
0297       halfLengthXminY, halfLengthXmaxY, halfLengthY);
0298   return std::make_tuple(std::move(tBounds), rAxes, thickness);
0299 }
0300
0301 std::tuple<std::shared_ptr<PlanarBounds>, std::array<int, 2u>, double>
0302 ActsPlugins::Geant4ShapeConverter::planarBounds(const G4VSolid& g4Solid) {
0303   const G4Box* box = dynamic_cast<const G4Box*>(&g4Solid);
0304   if (box != nullptr) {
0305     auto [rBounds, axes, thickness] = rectangleBounds(*box);
0306     return {std::move(rBounds), axes, thickness};
0307   }
0308
0309   const G4Trd* trd = dynamic_cast<const G4Trd*>(&g4Solid);
0310   if (trd != nullptr) {
0311     auto [tBounds, axes, thickness] = trapezoidBounds(*trd);
0312     return {std::move(tBounds), axes, thickness};
0313   }
0314
0315   std::shared_ptr<PlanarBounds> pBounds = nullptr;
0316   std::array<int, 2u> rAxes = {};
0317   double rThickness = 0.;
0318   return {std::move(pBounds), rAxes, rThickness};
0319 }
0320
0321 namespace {
0322 Transform3 axesOriented(const Transform3& toGlobalOriginal,
0323                         const std::array<int, 2u>& axes) {
0324   auto originalRotation = toGlobalOriginal.rotation();
0325   auto colX = originalRotation.col(std::abs(axes[0u]));
0326   auto colY = originalRotation.col(std::abs(axes[1u]));
0327   colX *= std::copysign(1, axes[0u]);
0328   colY *= std::copysign(1, axes[1u]);
0329   Vector3 colZ = colX.cross(colY);
0330
0331   Transform3 orientedTransform = Transform3::Identity();
0332   orientedTransform.matrix().block(0, 0, 3, 1) = colX;
0333   orientedTransform.matrix().block(0, 1, 3, 1) = colY;
0334   orientedTransform.matrix().block(0, 2, 3, 1) = colZ;
0335   orientedTransform.matrix().block(0, 3, 3, 1) = toGlobalOriginal.translation();
0336
0337   return orientedTransform;
0338 }
0339 }  // namespace
0340
0341 std::shared_ptr<Surface> ActsPlugins::Geant4PhysicalVolumeConverter::surface(
0342     const G4VPhysicalVolume& g4PhysVol, const Transform3& toGlobal,
0343     bool convertMaterial, double compressed) {
0344   // Get the logical volume
0345   auto g4LogVol = g4PhysVol.GetLogicalVolume();
0346   auto g4Solid = g4LogVol->GetSolid();
0347
0348   auto assignMaterial = [&](Surface& sf, double moriginal,
0349                             double mcompressed) -> void {
0350     auto g4Material = g4LogVol->GetMaterial();
0351     if (convertMaterial && g4Material != nullptr) {
0352       if (compressed < 0.) {
0353         mcompressed = moriginal;
0354       }
0355       auto surfaceMaterial = Geant4MaterialConverter{}.surfaceMaterial(
0356           *g4Material, moriginal, mcompressed);
0357       sf.assignSurfaceMaterial(std::move(surfaceMaterial));
0358     }
0359   };
0360
0361   // Dynamic cast chain & conversion
0362   std::shared_ptr<Surface> surface = nullptr;
0363
0364   // Into a rectangle
0365   auto g4Box = dynamic_cast<const G4Box*>(g4Solid);
0366   if (g4Box != nullptr) {
0367     if (forcedType == Surface::SurfaceType::Other ||
0368         forcedType == Surface::SurfaceType::Plane) {
0369       auto [bounds, axes, original] =
0370           Geant4ShapeConverter{}.rectangleBounds(*g4Box);
0371       auto orientedToGlobal = axesOriented(toGlobal, axes);
0372       surface = Surface::makeShared<PlaneSurface>(orientedToGlobal,
0373                                                   std::move(bounds));
0374       assignMaterial(*surface, original, compressed);
0375       return surface;
0376     } else {
0377       throw std::runtime_error("Can not convert 'G4Box' into forced shape.");
0378     }
0379   }
0380
0381   // Into a Trapezoid - from Trd
0382   auto g4Trd = dynamic_cast<const G4Trd*>(g4Solid);
0383   if (g4Trd != nullptr) {
0384     if (forcedType == Surface::SurfaceType::Other ||
0385         forcedType == Surface::SurfaceType::Plane) {
0386       auto [bounds, axes, original] =
0387           Geant4ShapeConverter{}.trapezoidBounds(*g4Trd);
0388       auto orientedToGlobal = axesOriented(toGlobal, axes);
0389       surface = Surface::makeShared<PlaneSurface>(orientedToGlobal,
0390                                                   std::move(bounds));
0391       assignMaterial(*surface, original, compressed);
0392       return surface;
0393     } else {
0394       throw std::runtime_error("Can not convert 'G4Trd' into forced shape.");
0395     }
0396   }
0397
0398   // Into a Trapezoid - from Trap
0399   auto g4Trap = dynamic_cast<const G4Trap*>(g4Solid);
0400   if (g4Trap != nullptr) {
0401     if (forcedType == Surface::SurfaceType::Other ||
0402         forcedType == Surface::SurfaceType::Plane) {
0403       auto [bounds, axes, original] =
0404           Geant4ShapeConverter{}.trapezoidBounds(*g4Trap);
0405       auto orientedToGlobal = axesOriented(toGlobal, axes);
0406       surface = Surface::makeShared<PlaneSurface>(orientedToGlobal,
0407                                                   std::move(bounds));
0408       assignMaterial(*surface, original, compressed);
0409       return surface;
0410     } else {
0411       throw std::runtime_error("Can not convert 'G4Trap' into forced shape.");
0412     }
0413   }
0414
0415   // Into a Cylinder, disc or line
0416   auto g4Tubs = dynamic_cast<const G4Tubs*>(g4Solid);
0417   if (g4Tubs != nullptr) {
0418     double diffR = g4Tubs->GetOuterRadius() - g4Tubs->GetInnerRadius();
0419     double diffZ = 2 * g4Tubs->GetZHalfLength();
0420     // Detect if cylinder or disc case
0421     double original = 0.;
0422     if (forcedType == Surface::SurfaceType::Cylinder ||
0423         (diffR < diffZ && forcedType == Surface::SurfaceType::Other)) {
0424       auto [bounds, originalT] = Geant4ShapeConverter{}.cylinderBounds(*g4Tubs);
0425       original = originalT;
0426       surface =
0427           Surface::makeShared<CylinderSurface>(toGlobal, std::move(bounds));
0428     } else if (forcedType == Surface::SurfaceType::Disc ||
0429                forcedType == Surface::SurfaceType::Other) {
0430       auto [bounds, originalT] = Geant4ShapeConverter{}.radialBounds(*g4Tubs);
0431       original = originalT;
0432       surface = Surface::makeShared<DiscSurface>(toGlobal, std::move(bounds));
0433     } else if (forcedType == Surface::SurfaceType::Straw) {
0434       auto bounds = Geant4ShapeConverter{}.lineBounds(*g4Tubs);
0435       surface = Surface::makeShared<StrawSurface>(toGlobal, std::move(bounds));
0436
0437     } else {
0438       throw std::runtime_error("Can not convert 'G4Tubs' into forced shape.");
0439     }
0440     assignMaterial(*surface, original, compressed);
0441     return surface;
0442   }
0443
0444   return nullptr;
0445 }
0446
0447 Material ActsPlugins::Geant4MaterialConverter::material(
0448     const G4Material& g4Material, double compression) {
0449   auto X0 = g4Material.GetRadlen();
0450   auto L0 = g4Material.GetNuclearInterLength();
0451   auto Rho = g4Material.GetDensity();
0452
0453   // Get{A,Z} is only meaningful for single-element materials (according to
0454   // the Geant4 docs). Need to compute average manually.
0455   auto g4Elements = g4Material.GetElementVector();
0456   auto g4Fraction = g4Material.GetFractionVector();
0457   auto g4NElements = g4Material.GetNumberOfElements();
0458   double Ar = 0;
0459   double Z = 0;
0460   if (g4NElements == 1) {
0461     Ar = g4Elements->at(0)->GetN();
0462     Z = g4Material.GetZ();
0463   } else {
0464     for (std::size_t i = 0; i < g4NElements; i++) {
0465       Ar += g4Elements->at(i)->GetN() * g4Fraction[i];
0466       Z += g4Elements->at(i)->GetZ() * g4Fraction[i];
0467     }
0468   }
0469
0470   return Material::fromMassDensity(X0 / compression, L0 / compression, Ar, Z,
0471                                    compression * Rho);
0472 }
0473
0474 std::shared_ptr<HomogeneousSurfaceMaterial>
0475 ActsPlugins::Geant4MaterialConverter::surfaceMaterial(
0476     const G4Material& g4Material, double original, double compressed) {
0477   double compression = original / compressed;
0478   return std::make_shared<HomogeneousSurfaceMaterial>(
0479       MaterialSlab(material(g4Material, compression), compressed));
0480 }
0481
0482 std::shared_ptr<CylinderVolumeBounds>
0483 ActsPlugins::Geant4VolumeConverter::cylinderBounds(const G4Tubs& g4Tubs) {
0484   using C = CylinderVolumeBounds;
0485
0486   std::array<double, C::eSize> tArray = {};
0487   tArray[C::eMinR] = static_cast<double>(g4Tubs.GetInnerRadius());
0488   tArray[C::eMaxR] = static_cast<double>(g4Tubs.GetOuterRadius());
0489   tArray[C::eHalfLengthZ] = static_cast<double>(g4Tubs.GetZHalfLength());
0490   tArray[C::eHalfPhiSector] =
0491       0.5 * static_cast<double>(g4Tubs.GetDeltaPhiAngle());
0492   tArray[C::eAveragePhi] = static_cast<double>(g4Tubs.GetStartPhiAngle());
0493
0494   return std::make_shared<CylinderVolumeBounds>(tArray);
0495 }