u4/tests/U4MaterialLegacyReemissionTest.cc

0001 /**
0002  * Unit tests for converting legacy optical reemission properties to Geant4
0003  * wavelength-shifting properties.
0004  *
0005  * Here, legacy refers specifically to the `REEMISSIONPROB` property schema
0006  * used by the removed local scintillation process. In that schema, wavelength
0007  * shifting is expressed with an absorption length, `ABSLENGTH`, followed by
0008  * a conditional reemission probability. The conversion translates this
0009  * representation into competing ordinary-absorption and wavelength-shifting
0010  * lengths suitable for Geant4's `G4OpWLS` process:
0011  *
0012  * @code
0013  * ABSLENGTH    -> ABSLENGTH / (1 - REEMISSIONPROB)
0014  * WLSABSLENGTH -> ABSLENGTH / REEMISSIONPROB
0015  * @endcode
0016  *
0017  * These tests cover interpolated inputs on different energy grids, the zero-
0018  * and unit-probability limits, rejection of unsafe varying endpoints,
0019  * incomplete legacy definitions, preservation of authoritative WLS
0020  * properties, emission-spectrum cloning, time-constant propagation, removal
0021  * of consumed legacy properties, and idempotence.
0022  *
0023  * Requirements use an always-active failure helper rather than `assert`, so
0024  * the conversion calls and checks execute in both Debug and Release builds.
0025  */
0026
0027 #include <cmath>
0028 #include <cstdlib>
0029 #include <iostream>
0030 #include <vector>
0031
0032 #include "G4Material.hh"
0033 #include "G4MaterialPropertiesTable.hh"
0034 #include "G4MaterialPropertyVector.hh"
0035 #include "G4SystemOfUnits.hh"
0036
0037 #include "OPTICKS_LOG.hh"
0038 #include "U4Material.hh"
0039
0040 namespace
0041 {
0042 /**
0043  * Terminates the test executable with a diagnostic when a requirement fails.
0044  *
0045  * @param condition requirement result
0046  * @param message diagnostic printed when `condition` is false
0047  */
0048 void Require(bool condition, const char* message)
0049 {
0050     if (condition)
0051         return;
0052     std::cerr << "U4MaterialLegacyReemissionTest FAILED: " << message << std::endl;
0053     std::exit(EXIT_FAILURE);
0054 }
0055
0056 /**
0057  * Compares length- or time-valued quantities using an absolute tolerance.
0058  *
0059  * @param lhs first value
0060  * @param rhs second value
0061  * @param tolerance maximum accepted absolute difference
0062  * @return `true` when the values differ by less than `tolerance`
0063  */
0064 bool Close(G4double lhs, G4double rhs, G4double tolerance = 1e-9 * m)
0065 {
0066     return std::abs(lhs - rhs) < tolerance;
0067 }
0068
0069 /**
0070  * Compares non-zero quantities using a relative tolerance.
0071  *
0072  * @param lhs first value
0073  * @param rhs non-zero reference value
0074  * @param tolerance maximum accepted relative difference
0075  * @return `true` when the relative difference does not exceed `tolerance`
0076  */
0077 bool RelativeClose(G4double lhs, G4double rhs, G4double tolerance = 1e-3)
0078 {
0079     return std::abs(lhs - rhs) / std::abs(rhs) <= tolerance;
0080 }
0081
0082 /**
0083  * Constructs a minimal material with an initially empty property table.
0084  *
0085  * @param name material name
0086  * @return newly allocated material
0087  */
0088 G4Material* MakeMaterial(const char* name)
0089 {
0090     G4Material* material = new G4Material(name, 1., 1.01 * g / mole, 1. * g / cm3);
0091     material->SetMaterialPropertiesTable(new G4MaterialPropertiesTable);
0092     return material;
0093 }
0094
0095 /**
0096  * Constructs a Geant4 material-property vector from parallel arrays.
0097  *
0098  * @param energies photon-energy coordinates
0099  * @param values property values corresponding to `energies`
0100  * @return newly allocated property vector
0101  */
0102 G4MaterialPropertyVector* MakeProperty(
0103     const std::vector<G4double>& energies,
0104     const std::vector<G4double>& values)
0105 {
0106     return new G4MaterialPropertyVector(energies, values);
0107 }
0108
0109 /**
0110  * Verify competing absorption lengths on the union of mismatched grids.
0111  *
0112  * Also verifies emission-component cloning, time-constant propagation,
0113  * consumption of `REEMISSIONPROB`, and idempotence of the conversion.
0114  */
0115 void TestCompetingLengthsAndMismatchedGrids()
0116 {
0117     G4Material*                material = MakeMaterial("U4MaterialLegacyReemissionTestCompetingLengths");
0118     G4MaterialPropertiesTable* mpt = material->GetMaterialPropertiesTable();
0119
0120     const std::vector<G4double> absorptionEnergies = {2. * eV, 4. * eV};
0121     const std::vector<G4double> absorptionValues = {10. * m, 20. * m};
0122     const std::vector<G4double> probabilityEnergies = {2. * eV, 3. * eV, 4. * eV};
0123     const std::vector<G4double> probabilities = {0.25, 0.5, 0.75};
0124     const std::vector<G4double> emissionValues = {1., 2.};
0125
0126     G4MaterialPropertyVector* originalAbsorption =
0127         MakeProperty(absorptionEnergies, absorptionValues);
0128     G4MaterialPropertyVector* scintComponent =
0129         MakeProperty(absorptionEnergies, emissionValues);
0130
0131     mpt->AddProperty("ABSLENGTH", originalAbsorption);
0132     mpt->AddProperty(
0133         "REEMISSIONPROB",
0134         MakeProperty(probabilityEnergies, probabilities),
0135         true);
0136     mpt->AddProperty("SCINTILLATIONCOMPONENT1", scintComponent);
0137     mpt->AddConstProperty("SCINTILLATIONTIMECONSTANT1", 7. * ns);
0138
0139     const bool converted = U4Material::ConvertLegacyReemissionToWLS(material);
0140     Require(converted, "non-zero legacy properties were not converted");
0141     Require(mpt->GetProperty("REEMISSIONPROB") == nullptr, "legacy probability was retained");
0142
0143     G4MaterialPropertyVector* absorption = mpt->GetProperty("ABSLENGTH");
0144     G4MaterialPropertyVector* wlsAbsorption = mpt->GetProperty("WLSABSLENGTH");
0145     G4MaterialPropertyVector* wlsComponent = mpt->GetProperty("WLSCOMPONENT");
0146
0147     Require(absorption != nullptr, "converted ABSLENGTH is missing");
0148     Require(wlsAbsorption != nullptr, "converted WLSABSLENGTH is missing");
0149     Require(wlsComponent != nullptr, "converted WLSCOMPONENT is missing");
0150     Require(wlsComponent != scintComponent, "emission component was not cloned");
0151     Require(absorption->GetVectorLength() > 3, "converted grid was not adaptively subdivided");
0152     Require(wlsAbsorption->GetVectorLength() == absorption->GetVectorLength(),
0153             "converted absorption grids differ");
0154     Require(mpt->ConstPropertyExists("WLSTIMECONSTANT"), "WLSTIMECONSTANT is missing");
0155     Require(Close(mpt->GetConstProperty("WLSTIMECONSTANT"), 7. * ns, 1e-12 * ns),
0156             "WLSTIMECONSTANT was not copied");
0157
0158     Require(Close(absorption->Value(2. * eV), 10. * m / 0.75), "ordinary absorption at 2 eV is wrong");
0159     Require(Close(wlsAbsorption->Value(2. * eV), 40. * m), "WLS absorption at 2 eV is wrong");
0160     Require(Close(absorption->Value(3. * eV), 30. * m), "interpolated ordinary absorption is wrong");
0161     Require(Close(wlsAbsorption->Value(3. * eV), 30. * m), "interpolated WLS absorption is wrong");
0162     Require(Close(absorption->Value(4. * eV), 80. * m), "ordinary absorption at 4 eV is wrong");
0163     Require(Close(wlsAbsorption->Value(4. * eV), 20. * m / 0.75), "WLS absorption at 4 eV is wrong");
0164
0165     for (unsigned i = 0; i < 997; ++i)
0166     {
0167         const G4double energy = (2. + 2. * (i + 0.5) / 997.) * eV;
0168         const G4double originalLength = (10. + 5. * (energy / eV - 2.)) * m;
0169         const G4double originalProbability = 0.25 + 0.25 * (energy / eV - 2.);
0170         const G4double ordinaryRate = 1. / absorption->Value(energy);
0171         const G4double wavelengthShiftRate = 1. / wlsAbsorption->Value(energy);
0172         const G4double totalRate = ordinaryRate + wavelengthShiftRate;
0173         const G4double branchingFraction = wavelengthShiftRate / totalRate;
0174
0175         Require(RelativeClose(totalRate, 1. / originalLength),
0176                 "interpolated total attenuation rate exceeds tolerance");
0177         Require(Close(branchingFraction, originalProbability, 1e-3),
0178                 "interpolated WLS branching fraction exceeds tolerance");
0179     }
0180
0181     const bool convertedAgain = U4Material::ConvertLegacyReemissionToWLS(material);
0182     Require(!convertedAgain, "conversion is not idempotent");
0183 }
0184
0185 /**
0186  * Verify that legacy optical-parent timing wins over scintillation fallbacks.
0187  *
0188  * The removed local process used the first energy coordinate in
0189  * `OpticalCONSTANT` for reemission. Migration must preserve that choice even
0190  * when ordinary scintillation time constants are also present.
0191  */
0192 void TestOpticalTimeConstantPrecedence()
0193 {
0194     G4Material*                 material = MakeMaterial("U4MaterialLegacyReemissionTestOpticalTime");
0195     G4MaterialPropertiesTable*  mpt = material->GetMaterialPropertiesTable();
0196     const std::vector<G4double> energies = {2. * eV, 3. * eV};
0197
0198     mpt->AddProperty("ABSLENGTH", MakeProperty(energies, {10. * m, 10. * m}));
0199     mpt->AddProperty(
0200         "REEMISSIONPROB",
0201         MakeProperty(energies, {0.5, 0.5}),
0202         true);
0203     mpt->AddProperty("SCINTILLATIONCOMPONENT1", MakeProperty(energies, {1., 1.}));
0204     mpt->AddProperty(
0205         "OpticalCONSTANT",
0206         MakeProperty({11. * ns, 12. * ns}, {1., 0.}),
0207         true);
0208     mpt->AddConstProperty("SCINTILLATIONTIMECONSTANT1", 7. * ns);
0209     mpt->AddConstProperty("FASTTIMECONSTANT", 5. * ns, true);
0210
0211     const bool converted = U4Material::ConvertLegacyReemissionToWLS(material);
0212     Require(converted, "legacy optical timing material was not converted");
0213     Require(mpt->ConstPropertyExists("WLSTIMECONSTANT"), "WLSTIMECONSTANT is missing");
0214     Require(Close(mpt->GetConstProperty("WLSTIMECONSTANT"), 11. * ns, 1e-12 * ns),
0215             "WLSTIMECONSTANT did not prefer OpticalCONSTANT");
0216 }
0217
0218 /**
0219  * Verify that an all-zero reemission probability is consumed as a no-op.
0220  *
0221  * Ordinary absorption must retain its original vector and no WLS absorption
0222  * property should be introduced.
0223  */
0224 void TestZeroProbabilityRemoval()
0225 {
0226     G4Material*                 material = MakeMaterial("U4MaterialLegacyReemissionTestZeroProbability");
0227     G4MaterialPropertiesTable*  mpt = material->GetMaterialPropertiesTable();
0228     const std::vector<G4double> energies = {2. * eV, 3. * eV};
0229     const std::vector<G4double> absorptionValues = {10. * m, 20. * m};
0230     const std::vector<G4double> probabilities = {0., 0.};
0231
0232     G4MaterialPropertyVector* absorption = MakeProperty(energies, absorptionValues);
0233     mpt->AddProperty("ABSLENGTH", absorption);
0234     mpt->AddProperty("REEMISSIONPROB", MakeProperty(energies, probabilities), true);
0235
0236     const bool converted = U4Material::ConvertLegacyReemissionToWLS(material);
0237     Require(converted, "zero probability was not handled");
0238     Require(mpt->GetProperty("REEMISSIONPROB") == nullptr, "zero probability was not removed");
0239     Require(mpt->GetProperty("ABSLENGTH") == absorption, "zero probability changed ABSLENGTH");
0240     Require(mpt->GetProperty("WLSABSLENGTH") == nullptr, "zero probability created WLS absorption");
0241 }
0242
0243 /**
0244  * Verify the unit-probability limit where every absorption becomes WLS.
0245  *
0246  * The WLS length must equal the original absorption length while the competing
0247  * ordinary-absorption length becomes effectively infinite.
0248  */
0249 void TestUnitProbabilityConversion()
0250 {
0251     G4Material*                 material = MakeMaterial("U4MaterialLegacyReemissionTestUnitProbability");
0252     G4MaterialPropertiesTable*  mpt = material->GetMaterialPropertiesTable();
0253     const std::vector<G4double> energies = {2. * eV, 3. * eV};
0254     const std::vector<G4double> absorptionValues = {10. * m, 20. * m};
0255     const std::vector<G4double> probabilities = {1., 1.};
0256     const std::vector<G4double> emissionValues = {1., 2.};
0257
0258     mpt->AddProperty("ABSLENGTH", MakeProperty(energies, absorptionValues));
0259     mpt->AddProperty("REEMISSIONPROB", MakeProperty(energies, probabilities), true);
0260     mpt->AddProperty("SCINTILLATIONCOMPONENT1", MakeProperty(energies, emissionValues));
0261
0262     const bool converted = U4Material::ConvertLegacyReemissionToWLS(material);
0263     Require(converted, "unit probability was not converted");
0264     Require(mpt->GetProperty("REEMISSIONPROB") == nullptr, "unit probability was retained");
0265
0266     G4MaterialPropertyVector* absorption = mpt->GetProperty("ABSLENGTH");
0267     G4MaterialPropertyVector* wlsAbsorption = mpt->GetProperty("WLSABSLENGTH");
0268     Require(absorption != nullptr, "unit probability ordinary absorption is missing");
0269     Require(wlsAbsorption != nullptr, "unit probability WLS absorption is missing");
0270     Require(absorption->Value(2. * eV) > 1e100 * m,
0271             "unit probability ordinary absorption is not effectively infinite");
0272     Require(Close(wlsAbsorption->Value(2. * eV), 10. * m),
0273             "unit probability WLS absorption at 2 eV is wrong");
0274     Require(Close(wlsAbsorption->Value(3. * eV), 20. * m),
0275             "unit probability WLS absorption at 3 eV is wrong");
0276 }
0277
0278 /**
0279  * Verify that an incomplete non-zero legacy definition is left intact.
0280  *
0281  * Without an emission spectrum the conversion must fail without partially
0282  * mutating or consuming the legacy material properties.
0283  */
0284 void TestIncompleteDefinitionIsRetained()
0285 {
0286     G4Material*                 material = MakeMaterial("U4MaterialLegacyReemissionTestIncomplete");
0287     G4MaterialPropertiesTable*  mpt = material->GetMaterialPropertiesTable();
0288     const std::vector<G4double> energies = {2. * eV, 3. * eV};
0289     const std::vector<G4double> absorptionValues = {10. * m, 20. * m};
0290     const std::vector<G4double> probabilities = {0.25, 0.5};
0291
0292     G4MaterialPropertyVector* absorption = MakeProperty(energies, absorptionValues);
0293     G4MaterialPropertyVector* probability = MakeProperty(energies, probabilities);
0294     mpt->AddProperty("ABSLENGTH", absorption);
0295     mpt->AddProperty("REEMISSIONPROB", probability, true);
0296
0297     const bool converted = U4Material::ConvertLegacyReemissionToWLS(material);
0298     Require(!converted, "incomplete non-zero definition was reported as converted");
0299     Require(mpt->GetProperty("REEMISSIONPROB") == probability, "failed conversion removed legacy probability");
0300     Require(mpt->GetProperty("ABSLENGTH") == absorption, "failed conversion changed ABSLENGTH");
0301     Require(mpt->GetProperty("WLSABSLENGTH") == nullptr, "failed conversion created WLS absorption");
0302 }
0303
0304 /**
0305  * Verify that a varying probability reaching zero or one is left unconverted.
0306  *
0307  * A Geant4 property vector linearly interpolates its stored interaction
0308  * lengths. Conversion must therefore avoid adjoining a finite length to the
0309  * effectively infinite length required at either probability endpoint.
0310  */
0311 void TestVaryingEndpointProbabilityIsRetained()
0312 {
0313     const std::vector<G4double>              energies = {2. * eV, 4. * eV};
0314     const std::vector<G4double>              absorptionValues = {10. * m, 10. * m};
0315     const std::vector<G4double>              emissionValues = {1., 1.};
0316     const std::vector<std::vector<G4double>> probabilities = {
0317         {0., 0.5},
0318         {0.5, 1.}};
0319     const char* names[] = {
0320         "U4MaterialLegacyReemissionTestVaryingFromZero",
0321         "U4MaterialLegacyReemissionTestVaryingToOne"};
0322
0323     for (std::size_t i = 0; i < probabilities.size(); ++i)
0324     {
0325         G4Material*                material = MakeMaterial(names[i]);
0326         G4MaterialPropertiesTable* mpt = material->GetMaterialPropertiesTable();
0327         G4MaterialPropertyVector*  absorption = MakeProperty(energies, absorptionValues);
0328         G4MaterialPropertyVector*  probability = MakeProperty(energies, probabilities[i]);
0329
0330         mpt->AddProperty("ABSLENGTH", absorption);
0331         mpt->AddProperty("REEMISSIONPROB", probability, true);
0332         mpt->AddProperty(
0333             "SCINTILLATIONCOMPONENT1",
0334             MakeProperty(energies, emissionValues));
0335
0336         const bool converted = U4Material::ConvertLegacyReemissionToWLS(material);
0337         Require(!converted, "varying endpoint probability was reported as converted");
0338         Require(mpt->GetProperty("REEMISSIONPROB") == probability,
0339                 "rejected endpoint conversion removed legacy probability");
0340         Require(mpt->GetProperty("ABSLENGTH") == absorption,
0341                 "rejected endpoint conversion changed ABSLENGTH");
0342         Require(mpt->GetProperty("WLSABSLENGTH") == nullptr,
0343                 "rejected endpoint conversion created WLSABSLENGTH");
0344         Require(mpt->GetProperty("WLSCOMPONENT") == nullptr,
0345                 "rejected endpoint conversion created WLSCOMPONENT");
0346     }
0347 }
0348
0349 /**
0350  * Verify that existing WLS absorption data remains authoritative.
0351  *
0352  * Conversion may complete the missing WLS emission component and time
0353  * constant, but must not replace existing absorption vectors.
0354  */
0355 void TestAuthoritativeWLSIsPreserved()
0356 {
0357     G4Material*                 material = MakeMaterial("U4MaterialLegacyReemissionTestAuthoritativeWLS");
0358     G4MaterialPropertiesTable*  mpt = material->GetMaterialPropertiesTable();
0359     const std::vector<G4double> energies = {2. * eV, 3. * eV};
0360     const std::vector<G4double> absorptionValues = {10. * m, 20. * m};
0361     const std::vector<G4double> probabilities = {0.25, 0.5};
0362     const std::vector<G4double> emissionValues = {1., 2.};
0363     const std::vector<G4double> wlsAbsorptionValues = {30. * m, 40. * m};
0364
0365     G4MaterialPropertyVector* absorption = MakeProperty(energies, absorptionValues);
0366     G4MaterialPropertyVector* wlsAbsorption = MakeProperty(energies, wlsAbsorptionValues);
0367     G4MaterialPropertyVector* scintComponent = MakeProperty(energies, emissionValues);
0368     mpt->AddProperty("ABSLENGTH", absorption);
0369     mpt->AddProperty("WLSABSLENGTH", wlsAbsorption);
0370     mpt->AddProperty("SCINTILLATIONCOMPONENT1", scintComponent);
0371     mpt->AddProperty("REEMISSIONPROB", MakeProperty(energies, probabilities), true);
0372
0373     const bool converted = U4Material::ConvertLegacyReemissionToWLS(material);
0374     Require(converted, "authoritative WLS definition was not completed");
0375     Require(mpt->GetProperty("REEMISSIONPROB") == nullptr, "legacy probability survived authoritative WLS");
0376     Require(mpt->GetProperty("ABSLENGTH") == absorption, "authoritative WLS changed ABSLENGTH");
0377     Require(mpt->GetProperty("WLSABSLENGTH") == wlsAbsorption, "authoritative WLS absorption was replaced");
0378     Require(mpt->GetProperty("WLSCOMPONENT") != nullptr, "authoritative WLS component was not completed");
0379     Require(mpt->GetProperty("WLSCOMPONENT") != scintComponent, "authoritative WLS component was not cloned");
0380     Require(mpt->ConstPropertyExists("WLSTIMECONSTANT"), "authoritative WLS time constant is missing");
0381     Require(Close(mpt->GetConstProperty("WLSTIMECONSTANT"), 0., 1e-12 * ns),
0382             "default WLS time constant is not zero");
0383 }
0384 } // namespace
0385
0386 int main(int argc, char** argv)
0387 {
0388     OPTICKS_LOG(argc, argv);
0389
0390     TestCompetingLengthsAndMismatchedGrids();
0391     TestOpticalTimeConstantPrecedence();
0392     TestZeroProbabilityRemoval();
0393     TestUnitProbabilityConversion();
0394     TestIncompleteDefinitionIsRetained();
0395     TestAuthoritativeWLSIsPreserved();
0396     TestVaryingEndpointProbabilityIsRetained();
0397     return 0;
0398 }