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