TRestGeant4AnalysisProcess.cxx

/*************************************************************************
 * This file is part of the REST software framework.                     *
 *                                                                       *
 * Copyright (C) 2016 GIFNA/TREX (University of Zaragoza)                *
 * For more information see http://gifna.unizar.es/trex                  *
 *                                                                       *
 * REST is free software: you can redistribute it and/or modify          *
 * it under the terms of the GNU General Public License as published by  *
 * the Free Software Foundation, either version 3 of the License, or     *
 * (at your option) any later version.                                   *
 *                                                                       *
 * REST is distributed in the hope that it will be useful,               *
 * but WITHOUT ANY WARRANTY; without even the implied warranty of        *
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the          *
 * GNU General Public License for more details.                          *
 *                                                                       *
 * You should have a copy of the GNU General Public License along with   *
 * REST in $REST_PATH/LICENSE.                                           *
 * If not, see http://www.gnu.org/licenses/.                             *
 * For the list of contributors see $REST_PATH/CREDITS.                  *
 *************************************************************************/
 
 
#include "TRestGeant4AnalysisProcess.h"
 
#include <algorithm>
#include <tuple>
 
using namespace std;
 
ClassImp(TRestGeant4AnalysisProcess);
 
TRestGeant4AnalysisProcess::TRestGeant4AnalysisProcess() { Initialize(); }
 
TRestGeant4AnalysisProcess::TRestGeant4AnalysisProcess(const char* configFilename) {
    Initialize();
 
    if (LoadConfigFromFile(configFilename)) LoadDefaultConfig();
}
 
TRestGeant4AnalysisProcess::~TRestGeant4AnalysisProcess() { delete fOutputG4Event; }
 
void TRestGeant4AnalysisProcess::LoadDefaultConfig() { SetTitle("Default config"); }
 
void TRestGeant4AnalysisProcess::Initialize() {
    fG4Metadata = nullptr;
    SetSectionName(this->ClassName());
    SetLibraryVersion(LIBRARY_VERSION);
 
    fInputG4Event = nullptr;
    fOutputG4Event = new TRestGeant4Event();
}
 
void TRestGeant4AnalysisProcess::LoadConfig(const string& configFilename, const string& name) {
    if (LoadConfigFromFile(configFilename, name)) LoadDefaultConfig();
}
 
namespace {
string CapitalizeFirst(string name) {
    if (!name.empty()) {
        name[0] = toupper(name[0]);
    }
    return name;
}
 
bool EndsWith(const string& text, const string& suffix) {
    return text.size() >= suffix.size() &&
           text.compare(text.size() - suffix.size(), suffix.size(), suffix) == 0;
}
 
bool StartsWith(const string& text, const string& prefix) { return text.rfind(prefix, 0) == 0; }
}  // namespace
 
const map<string, string>& TRestGeant4AnalysisProcess::GetRestToGeant4ProcessNameMap() {
    static const map<string, string> processNameMap = {
        // electromagnetic, photons
        {"Photoelectric", "phot"},
        {"Compton", "compt"},
        {"Rayleigh", "Rayl"},
        {"Conversion", "conv"},
        {"PhotonNuclear", "photonNuclear"},
        // electromagnetic, charged particles
        {"Bremsstrahlung", "eBrem"},
        {"Bremstralung", "eBrem"},  // legacy misspelling
        {"Annihilation", "annihil"},
        {"EIoni", "eIoni"},
        {"HIoni", "hIoni"},
        {"IonIoni", "ionIoni"},
        {"MuIoni", "muIoni"},
        {"Msc", "msc"},
        {"CoulombScat", "CoulombScat"},
        // hadronic
        {"HadElastic", "hadElastic"},
        {"NeutronInelastic", "neutronInelastic"},
        {"NInelastic", "neutronInelastic"},  // legacy name
        {"ProtonInelastic", "protonInelastic"},
        {"NCapture", "nCapture"},
        // decay
        {"RadioactiveDecay", "RadioactiveDecay"},
        {"RadiactiveDecay", "RadioactiveDecay"},  // legacy misspelling
        {"Decay", "Decay"},
    };
    return processNameMap;
}
 
string TRestGeant4AnalysisProcess::GetGeant4ProcessName(const string& name) const {
    if (name.empty()) {
        return "";
    }
 
    const auto& processNameMap = GetRestToGeant4ProcessNameMap();
    const auto it = processNameMap.find(name);
    if (it != processNameMap.end()) {
        return it->second;
    }
 
    for (const auto& [restName, geant4Name] : processNameMap) {
        if (CapitalizeFirst(geant4Name) == name) {
            return geant4Name;
        }
    }
 
    if (fG4Metadata != nullptr) {
        for (const auto& geant4Name : fG4Metadata->GetGeant4PhysicsInfo().GetAllProcesses()) {
            if (CapitalizeFirst(geant4Name.Data()) == name) {
                return geant4Name.Data();
            }
        }
    }
 
    return "";
}
 
void TRestGeant4AnalysisProcess::PrintNotActiveVolumeWarning(const TString& volumeName) {
    cout << endl;
    cout << "??????????????????????????????????????????????????" << endl;
    cout << "REST warning : TRestGeant4AnalysisProcess." << endl;
    cout << "------------------------------------------" << endl;
    cout << endl;
    cout << " Volume " << volumeName << " is not an active volume" << endl;
    cout << endl;
    cout << "List of active volumes : " << endl;
    cout << "------------------------ " << endl;
 
    for (unsigned int n = 0; n < fG4Metadata->GetNumberOfActiveVolumes(); n++)
        cout << "Volume " << n << " : " << fG4Metadata->GetActiveVolumeName(n) << endl;
    cout << "??????????????????????????????????????????????????" << endl;
    cout << endl;
}
 
void TRestGeant4AnalysisProcess::InitProcess() {
    fG4Metadata = GetMetadata<TRestGeant4Metadata>();
 
    // InitProcess may be called more than once
    fEnergyInObservables.clear();
    fVolumeID.clear();
    fVolumeName.clear();
    fMeanPosObservables.clear();
    fVolumeID2.clear();
    fDirID.clear();
    fProcessObservables.clear();
    fVolumeID3.clear();
    fProcessVolumeName.clear();
    fProcessName.clear();
    fTrackCounterObservables.clear();
    fParticleTrackCounter.clear();
    fTracksEDepObservables.clear();
    fParticleTrackEdep.clear();
    fContainsProcessObservables.clear();
    fContainsProcessG4Name.clear();
    fPerProcessObservables.clear();
    fPerProcessG4Names.clear();
    fPerProcessParticle.clear();
 
    std::vector<string> fObservables;
    fObservables = TRestEventProcess::ReadObservables();
    std::map<std::string, std::string> aliasObsToVol = GetAliasObservableNameToVolume();
 
    const auto& physicsInfo = fG4Metadata->GetGeant4PhysicsInfo();
    const auto availableProcesses = physicsInfo.GetAllProcesses();
 
    // containsProcessPhot and containsProcessCompt are always evaluated (historical behaviour)
    fContainsProcessObservables = {"containsProcessPhot", "containsProcessCompt"};
    fContainsProcessG4Name = {"phot", "compt"};
 
    if (fPerProcessSensitiveEnergy) {
        // {observable name, particle (empty for any), Geant4 processes}
        const vector<tuple<string, string, vector<string>>> perProcessDefinitions = {
            {"PerProcessPhotoelectric", "", {"phot"}},
            {"PerProcessCompton", "", {"compt"}},
            {"PerProcessElectronicIoni", "", {"eIoni"}},
            {"PerProcessAlphaIoni", "alpha", {"ionIoni", "alphaIoni"}},
            {"PerProcessIonIoni", "", {"ionIoni"}},
            {"PerProcessHadronicIoni", "", {"hIoni"}},
            {"PerProcessProtonIoni", "proton", {"hIoni"}},
            {"PerProcessMsc", "", {"msc"}},
            {"PerProcessHadronElastic", "", {"hadElastic"}},
            {"PerProcessNeutronElastic", "neutron", {"hadElastic"}},
        };
        for (const auto& [obsName, particle, processes] : perProcessDefinitions) {
            fPerProcessObservables.push_back(obsName);
            fPerProcessParticle.push_back(particle);
            fPerProcessG4Names.push_back(processes);
            if (fAnalysisTree != nullptr) {
                fAnalysisTree->AddObservable((string)GetName() + "_" + obsName, "double",
                                             "Energy deposited in the sensitive volume by a given process");
            }
        }
    }
 
    for (unsigned int i = 0; i < fObservables.size(); i++) {
        if (fObservables[i].find("VolumeEDep") != string::npos) {
            TString volName = fObservables[i].substr(0, fObservables[i].length() - 10).c_str();
            if (aliasObsToVol.find(fObservables[i]) != aliasObsToVol.end()) {
                volName = aliasObsToVol[fObservables[i]].c_str();
            }
 
            Int_t volId = fG4Metadata->GetActiveVolumeID(volName);
            if (volId >= 0) {
                fEnergyInObservables.push_back(fObservables[i]);
                fVolumeID.push_back(volId);
                fVolumeName.push_back(volName.Data());
            }
 
            if (volId == -1) {
                PrintNotActiveVolumeWarning(volName);
            }
        }
 
        if (fObservables[i].find("MeanPos") != string::npos) {
            TString volName2 = fObservables[i].substr(0, fObservables[i].length() - 8).c_str();
            if (aliasObsToVol.find(fObservables[i]) != aliasObsToVol.end()) {
                volName2 = aliasObsToVol[fObservables[i]].c_str();
            }
            std::string dirId = fObservables[i].substr(fObservables[i].length() - 1, 1);
 
            Int_t volId2 = fG4Metadata->GetActiveVolumeID(volName2);
            if (volId2 >= 0) {
                fMeanPosObservables.push_back(fObservables[i]);
                fVolumeID2.push_back(volId2);
                fDirID.push_back(dirId);
            }
 
            if (volId2 == -1) {
                PrintNotActiveVolumeWarning(volName2);
            }
 
            if ((dirId != "X") && (dirId != "Y") && (dirId != "Z")) {
                cout << endl;
                cout << "??????????????????????????????????????????????????" << endl;
                cout << "REST warning : TRestGeant4AnalysisProcess." << endl;
                cout << "------------------------------------------" << endl;
                cout << endl;
                cout << " Direction " << dirId << " is not valid" << endl;
                cout << " Only X, Y or Z accepted" << endl;
                cout << endl;
            }
        }
 
        if (StartsWith(fObservables[i], "containsProcess")) {
            const string processSuffix = fObservables[i].substr(string("containsProcess").length());
            if (std::find(fContainsProcessObservables.begin(), fContainsProcessObservables.end(),
                          fObservables[i]) != fContainsProcessObservables.end()) {
                continue;
            }
            const string geant4Name = GetGeant4ProcessName(processSuffix);
            if (geant4Name.empty()) {
                RESTWarning << "TRestGeant4AnalysisProcess: unknown process '" << processSuffix
                            << "' in observable '" << fObservables[i] << "'. It will not be filled."
                            << RESTendl;
                continue;
            }
            if (availableProcesses.count(geant4Name) == 0) {
                RESTWarning << "TRestGeant4AnalysisProcess: Geant4 process '" << geant4Name
                            << "' (observable '" << fObservables[i]
                            << "') is not registered in the physics info of this simulation. "
                            << "Its value will always be 0 (processes sharing type and subtype, e.g. "
                            << "eIoni and ionIoni, are registered only once)" << RESTendl;
            }
            fContainsProcessObservables.push_back(fObservables[i]);
            fContainsProcessG4Name.push_back(geant4Name);
            continue;
        }
 
        // <volume><Process>Process observables. Other observables ending in "Process" are excluded.
        if (EndsWith(fObservables[i], "Process") && !StartsWith(fObservables[i], "firstTrackInSensitive")) {
            const string nameWithoutSuffix =
                fObservables[i].substr(0, fObservables[i].length() - string("Process").length());
            const bool hasAlias = aliasObsToVol.find(fObservables[i]) != aliasObsToVol.end();
 
            // find the longest trailing substring which is a known process name
            string processName;
            string geant4Name;
            for (size_t length = 1; length <= nameWithoutSuffix.length(); length++) {
                // the volume name may only be empty if it is given through the `volume` attribute
                if (length == nameWithoutSuffix.length() && !hasAlias) {
                    break;
                }
                const string candidate = nameWithoutSuffix.substr(nameWithoutSuffix.length() - length);
                const string candidateGeant4Name = GetGeant4ProcessName(candidate);
                if (!candidateGeant4Name.empty()) {
                    processName = candidate;
                    geant4Name = candidateGeant4Name;
                }
            }
 
            if (geant4Name.empty()) {
                RESTWarning << "TRestGeant4AnalysisProcess: no known process name found in observable '"
                            << fObservables[i] << "'. It will not be filled." << RESTendl;
                RESTWarning << "Valid process names are:";
                for (const auto& [restName, g4Name] : GetRestToGeant4ProcessNameMap()) {
                    RESTWarning << " " << restName;
                }
                RESTWarning << " (or a Geant4 process name with the first letter capitalized)" << RESTendl;
                continue;
            }
 
            TString volName3 = nameWithoutSuffix.substr(0, nameWithoutSuffix.length() - processName.length());
            if (hasAlias) {
                volName3 = aliasObsToVol[fObservables[i]].c_str();
            }
            Int_t volId3 = fG4Metadata->GetActiveVolumeID(volName3);
 
            if (volId3 >= 0) {
                fProcessObservables.push_back(fObservables[i]);
                fVolumeID3.push_back(volId3);
                fProcessVolumeName.emplace_back(volName3.Data());
                fProcessName.push_back(geant4Name);
                if (availableProcesses.count(geant4Name) == 0) {
                    RESTWarning << "TRestGeant4AnalysisProcess: Geant4 process '" << geant4Name
                                << "' (observable '" << fObservables[i]
                                << "') is not registered in the physics info of this simulation. "
                                << "It will be 0 unless energy is deposited under this process name "
                                << "(processes sharing type and subtype, e.g. eIoni and ionIoni, are "
                                << "registered only once)" << RESTendl;
                }
            } else {
                PrintNotActiveVolumeWarning(volName3);
            }
        }
 
        if (fObservables[i].find("TracksCounter") != string::npos) {
            TString particleName = fObservables[i].substr(0, fObservables[i].length() - 13).c_str();
            fTrackCounterObservables.push_back(fObservables[i]);
            fParticleTrackCounter.emplace_back(particleName.Data());
        }
 
        if (fObservables[i].find("TracksEDep") != string::npos) {
            TString particleName = fObservables[i].substr(0, fObservables[i].length() - 10).c_str();
            fTracksEDepObservables.push_back(fObservables[i]);
            fParticleTrackEdep.emplace_back(particleName.Data());
        }
    }
}
 
TRestEvent* TRestGeant4AnalysisProcess::ProcessEvent(TRestEvent* inputEvent) {
    fInputG4Event = (TRestGeant4Event*)inputEvent;
    *fOutputG4Event = *((TRestGeant4Event*)inputEvent);
 
    const auto sensitiveVolumeName = fG4Metadata->GetSensitiveVolume();
 
    Double_t sensitiveVolumeEnergy = fOutputG4Event->GetEnergyInVolume(sensitiveVolumeName.Data());
    // Get time of the first hit in the sensitive volume
    double hitTime = std::numeric_limits<double>::infinity();
    std::set<int> trackIdsInSensitiveVolume;
 
    for (const auto& track : fOutputG4Event->GetTracks()) {
        const auto& hits = track.GetHits();
        for (size_t hitIndex = 0; hitIndex < hits.GetNumberOfHits(); hitIndex++) {
            const auto volumeName = hits.GetVolumeName(hitIndex);
            if (volumeName != sensitiveVolumeName) {
                continue;
            }
            const double energy = hits.GetEnergy(hitIndex);
            if (energy <= 0) {
                continue;
            }
 
            const double time = hits.GetTime(hitIndex);
            if (time < hitTime) {
                hitTime = time;
            }
 
            trackIdsInSensitiveVolume.insert(track.GetTrackID());
        }
    }
 
    SetObservableValue("sensitiveVolumeFirstHitTime", hitTime);
 
    std::set<int> trackParentsOfInterest;
    for (const auto& trackIdInSensitiveVolume : trackIdsInSensitiveVolume) {
        const auto track = fOutputG4Event->GetTrackByID(trackIdInSensitiveVolume);
        TRestGeant4Track* parent = track;
        bool found = false;
        while (parent != nullptr && !found) {
            // iterate over hits
            const auto& hits = parent->GetHits();
            if (hits.GetVolumeName(0) != sensitiveVolumeName) {
                for (size_t hitIndex = 0; hitIndex < hits.GetNumberOfHits(); hitIndex++) {
                    const auto volumeName = hits.GetVolumeName(hitIndex);
                    if (volumeName != sensitiveVolumeName) {
                        found = true;
                        trackParentsOfInterest.insert(parent->GetTrackID());
                        break;
                    }
                }
            }
 
            parent = parent->GetParentTrack();
        }
    }
 
    if (!trackParentsOfInterest.empty()) {
        const auto track = fOutputG4Event->GetTrackByID(*trackParentsOfInterest.begin());
 
        const auto position = track->GetInitialPosition();
        SetObservableValue("firstTrackInSensitivePositionX", position.X());
        SetObservableValue("firstTrackInSensitivePositionY", position.Y());
        SetObservableValue("firstTrackInSensitivePositionZ", position.Z());
 
        const auto energy = track->GetInitialKineticEnergy();
        SetObservableValue("firstTrackInSensitiveEnergy", energy);
 
        const string particleName = track->GetParticleName().Data();
        SetObservableValue("firstTrackInSensitiveParticle", particleName);
 
        string parentParticleName;
        if (track->GetParentTrack() != nullptr) {
            parentParticleName = track->GetParentTrack()->GetParticleName().Data();
        }
        SetObservableValue("firstTrackInSensitiveParentParticle", parentParticleName);
 
        const string creatorProcess = track->GetCreatorProcess().Data();
        SetObservableValue("firstTrackInSensitiveCreatorProcess", creatorProcess);
 
        const string volumeName = track->GetHits().GetVolumeName(0).Data();
        SetObservableValue("firstTrackInSensitiveVolumeName", volumeName);
    }
 
    SetObservableValue("firstTrackInSensitiveOk", trackParentsOfInterest.size() == 1);
 
    if (GetVerboseLevel() >= TRestStringOutput::REST_Verbose_Level::REST_Debug) {
        cout << "----------------------------" << endl;
        cout << "TRestGeant4Event : " << fOutputG4Event->GetID() << endl;
        cout << "Sensitive volume Energy : " << sensitiveVolumeEnergy << endl;
        cout << "Total energy : " << fOutputG4Event->GetTotalDepositedEnergy() << endl;
    }
 
    SetObservableValue("sensitiveVolumeEnergy", sensitiveVolumeEnergy);
 
    Double_t xOrigin = fOutputG4Event->GetPrimaryEventOrigin().X();
    SetObservableValue("xOriginPrimary", xOrigin);
 
    Double_t yOrigin = fOutputG4Event->GetPrimaryEventOrigin().Y();
    SetObservableValue("yOriginPrimary", yOrigin);
 
    Double_t zOrigin = fOutputG4Event->GetPrimaryEventOrigin().Z();
    SetObservableValue("zOriginPrimary", zOrigin);
 
    Double_t xDirection = fOutputG4Event->GetPrimaryEventDirection(0).X();
    SetObservableValue("xDirectionPrimary", xDirection);
 
    Double_t yDirection = fOutputG4Event->GetPrimaryEventDirection(0).Y();
    SetObservableValue("yDirectionPrimary", yDirection);
 
    Double_t zDirection = fOutputG4Event->GetPrimaryEventDirection(0).Z();
    SetObservableValue("zDirectionPrimary", zDirection);
 
    SetObservableValue("eventPrimaryParticleName", fOutputG4Event->GetPrimaryEventParticleName(0));
    string subEventPrimaryParticleName;
    if (fOutputG4Event->GetSubID() != 0) {
        subEventPrimaryParticleName = fOutputG4Event->GetSubEventPrimaryEventParticleName();
    } else {
        subEventPrimaryParticleName = "generator";
    }
    SetObservableValue("subEventPrimaryParticleName", subEventPrimaryParticleName);
 
    TVector3 primaryDirection(xDirection, yDirection, zDirection);
    primaryDirection = primaryDirection.Unit();
    SetObservableValue("thetaPrimary", primaryDirection.Theta());
    SetObservableValue("phiPrimary", primaryDirection.Phi());
    SetObservableValue("zenithYDegrees", TMath::ACos(-1 * yDirection) * TMath::RadToDeg());
    TVector3 sourceDirection = fG4Metadata->GetParticleSource(0)->GetDirection();
    sourceDirection = sourceDirection.Unit();
    SetObservableValue("zenithSourceDegrees",
                       TMath::ACos(primaryDirection.Dot(sourceDirection)) * TMath::RadToDeg());
 
    Double_t energyPrimary = fOutputG4Event->GetPrimaryEventEnergy(0);
    SetObservableValue("energyPrimary", energyPrimary);
 
    Double_t energyTotal = fOutputG4Event->GetTotalDepositedEnergy();
    SetObservableValue("totalEdep", energyTotal);
 
    Double_t size = fOutputG4Event->GetBoundingBoxSize();
    SetObservableValue("boundingSize", size);
 
    // containsProcessXxx observables: 1 if the event contains the Geant4 process, 0 otherwise
    const auto& physicsInfo = fG4Metadata->GetGeant4PhysicsInfo();
    for (unsigned int n = 0; n < fContainsProcessObservables.size(); n++) {
        const auto& processName = fContainsProcessG4Name[n];
        const Int_t processID = physicsInfo.GetProcessID(processName);
        // GetProcessID returns a default id (0) for unknown processes, which may be a valid id
        const bool processIsRegistered = physicsInfo.GetProcessName(processID) == processName;
        Int_t containsProcess = 0;
        if (processIsRegistered && fOutputG4Event->ContainsProcess(processID)) {
            containsProcess = 1;
        }
        SetObservableValue(fContainsProcessObservables[n], containsProcess);
    }
 
    if (!fProcessObservables.empty()) {
        const auto energyInVolumePerProcess = fOutputG4Event->GetEnergyInVolumePerProcessMap();
        for (unsigned int n = 0; n < fProcessObservables.size(); n++) {
            Double_t energy = 0;
            const auto volume = energyInVolumePerProcess.find(fProcessVolumeName[n]);
            if (volume != energyInVolumePerProcess.end()) {
                const auto process = volume->second.find(fProcessName[n]);
                if (process != volume->second.end()) {
                    energy = process->second;
                }
            }
            SetObservableValue(fProcessObservables[n], energy);
        }
    }
 
    if (!fPerProcessObservables.empty()) {
        const auto energyInVolumePerParticlePerProcess =
            fOutputG4Event->GetEnergyInVolumePerParticlePerProcessMap();
        const auto sensitiveVolume = energyInVolumePerParticlePerProcess.find(sensitiveVolumeName.Data());
        for (unsigned int n = 0; n < fPerProcessObservables.size(); n++) {
            Double_t energy = 0;
            if (sensitiveVolume != energyInVolumePerParticlePerProcess.end()) {
                for (const auto& [particle, energyPerProcess] : sensitiveVolume->second) {
                    if (!fPerProcessParticle[n].empty() && particle != fPerProcessParticle[n]) {
                        continue;
                    }
                    for (const auto& processName : fPerProcessG4Names[n]) {
                        const auto process = energyPerProcess.find(processName);
                        if (process != energyPerProcess.end()) {
                            energy += process->second;
                        }
                    }
                }
            }
            if (fPerProcessSensitiveEnergyNorm) {
                energy = sensitiveVolumeEnergy > 0 ? energy / sensitiveVolumeEnergy : 0;
            }
            SetObservableValue(fPerProcessObservables[n], energy);
        }
    }
 
    for (unsigned int n = 0; n < fParticleTrackCounter.size(); n++) {
        Int_t nT = fOutputG4Event->GetNumberOfTracksForParticle(fParticleTrackCounter[n]);
        string obsName = fTrackCounterObservables[n];
        SetObservableValue(obsName, nT);
    }
 
    for (unsigned int n = 0; n < fTracksEDepObservables.size(); n++) {
        Double_t energy = fOutputG4Event->GetEnergyDepositedByParticle(fParticleTrackEdep[n]);
        string obsName = fTracksEDepObservables[n];
        SetObservableValue(obsName, energy);
    }
 
    for (unsigned int n = 0; n < fEnergyInObservables.size(); n++) {
        Double_t en = fOutputG4Event->GetEnergyInVolume(fVolumeName[n]);
        string obsName = fEnergyInObservables[n];
        SetObservableValue(obsName, en);
    }
 
    for (unsigned int n = 0; n < fMeanPosObservables.size(); n++) {
        string obsName = fMeanPosObservables[n];
 
        Double_t mpos = 0;
        if (fDirID[n] == (TString) "X")
            mpos = fOutputG4Event->GetMeanPositionInVolume(fVolumeID2[n]).X();
 
        else if (fDirID[n] == (TString) "Y")
            mpos = fOutputG4Event->GetMeanPositionInVolume(fVolumeID2[n]).Y();
 
        else if (fDirID[n] == (TString) "Z")
            mpos = fOutputG4Event->GetMeanPositionInVolume(fVolumeID2[n]).Z();
 
        SetObservableValue(obsName, mpos);
    }
 
    if (GetVerboseLevel() >= TRestStringOutput::REST_Verbose_Level::REST_Debug) {
        cout << "G4 Tracks : " << fOutputG4Event->GetNumberOfTracks() << endl;
        cout << "----------------------------" << endl;
    }
 
    return fOutputG4Event;
}
 
std::map<std::string, std::string> TRestGeant4AnalysisProcess::GetAliasObservableNameToVolume() {
    TiXmlElement* e = GetElement("observable");
    std::map<std::string, std::string> obsNamesToVolumes;
 
    while (e != nullptr) {
        const char* obschr = e->Attribute("name");
        const char* _value = e->Attribute("value");
        const char* _vol = e->Attribute("volume");
 
        string value;
        if (_value == nullptr)
            value = "ON";
        else
            value = _value;
 
        if (ToUpper(value) == "ON") {
            if (obschr != nullptr) {
                string volume;
                if (_vol) {
                    volume = _vol;
                    obsNamesToVolumes[(string)obschr] = volume;
                }
            }
        }
 
        e = e->NextSiblingElement("observable");
    }
 
    return obsNamesToVolumes;
}
 
void TRestGeant4AnalysisProcess::EndProcess() {}