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