Plugins/Geant4/Geant4MaterialStepConverterTests.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/Definitions/Algebra.hpp"
0012 #include "Acts/Definitions/Units.hpp"
0013 #include "Acts/Material/MaterialInteraction.hpp"
0014 #include "ActsPlugins/Geant4/Geant4MaterialStepConverter.hpp"
0015 #include "ActsTests/CommonHelpers/FloatComparisons.hpp"
0016
0017 #include "G4Element.hh"
0018 #include "G4Material.hh"
0019 #include "G4Step.hh"
0020 #include "G4StepPoint.hh"
0021 #include "G4SystemOfUnits.hh"
0022 #include "G4ThreeVector.hh"
0023 #include "G4Track.hh"
0024
0025 using namespace Acts;
0026 using namespace ActsPlugins;
0027 using namespace Acts::UnitLiterals;
0028
0029 namespace {
0030
0031 // Silicon-like single element material and a two element mixture
0032 G4Material* g4Silicon =
0033     new G4Material("TestSilicon", 14., 28.0855 * g / mole, 2.329 * g / cm3);
0034
0035 G4Material* makeMixture() {
0036   auto* h = new G4Element("TestHydrogen", "TH", 1., 1.008 * g / mole);
0037   auto* o = new G4Element("TestOxygen", "TO", 8., 15.999 * g / mole);
0038   auto* water = new G4Material("TestWater", 1.0 * g / cm3, 2);
0039   water->AddElement(h, 2);
0040   water->AddElement(o, 1);
0041   return water;
0042 }
0043 G4Material* g4Water = makeMixture();
0044
0045 /// Make a step of @p length through @p material, starting at @p position
0046 G4Step makeStep(G4Material* material, const G4ThreeVector& position,
0047                 const G4ThreeVector& direction, double length) {
0048   G4Step step;
0049   G4StepPoint* pre = step.GetPreStepPoint();
0050   pre->SetMaterial(material);
0051   pre->SetPosition(position);
0052   pre->SetMomentumDirection(direction);
0053   pre->SetKineticEnergy(1 * GeV);
0054   pre->SetMass(0.);
0055   step.SetStepLength(length);
0056   return step;
0057 }
0058
0059 }  // namespace
0060
0061 namespace ActsTests {
0062
0063 BOOST_AUTO_TEST_SUITE(Geant4Suite)
0064
0065 BOOST_AUTO_TEST_CASE(Geant4MaterialStepConversion) {
0066   Geant4MaterialStepConverter converter(Geant4MaterialStepConverter::Config{});
0067
0068   G4ThreeVector position(10. * mm, 20. * mm, 30. * mm);
0069   G4ThreeVector direction = G4ThreeVector(1., 1., 0.).unit();
0070   G4Step step = makeStep(g4Silicon, position, direction, 2. * mm);
0071
0072   RecordedMaterialTrack track;
0073   BOOST_CHECK(converter.addStep(step, track));
0074   BOOST_REQUIRE_EQUAL(track.second.materialInteractions.size(), 1u);
0075   const MaterialInteraction& interaction =
0076       track.second.materialInteractions.front();
0077
0078   CHECK_CLOSE_ABS(interaction.position, Vector3(10_mm, 20_mm, 30_mm), 1e-9);
0079   CHECK_CLOSE_ABS(interaction.direction, Vector3(1., 1., 0.).normalized(),
0080                   1e-9);
0081   CHECK_CLOSE_REL(interaction.pathCorrection, 2_mm, 1e-9);
0082
0083   // Material properties are stored in single precision
0084   const auto& slab = interaction.materialSlab;
0085   CHECK_CLOSE_REL(slab.thickness(), 2_mm, 1e-9);
0086   CHECK_CLOSE_REL(slab.material().X0(), g4Silicon->GetRadlen() / mm * 1_mm,
0087                   1e-5);
0088   CHECK_CLOSE_REL(slab.material().L0(),
0089                   g4Silicon->GetNuclearInterLength() / mm * 1_mm, 1e-5);
0090   CHECK_CLOSE_REL(slab.material().Ar(), 28.0855, 1e-5);
0091   CHECK_CLOSE_REL(slab.material().Z(), 14., 1e-5);
0092   CHECK_CLOSE_REL(slab.material().massDensity(), 2.329_g / 1_cm3, 1e-5);
0093
0094   // No element fractions unless requested
0095   BOOST_CHECK(interaction.elementZ.empty());
0096   BOOST_CHECK(interaction.elementFrac.empty());
0097 }
0098
0099 BOOST_AUTO_TEST_CASE(Geant4MaterialStepConversionMixture) {
0100   Geant4MaterialStepConverter::Config cfg;
0101   cfg.recordElementFractions = true;
0102   Geant4MaterialStepConverter converter(cfg);
0103
0104   G4Step step =
0105       makeStep(g4Water, G4ThreeVector(), G4ThreeVector(0., 0., 1.), 1. * cm);
0106   RecordedMaterialTrack track;
0107   BOOST_CHECK(converter.addStep(step, track));
0108   BOOST_REQUIRE_EQUAL(track.second.materialInteractions.size(), 1u);
0109   const MaterialInteraction& interaction =
0110       track.second.materialInteractions.front();
0111
0112   // Mass fraction weighted averages
0113   const double* fractions = g4Water->GetFractionVector();
0114   double expectedAr = fractions[0] * 1.008 + fractions[1] * 15.999;
0115   double expectedZ = fractions[0] * 1. + fractions[1] * 8.;
0116   CHECK_CLOSE_REL(interaction.materialSlab.material().Ar(), expectedAr, 1e-5);
0117   CHECK_CLOSE_REL(interaction.materialSlab.material().Z(), expectedZ, 1e-5);
0118
0119   BOOST_REQUIRE_EQUAL(interaction.elementZ.size(), 2u);
0120   BOOST_CHECK_EQUAL(interaction.elementZ[0], 1u);
0121   BOOST_CHECK_EQUAL(interaction.elementZ[1], 8u);
0122   BOOST_REQUIRE_EQUAL(interaction.elementFrac.size(), 2u);
0123   CHECK_CLOSE_REL(interaction.elementFrac[0], fractions[0], 1e-5);
0124   CHECK_CLOSE_REL(interaction.elementFrac[1], fractions[1], 1e-5);
0125 }
0126
0127 BOOST_AUTO_TEST_CASE(Geant4MaterialStepConversionExcluded) {
0128   Geant4MaterialStepConverter::Config cfg;
0129   cfg.excludeMaterials = {"TestWater"};
0130   Geant4MaterialStepConverter converter(cfg);
0131
0132   G4ThreeVector origin;
0133   G4ThreeVector zAxis(0., 0., 1.);
0134   RecordedMaterialTrack track;
0135   BOOST_CHECK(
0136       !converter.addStep(makeStep(g4Water, origin, zAxis, 1. * mm), track));
0137   // A step without material is not recorded
0138   BOOST_CHECK(
0139       !converter.addStep(makeStep(nullptr, origin, zAxis, 1. * mm), track));
0140   BOOST_CHECK(track.second.materialInteractions.empty());
0141   BOOST_CHECK_EQUAL(track.second.materialInX0, 0.);
0142
0143   BOOST_CHECK(
0144       converter.addStep(makeStep(g4Silicon, origin, zAxis, 1. * mm), track));
0145   BOOST_CHECK_EQUAL(track.second.materialInteractions.size(), 1u);
0146 }
0147
0148 BOOST_AUTO_TEST_CASE(Geant4MaterialTrackAccumulation) {
0149   Geant4MaterialStepConverter converter(Geant4MaterialStepConverter::Config{});
0150
0151   G4ThreeVector zAxis(0., 0., 1.);
0152   RecordedMaterialTrack track;
0153   converter.addStep(makeStep(g4Silicon, G4ThreeVector(), zAxis, 1. * mm),
0154                     track);
0155   converter.addStep(
0156       makeStep(g4Water, G4ThreeVector(0., 0., 1. * mm), zAxis, 2. * mm), track);
0157
0158   BOOST_REQUIRE_EQUAL(track.second.materialInteractions.size(), 2u);
0159   double expectedX0 =
0160       1. * mm / g4Silicon->GetRadlen() + 2. * mm / g4Water->GetRadlen();
0161   double expectedL0 = 1. * mm / g4Silicon->GetNuclearInterLength() +
0162                       2. * mm / g4Water->GetNuclearInterLength();
0163   CHECK_CLOSE_REL(track.second.materialInX0, expectedX0, 1e-5);
0164   CHECK_CLOSE_REL(track.second.materialInL0, expectedL0, 1e-5);
0165 }
0166
0167 BOOST_AUTO_TEST_CASE(Geant4MaterialTrackStart) {
0168   Geant4MaterialStepConverter converter(Geant4MaterialStepConverter::Config{});
0169
0170   // The start is taken from the track vertex only. An empty track avoids
0171   // creating a particle, which would need the Geant4 particle tables and
0172   // with them the Geant4 data sets.
0173   G4ThreeVector vertex(1. * mm, 2. * mm, 3. * mm);
0174   G4ThreeVector direction = G4ThreeVector(0., 1., 1.).unit();
0175   G4Track track;
0176   track.SetVertexPosition(vertex);
0177   track.SetVertexMomentumDirection(direction);
0178
0179   RecordedMaterialTrack rmTrack = converter.startTrack(track);
0180   CHECK_CLOSE_ABS(rmTrack.first.first, Vector3(1_mm, 2_mm, 3_mm), 1e-9);
0181   CHECK_CLOSE_ABS(rmTrack.first.second, Vector3(0., 1., 1.).normalized(), 1e-9);
0182   BOOST_CHECK(rmTrack.second.materialInteractions.empty());
0183   BOOST_CHECK_EQUAL(rmTrack.second.materialInX0, 0.);
0184   BOOST_CHECK_EQUAL(rmTrack.second.materialInL0, 0.);
0185 }
0186
0187 BOOST_AUTO_TEST_SUITE_END()
0188
0189 }  // namespace ActsTests