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"
 
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();
}
 
void TRestGeant4AnalysisProcess::InitProcess() {
    fG4Metadata = GetMetadata<TRestGeant4Metadata>();
 
    std::vector<string> fObservables;
    fObservables = TRestEventProcess::ReadObservables();
 
    if (fPerProcessSensitiveEnergy) {
        fObservables.emplace_back("PerProcessPhotoelectric");
        fObservables.emplace_back("PerProcessCompton");
        fObservables.emplace_back("PerProcessElectronicIoni");
        fObservables.emplace_back("PerProcessAlphaIoni");
        fObservables.emplace_back("PerProcessIonIoni");
        fObservables.emplace_back("PerProcessHadronicIoni");
        fObservables.emplace_back("PerProcessProtonIoni");
        fObservables.emplace_back("PerProcessMsc");
        fObservables.emplace_back("PerProcessHadronElastic");
        fObservables.emplace_back("PerProcessNeutronElastic");
    }
    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();
 
            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) {
                cout << endl;
                cout << "??????????????????????????????????????????????????" << endl;
                cout << "REST warning : TRestGeant4AnalysisProcess." << endl;
                cout << "------------------------------------------" << endl;
                cout << endl;
                cout << " Volume " << volName << " 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;
            }
        }
 
        if (fObservables[i].find("MeanPos") != string::npos) {
            TString volName2 = fObservables[i].substr(0, fObservables[i].length() - 8).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) {
                cout << endl;
                cout << "??????????????????????????????????????????????????" << endl;
                cout << "REST warning : TRestGeant4AnalysisProcess." << endl;
                cout << "------------------------------------------" << endl;
                cout << endl;
                cout << " Volume " << volName2 << " 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;
            }
 
            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 (fObservables[i].find("Process") != string::npos) {
            Int_t ls = 0;
            if (fObservables[i].find("RadiactiveDecay") != string::npos) ls = 15;
            if (fObservables[i].find("Photoelectric") != string::npos) ls = 13;
            if (fObservables[i].find("PhotonNuclear") != string::npos) ls = 13;
            if (fObservables[i].find("Bremstralung") != string::npos) ls = 12;
            if (fObservables[i].find("NInelastic") != string::npos) ls = 10;
            if (fObservables[i].find("HadElastic") != string::npos) ls = 10;
            if (fObservables[i].find("NCapture") != string::npos) ls = 8;
            if (fObservables[i].find("Compton") != string::npos) ls = 7;
            if (fObservables[i].find("Neutron") != string::npos) ls = 7;
            if (fObservables[i].find("Alpha") != string::npos) ls = 5;
            if (fObservables[i].find("Argon") != string::npos) ls = 5;
            if (fObservables[i].find("Xenon") != string::npos) ls = 5;
            if (fObservables[i].find("Neon") != string::npos) ls = 4;
 
            TString processName = fObservables[i].substr(fObservables[i].length() - (ls + 7), ls).c_str();
            TString volName3 = fObservables[i].substr(0, fObservables[i].length() - (ls + 7)).c_str();
            Int_t volId3 = fG4Metadata->GetActiveVolumeID(volName3);
 
            if (volId3 >= 0) {
                fProcessObservables.push_back(fObservables[i]);
                fVolumeID3.push_back(volId3);
                fProcessName.emplace_back(processName.Data());
            }
        }
        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);
 
    // process names as named by Geant4
    // processes present here will be added to the list of observables which can be used to see if the event
    // contains the process of interest.
    vector<string> processNames = {"phot", "compt"};
    for (auto& processName : processNames) {
        Int_t containsProcess = 0;
        if (fOutputG4Event->ContainsProcess(fG4Metadata->GetGeant4PhysicsInfo().GetProcessID(processName))) {
            containsProcess = 1;
        }
 
        if (!processName.empty()) {
            processName[0] = toupper(processName[0]);
            SetObservableValue("containsProcess" + processName, containsProcess);
        }
    }
 
    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;
}
 
void TRestGeant4AnalysisProcess::EndProcess() {}