u4/tests/U4GenstepTest.cc

0001 #include <stdexcept>
0002 #include <vector>
0003
0004 #include "G4DynamicParticle.hh"
0005 #include "G4Electron.hh"
0006 #include "G4MaterialPropertiesTable.hh"
0007 #include "G4ProcessVector.hh"
0008 #include "G4Scintillation.hh"
0009 #include "G4StepStatus.hh"
0010 #include "G4SystemOfUnits.hh"
0011 #include "G4Track.hh"
0012
0013 #include "U4Genstep.h"
0014
0015 namespace
0016 {
0017 void Require(bool condition, const char* message)
0018 {
0019     if (!condition)
0020         throw std::runtime_error(message);
0021 }
0022
0023 void AddSpectrum(G4MaterialPropertiesTable& mpt, const char* key)
0024 {
0025     const std::vector<G4double> energies{2. * eV, 3. * eV};
0026     const std::vector<G4double> intensities{1., 1.};
0027     mpt.AddProperty(key, energies, intensities);
0028 }
0029
0030 void TestProcessSelection()
0031 {
0032     G4ProcessVector at_rest(2);
0033     G4ProcessVector post_step(5);
0034
0035     const U4Genstep::ProcessRange at_rest_range =
0036         U4Genstep::Processes(fAtRestDoItProc, &at_rest, &post_step);
0037     Require(at_rest_range.processes == &at_rest, "at-rest process vector was not selected");
0038     Require(at_rest_range.count == 2, "at-rest process count was not selected");
0039
0040     const U4Genstep::ProcessRange post_step_range =
0041         U4Genstep::Processes(fPostStepDoItProc, &at_rest, &post_step);
0042     Require(post_step_range.processes == &post_step, "post-step process vector was not selected");
0043     Require(post_step_range.count == 5, "post-step process count was not selected");
0044
0045     const U4Genstep::ProcessRange missing = U4Genstep::Processes(nullptr);
0046     Require(missing.processes == nullptr, "null stepping manager returned a process vector");
0047     Require(missing.count == 0, "null stepping manager returned a process count");
0048 }
0049
0050 void TestTwoComponents()
0051 {
0052     G4MaterialPropertiesTable mpt;
0053     AddSpectrum(mpt, "SCINTILLATIONCOMPONENT1");
0054     AddSpectrum(mpt, "SCINTILLATIONCOMPONENT2");
0055     mpt.AddConstProperty("SCINTILLATIONYIELD1", 0.25);
0056     mpt.AddConstProperty("SCINTILLATIONYIELD2", 0.75);
0057     mpt.AddConstProperty("SCINTILLATIONTIMECONSTANT1", 10. * ns);
0058     mpt.AddConstProperty("SCINTILLATIONTIMECONSTANT2", 20. * ns);
0059
0060     G4Scintillation                    scintillation;
0061     const U4Genstep::ScintillationPlan plan =
0062         U4Genstep::Scintillation(nullptr, &scintillation, &mpt, 11);
0063
0064     Require(plan.num_components == 2, "two-component spectrum was not detected");
0065     Require(plan.counts[0] == 2, "first two-component count is incorrect");
0066     Require(plan.counts[1] == 9, "second component did not receive the remainder");
0067     Require(plan.counts[2] == 0, "inactive third component received photons");
0068     Require(plan.times[0] == 10. * ns, "first time constant is incorrect");
0069     Require(plan.times[1] == 20. * ns, "second time constant is incorrect");
0070     Require(plan.photonCount() == 11, "two-component photon count was not preserved");
0071 }
0072
0073 void TestThreeComponents()
0074 {
0075     G4MaterialPropertiesTable mpt;
0076     AddSpectrum(mpt, "SCINTILLATIONCOMPONENT1");
0077     AddSpectrum(mpt, "SCINTILLATIONCOMPONENT2");
0078     AddSpectrum(mpt, "SCINTILLATIONCOMPONENT3");
0079     mpt.AddConstProperty("SCINTILLATIONYIELD1", 0.2);
0080     mpt.AddConstProperty("SCINTILLATIONYIELD2", 0.3);
0081     mpt.AddConstProperty("SCINTILLATIONYIELD3", 0.5);
0082     mpt.AddConstProperty("SCINTILLATIONTIMECONSTANT1", 10. * ns);
0083     mpt.AddConstProperty("SCINTILLATIONTIMECONSTANT2", 20. * ns);
0084     mpt.AddConstProperty("SCINTILLATIONTIMECONSTANT3", 30. * ns);
0085
0086     G4Scintillation                    scintillation;
0087     const U4Genstep::ScintillationPlan plan =
0088         U4Genstep::Scintillation(nullptr, &scintillation, &mpt, 11);
0089
0090     Require(plan.num_components == 3, "three-component spectrum was not detected");
0091     Require(plan.counts[0] == 2, "first three-component count is incorrect");
0092     Require(plan.counts[1] == 3, "second three-component count is incorrect");
0093     Require(plan.counts[2] == 5, "third three-component count is incorrect");
0094     Require(plan.times[0] == 10. * ns, "first time constant is incorrect");
0095     Require(plan.times[1] == 20. * ns, "second time constant is incorrect");
0096     Require(plan.times[2] == 30. * ns, "third time constant is incorrect");
0097     Require(plan.photonCount() == 10, "three-component truncation does not match Geant4");
0098 }
0099
0100 void TestParticleSpecificComponents()
0101 {
0102     G4MaterialPropertiesTable mpt;
0103     AddSpectrum(mpt, "SCINTILLATIONCOMPONENT1");
0104     AddSpectrum(mpt, "SCINTILLATIONCOMPONENT2");
0105     mpt.AddConstProperty("SCINTILLATIONTIMECONSTANT1", 10. * ns);
0106     mpt.AddConstProperty("SCINTILLATIONTIMECONSTANT2", 20. * ns);
0107     mpt.AddConstProperty("ELECTRONSCINTILLATIONYIELD1", 0.4);
0108     mpt.AddConstProperty("ELECTRONSCINTILLATIONYIELD2", 0.6);
0109     mpt.AddConstProperty("ELECTRONSCINTILLATIONTIMECONSTANT1", 15. * ns);
0110
0111     G4Scintillation scintillation;
0112     scintillation.SetScintillationByParticleType(true);
0113     G4Track track(
0114         new G4DynamicParticle(G4Electron::Definition(), G4ThreeVector(0., 0., 1.), 1. * MeV),
0115         0.,
0116         G4ThreeVector());
0117
0118     const U4Genstep::ScintillationPlan plan =
0119         U4Genstep::Scintillation(&track, &scintillation, &mpt, 10);
0120
0121     Require(plan.counts[0] == 4, "particle-specific first yield was not used");
0122     Require(plan.counts[1] == 6, "particle-specific second yield was not used");
0123     Require(plan.times[0] == 15. * ns, "particle-specific time constant was not used");
0124     Require(plan.times[1] == 20. * ns, "standard time fallback was not used");
0125 }
0126 } // namespace
0127
0128 int main()
0129 {
0130     TestProcessSelection();
0131     TestTwoComponents();
0132     TestThreeComponents();
0133     TestParticleSpecificComponents();
0134     return 0;
0135 }