Plugins/Geant4/Geant4SurfaceProviderTests.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 <boost/test/unit_test.hpp>
0010
0011 #include "Acts/Detector/LayerStructureBuilder.hpp"
0012 #include "Acts/Geometry/GeometryContext.hpp"
0013 #include "Acts/Utilities/BinningType.hpp"
0014 #include "Acts/Utilities/StringHelpers.hpp"
0015 #include "ActsPlugins/Geant4/Geant4SurfaceProvider.hpp"
0016 #include "ActsTests/CommonHelpers/FloatComparisons.hpp"
0017
0018 #include <filesystem>
0019 #include <memory>
0020 #include <string>
0021
0022 #include "G4Box.hh"
0023 #include "G4GDMLParser.hh"
0024 #include "G4LogicalVolume.hh"
0025 #include "G4NistManager.hh"
0026 #include "G4PVPlacement.hh"
0027 #include "G4RotationMatrix.hh"
0028 #include "G4SystemOfUnits.hh"
0029 #include "G4ThreeVector.hh"
0030 #include "G4Transform3D.hh"
0031
0032 using namespace Acts;
0033 using namespace Acts::Experimental;
0034 using namespace ActsPlugins;
0035
0036 /// @brief Convert Acts binning value to Geant4 axis
0037 /// as Geant4 uses a different axis convention
0038 /// @param bv the Acts binning value
0039 EAxis binToGeant4Axis(const AxisDirection& bv) {
0040   switch (bv) {
0041     case AxisDirection::AxisX:
0042       return EAxis::kXAxis;
0043     case AxisDirection::AxisY:
0044       return EAxis::kYAxis;
0045     case AxisDirection::AxisZ:
0046       return EAxis::kZAxis;
0047     case AxisDirection::AxisR:
0048       return EAxis::kRho;
0049     case AxisDirection::AxisPhi:
0050       return EAxis::kPhi;
0051     default:
0052       throw std::invalid_argument(
0053           "No Geant4 axis conversion for this binning value");
0054   }
0055 }
0056
0057 const int nLayers = 5;
0058 const int nChips = 5;
0059 const int nArms = 2;
0060
0061 const double cellDimX = 0.4 * cm;
0062 const double cellDimY = 0.3 * cm;
0063 const double cellDimZ = 0.1 * cm;
0064
0065 const double armOffset = 10 * cm;
0066
0067 const std::filesystem::path gdmlPath = "two-arms-telescope.gdml";
0068
0069 namespace ActsTests {
0070
0071 BOOST_AUTO_TEST_SUITE(Geant4Suite)
0072
0073 std::tuple<G4VPhysicalVolume*, std::vector<std::string>>
0074 ConstructGeant4World() {
0075   // Get nist material manager
0076   G4NistManager* nist = G4NistManager::Instance();
0077
0078   // Option to switch on/off checking of volumes overlaps
0079   //
0080   G4bool checkOverlaps = true;
0081
0082   //
0083   // World
0084   //
0085   G4double worldSizeXY = 50 * cm;
0086   G4double worldSizeZ = 50 * cm;
0087   G4Material* worldMat = nist->FindOrBuildMaterial("G4_Galactic");
0088
0089   auto solidWorld = new G4Box("World", 0.5 * worldSizeXY, 0.5 * worldSizeXY,
0090                               0.5 * worldSizeZ);
0091
0092   auto logicWorld = new G4LogicalVolume(solidWorld, worldMat, "World");
0093
0094   auto physWorld = new G4PVPlacement(nullptr, G4ThreeVector(), logicWorld,
0095                                      "World", nullptr, false, 0, checkOverlaps);
0096
0097   G4Material* material = nist->FindOrBuildMaterial("G4_He");
0098
0099   std::vector<std::string> names;
0100   for (int nArm = 0; nArm < nArms; nArm++) {
0101     for (int nLayer = 0; nLayer < nLayers; nLayer++) {
0102       for (int nChip = 0; nChip < nChips; nChip++) {
0103         int sign = (nArm == 0) ? 1 : -1;
0104         double posX = sign * (armOffset + nChip * cm);
0105         double posY = 0;
0106         double posZ = nLayer * cm;
0107         G4ThreeVector pos = G4ThreeVector(posX, posY, posZ);
0108         G4ThreeVector dims = G4ThreeVector(cellDimX, cellDimY, cellDimZ);
0109
0110         int cellId = nChips * nLayers * nArm + nChips * nLayer + nChip;
0111         std::string name = "cell" + std::to_string(cellId);
0112
0113         names.push_back(name);
0114
0115         // Box cell
0116         auto solidCell = new G4Box(name, dims[0], dims[1], dims[2]);
0117
0118         G4LogicalVolume* cellLogical =
0119             new G4LogicalVolume(solidCell, material, "cellSensitive");
0120
0121         new G4PVPlacement(nullptr, pos, cellLogical, name, logicWorld, false, 0,
0122                           true);
0123       }
0124     }
0125   }
0126
0127   // Write GDML file
0128   G4GDMLParser parser;
0129   parser.SetOutputFileOverwrite(true);
0130   parser.Write(gdmlPath.string(), physWorld);
0131
0132   return {physWorld, names};
0133 }
0134
0135 auto gctx = GeometryContext();
0136
0137 // Construct the world
0138 auto [physWorld, names] = ConstructGeant4World();
0139
0140 BOOST_AUTO_TEST_CASE(Geant4SurfaceProviderNames) {
0141   /// Read the gdml file and get the world volume
0142   G4GDMLParser parser;
0143   parser.Read(gdmlPath.string(), false);
0144   auto world = parser.GetWorldVolume();
0145
0146   // Default template parameters are fine
0147   // when using names as identifiers
0148   auto spFullCfg = Geant4SurfaceProvider<>::Config();
0149   spFullCfg.g4World = world;
0150   spFullCfg.surfacePreselector =
0151       std::make_shared<Geant4PhysicalVolumeSelectors::NameSelector>(names,
0152                                                                     true);
0153
0154   auto spFull = std::make_shared<Geant4SurfaceProvider<>>(
0155       spFullCfg, Geant4SurfaceProvider<>::kdtOptions());
0156
0157   auto lbFullCfg = LayerStructureBuilder::Config();
0158   lbFullCfg.surfacesProvider = spFull;
0159
0160   auto lbFull = std::make_shared<LayerStructureBuilder>(lbFullCfg);
0161
0162   auto [sFull, vFull, suFull, vuFull] = lbFull->construct(gctx);
0163
0164   BOOST_CHECK_EQUAL(sFull.size(), names.size());
0165   for (int nArm = 0; nArm < nArms; nArm++) {
0166     for (int nLayer = 0; nLayer < nLayers; nLayer++) {
0167       for (int nChip = 0; nChip < nChips; nChip++) {
0168         int sign = (nArm == 0) ? 1 : -1;
0169         double posX = sign * (armOffset + nChip * cm);
0170         double posY = 0;
0171         double posZ = nLayer * cm;
0172         Vector3 pos = Vector3(posX, posY, posZ);
0173
0174         BOOST_CHECK_EQUAL(
0175             sFull.at(nChips * nLayers * nArm + nChips * nLayer + nChip)
0176                 ->center(gctx),
0177             pos);
0178       }
0179     }
0180   }
0181
0182   // Now check that we can extract only
0183   // a subset of the surfaces
0184   std::vector<std::string> leftArmNames(names.begin(),
0185                                         names.begin() + names.size() / 2);
0186
0187   auto spLeftArmCfg = Geant4SurfaceProvider<>::Config();
0188   spLeftArmCfg.g4World = world;
0189   spLeftArmCfg.surfacePreselector =
0190       std::make_shared<Geant4PhysicalVolumeSelectors::NameSelector>(
0191           leftArmNames, true);
0192
0193   auto spLeftArm = std::make_shared<Geant4SurfaceProvider<>>(
0194       spLeftArmCfg, Geant4SurfaceProvider<>::kdtOptions());
0195
0196   auto lbCfg = LayerStructureBuilder::Config();
0197   lbCfg.surfacesProvider = spLeftArm;
0198
0199   auto lbLeftArm = std::make_shared<LayerStructureBuilder>(lbCfg);
0200
0201   auto [sLeftArm, vLeftArm, suLeftArm, vuLeftArm] = lbLeftArm->construct(gctx);
0202
0203   BOOST_CHECK_EQUAL(sLeftArm.size(), leftArmNames.size());
0204   for (int nLayer = 0; nLayer < nLayers; nLayer++) {
0205     for (int nChip = 0; nChip < nChips; nChip++) {
0206       double posX = armOffset + nChip * cm;
0207       double posY = 0;
0208       double posZ = nLayer * cm;
0209       Vector3 pos = Vector3(posX, posY, posZ);
0210
0211       BOOST_CHECK_EQUAL(sLeftArm.at(nChips * nLayer + nChip)->center(gctx),
0212                         pos);
0213     }
0214   }
0215 }
0216
0217 BOOST_AUTO_TEST_CASE(Geant4SurfaceProviderRanges) {
0218   /// Read the gdml file and get the world volume
0219   G4GDMLParser parser;
0220   parser.Read(gdmlPath.string(), false);
0221   auto world = parser.GetWorldVolume();
0222
0223   // 1D selection -- select only the second row
0224   auto sp1DCfg = Geant4SurfaceProvider<1>::Config();
0225   sp1DCfg.g4World = world;
0226
0227   auto kdt1DOpt = Geant4SurfaceProvider<1>::kdtOptions();
0228   kdt1DOpt.range = RangeXD<1, double>();
0229   kdt1DOpt.range[0].set(8, 12);
0230   kdt1DOpt.binningValues = {AxisDirection::AxisZ};
0231
0232   auto sp1D = std::make_shared<Geant4SurfaceProvider<1>>(sp1DCfg, kdt1DOpt);
0233
0234   auto lb1DCfg = LayerStructureBuilder::Config();
0235   lb1DCfg.surfacesProvider = sp1D;
0236
0237   auto lb1D = std::make_shared<LayerStructureBuilder>(lb1DCfg);
0238
0239   auto [s1D, v1D, su1D, vu1D] = lb1D->construct(gctx);
0240
0241   BOOST_CHECK_EQUAL(s1D.size(), nChips * nArms);
0242   for (int nArm = 0; nArm < nArms; nArm++) {
0243     for (int nChip = 0; nChip < nChips; nChip++) {
0244       int sign = (nArm == 0) ? 1 : -1;
0245       double posX = sign * (armOffset + nChip * cm);
0246       double posY = 0;
0247       double posZ = 10;
0248       Vector3 pos = Vector3(posX, posY, posZ);
0249
0250       BOOST_CHECK_EQUAL(s1D.at(nChips * nArm + nChip)->center(gctx), pos);
0251     }
0252   }
0253
0254   // 2D selection -- select only the second row
0255   // of the left arm
0256   auto sp2DCfg = Geant4SurfaceProvider<2>::Config();
0257   sp2DCfg.g4World = world;
0258
0259   auto kdt2DOpt = Geant4SurfaceProvider<2>::kdtOptions();
0260   kdt2DOpt.range = RangeXD<2, double>();
0261   kdt2DOpt.range[0].set(8, 12);
0262   kdt2DOpt.range[1].set(armOffset - 5, armOffset + 100);
0263   kdt2DOpt.binningValues = {AxisDirection::AxisZ};
0264
0265   auto sp2D = std::make_shared<Geant4SurfaceProvider<2>>(sp2DCfg, kdt2DOpt);
0266
0267   auto lb2DCfg = LayerStructureBuilder::Config();
0268   lb2DCfg.surfacesProvider = sp2D;
0269
0270   auto lb2D = std::make_shared<LayerStructureBuilder>(lb2DCfg);
0271
0272   auto [s2D, v2D, su2D, vu2D] = lb2D->construct(gctx);
0273
0274   BOOST_CHECK_EQUAL(s2D.size(), nChips);
0275   for (int nChip = 0; nChip < nChips; nChip++) {
0276     double posX = armOffset + nChip * cm;
0277     double posY = 0, posZ = 10;
0278     Vector3 pos = Vector3(posX, posY, posZ);
0279
0280     BOOST_CHECK_EQUAL(s2D.at(nChip)->center(gctx), pos);
0281   }
0282
0283   // Preselect the left arm based on the position
0284   // and select only the second row
0285   auto sp2DPosCfg = Geant4SurfaceProvider<1>::Config();
0286   sp2DPosCfg.g4World = world;
0287   std::map<unsigned int, std::tuple<double, double>> ranges;
0288
0289   std::array<unsigned int, 3> g4Axes{0};
0290   for (auto& bv :
0291        {AxisDirection::AxisX, AxisDirection::AxisY, AxisDirection::AxisZ}) {
0292     g4Axes[toUnderlying(bv)] = binToGeant4Axis(bv);
0293   }
0294
0295   ranges[g4Axes[0]] = std::make_tuple(armOffset - 5, armOffset + 100);
0296   ranges[g4Axes[1]] = std::make_tuple(-100, 100);
0297   ranges[g4Axes[2]] = std::make_tuple(-100, 100);
0298
0299   sp2DPosCfg.surfacePreselector =
0300       std::make_shared<Geant4PhysicalVolumeSelectors::PositionSelector>(ranges);
0301
0302   auto sp2DPos =
0303       std::make_shared<Geant4SurfaceProvider<1>>(sp2DPosCfg, kdt1DOpt);
0304
0305   auto lb2DPosCfg = LayerStructureBuilder::Config();
0306   lb2DPosCfg.surfacesProvider = sp2DPos;
0307
0308   auto lb2DPos = std::make_shared<LayerStructureBuilder>(lb2DPosCfg);
0309
0310   auto [s2DPos, v2DPos, su2DPos, vu2DPos] = lb2DPos->construct(gctx);
0311
0312   BOOST_CHECK_EQUAL(s2DPos.size(), nChips);
0313   for (int nChip = 0; nChip < nChips; nChip++) {
0314     double posX = armOffset + nChip * cm;
0315     double posY = 0;
0316     double posZ = 10;
0317     Vector3 pos = Vector3(posX, posY, posZ);
0318
0319     BOOST_CHECK_EQUAL(s2DPos.at(nChip)->center(gctx), pos);
0320   }
0321 }
0322
0323 const char* gdml_head_xml =
0324     R"""(<?xml version="1.0" ?>
0325          <gdml xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance" xsi:noNamespaceSchemaLocation="http://service-spi.web.cern.ch/service-spi/app/releases/GDML/schema/gdml.xsd">)""";
0326
0327 BOOST_AUTO_TEST_CASE(Geant4RectangleFromGDML) {
0328   std::ofstream bgdml;
0329   bgdml.open("Plane.gdml");
0330   bgdml << gdml_head_xml;
0331   bgdml << "<solids>" << '\n';
0332   bgdml << "<box name=\"wb\" x=\"100\" y=\"100\" z=\"500\" lunit=\"mm\"/>"
0333         << '\n';
0334   bgdml << "<box name=\"c\" x=\"35\" y=\"55\" z=\"90\" lunit=\"mm\"/>" << '\n';
0335   bgdml << "<box name=\"b\" x=\"25\" y=\"50\" z=\"1\" lunit=\"mm\"/>" << '\n';
0336   bgdml << "</solids>" << '\n';
0337   bgdml << "<structure>" << '\n';
0338   bgdml << "    <volume name=\"b\">" << '\n';
0339   bgdml << "     <materialref ref=\"G4_Fe\"/>" << '\n';
0340   bgdml << "         <solidref ref=\"b\"/>" << '\n';
0341   bgdml << "    </volume>" << '\n';
0342   bgdml << "    <volume name=\"cl\">" << '\n';
0343   bgdml << "         <materialref ref=\"G4_Galactic\"/>" << '\n';
0344   bgdml << "         <solidref ref=\"c\"/>" << '\n';
0345   bgdml << "             <physvol name=\"b_pv\">" << '\n';
0346   bgdml << "                    <volumeref ref=\"b\"/>" << '\n';
0347   bgdml << "                    <position name=\"b_pv_pos\" unit=\"mm\" "
0348            "x=\"0\" y=\"5.\" z=\"0\"/>"
0349         << '\n';
0350   bgdml << "                    <rotation name=\"b_pv_rot\" unit=\"deg\" "
0351            "x=\"-90\" y=\"0\" z=\"0\"/>"
0352         << '\n';
0353   bgdml << "              </physvol>" << '\n';
0354   bgdml << "    </volume>" << '\n';
0355   bgdml << "    <volume name=\"wl\">" << '\n';
0356   bgdml << "         <materialref ref=\"G4_Galactic\"/>" << '\n';
0357   bgdml << "         <solidref ref=\"wb\"/>" << '\n';
0358   bgdml << "             <physvol name=\"cl_pv\">" << '\n';
0359   bgdml << "                    <volumeref ref=\"cl\"/>" << '\n';
0360   bgdml << "                    <rotation name=\"cl_pv_rot\" unit=\"deg\" "
0361            "x=\"-90\" y=\"0\" z=\"0\"/>"
0362         << '\n';
0363   bgdml << "              </physvol>" << '\n';
0364   bgdml << "    </volume>" << '\n';
0365   bgdml << "</structure>" << '\n';
0366   bgdml << "<setup name=\"Default\" version=\"1.0\">" << '\n';
0367   bgdml << "    <world ref=\"wl\"/>" << '\n';
0368   bgdml << "</setup>" << '\n';
0369   bgdml << "</gdml>" << '\n';
0370
0371   bgdml.close();
0372
0373   /// Read the gdml file and get the world volume
0374   G4GDMLParser parser;
0375   parser.Read("Plane.gdml", false);
0376   auto world = parser.GetWorldVolume();
0377
0378   // 1D selection -- select only the second row
0379   auto planeFromGDMLCfg = Geant4SurfaceProvider<1>::Config();
0380   planeFromGDMLCfg.g4World = world;
0381   planeFromGDMLCfg.surfacePreselector =
0382       std::make_shared<Geant4PhysicalVolumeSelectors::NameSelector>(
0383           std::vector<std::string>{"b_pv"}, true);
0384
0385   auto kdt1DOpt = Geant4SurfaceProvider<1>::kdtOptions();
0386   kdt1DOpt.range = RangeXD<1, double>();
0387   kdt1DOpt.range[0].set(-100, 100);
0388   kdt1DOpt.binningValues = {AxisDirection::AxisZ};
0389
0390   auto tContext = GeometryContext();
0391
0392   auto planeProvider =
0393       std::make_shared<Geant4SurfaceProvider<1>>(planeFromGDMLCfg, kdt1DOpt);
0394
0395   auto planes = planeProvider->surfaces(tContext);
0396   BOOST_CHECK_EQUAL(planes.size(), 1u);
0397   CHECK_CLOSE_ABS(planes.front()->center(tContext).z(), 5., 1e-6);
0398 }
0399
0400 BOOST_AUTO_TEST_SUITE_END()
0401
0402 }  // namespace ActsTests