Update DecayLaw in EventGeneratorIsotropic

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();

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
   };