| 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 // 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. Truscott.
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::levelTolerance = 10.0*CLHEP::eV;
100 const G4ThreeVector G4VRadioactiveDecay::origin(0.,0.,0.);
101
102 DecayTableMap* G4VRadioactiveDecay::master_dkmap = nullptr;
103 std::map<G4int, G4String>* G4VRadioactiveDecay::theUserRDataFiles = nullptr;
104 G4String G4VRadioactiveDecay::dirPath = "";
105
106 namespace
107 {
108 G4Mutex radioactiveDecayMutex = G4MUTEX_INITIALIZER;
109 }
110
111 G4VRadioactiveDecay::G4VRadioactiveDecay(const G4String& processName,
112 const G4double timeThreshold)
113 : G4VRestDiscreteProcess(processName, fDecay),
114 fThresholdForVeryLongDecayTime( 1.0*CLHEP::year )
115 {
116 if (GetVerboseLevel() > 1) {
117 G4cout << "G4VRadioactiveDecay constructor: processName = " << processName
118 << G4endl;
119 }
120
121 SetProcessSubType(fRadioactiveDecay);
122
123 theRadioactiveDecayMessenger = new G4RadioactiveDecayMessenger(this);
124 pParticleChange = &fParticleChangeForRadDecay;
125
126 // Check data directory
127 if (dirPath.empty()) {
128 const char* path_var = G4FindDataDir("G4RADIOACTIVEDATA");
129 if (nullptr == path_var) {
130 G4Exception("G4VRadioactiveDecay()", "HAD_RDM_200", FatalException,
131 "Environment variable G4RADIOACTIVEDATA is not set");
132 } else {
133 dirPath = path_var; // convert to string
134 std::ostringstream os;
135 os << dirPath << "/z1.a3"; // used as a dummy
136 std::ifstream testFile;
137 testFile.open(os.str() );
138 if ( !testFile.is_open() )
139 G4Exception("G4VRadioactiveDecay()","HAD_RDM_201",FatalException,
140 "Environment variable G4RADIOACTIVEDATA is set, but does not point to correct directory");
141 }
142 }
143 // Set up photon evaporation for use in G4ITDecay
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, G4String>;
155 }
156
157 // RDM applies to all logical volumes by default
158 SelectAllVolumes();
159 G4HadronicProcessStore::Instance()->RegisterExtraProcess(this);
160
161 // The time threshold for radioactive decays can be set in 3 ways:
162 // 1. Via C++ interface: G4HadronicParameters::Instance()->SetTimeThresholdForRadioactiveDecay(value)
163 // 2. Via the second parameter of the G4VRadioactiveDecay constructor
164 // 3. Via UI command: /process/had/rdm/thresholdForVeryLongDecayTime value
165 // If both 1. and 2. are specified (at the moment when the G4VRadioactiveDecay constructor is called),
166 // then we take the larger value, to be conservative.
167 // If, later on (after invoking the G4VRadioactiveDecay constructor) 3. is specified,
168 // then this value is used (and the eventual values 1. and/or 2. are ignored).
169 G4double timeThresholdBis = G4HadronicParameters::Instance()->GetTimeThresholdForRadioactiveDecay();
170 if ( timeThreshold > 0.0 || timeThresholdBis > 0.0 ) {
171 if ( timeThreshold > timeThresholdBis ) timeThresholdBis = timeThreshold;
172 fThresholdForVeryLongDecayTime = timeThresholdBis;
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 G4ParticleDefinition& aParticle)
199 {
200 const G4String& pname = aParticle.GetParticleName();
201 if (pname == "GenericIon" || pname == "triton") { return true; }
202 // All particles other than G4Ions, are rejected by default
203 const G4Ions* p = dynamic_cast<const G4Ions*>(&aParticle);
204 if (nullptr == p) { return false; }
205
206 // excited isomere may decay via gamma evaporation
207 if (p->GetExcitationEnergy() > 0.0) { return true; }
208
209 // Check on life time
210 G4double lifeTime = p->GetPDGLifeTime();
211 if (lifeTime < 0.0 || lifeTime > fThresholdForVeryLongDecayTime) {
212 return false;
213 }
214
215 // Determine whether the nuclide falls into the correct A and Z range
216 G4int A = p->GetAtomicMass();
217 G4int Z = p->GetAtomicNumber();
218
219 if (A > theNucleusLimits.GetAMax() || A < theNucleusLimits.GetAMin() ||
220 Z > theNucleusLimits.GetZMax() || Z < theNucleusLimits.GetZMin()) {
221 return false;
222 }
223
224 return true;
225 }
226
227
228 G4VParticleChange* G4VRadioactiveDecay::AtRestDoIt(const G4Track& theTrack,
229 const G4Step& theStep)
230 {
231 return DecayIt(theTrack, theStep);
232 }
233
234
235 G4VParticleChange* G4VRadioactiveDecay::PostStepDoIt(const G4Track& theTrack,
236 const G4Step& theStep)
237 {
238 return DecayIt(theTrack, theStep);
239 }
240
241
242 void G4VRadioactiveDecay::ProcessDescription(std::ostream& outFile) const
243 {
244 outFile << "The radioactive decay process (G4VRadioactiveDecay) handles the\n"
245 << "alpha, beta+, beta-, electron capture and isomeric transition\n"
246 << "decays of nuclei (G4GenericIon) with masses A > 4.\n"
247 << "The required half-lives and decay schemes are retrieved from\n"
248 << "the RadioactiveDecay database which was derived from ENSDF.\n";
249 }
250
251
252
253
254 G4DecayTable* G4VRadioactiveDecay::GetDecayTable(const G4ParticleDefinition* aNucleus)
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 G4Ions*>(aNucleus);
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 G4String& aVolume)
274 {
275 G4LogicalVolumeStore* theLogicalVolumes = G4LogicalVolumeStore::GetInstance();
276 G4LogicalVolume* volume = nullptr;
277 volume = theLogicalVolumes->GetVolume(aVolume);
278 if (volume != nullptr)
279 {
280 ValidVolumes.push_back(aVolume);
281 std::sort(ValidVolumes.begin(), ValidVolumes.end());
282 // sort need for performing binary_search
283
284 if (GetVerboseLevel() > 0)
285 G4cout << " Radioactive decay applied to " << aVolume << G4endl;
286 }
287 else
288 {
289 G4ExceptionDescription ed;
290 ed << aVolume << " is not a valid logical volume name."
291 << " Decay not activated for it."
292 << G4endl;
293 G4Exception("G4VRadioactiveDecay::SelectAVolume()", "HAD_RDM_300",
294 JustWarning, ed);
295 }
296 }
297
298
299 void G4VRadioactiveDecay::DeselectAVolume(const G4String& aVolume)
300 {
301 G4LogicalVolumeStore* theLogicalVolumes = G4LogicalVolumeStore::GetInstance();
302 G4LogicalVolume* volume = nullptr;
303 volume = theLogicalVolumes->GetVolume(aVolume);
304 if (volume != nullptr)
305 {
306 auto location= std::find(ValidVolumes.cbegin(),ValidVolumes.cend(),aVolume);
307 if (location != ValidVolumes.cend() )
308 {
309 ValidVolumes.erase(location);
310 std::sort(ValidVolumes.begin(), ValidVolumes.end());
311 isAllVolumesMode = false;
312 if (GetVerboseLevel() > 0)
313 G4cout << " G4VRadioactiveDecay::DeselectAVolume: " << aVolume
314 << " is removed from list " << G4endl;
315 }
316 else
317 {
318 G4ExceptionDescription ed;
319 ed << aVolume << " is not in the list. No action taken." << G4endl;
320 G4Exception("G4VRadioactiveDecay::DeselectAVolume()", "HAD_RDM_300",
321 JustWarning, ed);
322 }
323 }
324 else
325 {
326 G4ExceptionDescription ed;
327 ed << aVolume << " is not a valid logical volume name. No action taken."
328 << G4endl;
329 G4Exception("G4VRadioactiveDecay::DeselectAVolume()", "HAD_RDM_300",
330 JustWarning, ed);
331 }
332 }
333
334
335 void G4VRadioactiveDecay::SelectAllVolumes()
336 {
337 G4LogicalVolumeStore* theLogicalVolumes = G4LogicalVolumeStore::GetInstance();
338 G4LogicalVolume* volume = nullptr;
339 ValidVolumes.clear();
340 #ifdef G4VERBOSE
341 if (GetVerboseLevel()>1)
342 G4cout << " RDM Applies to all Volumes" << G4endl;
343 #endif
344 for (std::size_t i = 0; i < theLogicalVolumes->size(); ++i){
345 volume = (*theLogicalVolumes)[i];
346 ValidVolumes.push_back(volume->GetName());
347 #ifdef G4VERBOSE
348 if (GetVerboseLevel()>1)
349 G4cout << " RDM Applies to Volume " << volume->GetName() << G4endl;
350 #endif
351 }
352 std::sort(ValidVolumes.begin(), ValidVolumes.end());
353 // sort needed in order to allow binary_search
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 removed from all volumes" << G4endl;
364 #endif
365 }
366
367
368 // GetMeanLifeTime (required by the base class)
369 G4double G4VRadioactiveDecay::GetMeanLifeTime(const G4Track& theTrack,
370 G4ForceCondition*)
371 {
372 G4double meanlife = DBL_MAX;
373 const G4ParticleDefinition* theParticleDef = theTrack.GetParticleDefinition();
374 if (!IsApplicable(*theParticleDef)) { return meanlife; }
375 G4double theLife = theParticleDef->GetPDGLifeTime();
376 #ifdef G4VERBOSE
377 if (GetVerboseLevel() > 2) {
378 G4cout << "G4VRadioactiveDecay::GetMeanLifeTime() for "
379 << theParticleDef->GetParticleName() << G4endl;
380 G4cout << "KineticEnergy(GeV)=" << theTrack.GetKineticEnergy()/CLHEP::GeV
381 << " Mass(GeV)=" << theParticleDef->GetPDGMass()/CLHEP::GeV
382 << " LifeTime(ns)=" << theLife/CLHEP::ns << G4endl;
383 }
384 #endif
385 if (theLife >= 0.0 && theLife <= fThresholdForVeryLongDecayTime) {
386 meanlife = theLife;
387 }
388
389 if (meanlife == DBL_MAX) {
390 const G4Ions* ion = dynamic_cast<const G4Ions*>(theParticleDef);
391 if (nullptr != ion && ion->GetExcitationEnergy() > 0.0) {
392 meanlife = 0.0;
393 }
394 }
395
396 #ifdef G4VERBOSE
397 if (GetVerboseLevel() > 2)
398 G4cout << "G4VRadioactiveDecay::GetMeanLifeTime: "
399 << meanlife/CLHEP::s << " second " << G4endl;
400 #endif
401
402 return meanlife;
403 }
404
405
406 ////////////////////////////////////////////////////////////////////////////////
407 // //
408 // GetMeanFreePath for decay in flight //
409 // //
410 ////////////////////////////////////////////////////////////////////////////////
411
412 G4double G4VRadioactiveDecay::GetMeanFreePath(const G4Track& aTrack, G4double,
413 G4ForceCondition* fc)
414 {
415 G4double res = DBL_MAX;
416 G4double lifeTime = GetMeanLifeTime(aTrack, fc);
417 if (lifeTime > 0.0 && lifeTime < DBL_MAX) {
418 auto dParticle = aTrack.GetDynamicParticle();
419 res = lifeTime*dParticle->GetTotalEnergy()*aTrack.GetVelocity()/dParticle->GetMass();
420 } else {
421 res = lifeTime;
422 }
423
424 #ifdef G4VERBOSE
425 if (GetVerboseLevel() > 2) {
426 G4cout << "G4VRadioactiveDecay::GetMeanFreePath() for "
427 << aTrack.GetDefinition()->GetParticleName() << G4endl;
428 G4cout << " kinEnergy(GeV)=" << aTrack.GetKineticEnergy()/CLHEP::GeV
429 << " lifeTime(ns)=" << lifeTime
430 << " mean free path(cm)=" << res/CLHEP::cm << G4endl;
431 }
432 #endif
433 return res;
434 }
435
436 ////////////////////////////////////////////////////////////////////////////////
437 // //
438 // BuildPhysicsTable - initialization of atomic de-excitation //
439 // //
440 ////////////////////////////////////////////////////////////////////////////////
441
442 void G4VRadioactiveDecay::BuildPhysicsTable(const G4ParticleDefinition& p)
443 {
444 if (isInitialised) { return; }
445 isInitialised = true;
446 if (G4HadronicParameters::Instance()->GetVerboseLevel() > 0 &&
447 G4Threading::IsMasterThread() && "GenericIon" == p.GetParticleName()) {
448 StreamInfo(G4cout, "\n");
449 }
450 photonEvaporation->Initialise();
451 photonEvaporation->RDMForced(true);
452 photonEvaporation->SetICM(true);
453 decayIT->SetARM(applyARM);
454
455 G4HadronicProcessStore::Instance()->RegisterParticleForExtraProcess(this, &p);
456 G4HadronicProcessStore::Instance()->PrintInfo(&p);
457 }
458
459 ////////////////////////////////////////////////////////////////////////////////
460 // //
461 // StreamInfo - stream out parameters //
462 // //
463 ////////////////////////////////////////////////////////////////////////////////
464
465 void
466 G4VRadioactiveDecay::StreamInfo(std::ostream& os, const G4String& endline)
467 {
468 G4DeexPrecoParameters* deex =
469 G4NuclearLevelData::GetInstance()->GetParameters();
470 G4EmParameters* emparam = G4EmParameters::Instance();
471 G4double minMeanLife = G4NuclideTable::GetInstance()->GetMeanLifeThreshold();
472
473 G4long prec = os.precision(5);
474 os << "======================================================================"
475 << endline;
476 os << "====== Radioactive Decay Physics Parameters ======="
477 << endline;
478 os << "======================================================================"
479 << endline;
480 os << "min MeanLife (from G4NuclideTable) "
481 << G4BestUnit(minMeanLife, "Time") << endline;
482 os << "Max life time (from G4DeexPrecoParameters) "
483 << G4BestUnit(deex->GetMaxLifeTime(), "Time") << endline;
484 os << "Internal e- conversion flag "
485 << deex->GetInternalConversionFlag() << endline;
486 os << "Stored internal conversion coefficients "
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 correlations "
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-excitation "
499 << emparam->DeexcitationIgnoreCut() << endline;
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 rest "
505 << G4BestUnit(fThresholdForVeryLongDecayTime, "Time") << endline;
506 os << "======================================================================"
507 << G4endl;
508 os.precision(prec);
509 }
510
511
512 ////////////////////////////////////////////////////////////////////////////////
513 // //
514 // LoadDecayTable loads the decay scheme from the RadioactiveDecay database //
515 // for the parent nucleus. //
516 // //
517 ////////////////////////////////////////////////////////////////////////////////
518
519 G4DecayTable* G4VRadioactiveDecay::LoadDecayTable(const G4Ions* theIon)
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 A of the parent nucleus
530 // file containing radioactive decay data.
531 G4int A = theIon->GetAtomicMass();
532 G4int Z = theIon->GetAtomicNumber();
533
534 //G4cout << "LoadDecayTable for " << key << " Z=" << Z << " A=" << A << G4endl;
535
536 G4double levelEnergy = theIon->GetExcitationEnergy();
537 G4Ions::G4FloatLevelBase floatingLevel = theIon->GetFloatLevelBase();
538
539 //Check if data have been provided by the user
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 << '\0';
548 file = os.str();
549 }
550
551 G4DecayTable* theDecayTable = new G4DecayTable();
552 G4bool found(false); // True if energy level matches one in table
553
554 std::ifstream DecaySchemeFile;
555 DecaySchemeFile.open(file);
556
557 if (DecaySchemeFile.good()) {
558 // Initialize variables used for reading in radioactive decay data
559 G4bool floatMatch(false);
560 const G4int nMode = G4RadioactiveDecayModeSize;
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 daughterFloatLevel;
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 until you identify the decay
583 // data relating to the nuclide of concern.
584
585 G4bool complete(false); // bool insures only one set of values read for any
586 // given parent energy level
587 G4int loop = 0;
588 /* Loop checking, 01.09.2015, D.Wright */
589 while (!complete && !DecaySchemeFile.getline(inputChars, 120).eof()) {
590 loop++;
591 if (loop > 100000) {
592 G4Exception("G4VRadioactiveDecay::LoadDecayTable()", "HAD_RDM_100",
593 JustWarning, "While loop count exceeded");
594 break;
595 }
596
597 inputLine = inputChars;
598 G4StrUtil::rstrip(inputLine);
599 if (inputChars[0] != '#' && inputLine.length() != 0) {
600 std::istringstream tmpStream(inputLine);
601
602 if (inputChars[0] == 'P') {
603 // Nucleus is a parent type. Check excitation level to see if it
604 // matches that of theParentNucleus
605 tmpStream >> recordType >> parentExcitation >> floatingFlag >> dummy;
606 // "dummy" takes the place of half-life
607 // Now read in from ENSDFSTATE in particle category
608
609 if (found) {
610 complete = true;
611 } else {
612 // Take first level which matches excitation energy regardless of floating level
613 found = (std::abs(parentExcitation*keV - levelEnergy) < levelTolerance);
614 if (floatingLevel != noFloat) {
615 // If floating level specificed, require match of both energy and floating level
616 floatMatch = (floatingLevel == G4Ions::FloatLevelBase(floatingFlag.back()) );
617 if (!floatMatch) found = false;
618 }
619 }
620
621 } else if (found) {
622 // The right part of the radioactive decay data file has been found. Search
623 // through it to determine the mode of decay of the subsequent records.
624
625 // Store for later the total decay probability for each decay mode
626 if (inputLine.length() < 72) {
627 tmpStream >> theDecayMode >> dummy >> decayModeTotal;
628
629 switch (theDecayMode) {
630 case IT:
631 {
632 G4ITDecay* anITChannel = new G4ITDecay(theIon, decayModeTotal, 0.0, 0.0);
633 theDecayTable->Insert(anITChannel);
634 }
635 break;
636 case BetaMinus:
637 modeTotalBR[BetaMinus] = decayModeTotal; break;
638 case BetaPlus:
639 modeTotalBR[BetaPlus] = decayModeTotal; break;
640 case KshellEC:
641 modeTotalBR[KshellEC] = decayModeTotal; break;
642 case LshellEC:
643 modeTotalBR[LshellEC] = decayModeTotal; break;
644 case MshellEC:
645 modeTotalBR[MshellEC] = decayModeTotal; break;
646 case NshellEC:
647 modeTotalBR[NshellEC] = decayModeTotal; break;
648 case Alpha:
649 modeTotalBR[Alpha] = decayModeTotal; break;
650 case Proton:
651 modeTotalBR[Proton] = decayModeTotal; break;
652 case Neutron:
653 modeTotalBR[Neutron] = decayModeTotal; break;
654 case SpFission:
655 modeTotalBR[SpFission] = decayModeTotal; break;
656 case BDProton:
657 /* Not yet implemented */ break;
658 case BDNeutron:
659 /* Not yet implemented */ break;
660 case Beta2Minus:
661 /* Not yet implemented */ break;
662 case Beta2Plus:
663 /* Not yet implemented */ break;
664 case Proton2:
665 /* Not yet implemented */ break;
666 case Neutron2:
667 /* Not yet implemented */ break;
668 case Triton:
669 modeTotalBR[Triton] = decayModeTotal; break;
670 case RDM_ERROR:
671
672 default:
673 G4Exception("G4VRadioactiveDecay::LoadDecayTable()", "HAD_RDM_000",
674 FatalException, "Selected decay mode does not exist");
675 } // switch
676
677 } else {
678 if (inputLine.length() < 84) {
679 tmpStream >> theDecayMode >> a >> daughterFloatFlag >> b >> c;
680 betaType = allowed;
681 } else {
682 tmpStream >> theDecayMode >> a >> daughterFloatFlag >> b >> c >> betaType;
683 }
684
685 // Allowed transitions are the default. Forbidden transitions are
686 // indicated in the last column.
687 a /= 1000.;
688 c /= 1000.;
689 b /= 100.;
690 daughterFloatLevel = G4Ions::FloatLevelBase(daughterFloatFlag.back());
691
692 switch (theDecayMode) {
693 case BetaMinus:
694 {
695 G4BetaMinusDecay* aBetaMinusChannel =
696 new G4BetaMinusDecay(theIon, b, c*MeV, a*MeV,
697 daughterFloatLevel, betaType);
698 //aBetaMinusChannel->DumpNuclearInfo();
699 theDecayTable->Insert(aBetaMinusChannel);
700 modeSumBR[BetaMinus] += b;
701 }
702 break;
703
704 case BetaPlus:
705 {
706 G4BetaPlusDecay* aBetaPlusChannel =
707 new G4BetaPlusDecay(theIon, b, c*MeV, a*MeV,
708 daughterFloatLevel, betaType);
709 //aBetaPlusChannel->DumpNuclearInfo();
710 theDecayTable->Insert(aBetaPlusChannel);
711 modeSumBR[BetaPlus] += b;
712 }
713 break;
714
715 case KshellEC: // K-shell electron capture
716 {
717 G4ECDecay* aKECChannel =
718 new G4ECDecay(theIon, b, c*MeV, a*MeV,
719 daughterFloatLevel, KshellEC);
720 //aKECChannel->DumpNuclearInfo();
721 aKECChannel->SetARM(applyARM);
722 theDecayTable->Insert(aKECChannel);
723 modeSumBR[KshellEC] += b;
724 }
725 break;
726
727 case LshellEC: // L-shell electron capture
728 {
729 G4ECDecay* aLECChannel =
730 new G4ECDecay(theIon, b, c*MeV, a*MeV,
731 daughterFloatLevel, LshellEC);
732 // aLECChannel->DumpNuclearInfo();
733 aLECChannel->SetARM(applyARM);
734 theDecayTable->Insert(aLECChannel);
735 modeSumBR[LshellEC] += b;
736 }
737 break;
738
739 case MshellEC: // M-shell electron capture
740 {
741 G4ECDecay* aMECChannel =
742 new G4ECDecay(theIon, b, c*MeV, a*MeV,
743 daughterFloatLevel, MshellEC);
744 // aMECChannel->DumpNuclearInfo();
745 aMECChannel->SetARM(applyARM);
746 theDecayTable->Insert(aMECChannel);
747 modeSumBR[MshellEC] += b;
748 }
749 break;
750
751 case NshellEC: // N-shell electron capture
752 {
753 G4ECDecay* aNECChannel =
754 new G4ECDecay(theIon, b, c*MeV, a*MeV,
755 daughterFloatLevel, NshellEC);
756 // aNECChannel->DumpNuclearInfo();
757 aNECChannel->SetARM(applyARM);
758 theDecayTable->Insert(aNECChannel);
759 modeSumBR[NshellEC] += b;
760 }
761 break;
762
763 case Alpha:
764 {
765 G4AlphaDecay* anAlphaChannel =
766 new G4AlphaDecay(theIon, b, c*MeV, a*MeV,
767 daughterFloatLevel);
768 // anAlphaChannel->DumpNuclearInfo();
769 theDecayTable->Insert(anAlphaChannel);
770 modeSumBR[Alpha] += b;
771 }
772 break;
773
774 case Proton:
775 {
776 G4ProtonDecay* aProtonChannel =
777 new G4ProtonDecay(theIon, b, c*MeV, a*MeV,
778 daughterFloatLevel);
779 // aProtonChannel->DumpNuclearInfo();
780 theDecayTable->Insert(aProtonChannel);
781 modeSumBR[Proton] += b;
782 }
783 break;
784
785 case Neutron:
786 {
787 G4NeutronDecay* aNeutronChannel =
788 new G4NeutronDecay(theIon, b, c*MeV, a*MeV,
789 daughterFloatLevel);
790 // aNeutronChannel->DumpNuclearInfo();
791 theDecayTable->Insert(aNeutronChannel);
792 modeSumBR[Neutron] += b;
793 }
794 break;
795
796 case SpFission:
797 {
798 G4SFDecay* aSpontFissChannel =
799 new G4SFDecay(theIon, b, c*MeV, a*MeV, daughterFloatLevel);
800 theDecayTable->Insert(aSpontFissChannel);
801 modeSumBR[SpFission] += b;
802 }
803 break;
804
805 case BDProton:
806 // Not yet implemented
807 // G4cout << " beta-delayed proton decay, a = " << a << ", b = " << b << ", c = " << c << G4endl;
808 break;
809
810 case BDNeutron:
811 // Not yet implemented
812 // G4cout << " beta-delayed neutron decay, a = " << a << ", b = " << b << ", c = " << c << G4endl;
813 break;
814
815 case Beta2Minus:
816 // Not yet implemented
817 // G4cout << " Double beta- decay, a = " << a << ", b = " << b << ", c = " << c << G4endl;
818 break;
819
820 case Beta2Plus:
821 // Not yet implemented
822 // G4cout << " Double beta+ decay, a = " << a << ", b = " << b << ", c = " << c << G4endl;
823 break;
824
825 case Proton2:
826 // Not yet implemented
827 // G4cout << " Double proton decay, a = " << a << ", b = " << b << ", c = " << c << G4endl;
828 break;
829
830 case Neutron2:
831 // Not yet implemented
832 // G4cout << " Double beta- decay, a = " << a << ", b = " << b << ", c = " << c << G4endl;
833 break;
834
835 case Triton:
836 {
837 G4TritonDecay* aTritonChannel =
838 new G4TritonDecay(theIon, b, c*MeV, a*MeV, daughterFloatLevel);
839 theDecayTable->Insert(aTritonChannel);
840 modeSumBR[Triton] += b;
841 }
842 break;
843
844 case RDM_ERROR:
845
846 default:
847 G4Exception("G4VRadioactiveDecay::LoadDecayTable()", "HAD_RDM_000",
848 FatalException, "Selected decay mode does not exist");
849 } // switch
850 } // line < 72
851 } // if char == P
852 } // if char != #
853 } // While
854
855 // Go through the decay table and make sure that the branching ratios are
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->entries(); ++i) {
864 theChannel = theDecayTable->GetDecayChannel(i);
865 theNuclearDecayChannel = static_cast<G4NuclearDecay*>(theChannel);
866 theDecayMode = theNuclearDecayChannel->GetDecayMode();
867
868 if (theDecayMode != IT) {
869 theBR = theChannel->GetBR();
870 theChannel->SetBR(theBR*modeTotalBR[theDecayMode]/modeSumBR[theDecayMode]);
871 }
872 }
873 } // decay file exists
874
875 DecaySchemeFile.close();
876
877 if (!found && levelEnergy > 0) {
878 // Case where IT cascade for excited isotopes has no entries in RDM database
879 // Decay mode is isomeric transition.
880 G4ITDecay* anITChannel = new G4ITDecay(theIon, 1.0, 0.0, 0.0);
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()] = theDecayTable;
890 lk.unlock();
891 return theDecayTable;
892 }
893
894 void G4VRadioactiveDecay::AddUserDecayDataFile(G4int Z, G4int A,
895 const G4String& filename)
896 {
897 if (Z < 1 || A < 2) G4cout << "Z and A not valid!" << G4endl;
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! " << G4endl;
906 G4Exception("G4VRadioactiveDecay::AddUserDecayDataFile()", "HAD_RDM_001",
907 FatalException, ed);
908 }
909 }
910
911 ////////////////////////////////////////////////////////////////////////////////
912 // //
913 // DecayIt //
914 // //
915 ////////////////////////////////////////////////////////////////////////////////
916
917 G4VParticleChange*
918 G4VRadioactiveDecay::DecayIt(const G4Track& theTrack, const G4Step&)
919 {
920 // Initialize G4ParticleChange object, get particle details and decay table
921 fParticleChangeForRadDecay.Initialize(theTrack);
922 fParticleChangeForRadDecay.ProposeWeight(theTrack.GetWeight());
923 const G4DynamicParticle* theParticle = theTrack.GetDynamicParticle();
924 const G4ParticleDefinition* theParticleDef = theParticle->GetDefinition();
925
926 // First check whether RDM applies to the current logical volume
927 if (!isAllVolumesMode) {
928 if (!std::binary_search(ValidVolumes.begin(), ValidVolumes.end(),
929 theTrack.GetVolume()->GetLogicalVolume()->GetName())) {
930 #ifdef G4VERBOSE
931 if (GetVerboseLevel()>1) {
932 G4cout <<"G4VRadioactiveDecay::DecayIt : "
933 << theTrack.GetVolume()->GetLogicalVolume()->GetName()
934 << " is not selected for the RDM"<< G4endl;
935 G4cout << " There are " << ValidVolumes.size() << " volumes" << G4endl;
936 G4cout << " The Valid volumes are: ";
937 for (auto const & vol : ValidVolumes) { G4cout << vol << " " << G4endl; }
938 G4cout << G4endl;
939 }
940 #endif
941 fParticleChangeForRadDecay.SetNumberOfSecondaries(0);
942
943 // Kill the parent particle.
944 fParticleChangeForRadDecay.ProposeTrackStatus(fStopAndKill) ;
945 fParticleChangeForRadDecay.ProposeLocalEnergyDeposit(0.0);
946 ClearNumberOfInteractionLengthLeft();
947 return &fParticleChangeForRadDecay;
948 }
949 }
950
951 // Now check if particle is valid for RDM
952 G4DecayTable* theDecayTable = GetDecayTable(theParticleDef);
953 if ( theDecayTable == nullptr || theDecayTable->entries() == 0) {
954 // Particle is not an ion or is outside the nucleuslimits for decay
955 #ifdef G4VERBOSE
956 if (GetVerboseLevel() > 1) {
957 G4cout << "G4VRadioactiveDecay::DecayIt : "
958 << theParticleDef->GetParticleName()
959 << " is outside (Z,A) limits set for the decay or has no decays."
960 << G4endl;
961 }
962 #endif
963 fParticleChangeForRadDecay.SetNumberOfSecondaries(0);
964
965 // Kill the parent particle
966 fParticleChangeForRadDecay.ProposeTrackStatus(fStopAndKill) ;
967 fParticleChangeForRadDecay.ProposeLocalEnergyDeposit(0.0);
968 ClearNumberOfInteractionLengthLeft();
969 return &fParticleChangeForRadDecay;
970 }
971 //G4cout << "DecayIt for " << theParticleDef->GetParticleName()
972 // << " isAllVolumesMode:" << isAllVolumesMode
973 // << " decayTable=" << theDecayTable << G4endl;
974
975 // Data found. Decay nucleus without variance reduction.
976 DecayAnalog(theTrack, theDecayTable);
977 return &fParticleChangeForRadDecay;
978 }
979
980 void G4VRadioactiveDecay::DecayAnalog(const G4Track& theTrack,
981 G4DecayTable* decayTable)
982 {
983 const G4DynamicParticle* theParticle = theTrack.GetDynamicParticle();
984 const G4ParticleDefinition* theParticleDef = theParticle->GetDefinition();
985 //G4cout << "DecayIt for " << theParticleDef->GetParticleName() << G4endl;
986 G4DecayProducts* products = DoDecay(*theParticleDef, decayTable);
987
988 // Check if the product is the same as input and kill the track if
989 // necessary to prevent infinite loop (11/05/10, F.Lei)
990 if (nullptr == products || products->entries() == 1) {
991 fParticleChangeForRadDecay.SetNumberOfSecondaries(0);
992 fParticleChangeForRadDecay.ProposeTrackStatus(fStopAndKill);
993 fParticleChangeForRadDecay.ProposeLocalEnergyDeposit(0.0);
994 ClearNumberOfInteractionLengthLeft();
995 delete products;
996 return;
997 }
998
999 G4double energyDeposit = 0.0;
1000 G4double finalGlobalTime = theTrack.GetGlobalTime();
1001 G4double finalLocalTime = theTrack.GetLocalTime();
1002
1003 // Get parent particle information and boost the decay products to the
1004 // laboratory frame
1005
1006 // ParentEnergy used for the boost should be the total energy of the nucleus
1007 // of the parent ion without the energy of the shell electrons
1008 // (correction for bug 1359 by L. Desorgher)
1009 G4double ParentEnergy = theParticle->GetKineticEnergy()
1010 + theParticle->GetParticleDefinition()->GetPDGMass();
1011 G4ThreeVector ParentDirection(theParticle->GetMomentumDirection());
1012
1013 if (theTrack.GetTrackStatus() == fStopButAlive) {
1014 // this condition seems to be always True, further investigation is needed (L.Desorgher)
1015
1016 // The particle is decayed at rest
1017 // Since the time is for the particle at rest, need to add additional time
1018 // lapsed between particle coming to rest and the actual decay. This time
1019 // is sampled with the mean-life of the particle. Need to protect the case
1020 // PDGTime < 0. (F.Lei 11/05/10)
1021 G4double temptime = -std::log(G4UniformRand() ) *
1022 theParticleDef->GetPDGLifeTime();
1023 if (temptime < 0.) temptime = 0.;
1024 finalGlobalTime += temptime;
1025 finalLocalTime += temptime;
1026 energyDeposit += theParticle->GetKineticEnergy();
1027
1028 // Kill the parent particle, and ignore its decay, if it decays later than the
1029 // threshold fThresholdForVeryLongDecayTime (whose default value corresponds
1030 // to more than twice the age of the universe).
1031 // This kind of cut has been introduced (in April 2021) in order to avoid to
1032 // account energy depositions happening after many billions of years in
1033 // ordinary materials used in calorimetry, in particular Tungsten and Lead
1034 // (via their natural unstable, but very long lived, isotopes, such as
1035 // W183, W180 and Pb204).
1036 // Note that the cut is not on the average, mean lifetime, but on the actual
1037 // sampled global decay time.
1038 if ( finalGlobalTime > fThresholdForVeryLongDecayTime ) {
1039 fParticleChangeForRadDecay.SetNumberOfSecondaries(0);
1040 fParticleChangeForRadDecay.ProposeTrackStatus(fStopAndKill) ;
1041 fParticleChangeForRadDecay.ProposeLocalEnergyDeposit(0.0);
1042 ClearNumberOfInteractionLengthLeft();
1043 delete products;
1044 return;
1045 }
1046 }
1047 products->Boost(ParentEnergy, ParentDirection);
1048
1049 // Add products in theParticleChangeForRadDecay.
1050 G4int numberOfSecondaries = products->entries();
1051 fParticleChangeForRadDecay.SetNumberOfSecondaries(numberOfSecondaries);
1052
1053 if (GetVerboseLevel() > 1) {
1054 G4cout << "G4VRadioactiveDecay::DecayAnalog: Decay vertex :";
1055 G4cout << " Time: " << finalGlobalTime/ns << "[ns]";
1056 G4cout << " X:" << (theTrack.GetPosition()).x() /cm << "[cm]";
1057 G4cout << " Y:" << (theTrack.GetPosition()).y() /cm << "[cm]";
1058 G4cout << " Z:" << (theTrack.GetPosition()).z() /cm << "[cm]";
1059 G4cout << G4endl;
1060 G4cout << "G4Decay::DecayIt : decay products in Lab. Frame" << G4endl;
1061 products->DumpInfo();
1062 products->IsChecked();
1063 }
1064
1065 const G4int modelID_forIT = G4PhysicsModelCatalog::GetModelID( "model_RDM_IT" );
1066 G4int modelID = modelID_forIT + 10*theRadDecayMode;
1067 const G4int modelID_forAtomicRelaxation =
1068 G4PhysicsModelCatalog::GetModelID( "model_RDM_AtomicRelaxation" );
1069 for ( G4int index = 0; index < numberOfSecondaries; ++index ) {
1070 G4Track* secondary = new G4Track( products->PopProducts(), finalGlobalTime,
1071 theTrack.GetPosition() );
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_forAtomicRelaxation) ;
1080 }
1081 } else if ( theRadDecayMode >= KshellEC && theRadDecayMode <= NshellEC &&
1082 index < numberOfSecondaries-1 ) {
1083 secondary->SetCreatorModelID( modelID_forAtomicRelaxation );
1084 }
1085 secondary->SetGoodForTrackingFlag();
1086 secondary->SetTouchableHandle( theTrack.GetTouchableHandle() );
1087 fParticleChangeForRadDecay.AddSecondary( secondary );
1088 }
1089
1090 delete products;
1091
1092 // Kill the parent particle
1093 fParticleChangeForRadDecay.ProposeTrackStatus(fStopAndKill) ;
1094 fParticleChangeForRadDecay.ProposeLocalEnergyDeposit(energyDeposit);
1095 fParticleChangeForRadDecay.ProposeLocalTime(finalLocalTime);
1096
1097 // Reset NumberOfInteractionLengthLeft.
1098 ClearNumberOfInteractionLengthLeft();
1099 }
1100
1101
1102 G4DecayProducts*
1103 G4VRadioactiveDecay::DoDecay(const G4ParticleDefinition& theParticleDef,
1104 G4DecayTable* theDecayTable)
1105 {
1106 G4DecayProducts* products = nullptr;
1107
1108 // Choose a decay channel.
1109 // G4DecayTable::SelectADecayChannel checks to see if sum of daughter masses
1110 // exceeds parent mass. Pass it the parent mass + maximum Q value to account
1111 // for difference in mass defect.
1112 G4double parentPlusQ = theParticleDef.GetPDGMass() + 30.*MeV;
1113 G4VDecayChannel* theDecayChannel = theDecayTable->SelectADecayChannel(parentPlusQ);
1114
1115 if (theDecayChannel == nullptr) {
1116 // Decay channel not found.
1117 G4ExceptionDescription ed;
1118 ed << " Cannot determine decay channel for " << theParticleDef.GetParticleName() << G4endl;
1119 G4Exception("G4VRadioactiveDecay::DoDecay", "HAD_RDM_013",
1120 FatalException, ed);
1121 } else {
1122 // A decay channel has been identified, so execute the DecayIt.
1123 #ifdef G4VERBOSE
1124 if (GetVerboseLevel() > 1) {
1125 G4cout << "G4VRadioactiveDecay::DoIt : selected decay channel addr: "
1126 << theDecayChannel << G4endl;
1127 }
1128 #endif
1129 theRadDecayMode = static_cast<G4NuclearDecay*>(theDecayChannel)->GetDecayMode();
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(theParticleDef.GetPDGMass());
1138 }
1139
1140 // Apply directional bias if requested by user
1141 CollimateDecay(products);
1142 }
1143
1144 return products;
1145 }
1146
1147
1148 // Apply directional bias for "visible" daughters (e+-, gamma, n, p, alpha)
1149 void G4VRadioactiveDecay::CollimateDecay(G4DecayProducts* products) {
1150
1151 if (origin == forceDecayDirection) return; // No collimation requested
1152 if (180.*deg == forceDecayHalfAngle) return;
1153 if (0 == products || 0 == products->entries()) return;
1154
1155 #ifdef G4VERBOSE
1156 if (GetVerboseLevel() > 1) G4cout << "Begin of CollimateDecay..." << G4endl;
1157 #endif
1158
1159 // Particles suitable for directional biasing (for if-blocks below)
1160 static const G4ParticleDefinition* electron = G4Electron::Definition();
1161 static const G4ParticleDefinition* positron = G4Positron::Definition();
1162 static const G4ParticleDefinition* neutron = G4Neutron::Definition();
1163 static const G4ParticleDefinition* gamma = G4Gamma::Definition();
1164 static const G4ParticleDefinition* alpha = G4Alpha::Definition();
1165 static const G4ParticleDefinition* triton = G4Triton::Definition();
1166 static const G4ParticleDefinition* proton = G4Proton::Definition();
1167
1168 G4ThreeVector newDirection; // Re-use to avoid memory churn
1169 for (G4int i=0; i<products->entries(); ++i) {
1170 G4DynamicParticle* daughter = (*products)[i];
1171 const G4ParticleDefinition* daughterType =
1172 daughter->GetParticleDefinition();
1173 if (daughterType == electron || daughterType == positron ||
1174 daughterType == neutron || daughterType == gamma ||
1175 daughterType == alpha || daughterType == triton || daughterType == proton) {
1176 CollimateDecayProduct(daughter);
1177 }
1178 }
1179 }
1180
1181 void G4VRadioactiveDecay::CollimateDecayProduct(G4DynamicParticle* daughter) {
1182 #ifdef G4VERBOSE
1183 if (GetVerboseLevel() > 1) {
1184 G4cout << "CollimateDecayProduct for daughter "
1185 << daughter->GetParticleDefinition()->GetParticleName() << G4endl;
1186 }
1187 #endif
1188
1189 G4ThreeVector collimate = ChooseCollimationDirection();
1190 if (origin != collimate) daughter->SetMomentumDirection(collimate);
1191 }
1192
1193
1194 // Choose random direction within collimation cone
1195
1196 G4ThreeVector G4VRadioactiveDecay::ChooseCollimationDirection() const {
1197 if (origin == forceDecayDirection) return origin; // Don't do collimation
1198 if (forceDecayHalfAngle == 180.*deg) return origin;
1199
1200 G4ThreeVector dir = forceDecayDirection;
1201
1202 // Return direction offset by random throw
1203 if (forceDecayHalfAngle > 0.) {
1204 // Generate uniform direction around central axis
1205 G4double phi = 2.*pi*G4UniformRand();
1206 G4double cosMin = std::cos(forceDecayHalfAngle);
1207 G4double cosTheta = (1.-cosMin)*G4UniformRand() + cosMin; // [cosMin,1.)
1208
1209 dir.setPhi(dir.phi()+phi);
1210 dir.setTheta(dir.theta()+std::acos(cosTheta));
1211 }
1212
1213 #ifdef G4VERBOSE
1214 if (GetVerboseLevel()>1)
1215 G4cout << " ChooseCollimationDirection returns " << dir << G4endl;
1216 #endif
1217
1218 return dir;
1219 }
1220
1221