From 9024f778cdbdfb7e24402e4550cfc319584ac960 Mon Sep 17 00:00:00 2001
From: "audrey.chatillon" <audrey.chatillon@gmail.com>
Date: Tue, 21 Jan 2025 16:57:44 +0100
Subject: [PATCH] Update DecayLaw in EventGeneratorIsotropic
---
.../EventGenerator/EventGeneratorIsotropic.cc | 23 ++++++++++++-------
.../EventGenerator/EventGeneratorIsotropic.hh | 2 +-
2 files changed, 16 insertions(+), 9 deletions(-)
diff --git a/NPSimulation/EventGenerator/EventGeneratorIsotropic.cc b/NPSimulation/EventGenerator/EventGeneratorIsotropic.cc
index 0788d3c16..beb9a1b5d 100644
--- a/NPSimulation/EventGenerator/EventGeneratorIsotropic.cc
+++ b/NPSimulation/EventGenerator/EventGeneratorIsotropic.cc
@@ -129,10 +129,10 @@ void EventGeneratorIsotropic::ReadConfiguration(NPL::InputParser parser){
if(blocks[i]->HasToken("SigmaZ"))
it->m_SigmaZ=blocks[i]->GetDouble("SigmaZ","mm");
if(blocks[i]->HasToken("Multiplicity"))
- it->m_Multiplicty=blocks[i]->GetVectorInt("Multiplicity");
+ it->m_Multiplicity=blocks[i]->GetVectorInt("Multiplicity");
if(blocks[i]->HasToken("DecayLaw"))
it->m_DecayLaw=blocks[i]->GetString("DecayLaw");
- if(blocks[i]->HasToken("Activity_Bq"))
+ if(blocks[i]->HasToken("ActivityBq"))
it->m_ActivityBq=blocks[i]->GetDouble("ActivityBq","Bq");
}
else{
@@ -143,8 +143,8 @@ void EventGeneratorIsotropic::ReadConfiguration(NPL::InputParser parser){
for(auto& par : m_Parameters) {
if(par.m_ExcitationEnergy.size()==0)
par.m_ExcitationEnergy.resize(par.m_particleName.size(),0);
- if(par.m_Multiplicty.size()==0)
- par.m_Multiplicty.resize(par.m_particleName.size(),1);
+ if(par.m_Multiplicity.size()==0)
+ par.m_Multiplicity.resize(par.m_particleName.size(),1);
if(par.m_EnergyDistribution!="flat" && par.m_EnergyDistribution!="FromHisto"){
if(par.m_EnergyDistribution=="Watt"){
@@ -162,9 +162,14 @@ void EventGeneratorIsotropic::ReadConfiguration(NPL::InputParser parser){
void EventGeneratorIsotropic::GenerateEvent(G4Event*){
for(auto& par : m_Parameters) {
- double delta_t = 0;
+ double delta_t = 0.;
+ bool DecayOn = true;
for(unsigned int p=0; p<par.m_particleName.size(); p++){
- for(int i=0; i<par.m_Multiplicty[p]; i++){
+ double step_t = 1.e-2 / par.m_ActivityBq; // [s]
+ if(par.m_DecayLaw=="on"){
+ par.m_Multiplicity[p] = ((int)(50.e-9 / step_t) == 0) ? 1 : (int)(50.e-9 / step_t);
+ }
+ for(int i=0; i<par.m_Multiplicity[p]; i++){
par.m_particle=NULL;
if(par.m_particle==NULL){
@@ -178,11 +183,13 @@ void EventGeneratorIsotropic::GenerateEvent(G4Event*){
par.m_particle = G4ParticleTable::GetParticleTable()->GetIonTable()->GetIon(N->GetZ(), N->GetA(),par.m_ExcitationEnergy[p]);
delete N;
if(par.m_DecayLaw == "on" && i>0){
- delta_t += G4RandExponential::shoot(1./par.m_ActivityBq) ;
+ delta_t = i*step_t;
+ DecayOn = !(par.m_DecayLaw == "on" && RandFlat::shoot() > 1.e-2);
}
}
}
- if(delta_t > 50.e-9) continue; //to do to parametrize (time_window)
+ if(!DecayOn) continue;
+ //cout << " Process decay at delta_t = " << delta_t << endl;
G4double cos_theta_min = cos(par.m_HalfOpenAngleMin);
G4double cos_theta_max = cos(par.m_HalfOpenAngleMax);
G4double cos_theta = cos_theta_min + (cos_theta_max - cos_theta_min) * RandFlat::shoot();
diff --git a/NPSimulation/EventGenerator/EventGeneratorIsotropic.hh b/NPSimulation/EventGenerator/EventGeneratorIsotropic.hh
index e15e00f2d..7ed442381 100644
--- a/NPSimulation/EventGenerator/EventGeneratorIsotropic.hh
+++ b/NPSimulation/EventGenerator/EventGeneratorIsotropic.hh
@@ -71,7 +71,7 @@ private: // Source parameter from input file
G4ParticleDefinition* m_particle ; // Kind of particle to shoot isotropically
vector<G4double> m_ExcitationEnergy ; // Excitation energy of the emitted particle
vector<string> m_particleName ;
- vector<G4int> m_Multiplicty ;
+ vector<G4int> m_Multiplicity ;
TString m_DecayLaw ; // [on]: mult follows the decay law controlled by T1/2 [off]: fix mult
G4double m_ActivityBq ; // if m_DecayLaw = "on" gives here the Activity in Bq of the sample
};
--
GitLab