TRestGeant4VetoAnalysisProcess.cxx

/*************************************************************************
 * This file is part of the REST software framework.                     *
 *                                                                       *
 * Copyright (C) 2020 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 "TRestGeant4VetoAnalysisProcess.h"
 
using namespace std;
 
ClassImp(TRestGeant4VetoAnalysisProcess);
 
TRestGeant4VetoAnalysisProcess::TRestGeant4VetoAnalysisProcess() { Initialize(); }
 
TRestGeant4VetoAnalysisProcess::TRestGeant4VetoAnalysisProcess(const char* configFilename) {
    TRestGeant4VetoAnalysisProcess();
    if (LoadConfigFromFile(configFilename)) {
        LoadDefaultConfig();
    }
}
 
TRestGeant4VetoAnalysisProcess::~TRestGeant4VetoAnalysisProcess() { delete fOutputEvent; }
 
void TRestGeant4VetoAnalysisProcess::LoadDefaultConfig() { SetTitle("Default config"); }
 
void TRestGeant4VetoAnalysisProcess::Initialize() {
    SetSectionName(this->ClassName());
    SetLibraryVersion(LIBRARY_VERSION);
 
    fOutputEvent = new TRestGeant4Event();
}
 
void TRestGeant4VetoAnalysisProcess::LoadConfig(const string& configFilename, const string& name) {
    if (LoadConfigFromFile(configFilename, name)) LoadDefaultConfig();
}
 
void TRestGeant4VetoAnalysisProcess::InitProcess() {
    // CAREFUL THIS METHOD IS CALLED TWICE!
    fVetoVolumes.clear();
    fVetoDetectorVolumes.clear();
    fVetoDetectorBoundaryDirection.clear();
    fVetoDetectorBoundaryPosition.clear();
 
    if (fGeant4Metadata == nullptr) {
        // maybe it was manually initialized
        fGeant4Metadata = GetMetadata<TRestGeant4Metadata>();
    }
    if (fGeant4Metadata == nullptr) {
        cerr << "TRestGeant4VetoAnalysisProcess::InitProcess: Geant4 metadata not found" << endl;
        exit(1);
    }
 
    const auto& geometryInfo = fGeant4Metadata->GetGeant4GeometryInfo();
 
    fVetoVolumes = geometryInfo.GetAllPhysicalVolumesMatchingExpression(fVetoVolumesExpression);
    if (fVetoVolumes.empty()) {
        const auto logicalVolumes =
            geometryInfo.GetAllLogicalVolumesMatchingExpression(fVetoVolumesExpression);
        for (const auto& logicalVolume : logicalVolumes) {
            for (const auto& physicalVolume : geometryInfo.GetAllPhysicalVolumesFromLogical(logicalVolume)) {
                fVetoVolumes.push_back(geometryInfo.GetAlternativeNameFromGeant4PhysicalName(physicalVolume));
            }
        }
    }
    if (fVetoVolumes.empty()) {
        cerr << "TRestGeant4VetoAnalysisProcess::InitProcess: No veto volumes found" << endl;
        exit(1);
    }
 
    // get detector volumes if requested
    if (!fVetoDetectorsExpression.IsNull()) {
        fVetoDetectorVolumes = geometryInfo.GetAllPhysicalVolumesMatchingExpression(fVetoDetectorsExpression);
        if (fVetoDetectorVolumes.empty()) {
            const auto logicalVolumes =
                geometryInfo.GetAllLogicalVolumesMatchingExpression(fVetoDetectorsExpression);
            for (const auto& logicalVolume : logicalVolumes) {
                for (const auto& physicalVolume :
                     geometryInfo.GetAllPhysicalVolumesFromLogical(logicalVolume)) {
                    fVetoDetectorVolumes.push_back(
                        geometryInfo.GetAlternativeNameFromGeant4PhysicalName(physicalVolume));
                }
            }
        }
        if (fVetoDetectorVolumes.empty()) {
            cerr << "TRestGeant4VetoAnalysisProcess::InitProcess: No detector volumes found" << endl;
            exit(1);
        }
        if (fVetoDetectorVolumes.size() != fVetoVolumes.size()) {
            cerr << "TRestGeant4VetoAnalysisProcess::InitProcess: Number of detector volumes "
                 << "does not match number of veto volumes" << endl;
            exit(1);
        }
    }
 
    for (unsigned int i = 0; i < fVetoDetectorVolumes.size(); i++) {
        const auto& vetoName = fVetoVolumes[i];
        const auto& vetoPosition = geometryInfo.GetPosition(vetoName);
 
        const auto& vetoDetectorName = fVetoDetectorVolumes[i];
        const auto& vetoDetectorPosition = geometryInfo.GetPosition(vetoDetectorName);
 
        const auto distance = vetoDetectorPosition - vetoPosition;
        const auto direction = distance.Unit();
 
        fVetoDetectorBoundaryDirection[vetoName] = direction;
        fVetoDetectorBoundaryPosition[vetoName] = vetoDetectorPosition - direction * fVetoDetectorOffsetSize;
    }
 
    PrintMetadata();
}
 
TRestEvent* TRestGeant4VetoAnalysisProcess::ProcessEvent(TRestEvent* inputEvent) {
    /*
    fInputEvent = (TRestGeant4Event*)inputEvent;
    *fOutputEvent = *((TRestGeant4Event*)inputEvent);
 
    std::map<string, Double_t> volume_energy_map;
 
    for (unsigned int i = 0; i < fVetoVolumeIds.size(); i++) {
        int id = fVetoVolumeIds[i];
        string volume_name = (string)fGeant4Metadata->GetActiveVolumeName(id);
 
        Double_t energy = fOutputG4Event->GetEnergyDepositedInVolume(id);
        volume_energy_map[volume_name] = energy;
    }
 
    Double_t energy_veto_max = 0;
    for (const auto& pair : volume_energy_map) {
        Double_t veto_energy = pair.second;
        SetObservableValue(pair.first + "VolumeEDep", veto_energy);
        if (veto_energy > energy_veto_max) {
            energy_veto_max = veto_energy;
        };
    }
    SetObservableValue("vetoAllEVetoMax", energy_veto_max);
 
    // veto groups
    for (const auto& pair : fVetoGroupVolumeNames) {
        Double_t energy_veto_max_group = 0;
        for (unsigned int i = 0; i < pair.second.size(); i++) {
            string volume_name = pair.second[i];
            Double_t veto_energy = volume_energy_map[volume_name];
            if (veto_energy > energy_veto_max_group) {
                energy_veto_max_group = veto_energy;
            };
        }
        // convert to Upper + lower case (VetoGroupTopEVetoMax, ...)
        string group_name;
        for (auto it = pair.first.cbegin(); it != pair.first.cend(); ++it) {
            if (it == pair.first.cbegin()) {
                group_name += std::toupper(*it);
            } else {
                group_name += std::tolower(*it);
            }
        }
        SetObservableValue("vetoGroup" + group_name + "EVetoMax", energy_veto_max_group);
    }
 
    //
 
    for (const auto& quenching_factor : fQuenchingFactors) {
        string quenching_factor_string = std::to_string(quenching_factor);
        // replace "." in string by "_" because its gives very strange problems
        quenching_factor_string =
            quenching_factor_string.replace(quenching_factor_string.find("."), sizeof(".") - 1, "_");
        volume_energy_map.clear();
        for (int i = 0; i < fOutputG4Event->GetNumberOfTracks(); i++) {
            auto track = fOutputG4Event->GetTrack(i);
            string particle_name = (string)track.GetParticleName();
            for (const auto& id : fVetoVolumeIds) {
                string volume_name = (string)fGeant4Metadata->GetActiveVolumeName(id);
 
                if (particle_name == "e-" || particle_name == "e+" || particle_name == "gamma") {
                    // no quenching factor
                    volume_energy_map[volume_name] += track.GetEnergyInVolume(id);
                } else {
                    // apply quenching factor
                    volume_energy_map[volume_name] += quenching_factor * track.GetEnergyInVolume(id);
                }
            }
        }
 
        Double_t energy_veto_max = 0;
        for (const auto& pair : volume_energy_map) {
            Double_t veto_energy = pair.second;
            SetObservableValue(pair.first + "VolumeEDep" + "Qf" + quenching_factor_string, veto_energy);
            if (veto_energy > energy_veto_max) {
                energy_veto_max = veto_energy;
            };
        }
        SetObservableValue(string("vetoAllEVetoMax") + "Qf" + quenching_factor_string, energy_veto_max);
 
        // veto groups
        for (const auto& pair : fVetoGroupVolumeNames) {
            Double_t energy_veto_max_group = 0;
            for (unsigned int i = 0; i < pair.second.size(); i++) {
                string volume_name = pair.second[i];
                Double_t veto_energy = volume_energy_map[volume_name];
                if (veto_energy > energy_veto_max_group) {
                    energy_veto_max_group = veto_energy;
                };
            }
            // convert to Upper + lower case (VetoGroupTopEVetoMax, ...)
            string group_name;
            for (auto it = pair.first.cbegin(); it != pair.first.cend(); ++it) {
                if (it == pair.first.cbegin()) {
                    group_name += std::toupper(*it);
                } else {
                    group_name += std::tolower(*it);
                }
            }
            SetObservableValue("vetoGroup" + group_name + "EVetoMax" + "Qf" + quenching_factor_string,
                               energy_veto_max_group);
        }
    }
 
    return fOutputG4Event;
     */
    return fInputEvent;
}
 
void TRestGeant4VetoAnalysisProcess::EndProcess() {
    // Function to be executed once at the end of the process
    // (after all events have been processed)
 
    // Start by calling the EndProcess function of the abstract class.
    // Comment this if you don't want it.
    // TRestEventProcess::EndProcess();
}
 
void TRestGeant4VetoAnalysisProcess::InitFromConfigFile() {
    // word to identify active volume as veto (default = "veto" e.g. "vetoTop")
    fVetoVolumesExpression = GetParameter("vetoVolumesExpression", fVetoVolumesExpression);
    fVetoDetectorsExpression = GetParameter("vetoDetectorsExpression", fVetoDetectorsExpression);
 
    fVetoDetectorOffsetSize = GetDblParameterWithUnits("vetoDetectorOffset", fVetoDetectorOffsetSize);
    fVetoLightAttenuation = GetDblParameterWithUnits("vetoLightAttenuation", fVetoLightAttenuation);
    fVetoQuenchingFactor = GetDblParameterWithUnits("quenchingFactor", fVetoQuenchingFactor);
}
 
void TRestGeant4VetoAnalysisProcess::PrintMetadata() {
    BeginPrintProcess();
 
    cout << "Veto volume expression: " << fVetoVolumesExpression << endl;
    if (!fVetoDetectorsExpression.IsNull()) {
        cout << "Veto detector expression: " << fVetoDetectorsExpression << endl;
        cout << "Veto detector offset: " << fVetoDetectorOffsetSize << endl;
        cout << "Veto light attenuation: " << fVetoLightAttenuation << endl;
    } else {
        cout << "Veto detector expression: not set" << endl;
    }
    cout << "Veto quenching factor: " << fVetoQuenchingFactor << endl;
 
    RESTDebug << RESTendl;
 
    if (fVetoVolumes.empty()) {
        cout << "Process not initialized yet" << endl;
        return;
    }
 
    cout << "Number of veto volumes: " << fVetoVolumes.size() << endl;
    cout << "Number of veto detector volumes: " << fVetoDetectorVolumes.size() << endl;
 
    const auto& geometryInfo = fGeant4Metadata->GetGeant4GeometryInfo();
    for (unsigned int i = 0; i < fVetoVolumes.size(); i++) {
        const auto& vetoName = fVetoVolumes[i];
        const auto& vetoPosition = geometryInfo.GetPosition(vetoName);
 
        cout << TString::Format(" - Veto volume: %d - name: '%s' - position: %s mm", i, vetoName.Data(),
                                VectorToString(vetoPosition).c_str())
             << endl;
 
        if (fVetoDetectorVolumes.empty()) {
            continue;
        }
 
        const auto& vetoDetectorName = fVetoDetectorVolumes[i];
        const auto& vetoDetectorPosition = geometryInfo.GetPosition(vetoDetectorName);
 
        cout << TString::Format("   Veto detector name: '%s' - position: %s mm", vetoDetectorName.Data(),
                                VectorToString(vetoDetectorPosition).c_str())
             << endl;
 
        cout << TString::Format("   Boundary position: %s mm - direction: %s",
                                VectorToString(fVetoDetectorBoundaryPosition.at(vetoName)).c_str(),
                                VectorToString(fVetoDetectorBoundaryDirection.at(vetoName)).c_str());
    }
}