gdetectorConstruction.cc

Go to the documentation of this file.

// gdetectorConstruction
#include "gdetectorConstruction.h"
#include "gdetector_options.h"
 
// gemc
#include "gdynamicdigitizationConventions.h"
#include "gtouchableConventions.h"
#include "gsystemConventions.h"
#include "gsystem_options.h"
#include "gDosimeterDigitization.h"
#include "gFluxDigitization.h"
#include "gPhotonDetectorDigitization.h"
#include "gParticleCounterDigitization.h"
#include "gfield_options.h"
 
// geant4
#include "G4SDManager.hh"
#include "G4RunManager.hh"
#include "G4UserLimits.hh"
#include "G4VVisManager.hh"
 
namespace {
bool is_unset_field_name(const std::string& name) {
    return name.empty() || name == guts::SERIALIZED_NULL_TOKEN || name == "not provided";
}
 
class GVisManagerGuard : public G4VVisManager {
public:
    static void set(G4VVisManager* visManager) { SetConcreteInstance(visManager); }
 
    void Draw(const G4Circle&, const G4Transform3D&) override {}
    void Draw(const G4Polyhedron&, const G4Transform3D&) override {}
    void Draw(const G4Polyline&, const G4Transform3D&) override {}
    void Draw(const G4Polymarker&, const G4Transform3D&) override {}
    void Draw(const G4Square&, const G4Transform3D&) override {}
    void Draw(const G4Text&, const G4Transform3D&) override {}
    void Draw2D(const G4Circle&, const G4Transform3D&) override {}
    void Draw2D(const G4Polyhedron&, const G4Transform3D&) override {}
    void Draw2D(const G4Polyline&, const G4Transform3D&) override {}
    void Draw2D(const G4Polymarker&, const G4Transform3D&) override {}
    void Draw2D(const G4Square&, const G4Transform3D&) override {}
    void Draw2D(const G4Text&, const G4Transform3D&) override {}
    void Draw(const G4VTrajectory&) override {}
    void Draw(const G4VHit&) override {}
    void Draw(const G4VDigi&) override {}
    void Draw(const G4LogicalVolume&, const G4VisAttributes&, const G4Transform3D&) override {}
    void Draw(const G4VPhysicalVolume&, const G4VisAttributes&, const G4Transform3D&) override {}
    void Draw(const G4VSolid&, const G4VisAttributes&, const G4Transform3D&) override {}
    void BeginDraw(const G4Transform3D&) override {}
    void EndDraw() override {}
    void BeginDraw2D(const G4Transform3D&) override {}
    void EndDraw2D() override {}
    void GeometryHasChanged() override {}
    void DispatchToModel(const G4VTrajectory&) override {}
    G4bool FilterTrajectory(const G4VTrajectory&) override { return true; }
    G4bool FilterHit(const G4VHit&) override { return true; }
    G4bool FilterDigi(const G4VDigi&) override { return true; }
};
}
 
G4ThreadLocal GMagneto *GDetectorConstruction::gmagneto = nullptr;
G4ThreadLocal std::map<std::string, GSensitiveDetector*>* GDetectorConstruction::tlSDMap = nullptr;
 
GDetectorConstruction::GDetectorConstruction(std::shared_ptr<GOptions> gopts)
    : GBase(gopts, GDETECTOR_LOGGER),
      G4VUserDetectorConstruction(), // Geant4 base class.
      gopt(gopts) {
    // The map is populated with the master geometry before worker SDs are constructed.
    digitization_routines_map = std::make_shared<gdynamicdigitization::dRoutinesMap>();
}
 
// Builds (or rebuilds) the GEMC world and then the Geant4 world.
G4VPhysicalVolume *GDetectorConstruction::Construct() {
    log->debug(NORMAL, "%1", FUNCTION_NAME);
 
    // Delete old geometry objects if they exist.
    // These shared_ptr resets guarantee we won't keep references to stale world objects.
    gworld.reset();
    g4world.reset();
 
    // - if no systems are provided, we just launched gemc: create from options
    // - otherwise, it's a geometry re-load. use existing systems.
    if (gsystems.empty()) {
        log->debug(NORMAL, "%1creating world from options", FUNCTION_NAME);
        gworld = std::make_shared<GWorld>(gopt);
    } else {
        log->debug(NORMAL, "%1creating world from a gsystem vector of size %2", FUNCTION_NAME, gsystems.size());
        gworld = std::make_shared<GWorld>(gopt, cloneSystemDescriptors(gsystems));
    }
 
    // Build Geant4 world (solids, logical and physical volumes) based on the GEMC world.
    g4world = std::make_shared<G4World>(gworld.get(), gopt);
 
    auto nsdetectors = gworld->getSensitiveDetectorsList().size();
 
    // tally with number :
    log->info(
        0, "Tally summary: \n - %1 volumes\n - %2 geant4 built volumes\n - %3 sensitive detectors\n",
        gworld->get_number_of_volumes() - 1, g4world->number_of_volumes(), nsdetectors);
 
    // Logical volumes are shared by all workers, so load their digitization metadata and
    // install user limits once while constructing the master geometry.
    if (digiplugins_need_reload) {
        loadDigitizationPlugins();
        digiplugins_need_reload = false;
    }
    assignUserLimits();
 
    // Return the physical volume for the ROOT world volume.
    return g4world->getG4Volume(gsystem::ROOTWORLDGVOLUMENAME)->getPhysical();
}
 
// Installs sensitive detectors and EM fields for the constructed geometry.
void GDetectorConstruction::ConstructSDandField() {
    auto sdManager = G4SDManager::GetSDMpointer();
 
    log->debug(NORMAL, "%1", FUNCTION_NAME);
 
    // Deactivate all SDs this thread registered in prior geometry loads.
    // G4SDManager retains SDs indefinitely across reloads; an active stale SD has
    // Initialize() called at the start of every event even when no volume uses it,
    // producing orphan hit collections whose names are absent from the digitization map.
    // We deactivate them here and reactivate only those needed for the current geometry.
    if (!tlSDMap) {
        tlSDMap = new std::map<std::string, GSensitiveDetector*>();
    }
    for (auto& [name, sd] : *tlSDMap) {
        sd->Activate(false);
    }
 
    // Local cache of sensitive detectors keyed by digitization name.
    // Multiple volumes can share the same digitization name and therefore reuse one SD instance.
    std::unordered_map<std::string, GSensitiveDetector *> sensitiveDetectorsMap;
 
    // --- Electromagnetic field reset (-no_field) --------------------------------------------------
    // Parse -no_field: a gvolume name, a whitespace/comma-separated list of gvolume names, or 'all'.
    // 'all' resets every per-volume field and the global field; a name resets only that volume's field.
    bool                  disable_all_fields = false;
    std::set<std::string> no_field_volumes;
    {
        auto no_field_value = gopt->getOptionalScalarString(gfields::NO_FIELD_OPTION);
        if (no_field_value && !no_field_value->empty()) {
            if (*no_field_value == gfields::NO_FIELD_ALL) { disable_all_fields = true; }
            else {
                for (auto& c : *no_field_value) { if (c == ',') { c = ' '; } }
                for (const auto& v : gutilities::getStringVectorFromString(*no_field_value)) {
                    no_field_volumes.insert(v);
                }
            }
        }
    }
 
    // Global field name, honoring -no_field=all.
    const bool no_system_selected = gsystems.empty() && gsystem::getSystems(gopt).empty();
    auto global_field_name = gopt->getOptionalScalarString(gfields::GLOBAL_FIELD_OPTION).value_or("");
    if (no_system_selected && !disable_all_fields && is_unset_field_name(global_field_name)) {
        std::vector<std::string> configured_fields;
        for (const auto& field_definition : gfields::get_GFieldDefinition(gopt)) {
            if (!is_unset_field_name(field_definition.name)) {
                configured_fields.push_back(field_definition.name);
            }
        }
        if (configured_fields.size() == 1) {
            global_field_name = configured_fields.front();
            log->info(
                1,
                "No explicit global field selected: using the only configured field <%1> on the ROOT world "
                "volume.", global_field_name);
        }
    }
    const bool global_field_set = !disable_all_fields && !is_unset_field_name(global_field_name);
 
    // First pass: collect the fields that are actually used so only their plugins and maps are loaded.
    // Volumes (or the global field) reset via -no_field are excluded here and never trigger a load.
    std::set<std::string> required_fields;
    std::set<std::string> matched_no_field;
    if (!disable_all_fields) {
        for (const auto &[systemName, gsystemPtr]: *gworld->getSystemsMap()) {
            for (const auto &[volumeName, gvolumePtr]: gsystemPtr->getGVolumesMap()) {
                const auto &g4name = gvolumePtr->getG4Name();
                const bool  reset  = no_field_volumes.count(volumeName) || no_field_volumes.count(g4name);
                if (no_field_volumes.count(volumeName)) { matched_no_field.insert(volumeName); }
                if (no_field_volumes.count(g4name)) { matched_no_field.insert(g4name); }
 
                const auto &field_name = gvolumePtr->getEMField();
                if (!field_name || reset) { continue; }
                required_fields.insert(*field_name);
            }
        }
        for (const auto &name: no_field_volumes) {
            if (!matched_no_field.count(name)) {
                log->warning("-%1: volume <%2> not found in the geometry.", gfields::NO_FIELD_OPTION, name);
            }
        }
    }
    if (global_field_set) { required_fields.insert(global_field_name); }
 
    // Build the magnetic-field registry only when at least one field is needed: if every field was reset
    // (e.g. -no_field=all) no plugin or field map is loaded at all.
    if (!required_fields.empty() && gmagneto == nullptr) { gmagneto = new GMagneto(gopt, required_fields); }
 
    // Loop over all systems and their volumes.
    for (const auto &[systemName, gsystemPtr]: *gworld->getSystemsMap()) {
        for (const auto &[volumeName, gvolumePtr]: gsystemPtr->getGVolumesMap()) {
            auto const &digitizationName = gvolumePtr->getDigitization();
            auto const &g4name = gvolumePtr->getG4Name();
 
            // Component volumes are boolean-operation building blocks: they have a
            // solid but no logical volume, no field, and no sensitivity.
            if (gvolumePtr->getMaterial() == gsystem::GSYSTEMCOMPONENTMATERIAL) { continue; }
 
            auto *g4volume = g4world->getG4Volume(g4name)->getLogical();
 
            // Ensure the Geant4 logical volume exists.
            // Some GEMC volumes can be "copy-of" another volume; in that case, reuse the
            // referenced Geant4 logical volume rather than failing.
            if (g4volume == nullptr) {
                const auto& copyOf = gvolumePtr->getCopyOf();
                if (copyOf) {
                    auto gsystem = gvolumePtr->getSystem();
                    auto volume_copy = gsystem + "/" + *copyOf;
                    auto copyG4Volume = g4world->getG4Volume(volume_copy)->getLogical();
                    if (copyG4Volume != nullptr) { g4volume = copyG4Volume; } else {
                        log->error(
                            gsystem::ERR_GVOLUMENOTFOUND, "%1 Logical volume copy <%2> not found.",
                            FUNCTION_NAME, volume_copy);
                    }
                }
            }
            if (g4volume == nullptr) {
                log->error(
                    gsystem::ERR_GVOLUMENOTFOUND, "%1 Logical volume <%2> not found.", FUNCTION_NAME, g4name);
            }
 
            // Skip volumes with no digitization.
            if (digitizationName) {
                const std::string& digitization = *digitizationName;
 
                // Obtain (or create) the sensitive detector for this digitization name.
                // We reuse an existing SD object already in G4SDManager rather than creating a
                // new one because AddNewDetector() keeps the OLD object on duplicate names (DET1010).
                // If we created a new SD, G4SDManager would call Initialize() on the stale object
                // while SetSensitiveDetector() pointed the logical volume to the new one — leaving
                // the new SD's gHitsCollection uninitialized when ProcessHits() is called.
                if (!sensitiveDetectorsMap.contains(digitization)) {
                    // Reuse a previously registered SD for this name if one exists on this thread.
                    // AddNewDetector() silently keeps the OLD object on duplicate names (DET1010);
                    // reusing the same pointer avoids that and ensures G4SDManager's Initialize()
                    // fires on the same object that SetSensitiveDetector() wires to the volumes.
                    auto tlIt = tlSDMap->find(digitization);
                    if (tlIt != tlSDMap->end()) {
                        log->info(
                            2, "Reusing existing sensitive detector <%1> for volume <%2>", digitization,
                            g4name);
                        tlIt->second->resetTouchableMap();
                        tlIt->second->Activate(true);
                        sensitiveDetectorsMap[digitization] = tlIt->second;
                    } else {
                        log->info(
                            2, "Creating new sensitive detector <%1> for volume <%2>", digitization, g4name);
                        auto* newSD = new GSensitiveDetector(digitization, gopt);
                        sdManager->AddNewDetector(newSD);
                        sensitiveDetectorsMap[digitization] = newSD;
                        (*tlSDMap)[digitization] = newSD;
                    }
                } else {
                    log->info(
                        2, "Sensitive detector <%1> is already created and available for volume <%2>",
                        digitization, g4name);
                }
 
                // Register the volume touchable with the sensitive detector.
                // The touchable encodes identity and dimension metadata needed by digitization.
                const auto &vdimensions = gvolumePtr->getDetectorDimensions();
                const auto &identity = gvolumePtr->getGIdentity();
                if (!identity) {
                    log->error(
                        gsystem::ERR_GWRONGNUMBEROFPARS,
                        "Sensitive volume <%1> has digitization <%2> but no identifier.", g4name, digitization);
                }
                const auto &mass = g4volume->GetMass();
                auto this_gtouchable = std::make_shared<
                    GTouchable>(gopt, digitization, *identity, vdimensions, mass);
                sensitiveDetectorsMap[digitization]->registerGVolumeTouchable(g4name, this_gtouchable);
 
                // Attach the SD to the logical volume (no AddNewDetector call needed for reused SDs).
                g4volume->SetSensitiveDetector(sensitiveDetectorsMap[digitization]);
 
                //  auto maxStep =
 
                //g4volume->SetUserLimits(new G4UserLimits(0.1*mm, 0.1*mm));
 
                log->info(
                    2, "Logical Volume  <%1> has been successfully assigned to SD.%2", g4name,
                    sensitiveDetectorsMap[digitization]);
            }
 
            // Process electromagnetic fields.
            // If a volume declares an EM field, install a per-volume field manager configured by the
            // named field map, unless that field was reset via -no_field (a matching name, or =all).
            const auto &field_name           = gvolumePtr->getEMField();
            const bool  volume_field_present = field_name.has_value();
            const bool  volume_field_reset   = disable_all_fields ||
                no_field_volumes.count(volumeName) || no_field_volumes.count(g4name);
            if (volume_field_present && volume_field_reset) {
                log->info(
                    2, "Volume <%1> field <%2> reset by -%3: no field installed.", volumeName, *field_name,
                    gfields::NO_FIELD_OPTION);
            } else if (volume_field_present) {
                log->info(
                    2, "Volume <%1> has field: <%2>. Looking into field map definitions.", volumeName,
                    *field_name);
                log->info(2, "Setting field manager for volume <%1> with field <%2>", g4name, *field_name);
                g4world->setFieldManagerForVolume(g4name, gmagneto->getFieldMgr(*field_name).get(), true);
            }
        }
    }
 
    // Process the global field, if requested and not reset via -no_field=all.
    // A global field is associated with the ROOT world volume and propagated to all daughters, so it
    // applies everywhere a more specific per-volume field has not been installed.
    if (global_field_set) {
        log->info(
            2, "Setting global field manager for the ROOT world volume <%1> with field <%2>",
            gsystem::ROOTWORLDGVOLUMENAME, global_field_name);
        g4world->setFieldManagerForVolume(gsystem::ROOTWORLDGVOLUMENAME,
                                          gmagneto->getFieldMgr(global_field_name).get(), true);
    } else if (disable_all_fields && global_field_name != "" &&
        !is_unset_field_name(global_field_name)) {
        log->info(
            2, "Global field <%1> reset by -%2=%3: none installed.", global_field_name,
            gfields::NO_FIELD_OPTION, gfields::NO_FIELD_ALL);
    }
 
    // Bind each digitization routine to its corresponding sensitive detector.
    const auto sdetectors = gworld->getSensitiveDetectorsList();
    for (auto &sdname: sdetectors) {
        auto digitization_routine = digitization_routines_map->at(sdname);
 
        sensitiveDetectorsMap[sdname]->assign_digi_routine(digitization_routine);
        log->info(
            1, "Digitization routine <%1> has been successfully assigned to SD.%2", sdname,
            sensitiveDetectorsMap[sdname]);
    }
}
 
 
// Installs readout step limits once on the master-owned logical volumes.
void GDetectorConstruction::assignUserLimits() {
    const auto sdetectors = gworld->getSensitiveDetectorsList();
    for (const auto& sdname : sdetectors) {
        const auto digitization_routine = digitization_routines_map->at(sdname);
        const auto max_step = digitization_routine->readoutSpecs->getMaxStep();
        auto limits = std::make_unique<G4UserLimits>(max_step);
        auto* limits_ptr = limits.get();
        user_limits.emplace_back(std::move(limits));
 
        for (const auto& [system_name, system] : *gworld->getSystemsMap()) {
            for (const auto& [volume_name, volume] : system->getGVolumesMap()) {
                const auto& digitization_name = volume->getDigitization();
                if (!digitization_name || *digitization_name != sdname) { continue; }
 
                const auto& g4name = volume->getG4Name();
                const auto* g4volume = g4world->getG4Volume(g4name);
                if (g4volume == nullptr || g4volume->getLogical() == nullptr) {
                    log->error(
                        gsystem::ERR_GVOLUMENOTFOUND,
                        "%1 Logical volume <%2> not found while assigning user limits.", FUNCTION_NAME, g4name);
                    continue;
                }
 
                g4volume->getLogical()->SetUserLimits(limits_ptr);
                log->info(1, "Setting G4UserLimits for volume <%1> with maxStep <%2>", g4name, max_step);
            }
        }
    }
}
 
// Loads (or dynamically resolves) the digitization routine for each sensitive detector.
void GDetectorConstruction::loadDigitizationPlugins() {
    // Clear stale entries so a reloaded geometry always gets fresh routines.
    // emplace() would silently skip existing keys, causing hits not to be recorded
    // when the geometry is reloaded with the same SD names.
    digitization_routines_map->clear();
 
    // Resolve the variation each routine uses to load constants / translation tables.
    // By default a routine follows the variation of the gsystem it belongs to; the
    // digitization_variation option, when set, overrides that for every routine.
    const auto digiVariationOverride = gopt->getOptionalScalarString("digitization_variation");
 
    // Map each digitization (sensitive-detector) name to its system's variation.
    std::map<std::string, std::string> systemVariationFor;
    for (const auto& [systemName, gsystemPtr] : *gworld->getSystemsMap()) {
        for (const auto& [volumeName, gvolumePtr] : gsystemPtr->getGVolumesMap()) {
            const auto& digiName = gvolumePtr->getDigitization();
            if (digiName) {
                systemVariationFor.emplace(*digiName, gsystemPtr->getVariation());
            }
        }
    }
 
    const auto sdetectors = gworld->getSensitiveDetectorsList();
 
    for (auto &sdname: sdetectors) {
        if (sdname == gtouchable::FLUXNAME) {
            log->info(1, "Loading flux digitization plugin for routine <%1>", sdname);
            digitization_routines_map->emplace(sdname, std::make_shared<GFluxDigitization>(gopt));
        } else if (sdname == gtouchable::GPHOTON_DETECTORNAME) {
            log->info(1, "Loading gPhotonDetector digitization plugin for routine <%1>", sdname);
            digitization_routines_map->emplace(sdname, std::make_shared<GPhotonDetectorDigitization>(gopt));
        } else if (sdname == gtouchable::COUNTERNAME) {
            log->info(1, "Loading particle counter digitization plugin for routine <%1>", sdname);
            digitization_routines_map->emplace(sdname, std::make_shared<GParticleCounterDigitization>(gopt));
        } else if (sdname == gtouchable::DOSIMETERNAME) {
            log->info(1, "Loading dosimeter digitization plugin for routine <%1>", sdname);
            digitization_routines_map->emplace(sdname, std::make_shared<GDosimeterDigitization>(gopt));
        } else {
            // if it's not in the map already, add it
            log->info(0, "Loading new digitization plugin for routine <%1>", sdname);
            digitization_routines_map->emplace(sdname, gdynamicdigitization::load_dynamicRoutine(sdname, gopt));
        }
 
        // Ensure each routine uses the correct logger and is configured for readout.
        digitization_routines_map->at(sdname)->set_loggers(gopt);
 
        // Resolve and store the variation used when loading this routine's constants/TT:
        // the digitization_variation option when set, otherwise the routine's gsystem variation.
        std::string variation = "default";
        if (const auto it = systemVariationFor.find(sdname); it != systemVariationFor.end()) {
            variation = it->second;
        }
        if (digiVariationOverride) { variation = *digiVariationOverride; }
        digitization_routines_map->at(sdname)->setDigitizationVariation(variation);
 
        // Resolve whether this system applies threshold / efficiency rejection, from the
        // global applyThresholds / applyInefficiencies options (default: neither).
        digitization_routines_map->at(sdname)->setHitRejectionPolicies(sdname);
 
        if (digitization_routines_map->at(sdname)->defineReadoutSpecs()) {
            log->info(1, "Digitization routine <%1> has been successfully defined.", sdname);
        } else { log->error(ERR_DEFINESPECFAIL, "defineReadoutSpecs failure for <%1>", sdname); }
    }
}
 
 
SystemList GDetectorConstruction::cloneSystemDescriptors(const SystemList& systems) const {
    SystemList descriptors;
    descriptors.reserve(systems.size());
 
    for (const auto& system: systems) {
        if (system != nullptr) {
            descriptors.emplace_back(system->descriptorClone(gopt));
        }
    }
 
    return descriptors;
}
 
 
void GDetectorConstruction::reload_geometry(SystemList sl) {
    // Geometry is changing: ensure Construct() reloads the digitization plugins.
    digiplugins_need_reload = true;
 
    // it could be empty for tests
    if (!sl.empty()) {
        gsystems = cloneSystemDescriptors(sl);
    }
 
    // Reconstruct the master geometry immediately so GUI pages can inspect the
    // new world without starting worker threads during a setup-tab reload.
    auto rm = G4RunManager::GetRunManager();
 
    if (rm) {
        // Null out the vis manager while we clean and rebuild the geometry stores.
        // G4*Store::Clean() would otherwise leave the ToolsSG scene graph with
        // dangling references to deleted Geant4 objects, causing a crash on the
        // next visualization flush.
        auto* visManager = G4VVisManager::GetConcreteInstance();
        GVisManagerGuard::set(nullptr);
        rm->DefineWorldVolume(Construct());
        GVisManagerGuard::set(visManager);
        ConstructSDandField();
    } else { log->error(1, "GDetectorConstruction::reload_geometryGeant4 Run manager not found."); }
}
 
void GDetectorConstruction::prepare_geometry_for_run() {
    auto rm = G4RunManager::GetRunManager();
 
    if (rm) {
        auto* visManager = G4VVisManager::GetConcreteInstance();
        GVisManagerGuard::set(nullptr);
        // Geometry is being rebuilt: ensure Construct() reloads the digitization plugins.
        digiplugins_need_reload = true;
        rm->ReinitializeGeometry(true, true);
        rm->GeometryHasBeenModified();
        // Optical processes such as G4Cerenkov cache material-property tables by
        // material index. A setup-tab reload can introduce a different variation's
        // optical material, so force physics tables to rebuild after the new
        // geometry/material set has been installed.
        rm->PhysicsHasBeenModified();
        rm->Initialize();
        GVisManagerGuard::set(visManager);
    } else { log->error(1, "GDetectorConstruction::prepare_geometry_for_runGeant4 Run manager not found."); }
}