diff --git a/inc/TRestGeant4QuenchingProcess.h b/inc/TRestGeant4QuenchingProcess.h index b76fe08..726e561 100644 --- a/inc/TRestGeant4QuenchingProcess.h +++ b/inc/TRestGeant4QuenchingProcess.h @@ -25,10 +25,14 @@ #include +#include +#include +#include + #include "TRestGeant4Event.h" #include "TRestGeant4Metadata.h" -//! Recomputes the energy of every hit based on quenching factor for each particle and volume +//! Recomputes visible hit energies using simple Lindhard and Birks quenching models class TRestGeant4QuenchingProcess : public TRestEventProcess { private: /// A pointer to the specific TRestGeant4Event input @@ -42,6 +46,17 @@ class TRestGeant4QuenchingProcess : public TRestEventProcess { std::set fUserVolumeExpressions; std::set fVolumes; + std::map fUserVolumeModels; + std::map fUserVolumeBirksConstants; + std::map fUserVolumeBirksFallbackStepLengths; + + std::map fVolumeModels; + std::map fVolumeBirksConstants; + std::map fVolumeBirksFallbackStepLengths; + + Double_t fBirksConstant = 0.000126; // mm/keV, equivalent to 0.126 mm/MeV + Double_t fBirksFallbackStepLength = 0.5; // mm + Bool_t fApplyToHitEnergies = true; void Initialize() override; void InitFromConfigFile() override; @@ -72,6 +87,6 @@ class TRestGeant4QuenchingProcess : public TRestEventProcess { explicit TRestGeant4QuenchingProcess(const char* configFilename); ~TRestGeant4QuenchingProcess() override; - ClassDefOverride(TRestGeant4QuenchingProcess, 2); + ClassDefOverride(TRestGeant4QuenchingProcess, 3); }; #endif diff --git a/macros/REST_Geant4_ViewEvent.C b/macros/REST_Geant4_ViewEvent.C index f187941..afb2493 100644 --- a/macros/REST_Geant4_ViewEvent.C +++ b/macros/REST_Geant4_ViewEvent.C @@ -1,8 +1,8 @@ +#include + #include "TRestBrowser.h" #include "TRestTask.h" -#include - #ifndef RestTask_ViewG4Events #define RestTask_ViewG4Events diff --git a/src/TRestGeant4QuenchingProcess.cxx b/src/TRestGeant4QuenchingProcess.cxx index 2a3bcd8..bbfbbf7 100644 --- a/src/TRestGeant4QuenchingProcess.cxx +++ b/src/TRestGeant4QuenchingProcess.cxx @@ -21,24 +21,19 @@ *************************************************************************/ ////////////////////////////////////////////////////////////////////////// -/// The TRestGeant4QuenchingProcess modifies the energy deposits of the simulation by multiplying by a -/// user-defined quenching factor. This quenching factor is specified by volume and by particle and has the -/// value 1.0 by default (no change). Physical volumes, logical volumes, or expressions (to be matched to -/// either physical or logical) can be used for volumes. Example usage: \code +/// The TRestGeant4QuenchingProcess converts Geant4 deposited energies into a simple visible-energy estimate. +/// It supports a Lindhard-style nuclear recoil model and a Birks-law scintillator model. By default the model +/// is inferred from the volume name: scintillator/VETO-like volumes use Birks and the rest use Lindhard. +/// Physical volumes, logical volumes, or expressions (to be matched to either physical or logical) can be +/// used. Example usage: \code ///< TRestGeant4QuenchingProcess> /// -/// -/// -/// -/// -/// -/// -/// +/// +/// +/// /// -/// -/// -/// -/// +/// +/// /// /// /// \endcode @@ -46,10 +41,108 @@ #include "TRestGeant4QuenchingProcess.h" +#include + +#include +#include +#include +#include + using namespace std; ClassImp(TRestGeant4QuenchingProcess); +namespace { +string ToLower(string value) { + transform(value.begin(), value.end(), value.begin(), + [](unsigned char c) { return static_cast(tolower(c)); }); + return value; +} + +string NormalizeQuenchingModel(const string& model) { + const auto normalized = ToLower(model); + if (normalized == "auto" || normalized == "lindhard" || normalized == "birks") { + return normalized; + } + + RESTWarning << "TRestGeant4QuenchingProcess: Unknown quenching model '" << model + << "'. Falling back to auto." << RESTendl; + return "auto"; +} + +string InferQuenchingModel(const string& requestedModel, const string& userVolumeExpression, + const string& volumeName) { + if (requestedModel != "auto") { + return requestedModel; + } + + const auto haystack = ToLower(userVolumeExpression + " " + volumeName); + if (haystack.find("scintillator") != string::npos || haystack.find("veto") != string::npos) { + return "birks"; + } + return "lindhard"; +} + +double Clamp01(double value) { + if (value < 0) { + return 0; + } + if (value > 1) { + return 1; + } + return value; +} + +double LindhardQuenchingFactor(double recoilEnergy, int A, int Z) { + if (recoilEnergy <= 0 || A <= 0 || Z <= 0) { + return 1.0; + } + + const double gamma = 11.5 * recoilEnergy * TMath::Power(Z, -7.0 / 3.0); + const double g = 3 * TMath::Power(gamma, 0.15) + 0.7 * TMath::Power(gamma, 0.6) + gamma; + const double k = 0.133 * TMath::Power(Z, 2.0 / 3.0) * TMath::Power(A, -1.0 / 2.0); + return Clamp01(k * g / (1 + k * g)); +} + +bool IsNeutralParticleWithoutBirksQuenching(const string& particleName) { + const auto name = ToLower(particleName); + return name == "gamma" || name == "opticalphoton" || name == "geantino" || name == "chargedgeantino" || + name.find("neutrino") != string::npos; +} + +double EstimateStepLength(const TRestGeant4Hits& hits, int hitIndex, double fallbackStepLength) { + const auto nHits = static_cast(hits.GetNumberOfHits()); + if (nHits <= 1) { + return fallbackStepLength; + } + + const auto volumeId = hits.GetVolumeId(hitIndex); + const auto position = hits.GetPosition(hitIndex); + double stepLength = numeric_limits::max(); + + if (hitIndex > 0 && hits.GetVolumeId(hitIndex - 1) == volumeId) { + stepLength = min(stepLength, (position - hits.GetPosition(hitIndex - 1)).Mag()); + } + if (hitIndex + 1 < nHits && hits.GetVolumeId(hitIndex + 1) == volumeId) { + stepLength = min(stepLength, (hits.GetPosition(hitIndex + 1) - position).Mag()); + } + + if (!std::isfinite(stepLength) || stepLength <= 0 || stepLength == numeric_limits::max()) { + return fallbackStepLength; + } + return stepLength; +} + +double BirksQuenchingFactor(double energy, double stepLength, double birksConstant) { + if (energy <= 0 || stepLength <= 0 || birksConstant <= 0) { + return 1.0; + } + + const double dEdx = energy / stepLength; + return Clamp01(1.0 / (1.0 + birksConstant * dEdx)); +} +} // namespace + /////////////////////////////////////////////// /// \brief Default constructor /// @@ -96,6 +189,10 @@ void TRestGeant4QuenchingProcess::Initialize() { fInputG4Event = nullptr; fOutputG4Event = new TRestGeant4Event(); + + fBirksConstant = 0.000126; // 0.126 mm/MeV in REST units, mm/keV + fBirksFallbackStepLength = 0.5; // mm + fApplyToHitEnergies = true; } /////////////////////////////////////////////// @@ -117,6 +214,15 @@ void TRestGeant4QuenchingProcess::LoadConfig(const string& configFilename, const } void TRestGeant4QuenchingProcess::InitFromConfigFile() { + fUserVolumeExpressions.clear(); + fUserVolumeModels.clear(); + fUserVolumeBirksConstants.clear(); + fUserVolumeBirksFallbackStepLengths.clear(); + + fApplyToHitEnergies = StringToBool(GetParameter("applyToHitEnergies", fApplyToHitEnergies)); + fBirksConstant = GetDblParameterWithUnits("birksConstant", fBirksConstant); + fBirksFallbackStepLength = GetDblParameterWithUnits("birksFallbackStepLength", fBirksFallbackStepLength); + TiXmlElement* volumeElement = GetElement("volume"); while (volumeElement) { const string volumeName = GetParameter("name", volumeElement, ""); @@ -125,6 +231,12 @@ void TRestGeant4QuenchingProcess::InitFromConfigFile() { exit(1); } fUserVolumeExpressions.insert(volumeName); + fUserVolumeModels[volumeName] = + NormalizeQuenchingModel(GetParameter("model", volumeElement, string("auto"))); + fUserVolumeBirksConstants[volumeName] = + GetDblParameterWithUnits("birksConstant", volumeElement, fBirksConstant); + fUserVolumeBirksFallbackStepLengths[volumeName] = + GetDblParameterWithUnits("birksFallbackStepLength", volumeElement, fBirksFallbackStepLength); volumeElement = GetNextElement(volumeElement); } } @@ -138,6 +250,11 @@ void TRestGeant4QuenchingProcess::InitProcess() { exit(1); } + fVolumes.clear(); + fVolumeModels.clear(); + fVolumeBirksConstants.clear(); + fVolumeBirksFallbackStepLengths.clear(); + const auto geometryInfo = fGeant4Metadata->GetGeant4GeometryInfo(); // check all the user volume expressions are valid and correspond to at least a volume for (const auto& userVolume : fUserVolumeExpressions) { @@ -145,40 +262,41 @@ void TRestGeant4QuenchingProcess::InitProcess() { for (const auto& volume : geometryInfo.GetAllPhysicalVolumesMatchingExpression(userVolume)) { physicalVolumes.insert(volume.Data()); } - if (!physicalVolumes.empty()) { - continue; - } - // maybe it refers to a logical volume - for (const auto& logicalVolume : geometryInfo.GetAllLogicalVolumesMatchingExpression(userVolume)) { - for (const auto& volume : geometryInfo.GetAllPhysicalVolumesFromLogical(logicalVolume.Data())) { - physicalVolumes.insert(volume.Data()); + if (physicalVolumes.empty()) { + // maybe it refers to a logical volume + for (const auto& logicalVolume : + geometryInfo.GetAllLogicalVolumesMatchingExpression(userVolume)) { + for (const auto& volume : + geometryInfo.GetAllPhysicalVolumesFromLogical(logicalVolume.Data())) { + physicalVolumes.insert(volume.Data()); + } } } if (physicalVolumes.empty()) { RESTWarning << "TRestGeant4QuenchingProcess: No volume found matching expression: " << userVolume << RESTendl; + continue; } for (const auto& physicalVolume : physicalVolumes) { const auto volumeName = geometryInfo.GetAlternativeNameFromGeant4PhysicalName(physicalVolume); - fVolumes.insert(volumeName.Data()); + const string volumeNameString = volumeName.Data(); + fVolumes.insert(volumeNameString); + fVolumeModels[volumeNameString] = + InferQuenchingModel(fUserVolumeModels[userVolume], userVolume, volumeNameString); + fVolumeBirksConstants[volumeNameString] = fUserVolumeBirksConstants[userVolume]; + fVolumeBirksFallbackStepLengths[volumeNameString] = + fUserVolumeBirksFallbackStepLengths[userVolume]; } } RESTDebug << "TRestGeant4QuenchingProcess initialized with volumes" << RESTendl; for (const auto& volume : fVolumes) { - RESTDebug << " " << volume << RESTendl; + RESTDebug << " " << volume << " (" << fVolumeModels[volume] << ")" << RESTendl; } } -double QuenchingFactor(double recoilEnergy, int A, int Z) { - // Lindhard formula - double gamma = 11.5 * recoilEnergy * TMath::Power(Z, -7.0 / 3.0); - double g = 3 * TMath::Power(gamma, 0.15) + 0.7 * TMath::Power(gamma, 0.6) + gamma; - double k = 0.133 * TMath::Power(Z, 2.0 / 3.0) * TMath::Power(A, -1.0 / 2.0); - return k * g / (1 + k * g); -} /////////////////////////////////////////////// /// \brief The main processing event function /// @@ -186,51 +304,80 @@ TRestEvent* TRestGeant4QuenchingProcess::ProcessEvent(TRestEvent* inputEvent) { fInputG4Event = (TRestGeant4Event*)inputEvent; *fOutputG4Event = *((TRestGeant4Event*)inputEvent); - const double sensitiveVolumeEnergyBefore = fOutputG4Event->GetSensitiveVolumeEnergy(); + const string sensitiveVolumeName = fGeant4Metadata->GetSensitiveVolume().Data(); + const double sensitiveVolumeEnergyBefore = fOutputG4Event->GetEnergyInVolume(sensitiveVolumeName); fOutputG4Event->InitializeReferences(GetRunInfo()); fOutputG4Event->fEnergyInVolumePerParticlePerProcess.clear(); + fOutputG4Event->fTotalDepositedEnergy = 0; + fOutputG4Event->SetSensitiveVolumeEnergy(0); + fill(fOutputG4Event->fVolumeDepositedEnergy.begin(), fOutputG4Event->fVolumeDepositedEnergy.end(), 0.0); + + bool missingHadronicInfoWarningPrinted = false; // loop over all tracks for (int trackIndex = 0; trackIndex < int(fOutputG4Event->GetNumberOfTracks()); trackIndex++) { // get the track TRestGeant4Track* track = fOutputG4Event->GetTrackPointer(trackIndex); - const auto& particleName = track->GetParticleName(); + const string particleName = track->GetParticleName().Data(); auto hits = track->GetHitsPointer(); - if (!hits->GetHadronicOk()) { - cerr << "TRestGeant4QuenchingProcess: Hadronic information not available. Use the " - "'storeHadronicTargetInfo' parameter in the restG4 configuration" - << endl; - exit(1); - } auto& energy = hits->GetEnergyRef(); for (int hitIndex = 0; hitIndex < int(hits->GetNumberOfHits()); hitIndex++) { - const auto& volumeName = hits->GetVolumeName(hitIndex); - - const string isotopeName = hits->GetHadronicTargetIsotopeName(hitIndex); - const int isotopeA = hits->GetHadronicTargetIsotopeA(hitIndex); - const int isotopeZ = hits->GetHadronicTargetIsotopeZ(hitIndex); + const string volumeName = hits->GetVolumeName(hitIndex).Data(); + const double depositedEnergy = hits->GetEnergy(hitIndex); + double quenchingFactor = 1.0; - double recoilEnergy = hits->GetEnergy(hitIndex); + if (fVolumes.count(volumeName) && depositedEnergy > 0) { + const auto model = + fVolumeModels.count(volumeName) ? fVolumeModels[volumeName] : string("lindhard"); + if (model == "birks") { + if (!IsNeutralParticleWithoutBirksQuenching(particleName)) { + const double stepLength = + EstimateStepLength(*hits, hitIndex, fVolumeBirksFallbackStepLengths[volumeName]); + quenchingFactor = BirksQuenchingFactor(depositedEnergy, stepLength, + fVolumeBirksConstants[volumeName]); + } + } else if (model == "lindhard") { + if (hits->GetHadronicOk()) { + const string isotopeName = hits->GetHadronicTargetIsotopeName(hitIndex); + const int isotopeA = hits->GetHadronicTargetIsotopeA(hitIndex); + const int isotopeZ = hits->GetHadronicTargetIsotopeZ(hitIndex); + if (!isotopeName.empty()) { + quenchingFactor = LindhardQuenchingFactor(depositedEnergy, isotopeA, isotopeZ); + } + } else if (!missingHadronicInfoWarningPrinted) { + RESTWarning + << "TRestGeant4QuenchingProcess: Lindhard quenching requested but " + "hadronic target information is not available. Unquenched hit energies " + "will be used for those hits. Enable storeHadronicTargetInfo in restG4 " + "to apply this model." + << RESTendl; + missingHadronicInfoWarningPrinted = true; + } + } + } - double quenchingFactor = 1.0; - if (fVolumes.count(volumeName.Data()) && recoilEnergy > 0 && !isotopeName.empty()) { - quenchingFactor = QuenchingFactor(recoilEnergy, isotopeA, isotopeZ); + const auto visibleEnergy = depositedEnergy * quenchingFactor; + if (fApplyToHitEnergies) { + energy[hitIndex] = visibleEnergy; } const auto processName = hits->GetProcessName(hitIndex); - if (energy[hitIndex] > 0) { - fOutputG4Event->fEnergyInVolumePerParticlePerProcess[volumeName.Data()][particleName.Data()] - [processName.Data()] += - energy[hitIndex] * quenchingFactor; + if (visibleEnergy > 0) { + fOutputG4Event->AddEnergyInVolumeForParticleForProcess(visibleEnergy, volumeName, + particleName, processName.Data()); + const auto volumeId = hits->GetVolumeId(hitIndex); + if (volumeId >= 0 && volumeId < int(fOutputG4Event->fVolumeDepositedEnergy.size())) { + fOutputG4Event->fVolumeDepositedEnergy[volumeId] += visibleEnergy; + } } } } - const double sensitiveVolumeEnergyAfter = - fOutputG4Event->GetEnergyInVolume(fGeant4Metadata->GetSensitiveVolume().Data()); + const double sensitiveVolumeEnergyAfter = fOutputG4Event->GetEnergyInVolume(sensitiveVolumeName); + fOutputG4Event->SetSensitiveVolumeEnergy(sensitiveVolumeEnergyAfter); bool sensitiveQuenched = TMath::Abs(sensitiveVolumeEnergyAfter - sensitiveVolumeEnergyBefore) > 1e-2; SetObservableValue("sensitiveQuenched", sensitiveQuenched); @@ -244,8 +391,11 @@ void TRestGeant4QuenchingProcess::EndProcess() {} void TRestGeant4QuenchingProcess::PrintMetadata() { BeginPrintProcess(); + RESTMetadata << "Apply to hit energies: " << (fApplyToHitEnergies ? "true" : "false") << RESTendl; + RESTMetadata << "Default Birks constant: " << fBirksConstant << " mm/keV" << RESTendl; + RESTMetadata << "Default Birks fallback step length: " << fBirksFallbackStepLength << " mm" << RESTendl; for (auto const& volume : fVolumes) { - RESTMetadata << "Volume: " << volume << RESTendl; + RESTMetadata << "Volume: " << volume << " model=" << fVolumeModels[volume] << RESTendl; } EndPrintProcess(); } diff --git a/test/files/TRestGeant4QuenchingProcessExample.rml b/test/files/TRestGeant4QuenchingProcessExample.rml index f3ca880..45d264b 100644 --- a/test/files/TRestGeant4QuenchingProcessExample.rml +++ b/test/files/TRestGeant4QuenchingProcessExample.rml @@ -4,18 +4,12 @@ - - - - - - - + + + - - - - + +