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 "gDosimeterDigitization.h"
#include "gFluxDigitization.h"
#include "gParticleCounterDigitization.h"
 
// geant4
#include "G4SDManager.hh"
#include "G4GeometryManager.hh"
#include "G4SolidStore.hh"
#include "G4LogicalVolumeStore.hh"
#include "G4PhysicalVolumeStore.hh"
#include "G4ReflectionFactory.hh"
#include "G4RunManager.hh"
#include "G4UserLimits.hh"
 
G4ThreadLocal GMagneto *GDetectorConstruction::gmagneto = nullptr;
 
GDetectorConstruction::GDetectorConstruction(std::shared_ptr<GOptions> gopts)
    : GBase(gopts, GDETECTOR_LOGGER),
      G4VUserDetectorConstruction(), // Geant4 base class.
      gopt(gopts) {
    // Map is populated after SDs exist, in the SD/field construction path.
    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, FUNCTION_NAME);
 
    // Clean any old geometry.
    // This is required when reloading geometry, to prevent stale stores and duplicated objects.
    G4GeometryManager::GetInstance()->OpenGeometry();
    G4PhysicalVolumeStore::Clean();
    G4LogicalVolumeStore::Clean();
    G4SolidStore::Clean();
    G4ReflectionFactory::Instance()->Clean();
 
    // 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, FUNCTION_NAME, "creating world from options");
        gworld = std::make_shared<GWorld>(gopt);
    } else {
        log->debug(NORMAL, FUNCTION_NAME, "creating world from a gsystem vector of size ", gsystems.size());
        gworld = std::make_shared<GWorld>(gopt, 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 - ", gworld->get_number_of_volumes() - 1, " volumes\n - ",
              g4world->number_of_volumes(), " geant4 built volumes\n - ",
              nsdetectors, " sensitive detectors\n");
 
 
    // Return the physical volume for the ROOT world volume.
    return g4world->getG4Volume(ROOTWORLDGVOLUMENAME)->getPhysical();
}
 
// Installs sensitive detectors and EM fields for the constructed geometry.
void GDetectorConstruction::ConstructSDandField() {
    auto sdManager = G4SDManager::GetSDMpointer();
 
    log->debug(NORMAL, FUNCTION_NAME);
 
    // 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;
 
    // 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();
            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) {
                std::string copyOf = gvolumePtr->getCopyOf();
                if (copyOf != "" && copyOf != UNINITIALIZEDSTRINGQUANTITY) {
                    auto gsystem = gvolumePtr->getSystem();
                    auto volume_copy = gsystem + "/" + copyOf;
                    auto copyG4Volume = g4world->getG4Volume(volume_copy)->getLogical();
                    if (copyG4Volume != nullptr) { g4volume = copyG4Volume; } else {
                        log->error(ERR_GVOLUMENOTFOUND, FUNCTION_NAME,
                                   " Logical volume copy <" + volume_copy + "> not found.");
                    }
                }
            }
            if (g4volume == nullptr) {
                log->error(ERR_GVOLUMENOTFOUND, FUNCTION_NAME, " Logical volume <" + g4name + "> not found.");
            }
 
            // Skip volumes with no digitization.
            if (digitizationName != "" && digitizationName != UNINITIALIZEDSTRINGQUANTITY) {
                // Create the sensitive detector if it does not exist yet.
                if (sensitiveDetectorsMap.find(digitizationName) == sensitiveDetectorsMap.end()) {
                    log->info(2, "Creating new sensitive detector <", digitizationName, "> for volume <", g4name, ">");
 
                    sensitiveDetectorsMap[digitizationName] = new GSensitiveDetector(digitizationName, gopt);
                } else {
                    log->info(2, "Sensitive detector <", digitizationName,
                              "> is already created and available for volume <", 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();
                const auto &mass = g4volume->GetMass();
                auto this_gtouchable = std::make_shared<
                    GTouchable>(gopt, digitizationName, identity, vdimensions, mass);
                sensitiveDetectorsMap[digitizationName]->registerGVolumeTouchable(g4name, this_gtouchable);
 
                // Register the detector with Geant4 and attach it to the logical volume.
                sdManager->AddNewDetector(sensitiveDetectorsMap[digitizationName]);
                g4volume->SetSensitiveDetector(sensitiveDetectorsMap[digitizationName]);
 
                //  auto maxStep =
 
                //g4volume->SetUserLimits(new G4UserLimits(0.1*mm, 0.1*mm));
 
                log->info(2, "Logical Volume  <" + g4name + "> has been successfully assigned to SD.",
                          sensitiveDetectorsMap[digitizationName]);
            }
 
            // Process electromagnetic fields.
            // If a volume declares an EM field, ensure the field container exists and install
            // a per-volume field manager configured by the named field map.
            const auto &field_name = gvolumePtr->getEMField();
            if (field_name != "" && field_name != UNINITIALIZEDSTRINGQUANTITY) {
                if (gmagneto == nullptr) { gmagneto = new GMagneto(gopt); }
                log->info(2, "Volume <", volumeName, "> has field: <", field_name,
                          ">. Looking into field map definitions.");
                log->info(2, "Setting field manager for volume <", g4name, "> with field <", field_name, ">");
                g4world->setFieldManagerForVolume(g4name, gmagneto->getFieldMgr(field_name).get(), true);
            }
        }
    }
 
    // Load digitization plugins after constructing sensitive detectors.
    // This populates digitization_routines_map for each sensitive detector name.
    loadDigitizationPlugins();
 
    // 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);
        double maxStep = digitization_routine->readoutSpecs->getMaxStep();
 
        sensitiveDetectorsMap[sdname]->assign_digi_routine(digitization_routine);
        log->info(1, "Digitization routine <" + sdname + "> has been successfully assigned to SD.",
                  sensitiveDetectorsMap[sdname]);
 
        // Loop over all systems and their volumes.
        // and assign max step to the corresponding logical volume
        for (const auto &[systemName, gsystemPtr]: *gworld->getSystemsMap()) {
            for (const auto &[volumeName, gvolumePtr]: gsystemPtr->getGVolumesMap()) {
                auto const &digitizationName = gvolumePtr->getDigitization();
                if (digitizationName == sdname) {
                    auto const &g4name = gvolumePtr->getG4Name();
                    auto *g4volume = g4world->getG4Volume(g4name)->getLogical();
 
                    // g4volume->SetUserLimits(new G4UserLimits(maxStep, maxStep)); // this will also kill track cause
                    // the second argument is max track length
                    g4volume->SetUserLimits(new G4UserLimits(maxStep));
 
                    log->info(1, "Setting G4UserLimits for volume <", g4name, "> with maxStep <", maxStep, ">");
                }
            }
        }
    }
}
 
// Loads (or dynamically resolves) the digitization routine for each sensitive detector.
void GDetectorConstruction::loadDigitizationPlugins() {
    const auto sdetectors = gworld->getSensitiveDetectorsList();
 
    for (auto &sdname: sdetectors) {
        if (sdname == FLUXNAME) {
            log->info(1, "Loading flux digitization plugin for routine <" + sdname + ">");
            digitization_routines_map->emplace(sdname, std::make_shared<GFluxDigitization>(gopt));
        } else if (sdname == COUNTERNAME) {
            log->info(1, "Loading particle counter digitization plugin for routine <" + sdname + ">");
            digitization_routines_map->emplace(sdname, std::make_shared<GParticleCounterDigitization>(gopt));
        } else if (sdname == DOSIMETERNAME) {
            log->info(1, "Loading dosimeter digitization plugin for routine <" + 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 <" + 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);
 
        if (digitization_routines_map->at(sdname)->defineReadoutSpecs()) {
            log->info(1, "Digitization routine <" + sdname + "> has been successfully defined.");
        } else { log->error(ERR_DEFINESPECFAIL, "defineReadoutSpecs failure for <" + sdname + ">"); }
    }
}
 
 
void GDetectorConstruction::reload_geometry(SystemList sl) {
    // it could be empty for tests
    if (!sl.empty()) {
        // Use vector assignment to update the local systems.
        gsystems = sl;
    }
 
    // Reconstruct the geometry and update the world volume - if the run manager exists
    auto rm = G4RunManager::GetRunManager();
 
    // TODO: not sure if DefineWorldVolume also call ConstructSDandField automatically?
    // is this all there is to do here to reload a geometry?
    if (rm) {
        rm->DefineWorldVolume(Construct());
        ConstructSDandField();
    } else { log->error(1, "GDetectorConstruction::reload_geometry", "Geant4 Run manager not found."); }
}