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