| Geant4 Cross Reference |
1 // 1
2 // *******************************************
3 // * License and Disclaimer
4 // *
5 // * The Geant4 software is copyright of th
6 // * the Geant4 Collaboration. It is provided
7 // * conditions of the Geant4 Software License
8 // * LICENSE and available at http://cern.ch/
9 // * include a list of copyright holders.
10 // *
11 // * Neither the authors of this software syst
12 // * institutes,nor the agencies providing fin
13 // * work make any representation or warran
14 // * regarding this software system or assum
15 // * use. Please see the license in the file
16 // * for the full disclaimer and the limitatio
17 // *
18 // * This code implementation is the result
19 // * technical work of the GEANT4 collaboratio
20 // * By using, copying, modifying or distri
21 // * any work based on the software) you ag
22 // * use in resulting scientific publicati
23 // * acceptance of all terms of the Geant4 Sof
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. Trusc
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 G
94 const G
95 : G4VRadioactiveDecay(processName, timeThres
96 {
97 #ifdef G4VERBOSE
98 if (GetVerboseLevel() > 1) {
99 G4cout << "G4RadioactiveDecay constructor:
100 << G4endl;
101 }
102 #endif
103
104 theRadioactivationMessenger = new G4Radioact
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 G4Radioac
119 theRadioactivityTables.push_back(rTable);
120 NSplit = 1;
121 AnalogueMC = true;
122 BRBias = true;
123 halflifethreshold = 1000.*nanosecond;
124 }
125
126 void G4RadioactiveDecay::ProcessDescription(st
127 {
128 outFile << "The G4RadioactiveDecay process p
129 << "nuclides (G4GenericIon) in biase
130 << "duplication, branching ratio bia
131 << "and detector time convolution.
132 << "activation physics.\n"
133 << "The required half-lives and deca
134 << "the RadioactiveDecay database wh
135 }
136
137
138 G4RadioactiveDecay::~G4RadioactiveDecay()
139 {
140 delete theRadioactivationMessenger;
141 }
142
143 G4bool
144 G4RadioactiveDecay::IsRateTableReady(const G4P
145 {
146 // Check whether the radioactive decay rates
147 // been calculated.
148 G4String aParticleName = aParticle.GetPartic
149 for (std::size_t i = 0; i < theParentChainTa
150 if (theParentChainTable[i].GetIonName() ==
151 }
152 return false;
153 }
154
155 void
156 G4RadioactiveDecay::GetChainsFromParent(const
157 {
158 // Retrieve the decay rate table for the spe
159 G4String aParticleName = aParticle.GetPartic
160
161 for (std::size_t i = 0; i < theParentChainTa
162 if (theParentChainTable[i].GetIonName() ==
163 theDecayRateVector = theParentChainTable
164 }
165 }
166 #ifdef G4VERBOSE
167 if (GetVerboseLevel() > 1) {
168 G4cout << "The DecayRate Table for " << aP
169 << G4endl;
170 }
171 #endif
172 }
173
174 // ConvolveSourceTimeProfile performs the conv
175 // function with a single exponential characte
176 // decay chain. The time profile is treated a
177 // convolution integral can be done bin-by-bin
178 // This implements Eq. 4.13 of DERA technical
179
180 G4double
181 G4RadioactiveDecay::ConvolveSourceTimeProfile(
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
192 loop++;
193 if (loop > 1000) {
194 G4Exception("G4RadioactiveDecay::Convo
195 "HAD_RDM_100", JustWarning
196 break;
197 }
198 ++nbin;
199 }
200 --nbin;
201 }
202
203 // Use expm1 wherever possible to avoid larg
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::ex
211 std::expm1(earg);
212 } else {
213 convolvedTime += SProfile[i] *
214 (std::exp(-(t-SBin[i+1])/tau)-std::e
215 }
216 }
217 }
218 convolvedTime -= SProfile[nbin] * std::expm1
219 // tau divided out of final result to provid
220
221 if (convolvedTime < 0.) {
222 G4cout << " Convolved time =: " << convolv
223 G4cout << " t = " << t << " tau = " << tau
224 G4cout << SBin[nbin] << " " << SBin[0] <<
225 convolvedTime = 0.;
226 }
227 #ifdef G4VERBOSE
228 if (GetVerboseLevel() > 2)
229 G4cout << " Convolved time: " << convolved
230 #endif
231 return convolvedTime;
232 }
233
234
235 //////////////////////////////////////////////
236 //
237 // GetDecayTime
238 // Randomly select a decay time for the dec
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 checki
251 // Entries in DProfile[i] are all between
252 // Comparison with rand chooses which time
253 ++i;
254 loop++;
255 if (loop > 100000) {
256 G4Exception("G4RadioactiveDecay::GetDeca
257 JustWarning, "While loop cou
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 <<"
267 #endif
268 return decaytime;
269 }
270
271
272 G4int G4RadioactiveDecay::GetDecayTimeBin(cons
273 {
274 G4int i = 0;
275
276 G4int loop = 0;
277 while (aDecayTime > DBin[i] ) { /* Loop ch
278 ++i;
279 loop++;
280 if (loop > 100000) {
281 G4Exception("G4RadioactiveDecay::GetDeca
282 JustWarning, "While loop cou
283 break;
284 }
285 }
286
287 return i;
288 }
289
290 //////////////////////////////////////////////
291 //
292 // GetMeanLifeTime (required by the base clas
293 //
294 //////////////////////////////////////////////
295
296 G4double G4RadioactiveDecay::GetMeanLifeTime(c
297 G4
298 {
299 // For variance reduction time is set to 0 s
300 // to decay immediately.
301 // In analogue mode it returns the particle'
302 G4double meanlife = 0.;
303 if (AnalogueMC) meanlife = G4VRadioactiveDec
304 return meanlife;
305 }
306
307
308 void
309 G4RadioactiveDecay::SetDecayRate(G4int theZ, G
310 G4int theG, std::vector<G4double>& the
311 std::vector<G4double>& theTaos)
312 // Why not make this a method of G4Radioactiv
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(theCoefficient
320 ratesToDaughter.SetTaos(theTaos);
321 }
322
323
324 void G4RadioactiveDecay::
325 CalculateChainsFromParent(const G4ParticleDefi
326 {
327 // Use extended Bateman equation to calculat
328 // progeny of theParentNucleus. The coeffic
329 // calculated using the method of P. Truscot
330 // DERA Technical Note DERA/CIS/CIS2/7/36/4/
331 // Coefficients are then added to the decay
332
333 // Create and initialise variables used in t
334 theDecayRateVector.clear();
335
336 G4int nGeneration = 0;
337
338 std::vector<G4double> taos;
339
340 // Dimensionless A coefficients of Eqs. 4.24
341 std::vector<G4double> Acoeffs;
342
343 // According to Eq. 4.26 the first coefficie
344 Acoeffs.push_back(-1.);
345
346 const G4Ions* ion = static_cast<const G4Ions
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 (G4Radioact
356 // isotope data
357 SetDecayRate(Z,A,E,nGeneration,Acoeffs,taos)
358
359 // store the decay rate in decay rate vector
360 theDecayRateVector.push_back(ratesToDaughter
361 ++nEntry;
362
363 // Now start treating the secondary generati
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 coefficient
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 fir
397 G4int nT = nEntry; // Total number of deca
398 const G4int nMode = G4RadioactiveDecayModeSi
399 G4double brs[nMode];
400 //
401 theIonTable = G4ParticleTable::GetParticleTa
402
403 G4int loop = 0;
404 while (!stable) { /* Loop checking, 01.09.
405 loop++;
406 if (loop > 10000) {
407 G4Exception("G4RadioactiveDecay::Calcula
408 JustWarning, "While loop cou
409 break;
410 }
411 nGeneration++;
412 for (j = nS; j < nT; ++j) {
413 // First time through, get data for pare
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::Calcula
421 << ZP << ", " << AP << ", " <<
422 << ") are being calculated, ge
423 << G4endl;
424 }
425
426 aParentNucleus = theIonTable->GetIon(ZP,
427 parentDecayTable = GetDecayTable(aParent
428 if (nullptr == parentDecayTable) { conti
429
430 G4DecayTable* summedDecayTable = new G4D
431 // This instance of G4DecayTable is for
432 // channels. It will contain one decay
433 // (alpha, beta, etc.); its branching ra
434 // branching ratios for that type of dec
435 // halflife of a particular channel is l
436 // that channel will be inserted specifi
437 // ratio will not be included in the abo
438 // This instance is not used to perform
439
440 for (G4int k = 0; k < nMode; ++k) brs[k]
441
442 // Go through the decay table and sum al
443 for (G4int i = 0; i < parentDecayTable->
444 theChannel = parentDecayTable->GetDeca
445 theNuclearDecayChannel = static_cast<G
446 theDecayMode = theNuclearDecayChannel-
447 daughterExcitation = theNuclearDecayCh
448 theDaughterNucleus = theNuclearDecayCh
449 AD = ((const G4Ions*)(theDaughterNucle
450 ZD = ((const G4Ions*)(theDaughterNucle
451 const G4LevelManager* levelManager =
452 G4NuclearLevelData::GetInstance()->G
453
454 // Check each nuclide to see if it is
455 // If so, add it to the decay chain by
456 // summedDecayTable. If not, just add
457 if (levelManager->NumberOfTransitions(
458 nearestEnergy = levelManager->Neares
459 if ((std::abs(daughterExcitation - n
460 (std::abs(daughterExcitation - nearest
461 // Level half-life is in ns and th
462 // by default, user can set it via
463 nearestLevelIndex = (G4int)levelMa
464 if (levelManager->LifeTime(nearest
465 // save the metastable decay cha
466 summedDecayTable->Insert(theChan
467 } else {
468 brs[theDecayMode] += theChannel-
469 }
470 } else {
471 brs[theDecayMode] += theChannel->G
472 }
473 } else {
474 brs[theDecayMode] += theChannel->Get
475 }
476
477 } // Combine decay channels (loop i)
478
479 brs[BetaPlus] = brs[BetaPlus]+brs[Kshell
480 brs[KshellEC] = brs[LshellEC] = brs[Mshe
481 for (G4int i = 0; i < nMode; ++i) {
482 if (brs[i] > 0.) {
483 switch (i) {
484 case IT:
485 // Decay mode is isomeric transiti
486 theITChannel = new G4ITDecay(aPare
487
488 summedDecayTable->Insert(theITChan
489 break;
490
491 case BetaMinus:
492 // Decay mode is beta-
493 theBetaMinusChannel = new G4BetaMi
494
495
496 summedDecayTable->Insert(theBetaMi
497 break;
498
499 case BetaPlus:
500 // Decay mode is beta+ + EC.
501 theBetaPlusChannel = new G4BetaPlu
502
503
504 summedDecayTable->Insert(theBetaPl
505 break;
506
507 case Alpha:
508 // Decay mode is alpha.
509 theAlphaChannel = new G4AlphaDecay
510
511 summedDecayTable->Insert(theAlphaC
512 break;
513
514 case Proton:
515 // Decay mode is proton.
516 theProtonChannel = new G4ProtonDec
517
518 summedDecayTable->Insert(theProton
519 break;
520
521 case Neutron:
522 // Decay mode is neutron.
523 theNeutronChannel = new G4NeutronD
524
525 summedDecayTable->Insert(theNeutro
526 break;
527
528 case SpFission:
529 // Decay mode is spontaneous fissi
530 theFissionChannel = new G4SFDecay(
531
532 summedDecayTable->Insert(theFissio
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 G4TritonDec
562
563 summedDecayTable->Insert(theTriton
564 break;
565
566 default:
567 break;
568 }
569 }
570 }
571
572 // loop over all branches in summedDecay
573 //
574 for (G4int i = 0; i < summedDecayTable->
575 theChannel = summedDecayTable->GetDeca
576 theNuclearDecayChannel = static_cast<G
577 theBR = theChannel->GetBR();
578 theDaughterNucleus = theNuclearDecayCh
579
580 // First check if the decay of the ori
581 // if true create a new ground-state n
582 if (theNuclearDecayChannel->GetDecayMo
583
584 A = ((const G4Ions*)(theDaughterNucl
585 Z = ((const G4Ions*)(theDaughterNucl
586 theDaughterNucleus=theIonTable->GetI
587 }
588 if (IsApplicable(*theDaughterNucleus)
589 aParentNucleus != theDaughterNucle
590 // need to make sure daughter has de
591 parentDecayTable = GetDecayTable(the
592 if (nullptr != parentDecayTable && p
593 A = ((const G4Ions*)(theDaughterNu
594 Z = ((const G4Ions*)(theDaughterNu
595 E = ((const G4Ions*)(theDaughterNu
596
597 TaoPlus = theDaughterNucleus->GetP
598 if (TaoPlus <= 0.) TaoPlus = 1e-1
599
600 // first set the taos, one simply
601 taos.clear();
602 taos = TP; // load lifetimes of
603 std::size_t k;
604 //check that TaoPlus differs from
605 //for (k = 0; k < TP.size(); ++k){
606 //if (std::abs((TaoPlus-TP[k])/TP[
607 //}
608 taos.push_back(TaoPlus); // add d
609 // now calculate the coefficiencie
610 //
611 // they are in two parts, first th
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]; //
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 coeffi
628 // Eq 4.25 of the TN. Note Yn+1 i
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)(th
641 aRate1 += aRate;
642 }
643 theRate = -aRate1;
644 Acoeffs.push_back(theRate);
645 SetDecayRate (Z,A,E,nGeneration,Ac
646 theDecayRateVector.push_back(rates
647 nEntry++;
648 } // there are entries in the table
649 } // nuclide is OK to decay
650 } // end of loop (i) over decay table br
651
652 delete summedDecayTable;
653
654 } // Getting contents of decay rate vector
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 tab
665 // which is the name of the original particl
666 chainsFromParent.SetIonName(theParentNucleus
667
668 // now fill the decay table with the newly c
669 chainsFromParent.SetItsRates(theDecayRateVec
670
671 // finally add the decayratetable to the tab
672 theParentChainTable.push_back(chainsFromPare
673 }
674
675 //////////////////////////////////////////////
676 //
677 // SetSourceTimeProfile
678 // read in the source time profile functio
679 //
680 //////////////////////////////////////////////
681
682 void G4RadioactiveDecay::SetSourceTimeProfile(
683 {
684 std::ifstream infile ( filename, std::ios::i
685 if (!infile) {
686 G4ExceptionDescription ed;
687 ed << " Could not open file " << filename
688 G4Exception("G4RadioactiveDecay::SetSource
689 FatalException, ed);
690 }
691
692 G4double bin, flux;
693 NSourceBin = -1;
694
695 G4int loop = 0;
696 while (infile >> bin >> flux) { /* Loop che
697 loop++;
698 if (loop > 10000) {
699 G4Exception("G4RadioactiveDecay::SetSour
700 JustWarning, "While loop cou
701 break;
702 }
703
704 NSourceBin++;
705 if (NSourceBin > 99) {
706 G4Exception("G4RadioactiveDecay::SetSour
707 FatalException, "Input sourc
708
709 } else {
710 SBin[NSourceBin] = bin * s; // Conver
711 SProfile[NSourceBin] = flux; // Dimens
712 }
713 }
714
715 AnalogueMC = false;
716 infile.close();
717
718 #ifdef G4VERBOSE
719 if (GetVerboseLevel() > 2)
720 G4cout <<" Source Timeprofile Nbin = " <<
721 #endif
722 }
723
724 //////////////////////////////////////////////
725 //
726 // SetDecayBiasProfile
727 // read in the decay bias scheme function (
728 //
729 //////////////////////////////////////////////
730
731 void G4RadioactiveDecay::SetDecayBias(const G4
732 {
733 std::ifstream infile(filename, std::ios::in)
734 if (!infile) G4Exception("G4RadioactiveDecay
735 FatalException, "Un
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 ch
747 NDecayBin++;
748 loop++;
749 if (loop > 10000) {
750 G4Exception("G4RadioactiveDecay::SetDeca
751 JustWarning, "While loop cou
752 break;
753 }
754
755 if (NDecayBin > 99) {
756 G4Exception("G4RadioactiveDecay::SetDeca
757 FatalException, "Input bias
758 } else {
759 DBin[NDecayBin] = bin * s; // Conver
760 DProfile[NDecayBin] = flux; // Dimens
761 if (flux > 0.) {
762 decayWindows[NDecayBin] = dWindows;
763 dWindows++;
764 G4RadioactivityTable *rTable = new G4R
765 theRadioactivityTables.push_back(rTabl
766 }
767 }
768 }
769 for ( i = 1; i<= NDecayBin; ++i) DProfile[i]
770 for ( i = 0; i<= NDecayBin; ++i) DProfile[i]
771 // Normalize so e
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 = " <<
780 #endif
781 }
782
783 //////////////////////////////////////////////
784 //
785 // DecayIt
786 //
787 //////////////////////////////////////////////
788
789 G4VParticleChange*
790 G4RadioactiveDecay::DecayIt(const G4Track& the
791 {
792 // Initialize G4ParticleChange object, get p
793 fParticleChangeForRadDecay.Initialize(theTra
794 fParticleChangeForRadDecay.ProposeWeight(the
795 const G4DynamicParticle* theParticle = theTr
796 const G4ParticleDefinition* theParticleDef =
797
798 // First check whether RDM applies to the cu
799 if (!isAllVolumesMode) {
800 if (!std::binary_search(ValidVolumes.begin
801 theTrack.GetVolume()->Get
802 #ifdef G4VERBOSE
803 if (GetVerboseLevel()>0) {
804 G4cout <<"G4RadioactiveDecay::DecayIt
805 << theTrack.GetVolume()->GetLog
806 << " is not selected for the RD
807 G4cout << " There are " << ValidVolume
808 G4cout << " The Valid volumes are " <<
809 for (std::size_t i = 0; i< ValidVolume
810 G4cout << Va
811 }
812 #endif
813 fParticleChangeForRadDecay.SetNumberOfSe
814
815 // Kill the parent particle.
816 fParticleChangeForRadDecay.ProposeTrackS
817 fParticleChangeForRadDecay.ProposeLocalE
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 th
826 #ifdef G4VERBOSE
827 if (GetVerboseLevel() > 1) {
828 G4cout << "G4RadioactiveDecay::DecayIt :
829 << theParticleDef->GetParticleNam
830 << " is not an ion or is outside
831 << " Set particle change accordin
832 << G4endl;
833 }
834 #endif
835 fParticleChangeForRadDecay.SetNumberOfSeco
836
837 // Kill the parent particle
838 fParticleChangeForRadDecay.ProposeTrackSta
839 fParticleChangeForRadDecay.ProposeLocalEne
840 ClearNumberOfInteractionLengthLeft();
841 return &fParticleChangeForRadDecay;
842 }
843
844 G4DecayTable* theDecayTable = GetDecayTable(
845
846 if (theDecayTable == nullptr || theDecayTabl
847 // No data in the decay table. Set partic
848 // to indicate this.
849 #ifdef G4VERBOSE
850 if (GetVerboseLevel() > 1) {
851 G4cout << "G4RadioactiveDecay::DecayIt :
852 << "decay table not defined for "
853 << theParticleDef->GetParticleNam
854 << ". Set particle change accordi
855 << G4endl;
856 }
857 #endif
858 fParticleChangeForRadDecay.SetNumberOfSeco
859
860 // Kill the parent particle.
861 fParticleChangeForRadDecay.ProposeTrackSta
862 fParticleChangeForRadDecay.ProposeLocalEne
863 ClearNumberOfInteractionLengthLeft();
864 return &fParticleChangeForRadDecay;
865
866 } else {
867 // Data found. Try to decay nucleus
868 if (AnalogueMC) {
869 G4VRadioactiveDecay::DecayAnalog(theTrac
870
871 } else {
872 // Proceed with decay using variance red
873 G4double energyDeposit = 0.0;
874 G4double finalGlobalTime = theTrack.GetG
875 G4double finalLocalTime = theTrack.GetLo
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 nuclid
884 if (!IsRateTableReady(*theParticleDef))
885 CalculateChainsFromParent(*theParticleDef);
886 GetChainsFromParent(*theParticleDef);
887
888 // Declare some of the variables require
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*> secondar
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 de
913 // supplied by user
914 G4double theDecayTime = GetDecayTime()
915 G4int nbin = GetDecayTimeBin(theDecayT
916
917 // calculate the first part of the wei
918 G4double weight1 = 1.;
919 if (nbin == 1) {
920 weight1 = 1./DProfile[nbin-1]
921 *(DBin[nbin]-DBin[nbin-1])
922 } else if (nbin > 1) {
923 // Go from nbin to nbin-2 because fl
924 weight1 = 1./(DProfile[nbin]-DProfil
925 *(DBin[nbin]-DBin[nbin-1])
926 // weight1 = (probability of choosin
927 }
928 // it should be calculated in seconds
929 weight1 /= s ;
930
931 // loop over all the possible secondar
932 // the first one is itself.
933 for (i = 0; i < theDecayRateVector.siz
934 PZ = theDecayRateVector[i].GetZ();
935 PA = theDecayRateVector[i].GetA();
936 PE = theDecayRateVector[i].GetE();
937 PT = theDecayRateVector[i].GetTaos()
938 PR = theDecayRateVector[i].GetDecayR
939
940 // The array of arrays theDecayRateV
941 // chains of a given parent nucleus
942 // nuclide (Z,A,E).
943 //
944 // theDecayRateVector[0] contains th
945 // nucleus
946 // PZ = ZP
947 // PA = AP
948 // PE = EP
949 // PT[] = {TP}
950 // PR[] = {RP}
951 //
952 // theDecayRateVector[1] contains th
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 th
961 // generation daughter (Z1,A1,E1) an
962 // generation daughter to the second
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
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 i
979 // radioactivity of isotope (PZ,PA,P
980 // it will be used to calculate the
981 // decay products of this isotope
982
983 // For each nuclide, calculate all t
984 // the parent nuclide
985 decayRate = 0.L;
986 for (G4int j = 0; j < G4int(PT.size(
987 taotime = ConvolveSourceTimeProfil
988 decayRate -= PR[j] * (long double)
989 // Eq.4.23 of of the TN
990 // note the negative here is requi
991 // equation is defined to be negat
992 // i.e. decay away, but we need po
993
994 // G4cout << j << "\t"<< PT[j]/s <
995 }
996
997 // At this point any negative decay
998 // (order 10**-30) that negative val
999 // errors. Set them to zero.
1000 if (decayRate < 0.0) decayRate = 0.
1001
1002 // G4cout <<theDecayTime/s <<"\t"<
1003 // G4cout << theTrack.GetWeight()
1004
1005 // Add isotope to the radioactivity
1006 // One table for each observation t
1007 // SetDecayBias(G4String filename)
1008
1009 theRadioactivityTables[decayWindows
1010 ->AddIsotope(PZ,PA,PE,weight1
1011
1012 // Now calculate the statistical we
1013 // One needs to fold the source bia
1014 // also need to include the track w
1015 G4double weight = weight1*decayRate
1016
1017 // decay the isotope
1018 theIonTable = (G4IonTable *)(G4Part
1019 parentNucleus = theIonTable->GetIon
1020
1021 // Create a temprary products buffe
1022 // Its contents to be transfered to
1023 G4DecayProducts* tempprods = nullpt
1024
1025 // Decide whether to apply branchin
1026 if (BRBias) {
1027 G4DecayTable* decayTable = GetDec
1028 G4VDecayChannel* theDecayChannel = null
1029 if (nullptr != decayTable) {
1030 ndecaych = G4int(decayTable->entries(
1031 theDecayChannel = decayTable->G
1032 }
1033
1034 if (theDecayChannel == nullptr) {
1035 // Decay channel not found.
1036
1037 if (GetVerboseLevel() > 0) {
1038 G4cout << " G4RadioactiveDeca
1039 G4cout << " for this nucleus;
1040 G4cout << G4endl;
1041 if (nullptr != decayTable) {
1042 }
1043 // DHW 6 Dec 2010 - do decay as if no
1044 tempprods = DoDecay(*parentNucl
1045 } else {
1046 // A decay channel has been ide
1047 G4double tempmass = parentNucle
1048 tempprods = theDecayChannel->De
1049 weight *= (theDecayChannel->Get
1050 }
1051 } else {
1052 tempprods = DoDecay(*parentNucleu
1053 }
1054
1055 // save the secondaries for buffers
1056 numberOfSecondaries = tempprods->en
1057 currentTime = finalGlobalTime + the
1058 for (index = 0; index < numberOfSec
1059 asecondaryparticle = tempprods->P
1060 if (asecondaryparticle->GetDefini
1061 pw.push_back(weight);
1062 ptime.push_back(currentTime);
1063 secondaryparticles.push_back(as
1064 }
1065 // Generate gammas and Xrays from
1066 else if (((const G4Ions*)(asecond
1067 ->GetExcitationEnergy()
1068 G4ParticleDefinition* apartDef
1069 AddDeexcitationSpectrumForBiasM
1070
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
1079 // and submmit them back to the trackin
1080 totalNumberOfSecondaries = (G4int)pw.si
1081 fParticleChangeForRadDecay.SetNumberOfS
1082 for (index=0; index < totalNumberOfSeco
1083 G4Track* secondary = new G4Track(seco
1084 ptim
1085 secondary->SetGoodForTrackingFlag();
1086 secondary->SetTouchableHandle(theTrac
1087 secondary->SetWeight(pw[index]);
1088 fParticleChangeForRadDecay.AddSeconda
1089 }
1090
1091 // Kill the parent particle
1092 fParticleChangeForRadDecay.ProposeTrack
1093 fParticleChangeForRadDecay.ProposeLocal
1094 fParticleChangeForRadDecay.ProposeLocal
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 ele
1105 void
1106 G4RadioactiveDecay::AddDeexcitationSpectrumFo
1107 G4dou
1108 std::
1109 std::
1110 std::
1111 {
1112 G4double elevel=((const G4Ions*)(apartDef))
1113 G4double life_time=apartDef->GetPDGLifeTime
1114 G4ITDecay* anITChannel = 0;
1115
1116 while (life_time < halflifethreshold && ele
1117 decayIT->SetupDecay(apartDef);
1118 G4DecayProducts* pevap_products = decayIT
1119 G4int nb_pevapSecondaries = pevap_product
1120
1121 G4DynamicParticle* a_pevap_secondary = 0;
1122 G4ParticleDefinition* secDef = 0;
1123 for (G4int ind = 0; ind < nb_pevapSeconda
1124 a_pevap_secondary= pevap_products->PopP
1125 secDef = a_pevap_secondary->GetDefiniti
1126
1127 if (secDef->GetBaryonNumber() > 4) {
1128 elevel = ((const G4Ions*)(secDef))->G
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_sec
1135 }
1136 } else {
1137 weights_v.push_back(weight);
1138 times_v.push_back(currentTime);
1139 secondaries_v.push_back(a_pevap_secon
1140 }
1141 }
1142
1143 delete anITChannel;
1144 delete pevap_products;
1145 }
1146 }
1147
1148