Plugins/Geant4/Geant4SurfaceProviderTests.cpp

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