g4materials.cc

Go to the documentation of this file.

#include "g4world.h"
 
// geant4
#include "G4Material.hh"
#include "G4Element.hh"
#include "G4Isotope.hh"
#include "gsystemConventions.h"
 
bool G4World::createG4Material(const std::shared_ptr<GMaterial> &gmaterial) {
    auto NISTman = G4NistManager::Instance(); // material G4 Manager
    auto materialName = gmaterial->getName();
 
    // Only build the material if it is not already available in Geant4.
    auto g4material = NISTman->FindMaterial(materialName);
    if (g4material != nullptr) {
        log->info(2, "Material <", materialName, "> already exists in G4NistManager");
        return true;
    }
 
    auto components = gmaterial->getComponents();
    auto amounts = gmaterial->getAmounts();
    bool isChemical = gmaterial->isChemicalFormula();
 
    // Scan material components:
    // return false if any component does not exist yet (caller will retry later).
    for (auto &componentName: components) {
        if (isChemical) {
            if (NISTman->FindOrBuildElement(componentName) == nullptr) {
                log->info(2, "Element <", componentName, ">, needed by ", materialName, ",  not found yet");
                return false;
            } else { log->info(2, "Element <", componentName, "> needed by ", materialName, " now found"); }
        } else {
            if (NISTman->FindOrBuildMaterial(componentName) == nullptr) {
                log->info(2, "Material <", componentName, ">, needed by ", materialName, ",  not found yet");
                return false;
            } else { log->info(2, "Material <", componentName, "> needed by ", materialName, " now found"); }
        }
    }
 
    // Build the composed material from its components.
    auto density = gmaterial->getDensity();
    g4materialsMap[materialName] = new G4Material(materialName, density * CLHEP::g / CLHEP::cm3,
                                                  static_cast<G4int>(components.size()));
 
    if (isChemical) {
        log->info(2, "Building material <", materialName, "> with components:");
        for (size_t i = 0; i < components.size(); i++) {
            log->info(2, "element <", components[i], "> with amount: ", amounts[i]);
        }
 
        for (size_t i = 0; i < components.size(); i++) {
            auto element = NISTman->FindOrBuildElement(components[i]);
            g4materialsMap[materialName]->AddElement(element, static_cast<G4int>(amounts[i]));
        }
    } else {
        log->info(2, "Building material <", materialName, "> with components:");
        for (size_t i = 0; i < components.size(); i++) {
            log->info(2, "material <", components[i], "> with fractional mass: ", amounts[i]);
        }
 
        for (size_t i = 0; i < components.size(); i++) {
            auto material = NISTman->FindOrBuildMaterial(components[i]);
            g4materialsMap[materialName]->AddMaterial(material, amounts[i]);
        }
    }
 
    // optical properties
    auto photonEnergy = gmaterial->getPhotonEnergy();
    if (!photonEnergy.empty()) {
        auto *materialPropertiesTable = new G4MaterialPropertiesTable();
        bool hasOpticalProperties = false;
 
        auto addProperty = [&](const char *propertyName, const std::vector<double> &values) {
            if (values.empty()) { return; }
            if (values.size() != photonEnergy.size()) {
                log->error(gsystem::ERR_GMATERIALOPTICALPROPERTYMISMATCH,
                           "material <", materialName, "> optical property <", propertyName, "> has ",
                           values.size(), " entries but photonEnergy has ", photonEnergy.size());
            }
            materialPropertiesTable->AddProperty(propertyName, photonEnergy, values);
            hasOpticalProperties = true;
        };
 
        auto addConstantProperty = [&](const char *propertyName, double value, bool isSet) {
            if (!isSet) { return; }
            materialPropertiesTable->AddConstProperty(propertyName, value);
            hasOpticalProperties = true;
        };
 
        // optical properties, these
        addProperty("RINDEX", gmaterial->getIndexOfRefraction());
        addProperty("ABSLENGTH", gmaterial->getAbsorptionLength());
        addProperty("REFLECTIVITY", gmaterial->getReflectivity());
        addProperty("EFFICIENCY", gmaterial->getEfficiency());
        // Geant4 11.x renamed the scintillation property keys (formerly FASTCOMPONENT etc.)
        addProperty("SCINTILLATIONCOMPONENT1", gmaterial->getFastComponent());
        addProperty("SCINTILLATIONCOMPONENT2", gmaterial->getSlowComponent());
 
        // scalar properties
        addConstantProperty("SCINTILLATIONYIELD", gmaterial->getScintillationYield(),
                            gmaterial->hasScintillationYield());
        addConstantProperty("RESOLUTIONSCALE", gmaterial->getResolutionScale(),
                            gmaterial->hasResolutionScale());
 
        addConstantProperty("SCINTILLATIONTIMECONSTANT1", gmaterial->getFasttimeConstant() * CLHEP::ns,
                            gmaterial->hasFasttimeConstant());
        addConstantProperty("SCINTILLATIONTIMECONSTANT2", gmaterial->getSlowtimeConstant() * CLHEP::ns,
                            gmaterial->hasSlowtimeConstant());
        addConstantProperty("SCINTILLATIONYIELD1", gmaterial->getYieldratio(),
                            gmaterial->hasYieldratio());
 
        double getBirksConstant = gmaterial->getBirksConstant();
        if (gmaterial->hasBirksConstant()) {
            g4materialsMap[materialName]->GetIonisation()->SetBirksConstant(getBirksConstant);
            hasOpticalProperties = true;
        }
 
        // leaving RAYLEIGH last for consistency checks in
        addProperty("RAYLEIGH", gmaterial->getRayleigh());
 
        if (hasOpticalProperties) {
            g4materialsMap[materialName]->SetMaterialPropertiesTable(materialPropertiesTable);
            log->info(2, "Attached optical material properties table to material <", materialName, ">");
        } else { delete materialPropertiesTable; }
    }
 
 
    return true;
}
 
void G4World::buildDefaultMaterialsElementsAndIsotopes() {
    // Create a small set of commonly-used isotopes/elements/materials if they are missing.
    // These are defined using Geant4 primitives and then registered in the local map for reference.
    int Z, N;
    double a, d, T;
 
    // ----  Hydrogen
 
    // Hydrogen gas material definition (Hydrogen element + state/gas parameters).
    if (G4NistManager::Instance()->FindMaterial(g4system::HGAS_MATERIAL) == nullptr) {
        Z = 1;
        a = 1.01 * CLHEP::g / CLHEP::mole;
        d = 0.00275 * CLHEP::g / CLHEP::cm3;
        T = 50.0 * CLHEP::kelvin;
        auto Hydrogen = new G4Element(g4system::HYDROGEN_ELEMENT, g4system::HYDROGEN_ELEMENT, Z, a);
        g4materialsMap[g4system::HGAS_MATERIAL] = new G4Material(g4system::HGAS_MATERIAL,
                                                       d,
                                                       1,
                                                       kStateGas,
                                                       T);
        g4materialsMap[g4system::HGAS_MATERIAL]->AddElement(Hydrogen, 1);
    }
    log->info(2, "G4World: Hydrogen gas material <", g4system::HGAS_MATERIAL, "> created with density <", d, ">");
 
    // ----  Deuterium
 
    // Deuteron isotope and Deuterium element definition.
    if (G4NistManager::Instance()->FindOrBuildElement(g4system::DEUTERIUM_ELEMENT) == nullptr) {
        Z = 1;
        N = 2;
        a = 2.0141018 * CLHEP::g / CLHEP::mole;
        auto Deuteron = new G4Isotope(g4system::DEUTERON_ISOTOPE, Z, N, a);
 
        // Deuterium element: isotope composition is explicitly set to the Deuteron isotope.
        Deuterium = new G4Element(g4system::DEUTERIUM_ELEMENT, g4system::DEUTERIUM_ELEMENT, 1);
        Deuterium->AddIsotope(Deuteron, 1);
    }
    log->info(2, "G4World: Deuterium element <", g4system::DEUTERIUM_ELEMENT, "> created with density <", d, ">");
 
    // Deuterium gas material.
    if (G4NistManager::Instance()->FindMaterial(g4system::DEUTERIUMGAS_MATERIAL) == nullptr) {
        d = 0.000452 * CLHEP::g / CLHEP::cm3;
        T = 294.25 * CLHEP::kelvin;
        g4materialsMap[g4system::DEUTERIUMGAS_MATERIAL] = new G4Material(g4system::DEUTERIUMGAS_MATERIAL,
                                                               d,
                                                               1,
                                                               kStateGas,
                                                               T);
        g4materialsMap[g4system::DEUTERIUMGAS_MATERIAL]->AddElement(Deuterium, 1);
    }
    log->info(2, "G4World: Deuterium gas material <", g4system::DEUTERIUMGAS_MATERIAL,
              "> created with density <", d, ">");
 
    // Liquid Deuterium material.
    if (G4NistManager::Instance()->FindMaterial(g4system::LD2_MATERIAL) == nullptr) {
        d = 0.169 * CLHEP::g / CLHEP::cm3;
        T = 22.0 * CLHEP::kelvin;
        g4materialsMap[g4system::LD2_MATERIAL] = new G4Material(g4system::LD2_MATERIAL,
                                                      d,
                                                      1,
                                                      kStateLiquid,
                                                      T);
        g4materialsMap[g4system::LD2_MATERIAL]->AddElement(Deuterium, 2);
    }
    log->info(2, "G4World: Liquid Deuterium material <", g4system::LD2_MATERIAL, "> created with density <", d,
              ">");
 
    // Ammonia (ND3) material definition.
    if (G4NistManager::Instance()->FindMaterial(g4system::ND3_MATERIAL) == nullptr) {
        Z = 7;
        a = 14.01 * CLHEP::g / CLHEP::mole;
        d = 1.007 * CLHEP::g / CLHEP::cm3;
        T = 1.0 * CLHEP::kelvin;
        auto Nitrogen = new G4Element(g4system::NITRO_ELEMENT, g4system::NITRO_ELEMENT, Z, a);
        g4materialsMap[g4system::ND3_MATERIAL] = new G4Material(g4system::ND3_MATERIAL,
                                                      d,
                                                      2,
                                                      kStateLiquid,
                                                      T);
        g4materialsMap[g4system::ND3_MATERIAL]->AddElement(Nitrogen, 1);
        g4materialsMap[g4system::ND3_MATERIAL]->AddElement(Deuterium, 3);
    }
    log->info(2, "G4World: Ammonia material <", g4system::ND3_MATERIAL, "> created with density <", d, ">");
 
    // ---- Helium 3
 
    // Helion isotope and Helium3 element definition.
    if (G4NistManager::Instance()->FindOrBuildElement(g4system::HELIUM3_ELEMENT) == nullptr) {
        Z = 2;
        N = 3;
        a = 3.0160293 * CLHEP::g / CLHEP::mole;
        auto Helion = new G4Isotope(g4system::HELION_ISOTOPE, Z, N, a);
 
        // Helium-3 element: isotope composition is explicitly set to the Helion isotope.
        Helium3 = new G4Element(g4system::HELIUM3_ELEMENT, g4system::HELIUM3_ELEMENT, 1);
        Helium3->AddIsotope(Helion, 1);
    }
    log->info(2, "G4World: Helium 3 element <", g4system::HELIUM3_ELEMENT, "> created with density <", d, ">");
 
    // Helium-3 gas material definition.
    if (G4NistManager::Instance()->FindMaterial(g4system::HELIUM3GAS_MATERIAL) == nullptr) {
        // Density at 21.1°C (70°F): 0.1650 kg/m3.
        d = 0.1650 * CLHEP::mg / CLHEP::cm3;
        T = 294.25 * CLHEP::kelvin;
        g4materialsMap[g4system::HELIUM3GAS_MATERIAL] = new G4Material(g4system::HELIUM3GAS_MATERIAL,
                                                             d,
                                                             1,
                                                             kStateGas,
                                                             T);
        g4materialsMap[g4system::HELIUM3GAS_MATERIAL]->AddElement(Helium3, 1);
    }
    log->info(2, "G4World: Helium 3 gas material <", g4system::HELIUM3GAS_MATERIAL, "> created with density <",
              d, ">");
 
    // ---- Tritium
 
    // Triton isotope and Tritium element definition.
    if (G4NistManager::Instance()->FindOrBuildElement(g4system::TRITIUM_ELEMENT) == nullptr) {
        Z = 1;
        N = 3;
        a = 3.0160492 * CLHEP::g / CLHEP::mole;
        auto Triton = new G4Isotope(g4system::TRITON_ISOTOPE, Z, N, a);
 
        Tritium = new G4Element(g4system::TRITIUM_ELEMENT, g4system::TRITIUM_ELEMENT, 1);
        Tritium->AddIsotope(Triton, 1);
    }
    log->info(2, "G4World: Tritium element <", g4system::TRITIUM_ELEMENT, "> created with density <", d, ">");
 
    // Tritium gas material definition.
    if (G4NistManager::Instance()->FindMaterial(g4system::TRITIUMGAS_MATERIAL) == nullptr) {
        d = 0.0034 * CLHEP::g / CLHEP::cm3;
        T = 40.0 * CLHEP::kelvin;
        g4materialsMap[g4system::TRITIUMGAS_MATERIAL] = new G4Material(g4system::TRITIUMGAS_MATERIAL,
                                                             d,
                                                             1,
                                                             kStateGas, T);
        g4materialsMap[g4system::TRITIUMGAS_MATERIAL]->AddElement(Tritium, 1);
    }
    log->info(2, "G4World: Tritium gas material <", g4system::TRITIUMGAS_MATERIAL, "> created with density <",
              d, ">");
 
    // ---- Optical Air
 
    // Air with a refractive index, used by optical volumes (GEMC2's built-in "Air_Opt").
    // Without RINDEX Geant4 kills optical photons at the material boundary.
    if (G4NistManager::Instance()->FindMaterial(g4system::AIROPTICAL_MATERIAL) == nullptr) {
        auto NISTman = G4NistManager::Instance();
        d = 1.29 * CLHEP::mg / CLHEP::cm3;
        g4materialsMap[g4system::AIROPTICAL_MATERIAL] = new G4Material(g4system::AIROPTICAL_MATERIAL, d, 2);
        g4materialsMap[g4system::AIROPTICAL_MATERIAL]->AddMaterial(NISTman->FindOrBuildMaterial("G4_N"),
                                                                   70 * CLHEP::perCent);
        g4materialsMap[g4system::AIROPTICAL_MATERIAL]->AddMaterial(NISTman->FindOrBuildMaterial("G4_O"),
                                                                   30 * CLHEP::perCent);
 
        std::vector<G4double> photonEnergy = {2.034 * CLHEP::eV, 4.136 * CLHEP::eV};
        std::vector<G4double> refractiveIndex = {1.00, 1.00};
        auto airOpticalMPT = new G4MaterialPropertiesTable();
        airOpticalMPT->AddProperty("RINDEX", photonEnergy, refractiveIndex);
        g4materialsMap[g4system::AIROPTICAL_MATERIAL]->SetMaterialPropertiesTable(airOpticalMPT);
    }
    log->info(2, "G4World: Optical air material <", g4system::AIROPTICAL_MATERIAL, "> created with density <",
              d, ">");
 
    // ---- Kryptonite
 
    // Material used to kill every track that touches it (GEMC2's built-in).
    if (G4NistManager::Instance()->FindMaterial(g4system::KRYPTONITE_MATERIAL) == nullptr) {
        d = 1.0e-8 * CLHEP::mg / CLHEP::cm3;
        g4materialsMap[g4system::KRYPTONITE_MATERIAL] = new G4Material(g4system::KRYPTONITE_MATERIAL, d, 1);
        g4materialsMap[g4system::KRYPTONITE_MATERIAL]->AddMaterial(
            G4NistManager::Instance()->FindOrBuildMaterial("G4_Ar"), 100 * CLHEP::perCent);
    }
    log->info(2, "G4World: Kryptonite material <", g4system::KRYPTONITE_MATERIAL, "> created with density <", d,
              ">");
}