gDosimeterDigitization.cc

Go to the documentation of this file.

 
#include "gDosimeterDigitization.h"
#include <fstream>
 
// See header for API docs.
bool GDosimeterDigitization::defineReadoutSpecsImpl() {
    // Time window is the width of one electronics time cell (unit follows the convention
    // used by the rest of the digitization chain; typically ns).
    double timeWindow    = 10;                  // electronic readout time-window of the detector
    double gridStartTime = 0;                   // defines the windows grid
    auto   hitBitSet     = HitBitSet("000001"); // defines what information to be stored in the hit
    double maxStep       = 1 * CLHEP::mm;
 
    // Readout specs are immutable after initialization and shared by all processed hits.
    readoutSpecs = std::make_shared<GReadoutSpecs>(timeWindow, gridStartTime, hitBitSet, maxStep, log);
 
    return true;
}
 
// See header for API docs.
std::unique_ptr<GDigitizedData> GDosimeterDigitization::digitizeHitImpl(GHit* ghit, [[maybe_unused]] size_t hitn) {
    auto etot = ghit->getTotalEnergyDeposited();
    auto mass = ghit->getMass();
    auto dose = etot / mass;
 
    auto gdata = std::make_unique<GDigitizedData>(gopts, ghit);
 
    // Store the detector identity and the total deposited energy.
    gdata->includeVariable("etot", etot / MeV);
    gdata->includeVariable("dose", dose / gemc_units::picogray);
 
    //log->info(0, FUNCTION_NAME, "Edep: ", etot, ", hitn:", ghit->getEdeps().size());
 
    // // Per-step information used to build the NIEL-weight.
    // auto pids      = ghit->getPids();
    // auto pEnergies = ghit->getEs();
    //
    // double nielWeight = 0;
    //
    // // Accumulate NIEL factor step-by-step. This treats each step independently and sums
    // // the interpolated factors for supported particle species.
    // for (size_t stepIndex = 0; stepIndex < pids.size(); stepIndex++) {
    //  // Use absolute PID so negative particle IDs (e.g. -11) are handled.
    //  int pid = std::abs(pids[stepIndex]);
    //
    //  // Only a few particle types are supported by the calibration data files.
    //  if (pid == 11 || pid == 211 || pid == 2212 || pid == 2112) {
    //      // Convert from total energy to kinetic-like quantity used by the NIEL tables
    //      // by subtracting the particle rest mass (in MeV).
    //      double E   = pEnergies[stepIndex] - pMassMeV[pid];
    //      nielWeight += getNielFactorForParticleAtEnergy(pid, E);
    //  }
    // }
 
    // gdata->includeVariable("nielWeight", nielWeight);
 
    return gdata;
}
 
// See header for API docs.
bool GDosimeterDigitization::loadConstantsImpl([[maybe_unused]] int                runno,
                                               [[maybe_unused]] std::string const& variation) {
    // Map from particle id to the calibration filename (text, 2 columns).
    std::map<int, std::string> nielDataFiles;
    nielDataFiles[11]   = "niel_electron.txt";
    nielDataFiles[211]  = "niel_pion.txt";
    nielDataFiles[2112] = "niel_neutron.txt";
    nielDataFiles[2212] = "niel_proton.txt";
 
    // GEMC installation root used to locate plugin data.
    std::filesystem::path gemcRoot = gutilities::gemc_root();
 
    for (const auto& [pid, filename] : nielDataFiles) {
        std::string dataFileWithPath = gemcRoot.string() + "/dosimeterData" + "/Niel/" + filename;
 
        std::ifstream inputfile(dataFileWithPath);
        if (!inputfile) {
            // On Linux, tests may run from the build directory, where plugin data lives under
            // gdynamicDigitization/...
            dataFileWithPath = gemcRoot.string() + "/gdynamicDigitization/dosimeterData" + "/Niel/" + filename;
            inputfile.open(dataFileWithPath);
            if (!inputfile) {
                log->error(EC__FILENOTFOUND, "Error loading dosimeter data for pid <", pid, "> from file ",
                           dataFileWithPath);
            }
        }
 
        log->info(1, " Loading dosimeter data for pid <", pid, "> from file ", dataFileWithPath);
 
        // Expected file format: repeated pairs (factor, energyMeV).
        double p0, p1;
        while (inputfile >> p0 >> p1) {
            nielfactorMap[pid].push_back(p0);
            E_nielfactorMap[pid].push_back(p1);
        }
        inputfile.close();
    }
 
    // Particle rest masses used by the interpolation routine (MeV).
    pMassMeV[11]   = 0.510;
    pMassMeV[211]  = 139.570;
    pMassMeV[2112] = 939.565;
    pMassMeV[2212] = 938.272;
 
    return true;
}
 
// See header for API docs.
double GDosimeterDigitization::getNielFactorForParticleAtEnergy(int pid, double energyMeV) {
    const auto niel_N = nielfactorMap[pid].size();
 
    // Guard against missing/empty calibration data.
    // (If this happens, it indicates loadConstantsImpl did not populate the maps as expected.)
    if (niel_N == 0) {
        log->error(EC__FILENOTFOUND, "NIEL tables are empty for pid <", pid, ">. Did loadConstantsImpl fail?");
        return 0.0;
    }
 
    // j is the first index for which energyMeV < E_table[j]. Initialize to "past the end".
    size_t j = niel_N;
 
    for (size_t i = 0; i < niel_N; i++) {
        if (energyMeV < E_nielfactorMap[pid][i]) {
            j = i;
            break;
        }
    }
 
    double value = 0.0;
 
    if (j > 0 && j < niel_N) {
        // Linear interpolation between (j-1) and j.
        const auto nielfactor_low  = nielfactorMap[pid][j - 1];
        const auto nielfactor_high = nielfactorMap[pid][j];
        const auto energy_low      = E_nielfactorMap[pid][j - 1];
        const auto energy_high     = E_nielfactorMap[pid][j];
 
        value = nielfactor_low +
                ( nielfactor_high - nielfactor_low ) / ( energy_high - energy_low ) * ( energyMeV - energy_low );
    }
    else if (j == 0) {
        // energy below the first threshold: clamp to first value.
        value = nielfactorMap[pid].front();
    }
    else {
        // energy beyond the last threshold: clamp to last value.
        value = nielfactorMap[pid].back();
    }
 
    log->debug(NORMAL, " pid: ", pid, ", j: ", j, ", value: ", value, ", energy: ", energyMeV);
 
    return value;
}