| 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 // --------------
>> 46 // 17 Oct 2011, L. Desorgher
>> 47 // -Add possibility for the user to load its own decay file.
>> 48 // -Set halflifethreshold negative by default to allow the tracking of all
>> 49 // excited nuclei resulting from a radioactive decay
34 // 50 //
35 // Based on the code of F. Lei and P.R. Trusc << 51 // 01 June 2011, M. Kelsey -- Add directional biasing interface to allow for
>> 52 // "collimation" of decay daughters.
>> 53 // 16 February 2006, V.Ivanchenko fix problem in IsApplicable connected with
>> 54 // 8.0 particle design
>> 55 // 18 October 2002, F. Lei
>> 56 // in the case of beta decay, added a check of the end-energy
>> 57 // to ensure it is > 0.
>> 58 // ENSDF occationally have beta decay entries with zero energies
>> 59 //
>> 60 // 27 Sepetember 2001, F. Lei
>> 61 // verboselevel(0) used in constructor
>> 62 //
>> 63 // 01 November 2000, F.Lei
>> 64 // added " ee = e0 +1. ;" as line 763
>> 65 // tagged as "radiative_decay-V02-00-02"
>> 66 // 28 October 2000, F Lei
>> 67 // added fast beta decay mode. Many files have been changed.
>> 68 // tagged as "radiative_decay-V02-00-01"
>> 69 //
>> 70 // 25 October 2000, F Lei, DERA UK
>> 71 // 1) line 1185 added 'const' to work with tag "Track-V02-00-00"
>> 72 // tagged as "radiative_decay-V02-00-00"
>> 73 // 14 April 2000, F Lei, DERA UK
>> 74 // 0.b.4 release. Changes are:
>> 75 // 1) Use PhotonEvaporation instead of DiscreteGammaDeexcitation
>> 76 // 2) VR: Significant efficiency improvement
>> 77 //
>> 78 // 29 February 2000, P R Truscott, DERA UK
>> 79 // 0.b.3 release.
>> 80 //
>> 81 ///////////////////////////////////////////////////////////////////////////////
36 // 82 //
37 ////////////////////////////////////////////// <<
38 <<
39 #include "G4RadioactiveDecay.hh" 83 #include "G4RadioactiveDecay.hh"
40 #include "G4RadioactivationMessenger.hh" << 84 #include "G4RadioactiveDecaymessenger.hh"
41 85
42 #include "G4SystemOfUnits.hh" <<
43 #include "G4DynamicParticle.hh" 86 #include "G4DynamicParticle.hh"
44 #include "G4DecayProducts.hh" 87 #include "G4DecayProducts.hh"
45 #include "G4DecayTable.hh" 88 #include "G4DecayTable.hh"
>> 89 #include "G4PhysicsLogVector.hh"
46 #include "G4ParticleChangeForRadDecay.hh" 90 #include "G4ParticleChangeForRadDecay.hh"
47 #include "G4ITDecay.hh" << 91 #include "G4ITDecayChannel.hh"
48 #include "G4BetaDecayType.hh" << 92 #include "G4BetaMinusDecayChannel.hh"
49 #include "G4BetaMinusDecay.hh" << 93 #include "G4BetaPlusDecayChannel.hh"
50 #include "G4BetaPlusDecay.hh" << 94 #include "G4KshellECDecayChannel.hh"
51 #include "G4ECDecay.hh" << 95 #include "G4LshellECDecayChannel.hh"
52 #include "G4AlphaDecay.hh" << 96 #include "G4MshellECDecayChannel.hh"
53 #include "G4TritonDecay.hh" << 97 #include "G4AlphaDecayChannel.hh"
54 #include "G4ProtonDecay.hh" <<
55 #include "G4NeutronDecay.hh" <<
56 #include "G4SFDecay.hh" <<
57 #include "G4VDecayChannel.hh" 98 #include "G4VDecayChannel.hh"
58 #include "G4NuclearDecay.hh" <<
59 #include "G4RadioactiveDecayMode.hh" 99 #include "G4RadioactiveDecayMode.hh"
60 #include "G4Fragment.hh" <<
61 #include "G4Ions.hh" 100 #include "G4Ions.hh"
62 #include "G4IonTable.hh" 101 #include "G4IonTable.hh"
>> 102 #include "G4RIsotopeTable.hh"
63 #include "G4BetaDecayType.hh" 103 #include "G4BetaDecayType.hh"
>> 104 #include "G4BetaDecayCorrections.hh"
64 #include "Randomize.hh" 105 #include "Randomize.hh"
65 #include "G4LogicalVolumeStore.hh" 106 #include "G4LogicalVolumeStore.hh"
66 #include "G4NuclearLevelData.hh" << 107 #include "G4NuclearLevelManager.hh"
67 #include "G4DeexPrecoParameters.hh" << 108 #include "G4NuclearLevelStore.hh"
68 #include "G4LevelManager.hh" <<
69 #include "G4ThreeVector.hh" 109 #include "G4ThreeVector.hh"
70 #include "G4Electron.hh" 110 #include "G4Electron.hh"
71 #include "G4Positron.hh" 111 #include "G4Positron.hh"
72 #include "G4Neutron.hh" 112 #include "G4Neutron.hh"
73 #include "G4Gamma.hh" 113 #include "G4Gamma.hh"
74 #include "G4Alpha.hh" 114 #include "G4Alpha.hh"
75 #include "G4Triton.hh" <<
76 #include "G4Proton.hh" <<
77 115
>> 116 #include "G4HadTmpUtil.hh"
78 #include "G4HadronicProcessType.hh" 117 #include "G4HadronicProcessType.hh"
79 #include "G4HadronicProcessStore.hh" <<
80 #include "G4HadronicException.hh" <<
81 #include "G4LossTableManager.hh" <<
82 #include "G4VAtomDeexcitation.hh" <<
83 #include "G4UAtomicDeexcitation.hh" <<
84 #include "G4PhotonEvaporation.hh" <<
85 118
86 #include <vector> 119 #include <vector>
87 #include <sstream> 120 #include <sstream>
88 #include <algorithm> 121 #include <algorithm>
89 #include <fstream> 122 #include <fstream>
90 123
91 using namespace CLHEP; 124 using namespace CLHEP;
92 125
93 G4RadioactiveDecay::G4RadioactiveDecay(const G << 126 const G4double G4RadioactiveDecay::levelTolerance =2.0*keV;
94 const G << 127 const G4ThreeVector G4RadioactiveDecay::origin(0.,0.,0.);
95 : G4VRadioactiveDecay(processName, timeThres << 128
>> 129 G4RadioactiveDecay::G4RadioactiveDecay(const G4String& processName)
>> 130 : G4VRestDiscreteProcess(processName, fDecay), HighestBinValue(10.0),
>> 131 LowestBinValue(1.0e-3), TotBin(200), forceDecayDirection(0.,0.,0.),
>> 132 forceDecayHalfAngle(0.*deg), verboseLevel(0)
96 { 133 {
97 #ifdef G4VERBOSE 134 #ifdef G4VERBOSE
98 if (GetVerboseLevel() > 1) { << 135 if (GetVerboseLevel()>1) {
99 G4cout << "G4RadioactiveDecay constructor: << 136 G4cout <<"G4RadioactiveDecay constructor Name: ";
100 << G4endl; << 137 G4cout <<processName << G4endl; }
101 } <<
102 #endif 138 #endif
103 139
104 theRadioactivationMessenger = new G4Radioact << 140 SetProcessSubType(fRadioactiveDecay);
105 141
>> 142 theRadioactiveDecaymessenger = new G4RadioactiveDecaymessenger(this);
>> 143 theIsotopeTable = new G4RIsotopeTable();
>> 144 aPhysicsTable = 0;
>> 145 pParticleChange = &fParticleChangeForRadDecay;
>> 146
>> 147 // Now register the Isotope table with G4IonTable.
>> 148 G4IonTable *theIonTable =
>> 149 (G4IonTable *)(G4ParticleTable::GetParticleTable()->GetIonTable());
>> 150 G4VIsotopeTable *aVirtualTable = theIsotopeTable;
>> 151 theIonTable->RegisterIsotopeTable(aVirtualTable);
>> 152
>> 153 // Reset the contents of the list of nuclei for which decay scheme data
>> 154 // have been loaded.
>> 155 LoadedNuclei.clear();
>> 156
>> 157 //
>> 158 //Reset the list of user define data file
>> 159 //
>> 160 theUserRadioactiveDataFiles.clear();
>> 161
>> 162 //
>> 163 //
106 // Apply default values. 164 // Apply default values.
>> 165 //
107 NSourceBin = 1; 166 NSourceBin = 1;
108 SBin[0] = 0.* s; 167 SBin[0] = 0.* s;
109 SBin[1] = 1.* s; // Convert to ns << 168 SBin[1] = 1.* s;
110 SProfile[0] = 1.; 169 SProfile[0] = 1.;
111 SProfile[1] = 0.; 170 SProfile[1] = 0.;
112 NDecayBin = 1; 171 NDecayBin = 1;
113 DBin[0] = 0. * s ; 172 DBin[0] = 0. * s ;
114 DBin[1] = 1. * s; 173 DBin[1] = 1. * s;
115 DProfile[0] = 1.; 174 DProfile[0] = 1.;
116 DProfile[1] = 0.; 175 DProfile[1] = 0.;
117 decayWindows[0] = 0; 176 decayWindows[0] = 0;
118 G4RadioactivityTable* rTable = new G4Radioac 177 G4RadioactivityTable* rTable = new G4RadioactivityTable() ;
119 theRadioactivityTables.push_back(rTable); 178 theRadioactivityTables.push_back(rTable);
120 NSplit = 1; 179 NSplit = 1;
121 AnalogueMC = true; << 180 AnalogueMC = true ;
122 BRBias = true; << 181 FBeta = false ;
123 halflifethreshold = 1000.*nanosecond; << 182 BRBias = true ;
>> 183 applyICM = true ;
>> 184 applyARM = true ;
>> 185 halflifethreshold = -1.*second;
>> 186 //
>> 187 // RDM applies to xall logical volumes as default
>> 188 SelectAllVolumes();
124 } 189 }
125 190
126 void G4RadioactiveDecay::ProcessDescription(st << 191
127 { << 192 G4RadioactiveDecay::~G4RadioactiveDecay()
128 outFile << "The G4RadioactiveDecay process p << 193 {
129 << "nuclides (G4GenericIon) in biase << 194 if (aPhysicsTable != 0) {
130 << "duplication, branching ratio bia << 195 aPhysicsTable->clearAndDestroy();
131 << "and detector time convolution. << 196 delete aPhysicsTable;
132 << "activation physics.\n" << 197 }
133 << "The required half-lives and deca << 198 delete theRadioactiveDecaymessenger;
134 << "the RadioactiveDecay database wh <<
135 } 199 }
136 200
137 201
138 G4RadioactiveDecay::~G4RadioactiveDecay() << 202 G4bool
>> 203 G4RadioactiveDecay::IsApplicable(const G4ParticleDefinition& aParticle)
>> 204 {
>> 205 // All particles, other than G4Ions, are rejected by default
>> 206 if (((const G4Ions*)(&aParticle))->GetExcitationEnergy() > 0.) {return true;}
>> 207 if (aParticle.GetParticleName() == "GenericIon") {
>> 208 return true;
>> 209 } else if (!(aParticle.GetParticleType() == "nucleus")
>> 210 || aParticle.GetPDGLifeTime() < 0. ) {
>> 211 return false;
>> 212 }
>> 213
>> 214 // Determine whether the nuclide falls into the correct A and Z range
>> 215
>> 216 G4int A = ((const G4Ions*) (&aParticle))->GetAtomicMass();
>> 217 G4int Z = ((const G4Ions*) (&aParticle))->GetAtomicNumber();
>> 218 if (A > theNucleusLimits.GetAMax() || A < theNucleusLimits.GetAMin())
>> 219 {return false;}
>> 220 else if (Z > theNucleusLimits.GetZMax() || Z < theNucleusLimits.GetZMin())
>> 221 {return false;}
>> 222 return true;
>> 223 }
>> 224
>> 225
>> 226 G4bool G4RadioactiveDecay::IsLoaded(const G4ParticleDefinition &aParticle)
>> 227 {
>> 228 // Check whether the radioactive decay data on the ion have already been
>> 229 // loaded
>> 230
>> 231 return std::binary_search(LoadedNuclei.begin(),
>> 232 LoadedNuclei.end(),
>> 233 aParticle.GetParticleName());
>> 234 }
>> 235
>> 236
>> 237 void G4RadioactiveDecay::SelectAVolume(const G4String aVolume)
>> 238 {
>> 239 G4LogicalVolumeStore* theLogicalVolumes;
>> 240 G4LogicalVolume *volume;
>> 241 theLogicalVolumes=G4LogicalVolumeStore::GetInstance();
>> 242 for (size_t i = 0; i < theLogicalVolumes->size(); i++){
>> 243 volume=(*theLogicalVolumes)[i];
>> 244 if (volume->GetName() == aVolume) {
>> 245 ValidVolumes.push_back(aVolume);
>> 246 std::sort(ValidVolumes.begin(), ValidVolumes.end());
>> 247 // sort need for performing binary_search
>> 248 #ifdef G4VERBOSE
>> 249 if (GetVerboseLevel()>0)
>> 250 G4cout << " RDM Applies to : " << aVolume << G4endl;
>> 251 #endif
>> 252 }else if(i == theLogicalVolumes->size())
>> 253 {
>> 254 G4cerr << "SelectAVolume: "<<aVolume << " is not a valid logical volume name"<< G4endl;
>> 255 }
>> 256 }
>> 257 }
>> 258
>> 259
>> 260 void G4RadioactiveDecay::DeselectAVolume(const G4String aVolume)
>> 261 {
>> 262 G4LogicalVolumeStore *theLogicalVolumes;
>> 263 G4LogicalVolume *volume;
>> 264 theLogicalVolumes=G4LogicalVolumeStore::GetInstance();
>> 265 for (size_t i = 0; i < theLogicalVolumes->size(); i++){
>> 266 volume=(*theLogicalVolumes)[i];
>> 267 if (volume->GetName() == aVolume) {
>> 268 std::vector<G4String>::iterator location;
>> 269 location = std::find(ValidVolumes.begin(),ValidVolumes.end(),aVolume);
>> 270 if (location != ValidVolumes.end()) {
>> 271 ValidVolumes.erase(location);
>> 272 std::sort(ValidVolumes.begin(), ValidVolumes.end());
>> 273 } else {
>> 274 G4cerr << " DeselectVolume:" << aVolume << " is not in the list"<< G4endl;
>> 275 }
>> 276 #ifdef G4VERBOSE
>> 277 if (GetVerboseLevel()>0)
>> 278 G4cout << " DeselectVolume: " << aVolume << " is removed from list"<<G4endl;
>> 279 #endif
>> 280 } else if (i == theLogicalVolumes->size()) {
>> 281 G4cerr << " DeselectVolume:" << aVolume
>> 282 << "is not a valid logical volume name"<< G4endl;
>> 283 }
>> 284 }
>> 285 }
>> 286
>> 287
>> 288 void G4RadioactiveDecay::SelectAllVolumes()
139 { 289 {
140 delete theRadioactivationMessenger; << 290 G4LogicalVolumeStore *theLogicalVolumes;
>> 291 G4LogicalVolume *volume;
>> 292 theLogicalVolumes=G4LogicalVolumeStore::GetInstance();
>> 293 ValidVolumes.clear();
>> 294 #ifdef G4VERBOSE
>> 295 if (GetVerboseLevel()>0)
>> 296 G4cout << " RDM Applies to all Volumes" << G4endl;
>> 297 #endif
>> 298 for (size_t i = 0; i < theLogicalVolumes->size(); i++){
>> 299 volume=(*theLogicalVolumes)[i];
>> 300 ValidVolumes.push_back(volume->GetName());
>> 301 #ifdef G4VERBOSE
>> 302 if (GetVerboseLevel()>0)
>> 303 G4cout << " RDM Applies to Volume " << volume->GetName() << G4endl;
>> 304 #endif
>> 305 }
>> 306 std::sort(ValidVolumes.begin(), ValidVolumes.end());
>> 307 // sort needed in order to allow binary_search
>> 308 }
>> 309
>> 310
>> 311 void G4RadioactiveDecay::DeselectAllVolumes()
>> 312 {
>> 313 ValidVolumes.clear();
>> 314 #ifdef G4VERBOSE
>> 315 if (GetVerboseLevel()>0)
>> 316 G4cout << " RDM removed from all volumes" << G4endl;
>> 317 #endif
141 } 318 }
142 319
>> 320
143 G4bool 321 G4bool
144 G4RadioactiveDecay::IsRateTableReady(const G4P 322 G4RadioactiveDecay::IsRateTableReady(const G4ParticleDefinition& aParticle)
145 { 323 {
146 // Check whether the radioactive decay rates 324 // Check whether the radioactive decay rates table for the ion has already
147 // been calculated. 325 // been calculated.
148 G4String aParticleName = aParticle.GetPartic 326 G4String aParticleName = aParticle.GetParticleName();
149 for (std::size_t i = 0; i < theParentChainTa << 327 for (size_t i = 0; i < theDecayRateTableVector.size(); i++) {
150 if (theParentChainTable[i].GetIonName() == << 328 if (theDecayRateTableVector[i].GetIonName() == aParticleName)
>> 329 return true;
151 } 330 }
152 return false; 331 return false;
153 } 332 }
154 333
>> 334 // GetDecayRateTable
>> 335 // retrieve the decayratetable for the specified aParticle
>> 336
155 void 337 void
156 G4RadioactiveDecay::GetChainsFromParent(const << 338 G4RadioactiveDecay::GetDecayRateTable(const G4ParticleDefinition& aParticle)
157 { 339 {
158 // Retrieve the decay rate table for the spe <<
159 G4String aParticleName = aParticle.GetPartic 340 G4String aParticleName = aParticle.GetParticleName();
160 341
161 for (std::size_t i = 0; i < theParentChainTa << 342 for (size_t i = 0; i < theDecayRateTableVector.size(); i++) {
162 if (theParentChainTable[i].GetIonName() == << 343 if (theDecayRateTableVector[i].GetIonName() == aParticleName) {
163 theDecayRateVector = theParentChainTable << 344 theDecayRateVector = theDecayRateTableVector[i].GetItsRates();
164 } 345 }
165 } 346 }
166 #ifdef G4VERBOSE 347 #ifdef G4VERBOSE
167 if (GetVerboseLevel() > 1) { << 348 if (GetVerboseLevel() > 0) {
168 G4cout << "The DecayRate Table for " << aP << 349 G4cout << "The DecayRate Table for "
169 << G4endl; << 350 << aParticleName << " is selected." << G4endl;
170 } 351 }
171 #endif 352 #endif
172 } 353 }
173 354
174 // ConvolveSourceTimeProfile performs the conv << 355 // GetTaoTime performs the convolution of the source time profile function
175 // function with a single exponential characte << 356 // with the decay constants in the decay chain.
176 // decay chain. The time profile is treated a <<
177 // convolution integral can be done bin-by-bin <<
178 // This implements Eq. 4.13 of DERA technical <<
179 357
180 G4double << 358 G4double G4RadioactiveDecay::GetTaoTime(G4double t, G4double tao)
181 G4RadioactiveDecay::ConvolveSourceTimeProfile( <<
182 { 359 {
183 G4double convolvedTime = 0.0; << 360 G4double taotime =0.;
184 G4int nbin; 361 G4int nbin;
185 if ( t > SBin[NSourceBin]) { 362 if ( t > SBin[NSourceBin]) {
186 nbin = NSourceBin; << 363 nbin = NSourceBin;}
187 } else { << 364 else {
188 nbin = 0; 365 nbin = 0;
189 << 366 while (t > SBin[nbin]) nbin++;
190 G4int loop = 0; << 367 nbin--;}
191 while (t > SBin[nbin]) { // Loop checking << 368 if (nbin > 0) {
192 loop++; << 369 for (G4int i = 0; i < nbin; i++)
193 if (loop > 1000) { << 370 {
194 G4Exception("G4RadioactiveDecay::Convo << 371 taotime += SProfile[i] * (std::exp(-(t-SBin[i+1])/tao)-std::exp(-(t-SBin[i])/tao));
195 "HAD_RDM_100", JustWarning <<
196 break; <<
197 } 372 }
198 ++nbin; <<
199 } <<
200 --nbin; <<
201 } 373 }
202 << 374 taotime += SProfile[nbin] * (1-std::exp(-(t-SBin[nbin])/tao));
203 // Use expm1 wherever possible to avoid larg << 375 if (taotime < 0.) {
204 // 1 - exp(x) for small x << 376 G4cout <<" Tao time: " <<taotime << " reset to zero!"<<G4endl;
205 G4double earg = 0.0; << 377 taotime = 0.;
206 if (nbin > 0) { <<
207 for (G4int i = 0; i < nbin; ++i) { <<
208 earg = (SBin[i+1] - SBin[i])/tau; <<
209 if (earg < 100.) { <<
210 convolvedTime += SProfile[i] * std::ex <<
211 std::expm1(earg); <<
212 } else { <<
213 convolvedTime += SProfile[i] * <<
214 (std::exp(-(t-SBin[i+1])/tau)-std::e <<
215 } <<
216 } <<
217 } 378 }
218 convolvedTime -= SProfile[nbin] * std::expm1 <<
219 // tau divided out of final result to provid <<
220 379
221 if (convolvedTime < 0.) { <<
222 G4cout << " Convolved time =: " << convolv <<
223 G4cout << " t = " << t << " tau = " << tau <<
224 G4cout << SBin[nbin] << " " << SBin[0] << <<
225 convolvedTime = 0.; <<
226 } <<
227 #ifdef G4VERBOSE 380 #ifdef G4VERBOSE
228 if (GetVerboseLevel() > 2) << 381 if (GetVerboseLevel()>1)
229 G4cout << " Convolved time: " << convolved << 382 {G4cout <<" Tao time: " <<taotime <<G4endl;}
230 #endif 383 #endif
231 return convolvedTime; << 384 return taotime ;
232 } 385 }
233 386
234 <<
235 ////////////////////////////////////////////// 387 ////////////////////////////////////////////////////////////////////////////////
236 // 388 // //
237 // GetDecayTime 389 // GetDecayTime //
238 // Randomly select a decay time for the dec 390 // Randomly select a decay time for the decay process, following the //
239 // supplied decay time bias scheme. 391 // supplied decay time bias scheme. //
240 // 392 // //
241 ////////////////////////////////////////////// 393 ////////////////////////////////////////////////////////////////////////////////
242 394
243 G4double G4RadioactiveDecay::GetDecayTime() 395 G4double G4RadioactiveDecay::GetDecayTime()
244 { 396 {
245 G4double decaytime = 0.; 397 G4double decaytime = 0.;
246 G4double rand = G4UniformRand(); 398 G4double rand = G4UniformRand();
247 G4int i = 0; 399 G4int i = 0;
248 << 400 while ( DProfile[i] < rand) i++;
249 G4int loop = 0; <<
250 while (DProfile[i] < rand) { /* Loop checki <<
251 // Entries in DProfile[i] are all between <<
252 // Comparison with rand chooses which time <<
253 ++i; <<
254 loop++; <<
255 if (loop > 100000) { <<
256 G4Exception("G4RadioactiveDecay::GetDeca <<
257 JustWarning, "While loop cou <<
258 break; <<
259 } <<
260 } <<
261 <<
262 rand = G4UniformRand(); 401 rand = G4UniformRand();
263 decaytime = DBin[i] + rand*(DBin[i+1]-DBin[i 402 decaytime = DBin[i] + rand*(DBin[i+1]-DBin[i]);
264 #ifdef G4VERBOSE 403 #ifdef G4VERBOSE
265 if (GetVerboseLevel() > 2) << 404 if (GetVerboseLevel()>1)
266 G4cout <<" Decay time: " <<decaytime/s <<" << 405 {G4cout <<" Decay time: " <<decaytime/s <<"[s]" <<G4endl;}
267 #endif 406 #endif
268 return decaytime; 407 return decaytime;
269 } 408 }
270 409
271 410
272 G4int G4RadioactiveDecay::GetDecayTimeBin(cons 411 G4int G4RadioactiveDecay::GetDecayTimeBin(const G4double aDecayTime)
273 { 412 {
274 G4int i = 0; 413 G4int i = 0;
275 << 414 while ( aDecayTime > DBin[i] ) i++;
276 G4int loop = 0; <<
277 while (aDecayTime > DBin[i] ) { /* Loop ch <<
278 ++i; <<
279 loop++; <<
280 if (loop > 100000) { <<
281 G4Exception("G4RadioactiveDecay::GetDeca <<
282 JustWarning, "While loop cou <<
283 break; <<
284 } <<
285 } <<
286 <<
287 return i; 415 return i;
288 } 416 }
289 417
290 ////////////////////////////////////////////// 418 ////////////////////////////////////////////////////////////////////////////////
291 // 419 // //
292 // GetMeanLifeTime (required by the base clas 420 // GetMeanLifeTime (required by the base class) //
293 // 421 // //
294 ////////////////////////////////////////////// 422 ////////////////////////////////////////////////////////////////////////////////
295 423
296 G4double G4RadioactiveDecay::GetMeanLifeTime(c 424 G4double G4RadioactiveDecay::GetMeanLifeTime(const G4Track& theTrack,
297 G4 << 425 G4ForceCondition* )
298 { 426 {
299 // For variance reduction time is set to 0 s << 427 // For varience reduction implementation the time is set to 0 so as to
300 // to decay immediately. << 428 // force the particle to decay immediately.
301 // In analogue mode it returns the particle' << 429 // in analogueMC mode it return the particles meanlife.
>> 430 //
302 G4double meanlife = 0.; 431 G4double meanlife = 0.;
303 if (AnalogueMC) meanlife = G4VRadioactiveDec << 432 if (AnalogueMC) {
304 return meanlife; << 433 const G4DynamicParticle* theParticle = theTrack.GetDynamicParticle();
>> 434 G4ParticleDefinition* theParticleDef = theParticle->GetDefinition();
>> 435 G4double theLife = theParticleDef->GetPDGLifeTime();
>> 436
>> 437 #ifdef G4VERBOSE
>> 438 if (GetVerboseLevel()>2)
>> 439 {
>> 440 G4cout <<"G4RadioactiveDecay::GetMeanLifeTime() " <<G4endl;
>> 441 G4cout <<"KineticEnergy:" <<theParticle->GetKineticEnergy()/GeV <<"[GeV]";
>> 442 G4cout <<"Mass:" <<theParticle->GetMass()/GeV <<"[GeV]";
>> 443 G4cout <<"Life time: " <<theLife/ns <<"[ns]" << G4endl;
>> 444 }
>> 445 #endif
>> 446 if (theParticleDef->GetPDGStable()) {meanlife = DBL_MAX;}
>> 447 else if (theLife < 0.0) {meanlife = DBL_MAX;}
>> 448 else {meanlife = theLife;}
>> 449 // set the meanlife to zero for excited istopes which is not in the RDM database
>> 450 if (((const G4Ions*)(theParticleDef))->GetExcitationEnergy() > 0. && meanlife == DBL_MAX) {meanlife = 0.;}
>> 451 }
>> 452 #ifdef G4VERBOSE
>> 453 if (GetVerboseLevel()>1)
>> 454 {G4cout <<"mean life time: " <<meanlife/s <<"[s]" <<G4endl;}
>> 455 #endif
>> 456
>> 457 return meanlife;
>> 458 }
>> 459
>> 460 ////////////////////////////////////////////////////////////////////////////////
>> 461 // //
>> 462 // GetMeanFreePath (similar in function to GetMeanFreeTime) //
>> 463 // //
>> 464 ////////////////////////////////////////////////////////////////////////////////
>> 465
>> 466 G4double G4RadioactiveDecay::GetMeanFreePath (const G4Track& aTrack,
>> 467 G4double, G4ForceCondition*)
>> 468 {
>> 469 // constants
>> 470 G4bool isOutRange ;
>> 471
>> 472 // get particle
>> 473 const G4DynamicParticle* aParticle = aTrack.GetDynamicParticle();
>> 474
>> 475 // returns the mean free path in GEANT4 internal units
>> 476 G4double pathlength;
>> 477 G4ParticleDefinition* aParticleDef = aParticle->GetDefinition();
>> 478 G4double aCtau = c_light * aParticleDef->GetPDGLifeTime();
>> 479 G4double aMass = aParticle->GetMass();
>> 480
>> 481 #ifdef G4VERBOSE
>> 482 if (GetVerboseLevel()>2) {
>> 483 G4cout << "G4RadioactiveDecay::GetMeanFreePath() "<< G4endl;
>> 484 G4cout << "KineticEnergy:" << aParticle->GetKineticEnergy()/GeV <<"[GeV]";
>> 485 G4cout << "Mass:" << aMass/GeV <<"[GeV]";
>> 486 G4cout << "c*Tau:" << aCtau/m <<"[m]" <<G4endl;
>> 487 }
>> 488 #endif
>> 489
>> 490 // check if the particle is stable?
>> 491 if (aParticleDef->GetPDGStable()) {
>> 492 pathlength = DBL_MAX;
>> 493
>> 494 } else if (aCtau < 0.0) {
>> 495 pathlength = DBL_MAX;
>> 496
>> 497 //check if the particle has very short life time ?
>> 498 } else if (aCtau < DBL_MIN) {
>> 499 pathlength = DBL_MIN;
>> 500
>> 501 //check if zero mass
>> 502 } else if (aMass < DBL_MIN) {
>> 503 pathlength = DBL_MAX;
>> 504 #ifdef G4VERBOSE
>> 505 if (GetVerboseLevel()>1) {
>> 506 G4cerr << " Zero Mass particle " << G4endl;
>> 507 }
>> 508 #endif
>> 509 } else {
>> 510 //calculate the mean free path
>> 511 // by using normalized kinetic energy (= Ekin/mass)
>> 512 G4double rKineticEnergy = aParticle->GetKineticEnergy()/aMass;
>> 513 if ( rKineticEnergy > HighestBinValue) {
>> 514 // beta >> 1
>> 515 pathlength = ( rKineticEnergy + 1.0)* aCtau;
>> 516 } else if ( rKineticEnergy > LowestBinValue) {
>> 517 // check if aPhysicsTable exists
>> 518 if (aPhysicsTable == 0) BuildPhysicsTable(*aParticleDef);
>> 519 // beta is in the range valid for PhysicsTable
>> 520 pathlength = aCtau *
>> 521 ((*aPhysicsTable)(0))-> GetValue(rKineticEnergy,isOutRange);
>> 522 } else if ( rKineticEnergy < DBL_MIN ) {
>> 523 // too slow particle
>> 524 #ifdef G4VERBOSE
>> 525 if (GetVerboseLevel()>2) {
>> 526 G4cout << "G4Decay::GetMeanFreePath() !!particle stops!!";
>> 527 G4cout << aParticleDef->GetParticleName() << G4endl;
>> 528 G4cout << "KineticEnergy:" << aParticle->GetKineticEnergy()/GeV <<"[GeV]";
>> 529 }
>> 530 #endif
>> 531 pathlength = DBL_MIN;
>> 532 } else {
>> 533 // beta << 1
>> 534 pathlength = (aParticle->GetTotalMomentum())/aMass*aCtau ;
>> 535 }
>> 536 }
>> 537 #ifdef G4VERBOSE
>> 538 if (GetVerboseLevel()>1) {
>> 539 G4cout << "mean free path: "<< pathlength/m << "[m]" << G4endl;
>> 540 }
>> 541 #endif
>> 542 return pathlength;
>> 543 }
>> 544
>> 545
>> 546 void G4RadioactiveDecay::BuildPhysicsTable(const G4ParticleDefinition&)
>> 547 {
>> 548 // if aPhysicsTableis has already been created, do nothing
>> 549 if (aPhysicsTable != 0) return;
>> 550
>> 551 // create aPhysicsTable
>> 552 if (GetVerboseLevel()>1) G4cerr <<" G4Decay::BuildPhysicsTable() "<< G4endl;
>> 553 aPhysicsTable = new G4PhysicsTable(1);
>> 554
>> 555 //create physics vector
>> 556 G4PhysicsLogVector* aVector = new G4PhysicsLogVector(
>> 557 LowestBinValue,
>> 558 HighestBinValue,
>> 559 TotBin);
>> 560
>> 561 G4double beta, gammainv;
>> 562 // fill physics Vector
>> 563 G4int i;
>> 564 for ( i = 0 ; i < TotBin ; i++ ) {
>> 565 gammainv = 1.0/(aVector->GetLowEdgeEnergy(i) + 1.0);
>> 566 beta = std::sqrt((1.0 - gammainv)*(1.0 +gammainv));
>> 567 aVector->PutValue(i, beta/gammainv);
>> 568 }
>> 569 aPhysicsTable->insert(aVector);
>> 570 }
>> 571
>> 572 ////////////////////////////////////////////////////////////////////////////////
>> 573 // //
>> 574 // LoadDecayTable loads the decay scheme from the RadioactiveDecay database //
>> 575 // for the parent nucleus. //
>> 576 // //
>> 577 ////////////////////////////////////////////////////////////////////////////////
>> 578
>> 579 G4DecayTable*
>> 580 G4RadioactiveDecay::LoadDecayTable(G4ParticleDefinition& theParentNucleus)
>> 581 {
>> 582 // Create and initialise variables used in the method.
>> 583 G4DecayTable* theDecayTable = new G4DecayTable();
>> 584
>> 585 // Generate input data file name using Z and A of the parent nucleus
>> 586 // file containing radioactive decay data.
>> 587 G4int A = ((const G4Ions*)(&theParentNucleus))->GetAtomicMass();
>> 588 G4int Z = ((const G4Ions*)(&theParentNucleus))->GetAtomicNumber();
>> 589 G4double E = ((const G4Ions*)(&theParentNucleus))->GetExcitationEnergy();
>> 590
>> 591
>> 592
>> 593 //Check if data have been provided by the user
>> 594 G4String file= theUserRadioactiveDataFiles[1000*A+Z];
>> 595
>> 596 if (file ==""){
>> 597 if ( !getenv("G4RADIOACTIVEDATA") ) {
>> 598 G4cout << "Please setenv G4RADIOACTIVEDATA to point to the radioactive decay data files." << G4endl;
>> 599 throw G4HadronicException(__FILE__, __LINE__,
>> 600 "Please setenv G4RADIOACTIVEDATA to point to the radioactive decay data files.");
>> 601 }
>> 602 G4String dirName = getenv("G4RADIOACTIVEDATA");
>> 603 LoadedNuclei.push_back(theParentNucleus.GetParticleName());
>> 604 std::sort( LoadedNuclei.begin(), LoadedNuclei.end() );
>> 605 // sort needed to allow binary_search
>> 606
>> 607 std::ostringstream os;
>> 608 os <<dirName <<"/z" <<Z <<".a" <<A <<'\0';
>> 609 file = os.str();
>> 610 }
>> 611
>> 612 std::ifstream DecaySchemeFile(file);
>> 613
>> 614 G4bool found(false);
>> 615 if (DecaySchemeFile) {
>> 616 // Initialise variables used for reading in radioactive decay data.
>> 617 G4int nMode = 7;
>> 618 G4bool modeFirstRecord[7];
>> 619 G4double modeTotalBR[7] = {0.0};
>> 620 G4double modeSumBR[7];
>> 621 for (G4int i = 0; i < nMode; i++) {
>> 622 modeFirstRecord[i] = true;
>> 623 modeSumBR[i] = 0.0;
>> 624 }
>> 625
>> 626 G4bool complete(false);
>> 627 char inputChars[100]={' '};
>> 628 G4String inputLine;
>> 629 G4String recordType("");
>> 630 G4RadioactiveDecayMode theDecayMode;
>> 631 G4double a(0.0);
>> 632 G4double b(0.0);
>> 633 G4double c(0.0);
>> 634 G4BetaDecayType betaType(allowed);
>> 635 G4double e0;
>> 636
>> 637 // Loop through each data file record until you identify the decay
>> 638 // data relating to the nuclide of concern.
>> 639
>> 640 while (!complete && !DecaySchemeFile.getline(inputChars, 100).eof()) {
>> 641 inputLine = inputChars;
>> 642 inputLine = inputLine.strip(1);
>> 643 if (inputChars[0] != '#' && inputLine.length() != 0) {
>> 644 std::istringstream tmpStream(inputLine);
>> 645
>> 646 if (inputChars[0] == 'P') {
>> 647 // Nucleus is a parent type. Check excitation level to see if it
>> 648 // matches that of theParentNucleus
>> 649 tmpStream >> recordType >> a >> b;
>> 650 if (found) {complete = true;}
>> 651 else {found = (std::abs(a*keV - E) < levelTolerance);}
>> 652
>> 653 } else if (found) {
>> 654 // The right part of the radioactive decay data file has been found. Search
>> 655 // through it to determine the mode of decay of the subsequent records.
>> 656 if (inputChars[0] == 'W') {
>> 657 #ifdef G4VERBOSE
>> 658 if (GetVerboseLevel() > 0) {
>> 659 // a comment line identified and print out the message
>> 660 //
>> 661 G4cout << " Warning in G4RadioactiveDecay::LoadDecayTable " << G4endl;
>> 662 G4cout << " In data file " << file << G4endl;
>> 663 G4cout << " " << inputLine << G4endl;
>> 664 }
>> 665 #endif
>> 666 } else {
>> 667 tmpStream >> theDecayMode >> a >> b >> c >> betaType;
>> 668
>> 669 // Allowed transitions are the default. Forbidden transitions are
>> 670 // indicated in the last column.
>> 671 if (inputLine.length() < 80) betaType = allowed;
>> 672 a /= 1000.;
>> 673 c /= 1000.;
>> 674
>> 675 switch (theDecayMode) {
>> 676
>> 677 case IT: // Isomeric transition
>> 678 {
>> 679 G4ITDecayChannel* anITChannel =
>> 680 new G4ITDecayChannel(GetVerboseLevel(),
>> 681 (const G4Ions*)& theParentNucleus, b);
>> 682 anITChannel->SetICM(applyICM);
>> 683 anITChannel->SetARM(applyARM);
>> 684 anITChannel->SetHLThreshold(halflifethreshold);
>> 685 theDecayTable->Insert(anITChannel);
>> 686 }
>> 687 break;
>> 688
>> 689 case BetaMinus:
>> 690 {
>> 691 if (modeFirstRecord[1]) {
>> 692 modeFirstRecord[1] = false;
>> 693 modeTotalBR[1] = b;
>> 694 } else {
>> 695 if (c > 0.) {
>> 696 e0 = c*MeV/0.511;
>> 697 G4BetaDecayCorrections corrections(Z+1, A);
>> 698
>> 699 // array to store function shape
>> 700 G4int npti = 100;
>> 701 G4double* pdf = new G4double[npti];
>> 702
>> 703 G4double e; // Total electron energy in units of electron mass
>> 704 G4double p; // Electron momentum in units of electron mass
>> 705 G4double f; // Spectral shape function value
>> 706 for (G4int ptn = 0; ptn < npti; ptn++) {
>> 707 // Calculate simple phase space spectrum
>> 708 e = 1. + e0*(ptn+0.5)/100.;
>> 709 p = std::sqrt(e*e - 1.);
>> 710 f = p*e*(e0-e+1)*(e0-e+1);
>> 711
>> 712 // Apply Fermi factor to get allowed shape
>> 713 f *= corrections.FermiFunction(e);
>> 714
>> 715 // Apply shape factor for forbidden transitions
>> 716 f *= corrections.ShapeFactor(betaType, p, e0-e+1.);
>> 717 pdf[ptn] = f;
>> 718 }
>> 719
>> 720 RandGeneral* aRandomEnergy = new RandGeneral( pdf, npti);
>> 721 G4BetaMinusDecayChannel *aBetaMinusChannel = new
>> 722 G4BetaMinusDecayChannel(GetVerboseLevel(), &theParentNucleus,
>> 723 b, c*MeV, a*MeV, 0, FBeta, aRandomEnergy);
>> 724 aBetaMinusChannel->SetICM(applyICM);
>> 725 aBetaMinusChannel->SetARM(applyARM);
>> 726 aBetaMinusChannel->SetHLThreshold(halflifethreshold);
>> 727 theDecayTable->Insert(aBetaMinusChannel);
>> 728 modeSumBR[1] += b;
>> 729 delete[] pdf;
>> 730 } // c > 0
>> 731 } // if not first record
>> 732 }
>> 733 break;
>> 734
>> 735 case BetaPlus:
>> 736 {
>> 737 if (modeFirstRecord[2]) {
>> 738 modeFirstRecord[2] = false;
>> 739 modeTotalBR[2] = b;
>> 740 } else {
>> 741 e0 = c*MeV/0.511 - 2.;
>> 742 // Need to test e0 for nuclei which have Q < 2Me in their
>> 743 // data files (e.g. z67.a162)
>> 744 if (e0 > 0.) {
>> 745 G4BetaDecayCorrections corrections(1-Z, A);
>> 746
>> 747 // array to store function shape
>> 748 G4int npti = 100;
>> 749 G4double* pdf = new G4double[npti];
>> 750
>> 751 G4double e; // Total positron energy in units of electron mass
>> 752 G4double p; // Positron momentum in units of electron mass
>> 753 G4double f; // Spectral shape function value
>> 754 for (G4int ptn = 0; ptn < npti; ptn++) {
>> 755 // Calculate simple phase space spectrum
>> 756 e = 1. + e0*(ptn+0.5)/100.;
>> 757 p = std::sqrt(e*e - 1.);
>> 758 f = p*e*(e0-e+1)*(e0-e+1);
>> 759
>> 760 // Apply Fermi factor to get allowed shape
>> 761 f *= corrections.FermiFunction(e);
>> 762
>> 763 // Apply shape factor for forbidden transitions
>> 764 f *= corrections.ShapeFactor(betaType, p, e0-e+1.);
>> 765 pdf[ptn] = f;
>> 766 }
>> 767 RandGeneral* aRandomEnergy = new RandGeneral( pdf, npti);
>> 768 G4BetaPlusDecayChannel *aBetaPlusChannel = new
>> 769 G4BetaPlusDecayChannel(GetVerboseLevel(), &theParentNucleus,
>> 770 b, c*MeV, a*MeV, 0, FBeta, aRandomEnergy);
>> 771 aBetaPlusChannel->SetICM(applyICM);
>> 772 aBetaPlusChannel->SetARM(applyARM);
>> 773 aBetaPlusChannel->SetHLThreshold(halflifethreshold);
>> 774 theDecayTable->Insert(aBetaPlusChannel);
>> 775 modeSumBR[2] += b;
>> 776 delete[] pdf;
>> 777 } // if e0 > 0
>> 778 } // if not first record
>> 779 }
>> 780 break;
>> 781
>> 782 case KshellEC: // K-shell electron capture
>> 783
>> 784 if (modeFirstRecord[3]) {
>> 785 modeFirstRecord[3] = false;
>> 786 modeTotalBR[3] = b;
>> 787 } else {
>> 788 G4KshellECDecayChannel* aKECChannel =
>> 789 new G4KshellECDecayChannel(GetVerboseLevel(),
>> 790 &theParentNucleus,
>> 791 b, c*MeV, a*MeV);
>> 792 aKECChannel->SetICM(applyICM);
>> 793 aKECChannel->SetARM(applyARM);
>> 794 aKECChannel->SetHLThreshold(halflifethreshold);
>> 795 theDecayTable->Insert(aKECChannel);
>> 796 modeSumBR[3] += b;
>> 797 }
>> 798 break;
>> 799
>> 800 case LshellEC: // L-shell electron capture
>> 801
>> 802 if (modeFirstRecord[4]) {
>> 803 modeFirstRecord[4] = false;
>> 804 modeTotalBR[4] = b;
>> 805 } else {
>> 806 G4LshellECDecayChannel *aLECChannel = new
>> 807 G4LshellECDecayChannel (GetVerboseLevel(), &theParentNucleus,
>> 808 b, c*MeV, a*MeV);
>> 809 aLECChannel->SetICM(applyICM);
>> 810 aLECChannel->SetARM(applyARM);
>> 811 aLECChannel->SetHLThreshold(halflifethreshold);
>> 812 theDecayTable->Insert(aLECChannel);
>> 813 modeSumBR[4] += b;
>> 814 }
>> 815 break;
>> 816
>> 817 case MshellEC: // M-shell electron capture
>> 818 // In this implementation it is added to L-shell case
>> 819 if (modeFirstRecord[5]) {
>> 820 modeFirstRecord[5] = false;
>> 821 modeTotalBR[5] = b;
>> 822 } else {
>> 823 G4MshellECDecayChannel* aMECChannel =
>> 824 new G4MshellECDecayChannel(GetVerboseLevel(),
>> 825 &theParentNucleus,
>> 826 b, c*MeV, a*MeV);
>> 827 aMECChannel->SetICM(applyICM);
>> 828 aMECChannel->SetARM(applyARM);
>> 829 aMECChannel->SetHLThreshold(halflifethreshold);
>> 830 theDecayTable->Insert(aMECChannel);
>> 831 modeSumBR[5] += b;
>> 832 }
>> 833 break;
>> 834
>> 835 case Alpha:
>> 836 //G4cout<<"Alpha channel"<<a<<'\t'<<b<<'\t'<<c<<std::endl;
>> 837
>> 838 if (modeFirstRecord[6]) {
>> 839 modeFirstRecord[6] = false;
>> 840 modeTotalBR[6] = b;
>> 841 } else {
>> 842 G4AlphaDecayChannel* anAlphaChannel =
>> 843 new G4AlphaDecayChannel(GetVerboseLevel(),
>> 844 &theParentNucleus,
>> 845 b, c*MeV, a*MeV);
>> 846 anAlphaChannel->SetICM(applyICM);
>> 847 anAlphaChannel->SetARM(applyARM);
>> 848 anAlphaChannel->SetHLThreshold(halflifethreshold);
>> 849 theDecayTable->Insert(anAlphaChannel);
>> 850 modeSumBR[6] += b;
>> 851 }
>> 852 break;
>> 853 case SpFission:
>> 854 //Still needed to be implemented
>> 855 //G4cout<<"Sp fission channel"<<a<<'\t'<<b<<'\t'<<c<<std::endl;
>> 856 break;
>> 857 case ERROR:
>> 858
>> 859 default:
>> 860 G4Exception("G4RadioactiveDecay::LoadDecayTable()", "HAD_RDM_000",
>> 861 FatalException, "Selected decay mode does not exist");
>> 862 } // switch
>> 863 } // if char == W
>> 864 } // if char == P
>> 865 } // if char != #
>> 866 } // While
>> 867
>> 868 // Go through the decay table and make sure that the branching ratios are
>> 869 // correctly normalised.
>> 870
>> 871 G4VDecayChannel* theChannel = 0;
>> 872 G4NuclearDecayChannel* theNuclearDecayChannel = 0;
>> 873 G4String mode = "";
>> 874
>> 875 G4double theBR = 0.0;
>> 876 for (G4int i = 0; i < theDecayTable->entries(); i++) {
>> 877 theChannel = theDecayTable->GetDecayChannel(i);
>> 878 theNuclearDecayChannel = static_cast<G4NuclearDecayChannel*>(theChannel);
>> 879 theDecayMode = theNuclearDecayChannel->GetDecayMode();
>> 880
>> 881 if (theDecayMode != IT) {
>> 882 theBR = theChannel->GetBR();
>> 883 theChannel->SetBR(theBR*modeTotalBR[theDecayMode]/modeSumBR[theDecayMode]);
>> 884 }
>> 885 }
>> 886 } // if (DecaySchemeFile)
>> 887 DecaySchemeFile.close();
>> 888
>> 889
>> 890 if (!found && E > 0.) {
>> 891 // cases where IT cascade for exited isotopes without entry in RDM database
>> 892 // Decay mode is isomeric transition.
>> 893 //
>> 894 G4ITDecayChannel *anITChannel = new G4ITDecayChannel
>> 895 (GetVerboseLevel(), (const G4Ions*) &theParentNucleus, 1);
>> 896 anITChannel->SetICM(applyICM);
>> 897 anITChannel->SetARM(applyARM);
>> 898 anITChannel->SetHLThreshold(halflifethreshold);
>> 899 theDecayTable->Insert(anITChannel);
>> 900 }
>> 901 if (!theDecayTable) {
>> 902 //
>> 903 // There is no radioactive decay data for this nucleus. Return a null
>> 904 // decay table.
>> 905 //
>> 906 G4cerr <<"G4RadoactiveDecay::LoadDecayTable() : cannot find ion radioactive decay file " <<G4endl;
>> 907 theDecayTable = 0;
>> 908 return theDecayTable;
>> 909 }
>> 910 if (theDecayTable && GetVerboseLevel()>1)
>> 911 {
>> 912 G4cout <<"G4RadioactiveDecay::LoadDecayTable()" << G4endl;
>> 913 G4cout << " No. of entries: "<< theDecayTable->entries() <<G4endl;
>> 914 theDecayTable ->DumpInfo();
>> 915 }
>> 916
>> 917 return theDecayTable;
305 } 918 }
>> 919 ////////////////////////////////////////////////////////////////////
>> 920 //
>> 921 void G4RadioactiveDecay::AddUserDecayDataFile(G4int Z, G4int A,G4String filename)
>> 922 { if (Z<1 || A<2) {
>> 923 G4cout<<"Z and A not valid!"<<G4endl;
>> 924 }
>> 925
>> 926 std::ifstream DecaySchemeFile(filename);
>> 927 if (DecaySchemeFile){
>> 928 G4int ID_ion=A*1000+Z;
>> 929 theUserRadioactiveDataFiles[ID_ion]=filename;
>> 930 }
>> 931 else {
>> 932 G4cout<<"The file "<<filename<<" does not exist!"<<G4endl;
>> 933 }
>> 934 }
>> 935 ////////////////////////////////////////////////////////////////////////
>> 936 //
306 937
307 938
308 void 939 void
309 G4RadioactiveDecay::SetDecayRate(G4int theZ, G 940 G4RadioactiveDecay::SetDecayRate(G4int theZ, G4int theA, G4double theE,
310 G4int theG, std::vector<G4double>& the << 941 G4int theG, std::vector<G4double> theRates,
311 std::vector<G4double>& theTaos) << 942 std::vector<G4double> theTaos)
312 // Why not make this a method of G4Radioactiv <<
313 { 943 {
314 //fill the decay rate vector 944 //fill the decay rate vector
315 ratesToDaughter.SetZ(theZ); << 945 theDecayRate.SetZ(theZ);
316 ratesToDaughter.SetA(theA); << 946 theDecayRate.SetA(theA);
317 ratesToDaughter.SetE(theE); << 947 theDecayRate.SetE(theE);
318 ratesToDaughter.SetGeneration(theG); << 948 theDecayRate.SetGeneration(theG);
319 ratesToDaughter.SetDecayRateC(theCoefficient << 949 theDecayRate.SetDecayRateC(theRates);
320 ratesToDaughter.SetTaos(theTaos); << 950 theDecayRate.SetTaos(theTaos);
321 } 951 }
322 952
>> 953 //////////////////////////////////////////////////////////////////////////
>> 954 //
>> 955 void
>> 956 G4RadioactiveDecay::AddDecayRateTable(const G4ParticleDefinition& theParentNucleus)
>> 957 {
>> 958 // 1) To calculate all the coefficiecies required to derive the
>> 959 // radioactivities for all progeny of theParentNucleus
>> 960 //
>> 961 // 2) Add the coefficiencies to the decay rate table vector
>> 962 //
323 963
324 void G4RadioactiveDecay:: << 964 //
325 CalculateChainsFromParent(const G4ParticleDefi <<
326 { <<
327 // Use extended Bateman equation to calculat <<
328 // progeny of theParentNucleus. The coeffic <<
329 // calculated using the method of P. Truscot <<
330 // DERA Technical Note DERA/CIS/CIS2/7/36/4/ <<
331 // Coefficients are then added to the decay <<
332 <<
333 // Create and initialise variables used in t 965 // Create and initialise variables used in the method.
>> 966 //
>> 967
334 theDecayRateVector.clear(); 968 theDecayRateVector.clear();
335 969
336 G4int nGeneration = 0; 970 G4int nGeneration = 0;
337 << 971 std::vector<G4double> rates;
338 std::vector<G4double> taos; 972 std::vector<G4double> taos;
339 973
340 // Dimensionless A coefficients of Eqs. 4.24 << 974 // start rate is -1.
341 std::vector<G4double> Acoeffs; << 975 // Eq.4.26 of the Technical Note
342 << 976 rates.push_back(-1.);
343 // According to Eq. 4.26 the first coefficie << 977 //
344 Acoeffs.push_back(-1.); << 978 //
345 << 979 G4int A = ((const G4Ions*)(&theParentNucleus))->GetAtomicMass();
346 const G4Ions* ion = static_cast<const G4Ions << 980 G4int Z = ((const G4Ions*)(&theParentNucleus))->GetAtomicNumber();
347 G4int A = ion->GetAtomicMass(); << 981 G4double E = ((const G4Ions*)(&theParentNucleus))->GetExcitationEnergy();
348 G4int Z = ion->GetAtomicNumber(); << 982 G4double tao = theParentNucleus.GetPDGLifeTime();
349 G4double E = ion->GetExcitationEnergy(); << 983 if (tao < 0.) tao = 1e-30;
350 G4double tao = ion->GetPDGLifeTime(); <<
351 if (tao < 0.) tao = 1e-100; <<
352 taos.push_back(tao); 984 taos.push_back(tao);
353 G4int nEntry = 0; 985 G4int nEntry = 0;
354 986
355 // Fill the decay rate container (G4Radioact << 987 //fill the decay rate with the intial isotope data
356 // isotope data << 988 SetDecayRate(Z,A,E,nGeneration,rates,taos);
357 SetDecayRate(Z,A,E,nGeneration,Acoeffs,taos) <<
358 989
359 // store the decay rate in decay rate vector 990 // store the decay rate in decay rate vector
360 theDecayRateVector.push_back(ratesToDaughter << 991 theDecayRateVector.push_back(theDecayRate);
361 ++nEntry; << 992 nEntry++;
>> 993
>> 994 // now start treating the sencondary generations..
362 995
363 // Now start treating the secondary generati <<
364 G4bool stable = false; 996 G4bool stable = false;
>> 997 G4int i;
365 G4int j; 998 G4int j;
366 G4VDecayChannel* theChannel = 0; 999 G4VDecayChannel* theChannel = 0;
367 G4NuclearDecay* theNuclearDecayChannel = 0; << 1000 G4NuclearDecayChannel* theNuclearDecayChannel = 0;
368 << 1001 G4ITDecayChannel* theITChannel = 0;
369 G4ITDecay* theITChannel = 0; << 1002 G4BetaMinusDecayChannel *theBetaMinusChannel = 0;
370 G4BetaMinusDecay* theBetaMinusChannel = 0; << 1003 G4BetaPlusDecayChannel *theBetaPlusChannel = 0;
371 G4BetaPlusDecay* theBetaPlusChannel = 0; << 1004 G4AlphaDecayChannel *theAlphaChannel = 0;
372 G4AlphaDecay* theAlphaChannel = 0; <<
373 G4ProtonDecay* theProtonChannel = 0; <<
374 G4TritonDecay* theTritonChannel = 0; <<
375 G4NeutronDecay* theNeutronChannel = 0; <<
376 G4SFDecay* theFissionChannel = 0; <<
377 <<
378 G4RadioactiveDecayMode theDecayMode; 1005 G4RadioactiveDecayMode theDecayMode;
379 G4double theBR = 0.0; 1006 G4double theBR = 0.0;
380 G4int AP = 0; 1007 G4int AP = 0;
381 G4int ZP = 0; 1008 G4int ZP = 0;
382 G4int AD = 0; 1009 G4int AD = 0;
383 G4int ZD = 0; 1010 G4int ZD = 0;
384 G4double EP = 0.; 1011 G4double EP = 0.;
385 std::vector<G4double> TP; 1012 std::vector<G4double> TP;
386 std::vector<G4double> RP; // A coefficient << 1013 std::vector<G4double> RP;
387 G4ParticleDefinition *theDaughterNucleus; 1014 G4ParticleDefinition *theDaughterNucleus;
388 G4double daughterExcitation; 1015 G4double daughterExcitation;
389 G4double nearestEnergy = 0.0; <<
390 G4int nearestLevelIndex = 0; <<
391 G4ParticleDefinition *aParentNucleus; 1016 G4ParticleDefinition *aParentNucleus;
392 G4IonTable* theIonTable; 1017 G4IonTable* theIonTable;
393 G4DecayTable* parentDecayTable; << 1018 G4DecayTable *aTempDecayTable;
394 G4double theRate; 1019 G4double theRate;
395 G4double TaoPlus; 1020 G4double TaoPlus;
396 G4int nS = 0; // Running index of fir << 1021 G4int nS = 0;
397 G4int nT = nEntry; // Total number of deca << 1022 G4int nT = nEntry;
398 const G4int nMode = G4RadioactiveDecayModeSi << 1023 G4double brs[7];
399 G4double brs[nMode]; <<
400 // 1024 //
401 theIonTable = G4ParticleTable::GetParticleTa << 1025 theIonTable =
>> 1026 (G4IonTable*)(G4ParticleTable::GetParticleTable()->GetIonTable());
402 1027
403 G4int loop = 0; << 1028 while (!stable) {
404 while (!stable) { /* Loop checking, 01.09. <<
405 loop++; <<
406 if (loop > 10000) { <<
407 G4Exception("G4RadioactiveDecay::Calcula <<
408 JustWarning, "While loop cou <<
409 break; <<
410 } <<
411 nGeneration++; 1029 nGeneration++;
412 for (j = nS; j < nT; ++j) { << 1030 for (j = nS; j< nT; j++) {
413 // First time through, get data for pare <<
414 ZP = theDecayRateVector[j].GetZ(); 1031 ZP = theDecayRateVector[j].GetZ();
415 AP = theDecayRateVector[j].GetA(); 1032 AP = theDecayRateVector[j].GetA();
416 EP = theDecayRateVector[j].GetE(); 1033 EP = theDecayRateVector[j].GetE();
417 RP = theDecayRateVector[j].GetDecayRateC 1034 RP = theDecayRateVector[j].GetDecayRateC();
418 TP = theDecayRateVector[j].GetTaos(); << 1035 TP = theDecayRateVector[j].GetTaos();
419 if (GetVerboseLevel() > 1) { << 1036 if (GetVerboseLevel()>0){
420 G4cout << "G4RadioactiveDecay::Calcula << 1037 G4cout <<"G4RadioactiveDecay::AddDecayRateTable : "
421 << ZP << ", " << AP << ", " << << 1038 << " daughters of ("<< ZP <<", "<<AP<<", "
422 << ") are being calculated, ge << 1039 << EP <<") "
423 << G4endl; << 1040 << " are being calculated. "
>> 1041 <<" generation = "
>> 1042 << nGeneration << G4endl;
424 } 1043 }
425 <<
426 aParentNucleus = theIonTable->GetIon(ZP, 1044 aParentNucleus = theIonTable->GetIon(ZP,AP,EP);
427 parentDecayTable = GetDecayTable(aParent << 1045 if (!IsLoaded(*aParentNucleus)){
428 if (nullptr == parentDecayTable) { conti << 1046 aParentNucleus->SetDecayTable(LoadDecayTable(*aParentNucleus));
429 <<
430 G4DecayTable* summedDecayTable = new G4D <<
431 // This instance of G4DecayTable is for <<
432 // channels. It will contain one decay <<
433 // (alpha, beta, etc.); its branching ra <<
434 // branching ratios for that type of dec <<
435 // halflife of a particular channel is l <<
436 // that channel will be inserted specifi <<
437 // ratio will not be included in the abo <<
438 // This instance is not used to perform <<
439 <<
440 for (G4int k = 0; k < nMode; ++k) brs[k] <<
441 <<
442 // Go through the decay table and sum al <<
443 for (G4int i = 0; i < parentDecayTable-> <<
444 theChannel = parentDecayTable->GetDeca <<
445 theNuclearDecayChannel = static_cast<G <<
446 theDecayMode = theNuclearDecayChannel- <<
447 daughterExcitation = theNuclearDecayCh <<
448 theDaughterNucleus = theNuclearDecayCh <<
449 AD = ((const G4Ions*)(theDaughterNucle <<
450 ZD = ((const G4Ions*)(theDaughterNucle <<
451 const G4LevelManager* levelManager = <<
452 G4NuclearLevelData::GetInstance()->G <<
453 <<
454 // Check each nuclide to see if it is <<
455 // If so, add it to the decay chain by <<
456 // summedDecayTable. If not, just add <<
457 if (levelManager->NumberOfTransitions( <<
458 nearestEnergy = levelManager->Neares <<
459 if ((std::abs(daughterExcitation - n <<
460 (std::abs(daughterExcitation - nearest <<
461 // Level half-life is in ns and th <<
462 // by default, user can set it via <<
463 nearestLevelIndex = (G4int)levelMa <<
464 if (levelManager->LifeTime(nearest <<
465 // save the metastable decay cha <<
466 summedDecayTable->Insert(theChan <<
467 } else { <<
468 brs[theDecayMode] += theChannel- <<
469 } <<
470 } else { <<
471 brs[theDecayMode] += theChannel->G <<
472 } <<
473 } else { <<
474 brs[theDecayMode] += theChannel->Get <<
475 } <<
476 <<
477 } // Combine decay channels (loop i) <<
478 <<
479 brs[BetaPlus] = brs[BetaPlus]+brs[Kshell <<
480 brs[KshellEC] = brs[LshellEC] = brs[Mshe <<
481 for (G4int i = 0; i < nMode; ++i) { <<
482 if (brs[i] > 0.) { <<
483 switch (i) { <<
484 case IT: <<
485 // Decay mode is isomeric transiti <<
486 theITChannel = new G4ITDecay(aPare <<
487 <<
488 summedDecayTable->Insert(theITChan <<
489 break; <<
490 <<
491 case BetaMinus: <<
492 // Decay mode is beta- <<
493 theBetaMinusChannel = new G4BetaMi <<
494 <<
495 <<
496 summedDecayTable->Insert(theBetaMi <<
497 break; <<
498 <<
499 case BetaPlus: <<
500 // Decay mode is beta+ + EC. <<
501 theBetaPlusChannel = new G4BetaPlu <<
502 <<
503 <<
504 summedDecayTable->Insert(theBetaPl <<
505 break; <<
506 <<
507 case Alpha: <<
508 // Decay mode is alpha. <<
509 theAlphaChannel = new G4AlphaDecay <<
510 <<
511 summedDecayTable->Insert(theAlphaC <<
512 break; <<
513 <<
514 case Proton: <<
515 // Decay mode is proton. <<
516 theProtonChannel = new G4ProtonDec <<
517 <<
518 summedDecayTable->Insert(theProton <<
519 break; <<
520 <<
521 case Neutron: <<
522 // Decay mode is neutron. <<
523 theNeutronChannel = new G4NeutronD <<
524 <<
525 summedDecayTable->Insert(theNeutro <<
526 break; <<
527 <<
528 case SpFission: <<
529 // Decay mode is spontaneous fissi <<
530 theFissionChannel = new G4SFDecay( <<
531 <<
532 summedDecayTable->Insert(theFissio <<
533 break; <<
534 <<
535 case BDProton: <<
536 // Not yet implemented <<
537 break; <<
538 <<
539 case BDNeutron: <<
540 // Not yet implemented <<
541 break; <<
542 <<
543 case Beta2Minus: <<
544 // Not yet implemented <<
545 break; <<
546 <<
547 case Beta2Plus: <<
548 // Not yet implemented <<
549 break; <<
550 <<
551 case Proton2: <<
552 // Not yet implemented <<
553 break; <<
554 <<
555 case Neutron2: <<
556 // Not yet implemented <<
557 break; <<
558 <<
559 case Triton: <<
560 // Decay mode is Triton. <<
561 theTritonChannel = new G4TritonDec <<
562 <<
563 summedDecayTable->Insert(theTriton <<
564 break; <<
565 <<
566 default: <<
567 break; <<
568 } <<
569 } <<
570 } 1047 }
>> 1048
>> 1049 G4DecayTable* theDecayTable = new G4DecayTable();
>> 1050 aTempDecayTable = aParentNucleus->GetDecayTable();
>> 1051 for (i=0; i< 7; i++) brs[i] = 0.0;
571 1052
572 // loop over all branches in summedDecay <<
573 // 1053 //
574 for (G4int i = 0; i < summedDecayTable-> << 1054 // Go through the decay table and to combine the same decay channels
575 theChannel = summedDecayTable->GetDeca << 1055 //
576 theNuclearDecayChannel = static_cast<G << 1056 for (i=0; i<aTempDecayTable->entries(); i++) {
577 theBR = theChannel->GetBR(); << 1057 theChannel = aTempDecayTable->GetDecayChannel(i);
578 theDaughterNucleus = theNuclearDecayCh << 1058 theNuclearDecayChannel = static_cast<G4NuclearDecayChannel *>(theChannel);
579 << 1059 theDecayMode = theNuclearDecayChannel->GetDecayMode();
580 // First check if the decay of the ori << 1060 daughterExcitation = theNuclearDecayChannel->GetDaughterExcitation ();
581 // if true create a new ground-state n << 1061 theDaughterNucleus = theNuclearDecayChannel->GetDaughterNucleus () ;
582 if (theNuclearDecayChannel->GetDecayMo << 1062 AD = ((const G4Ions*)(theDaughterNucleus))->GetAtomicMass();
583 << 1063 ZD = ((const G4Ions*)(theDaughterNucleus))->GetAtomicNumber();
584 A = ((const G4Ions*)(theDaughterNucl << 1064 G4NuclearLevelManager* levelManager =
585 Z = ((const G4Ions*)(theDaughterNucl << 1065 G4NuclearLevelStore::GetInstance()->GetManager(ZD, AD);
586 theDaughterNucleus=theIonTable->GetI << 1066 if (levelManager->NumberOfLevels() ) {
>> 1067 const G4NuclearLevel* level =
>> 1068 levelManager->NearestLevel (daughterExcitation);
>> 1069
>> 1070 if (std::abs(daughterExcitation - level->Energy()) < levelTolerance) {
>> 1071 // Level half-life is in ns and the threshold is set to 1 micros by default, user can set it via the UI command
>> 1072 if (level->HalfLife()*ns >= halflifethreshold ){
>> 1073 // save the metastable nucleus
>> 1074 theDecayTable->Insert(theChannel);
>> 1075 }
>> 1076 else{
>> 1077 brs[theDecayMode] += theChannel->GetBR();
>> 1078 }
>> 1079 }
>> 1080 else {
>> 1081 brs[theDecayMode] += theChannel->GetBR();
>> 1082 }
>> 1083 }
>> 1084 else{
>> 1085 brs[theDecayMode] += theChannel->GetBR();
>> 1086 }
>> 1087 }
>> 1088 brs[2] = brs[2]+brs[3]+brs[4]+brs[5];
>> 1089 brs[3] = brs[4] =brs[5] = 0.0;
>> 1090 for (i= 0; i<7; i++){
>> 1091 if (brs[i] > 0.) {
>> 1092 switch ( i ) {
>> 1093 case 0:
>> 1094 //
>> 1095 //
>> 1096 // Decay mode is isomeric transition.
>> 1097 //
>> 1098
>> 1099 theITChannel = new G4ITDecayChannel
>> 1100 (0, (const G4Ions*) aParentNucleus, brs[0]);
>> 1101
>> 1102 theDecayTable->Insert(theITChannel);
>> 1103 break;
>> 1104
>> 1105 case 1:
>> 1106 //
>> 1107 //
>> 1108 // Decay mode is beta-.
>> 1109 //
>> 1110 theBetaMinusChannel = new G4BetaMinusDecayChannel (0, aParentNucleus,
>> 1111 brs[1], 0.*MeV, 0.*MeV, 1, false, 0);
>> 1112 theDecayTable->Insert(theBetaMinusChannel);
>> 1113
>> 1114 break;
>> 1115
>> 1116 case 2:
>> 1117 //
>> 1118 //
>> 1119 // Decay mode is beta+ + EC.
>> 1120 //
>> 1121 theBetaPlusChannel = new G4BetaPlusDecayChannel (GetVerboseLevel(), aParentNucleus,
>> 1122 brs[2], 0.*MeV, 0.*MeV, 1, false, 0);
>> 1123 theDecayTable->Insert(theBetaPlusChannel);
>> 1124 break;
>> 1125
>> 1126 case 6:
>> 1127 //
>> 1128 //
>> 1129 // Decay mode is alpha.
>> 1130 //
>> 1131 theAlphaChannel = new G4AlphaDecayChannel(GetVerboseLevel(),
>> 1132 aParentNucleus,
>> 1133 brs[6], 0.*MeV, 0.*MeV);
>> 1134 theDecayTable->Insert(theAlphaChannel);
>> 1135 break;
>> 1136
>> 1137 default:
>> 1138 break;
>> 1139 }
>> 1140 }
>> 1141 }
>> 1142 //
>> 1143 // loop over all branches in theDecayTable
>> 1144 //
>> 1145 for (i = 0; i < theDecayTable->entries(); i++){
>> 1146 theChannel = theDecayTable->GetDecayChannel(i);
>> 1147 theNuclearDecayChannel = static_cast<G4NuclearDecayChannel*>(theChannel);
>> 1148 theBR = theChannel->GetBR();
>> 1149 theDaughterNucleus = theNuclearDecayChannel->GetDaughterNucleus();
>> 1150 // first check if the decay of the original nucleus is an IT channel, if true create a new groud-level nucleus
>> 1151 if (theNuclearDecayChannel->GetDecayMode() == IT && nGeneration == 1) {
>> 1152 A = ((const G4Ions*)(theDaughterNucleus))->GetAtomicMass();
>> 1153 Z = ((const G4Ions*)(theDaughterNucleus))->GetAtomicNumber();
>> 1154 theDaughterNucleus=theIonTable->GetIon(Z,A,0.);
587 } 1155 }
588 if (IsApplicable(*theDaughterNucleus) << 1156 if (IsApplicable(*theDaughterNucleus) &&
>> 1157 theBR &&
589 aParentNucleus != theDaughterNucle 1158 aParentNucleus != theDaughterNucleus) {
590 // need to make sure daughter has de << 1159 // need to make sure daugher has decaytable
591 parentDecayTable = GetDecayTable(the << 1160 if (!IsLoaded(*theDaughterNucleus)){
592 if (nullptr != parentDecayTable && p << 1161 theDaughterNucleus->SetDecayTable(LoadDecayTable(*theDaughterNucleus));
593 A = ((const G4Ions*)(theDaughterNu << 1162 }
594 Z = ((const G4Ions*)(theDaughterNu << 1163 if (theDaughterNucleus->GetDecayTable()->entries() ) {
595 E = ((const G4Ions*)(theDaughterNu << 1164 //
596 << 1165 A = ((const G4Ions*)(theDaughterNucleus))->GetAtomicMass();
597 TaoPlus = theDaughterNucleus->GetP << 1166 Z = ((const G4Ions*)(theDaughterNucleus))->GetAtomicNumber();
598 if (TaoPlus <= 0.) TaoPlus = 1e-1 << 1167 E = ((const G4Ions*)(theDaughterNucleus))->GetExcitationEnergy();
599 << 1168
600 // first set the taos, one simply << 1169 TaoPlus = theDaughterNucleus->GetPDGLifeTime();
601 taos.clear(); << 1170 // cout << TaoPlus <<G4endl;
602 taos = TP; // load lifetimes of << 1171 if (TaoPlus <= 0.) TaoPlus = 1e-30;
603 std::size_t k; << 1172
604 //check that TaoPlus differs from << 1173 // first set the taos, one simply need to add to the parent ones
605 //for (k = 0; k < TP.size(); ++k){ << 1174 taos.clear();
606 //if (std::abs((TaoPlus-TP[k])/TP[ << 1175 taos = TP;
607 //} << 1176 taos.push_back(TaoPlus);
608 taos.push_back(TaoPlus); // add d << 1177 // now calculate the coefficiencies
609 // now calculate the coefficiencie << 1178 //
610 // << 1179 // they are in two parts, first the less than n ones
611 // they are in two parts, first th << 1180 // Eq 4.24 of the TN
612 // Eq 4.24 of the TN << 1181 rates.clear();
613 Acoeffs.clear(); << 1182 size_t k;
614 long double ta1,ta2; << 1183 for (k = 0; k < RP.size(); k++){
615 ta2 = (long double)TaoPlus; << 1184 if ((TP[k]-TaoPlus) == 0.) {
616 for (k = 0; k < RP.size(); ++k){ << 1185 theRate = 1e30;
617 ta1 = (long double)TP[k]; // << 1186 } else {
618 if (ta1 == ta2) { << 1187 theRate = TP[k]/(TP[k]-TaoPlus) * theBR * RP[k];
619 theRate = 1.e100; << 1188 }
620 } else { << 1189 rates.push_back(theRate);
621 theRate = ta1/(ta1-ta2); << 1190 }
622 } << 1191 //
623 theRate = theRate * theBR * RP[k << 1192 // the sencond part: the n:n coefficiency
624 Acoeffs.push_back(theRate); << 1193 // Eq 4.25 of the TN. Note Yn+1 is zero apart from Y1 which is -1 as treated at line 1013
625 } << 1194 //
626 << 1195 theRate = 0.;
627 // the second part: the n:n coeffi << 1196 G4double aRate;
628 // Eq 4.25 of the TN. Note Yn+1 i << 1197 for (k = 0; k < RP.size(); k++){
629 // as treated at line 1013 << 1198 if ((TP[k]-TaoPlus) == 0.) {
630 theRate = 0.; << 1199 aRate = 1e30;
631 long double aRate, aRate1; << 1200 } else {
632 aRate1 = 0.L; << 1201 aRate = TaoPlus/(TP[k]-TaoPlus) * theBR * RP[k];
633 for (k = 0; k < RP.size(); ++k){ << 1202 }
634 ta1 = (long double)TP[k]; << 1203 theRate -= aRate;
635 if (ta1 == ta2 ) { << 1204 }
636 aRate = 1.e100; << 1205 rates.push_back(theRate);
637 } else { << 1206 SetDecayRate (Z,A,E,nGeneration,rates,taos);
638 aRate = ta2/(ta1-ta2); << 1207 theDecayRateVector.push_back(theDecayRate);
639 } << 1208 nEntry++;
640 aRate = aRate * (long double)(th << 1209 }
641 aRate1 += aRate; << 1210 } // end of testing daughter nucleus
642 } << 1211 } // end of i loop( the branches)
643 theRate = -aRate1; << 1212 delete theDecayTable;
644 Acoeffs.push_back(theRate); <<
645 SetDecayRate (Z,A,E,nGeneration,Ac <<
646 theDecayRateVector.push_back(rates <<
647 nEntry++; <<
648 } // there are entries in the table <<
649 } // nuclide is OK to decay <<
650 } // end of loop (i) over decay table br <<
651 <<
652 delete summedDecayTable; <<
653 1213
654 } // Getting contents of decay rate vector << 1214 } //end of for j loop
655 nS = nT; 1215 nS = nT;
656 nT = nEntry; 1216 nT = nEntry;
657 if (nS == nT) stable = true; 1217 if (nS == nT) stable = true;
658 } // while nuclide is not stable << 1218 }
659 1219
660 // end of while loop << 1220 //end of while loop
661 // the calculation completed here 1221 // the calculation completed here
662 1222
663 1223
664 // fill the first part of the decay rate tab 1224 // fill the first part of the decay rate table
665 // which is the name of the original particl << 1225 // which is the name of the original particle (isotope)
666 chainsFromParent.SetIonName(theParentNucleus << 1226 //
667 << 1227 theDecayRateTable.SetIonName(theParentNucleus.GetParticleName());
>> 1228 //
>> 1229 //
668 // now fill the decay table with the newly c 1230 // now fill the decay table with the newly completed decay rate vector
669 chainsFromParent.SetItsRates(theDecayRateVec << 1231 //
>> 1232
>> 1233 theDecayRateTable.SetItsRates(theDecayRateVector);
670 1234
>> 1235 //
671 // finally add the decayratetable to the tab 1236 // finally add the decayratetable to the tablevector
672 theParentChainTable.push_back(chainsFromPare << 1237 //
>> 1238 theDecayRateTableVector.push_back(theDecayRateTable);
673 } 1239 }
674 1240
675 ////////////////////////////////////////////// 1241 ////////////////////////////////////////////////////////////////////////////////
676 // 1242 // //
677 // SetSourceTimeProfile 1243 // SetSourceTimeProfile //
678 // read in the source time profile functio 1244 // read in the source time profile function (histogram) //
679 // 1245 // //
680 ////////////////////////////////////////////// 1246 ////////////////////////////////////////////////////////////////////////////////
681 1247
682 void G4RadioactiveDecay::SetSourceTimeProfile( << 1248 void G4RadioactiveDecay::SetSourceTimeProfile(G4String filename)
683 { 1249 {
684 std::ifstream infile ( filename, std::ios::i 1250 std::ifstream infile ( filename, std::ios::in );
685 if (!infile) { 1251 if (!infile) {
686 G4ExceptionDescription ed; 1252 G4ExceptionDescription ed;
687 ed << " Could not open file " << filename 1253 ed << " Could not open file " << filename << G4endl;
688 G4Exception("G4RadioactiveDecay::SetSource 1254 G4Exception("G4RadioactiveDecay::SetSourceTimeProfile()", "HAD_RDM_001",
689 FatalException, ed); 1255 FatalException, ed);
690 } 1256 }
691 1257
>> 1258 // if (!infile) G4Exception(__FILE__, G4inttostring(__LINE__), FatalException,
>> 1259 // "Unable to open source data file");
>> 1260
692 G4double bin, flux; 1261 G4double bin, flux;
693 NSourceBin = -1; 1262 NSourceBin = -1;
694 << 1263 while (infile >> bin >> flux ) {
695 G4int loop = 0; <<
696 while (infile >> bin >> flux) { /* Loop che <<
697 loop++; <<
698 if (loop > 10000) { <<
699 G4Exception("G4RadioactiveDecay::SetSour <<
700 JustWarning, "While loop cou <<
701 break; <<
702 } <<
703 <<
704 NSourceBin++; 1264 NSourceBin++;
705 if (NSourceBin > 99) { 1265 if (NSourceBin > 99) {
706 G4Exception("G4RadioactiveDecay::SetSour 1266 G4Exception("G4RadioactiveDecay::SetSourceTimeProfile()", "HAD_RDM_002",
707 FatalException, "Input sourc 1267 FatalException, "Input source time file too big (>100 rows)");
708 1268
>> 1269 // G4Exception(__FILE__, G4inttostring(__LINE__), FatalException,
>> 1270 // "input source data file too big (>100 rows)" );
709 } else { 1271 } else {
710 SBin[NSourceBin] = bin * s; // Conver << 1272 SBin[NSourceBin] = bin * s;
711 SProfile[NSourceBin] = flux; // Dimens << 1273 SProfile[NSourceBin] = flux;
712 } 1274 }
713 } 1275 }
714 << 1276 SetAnalogueMonteCarlo(0);
715 AnalogueMC = false; <<
716 infile.close(); 1277 infile.close();
717 1278
718 #ifdef G4VERBOSE 1279 #ifdef G4VERBOSE
719 if (GetVerboseLevel() > 2) << 1280 if (GetVerboseLevel()>1)
720 G4cout <<" Source Timeprofile Nbin = " << << 1281 {G4cout <<" Source Timeprofile Nbin = " << NSourceBin <<G4endl;}
721 #endif 1282 #endif
722 } 1283 }
723 1284
724 ////////////////////////////////////////////// 1285 ////////////////////////////////////////////////////////////////////////////////
725 // 1286 // //
726 // SetDecayBiasProfile 1287 // SetDecayBiasProfile //
727 // read in the decay bias scheme function ( 1288 // read in the decay bias scheme function (histogram) //
728 // 1289 // //
729 ////////////////////////////////////////////// 1290 ////////////////////////////////////////////////////////////////////////////////
730 1291
731 void G4RadioactiveDecay::SetDecayBias(const G4 << 1292 void G4RadioactiveDecay::SetDecayBias(G4String filename)
732 { 1293 {
>> 1294
733 std::ifstream infile(filename, std::ios::in) 1295 std::ifstream infile(filename, std::ios::in);
734 if (!infile) G4Exception("G4RadioactiveDecay << 1296 if (!infile) G4Exception("G4RadioactiveDecay::SetDecayBias()", "HAD_RDM_003",
735 FatalException, "Un 1297 FatalException, "Unable to open bias data file" );
>> 1298 // if (!infile) G4Exception(__FILE__, G4inttostring(__LINE__), FatalException,
>> 1299 // "Unable to open bias data file" );
736 1300
737 G4double bin, flux; 1301 G4double bin, flux;
738 G4int dWindows = 0; 1302 G4int dWindows = 0;
739 G4int i ; 1303 G4int i ;
740 1304
741 theRadioactivityTables.clear(); 1305 theRadioactivityTables.clear();
742 1306
743 NDecayBin = -1; 1307 NDecayBin = -1;
744 << 1308 while (infile >> bin >> flux ) {
745 G4int loop = 0; <<
746 while (infile >> bin >> flux ) { /* Loop ch <<
747 NDecayBin++; 1309 NDecayBin++;
748 loop++; <<
749 if (loop > 10000) { <<
750 G4Exception("G4RadioactiveDecay::SetDeca <<
751 JustWarning, "While loop cou <<
752 break; <<
753 } <<
754 <<
755 if (NDecayBin > 99) { 1310 if (NDecayBin > 99) {
756 G4Exception("G4RadioactiveDecay::SetDeca << 1311 G4Exception("G4RadioactiveDecay::SetDecayBias()", "HAD_RDM_004",
757 FatalException, "Input bias 1312 FatalException, "Input bias file too big (>100 rows)" );
>> 1313 // G4Exception(__FILE__, G4inttostring(__LINE__), FatalException,
>> 1314 // "input bias data file too big (>100 rows)" );
758 } else { 1315 } else {
759 DBin[NDecayBin] = bin * s; // Conver << 1316 DBin[NDecayBin] = bin * s;
760 DProfile[NDecayBin] = flux; // Dimens << 1317 DProfile[NDecayBin] = flux;
761 if (flux > 0.) { 1318 if (flux > 0.) {
762 decayWindows[NDecayBin] = dWindows; 1319 decayWindows[NDecayBin] = dWindows;
763 dWindows++; 1320 dWindows++;
764 G4RadioactivityTable *rTable = new G4R 1321 G4RadioactivityTable *rTable = new G4RadioactivityTable() ;
765 theRadioactivityTables.push_back(rTabl 1322 theRadioactivityTables.push_back(rTable);
766 } 1323 }
767 } 1324 }
768 } 1325 }
769 for ( i = 1; i<= NDecayBin; ++i) DProfile[i] << 1326 for ( i = 1; i<= NDecayBin; i++) DProfile[i] += DProfile[i-1];
770 for ( i = 0; i<= NDecayBin; ++i) DProfile[i] << 1327 for ( i = 0; i<= NDecayBin; i++) DProfile[i] /= DProfile[NDecayBin];
771 // Normalize so e <<
772 // converted to accumulated probabilities 1328 // converted to accumulated probabilities
773 1329
774 AnalogueMC = false; << 1330 SetAnalogueMonteCarlo(0);
775 infile.close(); 1331 infile.close();
776 1332
777 #ifdef G4VERBOSE 1333 #ifdef G4VERBOSE
778 if (GetVerboseLevel() > 2) << 1334 if (GetVerboseLevel()>1)
779 G4cout <<" Decay Bias Profile Nbin = " << << 1335 {G4cout <<" Decay Bias Profile Nbin = " << NDecayBin <<G4endl;}
780 #endif 1336 #endif
781 } 1337 }
782 1338
783 ////////////////////////////////////////////// 1339 ////////////////////////////////////////////////////////////////////////////////
784 // 1340 // //
785 // DecayIt 1341 // DecayIt //
786 // 1342 // //
787 ////////////////////////////////////////////// 1343 ////////////////////////////////////////////////////////////////////////////////
788 1344
789 G4VParticleChange* 1345 G4VParticleChange*
790 G4RadioactiveDecay::DecayIt(const G4Track& the 1346 G4RadioactiveDecay::DecayIt(const G4Track& theTrack, const G4Step&)
791 { 1347 {
792 // Initialize G4ParticleChange object, get p << 1348 // Initialize the G4ParticleChange object. Get the particle details and the
>> 1349 // decay table.
>> 1350
793 fParticleChangeForRadDecay.Initialize(theTra 1351 fParticleChangeForRadDecay.Initialize(theTrack);
794 fParticleChangeForRadDecay.ProposeWeight(the <<
795 const G4DynamicParticle* theParticle = theTr 1352 const G4DynamicParticle* theParticle = theTrack.GetDynamicParticle();
796 const G4ParticleDefinition* theParticleDef = << 1353 G4ParticleDefinition *theParticleDef = theParticle->GetDefinition();
797 1354
798 // First check whether RDM applies to the cu 1355 // First check whether RDM applies to the current logical volume
799 if (!isAllVolumesMode) { <<
800 if (!std::binary_search(ValidVolumes.begin <<
801 theTrack.GetVolume()->Get <<
802 #ifdef G4VERBOSE <<
803 if (GetVerboseLevel()>0) { <<
804 G4cout <<"G4RadioactiveDecay::DecayIt <<
805 << theTrack.GetVolume()->GetLog <<
806 << " is not selected for the RD <<
807 G4cout << " There are " << ValidVolume <<
808 G4cout << " The Valid volumes are " << <<
809 for (std::size_t i = 0; i< ValidVolume <<
810 G4cout << Va <<
811 } <<
812 #endif <<
813 fParticleChangeForRadDecay.SetNumberOfSe <<
814 1356
815 // Kill the parent particle. << 1357 if (!std::binary_search(ValidVolumes.begin(), ValidVolumes.end(),
816 fParticleChangeForRadDecay.ProposeTrackS << 1358 theTrack.GetVolume()->GetLogicalVolume()->GetName())) {
817 fParticleChangeForRadDecay.ProposeLocalE << 1359 #ifdef G4VERBOSE
818 ClearNumberOfInteractionLengthLeft(); << 1360 if (GetVerboseLevel()>0) {
819 return &fParticleChangeForRadDecay; << 1361 G4cout <<"G4RadioactiveDecay::DecayIt : "
820 } << 1362 << theTrack.GetVolume()->GetLogicalVolume()->GetName()
821 } << 1363 << " is not selected for the RDM"<< G4endl;
822 << 1364 G4cout << " There are " << ValidVolumes.size() << " volumes" << G4endl;
823 // Now check if particle is valid for RDM << 1365 G4cout << " The Valid volumes are " << G4endl;
824 if (!(IsApplicable(*theParticleDef) ) ) { << 1366 for (size_t i = 0; i< ValidVolumes.size(); i++) G4cout << ValidVolumes[i] << G4endl;
825 // Particle is not an ion or is outside th <<
826 #ifdef G4VERBOSE <<
827 if (GetVerboseLevel() > 1) { <<
828 G4cout << "G4RadioactiveDecay::DecayIt : <<
829 << theParticleDef->GetParticleNam <<
830 << " is not an ion or is outside <<
831 << " Set particle change accordin <<
832 << G4endl; <<
833 } 1367 }
834 #endif 1368 #endif
835 fParticleChangeForRadDecay.SetNumberOfSeco 1369 fParticleChangeForRadDecay.SetNumberOfSecondaries(0);
836 1370
837 // Kill the parent particle << 1371 // Kill the parent particle.
838 fParticleChangeForRadDecay.ProposeTrackSta << 1372
>> 1373 fParticleChangeForRadDecay.ProposeTrackStatus( fStopAndKill ) ;
839 fParticleChangeForRadDecay.ProposeLocalEne 1374 fParticleChangeForRadDecay.ProposeLocalEnergyDeposit(0.0);
840 ClearNumberOfInteractionLengthLeft(); 1375 ClearNumberOfInteractionLengthLeft();
841 return &fParticleChangeForRadDecay; 1376 return &fParticleChangeForRadDecay;
842 } 1377 }
843 1378
844 G4DecayTable* theDecayTable = GetDecayTable( << 1379 // now check is the particle is valid for RDM
845 1380
846 if (theDecayTable == nullptr || theDecayTabl << 1381 if (!(IsApplicable(*theParticleDef))) {
847 // No data in the decay table. Set partic << 1382 //
848 // to indicate this. << 1383 // The particle is not a Ion or outside the nucleuslimits for decay
849 #ifdef G4VERBOSE << 1384 //
850 if (GetVerboseLevel() > 1) { << 1385 #ifdef G4VERBOSE
851 G4cout << "G4RadioactiveDecay::DecayIt : << 1386 if (GetVerboseLevel()>0) {
852 << "decay table not defined for " << 1387 G4cerr <<"G4RadioactiveDecay::DecayIt : "
853 << theParticleDef->GetParticleNam << 1388 <<theParticleDef->GetParticleName()
854 << ". Set particle change accordi << 1389 << " is not a valid nucleus for the RDM"<< G4endl;
855 << G4endl; <<
856 } 1390 }
857 #endif 1391 #endif
858 fParticleChangeForRadDecay.SetNumberOfSeco 1392 fParticleChangeForRadDecay.SetNumberOfSecondaries(0);
859 1393
>> 1394 //
860 // Kill the parent particle. 1395 // Kill the parent particle.
861 fParticleChangeForRadDecay.ProposeTrackSta << 1396 //
>> 1397 fParticleChangeForRadDecay.ProposeTrackStatus( fStopAndKill ) ;
862 fParticleChangeForRadDecay.ProposeLocalEne 1398 fParticleChangeForRadDecay.ProposeLocalEnergyDeposit(0.0);
863 ClearNumberOfInteractionLengthLeft(); 1399 ClearNumberOfInteractionLengthLeft();
864 return &fParticleChangeForRadDecay; 1400 return &fParticleChangeForRadDecay;
>> 1401 }
>> 1402
>> 1403 if (!IsLoaded(*theParticleDef))
>> 1404 theParticleDef->SetDecayTable(LoadDecayTable(*theParticleDef));
>> 1405
>> 1406 G4DecayTable* theDecayTable = theParticleDef->GetDecayTable();
>> 1407
>> 1408 if (theDecayTable == 0 || theDecayTable->entries() == 0) {
>> 1409 // There are no data in the decay table. Set the particle change parameters
>> 1410 // to indicate this.
>> 1411 #ifdef G4VERBOSE
>> 1412 if (GetVerboseLevel()>0)
>> 1413 {
>> 1414 G4cerr <<"G4RadioactiveDecay::DecayIt : decay table not defined for ";
>> 1415 G4cerr <<theParticleDef->GetParticleName() <<G4endl;
>> 1416 }
>> 1417 #endif
>> 1418 fParticleChangeForRadDecay.SetNumberOfSecondaries(0);
>> 1419
>> 1420 // Kill the parent particle.
>> 1421 fParticleChangeForRadDecay.ProposeTrackStatus( fStopAndKill ) ;
>> 1422 fParticleChangeForRadDecay.ProposeLocalEnergyDeposit(0.0);
>> 1423 ClearNumberOfInteractionLengthLeft();
>> 1424 return &fParticleChangeForRadDecay;
865 1425
866 } else { 1426 } else {
867 // Data found. Try to decay nucleus 1427 // Data found. Try to decay nucleus
>> 1428 G4double energyDeposit = 0.0;
>> 1429 G4double finalGlobalTime = theTrack.GetGlobalTime();
>> 1430 G4double finalLocalTime = theTrack.GetLocalTime();
>> 1431 G4int index;
>> 1432 G4ThreeVector currentPosition;
>> 1433 currentPosition = theTrack.GetPosition();
>> 1434
>> 1435 // Check whether use Analogue or VR implementation
868 if (AnalogueMC) { 1436 if (AnalogueMC) {
869 G4VRadioactiveDecay::DecayAnalog(theTrac << 1437 #ifdef G4VERBOSE
>> 1438 if (GetVerboseLevel() > 0)
>> 1439 G4cout <<"DecayIt: Analogue MC version " << G4endl;
>> 1440 #endif
870 1441
871 } else { << 1442 G4DecayProducts* products = DoDecay(*theParticleDef);
872 // Proceed with decay using variance red <<
873 G4double energyDeposit = 0.0; <<
874 G4double finalGlobalTime = theTrack.GetG <<
875 G4double finalLocalTime = theTrack.GetLo <<
876 G4int index; <<
877 G4ThreeVector currentPosition; <<
878 currentPosition = theTrack.GetPosition() <<
879 1443
880 G4IonTable* theIonTable; << 1444 // Check if the product is the same as input and kill the track if
881 G4ParticleDefinition* parentNucleus; << 1445 // necessary to prevent infinite loop (11/05/10, F.Lei)
>> 1446 if ( products->entries() == 1) {
>> 1447 fParticleChangeForRadDecay.SetNumberOfSecondaries(0);
>> 1448 fParticleChangeForRadDecay.ProposeTrackStatus( fStopAndKill ) ;
>> 1449 fParticleChangeForRadDecay.ProposeLocalEnergyDeposit(0.0);
>> 1450 ClearNumberOfInteractionLengthLeft();
>> 1451 return &fParticleChangeForRadDecay;
>> 1452 }
>> 1453
>> 1454 // Get parent particle information and boost the decay products to the
>> 1455 // laboratory frame based on this information.
>> 1456 G4double ParentEnergy = theParticle->GetTotalEnergy();
>> 1457 G4ThreeVector ParentDirection(theParticle->GetMomentumDirection());
>> 1458
>> 1459 if (theTrack.GetTrackStatus() == fStopButAlive) {
>> 1460 // The particle is decayed at rest.
>> 1461 // since the time is still for rest particle in G4 we need to add the
>> 1462 // additional time lapsed between the particle come to rest and the
>> 1463 // actual decay. This time is simply sampled with the mean-life of
>> 1464 // the particle. But we need to protect the case PDGTime < 0.
>> 1465 // (F.Lei 11/05/10)
>> 1466 G4double temptime = -std::log( G4UniformRand())
>> 1467 *theParticleDef->GetPDGLifeTime();
>> 1468 if (temptime < 0.) temptime = 0.;
>> 1469 finalGlobalTime += temptime;
>> 1470 finalLocalTime += temptime;
>> 1471 energyDeposit += theParticle->GetKineticEnergy();
>> 1472 } else {
>> 1473 // The particle is decayed in flight (PostStep case).
>> 1474 products->Boost( ParentEnergy, ParentDirection);
>> 1475 }
>> 1476
>> 1477 // Add products in theParticleChangeForRadDecay.
>> 1478 G4int numberOfSecondaries = products->entries();
>> 1479 fParticleChangeForRadDecay.SetNumberOfSecondaries(numberOfSecondaries);
>> 1480 #ifdef G4VERBOSE
>> 1481 if (GetVerboseLevel()>1) {
>> 1482 G4cout <<"G4RadioactiveDecay::DecayIt : Decay vertex :";
>> 1483 G4cout <<" Time: " <<finalGlobalTime/ns <<"[ns]";
>> 1484 G4cout <<" X:" <<(theTrack.GetPosition()).x() /cm <<"[cm]";
>> 1485 G4cout <<" Y:" <<(theTrack.GetPosition()).y() /cm <<"[cm]";
>> 1486 G4cout <<" Z:" <<(theTrack.GetPosition()).z() /cm <<"[cm]";
>> 1487 G4cout << G4endl;
>> 1488 G4cout <<"G4Decay::DecayIt : decay products in Lab. Frame" <<G4endl;
>> 1489 products->DumpInfo();
>> 1490 }
>> 1491 #endif
>> 1492 for (index=0; index < numberOfSecondaries; index++) {
>> 1493 G4Track* secondary = new G4Track(products->PopProducts(),
>> 1494 finalGlobalTime, currentPosition);
>> 1495 secondary->SetGoodForTrackingFlag();
>> 1496 secondary->SetTouchableHandle(theTrack.GetTouchableHandle());
>> 1497 fParticleChangeForRadDecay.AddSecondary(secondary);
>> 1498 }
>> 1499 delete products;
>> 1500 // end of analogue MC algarithm
882 1501
883 // Get decay chains for the given nuclid << 1502 } else {
884 if (!IsRateTableReady(*theParticleDef)) << 1503 // Variance Reduction Method
885 CalculateChainsFromParent(*theParticleDef); << 1504 #ifdef G4VERBOSE
886 GetChainsFromParent(*theParticleDef); << 1505 if (GetVerboseLevel()>0)
>> 1506 G4cout << "DecayIt: Variance Reduction version " << G4endl;
>> 1507 #endif
>> 1508 if (!IsRateTableReady(*theParticleDef)) {
>> 1509 // if the decayrates are not ready, calculate them and
>> 1510 // add to the rate table vector
>> 1511 AddDecayRateTable(*theParticleDef);
>> 1512 }
>> 1513 //retrieve the rates
>> 1514 GetDecayRateTable(*theParticleDef);
887 1515
888 // Declare some of the variables require << 1516 // declare some of the variables required in the implementation
>> 1517 G4ParticleDefinition* parentNucleus;
>> 1518 G4IonTable* theIonTable;
889 G4int PZ; 1519 G4int PZ;
890 G4int PA; 1520 G4int PA;
891 G4double PE; 1521 G4double PE;
892 G4String keyName; <<
893 std::vector<G4double> PT; 1522 std::vector<G4double> PT;
894 std::vector<G4double> PR; 1523 std::vector<G4double> PR;
895 G4double taotime; 1524 G4double taotime;
896 long double decayRate; << 1525 G4double decayRate;
897 1526
898 std::size_t i; << 1527 size_t i;
>> 1528 size_t j;
899 G4int numberOfSecondaries; 1529 G4int numberOfSecondaries;
900 G4int totalNumberOfSecondaries = 0; 1530 G4int totalNumberOfSecondaries = 0;
901 G4double currentTime = 0.; 1531 G4double currentTime = 0.;
902 G4int ndecaych; 1532 G4int ndecaych;
903 G4DynamicParticle* asecondaryparticle; 1533 G4DynamicParticle* asecondaryparticle;
904 std::vector<G4DynamicParticle*> secondar 1534 std::vector<G4DynamicParticle*> secondaryparticles;
905 std::vector<G4double> pw; 1535 std::vector<G4double> pw;
906 std::vector<G4double> ptime; 1536 std::vector<G4double> ptime;
907 pw.clear(); 1537 pw.clear();
908 ptime.clear(); 1538 ptime.clear();
909 1539
910 // Now apply the nucleus splitting << 1540 //now apply the nucleus splitting
911 for (G4int n = 0; n < NSplit; ++n) { << 1541 for (G4int n = 0; n < NSplit; n++) {
912 // Get the decay time following the de 1542 // Get the decay time following the decay probability function
913 // supplied by user << 1543 // suppllied by user
914 G4double theDecayTime = GetDecayTime() 1544 G4double theDecayTime = GetDecayTime();
915 G4int nbin = GetDecayTimeBin(theDecayT 1545 G4int nbin = GetDecayTimeBin(theDecayTime);
916 << 1546
917 // calculate the first part of the wei 1547 // calculate the first part of the weight function
918 G4double weight1 = 1.; 1548 G4double weight1 = 1.;
919 if (nbin == 1) { 1549 if (nbin == 1) {
920 weight1 = 1./DProfile[nbin-1] 1550 weight1 = 1./DProfile[nbin-1]
921 *(DBin[nbin]-DBin[nbin-1]) << 1551 *(DBin[nbin]-DBin[nbin-1])/NSplit;
922 } else if (nbin > 1) { 1552 } else if (nbin > 1) {
923 // Go from nbin to nbin-2 because fl <<
924 weight1 = 1./(DProfile[nbin]-DProfil 1553 weight1 = 1./(DProfile[nbin]-DProfile[nbin-2])
925 *(DBin[nbin]-DBin[nbin-1]) 1554 *(DBin[nbin]-DBin[nbin-1])/NSplit;
926 // weight1 = (probability of choosin <<
927 } 1555 }
>> 1556
928 // it should be calculated in seconds 1557 // it should be calculated in seconds
929 weight1 /= s ; 1558 weight1 /= s ;
930 1559
931 // loop over all the possible secondar 1560 // loop over all the possible secondaries of the nucleus
932 // the first one is itself. 1561 // the first one is itself.
933 for (i = 0; i < theDecayRateVector.siz << 1562 for (i = 0; i < theDecayRateVector.size(); i++) {
934 PZ = theDecayRateVector[i].GetZ(); 1563 PZ = theDecayRateVector[i].GetZ();
935 PA = theDecayRateVector[i].GetA(); 1564 PA = theDecayRateVector[i].GetA();
936 PE = theDecayRateVector[i].GetE(); 1565 PE = theDecayRateVector[i].GetE();
937 PT = theDecayRateVector[i].GetTaos() 1566 PT = theDecayRateVector[i].GetTaos();
938 PR = theDecayRateVector[i].GetDecayR 1567 PR = theDecayRateVector[i].GetDecayRateC();
939 1568
940 // The array of arrays theDecayRateV << 1569 // Calculate the decay rate of the isotope
941 // chains of a given parent nucleus << 1570 // decayRate is the radioactivity of isotope (PZ,PA,PE) at the
942 // nuclide (Z,A,E). << 1571 // time 'theDecayTime'
943 // <<
944 // theDecayRateVector[0] contains th <<
945 // nucleus <<
946 // PZ = ZP <<
947 // PA = AP <<
948 // PE = EP <<
949 // PT[] = {TP} <<
950 // PR[] = {RP} <<
951 // <<
952 // theDecayRateVector[1] contains th <<
953 // generation daughter (Z1,A1,E1). <<
954 // PZ = Z1 <<
955 // PA = A1 <<
956 // PE = E1 <<
957 // PT[] = {TP, T1} <<
958 // PR[] = {RP, R1} <<
959 // <<
960 // theDecayRateVector[2] contains th <<
961 // generation daughter (Z1,A1,E1) an <<
962 // generation daughter to the second <<
963 // PZ = Z2 <<
964 // PA = A2 <<
965 // PE = E2 <<
966 // PT[] = {TP, T1, T2} <<
967 // PR[] = {RP, R1, R2} <<
968 // <<
969 // theDecayRateVector[3] may contain <<
970 // PZ = Z2a <<
971 // PA = A2a <<
972 // PE = E2a <<
973 // PT[] = {TP, T1, T2a} <<
974 // PR[] = {RP, R1, R2a} <<
975 // <<
976 // and so on. <<
977 <<
978 // Calculate the decay rate of the i <<
979 // radioactivity of isotope (PZ,PA,P <<
980 // it will be used to calculate the 1572 // it will be used to calculate the statistical weight of the
981 // decay products of this isotope 1573 // decay products of this isotope
982 1574
983 // For each nuclide, calculate all t << 1575 // G4cout <<"PA= "<< PA << " PZ= " << PZ << " PE= "<< PE <<G4endl;
984 // the parent nuclide << 1576 decayRate = 0.;
985 decayRate = 0.L; << 1577 for (j = 0; j < PT.size(); j++) {
986 for (G4int j = 0; j < G4int(PT.size( << 1578 taotime = GetTaoTime(theDecayTime,PT[j]);
987 taotime = ConvolveSourceTimeProfil << 1579 decayRate -= PR[j] * taotime;
988 decayRate -= PR[j] * (long double) <<
989 // Eq.4.23 of of the TN 1580 // Eq.4.23 of of the TN
990 // note the negative here is requi 1581 // note the negative here is required as the rate in the
991 // equation is defined to be negat 1582 // equation is defined to be negative,
992 // i.e. decay away, but we need po << 1583 // i.e. decay away, but we need pasitive value here.
993 1584
994 // G4cout << j << "\t"<< PT[j]/s < << 1585 // G4cout << j << "\t"<< PT[j]/s <<"\t"<<PR[j]<< "\t"
>> 1586 // << decayRate << G4endl;
995 } 1587 }
996 1588
997 // At this point any negative decay <<
998 // (order 10**-30) that negative val <<
999 // errors. Set them to zero. <<
1000 if (decayRate < 0.0) decayRate = 0. <<
1001 <<
1002 // G4cout <<theDecayTime/s <<"\t"< <<
1003 // G4cout << theTrack.GetWeight() <<
1004 <<
1005 // Add isotope to the radioactivity <<
1006 // One table for each observation t <<
1007 // SetDecayBias(G4String filename) <<
1008 <<
1009 theRadioactivityTables[decayWindows <<
1010 ->AddIsotope(PZ,PA,PE,weight1 <<
1011 <<
1012 // Now calculate the statistical we 1589 // Now calculate the statistical weight
1013 // One needs to fold the source bia 1590 // One needs to fold the source bias function with the decaytime
1014 // also need to include the track w 1591 // also need to include the track weight! (F.Lei, 28/10/10)
1015 G4double weight = weight1*decayRate << 1592 G4double weight = weight1*decayRate*theTrack.GetWeight();
1016 << 1593
>> 1594 // add the isotope to the radioactivity tables
>> 1595 // G4cout <<theDecayTime/s <<"\t"<<nbin<<G4endl;
>> 1596 // G4cout << theTrack.GetWeight() <<"\t"<<weight1<<"\t"<<decayRate<< G4endl;
>> 1597 theRadioactivityTables[decayWindows[nbin-1]]->AddIsotope(PZ,PA,PE,weight);
>> 1598
1017 // decay the isotope 1599 // decay the isotope
1018 theIonTable = (G4IonTable *)(G4Part 1600 theIonTable = (G4IonTable *)(G4ParticleTable::GetParticleTable()->GetIonTable());
1019 parentNucleus = theIonTable->GetIon 1601 parentNucleus = theIonTable->GetIon(PZ,PA,PE);
1020 1602
1021 // Create a temprary products buffe 1603 // Create a temprary products buffer.
1022 // Its contents to be transfered to 1604 // Its contents to be transfered to the products at the end of the loop
1023 G4DecayProducts* tempprods = nullpt << 1605 G4DecayProducts* tempprods = 0;
1024 1606
1025 // Decide whether to apply branchin << 1607 // Decide whether to apply branching ratio bias or not
1026 if (BRBias) { 1608 if (BRBias) {
1027 G4DecayTable* decayTable = GetDec << 1609 G4DecayTable* theDecayTable = parentNucleus->GetDecayTable();
1028 G4VDecayChannel* theDecayChannel = null << 1610 ndecaych = G4int(theDecayTable->entries()*G4UniformRand());
1029 if (nullptr != decayTable) { << 1611 G4VDecayChannel* theDecayChannel = theDecayTable->GetDecayChannel(ndecaych);
1030 ndecaych = G4int(decayTable->entries( << 1612 if (theDecayChannel == 0) {
1031 theDecayChannel = decayTable->G <<
1032 } <<
1033 <<
1034 if (theDecayChannel == nullptr) { <<
1035 // Decay channel not found. 1613 // Decay channel not found.
1036 << 1614 #ifdef G4VERBOSE
1037 if (GetVerboseLevel() > 0) { << 1615 if (GetVerboseLevel()>0) {
1038 G4cout << " G4RadioactiveDeca << 1616 G4cerr << " G4RadioactiveDecay::DoIt : cannot determine decay channel ";
1039 G4cout << " for this nucleus; << 1617 G4cerr << " for this nucleus; decay as if no biasing active ";
1040 G4cout << G4endl; << 1618 G4cerr << G4endl;
1041 if (nullptr != decayTable) { << 1619 theDecayTable ->DumpInfo();
1042 } 1620 }
1043 // DHW 6 Dec 2010 - do decay as if no << 1621 #endif
1044 tempprods = DoDecay(*parentNucl << 1622 tempprods = DoDecay(*parentNucleus); // DHW 6 Dec 2010 - do decay as if no biasing
>> 1623 // to avoid deref of temppprods = 0
1045 } else { 1624 } else {
1046 // A decay channel has been ide 1625 // A decay channel has been identified, so execute the DecayIt.
1047 G4double tempmass = parentNucle << 1626 G4double tempmass = parentNucleus->GetPDGMass();
1048 tempprods = theDecayChannel->De 1627 tempprods = theDecayChannel->DecayIt(tempmass);
1049 weight *= (theDecayChannel->Get << 1628 weight *= (theDecayChannel->GetBR())*(theDecayTable->entries());
1050 } 1629 }
1051 } else { 1630 } else {
1052 tempprods = DoDecay(*parentNucleu << 1631 tempprods = DoDecay(*parentNucleus);
1053 } 1632 }
1054 1633
1055 // save the secondaries for buffers 1634 // save the secondaries for buffers
1056 numberOfSecondaries = tempprods->en 1635 numberOfSecondaries = tempprods->entries();
1057 currentTime = finalGlobalTime + the 1636 currentTime = finalGlobalTime + theDecayTime;
1058 for (index = 0; index < numberOfSec << 1637 for (index = 0; index < numberOfSecondaries; index++) {
1059 asecondaryparticle = tempprods->P 1638 asecondaryparticle = tempprods->PopProducts();
1060 if (asecondaryparticle->GetDefini << 1639 if (asecondaryparticle->GetDefinition()->GetBaryonNumber() < 5) {
1061 pw.push_back(weight); 1640 pw.push_back(weight);
1062 ptime.push_back(currentTime); 1641 ptime.push_back(currentTime);
1063 secondaryparticles.push_back(as 1642 secondaryparticles.push_back(asecondaryparticle);
1064 } 1643 }
1065 // Generate gammas and Xrays from <<
1066 else if (((const G4Ions*)(asecond <<
1067 ->GetExcitationEnergy() <<
1068 G4ParticleDefinition* apartDef <<
1069 AddDeexcitationSpectrumForBiasM <<
1070 <<
1071 } <<
1072 } 1644 }
1073 <<
1074 delete tempprods; 1645 delete tempprods;
>> 1646
1075 } // end of i loop 1647 } // end of i loop
1076 } // end of n loop 1648 } // end of n loop
1077 1649
1078 // now deal with the secondaries in the 1650 // now deal with the secondaries in the two stl containers
1079 // and submmit them back to the trackin 1651 // and submmit them back to the tracking manager
1080 totalNumberOfSecondaries = (G4int)pw.si << 1652 totalNumberOfSecondaries = pw.size();
1081 fParticleChangeForRadDecay.SetNumberOfS 1653 fParticleChangeForRadDecay.SetNumberOfSecondaries(totalNumberOfSecondaries);
1082 for (index=0; index < totalNumberOfSeco << 1654 for (index=0; index < totalNumberOfSecondaries; index++) {
1083 G4Track* secondary = new G4Track(seco 1655 G4Track* secondary = new G4Track(secondaryparticles[index],
1084 ptim 1656 ptime[index], currentPosition);
1085 secondary->SetGoodForTrackingFlag(); 1657 secondary->SetGoodForTrackingFlag();
1086 secondary->SetTouchableHandle(theTrac 1658 secondary->SetTouchableHandle(theTrack.GetTouchableHandle());
1087 secondary->SetWeight(pw[index]); 1659 secondary->SetWeight(pw[index]);
1088 fParticleChangeForRadDecay.AddSeconda 1660 fParticleChangeForRadDecay.AddSecondary(secondary);
1089 } 1661 }
>> 1662 // make sure the original track is set to stop and its kinematic energy collected
>> 1663 //
>> 1664 //theTrack.SetTrackStatus(fStopButAlive);
>> 1665 //energyDeposit += theParticle->GetKineticEnergy();
1090 1666
1091 // Kill the parent particle << 1667 } // End of Variance Reduction
1092 fParticleChangeForRadDecay.ProposeTrack <<
1093 fParticleChangeForRadDecay.ProposeLocal <<
1094 fParticleChangeForRadDecay.ProposeLocal <<
1095 // Reset NumberOfInteractionLengthLeft. <<
1096 ClearNumberOfInteractionLengthLeft(); <<
1097 } // end VR decay <<
1098 1668
1099 return &fParticleChangeForRadDecay; << 1669 // Kill the parent particle
1100 } // end of data found branch << 1670 fParticleChangeForRadDecay.ProposeTrackStatus(fStopAndKill) ;
>> 1671 fParticleChangeForRadDecay.ProposeLocalEnergyDeposit(energyDeposit);
>> 1672 fParticleChangeForRadDecay.ProposeLocalTime(finalLocalTime);
>> 1673 // Reset NumberOfInteractionLengthLeft.
>> 1674 ClearNumberOfInteractionLengthLeft();
>> 1675
>> 1676 return &fParticleChangeForRadDecay ;
>> 1677 }
1101 } 1678 }
1102 1679
1103 1680
1104 // Add gamma, X-ray, conversion and auger ele << 1681 G4DecayProducts*
1105 void << 1682 G4RadioactiveDecay::DoDecay(G4ParticleDefinition& theParticleDef)
1106 G4RadioactiveDecay::AddDeexcitationSpectrumFo << 1683 {
1107 G4dou << 1684 G4DecayProducts* products = 0;
1108 std:: << 1685
1109 std:: << 1686 // follow the decaytable and generate the secondaries...
1110 std:: << 1687 #ifdef G4VERBOSE
1111 { << 1688 if (GetVerboseLevel()>0) G4cout<<"Begin of DoDecay..."<<G4endl;
1112 G4double elevel=((const G4Ions*)(apartDef)) << 1689 #endif
1113 G4double life_time=apartDef->GetPDGLifeTime << 1690
1114 G4ITDecay* anITChannel = 0; << 1691 G4DecayTable* theDecayTable = theParticleDef.GetDecayTable();
1115 << 1692
1116 while (life_time < halflifethreshold && ele << 1693 // Choose a decay channel.
1117 decayIT->SetupDecay(apartDef); << 1694 #ifdef G4VERBOSE
1118 G4DecayProducts* pevap_products = decayIT << 1695 if (GetVerboseLevel()>0) G4cout <<"Selecte a channel..."<<G4endl;
1119 G4int nb_pevapSecondaries = pevap_product << 1696 #endif
1120 << 1697
1121 G4DynamicParticle* a_pevap_secondary = 0; << 1698 G4VDecayChannel* theDecayChannel = theDecayTable->SelectADecayChannel();
1122 G4ParticleDefinition* secDef = 0; << 1699 if (theDecayChannel == 0) {
1123 for (G4int ind = 0; ind < nb_pevapSeconda << 1700 // Decay channel not found.
1124 a_pevap_secondary= pevap_products->PopP << 1701 G4cerr <<"G4RadioactiveDecay::DoIt : can not determine decay channel";
1125 secDef = a_pevap_secondary->GetDefiniti << 1702 G4cerr <<G4endl;
1126 << 1703 theDecayTable ->DumpInfo();
1127 if (secDef->GetBaryonNumber() > 4) { << 1704 } else {
1128 elevel = ((const G4Ions*)(secDef))->G << 1705 // A decay channel has been identified, so execute the DecayIt.
1129 life_time = secDef->GetPDGLifeTime(); << 1706 #ifdef G4VERBOSE
1130 apartDef = secDef; << 1707 if (GetVerboseLevel()>1) {
1131 if (secDef->GetPDGStable() ) { << 1708 G4cerr <<"G4RadioactiveDecay::DoIt : selected decay channel addr:";
1132 weights_v.push_back(weight); << 1709 G4cerr <<theDecayChannel <<G4endl;
1133 times_v.push_back(currentTime); <<
1134 secondaries_v.push_back(a_pevap_sec <<
1135 } <<
1136 } else { <<
1137 weights_v.push_back(weight); <<
1138 times_v.push_back(currentTime); <<
1139 secondaries_v.push_back(a_pevap_secon <<
1140 } <<
1141 } 1710 }
>> 1711 #endif
>> 1712 G4double tempmass = theParticleDef.GetPDGMass();
>> 1713 products = theDecayChannel->DecayIt(tempmass);
1142 1714
1143 delete anITChannel; << 1715 // Apply directional bias if requested by user
1144 delete pevap_products; << 1716 CollimateDecay(products);
1145 } 1717 }
>> 1718
>> 1719 return products;
1146 } 1720 }
1147 1721
>> 1722
>> 1723 // Apply directional bias for "visible" daughters (e+-, gamma, n, alpha)
>> 1724
>> 1725 void G4RadioactiveDecay::CollimateDecay(G4DecayProducts* products) {
>> 1726 if (origin == forceDecayDirection) return; // No collimation requested
>> 1727 if (180.*deg == forceDecayHalfAngle) return;
>> 1728 if (0 == products || 0 == products->entries()) return;
>> 1729
>> 1730 #ifdef G4VERBOSE
>> 1731 if (GetVerboseLevel()>0) G4cout<<"Begin of CollimateDecay..."<<G4endl;
>> 1732 #endif
>> 1733
>> 1734 // Particles suitable for directional biasing (for if-blocks below)
>> 1735 static const G4ParticleDefinition* electron = G4Electron::Definition();
>> 1736 static const G4ParticleDefinition* positron = G4Positron::Definition();
>> 1737 static const G4ParticleDefinition* neutron = G4Neutron::Definition();
>> 1738 static const G4ParticleDefinition* gamma = G4Gamma::Definition();
>> 1739 static const G4ParticleDefinition* alpha = G4Alpha::Definition();
>> 1740
>> 1741 G4ThreeVector newDirection; // Re-use to avoid memory churn
>> 1742 for (G4int i=0; i<products->entries(); i++) {
>> 1743 G4DynamicParticle* daughter = (*products)[i];
>> 1744 const G4ParticleDefinition* daughterType = daughter->GetParticleDefinition();
>> 1745
>> 1746 if (daughterType == electron || daughterType == positron ||
>> 1747 daughterType == neutron || daughterType == gamma ||
>> 1748 daughterType == alpha) CollimateDecayProduct(daughter);
>> 1749 }
>> 1750 }
>> 1751
>> 1752 void G4RadioactiveDecay::CollimateDecayProduct(G4DynamicParticle* daughter) {
>> 1753 #ifdef G4VERBOSE
>> 1754 if (GetVerboseLevel() > 1) {
>> 1755 G4cout << "CollimateDecayProduct for daughter "
>> 1756 << daughter->GetParticleDefinition()->GetParticleName() << G4endl;
>> 1757 }
>> 1758 #endif
>> 1759
>> 1760 G4ThreeVector collimate = ChooseCollimationDirection();
>> 1761 if (origin != collimate) daughter->SetMomentumDirection(collimate);
>> 1762 }
>> 1763
>> 1764
>> 1765 // Choose random direction within collimation cone
>> 1766
>> 1767 G4ThreeVector G4RadioactiveDecay::ChooseCollimationDirection() const {
>> 1768 if (origin == forceDecayDirection) return origin; // Don't do collimation
>> 1769 if (forceDecayHalfAngle == 180.*deg) return origin;
>> 1770
>> 1771 G4ThreeVector dir = forceDecayDirection;
>> 1772
>> 1773 // Return direction offset by random throw
>> 1774 if (forceDecayHalfAngle > 0.) {
>> 1775 // Generate uniform direction around central axis
>> 1776 G4double phi = 2.*pi*G4UniformRand();
>> 1777 G4double cosMin = std::cos(forceDecayHalfAngle);
>> 1778 G4double cosTheta = (1.-cosMin)*G4UniformRand() + cosMin; // [cosMin,1.)
>> 1779
>> 1780 dir.setPhi(dir.phi()+phi);
>> 1781 dir.setTheta(dir.theta()+std::acos(cosTheta));
>> 1782 }
>> 1783
>> 1784 #ifdef G4VERBOSE
>> 1785 if (GetVerboseLevel()>1)
>> 1786 G4cout << " ChooseCollimationDirection returns " << dir << G4endl;
>> 1787 #endif
>> 1788
>> 1789 return dir;
>> 1790 }
1148 1791