735#include "TRestGeant4Metadata.h"
738#include <TGeoManager.h>
740#include <TRestGDMLParser.h>
743#include "TRestGeant4ParticleSourceCosmics.h"
744#include "TRestGeant4PrimaryGeneratorInfo.h"
747using namespace TRestGeant4PrimaryGeneratorTypes;
802 fMagneticField = Get3DVectorParameterWithUnits(
"magneticField", TVector3(0, 0, 0));
810 cout <<
"\"gdmlFile\" parameter is not defined!" << endl;
817 cout <<
"Error downloading geometry file from URL: " <<
fGdmlFilename << endl;
827 if (
ToUpper(seedString) ==
"RANDOM" ||
ToUpper(seedString) ==
"RAND" ||
ToUpper(seedString) ==
"AUTO" ||
829 auto dd =
new double();
830 fSeed = (uintptr_t)dd + (uintptr_t)
this;
835 gRandom->SetSeed(
fSeed);
839 if ((((
string)
fGdmlFilename).find_first_not_of(
"./~") == 0 ||
849 Double_t defaultTime = 1. / REST_Units::s;
853 if (nEventsString == PARAMETER_NOT_FOUND_STR) {
858 const auto fNRequestedEntriesString =
GetParameter(
"nRequestedEntries");
860 fNRequestedEntriesString == PARAMETER_NOT_FOUND_STR ? 0 :
StringToInteger(fNRequestedEntriesString);
879 RESTWarning <<
"IMPORTANT: *maxTargetStepSize* parameter is now obsolete!" <<
RESTendl;
880 RESTWarning <<
"The sensitive volume will not define any integration step limit" <<
RESTendl;
882 RESTWarning <<
"In order to avoid this warning REMOVE the *maxTargetStepSize* definition,"
884 RESTWarning <<
"and replace it by the following statement at the <storage> section" <<
RESTendl;
886 RESTWarning <<
"<activeVolume name=\"" << this->
GetSensitiveVolume() <<
"\" maxStepSize=\""
889 RESTInfo <<
"Now, any active volume is allowed to define a maxStepSize" <<
RESTendl;
891 RESTInfo <<
"It is also possible to define a default step size for all active volumes," <<
RESTendl;
892 RESTInfo <<
"so that in case no step is defined, the default will be used." <<
RESTendl;
894 RESTInfo <<
"<storage sensitiveVolume=\"" << this->
GetSensitiveVolume() <<
"\" maxStepSize=\""
923 if (TRestGeant4PrimaryGeneratorTypes::StringToSpatialGeneratorTypes(
925 TRestGeant4PrimaryGeneratorTypes::SpatialGeneratorTypes::COSMIC) {
927 <<
"TRestGeant4Metadata::GetEquivalentSimulatedTime can only be called for 'cosmic' generator"
932 if (TRestGeant4PrimaryGeneratorTypes::StringToEnergyDistributionTypes(
933 source->GetEnergyDistributionType().Data()) !=
934 TRestGeant4PrimaryGeneratorTypes::EnergyDistributionTypes::FORMULA2) {
935 RESTError <<
"TRestGeant4Metadata::GetEquivalentSimulatedTime can only be called for 'formula2' "
936 "energy distribution"
940 if (TRestGeant4PrimaryGeneratorTypes::StringToAngularDistributionTypes(
941 source->GetAngularDistributionType().Data()) !=
942 TRestGeant4PrimaryGeneratorTypes::AngularDistributionTypes::FORMULA2) {
943 RESTError <<
"TRestGeant4Metadata::GetEquivalentSimulatedTime can only be called for 'formula' "
944 "angular distribution"
949 const auto energyRange = source->GetEnergyDistributionRange();
950 const auto angularRange = source->GetAngularDistributionRange();
951 auto function = (TF2*)source->GetEnergyAndAngularDistributionFunction()->Clone();
953 const auto countsPerSecondPerCm2 =
954 function->Integral(energyRange.X(), energyRange.Y(), angularRange.X(), angularRange.Y(), 1E-9);
955 return countsPerSecondPerCm2;
958Double_t TRestGeant4Metadata::GetCosmicIntensityInCountsPerSecond()
const {
961 const auto countsPerSecond = countsPerSecondPerCm2 * surface;
962 return countsPerSecond;
965Double_t TRestGeant4Metadata::GetEquivalentSimulatedTime()
const {
966 const auto countsPerSecond = GetCosmicIntensityInCountsPerSecond();
967 const auto seconds = double(
fNEvents) / countsPerSecond;
971void TRestGeant4Metadata::ReadBiasing() {
972 TiXmlElement* biasingDefinition =
GetElement(
"biasing");
973 if (biasingDefinition ==
nullptr) {
977 TString biasEnabled =
GetFieldValue(
"value", biasingDefinition);
980 RESTDebug <<
"Bias: " << biasEnabled <<
RESTendl;
982 if (biasEnabled ==
"on" || biasEnabled ==
"ON" || biasEnabled ==
"On" || biasEnabled ==
"oN") {
983 RESTInfo <<
"Biasing is enabled" <<
RESTendl;
985 TiXmlElement* biasVolumeDefinition =
GetElement(
"biasingVolume", biasingDefinition);
987 while (biasVolumeDefinition) {
989 RESTDebug <<
"Def: " << biasVolumeDefinition <<
RESTendl;
991 biasVolume.SetBiasingVolumePosition(
992 Get3DVectorParameterWithUnits(
"position", biasVolumeDefinition));
994 biasVolume.SetBiasingVolumeSize(GetDblParameterWithUnits(
"size", biasVolumeDefinition));
995 biasVolume.SetEnergyRange(Get2DVectorParameterWithUnits(
"energyRange", biasVolumeDefinition));
996 biasVolume.SetBiasingVolumeType(biasType);
1044 TiXmlElement* generatorDefinition =
GetElement(
"generator");
1047 SpatialGeneratorTypesToString(StringToSpatialGeneratorTypes(
GetParameter(
1048 "type", generatorDefinition, SpatialGeneratorTypesToString(SpatialGeneratorTypes::VOLUME))));
1050 SpatialGeneratorShapesToString(StringToSpatialGeneratorShapes(
GetParameter(
1051 "shape", generatorDefinition, SpatialGeneratorShapesToString(SpatialGeneratorShapes::BOX))));
1055 SpatialGeneratorShapesToString(SpatialGeneratorShapes::GDML);
1058 Get3DVectorParameterWithUnits(
"size", generatorDefinition, {0, 0, 0});
1060 Get3DVectorParameterWithUnits(
"position", generatorDefinition, {0, 0, 0});
1062 Get3DVectorParameterWithUnits(
"rotationAxis", generatorDefinition, {0, 0, 1});
1064 GetDblParameterWithUnits(
"rotationAngle", generatorDefinition, 0);
1068 TiXmlElement* sourceDefinition =
GetElement(
"source", generatorDefinition);
1069 while (sourceDefinition) {
1070 string use =
GetParameter(
"use", sourceDefinition,
"");
1076 if (
string(source->GetName()) ==
"TRestGeant4ParticleSourceCosmics") {
1077 TRestGeant4ParticleSourceCosmics::SetSeed(
GetSeed());
1085 RESTError <<
"No sources are added inside generator!" <<
RESTendl;
1094 TiXmlElement* sourceDefinition = definition;
1096 source->SetParticleName(
GetParameter(
"particle", sourceDefinition));
1100 TString fullChain =
GetParameter(
"fullChain", sourceDefinition);
1102 if (fullChain.EqualTo(
"on", TString::ECaseCompare::kIgnoreCase)) {
1107 TiXmlElement* angularDefinition =
GetElement(
"angular", sourceDefinition);
1108 if (angularDefinition ==
nullptr) {
1109 angularDefinition =
GetElement(
"angularDist", sourceDefinition);
1112 "type", angularDefinition, AngularDistributionTypesToString(AngularDistributionTypes::FLUX)));
1113 if ((StringToAngularDistributionTypes(source->GetAngularDistributionType().Data()) ==
1114 AngularDistributionTypes::TH1D) ||
1115 (StringToAngularDistributionTypes(source->GetAngularDistributionType().Data()) ==
1116 AngularDistributionTypes::TH2D)) {
1119 if (name ==
"NO_SUCH_PARA") {
1122 source->SetAngularDistributionNameInFile(name);
1124 source->SetAngularDistributionRange(
1125 Get2DVectorParameterWithUnits(
"range", angularDefinition, {0, TMath::Pi()}));
1126 if (StringToAngularDistributionTypes(source->GetAngularDistributionType().Data()) ==
1127 AngularDistributionTypes::FORMULA) {
1128 source->SetAngularDistributionFormula(
GetParameter(
"name", angularDefinition));
1130 const auto function = source->GetAngularDistributionFunction();
1131 if (source->GetAngularDistributionRangeMin() < function->GetXaxis()->GetXmin()) {
1132 source->SetAngularDistributionRange(
1133 {function->GetXaxis()->GetXmin(), source->GetAngularDistributionRangeMax()});
1135 if (source->GetAngularDistributionRangeMax() > function->GetXaxis()->GetXmax()) {
1136 source->SetAngularDistributionRange(
1137 {source->GetAngularDistributionRangeMin(), function->GetXaxis()->GetXmax()});
1140 if ((StringToAngularDistributionTypes(source->GetAngularDistributionType().Data()) ==
1141 AngularDistributionTypes::FORMULA) ||
1142 (StringToAngularDistributionTypes(source->GetAngularDistributionType().Data()) ==
1143 AngularDistributionTypes::FORMULA2)) {
1144 source->SetAngularDistributionFormulaNPoints(
static_cast<size_t>(GetDblParameterWithUnits(
1145 "nPoints", angularDefinition, source->GetAngularDistributionFormulaNPoints())));
1148 StringToAngularDistributionTypes(source->GetAngularDistributionType().Data()) ==
1149 AngularDistributionTypes::BACK_TO_BACK) {
1150 cout <<
"WARNING: First source cannot be backtoback. Setting it to isotropic" << endl;
1151 source->SetAngularDistributionType(
1152 AngularDistributionTypesToString(AngularDistributionTypes::ISOTROPIC));
1156 if ((StringToAngularDistributionTypes(source->GetAngularDistributionType().Data()) ==
1157 AngularDistributionTypes::ISOTROPIC)) {
1158 source->SetAngularDistributionRange({0, TMath::Pi()});
1159 const double coneHalfAngle = GetDblParameterWithUnits(
"coneHalfAngle", angularDefinition, 0);
1160 source->SetAngularDistributionIsotropicConeHalfAngle(coneHalfAngle);
1163 TiXmlElement* energyDefinition =
GetElement(
"energy", sourceDefinition);
1164 if (energyDefinition ==
nullptr) {
1165 energyDefinition =
GetElement(
"energyDist", sourceDefinition);
1168 "type", energyDefinition, EnergyDistributionTypesToString(EnergyDistributionTypes::MONO)));
1169 if ((StringToEnergyDistributionTypes(source->GetEnergyDistributionType().Data()) ==
1170 EnergyDistributionTypes::TH1D) ||
1171 (StringToEnergyDistributionTypes(source->GetEnergyDistributionType().Data()) ==
1172 EnergyDistributionTypes::TH2D)) {
1175 if (name ==
"NO_SUCH_PARA") {
1178 source->SetEnergyDistributionNameInFile(name);
1180 source->SetEnergyDistributionRange(Get2DVectorParameterWithUnits(
"range", energyDefinition, {0, 1E20}));
1181 if (StringToEnergyDistributionTypes(source->GetEnergyDistributionType().Data()) ==
1182 EnergyDistributionTypes::MONO) {
1184 energy = GetDblParameterWithUnits(
"energy", energyDefinition, 0);
1185 source->SetEnergyDistributionRange({energy, energy});
1186 source->SetEnergy(energy);
1188 if (StringToEnergyDistributionTypes(source->GetEnergyDistributionType().Data()) ==
1189 EnergyDistributionTypes::FORMULA) {
1190 source->SetEnergyDistributionFormula(
GetParameter(
"name", energyDefinition));
1192 const auto function = source->GetEnergyDistributionFunction();
1193 if (source->GetEnergyDistributionRangeMin() < function->GetXaxis()->GetXmin()) {
1194 source->SetEnergyDistributionRange(
1195 {function->GetXaxis()->GetXmin(), source->GetEnergyDistributionRangeMax()});
1197 if (source->GetEnergyDistributionRangeMax() > function->GetXaxis()->GetXmax()) {
1198 source->SetEnergyDistributionRange(
1199 {source->GetEnergyDistributionRangeMin(), function->GetXaxis()->GetXmax()});
1202 if ((StringToEnergyDistributionTypes(source->GetEnergyDistributionType().Data()) ==
1203 EnergyDistributionTypes::FORMULA) ||
1204 (StringToEnergyDistributionTypes(source->GetEnergyDistributionType().Data()) ==
1205 EnergyDistributionTypes::FORMULA2)) {
1206 source->SetEnergyDistributionFormulaNPoints(
static_cast<size_t>(GetDblParameterWithUnits(
1207 "nPoints", energyDefinition, source->GetEnergyDistributionFormulaNPoints())));
1209 if (StringToEnergyDistributionTypes(source->GetEnergyDistributionType().Data()) ==
1210 EnergyDistributionTypes::FORMULA2 &&
1211 StringToAngularDistributionTypes(source->GetAngularDistributionType().Data()) ==
1212 AngularDistributionTypes::FORMULA2) {
1213 const auto empty = TString(
"");
1214 auto energyDistName =
GetParameter(
"name", energyDefinition, empty);
1215 auto angularDistName =
GetParameter(
"name", energyDefinition, empty);
1217 if (energyDistName == empty && angularDistName == empty) {
1218 RESTError <<
"No name specified for energy and angular distribution" <<
RESTendl;
1221 if (energyDistName == empty) {
1222 angularDistName = energyDistName;
1223 RESTWarning <<
"No name specified for energy distribution. Using angular distribution name: "
1226 if (angularDistName == empty) {
1227 energyDistName = energyDistName;
1228 RESTWarning <<
"No name specified for angular distribution. Using energy distribution name: "
1232 if (energyDistName == empty || angularDistName == empty) {
1234 RESTError <<
"When using 'formula2' the name of energy and angular dist must match" <<
RESTendl;
1237 const auto energyAndAngularDistName = energyDistName;
1238 source->SetEnergyAndAngularDistributionFormula(energyAndAngularDistName);
1240 const auto function = source->GetEnergyAndAngularDistributionFunction();
1243 if (source->GetEnergyDistributionRangeMin() < function->GetXaxis()->GetXmin()) {
1244 source->SetEnergyDistributionRange(
1245 {function->GetXaxis()->GetXmin(), source->GetEnergyDistributionRangeMax()});
1247 if (source->GetEnergyDistributionRangeMax() > function->GetXaxis()->GetXmax()) {
1248 source->SetEnergyDistributionRange(
1249 {source->GetEnergyDistributionRangeMin(), function->GetXaxis()->GetXmax()});
1253 if (source->GetAngularDistributionRangeMin() < function->GetYaxis()->GetXmin()) {
1254 source->SetAngularDistributionRange(
1255 {function->GetYaxis()->GetXmin(), source->GetAngularDistributionRangeMax()});
1257 if (source->GetAngularDistributionRangeMax() > function->GetYaxis()->GetXmax()) {
1258 source->SetAngularDistributionRange(
1259 {source->GetAngularDistributionRangeMin(), function->GetYaxis()->GetXmax()});
1293 RESTInfo <<
"TRestGeant4Metadata: Processing detector section" <<
RESTendl;
1294 TiXmlElement* detectorDefinition =
GetElement(
"detector");
1295 if (detectorDefinition ==
nullptr) {
1296 RESTError <<
"Detector section (<detector>) is not present in metadata!" <<
RESTendl;
1300 const bool activateAllVolumes =
1301 StringToBool(
GetParameter(
"activateAllVolumes", detectorDefinition,
"true"));
1305 TiXmlElement* removeUnwantedTracksDefinition =
GetElement(
"removeUnwantedTracks", detectorDefinition);
1306 if (removeUnwantedTracksDefinition !=
nullptr) {
1309 StringToBool(
GetParameter(
"keepZeroEnergyTracks", removeUnwantedTracksDefinition,
"false"));
1315 Double_t defaultStep = GetDblParameterWithUnits(
"maxStepSize", detectorDefinition);
1316 if (defaultStep < 0) {
1320 TiXmlElement* volumeDefinition =
GetElement(
"volume", detectorDefinition);
1321 while (volumeDefinition !=
nullptr) {
1323 if (name.empty() || name ==
"Not defined") {
1324 RESTError <<
"volume inside detector section defined without name" <<
RESTendl;
1327 vector<string> physicalVolumes;
1331 physicalVolumes.emplace_back(
1336 if (physicalVolumes.empty()) {
1337 for (
const auto& physical :
1339 physicalVolumes.emplace_back(
1343 if (physicalVolumes.empty()) {
1344 for (
const auto& logical :
fGeant4GeometryInfo.GetAllLogicalVolumesMatchingExpression(name)) {
1345 for (
const auto& physical :
1347 physicalVolumes.emplace_back(
1353 physicalVolumes.push_back(name);
1356 if (physicalVolumes.empty()) {
1358 <<
"volume '" << name
1359 <<
"' is not defined in the geometry. Could not match to a physical, logical volume or regex"
1364 bool isSensitive =
false;
1365 const string isSensitiveValue =
GetFieldValue(
"sensitive", volumeDefinition);
1366 if (isSensitiveValue !=
"Not defined") {
1367 isSensitive = StringToBool(isSensitiveValue);
1370 bool isKeepTracks = isSensitive;
1371 const string isKeepTracksValue =
GetFieldValue(
"keepTracks", volumeDefinition);
1372 if (isKeepTracksValue !=
"Not defined") {
1373 isKeepTracks = StringToBool(isKeepTracksValue);
1376 bool isKill =
false;
1377 const string isKillValue =
GetFieldValue(
"kill", volumeDefinition);
1378 if (isKillValue !=
"Not defined") {
1379 isKill = StringToBool(isKillValue);
1383 if (chance == -1 || chance < 0) {
1386 Double_t maxStep = GetDblParameterWithUnits(
"maxStepSize", volumeDefinition);
1388 maxStep = defaultStep;
1391 for (
const auto& physical : physicalVolumes) {
1392 RESTInfo <<
"Adding " << (isSensitive ?
"sensitive" :
"active") <<
" volume from RML: '"
1393 << physical << (chance != 1 ?
" with change: " + to_string(chance) :
" ") <<
RESTendl;
1402 fKillVolumes.insert(physical);
1409 cout <<
"No sensitive volumes defined in detector section. Adding 'gas' as sensitive volume" << endl;
1439 if (!IsActiveVolume(name)) {
1441 RESTInfo <<
"Automatically adding active volume: '" << name <<
"' with chance: " << 1
1448 RESTError <<
"No active volumes defined. Please check the detector section" <<
RESTendl;
1468 RESTMetadata <<
"Magnetic field: (" << mx <<
", " << my <<
", " << mz <<
") T" <<
RESTendl;
1471 RESTMetadata <<
"Register empty tracks was enabled" <<
RESTendl;
1473 RESTMetadata <<
"Register empty tracks was NOT enabled" <<
RESTendl;
1477 RESTMetadata <<
" ++++++++++ Generator +++++++++++ " <<
RESTendl;
1484 RESTMetadata <<
" ++++++++++ Detector +++++++++++ " <<
RESTendl;
1487 RESTMetadata <<
"Number of sensitive volumes: " << GetNumberOfSensitiveVolumes() <<
RESTendl;
1489 RESTMetadata <<
"Sensitive volume: " << sensitiveVolume <<
RESTendl;
1503 if (GetRemoveUnwantedTracks()) {
1504 RESTMetadata <<
"removeUnwantedTracks is enabled "
1506 <<
" keeping zero energy tracks" <<
RESTendl;
1516 RESTMetadata <<
"+++++" <<
RESTendl;
1546 if (name == activeVolumeName) {
1555 fActiveVolumesSet.insert(name.Data());
1577 RESTWarning <<
"TRestGeant4Metadata::GetStorageChance. Volume " << volume <<
" not found" <<
RESTendl;
1591 RESTWarning <<
"TRestGeant4Metadata::GetMaxStepSize. Volume " << volume <<
" not found" <<
RESTendl;
1596size_t TRestGeant4Metadata::GetGeant4VersionMajor()
const {
1597 TString majorVersion =
"";
1605 return std::stoi(majorVersion.Data());
1620 fIsMerge = metadata.fIsMerge;
1646 fKillVolumes = metadata.fKillVolumes;
virtual void PrintMetadata() override
Implemented it in the derived metadata class to print out specific metadata information.
TVector3 fSpatialGeneratorSize
The size of the generator. Stores up to three dimensions according to the shape box: length,...
TString GetSpatialGeneratorType() const
Returns a std::string specifying the generator type (volume, surface, point, virtualWall,...
TString fSpatialGeneratorFrom
The volume name where the events are generated, in case of volume or surface generator types.
TString fSpatialGeneratorSpatialDensityFunction
Defines density distribution of the generator shape. rho=F(x,y,z), in range 0~1.
Double_t fSpatialGeneratorRotationValue
degrees of rotation for generator virtual shape along the axis
TVector3 fSpatialGeneratorPosition
The position of the generator for virtual, and point, generator types.
Double_t GetSpatialGeneratorCosmicSurfaceTermCm2() const
Returns cosmic surface term (cm2) for simulation time computation.
TVector3 fSpatialGeneratorRotationAxis
A 3d-std::vector with the angles, measured in degrees, of a XYZ rotation applied to the virtual gener...
TString fSpatialGeneratorShape
Shape of spatial generator (wall, GDML, sphere, etc)
TString fSpatialGeneratorType
Type of spatial generator (point, surface, volume, custom)
Int_t GetChar(std::string hint="Press a KEY to continue ...")
Helps to pause the program, printing a message before pausing.
Double_t StringToDouble(std::string in)
Gets a double from a string.
std::string ToUpper(std::string in)
Convert string to its upper case. Alternative of TString::ToUpper.
Int_t StringToInteger(std::string in)
Gets an integer from a string.
TVector3 StringTo3DVector(std::string in)
Gets a 3D-vector from a string. Format should be : (X,Y,Z).