TRestGeant4ToDetectorHitsProcess.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 "TRestGeant4ToDetectorHitsProcess.h"
 
using namespace std;
 
ClassImp(TRestGeant4ToDetectorHitsProcess);
 
TRestGeant4ToDetectorHitsProcess::TRestGeant4ToDetectorHitsProcess() { Initialize(); }
 
TRestGeant4ToDetectorHitsProcess::TRestGeant4ToDetectorHitsProcess(const char* configFilename) {
    Initialize();
 
    if (LoadConfigFromFile(configFilename)) {
        LoadDefaultConfig();
    }
}
 
TRestGeant4ToDetectorHitsProcess::~TRestGeant4ToDetectorHitsProcess() { delete fHitsEvent; }
 
void TRestGeant4ToDetectorHitsProcess::LoadDefaultConfig() {
    SetTitle("Default config");
 
    cout << "Geant4 to hits metadata not found. Loading default values" << endl;
}
 
void TRestGeant4ToDetectorHitsProcess::Initialize() {
    SetSectionName(this->ClassName());
    SetLibraryVersion(LIBRARY_VERSION);
 
    fGeant4Event = nullptr;
    fHitsEvent = new TRestDetectorHitsEvent();
}
 
void TRestGeant4ToDetectorHitsProcess::LoadConfig(const string& configFilename, const string& name) {
    if (LoadConfigFromFile(configFilename, name)) LoadDefaultConfig();
}
 
void TRestGeant4ToDetectorHitsProcess::InitProcess() {
    fGeant4Metadata = GetMetadata<TRestGeant4Metadata>();
 
    for (const auto& userVolume : fVolumeSelection) {
        if (fGeant4Metadata->GetActiveVolumeID(userVolume) >= 0) {
            fVolumeId.push_back(fGeant4Metadata->GetActiveVolumeID(userVolume));
        } else if (GetVerboseLevel() >= TRestStringOutput::REST_Verbose_Level::REST_Warning)
            cout << "TRestGeant4ToDetectorHitsProcess. volume name : " << userVolume
                 << " not found and will not be added." << endl;
    }
 
    sort(fVolumeId.begin(), fVolumeId.end());
    fVolumeId.erase(unique(fVolumeId.begin(), fVolumeId.end()), fVolumeId.end());
 
    for (size_t i = 0; i < fVolumeId.size(); i++) {
        RESTDebug << "TRestGeant4ToDetectorHitsProcess. Volume id : " << fVolumeId[i]
                  << " name : " << fGeant4Metadata->GetActiveVolumeName(fVolumeId[i]) << RESTendl;
    }
 
    RESTDebug << "Active volumes available in TRestGeant4Metadata" << RESTendl;
    RESTDebug << "-------------------------------------------" << RESTendl;
    for (unsigned int n = 0; n < fGeant4Metadata->GetNumberOfActiveVolumes(); n++) {
        RESTDebug << "Volume id : " << n << " name : " << fGeant4Metadata->GetActiveVolumeName(n) << RESTendl;
    }
    RESTDebug << RESTendl;
 
    RESTDebug << "TRestGeant4HitsProcess volumes enabled in RML : ";
    RESTDebug << "-------------------------------------------" << RESTendl;
    if (fVolumeSelection.empty())
        RESTDebug << "all" << RESTendl;
    else {
        for (const auto& volume : fVolumeSelection) {
            RESTDebug << "" << RESTendl;
            RESTDebug << " - " << volume << RESTendl;
        }
        RESTDebug << " " << RESTendl;
    }
 
    if (!fVolumeSelection.empty() && fVolumeSelection.size() != fVolumeId.size())
        RESTWarning << "TRestGeant4ToDetectorHitsProcess. Not all volumes were properly identified!"
                    << RESTendl;
 
    if (!fVolumeId.empty()) {
        RESTDebug << "TRestGeant4HitsProcess volumes identified : ";
        RESTDebug << "---------------------------------------" << RESTendl;
        if (fVolumeSelection.empty())
            RESTDebug << "all" << RESTendl;
        else
            for (const auto& volume : fVolumeSelection) {
                RESTDebug << "" << RESTendl;
                RESTDebug << " - " << volume << RESTendl;
            }
        RESTDebug << " " << RESTendl;
    }
}
 
TRestEvent* TRestGeant4ToDetectorHitsProcess::ProcessEvent(TRestEvent* inputEvent) {
    fGeant4Event = (TRestGeant4Event*)inputEvent;
 
    fGeant4Event->InitializeReferences(GetRunInfo());
 
    if (GetVerboseLevel() >= TRestStringOutput::REST_Verbose_Level::REST_Extreme) {
        cout << "------ TRestGeant4ToDetectorHitsProcess --- Printing Input Event --- START ----" << endl;
        fGeant4Event->PrintEvent();
        cout << "------ TRestGeant4ToDetectorHitsProcess --- Printing Input Event ---- END ----" << endl;
        GetChar();
    }
 
    fHitsEvent->SetRunOrigin(fGeant4Event->GetRunOrigin());
    fHitsEvent->SetSubRunOrigin(fGeant4Event->GetSubRunOrigin());
    fHitsEvent->SetID(fGeant4Event->GetID());
    fHitsEvent->SetSubID(fGeant4Event->GetSubID());
    fHitsEvent->SetSubEventTag(fGeant4Event->GetSubEventTag());
    fHitsEvent->SetTimeStamp(fGeant4Event->GetTimeStamp());
    fHitsEvent->SetState(fGeant4Event->isOk());
 
    for (const auto& track : fGeant4Event->GetTracks()) {
        const auto& hits = track.GetHits();
        for (unsigned int i = 0; i < track.GetNumberOfHits(); i++) {
            const auto energy = hits.GetEnergy(i);
            if (energy <= 0) {
                continue;
            }
            const TVector3& position = hits.GetPosition(i);
            const double time = hits.GetTime(i);
 
            if (time < 0) {
                cerr << "TRestGeant4ToDetectorHitsProcess. Negative time found. This should never happen"
                     << endl;
                exit(1);
            }
            if (fVolumeId.empty()) {
                // if no volume is selected, all hits are added
                fHitsEvent->AddHit(position.X(), position.Y(), position.Z(), energy, time);
            } else {
                const string volumeName = hits.GetVolumeName(i).Data();
                const auto volumeId = fGeant4Metadata->GetActiveVolumeID(volumeName);
 
                // cout << "volumeName : " << volumeName << " volumeId : " << volumeId << endl;
                if (find(fVolumeId.begin(), fVolumeId.end(), volumeId) != fVolumeId.end()) {
                    const REST_HitType type = fHitTypes.at(volumeName);
 
                    fHitsEvent->AddHit(position, energy, time, type);
                }
            }
        }
    }
 
    if (GetVerboseLevel() >= TRestStringOutput::REST_Verbose_Level::REST_Debug) {
        cout << "TRestGeant4ToDetectorHitsProcess. Hits added : " << fHitsEvent->GetNumberOfHits() << endl;
        cout << "TRestGeant4ToDetectorHitsProcess. Hits total energy : " << fHitsEvent->GetTotalEnergy()
             << endl;
    }
 
    return fHitsEvent;
}
 
void TRestGeant4ToDetectorHitsProcess::InitFromConfigFile() {
    // Attempt to access TRestGeant4Metadata
    fGeant4Metadata = GetMetadata<TRestGeant4Metadata>();
    if (fGeant4Metadata == nullptr) {
        RESTWarning << "TRestGeant4ToDetectorHitsProcess. No TRestGeant4Metadata found in the input file"
                    << RESTendl;
    }
 
    set<string> volumesToAdd;
    TiXmlElement* volumeDefinition = GetElement("volume");
    if (volumeDefinition == nullptr) {
        volumeDefinition = GetElement("addVolume");
        if (volumeDefinition != nullptr) {
            RESTWarning << "TRestGeant4ToDetectorHitsProcess. 'addVolume' tag is deprecated. Please use "
                           "'volume' instead."
                        << RESTendl;
        }
    }
    while (volumeDefinition != nullptr) {
        const auto userVolume = GetFieldValue("name", volumeDefinition);
        const auto typeName = GetFieldValue("type", volumeDefinition);
        REST_HitType type = XYZ;
        if (typeName == "veto") {
            type = VETO;
        }
 
        if (userVolume == "Not defined") {
            RESTError << "TRestGeant4ToDetectorHitsProcess. No name defined for volume" << RESTendl;
        }
        if (fGeant4Metadata != nullptr) {
            const auto& geometryInfo = fGeant4Metadata->GetGeant4GeometryInfo();
 
            auto physicalVolumes = geometryInfo.GetAllPhysicalVolumesMatchingExpression(userVolume);
            if (physicalVolumes.empty()) {
                const auto logicalVolumes = geometryInfo.GetAllLogicalVolumesMatchingExpression(userVolume);
                for (const auto& logicalVolume : logicalVolumes) {
                    for (const auto& physicalVolume :
                         geometryInfo.GetAllPhysicalVolumesFromLogical(logicalVolume)) {
                        physicalVolumes.push_back(
                            geometryInfo.GetAlternativeNameFromGeant4PhysicalName(physicalVolume));
                    }
                }
            }
            for (const auto& physicalVolume : physicalVolumes) {
                volumesToAdd.insert(physicalVolume.Data());
                fHitTypes[physicalVolume.Data()] = type;
            }
        } else {
            volumesToAdd.insert(userVolume);
            fHitTypes[userVolume] = type;
        }
 
        volumeDefinition = GetNextElement(volumeDefinition);
    }
 
    for (const auto& volume : volumesToAdd) {
        if (find(fVolumeSelection.begin(), fVolumeSelection.end(), volume) == fVolumeSelection.end()) {
            fVolumeSelection.emplace_back(volume);
        }
    }
}
 
void TRestGeant4ToDetectorHitsProcess::PrintMetadata() {
    BeginPrintProcess();
 
    for (const auto& volume : fVolumeSelection) {
        RESTMetadata << "Volume added : " << volume << RESTendl;
    }
 
    EndPrintProcess();
}