src/TRestGeant4Metadata.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 "TRestGeant4Metadata.h"
 
#include <TAxis.h>
#include <TGeoManager.h>
#include <TObjString.h>
#include <TRandom.h>
#include <TRestGDMLParser.h>
#include <TRestRun.h>
 
#include "TRestGeant4ParticleSourceCosmics.h"
#include "TRestGeant4PrimaryGeneratorInfo.h"
 
using namespace std;
using namespace TRestGeant4PrimaryGeneratorTypes;
 
ClassImp(TRestGeant4Metadata);
 
TRestGeant4Metadata::TRestGeant4Metadata() : TRestMetadata() { Initialize(); }
 
TRestGeant4Metadata::TRestGeant4Metadata(const char* configFilename, const string& name)
    : TRestMetadata(configFilename) {
    Initialize();
 
    LoadConfigFromFile(fConfigFileName, name);
}
 
TRestGeant4Metadata::~TRestGeant4Metadata() { RemoveParticleSources(); }
 
void TRestGeant4Metadata::Initialize() {
    SetSectionName(this->ClassName());
    SetLibraryVersion(LIBRARY_VERSION);
 
    fChance.clear();
    fActiveVolumes.clear();
    fBiasingVolumes.clear();
 
    RemoveParticleSources();
}
 
void TRestGeant4Metadata::InitFromConfigFile() {
    this->Initialize();
 
    fMagneticField = Get3DVectorParameterWithUnits("magneticField", TVector3(0, 0, 0));
 
    // Initialize the metadata members from a configfile
    fGdmlFilename = GetParameter("gdmlFile");
    if (fGdmlFilename == PARAMETER_NOT_FOUND_STR) {
        fGdmlFilename = GetParameter("gdml_file");  // old name
    }
    if (fGdmlFilename == PARAMETER_NOT_FOUND_STR) {
        cout << "\"gdmlFile\" parameter is not defined!" << endl;
        exit(1);
    }
    if (TRestTools::isURL(fGdmlFilename.Data())) {
        fGeometryPath = "";
        fGdmlFilename = TRestTools::DownloadRemoteFile(fGdmlFilename.Data());
        if (fGdmlFilename == "") {
            cout << "Error downloading geometry file from URL: " << fGdmlFilename << endl;
            exit(1);
        }
    }
 
    fGeometryPath = GetParameter("geometryPath", "");
 
    fStoreHadronicTargetInfo = StringToBool(GetParameter("storeHadronicTargetInfo", "false"));
 
    string seedString = GetParameter("seed", "0");
    if (ToUpper(seedString) == "RANDOM" || ToUpper(seedString) == "RAND" || ToUpper(seedString) == "AUTO" ||
        seedString == "0") {
        auto dd = new double();
        fSeed = (uintptr_t)dd + (uintptr_t)this;
        delete dd;
    } else {
        fSeed = (Long_t)StringToInteger(seedString);
    }
    gRandom->SetSeed(fSeed);
 
    // if "gdmlFile" is purely a file (without any path) and "geometryPath" is
    // defined, we recombine them together
    if ((((string)fGdmlFilename).find_first_not_of("./~") == 0 ||
         ((string)fGdmlFilename).find('/') == string::npos) &&
        fGeometryPath != "") {
        if (fGeometryPath[fGeometryPath.Length() - 1] == '/') {
            fGdmlFilename = fGeometryPath + GetParameter("gdmlFile");
        } else {
            fGdmlFilename = fGeometryPath + "/" + GetParameter("gdmlFile");
        }
    }
 
    Double_t defaultTime = 1. / REST_Units::s;
    fSubEventTimeDelay = GetDblParameterWithUnits("subEventTimeDelay", defaultTime);
 
    auto nEventsString = GetParameter("nEvents");
    if (nEventsString == PARAMETER_NOT_FOUND_STR) {
        nEventsString = GetParameter("Nevents");  // old name
    }
    fNEvents = nEventsString == PARAMETER_NOT_FOUND_STR ? 0 : StringToInteger(nEventsString);
 
    fStoreTracks = ToUpper(GetParameter("storeTracks", "true")) == "TRUE" ||
                   ToUpper(GetParameter("storeTracks", "on")) == "ON";
 
    const auto fNRequestedEntriesString = GetParameter("nRequestedEntries");
    fNRequestedEntries =
        fNRequestedEntriesString == PARAMETER_NOT_FOUND_STR ? 0 : StringToInteger(fNRequestedEntriesString);
 
    fSaveAllEvents = ToUpper(GetParameter("saveAllEvents", "false")) == "TRUE" ||
                     ToUpper(GetParameter("saveAllEvents", "off")) == "ON";
 
    fPrintProgress = ToUpper(GetParameter("printProgress", "true")) == "TRUE" ||
                     ToUpper(GetParameter("printProgress", "on")) == "ON";
 
    fRegisterEmptyTracks = ToUpper(GetParameter("registerEmptyTracks", "false")) == "TRUE" ||
                           ToUpper(GetParameter("registerEmptyTracks", "off")) == "ON";
 
    ReadGenerator();
    // Detector (old storage) section is processed after initializing geometry info in Detector Construction
    // This allows to use regular expression to match logical or physical volumes etc.
    ReadBiasing();
 
    fMaxTargetStepSize = GetDblParameterWithUnits("maxTargetStepSize", -1);
    if (fMaxTargetStepSize > 0) {
        cout << " " << endl;
        RESTWarning << "IMPORTANT: *maxTargetStepSize* parameter is now obsolete!" << RESTendl;
        RESTWarning << "The sensitive volume will not define any integration step limit" << RESTendl;
        cout << " " << endl;
        RESTWarning << "In order to avoid this warning REMOVE the *maxTargetStepSize* definition,"
                    << RESTendl;
        RESTWarning << "and replace it by the following statement at the <storage> section" << RESTendl;
        cout << " " << endl;
        RESTWarning << "<activeVolume name=\"" << this->GetSensitiveVolume() << "\" maxStepSize=\""
                    << fMaxTargetStepSize << "mm\" />" << RESTendl;
        cout << " " << endl;
        RESTInfo << "Now, any active volume is allowed to define a maxStepSize" << RESTendl;
        cout << " " << endl;
        RESTInfo << "It is also possible to define a default step size for all active volumes," << RESTendl;
        RESTInfo << "so that in case no step is defined, the default will be used." << RESTendl;
        cout << " " << endl;
        RESTInfo << "<storage sensitiveVolume=\"" << this->GetSensitiveVolume() << "\" maxStepSize=\""
                 << fMaxTargetStepSize << "mm\" />" << RESTendl;
        cout << " " << endl;
        GetChar();
    }
 
    // should return success or fail
}
 
 
Double_t TRestGeant4Metadata::GetCosmicFluxInCountsPerCm2PerSecond() const {
    const auto source = GetParticleSource();
 
    double countsPerSecondPerCm2 = 0;
    if (TRestGeant4PrimaryGeneratorTypes::StringToEnergyDistributionTypes(
            source->GetEnergyDistributionType().Data()) ==
            TRestGeant4PrimaryGeneratorTypes::EnergyDistributionTypes::FORMULA2 &&
        TRestGeant4PrimaryGeneratorTypes::StringToAngularDistributionTypes(
            source->GetAngularDistributionType().Data()) ==
            TRestGeant4PrimaryGeneratorTypes::AngularDistributionTypes::FORMULA2) {
        const auto energyRange = source->GetEnergyDistributionRange();
        const auto angularRange = source->GetAngularDistributionRange();
        auto function = (TF2*)source->GetEnergyAndAngularDistributionFunction()->Clone();
        // counts per second per cm2 (distribution is already integrated over uniform phi)
        countsPerSecondPerCm2 =
            function->Integral(energyRange.X(), energyRange.Y(), angularRange.X(), angularRange.Y(), 1E-9);
    }
 
    else if (std::string(source->GetName()) == "TRestGeant4ParticleSourceCosmics") {
        auto cosmicSource = dynamic_cast<TRestGeant4ParticleSourceCosmics*>(source);
        const auto names = cosmicSource->GetParticleNames();
        const auto histograms = cosmicSource->GetHistogramsTransformed();
        for (const auto& name : names) {
            countsPerSecondPerCm2 += histograms.at(name)->Integral();
        }
    } else if (TRestGeant4PrimaryGeneratorTypes::StringToEnergyDistributionTypes(
                   source->GetEnergyDistributionType().Data()) ==
                   TRestGeant4PrimaryGeneratorTypes::EnergyDistributionTypes::TH2D &&
               TRestGeant4PrimaryGeneratorTypes::StringToAngularDistributionTypes(
                   source->GetAngularDistributionType().Data()) ==
                   TRestGeant4PrimaryGeneratorTypes::AngularDistributionTypes::TH2D) {
        const auto filename = source->GetEnergyDistributionFilename();
        const auto name = source->GetEnergyDistributionNameInFile();
        TFile* file = TFile::Open(filename);
        auto energyAndAngularDistributionHistogram =
            file->Get<TH2D>(name);  // it's the same for both angular and energy
        if (energyAndAngularDistributionHistogram == nullptr) {
            throw std::runtime_error("Histogram not found in file");
        }
        countsPerSecondPerCm2 = energyAndAngularDistributionHistogram->Integral();
        file->Close();
        delete file;
    }
 
    else {
        throw std::runtime_error("Cosmic flux calculation is only supported for TFormula2 or TH2D sources");
    }
 
    return countsPerSecondPerCm2;
}
 
Double_t TRestGeant4Metadata::GetCosmicIntensityInCountsPerSecond() const {
    const auto countsPerSecondPerCm2 = GetCosmicFluxInCountsPerCm2PerSecond();
    const auto surface = GetGeneratorSurfaceCm2();
    const auto countsPerSecond = countsPerSecondPerCm2 * surface;
    return countsPerSecond;
}
 
Double_t TRestGeant4Metadata::GetEquivalentSimulatedTime() const {
    const auto type = ToLower(fGeant4PrimaryGeneratorInfo.GetSpatialGeneratorType());
 
    double scalingFactor = 1.0;
    if (type == "cosmic") {
        // get the cosmic generator
        auto cosmicSource = dynamic_cast<TRestGeant4ParticleSourceCosmics*>(GetParticleSource());
        if (cosmicSource != nullptr) {
            scalingFactor =
                cosmicSource
                    ->GetEnergyRangeScalingFactor();  // number less than 1, to account for energy range
            if (scalingFactor < 0 || scalingFactor > 1) {
                throw std::runtime_error("Energy range scaling factor must be between 0 and 1");
            }
        }
    }
 
    // counts per seconds should be reduced proportionally to the range we are sampling
    const auto countsPerSecond = GetCosmicIntensityInCountsPerSecond() * scalingFactor;
    const auto seconds = double(fNEvents) / countsPerSecond;
 
    return seconds;
}
 
void TRestGeant4Metadata::ReadBiasing() {
    TiXmlElement* biasingDefinition = GetElement("biasing");
    if (biasingDefinition == nullptr) {
        fNBiasingVolumes = 0;
        return;
    }
    TString biasEnabled = GetFieldValue("value", biasingDefinition);
    TString biasType = GetFieldValue("type", biasingDefinition);
 
    RESTDebug << "Bias: " << biasEnabled << RESTendl;
 
    if (biasEnabled == "on" || biasEnabled == "ON" || biasEnabled == "On" || biasEnabled == "oN") {
        RESTInfo << "Biasing is enabled" << RESTendl;
 
        TiXmlElement* biasVolumeDefinition = GetElement("biasingVolume", biasingDefinition);
        Int_t n = 0;
        while (biasVolumeDefinition) {
            TRestGeant4BiasingVolume biasVolume;
            RESTDebug << "Def: " << biasVolumeDefinition << RESTendl;
 
            biasVolume.SetBiasingVolumePosition(
                Get3DVectorParameterWithUnits("position", biasVolumeDefinition));
            biasVolume.SetBiasingFactor(StringToDouble(GetParameter("factor", biasVolumeDefinition)));
            biasVolume.SetBiasingVolumeSize(GetDblParameterWithUnits("size", biasVolumeDefinition));
            biasVolume.SetEnergyRange(Get2DVectorParameterWithUnits("energyRange", biasVolumeDefinition));
            biasVolume.SetBiasingVolumeType(biasType);  // For the moment all the volumes should be same type
 
            /* TODO check that values are right if not printBiasingVolume with
            getchar() biasVolume.PrintBiasingVolume(); getchar();
            */
 
            fBiasingVolumes.push_back(biasVolume);
            biasVolumeDefinition = GetNextElement(biasVolumeDefinition);
            n++;
        }
        fNBiasingVolumes = n;
    }
}
 
void TRestGeant4Metadata::ReadGenerator() {
    // TODO if some fields are defined in the generator but not finally used
    // i.e. <generator type="volume" from="gasTarget" position="(100,0,-100)
    // here position is irrelevant since the event will be generated from the
    // volume defined in the geometry we should take care of these values so they
    // are not stored in metadata (to avoid future confusion) In the particular
    // case of position, the value is overwritten in DetectorConstruction by the
    // center of the volume (i.e. gasTarget) but if i.e rotation or side are
    // defined and not relevant we should set it to -1
 
    TiXmlElement* generatorDefinition = GetElement("generator");
 
    fGeant4PrimaryGeneratorInfo.fSpatialGeneratorType =
        SpatialGeneratorTypesToString(StringToSpatialGeneratorTypes(GetParameter(
            "type", generatorDefinition, SpatialGeneratorTypesToString(SpatialGeneratorTypes::VOLUME))));
    fGeant4PrimaryGeneratorInfo.fSpatialGeneratorShape =
        SpatialGeneratorShapesToString(StringToSpatialGeneratorShapes(GetParameter(
            "shape", generatorDefinition, SpatialGeneratorShapesToString(SpatialGeneratorShapes::BOX))));
    fGeant4PrimaryGeneratorInfo.fSpatialGeneratorFrom = GetParameter("from", generatorDefinition);
    if (fGeant4PrimaryGeneratorInfo.fSpatialGeneratorFrom != PARAMETER_NOT_FOUND_STR) {
        fGeant4PrimaryGeneratorInfo.fSpatialGeneratorShape =
            SpatialGeneratorShapesToString(SpatialGeneratorShapes::GDML);
    }
    fGeant4PrimaryGeneratorInfo.fSpatialGeneratorSize =
        Get3DVectorParameterWithUnits("size", generatorDefinition, {0, 0, 0});
    fGeant4PrimaryGeneratorInfo.fSpatialGeneratorPosition =
        Get3DVectorParameterWithUnits("position", generatorDefinition, {0, 0, 0});
    fGeant4PrimaryGeneratorInfo.fSpatialGeneratorRotationAxis =
        Get3DVectorParameterWithUnits("rotationAxis", generatorDefinition, {0, 0, 1});
    fGeant4PrimaryGeneratorInfo.fSpatialGeneratorRotationValue =
        GetDblParameterWithUnits("rotationAngle", generatorDefinition, 0);
    fGeant4PrimaryGeneratorInfo.fSpatialGeneratorSpatialDensityFunction =
        GetParameter("densityFunc", generatorDefinition, "1");
 
    TiXmlElement* sourceDefinition = GetElement("source", generatorDefinition);
    while (sourceDefinition) {
        string use = GetParameter("use", sourceDefinition, "");
 
        TRestGeant4ParticleSource* source = TRestGeant4ParticleSource::instantiate(use);
        ReadParticleSource(source, sourceDefinition);
        AddParticleSource(source);
 
        if (string(source->GetName()) == "TRestGeant4ParticleSourceCosmics") {
            TRestGeant4ParticleSourceCosmics::SetSeed(GetSeed());
        }
 
        sourceDefinition = GetNextElement(sourceDefinition);
    }
 
    // check if the generator is valid.
    if (GetNumberOfSources() == 0) {
        RESTError << "No sources are added inside generator!" << RESTendl;
        exit(1);
    }
}
 
void TRestGeant4Metadata::ReadParticleSource(TRestGeant4ParticleSource* source, TiXmlElement* definition) {
    TiXmlElement* sourceDefinition = definition;
 
    source->SetParticleName(GetParameter("particle", sourceDefinition));
    source->SetExcitationLevel(StringToDouble(GetParameter("excitedLevel", sourceDefinition)));
    source->SetParticleCharge(StringToInteger(GetParameter("charge", sourceDefinition, "0")));
 
    TString fullChain = GetParameter("fullChain", sourceDefinition);
    SetFullChain(false);
    if (fullChain.EqualTo("on", TString::ECaseCompare::kIgnoreCase)) {
        SetFullChain(true);
    }
 
    const TString fullChainStopIsotopesString = GetParameter("fullChainStopIsotopes", sourceDefinition, "");
    if (fullChainStopIsotopesString != "") {
        TObjArray* fullChainStopIsotopesArray = fullChainStopIsotopesString.Tokenize(",");
        for (int i = 0; i < fullChainStopIsotopesArray->GetEntries(); i++) {
            TString isotope = ((TObjString*)fullChainStopIsotopesArray->At(i))->String();
            fFullChainStopIsotopes.insert(isotope.Data());
        }
    }
 
    // Angular distribution parameters
    TiXmlElement* angularDefinition = GetElement("angular", sourceDefinition);
    if (angularDefinition == nullptr) {
        angularDefinition = GetElement("angularDist", sourceDefinition);  // old name
    }
    source->SetAngularDistributionType(GetParameter(
        "type", angularDefinition, AngularDistributionTypesToString(AngularDistributionTypes::FLUX)));
    if ((StringToAngularDistributionTypes(source->GetAngularDistributionType().Data()) ==
         AngularDistributionTypes::TH1D) ||
        (StringToAngularDistributionTypes(source->GetAngularDistributionType().Data()) ==
         AngularDistributionTypes::TH2D)) {
        source->SetAngularDistributionFilename(SearchFile(GetParameter("file", angularDefinition)));
        auto name = GetParameter("name", angularDefinition);
        if (name == "NO_SUCH_PARA") {
            name = GetParameter("spctName", angularDefinition);  // old name
        }
        source->SetAngularDistributionNameInFile(name);
    }
    source->SetAngularDistributionRange(
        Get2DVectorParameterWithUnits("range", angularDefinition, {0, TMath::Pi()}));
    if (StringToAngularDistributionTypes(source->GetAngularDistributionType().Data()) ==
        AngularDistributionTypes::FORMULA) {
        source->SetAngularDistributionFormula(GetParameter("name", angularDefinition));
        // We cannot use an angular range bigger than the range of the formula
        const auto function = source->GetAngularDistributionFunction();
        if (source->GetAngularDistributionRangeMin() < function->GetXaxis()->GetXmin()) {
            source->SetAngularDistributionRange(
                {function->GetXaxis()->GetXmin(), source->GetAngularDistributionRangeMax()});
        }
        if (source->GetAngularDistributionRangeMax() > function->GetXaxis()->GetXmax()) {
            source->SetAngularDistributionRange(
                {source->GetAngularDistributionRangeMin(), function->GetXaxis()->GetXmax()});
        }
    }
    if ((StringToAngularDistributionTypes(source->GetAngularDistributionType().Data()) ==
         AngularDistributionTypes::FORMULA) ||
        (StringToAngularDistributionTypes(source->GetAngularDistributionType().Data()) ==
         AngularDistributionTypes::FORMULA2)) {
        source->SetAngularDistributionFormulaNPoints(static_cast<size_t>(GetDblParameterWithUnits(
            "nPoints", angularDefinition, source->GetAngularDistributionFormulaNPoints())));
    }
    if (GetNumberOfSources() == 0 &&
        StringToAngularDistributionTypes(source->GetAngularDistributionType().Data()) ==
            AngularDistributionTypes::BACK_TO_BACK) {
        cout << "WARNING: First source cannot be backtoback. Setting it to isotropic" << endl;
        source->SetAngularDistributionType(
            AngularDistributionTypesToString(AngularDistributionTypes::ISOTROPIC));
    }
    source->SetDirection(StringTo3DVector(GetParameter("direction", angularDefinition, "(1,0,0)")));
 
    if ((StringToAngularDistributionTypes(source->GetAngularDistributionType().Data()) ==
         AngularDistributionTypes::ISOTROPIC)) {
        source->SetAngularDistributionRange({0, TMath::Pi()});
        const double coneHalfAngle = GetDblParameterWithUnits("coneHalfAngle", angularDefinition, 0);
        source->SetAngularDistributionIsotropicConeHalfAngle(coneHalfAngle);
    }
    // Energy distribution parameters
    TiXmlElement* energyDefinition = GetElement("energy", sourceDefinition);
    if (energyDefinition == nullptr) {
        energyDefinition = GetElement("energyDist", sourceDefinition);  // old name
    }
    source->SetEnergyDistributionType(GetParameter(
        "type", energyDefinition, EnergyDistributionTypesToString(EnergyDistributionTypes::MONO)));
    if ((StringToEnergyDistributionTypes(source->GetEnergyDistributionType().Data()) ==
         EnergyDistributionTypes::TH1D) ||
        (StringToEnergyDistributionTypes(source->GetEnergyDistributionType().Data()) ==
         EnergyDistributionTypes::TH2D)) {
        source->SetEnergyDistributionFilename(SearchFile(GetParameter("file", energyDefinition)));
        auto name = GetParameter("name", energyDefinition);
        if (name == "NO_SUCH_PARA") {
            name = GetParameter("spctName", energyDefinition);  // old name
        }
        source->SetEnergyDistributionNameInFile(name);
    }
    source->SetEnergyDistributionRange(Get2DVectorParameterWithUnits("range", energyDefinition, {0, 1E20}));
    if (StringToEnergyDistributionTypes(source->GetEnergyDistributionType().Data()) ==
        EnergyDistributionTypes::MONO) {
        Double_t energy;
        energy = GetDblParameterWithUnits("energy", energyDefinition, 0);
        source->SetEnergyDistributionRange({energy, energy});
        source->SetEnergy(energy);
    }
    if (StringToEnergyDistributionTypes(source->GetEnergyDistributionType().Data()) ==
        EnergyDistributionTypes::FORMULA) {
        source->SetEnergyDistributionFormula(GetParameter("name", energyDefinition));
        // We cannot use an energy range bigger than the range of the formula
        const auto function = source->GetEnergyDistributionFunction();
        if (source->GetEnergyDistributionRangeMin() < function->GetXaxis()->GetXmin()) {
            source->SetEnergyDistributionRange(
                {function->GetXaxis()->GetXmin(), source->GetEnergyDistributionRangeMax()});
        }
        if (source->GetEnergyDistributionRangeMax() > function->GetXaxis()->GetXmax()) {
            source->SetEnergyDistributionRange(
                {source->GetEnergyDistributionRangeMin(), function->GetXaxis()->GetXmax()});
        }
    }
    if ((StringToEnergyDistributionTypes(source->GetEnergyDistributionType().Data()) ==
         EnergyDistributionTypes::FORMULA) ||
        (StringToEnergyDistributionTypes(source->GetEnergyDistributionType().Data()) ==
         EnergyDistributionTypes::FORMULA2)) {
        source->SetEnergyDistributionFormulaNPoints(static_cast<size_t>(GetDblParameterWithUnits(
            "nPoints", energyDefinition, source->GetEnergyDistributionFormulaNPoints())));
    }
    if (StringToEnergyDistributionTypes(source->GetEnergyDistributionType().Data()) ==
            EnergyDistributionTypes::FORMULA2 &&
        StringToAngularDistributionTypes(source->GetAngularDistributionType().Data()) ==
            AngularDistributionTypes::FORMULA2) {
        const auto empty = TString("");
        auto energyDistName = GetParameter("name", energyDefinition, empty);
        auto angularDistName = GetParameter("name", energyDefinition, empty);
 
        if (energyDistName == empty && angularDistName == empty) {
            RESTError << "No name specified for energy and angular distribution" << RESTendl;
            exit(1);
        }
        if (energyDistName == empty) {
            angularDistName = energyDistName;
            RESTWarning << "No name specified for energy distribution. Using angular distribution name: "
                        << angularDistName << RESTendl;
        }
        if (angularDistName == empty) {
            energyDistName = energyDistName;
            RESTWarning << "No name specified for angular distribution. Using energy distribution name: "
                        << energyDistName << RESTendl;
        }
 
        if (energyDistName == empty || angularDistName == empty) {
            // we should never enter here, just leave this as a check
            RESTError << "When using 'formula2' the name of energy and angular dist must match" << RESTendl;
            exit(1);
        }
        const auto energyAndAngularDistName = energyDistName;
        source->SetEnergyAndAngularDistributionFormula(energyAndAngularDistName);
 
        const auto function = source->GetEnergyAndAngularDistributionFunction();
 
        // Set energy range
        if (source->GetEnergyDistributionRangeMin() < function->GetXaxis()->GetXmin()) {
            source->SetEnergyDistributionRange(
                {function->GetXaxis()->GetXmin(), source->GetEnergyDistributionRangeMax()});
        }
        if (source->GetEnergyDistributionRangeMax() > function->GetXaxis()->GetXmax()) {
            source->SetEnergyDistributionRange(
                {source->GetEnergyDistributionRangeMin(), function->GetXaxis()->GetXmax()});
        }
 
        // Set angular range
        if (source->GetAngularDistributionRangeMin() < function->GetYaxis()->GetXmin()) {
            source->SetAngularDistributionRange(
                {function->GetYaxis()->GetXmin(), source->GetAngularDistributionRangeMax()});
        }
        if (source->GetAngularDistributionRangeMax() > function->GetYaxis()->GetXmax()) {
            source->SetAngularDistributionRange(
                {source->GetAngularDistributionRangeMin(), function->GetYaxis()->GetXmax()});
        }
    }
    // allow custom configuration from the class
    source->LoadConfigFromElement(sourceDefinition, fElementGlobal, fVariables);
    // AddParticleSource(source);
}
 
void TRestGeant4Metadata::RemoveParticleSources() {
    for (auto source : fParticleSource) {
        delete source;
    }
    fParticleSource.clear();
}
 
void TRestGeant4Metadata::AddParticleSource(TRestGeant4ParticleSource* src) {
    fParticleSource.push_back(src);
}
 
void TRestGeant4Metadata::ReadDetector() {
    RESTInfo << "TRestGeant4Metadata: Processing detector section" << RESTendl;
    TiXmlElement* detectorDefinition = GetElement("detector");
    if (detectorDefinition == nullptr) {
        RESTError << "Detector section (<detector>) is not present in metadata!" << RESTendl;
        exit(1);
    }
 
    const bool activateAllVolumes =
        StringToBool(GetParameter("activateAllVolumes", detectorDefinition, "true"));
    RESTInfo << "TRestGeant4Metadata: Setting 'fActivateAllVolumes' to " << fActivateAllVolumes << RESTendl;
    fActivateAllVolumes = activateAllVolumes;
 
    TiXmlElement* removeUnwantedTracksDefinition = GetElement("removeUnwantedTracks", detectorDefinition);
    if (removeUnwantedTracksDefinition != nullptr) {
        fRemoveUnwantedTracks = StringToBool(GetParameter("enabled", removeUnwantedTracksDefinition, "true"));
        fRemoveUnwantedTracksKeepZeroEnergyTracks =
            StringToBool(GetParameter("keepZeroEnergyTracks", removeUnwantedTracksDefinition, "false"));
        RESTInfo << "TRestGeant4Metadata: Setting 'fRemoveUnwantedTracks' to " << fRemoveUnwantedTracks
                 << " with 'keepZeroEnergyTracks' option set to " << fRemoveUnwantedTracksKeepZeroEnergyTracks
                 << RESTendl;
    }
 
    Double_t defaultStep = GetDblParameterWithUnits("maxStepSize", detectorDefinition);
    if (defaultStep < 0) {
        defaultStep = 0;
    }
 
    TiXmlElement* volumeDefinition = GetElement("volume", detectorDefinition);
    while (volumeDefinition != nullptr) {
        string name = GetFieldValue("name", volumeDefinition);
        if (name.empty() || name == "Not defined") {
            RESTError << "volume inside detector section defined without name" << RESTendl;
            exit(1);
        }
        vector<string> physicalVolumes;
        if (!fGeant4GeometryInfo.IsValidGdmlName(name)) {
            if (fGeant4GeometryInfo.IsValidLogicalVolume(name)) {
                for (const auto& physical : fGeant4GeometryInfo.GetAllPhysicalVolumesFromLogical(name)) {
                    physicalVolumes.emplace_back(
                        fGeant4GeometryInfo.GetAlternativeNameFromGeant4PhysicalName(physical));
                }
            }
            // does not match a explicit (gdml) physical name or a logical name, perhaps its a regular exp
            if (physicalVolumes.empty()) {
                for (const auto& physical :
                     fGeant4GeometryInfo.GetAllPhysicalVolumesMatchingExpression(name)) {
                    physicalVolumes.emplace_back(
                        fGeant4GeometryInfo.GetAlternativeNameFromGeant4PhysicalName(physical));
                }
            }
            if (physicalVolumes.empty()) {
                for (const auto& logical : fGeant4GeometryInfo.GetAllLogicalVolumesMatchingExpression(name)) {
                    for (const auto& physical :
                         fGeant4GeometryInfo.GetAllPhysicalVolumesFromLogical(logical)) {
                        physicalVolumes.emplace_back(
                            fGeant4GeometryInfo.GetAlternativeNameFromGeant4PhysicalName(physical));
                    }
                }
            }
        } else {
            physicalVolumes.push_back(name);
        }
 
        if (physicalVolumes.empty()) {
            RESTError
                << "volume '" << name
                << "' is not defined in the geometry. Could not match to a physical, logical volume or regex"
                << RESTendl;
            exit(1);
        }
 
        bool isSensitive = false;
        const string isSensitiveValue = GetFieldValue("sensitive", volumeDefinition);
        if (isSensitiveValue != "Not defined") {
            isSensitive = StringToBool(isSensitiveValue);
        }
 
        bool isKeepTracks = isSensitive;
        const string isKeepTracksValue = GetFieldValue("keepTracks", volumeDefinition);
        if (isKeepTracksValue != "Not defined") {
            isKeepTracks = StringToBool(isKeepTracksValue);
        }
 
        bool isKill = false;
        const string isKillValue = GetFieldValue("kill", volumeDefinition);
        if (isKillValue != "Not defined") {
            isKill = StringToBool(isKillValue);
        }
 
        Double_t chance = StringToDouble(GetFieldValue("chance", volumeDefinition));
        if (chance == -1 || chance < 0) {
            chance = 1.0;
        }
        Double_t maxStep = GetDblParameterWithUnits("maxStepSize", volumeDefinition);
        if (maxStep < 0) {
            maxStep = defaultStep;
        }
 
        for (const auto& physical : physicalVolumes) {
            RESTInfo << "Adding " << (isSensitive ? "sensitive" : "active") << " volume from RML: '"
                     << physical << (chance != 1 ? " with change: " + to_string(chance) : " ") << RESTendl;
            SetActiveVolume(physical, chance, maxStep);
            if (isSensitive) {
                InsertSensitiveVolume(physical);
            }
            if (fRemoveUnwantedTracks && isKeepTracks) {
                fRemoveUnwantedTracksVolumesToKeep.insert(physical);
            }
            if (isKill) {
                fKillVolumes.insert(physical);
            }
        }
        volumeDefinition = GetNextElement(volumeDefinition);
    }
 
    if (fSensitiveVolumes.empty()) {
        cout << "No sensitive volumes defined in detector section. Adding 'gas' as sensitive volume" << endl;
        InsertSensitiveVolume("gas");
    }
 
    fEnergyRangeStored = Get2DVectorParameterWithUnits("energyRange", detectorDefinition, fEnergyRangeStored);
    // TODO: Place this as a generic method
    if (fEnergyRangeStored.X() < 0) {
        fEnergyRangeStored.Set(0, fEnergyRangeStored.Y());
    }
    if (fEnergyRangeStored.Y() < 0) {
        fEnergyRangeStored.Set(fEnergyRangeStored.X(), 0);
    }
    if (fEnergyRangeStored.Y() > fEnergyRangeStored.Y()) {
        fEnergyRangeStored.Set(fEnergyRangeStored.Y(), fEnergyRangeStored.X());
    }
    if (fEnergyRangeStored.Y() <= 0) {
        RESTError << "Energy range is not valid: (" << fEnergyRangeStored.X() << ", "
                  << fEnergyRangeStored.Y() << ") keV" << RESTendl;
        exit(1);
    }
 
    auto gdml = new TRestGDMLParser();
    gdml->Load(GetGdmlFilename().Data());
 
    fGeant4GeometryInfo.PopulateFromGdml(gdml->GetOutputGDMLFile());
 
    // If the user did not add explicitly any volume to the storage section we understand
    // the user wants to register all the volumes
    if (fActivateAllVolumes) {
        for (const auto& name : fGeant4GeometryInfo.GetAllPhysicalVolumes()) {
            if (!IsActiveVolume(name)) {
                SetActiveVolume(name, 1, defaultStep);
                RESTInfo << "Automatically adding active volume: '" << name << "' with chance: " << 1
                         << RESTendl;
            }
        }
    }
 
    if (GetNumberOfActiveVolumes() == 0) {
        RESTError << "No active volumes defined. Please check the detector section" << RESTendl;
        exit(1);
    }
}
 
void TRestGeant4Metadata::PrintMetadata() {
    TRestMetadata::PrintMetadata();
 
    RESTMetadata << "Geant4 version: " << GetGeant4Version() << RESTendl;
    RESTMetadata << "Random seed: " << GetSeed() << RESTendl;
    RESTMetadata << "GDML geometry: " << GetGdmlReference() << RESTendl;
    RESTMetadata << "GDML materials reference: " << GetMaterialsReference() << RESTendl;
    RESTMetadata << "Sub-event time delay: " << GetSubEventTimeDelay() << " us" << RESTendl;
    Double_t mx = GetMagneticField().X();
    Double_t my = GetMagneticField().Y();
    Double_t mz = GetMagneticField().Z();
    RESTMetadata << "Magnetic field: (" << mx << ", " << my << ", " << mz << ") T" << RESTendl;
    if (fSaveAllEvents) RESTMetadata << "Save all events was enabled!" << RESTendl;
    if (fRegisterEmptyTracks) {
        RESTMetadata << "Register empty tracks was enabled" << RESTendl;
    } else {
        RESTMetadata << "Register empty tracks was NOT enabled" << RESTendl;
    }
 
    RESTMetadata << " " << RESTendl;
    RESTMetadata << "   ++++++++++ Generator +++++++++++   " << RESTendl;
    RESTMetadata << "Number of generated events: " << GetNumberOfEvents() << RESTendl;
    fGeant4PrimaryGeneratorInfo.Print();
 
    for (int i = 0; i < GetNumberOfSources(); i++) GetParticleSource(i)->PrintMetadata();
 
    RESTMetadata << " " << RESTendl;
    RESTMetadata << "   ++++++++++ Detector +++++++++++   " << RESTendl;
    RESTMetadata << "Energy range (keV): (" << GetMinimumEnergyStored() << ", " << GetMaximumEnergyStored()
                 << ")" << RESTendl;
    RESTMetadata << "Number of sensitive volumes: " << GetNumberOfSensitiveVolumes() << RESTendl;
    for (const auto& sensitiveVolume : fSensitiveVolumes) {
        RESTMetadata << "Sensitive volume: " << sensitiveVolume << RESTendl;
    }
    RESTMetadata << "Number of active volumes: " << GetNumberOfActiveVolumes() << RESTendl;
    for (unsigned int n = 0; n < GetNumberOfActiveVolumes(); n++) {
        const auto name = GetActiveVolumeName(n);
        RESTMetadata << "Name: " << name << ", ID: " << fGeant4GeometryInfo.GetIDFromVolume(name)
                     << ", maxStep: " << GetMaxStepSize(name) << "mm "
                     << ", chance: " << GetStorageChance(name) << RESTendl;
    }
 
    for (unsigned int n = 0; n < GetNumberOfBiasingVolumes(); n++) {
        GetBiasingVolume(n).PrintBiasingVolume();
    }
 
    if (GetRemoveUnwantedTracks()) {
        RESTMetadata << "removeUnwantedTracks is enabled "
                     << (fRemoveUnwantedTracksKeepZeroEnergyTracks ? "with" : "without")
                     << " keeping zero energy tracks" << RESTendl;
 
        /*
        RESTMetadata << "keep volumes for removeUnwantedTracks:" << RESTendl;
        for (const auto& volume : fRemoveUnwantedTracksVolumesToKeep) {
            RESTMetadata << "\t" << volume << RESTendl;
        }
         */
    }
 
    RESTMetadata << "+++++" << RESTendl;
}
 
Int_t TRestGeant4Metadata::GetActiveVolumeID(const TString& name) {
    Int_t id;
    for (id = 0; id < (Int_t)fActiveVolumes.size(); id++) {
        if (fActiveVolumes[id] == name) return id;
    }
    return -1;
}
 
void TRestGeant4Metadata::SetActiveVolume(const TString& name, Double_t chance, Double_t maxStep) {
    for (size_t i = 0; i < fActiveVolumes.size(); i++) {
        const auto activeVolumeName = fActiveVolumes[i];
        if (name == activeVolumeName) {
            fChance[i] = chance;
            fMaxStepSize[i] = maxStep;
            return;
        }
    }
    fActiveVolumes.push_back(name);
    fChance.push_back(chance);
    fMaxStepSize.push_back(maxStep);
    fActiveVolumesSet.insert(name.Data());
}
 
double TRestGeant4Metadata::GetGeneratorSurfaceCm2() const {
    const auto type = ToLower(fGeant4PrimaryGeneratorInfo.GetSpatialGeneratorType());
    const auto shape = ToLower(fGeant4PrimaryGeneratorInfo.GetSpatialGeneratorShape());
 
    if (type == "surface" && shape == "circle") {
        const auto radius = fGeant4PrimaryGeneratorInfo.GetSpatialGeneratorSize().X();
        return TMath::Pi() * radius * radius * 0.01;  // cm2
    } else if (type == "cosmic") {
        return fGeant4PrimaryGeneratorInfo.GetSpatialGeneratorCosmicSurfaceTermCm2();
    }
 
    throw std::runtime_error(
        "TRestGeant4Metadata::GetGeneratorSurfaceCm2: Can only be called for 'cosmic' and 'surface/circle' "
        "generators");
}
 
Bool_t TRestGeant4Metadata::isVolumeStored(const TString& volume) const {
    for (unsigned int n = 0; n < GetNumberOfActiveVolumes(); n++) {
        if (GetActiveVolumeName(n) == volume) {
            return true;
        }
    }
 
    return false;
}
 
Double_t TRestGeant4Metadata::GetStorageChance(const TString& volume) {
    for (auto id = 0; id < (Int_t)fActiveVolumes.size(); id++) {
        {
            if (fActiveVolumes[id] == volume) {
                return fChance[id];
            }
        }
    }
    RESTWarning << "TRestGeant4Metadata::GetStorageChance. Volume " << volume << " not found" << RESTendl;
 
    return 0;
}
 
Double_t TRestGeant4Metadata::GetMaxStepSize(const TString& volume) {
    for (Int_t i = 0; i < (Int_t)fActiveVolumes.size(); i++) {
        if (volume.EqualTo(fActiveVolumes[i], TString::kIgnoreCase)) {
            return fMaxStepSize[i];
        }
    }
    RESTWarning << "TRestGeant4Metadata::GetMaxStepSize. Volume " << volume << " not found" << RESTendl;
 
    return 0;
}
 
size_t TRestGeant4Metadata::GetGeant4VersionMajor() const {
    TString majorVersion = "";
    for (int i = 0; i < fGeant4Version.Length(); i++) {
        auto c = fGeant4Version[i];
        if (c == '.') {
            break;
        }
        majorVersion += c;
    }
    return std::stoi(majorVersion.Data());
}
 
void TRestGeant4Metadata::Merge(const TRestGeant4Metadata& metadata) {
    fIsMerge = true;
    fSeed = 0;  // seed makes no sense in a merged file
 
    fNEvents += metadata.fNEvents;
    fNRequestedEntries += metadata.fNRequestedEntries;
    fSimulationTime += metadata.fSimulationTime;
}
 
TRestGeant4Metadata::TRestGeant4Metadata(const TRestGeant4Metadata& metadata) : TRestMetadata(metadata) {
    *this = metadata;
}
 
TRestGeant4Metadata& TRestGeant4Metadata::operator=(const TRestGeant4Metadata& metadata) {
    fIsMerge = metadata.fIsMerge;
    fGeant4GeometryInfo = metadata.fGeant4GeometryInfo;
    fGeant4PhysicsInfo = metadata.fGeant4PhysicsInfo;
    fGeant4PrimaryGeneratorInfo = metadata.fGeant4PrimaryGeneratorInfo;
    fGeant4Version = metadata.fGeant4Version;
    fGdmlReference = metadata.fGdmlReference;
    fMaterialsReference = metadata.fMaterialsReference;
    fEnergyRangeStored = metadata.fEnergyRangeStored;
    fActiveVolumes = metadata.fActiveVolumes;
    fChance = metadata.fChance;
    fMaxStepSize = metadata.fMaxStepSize;
    // fParticleSource = metadata.fParticleSource; // segfaults (pointers!)
    fNBiasingVolumes = metadata.fNBiasingVolumes;
    fBiasingVolumes = metadata.fBiasingVolumes;
    fMaxTargetStepSize = metadata.fMaxTargetStepSize;
    fSubEventTimeDelay = metadata.fSubEventTimeDelay;
    fFullChain = metadata.fFullChain;
    fSensitiveVolumes = metadata.fSensitiveVolumes;
    fNEvents = metadata.fNEvents;
    fNRequestedEntries = metadata.fNRequestedEntries;
    fStoreTracks = metadata.fStoreTracks;
    fSimulationMaxTimeSeconds = metadata.fSimulationMaxTimeSeconds;
    fSeed = metadata.fSeed;
    fSaveAllEvents = metadata.fSaveAllEvents;
    fRemoveUnwantedTracks = metadata.fRemoveUnwantedTracks;
    fRemoveUnwantedTracksKeepZeroEnergyTracks = metadata.fRemoveUnwantedTracksKeepZeroEnergyTracks;
    fRemoveUnwantedTracksVolumesToKeep = metadata.fRemoveUnwantedTracksVolumesToKeep;
    fKillVolumes = metadata.fKillVolumes;
    fRegisterEmptyTracks = metadata.fRegisterEmptyTracks;
    fMagneticField = metadata.fMagneticField;
    return *this;
}