| Geant4 Cross Reference |
1 // 1 //
2 // ******************************************* 2 // ********************************************************************
3 // * License and Disclaimer 3 // * License and Disclaimer *
4 // * 4 // * *
5 // * The Geant4 software is copyright of th 5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided 6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License 7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/ 8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. 9 // * include a list of copyright holders. *
10 // * 10 // * *
11 // * Neither the authors of this software syst 11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing fin 12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warran 13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assum 14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file 15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitatio 16 // * for the full disclaimer and the limitation of liability. *
17 // * 17 // * *
18 // * This code implementation is the result 18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboratio 19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distri 20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you ag 21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publicati 22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Sof 23 // * acceptance of all terms of the Geant4 Software license. *
24 // ******************************************* 24 // ********************************************************************
25 // 25 //
26 ////////////////////////////////////////////// << 26 // MODULE: G4RadioactiveDecay.cc
27 // 27 //
28 // GEANT4 Class source file << 28 // Author: F Lei & P R Truscott
>> 29 // Organisation: DERA UK
>> 30 // Customer: ESA/ESTEC, NOORDWIJK
>> 31 // Contract: 12115/96/JG/NL Work Order No. 3
29 // 32 //
30 // G4RadioactiveDecay << 33 // Documentation avaialable at http://www.space.dera.gov.uk/space_env/rdm.html
>> 34 // These include:
>> 35 // User Requirement Document (URD)
>> 36 // Software Specification Documents (SSD)
>> 37 // Software User Manual (SUM)
>> 38 // Technical Note (TN) on the physics and algorithms
31 // 39 //
32 // Author: D.H. Wright (SLAC) << 40 // The test and example programs are not included in the public release of
33 // Date: 29 August 2017 << 41 // G4 but they can be downloaded from
>> 42 // http://www.space.qinetiq.com/space_env/rdm.html
>> 43 //
>> 44 // CHANGE HISTORY
>> 45 // --------------
34 // 46 //
35 // Based on the code of F. Lei and P.R. Trusc << 47 // 13 Oct 2015, L.G. Sarmiento Neutron emission added
>> 48 //
>> 49 // 06 Aug 2014, L.G. Sarmiento Proton decay mode added mimicking the alpha decay
>> 50 //
>> 51 // 03 Oct 2012, V. Ivanchenko removed internal table for mean free path
>> 52 // similar to what is done for as G4Decay
>> 53 // 10 July 2012, L. Desorgher
>> 54 // -In LoadDecayTable:
>> 55 // Add LoadedNuclei.push_back(theParentNucleus.GetParticleName());
>> 56 // also for the case where user data files are used. Correction for bug
>> 57 // 1324. Changes proposed by Joa L.
>> 58 //
>> 59 //
>> 60 // 01 May 2012, L. Desorgher
>> 61 // -Force the reading of user file to theIsotopeTable
>> 62 // -Merge the development by Fan Lei for activation computation
>> 63 //
>> 64 // 17 Oct 2011, L. Desorgher
>> 65 // -Add possibility for the user to load its own decay file.
>> 66 // -Set halflifethreshold negative by default to allow the tracking of all
>> 67 // excited nuclei resulting from a radioactive decay
>> 68 //
>> 69 // 01 June 2011, M. Kelsey -- Add directional biasing interface to allow for
>> 70 // "collimation" of decay daughters.
>> 71 // 16 February 2006, V.Ivanchenko fix problem in IsApplicable connected with
>> 72 // 8.0 particle design
>> 73 // 18 October 2002, F. Lei
>> 74 // in the case of beta decay, added a check of the end-energy
>> 75 // to ensure it is > 0.
>> 76 // ENSDF occationally have beta decay entries with zero energies
>> 77 //
>> 78 // 27 Sepetember 2001, F. Lei
>> 79 // verboselevel(0) used in constructor
>> 80 //
>> 81 // 01 November 2000, F.Lei
>> 82 // added " ee = e0 +1. ;" as line 763
>> 83 // tagged as "radiative_decay-V02-00-02"
>> 84 // 28 October 2000, F Lei
>> 85 // added fast beta decay mode. Many files have been changed.
>> 86 // tagged as "radiative_decay-V02-00-01"
>> 87 //
>> 88 // 25 October 2000, F Lei, DERA UK
>> 89 // 1) line 1185 added 'const' to work with tag "Track-V02-00-00"
>> 90 // tagged as "radiative_decay-V02-00-00"
>> 91 // 14 April 2000, F Lei, DERA UK
>> 92 // 0.b.4 release. Changes are:
>> 93 // 1) Use PhotonEvaporation instead of DiscreteGammaDeexcitation
>> 94 // 2) VR: Significant efficiency improvement
>> 95 //
>> 96 // 29 February 2000, P R Truscott, DERA UK
>> 97 // 0.b.3 release.
>> 98 //
>> 99 ///////////////////////////////////////////////////////////////////////////////
36 // 100 //
37 ////////////////////////////////////////////// <<
38 <<
39 #include "G4RadioactiveDecay.hh" 101 #include "G4RadioactiveDecay.hh"
40 #include "G4RadioactivationMessenger.hh" << 102 #include "G4RadioactiveDecaymessenger.hh"
41 103
42 #include "G4SystemOfUnits.hh" 104 #include "G4SystemOfUnits.hh"
43 #include "G4DynamicParticle.hh" 105 #include "G4DynamicParticle.hh"
44 #include "G4DecayProducts.hh" 106 #include "G4DecayProducts.hh"
45 #include "G4DecayTable.hh" 107 #include "G4DecayTable.hh"
46 #include "G4ParticleChangeForRadDecay.hh" 108 #include "G4ParticleChangeForRadDecay.hh"
47 #include "G4ITDecay.hh" 109 #include "G4ITDecay.hh"
48 #include "G4BetaDecayType.hh" 110 #include "G4BetaDecayType.hh"
49 #include "G4BetaMinusDecay.hh" 111 #include "G4BetaMinusDecay.hh"
50 #include "G4BetaPlusDecay.hh" 112 #include "G4BetaPlusDecay.hh"
51 #include "G4ECDecay.hh" 113 #include "G4ECDecay.hh"
52 #include "G4AlphaDecay.hh" 114 #include "G4AlphaDecay.hh"
53 #include "G4TritonDecay.hh" 115 #include "G4TritonDecay.hh"
54 #include "G4ProtonDecay.hh" 116 #include "G4ProtonDecay.hh"
55 #include "G4NeutronDecay.hh" 117 #include "G4NeutronDecay.hh"
56 #include "G4SFDecay.hh" 118 #include "G4SFDecay.hh"
57 #include "G4VDecayChannel.hh" 119 #include "G4VDecayChannel.hh"
58 #include "G4NuclearDecay.hh" 120 #include "G4NuclearDecay.hh"
59 #include "G4RadioactiveDecayMode.hh" 121 #include "G4RadioactiveDecayMode.hh"
60 #include "G4Fragment.hh" 122 #include "G4Fragment.hh"
61 #include "G4Ions.hh" 123 #include "G4Ions.hh"
62 #include "G4IonTable.hh" 124 #include "G4IonTable.hh"
63 #include "G4BetaDecayType.hh" 125 #include "G4BetaDecayType.hh"
64 #include "Randomize.hh" 126 #include "Randomize.hh"
65 #include "G4LogicalVolumeStore.hh" 127 #include "G4LogicalVolumeStore.hh"
66 #include "G4NuclearLevelData.hh" 128 #include "G4NuclearLevelData.hh"
67 #include "G4DeexPrecoParameters.hh" 129 #include "G4DeexPrecoParameters.hh"
>> 130 #include "G4EmParameters.hh"
68 #include "G4LevelManager.hh" 131 #include "G4LevelManager.hh"
69 #include "G4ThreeVector.hh" 132 #include "G4ThreeVector.hh"
70 #include "G4Electron.hh" 133 #include "G4Electron.hh"
71 #include "G4Positron.hh" 134 #include "G4Positron.hh"
72 #include "G4Neutron.hh" 135 #include "G4Neutron.hh"
73 #include "G4Gamma.hh" 136 #include "G4Gamma.hh"
74 #include "G4Alpha.hh" 137 #include "G4Alpha.hh"
75 #include "G4Triton.hh" 138 #include "G4Triton.hh"
76 #include "G4Proton.hh" 139 #include "G4Proton.hh"
77 140
78 #include "G4HadronicProcessType.hh" 141 #include "G4HadronicProcessType.hh"
79 #include "G4HadronicProcessStore.hh" 142 #include "G4HadronicProcessStore.hh"
80 #include "G4HadronicException.hh" 143 #include "G4HadronicException.hh"
81 #include "G4LossTableManager.hh" <<
82 #include "G4VAtomDeexcitation.hh" <<
83 #include "G4UAtomicDeexcitation.hh" <<
84 #include "G4PhotonEvaporation.hh" 144 #include "G4PhotonEvaporation.hh"
85 145
86 #include <vector> 146 #include <vector>
87 #include <sstream> 147 #include <sstream>
88 #include <algorithm> 148 #include <algorithm>
89 #include <fstream> 149 #include <fstream>
90 150
91 using namespace CLHEP; 151 using namespace CLHEP;
92 152
93 G4RadioactiveDecay::G4RadioactiveDecay(const G << 153 const G4double G4RadioactiveDecay::levelTolerance = 10.0*eV;
94 const G << 154 const G4ThreeVector G4RadioactiveDecay::origin(0.,0.,0.);
95 : G4VRadioactiveDecay(processName, timeThres << 155
>> 156 #ifdef G4MULTITHREADED
>> 157 #include "G4AutoLock.hh"
>> 158 G4Mutex G4RadioactiveDecay::radioactiveDecayMutex = G4MUTEX_INITIALIZER;
>> 159 DecayTableMap* G4RadioactiveDecay::master_dkmap = 0;
>> 160
>> 161 G4int& G4RadioactiveDecay::NumberOfInstances()
>> 162 {
>> 163 static G4int numberOfInstances = 0;
>> 164 return numberOfInstances;
>> 165 }
>> 166 #endif
>> 167
>> 168 G4RadioactiveDecay::G4RadioactiveDecay(const G4String& processName)
>> 169 : G4VRestDiscreteProcess(processName, fDecay), isInitialised(false),
>> 170 forceDecayDirection(0.,0.,0.), forceDecayHalfAngle(0.*deg), dirPath(""),
>> 171 verboseLevel(0)
96 { 172 {
97 #ifdef G4VERBOSE 173 #ifdef G4VERBOSE
98 if (GetVerboseLevel() > 1) { 174 if (GetVerboseLevel() > 1) {
99 G4cout << "G4RadioactiveDecay constructor: 175 G4cout << "G4RadioactiveDecay constructor: processName = " << processName
100 << G4endl; 176 << G4endl;
101 } 177 }
102 #endif 178 #endif
103 179
104 theRadioactivationMessenger = new G4Radioact << 180 SetProcessSubType(fRadioactiveDecay);
>> 181
>> 182 theRadioactiveDecaymessenger = new G4RadioactiveDecaymessenger(this);
>> 183 pParticleChange = &fParticleChangeForRadDecay;
>> 184
>> 185 // Set up photon evaporation for use in G4ITDecay
>> 186 photonEvaporation = new G4PhotonEvaporation();
>> 187 // photonEvaporation->SetVerboseLevel(2);
>> 188 photonEvaporation->RDMForced(true);
>> 189 photonEvaporation->SetICM(true);
>> 190
>> 191 // Check data directory
>> 192 char* path_var = std::getenv("G4RADIOACTIVEDATA");
>> 193 if (!path_var) {
>> 194 G4Exception("G4RadioactiveDecay()","HAD_RDM_200",FatalException,
>> 195 "Environment variable G4RADIOACTIVEDATA is not set");
>> 196 } else {
>> 197 dirPath = path_var; // convert to string
>> 198 std::ostringstream os;
>> 199 os << dirPath << "/z1.a3"; // used as a dummy
>> 200 std::ifstream testFile;
>> 201 testFile.open(os.str() );
>> 202 if (!testFile.is_open() )
>> 203 G4Exception("G4RadioactiveDecay()","HAD_RDM_201",FatalException,
>> 204 "Environment variable G4RADIOACTIVEDATA is set, but does not point to correct directory");
>> 205 }
>> 206
>> 207 // Reset the list of user defined data files
>> 208 theUserRadioactiveDataFiles.clear();
>> 209
>> 210 // Instantiate the map of decay tables
>> 211 #ifdef G4MULTITHREADED
>> 212 G4AutoLock lk(&G4RadioactiveDecay::radioactiveDecayMutex);
>> 213 NumberOfInstances()++;
>> 214 if(!master_dkmap) master_dkmap = new DecayTableMap;
>> 215 #endif
>> 216 dkmap = new DecayTableMap;
105 217
106 // Apply default values. 218 // Apply default values.
107 NSourceBin = 1; 219 NSourceBin = 1;
108 SBin[0] = 0.* s; 220 SBin[0] = 0.* s;
109 SBin[1] = 1.* s; // Convert to ns << 221 SBin[1] = 1.* s;
110 SProfile[0] = 1.; 222 SProfile[0] = 1.;
111 SProfile[1] = 0.; 223 SProfile[1] = 0.;
112 NDecayBin = 1; 224 NDecayBin = 1;
113 DBin[0] = 0. * s ; 225 DBin[0] = 0. * s ;
114 DBin[1] = 1. * s; 226 DBin[1] = 1. * s;
115 DProfile[0] = 1.; 227 DProfile[0] = 1.;
116 DProfile[1] = 0.; 228 DProfile[1] = 0.;
117 decayWindows[0] = 0; 229 decayWindows[0] = 0;
118 G4RadioactivityTable* rTable = new G4Radioac 230 G4RadioactivityTable* rTable = new G4RadioactivityTable() ;
119 theRadioactivityTables.push_back(rTable); 231 theRadioactivityTables.push_back(rTable);
120 NSplit = 1; 232 NSplit = 1;
121 AnalogueMC = true; << 233 AnalogueMC = true ;
122 BRBias = true; << 234 FBeta = false ;
123 halflifethreshold = 1000.*nanosecond; << 235 BRBias = true ;
>> 236 applyICM = true ;
>> 237 applyARM = true ;
>> 238 halflifethreshold = nanosecond;
>> 239
>> 240 // RDM applies to all logical volumes by default
>> 241 isAllVolumesMode = true;
>> 242 SelectAllVolumes();
>> 243 G4HadronicProcessStore::Instance()->RegisterExtraProcess(this);
124 } 244 }
125 245
>> 246
126 void G4RadioactiveDecay::ProcessDescription(st 247 void G4RadioactiveDecay::ProcessDescription(std::ostream& outFile) const
127 { 248 {
128 outFile << "The G4RadioactiveDecay process p << 249 outFile << "The radioactive decay process (G4RadioactiveDecay) handles the\n"
129 << "nuclides (G4GenericIon) in biase << 250 << "alpha, beta+, beta-, electron capture and isomeric transition\n"
130 << "duplication, branching ratio bia << 251 << "decays of nuclei (G4GenericIon) with masses A > 4.\n"
131 << "and detector time convolution. <<
132 << "activation physics.\n" <<
133 << "The required half-lives and deca 252 << "The required half-lives and decay schemes are retrieved from\n"
134 << "the RadioactiveDecay database wh 253 << "the RadioactiveDecay database which was derived from ENSDF.\n";
135 } 254 }
136 255
137 256
138 G4RadioactiveDecay::~G4RadioactiveDecay() 257 G4RadioactiveDecay::~G4RadioactiveDecay()
139 { 258 {
140 delete theRadioactivationMessenger; << 259 delete theRadioactiveDecaymessenger;
>> 260 delete photonEvaporation;
>> 261 for (DecayTableMap::iterator i = dkmap->begin(); i != dkmap->end(); i++) {
>> 262 delete i->second;
>> 263 }
>> 264 dkmap->clear();
>> 265 delete dkmap;
>> 266 #ifdef G4MULTITHREADED
>> 267 G4AutoLock lk(&G4RadioactiveDecay::radioactiveDecayMutex);
>> 268 --NumberOfInstances();
>> 269 if(NumberOfInstances()==0)
>> 270 {
>> 271 for (DecayTableMap::iterator i = master_dkmap->begin(); i != master_dkmap->end(); i++) {
>> 272 delete i->second;
>> 273 }
>> 274 master_dkmap->clear();
>> 275 delete master_dkmap;
>> 276 }
>> 277 #endif
>> 278 }
>> 279
>> 280
>> 281 G4bool G4RadioactiveDecay::IsApplicable(const G4ParticleDefinition& aParticle)
>> 282 {
>> 283 // All particles other than G4Ions are rejected by default
>> 284 if (((const G4Ions*)(&aParticle))->GetExcitationEnergy() > 0.) {
>> 285 return true; // Not ground state - decay
>> 286 }
>> 287
>> 288 if (aParticle.GetParticleName() == "GenericIon") {
>> 289 return true;
>> 290 } else if (!(aParticle.GetParticleType() == "nucleus")
>> 291 || aParticle.GetPDGLifeTime() < 0. ) {
>> 292 return false; // Nuclide is stable - no decay
>> 293 }
>> 294
>> 295 // At this point nuclide must be an unstable ground state
>> 296 // Determine whether it falls into the correct A and Z range
>> 297 G4int A = ((const G4Ions*) (&aParticle))->GetAtomicMass();
>> 298 G4int Z = ((const G4Ions*) (&aParticle))->GetAtomicNumber();
>> 299
>> 300 if (A > theNucleusLimits.GetAMax() || A < theNucleusLimits.GetAMin())
>> 301 {return false;}
>> 302 else if (Z > theNucleusLimits.GetZMax() || Z < theNucleusLimits.GetZMin())
>> 303 {return false;}
>> 304 return true;
141 } 305 }
142 306
>> 307 G4DecayTable* G4RadioactiveDecay::GetDecayTable(const G4ParticleDefinition* aNucleus)
>> 308 {
>> 309 G4String key = aNucleus->GetParticleName();
>> 310 DecayTableMap::iterator table_ptr = dkmap->find(key);
>> 311
>> 312 G4DecayTable* theDecayTable = 0;
>> 313 if (table_ptr == dkmap->end() ) { // If table not there,
>> 314 theDecayTable = LoadDecayTable(*aNucleus); // load from file and
>> 315 if(theDecayTable) (*dkmap)[key] = theDecayTable; // store in library
>> 316 } else {
>> 317 theDecayTable = table_ptr->second;
>> 318 }
>> 319
>> 320 return theDecayTable;
>> 321 }
>> 322
>> 323
>> 324 void G4RadioactiveDecay::SelectAVolume(const G4String aVolume)
>> 325 {
>> 326 G4LogicalVolumeStore* theLogicalVolumes;
>> 327 G4LogicalVolume* volume;
>> 328 theLogicalVolumes = G4LogicalVolumeStore::GetInstance();
>> 329 for (size_t i = 0; i < theLogicalVolumes->size(); i++) {
>> 330 volume=(*theLogicalVolumes)[i];
>> 331 if (volume->GetName() == aVolume) {
>> 332 ValidVolumes.push_back(aVolume);
>> 333 std::sort(ValidVolumes.begin(), ValidVolumes.end());
>> 334 // sort need for performing binary_search
>> 335 #ifdef G4VERBOSE
>> 336 if (GetVerboseLevel()>0)
>> 337 G4cout << " RDM Applies to : " << aVolume << G4endl;
>> 338 #endif
>> 339 } else if(i == theLogicalVolumes->size()) {
>> 340 G4cerr << "SelectAVolume: "<< aVolume
>> 341 << " is not a valid logical volume name" << G4endl;
>> 342 }
>> 343 }
>> 344 }
>> 345
>> 346
>> 347 void G4RadioactiveDecay::DeselectAVolume(const G4String aVolume)
>> 348 {
>> 349 G4LogicalVolumeStore* theLogicalVolumes;
>> 350 G4LogicalVolume* volume;
>> 351 theLogicalVolumes = G4LogicalVolumeStore::GetInstance();
>> 352 for (size_t i = 0; i < theLogicalVolumes->size(); i++){
>> 353 volume=(*theLogicalVolumes)[i];
>> 354 if (volume->GetName() == aVolume) {
>> 355 std::vector<G4String>::iterator location;
>> 356 location = std::find(ValidVolumes.begin(),ValidVolumes.end(),aVolume);
>> 357 if (location != ValidVolumes.end()) {
>> 358 ValidVolumes.erase(location);
>> 359 std::sort(ValidVolumes.begin(), ValidVolumes.end());
>> 360 isAllVolumesMode =false;
>> 361 } else {
>> 362 G4cerr << " DeselectVolume:" << aVolume << " is not in the list "
>> 363 << G4endl;
>> 364 }
>> 365 #ifdef G4VERBOSE
>> 366 if (GetVerboseLevel() > 0)
>> 367 G4cout << " DeselectVolume: " << aVolume << " is removed from list "
>> 368 << G4endl;
>> 369 #endif
>> 370 } else if (i == theLogicalVolumes->size()) {
>> 371 G4cerr << " DeselectVolume:" << aVolume
>> 372 << "is not a valid logical volume name" << G4endl;
>> 373 }
>> 374 }
>> 375 }
>> 376
>> 377
>> 378 void G4RadioactiveDecay::SelectAllVolumes()
>> 379 {
>> 380 G4LogicalVolumeStore* theLogicalVolumes;
>> 381 G4LogicalVolume* volume;
>> 382 theLogicalVolumes = G4LogicalVolumeStore::GetInstance();
>> 383 ValidVolumes.clear();
>> 384 #ifdef G4VERBOSE
>> 385 if (GetVerboseLevel()>0)
>> 386 G4cout << " RDM Applies to all Volumes" << G4endl;
>> 387 #endif
>> 388 for (size_t i = 0; i < theLogicalVolumes->size(); i++){
>> 389 volume = (*theLogicalVolumes)[i];
>> 390 ValidVolumes.push_back(volume->GetName());
>> 391 #ifdef G4VERBOSE
>> 392 if (GetVerboseLevel()>0)
>> 393 G4cout << " RDM Applies to Volume " << volume->GetName() << G4endl;
>> 394 #endif
>> 395 }
>> 396 std::sort(ValidVolumes.begin(), ValidVolumes.end());
>> 397 // sort needed in order to allow binary_search
>> 398 isAllVolumesMode=true;
>> 399 }
>> 400
>> 401
>> 402 void G4RadioactiveDecay::DeselectAllVolumes()
>> 403 {
>> 404 ValidVolumes.clear();
>> 405 isAllVolumesMode=false;
>> 406 #ifdef G4VERBOSE
>> 407 if (GetVerboseLevel() > 0) G4cout << "RDM removed from all volumes" << G4endl;
>> 408 #endif
>> 409 }
>> 410
>> 411
143 G4bool 412 G4bool
144 G4RadioactiveDecay::IsRateTableReady(const G4P 413 G4RadioactiveDecay::IsRateTableReady(const G4ParticleDefinition& aParticle)
145 { 414 {
146 // Check whether the radioactive decay rates 415 // Check whether the radioactive decay rates table for the ion has already
147 // been calculated. 416 // been calculated.
148 G4String aParticleName = aParticle.GetPartic 417 G4String aParticleName = aParticle.GetParticleName();
149 for (std::size_t i = 0; i < theParentChainTa << 418 for (size_t i = 0; i < theParentChainTable.size(); i++) {
150 if (theParentChainTable[i].GetIonName() == 419 if (theParentChainTable[i].GetIonName() == aParticleName) return true;
151 } 420 }
152 return false; 421 return false;
153 } 422 }
154 423
>> 424 // retrieve the decayratetable for the specified aParticle
155 void 425 void
156 G4RadioactiveDecay::GetChainsFromParent(const 426 G4RadioactiveDecay::GetChainsFromParent(const G4ParticleDefinition& aParticle)
157 { 427 {
158 // Retrieve the decay rate table for the spe <<
159 G4String aParticleName = aParticle.GetPartic 428 G4String aParticleName = aParticle.GetParticleName();
160 429
161 for (std::size_t i = 0; i < theParentChainTa << 430 for (size_t i = 0; i < theParentChainTable.size(); i++) {
162 if (theParentChainTable[i].GetIonName() == 431 if (theParentChainTable[i].GetIonName() == aParticleName) {
163 theDecayRateVector = theParentChainTable 432 theDecayRateVector = theParentChainTable[i].GetItsRates();
164 } 433 }
165 } 434 }
166 #ifdef G4VERBOSE 435 #ifdef G4VERBOSE
167 if (GetVerboseLevel() > 1) { << 436 if (GetVerboseLevel() > 0) {
168 G4cout << "The DecayRate Table for " << aP 437 G4cout << "The DecayRate Table for " << aParticleName << " is selected."
169 << G4endl; 438 << G4endl;
170 } 439 }
171 #endif 440 #endif
172 } 441 }
173 442
>> 443 /* DHW: long double version - only few % improvement, but don't delete yet
>> 444 G4double
>> 445 G4RadioactiveDecay::ConvolveSourceTimeProfile(const G4double t, const G4double tau)
>> 446 {
>> 447 long double convolvedTime = 0.L;
>> 448 G4int nbin;
>> 449 if ( t > SBin[NSourceBin]) {
>> 450 nbin = NSourceBin;
>> 451 } else {
>> 452 nbin = 0;
>> 453
>> 454 G4int loop = 0;
>> 455 while (t > SBin[nbin]) {
>> 456 loop++;
>> 457 if (loop > 1000) {
>> 458 G4Exception("G4RadioactiveDecay::ConvolveSourceTimeProfile()",
>> 459 "HAD_RDM_100", JustWarning, "While loop count exceeded");
>> 460 break;
>> 461 }
>> 462
>> 463 nbin++;
>> 464 }
>> 465 nbin--;
>> 466 }
>> 467
>> 468 long double lt = t ;
>> 469 long double ltau = tau;
>> 470 long double earg = 0.L;
>> 471 if (nbin > 0) {
>> 472 for (G4int i = 0; i < nbin; i++) {
>> 473 earg = (long double)(SBin[i+1] - SBin[i])/ltau;
>> 474 if (earg < 100.) {
>> 475 convolvedTime += (long double)SProfile[i] *
>> 476 std::exp(((long double)SBin[i] - lt)/ltau) *
>> 477 std::expm1(earg);
>> 478 } else {
>> 479 convolvedTime += (long double)SProfile[i] *
>> 480 (std::exp(-(lt-(long double)SBin[i+1])/ltau)-std::exp(-(lt-(long double)SBin[i])/ltau));
>> 481 }
>> 482 }
>> 483 }
>> 484 // Use -expm1 instead of 1 - exp
>> 485 convolvedTime -= (long double)SProfile[nbin] * std::expm1(((long double)SBin[nbin] - lt)/ltau);
>> 486
>> 487 if (convolvedTime < 0.) {
>> 488 G4cout << " Convolved time =: " << convolvedTime << " reset to zero! " << G4endl;
>> 489 G4cout << " t = " << t << " tau = " << tau << G4endl;
>> 490 G4cout << SBin[nbin] << " " << SBin[0] << G4endl;
>> 491 convolvedTime = 0.;
>> 492 }
>> 493 #ifdef G4VERBOSE
>> 494 if (GetVerboseLevel() > 1)
>> 495 G4cout << " Convolved time: " << convolvedTime << G4endl;
>> 496 #endif
>> 497 return (G4double)convolvedTime;
>> 498 }
>> 499 */
>> 500
174 // ConvolveSourceTimeProfile performs the conv 501 // ConvolveSourceTimeProfile performs the convolution of the source time profile
175 // function with a single exponential characte 502 // function with a single exponential characterized by a decay constant in the
176 // decay chain. The time profile is treated a << 503 // decay chain. The time profile is treated as a set of step functions so that
177 // convolution integral can be done bin-by-bin << 504 // the convolution integral can be done bin-by-bin.
178 // This implements Eq. 4.13 of DERA technical 505 // This implements Eq. 4.13 of DERA technical note, with SProfile[i] = F(t')
179 506
180 G4double 507 G4double
181 G4RadioactiveDecay::ConvolveSourceTimeProfile( 508 G4RadioactiveDecay::ConvolveSourceTimeProfile(const G4double t, const G4double tau)
182 { 509 {
183 G4double convolvedTime = 0.0; 510 G4double convolvedTime = 0.0;
184 G4int nbin; 511 G4int nbin;
185 if ( t > SBin[NSourceBin]) { 512 if ( t > SBin[NSourceBin]) {
186 nbin = NSourceBin; 513 nbin = NSourceBin;
187 } else { 514 } else {
188 nbin = 0; 515 nbin = 0;
189 516
190 G4int loop = 0; 517 G4int loop = 0;
191 while (t > SBin[nbin]) { // Loop checking 518 while (t > SBin[nbin]) { // Loop checking, 01.09.2015, D.Wright
192 loop++; 519 loop++;
193 if (loop > 1000) { 520 if (loop > 1000) {
194 G4Exception("G4RadioactiveDecay::Convo 521 G4Exception("G4RadioactiveDecay::ConvolveSourceTimeProfile()",
195 "HAD_RDM_100", JustWarning << 522 "HAD_RDM_100", JustWarning,"While loop count exceeded");
196 break; 523 break;
197 } 524 }
198 ++nbin; << 525 nbin++;
199 } 526 }
200 --nbin; << 527 nbin--;
201 } 528 }
202 529
203 // Use expm1 wherever possible to avoid larg 530 // Use expm1 wherever possible to avoid large cancellation errors in
204 // 1 - exp(x) for small x << 531 // 1 - exp(x) for small x
205 G4double earg = 0.0; 532 G4double earg = 0.0;
206 if (nbin > 0) { 533 if (nbin > 0) {
207 for (G4int i = 0; i < nbin; ++i) { << 534 for (G4int i = 0; i < nbin; i++) {
208 earg = (SBin[i+1] - SBin[i])/tau; 535 earg = (SBin[i+1] - SBin[i])/tau;
209 if (earg < 100.) { 536 if (earg < 100.) {
210 convolvedTime += SProfile[i] * std::ex 537 convolvedTime += SProfile[i] * std::exp((SBin[i] - t)/tau) *
211 std::expm1(earg); 538 std::expm1(earg);
212 } else { 539 } else {
213 convolvedTime += SProfile[i] * 540 convolvedTime += SProfile[i] *
214 (std::exp(-(t-SBin[i+1])/tau)-std::e 541 (std::exp(-(t-SBin[i+1])/tau)-std::exp(-(t-SBin[i])/tau));
215 } 542 }
216 } 543 }
217 } 544 }
218 convolvedTime -= SProfile[nbin] * std::expm1 545 convolvedTime -= SProfile[nbin] * std::expm1((SBin[nbin] - t)/tau);
219 // tau divided out of final result to provid 546 // tau divided out of final result to provide probability of decay in window
220 547
221 if (convolvedTime < 0.) { 548 if (convolvedTime < 0.) {
222 G4cout << " Convolved time =: " << convolv 549 G4cout << " Convolved time =: " << convolvedTime << " reset to zero! " << G4endl;
223 G4cout << " t = " << t << " tau = " << tau 550 G4cout << " t = " << t << " tau = " << tau << G4endl;
224 G4cout << SBin[nbin] << " " << SBin[0] << 551 G4cout << SBin[nbin] << " " << SBin[0] << G4endl;
225 convolvedTime = 0.; 552 convolvedTime = 0.;
226 } 553 }
227 #ifdef G4VERBOSE 554 #ifdef G4VERBOSE
228 if (GetVerboseLevel() > 2) << 555 if (GetVerboseLevel() > 1)
229 G4cout << " Convolved time: " << convolved 556 G4cout << " Convolved time: " << convolvedTime << G4endl;
230 #endif 557 #endif
231 return convolvedTime; 558 return convolvedTime;
232 } 559 }
233 560
234 561
235 ////////////////////////////////////////////// 562 ////////////////////////////////////////////////////////////////////////////////
236 // 563 // //
237 // GetDecayTime 564 // GetDecayTime //
238 // Randomly select a decay time for the dec 565 // Randomly select a decay time for the decay process, following the //
239 // supplied decay time bias scheme. 566 // supplied decay time bias scheme. //
240 // 567 // //
241 ////////////////////////////////////////////// 568 ////////////////////////////////////////////////////////////////////////////////
242 569
243 G4double G4RadioactiveDecay::GetDecayTime() 570 G4double G4RadioactiveDecay::GetDecayTime()
244 { 571 {
245 G4double decaytime = 0.; 572 G4double decaytime = 0.;
246 G4double rand = G4UniformRand(); 573 G4double rand = G4UniformRand();
247 G4int i = 0; 574 G4int i = 0;
248 575
249 G4int loop = 0; 576 G4int loop = 0;
250 while (DProfile[i] < rand) { /* Loop checki << 577 while ( DProfile[i] < rand) { /* Loop checking, 01.09.2015, D.Wright */
251 // Entries in DProfile[i] are all between << 578 i++;
252 // Comparison with rand chooses which time <<
253 ++i; <<
254 loop++; 579 loop++;
255 if (loop > 100000) { 580 if (loop > 100000) {
256 G4Exception("G4RadioactiveDecay::GetDeca 581 G4Exception("G4RadioactiveDecay::GetDecayTime()", "HAD_RDM_100",
257 JustWarning, "While loop cou 582 JustWarning, "While loop count exceeded");
258 break; 583 break;
259 } 584 }
260 } 585 }
261 586
262 rand = G4UniformRand(); 587 rand = G4UniformRand();
263 decaytime = DBin[i] + rand*(DBin[i+1]-DBin[i 588 decaytime = DBin[i] + rand*(DBin[i+1]-DBin[i]);
264 #ifdef G4VERBOSE 589 #ifdef G4VERBOSE
265 if (GetVerboseLevel() > 2) << 590 if (GetVerboseLevel() > 1)
266 G4cout <<" Decay time: " <<decaytime/s <<" 591 G4cout <<" Decay time: " <<decaytime/s <<"[s]" <<G4endl;
267 #endif 592 #endif
268 return decaytime; 593 return decaytime;
269 } 594 }
270 595
271 596
272 G4int G4RadioactiveDecay::GetDecayTimeBin(cons 597 G4int G4RadioactiveDecay::GetDecayTimeBin(const G4double aDecayTime)
273 { 598 {
274 G4int i = 0; 599 G4int i = 0;
275 600
276 G4int loop = 0; 601 G4int loop = 0;
277 while (aDecayTime > DBin[i] ) { /* Loop ch << 602 while ( aDecayTime > DBin[i] ) { /* Loop checking, 01.09.2015, D.Wright */
278 ++i; << 603 i++;
279 loop++; 604 loop++;
280 if (loop > 100000) { 605 if (loop > 100000) {
281 G4Exception("G4RadioactiveDecay::GetDeca << 606 G4Exception("G4RadioactiveDecay::GetDecayTimeBin()", "HAD_RDM_100",
282 JustWarning, "While loop cou 607 JustWarning, "While loop count exceeded");
283 break; 608 break;
284 } 609 }
285 } 610 }
286 611
287 return i; 612 return i;
288 } 613 }
289 614
290 ////////////////////////////////////////////// 615 ////////////////////////////////////////////////////////////////////////////////
291 // 616 // //
292 // GetMeanLifeTime (required by the base clas 617 // GetMeanLifeTime (required by the base class) //
293 // 618 // //
294 ////////////////////////////////////////////// 619 ////////////////////////////////////////////////////////////////////////////////
295 620
296 G4double G4RadioactiveDecay::GetMeanLifeTime(c 621 G4double G4RadioactiveDecay::GetMeanLifeTime(const G4Track& theTrack,
297 G4 << 622 G4ForceCondition*)
298 { 623 {
299 // For variance reduction time is set to 0 s << 624 // For variance reduction the time is set to 0 so as to force the particle
300 // to decay immediately. 625 // to decay immediately.
301 // In analogue mode it returns the particle' << 626 // In analogueMC mode it returns the particle's mean-life.
>> 627
302 G4double meanlife = 0.; 628 G4double meanlife = 0.;
303 if (AnalogueMC) meanlife = G4VRadioactiveDec << 629 if (AnalogueMC) {
304 return meanlife; << 630 const G4DynamicParticle* theParticle = theTrack.GetDynamicParticle();
>> 631 const G4ParticleDefinition* theParticleDef = theParticle->GetDefinition();
>> 632 G4double theLife = theParticleDef->GetPDGLifeTime();
>> 633 #ifdef G4VERBOSE
>> 634 if (GetVerboseLevel() > 2) {
>> 635 G4cout << "G4RadioactiveDecay::GetMeanLifeTime() " << G4endl;
>> 636 G4cout << "KineticEnergy: " << theParticle->GetKineticEnergy()/GeV
>> 637 << " GeV, Mass: " << theParticle->GetMass()/GeV
>> 638 << " GeV, Life time: " << theLife/ns << " ns " << G4endl;
>> 639 }
>> 640 #endif
>> 641 if (theParticleDef->GetPDGStable()) {meanlife = DBL_MAX;}
>> 642 else if (theLife < 0.0) {meanlife = DBL_MAX;}
>> 643 else {meanlife = theLife;}
>> 644 // Set meanlife to zero for excited istopes which are not in the
>> 645 // RDM database
>> 646 if (((const G4Ions*)(theParticleDef))->GetExcitationEnergy() > 0. &&
>> 647 meanlife == DBL_MAX) {meanlife = 0.;}
>> 648 }
>> 649 #ifdef G4VERBOSE
>> 650 if (GetVerboseLevel() > 1)
>> 651 G4cout << " mean life time: " << meanlife/s << " s " << G4endl;
>> 652 #endif
>> 653
>> 654 return meanlife;
>> 655 }
>> 656
>> 657 ////////////////////////////////////////////////////////////////////////
>> 658 // //
>> 659 // GetMeanFreePath for decay in flight //
>> 660 // //
>> 661 ////////////////////////////////////////////////////////////////////////
>> 662
>> 663 G4double G4RadioactiveDecay::GetMeanFreePath (const G4Track& aTrack, G4double,
>> 664 G4ForceCondition*)
>> 665 {
>> 666 const G4DynamicParticle* aParticle = aTrack.GetDynamicParticle();
>> 667 const G4ParticleDefinition* aParticleDef = aParticle->GetDefinition();
>> 668 G4double tau = aParticleDef->GetPDGLifeTime();
>> 669 G4double aMass = aParticle->GetMass();
>> 670
>> 671 #ifdef G4VERBOSE
>> 672 if (GetVerboseLevel() > 2) {
>> 673 G4cout << "G4RadioactiveDecay::GetMeanFreePath() " << G4endl;
>> 674 G4cout << " KineticEnergy: " << aParticle->GetKineticEnergy()/GeV
>> 675 << " GeV, Mass: " << aMass/GeV << " GeV, tau: " << tau << " ns "
>> 676 << G4endl;
>> 677 }
>> 678 #endif
>> 679 G4double pathlength = DBL_MAX;
>> 680 if (tau != -1) {
>> 681 // Ion can decay
>> 682
>> 683 if (tau < -1000.0) {
>> 684 pathlength = DBL_MIN; // nuclide had very short lifetime or wasn't in table
>> 685
>> 686 } else if (tau < 0.0) {
>> 687 G4cout << aParticleDef->GetParticleName() << " has lifetime " << tau << G4endl;
>> 688 G4ExceptionDescription ed;
>> 689 ed << "Ion has negative lifetime " << tau
>> 690 << " but is not stable. Setting mean free path to DBL_MAX" << G4endl;
>> 691 G4Exception("G4RadioactiveDecay::GetMeanFreePath()", "HAD_RDM_011",
>> 692 JustWarning, ed);
>> 693 pathlength = DBL_MAX;
>> 694
>> 695 } else {
>> 696 // Calculate mean free path
>> 697 G4double betaGamma = aParticle->GetTotalMomentum()/aMass;
>> 698 pathlength = c_light*tau*betaGamma;
>> 699
>> 700 if (pathlength < DBL_MIN) {
>> 701 pathlength = DBL_MIN;
>> 702 #ifdef G4VERBOSE
>> 703 if (GetVerboseLevel() > 2) {
>> 704 G4cout << "G4Decay::GetMeanFreePath: "
>> 705 << aParticleDef->GetParticleName()
>> 706 << " stops, kinetic energy = "
>> 707 << aParticle->GetKineticEnergy()/keV <<" keV " << G4endl;
>> 708 }
>> 709 #endif
>> 710 }
>> 711 }
>> 712 }
>> 713
>> 714 #ifdef G4VERBOSE
>> 715 if (GetVerboseLevel() > 1) {
>> 716 G4cout << "mean free path: "<< pathlength/m << " m" << G4endl;
>> 717 }
>> 718 #endif
>> 719 return pathlength;
>> 720 }
>> 721
>> 722 ////////////////////////////////////////////////////////////////////////
>> 723 // //
>> 724 // BuildPhysicsTable - enable print of parameters //
>> 725 // //
>> 726 ////////////////////////////////////////////////////////////////////////
>> 727
>> 728 void G4RadioactiveDecay::BuildPhysicsTable(const G4ParticleDefinition&)
>> 729 {
>> 730 if (!isInitialised) {
>> 731 isInitialised = true;
>> 732 #ifdef G4VERBOSE
>> 733 if(G4Threading::IsMasterThread()) { StreamInfo(G4cout, "\n"); }
>> 734 #endif
>> 735 }
>> 736 G4HadronicProcessStore::
>> 737 Instance()->RegisterParticleForExtraProcess(this,G4GenericIon::GenericIon());
>> 738 }
>> 739
>> 740 ////////////////////////////////////////////////////////////////////////
>> 741 // //
>> 742 // StreamInfo - stream out parameters //
>> 743 // //
>> 744 ////////////////////////////////////////////////////////////////////////
>> 745
>> 746 void G4RadioactiveDecay::StreamInfo(std::ostream& os, const G4String& endline)
>> 747 {
>> 748 G4DeexPrecoParameters* deex =
>> 749 G4NuclearLevelData::GetInstance()->GetParameters();
>> 750 G4EmParameters* emparam = G4EmParameters::Instance();
>> 751
>> 752 G4int prec = os.precision(5);
>> 753 os << "======================================================================="
>> 754 << endline;
>> 755 os << "====== Radioactive Decay Physics Parameters ========"
>> 756 << endline;
>> 757 os << "======================================================================="
>> 758 << endline;
>> 759 os << "Max life time "
>> 760 << deex->GetMaxLifeTime()/CLHEP::ps << " ps" << endline;
>> 761 os << "Internal e- conversion flag "
>> 762 << deex->GetInternalConversionFlag() << endline;
>> 763 os << "Stored internal conversion coefficients "
>> 764 << deex->StoreICLevelData() << endline;
>> 765 os << "Enable correlated gamma emission "
>> 766 << deex->CorrelatedGamma() << endline;
>> 767 os << "Max 2J for sampling of angular correlations "
>> 768 << deex->GetTwoJMAX() << endline;
>> 769 os << "Atomic de-excitation enabled "
>> 770 << emparam->Fluo() << endline;
>> 771 os << "Auger electron emission enabled "
>> 772 << emparam->Auger() << endline;
>> 773 os << "Auger cascade enabled "
>> 774 << emparam->AugerCascade() << endline;
>> 775 os << "Check EM cuts disabled for atomic de-excitation "
>> 776 << emparam->DeexcitationIgnoreCut() << endline;
>> 777 os << "Use Bearden atomic level energies "
>> 778 << emparam->BeardenFluoDir() << endline;
>> 779 os << "======================================================================="
>> 780 << endline;
>> 781 os.precision(prec);
>> 782 }
>> 783
>> 784 ////////////////////////////////////////////////////////////////////////////////
>> 785 // //
>> 786 // LoadDecayTable loads the decay scheme from the RadioactiveDecay database //
>> 787 // for the parent nucleus. //
>> 788 // //
>> 789 ////////////////////////////////////////////////////////////////////////////////
>> 790
>> 791 G4DecayTable*
>> 792 G4RadioactiveDecay::LoadDecayTable(const G4ParticleDefinition& theParentNucleus)
>> 793 {
>> 794 // Generate input data file name using Z and A of the parent nucleus
>> 795 // file containing radioactive decay data.
>> 796 G4int A = ((const G4Ions*)(&theParentNucleus))->GetAtomicMass();
>> 797 G4int Z = ((const G4Ions*)(&theParentNucleus))->GetAtomicNumber();
>> 798
>> 799 G4double levelEnergy = ((const G4Ions*)(&theParentNucleus))->GetExcitationEnergy();
>> 800 G4Ions::G4FloatLevelBase floatingLevel =
>> 801 ((const G4Ions*)(&theParentNucleus))->GetFloatLevelBase();
>> 802
>> 803 #ifdef G4MULTITHREADED
>> 804 G4AutoLock lk(&G4RadioactiveDecay::radioactiveDecayMutex);
>> 805
>> 806 G4String key = theParentNucleus.GetParticleName();
>> 807 DecayTableMap::iterator master_table_ptr = master_dkmap->find(key);
>> 808
>> 809 if (master_table_ptr != master_dkmap->end() ) { // If table is there
>> 810 return master_table_ptr->second;
>> 811 }
>> 812 #endif
>> 813
>> 814 //Check if data have been provided by the user
>> 815 G4String file = theUserRadioactiveDataFiles[1000*A+Z];
>> 816
>> 817 if (file == "") {
>> 818 std::ostringstream os;
>> 819 os << dirPath << "/z" << Z << ".a" << A << '\0';
>> 820 file = os.str();
>> 821 }
>> 822
>> 823 G4DecayTable* theDecayTable = new G4DecayTable();
>> 824 G4bool found(false); // True if energy level matches one in table
>> 825
>> 826 std::ifstream DecaySchemeFile;
>> 827 DecaySchemeFile.open(file);
>> 828
>> 829 if (DecaySchemeFile.good()) {
>> 830 // Initialize variables used for reading in radioactive decay data
>> 831 G4bool floatMatch(false);
>> 832 const G4int nMode = 12;
>> 833 G4double modeTotalBR[nMode] = {0.0};
>> 834 G4double modeSumBR[nMode];
>> 835 for (G4int i = 0; i < nMode; i++) {
>> 836 modeSumBR[i] = 0.0;
>> 837 }
>> 838
>> 839 char inputChars[120]={' '};
>> 840 G4String inputLine;
>> 841 G4String recordType("");
>> 842 G4String floatingFlag("");
>> 843 G4String daughterFloatFlag("");
>> 844 G4Ions::G4FloatLevelBase daughterFloatLevel;
>> 845 G4RadioactiveDecayMode theDecayMode;
>> 846 G4double decayModeTotal(0.0);
>> 847 G4double parentExcitation(0.0);
>> 848 G4double a(0.0);
>> 849 G4double b(0.0);
>> 850 G4double c(0.0);
>> 851 G4double dummy(0.0);
>> 852 G4BetaDecayType betaType(allowed);
>> 853
>> 854 // Loop through each data file record until you identify the decay
>> 855 // data relating to the nuclide of concern.
>> 856
>> 857 G4bool complete(false); // bool insures only one set of values read for any
>> 858 // given parent energy level
>> 859 G4int loop = 0;
>> 860 while (!complete && !DecaySchemeFile.getline(inputChars, 120).eof()) { /* Loop checking, 01.09.2015, D.Wright */
>> 861 loop++;
>> 862 if (loop > 100000) {
>> 863 G4Exception("G4RadioactiveDecay::LoadDecayTable()", "HAD_RDM_100",
>> 864 JustWarning, "While loop count exceeded");
>> 865 break;
>> 866 }
>> 867
>> 868 inputLine = inputChars;
>> 869 inputLine = inputLine.strip(1);
>> 870 if (inputChars[0] != '#' && inputLine.length() != 0) {
>> 871 std::istringstream tmpStream(inputLine);
>> 872
>> 873 if (inputChars[0] == 'P') {
>> 874 // Nucleus is a parent type. Check excitation level to see if it
>> 875 // matches that of theParentNucleus
>> 876 tmpStream >> recordType >> parentExcitation >> floatingFlag >> dummy;
>> 877 // "dummy" takes the place of half-life
>> 878 // Now read in from ENSDFSTATE in particle category
>> 879
>> 880 if (found) {
>> 881 complete = true;
>> 882 } else {
>> 883 // Take first level which matches excitation energy regardless of floating level
>> 884 found = (std::abs(parentExcitation*keV - levelEnergy) < levelTolerance);
>> 885 if (floatingLevel != noFloat) {
>> 886 // If floating level specificed, require match of both energy and floating level
>> 887 floatMatch = (floatingLevel == G4Ions::FloatLevelBase(floatingFlag.back()) );
>> 888 if (!floatMatch) found = false;
>> 889 }
>> 890 }
>> 891
>> 892 } else if (found) {
>> 893 // The right part of the radioactive decay data file has been found. Search
>> 894 // through it to determine the mode of decay of the subsequent records.
>> 895
>> 896 // Store for later the total decay probability for each decay mode
>> 897 if (inputLine.length() < 72) {
>> 898 tmpStream >> theDecayMode >> dummy >> decayModeTotal;
>> 899 switch (theDecayMode) {
>> 900 case IT:
>> 901 {
>> 902 G4ITDecay* anITChannel = new G4ITDecay(&theParentNucleus, decayModeTotal,
>> 903 0.0, 0.0, photonEvaporation);
>> 904 anITChannel->SetHLThreshold(halflifethreshold);
>> 905 anITChannel->SetARM(applyARM);
>> 906 theDecayTable->Insert(anITChannel);
>> 907 // anITChannel->DumpNuclearInfo();
>> 908 }
>> 909 break;
>> 910 case BetaMinus:
>> 911 modeTotalBR[1] = decayModeTotal; break;
>> 912 case BetaPlus:
>> 913 modeTotalBR[2] = decayModeTotal; break;
>> 914 case KshellEC:
>> 915 modeTotalBR[3] = decayModeTotal; break;
>> 916 case LshellEC:
>> 917 modeTotalBR[4] = decayModeTotal; break;
>> 918 case MshellEC:
>> 919 modeTotalBR[5] = decayModeTotal; break;
>> 920 case NshellEC:
>> 921 modeTotalBR[6] = decayModeTotal; break;
>> 922 case Alpha:
>> 923 modeTotalBR[7] = decayModeTotal; break;
>> 924 case Proton:
>> 925 modeTotalBR[8] = decayModeTotal; break;
>> 926 case Neutron:
>> 927 modeTotalBR[9] = decayModeTotal; break;
>> 928 case BDProton:
>> 929 break;
>> 930 case BDNeutron:
>> 931 break;
>> 932 case Beta2Minus:
>> 933 break;
>> 934 case Beta2Plus:
>> 935 break;
>> 936 case Proton2:
>> 937 break;
>> 938 case Neutron2:
>> 939 break;
>> 940 case SpFission:
>> 941 modeTotalBR[10] = decayModeTotal; break;
>> 942 case Triton:
>> 943 modeTotalBR[11] = decayModeTotal; break;
>> 944 case RDM_ERROR:
>> 945
>> 946 default:
>> 947 G4Exception("G4RadioactiveDecay::LoadDecayTable()", "HAD_RDM_000",
>> 948 FatalException, "Selected decay mode does not exist");
>> 949 } // switch
>> 950
>> 951 } else {
>> 952 if (inputLine.length() < 84) {
>> 953 tmpStream >> theDecayMode >> a >> daughterFloatFlag >> b >> c;
>> 954 betaType = allowed;
>> 955 } else {
>> 956 tmpStream >> theDecayMode >> a >> daughterFloatFlag >> b >> c >> betaType;
>> 957 }
>> 958
>> 959 // Allowed transitions are the default. Forbidden transitions are
>> 960 // indicated in the last column.
>> 961 a /= 1000.;
>> 962 c /= 1000.;
>> 963 b/=100.;
>> 964 daughterFloatLevel = G4Ions::FloatLevelBase(daughterFloatFlag.back());
>> 965
>> 966 switch (theDecayMode) {
>> 967 case BetaMinus:
>> 968 {
>> 969 G4BetaMinusDecay* aBetaMinusChannel =
>> 970 new G4BetaMinusDecay(&theParentNucleus, b, c*MeV, a*MeV,
>> 971 daughterFloatLevel, betaType);
>> 972 // aBetaMinusChannel->DumpNuclearInfo();
>> 973 aBetaMinusChannel->SetHLThreshold(halflifethreshold);
>> 974 theDecayTable->Insert(aBetaMinusChannel);
>> 975 modeSumBR[1] += b;
>> 976 }
>> 977 break;
>> 978
>> 979 case BetaPlus:
>> 980 {
>> 981 G4BetaPlusDecay* aBetaPlusChannel =
>> 982 new G4BetaPlusDecay(&theParentNucleus, b, c*MeV, a*MeV,
>> 983 daughterFloatLevel, betaType);
>> 984 // aBetaPlusChannel->DumpNuclearInfo();
>> 985 aBetaPlusChannel->SetHLThreshold(halflifethreshold);
>> 986 theDecayTable->Insert(aBetaPlusChannel);
>> 987 modeSumBR[2] += b;
>> 988 }
>> 989 break;
>> 990
>> 991 case KshellEC: // K-shell electron capture
>> 992 {
>> 993 G4ECDecay* aKECChannel =
>> 994 new G4ECDecay(&theParentNucleus, b, c*MeV, a*MeV,
>> 995 daughterFloatLevel, KshellEC);
>> 996 // aKECChannel->DumpNuclearInfo();
>> 997 aKECChannel->SetHLThreshold(halflifethreshold);
>> 998 aKECChannel->SetARM(applyARM);
>> 999 theDecayTable->Insert(aKECChannel);
>> 1000 modeSumBR[3] += b;
>> 1001 }
>> 1002 break;
>> 1003
>> 1004 case LshellEC: // L-shell electron capture
>> 1005 {
>> 1006 G4ECDecay* aLECChannel =
>> 1007 new G4ECDecay(&theParentNucleus, b, c*MeV, a*MeV,
>> 1008 daughterFloatLevel, LshellEC);
>> 1009 // aLECChannel->DumpNuclearInfo();
>> 1010 aLECChannel->SetHLThreshold(halflifethreshold);
>> 1011 aLECChannel->SetARM(applyARM);
>> 1012 theDecayTable->Insert(aLECChannel);
>> 1013 modeSumBR[4] += b;
>> 1014 }
>> 1015 break;
>> 1016
>> 1017 case MshellEC: // M-shell electron capture
>> 1018 {
>> 1019 G4ECDecay* aMECChannel =
>> 1020 new G4ECDecay(&theParentNucleus, b, c*MeV, a*MeV,
>> 1021 daughterFloatLevel, MshellEC);
>> 1022 // aMECChannel->DumpNuclearInfo();
>> 1023 aMECChannel->SetHLThreshold(halflifethreshold);
>> 1024 aMECChannel->SetARM(applyARM);
>> 1025 theDecayTable->Insert(aMECChannel);
>> 1026 modeSumBR[5] += b;
>> 1027 }
>> 1028 break;
>> 1029
>> 1030 case NshellEC: // N-shell electron capture
>> 1031 {
>> 1032 G4ECDecay* aNECChannel =
>> 1033 new G4ECDecay(&theParentNucleus, b, c*MeV, a*MeV,
>> 1034 daughterFloatLevel, NshellEC);
>> 1035 // aNECChannel->DumpNuclearInfo();
>> 1036 aNECChannel->SetHLThreshold(halflifethreshold);
>> 1037 aNECChannel->SetARM(applyARM);
>> 1038 theDecayTable->Insert(aNECChannel);
>> 1039 modeSumBR[6] += b;
>> 1040 }
>> 1041 break;
>> 1042
>> 1043 case Alpha:
>> 1044 {
>> 1045 G4AlphaDecay* anAlphaChannel =
>> 1046 new G4AlphaDecay(&theParentNucleus, b, c*MeV, a*MeV,
>> 1047 daughterFloatLevel);
>> 1048 // anAlphaChannel->DumpNuclearInfo();
>> 1049 anAlphaChannel->SetHLThreshold(halflifethreshold);
>> 1050 theDecayTable->Insert(anAlphaChannel);
>> 1051 modeSumBR[7] += b;
>> 1052 }
>> 1053 break;
>> 1054
>> 1055 case Proton:
>> 1056 {
>> 1057 G4ProtonDecay* aProtonChannel =
>> 1058 new G4ProtonDecay(&theParentNucleus, b, c*MeV, a*MeV,
>> 1059 daughterFloatLevel);
>> 1060 // aProtonChannel->DumpNuclearInfo();
>> 1061 aProtonChannel->SetHLThreshold(halflifethreshold);
>> 1062 theDecayTable->Insert(aProtonChannel);
>> 1063 modeSumBR[8] += b;
>> 1064 }
>> 1065 break;
>> 1066
>> 1067 case Neutron:
>> 1068 {
>> 1069 G4NeutronDecay* aNeutronChannel =
>> 1070 new G4NeutronDecay(&theParentNucleus, b, c*MeV, a*MeV,
>> 1071 daughterFloatLevel);
>> 1072 // aNeutronChannel->DumpNuclearInfo();
>> 1073 aNeutronChannel->SetHLThreshold(halflifethreshold);
>> 1074 theDecayTable->Insert(aNeutronChannel);
>> 1075 modeSumBR[9] += b;
>> 1076 }
>> 1077 break;
>> 1078
>> 1079 case BDProton:
>> 1080 // Not yet implemented
>> 1081 // G4cout << " beta-delayed proton decay, a = " << a << ", b = " << b << ", c = " << c << G4endl;
>> 1082 break;
>> 1083 case BDNeutron:
>> 1084 // Not yet implemented
>> 1085 // G4cout << " beta-delayed neutron decay, a = " << a << ", b = " << b << ", c = " << c << G4endl;
>> 1086 break;
>> 1087 case Beta2Minus:
>> 1088 // Not yet implemented
>> 1089 // G4cout << " Double beta- decay, a = " << a << ", b = " << b << ", c = " << c << G4endl;
>> 1090 break;
>> 1091 case Beta2Plus:
>> 1092 // Not yet implemented
>> 1093 // G4cout << " Double beta+ decay, a = " << a << ", b = " << b << ", c = " << c << G4endl;
>> 1094 break;
>> 1095 case Proton2:
>> 1096 // Not yet implemented
>> 1097 // G4cout << " Double proton decay, a = " << a << ", b = " << b << ", c = " << c << G4endl;
>> 1098 break;
>> 1099 case Neutron2:
>> 1100 // Not yet implemented
>> 1101 // G4cout << " Double beta- decay, a = " << a << ", b = " << b << ", c = " << c << G4endl;
>> 1102 break;
>> 1103 case SpFission:
>> 1104 {
>> 1105 G4SFDecay* aSpontFissChannel =
>> 1106 // new G4SFDecay(&theParentNucleus, decayModeTotal, 0.0, 0.0);
>> 1107 new G4SFDecay(&theParentNucleus, b, c*MeV, a*MeV,
>> 1108 daughterFloatLevel);
>> 1109 theDecayTable->Insert(aSpontFissChannel);
>> 1110 modeSumBR[10] += b;
>> 1111 }
>> 1112 break;
>> 1113 case Triton:
>> 1114 {
>> 1115 G4TritonDecay* aTritonChannel =
>> 1116 new G4TritonDecay(&theParentNucleus, b, c*MeV, a*MeV,
>> 1117 daughterFloatLevel);
>> 1118 // anTritonChannel->DumpNuclearInfo();
>> 1119 aTritonChannel->SetHLThreshold(halflifethreshold);
>> 1120 theDecayTable->Insert(aTritonChannel);
>> 1121 modeSumBR[11] += b;
>> 1122 }
>> 1123 break;
>> 1124
>> 1125 case RDM_ERROR:
>> 1126
>> 1127 default:
>> 1128 G4Exception("G4RadioactiveDecay::LoadDecayTable()", "HAD_RDM_000",
>> 1129 FatalException, "Selected decay mode does not exist");
>> 1130 } // switch
>> 1131 } // line < 72
>> 1132 } // if char == P
>> 1133 } // if char != #
>> 1134 } // While
>> 1135
>> 1136 // Go through the decay table and make sure that the branching ratios are
>> 1137 // correctly normalised.
>> 1138
>> 1139 G4VDecayChannel* theChannel = 0;
>> 1140 G4NuclearDecay* theNuclearDecayChannel = 0;
>> 1141 G4String mode = "";
>> 1142
>> 1143 G4double theBR = 0.0;
>> 1144 for (G4int i = 0; i < theDecayTable->entries(); i++) {
>> 1145 theChannel = theDecayTable->GetDecayChannel(i);
>> 1146 theNuclearDecayChannel = static_cast<G4NuclearDecay*>(theChannel);
>> 1147 theDecayMode = theNuclearDecayChannel->GetDecayMode();
>> 1148
>> 1149 if (theDecayMode != IT) {
>> 1150 theBR = theChannel->GetBR();
>> 1151 theChannel->SetBR(theBR*modeTotalBR[theDecayMode]/modeSumBR[theDecayMode]);
>> 1152 }
>> 1153 }
>> 1154 } // decay file exists
>> 1155
>> 1156 DecaySchemeFile.close();
>> 1157
>> 1158 if (!found && levelEnergy > 0) {
>> 1159 // Case where IT cascade for excited isotopes has no entries in RDM database
>> 1160 // Decay mode is isomeric transition.
>> 1161 G4ITDecay* anITChannel = new G4ITDecay(&theParentNucleus, 1.0, 0.0, 0.0,
>> 1162 photonEvaporation);
>> 1163 anITChannel->SetHLThreshold(halflifethreshold);
>> 1164 anITChannel->SetARM(applyARM);
>> 1165 theDecayTable->Insert(anITChannel);
>> 1166 }
>> 1167
>> 1168 if (theDecayTable && GetVerboseLevel() > 1) {
>> 1169 theDecayTable->DumpInfo();
>> 1170 }
>> 1171
>> 1172 #ifdef G4MULTITHREADED
>> 1173 //(*master_dkmap)[key] = theDecayTable; // store in master library
>> 1174 #endif
>> 1175 return theDecayTable;
>> 1176 }
>> 1177
>> 1178 void
>> 1179 G4RadioactiveDecay::AddUserDecayDataFile(G4int Z, G4int A, G4String filename)
>> 1180 {
>> 1181 if (Z < 1 || A < 2) G4cout << "Z and A not valid!" << G4endl;
>> 1182
>> 1183 std::ifstream DecaySchemeFile(filename);
>> 1184 if (DecaySchemeFile) {
>> 1185 G4int ID_ion = A*1000 + Z;
>> 1186 theUserRadioactiveDataFiles[ID_ion] = filename;
>> 1187 } else {
>> 1188 G4cout << "The file " << filename << " does not exist!" << G4endl;
>> 1189 }
305 } 1190 }
306 1191
307 1192
308 void 1193 void
309 G4RadioactiveDecay::SetDecayRate(G4int theZ, G 1194 G4RadioactiveDecay::SetDecayRate(G4int theZ, G4int theA, G4double theE,
310 G4int theG, std::vector<G4double>& the << 1195 G4int theG, std::vector<G4double> theCoefficients,
311 std::vector<G4double>& theTaos) << 1196 std::vector<G4double> theTaos)
312 // Why not make this a method of G4Radioactiv 1197 // Why not make this a method of G4RadioactiveDecayRate? (e.g. SetParameters)
313 { 1198 {
314 //fill the decay rate vector 1199 //fill the decay rate vector
315 ratesToDaughter.SetZ(theZ); 1200 ratesToDaughter.SetZ(theZ);
316 ratesToDaughter.SetA(theA); 1201 ratesToDaughter.SetA(theA);
317 ratesToDaughter.SetE(theE); 1202 ratesToDaughter.SetE(theE);
318 ratesToDaughter.SetGeneration(theG); 1203 ratesToDaughter.SetGeneration(theG);
319 ratesToDaughter.SetDecayRateC(theCoefficient 1204 ratesToDaughter.SetDecayRateC(theCoefficients);
320 ratesToDaughter.SetTaos(theTaos); 1205 ratesToDaughter.SetTaos(theTaos);
321 } 1206 }
322 1207
323 1208
324 void G4RadioactiveDecay:: 1209 void G4RadioactiveDecay::
325 CalculateChainsFromParent(const G4ParticleDefi 1210 CalculateChainsFromParent(const G4ParticleDefinition& theParentNucleus)
326 { 1211 {
327 // Use extended Bateman equation to calculat 1212 // Use extended Bateman equation to calculate the radioactivities of all
328 // progeny of theParentNucleus. The coeffic 1213 // progeny of theParentNucleus. The coefficients required to do this are
329 // calculated using the method of P. Truscot 1214 // calculated using the method of P. Truscott (Ph.D. thesis and
330 // DERA Technical Note DERA/CIS/CIS2/7/36/4/ 1215 // DERA Technical Note DERA/CIS/CIS2/7/36/4/10) 11 January 2000.
331 // Coefficients are then added to the decay 1216 // Coefficients are then added to the decay rate table vector
332 1217
333 // Create and initialise variables used in t 1218 // Create and initialise variables used in the method.
334 theDecayRateVector.clear(); 1219 theDecayRateVector.clear();
335 1220
336 G4int nGeneration = 0; 1221 G4int nGeneration = 0;
337 1222
338 std::vector<G4double> taos; 1223 std::vector<G4double> taos;
339 1224
340 // Dimensionless A coefficients of Eqs. 4.24 1225 // Dimensionless A coefficients of Eqs. 4.24 and 4.25 of the TN
341 std::vector<G4double> Acoeffs; 1226 std::vector<G4double> Acoeffs;
342 1227
343 // According to Eq. 4.26 the first coefficie 1228 // According to Eq. 4.26 the first coefficient (A_1:1) is -1
344 Acoeffs.push_back(-1.); 1229 Acoeffs.push_back(-1.);
345 1230
346 const G4Ions* ion = static_cast<const G4Ions << 1231 G4int A = ((const G4Ions*)(&theParentNucleus))->GetAtomicMass();
347 G4int A = ion->GetAtomicMass(); << 1232 G4int Z = ((const G4Ions*)(&theParentNucleus))->GetAtomicNumber();
348 G4int Z = ion->GetAtomicNumber(); << 1233 G4double E = ((const G4Ions*)(&theParentNucleus))->GetExcitationEnergy();
349 G4double E = ion->GetExcitationEnergy(); << 1234 G4double tao = theParentNucleus.GetPDGLifeTime();
350 G4double tao = ion->GetPDGLifeTime(); <<
351 if (tao < 0.) tao = 1e-100; 1235 if (tao < 0.) tao = 1e-100;
352 taos.push_back(tao); 1236 taos.push_back(tao);
353 G4int nEntry = 0; 1237 G4int nEntry = 0;
354 1238
355 // Fill the decay rate container (G4Radioact << 1239 // Fill the decay rate container (G4RadioactiveDecayRatesToDaughter) with
356 // isotope data << 1240 // the parent isotope data
357 SetDecayRate(Z,A,E,nGeneration,Acoeffs,taos) << 1241 SetDecayRate(Z,A,E,nGeneration,Acoeffs,taos);
358 1242
359 // store the decay rate in decay rate vector 1243 // store the decay rate in decay rate vector
360 theDecayRateVector.push_back(ratesToDaughter 1244 theDecayRateVector.push_back(ratesToDaughter);
361 ++nEntry; << 1245 nEntry++;
362 1246
363 // Now start treating the secondary generati 1247 // Now start treating the secondary generations.
364 G4bool stable = false; 1248 G4bool stable = false;
>> 1249 G4int i;
365 G4int j; 1250 G4int j;
366 G4VDecayChannel* theChannel = 0; 1251 G4VDecayChannel* theChannel = 0;
367 G4NuclearDecay* theNuclearDecayChannel = 0; 1252 G4NuclearDecay* theNuclearDecayChannel = 0;
368 1253
369 G4ITDecay* theITChannel = 0; 1254 G4ITDecay* theITChannel = 0;
370 G4BetaMinusDecay* theBetaMinusChannel = 0; 1255 G4BetaMinusDecay* theBetaMinusChannel = 0;
371 G4BetaPlusDecay* theBetaPlusChannel = 0; 1256 G4BetaPlusDecay* theBetaPlusChannel = 0;
372 G4AlphaDecay* theAlphaChannel = 0; 1257 G4AlphaDecay* theAlphaChannel = 0;
373 G4ProtonDecay* theProtonChannel = 0; 1258 G4ProtonDecay* theProtonChannel = 0;
374 G4TritonDecay* theTritonChannel = 0; 1259 G4TritonDecay* theTritonChannel = 0;
375 G4NeutronDecay* theNeutronChannel = 0; 1260 G4NeutronDecay* theNeutronChannel = 0;
376 G4SFDecay* theFissionChannel = 0; 1261 G4SFDecay* theFissionChannel = 0;
377 1262
378 G4RadioactiveDecayMode theDecayMode; 1263 G4RadioactiveDecayMode theDecayMode;
379 G4double theBR = 0.0; 1264 G4double theBR = 0.0;
380 G4int AP = 0; 1265 G4int AP = 0;
381 G4int ZP = 0; 1266 G4int ZP = 0;
382 G4int AD = 0; 1267 G4int AD = 0;
383 G4int ZD = 0; 1268 G4int ZD = 0;
384 G4double EP = 0.; 1269 G4double EP = 0.;
385 std::vector<G4double> TP; 1270 std::vector<G4double> TP;
386 std::vector<G4double> RP; // A coefficient 1271 std::vector<G4double> RP; // A coefficients of the previous generation
387 G4ParticleDefinition *theDaughterNucleus; 1272 G4ParticleDefinition *theDaughterNucleus;
388 G4double daughterExcitation; 1273 G4double daughterExcitation;
389 G4double nearestEnergy = 0.0; 1274 G4double nearestEnergy = 0.0;
390 G4int nearestLevelIndex = 0; 1275 G4int nearestLevelIndex = 0;
391 G4ParticleDefinition *aParentNucleus; 1276 G4ParticleDefinition *aParentNucleus;
392 G4IonTable* theIonTable; 1277 G4IonTable* theIonTable;
393 G4DecayTable* parentDecayTable; 1278 G4DecayTable* parentDecayTable;
394 G4double theRate; 1279 G4double theRate;
395 G4double TaoPlus; 1280 G4double TaoPlus;
396 G4int nS = 0; // Running index of fir 1281 G4int nS = 0; // Running index of first decay in a given generation
397 G4int nT = nEntry; // Total number of deca 1282 G4int nT = nEntry; // Total number of decays accumulated over entire history
398 const G4int nMode = G4RadioactiveDecayModeSi << 1283 const G4int nMode = 12;
399 G4double brs[nMode]; 1284 G4double brs[nMode];
>> 1285
400 // 1286 //
401 theIonTable = G4ParticleTable::GetParticleTa << 1287 theIonTable =
>> 1288 (G4IonTable*)(G4ParticleTable::GetParticleTable()->GetIonTable());
402 1289
403 G4int loop = 0; 1290 G4int loop = 0;
>> 1291
404 while (!stable) { /* Loop checking, 01.09. 1292 while (!stable) { /* Loop checking, 01.09.2015, D.Wright */
405 loop++; 1293 loop++;
406 if (loop > 10000) { 1294 if (loop > 10000) {
407 G4Exception("G4RadioactiveDecay::Calcula 1295 G4Exception("G4RadioactiveDecay::CalculateChainsFromParent()", "HAD_RDM_100",
408 JustWarning, "While loop cou 1296 JustWarning, "While loop count exceeded");
409 break; 1297 break;
410 } 1298 }
411 nGeneration++; 1299 nGeneration++;
412 for (j = nS; j < nT; ++j) { << 1300 for (j = nS; j < nT; j++) {
413 // First time through, get data for pare 1301 // First time through, get data for parent nuclide
414 ZP = theDecayRateVector[j].GetZ(); 1302 ZP = theDecayRateVector[j].GetZ();
415 AP = theDecayRateVector[j].GetA(); 1303 AP = theDecayRateVector[j].GetA();
416 EP = theDecayRateVector[j].GetE(); 1304 EP = theDecayRateVector[j].GetE();
417 RP = theDecayRateVector[j].GetDecayRateC 1305 RP = theDecayRateVector[j].GetDecayRateC();
418 TP = theDecayRateVector[j].GetTaos(); 1306 TP = theDecayRateVector[j].GetTaos();
419 if (GetVerboseLevel() > 1) { << 1307 if (GetVerboseLevel() > 0) {
420 G4cout << "G4RadioactiveDecay::Calcula 1308 G4cout << "G4RadioactiveDecay::CalculateChainsFromParent: daughters of ("
421 << ZP << ", " << AP << ", " << 1309 << ZP << ", " << AP << ", " << EP
422 << ") are being calculated, ge 1310 << ") are being calculated, generation = " << nGeneration
423 << G4endl; 1311 << G4endl;
424 } 1312 }
>> 1313 // G4cout << " Taus = " << G4endl;
>> 1314 // for (G4int ii = 0; ii < TP.size(); ii++) G4cout << TP[ii] << ", " ;
>> 1315 // G4cout << G4endl;
425 1316
426 aParentNucleus = theIonTable->GetIon(ZP, 1317 aParentNucleus = theIonTable->GetIon(ZP,AP,EP);
427 parentDecayTable = GetDecayTable(aParent 1318 parentDecayTable = GetDecayTable(aParentNucleus);
428 if (nullptr == parentDecayTable) { conti <<
429 1319
430 G4DecayTable* summedDecayTable = new G4D 1320 G4DecayTable* summedDecayTable = new G4DecayTable();
431 // This instance of G4DecayTable is for 1321 // This instance of G4DecayTable is for accumulating BRs and decay
432 // channels. It will contain one decay 1322 // channels. It will contain one decay channel per type of decay
433 // (alpha, beta, etc.); its branching ra 1323 // (alpha, beta, etc.); its branching ratio will be the sum of all
434 // branching ratios for that type of dec 1324 // branching ratios for that type of decay of the parent. If the
435 // halflife of a particular channel is l 1325 // halflife of a particular channel is longer than some threshold,
436 // that channel will be inserted specifi 1326 // that channel will be inserted specifically and its branching
437 // ratio will not be included in the abo 1327 // ratio will not be included in the above sums.
438 // This instance is not used to perform 1328 // This instance is not used to perform actual decays.
439 1329
440 for (G4int k = 0; k < nMode; ++k) brs[k] << 1330 for (G4int k = 0; k < nMode; k++) brs[k] = 0.0;
441 1331
442 // Go through the decay table and sum al 1332 // Go through the decay table and sum all channels having the same decay mode
443 for (G4int i = 0; i < parentDecayTable-> << 1333 for (i = 0; i < parentDecayTable->entries(); i++) {
444 theChannel = parentDecayTable->GetDeca 1334 theChannel = parentDecayTable->GetDecayChannel(i);
445 theNuclearDecayChannel = static_cast<G 1335 theNuclearDecayChannel = static_cast<G4NuclearDecay*>(theChannel);
446 theDecayMode = theNuclearDecayChannel- 1336 theDecayMode = theNuclearDecayChannel->GetDecayMode();
447 daughterExcitation = theNuclearDecayCh 1337 daughterExcitation = theNuclearDecayChannel->GetDaughterExcitation();
448 theDaughterNucleus = theNuclearDecayCh << 1338 theDaughterNucleus = theNuclearDecayChannel->GetDaughterNucleus() ;
>> 1339
449 AD = ((const G4Ions*)(theDaughterNucle 1340 AD = ((const G4Ions*)(theDaughterNucleus))->GetAtomicMass();
450 ZD = ((const G4Ions*)(theDaughterNucle 1341 ZD = ((const G4Ions*)(theDaughterNucleus))->GetAtomicNumber();
451 const G4LevelManager* levelManager = 1342 const G4LevelManager* levelManager =
452 G4NuclearLevelData::GetInstance()->G 1343 G4NuclearLevelData::GetInstance()->GetLevelManager(ZD,AD);
453 1344
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( 1345 if (levelManager->NumberOfTransitions() ) {
458 nearestEnergy = levelManager->Neares 1346 nearestEnergy = levelManager->NearestLevelEnergy(daughterExcitation);
459 if ((std::abs(daughterExcitation - n << 1347 if (std::abs(daughterExcitation - nearestEnergy) < levelTolerance) {
460 (std::abs(daughterExcitation - nearest <<
461 // Level half-life is in ns and th 1348 // Level half-life is in ns and the threshold is set to 1 micros
462 // by default, user can set it via 1349 // by default, user can set it via the UI command
463 nearestLevelIndex = (G4int)levelMa << 1350 nearestLevelIndex = levelManager->NearestLevelIndex(daughterExcitation);
464 if (levelManager->LifeTime(nearest 1351 if (levelManager->LifeTime(nearestLevelIndex)*ns >= halflifethreshold){
465 // save the metastable decay cha << 1352 // save the metastable nucleus
466 summedDecayTable->Insert(theChan 1353 summedDecayTable->Insert(theChannel);
467 } else { 1354 } else {
468 brs[theDecayMode] += theChannel- 1355 brs[theDecayMode] += theChannel->GetBR();
469 } 1356 }
470 } else { 1357 } else {
471 brs[theDecayMode] += theChannel->G 1358 brs[theDecayMode] += theChannel->GetBR();
472 } 1359 }
473 } else { 1360 } else {
474 brs[theDecayMode] += theChannel->Get 1361 brs[theDecayMode] += theChannel->GetBR();
475 } 1362 }
476 <<
477 } // Combine decay channels (loop i) 1363 } // Combine decay channels (loop i)
478 << 1364
479 brs[BetaPlus] = brs[BetaPlus]+brs[Kshell << 1365 brs[2] = brs[2]+brs[3]+brs[4]+brs[5]+brs[6]; // Combine beta+ and EC
480 brs[KshellEC] = brs[LshellEC] = brs[Mshe << 1366 brs[3] = brs[4] = brs[5] = brs[6] = 0.0;
481 for (G4int i = 0; i < nMode; ++i) { << 1367 for (i = 0; i < nMode; i++) { // loop over decay modes
482 if (brs[i] > 0.) { 1368 if (brs[i] > 0.) {
483 switch (i) { 1369 switch (i) {
484 case IT: << 1370 case 0:
485 // Decay mode is isomeric transiti 1371 // Decay mode is isomeric transition
486 theITChannel = new G4ITDecay(aPare << 1372 theITChannel = new G4ITDecay(aParentNucleus, brs[0], 0.0, 0.0,
>> 1373 photonEvaporation);
487 1374
488 summedDecayTable->Insert(theITChan 1375 summedDecayTable->Insert(theITChannel);
489 break; 1376 break;
490 1377
491 case BetaMinus: << 1378 case 1:
492 // Decay mode is beta- 1379 // Decay mode is beta-
493 theBetaMinusChannel = new G4BetaMi << 1380 theBetaMinusChannel = new G4BetaMinusDecay(aParentNucleus, brs[1],
494 1381 0.*MeV, 0.*MeV,
495 1382 noFloat, allowed);
496 summedDecayTable->Insert(theBetaMi 1383 summedDecayTable->Insert(theBetaMinusChannel);
497 break; 1384 break;
498 1385
499 case BetaPlus: << 1386 case 2:
500 // Decay mode is beta+ + EC. 1387 // Decay mode is beta+ + EC.
501 theBetaPlusChannel = new G4BetaPlu << 1388 theBetaPlusChannel = new G4BetaPlusDecay(aParentNucleus, brs[2], // DHW: April 2015
502 1389 0.*MeV, 0.*MeV,
503 1390 noFloat, allowed);
504 summedDecayTable->Insert(theBetaPl 1391 summedDecayTable->Insert(theBetaPlusChannel);
505 break; 1392 break;
506 1393
507 case Alpha: << 1394 case 7:
508 // Decay mode is alpha. 1395 // Decay mode is alpha.
509 theAlphaChannel = new G4AlphaDecay << 1396 theAlphaChannel = new G4AlphaDecay(aParentNucleus, brs[7], 0.*MeV,
510 1397 0.*MeV, noFloat);
511 summedDecayTable->Insert(theAlphaC 1398 summedDecayTable->Insert(theAlphaChannel);
512 break; 1399 break;
513 1400
514 case Proton: << 1401 case 8:
515 // Decay mode is proton. 1402 // Decay mode is proton.
516 theProtonChannel = new G4ProtonDec << 1403 theProtonChannel = new G4ProtonDecay(aParentNucleus, brs[8], 0.*MeV,
517 1404 0.*MeV, noFloat);
518 summedDecayTable->Insert(theProton 1405 summedDecayTable->Insert(theProtonChannel);
519 break; 1406 break;
520 1407
521 case Neutron: << 1408 case 9:
522 // Decay mode is neutron. 1409 // Decay mode is neutron.
523 theNeutronChannel = new G4NeutronD << 1410 theNeutronChannel = new G4NeutronDecay(aParentNucleus, brs[9], 0.*MeV,
524 1411 0.*MeV, noFloat);
525 summedDecayTable->Insert(theNeutro 1412 summedDecayTable->Insert(theNeutronChannel);
526 break; 1413 break;
527 1414
528 case SpFission: << 1415 case 10:
529 // Decay mode is spontaneous fissi 1416 // Decay mode is spontaneous fission
530 theFissionChannel = new G4SFDecay( << 1417 theFissionChannel = new G4SFDecay(aParentNucleus, brs[10], 0.*MeV,
531 1418 0.*MeV, noFloat);
532 summedDecayTable->Insert(theFissio 1419 summedDecayTable->Insert(theFissionChannel);
533 break; 1420 break;
534 << 1421 case 11:
535 case BDProton: << 1422 // Decay mode is triton.
536 // Not yet implemented << 1423 theTritonChannel = new G4TritonDecay(aParentNucleus, brs[11], 0.*MeV,
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 1424 0.*MeV, noFloat);
563 summedDecayTable->Insert(theTriton 1425 summedDecayTable->Insert(theTritonChannel);
564 break; 1426 break;
565 << 1427
566 default: 1428 default:
567 break; 1429 break;
568 } 1430 }
569 } 1431 }
570 } 1432 }
571 <<
572 // loop over all branches in summedDecay 1433 // loop over all branches in summedDecayTable
573 // 1434 //
574 for (G4int i = 0; i < summedDecayTable-> << 1435 for (i = 0; i < summedDecayTable->entries(); i++){
575 theChannel = summedDecayTable->GetDeca 1436 theChannel = summedDecayTable->GetDecayChannel(i);
576 theNuclearDecayChannel = static_cast<G 1437 theNuclearDecayChannel = static_cast<G4NuclearDecay*>(theChannel);
577 theBR = theChannel->GetBR(); 1438 theBR = theChannel->GetBR();
578 theDaughterNucleus = theNuclearDecayCh 1439 theDaughterNucleus = theNuclearDecayChannel->GetDaughterNucleus();
579 1440
580 // First check if the decay of the ori 1441 // First check if the decay of the original nucleus is an IT channel,
581 // if true create a new ground-state n 1442 // if true create a new ground-state nucleus
582 if (theNuclearDecayChannel->GetDecayMo 1443 if (theNuclearDecayChannel->GetDecayMode() == IT && nGeneration == 1) {
583 1444
584 A = ((const G4Ions*)(theDaughterNucl 1445 A = ((const G4Ions*)(theDaughterNucleus))->GetAtomicMass();
585 Z = ((const G4Ions*)(theDaughterNucl 1446 Z = ((const G4Ions*)(theDaughterNucleus))->GetAtomicNumber();
586 theDaughterNucleus=theIonTable->GetI 1447 theDaughterNucleus=theIonTable->GetIon(Z,A,0.);
587 } 1448 }
588 if (IsApplicable(*theDaughterNucleus) << 1449 if (IsApplicable(*theDaughterNucleus) && theBR &&
589 aParentNucleus != theDaughterNucle 1450 aParentNucleus != theDaughterNucleus) {
590 // need to make sure daughter has de 1451 // need to make sure daughter has decay table
591 parentDecayTable = GetDecayTable(the 1452 parentDecayTable = GetDecayTable(theDaughterNucleus);
592 if (nullptr != parentDecayTable && p << 1453
>> 1454 if (parentDecayTable->entries() ) {
593 A = ((const G4Ions*)(theDaughterNu 1455 A = ((const G4Ions*)(theDaughterNucleus))->GetAtomicMass();
594 Z = ((const G4Ions*)(theDaughterNu 1456 Z = ((const G4Ions*)(theDaughterNucleus))->GetAtomicNumber();
595 E = ((const G4Ions*)(theDaughterNu 1457 E = ((const G4Ions*)(theDaughterNucleus))->GetExcitationEnergy();
596 1458
597 TaoPlus = theDaughterNucleus->GetP 1459 TaoPlus = theDaughterNucleus->GetPDGLifeTime();
598 if (TaoPlus <= 0.) TaoPlus = 1e-1 1460 if (TaoPlus <= 0.) TaoPlus = 1e-100;
599 1461
600 // first set the taos, one simply 1462 // first set the taos, one simply need to add to the parent ones
601 taos.clear(); 1463 taos.clear();
602 taos = TP; // load lifetimes of 1464 taos = TP; // load lifetimes of all previous generations
603 std::size_t k; << 1465 size_t k;
604 //check that TaoPlus differs from 1466 //check that TaoPlus differs from other taos from at least 1.e5 relative difference
605 //for (k = 0; k < TP.size(); ++k){ << 1467 //for (k = 0; k < TP.size(); k++){
606 //if (std::abs((TaoPlus-TP[k])/TP[ 1468 //if (std::abs((TaoPlus-TP[k])/TP[k])<1.e-5 ) TaoPlus=1.00001*TP[k];
607 //} 1469 //}
608 taos.push_back(TaoPlus); // add d 1470 taos.push_back(TaoPlus); // add daughter lifetime to list
609 // now calculate the coefficiencie 1471 // now calculate the coefficiencies
610 // 1472 //
611 // they are in two parts, first th 1473 // they are in two parts, first the less than n ones
612 // Eq 4.24 of the TN 1474 // Eq 4.24 of the TN
613 Acoeffs.clear(); 1475 Acoeffs.clear();
614 long double ta1,ta2; 1476 long double ta1,ta2;
615 ta2 = (long double)TaoPlus; 1477 ta2 = (long double)TaoPlus;
616 for (k = 0; k < RP.size(); ++k){ << 1478 for (k = 0; k < RP.size(); k++){
617 ta1 = (long double)TP[k]; // 1479 ta1 = (long double)TP[k]; // loop over lifetimes of all previous generations
618 if (ta1 == ta2) { 1480 if (ta1 == ta2) {
619 theRate = 1.e100; 1481 theRate = 1.e100;
620 } else { 1482 } else {
621 theRate = ta1/(ta1-ta2); 1483 theRate = ta1/(ta1-ta2);
622 } 1484 }
623 theRate = theRate * theBR * RP[k 1485 theRate = theRate * theBR * RP[k];
624 Acoeffs.push_back(theRate); 1486 Acoeffs.push_back(theRate);
625 } 1487 }
626 1488
627 // the second part: the n:n coeffi 1489 // the second part: the n:n coefficiency
628 // Eq 4.25 of the TN. Note Yn+1 i 1490 // Eq 4.25 of the TN. Note Yn+1 is zero apart from Y1 which is -1
629 // as treated at line 1013 1491 // as treated at line 1013
630 theRate = 0.; 1492 theRate = 0.;
631 long double aRate, aRate1; 1493 long double aRate, aRate1;
632 aRate1 = 0.L; 1494 aRate1 = 0.L;
633 for (k = 0; k < RP.size(); ++k){ << 1495 for (k = 0; k < RP.size(); k++){
634 ta1 = (long double)TP[k]; 1496 ta1 = (long double)TP[k];
635 if (ta1 == ta2 ) { 1497 if (ta1 == ta2 ) {
636 aRate = 1.e100; 1498 aRate = 1.e100;
637 } else { 1499 } else {
638 aRate = ta2/(ta1-ta2); 1500 aRate = ta2/(ta1-ta2);
639 } 1501 }
640 aRate = aRate * (long double)(th 1502 aRate = aRate * (long double)(theBR * RP[k]);
641 aRate1 += aRate; 1503 aRate1 += aRate;
642 } 1504 }
643 theRate = -aRate1; 1505 theRate = -aRate1;
644 Acoeffs.push_back(theRate); 1506 Acoeffs.push_back(theRate);
645 SetDecayRate (Z,A,E,nGeneration,Ac 1507 SetDecayRate (Z,A,E,nGeneration,Acoeffs,taos);
646 theDecayRateVector.push_back(rates 1508 theDecayRateVector.push_back(ratesToDaughter);
647 nEntry++; 1509 nEntry++;
648 } // there are entries in the table 1510 } // there are entries in the table
649 } // nuclide is OK to decay 1511 } // nuclide is OK to decay
650 } // end of loop (i) over decay table br 1512 } // end of loop (i) over decay table branches
651 1513
652 delete summedDecayTable; 1514 delete summedDecayTable;
653 1515
654 } // Getting contents of decay rate vector 1516 } // Getting contents of decay rate vector (end loop on j)
655 nS = nT; 1517 nS = nT;
656 nT = nEntry; 1518 nT = nEntry;
657 if (nS == nT) stable = true; 1519 if (nS == nT) stable = true;
658 } // while nuclide is not stable 1520 } // while nuclide is not stable
659 1521
660 // end of while loop 1522 // end of while loop
661 // the calculation completed here 1523 // the calculation completed here
662 1524
663 1525
664 // fill the first part of the decay rate tab 1526 // fill the first part of the decay rate table
665 // which is the name of the original particl 1527 // which is the name of the original particle (isotope)
666 chainsFromParent.SetIonName(theParentNucleus << 1528 chainsFromParent.SetIonName(theParentNucleus.GetParticleName());
667 1529
668 // now fill the decay table with the newly c 1530 // now fill the decay table with the newly completed decay rate vector
669 chainsFromParent.SetItsRates(theDecayRateVec 1531 chainsFromParent.SetItsRates(theDecayRateVector);
670 1532
671 // finally add the decayratetable to the tab 1533 // finally add the decayratetable to the tablevector
672 theParentChainTable.push_back(chainsFromPare 1534 theParentChainTable.push_back(chainsFromParent);
673 } 1535 }
674 1536
675 ////////////////////////////////////////////// 1537 ////////////////////////////////////////////////////////////////////////////////
676 // 1538 // //
677 // SetSourceTimeProfile 1539 // SetSourceTimeProfile //
678 // read in the source time profile functio 1540 // read in the source time profile function (histogram) //
679 // 1541 // //
680 ////////////////////////////////////////////// 1542 ////////////////////////////////////////////////////////////////////////////////
681 1543
682 void G4RadioactiveDecay::SetSourceTimeProfile( << 1544 void G4RadioactiveDecay::SetSourceTimeProfile(G4String filename)
683 { 1545 {
684 std::ifstream infile ( filename, std::ios::i 1546 std::ifstream infile ( filename, std::ios::in );
685 if (!infile) { 1547 if (!infile) {
686 G4ExceptionDescription ed; 1548 G4ExceptionDescription ed;
687 ed << " Could not open file " << filename 1549 ed << " Could not open file " << filename << G4endl;
688 G4Exception("G4RadioactiveDecay::SetSource 1550 G4Exception("G4RadioactiveDecay::SetSourceTimeProfile()", "HAD_RDM_001",
689 FatalException, ed); 1551 FatalException, ed);
690 } 1552 }
691 1553
692 G4double bin, flux; 1554 G4double bin, flux;
693 NSourceBin = -1; 1555 NSourceBin = -1;
694 1556
695 G4int loop = 0; 1557 G4int loop = 0;
696 while (infile >> bin >> flux) { /* Loop che 1558 while (infile >> bin >> flux) { /* Loop checking, 01.09.2015, D.Wright */
697 loop++; 1559 loop++;
698 if (loop > 10000) { 1560 if (loop > 10000) {
699 G4Exception("G4RadioactiveDecay::SetSour 1561 G4Exception("G4RadioactiveDecay::SetSourceTimeProfile()", "HAD_RDM_100",
700 JustWarning, "While loop cou 1562 JustWarning, "While loop count exceeded");
701 break; 1563 break;
702 } 1564 }
703 1565
704 NSourceBin++; 1566 NSourceBin++;
705 if (NSourceBin > 99) { 1567 if (NSourceBin > 99) {
706 G4Exception("G4RadioactiveDecay::SetSour 1568 G4Exception("G4RadioactiveDecay::SetSourceTimeProfile()", "HAD_RDM_002",
707 FatalException, "Input sourc 1569 FatalException, "Input source time file too big (>100 rows)");
708 1570
709 } else { 1571 } else {
710 SBin[NSourceBin] = bin * s; // Conver << 1572 SBin[NSourceBin] = bin * s;
711 SProfile[NSourceBin] = flux; // Dimens << 1573 SProfile[NSourceBin] = flux;
712 } 1574 }
713 } 1575 }
714 << 1576 SetAnalogueMonteCarlo(0);
715 AnalogueMC = false; <<
716 infile.close(); 1577 infile.close();
717 1578
718 #ifdef G4VERBOSE 1579 #ifdef G4VERBOSE
719 if (GetVerboseLevel() > 2) << 1580 if (GetVerboseLevel() > 1)
720 G4cout <<" Source Timeprofile Nbin = " << 1581 G4cout <<" Source Timeprofile Nbin = " << NSourceBin <<G4endl;
721 #endif 1582 #endif
722 } 1583 }
723 1584
724 ////////////////////////////////////////////// 1585 ////////////////////////////////////////////////////////////////////////////////
725 // 1586 // //
726 // SetDecayBiasProfile 1587 // SetDecayBiasProfile //
727 // read in the decay bias scheme function ( 1588 // read in the decay bias scheme function (histogram) //
728 // 1589 // //
729 ////////////////////////////////////////////// 1590 ////////////////////////////////////////////////////////////////////////////////
730 1591
731 void G4RadioactiveDecay::SetDecayBias(const G4 << 1592 void G4RadioactiveDecay::SetDecayBias(G4String filename)
732 { 1593 {
>> 1594
733 std::ifstream infile(filename, std::ios::in) 1595 std::ifstream infile(filename, std::ios::in);
734 if (!infile) G4Exception("G4RadioactiveDecay << 1596 if (!infile) G4Exception("G4RadioactiveDecay::SetDecayBias()", "HAD_RDM_003",
735 FatalException, "Un 1597 FatalException, "Unable to open bias data file" );
736 1598
737 G4double bin, flux; 1599 G4double bin, flux;
738 G4int dWindows = 0; 1600 G4int dWindows = 0;
739 G4int i ; 1601 G4int i ;
740 1602
741 theRadioactivityTables.clear(); 1603 theRadioactivityTables.clear();
742 1604
743 NDecayBin = -1; 1605 NDecayBin = -1;
744 1606
745 G4int loop = 0; 1607 G4int loop = 0;
746 while (infile >> bin >> flux ) { /* Loop ch 1608 while (infile >> bin >> flux ) { /* Loop checking, 01.09.2015, D.Wright */
747 NDecayBin++; 1609 NDecayBin++;
748 loop++; 1610 loop++;
749 if (loop > 10000) { 1611 if (loop > 10000) {
750 G4Exception("G4RadioactiveDecay::SetDeca 1612 G4Exception("G4RadioactiveDecay::SetDecayBias()", "HAD_RDM_100",
751 JustWarning, "While loop cou 1613 JustWarning, "While loop count exceeded");
752 break; 1614 break;
753 } 1615 }
754 1616
755 if (NDecayBin > 99) { 1617 if (NDecayBin > 99) {
756 G4Exception("G4RadioactiveDecay::SetDeca << 1618 G4Exception("G4RadioactiveDecay::SetDecayBias()", "HAD_RDM_004",
757 FatalException, "Input bias 1619 FatalException, "Input bias file too big (>100 rows)" );
758 } else { 1620 } else {
759 DBin[NDecayBin] = bin * s; // Conver << 1621 DBin[NDecayBin] = bin * s;
760 DProfile[NDecayBin] = flux; // Dimens << 1622 DProfile[NDecayBin] = flux;
761 if (flux > 0.) { 1623 if (flux > 0.) {
762 decayWindows[NDecayBin] = dWindows; 1624 decayWindows[NDecayBin] = dWindows;
763 dWindows++; 1625 dWindows++;
764 G4RadioactivityTable *rTable = new G4R 1626 G4RadioactivityTable *rTable = new G4RadioactivityTable() ;
765 theRadioactivityTables.push_back(rTabl 1627 theRadioactivityTables.push_back(rTable);
766 } 1628 }
767 } 1629 }
768 } 1630 }
769 for ( i = 1; i<= NDecayBin; ++i) DProfile[i] << 1631 for ( i = 1; i<= NDecayBin; i++) DProfile[i] += DProfile[i-1];
770 for ( i = 0; i<= NDecayBin; ++i) DProfile[i] << 1632 for ( i = 0; i<= NDecayBin; i++) DProfile[i] /= DProfile[NDecayBin];
771 // Normalize so e <<
772 // converted to accumulated probabilities 1633 // converted to accumulated probabilities
773 1634
774 AnalogueMC = false; << 1635 SetAnalogueMonteCarlo(0);
775 infile.close(); 1636 infile.close();
776 1637
777 #ifdef G4VERBOSE 1638 #ifdef G4VERBOSE
778 if (GetVerboseLevel() > 2) << 1639 if (GetVerboseLevel() > 1)
779 G4cout <<" Decay Bias Profile Nbin = " << 1640 G4cout <<" Decay Bias Profile Nbin = " << NDecayBin <<G4endl;
780 #endif 1641 #endif
781 } 1642 }
782 1643
783 ////////////////////////////////////////////// 1644 ////////////////////////////////////////////////////////////////////////////////
784 // 1645 // //
785 // DecayIt 1646 // DecayIt //
786 // 1647 // //
787 ////////////////////////////////////////////// 1648 ////////////////////////////////////////////////////////////////////////////////
788 1649
789 G4VParticleChange* 1650 G4VParticleChange*
790 G4RadioactiveDecay::DecayIt(const G4Track& the 1651 G4RadioactiveDecay::DecayIt(const G4Track& theTrack, const G4Step&)
791 { 1652 {
792 // Initialize G4ParticleChange object, get p 1653 // Initialize G4ParticleChange object, get particle details and decay table
>> 1654
793 fParticleChangeForRadDecay.Initialize(theTra 1655 fParticleChangeForRadDecay.Initialize(theTrack);
794 fParticleChangeForRadDecay.ProposeWeight(the 1656 fParticleChangeForRadDecay.ProposeWeight(theTrack.GetWeight());
795 const G4DynamicParticle* theParticle = theTr 1657 const G4DynamicParticle* theParticle = theTrack.GetDynamicParticle();
796 const G4ParticleDefinition* theParticleDef = 1658 const G4ParticleDefinition* theParticleDef = theParticle->GetDefinition();
797 1659
798 // First check whether RDM applies to the cu 1660 // First check whether RDM applies to the current logical volume
799 if (!isAllVolumesMode) { 1661 if (!isAllVolumesMode) {
800 if (!std::binary_search(ValidVolumes.begin 1662 if (!std::binary_search(ValidVolumes.begin(), ValidVolumes.end(),
801 theTrack.GetVolume()->Get 1663 theTrack.GetVolume()->GetLogicalVolume()->GetName())) {
802 #ifdef G4VERBOSE 1664 #ifdef G4VERBOSE
803 if (GetVerboseLevel()>0) { 1665 if (GetVerboseLevel()>0) {
804 G4cout <<"G4RadioactiveDecay::DecayIt 1666 G4cout <<"G4RadioactiveDecay::DecayIt : "
805 << theTrack.GetVolume()->GetLog 1667 << theTrack.GetVolume()->GetLogicalVolume()->GetName()
806 << " is not selected for the RD 1668 << " is not selected for the RDM"<< G4endl;
807 G4cout << " There are " << ValidVolume 1669 G4cout << " There are " << ValidVolumes.size() << " volumes" << G4endl;
808 G4cout << " The Valid volumes are " << 1670 G4cout << " The Valid volumes are " << G4endl;
809 for (std::size_t i = 0; i< ValidVolume << 1671 for (size_t i = 0; i< ValidVolumes.size(); i++)
810 G4cout << Va 1672 G4cout << ValidVolumes[i] << G4endl;
811 } 1673 }
812 #endif 1674 #endif
813 fParticleChangeForRadDecay.SetNumberOfSe 1675 fParticleChangeForRadDecay.SetNumberOfSecondaries(0);
814 1676
815 // Kill the parent particle. 1677 // Kill the parent particle.
816 fParticleChangeForRadDecay.ProposeTrackS 1678 fParticleChangeForRadDecay.ProposeTrackStatus(fStopAndKill) ;
817 fParticleChangeForRadDecay.ProposeLocalE 1679 fParticleChangeForRadDecay.ProposeLocalEnergyDeposit(0.0);
818 ClearNumberOfInteractionLengthLeft(); 1680 ClearNumberOfInteractionLengthLeft();
819 return &fParticleChangeForRadDecay; 1681 return &fParticleChangeForRadDecay;
820 } 1682 }
821 } 1683 }
822 1684
823 // Now check if particle is valid for RDM 1685 // Now check if particle is valid for RDM
824 if (!(IsApplicable(*theParticleDef) ) ) { 1686 if (!(IsApplicable(*theParticleDef) ) ) {
825 // Particle is not an ion or is outside th 1687 // Particle is not an ion or is outside the nucleuslimits for decay
>> 1688
826 #ifdef G4VERBOSE 1689 #ifdef G4VERBOSE
827 if (GetVerboseLevel() > 1) { << 1690 if (GetVerboseLevel()>0) {
828 G4cout << "G4RadioactiveDecay::DecayIt : << 1691 G4cerr << "G4RadioactiveDecay::DecayIt : "
829 << theParticleDef->GetParticleNam << 1692 << theParticleDef->GetParticleName()
830 << " is not an ion or is outside << 1693 << " is not a valid nucleus for the RDM"<< G4endl;
831 << " Set particle change accordin <<
832 << G4endl; <<
833 } 1694 }
834 #endif 1695 #endif
835 fParticleChangeForRadDecay.SetNumberOfSeco 1696 fParticleChangeForRadDecay.SetNumberOfSecondaries(0);
836 1697
837 // Kill the parent particle 1698 // Kill the parent particle
838 fParticleChangeForRadDecay.ProposeTrackSta 1699 fParticleChangeForRadDecay.ProposeTrackStatus(fStopAndKill) ;
839 fParticleChangeForRadDecay.ProposeLocalEne 1700 fParticleChangeForRadDecay.ProposeLocalEnergyDeposit(0.0);
840 ClearNumberOfInteractionLengthLeft(); 1701 ClearNumberOfInteractionLengthLeft();
841 return &fParticleChangeForRadDecay; 1702 return &fParticleChangeForRadDecay;
842 } 1703 }
843 <<
844 G4DecayTable* theDecayTable = GetDecayTable( 1704 G4DecayTable* theDecayTable = GetDecayTable(theParticleDef);
845 1705
846 if (theDecayTable == nullptr || theDecayTabl << 1706 if (theDecayTable == 0 || theDecayTable->entries() == 0) {
847 // No data in the decay table. Set partic 1707 // No data in the decay table. Set particle change parameters
848 // to indicate this. 1708 // to indicate this.
849 #ifdef G4VERBOSE 1709 #ifdef G4VERBOSE
850 if (GetVerboseLevel() > 1) { << 1710 if (GetVerboseLevel() > 0) {
851 G4cout << "G4RadioactiveDecay::DecayIt : << 1711 G4cerr <<"G4RadioactiveDecay::DecayIt : decay table not defined for ";
852 << "decay table not defined for " << 1712 G4cerr <<theParticleDef->GetParticleName() <<G4endl;
853 << theParticleDef->GetParticleNam <<
854 << ". Set particle change accordi <<
855 << G4endl; <<
856 } 1713 }
857 #endif 1714 #endif
858 fParticleChangeForRadDecay.SetNumberOfSeco 1715 fParticleChangeForRadDecay.SetNumberOfSecondaries(0);
859 1716
860 // Kill the parent particle. 1717 // Kill the parent particle.
861 fParticleChangeForRadDecay.ProposeTrackSta 1718 fParticleChangeForRadDecay.ProposeTrackStatus(fStopAndKill) ;
862 fParticleChangeForRadDecay.ProposeLocalEne 1719 fParticleChangeForRadDecay.ProposeLocalEnergyDeposit(0.0);
863 ClearNumberOfInteractionLengthLeft(); 1720 ClearNumberOfInteractionLengthLeft();
864 return &fParticleChangeForRadDecay; 1721 return &fParticleChangeForRadDecay;
865 1722
866 } else { 1723 } else {
867 // Data found. Try to decay nucleus 1724 // Data found. Try to decay nucleus
>> 1725 G4double energyDeposit = 0.0;
>> 1726 G4double finalGlobalTime = theTrack.GetGlobalTime();
>> 1727 G4double finalLocalTime = theTrack.GetLocalTime();
>> 1728 G4int index;
>> 1729 G4ThreeVector currentPosition;
>> 1730 currentPosition = theTrack.GetPosition();
>> 1731
>> 1732 // Check whether use Analogue or VR implementation
868 if (AnalogueMC) { 1733 if (AnalogueMC) {
869 G4VRadioactiveDecay::DecayAnalog(theTrac << 1734 #ifdef G4VERBOSE
>> 1735 if (GetVerboseLevel() > 0)
>> 1736 G4cout <<"DecayIt: Analogue MC version " << G4endl;
>> 1737 # endif
>> 1738
>> 1739 G4DecayProducts* products = DoDecay(*theParticleDef);
>> 1740
>> 1741 // Check if the product is the same as input and kill the track if
>> 1742 // necessary to prevent infinite loop (11/05/10, F.Lei)
>> 1743 if ( products->entries() == 1) {
>> 1744 fParticleChangeForRadDecay.SetNumberOfSecondaries(0);
>> 1745 fParticleChangeForRadDecay.ProposeTrackStatus(fStopAndKill);
>> 1746 fParticleChangeForRadDecay.ProposeLocalEnergyDeposit(0.0);
>> 1747 ClearNumberOfInteractionLengthLeft();
>> 1748 return &fParticleChangeForRadDecay;
>> 1749 }
870 1750
871 } else { << 1751 // Get parent particle information and boost the decay products to the
872 // Proceed with decay using variance red << 1752 // laboratory frame based on this information.
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 1753
880 G4IonTable* theIonTable; << 1754 //The Parent Energy used for the boost should be the total energy of
881 G4ParticleDefinition* parentNucleus; << 1755 // the nucleus of the parent ion without the energy of the shell electrons
>> 1756 // (correction for bug 1359 by L. Desorgher)
>> 1757 G4double ParentEnergy = theParticle->GetKineticEnergy()
>> 1758 + theParticle->GetParticleDefinition()->GetPDGMass();
>> 1759 G4ThreeVector ParentDirection(theParticle->GetMomentumDirection());
>> 1760
>> 1761 if (theTrack.GetTrackStatus() == fStopButAlive) {
>> 1762 //this condition seems to be always True, further investigation is needed (L.Desorgher)
>> 1763
>> 1764 // The particle is decayed at rest.
>> 1765 // since the time is still for rest particle in G4 we need to add the
>> 1766 // additional time lapsed between the particle come to rest and the
>> 1767 // actual decay. This time is simply sampled with the mean-life of
>> 1768 // the particle. But we need to protect the case PDGTime < 0.
>> 1769 // (F.Lei 11/05/10)
>> 1770 G4double temptime = -std::log( G4UniformRand())
>> 1771 *theParticleDef->GetPDGLifeTime();
>> 1772 if (temptime < 0.) temptime = 0.;
>> 1773 finalGlobalTime += temptime;
>> 1774 finalLocalTime += temptime;
>> 1775 energyDeposit += theParticle->GetKineticEnergy();
>> 1776 }
>> 1777 products->Boost(ParentEnergy, ParentDirection);
882 1778
883 // Get decay chains for the given nuclid << 1779 // Add products in theParticleChangeForRadDecay.
884 if (!IsRateTableReady(*theParticleDef)) << 1780 G4int numberOfSecondaries = products->entries();
885 CalculateChainsFromParent(*theParticleDef); << 1781 fParticleChangeForRadDecay.SetNumberOfSecondaries(numberOfSecondaries);
>> 1782 #ifdef G4VERBOSE
>> 1783 if (GetVerboseLevel()>1) {
>> 1784 G4cout <<"G4RadioactiveDecay::DecayIt : Decay vertex :";
>> 1785 G4cout <<" Time: " <<finalGlobalTime/ns <<"[ns]";
>> 1786 G4cout <<" X:" <<(theTrack.GetPosition()).x() /cm <<"[cm]";
>> 1787 G4cout <<" Y:" <<(theTrack.GetPosition()).y() /cm <<"[cm]";
>> 1788 G4cout <<" Z:" <<(theTrack.GetPosition()).z() /cm <<"[cm]";
>> 1789 G4cout << G4endl;
>> 1790 G4cout <<"G4Decay::DecayIt : decay products in Lab. Frame" <<G4endl;
>> 1791 products->DumpInfo();
>> 1792 products->IsChecked();
>> 1793 }
>> 1794 #endif
>> 1795 for (index=0; index < numberOfSecondaries; index++) {
>> 1796 G4Track* secondary = new G4Track(products->PopProducts(),
>> 1797 finalGlobalTime, currentPosition);
>> 1798
>> 1799 secondary->SetCreatorModelIndex(theRadDecayMode);
>> 1800 //Change for atomics relaxation
>> 1801 if (theRadDecayMode == IT && index>0){
>> 1802 if (index == numberOfSecondaries-1) secondary->SetCreatorModelIndex(IT);
>> 1803 else secondary->SetCreatorModelIndex(30);
>> 1804 }
>> 1805 else if (theRadDecayMode >= KshellEC && theRadDecayMode <= NshellEC
>> 1806 && index <numberOfSecondaries-1){
>> 1807 secondary->SetCreatorModelIndex(30);
>> 1808 }
>> 1809 secondary->SetGoodForTrackingFlag();
>> 1810 secondary->SetTouchableHandle(theTrack.GetTouchableHandle());
>> 1811 fParticleChangeForRadDecay.AddSecondary(secondary);
>> 1812 }
>> 1813 delete products;
>> 1814 // end of analogue MC algorithm
>> 1815
>> 1816 } else {
>> 1817 // Variance Reduction Method
>> 1818 #ifdef G4VERBOSE
>> 1819 if (GetVerboseLevel()>0)
>> 1820 G4cout << "DecayIt: Variance Reduction version " << G4endl;
>> 1821 #endif
>> 1822 // Get decay chains for the given nuclide
>> 1823 if (!IsRateTableReady(*theParticleDef)) CalculateChainsFromParent(*theParticleDef);
886 GetChainsFromParent(*theParticleDef); 1824 GetChainsFromParent(*theParticleDef);
887 1825
888 // Declare some of the variables require << 1826 // declare some of the variables required in the implementation
>> 1827 G4ParticleDefinition* parentNucleus;
>> 1828 G4IonTable* theIonTable;
889 G4int PZ; 1829 G4int PZ;
890 G4int PA; 1830 G4int PA;
891 G4double PE; 1831 G4double PE;
892 G4String keyName; 1832 G4String keyName;
893 std::vector<G4double> PT; 1833 std::vector<G4double> PT;
894 std::vector<G4double> PR; 1834 std::vector<G4double> PR;
895 G4double taotime; << 1835 G4double tauprob;
896 long double decayRate; 1836 long double decayRate;
897 1837
898 std::size_t i; << 1838 size_t i;
>> 1839 // size_t j;
899 G4int numberOfSecondaries; 1840 G4int numberOfSecondaries;
900 G4int totalNumberOfSecondaries = 0; 1841 G4int totalNumberOfSecondaries = 0;
901 G4double currentTime = 0.; 1842 G4double currentTime = 0.;
902 G4int ndecaych; 1843 G4int ndecaych;
903 G4DynamicParticle* asecondaryparticle; 1844 G4DynamicParticle* asecondaryparticle;
904 std::vector<G4DynamicParticle*> secondar 1845 std::vector<G4DynamicParticle*> secondaryparticles;
905 std::vector<G4double> pw; 1846 std::vector<G4double> pw;
906 std::vector<G4double> ptime; 1847 std::vector<G4double> ptime;
907 pw.clear(); 1848 pw.clear();
908 ptime.clear(); 1849 ptime.clear();
909 1850
910 // Now apply the nucleus splitting << 1851 //now apply the nucleus splitting
911 for (G4int n = 0; n < NSplit; ++n) { << 1852 for (G4int n = 0; n < NSplit; n++) {
912 // Get the decay time following the de 1853 // Get the decay time following the decay probability function
913 // supplied by user << 1854 // suppllied by user
914 G4double theDecayTime = GetDecayTime() 1855 G4double theDecayTime = GetDecayTime();
915 G4int nbin = GetDecayTimeBin(theDecayT 1856 G4int nbin = GetDecayTimeBin(theDecayTime);
916 1857
917 // calculate the first part of the wei 1858 // calculate the first part of the weight function
918 G4double weight1 = 1.; 1859 G4double weight1 = 1.;
919 if (nbin == 1) { 1860 if (nbin == 1) {
920 weight1 = 1./DProfile[nbin-1] 1861 weight1 = 1./DProfile[nbin-1]
921 *(DBin[nbin]-DBin[nbin-1]) << 1862 *(DBin[nbin]-DBin[nbin-1])/NSplit;
922 } else if (nbin > 1) { 1863 } else if (nbin > 1) {
923 // Go from nbin to nbin-2 because fl <<
924 weight1 = 1./(DProfile[nbin]-DProfil 1864 weight1 = 1./(DProfile[nbin]-DProfile[nbin-2])
925 *(DBin[nbin]-DBin[nbin-1]) 1865 *(DBin[nbin]-DBin[nbin-1])/NSplit;
926 // weight1 = (probability of choosin <<
927 } 1866 }
>> 1867
928 // it should be calculated in seconds 1868 // it should be calculated in seconds
929 weight1 /= s ; 1869 weight1 /= s ;
930 1870
931 // loop over all the possible secondar << 1871 // Loop over all the possible secondaries of the nucleus
932 // the first one is itself. 1872 // the first one is itself.
933 for (i = 0; i < theDecayRateVector.siz << 1873 for (i = 0; i < theDecayRateVector.size(); i++) {
934 PZ = theDecayRateVector[i].GetZ(); 1874 PZ = theDecayRateVector[i].GetZ();
935 PA = theDecayRateVector[i].GetA(); 1875 PA = theDecayRateVector[i].GetA();
936 PE = theDecayRateVector[i].GetE(); 1876 PE = theDecayRateVector[i].GetE();
937 PT = theDecayRateVector[i].GetTaos() 1877 PT = theDecayRateVector[i].GetTaos();
938 PR = theDecayRateVector[i].GetDecayR 1878 PR = theDecayRateVector[i].GetDecayRateC();
939 1879
940 // The array of arrays theDecayRateV 1880 // The array of arrays theDecayRateVector contains all possible decay
941 // chains of a given parent nucleus 1881 // chains of a given parent nucleus (ZP,AP,EP) to a given descendant
942 // nuclide (Z,A,E). 1882 // nuclide (Z,A,E).
943 // 1883 //
944 // theDecayRateVector[0] contains th 1884 // theDecayRateVector[0] contains the decay parameters of the parent
945 // nucleus 1885 // nucleus
946 // PZ = ZP 1886 // PZ = ZP
947 // PA = AP 1887 // PA = AP
948 // PE = EP 1888 // PE = EP
949 // PT[] = {TP} 1889 // PT[] = {TP}
950 // PR[] = {RP} 1890 // PR[] = {RP}
951 // 1891 //
952 // theDecayRateVector[1] contains th 1892 // theDecayRateVector[1] contains the decay of the parent to the first
953 // generation daughter (Z1,A1,E1). 1893 // generation daughter (Z1,A1,E1).
954 // PZ = Z1 1894 // PZ = Z1
955 // PA = A1 1895 // PA = A1
956 // PE = E1 1896 // PE = E1
957 // PT[] = {TP, T1} 1897 // PT[] = {TP, T1}
958 // PR[] = {RP, R1} 1898 // PR[] = {RP, R1}
959 // 1899 //
960 // theDecayRateVector[2] contains th 1900 // theDecayRateVector[2] contains the decay of the parent to the first
961 // generation daughter (Z1,A1,E1) an 1901 // generation daughter (Z1,A1,E1) and the decay of the first
962 // generation daughter to the second << 1902 // generation daughter to the second generation daughter (Z2,A2,E2).
963 // PZ = Z2 1903 // PZ = Z2
964 // PA = A2 1904 // PA = A2
965 // PE = E2 1905 // PE = E2
966 // PT[] = {TP, T1, T2} 1906 // PT[] = {TP, T1, T2}
967 // PR[] = {RP, R1, R2} 1907 // PR[] = {RP, R1, R2}
968 // 1908 //
969 // theDecayRateVector[3] may contain 1909 // theDecayRateVector[3] may contain a branch chain
970 // PZ = Z2a 1910 // PZ = Z2a
971 // PA = A2a 1911 // PA = A2a
972 // PE = E2a 1912 // PE = E2a
973 // PT[] = {TP, T1, T2a} 1913 // PT[] = {TP, T1, T2a}
974 // PR[] = {RP, R1, R2a} 1914 // PR[] = {RP, R1, R2a}
975 // 1915 //
976 // and so on. 1916 // and so on.
977 1917
978 // Calculate the decay rate of the i << 1918 // Calculate the decay rate of the isotope. decayRate is the
979 // radioactivity of isotope (PZ,PA,P << 1919 // radioactivity of isotope (PZ,PA,PE) at 'theDecayTime'.
980 // it will be used to calculate the << 1920 // It will be used to calculate the statistical weight of the
981 // decay products of this isotope 1921 // decay products of this isotope
982 1922
983 // For each nuclide, calculate all t 1923 // For each nuclide, calculate all the decay chains which can reach
984 // the parent nuclide 1924 // the parent nuclide
985 decayRate = 0.L; 1925 decayRate = 0.L;
986 for (G4int j = 0; j < G4int(PT.size( << 1926 for (G4int j = 0; j < G4int(PT.size()); j++) {
987 taotime = ConvolveSourceTimeProfil << 1927 tauprob = ConvolveSourceTimeProfile(theDecayTime,PT[j]);
988 decayRate -= PR[j] * (long double) << 1928 // tauprob is dimensionless, PR has units of s-1
>> 1929 decayRate -= PR[j] * (long double)tauprob;
989 // Eq.4.23 of of the TN 1930 // Eq.4.23 of of the TN
990 // note the negative here is requi 1931 // note the negative here is required as the rate in the
991 // equation is defined to be negat 1932 // equation is defined to be negative,
992 // i.e. decay away, but we need po 1933 // i.e. decay away, but we need positive value here.
993 1934
994 // G4cout << j << "\t"<< PT[j]/s < << 1935 // G4cout << j << "\t" << PT[j]/s << "\t" << PR[j] << "\t" << decayRate << G4endl;
995 } 1936 }
996 1937
997 // At this point any negative decay 1938 // At this point any negative decay rates are probably small enough
998 // (order 10**-30) that negative val 1939 // (order 10**-30) that negative values are likely due to cancellation
999 // errors. Set them to zero. 1940 // errors. Set them to zero.
1000 if (decayRate < 0.0) decayRate = 0. 1941 if (decayRate < 0.0) decayRate = 0.0;
>> 1942 /*
>> 1943 if (decayRate < 0.0) {
>> 1944 if (-decayRate > 1.0e-30) {
>> 1945 G4ExceptionDescription ed;
>> 1946 ed << " Negative decay probability (magnitude > 1e-30) \n"
>> 1947 << " in variance reduction branch " << G4endl;
>> 1948 G4Exception("G4RadioactiveDecay::DecayIt()",
>> 1949 "HAD_RDM_200", JustWarning, ed);
>> 1950 } else {
>> 1951 // Decay probability is small enough that negative value is likely
>> 1952 // due to cancellation errors. Set it to zero.
>> 1953 decayRate = 0.0;
>> 1954 }
>> 1955 }
1001 1956
1002 // G4cout <<theDecayTime/s <<"\t"< << 1957 if (decayRate < 0.0) G4cout << " NEGATIVE decay rate = " << decayRate << G4endl;
1003 // G4cout << theTrack.GetWeight() << 1958 */
>> 1959 // G4cout << theDecayTime/s << "\t" << nbin << G4endl;
>> 1960 // G4cout << theTrack.GetWeight() << "\t" << weight1 << "\t" << decayRate << G4endl;
1004 1961
1005 // Add isotope to the radioactivity 1962 // Add isotope to the radioactivity tables
1006 // One table for each observation t << 1963 // One table for each observation time window specifed in
1007 // SetDecayBias(G4String filename) 1964 // SetDecayBias(G4String filename)
1008 <<
1009 theRadioactivityTables[decayWindows 1965 theRadioactivityTables[decayWindows[nbin-1]]
1010 ->AddIsotope(PZ,PA,PE,weight1 1966 ->AddIsotope(PZ,PA,PE,weight1*decayRate,theTrack.GetWeight());
1011 1967
1012 // Now calculate the statistical we 1968 // Now calculate the statistical weight
1013 // One needs to fold the source bia 1969 // One needs to fold the source bias function with the decaytime
1014 // also need to include the track w 1970 // also need to include the track weight! (F.Lei, 28/10/10)
1015 G4double weight = weight1*decayRate 1971 G4double weight = weight1*decayRate*theTrack.GetWeight();
1016 1972
1017 // decay the isotope << 1973 // Decay the isotope
1018 theIonTable = (G4IonTable *)(G4Part 1974 theIonTable = (G4IonTable *)(G4ParticleTable::GetParticleTable()->GetIonTable());
1019 parentNucleus = theIonTable->GetIon 1975 parentNucleus = theIonTable->GetIon(PZ,PA,PE);
1020 1976
1021 // Create a temprary products buffe 1977 // Create a temprary products buffer.
1022 // Its contents to be transfered to 1978 // Its contents to be transfered to the products at the end of the loop
1023 G4DecayProducts* tempprods = nullpt << 1979 G4DecayProducts* tempprods = 0;
1024 1980
1025 // Decide whether to apply branchin << 1981 // Decide whether to apply branching ratio bias or not
1026 if (BRBias) { 1982 if (BRBias) {
1027 G4DecayTable* decayTable = GetDec 1983 G4DecayTable* decayTable = GetDecayTable(parentNucleus);
1028 G4VDecayChannel* theDecayChannel = null <<
1029 if (nullptr != decayTable) { <<
1030 ndecaych = G4int(decayTable->entries( <<
1031 theDecayChannel = decayTable->G <<
1032 } <<
1033 1984
1034 if (theDecayChannel == nullptr) { << 1985 ndecaych = G4int(decayTable->entries()*G4UniformRand());
>> 1986 G4VDecayChannel* theDecayChannel = decayTable->GetDecayChannel(ndecaych);
>> 1987 if (theDecayChannel == 0) {
1035 // Decay channel not found. 1988 // Decay channel not found.
1036 << 1989 #ifdef G4VERBOSE
1037 if (GetVerboseLevel() > 0) { << 1990 if (GetVerboseLevel()>0) {
1038 G4cout << " G4RadioactiveDeca << 1991 G4cerr << " G4RadioactiveDecay::DoIt : cannot determine decay channel ";
1039 G4cout << " for this nucleus; << 1992 G4cerr << " for this nucleus; decay as if no biasing active ";
1040 G4cout << G4endl; << 1993 G4cerr << G4endl;
1041 if (nullptr != decayTable) { << 1994 decayTable ->DumpInfo();
1042 } 1995 }
1043 // DHW 6 Dec 2010 - do decay as if no << 1996 #endif
1044 tempprods = DoDecay(*parentNucl << 1997 tempprods = DoDecay(*parentNucleus); // DHW 6 Dec 2010 - do decay as if no biasing
>> 1998 // to avoid deref of temppprods = 0
1045 } else { 1999 } else {
1046 // A decay channel has been ide 2000 // A decay channel has been identified, so execute the DecayIt.
1047 G4double tempmass = parentNucle 2001 G4double tempmass = parentNucleus->GetPDGMass();
1048 tempprods = theDecayChannel->De 2002 tempprods = theDecayChannel->DecayIt(tempmass);
1049 weight *= (theDecayChannel->Get 2003 weight *= (theDecayChannel->GetBR())*(decayTable->entries());
1050 } 2004 }
1051 } else { 2005 } else {
1052 tempprods = DoDecay(*parentNucleu << 2006 tempprods = DoDecay(*parentNucleus);
1053 } 2007 }
1054 2008
1055 // save the secondaries for buffers << 2009 // Save the secondaries for buffers
1056 numberOfSecondaries = tempprods->en 2010 numberOfSecondaries = tempprods->entries();
1057 currentTime = finalGlobalTime + the 2011 currentTime = finalGlobalTime + theDecayTime;
1058 for (index = 0; index < numberOfSec << 2012 for (index = 0; index < numberOfSecondaries; index++) {
1059 asecondaryparticle = tempprods->P 2013 asecondaryparticle = tempprods->PopProducts();
1060 if (asecondaryparticle->GetDefini 2014 if (asecondaryparticle->GetDefinition()->GetPDGStable() ) {
1061 pw.push_back(weight); 2015 pw.push_back(weight);
1062 ptime.push_back(currentTime); 2016 ptime.push_back(currentTime);
1063 secondaryparticles.push_back(as 2017 secondaryparticles.push_back(asecondaryparticle);
1064 } << 2018 } else if (((const G4Ions*)(asecondaryparticle->GetDefinition()))->GetExcitationEnergy() > 0.
1065 // Generate gammas and Xrays from << 2019 && weight > 0.) {
1066 else if (((const G4Ions*)(asecond << 2020 // Compute the gamma
1067 ->GetExcitationEnergy() << 2021 // Generate gammas and XRays from excited nucleus, added by L.Desorgher
1068 G4ParticleDefinition* apartDef << 2022 G4ParticleDefinition* apartDef =asecondaryparticle->GetDefinition();
1069 AddDeexcitationSpectrumForBiasM << 2023 AddDeexcitationSpectrumForBiasMode(apartDef,weight,currentTime,pw,ptime,secondaryparticles);
1070 <<
1071 } 2024 }
1072 } 2025 }
1073 <<
1074 delete tempprods; 2026 delete tempprods;
>> 2027
1075 } // end of i loop 2028 } // end of i loop
1076 } // end of n loop 2029 } // end of n loop
1077 2030
1078 // now deal with the secondaries in the 2031 // now deal with the secondaries in the two stl containers
1079 // and submmit them back to the trackin 2032 // and submmit them back to the tracking manager
1080 totalNumberOfSecondaries = (G4int)pw.si << 2033 totalNumberOfSecondaries = pw.size();
1081 fParticleChangeForRadDecay.SetNumberOfS 2034 fParticleChangeForRadDecay.SetNumberOfSecondaries(totalNumberOfSecondaries);
1082 for (index=0; index < totalNumberOfSeco << 2035 for (index=0; index < totalNumberOfSecondaries; index++) {
1083 G4Track* secondary = new G4Track(seco 2036 G4Track* secondary = new G4Track(secondaryparticles[index],
1084 ptim 2037 ptime[index], currentPosition);
1085 secondary->SetGoodForTrackingFlag(); 2038 secondary->SetGoodForTrackingFlag();
1086 secondary->SetTouchableHandle(theTrac 2039 secondary->SetTouchableHandle(theTrack.GetTouchableHandle());
1087 secondary->SetWeight(pw[index]); 2040 secondary->SetWeight(pw[index]);
1088 fParticleChangeForRadDecay.AddSeconda 2041 fParticleChangeForRadDecay.AddSecondary(secondary);
1089 } 2042 }
>> 2043 // make sure the original track is set to stop and its kinematic energy collected
>> 2044 //
>> 2045 //theTrack.SetTrackStatus(fStopButAlive);
>> 2046 //energyDeposit += theParticle->GetKineticEnergy();
1090 2047
1091 // Kill the parent particle << 2048 } // End of Variance Reduction
1092 fParticleChangeForRadDecay.ProposeTrack << 2049
1093 fParticleChangeForRadDecay.ProposeLocal << 2050 // Kill the parent particle
1094 fParticleChangeForRadDecay.ProposeLocal << 2051 fParticleChangeForRadDecay.ProposeTrackStatus(fStopAndKill) ;
1095 // Reset NumberOfInteractionLengthLeft. << 2052 fParticleChangeForRadDecay.ProposeLocalEnergyDeposit(energyDeposit);
1096 ClearNumberOfInteractionLengthLeft(); << 2053 fParticleChangeForRadDecay.ProposeLocalTime(finalLocalTime);
1097 } // end VR decay << 2054 // Reset NumberOfInteractionLengthLeft.
>> 2055 ClearNumberOfInteractionLengthLeft();
1098 2056
1099 return &fParticleChangeForRadDecay; 2057 return &fParticleChangeForRadDecay;
1100 } // end of data found branch << 2058 }
1101 } 2059 }
1102 2060
1103 2061
>> 2062 G4DecayProducts*
>> 2063 G4RadioactiveDecay::DoDecay(const G4ParticleDefinition& theParticleDef)
>> 2064 {
>> 2065 G4DecayProducts* products = 0;
>> 2066 G4DecayTable* theDecayTable = GetDecayTable(&theParticleDef);
>> 2067
>> 2068 // Choose a decay channel.
>> 2069 #ifdef G4VERBOSE
>> 2070 if (GetVerboseLevel() > 0) G4cout << "Select a channel..." << G4endl;
>> 2071 #endif
>> 2072
>> 2073 // G4DecayTable::SelectADecayChannel checks to see if sum of daughter masses
>> 2074 // exceeds parent mass. Pass it the parent mass + maximum Q value to account
>> 2075 // for difference in mass defect.
>> 2076 G4double parentPlusQ = theParticleDef.GetPDGMass() + 30.*MeV;
>> 2077 G4VDecayChannel* theDecayChannel = theDecayTable->SelectADecayChannel(parentPlusQ);
>> 2078 theRadDecayMode = (static_cast<G4NuclearDecay*>(theDecayChannel))->GetDecayMode();
>> 2079 if (theDecayChannel == 0) {
>> 2080 // Decay channel not found.
>> 2081 G4ExceptionDescription ed;
>> 2082 ed << " Cannot determine decay channel for " << theParticleDef.GetParticleName() << G4endl;
>> 2083 G4Exception("G4RadioactiveDecay::DoDecay", "HAD_RDM_013",
>> 2084 FatalException, ed);
>> 2085 } else {
>> 2086 // A decay channel has been identified, so execute the DecayIt.
>> 2087 #ifdef G4VERBOSE
>> 2088 if (GetVerboseLevel() > 1) {
>> 2089 G4cerr << "G4RadioactiveDecay::DoIt : selected decay channel addr:";
>> 2090 G4cerr << theDecayChannel << G4endl;
>> 2091 }
>> 2092 #endif
>> 2093 products = theDecayChannel->DecayIt(theParticleDef.GetPDGMass() );
>> 2094
>> 2095 // Apply directional bias if requested by user
>> 2096 CollimateDecay(products);
>> 2097 }
>> 2098
>> 2099 return products;
>> 2100 }
>> 2101
>> 2102
>> 2103 // Apply directional bias for "visible" daughters (e+-, gamma, n, p, alpha)
>> 2104
>> 2105 void G4RadioactiveDecay::CollimateDecay(G4DecayProducts* products) {
>> 2106 if (origin == forceDecayDirection) return; // No collimation requested
>> 2107 if (180.*deg == forceDecayHalfAngle) return;
>> 2108 if (0 == products || 0 == products->entries()) return;
>> 2109
>> 2110 #ifdef G4VERBOSE
>> 2111 if (GetVerboseLevel() > 0) G4cout << "Begin of CollimateDecay..." << G4endl;
>> 2112 #endif
>> 2113
>> 2114 // Particles suitable for directional biasing (for if-blocks below)
>> 2115 static const G4ParticleDefinition* electron = G4Electron::Definition();
>> 2116 static const G4ParticleDefinition* positron = G4Positron::Definition();
>> 2117 static const G4ParticleDefinition* neutron = G4Neutron::Definition();
>> 2118 static const G4ParticleDefinition* gamma = G4Gamma::Definition();
>> 2119 static const G4ParticleDefinition* alpha = G4Alpha::Definition();
>> 2120 static const G4ParticleDefinition* triton = G4Triton::Definition();
>> 2121 static const G4ParticleDefinition* proton = G4Proton::Definition();
>> 2122
>> 2123 G4ThreeVector newDirection; // Re-use to avoid memory churn
>> 2124 for (G4int i=0; i<products->entries(); i++) {
>> 2125 G4DynamicParticle* daughter = (*products)[i];
>> 2126 const G4ParticleDefinition* daughterType =
>> 2127 daughter->GetParticleDefinition();
>> 2128 if (daughterType == electron || daughterType == positron ||
>> 2129 daughterType == neutron || daughterType == gamma ||
>> 2130 daughterType == alpha || daughterType == triton || daughterType == proton) CollimateDecayProduct(daughter);
>> 2131 }
>> 2132 }
>> 2133
>> 2134 void G4RadioactiveDecay::CollimateDecayProduct(G4DynamicParticle* daughter) {
>> 2135 #ifdef G4VERBOSE
>> 2136 if (GetVerboseLevel() > 1) {
>> 2137 G4cout << "CollimateDecayProduct for daughter "
>> 2138 << daughter->GetParticleDefinition()->GetParticleName() << G4endl;
>> 2139 }
>> 2140 #endif
>> 2141
>> 2142 G4ThreeVector collimate = ChooseCollimationDirection();
>> 2143 if (origin != collimate) daughter->SetMomentumDirection(collimate);
>> 2144 }
>> 2145
>> 2146
>> 2147 // Choose random direction within collimation cone
>> 2148
>> 2149 G4ThreeVector G4RadioactiveDecay::ChooseCollimationDirection() const {
>> 2150 if (origin == forceDecayDirection) return origin; // Don't do collimation
>> 2151 if (forceDecayHalfAngle == 180.*deg) return origin;
>> 2152
>> 2153 G4ThreeVector dir = forceDecayDirection;
>> 2154
>> 2155 // Return direction offset by random throw
>> 2156 if (forceDecayHalfAngle > 0.) {
>> 2157 // Generate uniform direction around central axis
>> 2158 G4double phi = 2.*pi*G4UniformRand();
>> 2159 G4double cosMin = std::cos(forceDecayHalfAngle);
>> 2160 G4double cosTheta = (1.-cosMin)*G4UniformRand() + cosMin; // [cosMin,1.)
>> 2161
>> 2162 dir.setPhi(dir.phi()+phi);
>> 2163 dir.setTheta(dir.theta()+std::acos(cosTheta));
>> 2164 }
>> 2165
>> 2166 #ifdef G4VERBOSE
>> 2167 if (GetVerboseLevel()>1)
>> 2168 G4cout << " ChooseCollimationDirection returns " << dir << G4endl;
>> 2169 #endif
>> 2170
>> 2171 return dir;
>> 2172 }
>> 2173
1104 // Add gamma, X-ray, conversion and auger ele 2174 // Add gamma, X-ray, conversion and auger electrons for bias mode
1105 void << 2175 void
1106 G4RadioactiveDecay::AddDeexcitationSpectrumFo 2176 G4RadioactiveDecay::AddDeexcitationSpectrumForBiasMode(G4ParticleDefinition* apartDef,
1107 G4dou 2177 G4double weight,G4double currentTime,
1108 std:: 2178 std::vector<double>& weights_v,
1109 std:: 2179 std::vector<double>& times_v,
1110 std:: 2180 std::vector<G4DynamicParticle*>& secondaries_v)
1111 { 2181 {
1112 G4double elevel=((const G4Ions*)(apartDef)) << 2182 G4double elevel = ((const G4Ions*)(apartDef))->GetExcitationEnergy();
1113 G4double life_time=apartDef->GetPDGLifeTime << 2183 G4double life_time = apartDef->GetPDGLifeTime();
1114 G4ITDecay* anITChannel = 0; 2184 G4ITDecay* anITChannel = 0;
1115 2185
1116 while (life_time < halflifethreshold && ele 2186 while (life_time < halflifethreshold && elevel > 0.) {
1117 decayIT->SetupDecay(apartDef); << 2187 anITChannel = new G4ITDecay(apartDef, 100., elevel, elevel, photonEvaporation);
1118 G4DecayProducts* pevap_products = decayIT << 2188 G4DecayProducts* pevap_products = anITChannel->DecayIt(0.);
1119 G4int nb_pevapSecondaries = pevap_product 2189 G4int nb_pevapSecondaries = pevap_products->entries();
1120 2190
1121 G4DynamicParticle* a_pevap_secondary = 0; 2191 G4DynamicParticle* a_pevap_secondary = 0;
1122 G4ParticleDefinition* secDef = 0; 2192 G4ParticleDefinition* secDef = 0;
1123 for (G4int ind = 0; ind < nb_pevapSeconda 2193 for (G4int ind = 0; ind < nb_pevapSecondaries; ind++) {
1124 a_pevap_secondary= pevap_products->PopP << 2194 a_pevap_secondary = pevap_products->PopProducts();
1125 secDef = a_pevap_secondary->GetDefiniti 2195 secDef = a_pevap_secondary->GetDefinition();
1126 2196
1127 if (secDef->GetBaryonNumber() > 4) { 2197 if (secDef->GetBaryonNumber() > 4) {
1128 elevel = ((const G4Ions*)(secDef))->G 2198 elevel = ((const G4Ions*)(secDef))->GetExcitationEnergy();
1129 life_time = secDef->GetPDGLifeTime(); 2199 life_time = secDef->GetPDGLifeTime();
1130 apartDef = secDef; 2200 apartDef = secDef;
1131 if (secDef->GetPDGStable() ) { 2201 if (secDef->GetPDGStable() ) {
1132 weights_v.push_back(weight); 2202 weights_v.push_back(weight);
1133 times_v.push_back(currentTime); 2203 times_v.push_back(currentTime);
1134 secondaries_v.push_back(a_pevap_sec 2204 secondaries_v.push_back(a_pevap_secondary);
1135 } 2205 }
1136 } else { 2206 } else {
1137 weights_v.push_back(weight); 2207 weights_v.push_back(weight);
1138 times_v.push_back(currentTime); 2208 times_v.push_back(currentTime);
1139 secondaries_v.push_back(a_pevap_secon 2209 secondaries_v.push_back(a_pevap_secondary);
1140 } 2210 }
1141 } 2211 }
1142 2212
1143 delete anITChannel; 2213 delete anITChannel;
1144 delete pevap_products; <<
1145 } 2214 }
1146 } 2215 }
>> 2216
1147 2217
1148 2218