Geant4 10.7.0
Toolkit for the simulation of the passage of particles through matter
Loading...
Searching...
No Matches
Public Member Functions | Protected Member Functions | List of all members
G4RadioactiveDecay Class Reference

#include <G4RadioactiveDecay.hh>

+ Inheritance diagram for G4RadioactiveDecay:

Public Member Functions

 G4RadioactiveDecay (const G4String &processName="RadioactiveDecay")
 
 ~G4RadioactiveDecay ()
 
virtual void ProcessDescription (std::ostream &outFile) const
 
G4bool IsApplicable (const G4ParticleDefinition &)
 
G4DecayTable * GetDecayTable (const G4ParticleDefinition *)
 
void SelectAVolume (const G4String aVolume)
 
void DeselectAVolume (const G4String aVolume)
 
void SelectAllVolumes ()
 
void DeselectAllVolumes ()
 
void SetDecayBias (G4String filename)
 
void SetHLThreshold (G4double hl)
 
void SetICM (G4bool icm)
 
void SetARM (G4bool arm)
 
void SetSourceTimeProfile (G4String filename)
 
G4bool IsRateTableReady (const G4ParticleDefinition &)
 
void CalculateChainsFromParent (const G4ParticleDefinition &)
 
void GetChainsFromParent (const G4ParticleDefinition &)
 
void SetDecayRate (G4int, G4int, G4double, G4int, std::vector< G4double >, std::vector< G4double >)
 
std::vector< G4RadioactivityTable * > GetTheRadioactivityTables ()
 
G4DecayTable * LoadDecayTable (const G4ParticleDefinition &theParentNucleus)
 
void AddUserDecayDataFile (G4int Z, G4int A, G4String filename)
 
void SetVerboseLevel (G4int value)
 
G4int GetVerboseLevel () const
 
void SetNucleusLimits (G4NucleusLimits theNucleusLimits1)
 
G4NucleusLimits GetNucleusLimits () const
 
void SetAnalogueMonteCarlo (G4bool r)
 
void SetFBeta (G4bool r)
 
G4bool IsAnalogueMonteCarlo ()
 
void SetBRBias (G4bool r)
 
void SetSplitNuclei (G4int r)
 
G4int GetSplitNuclei ()
 
void SetDecayDirection (const G4ThreeVector &theDir)
 
const G4ThreeVector & GetDecayDirection () const
 
void SetDecayHalfAngle (G4double halfAngle=0.*CLHEP::deg)
 
G4double GetDecayHalfAngle () const
 
void SetDecayCollimation (const G4ThreeVector &theDir, G4double halfAngle=0.*CLHEP::deg)
 
void BuildPhysicsTable (const G4ParticleDefinition &)
 
G4VParticleChange * DecayIt (const G4Track &theTrack, const G4Step &theStep)
 
- Public Member Functions inherited from G4VRestDiscreteProcess
 G4VRestDiscreteProcess (const G4String &, G4ProcessType aType=fNotDefined)
 
 G4VRestDiscreteProcess (G4VRestDiscreteProcess &)
 
virtual ~G4VRestDiscreteProcess ()
 
G4VRestDiscreteProcess & operator= (const G4VRestDiscreteProcess &)=delete
 
virtual G4double PostStepGetPhysicalInteractionLength (const G4Track &track, G4double previousStepSize, G4ForceCondition *condition)
 
virtual G4VParticleChange * PostStepDoIt (const G4Track &, const G4Step &)
 
virtual G4double AtRestGetPhysicalInteractionLength (const G4Track &, G4ForceCondition *)
 
virtual G4VParticleChange * AtRestDoIt (const G4Track &, const G4Step &)
 
virtual G4double AlongStepGetPhysicalInteractionLength (const G4Track &, G4double, G4double, G4double &, G4GPILSelection *)
 
virtual G4VParticleChange * AlongStepDoIt (const G4Track &, const G4Step &)
 
- Public Member Functions inherited from G4VProcess
 G4VProcess (const G4String &aName="NoName", G4ProcessType aType=fNotDefined)
 
 G4VProcess (const G4VProcess &right)
 
virtual ~G4VProcess ()
 
G4VProcess & operator= (const G4VProcess &)=delete
 
G4bool operator== (const G4VProcess &right) const
 
G4bool operator!= (const G4VProcess &right) const
 
virtual G4VParticleChange * PostStepDoIt (const G4Track &track, const G4Step &stepData)=0
 
virtual G4VParticleChange * AlongStepDoIt (const G4Track &track, const G4Step &stepData)=0
 
virtual G4VParticleChange * AtRestDoIt (const G4Track &track, const G4Step &stepData)=0
 
virtual G4double AlongStepGetPhysicalInteractionLength (const G4Track &track, G4double previousStepSize, G4double currentMinimumStep, G4double &proposedSafety, G4GPILSelection *selection)=0
 
virtual G4double AtRestGetPhysicalInteractionLength (const G4Track &track, G4ForceCondition *condition)=0
 
virtual G4double PostStepGetPhysicalInteractionLength (const G4Track &track, G4double previousStepSize, G4ForceCondition *condition)=0
 
G4double GetCurrentInteractionLength () const
 
void SetPILfactor (G4double value)
 
G4double GetPILfactor () const
 
G4double AlongStepGPIL (const G4Track &track, G4double previousStepSize, G4double currentMinimumStep, G4double &proposedSafety, G4GPILSelection *selection)
 
G4double AtRestGPIL (const G4Track &track, G4ForceCondition *condition)
 
G4double PostStepGPIL (const G4Track &track, G4double previousStepSize, G4ForceCondition *condition)
 
virtual G4bool IsApplicable (const G4ParticleDefinition &)
 
virtual void BuildPhysicsTable (const G4ParticleDefinition &)
 
virtual void PreparePhysicsTable (const G4ParticleDefinition &)
 
virtual G4bool StorePhysicsTable (const G4ParticleDefinition *, const G4String &, G4bool)
 
virtual G4bool RetrievePhysicsTable (const G4ParticleDefinition *, const G4String &, G4bool)
 
const G4String & GetPhysicsTableFileName (const G4ParticleDefinition *, const G4String &directory, const G4String &tableName, G4bool ascii=false)
 
const G4String & GetProcessName () const
 
G4ProcessType GetProcessType () const
 
void SetProcessType (G4ProcessType)
 
G4int GetProcessSubType () const
 
void SetProcessSubType (G4int)
 
virtual void StartTracking (G4Track *)
 
virtual void EndTracking ()
 
virtual void SetProcessManager (const G4ProcessManager *)
 
virtual const G4ProcessManager * GetProcessManager ()
 
virtual void ResetNumberOfInteractionLengthLeft ()
 
G4double GetNumberOfInteractionLengthLeft () const
 
G4double GetTotalNumberOfInteractionLengthTraversed () const
 
G4bool isAtRestDoItIsEnabled () const
 
G4bool isAlongStepDoItIsEnabled () const
 
G4bool isPostStepDoItIsEnabled () const
 
virtual void DumpInfo () const
 
virtual void ProcessDescription (std::ostream &outfile) const
 
void SetVerboseLevel (G4int value)
 
G4int GetVerboseLevel () const
 
virtual void SetMasterProcess (G4VProcess *masterP)
 
const G4VProcess * GetMasterProcess () const
 
virtual void BuildWorkerPhysicsTable (const G4ParticleDefinition &part)
 
virtual void PrepareWorkerPhysicsTable (const G4ParticleDefinition &)
 

Protected Member Functions

G4DecayProducts * DoDecay (const G4ParticleDefinition &theParticleDef)
 
void CollimateDecay (G4DecayProducts *products)
 
void CollimateDecayProduct (G4DynamicParticle *product)
 
G4ThreeVector ChooseCollimationDirection () const
 
G4double GetMeanFreePath (const G4Track &theTrack, G4double previousStepSize, G4ForceCondition *condition)
 
G4double GetMeanLifeTime (const G4Track &theTrack, G4ForceCondition *condition)
 
G4double ConvolveSourceTimeProfile (const G4double, const G4double)
 
G4double GetDecayTime ()
 
G4int GetDecayTimeBin (const G4double aDecayTime)
 
void AddDeexcitationSpectrumForBiasMode (G4ParticleDefinition *apartDef, G4double weight, G4double currenTime, std::vector< double > &weights_v, std::vector< double > &times_v, std::vector< G4DynamicParticle * > &secondaries_v)
 
virtual G4double GetMeanFreePath (const G4Track &aTrack, G4double previousStepSize, G4ForceCondition *condition)=0
 
virtual G4double GetMeanLifeTime (const G4Track &aTrack, G4ForceCondition *condition)=0
 
- Protected Member Functions inherited from G4VProcess
void SubtractNumberOfInteractionLengthLeft (G4double prevStepSize)
 
void ClearNumberOfInteractionLengthLeft ()
 

Additional Inherited Members

- Static Public Member Functions inherited from G4VProcess
static const G4String & GetProcessTypeName (G4ProcessType)
 
- Protected Attributes inherited from G4VProcess
const G4ProcessManager * aProcessManager = nullptr
 
G4VParticleChange * pParticleChange = nullptr
 
G4ParticleChange aParticleChange
 
G4double theNumberOfInteractionLengthLeft = -1.0
 
G4double currentInteractionLength = -1.0
 
G4double theInitialNumberOfInteractionLength = -1.0
 
G4String theProcessName
 
G4String thePhysicsTableFileName
 
G4ProcessType theProcessType = fNotDefined
 
G4int theProcessSubType = -1
 
G4double thePILfactor = 1.0
 
G4int verboseLevel = 0
 
G4bool enableAtRestDoIt = true
 
G4bool enableAlongStepDoIt = true
 
G4bool enablePostStepDoIt = true
 

Detailed Description

Definition at line 83 of file G4RadioactiveDecay.hh.

Constructor & Destructor Documentation

◆ G4RadioactiveDecay()

G4RadioactiveDecay::G4RadioactiveDecay ( const G4String &  processName = "RadioactiveDecay")

Definition at line 169 of file G4RadioactiveDecay.cc.

170 : G4VRestDiscreteProcess(processName, fDecay), isInitialised(false),
171 forceDecayDirection(0.,0.,0.), forceDecayHalfAngle(0.*deg), dirPath(""),
172 verboseLevel(1)
173{
174#ifdef G4VERBOSE
175 if (GetVerboseLevel() > 1) {
176 G4cout << "G4RadioactiveDecay constructor: processName = " << processName
177 << G4endl;
178 }
179#endif
180
181 G4cout << " G4RadioactiveDecay is deprecated and will be removed in Geant4 version 11. " << G4endl;
182 G4cout << " Please replace it with G4RadioactiveDecayBase if you want the unbiased radioactive deacy process." << G4endl;
183 G4cout << " If you want the general process, with optional biasing, use G4Radioactivation. " << G4endl;
184
185 SetProcessSubType(fRadioactiveDecay);
186
187 theRadioactiveDecaymessenger = new G4RadioactiveDecaymessenger(this);
188 pParticleChange = &fParticleChangeForRadDecay;
189
190 // Set up photon evaporation for use in G4ITDecay
191 photonEvaporation = new G4PhotonEvaporation();
192 // photonEvaporation->SetVerboseLevel(2);
193 photonEvaporation->RDMForced(true);
194 photonEvaporation->SetICM(true);
195
196 // Check data directory
197 char* path_var = std::getenv("G4RADIOACTIVEDATA");
198 if (!path_var) {
199 G4Exception("G4RadioactiveDecay()","HAD_RDM_200",FatalException,
200 "Environment variable G4RADIOACTIVEDATA is not set");
201 } else {
202 dirPath = path_var; // convert to string
203 std::ostringstream os;
204 os << dirPath << "/z1.a3"; // used as a dummy
205 std::ifstream testFile;
206 testFile.open(os.str() );
207 if (!testFile.is_open() )
208 G4Exception("G4RadioactiveDecay()","HAD_RDM_201",FatalException,
209 "Environment variable G4RADIOACTIVEDATA is set, but does not point to correct directory");
210 }
211
212 // Reset the list of user defined data files
213 theUserRadioactiveDataFiles.clear();
214
215 // Instantiate the map of decay tables
216#ifdef G4MULTITHREADED
217 G4AutoLock lk(&G4RadioactiveDecay::radioactiveDecayMutex);
218 NumberOfInstances()++;
219 if(!master_dkmap) master_dkmap = new DecayTableMap;
220#endif
221 dkmap = new DecayTableMap;
222
223 // Apply default values.
224 NSourceBin = 1;
225 SBin[0] = 0.* s;
226 SBin[1] = 1.* s;
227 SProfile[0] = 1.;
228 SProfile[1] = 0.;
229 NDecayBin = 1;
230 DBin[0] = 0. * s ;
231 DBin[1] = 1. * s;
232 DProfile[0] = 1.;
233 DProfile[1] = 0.;
234 decayWindows[0] = 0;
235 G4RadioactivityTable* rTable = new G4RadioactivityTable() ;
236 theRadioactivityTables.push_back(rTable);
237 NSplit = 1;
238 AnalogueMC = true ;
239 FBeta = false ;
240 BRBias = true ;
241 applyICM = true ;
242 applyARM = true ;
243 halflifethreshold = nanosecond;
244
245 // RDM applies to all logical volumes by default
246 isAllVolumesMode = true;
247 SelectAllVolumes();
248 G4HadronicProcessStore::Instance()->RegisterExtraProcess(this);
249}
FatalException
@ FatalException
Definition: G4ExceptionSeverity.hh:61
G4Exception
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *description)
Definition: G4Exception.cc:35
fRadioactiveDecay
@ fRadioactiveDecay
Definition: G4HadronicProcessType.hh:51
fDecay
@ fDecay
Definition: G4ProcessType.hh:45
DecayTableMap
std::map< G4String, G4DecayTable * > DecayTableMap
Definition: G4RadioactiveDecay.hh:80
G4endl
#define G4endl
Definition: G4ios.hh:57
G4cout
G4GLOB_DLL std::ostream G4cout
G4HadronicProcessStore::RegisterExtraProcess
void RegisterExtraProcess(G4VProcess *)
Definition: G4HadronicProcessStore.cc:487
G4HadronicProcessStore::Instance
static G4HadronicProcessStore * Instance()
Definition: G4HadronicProcessStore.cc:70
G4PhotonEvaporation::SetICM
virtual void SetICM(G4bool)
Definition: G4PhotonEvaporation.cc:555
G4PhotonEvaporation::RDMForced
virtual void RDMForced(G4bool)
Definition: G4PhotonEvaporation.cc:560
G4RadioactiveDecay::GetVerboseLevel
G4int GetVerboseLevel() const
Definition: G4RadioactiveDecay.hh:173
G4VProcess::SetProcessSubType
void SetProcessSubType(G4int)
Definition: G4VProcess.hh:406
G4VProcess::pParticleChange
G4VParticleChange * pParticleChange
Definition: G4VProcess.hh:321

◆ ~G4RadioactiveDecay()

G4RadioactiveDecay::~G4RadioactiveDecay ( )

Definition at line 262 of file G4RadioactiveDecay.cc.

263{
264 delete theRadioactiveDecaymessenger;
265 delete photonEvaporation;
266 for (DecayTableMap::iterator i = dkmap->begin(); i != dkmap->end(); i++) {
267 delete i->second;
268 }
269 dkmap->clear();
270 delete dkmap;
271#ifdef G4MULTITHREADED
272 G4AutoLock lk(&G4RadioactiveDecay::radioactiveDecayMutex);
273 --NumberOfInstances();
274 if(NumberOfInstances()==0)
275 {
276 for (DecayTableMap::iterator i = master_dkmap->begin(); i != master_dkmap->end(); i++) {
277 delete i->second;
278 }
279 master_dkmap->clear();
280 delete master_dkmap;
281 }
282#endif
283}

Member Function Documentation

◆ AddDeexcitationSpectrumForBiasMode()

void G4RadioactiveDecay::AddDeexcitationSpectrumForBiasMode ( G4ParticleDefinition *  apartDef,
G4double  weight,
G4double  currenTime,
std::vector< double > &  weights_v,
std::vector< double > &  times_v,
std::vector< G4DynamicParticle * > &  secondaries_v 
)
protected

Definition at line 2186 of file G4RadioactiveDecay.cc.

2191{
2192 G4double elevel = ((const G4Ions*)(apartDef))->GetExcitationEnergy();
2193 G4double life_time = apartDef->GetPDGLifeTime();
2194 G4ITDecay* anITChannel = 0;
2195
2196 while (life_time < halflifethreshold && elevel > 0.) {
2197 anITChannel = new G4ITDecay(apartDef, 100., elevel, elevel, photonEvaporation);
2198 G4DecayProducts* pevap_products = anITChannel->DecayIt(0.);
2199 G4int nb_pevapSecondaries = pevap_products->entries();
2200
2201 G4DynamicParticle* a_pevap_secondary = 0;
2202 G4ParticleDefinition* secDef = 0;
2203 for (G4int ind = 0; ind < nb_pevapSecondaries; ind++) {
2204 a_pevap_secondary = pevap_products->PopProducts();
2205 secDef = a_pevap_secondary->GetDefinition();
2206
2207 if (secDef->GetBaryonNumber() > 4) {
2208 elevel = ((const G4Ions*)(secDef))->GetExcitationEnergy();
2209 life_time = secDef->GetPDGLifeTime();
2210 apartDef = secDef;
2211 if (secDef->GetPDGStable() ) {
2212 weights_v.push_back(weight);
2213 times_v.push_back(currentTime);
2214 secondaries_v.push_back(a_pevap_secondary);
2215 }
2216 } else {
2217 weights_v.push_back(weight);
2218 times_v.push_back(currentTime);
2219 secondaries_v.push_back(a_pevap_secondary);
2220 }
2221 }
2222
2223 delete anITChannel;
2224 }
2225}
G4double
double G4double
Definition: G4Types.hh:83
G4int
int G4int
Definition: G4Types.hh:85
G4DecayProducts::entries
G4int entries() const
Definition: G4DecayProducts.hh:79
G4DecayProducts::PopProducts
G4DynamicParticle * PopProducts()
Definition: G4DecayProducts.cc:136
G4DynamicParticle::GetDefinition
G4ParticleDefinition * GetDefinition() const
G4ITDecay::DecayIt
virtual G4DecayProducts * DecayIt(G4double)
Definition: G4ITDecay.cc:74
G4Ions
Definition: G4Ions.hh:52
G4ParticleDefinition::GetPDGStable
G4bool GetPDGStable() const
G4ParticleDefinition::GetPDGLifeTime
G4double GetPDGLifeTime() const

Referenced by DecayIt().

◆ AddUserDecayDataFile()

void G4RadioactiveDecay::AddUserDecayDataFile ( G4int  Z,
G4int  A,
G4String  filename 
)

Definition at line 1192 of file G4RadioactiveDecay.cc.

1193{
1194 if (Z < 1 || A < 2) G4cout << "Z and A not valid!" << G4endl;
1195
1196 std::ifstream DecaySchemeFile(filename);
1197 if (DecaySchemeFile) {
1198 G4int ID_ion = A*1000 + Z;
1199 theUserRadioactiveDataFiles[ID_ion] = filename;
1200 } else {
1201 G4cout << "The file " << filename << " does not exist!" << G4endl;
1202 }
1203}
A
double A(double temperature)
Definition: G4DNAElectronHoleRecombination.cc:60

Referenced by G4RadioactiveDecaymessenger::SetNewValue().

◆ BuildPhysicsTable()

void G4RadioactiveDecay::BuildPhysicsTable ( const G4ParticleDefinition &  )
virtual

Reimplemented from G4VProcess.

Definition at line 740 of file G4RadioactiveDecay.cc.

741{
742 if (!isInitialised) {
743 isInitialised = true;
744#ifdef G4VERBOSE
745 if(G4HadronicParameters::Instance()->GetVerboseLevel() > 0 &&
746 G4Threading::IsMasterThread()) { StreamInfo(G4cout, "\n"); }
747#endif
748 }
749 G4HadronicProcessStore::
750 Instance()->RegisterParticleForExtraProcess(this,G4GenericIon::GenericIon());
751}
G4GenericIon::GenericIon
static G4GenericIon * GenericIon()
Definition: G4GenericIon.cc:92
G4HadronicParameters::Instance
static G4HadronicParameters * Instance()
Definition: G4HadronicParameters.cc:52
G4HadronicProcessStore::RegisterParticleForExtraProcess
void RegisterParticleForExtraProcess(G4VProcess *, const G4ParticleDefinition *)
Definition: G4HadronicProcessStore.cc:508
G4Threading::IsMasterThread
G4bool IsMasterThread()
Definition: G4Threading.cc:124

◆ CalculateChainsFromParent()

void G4RadioactiveDecay::CalculateChainsFromParent ( const G4ParticleDefinition &  theParentNucleus)

Definition at line 1222 of file G4RadioactiveDecay.cc.

1224{
1225 // Use extended Bateman equation to calculate the radioactivities of all
1226 // progeny of theParentNucleus. The coefficients required to do this are
1227 // calculated using the method of P. Truscott (Ph.D. thesis and
1228 // DERA Technical Note DERA/CIS/CIS2/7/36/4/10) 11 January 2000.
1229 // Coefficients are then added to the decay rate table vector
1230
1231 // Create and initialise variables used in the method.
1232 theDecayRateVector.clear();
1233
1234 G4int nGeneration = 0;
1235
1236 std::vector<G4double> taos;
1237
1238 // Dimensionless A coefficients of Eqs. 4.24 and 4.25 of the TN
1239 std::vector<G4double> Acoeffs;
1240
1241 // According to Eq. 4.26 the first coefficient (A_1:1) is -1
1242 Acoeffs.push_back(-1.);
1243
1244 G4int A = ((const G4Ions*)(&theParentNucleus))->GetAtomicMass();
1245 G4int Z = ((const G4Ions*)(&theParentNucleus))->GetAtomicNumber();
1246 G4double E = ((const G4Ions*)(&theParentNucleus))->GetExcitationEnergy();
1247 G4double tao = theParentNucleus.GetPDGLifeTime();
1248 if (tao < 0.) tao = 1e-100;
1249 taos.push_back(tao);
1250 G4int nEntry = 0;
1251
1252 // Fill the decay rate container (G4RadioactiveDecayRatesToDaughter) with
1253 // the parent isotope data
1254 SetDecayRate(Z,A,E,nGeneration,Acoeffs,taos);
1255
1256 // store the decay rate in decay rate vector
1257 theDecayRateVector.push_back(ratesToDaughter);
1258 nEntry++;
1259
1260 // Now start treating the secondary generations.
1261 G4bool stable = false;
1262 G4int i;
1263 G4int j;
1264 G4VDecayChannel* theChannel = 0;
1265 G4NuclearDecay* theNuclearDecayChannel = 0;
1266
1267 G4ITDecay* theITChannel = 0;
1268 G4BetaMinusDecay* theBetaMinusChannel = 0;
1269 G4BetaPlusDecay* theBetaPlusChannel = 0;
1270 G4AlphaDecay* theAlphaChannel = 0;
1271 G4ProtonDecay* theProtonChannel = 0;
1272 G4TritonDecay* theTritonChannel = 0;
1273 G4NeutronDecay* theNeutronChannel = 0;
1274 G4SFDecay* theFissionChannel = 0;
1275
1276 G4RadioactiveDecayMode theDecayMode;
1277 G4double theBR = 0.0;
1278 G4int AP = 0;
1279 G4int ZP = 0;
1280 G4int AD = 0;
1281 G4int ZD = 0;
1282 G4double EP = 0.;
1283 std::vector<G4double> TP;
1284 std::vector<G4double> RP; // A coefficients of the previous generation
1285 G4ParticleDefinition *theDaughterNucleus;
1286 G4double daughterExcitation;
1287 G4double nearestEnergy = 0.0;
1288 G4int nearestLevelIndex = 0;
1289 G4ParticleDefinition *aParentNucleus;
1290 G4IonTable* theIonTable;
1291 G4DecayTable* parentDecayTable;
1292 G4double theRate;
1293 G4double TaoPlus;
1294 G4int nS = 0; // Running index of first decay in a given generation
1295 G4int nT = nEntry; // Total number of decays accumulated over entire history
1296 const G4int nMode = 12;
1297 G4double brs[nMode];
1298
1299 //
1300 theIonTable =
1301 (G4IonTable*)(G4ParticleTable::GetParticleTable()->GetIonTable());
1302
1303 G4int loop = 0;
1304
1305 while (!stable) { /* Loop checking, 01.09.2015, D.Wright */
1306 loop++;
1307 if (loop > 10000) {
1308 G4Exception("G4RadioactiveDecay::CalculateChainsFromParent()", "HAD_RDM_100",
1309 JustWarning, "While loop count exceeded");
1310 break;
1311 }
1312 nGeneration++;
1313 for (j = nS; j < nT; j++) {
1314 // First time through, get data for parent nuclide
1315 ZP = theDecayRateVector[j].GetZ();
1316 AP = theDecayRateVector[j].GetA();
1317 EP = theDecayRateVector[j].GetE();
1318 RP = theDecayRateVector[j].GetDecayRateC();
1319 TP = theDecayRateVector[j].GetTaos();
1320 if (GetVerboseLevel() > 1) {
1321 G4cout << "G4RadioactiveDecay::CalculateChainsFromParent: daughters of ("
1322 << ZP << ", " << AP << ", " << EP
1323 << ") are being calculated, generation = " << nGeneration
1324 << G4endl;
1325 }
1326// G4cout << " Taus = " << G4endl;
1327// for (G4int ii = 0; ii < TP.size(); ii++) G4cout << TP[ii] << ", " ;
1328// G4cout << G4endl;
1329
1330 aParentNucleus = theIonTable->GetIon(ZP,AP,EP);
1331 parentDecayTable = GetDecayTable(aParentNucleus);
1332
1333 G4DecayTable* summedDecayTable = new G4DecayTable();
1334 // This instance of G4DecayTable is for accumulating BRs and decay
1335 // channels. It will contain one decay channel per type of decay
1336 // (alpha, beta, etc.); its branching ratio will be the sum of all
1337 // branching ratios for that type of decay of the parent. If the
1338 // halflife of a particular channel is longer than some threshold,
1339 // that channel will be inserted specifically and its branching
1340 // ratio will not be included in the above sums.
1341 // This instance is not used to perform actual decays.
1342
1343 for (G4int k = 0; k < nMode; k++) brs[k] = 0.0;
1344
1345 // Go through the decay table and sum all channels having the same decay mode
1346 for (i = 0; i < parentDecayTable->entries(); i++) {
1347 theChannel = parentDecayTable->GetDecayChannel(i);
1348 theNuclearDecayChannel = static_cast<G4NuclearDecay*>(theChannel);
1349 theDecayMode = theNuclearDecayChannel->GetDecayMode();
1350 daughterExcitation = theNuclearDecayChannel->GetDaughterExcitation();
1351 theDaughterNucleus = theNuclearDecayChannel->GetDaughterNucleus() ;
1352
1353 AD = ((const G4Ions*)(theDaughterNucleus))->GetAtomicMass();
1354 ZD = ((const G4Ions*)(theDaughterNucleus))->GetAtomicNumber();
1355 const G4LevelManager* levelManager =
1356 G4NuclearLevelData::GetInstance()->GetLevelManager(ZD,AD);
1357
1358 if (levelManager->NumberOfTransitions() ) {
1359 nearestEnergy = levelManager->NearestLevelEnergy(daughterExcitation);
1360 if (std::abs(daughterExcitation - nearestEnergy) < levelTolerance) {
1361 // Level half-life is in ns and the threshold is set to 1 micros
1362 // by default, user can set it via the UI command
1363 nearestLevelIndex = levelManager->NearestLevelIndex(daughterExcitation);
1364 if (levelManager->LifeTime(nearestLevelIndex)*ns >= halflifethreshold){
1365 // save the metastable nucleus
1366 summedDecayTable->Insert(theChannel);
1367 } else {
1368 brs[theDecayMode] += theChannel->GetBR();
1369 }
1370 } else {
1371 brs[theDecayMode] += theChannel->GetBR();
1372 }
1373 } else {
1374 brs[theDecayMode] += theChannel->GetBR();
1375 }
1376 } // Combine decay channels (loop i)
1377
1378 brs[2] = brs[2]+brs[3]+brs[4]+brs[5]+brs[6]; // Combine beta+ and EC
1379 brs[3] = brs[4] = brs[5] = brs[6] = 0.0;
1380 for (i = 0; i < nMode; i++) { // loop over decay modes
1381 if (brs[i] > 0.) {
1382 switch (i) {
1383 case 0:
1384 // Decay mode is isomeric transition
1385 theITChannel = new G4ITDecay(aParentNucleus, brs[0], 0.0, 0.0,
1386 photonEvaporation);
1387
1388 summedDecayTable->Insert(theITChannel);
1389 break;
1390
1391 case 1:
1392 // Decay mode is beta-
1393 theBetaMinusChannel = new G4BetaMinusDecay(aParentNucleus, brs[1],
1394 0.*MeV, 0.*MeV,
1395 noFloat, allowed);
1396 summedDecayTable->Insert(theBetaMinusChannel);
1397 break;
1398
1399 case 2:
1400 // Decay mode is beta+ + EC.
1401 theBetaPlusChannel = new G4BetaPlusDecay(aParentNucleus, brs[2], // DHW: April 2015
1402 0.*MeV, 0.*MeV,
1403 noFloat, allowed);
1404 summedDecayTable->Insert(theBetaPlusChannel);
1405 break;
1406
1407 case 7:
1408 // Decay mode is alpha.
1409 theAlphaChannel = new G4AlphaDecay(aParentNucleus, brs[7], 0.*MeV,
1410 0.*MeV, noFloat);
1411 summedDecayTable->Insert(theAlphaChannel);
1412 break;
1413
1414 case 8:
1415 // Decay mode is proton.
1416 theProtonChannel = new G4ProtonDecay(aParentNucleus, brs[8], 0.*MeV,
1417 0.*MeV, noFloat);
1418 summedDecayTable->Insert(theProtonChannel);
1419 break;
1420
1421 case 9:
1422 // Decay mode is neutron.
1423 theNeutronChannel = new G4NeutronDecay(aParentNucleus, brs[9], 0.*MeV,
1424 0.*MeV, noFloat);
1425 summedDecayTable->Insert(theNeutronChannel);
1426 break;
1427
1428 case 10:
1429 // Decay mode is spontaneous fission
1430 theFissionChannel = new G4SFDecay(aParentNucleus, brs[10], 0.*MeV,
1431 0.*MeV, noFloat);
1432 summedDecayTable->Insert(theFissionChannel);
1433 break;
1434 case 11:
1435 // Decay mode is triton.
1436 theTritonChannel = new G4TritonDecay(aParentNucleus, brs[11], 0.*MeV,
1437 0.*MeV, noFloat);
1438 summedDecayTable->Insert(theTritonChannel);
1439 break;
1440
1441 default:
1442 break;
1443 }
1444 }
1445 }
1446 // loop over all branches in summedDecayTable
1447 //
1448 for (i = 0; i < summedDecayTable->entries(); i++){
1449 theChannel = summedDecayTable->GetDecayChannel(i);
1450 theNuclearDecayChannel = static_cast<G4NuclearDecay*>(theChannel);
1451 theBR = theChannel->GetBR();
1452 theDaughterNucleus = theNuclearDecayChannel->GetDaughterNucleus();
1453
1454 // First check if the decay of the original nucleus is an IT channel,
1455 // if true create a new ground-state nucleus
1456 if (theNuclearDecayChannel->GetDecayMode() == IT && nGeneration == 1) {
1457
1458 A = ((const G4Ions*)(theDaughterNucleus))->GetAtomicMass();
1459 Z = ((const G4Ions*)(theDaughterNucleus))->GetAtomicNumber();
1460 theDaughterNucleus=theIonTable->GetIon(Z,A,0.);
1461 }
1462 if (IsApplicable(*theDaughterNucleus) && theBR &&
1463 aParentNucleus != theDaughterNucleus) {
1464 // need to make sure daughter has decay table
1465 parentDecayTable = GetDecayTable(theDaughterNucleus);
1466
1467 if (parentDecayTable->entries() ) {
1468 A = ((const G4Ions*)(theDaughterNucleus))->GetAtomicMass();
1469 Z = ((const G4Ions*)(theDaughterNucleus))->GetAtomicNumber();
1470 E = ((const G4Ions*)(theDaughterNucleus))->GetExcitationEnergy();
1471
1472 TaoPlus = theDaughterNucleus->GetPDGLifeTime();
1473 if (TaoPlus <= 0.) TaoPlus = 1e-100;
1474
1475 // first set the taos, one simply need to add to the parent ones
1476 taos.clear();
1477 taos = TP; // load lifetimes of all previous generations
1478 size_t k;
1479 //check that TaoPlus differs from other taos from at least 1.e5 relative difference
1480 //for (k = 0; k < TP.size(); k++){
1481 //if (std::abs((TaoPlus-TP[k])/TP[k])<1.e-5 ) TaoPlus=1.00001*TP[k];
1482 //}
1483 taos.push_back(TaoPlus); // add daughter lifetime to list
1484 // now calculate the coefficiencies
1485 //
1486 // they are in two parts, first the less than n ones
1487 // Eq 4.24 of the TN
1488 Acoeffs.clear();
1489 long double ta1,ta2;
1490 ta2 = (long double)TaoPlus;
1491 for (k = 0; k < RP.size(); k++){
1492 ta1 = (long double)TP[k]; // loop over lifetimes of all previous generations
1493 if (ta1 == ta2) {
1494 theRate = 1.e100;
1495 } else {
1496 theRate = ta1/(ta1-ta2);
1497 }
1498 theRate = theRate * theBR * RP[k];
1499 Acoeffs.push_back(theRate);
1500 }
1501
1502 // the second part: the n:n coefficiency
1503 // Eq 4.25 of the TN. Note Yn+1 is zero apart from Y1 which is -1
1504 // as treated at line 1013
1505 theRate = 0.;
1506 long double aRate, aRate1;
1507 aRate1 = 0.L;
1508 for (k = 0; k < RP.size(); k++){
1509 ta1 = (long double)TP[k];
1510 if (ta1 == ta2 ) {
1511 aRate = 1.e100;
1512 } else {
1513 aRate = ta2/(ta1-ta2);
1514 }
1515 aRate = aRate * (long double)(theBR * RP[k]);
1516 aRate1 += aRate;
1517 }
1518 theRate = -aRate1;
1519 Acoeffs.push_back(theRate);
1520 SetDecayRate (Z,A,E,nGeneration,Acoeffs,taos);
1521 theDecayRateVector.push_back(ratesToDaughter);
1522 nEntry++;
1523 } // there are entries in the table
1524 } // nuclide is OK to decay
1525 } // end of loop (i) over decay table branches
1526
1527 delete summedDecayTable;
1528
1529 } // Getting contents of decay rate vector (end loop on j)
1530 nS = nT;
1531 nT = nEntry;
1532 if (nS == nT) stable = true;
1533 } // while nuclide is not stable
1534
1535 // end of while loop
1536 // the calculation completed here
1537
1538
1539 // fill the first part of the decay rate table
1540 // which is the name of the original particle (isotope)
1541 chainsFromParent.SetIonName(theParentNucleus.GetParticleName());
1542
1543 // now fill the decay table with the newly completed decay rate vector
1544 chainsFromParent.SetItsRates(theDecayRateVector);
1545
1546 // finally add the decayratetable to the tablevector
1547 theParentChainTable.push_back(chainsFromParent);
1548}
allowed
@ allowed
Definition: G4BetaDecayType.hh:40
JustWarning
@ JustWarning
Definition: G4ExceptionSeverity.hh:65
noFloat
#define noFloat
Definition: G4Ions.hh:112
G4RadioactiveDecayMode
G4RadioactiveDecayMode
Definition: G4RadioactiveDecayMode.hh:67
G4bool
bool G4bool
Definition: G4Types.hh:86
G4DecayTable::GetDecayChannel
G4VDecayChannel * GetDecayChannel(G4int index) const
Definition: G4DecayTable.hh:112
G4DecayTable::Insert
void Insert(G4VDecayChannel *aChannel)
Definition: G4DecayTable.cc:53
G4DecayTable::entries
G4int entries() const
Definition: G4DecayTable.hh:100
G4IonTable::GetIon
G4ParticleDefinition * GetIon(G4int Z, G4int A, G4int lvl=0)
Definition: G4IonTable.cc:522
G4LevelManager::LifeTime
G4double LifeTime(size_t i) const
Definition: G4LevelManager.hh:186
G4LevelManager::NearestLevelIndex
size_t NearestLevelIndex(G4double energy, size_t index=0) const
Definition: G4LevelManager.cc:95
G4LevelManager::NearestLevelEnergy
G4double NearestLevelEnergy(G4double energy, size_t index=0) const
Definition: G4LevelManager.hh:176
G4LevelManager::NumberOfTransitions
size_t NumberOfTransitions() const
Definition: G4LevelManager.hh:133
G4NuclearDecay::GetDaughterExcitation
G4double GetDaughterExcitation()
Definition: G4NuclearDecay.hh:55
G4NuclearDecay::GetDaughterNucleus
G4ParticleDefinition * GetDaughterNucleus()
Definition: G4NuclearDecay.hh:59
G4NuclearDecay::GetDecayMode
G4RadioactiveDecayMode GetDecayMode()
Definition: G4NuclearDecay.hh:53
G4NuclearLevelData::GetLevelManager
const G4LevelManager * GetLevelManager(G4int Z, G4int A)
Definition: G4NuclearLevelData.cc:467
G4NuclearLevelData::GetInstance
static G4NuclearLevelData * GetInstance()
Definition: G4NuclearLevelData.cc:421
G4ParticleDefinition::GetParticleName
const G4String & GetParticleName() const
Definition: G4ParticleDefinition.hh:108
G4ParticleTable::GetIonTable
G4IonTable * GetIonTable() const
Definition: G4ParticleTable.cc:695
G4ParticleTable::GetParticleTable
static G4ParticleTable * GetParticleTable()
Definition: G4ParticleTable.cc:87
G4RadioactiveDecayChainsFromParent::SetItsRates
void SetItsRates(G4RadioactiveDecayRates arate)
Definition: G4RadioactiveDecayChainsFromParent.hh:90
G4RadioactiveDecay::GetDecayTable
G4DecayTable * GetDecayTable(const G4ParticleDefinition *)
Definition: G4RadioactiveDecay.cc:312
G4RadioactiveDecay::SetDecayRate
void SetDecayRate(G4int, G4int, G4double, G4int, std::vector< G4double >, std::vector< G4double >)
Definition: G4RadioactiveDecay.cc:1207
G4RadioactiveDecay::IsApplicable
G4bool IsApplicable(const G4ParticleDefinition &)
Definition: G4RadioactiveDecay.cc:286
G4VDecayChannel::GetBR
G4double GetBR() const
Definition: G4VDecayChannel.hh:297

Referenced by DecayIt().

◆ ChooseCollimationDirection()

G4ThreeVector G4RadioactiveDecay::ChooseCollimationDirection ( ) const
protected

Definition at line 2159 of file G4RadioactiveDecay.cc.

2159 {
2160 if (origin == forceDecayDirection) return origin; // Don't do collimation
2161 if (forceDecayHalfAngle == 180.*deg) return origin;
2162
2163 G4ThreeVector dir = forceDecayDirection;
2164
2165 // Return direction offset by random throw
2166 if (forceDecayHalfAngle > 0.) {
2167 // Generate uniform direction around central axis
2168 G4double phi = 2.*pi*G4UniformRand();
2169 G4double cosMin = std::cos(forceDecayHalfAngle);
2170 G4double cosTheta = (1.-cosMin)*G4UniformRand() + cosMin; // [cosMin,1.)
2171
2172 dir.setPhi(dir.phi()+phi);
2173 dir.setTheta(dir.theta()+std::acos(cosTheta));
2174 }
2175
2176#ifdef G4VERBOSE
2177 if (GetVerboseLevel()>1)
2178 G4cout << " ChooseCollimationDirection returns " << dir << G4endl;
2179#endif
2180
2181 return dir;
2182}
G4UniformRand
#define G4UniformRand()
Definition: Randomize.hh:52
CLHEP::Hep3Vector::phi
double phi() const
CLHEP::Hep3Vector::theta
double theta() const
CLHEP::Hep3Vector::setTheta
void setTheta(double)
CLHEP::Hep3Vector::setPhi
void setPhi(double)
G4INCL::Math::pi
const G4double pi
Definition: G4INCLGlobals.hh:68

Referenced by CollimateDecayProduct().

◆ CollimateDecay()

void G4RadioactiveDecay::CollimateDecay ( G4DecayProducts *  products)
protected

Definition at line 2115 of file G4RadioactiveDecay.cc.

2115 {
2116 if (origin == forceDecayDirection) return; // No collimation requested
2117 if (180.*deg == forceDecayHalfAngle) return;
2118 if (0 == products || 0 == products->entries()) return;
2119
2120#ifdef G4VERBOSE
2121 if (GetVerboseLevel() > 1) G4cout << "Begin decay collimation " << G4endl;
2122#endif
2123
2124 // Particles suitable for directional biasing (for if-blocks below)
2125 static const G4ParticleDefinition* electron = G4Electron::Definition();
2126 static const G4ParticleDefinition* positron = G4Positron::Definition();
2127 static const G4ParticleDefinition* neutron = G4Neutron::Definition();
2128 static const G4ParticleDefinition* gamma = G4Gamma::Definition();
2129 static const G4ParticleDefinition* alpha = G4Alpha::Definition();
2130 static const G4ParticleDefinition* triton = G4Triton::Definition();
2131 static const G4ParticleDefinition* proton = G4Proton::Definition();
2132
2133 G4ThreeVector newDirection; // Re-use to avoid memory churn
2134 for (G4int i=0; i<products->entries(); i++) {
2135 G4DynamicParticle* daughter = (*products)[i];
2136 const G4ParticleDefinition* daughterType =
2137 daughter->GetParticleDefinition();
2138 if (daughterType == electron || daughterType == positron ||
2139 daughterType == neutron || daughterType == gamma ||
2140 daughterType == alpha || daughterType == triton || daughterType == proton) CollimateDecayProduct(daughter);
2141 }
2142}
G4Alpha::Definition
static G4Alpha * Definition()
Definition: G4Alpha.cc:48
G4DynamicParticle::GetParticleDefinition
const G4ParticleDefinition * GetParticleDefinition() const
G4Electron::Definition
static G4Electron * Definition()
Definition: G4Electron.cc:48
G4Gamma::Definition
static G4Gamma * Definition()
Definition: G4Gamma.cc:48
G4Neutron::Definition
static G4Neutron * Definition()
Definition: G4Neutron.cc:53
G4Positron::Definition
static G4Positron * Definition()
Definition: G4Positron.cc:48
G4Proton::Definition
static G4Proton * Definition()
Definition: G4Proton.cc:48
G4RadioactiveDecay::CollimateDecayProduct
void CollimateDecayProduct(G4DynamicParticle *product)
Definition: G4RadioactiveDecay.cc:2144
G4Triton::Definition
static G4Triton * Definition()
Definition: G4Triton.cc:49

Referenced by DoDecay().

◆ CollimateDecayProduct()

void G4RadioactiveDecay::CollimateDecayProduct ( G4DynamicParticle *  product)
protected

Definition at line 2144 of file G4RadioactiveDecay.cc.

2144 {
2145#ifdef G4VERBOSE
2146 if (GetVerboseLevel() > 1) {
2147 G4cout << "CollimateDecayProduct for daughter "
2148 << daughter->GetParticleDefinition()->GetParticleName() << G4endl;
2149 }
2150#endif
2151
2152 G4ThreeVector collimate = ChooseCollimationDirection();
2153 if (origin != collimate) daughter->SetMomentumDirection(collimate);
2154}
G4RadioactiveDecay::ChooseCollimationDirection
G4ThreeVector ChooseCollimationDirection() const
Definition: G4RadioactiveDecay.cc:2159

Referenced by CollimateDecay().

◆ ConvolveSourceTimeProfile()

G4double G4RadioactiveDecay::ConvolveSourceTimeProfile ( const G4double  t,
const G4double  tau 
)
protected

Definition at line 520 of file G4RadioactiveDecay.cc.

521{
522 G4double convolvedTime = 0.0;
523 G4int nbin;
524 if ( t > SBin[NSourceBin]) {
525 nbin = NSourceBin;
526 } else {
527 nbin = 0;
528
529 G4int loop = 0;
530 while (t > SBin[nbin]) { // Loop checking, 01.09.2015, D.Wright
531 loop++;
532 if (loop > 1000) {
533 G4Exception("G4RadioactiveDecay::ConvolveSourceTimeProfile()",
534 "HAD_RDM_100", JustWarning,"While loop count exceeded");
535 break;
536 }
537 nbin++;
538 }
539 nbin--;
540 }
541
542 // Use expm1 wherever possible to avoid large cancellation errors in
543 // 1 - exp(x) for small x
544 G4double earg = 0.0;
545 if (nbin > 0) {
546 for (G4int i = 0; i < nbin; i++) {
547 earg = (SBin[i+1] - SBin[i])/tau;
548 if (earg < 100.) {
549 convolvedTime += SProfile[i] * std::exp((SBin[i] - t)/tau) *
550 std::expm1(earg);
551 } else {
552 convolvedTime += SProfile[i] *
553 (std::exp(-(t-SBin[i+1])/tau)-std::exp(-(t-SBin[i])/tau));
554 }
555 }
556 }
557 convolvedTime -= SProfile[nbin] * std::expm1((SBin[nbin] - t)/tau);
558 // tau divided out of final result to provide probability of decay in window
559
560 if (convolvedTime < 0.) {
561 G4cout << " Convolved time =: " << convolvedTime << " reset to zero! " << G4endl;
562 G4cout << " t = " << t << " tau = " << tau << G4endl;
563 G4cout << SBin[nbin] << " " << SBin[0] << G4endl;
564 convolvedTime = 0.;
565 }
566#ifdef G4VERBOSE
567 if (GetVerboseLevel() > 1)
568 G4cout << " Convolved time: " << convolvedTime << G4endl;
569#endif
570 return convolvedTime;
571}

Referenced by DecayIt().

◆ DecayIt()

G4VParticleChange * G4RadioactiveDecay::DecayIt ( const G4Track &  theTrack,
const G4Step &  theStep 
)

Definition at line 1664 of file G4RadioactiveDecay.cc.

1665{
1666 // Initialize G4ParticleChange object, get particle details and decay table
1667
1668 fParticleChangeForRadDecay.Initialize(theTrack);
1669 fParticleChangeForRadDecay.ProposeWeight(theTrack.GetWeight());
1670 const G4DynamicParticle* theParticle = theTrack.GetDynamicParticle();
1671 const G4ParticleDefinition* theParticleDef = theParticle->GetDefinition();
1672
1673 // First check whether RDM applies to the current logical volume
1674 if (!isAllVolumesMode) {
1675 if (!std::binary_search(ValidVolumes.begin(), ValidVolumes.end(),
1676 theTrack.GetVolume()->GetLogicalVolume()->GetName())) {
1677#ifdef G4VERBOSE
1678 if (GetVerboseLevel()>1) {
1679 G4cout <<"G4RadioactiveDecay::DecayIt : "
1680 << theTrack.GetVolume()->GetLogicalVolume()->GetName()
1681 << " is not selected for the RDM"<< G4endl;
1682 G4cout << " There are " << ValidVolumes.size() << " volumes" << G4endl;
1683 G4cout << " The Valid volumes are " << G4endl;
1684 for (size_t i = 0; i< ValidVolumes.size(); i++)
1685 G4cout << ValidVolumes[i] << G4endl;
1686 }
1687#endif
1688 fParticleChangeForRadDecay.SetNumberOfSecondaries(0);
1689
1690 // Kill the parent particle.
1691 fParticleChangeForRadDecay.ProposeTrackStatus(fStopAndKill) ;
1692 fParticleChangeForRadDecay.ProposeLocalEnergyDeposit(0.0);
1693 ClearNumberOfInteractionLengthLeft();
1694 return &fParticleChangeForRadDecay;
1695 }
1696 }
1697
1698 // Now check if particle is valid for RDM
1699 if (!(IsApplicable(*theParticleDef) ) ) {
1700 // Particle is not an ion or is outside the nucleuslimits for decay
1701#ifdef G4VERBOSE
1702 if (GetVerboseLevel() > 1) {
1703 G4cout << " G4RadioactiveDecay::DecayIt : "
1704 << theParticleDef->GetParticleName()
1705 << " is not a valid nucleus for the RDM. "<< G4endl;
1706 G4cout << " Set particle change accordingly. " << G4endl;
1707 }
1708#endif
1709 fParticleChangeForRadDecay.SetNumberOfSecondaries(0);
1710
1711 // Kill the parent particle
1712 fParticleChangeForRadDecay.ProposeTrackStatus(fStopAndKill) ;
1713 fParticleChangeForRadDecay.ProposeLocalEnergyDeposit(0.0);
1714 ClearNumberOfInteractionLengthLeft();
1715 return &fParticleChangeForRadDecay;
1716 }
1717 G4DecayTable* theDecayTable = GetDecayTable(theParticleDef);
1718
1719 if (theDecayTable == 0 || theDecayTable->entries() == 0) {
1720 // No data in the decay table. Set particle change parameters
1721 // to indicate this.
1722#ifdef G4VERBOSE
1723 if (GetVerboseLevel() > 1) {
1724 G4cout <<" G4RadioactiveDecay::DecayIt : decay table not defined for ";
1725 G4cout << theParticleDef->GetParticleName() << G4endl;
1726 G4cout << " Set particle change to indicate this. " << G4endl;
1727 }
1728#endif
1729 fParticleChangeForRadDecay.SetNumberOfSecondaries(0);
1730
1731 // Kill the parent particle.
1732 fParticleChangeForRadDecay.ProposeTrackStatus(fStopAndKill) ;
1733 fParticleChangeForRadDecay.ProposeLocalEnergyDeposit(0.0);
1734 ClearNumberOfInteractionLengthLeft();
1735 return &fParticleChangeForRadDecay;
1736
1737 } else {
1738 // Data found. Try to decay nucleus
1739 G4double energyDeposit = 0.0;
1740 G4double finalGlobalTime = theTrack.GetGlobalTime();
1741 G4double finalLocalTime = theTrack.GetLocalTime();
1742 G4int index;
1743 G4ThreeVector currentPosition;
1744 currentPosition = theTrack.GetPosition();
1745
1746 // Check whether use Analogue or VR implementation
1747 if (AnalogueMC) {
1748#ifdef G4VERBOSE
1749 if (GetVerboseLevel() > 1)
1750 G4cout <<"DecayIt: Analogue MC version " << G4endl;
1751# endif
1752
1753 G4DecayProducts* products = DoDecay(*theParticleDef);
1754
1755 // Check if the product is the same as input and kill the track if
1756 // necessary to prevent infinite loop (11/05/10, F.Lei)
1757 if ( products->entries() == 1) {
1758 fParticleChangeForRadDecay.SetNumberOfSecondaries(0);
1759 fParticleChangeForRadDecay.ProposeTrackStatus(fStopAndKill);
1760 fParticleChangeForRadDecay.ProposeLocalEnergyDeposit(0.0);
1761 ClearNumberOfInteractionLengthLeft();
1762 return &fParticleChangeForRadDecay;
1763 }
1764
1765 // Get parent particle information and boost the decay products to the
1766 // laboratory frame based on this information.
1767
1768 //The Parent Energy used for the boost should be the total energy of
1769 // the nucleus of the parent ion without the energy of the shell electrons
1770 // (correction for bug 1359 by L. Desorgher)
1771 G4double ParentEnergy = theParticle->GetKineticEnergy()
1772 + theParticle->GetParticleDefinition()->GetPDGMass();
1773 G4ThreeVector ParentDirection(theParticle->GetMomentumDirection());
1774
1775 if (theTrack.GetTrackStatus() == fStopButAlive) {
1776 //this condition seems to be always True, further investigation is needed (L.Desorgher)
1777
1778 // The particle is decayed at rest.
1779 // since the time is still for rest particle in G4 we need to add the
1780 // additional time lapsed between the particle come to rest and the
1781 // actual decay. This time is simply sampled with the mean-life of
1782 // the particle. But we need to protect the case PDGTime < 0.
1783 // (F.Lei 11/05/10)
1784 G4double temptime = -std::log( G4UniformRand())
1785 *theParticleDef->GetPDGLifeTime();
1786 if (temptime < 0.) temptime = 0.;
1787 finalGlobalTime += temptime;
1788 finalLocalTime += temptime;
1789 energyDeposit += theParticle->GetKineticEnergy();
1790 }
1791 products->Boost(ParentEnergy, ParentDirection);
1792
1793 // Add products in theParticleChangeForRadDecay.
1794 G4int numberOfSecondaries = products->entries();
1795 fParticleChangeForRadDecay.SetNumberOfSecondaries(numberOfSecondaries);
1796#ifdef G4VERBOSE
1797 if (GetVerboseLevel()>1) {
1798 G4cout <<"G4RadioactiveDecay::DecayIt : Decay vertex :";
1799 G4cout <<" Time: " <<finalGlobalTime/ns <<"[ns]";
1800 G4cout <<" X:" <<(theTrack.GetPosition()).x() /cm <<"[cm]";
1801 G4cout <<" Y:" <<(theTrack.GetPosition()).y() /cm <<"[cm]";
1802 G4cout <<" Z:" <<(theTrack.GetPosition()).z() /cm <<"[cm]";
1803 G4cout << G4endl;
1804 G4cout <<"G4Decay::DecayIt : decay products in Lab. Frame" <<G4endl;
1805 products->DumpInfo();
1806 products->IsChecked();
1807 }
1808#endif
1809 for (index=0; index < numberOfSecondaries; index++) {
1810 G4Track* secondary = new G4Track(products->PopProducts(),
1811 finalGlobalTime, currentPosition);
1812
1813 secondary->SetCreatorModelIndex(theRadDecayMode);
1814 //Change for atomics relaxation
1815 if (theRadDecayMode == IT && index>0){
1816 if (index == numberOfSecondaries-1) secondary->SetCreatorModelIndex(IT);
1817 else secondary->SetCreatorModelIndex(30);
1818 }
1819 else if (theRadDecayMode >= KshellEC && theRadDecayMode <= NshellEC
1820 && index <numberOfSecondaries-1){
1821 secondary->SetCreatorModelIndex(30);
1822 }
1823 secondary->SetGoodForTrackingFlag();
1824 secondary->SetTouchableHandle(theTrack.GetTouchableHandle());
1825 fParticleChangeForRadDecay.AddSecondary(secondary);
1826 }
1827 delete products;
1828 // end of analogue MC algorithm
1829
1830 } else {
1831 // Variance Reduction Method
1832#ifdef G4VERBOSE
1833 if (GetVerboseLevel()>0)
1834 G4cout << "DecayIt: Variance Reduction version " << G4endl;
1835#endif
1836 // Get decay chains for the given nuclide
1837 if (!IsRateTableReady(*theParticleDef)) CalculateChainsFromParent(*theParticleDef);
1838 GetChainsFromParent(*theParticleDef);
1839
1840 // declare some of the variables required in the implementation
1841 G4ParticleDefinition* parentNucleus;
1842 G4IonTable* theIonTable;
1843 G4int PZ;
1844 G4int PA;
1845 G4double PE;
1846 G4String keyName;
1847 std::vector<G4double> PT;
1848 std::vector<G4double> PR;
1849 G4double tauprob;
1850 long double decayRate;
1851
1852 size_t i;
1853// size_t j;
1854 G4int numberOfSecondaries;
1855 G4int totalNumberOfSecondaries = 0;
1856 G4double currentTime = 0.;
1857 G4int ndecaych;
1858 G4DynamicParticle* asecondaryparticle;
1859 std::vector<G4DynamicParticle*> secondaryparticles;
1860 std::vector<G4double> pw;
1861 std::vector<G4double> ptime;
1862 pw.clear();
1863 ptime.clear();
1864
1865 //now apply the nucleus splitting
1866 for (G4int n = 0; n < NSplit; n++) {
1867 // Get the decay time following the decay probability function
1868 // suppllied by user
1869 G4double theDecayTime = GetDecayTime();
1870 G4int nbin = GetDecayTimeBin(theDecayTime);
1871
1872 // calculate the first part of the weight function
1873 G4double weight1 = 1.;
1874 if (nbin == 1) {
1875 weight1 = 1./DProfile[nbin-1]
1876 *(DBin[nbin]-DBin[nbin-1])/NSplit;
1877 } else if (nbin > 1) {
1878 weight1 = 1./(DProfile[nbin]-DProfile[nbin-2])
1879 *(DBin[nbin]-DBin[nbin-1])/NSplit;
1880 }
1881
1882 // it should be calculated in seconds
1883 weight1 /= s ;
1884
1885 // Loop over all the possible secondaries of the nucleus
1886 // the first one is itself.
1887 for (i = 0; i < theDecayRateVector.size(); i++) {
1888 PZ = theDecayRateVector[i].GetZ();
1889 PA = theDecayRateVector[i].GetA();
1890 PE = theDecayRateVector[i].GetE();
1891 PT = theDecayRateVector[i].GetTaos();
1892 PR = theDecayRateVector[i].GetDecayRateC();
1893
1894 // The array of arrays theDecayRateVector contains all possible decay
1895 // chains of a given parent nucleus (ZP,AP,EP) to a given descendant
1896 // nuclide (Z,A,E).
1897 //
1898 // theDecayRateVector[0] contains the decay parameters of the parent
1899 // nucleus
1900 // PZ = ZP
1901 // PA = AP
1902 // PE = EP
1903 // PT[] = {TP}
1904 // PR[] = {RP}
1905 //
1906 // theDecayRateVector[1] contains the decay of the parent to the first
1907 // generation daughter (Z1,A1,E1).
1908 // PZ = Z1
1909 // PA = A1
1910 // PE = E1
1911 // PT[] = {TP, T1}
1912 // PR[] = {RP, R1}
1913 //
1914 // theDecayRateVector[2] contains the decay of the parent to the first
1915 // generation daughter (Z1,A1,E1) and the decay of the first
1916 // generation daughter to the second generation daughter (Z2,A2,E2).
1917 // PZ = Z2
1918 // PA = A2
1919 // PE = E2
1920 // PT[] = {TP, T1, T2}
1921 // PR[] = {RP, R1, R2}
1922 //
1923 // theDecayRateVector[3] may contain a branch chain
1924 // PZ = Z2a
1925 // PA = A2a
1926 // PE = E2a
1927 // PT[] = {TP, T1, T2a}
1928 // PR[] = {RP, R1, R2a}
1929 //
1930 // and so on.
1931
1932 // Calculate the decay rate of the isotope. decayRate is the
1933 // radioactivity of isotope (PZ,PA,PE) at 'theDecayTime'.
1934 // It will be used to calculate the statistical weight of the
1935 // decay products of this isotope
1936
1937 // For each nuclide, calculate all the decay chains which can reach
1938 // the parent nuclide
1939 decayRate = 0.L;
1940 for (G4int j = 0; j < G4int(PT.size()); j++) {
1941 tauprob = ConvolveSourceTimeProfile(theDecayTime,PT[j]);
1942 // tauprob is dimensionless, PR has units of s-1
1943 decayRate -= PR[j] * (long double)tauprob;
1944 // Eq.4.23 of of the TN
1945 // note the negative here is required as the rate in the
1946 // equation is defined to be negative,
1947 // i.e. decay away, but we need positive value here.
1948
1949 // G4cout << j << "\t" << PT[j]/s << "\t" << PR[j] << "\t" << decayRate << G4endl;
1950 }
1951
1952 // At this point any negative decay rates are probably small enough
1953 // (order 10**-30) that negative values are likely due to cancellation
1954 // errors. Set them to zero.
1955 if (decayRate < 0.0) decayRate = 0.0;
1956/*
1957 if (decayRate < 0.0) {
1958 if (-decayRate > 1.0e-30) {
1959 G4ExceptionDescription ed;
1960 ed << " Negative decay probability (magnitude > 1e-30) \n"
1961 << " in variance reduction branch " << G4endl;
1962 G4Exception("G4RadioactiveDecay::DecayIt()",
1963 "HAD_RDM_200", JustWarning, ed);
1964 } else {
1965 // Decay probability is small enough that negative value is likely
1966 // due to cancellation errors. Set it to zero.
1967 decayRate = 0.0;
1968 }
1969 }
1970
1971 if (decayRate < 0.0) G4cout << " NEGATIVE decay rate = " << decayRate << G4endl;
1972*/
1973 // G4cout << theDecayTime/s << "\t" << nbin << G4endl;
1974 // G4cout << theTrack.GetWeight() << "\t" << weight1 << "\t" << decayRate << G4endl;
1975
1976 // Add isotope to the radioactivity tables
1977 // One table for each observation time window specified in
1978 // SetDecayBias(G4String filename)
1979 theRadioactivityTables[decayWindows[nbin-1]]
1980 ->AddIsotope(PZ,PA,PE,weight1*decayRate,theTrack.GetWeight());
1981
1982 // Now calculate the statistical weight
1983 // One needs to fold the source bias function with the decaytime
1984 // also need to include the track weight! (F.Lei, 28/10/10)
1985 G4double weight = weight1*decayRate*theTrack.GetWeight();
1986
1987 // Decay the isotope
1988 theIonTable = (G4IonTable *)(G4ParticleTable::GetParticleTable()->GetIonTable());
1989 parentNucleus = theIonTable->GetIon(PZ,PA,PE);
1990
1991 // Create a temprary products buffer.
1992 // Its contents to be transfered to the products at the end of the loop
1993 G4DecayProducts* tempprods = 0;
1994
1995 // Decide whether to apply branching ratio bias or not
1996 if (BRBias) {
1997 G4DecayTable* decayTable = GetDecayTable(parentNucleus);
1998
1999 ndecaych = G4int(decayTable->entries()*G4UniformRand());
2000 G4VDecayChannel* theDecayChannel = decayTable->GetDecayChannel(ndecaych);
2001 if (theDecayChannel == 0) {
2002 // Decay channel not found.
2003#ifdef G4VERBOSE
2004 if (GetVerboseLevel() > 0) {
2005 G4cout << " G4RadioactiveDecay::DoIt : cannot determine decay ";
2006 G4cout << " channel for this nucleus; decay as if no biasing active ";
2007 G4cout << G4endl;
2008 decayTable ->DumpInfo();
2009 }
2010#endif
2011 tempprods = DoDecay(*parentNucleus); // DHW 6 Dec 2010 - do decay as if no biasing
2012 // to avoid deref of temppprods = 0
2013 } else {
2014 // A decay channel has been identified, so execute the DecayIt.
2015 G4double tempmass = parentNucleus->GetPDGMass();
2016 tempprods = theDecayChannel->DecayIt(tempmass);
2017 weight *= (theDecayChannel->GetBR())*(decayTable->entries());
2018 }
2019 } else {
2020 tempprods = DoDecay(*parentNucleus);
2021 }
2022
2023 // Save the secondaries for buffers
2024 numberOfSecondaries = tempprods->entries();
2025 currentTime = finalGlobalTime + theDecayTime;
2026 for (index = 0; index < numberOfSecondaries; index++) {
2027 asecondaryparticle = tempprods->PopProducts();
2028 if (asecondaryparticle->GetDefinition()->GetPDGStable() ) {
2029 pw.push_back(weight);
2030 ptime.push_back(currentTime);
2031 secondaryparticles.push_back(asecondaryparticle);
2032 } else if (((const G4Ions*)(asecondaryparticle->GetDefinition()))->GetExcitationEnergy() > 0.
2033 && weight > 0.) {
2034 // Compute the gamma
2035 // Generate gammas and XRays from excited nucleus, added by L.Desorgher
2036 G4ParticleDefinition* apartDef =asecondaryparticle->GetDefinition();
2037 AddDeexcitationSpectrumForBiasMode(apartDef,weight,currentTime,pw,ptime,secondaryparticles);
2038 }
2039 }
2040 delete tempprods;
2041
2042 } // end of i loop
2043 } // end of n loop
2044
2045 // now deal with the secondaries in the two stl containers
2046 // and submmit them back to the tracking manager
2047 totalNumberOfSecondaries = pw.size();
2048 fParticleChangeForRadDecay.SetNumberOfSecondaries(totalNumberOfSecondaries);
2049 for (index=0; index < totalNumberOfSecondaries; index++) {
2050 G4Track* secondary = new G4Track(secondaryparticles[index],
2051 ptime[index], currentPosition);
2052 secondary->SetGoodForTrackingFlag();
2053 secondary->SetTouchableHandle(theTrack.GetTouchableHandle());
2054 secondary->SetWeight(pw[index]);
2055 fParticleChangeForRadDecay.AddSecondary(secondary);
2056 }
2057 // make sure the original track is set to stop and its kinematic energy collected
2058 //
2059 //theTrack.SetTrackStatus(fStopButAlive);
2060 //energyDeposit += theParticle->GetKineticEnergy();
2061
2062 } // End of Variance Reduction
2063
2064 // Kill the parent particle
2065 fParticleChangeForRadDecay.ProposeTrackStatus(fStopAndKill) ;
2066 fParticleChangeForRadDecay.ProposeLocalEnergyDeposit(energyDeposit);
2067 fParticleChangeForRadDecay.ProposeLocalTime(finalLocalTime);
2068 // Reset NumberOfInteractionLengthLeft.
2069 ClearNumberOfInteractionLengthLeft();
2070
2071 return &fParticleChangeForRadDecay;
2072 }
2073}
fStopAndKill
@ fStopAndKill
Definition: G4TrackStatus.hh:46
fStopButAlive
@ fStopButAlive
Definition: G4TrackStatus.hh:44
G4DecayProducts::DumpInfo
void DumpInfo() const
Definition: G4DecayProducts.cc:334
G4DecayProducts::IsChecked
G4bool IsChecked() const
Definition: G4DecayProducts.cc:254
G4DecayProducts::Boost
void Boost(G4double totalEnergy, const G4ThreeVector &momentumDirection)
Definition: G4DecayProducts.cc:177
G4DecayTable::DumpInfo
void DumpInfo() const
Definition: G4DecayTable.cc:120
G4DynamicParticle::GetMomentumDirection
const G4ThreeVector & GetMomentumDirection() const
G4DynamicParticle::GetKineticEnergy
G4double GetKineticEnergy() const
G4LogicalVolume::GetName
const G4String & GetName() const
G4ParticleChangeForDecay::Initialize
virtual void Initialize(const G4Track &)
Definition: G4ParticleChangeForDecay.cc:112
G4ParticleChangeForRadDecay::AddSecondary
void AddSecondary(G4Track *aSecondary)
Definition: G4ParticleChangeForRadDecay.hh:79
G4RadioactiveDecay::DoDecay
G4DecayProducts * DoDecay(const G4ParticleDefinition &theParticleDef)
Definition: G4RadioactiveDecay.cc:2077
G4RadioactiveDecay::GetDecayTimeBin
G4int GetDecayTimeBin(const G4double aDecayTime)
Definition: G4RadioactiveDecay.cc:609
G4RadioactiveDecay::CalculateChainsFromParent
void CalculateChainsFromParent(const G4ParticleDefinition &)
Definition: G4RadioactiveDecay.cc:1223
G4RadioactiveDecay::IsRateTableReady
G4bool IsRateTableReady(const G4ParticleDefinition &)
Definition: G4RadioactiveDecay.cc:425
G4RadioactiveDecay::ConvolveSourceTimeProfile
G4double ConvolveSourceTimeProfile(const G4double, const G4double)
Definition: G4RadioactiveDecay.cc:520
G4RadioactiveDecay::AddDeexcitationSpectrumForBiasMode
void AddDeexcitationSpectrumForBiasMode(G4ParticleDefinition *apartDef, G4double weight, G4double currenTime, std::vector< double > &weights_v, std::vector< double > &times_v, std::vector< G4DynamicParticle * > &secondaries_v)
Definition: G4RadioactiveDecay.cc:2186
G4RadioactiveDecay::GetChainsFromParent
void GetChainsFromParent(const G4ParticleDefinition &)
Definition: G4RadioactiveDecay.cc:438
G4Track::GetTrackStatus
G4TrackStatus GetTrackStatus() const
G4Track::GetVolume
G4VPhysicalVolume * GetVolume() const
G4Track::GetWeight
G4double GetWeight() const
G4Track::SetWeight
void SetWeight(G4double aValue)
G4Track::GetPosition
const G4ThreeVector & GetPosition() const
G4Track::SetTouchableHandle
void SetTouchableHandle(const G4TouchableHandle &apValue)
G4Track::GetGlobalTime
G4double GetGlobalTime() const
G4Track::GetLocalTime
G4double GetLocalTime() const
G4Track::GetDynamicParticle
const G4DynamicParticle * GetDynamicParticle() const
G4Track::GetTouchableHandle
const G4TouchableHandle & GetTouchableHandle() const
G4Track::SetCreatorModelIndex
void SetCreatorModelIndex(G4int idx)
G4Track::SetGoodForTrackingFlag
void SetGoodForTrackingFlag(G4bool value=true)
G4VDecayChannel::DecayIt
virtual G4DecayProducts * DecayIt(G4double parentMass=-1.0)=0
G4VParticleChange::ProposeTrackStatus
void ProposeTrackStatus(G4TrackStatus status)
G4VParticleChange::ProposeWeight
void ProposeWeight(G4double finalWeight)
G4VParticleChange::ProposeLocalEnergyDeposit
void ProposeLocalEnergyDeposit(G4double anEnergyPart)
G4VParticleChange::SetNumberOfSecondaries
void SetNumberOfSecondaries(G4int totSecondaries)
G4VPhysicalVolume::GetLogicalVolume
G4LogicalVolume * GetLogicalVolume() const
G4VProcess::ClearNumberOfInteractionLengthLeft
void ClearNumberOfInteractionLengthLeft()
Definition: G4VProcess.hh:424
ns
#define ns
Definition: xmlparse.cc:614

◆ DeselectAllVolumes()

void G4RadioactiveDecay::DeselectAllVolumes ( )

Definition at line 414 of file G4RadioactiveDecay.cc.

415{
416 ValidVolumes.clear();
417 isAllVolumesMode=false;
418#ifdef G4VERBOSE
419 if (GetVerboseLevel() > 1) G4cout << "RDM removed from all volumes" << G4endl;
420#endif
421}

Referenced by G4RadioactiveDecaymessenger::SetNewValue().

◆ DeselectAVolume()

void G4RadioactiveDecay::DeselectAVolume ( const G4String  aVolume)

Definition at line 355 of file G4RadioactiveDecay.cc.

356{
357 G4LogicalVolumeStore* theLogicalVolumes;
358 G4LogicalVolume* volume;
359 theLogicalVolumes = G4LogicalVolumeStore::GetInstance();
360 for (size_t i = 0; i < theLogicalVolumes->size(); i++){
361 volume=(*theLogicalVolumes)[i];
362 if (volume->GetName() == aVolume) {
363 std::vector<G4String>::iterator location;
364 location = std::find(ValidVolumes.begin(),ValidVolumes.end(),aVolume);
365 if (location != ValidVolumes.end()) {
366 ValidVolumes.erase(location);
367 std::sort(ValidVolumes.begin(), ValidVolumes.end());
368 isAllVolumesMode =false;
369 } else {
370 G4ExceptionDescription ed;
371 ed << aVolume << " is not in the list " << G4endl;
372 G4Exception("G4RadioactiveDecayBase::DeselectAVolume()", "HAD_RDM_300",
373 JustWarning, ed);
374 }
375#ifdef G4VERBOSE
376 if (GetVerboseLevel() > 0)
377 G4cout << " DeselectVolume: " << aVolume << " is removed from list "
378 << G4endl;
379#endif
380 } else if (i == theLogicalVolumes->size()) {
381 G4ExceptionDescription ed;
382 ed << aVolume << " is not a valid logical volume name " << G4endl;
383 G4Exception("G4RadioactiveDecayBase::SelectAVolume()", "HAD_RDM_300",
384 JustWarning, ed);
385 }
386 }
387}
G4ExceptionDescription
std::ostringstream G4ExceptionDescription
Definition: G4Exception.hh:40
G4LogicalVolumeStore::GetInstance
static G4LogicalVolumeStore * GetInstance()
Definition: G4LogicalVolumeStore.cc:176

Referenced by G4RadioactiveDecaymessenger::SetNewValue().

◆ DoDecay()

G4DecayProducts * G4RadioactiveDecay::DoDecay ( const G4ParticleDefinition &  theParticleDef)
protected

Definition at line 2077 of file G4RadioactiveDecay.cc.

2078{
2079 G4DecayProducts* products = 0;
2080 G4DecayTable* theDecayTable = GetDecayTable(&theParticleDef);
2081
2082 // Choose a decay channel.
2083 // G4DecayTable::SelectADecayChannel checks to see if sum of daughter masses
2084 // exceeds parent mass. Pass it the parent mass + maximum Q value to account
2085 // for difference in mass defect.
2086 G4double parentPlusQ = theParticleDef.GetPDGMass() + 30.*MeV;
2087 G4VDecayChannel* theDecayChannel = theDecayTable->SelectADecayChannel(parentPlusQ);
2088 if (theDecayChannel == 0) {
2089 // Decay channel not found.
2090 G4ExceptionDescription ed;
2091 ed << " Cannot determine decay channel for " << theParticleDef.GetParticleName() << G4endl;
2092 G4Exception("G4RadioactiveDecay::DoDecay", "HAD_RDM_013",
2093 FatalException, ed);
2094 } else {
2095 // A decay channel has been identified, so execute the DecayIt.
2096#ifdef G4VERBOSE
2097 if (GetVerboseLevel() > 1) {
2098 G4cout << "G4RadioactiveDecay::DoIt : selected decay channel address:";
2099 G4cout << theDecayChannel << G4endl;
2100 }
2101#endif
2102 theRadDecayMode = (static_cast<G4NuclearDecay*>(theDecayChannel))->GetDecayMode();
2103 products = theDecayChannel->DecayIt(theParticleDef.GetPDGMass() );
2104
2105 // Apply directional bias if requested by user
2106 CollimateDecay(products);
2107 }
2108
2109 return products;
2110}
G4DecayTable::SelectADecayChannel
G4VDecayChannel * SelectADecayChannel(G4double parentMass=-1.)
Definition: G4DecayTable.cc:82
G4RadioactiveDecay::CollimateDecay
void CollimateDecay(G4DecayProducts *products)
Definition: G4RadioactiveDecay.cc:2115

Referenced by DecayIt().

◆ GetChainsFromParent()

void G4RadioactiveDecay::GetChainsFromParent ( const G4ParticleDefinition &  aParticle)

Definition at line 438 of file G4RadioactiveDecay.cc.

439{
440 G4String aParticleName = aParticle.GetParticleName();
441
442 for (size_t i = 0; i < theParentChainTable.size(); i++) {
443 if (theParentChainTable[i].GetIonName() == aParticleName) {
444 theDecayRateVector = theParentChainTable[i].GetItsRates();
445 }
446 }
447#ifdef G4VERBOSE
448 if (GetVerboseLevel() > 1) {
449 G4cout << "The DecayRate Table for " << aParticleName << " is selected."
450 << G4endl;
451 }
452#endif
453}

Referenced by DecayIt().

◆ GetDecayDirection()

const G4ThreeVector & G4RadioactiveDecay::GetDecayDirection ( ) const
inline

Definition at line 221 of file G4RadioactiveDecay.hh.

221 {
222 return forceDecayDirection;
223 }

◆ GetDecayHalfAngle()

G4double G4RadioactiveDecay::GetDecayHalfAngle ( ) const
inline

Definition at line 229 of file G4RadioactiveDecay.hh.

229{return forceDecayHalfAngle;}

◆ GetDecayTable()

G4DecayTable * G4RadioactiveDecay::GetDecayTable ( const G4ParticleDefinition *  aNucleus)

Definition at line 312 of file G4RadioactiveDecay.cc.

313{
314 G4String key = aNucleus->GetParticleName();
315 DecayTableMap::iterator table_ptr = dkmap->find(key);
316
317 G4DecayTable* theDecayTable = 0;
318 if (table_ptr == dkmap->end() ) { // If table not there,
319 theDecayTable = LoadDecayTable(*aNucleus); // load from file and
320 if(theDecayTable) (*dkmap)[key] = theDecayTable; // store in library
321 } else {
322 theDecayTable = table_ptr->second;
323 }
324
325 return theDecayTable;
326}
G4RadioactiveDecay::LoadDecayTable
G4DecayTable * LoadDecayTable(const G4ParticleDefinition &theParentNucleus)
Definition: G4RadioactiveDecay.cc:805

Referenced by CalculateChainsFromParent(), DecayIt(), and DoDecay().

◆ GetDecayTime()

G4double G4RadioactiveDecay::GetDecayTime ( )
protected

Definition at line 582 of file G4RadioactiveDecay.cc.

583{
584 G4double decaytime = 0.;
585 G4double rand = G4UniformRand();
586 G4int i = 0;
587
588 G4int loop = 0;
589 while ( DProfile[i] < rand) { /* Loop checking, 01.09.2015, D.Wright */
590 i++;
591 loop++;
592 if (loop > 100000) {
593 G4Exception("G4RadioactiveDecay::GetDecayTime()", "HAD_RDM_100",
594 JustWarning, "While loop count exceeded");
595 break;
596 }
597 }
598
599 rand = G4UniformRand();
600 decaytime = DBin[i] + rand*(DBin[i+1]-DBin[i]);
601#ifdef G4VERBOSE
602 if (GetVerboseLevel() > 1)
603 G4cout <<" Decay time: " <<decaytime/s <<"[s]" <<G4endl;
604#endif
605 return decaytime;
606}

Referenced by DecayIt().

◆ GetDecayTimeBin()

G4int G4RadioactiveDecay::GetDecayTimeBin ( const G4double  aDecayTime)
protected

Definition at line 609 of file G4RadioactiveDecay.cc.

610{
611 G4int i = 0;
612
613 G4int loop = 0;
614 while ( aDecayTime > DBin[i] ) { /* Loop checking, 01.09.2015, D.Wright */
615 i++;
616 loop++;
617 if (loop > 100000) {
618 G4Exception("G4RadioactiveDecay::GetDecayTimeBin()", "HAD_RDM_100",
619 JustWarning, "While loop count exceeded");
620 break;
621 }
622 }
623
624 return i;
625}

Referenced by DecayIt().

◆ GetMeanFreePath()

G4double G4RadioactiveDecay::GetMeanFreePath ( const G4Track &  theTrack,
G4double  previousStepSize,
G4ForceCondition *  condition 
)
protectedvirtual

Implements G4VRestDiscreteProcess.

Definition at line 675 of file G4RadioactiveDecay.cc.

677{
678 const G4DynamicParticle* aParticle = aTrack.GetDynamicParticle();
679 const G4ParticleDefinition* aParticleDef = aParticle->GetDefinition();
680 G4double tau = aParticleDef->GetPDGLifeTime();
681 G4double aMass = aParticle->GetMass();
682
683#ifdef G4VERBOSE
684 if (GetVerboseLevel() > 2) {
685 G4cout << "G4RadioactiveDecay::GetMeanFreePath() " << G4endl;
686 G4cout << " KineticEnergy: " << aParticle->GetKineticEnergy()/GeV
687 << " GeV, Mass: " << aMass/GeV << " GeV, tau: " << tau << " ns "
688 << G4endl;
689 }
690#endif
691 G4double pathlength = DBL_MAX;
692 if (tau != -1) {
693 // Ion can decay
694
695 if (tau < -1000.0) {
696 pathlength = DBL_MIN; // nuclide had very short lifetime or wasn't in table
697
698 } else if (tau < 0.0) {
699 G4cout << aParticleDef->GetParticleName() << " has lifetime " << tau << G4endl;
700 G4ExceptionDescription ed;
701 ed << "Ion has negative lifetime " << tau
702 << " but is not stable. Setting mean free path to DBL_MAX" << G4endl;
703 G4Exception("G4RadioactiveDecay::GetMeanFreePath()", "HAD_RDM_011",
704 JustWarning, ed);
705 pathlength = DBL_MAX;
706
707 } else {
708 // Calculate mean free path
709 G4double betaGamma = aParticle->GetTotalMomentum()/aMass;
710 pathlength = c_light*tau*betaGamma;
711
712 if (pathlength < DBL_MIN) {
713 pathlength = DBL_MIN;
714#ifdef G4VERBOSE
715 if (GetVerboseLevel() > 2) {
716 G4cout << "G4Decay::GetMeanFreePath: "
717 << aParticleDef->GetParticleName()
718 << " stops, kinetic energy = "
719 << aParticle->GetKineticEnergy()/keV <<" keV " << G4endl;
720 }
721#endif
722 }
723 }
724 }
725
726#ifdef G4VERBOSE
727 if (GetVerboseLevel() > 1) {
728 G4cout << "mean free path: "<< pathlength/m << " m" << G4endl;
729 }
730#endif
731 return pathlength;
732}
G4DynamicParticle::GetMass
G4double GetMass() const
G4DynamicParticle::GetTotalMomentum
G4double GetTotalMomentum() const
DBL_MIN
#define DBL_MIN
Definition: templates.hh:54
DBL_MAX
#define DBL_MAX
Definition: templates.hh:62

◆ GetMeanLifeTime()

G4double G4RadioactiveDecay::GetMeanLifeTime ( const G4Track &  theTrack,
G4ForceCondition *  condition 
)
protectedvirtual

Implements G4VRestDiscreteProcess.

Definition at line 633 of file G4RadioactiveDecay.cc.

635{
636 // For variance reduction the time is set to 0 so as to force the particle
637 // to decay immediately.
638 // In analogueMC mode it returns the particle's mean-life.
639
640 G4double meanlife = 0.;
641 if (AnalogueMC) {
642 const G4DynamicParticle* theParticle = theTrack.GetDynamicParticle();
643 const G4ParticleDefinition* theParticleDef = theParticle->GetDefinition();
644 G4double theLife = theParticleDef->GetPDGLifeTime();
645#ifdef G4VERBOSE
646 if (GetVerboseLevel() > 2) {
647 G4cout << "G4RadioactiveDecay::GetMeanLifeTime() " << G4endl;
648 G4cout << "KineticEnergy: " << theParticle->GetKineticEnergy()/GeV
649 << " GeV, Mass: " << theParticle->GetMass()/GeV
650 << " GeV, Life time: " << theLife/ns << " ns " << G4endl;
651 }
652#endif
653 if (theParticleDef->GetPDGStable()) {meanlife = DBL_MAX;}
654 else if (theLife < 0.0) {meanlife = DBL_MAX;}
655 else {meanlife = theLife;}
656 // Set meanlife to zero for excited istopes which are not in the
657 // RDM database
658 if (((const G4Ions*)(theParticleDef))->GetExcitationEnergy() > 0. &&
659 meanlife == DBL_MAX) {meanlife = 0.;}
660 }
661#ifdef G4VERBOSE
662 if (GetVerboseLevel() > 2)
663 G4cout << " mean life time: " << meanlife/s << " s " << G4endl;
664#endif
665
666 return meanlife;
667}

◆ GetNucleusLimits()

G4NucleusLimits G4RadioactiveDecay::GetNucleusLimits ( ) const
inline

Definition at line 181 of file G4RadioactiveDecay.hh.

182 {return theNucleusLimits;}

◆ GetSplitNuclei()

G4int G4RadioactiveDecay::GetSplitNuclei ( )
inline

Definition at line 215 of file G4RadioactiveDecay.hh.

215{return NSplit;}

◆ GetTheRadioactivityTables()

std::vector< G4RadioactivityTable * > G4RadioactiveDecay::GetTheRadioactivityTables ( )
inline

Definition at line 159 of file G4RadioactiveDecay.hh.

160 {return theRadioactivityTables;}

◆ GetVerboseLevel()

G4int G4RadioactiveDecay::GetVerboseLevel ( ) const
inline

Definition at line 173 of file G4RadioactiveDecay.hh.

173{return verboseLevel;}

Referenced by BuildPhysicsTable(), CalculateChainsFromParent(), ChooseCollimationDirection(), CollimateDecay(), CollimateDecayProduct(), ConvolveSourceTimeProfile(), DecayIt(), DeselectAllVolumes(), DeselectAVolume(), DoDecay(), G4RadioactiveDecay(), GetChainsFromParent(), GetDecayTime(), GetMeanFreePath(), GetMeanLifeTime(), LoadDecayTable(), SelectAllVolumes(), SelectAVolume(), SetDecayBias(), and SetSourceTimeProfile().

◆ IsAnalogueMonteCarlo()

G4bool G4RadioactiveDecay::IsAnalogueMonteCarlo ( )
inline

Definition at line 200 of file G4RadioactiveDecay.hh.

200{return AnalogueMC;}

◆ IsApplicable()

G4bool G4RadioactiveDecay::IsApplicable ( const G4ParticleDefinition &  aParticle)
virtual

Reimplemented from G4VProcess.

Definition at line 286 of file G4RadioactiveDecay.cc.

287{
288 // All particles other than G4Ions are rejected by default
289 if (((const G4Ions*)(&aParticle))->GetExcitationEnergy() > 0.) {
290 return true; // Not ground state - decay
291 }
292
293 if (aParticle.GetParticleName() == "GenericIon") {
294 return true;
295 } else if (!(aParticle.GetParticleType() == "nucleus")
296 || aParticle.GetPDGLifeTime() < 0. ) {
297 return false; // Nuclide is stable - no decay
298 }
299
300 // At this point nuclide must be an unstable ground state
301 // Determine whether it falls into the correct A and Z range
302 G4int A = ((const G4Ions*) (&aParticle))->GetAtomicMass();
303 G4int Z = ((const G4Ions*) (&aParticle))->GetAtomicNumber();
304
305 if (A > theNucleusLimits.GetAMax() || A < theNucleusLimits.GetAMin())
306 {return false;}
307 else if (Z > theNucleusLimits.GetZMax() || Z < theNucleusLimits.GetZMin())
308 {return false;}
309 return true;
310}
G4NucleusLimits::GetZMax
G4int GetZMax() const
Definition: G4NucleusLimits.hh:100
G4NucleusLimits::GetZMin
G4int GetZMin() const
Definition: G4NucleusLimits.hh:97
G4NucleusLimits::GetAMin
G4int GetAMin() const
Definition: G4NucleusLimits.hh:91
G4NucleusLimits::GetAMax
G4int GetAMax() const
Definition: G4NucleusLimits.hh:94
G4ParticleDefinition::GetParticleType
const G4String & GetParticleType() const
Definition: G4ParticleDefinition.hh:129

Referenced by CalculateChainsFromParent(), and DecayIt().

◆ IsRateTableReady()

G4bool G4RadioactiveDecay::IsRateTableReady ( const G4ParticleDefinition &  aParticle)

Definition at line 425 of file G4RadioactiveDecay.cc.

426{
427 // Check whether the radioactive decay rates table for the ion has already
428 // been calculated.
429 G4String aParticleName = aParticle.GetParticleName();
430 for (size_t i = 0; i < theParentChainTable.size(); i++) {
431 if (theParentChainTable[i].GetIonName() == aParticleName) return true;
432 }
433 return false;
434}

Referenced by DecayIt().

◆ LoadDecayTable()

G4DecayTable * G4RadioactiveDecay::LoadDecayTable ( const G4ParticleDefinition &  theParentNucleus)

Definition at line 805 of file G4RadioactiveDecay.cc.

806{
807 // Generate input data file name using Z and A of the parent nucleus
808 // file containing radioactive decay data.
809 G4int A = ((const G4Ions*)(&theParentNucleus))->GetAtomicMass();
810 G4int Z = ((const G4Ions*)(&theParentNucleus))->GetAtomicNumber();
811
812 G4double levelEnergy = ((const G4Ions*)(&theParentNucleus))->GetExcitationEnergy();
813 G4Ions::G4FloatLevelBase floatingLevel =
814 ((const G4Ions*)(&theParentNucleus))->GetFloatLevelBase();
815
816#ifdef G4MULTITHREADED
817 G4AutoLock lk(&G4RadioactiveDecay::radioactiveDecayMutex);
818
819 G4String key = theParentNucleus.GetParticleName();
820 DecayTableMap::iterator master_table_ptr = master_dkmap->find(key);
821
822 if (master_table_ptr != master_dkmap->end() ) { // If table is there
823 return master_table_ptr->second;
824 }
825#endif
826
827 //Check if data have been provided by the user
828 G4String file = theUserRadioactiveDataFiles[1000*A+Z];
829
830 if (file == "") {
831 std::ostringstream os;
832 os << dirPath << "/z" << Z << ".a" << A << '\0';
833 file = os.str();
834 }
835
836 G4DecayTable* theDecayTable = new G4DecayTable();
837 G4bool found(false); // True if energy level matches one in table
838
839 std::ifstream DecaySchemeFile;
840 DecaySchemeFile.open(file);
841
842 if (DecaySchemeFile.good()) {
843 // Initialize variables used for reading in radioactive decay data
844 G4bool floatMatch(false);
845 const G4int nMode = 12;
846 G4double modeTotalBR[nMode] = {0.0};
847 G4double modeSumBR[nMode];
848 for (G4int i = 0; i < nMode; i++) {
849 modeSumBR[i] = 0.0;
850 }
851
852 char inputChars[120]={' '};
853 G4String inputLine;
854 G4String recordType("");
855 G4String floatingFlag("");
856 G4String daughterFloatFlag("");
857 G4Ions::G4FloatLevelBase daughterFloatLevel;
858 G4RadioactiveDecayMode theDecayMode;
859 G4double decayModeTotal(0.0);
860 G4double parentExcitation(0.0);
861 G4double a(0.0);
862 G4double b(0.0);
863 G4double c(0.0);
864 G4double dummy(0.0);
865 G4BetaDecayType betaType(allowed);
866
867 // Loop through each data file record until you identify the decay
868 // data relating to the nuclide of concern.
869
870 G4bool complete(false); // bool insures only one set of values read for any
871 // given parent energy level
872 G4int loop = 0;
873 while (!complete && !DecaySchemeFile.getline(inputChars, 120).eof()) { /* Loop checking, 01.09.2015, D.Wright */
874 loop++;
875 if (loop > 100000) {
876 G4Exception("G4RadioactiveDecay::LoadDecayTable()", "HAD_RDM_100",
877 JustWarning, "While loop count exceeded");
878 break;
879 }
880
881 inputLine = inputChars;
882 inputLine = inputLine.strip(1);
883 if (inputChars[0] != '#' && inputLine.length() != 0) {
884 std::istringstream tmpStream(inputLine);
885
886 if (inputChars[0] == 'P') {
887 // Nucleus is a parent type. Check excitation level to see if it
888 // matches that of theParentNucleus
889 tmpStream >> recordType >> parentExcitation >> floatingFlag >> dummy;
890 // "dummy" takes the place of half-life
891 // Now read in from ENSDFSTATE in particle category
892
893 if (found) {
894 complete = true;
895 } else {
896 // Take first level which matches excitation energy regardless of floating level
897 found = (std::abs(parentExcitation*keV - levelEnergy) < levelTolerance);
898 if (floatingLevel != noFloat) {
899 // If floating level specificed, require match of both energy and floating level
900 floatMatch = (floatingLevel == G4Ions::FloatLevelBase(floatingFlag.back()) );
901 if (!floatMatch) found = false;
902 }
903 }
904
905 } else if (found) {
906 // The right part of the radioactive decay data file has been found. Search
907 // through it to determine the mode of decay of the subsequent records.
908
909 // Store for later the total decay probability for each decay mode
910 if (inputLine.length() < 72) {
911 tmpStream >> theDecayMode >> dummy >> decayModeTotal;
912 switch (theDecayMode) {
913 case IT:
914 {
915 G4ITDecay* anITChannel = new G4ITDecay(&theParentNucleus, decayModeTotal,
916 0.0, 0.0, photonEvaporation);
917 anITChannel->SetHLThreshold(halflifethreshold);
918 anITChannel->SetARM(applyARM);
919 theDecayTable->Insert(anITChannel);
920// anITChannel->DumpNuclearInfo();
921 }
922 break;
923 case BetaMinus:
924 modeTotalBR[1] = decayModeTotal; break;
925 case BetaPlus:
926 modeTotalBR[2] = decayModeTotal; break;
927 case KshellEC:
928 modeTotalBR[3] = decayModeTotal; break;
929 case LshellEC:
930 modeTotalBR[4] = decayModeTotal; break;
931 case MshellEC:
932 modeTotalBR[5] = decayModeTotal; break;
933 case NshellEC:
934 modeTotalBR[6] = decayModeTotal; break;
935 case Alpha:
936 modeTotalBR[7] = decayModeTotal; break;
937 case Proton:
938 modeTotalBR[8] = decayModeTotal; break;
939 case Neutron:
940 modeTotalBR[9] = decayModeTotal; break;
941 case BDProton:
942 break;
943 case BDNeutron:
944 break;
945 case Beta2Minus:
946 break;
947 case Beta2Plus:
948 break;
949 case Proton2:
950 break;
951 case Neutron2:
952 break;
953 case SpFission:
954 modeTotalBR[10] = decayModeTotal; break;
955 case Triton:
956 modeTotalBR[11] = decayModeTotal; break;
957 case RDM_ERROR:
958
959 default:
960 G4Exception("G4RadioactiveDecay::LoadDecayTable()", "HAD_RDM_000",
961 FatalException, "Selected decay mode does not exist");
962 } // switch
963
964 } else {
965 if (inputLine.length() < 84) {
966 tmpStream >> theDecayMode >> a >> daughterFloatFlag >> b >> c;
967 betaType = allowed;
968 } else {
969 tmpStream >> theDecayMode >> a >> daughterFloatFlag >> b >> c >> betaType;
970 }
971
972 // Allowed transitions are the default. Forbidden transitions are
973 // indicated in the last column.
974 a /= 1000.;
975 c /= 1000.;
976 b/=100.;
977 daughterFloatLevel = G4Ions::FloatLevelBase(daughterFloatFlag.back());
978
979 switch (theDecayMode) {
980 case BetaMinus:
981 {
982 G4BetaMinusDecay* aBetaMinusChannel =
983 new G4BetaMinusDecay(&theParentNucleus, b, c*MeV, a*MeV,
984 daughterFloatLevel, betaType);
985// aBetaMinusChannel->DumpNuclearInfo();
986 aBetaMinusChannel->SetHLThreshold(halflifethreshold);
987 theDecayTable->Insert(aBetaMinusChannel);
988 modeSumBR[1] += b;
989 }
990 break;
991
992 case BetaPlus:
993 {
994 G4BetaPlusDecay* aBetaPlusChannel =
995 new G4BetaPlusDecay(&theParentNucleus, b, c*MeV, a*MeV,
996 daughterFloatLevel, betaType);
997// aBetaPlusChannel->DumpNuclearInfo();
998 aBetaPlusChannel->SetHLThreshold(halflifethreshold);
999 theDecayTable->Insert(aBetaPlusChannel);
1000 modeSumBR[2] += b;
1001 }
1002 break;
1003
1004 case KshellEC: // K-shell electron capture
1005 {
1006 G4ECDecay* aKECChannel =
1007 new G4ECDecay(&theParentNucleus, b, c*MeV, a*MeV,
1008 daughterFloatLevel, KshellEC);
1009// aKECChannel->DumpNuclearInfo();
1010 aKECChannel->SetHLThreshold(halflifethreshold);
1011 aKECChannel->SetARM(applyARM);
1012 theDecayTable->Insert(aKECChannel);
1013 modeSumBR[3] += b;
1014 }
1015 break;
1016
1017 case LshellEC: // L-shell electron capture
1018 {
1019 G4ECDecay* aLECChannel =
1020 new G4ECDecay(&theParentNucleus, b, c*MeV, a*MeV,
1021 daughterFloatLevel, LshellEC);
1022// aLECChannel->DumpNuclearInfo();
1023 aLECChannel->SetHLThreshold(halflifethreshold);
1024 aLECChannel->SetARM(applyARM);
1025 theDecayTable->Insert(aLECChannel);
1026 modeSumBR[4] += b;
1027 }
1028 break;
1029
1030 case MshellEC: // M-shell electron capture
1031 {
1032 G4ECDecay* aMECChannel =
1033 new G4ECDecay(&theParentNucleus, b, c*MeV, a*MeV,
1034 daughterFloatLevel, MshellEC);
1035// aMECChannel->DumpNuclearInfo();
1036 aMECChannel->SetHLThreshold(halflifethreshold);
1037 aMECChannel->SetARM(applyARM);
1038 theDecayTable->Insert(aMECChannel);
1039 modeSumBR[5] += b;
1040 }
1041 break;
1042
1043 case NshellEC: // N-shell electron capture
1044 {
1045 G4ECDecay* aNECChannel =
1046 new G4ECDecay(&theParentNucleus, b, c*MeV, a*MeV,
1047 daughterFloatLevel, NshellEC);
1048// aNECChannel->DumpNuclearInfo();
1049 aNECChannel->SetHLThreshold(halflifethreshold);
1050 aNECChannel->SetARM(applyARM);
1051 theDecayTable->Insert(aNECChannel);
1052 modeSumBR[6] += b;
1053 }
1054 break;
1055
1056 case Alpha:
1057 {
1058 G4AlphaDecay* anAlphaChannel =
1059 new G4AlphaDecay(&theParentNucleus, b, c*MeV, a*MeV,
1060 daughterFloatLevel);
1061// anAlphaChannel->DumpNuclearInfo();
1062 anAlphaChannel->SetHLThreshold(halflifethreshold);
1063 theDecayTable->Insert(anAlphaChannel);
1064 modeSumBR[7] += b;
1065 }
1066 break;
1067
1068 case Proton:
1069 {
1070 G4ProtonDecay* aProtonChannel =
1071 new G4ProtonDecay(&theParentNucleus, b, c*MeV, a*MeV,
1072 daughterFloatLevel);
1073// aProtonChannel->DumpNuclearInfo();
1074 aProtonChannel->SetHLThreshold(halflifethreshold);
1075 theDecayTable->Insert(aProtonChannel);
1076 modeSumBR[8] += b;
1077 }
1078 break;
1079
1080 case Neutron:
1081 {
1082 G4NeutronDecay* aNeutronChannel =
1083 new G4NeutronDecay(&theParentNucleus, b, c*MeV, a*MeV,
1084 daughterFloatLevel);
1085// aNeutronChannel->DumpNuclearInfo();
1086 aNeutronChannel->SetHLThreshold(halflifethreshold);
1087 theDecayTable->Insert(aNeutronChannel);
1088 modeSumBR[9] += b;
1089 }
1090 break;
1091
1092 case BDProton:
1093 // Not yet implemented
1094 // G4cout << " beta-delayed proton decay, a = " << a << ", b = " << b << ", c = " << c << G4endl;
1095 break;
1096 case BDNeutron:
1097 // Not yet implemented
1098 // G4cout << " beta-delayed neutron decay, a = " << a << ", b = " << b << ", c = " << c << G4endl;
1099 break;
1100 case Beta2Minus:
1101 // Not yet implemented
1102 // G4cout << " Double beta- decay, a = " << a << ", b = " << b << ", c = " << c << G4endl;
1103 break;
1104 case Beta2Plus:
1105 // Not yet implemented
1106 // G4cout << " Double beta+ decay, a = " << a << ", b = " << b << ", c = " << c << G4endl;
1107 break;
1108 case Proton2:
1109 // Not yet implemented
1110 // G4cout << " Double proton decay, a = " << a << ", b = " << b << ", c = " << c << G4endl;
1111 break;
1112 case Neutron2:
1113 // Not yet implemented
1114 // G4cout << " Double beta- decay, a = " << a << ", b = " << b << ", c = " << c << G4endl;
1115 break;
1116 case SpFission:
1117 {
1118 G4SFDecay* aSpontFissChannel =
1119// new G4SFDecay(&theParentNucleus, decayModeTotal, 0.0, 0.0);
1120 new G4SFDecay(&theParentNucleus, b, c*MeV, a*MeV,
1121 daughterFloatLevel);
1122 theDecayTable->Insert(aSpontFissChannel);
1123 modeSumBR[10] += b;
1124 }
1125 break;
1126 case Triton:
1127 {
1128 G4TritonDecay* aTritonChannel =
1129 new G4TritonDecay(&theParentNucleus, b, c*MeV, a*MeV,
1130 daughterFloatLevel);
1131 // anTritonChannel->DumpNuclearInfo();
1132 aTritonChannel->SetHLThreshold(halflifethreshold);
1133 theDecayTable->Insert(aTritonChannel);
1134 modeSumBR[11] += b;
1135 }
1136 break;
1137
1138 case RDM_ERROR:
1139
1140 default:
1141 G4Exception("G4RadioactiveDecay::LoadDecayTable()", "HAD_RDM_000",
1142 FatalException, "Selected decay mode does not exist");
1143 } // switch
1144 } // line < 72
1145 } // if char == P
1146 } // if char != #
1147 } // While
1148
1149 // Go through the decay table and make sure that the branching ratios are
1150 // correctly normalised.
1151
1152 G4VDecayChannel* theChannel = 0;
1153 G4NuclearDecay* theNuclearDecayChannel = 0;
1154 G4String mode = "";
1155
1156 G4double theBR = 0.0;
1157 for (G4int i = 0; i < theDecayTable->entries(); i++) {
1158 theChannel = theDecayTable->GetDecayChannel(i);
1159 theNuclearDecayChannel = static_cast<G4NuclearDecay*>(theChannel);
1160 theDecayMode = theNuclearDecayChannel->GetDecayMode();
1161
1162 if (theDecayMode != IT) {
1163 theBR = theChannel->GetBR();
1164 theChannel->SetBR(theBR*modeTotalBR[theDecayMode]/modeSumBR[theDecayMode]);
1165 }
1166 }
1167 } // decay file exists
1168
1169 DecaySchemeFile.close();
1170
1171 if (!found && levelEnergy > 0) {
1172 // Case where IT cascade for excited isotopes has no entries in RDM database
1173 // Decay mode is isomeric transition.
1174 G4ITDecay* anITChannel = new G4ITDecay(&theParentNucleus, 1.0, 0.0, 0.0,
1175 photonEvaporation);
1176 anITChannel->SetHLThreshold(halflifethreshold);
1177 anITChannel->SetARM(applyARM);
1178 theDecayTable->Insert(anITChannel);
1179 }
1180
1181 if (theDecayTable && GetVerboseLevel() > 1) {
1182 theDecayTable->DumpInfo();
1183 }
1184
1185#ifdef G4MULTITHREADED
1186 //(*master_dkmap)[key] = theDecayTable; // store in master library
1187#endif
1188 return theDecayTable;
1189}
G4BetaDecayType
G4BetaDecayType
Definition: G4BetaDecayType.hh:39
G4ECDecay::SetARM
void SetARM(G4bool onoff)
Definition: G4ECDecay.hh:56
G4ITDecay::SetARM
void SetARM(G4bool onoff)
Definition: G4ITDecay.hh:59
G4Ions::FloatLevelBase
static G4Ions::G4FloatLevelBase FloatLevelBase(char flbChar)
Definition: G4Ions.cc:103
G4Ions::G4FloatLevelBase
G4FloatLevelBase
Definition: G4Ions.hh:83
G4NuclearDecay::SetHLThreshold
void SetHLThreshold(G4double HLT)
Definition: G4NuclearDecay.hh:61
G4String::strip
G4String strip(G4int strip_Type=trailing, char c=' ')
G4VDecayChannel::SetBR
void SetBR(G4double value)
Definition: G4VDecayChannel.cc:622

Referenced by GetDecayTable().

◆ ProcessDescription()

void G4RadioactiveDecay::ProcessDescription ( std::ostream &  outFile) const
virtual

Reimplemented from G4VProcess.

Definition at line 252 of file G4RadioactiveDecay.cc.

253{
254 outFile << "The radioactive decay process (G4RadioactiveDecay) handles the\n"
255 << "alpha, beta+, beta-, electron capture and isomeric transition\n"
256 << "decays of nuclei (G4GenericIon) with masses A > 4.\n"
257 << "The required half-lives and decay schemes are retrieved from\n"
258 << "the RadioactiveDecay database which was derived from ENSDF.\n";
259}

◆ SelectAllVolumes()

void G4RadioactiveDecay::SelectAllVolumes ( )

Definition at line 390 of file G4RadioactiveDecay.cc.

391{
392 G4LogicalVolumeStore* theLogicalVolumes;
393 G4LogicalVolume* volume;
394 theLogicalVolumes = G4LogicalVolumeStore::GetInstance();
395 ValidVolumes.clear();
396#ifdef G4VERBOSE
397 if (GetVerboseLevel()>1)
398 G4cout << " RDM Applies to all Volumes" << G4endl;
399#endif
400 for (size_t i = 0; i < theLogicalVolumes->size(); i++){
401 volume = (*theLogicalVolumes)[i];
402 ValidVolumes.push_back(volume->GetName());
403#ifdef G4VERBOSE
404 if (GetVerboseLevel()>1)
405 G4cout << " RDM Applies to Volume " << volume->GetName() << G4endl;
406#endif
407 }
408 std::sort(ValidVolumes.begin(), ValidVolumes.end());
409 // sort needed in order to allow binary_search
410 isAllVolumesMode=true;
411}

Referenced by G4RadioactiveDecay(), and G4RadioactiveDecaymessenger::SetNewValue().

◆ SelectAVolume()

void G4RadioactiveDecay::SelectAVolume ( const G4String  aVolume)

Definition at line 329 of file G4RadioactiveDecay.cc.

330{
331 G4LogicalVolumeStore* theLogicalVolumes;
332 G4LogicalVolume* volume;
333 theLogicalVolumes = G4LogicalVolumeStore::GetInstance();
334 for (size_t i = 0; i < theLogicalVolumes->size(); i++) {
335 volume=(*theLogicalVolumes)[i];
336 if (volume->GetName() == aVolume) {
337 ValidVolumes.push_back(aVolume);
338 std::sort(ValidVolumes.begin(), ValidVolumes.end());
339 // sort need for performing binary_search
340#ifdef G4VERBOSE
341 if (GetVerboseLevel()>1)
342 G4cout << " RDM Applies to : " << aVolume << G4endl;
343#endif
344 } else if(i == theLogicalVolumes->size()) {
345 G4ExceptionDescription ed;
346 ed << aVolume << " is not a valid logical volume name. Decay not activated for it."
347 << G4endl;
348 G4Exception("G4RadioactiveDecayBase::SelectAVolume()", "HAD_RDM_300",
349 JustWarning, ed);
350 }
351 }
352}

Referenced by G4RadioactiveDecaymessenger::SetNewValue().

◆ SetAnalogueMonteCarlo()

void G4RadioactiveDecay::SetAnalogueMonteCarlo ( G4bool  r)
inline

Definition at line 189 of file G4RadioactiveDecay.hh.

189 {
190 AnalogueMC = r;
191 if (!AnalogueMC) halflifethreshold = 1e-6*CLHEP::s;
192 }

Referenced by SetBRBias(), SetDecayBias(), SetSourceTimeProfile(), and SetSplitNuclei().

◆ SetARM()

void G4RadioactiveDecay::SetARM ( G4bool  arm)
inline

Definition at line 132 of file G4RadioactiveDecay.hh.

132{applyARM = arm;}

Referenced by G4RadioactiveDecaymessenger::SetNewValue().

◆ SetBRBias()

void G4RadioactiveDecay::SetBRBias ( G4bool  r)
inline

Definition at line 203 of file G4RadioactiveDecay.hh.

203 {
204 BRBias = r;
205 SetAnalogueMonteCarlo(0);
206 }
G4RadioactiveDecay::SetAnalogueMonteCarlo
void SetAnalogueMonteCarlo(G4bool r)
Definition: G4RadioactiveDecay.hh:189

Referenced by G4RadioactiveDecaymessenger::SetNewValue().

◆ SetDecayBias()

void G4RadioactiveDecay::SetDecayBias ( G4String  filename)

Definition at line 1605 of file G4RadioactiveDecay.cc.

1606{
1607
1608 std::ifstream infile(filename, std::ios::in);
1609 if (!infile) G4Exception("G4RadioactiveDecay::SetDecayBias()", "HAD_RDM_003",
1610 FatalException, "Unable to open bias data file" );
1611
1612 G4double bin, flux;
1613 G4int dWindows = 0;
1614 G4int i ;
1615
1616 theRadioactivityTables.clear();
1617
1618 NDecayBin = -1;
1619
1620 G4int loop = 0;
1621 while (infile >> bin >> flux ) { /* Loop checking, 01.09.2015, D.Wright */
1622 NDecayBin++;
1623 loop++;
1624 if (loop > 10000) {
1625 G4Exception("G4RadioactiveDecay::SetDecayBias()", "HAD_RDM_100",
1626 JustWarning, "While loop count exceeded");
1627 break;
1628 }
1629
1630 if (NDecayBin > 99) {
1631 G4Exception("G4RadioactiveDecay::SetDecayBias()", "HAD_RDM_004",
1632 FatalException, "Input bias file too big (>100 rows)" );
1633 } else {
1634 DBin[NDecayBin] = bin * s;
1635 DProfile[NDecayBin] = flux;
1636 if (flux > 0.) {
1637 decayWindows[NDecayBin] = dWindows;
1638 dWindows++;
1639 G4RadioactivityTable *rTable = new G4RadioactivityTable() ;
1640 theRadioactivityTables.push_back(rTable);
1641 }
1642 }
1643 }
1644 for ( i = 1; i<= NDecayBin; i++) DProfile[i] += DProfile[i-1];
1645 for ( i = 0; i<= NDecayBin; i++) DProfile[i] /= DProfile[NDecayBin];
1646 // converted to accumulated probabilities
1647
1648 SetAnalogueMonteCarlo(0);
1649 infile.close();
1650
1651#ifdef G4VERBOSE
1652 if (GetVerboseLevel() > 1)
1653 G4cout <<" Decay Bias Profile Nbin = " << NDecayBin <<G4endl;
1654#endif
1655}

Referenced by G4RadioactiveDecaymessenger::SetNewValue().

◆ SetDecayCollimation()

void G4RadioactiveDecay::SetDecayCollimation ( const G4ThreeVector &  theDir,
G4double  halfAngle = 0.*CLHEP::deg 
)
inline

Definition at line 233 of file G4RadioactiveDecay.hh.

234 {
235 SetDecayDirection(theDir);
236 SetDecayHalfAngle(halfAngle);
237 }
G4RadioactiveDecay::SetDecayHalfAngle
void SetDecayHalfAngle(G4double halfAngle=0.*CLHEP::deg)
Definition: G4RadioactiveDecay.hh:225
G4RadioactiveDecay::SetDecayDirection
void SetDecayDirection(const G4ThreeVector &theDir)
Definition: G4RadioactiveDecay.hh:217

◆ SetDecayDirection()

void G4RadioactiveDecay::SetDecayDirection ( const G4ThreeVector &  theDir)
inline

Definition at line 217 of file G4RadioactiveDecay.hh.

217 {
218 forceDecayDirection = theDir.unit();
219 }
CLHEP::Hep3Vector::unit
Hep3Vector unit() const

Referenced by SetDecayCollimation(), and G4RadioactiveDecaymessenger::SetNewValue().

◆ SetDecayHalfAngle()

void G4RadioactiveDecay::SetDecayHalfAngle ( G4double  halfAngle = 0.*CLHEP::deg)
inline

Definition at line 225 of file G4RadioactiveDecay.hh.

225 {
226 forceDecayHalfAngle = std::min(std::max(0.*CLHEP::deg,halfAngle),180.*CLHEP::deg);
227 }

Referenced by SetDecayCollimation(), and G4RadioactiveDecaymessenger::SetNewValue().

◆ SetDecayRate()

void G4RadioactiveDecay::SetDecayRate ( G4int  theZ,
G4int  theA,
G4double  theE,
G4int  theG,
std::vector< G4double >  theCoefficients,
std::vector< G4double >  theTaos 
)

Definition at line 1207 of file G4RadioactiveDecay.cc.

1211{
1212 //fill the decay rate vector
1213 ratesToDaughter.SetZ(theZ);
1214 ratesToDaughter.SetA(theA);
1215 ratesToDaughter.SetE(theE);
1216 ratesToDaughter.SetGeneration(theG);
1217 ratesToDaughter.SetDecayRateC(theCoefficients);
1218 ratesToDaughter.SetTaos(theTaos);
1219}
G4RadioactiveDecayRatesToDaughter::SetTaos
void SetTaos(std::vector< G4double > value)
Definition: G4RadioactiveDecayRatesToDaughter.hh:94
G4RadioactiveDecayRatesToDaughter::SetDecayRateC
void SetDecayRateC(std::vector< G4double > value)
Definition: G4RadioactiveDecayRatesToDaughter.hh:92

Referenced by CalculateChainsFromParent().

◆ SetFBeta()

void G4RadioactiveDecay::SetFBeta ( G4bool  r)
inline

Definition at line 196 of file G4RadioactiveDecay.hh.

196{ FBeta = r; }

Referenced by G4RadioactiveDecaymessenger::SetNewValue().

◆ SetHLThreshold()

void G4RadioactiveDecay::SetHLThreshold ( G4double  hl)
inline

Definition at line 126 of file G4RadioactiveDecay.hh.

126{halflifethreshold = hl;}

Referenced by G4RadioactiveDecaymessenger::SetNewValue().

◆ SetICM()

void G4RadioactiveDecay::SetICM ( G4bool  icm)
inline

Definition at line 129 of file G4RadioactiveDecay.hh.

129{applyICM = icm;}

Referenced by G4RadioactiveDecaymessenger::SetNewValue().

◆ SetNucleusLimits()

void G4RadioactiveDecay::SetNucleusLimits ( G4NucleusLimits  theNucleusLimits1)
inline

Definition at line 176 of file G4RadioactiveDecay.hh.

177 {theNucleusLimits = theNucleusLimits1 ;}

◆ SetSourceTimeProfile()

void G4RadioactiveDecay::SetSourceTimeProfile ( G4String  filename)

Definition at line 1557 of file G4RadioactiveDecay.cc.

1558{
1559 std::ifstream infile ( filename, std::ios::in );
1560 if (!infile) {
1561 G4ExceptionDescription ed;
1562 ed << " Could not open file " << filename << G4endl;
1563 G4Exception("G4RadioactiveDecay::SetSourceTimeProfile()", "HAD_RDM_001",
1564 FatalException, ed);
1565 }
1566
1567 G4double bin, flux;
1568 NSourceBin = -1;
1569
1570 G4int loop = 0;
1571 while (infile >> bin >> flux) { /* Loop checking, 01.09.2015, D.Wright */
1572 loop++;
1573 if (loop > 10000) {
1574 G4Exception("G4RadioactiveDecay::SetSourceTimeProfile()", "HAD_RDM_100",
1575 JustWarning, "While loop count exceeded");
1576 break;
1577 }
1578
1579 NSourceBin++;
1580 if (NSourceBin > 99) {
1581 G4Exception("G4RadioactiveDecay::SetSourceTimeProfile()", "HAD_RDM_002",
1582 FatalException, "Input source time file too big (>100 rows)");
1583
1584 } else {
1585 SBin[NSourceBin] = bin * s;
1586 SProfile[NSourceBin] = flux;
1587 }
1588 }
1589 SetAnalogueMonteCarlo(0);
1590 infile.close();
1591
1592#ifdef G4VERBOSE
1593 if (GetVerboseLevel() > 1)
1594 G4cout <<" Source Timeprofile Nbin = " << NSourceBin <<G4endl;
1595#endif
1596}

Referenced by G4RadioactiveDecaymessenger::SetNewValue().

◆ SetSplitNuclei()

void G4RadioactiveDecay::SetSplitNuclei ( G4int  r)
inline

Definition at line 209 of file G4RadioactiveDecay.hh.

209 {
210 NSplit = r;
211 SetAnalogueMonteCarlo(0);
212 }

Referenced by G4RadioactiveDecaymessenger::SetNewValue().

◆ SetVerboseLevel()

void G4RadioactiveDecay::SetVerboseLevel ( G4int  value)
inline

Definition at line 170 of file G4RadioactiveDecay.hh.

170{verboseLevel = value;}

Referenced by G4RadioactiveDecaymessenger::SetNewValue().

The documentation for this class was generated from the following files: