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