Skip to content

Commit

Permalink
Merge pull request #121 from rest-for-physics/lobis-quenching
Browse files Browse the repository at this point in the history 
Quenching
  • Loading branch information
@lobis
lobis authored Mar 7, 2024
2 parents 4df11fb + 420fe4f commit b67ace2
Show file tree
Hide file tree
Showing 6 changed files with 72 additions and 91 deletions.
11 changes: 10 additions & 1 deletion inc/TRestGeant4Hits.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -35,6 +35,10 @@ class TRestGeant4Hits : public TRestHits {
std::vector<Float_t> fKineticEnergy = {};
std::vector<TVector3> fMomentumDirection = {};

std::vector<std::string> fHadronicTargetIsotopeName = {};
std::vector<int> fHadronicTargetIsotopeA = {};
std::vector<int> fHadronicTargetIsotopeZ = {};

TRestGeant4Track* fTrack = nullptr; //!
TRestGeant4Event* fEvent = nullptr; //!

Expand All @@ -58,6 +62,11 @@ class TRestGeant4Hits : public TRestHits {
inline Int_t GetHitVolume(size_t n) const { return GetVolumeId(n); }
TString GetVolumeName(size_t n) const;

inline bool GetHadronicOk() const { return fHadronicTargetIsotopeName.size() > 0; }
inline std::string GetHadronicTargetIsotopeName(size_t n) const { return fHadronicTargetIsotopeName[n]; }
inline int GetHadronicTargetIsotopeA(size_t n) const { return fHadronicTargetIsotopeA[n]; }
inline int GetHadronicTargetIsotopeZ(size_t n) const { return fHadronicTargetIsotopeZ[n]; }

void RemoveG4Hits();

inline Double_t GetKineticEnergy(size_t n) const { return fKineticEnergy[n]; }
Expand All @@ -78,7 +87,7 @@ class TRestGeant4Hits : public TRestHits {
// Destructor
virtual ~TRestGeant4Hits();

ClassDef(TRestGeant4Hits, 7); // REST event superclass
ClassDef(TRestGeant4Hits, 8); // REST event superclass

// restG4
public:
Expand Down
7 changes: 6 additions & 1 deletion inc/TRestGeant4Metadata.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -141,6 +141,9 @@ class TRestGeant4Metadata : public TRestMetadata {
/// sensitive volume (used for debugging purposes). It is set to 'false' by default.
Bool_t fSaveAllEvents = false;

/// \brief Option to store target isotope information on hadronic processes (disabled by default)
Bool_t fStoreHadronicTargetInfo = false;

/// \brief Sets all volume as active without having to explicitly list them.
Bool_t fActivateAllVolumes = false; //!

Expand Down Expand Up @@ -190,6 +193,8 @@ class TRestGeant4Metadata : public TRestMetadata {

size_t GetGeant4VersionMajor() const;

bool GetStoreHadronicTargetInfo() const { return fStoreHadronicTargetInfo; }

/// Returns the local path to the GDML geometry
inline TString GetGeometryPath() const { return fGeometryPath; }

Expand Down Expand Up @@ -393,7 +398,7 @@ class TRestGeant4Metadata : public TRestMetadata {
TRestGeant4Metadata(const TRestGeant4Metadata& metadata);
TRestGeant4Metadata& operator=(const TRestGeant4Metadata& metadata);

ClassDefOverride(TRestGeant4Metadata, 15);
ClassDefOverride(TRestGeant4Metadata, 16);

// Allow modification of otherwise inaccessible / immutable members that shouldn't be modified by the user
friend class SteppingAction;
Expand Down
15 changes: 5 additions & 10 deletions inc/TRestGeant4QuenchingProcess.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -40,22 +40,17 @@ class TRestGeant4QuenchingProcess : public TRestEventProcess {
/// A pointer to the simulation metadata information accessible to TRestRun
TRestGeant4Metadata* fGeant4Metadata{}; //!

/// It stores the quenching factor for each particle for each user volume expression
std::map<std::string, std::map<std::string, float_t>> fUserVolumeParticleQuenchingFactor = {};

/// It stores the quenching factor for each particle for each volume of the simulation
std::map<std::string, std::map<std::string, float_t>> fVolumeParticleQuenchingFactor = {};
std::set<std::string> fUserVolumeExpressions;
std::set<std::string> fVolumes;

void Initialize() override;
void InitFromConfigFile() override;
void LoadDefaultConfig();

public:
std::vector<std::string> GetUserVolumeExpressions() const;
float_t GetQuenchingFactorForVolumeAndParticle(const std::string& volumeName,
const std::string& particleName) const;
std::set<std::string> GetUserVolumeExpressions() const;
std::set<std::string> GetVolumes() const;

//
RESTValue GetInputEvent() const override { return fInputG4Event; }
RESTValue GetOutputEvent() const override { return fOutputG4Event; }

Expand All @@ -77,6 +72,6 @@ class TRestGeant4QuenchingProcess : public TRestEventProcess {
explicit TRestGeant4QuenchingProcess(const char* configFilename);
~TRestGeant4QuenchingProcess() override;

ClassDefOverride(TRestGeant4QuenchingProcess, 1);
ClassDefOverride(TRestGeant4QuenchingProcess, 2);
};
#endif
2 changes: 2 additions & 0 deletions src/TRestGeant4Metadata.cxx
View file
Original file line number Diff line number Diff line change
Expand Up @@ -816,6 +816,8 @@ void TRestGeant4Metadata::InitFromConfigFile() {

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") {
Expand Down
120 changes: 48 additions & 72 deletions src/TRestGeant4QuenchingProcess.cxx
View file
Original file line number Diff line number Diff line change
Expand Up @@ -119,27 +119,12 @@ void TRestGeant4QuenchingProcess::LoadConfig(const string& configFilename, const
void TRestGeant4QuenchingProcess::InitFromConfigFile() {
TiXmlElement* volumeElement = GetElement("volume");
while (volumeElement) {
TiXmlElement* particleElement = GetElement("particle", volumeElement);
const string volumeName = GetParameter("name", volumeElement, "");
if (volumeName.empty()) {
cerr << "TRestGeant4QuenchingProcess: No volume expression specified" << endl;
exit(1);
}
while (particleElement) {
const string particleName = GetParameter("name", particleElement, "");
if (particleName.empty()) {
cerr << "TRestGeant4QuenchingProcess: No particle name specified" << endl;
exit(1);
}
const auto quenchingFactor = stof(GetParameter("quenchingFactor", particleElement, -1.0));
// if no value is specified, give error
if (quenchingFactor < 0.0 || quenchingFactor > 1.0) {
cerr << "TRestGeant4QuenchingProcess: Quenching factor must be between 0 and 1" << endl;
exit(1);
}
fUserVolumeParticleQuenchingFactor[volumeName][particleName] = quenchingFactor;
particleElement = GetNextElement(particleElement);
}
fUserVolumeExpressions.insert(volumeName);
volumeElement = GetNextElement(volumeElement);
}
}
Expand All @@ -155,7 +140,7 @@ void TRestGeant4QuenchingProcess::InitProcess() {

const auto geometryInfo = fGeant4Metadata->GetGeant4GeometryInfo();
// check all the user volume expressions are valid and correspond to at least a volume
for (const auto& [userVolume, particleToQuenchingMap] : fUserVolumeParticleQuenchingFactor) {
for (const auto& userVolume : fUserVolumeExpressions) {
set<string> physicalVolumes = {};
for (const auto& volume : geometryInfo.GetAllPhysicalVolumesMatchingExpression(userVolume)) {
physicalVolumes.insert(volume.Data());
Expand All @@ -177,26 +162,32 @@ void TRestGeant4QuenchingProcess::InitProcess() {

for (const auto& physicalVolume : physicalVolumes) {
const auto volumeName = geometryInfo.GetAlternativeNameFromGeant4PhysicalName(physicalVolume);
fVolumeParticleQuenchingFactor[volumeName.Data()] = particleToQuenchingMap;
fVolumes.insert(volumeName.Data());
}
}

RESTDebug << "TRestGeant4QuenchingProcess initialized with volumes" << RESTendl;
for (const auto& [volume, particleToQuenchingMap] : fVolumeParticleQuenchingFactor) {
RESTDebug << volume << RESTendl;
for (const auto& [particle, quenchingFactor] : particleToQuenchingMap) {
RESTDebug << "\t" << particle << " " << quenchingFactor << RESTendl;
}
cout << "TRestGeant4QuenchingProcess initialized with volumes" << endl;
for (const auto& volume : fVolumes) {
cout << " " << volume << endl;
}
}

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
///
TRestEvent* TRestGeant4QuenchingProcess::ProcessEvent(TRestEvent* inputEvent) {
fInputG4Event = (TRestGeant4Event*)inputEvent;
*fOutputG4Event = *((TRestGeant4Event*)inputEvent);

const double sensitiveVolumeEnergyBefore = fOutputG4Event->GetSensitiveVolumeEnergy();

fOutputG4Event->InitializeReferences(GetRunInfo());
fOutputG4Event->fEnergyInVolumePerParticlePerProcess.clear();

Expand All @@ -207,76 +198,61 @@ TRestEvent* TRestGeant4QuenchingProcess::ProcessEvent(TRestEvent* inputEvent) {
const auto& particleName = track->GetParticleName();

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 float_t quenchingFactor =
GetQuenchingFactorForVolumeAndParticle(volumeName.Data(), particleName.Data());
// default is 1.0, so no quenching
energy[hitIndex] *= quenchingFactor;

// cout << "TRestGeant4QuenchingProcess: Quenching " << particleName << " in " << volumeName
// << " by " << quenchingFactor << endl;
const string isotopeName = hits->GetHadronicTargetIsotopeName(hitIndex);
const int isotopeA = hits->GetHadronicTargetIsotopeA(hitIndex);
const int isotopeZ = hits->GetHadronicTargetIsotopeZ(hitIndex);

double recoilEnergy = hits->GetEnergy(hitIndex);

double quenchingFactor = 1.0;
if (fVolumes.count(volumeName.Data()) && recoilEnergy > 0 && !isotopeName.empty()) {
quenchingFactor = QuenchingFactor(recoilEnergy, isotopeA, isotopeZ);
}

const string processName = hits->GetProcessName(hitIndex).Data();
const auto processName = hits->GetProcessName(hitIndex);

if (energy[hitIndex] > 0) {
fOutputG4Event->fEnergyInVolumePerParticlePerProcess[volumeName.Data()][particleName.Data()]
[processName] += energy[hitIndex];
[processName.Data()] +=
energy[hitIndex] * quenchingFactor;
}
}
}

// Update the stores for energy in volumes (this should be automatic and not duplicated)
const double sensitiveVolumeEnergyAfter =
fOutputG4Event->GetEnergyInVolume(fGeant4Metadata->GetSensitiveVolume().Data());

bool sensitiveQuenched = TMath::Abs(sensitiveVolumeEnergyAfter - sensitiveVolumeEnergyBefore) > 1e-2;
SetObservableValue("sensitiveQuenched", sensitiveQuenched);
SetObservableValue("sensitiveVolumeEnergyBefore", sensitiveVolumeEnergyBefore);
SetObservableValue("sensitiveVolumeEnergyAfter", sensitiveVolumeEnergyAfter);

return fOutputG4Event;
}

///////////////////////////////////////////////
/// \brief Function to include required actions after all events have been processed.
void TRestGeant4QuenchingProcess::EndProcess() {}

vector<string> TRestGeant4QuenchingProcess::GetUserVolumeExpressions() const {
vector<string> userVolumeExpressions;
userVolumeExpressions.reserve(fUserVolumeParticleQuenchingFactor.size());
for (auto const& x : fUserVolumeParticleQuenchingFactor) {
userVolumeExpressions.push_back(x.first);
}
return userVolumeExpressions;
}

float_t TRestGeant4QuenchingProcess::GetQuenchingFactorForVolumeAndParticle(
const string& volumeName, const string& particleName) const {
float_t quenchingFactor = 1.0;
// check if the volume is in the map
if (fVolumeParticleQuenchingFactor.find(volumeName) != fVolumeParticleQuenchingFactor.end()) {
// check if the particle is in the map
if (fVolumeParticleQuenchingFactor.at(volumeName).find(particleName) !=
fVolumeParticleQuenchingFactor.at(volumeName).end()) {
quenchingFactor = fVolumeParticleQuenchingFactor.at(volumeName).at(particleName);
}
}
return quenchingFactor;
}

void TRestGeant4QuenchingProcess::PrintMetadata() {
BeginPrintProcess();
cout << "Printing TRestGeant4QuenchingProcess user configuration" << endl;
for (auto const& [volume, particleQuenchingFactorMap] : fUserVolumeParticleQuenchingFactor) {
for (auto const& volume : fVolumes) {
cout << "Volume: " << volume << endl;
for (auto const& [particle, quenchingFactor] : particleQuenchingFactorMap) {
cout << " - Particle: " << particle << " Quenching factor: " << quenchingFactor << endl;
}
}

if (!fVolumeParticleQuenchingFactor.empty()) {
cout << "Printing TRestGeant4QuenchingProcess configuration with actual volumes" << endl;
for (auto const& [volume, particleQuenchingFactorMap] : fVolumeParticleQuenchingFactor) {
cout << "Volume: " << volume << endl;
for (auto const& [particle, quenchingFactor] : particleQuenchingFactorMap) {
cout << " - Particle: " << particle << " Quenching factor: " << quenchingFactor << endl;
}
}
}
EndPrintProcess();
}

std::set<std::string> TRestGeant4QuenchingProcess::GetVolumes() const { return fVolumes; }

std::set<std::string> TRestGeant4QuenchingProcess::GetUserVolumeExpressions() const {
return fUserVolumeExpressions;
}
8 changes: 1 addition & 7 deletions test/src/Geant4QuenchingProcesses.cxx
View file
Original file line number Diff line number Diff line change
Expand Up @@ -10,7 +10,7 @@ using namespace std;

const auto filesPath = fs::path(__FILE__).parent_path().parent_path() / "files";
const auto processRmlFile = filesPath / "TRestGeant4QuenchingProcessExample.rml";
const auto simulationFile = "/tmp/tmp.CXV7dwGPdN/source/packages/restG4/examples/13.IAXO/Neutrons.root";
// const auto simulationFile = "/tmp/tmp.CXV7dwGPdN/source/packages/restG4/examples/13.IAXO/Neutrons.root";

TEST(TRestGeant4QuenchingProcess, TestFiles) {
cout << "Test files path: " << filesPath << endl;
Expand Down Expand Up @@ -41,12 +41,6 @@ TEST(TRestGeant4QuenchingProcess, Default) {

// verify default quenching factor is 1.0
vector<string> randomVolumes = {"scintillatorVolume", "gasVolume", "shieldingVolume"};
vector<string> particles = {"gamma", "e-", "e+"};
for (const auto& volume : randomVolumes) {
for (const auto& particle : particles) {
EXPECT_TRUE(process.GetQuenchingFactorForVolumeAndParticle(volume, particle) == 1.0);
}
}
}

TEST(TRestGeant4QuenchingProcess, FromRml) {
Expand Down

0 comments on commit b67ace2

Please sign in to comment.