Geant4/src/Geant4Converters.cpp

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