| Geant4 Cross Reference |
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25 //
26 ////////////////////////////////////////////////////////////////////////////////
27 //
28 // GEANT4 Class source file
29 //
30 // G4RadioactiveDecay
31 //
32 // Author: D.H. Wright (SLAC)
33 // Date: 29 August 2017
34 //
35 // Based on the code of F. Lei and P.R. Truscott.
36 //
37 ////////////////////////////////////////////////////////////////////////////////
38
39 #include "G4RadioactiveDecay.hh"
40 #include "G4RadioactivationMessenger.hh"
41
42 #include "G4SystemOfUnits.hh"
43 #include "G4DynamicParticle.hh"
44 #include "G4DecayProducts.hh"
45 #include "G4DecayTable.hh"
46 #include "G4ParticleChangeForRadDecay.hh"
47 #include "G4ITDecay.hh"
48 #include "G4BetaDecayType.hh"
49 #include "G4BetaMinusDecay.hh"
50 #include "G4BetaPlusDecay.hh"
51 #include "G4ECDecay.hh"
52 #include "G4AlphaDecay.hh"
53 #include "G4TritonDecay.hh"
54 #include "G4ProtonDecay.hh"
55 #include "G4NeutronDecay.hh"
56 #include "G4SFDecay.hh"
57 #include "G4VDecayChannel.hh"
58 #include "G4NuclearDecay.hh"
59 #include "G4RadioactiveDecayMode.hh"
60 #include "G4Fragment.hh"
61 #include "G4Ions.hh"
62 #include "G4IonTable.hh"
63 #include "G4BetaDecayType.hh"
64 #include "Randomize.hh"
65 #include "G4LogicalVolumeStore.hh"
66 #include "G4NuclearLevelData.hh"
67 #include "G4DeexPrecoParameters.hh"
68 #include "G4LevelManager.hh"
69 #include "G4ThreeVector.hh"
70 #include "G4Electron.hh"
71 #include "G4Positron.hh"
72 #include "G4Neutron.hh"
73 #include "G4Gamma.hh"
74 #include "G4Alpha.hh"
75 #include "G4Triton.hh"
76 #include "G4Proton.hh"
77
78 #include "G4HadronicProcessType.hh"
79 #include "G4HadronicProcessStore.hh"
80 #include "G4HadronicException.hh"
81 #include "G4LossTableManager.hh"
82 #include "G4VAtomDeexcitation.hh"
83 #include "G4UAtomicDeexcitation.hh"
84 #include "G4PhotonEvaporation.hh"
85
86 #include <vector>
87 #include <sstream>
88 #include <algorithm>
89 #include <fstream>
90
91 using namespace CLHEP;
92
93 G4RadioactiveDecay::G4RadioactiveDecay(const G4String& processName,
94 const G4double timeThresholdForRadioactiveDecays)
95 : G4VRadioactiveDecay(processName, timeThresholdForRadioactiveDecays)
96 {
97 #ifdef G4VERBOSE
98 if (GetVerboseLevel() > 1) {
99 G4cout << "G4RadioactiveDecay constructor: processName = " << processName
100 << G4endl;
101 }
102 #endif
103
104 theRadioactivationMessenger = new G4RadioactivationMessenger(this);
105
106 // Apply default values.
107 NSourceBin = 1;
108 SBin[0] = 0.* s;
109 SBin[1] = 1.* s; // Convert to ns
110 SProfile[0] = 1.;
111 SProfile[1] = 0.;
112 NDecayBin = 1;
113 DBin[0] = 0. * s ;
114 DBin[1] = 1. * s;
115 DProfile[0] = 1.;
116 DProfile[1] = 0.;
117 decayWindows[0] = 0;
118 G4RadioactivityTable* rTable = new G4RadioactivityTable() ;
119 theRadioactivityTables.push_back(rTable);
120 NSplit = 1;
121 AnalogueMC = true;
122 BRBias = true;
123 halflifethreshold = 1000.*nanosecond;
124 }
125
126 void G4RadioactiveDecay::ProcessDescription(std::ostream& outFile) const
127 {
128 outFile << "The G4RadioactiveDecay process performs radioactive decay of\n"
129 << "nuclides (G4GenericIon) in biased mode which includes nucleus\n"
130 << "duplication, branching ratio biasing, source time convolution\n"
131 << "and detector time convolution. It is designed for use in\n"
132 << "activation physics.\n"
133 << "The required half-lives and decay schemes are retrieved from\n"
134 << "the RadioactiveDecay database which was derived from ENSDF.\n";
135 }
136
137
138 G4RadioactiveDecay::~G4RadioactiveDecay()
139 {
140 delete theRadioactivationMessenger;
141 }
142
143 G4bool
144 G4RadioactiveDecay::IsRateTableReady(const G4ParticleDefinition& aParticle)
145 {
146 // Check whether the radioactive decay rates table for the ion has already
147 // been calculated.
148 G4String aParticleName = aParticle.GetParticleName();
149 for (std::size_t i = 0; i < theParentChainTable.size(); ++i) {
150 if (theParentChainTable[i].GetIonName() == aParticleName) return true;
151 }
152 return false;
153 }
154
155 void
156 G4RadioactiveDecay::GetChainsFromParent(const G4ParticleDefinition& aParticle)
157 {
158 // Retrieve the decay rate table for the specified aParticle
159 G4String aParticleName = aParticle.GetParticleName();
160
161 for (std::size_t i = 0; i < theParentChainTable.size(); ++i) {
162 if (theParentChainTable[i].GetIonName() == aParticleName) {
163 theDecayRateVector = theParentChainTable[i].GetItsRates();
164 }
165 }
166 #ifdef G4VERBOSE
167 if (GetVerboseLevel() > 1) {
168 G4cout << "The DecayRate Table for " << aParticleName << " is selected."
169 << G4endl;
170 }
171 #endif
172 }
173
174 // ConvolveSourceTimeProfile performs the convolution of the source time profile
175 // function with a single exponential characterized by a decay constant in the
176 // decay chain. The time profile is treated as a step function so that the
177 // convolution integral can be done bin-by-bin.
178 // This implements Eq. 4.13 of DERA technical note, with SProfile[i] = F(t')
179
180 G4double
181 G4RadioactiveDecay::ConvolveSourceTimeProfile(const G4double t, const G4double tau)
182 {
183 G4double convolvedTime = 0.0;
184 G4int nbin;
185 if ( t > SBin[NSourceBin]) {
186 nbin = NSourceBin;
187 } else {
188 nbin = 0;
189
190 G4int loop = 0;
191 while (t > SBin[nbin]) { // Loop checking, 01.09.2015, D.Wright
192 loop++;
193 if (loop > 1000) {
194 G4Exception("G4RadioactiveDecay::ConvolveSourceTimeProfile()",
195 "HAD_RDM_100", JustWarning, "While loop count exceeded");
196 break;
197 }
198 ++nbin;
199 }
200 --nbin;
201 }
202
203 // Use expm1 wherever possible to avoid large cancellation errors in
204 // 1 - exp(x) for small x
205 G4double earg = 0.0;
206 if (nbin > 0) {
207 for (G4int i = 0; i < nbin; ++i) {
208 earg = (SBin[i+1] - SBin[i])/tau;
209 if (earg < 100.) {
210 convolvedTime += SProfile[i] * std::exp((SBin[i] - t)/tau) *
211 std::expm1(earg);
212 } else {
213 convolvedTime += SProfile[i] *
214 (std::exp(-(t-SBin[i+1])/tau)-std::exp(-(t-SBin[i])/tau));
215 }
216 }
217 }
218 convolvedTime -= SProfile[nbin] * std::expm1((SBin[nbin] - t)/tau);
219 // tau divided out of final result to provide probability of decay in window
220
221 if (convolvedTime < 0.) {
222 G4cout << " Convolved time =: " << convolvedTime << " reset to zero! " << G4endl;
223 G4cout << " t = " << t << " tau = " << tau << G4endl;
224 G4cout << SBin[nbin] << " " << SBin[0] << G4endl;
225 convolvedTime = 0.;
226 }
227 #ifdef G4VERBOSE
228 if (GetVerboseLevel() > 2)
229 G4cout << " Convolved time: " << convolvedTime << G4endl;
230 #endif
231 return convolvedTime;
232 }
233
234
235 ////////////////////////////////////////////////////////////////////////////////
236 // //
237 // GetDecayTime //
238 // Randomly select a decay time for the decay process, following the //
239 // supplied decay time bias scheme. //
240 // //
241 ////////////////////////////////////////////////////////////////////////////////
242
243 G4double G4RadioactiveDecay::GetDecayTime()
244 {
245 G4double decaytime = 0.;
246 G4double rand = G4UniformRand();
247 G4int i = 0;
248
249 G4int loop = 0;
250 while (DProfile[i] < rand) { /* Loop checking, 01.09.2015, D.Wright */
251 // Entries in DProfile[i] are all between 0 and 1 and arranged in inreaseing order
252 // Comparison with rand chooses which time bin to sample
253 ++i;
254 loop++;
255 if (loop > 100000) {
256 G4Exception("G4RadioactiveDecay::GetDecayTime()", "HAD_RDM_100",
257 JustWarning, "While loop count exceeded");
258 break;
259 }
260 }
261
262 rand = G4UniformRand();
263 decaytime = DBin[i] + rand*(DBin[i+1]-DBin[i]);
264 #ifdef G4VERBOSE
265 if (GetVerboseLevel() > 2)
266 G4cout <<" Decay time: " <<decaytime/s <<"[s]" <<G4endl;
267 #endif
268 return decaytime;
269 }
270
271
272 G4int G4RadioactiveDecay::GetDecayTimeBin(const G4double aDecayTime)
273 {
274 G4int i = 0;
275
276 G4int loop = 0;
277 while (aDecayTime > DBin[i] ) { /* Loop checking, 01.09.2015, D.Wright */
278 ++i;
279 loop++;
280 if (loop > 100000) {
281 G4Exception("G4RadioactiveDecay::GetDecayTimeBin()", "HAD_RDM_100",
282 JustWarning, "While loop count exceeded");
283 break;
284 }
285 }
286
287 return i;
288 }
289
290 ////////////////////////////////////////////////////////////////////////////////
291 // //
292 // GetMeanLifeTime (required by the base class) //
293 // //
294 ////////////////////////////////////////////////////////////////////////////////
295
296 G4double G4RadioactiveDecay::GetMeanLifeTime(const G4Track& theTrack,
297 G4ForceCondition* fc)
298 {
299 // For variance reduction time is set to 0 so as to force the particle
300 // to decay immediately.
301 // In analogue mode it returns the particle's mean-life.
302 G4double meanlife = 0.;
303 if (AnalogueMC) meanlife = G4VRadioactiveDecay::GetMeanLifeTime(theTrack, fc);
304 return meanlife;
305 }
306
307
308 void
309 G4RadioactiveDecay::SetDecayRate(G4int theZ, G4int theA, G4double theE,
310 G4int theG, std::vector<G4double>& theCoefficients,
311 std::vector<G4double>& theTaos)
312 // Why not make this a method of G4RadioactiveDecayRate? (e.g. SetParameters)
313 {
314 //fill the decay rate vector
315 ratesToDaughter.SetZ(theZ);
316 ratesToDaughter.SetA(theA);
317 ratesToDaughter.SetE(theE);
318 ratesToDaughter.SetGeneration(theG);
319 ratesToDaughter.SetDecayRateC(theCoefficients);
320 ratesToDaughter.SetTaos(theTaos);
321 }
322
323
324 void G4RadioactiveDecay::
325 CalculateChainsFromParent(const G4ParticleDefinition& theParentNucleus)
326 {
327 // Use extended Bateman equation to calculate the radioactivities of all
328 // progeny of theParentNucleus. The coefficients required to do this are
329 // calculated using the method of P. Truscott (Ph.D. thesis and
330 // DERA Technical Note DERA/CIS/CIS2/7/36/4/10) 11 January 2000.
331 // Coefficients are then added to the decay rate table vector
332
333 // Create and initialise variables used in the method.
334 theDecayRateVector.clear();
335
336 G4int nGeneration = 0;
337
338 std::vector<G4double> taos;
339
340 // Dimensionless A coefficients of Eqs. 4.24 and 4.25 of the TN
341 std::vector<G4double> Acoeffs;
342
343 // According to Eq. 4.26 the first coefficient (A_1:1) is -1
344 Acoeffs.push_back(-1.);
345
346 const G4Ions* ion = static_cast<const G4Ions*>(&theParentNucleus);
347 G4int A = ion->GetAtomicMass();
348 G4int Z = ion->GetAtomicNumber();
349 G4double E = ion->GetExcitationEnergy();
350 G4double tao = ion->GetPDGLifeTime();
351 if (tao < 0.) tao = 1e-100;
352 taos.push_back(tao);
353 G4int nEntry = 0;
354
355 // Fill the decay rate container (G4RadioactiveDecayRate) with the parent
356 // isotope data
357 SetDecayRate(Z,A,E,nGeneration,Acoeffs,taos); // Fill TP with parent lifetime
358
359 // store the decay rate in decay rate vector
360 theDecayRateVector.push_back(ratesToDaughter);
361 ++nEntry;
362
363 // Now start treating the secondary generations.
364 G4bool stable = false;
365 G4int j;
366 G4VDecayChannel* theChannel = 0;
367 G4NuclearDecay* theNuclearDecayChannel = 0;
368
369 G4ITDecay* theITChannel = 0;
370 G4BetaMinusDecay* theBetaMinusChannel = 0;
371 G4BetaPlusDecay* theBetaPlusChannel = 0;
372 G4AlphaDecay* theAlphaChannel = 0;
373 G4ProtonDecay* theProtonChannel = 0;
374 G4TritonDecay* theTritonChannel = 0;
375 G4NeutronDecay* theNeutronChannel = 0;
376 G4SFDecay* theFissionChannel = 0;
377
378 G4RadioactiveDecayMode theDecayMode;
379 G4double theBR = 0.0;
380 G4int AP = 0;
381 G4int ZP = 0;
382 G4int AD = 0;
383 G4int ZD = 0;
384 G4double EP = 0.;
385 std::vector<G4double> TP;
386 std::vector<G4double> RP; // A coefficients of the previous generation
387 G4ParticleDefinition *theDaughterNucleus;
388 G4double daughterExcitation;
389 G4double nearestEnergy = 0.0;
390 G4int nearestLevelIndex = 0;
391 G4ParticleDefinition *aParentNucleus;
392 G4IonTable* theIonTable;
393 G4DecayTable* parentDecayTable;
394 G4double theRate;
395 G4double TaoPlus;
396 G4int nS = 0; // Running index of first decay in a given generation
397 G4int nT = nEntry; // Total number of decays accumulated over entire history
398 const G4int nMode = G4RadioactiveDecayModeSize;
399 G4double brs[nMode];
400 //
401 theIonTable = G4ParticleTable::GetParticleTable()->GetIonTable();
402
403 G4int loop = 0;
404 while (!stable) { /* Loop checking, 01.09.2015, D.Wright */
405 loop++;
406 if (loop > 10000) {
407 G4Exception("G4RadioactiveDecay::CalculateChainsFromParent()", "HAD_RDM_100",
408 JustWarning, "While loop count exceeded");
409 break;
410 }
411 nGeneration++;
412 for (j = nS; j < nT; ++j) {
413 // First time through, get data for parent nuclide
414 ZP = theDecayRateVector[j].GetZ();
415 AP = theDecayRateVector[j].GetA();
416 EP = theDecayRateVector[j].GetE();
417 RP = theDecayRateVector[j].GetDecayRateC();
418 TP = theDecayRateVector[j].GetTaos();
419 if (GetVerboseLevel() > 1) {
420 G4cout << "G4RadioactiveDecay::CalculateChainsFromParent: daughters of ("
421 << ZP << ", " << AP << ", " << EP
422 << ") are being calculated, generation = " << nGeneration
423 << G4endl;
424 }
425
426 aParentNucleus = theIonTable->GetIon(ZP,AP,EP);
427 parentDecayTable = GetDecayTable(aParentNucleus);
428 if (nullptr == parentDecayTable) { continue; }
429
430 G4DecayTable* summedDecayTable = new G4DecayTable();
431 // This instance of G4DecayTable is for accumulating BRs and decay
432 // channels. It will contain one decay channel per type of decay
433 // (alpha, beta, etc.); its branching ratio will be the sum of all
434 // branching ratios for that type of decay of the parent. If the
435 // halflife of a particular channel is longer than some threshold,
436 // that channel will be inserted specifically and its branching
437 // ratio will not be included in the above sums.
438 // This instance is not used to perform actual decays.
439
440 for (G4int k = 0; k < nMode; ++k) brs[k] = 0.0;
441
442 // Go through the decay table and sum all channels having the same decay mode
443 for (G4int i = 0; i < parentDecayTable->entries(); ++i) {
444 theChannel = parentDecayTable->GetDecayChannel(i);
445 theNuclearDecayChannel = static_cast<G4NuclearDecay*>(theChannel);
446 theDecayMode = theNuclearDecayChannel->GetDecayMode();
447 daughterExcitation = theNuclearDecayChannel->GetDaughterExcitation();
448 theDaughterNucleus = theNuclearDecayChannel->GetDaughterNucleus();
449 AD = ((const G4Ions*)(theDaughterNucleus))->GetAtomicMass();
450 ZD = ((const G4Ions*)(theDaughterNucleus))->GetAtomicNumber();
451 const G4LevelManager* levelManager =
452 G4NuclearLevelData::GetInstance()->GetLevelManager(ZD,AD);
453
454 // Check each nuclide to see if it is metastable (lifetime > 1 usec)
455 // If so, add it to the decay chain by inserting its decay channel in
456 // summedDecayTable. If not, just add its BR to sum for that decay mode.
457 if (levelManager->NumberOfTransitions() ) {
458 nearestEnergy = levelManager->NearestLevelEnergy(daughterExcitation);
459 if ((std::abs(daughterExcitation - nearestEnergy) < levelTolerance) &&
460 (std::abs(daughterExcitation - nearestEnergy) > DBL_EPSILON)) {
461 // Level half-life is in ns and the threshold is set to 1 micros
462 // by default, user can set it via the UI command
463 nearestLevelIndex = (G4int)levelManager->NearestLevelIndex(daughterExcitation);
464 if (levelManager->LifeTime(nearestLevelIndex)*ns >= halflifethreshold){
465 // save the metastable decay channel
466 summedDecayTable->Insert(theChannel);
467 } else {
468 brs[theDecayMode] += theChannel->GetBR();
469 }
470 } else {
471 brs[theDecayMode] += theChannel->GetBR();
472 }
473 } else {
474 brs[theDecayMode] += theChannel->GetBR();
475 }
476
477 } // Combine decay channels (loop i)
478
479 brs[BetaPlus] = brs[BetaPlus]+brs[KshellEC]+brs[LshellEC]+brs[MshellEC]+brs[NshellEC]; // Combine beta+ and EC
480 brs[KshellEC] = brs[LshellEC] = brs[MshellEC] = brs[NshellEC] = 0.0;
481 for (G4int i = 0; i < nMode; ++i) { // loop over decay modes
482 if (brs[i] > 0.) {
483 switch (i) {
484 case IT:
485 // Decay mode is isomeric transition
486 theITChannel = new G4ITDecay(aParentNucleus, brs[IT], 0.0, 0.0);
487
488 summedDecayTable->Insert(theITChannel);
489 break;
490
491 case BetaMinus:
492 // Decay mode is beta-
493 theBetaMinusChannel = new G4BetaMinusDecay(aParentNucleus, brs[BetaMinus],
494 0.*MeV, 0.*MeV,
495 noFloat, allowed);
496 summedDecayTable->Insert(theBetaMinusChannel);
497 break;
498
499 case BetaPlus:
500 // Decay mode is beta+ + EC.
501 theBetaPlusChannel = new G4BetaPlusDecay(aParentNucleus, brs[BetaPlus],
502 0.*MeV, 0.*MeV,
503 noFloat, allowed);
504 summedDecayTable->Insert(theBetaPlusChannel);
505 break;
506
507 case Alpha:
508 // Decay mode is alpha.
509 theAlphaChannel = new G4AlphaDecay(aParentNucleus, brs[Alpha], 0.*MeV,
510 0.*MeV, noFloat);
511 summedDecayTable->Insert(theAlphaChannel);
512 break;
513
514 case Proton:
515 // Decay mode is proton.
516 theProtonChannel = new G4ProtonDecay(aParentNucleus, brs[Proton], 0.*MeV,
517 0.*MeV, noFloat);
518 summedDecayTable->Insert(theProtonChannel);
519 break;
520
521 case Neutron:
522 // Decay mode is neutron.
523 theNeutronChannel = new G4NeutronDecay(aParentNucleus, brs[Neutron], 0.*MeV,
524 0.*MeV, noFloat);
525 summedDecayTable->Insert(theNeutronChannel);
526 break;
527
528 case SpFission:
529 // Decay mode is spontaneous fission
530 theFissionChannel = new G4SFDecay(aParentNucleus, brs[SpFission], 0.*MeV,
531 0.*MeV, noFloat);
532 summedDecayTable->Insert(theFissionChannel);
533 break;
534
535 case BDProton:
536 // Not yet implemented
537 break;
538
539 case BDNeutron:
540 // Not yet implemented
541 break;
542
543 case Beta2Minus:
544 // Not yet implemented
545 break;
546
547 case Beta2Plus:
548 // Not yet implemented
549 break;
550
551 case Proton2:
552 // Not yet implemented
553 break;
554
555 case Neutron2:
556 // Not yet implemented
557 break;
558
559 case Triton:
560 // Decay mode is Triton.
561 theTritonChannel = new G4TritonDecay(aParentNucleus, brs[Triton], 0.*MeV,
562 0.*MeV, noFloat);
563 summedDecayTable->Insert(theTritonChannel);
564 break;
565
566 default:
567 break;
568 }
569 }
570 }
571
572 // loop over all branches in summedDecayTable
573 //
574 for (G4int i = 0; i < summedDecayTable->entries(); ++i){
575 theChannel = summedDecayTable->GetDecayChannel(i);
576 theNuclearDecayChannel = static_cast<G4NuclearDecay*>(theChannel);
577 theBR = theChannel->GetBR();
578 theDaughterNucleus = theNuclearDecayChannel->GetDaughterNucleus();
579
580 // First check if the decay of the original nucleus is an IT channel,
581 // if true create a new ground-state nucleus
582 if (theNuclearDecayChannel->GetDecayMode() == IT && nGeneration == 1) {
583
584 A = ((const G4Ions*)(theDaughterNucleus))->GetAtomicMass();
585 Z = ((const G4Ions*)(theDaughterNucleus))->GetAtomicNumber();
586 theDaughterNucleus=theIonTable->GetIon(Z,A,0.);
587 }
588 if (IsApplicable(*theDaughterNucleus) && theBR > 0.0 &&
589 aParentNucleus != theDaughterNucleus) {
590 // need to make sure daughter has decay table
591 parentDecayTable = GetDecayTable(theDaughterNucleus);
592 if (nullptr != parentDecayTable && parentDecayTable->entries() > 0) {
593 A = ((const G4Ions*)(theDaughterNucleus))->GetAtomicMass();
594 Z = ((const G4Ions*)(theDaughterNucleus))->GetAtomicNumber();
595 E = ((const G4Ions*)(theDaughterNucleus))->GetExcitationEnergy();
596
597 TaoPlus = theDaughterNucleus->GetPDGLifeTime();
598 if (TaoPlus <= 0.) TaoPlus = 1e-100;
599
600 // first set the taos, one simply need to add to the parent ones
601 taos.clear();
602 taos = TP; // load lifetimes of all previous generations
603 std::size_t k;
604 //check that TaoPlus differs from other taos from at least 1.e5 relative difference
605 //for (k = 0; k < TP.size(); ++k){
606 //if (std::abs((TaoPlus-TP[k])/TP[k])<1.e-5 ) TaoPlus=1.00001*TP[k];
607 //}
608 taos.push_back(TaoPlus); // add daughter lifetime to list
609 // now calculate the coefficiencies
610 //
611 // they are in two parts, first the less than n ones
612 // Eq 4.24 of the TN
613 Acoeffs.clear();
614 long double ta1,ta2;
615 ta2 = (long double)TaoPlus;
616 for (k = 0; k < RP.size(); ++k){
617 ta1 = (long double)TP[k]; // loop over lifetimes of all previous generations
618 if (ta1 == ta2) {
619 theRate = 1.e100;
620 } else {
621 theRate = ta1/(ta1-ta2);
622 }
623 theRate = theRate * theBR * RP[k];
624 Acoeffs.push_back(theRate);
625 }
626
627 // the second part: the n:n coefficiency
628 // Eq 4.25 of the TN. Note Yn+1 is zero apart from Y1 which is -1
629 // as treated at line 1013
630 theRate = 0.;
631 long double aRate, aRate1;
632 aRate1 = 0.L;
633 for (k = 0; k < RP.size(); ++k){
634 ta1 = (long double)TP[k];
635 if (ta1 == ta2 ) {
636 aRate = 1.e100;
637 } else {
638 aRate = ta2/(ta1-ta2);
639 }
640 aRate = aRate * (long double)(theBR * RP[k]);
641 aRate1 += aRate;
642 }
643 theRate = -aRate1;
644 Acoeffs.push_back(theRate);
645 SetDecayRate (Z,A,E,nGeneration,Acoeffs,taos);
646 theDecayRateVector.push_back(ratesToDaughter);
647 nEntry++;
648 } // there are entries in the table
649 } // nuclide is OK to decay
650 } // end of loop (i) over decay table branches
651
652 delete summedDecayTable;
653
654 } // Getting contents of decay rate vector (end loop on j)
655 nS = nT;
656 nT = nEntry;
657 if (nS == nT) stable = true;
658 } // while nuclide is not stable
659
660 // end of while loop
661 // the calculation completed here
662
663
664 // fill the first part of the decay rate table
665 // which is the name of the original particle (isotope)
666 chainsFromParent.SetIonName(theParentNucleus.GetParticleName());
667
668 // now fill the decay table with the newly completed decay rate vector
669 chainsFromParent.SetItsRates(theDecayRateVector);
670
671 // finally add the decayratetable to the tablevector
672 theParentChainTable.push_back(chainsFromParent);
673 }
674
675 ////////////////////////////////////////////////////////////////////////////////
676 // //
677 // SetSourceTimeProfile //
678 // read in the source time profile function (histogram) //
679 // //
680 ////////////////////////////////////////////////////////////////////////////////
681
682 void G4RadioactiveDecay::SetSourceTimeProfile(const G4String& filename)
683 {
684 std::ifstream infile ( filename, std::ios::in );
685 if (!infile) {
686 G4ExceptionDescription ed;
687 ed << " Could not open file " << filename << G4endl;
688 G4Exception("G4RadioactiveDecay::SetSourceTimeProfile()", "HAD_RDM_001",
689 FatalException, ed);
690 }
691
692 G4double bin, flux;
693 NSourceBin = -1;
694
695 G4int loop = 0;
696 while (infile >> bin >> flux) { /* Loop checking, 01.09.2015, D.Wright */
697 loop++;
698 if (loop > 10000) {
699 G4Exception("G4RadioactiveDecay::SetSourceTimeProfile()", "HAD_RDM_100",
700 JustWarning, "While loop count exceeded");
701 break;
702 }
703
704 NSourceBin++;
705 if (NSourceBin > 99) {
706 G4Exception("G4RadioactiveDecay::SetSourceTimeProfile()", "HAD_RDM_002",
707 FatalException, "Input source time file too big (>100 rows)");
708
709 } else {
710 SBin[NSourceBin] = bin * s; // Convert read-in time to ns
711 SProfile[NSourceBin] = flux; // Dimensionless
712 }
713 }
714
715 AnalogueMC = false;
716 infile.close();
717
718 #ifdef G4VERBOSE
719 if (GetVerboseLevel() > 2)
720 G4cout <<" Source Timeprofile Nbin = " << NSourceBin <<G4endl;
721 #endif
722 }
723
724 ////////////////////////////////////////////////////////////////////////////////
725 // //
726 // SetDecayBiasProfile //
727 // read in the decay bias scheme function (histogram) //
728 // //
729 ////////////////////////////////////////////////////////////////////////////////
730
731 void G4RadioactiveDecay::SetDecayBias(const G4String& filename)
732 {
733 std::ifstream infile(filename, std::ios::in);
734 if (!infile) G4Exception("G4RadioactiveDecay::SetDecayBias()", "HAD_RDM_001",
735 FatalException, "Unable to open bias data file" );
736
737 G4double bin, flux;
738 G4int dWindows = 0;
739 G4int i ;
740
741 theRadioactivityTables.clear();
742
743 NDecayBin = -1;
744
745 G4int loop = 0;
746 while (infile >> bin >> flux ) { /* Loop checking, 01.09.2015, D.Wright */
747 NDecayBin++;
748 loop++;
749 if (loop > 10000) {
750 G4Exception("G4RadioactiveDecay::SetDecayBias()", "HAD_RDM_100",
751 JustWarning, "While loop count exceeded");
752 break;
753 }
754
755 if (NDecayBin > 99) {
756 G4Exception("G4RadioactiveDecay::SetDecayBias()", "HAD_RDM_002",
757 FatalException, "Input bias file too big (>100 rows)" );
758 } else {
759 DBin[NDecayBin] = bin * s; // Convert read-in time to ns
760 DProfile[NDecayBin] = flux; // Dimensionless
761 if (flux > 0.) {
762 decayWindows[NDecayBin] = dWindows;
763 dWindows++;
764 G4RadioactivityTable *rTable = new G4RadioactivityTable() ;
765 theRadioactivityTables.push_back(rTable);
766 }
767 }
768 }
769 for ( i = 1; i<= NDecayBin; ++i) DProfile[i] += DProfile[i-1]; // Cumulative flux vs i
770 for ( i = 0; i<= NDecayBin; ++i) DProfile[i] /= DProfile[NDecayBin];
771 // Normalize so entries increase from 0 to 1
772 // converted to accumulated probabilities
773
774 AnalogueMC = false;
775 infile.close();
776
777 #ifdef G4VERBOSE
778 if (GetVerboseLevel() > 2)
779 G4cout <<" Decay Bias Profile Nbin = " << NDecayBin <<G4endl;
780 #endif
781 }
782
783 ////////////////////////////////////////////////////////////////////////////////
784 // //
785 // DecayIt //
786 // //
787 ////////////////////////////////////////////////////////////////////////////////
788
789 G4VParticleChange*
790 G4RadioactiveDecay::DecayIt(const G4Track& theTrack, const G4Step&)
791 {
792 // Initialize G4ParticleChange object, get particle details and decay table
793 fParticleChangeForRadDecay.Initialize(theTrack);
794 fParticleChangeForRadDecay.ProposeWeight(theTrack.GetWeight());
795 const G4DynamicParticle* theParticle = theTrack.GetDynamicParticle();
796 const G4ParticleDefinition* theParticleDef = theParticle->GetDefinition();
797
798 // First check whether RDM applies to the current logical volume
799 if (!isAllVolumesMode) {
800 if (!std::binary_search(ValidVolumes.begin(), ValidVolumes.end(),
801 theTrack.GetVolume()->GetLogicalVolume()->GetName())) {
802 #ifdef G4VERBOSE
803 if (GetVerboseLevel()>0) {
804 G4cout <<"G4RadioactiveDecay::DecayIt : "
805 << theTrack.GetVolume()->GetLogicalVolume()->GetName()
806 << " is not selected for the RDM"<< G4endl;
807 G4cout << " There are " << ValidVolumes.size() << " volumes" << G4endl;
808 G4cout << " The Valid volumes are " << G4endl;
809 for (std::size_t i = 0; i< ValidVolumes.size(); ++i)
810 G4cout << ValidVolumes[i] << G4endl;
811 }
812 #endif
813 fParticleChangeForRadDecay.SetNumberOfSecondaries(0);
814
815 // Kill the parent particle.
816 fParticleChangeForRadDecay.ProposeTrackStatus(fStopAndKill) ;
817 fParticleChangeForRadDecay.ProposeLocalEnergyDeposit(0.0);
818 ClearNumberOfInteractionLengthLeft();
819 return &fParticleChangeForRadDecay;
820 }
821 }
822
823 // Now check if particle is valid for RDM
824 if (!(IsApplicable(*theParticleDef) ) ) {
825 // Particle is not an ion or is outside the nucleuslimits for decay
826 #ifdef G4VERBOSE
827 if (GetVerboseLevel() > 1) {
828 G4cout << "G4RadioactiveDecay::DecayIt : "
829 << theParticleDef->GetParticleName()
830 << " is not an ion or is outside (Z,A) limits set for the decay. "
831 << " Set particle change accordingly. "
832 << G4endl;
833 }
834 #endif
835 fParticleChangeForRadDecay.SetNumberOfSecondaries(0);
836
837 // Kill the parent particle
838 fParticleChangeForRadDecay.ProposeTrackStatus(fStopAndKill) ;
839 fParticleChangeForRadDecay.ProposeLocalEnergyDeposit(0.0);
840 ClearNumberOfInteractionLengthLeft();
841 return &fParticleChangeForRadDecay;
842 }
843
844 G4DecayTable* theDecayTable = GetDecayTable(theParticleDef);
845
846 if (theDecayTable == nullptr || theDecayTable->entries() == 0) {
847 // No data in the decay table. Set particle change parameters
848 // to indicate this.
849 #ifdef G4VERBOSE
850 if (GetVerboseLevel() > 1) {
851 G4cout << "G4RadioactiveDecay::DecayIt : "
852 << "decay table not defined for "
853 << theParticleDef->GetParticleName()
854 << ". Set particle change accordingly. "
855 << G4endl;
856 }
857 #endif
858 fParticleChangeForRadDecay.SetNumberOfSecondaries(0);
859
860 // Kill the parent particle.
861 fParticleChangeForRadDecay.ProposeTrackStatus(fStopAndKill) ;
862 fParticleChangeForRadDecay.ProposeLocalEnergyDeposit(0.0);
863 ClearNumberOfInteractionLengthLeft();
864 return &fParticleChangeForRadDecay;
865
866 } else {
867 // Data found. Try to decay nucleus
868 if (AnalogueMC) {
869 G4VRadioactiveDecay::DecayAnalog(theTrack, theDecayTable);
870
871 } else {
872 // Proceed with decay using variance reduction
873 G4double energyDeposit = 0.0;
874 G4double finalGlobalTime = theTrack.GetGlobalTime();
875 G4double finalLocalTime = theTrack.GetLocalTime();
876 G4int index;
877 G4ThreeVector currentPosition;
878 currentPosition = theTrack.GetPosition();
879
880 G4IonTable* theIonTable;
881 G4ParticleDefinition* parentNucleus;
882
883 // Get decay chains for the given nuclide
884 if (!IsRateTableReady(*theParticleDef))
885 CalculateChainsFromParent(*theParticleDef);
886 GetChainsFromParent(*theParticleDef);
887
888 // Declare some of the variables required in the implementation
889 G4int PZ;
890 G4int PA;
891 G4double PE;
892 G4String keyName;
893 std::vector<G4double> PT;
894 std::vector<G4double> PR;
895 G4double taotime;
896 long double decayRate;
897
898 std::size_t i;
899 G4int numberOfSecondaries;
900 G4int totalNumberOfSecondaries = 0;
901 G4double currentTime = 0.;
902 G4int ndecaych;
903 G4DynamicParticle* asecondaryparticle;
904 std::vector<G4DynamicParticle*> secondaryparticles;
905 std::vector<G4double> pw;
906 std::vector<G4double> ptime;
907 pw.clear();
908 ptime.clear();
909
910 // Now apply the nucleus splitting
911 for (G4int n = 0; n < NSplit; ++n) {
912 // Get the decay time following the decay probability function
913 // supplied by user
914 G4double theDecayTime = GetDecayTime();
915 G4int nbin = GetDecayTimeBin(theDecayTime);
916
917 // calculate the first part of the weight function
918 G4double weight1 = 1.;
919 if (nbin == 1) {
920 weight1 = 1./DProfile[nbin-1]
921 *(DBin[nbin]-DBin[nbin-1])/NSplit; // width of window in ns
922 } else if (nbin > 1) {
923 // Go from nbin to nbin-2 because flux entries in file alternate between 0 and 1
924 weight1 = 1./(DProfile[nbin]-DProfile[nbin-2])
925 *(DBin[nbin]-DBin[nbin-1])/NSplit;
926 // weight1 = (probability of choosing one of the bins)*(time width of bin)/NSplit
927 }
928 // it should be calculated in seconds
929 weight1 /= s ;
930
931 // loop over all the possible secondaries of the nucleus
932 // the first one is itself.
933 for (i = 0; i < theDecayRateVector.size(); ++i) {
934 PZ = theDecayRateVector[i].GetZ();
935 PA = theDecayRateVector[i].GetA();
936 PE = theDecayRateVector[i].GetE();
937 PT = theDecayRateVector[i].GetTaos();
938 PR = theDecayRateVector[i].GetDecayRateC();
939
940 // The array of arrays theDecayRateVector contains all possible decay
941 // chains of a given parent nucleus (ZP,AP,EP) to a given descendant
942 // nuclide (Z,A,E).
943 //
944 // theDecayRateVector[0] contains the decay parameters of the parent
945 // nucleus
946 // PZ = ZP
947 // PA = AP
948 // PE = EP
949 // PT[] = {TP}
950 // PR[] = {RP}
951 //
952 // theDecayRateVector[1] contains the decay of the parent to the first
953 // generation daughter (Z1,A1,E1).
954 // PZ = Z1
955 // PA = A1
956 // PE = E1
957 // PT[] = {TP, T1}
958 // PR[] = {RP, R1}
959 //
960 // theDecayRateVector[2] contains the decay of the parent to the first
961 // generation daughter (Z1,A1,E1) and the decay of the first
962 // generation daughter to the second generation daughter (Z2,A2,E2).
963 // PZ = Z2
964 // PA = A2
965 // PE = E2
966 // PT[] = {TP, T1, T2}
967 // PR[] = {RP, R1, R2}
968 //
969 // theDecayRateVector[3] may contain a branch chain
970 // PZ = Z2a
971 // PA = A2a
972 // PE = E2a
973 // PT[] = {TP, T1, T2a}
974 // PR[] = {RP, R1, R2a}
975 //
976 // and so on.
977
978 // Calculate the decay rate of the isotope. decayRate is the
979 // radioactivity of isotope (PZ,PA,PE) at 'theDecayTime'
980 // it will be used to calculate the statistical weight of the
981 // decay products of this isotope
982
983 // For each nuclide, calculate all the decay chains which can reach
984 // the parent nuclide
985 decayRate = 0.L;
986 for (G4int j = 0; j < G4int(PT.size() ); ++j) {
987 taotime = ConvolveSourceTimeProfile(theDecayTime,PT[j]);
988 decayRate -= PR[j] * (long double)taotime; // PRs are Acoeffs, taotime is inverse time
989 // Eq.4.23 of of the TN
990 // note the negative here is required as the rate in the
991 // equation is defined to be negative,
992 // i.e. decay away, but we need positive value here.
993
994 // G4cout << j << "\t"<< PT[j]/s << "\t" << PR[j] << "\t" << decayRate << G4endl;
995 }
996
997 // At this point any negative decay rates are probably small enough
998 // (order 10**-30) that negative values are likely due to cancellation
999 // errors. Set them to zero.
1000 if (decayRate < 0.0) decayRate = 0.0;
1001
1002 // G4cout <<theDecayTime/s <<"\t"<<nbin<<G4endl;
1003 // G4cout << theTrack.GetWeight() <<"\t"<<weight1<<"\t"<<decayRate<< G4endl;
1004
1005 // Add isotope to the radioactivity tables
1006 // One table for each observation time window specified in
1007 // SetDecayBias(G4String filename)
1008
1009 theRadioactivityTables[decayWindows[nbin-1]]
1010 ->AddIsotope(PZ,PA,PE,weight1*decayRate,theTrack.GetWeight());
1011
1012 // Now calculate the statistical weight
1013 // One needs to fold the source bias function with the decaytime
1014 // also need to include the track weight! (F.Lei, 28/10/10)
1015 G4double weight = weight1*decayRate*theTrack.GetWeight();
1016
1017 // decay the isotope
1018 theIonTable = (G4IonTable *)(G4ParticleTable::GetParticleTable()->GetIonTable());
1019 parentNucleus = theIonTable->GetIon(PZ,PA,PE);
1020
1021 // Create a temprary products buffer.
1022 // Its contents to be transfered to the products at the end of the loop
1023 G4DecayProducts* tempprods = nullptr;
1024
1025 // Decide whether to apply branching ratio bias or not
1026 if (BRBias) {
1027 G4DecayTable* decayTable = GetDecayTable(parentNucleus);
1028 G4VDecayChannel* theDecayChannel = nullptr;
1029 if (nullptr != decayTable) {
1030 ndecaych = G4int(decayTable->entries()*G4UniformRand());
1031 theDecayChannel = decayTable->GetDecayChannel(ndecaych);
1032 }
1033
1034 if (theDecayChannel == nullptr) {
1035 // Decay channel not found.
1036
1037 if (GetVerboseLevel() > 0) {
1038 G4cout << " G4RadioactiveDecay::DoIt : cannot determine decay channel ";
1039 G4cout << " for this nucleus; decay as if no biasing active. ";
1040 G4cout << G4endl;
1041 if (nullptr != decayTable) { decayTable ->DumpInfo(); }
1042 }
1043 // DHW 6 Dec 2010 - do decay as if no biasing to avoid deref of temppprods
1044 tempprods = DoDecay(*parentNucleus, theDecayTable);
1045 } else {
1046 // A decay channel has been identified, so execute the DecayIt.
1047 G4double tempmass = parentNucleus->GetPDGMass();
1048 tempprods = theDecayChannel->DecayIt(tempmass);
1049 weight *= (theDecayChannel->GetBR())*(decayTable->entries());
1050 }
1051 } else {
1052 tempprods = DoDecay(*parentNucleus, theDecayTable);
1053 }
1054
1055 // save the secondaries for buffers
1056 numberOfSecondaries = tempprods->entries();
1057 currentTime = finalGlobalTime + theDecayTime;
1058 for (index = 0; index < numberOfSecondaries; ++index) {
1059 asecondaryparticle = tempprods->PopProducts();
1060 if (asecondaryparticle->GetDefinition()->GetPDGStable() ) {
1061 pw.push_back(weight);
1062 ptime.push_back(currentTime);
1063 secondaryparticles.push_back(asecondaryparticle);
1064 }
1065 // Generate gammas and Xrays from excited nucleus, added by L.Desorgher
1066 else if (((const G4Ions*)(asecondaryparticle->GetDefinition()))
1067 ->GetExcitationEnergy() > 0. && weight > 0.) { //Compute the gamma
1068 G4ParticleDefinition* apartDef = asecondaryparticle->GetDefinition();
1069 AddDeexcitationSpectrumForBiasMode(apartDef,weight,currentTime,pw,
1070 ptime,secondaryparticles);
1071 }
1072 }
1073
1074 delete tempprods;
1075 } // end of i loop
1076 } // end of n loop
1077
1078 // now deal with the secondaries in the two stl containers
1079 // and submmit them back to the tracking manager
1080 totalNumberOfSecondaries = (G4int)pw.size();
1081 fParticleChangeForRadDecay.SetNumberOfSecondaries(totalNumberOfSecondaries);
1082 for (index=0; index < totalNumberOfSecondaries; ++index) {
1083 G4Track* secondary = new G4Track(secondaryparticles[index],
1084 ptime[index], currentPosition);
1085 secondary->SetGoodForTrackingFlag();
1086 secondary->SetTouchableHandle(theTrack.GetTouchableHandle());
1087 secondary->SetWeight(pw[index]);
1088 fParticleChangeForRadDecay.AddSecondary(secondary);
1089 }
1090
1091 // Kill the parent particle
1092 fParticleChangeForRadDecay.ProposeTrackStatus(fStopAndKill) ;
1093 fParticleChangeForRadDecay.ProposeLocalEnergyDeposit(energyDeposit);
1094 fParticleChangeForRadDecay.ProposeLocalTime(finalLocalTime);
1095 // Reset NumberOfInteractionLengthLeft.
1096 ClearNumberOfInteractionLengthLeft();
1097 } // end VR decay
1098
1099 return &fParticleChangeForRadDecay;
1100 } // end of data found branch
1101 }
1102
1103
1104 // Add gamma, X-ray, conversion and auger electrons for bias mode
1105 void
1106 G4RadioactiveDecay::AddDeexcitationSpectrumForBiasMode(G4ParticleDefinition* apartDef,
1107 G4double weight,G4double currentTime,
1108 std::vector<double>& weights_v,
1109 std::vector<double>& times_v,
1110 std::vector<G4DynamicParticle*>& secondaries_v)
1111 {
1112 G4double elevel=((const G4Ions*)(apartDef))->GetExcitationEnergy();
1113 G4double life_time=apartDef->GetPDGLifeTime();
1114 G4ITDecay* anITChannel = 0;
1115
1116 while (life_time < halflifethreshold && elevel > 0.) {
1117 decayIT->SetupDecay(apartDef);
1118 G4DecayProducts* pevap_products = decayIT->DecayIt(0.);
1119 G4int nb_pevapSecondaries = pevap_products->entries();
1120
1121 G4DynamicParticle* a_pevap_secondary = 0;
1122 G4ParticleDefinition* secDef = 0;
1123 for (G4int ind = 0; ind < nb_pevapSecondaries; ind++) {
1124 a_pevap_secondary= pevap_products->PopProducts();
1125 secDef = a_pevap_secondary->GetDefinition();
1126
1127 if (secDef->GetBaryonNumber() > 4) {
1128 elevel = ((const G4Ions*)(secDef))->GetExcitationEnergy();
1129 life_time = secDef->GetPDGLifeTime();
1130 apartDef = secDef;
1131 if (secDef->GetPDGStable() ) {
1132 weights_v.push_back(weight);
1133 times_v.push_back(currentTime);
1134 secondaries_v.push_back(a_pevap_secondary);
1135 }
1136 } else {
1137 weights_v.push_back(weight);
1138 times_v.push_back(currentTime);
1139 secondaries_v.push_back(a_pevap_secondary);
1140 }
1141 }
1142
1143 delete anITChannel;
1144 delete pevap_products;
1145 }
1146 }
1147
1148