| 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 // G4VRadioactiveDecay
31 //
32 // Author: D.H. Wright (SLAC)
33 // Date: 9 August 2017
34 //
35 // Based on the code of F. Lei and P.R. Trusc
36 //
37 //////////////////////////////////////////////
38
39 #include "G4VRadioactiveDecay.hh"
40 #include "G4RadioactiveDecayMessenger.hh"
41
42 #include "G4SystemOfUnits.hh"
43 #include "G4UnitsTable.hh"
44 #include "G4NuclideTable.hh"
45 #include "G4DynamicParticle.hh"
46 #include "G4DecayProducts.hh"
47 #include "G4DecayTable.hh"
48 #include "G4ParticleChangeForRadDecay.hh"
49 #include "G4ITDecay.hh"
50 #include "G4BetaDecayType.hh"
51 #include "G4BetaMinusDecay.hh"
52 #include "G4BetaPlusDecay.hh"
53 #include "G4ECDecay.hh"
54 #include "G4AlphaDecay.hh"
55 #include "G4TritonDecay.hh"
56 #include "G4ProtonDecay.hh"
57 #include "G4NeutronDecay.hh"
58 #include "G4SFDecay.hh"
59 #include "G4VDecayChannel.hh"
60 #include "G4NuclearDecay.hh"
61 #include "G4RadioactiveDecayMode.hh"
62 #include "G4Fragment.hh"
63 #include "G4Ions.hh"
64 #include "G4IonTable.hh"
65 #include "G4BetaDecayType.hh"
66 #include "Randomize.hh"
67 #include "G4LogicalVolumeStore.hh"
68 #include "G4NuclearLevelData.hh"
69 #include "G4DeexPrecoParameters.hh"
70 #include "G4LevelManager.hh"
71 #include "G4ThreeVector.hh"
72 #include "G4Electron.hh"
73 #include "G4Positron.hh"
74 #include "G4Neutron.hh"
75 #include "G4Gamma.hh"
76 #include "G4Alpha.hh"
77 #include "G4Triton.hh"
78 #include "G4Proton.hh"
79
80 #include "G4HadronicProcessType.hh"
81 #include "G4HadronicProcessStore.hh"
82 #include "G4HadronicException.hh"
83 #include "G4LossTableManager.hh"
84 #include "G4VAtomDeexcitation.hh"
85 #include "G4UAtomicDeexcitation.hh"
86 #include "G4PhotonEvaporation.hh"
87 #include "G4HadronicParameters.hh"
88
89 #include "G4PhysicsModelCatalog.hh"
90 #include "G4AutoLock.hh"
91
92 #include <vector>
93 #include <sstream>
94 #include <algorithm>
95 #include <fstream>
96
97 using namespace CLHEP;
98
99 const G4double G4VRadioactiveDecay::levelToler
100 const G4ThreeVector G4VRadioactiveDecay::origi
101
102 DecayTableMap* G4VRadioactiveDecay::master_dkm
103 std::map<G4int, G4String>* G4VRadioactiveDecay
104 G4String G4VRadioactiveDecay::dirPath = "";
105
106 namespace
107 {
108 G4Mutex radioactiveDecayMutex = G4MUTEX_INIT
109 }
110
111 G4VRadioactiveDecay::G4VRadioactiveDecay(const
112 const
113 : G4VRestDiscreteProcess(processName, fDecay)
114 fThresholdForVeryLongDecayTime( 1.0*CLHEP::
115 {
116 if (GetVerboseLevel() > 1) {
117 G4cout << "G4VRadioactiveDecay constructor
118 << G4endl;
119 }
120
121 SetProcessSubType(fRadioactiveDecay);
122
123 theRadioactiveDecayMessenger = new G4Radioac
124 pParticleChange = &fParticleChangeForRadDeca
125
126 // Check data directory
127 if (dirPath.empty()) {
128 const char* path_var = G4FindDataDir("G4RA
129 if (nullptr == path_var) {
130 G4Exception("G4VRadioactiveDecay()", "HA
131 "Environment variable G4RADI
132 } else {
133 dirPath = path_var; // convert to stri
134 std::ostringstream os;
135 os << dirPath << "/z1.a3"; // used as
136 std::ifstream testFile;
137 testFile.open(os.str() );
138 if ( !testFile.is_open() )
139 G4Exception("G4VRadioactiveDecay()","H
140 "Environment variable G4RA
141 }
142 }
143 // Set up photon evaporation for use in G4IT
144 photonEvaporation = new G4PhotonEvaporation(
145 photonEvaporation->RDMForced(true);
146 photonEvaporation->SetICM(true);
147 decayIT = new G4ITDecay(photonEvaporation);
148
149 // Instantiate the map of decay tables
150 if (nullptr == master_dkmap) {
151 master_dkmap = new DecayTableMap();
152 }
153 if (nullptr == theUserRDataFiles) {
154 theUserRDataFiles = new std::map<G4int, G4
155 }
156
157 // RDM applies to all logical volumes by def
158 SelectAllVolumes();
159 G4HadronicProcessStore::Instance()->Register
160
161 // The time threshold for radioactive decays
162 // 1. Via C++ interface: G4HadronicParameter
163 // 2. Via the second parameter of the G4VRad
164 // 3. Via UI command: /process/had/rdm/thres
165 // If both 1. and 2. are specified (at the m
166 // then we take the larger value, to be cons
167 // If, later on (after invoking the G4VRadio
168 // then this value is used (and the eventual
169 G4double timeThresholdBis = G4HadronicParame
170 if ( timeThreshold > 0.0 || timeThresholdBis
171 if ( timeThreshold > timeThresholdBis ) ti
172 fThresholdForVeryLongDecayTime = timeThres
173 }
174 }
175
176 G4VRadioactiveDecay::~G4VRadioactiveDecay()
177 {
178 delete theRadioactiveDecayMessenger;
179 delete photonEvaporation;
180 delete decayIT;
181 if (nullptr != master_dkmap) {
182 G4AutoLock lk(&radioactiveDecayMutex);
183 if (nullptr != master_dkmap) {
184 for (auto const & i : *master_dkmap) {
185 delete i.second;
186 }
187 master_dkmap->clear();
188 delete master_dkmap;
189 master_dkmap = nullptr;
190 }
191 delete theUserRDataFiles;
192 theUserRDataFiles = nullptr;
193 lk.unlock();
194 }
195 }
196
197
198 G4bool G4VRadioactiveDecay::IsApplicable(const
199 {
200 const G4String& pname = aParticle.GetParticl
201 if (pname == "GenericIon" || pname == "trito
202 // All particles other than G4Ions, are reje
203 const G4Ions* p = dynamic_cast<const G4Ions*
204 if (nullptr == p) { return false; }
205
206 // excited isomere may decay via gamma evapo
207 if (p->GetExcitationEnergy() > 0.0) { return
208
209 // Check on life time
210 G4double lifeTime = p->GetPDGLifeTime();
211 if (lifeTime < 0.0 || lifeTime > fThresholdF
212 return false;
213 }
214
215 // Determine whether the nuclide falls into
216 G4int A = p->GetAtomicMass();
217 G4int Z = p->GetAtomicNumber();
218
219 if (A > theNucleusLimits.GetAMax() || A < th
220 Z > theNucleusLimits.GetZMax() || Z < th
221 return false;
222 }
223
224 return true;
225 }
226
227
228 G4VParticleChange* G4VRadioactiveDecay::AtRest
229
230 {
231 return DecayIt(theTrack, theStep);
232 }
233
234
235 G4VParticleChange* G4VRadioactiveDecay::PostSt
236
237 {
238 return DecayIt(theTrack, theStep);
239 }
240
241
242 void G4VRadioactiveDecay::ProcessDescription(s
243 {
244 outFile << "The radioactive decay process (G
245 << "alpha, beta+, beta-, electron ca
246 << "decays of nuclei (G4GenericIon)
247 << "The required half-lives and deca
248 << "the RadioactiveDecay database wh
249 }
250
251
252
253
254 G4DecayTable* G4VRadioactiveDecay::GetDecayTab
255 {
256 G4String key = aNucleus->GetParticleName();
257 auto ptr = master_dkmap->find(key);
258
259 G4DecayTable* theDecayTable = nullptr;
260 if ( ptr == master_dkmap->end() ) {
261 // Load new file if table not there
262 const G4Ions* ion = dynamic_cast<const G4I
263 if (nullptr != ion) {
264 theDecayTable = LoadDecayTable(ion);
265 }
266 } else {
267 theDecayTable = ptr->second;
268 }
269 return theDecayTable;
270 }
271
272
273 void G4VRadioactiveDecay::SelectAVolume(const
274 {
275 G4LogicalVolumeStore* theLogicalVolumes = G4
276 G4LogicalVolume* volume = nullptr;
277 volume = theLogicalVolumes->GetVolume(aVolum
278 if (volume != nullptr)
279 {
280 ValidVolumes.push_back(aVolume);
281 std::sort(ValidVolumes.begin(), ValidVolum
282 // sort need for performing binary_search
283
284 if (GetVerboseLevel() > 0)
285 G4cout << " Radioactive decay applied to
286 }
287 else
288 {
289 G4ExceptionDescription ed;
290 ed << aVolume << " is not a valid logical
291 << " Decay not activated for it."
292 << G4endl;
293 G4Exception("G4VRadioactiveDecay::SelectAV
294 JustWarning, ed);
295 }
296 }
297
298
299 void G4VRadioactiveDecay::DeselectAVolume(cons
300 {
301 G4LogicalVolumeStore* theLogicalVolumes = G4
302 G4LogicalVolume* volume = nullptr;
303 volume = theLogicalVolumes->GetVolume(aVolum
304 if (volume != nullptr)
305 {
306 auto location= std::find(ValidVolumes.cbeg
307 if (location != ValidVolumes.cend() )
308 {
309 ValidVolumes.erase(location);
310 std::sort(ValidVolumes.begin(), ValidVol
311 isAllVolumesMode = false;
312 if (GetVerboseLevel() > 0)
313 G4cout << " G4VRadioactiveDecay::Desel
314 << " is removed from list " <<
315 }
316 else
317 {
318 G4ExceptionDescription ed;
319 ed << aVolume << " is not in the list.
320 G4Exception("G4VRadioactiveDecay::Desele
321 JustWarning, ed);
322 }
323 }
324 else
325 {
326 G4ExceptionDescription ed;
327 ed << aVolume << " is not a valid logical
328 << G4endl;
329 G4Exception("G4VRadioactiveDecay::Deselect
330 JustWarning, ed);
331 }
332 }
333
334
335 void G4VRadioactiveDecay::SelectAllVolumes()
336 {
337 G4LogicalVolumeStore* theLogicalVolumes = G4
338 G4LogicalVolume* volume = nullptr;
339 ValidVolumes.clear();
340 #ifdef G4VERBOSE
341 if (GetVerboseLevel()>1)
342 G4cout << " RDM Applies to all Volumes" <
343 #endif
344 for (std::size_t i = 0; i < theLogicalVolume
345 volume = (*theLogicalVolumes)[i];
346 ValidVolumes.push_back(volume->GetName());
347 #ifdef G4VERBOSE
348 if (GetVerboseLevel()>1)
349 G4cout << " RDM Applies to Volume
350 #endif
351 }
352 std::sort(ValidVolumes.begin(), ValidVolumes
353 // sort needed in order to allow binary_sear
354 isAllVolumesMode=true;
355 }
356
357
358 void G4VRadioactiveDecay::DeselectAllVolumes()
359 {
360 ValidVolumes.clear();
361 isAllVolumesMode=false;
362 #ifdef G4VERBOSE
363 if (GetVerboseLevel() > 1) G4cout << "RDM re
364 #endif
365 }
366
367
368 // GetMeanLifeTime (required by the base clas
369 G4double G4VRadioactiveDecay::GetMeanLifeTime(
370 G
371 {
372 G4double meanlife = DBL_MAX;
373 const G4ParticleDefinition* theParticleDef =
374 if (!IsApplicable(*theParticleDef)) { return
375 G4double theLife = theParticleDef->GetPDGLif
376 #ifdef G4VERBOSE
377 if (GetVerboseLevel() > 2) {
378 G4cout << "G4VRadioactiveDecay::GetMeanLif
379 << theParticleDef->GetParticleName() << G
380 G4cout << "KineticEnergy(GeV)=" << theTrac
381 << " Mass(GeV)=" << theParticleDef-
382 << " LifeTime(ns)=" << theLife/CLHE
383 }
384 #endif
385 if (theLife >= 0.0 && theLife <= fThresholdF
386 meanlife = theLife;
387 }
388
389 if (meanlife == DBL_MAX) {
390 const G4Ions* ion = dynamic_cast<const G4I
391 if (nullptr != ion && ion->GetExcitationEn
392 meanlife = 0.0;
393 }
394 }
395
396 #ifdef G4VERBOSE
397 if (GetVerboseLevel() > 2)
398 G4cout << "G4VRadioactiveDecay::GetMeanLif
399 << meanlife/CLHEP::s << " second " << G4e
400 #endif
401
402 return meanlife;
403 }
404
405
406 //////////////////////////////////////////////
407 //
408 // GetMeanFreePath for decay in flight
409 //
410 //////////////////////////////////////////////
411
412 G4double G4VRadioactiveDecay::GetMeanFreePath(
413 G
414 {
415 G4double res = DBL_MAX;
416 G4double lifeTime = GetMeanLifeTime(aTrack,
417 if (lifeTime > 0.0 && lifeTime < DBL_MAX) {
418 auto dParticle = aTrack.GetDynamicParticle
419 res = lifeTime*dParticle->GetTotalEnergy()
420 } else {
421 res = lifeTime;
422 }
423
424 #ifdef G4VERBOSE
425 if (GetVerboseLevel() > 2) {
426 G4cout << "G4VRadioactiveDecay::GetMeanFre
427 << aTrack.GetDefinition()->GetParticleNam
428 G4cout << " kinEnergy(GeV)=" << aTrack.Ge
429 << " lifeTime(ns)=" << lifeTime
430 << " mean free path(cm)=" << res/CL
431 }
432 #endif
433 return res;
434 }
435
436 //////////////////////////////////////////////
437 //
438 // BuildPhysicsTable - initialization of atom
439 //
440 //////////////////////////////////////////////
441
442 void G4VRadioactiveDecay::BuildPhysicsTable(co
443 {
444 if (isInitialised) { return; }
445 isInitialised = true;
446 if (G4HadronicParameters::Instance()->GetVer
447 G4Threading::IsMasterThread() && "Generi
448 StreamInfo(G4cout, "\n");
449 }
450 photonEvaporation->Initialise();
451 photonEvaporation->RDMForced(true);
452 photonEvaporation->SetICM(true);
453 decayIT->SetARM(applyARM);
454
455 G4HadronicProcessStore::Instance()->Register
456 G4HadronicProcessStore::Instance()->PrintInf
457 }
458
459 //////////////////////////////////////////////
460 //
461 // StreamInfo - stream out parameters
462 //
463 //////////////////////////////////////////////
464
465 void
466 G4VRadioactiveDecay::StreamInfo(std::ostream&
467 {
468 G4DeexPrecoParameters* deex =
469 G4NuclearLevelData::GetInstance()->GetPara
470 G4EmParameters* emparam = G4EmParameters::In
471 G4double minMeanLife = G4NuclideTable::GetIn
472
473 G4long prec = os.precision(5);
474 os << "=====================================
475 << endline;
476 os << "====== Radioactive Decay Phy
477 << endline;
478 os << "=====================================
479 << endline;
480 os << "min MeanLife (from G4NuclideTable)
481 << G4BestUnit(minMeanLife, "Time") << end
482 os << "Max life time (from G4DeexPrecoParame
483 << G4BestUnit(deex->GetMaxLifeTime(), "Ti
484 os << "Internal e- conversion flag
485 << deex->GetInternalConversionFlag() << e
486 os << "Stored internal conversion coefficien
487 << deex->StoreICLevelData() << endline;
488 os << "Enabled atomic relaxation mode
489 << applyARM << endline;
490 os << "Enable correlated gamma emission
491 << deex->CorrelatedGamma() << endline;
492 os << "Max 2J for sampling of angular correl
493 << deex->GetTwoJMAX() << endline;
494 os << "Atomic de-excitation enabled
495 << emparam->Fluo() << endline;
496 os << "Auger electron emission enabled
497 << emparam->Auger() << endline;
498 os << "Check EM cuts disabled for atomic de-
499 << emparam->DeexcitationIgnoreCut() << en
500 os << "Use Bearden atomic level energies
501 << emparam->BeardenFluoDir() << endline;
502 os << "Use ANSTO fluorescence model
503 << emparam->ANSTOFluoDir() << endline;
504 os << "Threshold for very long decay time at
505 << G4BestUnit(fThresholdForVeryLongDecayT
506 os << "=====================================
507 << G4endl;
508 os.precision(prec);
509 }
510
511
512 //////////////////////////////////////////////
513 //
514 // LoadDecayTable loads the decay scheme from
515 // for the parent nucleus.
516 //
517 //////////////////////////////////////////////
518
519 G4DecayTable* G4VRadioactiveDecay::LoadDecayTa
520 {
521 G4AutoLock lk(&radioactiveDecayMutex);
522 const G4String key = theIon->GetParticleName
523 auto dtptr = master_dkmap->find(key);
524 if (dtptr != master_dkmap->end()) {
525 lk.unlock();
526 return dtptr->second;
527 }
528
529 // Generate input data file name using Z and
530 // file containing radioactive decay data.
531 G4int A = theIon->GetAtomicMass();
532 G4int Z = theIon->GetAtomicNumber();
533
534 //G4cout << "LoadDecayTable for " << key <<
535
536 G4double levelEnergy = theIon->GetExcitation
537 G4Ions::G4FloatLevelBase floatingLevel = the
538
539 //Check if data have been provided by the us
540 G4String file;
541 G4int ke = 1000*A + Z;
542 auto ptr = theUserRDataFiles->find(ke);
543 if (ptr != theUserRDataFiles->end()) {
544 file = ptr->second;
545 } else {
546 std::ostringstream os;
547 os << dirPath << "/z" << Z << ".a" << A <<
548 file = os.str();
549 }
550
551 G4DecayTable* theDecayTable = new G4DecayTab
552 G4bool found(false); // True if energy l
553
554 std::ifstream DecaySchemeFile;
555 DecaySchemeFile.open(file);
556
557 if (DecaySchemeFile.good()) {
558 // Initialize variables used for reading i
559 G4bool floatMatch(false);
560 const G4int nMode = G4RadioactiveDecayMode
561 G4double modeTotalBR[nMode] = {0.0};
562 G4double modeSumBR[nMode];
563 for (G4int i = 0; i < nMode; ++i) {
564 modeSumBR[i] = 0.0;
565 }
566
567 char inputChars[120]={' '};
568 G4String inputLine;
569 G4String recordType("");
570 G4String floatingFlag("");
571 G4String daughterFloatFlag("");
572 G4Ions::G4FloatLevelBase daughterFloatLeve
573 G4RadioactiveDecayMode theDecayMode;
574 G4double decayModeTotal(0.0);
575 G4double parentExcitation(0.0);
576 G4double a(0.0);
577 G4double b(0.0);
578 G4double c(0.0);
579 G4double dummy(0.0);
580 G4BetaDecayType betaType(allowed);
581
582 // Loop through each data file record unti
583 // data relating to the nuclide of concern
584
585 G4bool complete(false); // bool insures o
586 // given parent e
587 G4int loop = 0;
588 /* Loop checking, 01.09.2015, D.Wright */
589 while (!complete && !DecaySchemeFile.getli
590 loop++;
591 if (loop > 100000) {
592 G4Exception("G4VRadioactiveDecay::Load
593 JustWarning, "While loop c
594 break;
595 }
596
597 inputLine = inputChars;
598 G4StrUtil::rstrip(inputLine);
599 if (inputChars[0] != '#' && inputLine.le
600 std::istringstream tmpStream(inputLine
601
602 if (inputChars[0] == 'P') {
603 // Nucleus is a parent type. Check
604 // matches that of theParentNucleus
605 tmpStream >> recordType >> parentExc
606 // "dummy" takes the place of half-l
607 // Now read in from ENSDFSTATE in p
608
609 if (found) {
610 complete = true;
611 } else {
612 // Take first level which matches
613 found = (std::abs(parentExcitation
614 if (floatingLevel != noFloat) {
615 // If floating level specificed,
616 floatMatch = (floatingLevel == G
617 if (!floatMatch) found = false;
618 }
619 }
620
621 } else if (found) {
622 // The right part of the radioactive
623 // through it to determine the mode
624
625 // Store for later the total decay p
626 if (inputLine.length() < 72) {
627 tmpStream >> theDecayMode >> dummy
628
629 switch (theDecayMode) {
630 case IT:
631 {
632 G4ITDecay* anITChannel = new G4ITDecay(t
633 theDecayTable->Insert(anITChannel);
634 }
635 break;
636 case BetaMinus:
637 modeTotalBR[BetaMinus] = decay
638 case BetaPlus:
639 modeTotalBR[BetaPlus] = decayM
640 case KshellEC:
641 modeTotalBR[KshellEC] = decayM
642 case LshellEC:
643 modeTotalBR[LshellEC] = decayM
644 case MshellEC:
645 modeTotalBR[MshellEC] = decayM
646 case NshellEC:
647 modeTotalBR[NshellEC] = decayM
648 case Alpha:
649 modeTotalBR[Alpha] = decayMode
650 case Proton:
651 modeTotalBR[Proton] = decayMod
652 case Neutron:
653 modeTotalBR[Neutron] = decayMo
654 case SpFission:
655 modeTotalBR[SpFission] = decay
656 case BDProton:
657 /* Not yet implemented */ bre
658 case BDNeutron:
659 /* Not yet implemented */ bre
660 case Beta2Minus:
661 /* Not yet implemented */ bre
662 case Beta2Plus:
663 /* Not yet implemented */ bre
664 case Proton2:
665 /* Not yet implemented */ bre
666 case Neutron2:
667 /* Not yet implemented */ bre
668 case Triton:
669 modeTotalBR[Triton] = decayMod
670 case RDM_ERROR:
671
672 default:
673 G4Exception("G4VRadioactiveDec
674 FatalException, "S
675 } // switch
676
677 } else {
678 if (inputLine.length() < 84) {
679 tmpStream >> theDecayMode >> a >
680 betaType = allowed;
681 } else {
682 tmpStream >> theDecayMode >> a >
683 }
684
685 // Allowed transitions are the def
686 // indicated in the last column.
687 a /= 1000.;
688 c /= 1000.;
689 b /= 100.;
690 daughterFloatLevel = G4Ions::Float
691
692 switch (theDecayMode) {
693 case BetaMinus:
694 {
695 G4BetaMinusDecay* aBetaMinusCh
696 new G4BetaMinusDecay(theIon,
697 daughte
698 //aBetaMinusChannel->DumpNuclearInfo();
699 theDecayTable->Insert(aBetaMin
700 modeSumBR[BetaMinus] += b;
701 }
702 break;
703
704 case BetaPlus:
705 {
706 G4BetaPlusDecay* aBetaPlusChan
707 new G4BetaPlusDecay(theIon,
708 daughter
709 //aBetaPlusChannel->DumpNuclearInfo();
710 theDecayTable->Insert(aBetaPlu
711 modeSumBR[BetaPlus] += b;
712 }
713 break;
714
715 case KshellEC: // K-shell elect
716 {
717 G4ECDecay* aKECChannel =
718 new G4ECDecay(theIon, b, c*M
719 daughterFloatL
720 //aKECChannel->DumpNuclearInfo();
721 aKECChannel->SetARM(applyARM);
722 theDecayTable->Insert(aKECChan
723 modeSumBR[KshellEC] += b;
724 }
725 break;
726
727 case LshellEC: // L-shell elect
728 {
729 G4ECDecay* aLECChannel =
730 new G4ECDecay(theIon, b, c*M
731 daughterFloatL
732 // aLECChannel->DumpNuclearInfo()
733 aLECChannel->SetARM(applyARM);
734 theDecayTable->Insert(aLECChan
735 modeSumBR[LshellEC] += b;
736 }
737 break;
738
739 case MshellEC: // M-shell elect
740 {
741 G4ECDecay* aMECChannel =
742 new G4ECDecay(theIon, b, c*M
743 daughterFloatL
744 // aMECChannel->DumpNuclearInfo()
745 aMECChannel->SetARM(applyARM);
746 theDecayTable->Insert(aMECChan
747 modeSumBR[MshellEC] += b;
748 }
749 break;
750
751 case NshellEC: // N-shell elect
752 {
753 G4ECDecay* aNECChannel =
754 new G4ECDecay(theIon, b, c*M
755 daughterFloatL
756 // aNECChannel->DumpNuclearInfo()
757 aNECChannel->SetARM(applyARM);
758 theDecayTable->Insert(aNECChan
759 modeSumBR[NshellEC] += b;
760 }
761 break;
762
763 case Alpha:
764 {
765 G4AlphaDecay* anAlphaChannel =
766 new G4AlphaDecay(theIon, b,
767 daughterFlo
768 // anAlphaChannel->DumpNuclearInf
769 theDecayTable->Insert(anAlphaC
770 modeSumBR[Alpha] += b;
771 }
772 break;
773
774 case Proton:
775 {
776 G4ProtonDecay* aProtonChannel
777 new G4ProtonDecay(theIon, b,
778 daughterFl
779 // aProtonChannel->DumpNuclearInf
780 theDecayTable->Insert(aProtonC
781 modeSumBR[Proton] += b;
782 }
783 break;
784
785 case Neutron:
786 {
787 G4NeutronDecay* aNeutronChanne
788 new G4NeutronDecay(theIon, b
789 daughterF
790 // aNeutronChannel->DumpNuclearIn
791 theDecayTable->Insert(aNeutron
792 modeSumBR[Neutron] += b;
793 }
794 break;
795
796 case SpFission:
797 {
798 G4SFDecay* aSpontFissChannel =
799 new G4SFDecay(theIon, b, c*M
800 theDecayTable->Insert(aSpontFi
801 modeSumBR[SpFission] += b;
802 }
803 break;
804
805 case BDProton:
806 // Not yet implemented
807 // G4cout << " beta-delayed
808 break;
809
810 case BDNeutron:
811 // Not yet implemented
812 // G4cout << " beta-delayed
813 break;
814
815 case Beta2Minus:
816 // Not yet implemented
817 // G4cout << " Double beta-
818 break;
819
820 case Beta2Plus:
821 // Not yet implemented
822 // G4cout << " Double beta+
823 break;
824
825 case Proton2:
826 // Not yet implemented
827 // G4cout << " Double proton
828 break;
829
830 case Neutron2:
831 // Not yet implemented
832 // G4cout << " Double beta-
833 break;
834
835 case Triton:
836 {
837 G4TritonDecay* aTritonChannel =
838 new G4TritonDecay(theIon,
839 theDecayTable->Insert(aTritonChannel);
840 modeSumBR[Triton] += b;
841 }
842 break;
843
844 case RDM_ERROR:
845
846 default:
847 G4Exception("G4VRadioactiveDec
848 FatalException, "S
849 } // switch
850 } // line < 72
851 } // if char == P
852 } // if char != #
853 } // While
854
855 // Go through the decay table and make sur
856 // correctly normalised.
857
858 G4VDecayChannel* theChannel = 0;
859 G4NuclearDecay* theNuclearDecayChannel = 0
860 G4String mode = "";
861
862 G4double theBR = 0.0;
863 for (G4int i = 0; i < theDecayTable->entri
864 theChannel = theDecayTable->GetDecayChan
865 theNuclearDecayChannel = static_cast<G4N
866 theDecayMode = theNuclearDecayChannel->G
867
868 if (theDecayMode != IT) {
869 theBR = theChannel->GetBR();
870 theChannel->SetBR(theBR*modeTotalBR[theDecay
871 }
872 }
873 } // decay file exists
874
875 DecaySchemeFile.close();
876
877 if (!found && levelEnergy > 0) {
878 // Case where IT cascade for excited isoto
879 // Decay mode is isomeric transition.
880 G4ITDecay* anITChannel = new G4ITDecay(the
881 theDecayTable->Insert(anITChannel);
882 }
883
884 if (GetVerboseLevel() > 1) {
885 theDecayTable->DumpInfo();
886 }
887
888 // store in master library
889 (*master_dkmap)[theIon->GetParticleName()] =
890 lk.unlock();
891 return theDecayTable;
892 }
893
894 void G4VRadioactiveDecay::AddUserDecayDataFile
895
896 {
897 if (Z < 1 || A < 2) G4cout << "Z and A not v
898
899 std::ifstream DecaySchemeFile(filename);
900 if (DecaySchemeFile) {
901 G4int ID_ion = A*1000 + Z;
902 (*theUserRDataFiles)[ID_ion] = filename;
903 } else {
904 G4ExceptionDescription ed;
905 ed << filename << " does not exist! " << G
906 G4Exception("G4VRadioactiveDecay::AddUserD
907 FatalException, ed);
908 }
909 }
910
911 //////////////////////////////////////////////
912 //
913 // DecayIt
914 //
915 //////////////////////////////////////////////
916
917 G4VParticleChange*
918 G4VRadioactiveDecay::DecayIt(const G4Track& th
919 {
920 // Initialize G4ParticleChange object, get p
921 fParticleChangeForRadDecay.Initialize(theTra
922 fParticleChangeForRadDecay.ProposeWeight(the
923 const G4DynamicParticle* theParticle = theTr
924 const G4ParticleDefinition* theParticleDef =
925
926 // First check whether RDM applies to the cu
927 if (!isAllVolumesMode) {
928 if (!std::binary_search(ValidVolumes.begin
929 theTrack.GetVolume()->Get
930 #ifdef G4VERBOSE
931 if (GetVerboseLevel()>1) {
932 G4cout <<"G4VRadioactiveDecay::DecayIt
933 << theTrack.GetVolume()->GetLog
934 << " is not selected for the RD
935 G4cout << " There are " << ValidVolume
936 G4cout << " The Valid volumes are: ";
937 for (auto const & vol : ValidVolumes) { G4co
938 G4cout << G4endl;
939 }
940 #endif
941 fParticleChangeForRadDecay.SetNumberOfSe
942
943 // Kill the parent particle.
944 fParticleChangeForRadDecay.ProposeTrackS
945 fParticleChangeForRadDecay.ProposeLocalE
946 ClearNumberOfInteractionLengthLeft();
947 return &fParticleChangeForRadDecay;
948 }
949 }
950
951 // Now check if particle is valid for RDM
952 G4DecayTable* theDecayTable = GetDecayTable(
953 if ( theDecayTable == nullptr || theDecayTab
954 // Particle is not an ion or is outside th
955 #ifdef G4VERBOSE
956 if (GetVerboseLevel() > 1) {
957 G4cout << "G4VRadioactiveDecay::DecayIt
958 << theParticleDef->GetParticleNam
959 << " is outside (Z,A) limits set
960 << G4endl;
961 }
962 #endif
963 fParticleChangeForRadDecay.SetNumberOfSeco
964
965 // Kill the parent particle
966 fParticleChangeForRadDecay.ProposeTrackSta
967 fParticleChangeForRadDecay.ProposeLocalEne
968 ClearNumberOfInteractionLengthLeft();
969 return &fParticleChangeForRadDecay;
970 }
971 //G4cout << "DecayIt for " << theParticleDef
972 // << " isAllVolumesMode:" << isAllVolumes
973 // << " decayTable=" << theDecayTable << G
974
975 // Data found. Decay nucleus without varianc
976 DecayAnalog(theTrack, theDecayTable);
977 return &fParticleChangeForRadDecay;
978 }
979
980 void G4VRadioactiveDecay::DecayAnalog(const G4
981 G4DecayT
982 {
983 const G4DynamicParticle* theParticle = theTr
984 const G4ParticleDefinition* theParticleDef =
985 //G4cout << "DecayIt for " << theParticleDef
986 G4DecayProducts* products = DoDecay(*thePart
987
988 // Check if the product is the same as input
989 // necessary to prevent infinite loop (11/05
990 if (nullptr == products || products->entries
991 fParticleChangeForRadDecay.SetNumberOfSeco
992 fParticleChangeForRadDecay.ProposeTrackSta
993 fParticleChangeForRadDecay.ProposeLocalEne
994 ClearNumberOfInteractionLengthLeft();
995 delete products;
996 return;
997 }
998
999 G4double energyDeposit = 0.0;
1000 G4double finalGlobalTime = theTrack.GetGlob
1001 G4double finalLocalTime = theTrack.GetLocal
1002
1003 // Get parent particle information and boos
1004 // laboratory frame
1005
1006 // ParentEnergy used for the boost should b
1007 // of the parent ion without the energy of
1008 // (correction for bug 1359 by L. Desorgher
1009 G4double ParentEnergy = theParticle->GetKin
1010 + theParticle->GetPar
1011 G4ThreeVector ParentDirection(theParticle->
1012
1013 if (theTrack.GetTrackStatus() == fStopButAl
1014 // this condition seems to be always True
1015
1016 // The particle is decayed at rest
1017 // Since the time is for the particle at
1018 // lapsed between particle coming to rest
1019 // is sampled with the mean-life of the p
1020 // PDGTime < 0. (F.Lei 11/05/10)
1021 G4double temptime = -std::log(G4UniformRa
1022 theParticleDef->GetPD
1023 if (temptime < 0.) temptime = 0.;
1024 finalGlobalTime += temptime;
1025 finalLocalTime += temptime;
1026 energyDeposit += theParticle->GetKineticE
1027
1028 // Kill the parent particle, and ignore i
1029 // threshold fThresholdForVeryLongDecayTi
1030 // to more than twice the age of the univ
1031 // This kind of cut has been introduced (
1032 // account energy depositions happening a
1033 // ordinary materials used in calorimetry
1034 // (via their natural unstable, but very
1035 // W183, W180 and Pb204).
1036 // Note that the cut is not on the averag
1037 // sampled global decay time.
1038 if ( finalGlobalTime > fThresholdForVeryL
1039 fParticleChangeForRadDecay.SetNumberOfS
1040 fParticleChangeForRadDecay.ProposeTrack
1041 fParticleChangeForRadDecay.ProposeLocal
1042 ClearNumberOfInteractionLengthLeft();
1043 delete products;
1044 return;
1045 }
1046 }
1047 products->Boost(ParentEnergy, ParentDirecti
1048
1049 // Add products in theParticleChangeForRadD
1050 G4int numberOfSecondaries = products->entri
1051 fParticleChangeForRadDecay.SetNumberOfSecon
1052
1053 if (GetVerboseLevel() > 1) {
1054 G4cout << "G4VRadioactiveDecay::DecayAnal
1055 G4cout << " Time: " << finalGlobalTime/ns
1056 G4cout << " X:" << (theTrack.GetPosition(
1057 G4cout << " Y:" << (theTrack.GetPosition(
1058 G4cout << " Z:" << (theTrack.GetPosition(
1059 G4cout << G4endl;
1060 G4cout << "G4Decay::DecayIt : decay produ
1061 products->DumpInfo();
1062 products->IsChecked();
1063 }
1064
1065 const G4int modelID_forIT = G4PhysicsModelC
1066 G4int modelID = modelID_forIT + 10*theRadDe
1067 const G4int modelID_forAtomicRelaxation =
1068 G4PhysicsModelCatalog::GetModelID( "model
1069 for ( G4int index = 0; index < numberOfSeco
1070 G4Track* secondary = new G4Track( product
1071 theTrac
1072 secondary->SetWeight( theTrack.GetWeight(
1073 secondary->SetCreatorModelID( modelID );
1074 // Change for atomics relaxation
1075 if ( theRadDecayMode == IT && index > 0
1076 if ( index == numberOfSecondaries-1 ) {
1077 secondary->SetCreatorModelID( modelID_forIT
1078 } else {
1079 secondary->SetCreatorModelID( modelID_forAt
1080 }
1081 } else if ( theRadDecayMode >= KshellEC
1082 index < numberOfSecondaries-1 ) {
1083 secondary->SetCreatorModelID( modelID_f
1084 }
1085 secondary->SetGoodForTrackingFlag();
1086 secondary->SetTouchableHandle( theTrack.G
1087 fParticleChangeForRadDecay.AddSecondary(
1088 }
1089
1090 delete products;
1091
1092 // Kill the parent particle
1093 fParticleChangeForRadDecay.ProposeTrackStat
1094 fParticleChangeForRadDecay.ProposeLocalEner
1095 fParticleChangeForRadDecay.ProposeLocalTime
1096
1097 // Reset NumberOfInteractionLengthLeft.
1098 ClearNumberOfInteractionLengthLeft();
1099 }
1100
1101
1102 G4DecayProducts*
1103 G4VRadioactiveDecay::DoDecay(const G4Particle
1104 G4DecayTable* th
1105 {
1106 G4DecayProducts* products = nullptr;
1107
1108 // Choose a decay channel.
1109 // G4DecayTable::SelectADecayChannel checks
1110 // exceeds parent mass. Pass it the parent
1111 // for difference in mass defect.
1112 G4double parentPlusQ = theParticleDef.GetPD
1113 G4VDecayChannel* theDecayChannel = theDecay
1114
1115 if (theDecayChannel == nullptr) {
1116 // Decay channel not found.
1117 G4ExceptionDescription ed;
1118 ed << " Cannot determine decay channel fo
1119 G4Exception("G4VRadioactiveDecay::DoDecay
1120 FatalException, ed);
1121 } else {
1122 // A decay channel has been identified, s
1123 #ifdef G4VERBOSE
1124 if (GetVerboseLevel() > 1) {
1125 G4cout << "G4VRadioactiveDecay::DoIt :
1126 << theDecayChannel << G4endl;
1127 }
1128 #endif
1129 theRadDecayMode = static_cast<G4NuclearDe
1130
1131 // for IT decay use local G4ITDecay class
1132 if (theRadDecayMode == IT) {
1133 decayIT->SetupDecay(&theParticleDef);
1134 products = decayIT->DecayIt(0.0);
1135 } else {
1136 // for others decayes use shared class
1137 products = theDecayChannel->DecayIt(the
1138 }
1139
1140 // Apply directional bias if requested by
1141 CollimateDecay(products);
1142 }
1143
1144 return products;
1145 }
1146
1147
1148 // Apply directional bias for "visible" daugh
1149 void G4VRadioactiveDecay::CollimateDecay(G4De
1150
1151 if (origin == forceDecayDirection) return;
1152 if (180.*deg == forceDecayHalfAngle) return
1153 if (0 == products || 0 == products->entries
1154
1155 #ifdef G4VERBOSE
1156 if (GetVerboseLevel() > 1) G4cout << "Begin
1157 #endif
1158
1159 // Particles suitable for directional biasi
1160 static const G4ParticleDefinition* electron
1161 static const G4ParticleDefinition* positron
1162 static const G4ParticleDefinition* neutron
1163 static const G4ParticleDefinition* gamma
1164 static const G4ParticleDefinition* alpha
1165 static const G4ParticleDefinition* triton
1166 static const G4ParticleDefinition* proton
1167
1168 G4ThreeVector newDirection; // Re-use to
1169 for (G4int i=0; i<products->entries(); ++i)
1170 G4DynamicParticle* daughter = (*products)
1171 const G4ParticleDefinition* daughterType
1172 daughter->G
1173 if (daughterType == electron || daughterT
1174 daughterType == neutron || daughterType ==
1175 daughterType == alpha || daughterType == tr
1176 CollimateDecayProduct(daughter);
1177 }
1178 }
1179 }
1180
1181 void G4VRadioactiveDecay::CollimateDecayProdu
1182 #ifdef G4VERBOSE
1183 if (GetVerboseLevel() > 1) {
1184 G4cout << "CollimateDecayProduct for daug
1185 << daughter->GetParticleDefinition()->Ge
1186 }
1187 #endif
1188
1189 G4ThreeVector collimate = ChooseCollimation
1190 if (origin != collimate) daughter->SetMomen
1191 }
1192
1193
1194 // Choose random direction within collimation
1195
1196 G4ThreeVector G4VRadioactiveDecay::ChooseColl
1197 if (origin == forceDecayDirection) return o
1198 if (forceDecayHalfAngle == 180.*deg) return
1199
1200 G4ThreeVector dir = forceDecayDirection;
1201
1202 // Return direction offset by random throw
1203 if (forceDecayHalfAngle > 0.) {
1204 // Generate uniform direction around cent
1205 G4double phi = 2.*pi*G4UniformRand();
1206 G4double cosMin = std::cos(forceDecayHalf
1207 G4double cosTheta = (1.-cosMin)*G4Uniform
1208
1209 dir.setPhi(dir.phi()+phi);
1210 dir.setTheta(dir.theta()+std::acos(cosThe
1211 }
1212
1213 #ifdef G4VERBOSE
1214 if (GetVerboseLevel()>1)
1215 G4cout << " ChooseCollimationDirection re
1216 #endif
1217
1218 return dir;
1219 }
1220
1221