| 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 ////////////////////////////////////////////////////////////////////////////////
27 // << 27 // //
28 // GEANT4 Class source file << 28 // File: G4RadioactiveDecay.cc //
29 // << 29 // Author: D.H. Wright (SLAC) //
30 // G4RadioactiveDecay << 30 // Date: 9 August 2017 //
31 // << 31 // Description: version the G4RadioactiveDecay process by F. Lei and //
32 // Author: D.H. Wright (SLAC) << 32 // P.R. Truscott with biasing and activation calculations //
33 // Date: 29 August 2017 << 33 // removed to a derived class. It performs alpha, beta, //
34 // << 34 // electron capture and isomeric transition decays of //
35 // Based on the code of F. Lei and P.R. Trusc << 35 // radioactive nuclei. //
36 // << 36 // //
37 ////////////////////////////////////////////// 37 ////////////////////////////////////////////////////////////////////////////////
38 38
39 #include "G4RadioactiveDecay.hh" 39 #include "G4RadioactiveDecay.hh"
40 #include "G4RadioactivationMessenger.hh" << 40 #include "G4RadioactiveDecayMessenger.hh"
41 41
42 #include "G4SystemOfUnits.hh" 42 #include "G4SystemOfUnits.hh"
43 #include "G4DynamicParticle.hh" 43 #include "G4DynamicParticle.hh"
44 #include "G4DecayProducts.hh" 44 #include "G4DecayProducts.hh"
45 #include "G4DecayTable.hh" 45 #include "G4DecayTable.hh"
46 #include "G4ParticleChangeForRadDecay.hh" 46 #include "G4ParticleChangeForRadDecay.hh"
47 #include "G4ITDecay.hh" 47 #include "G4ITDecay.hh"
48 #include "G4BetaDecayType.hh" 48 #include "G4BetaDecayType.hh"
49 #include "G4BetaMinusDecay.hh" 49 #include "G4BetaMinusDecay.hh"
50 #include "G4BetaPlusDecay.hh" 50 #include "G4BetaPlusDecay.hh"
51 #include "G4ECDecay.hh" 51 #include "G4ECDecay.hh"
52 #include "G4AlphaDecay.hh" 52 #include "G4AlphaDecay.hh"
53 #include "G4TritonDecay.hh" 53 #include "G4TritonDecay.hh"
54 #include "G4ProtonDecay.hh" 54 #include "G4ProtonDecay.hh"
55 #include "G4NeutronDecay.hh" 55 #include "G4NeutronDecay.hh"
56 #include "G4SFDecay.hh" 56 #include "G4SFDecay.hh"
57 #include "G4VDecayChannel.hh" 57 #include "G4VDecayChannel.hh"
58 #include "G4NuclearDecay.hh" 58 #include "G4NuclearDecay.hh"
59 #include "G4RadioactiveDecayMode.hh" 59 #include "G4RadioactiveDecayMode.hh"
60 #include "G4Fragment.hh" 60 #include "G4Fragment.hh"
61 #include "G4Ions.hh" 61 #include "G4Ions.hh"
62 #include "G4IonTable.hh" 62 #include "G4IonTable.hh"
63 #include "G4BetaDecayType.hh" 63 #include "G4BetaDecayType.hh"
64 #include "Randomize.hh" 64 #include "Randomize.hh"
65 #include "G4LogicalVolumeStore.hh" 65 #include "G4LogicalVolumeStore.hh"
66 #include "G4NuclearLevelData.hh" 66 #include "G4NuclearLevelData.hh"
67 #include "G4DeexPrecoParameters.hh" 67 #include "G4DeexPrecoParameters.hh"
68 #include "G4LevelManager.hh" 68 #include "G4LevelManager.hh"
69 #include "G4ThreeVector.hh" 69 #include "G4ThreeVector.hh"
70 #include "G4Electron.hh" 70 #include "G4Electron.hh"
71 #include "G4Positron.hh" 71 #include "G4Positron.hh"
72 #include "G4Neutron.hh" 72 #include "G4Neutron.hh"
73 #include "G4Gamma.hh" 73 #include "G4Gamma.hh"
74 #include "G4Alpha.hh" 74 #include "G4Alpha.hh"
75 #include "G4Triton.hh" 75 #include "G4Triton.hh"
76 #include "G4Proton.hh" 76 #include "G4Proton.hh"
77 77
78 #include "G4HadronicProcessType.hh" 78 #include "G4HadronicProcessType.hh"
79 #include "G4HadronicProcessStore.hh" 79 #include "G4HadronicProcessStore.hh"
80 #include "G4HadronicException.hh" 80 #include "G4HadronicException.hh"
81 #include "G4LossTableManager.hh" 81 #include "G4LossTableManager.hh"
82 #include "G4VAtomDeexcitation.hh" 82 #include "G4VAtomDeexcitation.hh"
83 #include "G4UAtomicDeexcitation.hh" 83 #include "G4UAtomicDeexcitation.hh"
84 #include "G4PhotonEvaporation.hh" 84 #include "G4PhotonEvaporation.hh"
>> 85 #include "G4HadronicParameters.hh"
85 86
86 #include <vector> 87 #include <vector>
87 #include <sstream> 88 #include <sstream>
88 #include <algorithm> 89 #include <algorithm>
89 #include <fstream> 90 #include <fstream>
90 91
>> 92 #include "G4PhysicsModelCatalog.hh"
>> 93
91 using namespace CLHEP; 94 using namespace CLHEP;
92 95
93 G4RadioactiveDecay::G4RadioactiveDecay(const G << 96 const G4double G4RadioactiveDecay::levelTolerance = 10.0*eV;
94 const G << 97 const G4ThreeVector G4RadioactiveDecay::origin(0.,0.,0.);
95 : G4VRadioactiveDecay(processName, timeThres << 98
>> 99 #ifdef G4MULTITHREADED
>> 100 #include "G4AutoLock.hh"
>> 101 G4Mutex G4RadioactiveDecay::radioactiveDecayMutex = G4MUTEX_INITIALIZER;
>> 102 DecayTableMap* G4RadioactiveDecay::master_dkmap = 0;
>> 103
>> 104 G4int& G4RadioactiveDecay::NumberOfInstances()
>> 105 {
>> 106 static G4int numberOfInstances = 0;
>> 107 return numberOfInstances;
>> 108 }
>> 109 #endif
>> 110
>> 111 G4RadioactiveDecay::G4RadioactiveDecay(const G4String& processName)
>> 112 : G4VRestDiscreteProcess(processName, fDecay), isInitialised(false),
>> 113 forceDecayDirection(0.,0.,0.), forceDecayHalfAngle(0.*deg), dirPath(""),
>> 114 verboseLevel(1),
>> 115 fThresholdForVeryLongDecayTime( 1.0e+27*CLHEP::nanosecond ) // Longer than twice Universe's age
96 { 116 {
97 #ifdef G4VERBOSE 117 #ifdef G4VERBOSE
98 if (GetVerboseLevel() > 1) { 118 if (GetVerboseLevel() > 1) {
99 G4cout << "G4RadioactiveDecay constructor: 119 G4cout << "G4RadioactiveDecay constructor: processName = " << processName
100 << G4endl; 120 << G4endl;
101 } 121 }
102 #endif 122 #endif
103 123
104 theRadioactivationMessenger = new G4Radioact << 124 SetProcessSubType(fRadioactiveDecay);
105 125
106 // Apply default values. << 126 theRadioactiveDecayMessenger = new G4RadioactiveDecayMessenger(this);
107 NSourceBin = 1; << 127 pParticleChange = &fParticleChangeForRadDecay;
108 SBin[0] = 0.* s; << 128
109 SBin[1] = 1.* s; // Convert to ns << 129 // Set up photon evaporation for use in G4ITDecay
110 SProfile[0] = 1.; << 130 photonEvaporation = new G4PhotonEvaporation();
111 SProfile[1] = 0.; << 131 photonEvaporation->RDMForced(true);
112 NDecayBin = 1; << 132 photonEvaporation->SetICM(true);
113 DBin[0] = 0. * s ; << 133
114 DBin[1] = 1. * s; << 134 // DHW G4DeexPrecoParameters* deex = G4NuclearLevelData::GetInstance()->GetParameters();
115 DProfile[0] = 1.; << 135 // DHW deex->SetCorrelatedGamma(true);
116 DProfile[1] = 0.; << 136
117 decayWindows[0] = 0; << 137 // Check data directory
118 G4RadioactivityTable* rTable = new G4Radioac << 138 char* path_var = std::getenv("G4RADIOACTIVEDATA");
119 theRadioactivityTables.push_back(rTable); << 139 if (!path_var) {
120 NSplit = 1; << 140 G4Exception("G4RadioactiveDecay()","HAD_RDM_200",FatalException,
121 AnalogueMC = true; << 141 "Environment variable G4RADIOACTIVEDATA is not set");
122 BRBias = true; << 142 } else {
123 halflifethreshold = 1000.*nanosecond; << 143 dirPath = path_var; // convert to string
>> 144 std::ostringstream os;
>> 145 os << dirPath << "/z1.a3"; // used as a dummy
>> 146 std::ifstream testFile;
>> 147 testFile.open(os.str() );
>> 148 if (!testFile.is_open() )
>> 149 G4Exception("G4RadioactiveDecay()","HAD_RDM_201",FatalException,
>> 150 "Environment variable G4RADIOACTIVEDATA is set, but does not point to correct directory");
>> 151 }
>> 152
>> 153 // Reset the list of user defined data files
>> 154 theUserRadioactiveDataFiles.clear();
>> 155
>> 156 // Instantiate the map of decay tables
>> 157 #ifdef G4MULTITHREADED
>> 158 G4AutoLock lk(&G4RadioactiveDecay::radioactiveDecayMutex);
>> 159 NumberOfInstances()++;
>> 160 if(!master_dkmap) master_dkmap = new DecayTableMap;
>> 161 #endif
>> 162 dkmap = new DecayTableMap;
>> 163
>> 164 // Apply default values
>> 165 applyARM = true;
>> 166 applyICM = true; // Always on; keep only for backward compatibility
>> 167
>> 168 // RDM applies to all logical volumes by default
>> 169 isAllVolumesMode = true;
>> 170 SelectAllVolumes();
>> 171 G4HadronicProcessStore::Instance()->RegisterExtraProcess(this);
124 } 172 }
125 173
126 void G4RadioactiveDecay::ProcessDescription(st 174 void G4RadioactiveDecay::ProcessDescription(std::ostream& outFile) const
127 { 175 {
128 outFile << "The G4RadioactiveDecay process p << 176 outFile << "The radioactive decay process (G4RadioactiveDecay) handles the\n"
129 << "nuclides (G4GenericIon) in biase << 177 << "alpha, beta+, beta-, electron capture and isomeric transition\n"
130 << "duplication, branching ratio bia << 178 << "decays of nuclei (G4GenericIon) with masses A > 4.\n"
131 << "and detector time convolution. <<
132 << "activation physics.\n" <<
133 << "The required half-lives and deca 179 << "The required half-lives and decay schemes are retrieved from\n"
134 << "the RadioactiveDecay database wh 180 << "the RadioactiveDecay database which was derived from ENSDF.\n";
135 } 181 }
136 182
137 183
138 G4RadioactiveDecay::~G4RadioactiveDecay() 184 G4RadioactiveDecay::~G4RadioactiveDecay()
139 { 185 {
140 delete theRadioactivationMessenger; << 186 delete theRadioactiveDecayMessenger;
>> 187 delete photonEvaporation;
>> 188 for (DecayTableMap::iterator i = dkmap->begin(); i != dkmap->end(); i++) {
>> 189 delete i->second;
>> 190 }
>> 191 dkmap->clear();
>> 192 delete dkmap;
>> 193 #ifdef G4MULTITHREADED
>> 194 G4AutoLock lk(&G4RadioactiveDecay::radioactiveDecayMutex);
>> 195 --NumberOfInstances();
>> 196 if(NumberOfInstances()==0)
>> 197 {
>> 198 for (DecayTableMap::iterator i = master_dkmap->begin(); i != master_dkmap->end(); i++) {
>> 199 delete i->second;
>> 200 }
>> 201 master_dkmap->clear();
>> 202 delete master_dkmap;
>> 203 }
>> 204 #endif
141 } 205 }
142 206
143 G4bool << 207
144 G4RadioactiveDecay::IsRateTableReady(const G4P << 208 G4bool G4RadioactiveDecay::IsApplicable(const G4ParticleDefinition& aParticle)
145 { 209 {
146 // Check whether the radioactive decay rates << 210 // All particles other than G4Ions, are rejected by default
147 // been calculated. << 211 if (((const G4Ions*)(&aParticle))->GetExcitationEnergy() > 0.) {return true;}
148 G4String aParticleName = aParticle.GetPartic << 212 if (aParticle.GetParticleName() == "GenericIon") {
149 for (std::size_t i = 0; i < theParentChainTa << 213 return true;
150 if (theParentChainTable[i].GetIonName() == << 214 } else if (!(aParticle.GetParticleType() == "nucleus")
>> 215 || aParticle.GetPDGLifeTime() < 0. ) {
>> 216 return false;
151 } 217 }
152 return false; << 218
>> 219 // Determine whether the nuclide falls into the correct A and Z range
>> 220 G4int A = ((const G4Ions*) (&aParticle))->GetAtomicMass();
>> 221 G4int Z = ((const G4Ions*) (&aParticle))->GetAtomicNumber();
>> 222
>> 223 if (A > theNucleusLimits.GetAMax() || A < theNucleusLimits.GetAMin())
>> 224 {return false;}
>> 225 else if (Z > theNucleusLimits.GetZMax() || Z < theNucleusLimits.GetZMin())
>> 226 {return false;}
>> 227 return true;
153 } 228 }
154 229
155 void << 230 G4DecayTable* G4RadioactiveDecay::GetDecayTable(const G4ParticleDefinition* aNucleus)
156 G4RadioactiveDecay::GetChainsFromParent(const <<
157 { 231 {
158 // Retrieve the decay rate table for the spe << 232 G4String key = aNucleus->GetParticleName();
159 G4String aParticleName = aParticle.GetPartic << 233 DecayTableMap::iterator table_ptr = dkmap->find(key);
160 234
161 for (std::size_t i = 0; i < theParentChainTa << 235 G4DecayTable* theDecayTable = 0;
162 if (theParentChainTable[i].GetIonName() == << 236 if (table_ptr == dkmap->end() ) { // If table not there,
163 theDecayRateVector = theParentChainTable << 237 theDecayTable = LoadDecayTable(*aNucleus); // load from file and
164 } << 238 if(theDecayTable) (*dkmap)[key] = theDecayTable; // store in library
165 } << 239 } else {
166 #ifdef G4VERBOSE << 240 theDecayTable = table_ptr->second;
167 if (GetVerboseLevel() > 1) { <<
168 G4cout << "The DecayRate Table for " << aP <<
169 << G4endl; <<
170 } 241 }
171 #endif << 242 return theDecayTable;
172 } 243 }
173 244
174 // ConvolveSourceTimeProfile performs the conv <<
175 // function with a single exponential characte <<
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 245
180 G4double << 246 void G4RadioactiveDecay::SelectAVolume(const G4String& aVolume)
181 G4RadioactiveDecay::ConvolveSourceTimeProfile( <<
182 { 247 {
183 G4double convolvedTime = 0.0; << 248 G4LogicalVolumeStore* theLogicalVolumes = G4LogicalVolumeStore::GetInstance();
184 G4int nbin; << 249 G4LogicalVolume* volume = nullptr;
185 if ( t > SBin[NSourceBin]) { << 250 volume = theLogicalVolumes->GetVolume(aVolume);
186 nbin = NSourceBin; << 251 if (volume != nullptr)
187 } else { << 252 {
188 nbin = 0; << 253 ValidVolumes.push_back(aVolume);
>> 254 std::sort(ValidVolumes.begin(), ValidVolumes.end());
>> 255 // sort need for performing binary_search
189 256
190 G4int loop = 0; << 257 if (GetVerboseLevel() > 0)
191 while (t > SBin[nbin]) { // Loop checking << 258 G4cout << " Radioactive decay applied to " << aVolume << G4endl;
192 loop++; << 259 }
193 if (loop > 1000) { << 260 else
194 G4Exception("G4RadioactiveDecay::Convo << 261 {
195 "HAD_RDM_100", JustWarning << 262 G4ExceptionDescription ed;
196 break; << 263 ed << aVolume << " is not a valid logical volume name."
197 } << 264 << " Decay not activated for it."
198 ++nbin; << 265 << G4endl;
199 } << 266 G4Exception("G4RadioactiveDecay::SelectAVolume()", "HAD_RDM_300",
200 --nbin; << 267 JustWarning, ed);
201 } 268 }
>> 269 }
202 270
203 // Use expm1 wherever possible to avoid larg << 271
204 // 1 - exp(x) for small x << 272 void G4RadioactiveDecay::DeselectAVolume(const G4String& aVolume)
205 G4double earg = 0.0; << 273 {
206 if (nbin > 0) { << 274 G4LogicalVolumeStore* theLogicalVolumes = G4LogicalVolumeStore::GetInstance();
207 for (G4int i = 0; i < nbin; ++i) { << 275 G4LogicalVolume* volume = nullptr;
208 earg = (SBin[i+1] - SBin[i])/tau; << 276 volume = theLogicalVolumes->GetVolume(aVolume);
209 if (earg < 100.) { << 277 if (volume != nullptr)
210 convolvedTime += SProfile[i] * std::ex << 278 {
211 std::expm1(earg); << 279 auto location= std::find(ValidVolumes.cbegin(),ValidVolumes.cend(),aVolume);
212 } else { << 280 if (location != ValidVolumes.cend() )
213 convolvedTime += SProfile[i] * << 281 {
214 (std::exp(-(t-SBin[i+1])/tau)-std::e << 282 ValidVolumes.erase(location);
215 } << 283 std::sort(ValidVolumes.begin(), ValidVolumes.end());
>> 284 isAllVolumesMode = false;
>> 285 if (GetVerboseLevel() > 0)
>> 286 G4cout << " G4RadioactiveDecay::DeselectAVolume: " << aVolume
>> 287 << " is removed from list " << G4endl;
>> 288 }
>> 289 else
>> 290 {
>> 291 G4ExceptionDescription ed;
>> 292 ed << aVolume << " is not in the list. No action taken." << G4endl;
>> 293 G4Exception("G4RadioactiveDecay::DeselectAVolume()", "HAD_RDM_300",
>> 294 JustWarning, ed);
216 } 295 }
217 } 296 }
218 convolvedTime -= SProfile[nbin] * std::expm1 << 297 else
219 // tau divided out of final result to provid << 298 {
>> 299 G4ExceptionDescription ed;
>> 300 ed << aVolume << " is not a valid logical volume name. No action taken."
>> 301 << G4endl;
>> 302 G4Exception("G4RadioactiveDecay::DeselectAVolume()", "HAD_RDM_300",
>> 303 JustWarning, ed);
>> 304 }
>> 305 }
>> 306
220 307
221 if (convolvedTime < 0.) { << 308 void G4RadioactiveDecay::SelectAllVolumes()
222 G4cout << " Convolved time =: " << convolv << 309 {
223 G4cout << " t = " << t << " tau = " << tau << 310 G4LogicalVolumeStore* theLogicalVolumes = G4LogicalVolumeStore::GetInstance();
224 G4cout << SBin[nbin] << " " << SBin[0] << << 311 G4LogicalVolume* volume = nullptr;
225 convolvedTime = 0.; << 312 ValidVolumes.clear();
>> 313 #ifdef G4VERBOSE
>> 314 if (GetVerboseLevel()>1)
>> 315 G4cout << " RDM Applies to all Volumes" << G4endl;
>> 316 #endif
>> 317 for (std::size_t i = 0; i < theLogicalVolumes->size(); ++i){
>> 318 volume = (*theLogicalVolumes)[i];
>> 319 ValidVolumes.push_back(volume->GetName());
>> 320 #ifdef G4VERBOSE
>> 321 if (GetVerboseLevel()>1)
>> 322 G4cout << " RDM Applies to Volume " << volume->GetName() << G4endl;
>> 323 #endif
226 } 324 }
>> 325 std::sort(ValidVolumes.begin(), ValidVolumes.end());
>> 326 // sort needed in order to allow binary_search
>> 327 isAllVolumesMode=true;
>> 328 }
>> 329
>> 330
>> 331 void G4RadioactiveDecay::DeselectAllVolumes()
>> 332 {
>> 333 ValidVolumes.clear();
>> 334 isAllVolumesMode=false;
227 #ifdef G4VERBOSE 335 #ifdef G4VERBOSE
228 if (GetVerboseLevel() > 2) << 336 if (GetVerboseLevel() > 1) G4cout << "RDM removed from all volumes" << G4endl;
229 G4cout << " Convolved time: " << convolved <<
230 #endif 337 #endif
231 return convolvedTime; <<
232 } 338 }
233 339
234 340
235 ////////////////////////////////////////////// 341 ////////////////////////////////////////////////////////////////////////////////
236 // 342 // //
237 // GetDecayTime << 343 // GetMeanLifeTime (required by the base class) //
238 // Randomly select a decay time for the dec <<
239 // supplied decay time bias scheme. <<
240 // 344 // //
241 ////////////////////////////////////////////// 345 ////////////////////////////////////////////////////////////////////////////////
242 346
243 G4double G4RadioactiveDecay::GetDecayTime() << 347 G4double G4RadioactiveDecay::GetMeanLifeTime(const G4Track& theTrack,
>> 348 G4ForceCondition*)
244 { 349 {
245 G4double decaytime = 0.; << 350 G4double meanlife = 0.;
246 G4double rand = G4UniformRand(); << 351 const G4DynamicParticle* theParticle = theTrack.GetDynamicParticle();
247 G4int i = 0; << 352 const G4ParticleDefinition* theParticleDef = theParticle->GetDefinition();
248 << 353 G4double theLife = theParticleDef->GetPDGLifeTime();
249 G4int loop = 0; << 354 #ifdef G4VERBOSE
250 while (DProfile[i] < rand) { /* Loop checki << 355 if (GetVerboseLevel() > 2) {
251 // Entries in DProfile[i] are all between << 356 G4cout << "G4RadioactiveDecay::GetMeanLifeTime() " << G4endl;
252 // Comparison with rand chooses which time << 357 G4cout << "KineticEnergy: " << theParticle->GetKineticEnergy()/GeV
253 ++i; << 358 << " GeV, Mass: " << theParticle->GetMass()/GeV
254 loop++; << 359 << " GeV, Life time: " << theLife/ns << " ns " << G4endl;
255 if (loop > 100000) { <<
256 G4Exception("G4RadioactiveDecay::GetDeca <<
257 JustWarning, "While loop cou <<
258 break; <<
259 } <<
260 } 360 }
261 << 361 #endif
262 rand = G4UniformRand(); << 362 if (theParticleDef->GetPDGStable()) {meanlife = DBL_MAX;}
263 decaytime = DBin[i] + rand*(DBin[i+1]-DBin[i << 363 else if (theLife < 0.0) {meanlife = DBL_MAX;}
>> 364 else {meanlife = theLife;}
>> 365 // Set meanlife to zero for excited istopes which are not in the
>> 366 // RDM database
>> 367 if (((const G4Ions*)(theParticleDef))->GetExcitationEnergy() > 0. &&
>> 368 meanlife == DBL_MAX) {meanlife = 0.;}
264 #ifdef G4VERBOSE 369 #ifdef G4VERBOSE
265 if (GetVerboseLevel() > 2) 370 if (GetVerboseLevel() > 2)
266 G4cout <<" Decay time: " <<decaytime/s <<" << 371 G4cout << " mean life time: " << meanlife/s << " s " << G4endl;
267 #endif 372 #endif
268 return decaytime; << 373
>> 374 return meanlife;
269 } 375 }
270 376
>> 377 ////////////////////////////////////////////////////////////////////////////////
>> 378 // //
>> 379 // GetMeanFreePath for decay in flight //
>> 380 // //
>> 381 ////////////////////////////////////////////////////////////////////////////////
271 382
272 G4int G4RadioactiveDecay::GetDecayTimeBin(cons << 383 G4double G4RadioactiveDecay::GetMeanFreePath(const G4Track& aTrack, G4double,
>> 384 G4ForceCondition*)
273 { 385 {
274 G4int i = 0; << 386 const G4DynamicParticle* aParticle = aTrack.GetDynamicParticle();
>> 387 const G4ParticleDefinition* aParticleDef = aParticle->GetDefinition();
>> 388 G4double tau = aParticleDef->GetPDGLifeTime();
>> 389 G4double aMass = aParticle->GetMass();
>> 390
>> 391 #ifdef G4VERBOSE
>> 392 if (GetVerboseLevel() > 2) {
>> 393 G4cout << "G4RadioactiveDecay::GetMeanFreePath() " << G4endl;
>> 394 G4cout << " KineticEnergy: " << aParticle->GetKineticEnergy()/GeV
>> 395 << " GeV, Mass: " << aMass/GeV << " GeV, tau: " << tau << " ns "
>> 396 << G4endl;
>> 397 }
>> 398 #endif
>> 399 G4double pathlength = DBL_MAX;
>> 400 if (tau != -1) {
>> 401 // Ion can decay
>> 402
>> 403 if (tau < -1000.0) {
>> 404 pathlength = DBL_MIN; // nuclide had very short lifetime or wasn't in table
>> 405
>> 406 } else if (tau < 0.0) {
>> 407 G4cout << aParticleDef->GetParticleName() << " has lifetime " << tau << G4endl;
>> 408 G4ExceptionDescription ed;
>> 409 ed << "Ion has negative lifetime " << tau
>> 410 << " but is not stable. Setting mean free path to DBL_MAX" << G4endl;
>> 411 G4Exception("G4RadioactiveDecay::GetMeanFreePath()", "HAD_RDM_011",
>> 412 JustWarning, ed);
>> 413 pathlength = DBL_MAX;
275 414
276 G4int loop = 0; << 415 } else {
277 while (aDecayTime > DBin[i] ) { /* Loop ch << 416 // Calculate mean free path
278 ++i; << 417 G4double betaGamma = aParticle->GetTotalMomentum()/aMass;
279 loop++; << 418 pathlength = c_light*tau*betaGamma;
280 if (loop > 100000) { << 419
281 G4Exception("G4RadioactiveDecay::GetDeca << 420 if (pathlength < DBL_MIN) {
282 JustWarning, "While loop cou << 421 pathlength = DBL_MIN;
283 break; << 422 #ifdef G4VERBOSE
>> 423 if (GetVerboseLevel() > 2) {
>> 424 G4cout << "G4Decay::GetMeanFreePath: "
>> 425 << aParticleDef->GetParticleName()
>> 426 << " stops, kinetic energy = "
>> 427 << aParticle->GetKineticEnergy()/keV <<" keV " << G4endl;
>> 428 }
>> 429 #endif
>> 430 }
284 } 431 }
285 } 432 }
286 433
287 return i; << 434 #ifdef G4VERBOSE
>> 435 if (GetVerboseLevel() > 2) {
>> 436 G4cout << "mean free path: "<< pathlength/m << " m" << G4endl;
>> 437 }
>> 438 #endif
>> 439 return pathlength;
288 } 440 }
289 441
290 ////////////////////////////////////////////// 442 ////////////////////////////////////////////////////////////////////////////////
291 // 443 // //
292 // GetMeanLifeTime (required by the base clas << 444 // BuildPhysicsTable - initialization of atomic de-excitation //
293 // 445 // //
294 ////////////////////////////////////////////// 446 ////////////////////////////////////////////////////////////////////////////////
295 447
296 G4double G4RadioactiveDecay::GetMeanLifeTime(c << 448 void G4RadioactiveDecay::BuildPhysicsTable(const G4ParticleDefinition&)
297 G4 <<
298 { 449 {
299 // For variance reduction time is set to 0 s << 450 if (!isInitialised) {
300 // to decay immediately. << 451 isInitialised = true;
301 // In analogue mode it returns the particle' << 452 #ifdef G4VERBOSE
302 G4double meanlife = 0.; << 453 if(G4HadronicParameters::Instance()->GetVerboseLevel() > 0 &&
303 if (AnalogueMC) meanlife = G4VRadioactiveDec << 454 G4Threading::IsMasterThread()) { StreamInfo(G4cout, "\n"); }
304 return meanlife; << 455 #endif
>> 456 }
>> 457 G4HadronicProcessStore::
>> 458 Instance()->RegisterParticleForExtraProcess(this,G4GenericIon::GenericIon());
305 } 459 }
306 460
>> 461 ////////////////////////////////////////////////////////////////////////////////
>> 462 // //
>> 463 // StreamInfo - stream out parameters //
>> 464 // //
>> 465 ////////////////////////////////////////////////////////////////////////////////
307 466
308 void 467 void
309 G4RadioactiveDecay::SetDecayRate(G4int theZ, G << 468 G4RadioactiveDecay::StreamInfo(std::ostream& os, const G4String& endline)
310 G4int theG, std::vector<G4double>& the << 469 {
311 std::vector<G4double>& theTaos) << 470 G4DeexPrecoParameters* deex =
312 // Why not make this a method of G4Radioactiv << 471 G4NuclearLevelData::GetInstance()->GetParameters();
313 { << 472 G4EmParameters* emparam = G4EmParameters::Instance();
314 //fill the decay rate vector << 473
315 ratesToDaughter.SetZ(theZ); << 474 G4int prec = os.precision(5);
316 ratesToDaughter.SetA(theA); << 475 os << "======================================================================"
317 ratesToDaughter.SetE(theE); << 476 << endline;
318 ratesToDaughter.SetGeneration(theG); << 477 os << "====== Radioactive Decay Physics Parameters ======="
319 ratesToDaughter.SetDecayRateC(theCoefficient << 478 << endline;
320 ratesToDaughter.SetTaos(theTaos); << 479 os << "======================================================================"
>> 480 << endline;
>> 481 os << "Max life time "
>> 482 << deex->GetMaxLifeTime()/CLHEP::ps << " ps" << endline;
>> 483 os << "Internal e- conversion flag "
>> 484 << deex->GetInternalConversionFlag() << endline;
>> 485 os << "Stored internal conversion coefficients "
>> 486 << deex->StoreICLevelData() << endline;
>> 487 os << "Enable correlated gamma emission "
>> 488 << deex->CorrelatedGamma() << endline;
>> 489 os << "Max 2J for sampling of angular correlations "
>> 490 << deex->GetTwoJMAX() << endline;
>> 491 os << "Atomic de-excitation enabled "
>> 492 << emparam->Fluo() << endline;
>> 493 os << "Auger electron emission enabled "
>> 494 << emparam->Auger() << endline;
>> 495 os << "Check EM cuts disabled for atomic de-excitation "
>> 496 << emparam->DeexcitationIgnoreCut() << endline;
>> 497 os << "Use Bearden atomic level energies "
>> 498 << emparam->BeardenFluoDir() << endline;
>> 499 os << "Use ANSTO fluorescence model "
>> 500 << emparam->ANSTOFluoDir() << endline;
>> 501 os << "Threshold for very long decay time at rest "
>> 502 << fThresholdForVeryLongDecayTime/CLHEP::ns << " ns" << endline;
>> 503 os << "======================================================================"
>> 504 << G4endl;
>> 505 os.precision(prec);
321 } 506 }
322 507
323 508
324 void G4RadioactiveDecay:: << 509 ////////////////////////////////////////////////////////////////////////////////
325 CalculateChainsFromParent(const G4ParticleDefi << 510 // //
>> 511 // LoadDecayTable loads the decay scheme from the RadioactiveDecay database //
>> 512 // for the parent nucleus. //
>> 513 // //
>> 514 ////////////////////////////////////////////////////////////////////////////////
>> 515
>> 516 G4DecayTable*
>> 517 G4RadioactiveDecay::LoadDecayTable(const G4ParticleDefinition& theParentNucleus)
326 { 518 {
327 // Use extended Bateman equation to calculat << 519 // Generate input data file name using Z and A of the parent nucleus
328 // progeny of theParentNucleus. The coeffic << 520 // file containing radioactive decay data.
329 // calculated using the method of P. Truscot << 521 G4int A = ((const G4Ions*)(&theParentNucleus))->GetAtomicMass();
330 // DERA Technical Note DERA/CIS/CIS2/7/36/4/ << 522 G4int Z = ((const G4Ions*)(&theParentNucleus))->GetAtomicNumber();
331 // Coefficients are then added to the decay << 523
332 << 524 G4double levelEnergy = ((const G4Ions*)(&theParentNucleus))->GetExcitationEnergy();
333 // Create and initialise variables used in t << 525 G4Ions::G4FloatLevelBase floatingLevel =
334 theDecayRateVector.clear(); << 526 ((const G4Ions*)(&theParentNucleus))->GetFloatLevelBase();
335 << 527
336 G4int nGeneration = 0; << 528 #ifdef G4MULTITHREADED
337 << 529 G4AutoLock lk(&G4RadioactiveDecay::radioactiveDecayMutex);
338 std::vector<G4double> taos; << 530
339 << 531 G4String key = theParentNucleus.GetParticleName();
340 // Dimensionless A coefficients of Eqs. 4.24 << 532 DecayTableMap::iterator master_table_ptr = master_dkmap->find(key);
341 std::vector<G4double> Acoeffs; << 533
342 << 534 if (master_table_ptr != master_dkmap->end() ) { // If table is there
343 // According to Eq. 4.26 the first coefficie << 535 return master_table_ptr->second;
344 Acoeffs.push_back(-1.); << 536 }
345 << 537 #endif
346 const G4Ions* ion = static_cast<const G4Ions << 538
347 G4int A = ion->GetAtomicMass(); << 539 //Check if data have been provided by the user
348 G4int Z = ion->GetAtomicNumber(); << 540 G4String file = theUserRadioactiveDataFiles[1000*A+Z];
349 G4double E = ion->GetExcitationEnergy(); << 541
350 G4double tao = ion->GetPDGLifeTime(); << 542 if (file == "") {
351 if (tao < 0.) tao = 1e-100; << 543 std::ostringstream os;
352 taos.push_back(tao); << 544 os << dirPath << "/z" << Z << ".a" << A << '\0';
353 G4int nEntry = 0; << 545 file = os.str();
354 << 546 }
355 // Fill the decay rate container (G4Radioact << 547
356 // isotope data << 548 G4DecayTable* theDecayTable = new G4DecayTable();
357 SetDecayRate(Z,A,E,nGeneration,Acoeffs,taos) << 549 G4bool found(false); // True if energy level matches one in table
358 << 550
359 // store the decay rate in decay rate vector << 551 std::ifstream DecaySchemeFile;
360 theDecayRateVector.push_back(ratesToDaughter << 552 DecaySchemeFile.open(file);
361 ++nEntry; << 553
362 << 554 if (DecaySchemeFile.good()) {
363 // Now start treating the secondary generati << 555 // Initialize variables used for reading in radioactive decay data
364 G4bool stable = false; << 556 G4bool floatMatch(false);
365 G4int j; << 557 const G4int nMode = G4RadioactiveDecayModeSize;
366 G4VDecayChannel* theChannel = 0; << 558 G4double modeTotalBR[nMode] = {0.0};
367 G4NuclearDecay* theNuclearDecayChannel = 0; << 559 G4double modeSumBR[nMode];
368 << 560 for (G4int i = 0; i < nMode; i++) {
369 G4ITDecay* theITChannel = 0; << 561 modeSumBR[i] = 0.0;
370 G4BetaMinusDecay* theBetaMinusChannel = 0; <<
371 G4BetaPlusDecay* theBetaPlusChannel = 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; <<
379 G4double theBR = 0.0; <<
380 G4int AP = 0; <<
381 G4int ZP = 0; <<
382 G4int AD = 0; <<
383 G4int ZD = 0; <<
384 G4double EP = 0.; <<
385 std::vector<G4double> TP; <<
386 std::vector<G4double> RP; // A coefficient <<
387 G4ParticleDefinition *theDaughterNucleus; <<
388 G4double daughterExcitation; <<
389 G4double nearestEnergy = 0.0; <<
390 G4int nearestLevelIndex = 0; <<
391 G4ParticleDefinition *aParentNucleus; <<
392 G4IonTable* theIonTable; <<
393 G4DecayTable* parentDecayTable; <<
394 G4double theRate; <<
395 G4double TaoPlus; <<
396 G4int nS = 0; // Running index of fir <<
397 G4int nT = nEntry; // Total number of deca <<
398 const G4int nMode = G4RadioactiveDecayModeSi <<
399 G4double brs[nMode]; <<
400 // <<
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 } 562 }
411 nGeneration++; << 563
412 for (j = nS; j < nT; ++j) { << 564 char inputChars[120]={' '};
413 // First time through, get data for pare << 565 G4String inputLine;
414 ZP = theDecayRateVector[j].GetZ(); << 566 G4String recordType("");
415 AP = theDecayRateVector[j].GetA(); << 567 G4String floatingFlag("");
416 EP = theDecayRateVector[j].GetE(); << 568 G4String daughterFloatFlag("");
417 RP = theDecayRateVector[j].GetDecayRateC << 569 G4Ions::G4FloatLevelBase daughterFloatLevel;
418 TP = theDecayRateVector[j].GetTaos(); << 570 G4RadioactiveDecayMode theDecayMode;
419 if (GetVerboseLevel() > 1) { << 571 G4double decayModeTotal(0.0);
420 G4cout << "G4RadioactiveDecay::Calcula << 572 G4double parentExcitation(0.0);
421 << ZP << ", " << AP << ", " << << 573 G4double a(0.0);
422 << ") are being calculated, ge << 574 G4double b(0.0);
423 << G4endl; << 575 G4double c(0.0);
>> 576 G4double dummy(0.0);
>> 577 G4BetaDecayType betaType(allowed);
>> 578
>> 579 // Loop through each data file record until you identify the decay
>> 580 // data relating to the nuclide of concern.
>> 581
>> 582 G4bool complete(false); // bool insures only one set of values read for any
>> 583 // given parent energy level
>> 584 G4int loop = 0;
>> 585 while (!complete && !DecaySchemeFile.getline(inputChars, 120).eof()) { /* Loop checking, 01.09.2015, D.Wright */
>> 586 loop++;
>> 587 if (loop > 100000) {
>> 588 G4Exception("G4RadioactiveDecay::LoadDecayTable()", "HAD_RDM_100",
>> 589 JustWarning, "While loop count exceeded");
>> 590 break;
424 } 591 }
>> 592
>> 593 inputLine = inputChars;
>> 594 G4StrUtil::rstrip(inputLine);
>> 595 if (inputChars[0] != '#' && inputLine.length() != 0) {
>> 596 std::istringstream tmpStream(inputLine);
>> 597
>> 598 if (inputChars[0] == 'P') {
>> 599 // Nucleus is a parent type. Check excitation level to see if it
>> 600 // matches that of theParentNucleus
>> 601 tmpStream >> recordType >> parentExcitation >> floatingFlag >> dummy;
>> 602 // "dummy" takes the place of half-life
>> 603 // Now read in from ENSDFSTATE in particle category
425 604
426 aParentNucleus = theIonTable->GetIon(ZP, << 605 if (found) {
427 parentDecayTable = GetDecayTable(aParent << 606 complete = true;
428 if (nullptr == parentDecayTable) { conti <<
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 { 607 } else {
471 brs[theDecayMode] += theChannel->G << 608 // Take first level which matches excitation energy regardless of floating level
>> 609 found = (std::abs(parentExcitation*keV - levelEnergy) < levelTolerance);
>> 610 if (floatingLevel != noFloat) {
>> 611 // If floating level specificed, require match of both energy and floating level
>> 612 floatMatch = (floatingLevel == G4Ions::FloatLevelBase(floatingFlag.back()) );
>> 613 if (!floatMatch) found = false;
>> 614 }
472 } 615 }
473 } else { <<
474 brs[theDecayMode] += theChannel->Get <<
475 } <<
476 616
477 } // Combine decay channels (loop i) << 617 } else if (found) {
>> 618 // The right part of the radioactive decay data file has been found. Search
>> 619 // through it to determine the mode of decay of the subsequent records.
>> 620
>> 621 // Store for later the total decay probability for each decay mode
>> 622 if (inputLine.length() < 72) {
>> 623 tmpStream >> theDecayMode >> dummy >> decayModeTotal;
>> 624 switch (theDecayMode) {
>> 625 case IT:
>> 626 {
>> 627 G4ITDecay* anITChannel = new G4ITDecay(&theParentNucleus, decayModeTotal,
>> 628 0.0, 0.0, photonEvaporation);
>> 629 // anITChannel->SetHLThreshold(halflifethreshold);
>> 630 anITChannel->SetARM(applyARM);
>> 631 theDecayTable->Insert(anITChannel);
>> 632 // anITChannel->DumpNuclearInfo();
>> 633 }
>> 634 break;
>> 635 case BetaMinus:
>> 636 modeTotalBR[BetaMinus] = decayModeTotal; break;
>> 637 case BetaPlus:
>> 638 modeTotalBR[BetaPlus] = decayModeTotal; break;
>> 639 case KshellEC:
>> 640 modeTotalBR[KshellEC] = decayModeTotal; break;
>> 641 case LshellEC:
>> 642 modeTotalBR[LshellEC] = decayModeTotal; break;
>> 643 case MshellEC:
>> 644 modeTotalBR[MshellEC] = decayModeTotal; break;
>> 645 case NshellEC:
>> 646 modeTotalBR[NshellEC] = decayModeTotal; break;
>> 647 case Alpha:
>> 648 modeTotalBR[Alpha] = decayModeTotal; break;
>> 649 case Proton:
>> 650 modeTotalBR[Proton] = decayModeTotal; break;
>> 651 case Neutron:
>> 652 modeTotalBR[Neutron] = decayModeTotal; break;
>> 653 case SpFission:
>> 654 modeTotalBR[SpFission] = decayModeTotal; break;
>> 655 case BDProton:
>> 656 /* Not yet implemented */ break;
>> 657 case BDNeutron:
>> 658 /* Not yet implemented */ break;
>> 659 case Beta2Minus:
>> 660 /* Not yet implemented */ break;
>> 661 case Beta2Plus:
>> 662 /* Not yet implemented */ break;
>> 663 case Proton2:
>> 664 /* Not yet implemented */ break;
>> 665 case Neutron2:
>> 666 /* Not yet implemented */ break;
>> 667 case Triton:
>> 668 modeTotalBR[Triton] = decayModeTotal; break;
>> 669 case RDM_ERROR:
>> 670
>> 671 default:
>> 672 G4Exception("G4RadioactiveDecay::LoadDecayTable()", "HAD_RDM_000",
>> 673 FatalException, "Selected decay mode does not exist");
>> 674 } // switch
478 675
479 brs[BetaPlus] = brs[BetaPlus]+brs[Kshell << 676 } else {
480 brs[KshellEC] = brs[LshellEC] = brs[Mshe << 677 if (inputLine.length() < 84) {
481 for (G4int i = 0; i < nMode; ++i) { << 678 tmpStream >> theDecayMode >> a >> daughterFloatFlag >> b >> c;
482 if (brs[i] > 0.) { << 679 betaType = allowed;
483 switch (i) { << 680 } else {
484 case IT: << 681 tmpStream >> theDecayMode >> a >> daughterFloatFlag >> b >> c >> betaType;
485 // Decay mode is isomeric transiti << 682 }
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 683
566 default: << 684 // Allowed transitions are the default. Forbidden transitions are
567 break; << 685 // indicated in the last column.
568 } << 686 a /= 1000.;
569 } << 687 c /= 1000.;
570 } << 688 b /= 100.;
>> 689 daughterFloatLevel = G4Ions::FloatLevelBase(daughterFloatFlag.back());
>> 690
>> 691 switch (theDecayMode) {
>> 692 case BetaMinus:
>> 693 {
>> 694 G4BetaMinusDecay* aBetaMinusChannel =
>> 695 new G4BetaMinusDecay(&theParentNucleus, b, c*MeV, a*MeV,
>> 696 daughterFloatLevel, betaType);
>> 697 // aBetaMinusChannel->DumpNuclearInfo();
>> 698 // aBetaMinusChannel->SetHLThreshold(halflifethreshold);
>> 699 theDecayTable->Insert(aBetaMinusChannel);
>> 700 modeSumBR[BetaMinus] += b;
>> 701 }
>> 702 break;
571 703
572 // loop over all branches in summedDecay << 704 case BetaPlus:
573 // << 705 {
574 for (G4int i = 0; i < summedDecayTable-> << 706 G4BetaPlusDecay* aBetaPlusChannel =
575 theChannel = summedDecayTable->GetDeca << 707 new G4BetaPlusDecay(&theParentNucleus, b, c*MeV, a*MeV,
576 theNuclearDecayChannel = static_cast<G << 708 daughterFloatLevel, betaType);
577 theBR = theChannel->GetBR(); << 709 // aBetaPlusChannel->DumpNuclearInfo();
578 theDaughterNucleus = theNuclearDecayCh << 710 // aBetaPlusChannel->SetHLThreshold(halflifethreshold);
579 << 711 theDecayTable->Insert(aBetaPlusChannel);
580 // First check if the decay of the ori << 712 modeSumBR[BetaPlus] += b;
581 // if true create a new ground-state n <<
582 if (theNuclearDecayChannel->GetDecayMo <<
583 <<
584 A = ((const G4Ions*)(theDaughterNucl <<
585 Z = ((const G4Ions*)(theDaughterNucl <<
586 theDaughterNucleus=theIonTable->GetI <<
587 } <<
588 if (IsApplicable(*theDaughterNucleus) <<
589 aParentNucleus != theDaughterNucle <<
590 // need to make sure daughter has de <<
591 parentDecayTable = GetDecayTable(the <<
592 if (nullptr != parentDecayTable && p <<
593 A = ((const G4Ions*)(theDaughterNu <<
594 Z = ((const G4Ions*)(theDaughterNu <<
595 E = ((const G4Ions*)(theDaughterNu <<
596 <<
597 TaoPlus = theDaughterNucleus->GetP <<
598 if (TaoPlus <= 0.) TaoPlus = 1e-1 <<
599 <<
600 // first set the taos, one simply <<
601 taos.clear(); <<
602 taos = TP; // load lifetimes of <<
603 std::size_t k; <<
604 //check that TaoPlus differs from <<
605 //for (k = 0; k < TP.size(); ++k){ <<
606 //if (std::abs((TaoPlus-TP[k])/TP[ <<
607 //} <<
608 taos.push_back(TaoPlus); // add d <<
609 // now calculate the coefficiencie <<
610 // <<
611 // they are in two parts, first th <<
612 // Eq 4.24 of the TN <<
613 Acoeffs.clear(); <<
614 long double ta1,ta2; <<
615 ta2 = (long double)TaoPlus; <<
616 for (k = 0; k < RP.size(); ++k){ <<
617 ta1 = (long double)TP[k]; // <<
618 if (ta1 == ta2) { <<
619 theRate = 1.e100; <<
620 } else { <<
621 theRate = ta1/(ta1-ta2); <<
622 } 713 }
623 theRate = theRate * theBR * RP[k << 714 break;
624 Acoeffs.push_back(theRate); <<
625 } <<
626 715
627 // the second part: the n:n coeffi << 716 case KshellEC: // K-shell electron capture
628 // Eq 4.25 of the TN. Note Yn+1 i << 717 {
629 // as treated at line 1013 << 718 G4ECDecay* aKECChannel =
630 theRate = 0.; << 719 new G4ECDecay(&theParentNucleus, b, c*MeV, a*MeV,
631 long double aRate, aRate1; << 720 daughterFloatLevel, KshellEC);
632 aRate1 = 0.L; << 721 // aKECChannel->DumpNuclearInfo();
633 for (k = 0; k < RP.size(); ++k){ << 722 // aKECChannel->SetHLThreshold(halflifethreshold);
634 ta1 = (long double)TP[k]; << 723 aKECChannel->SetARM(applyARM);
635 if (ta1 == ta2 ) { << 724 theDecayTable->Insert(aKECChannel);
636 aRate = 1.e100; << 725 modeSumBR[KshellEC] += b;
637 } else { <<
638 aRate = ta2/(ta1-ta2); <<
639 } 726 }
640 aRate = aRate * (long double)(th << 727 break;
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; <<
656 nT = nEntry; <<
657 if (nS == nT) stable = true; <<
658 } // while nuclide is not stable <<
659 <<
660 // end of while loop <<
661 // the calculation completed here <<
662 <<
663 <<
664 // fill the first part of the decay rate tab <<
665 // which is the name of the original particl <<
666 chainsFromParent.SetIonName(theParentNucleus <<
667 728
668 // now fill the decay table with the newly c << 729 case LshellEC: // L-shell electron capture
669 chainsFromParent.SetItsRates(theDecayRateVec << 730 {
>> 731 G4ECDecay* aLECChannel =
>> 732 new G4ECDecay(&theParentNucleus, b, c*MeV, a*MeV,
>> 733 daughterFloatLevel, LshellEC);
>> 734 // aLECChannel->DumpNuclearInfo();
>> 735 // aLECChannel->SetHLThreshold(halflifethreshold);
>> 736 aLECChannel->SetARM(applyARM);
>> 737 theDecayTable->Insert(aLECChannel);
>> 738 modeSumBR[LshellEC] += b;
>> 739 }
>> 740 break;
670 741
671 // finally add the decayratetable to the tab << 742 case MshellEC: // M-shell electron capture
672 theParentChainTable.push_back(chainsFromPare << 743 {
673 } << 744 G4ECDecay* aMECChannel =
>> 745 new G4ECDecay(&theParentNucleus, b, c*MeV, a*MeV,
>> 746 daughterFloatLevel, MshellEC);
>> 747 // aMECChannel->DumpNuclearInfo();
>> 748 // aMECChannel->SetHLThreshold(halflifethreshold);
>> 749 aMECChannel->SetARM(applyARM);
>> 750 theDecayTable->Insert(aMECChannel);
>> 751 modeSumBR[MshellEC] += b;
>> 752 }
>> 753 break;
674 754
675 ////////////////////////////////////////////// << 755 case NshellEC: // N-shell electron capture
676 // << 756 {
677 // SetSourceTimeProfile << 757 G4ECDecay* aNECChannel =
678 // read in the source time profile functio << 758 new G4ECDecay(&theParentNucleus, b, c*MeV, a*MeV,
679 // << 759 daughterFloatLevel, NshellEC);
680 ////////////////////////////////////////////// << 760 // aNECChannel->DumpNuclearInfo();
>> 761 // aNECChannel->SetHLThreshold(halflifethreshold);
>> 762 aNECChannel->SetARM(applyARM);
>> 763 theDecayTable->Insert(aNECChannel);
>> 764 modeSumBR[NshellEC] += b;
>> 765 }
>> 766 break;
681 767
682 void G4RadioactiveDecay::SetSourceTimeProfile( << 768 case Alpha:
683 { << 769 {
684 std::ifstream infile ( filename, std::ios::i << 770 G4AlphaDecay* anAlphaChannel =
685 if (!infile) { << 771 new G4AlphaDecay(&theParentNucleus, b, c*MeV, a*MeV,
686 G4ExceptionDescription ed; << 772 daughterFloatLevel);
687 ed << " Could not open file " << filename << 773 // anAlphaChannel->DumpNuclearInfo();
688 G4Exception("G4RadioactiveDecay::SetSource << 774 // anAlphaChannel->SetHLThreshold(halflifethreshold);
689 FatalException, ed); << 775 theDecayTable->Insert(anAlphaChannel);
690 } << 776 modeSumBR[Alpha] += b;
>> 777 }
>> 778 break;
691 779
692 G4double bin, flux; << 780 case Proton:
693 NSourceBin = -1; << 781 {
>> 782 G4ProtonDecay* aProtonChannel =
>> 783 new G4ProtonDecay(&theParentNucleus, b, c*MeV, a*MeV,
>> 784 daughterFloatLevel);
>> 785 // aProtonChannel->DumpNuclearInfo();
>> 786 // aProtonChannel->SetHLThreshold(halflifethreshold);
>> 787 theDecayTable->Insert(aProtonChannel);
>> 788 modeSumBR[Proton] += b;
>> 789 }
>> 790 break;
694 791
695 G4int loop = 0; << 792 case Neutron:
696 while (infile >> bin >> flux) { /* Loop che << 793 {
697 loop++; << 794 G4NeutronDecay* aNeutronChannel =
698 if (loop > 10000) { << 795 new G4NeutronDecay(&theParentNucleus, b, c*MeV, a*MeV,
699 G4Exception("G4RadioactiveDecay::SetSour << 796 daughterFloatLevel);
700 JustWarning, "While loop cou << 797 // aNeutronChannel->DumpNuclearInfo();
701 break; << 798 // aNeutronChannel->SetHLThreshold(halflifethreshold);
702 } << 799 theDecayTable->Insert(aNeutronChannel);
703 << 800 modeSumBR[Neutron] += b;
704 NSourceBin++; << 801 }
705 if (NSourceBin > 99) { << 802 break;
706 G4Exception("G4RadioactiveDecay::SetSour <<
707 FatalException, "Input sourc <<
708 803
709 } else { << 804 case SpFission:
710 SBin[NSourceBin] = bin * s; // Conver << 805 {
711 SProfile[NSourceBin] = flux; // Dimens << 806 G4SFDecay* aSpontFissChannel =
>> 807 // new G4SFDecay(&theParentNucleus, decayModeTotal, 0.0, 0.0);
>> 808 new G4SFDecay(&theParentNucleus, b, c*MeV, a*MeV,
>> 809 daughterFloatLevel);
>> 810 theDecayTable->Insert(aSpontFissChannel);
>> 811 modeSumBR[SpFission] += b;
>> 812 }
>> 813 break;
>> 814
>> 815 case BDProton:
>> 816 // Not yet implemented
>> 817 // G4cout << " beta-delayed proton decay, a = " << a << ", b = " << b << ", c = " << c << G4endl;
>> 818 break;
>> 819
>> 820 case BDNeutron:
>> 821 // Not yet implemented
>> 822 // G4cout << " beta-delayed neutron decay, a = " << a << ", b = " << b << ", c = " << c << G4endl;
>> 823 break;
>> 824
>> 825 case Beta2Minus:
>> 826 // Not yet implemented
>> 827 // G4cout << " Double beta- decay, a = " << a << ", b = " << b << ", c = " << c << G4endl;
>> 828 break;
>> 829
>> 830 case Beta2Plus:
>> 831 // Not yet implemented
>> 832 // G4cout << " Double beta+ decay, a = " << a << ", b = " << b << ", c = " << c << G4endl;
>> 833 break;
>> 834
>> 835 case Proton2:
>> 836 // Not yet implemented
>> 837 // G4cout << " Double proton decay, a = " << a << ", b = " << b << ", c = " << c << G4endl;
>> 838 break;
>> 839
>> 840 case Neutron2:
>> 841 // Not yet implemented
>> 842 // G4cout << " Double beta- decay, a = " << a << ", b = " << b << ", c = " << c << G4endl;
>> 843 break;
>> 844
>> 845 case Triton:
>> 846 {
>> 847 G4TritonDecay* aTritonChannel =
>> 848 new G4TritonDecay(&theParentNucleus, b, c*MeV, a*MeV,
>> 849 daughterFloatLevel);
>> 850 // anAlphaChannel->DumpNuclearInfo();
>> 851 // anAlphaChannel->SetHLThreshold(halflifethreshold);
>> 852 theDecayTable->Insert(aTritonChannel);
>> 853 modeSumBR[Triton] += b;
>> 854 }
>> 855 break;
>> 856
>> 857 case RDM_ERROR:
>> 858
>> 859 default:
>> 860 G4Exception("G4RadioactiveDecay::LoadDecayTable()", "HAD_RDM_000",
>> 861 FatalException, "Selected decay mode does not exist");
>> 862 } // switch
>> 863 } // line < 72
>> 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 G4NuclearDecay* 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<G4NuclearDecay*>(theChannel);
>> 879 theDecayMode = theNuclearDecayChannel->GetDecayMode();
>> 880
>> 881 if (theDecayMode != IT) {
>> 882 theBR = theChannel->GetBR();
>> 883 theChannel->SetBR(theBR*modeTotalBR[theDecayMode]/modeSumBR[theDecayMode]);
>> 884 }
712 } 885 }
>> 886 } // decay file exists
>> 887
>> 888 DecaySchemeFile.close();
>> 889
>> 890 if (!found && levelEnergy > 0) {
>> 891 // Case where IT cascade for excited isotopes has no entries in RDM database
>> 892 // Decay mode is isomeric transition.
>> 893 G4ITDecay* anITChannel = new G4ITDecay(&theParentNucleus, 1.0, 0.0, 0.0,
>> 894 photonEvaporation);
>> 895 // anITChannel->SetHLThreshold(halflifethreshold);
>> 896 anITChannel->SetARM(applyARM);
>> 897 theDecayTable->Insert(anITChannel);
713 } 898 }
714 899
715 AnalogueMC = false; << 900 if (theDecayTable && GetVerboseLevel() > 1) {
716 infile.close(); << 901 theDecayTable->DumpInfo();
>> 902 }
717 903
718 #ifdef G4VERBOSE << 904 #ifdef G4MULTITHREADED
719 if (GetVerboseLevel() > 2) << 905 //(*master_dkmap)[key] = theDecayTable; // store in master library
720 G4cout <<" Source Timeprofile Nbin = " << <<
721 #endif 906 #endif
>> 907 return theDecayTable;
722 } 908 }
723 909
724 ////////////////////////////////////////////// << 910 void
725 // << 911 G4RadioactiveDecay::AddUserDecayDataFile(G4int Z, G4int A, G4String filename)
726 // SetDecayBiasProfile <<
727 // read in the decay bias scheme function ( <<
728 // <<
729 ////////////////////////////////////////////// <<
730 <<
731 void G4RadioactiveDecay::SetDecayBias(const G4 <<
732 { 912 {
733 std::ifstream infile(filename, std::ios::in) << 913 if (Z < 1 || A < 2) G4cout << "Z and A not valid!" << G4endl;
734 if (!infile) G4Exception("G4RadioactiveDecay <<
735 FatalException, "Un <<
736 <<
737 G4double bin, flux; <<
738 G4int dWindows = 0; <<
739 G4int i ; <<
740 <<
741 theRadioactivityTables.clear(); <<
742 <<
743 NDecayBin = -1; <<
744 <<
745 G4int loop = 0; <<
746 while (infile >> bin >> flux ) { /* Loop ch <<
747 NDecayBin++; <<
748 loop++; <<
749 if (loop > 10000) { <<
750 G4Exception("G4RadioactiveDecay::SetDeca <<
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 } <<
769 for ( i = 1; i<= NDecayBin; ++i) DProfile[i] <<
770 for ( i = 0; i<= NDecayBin; ++i) DProfile[i] <<
771 // Normalize so e <<
772 // converted to accumulated probabilities <<
773 914
774 AnalogueMC = false; << 915 std::ifstream DecaySchemeFile(filename);
775 infile.close(); << 916 if (DecaySchemeFile) {
776 << 917 G4int ID_ion = A*1000 + Z;
777 #ifdef G4VERBOSE << 918 theUserRadioactiveDataFiles[ID_ion] = filename;
778 if (GetVerboseLevel() > 2) << 919 } else {
779 G4cout <<" Decay Bias Profile Nbin = " << << 920 G4ExceptionDescription ed;
780 #endif << 921 ed << filename << " does not exist! " << G4endl;
>> 922 G4Exception("G4RadioactiveDecay::AddUserDecayDataFile()", "HAD_RDM_001",
>> 923 FatalException, ed);
>> 924 }
781 } 925 }
782 926
783 ////////////////////////////////////////////// 927 ////////////////////////////////////////////////////////////////////////////////
784 // 928 // //
785 // DecayIt 929 // DecayIt //
786 // 930 // //
787 ////////////////////////////////////////////// 931 ////////////////////////////////////////////////////////////////////////////////
788 932
789 G4VParticleChange* 933 G4VParticleChange*
790 G4RadioactiveDecay::DecayIt(const G4Track& the 934 G4RadioactiveDecay::DecayIt(const G4Track& theTrack, const G4Step&)
791 { 935 {
792 // Initialize G4ParticleChange object, get p 936 // Initialize G4ParticleChange object, get particle details and decay table
793 fParticleChangeForRadDecay.Initialize(theTra 937 fParticleChangeForRadDecay.Initialize(theTrack);
794 fParticleChangeForRadDecay.ProposeWeight(the 938 fParticleChangeForRadDecay.ProposeWeight(theTrack.GetWeight());
795 const G4DynamicParticle* theParticle = theTr 939 const G4DynamicParticle* theParticle = theTrack.GetDynamicParticle();
796 const G4ParticleDefinition* theParticleDef = 940 const G4ParticleDefinition* theParticleDef = theParticle->GetDefinition();
797 941
798 // First check whether RDM applies to the cu 942 // First check whether RDM applies to the current logical volume
799 if (!isAllVolumesMode) { 943 if (!isAllVolumesMode) {
800 if (!std::binary_search(ValidVolumes.begin 944 if (!std::binary_search(ValidVolumes.begin(), ValidVolumes.end(),
801 theTrack.GetVolume()->Get 945 theTrack.GetVolume()->GetLogicalVolume()->GetName())) {
802 #ifdef G4VERBOSE 946 #ifdef G4VERBOSE
803 if (GetVerboseLevel()>0) { << 947 if (GetVerboseLevel()>1) {
804 G4cout <<"G4RadioactiveDecay::DecayIt 948 G4cout <<"G4RadioactiveDecay::DecayIt : "
805 << theTrack.GetVolume()->GetLog 949 << theTrack.GetVolume()->GetLogicalVolume()->GetName()
806 << " is not selected for the RD 950 << " is not selected for the RDM"<< G4endl;
807 G4cout << " There are " << ValidVolume 951 G4cout << " There are " << ValidVolumes.size() << " volumes" << G4endl;
808 G4cout << " The Valid volumes are " << 952 G4cout << " The Valid volumes are " << G4endl;
809 for (std::size_t i = 0; i< ValidVolume << 953 for (size_t i = 0; i< ValidVolumes.size(); i++)
810 G4cout << Va 954 G4cout << ValidVolumes[i] << G4endl;
811 } 955 }
812 #endif 956 #endif
813 fParticleChangeForRadDecay.SetNumberOfSe 957 fParticleChangeForRadDecay.SetNumberOfSecondaries(0);
814 958
815 // Kill the parent particle. 959 // Kill the parent particle.
816 fParticleChangeForRadDecay.ProposeTrackS 960 fParticleChangeForRadDecay.ProposeTrackStatus(fStopAndKill) ;
817 fParticleChangeForRadDecay.ProposeLocalE 961 fParticleChangeForRadDecay.ProposeLocalEnergyDeposit(0.0);
818 ClearNumberOfInteractionLengthLeft(); 962 ClearNumberOfInteractionLengthLeft();
819 return &fParticleChangeForRadDecay; 963 return &fParticleChangeForRadDecay;
820 } 964 }
821 } 965 }
822 966
823 // Now check if particle is valid for RDM 967 // Now check if particle is valid for RDM
824 if (!(IsApplicable(*theParticleDef) ) ) { 968 if (!(IsApplicable(*theParticleDef) ) ) {
825 // Particle is not an ion or is outside th 969 // Particle is not an ion or is outside the nucleuslimits for decay
826 #ifdef G4VERBOSE 970 #ifdef G4VERBOSE
827 if (GetVerboseLevel() > 1) { 971 if (GetVerboseLevel() > 1) {
828 G4cout << "G4RadioactiveDecay::DecayIt : 972 G4cout << "G4RadioactiveDecay::DecayIt : "
829 << theParticleDef->GetParticleNam << 973 << theParticleDef->GetParticleName()
830 << " is not an ion or is outside << 974 << " is not an ion or is outside (Z,A) limits set for the decay. "
831 << " Set particle change accordin 975 << " Set particle change accordingly. "
832 << G4endl; 976 << G4endl;
833 } 977 }
834 #endif 978 #endif
835 fParticleChangeForRadDecay.SetNumberOfSeco 979 fParticleChangeForRadDecay.SetNumberOfSecondaries(0);
836 980
837 // Kill the parent particle 981 // Kill the parent particle
838 fParticleChangeForRadDecay.ProposeTrackSta 982 fParticleChangeForRadDecay.ProposeTrackStatus(fStopAndKill) ;
839 fParticleChangeForRadDecay.ProposeLocalEne 983 fParticleChangeForRadDecay.ProposeLocalEnergyDeposit(0.0);
840 ClearNumberOfInteractionLengthLeft(); 984 ClearNumberOfInteractionLengthLeft();
841 return &fParticleChangeForRadDecay; 985 return &fParticleChangeForRadDecay;
842 } 986 }
843 987
844 G4DecayTable* theDecayTable = GetDecayTable( 988 G4DecayTable* theDecayTable = GetDecayTable(theParticleDef);
845 989
846 if (theDecayTable == nullptr || theDecayTabl << 990 if (theDecayTable == 0 || theDecayTable->entries() == 0) {
847 // No data in the decay table. Set partic 991 // No data in the decay table. Set particle change parameters
848 // to indicate this. 992 // to indicate this.
849 #ifdef G4VERBOSE 993 #ifdef G4VERBOSE
850 if (GetVerboseLevel() > 1) { 994 if (GetVerboseLevel() > 1) {
851 G4cout << "G4RadioactiveDecay::DecayIt : 995 G4cout << "G4RadioactiveDecay::DecayIt : "
852 << "decay table not defined for " 996 << "decay table not defined for "
853 << theParticleDef->GetParticleNam << 997 << theParticleDef->GetParticleName()
854 << ". Set particle change accordi 998 << ". Set particle change accordingly. "
855 << G4endl; 999 << G4endl;
856 } 1000 }
857 #endif 1001 #endif
858 fParticleChangeForRadDecay.SetNumberOfSeco 1002 fParticleChangeForRadDecay.SetNumberOfSecondaries(0);
859 1003
860 // Kill the parent particle. 1004 // Kill the parent particle.
861 fParticleChangeForRadDecay.ProposeTrackSta 1005 fParticleChangeForRadDecay.ProposeTrackStatus(fStopAndKill) ;
862 fParticleChangeForRadDecay.ProposeLocalEne 1006 fParticleChangeForRadDecay.ProposeLocalEnergyDeposit(0.0);
863 ClearNumberOfInteractionLengthLeft(); 1007 ClearNumberOfInteractionLengthLeft();
864 return &fParticleChangeForRadDecay; 1008 return &fParticleChangeForRadDecay;
865 1009
866 } else { 1010 } else {
867 // Data found. Try to decay nucleus 1011 // Data found. Try to decay nucleus
868 if (AnalogueMC) { <<
869 G4VRadioactiveDecay::DecayAnalog(theTrac <<
870 1012
871 } else { << 1013 /*
872 // Proceed with decay using variance red << 1014 G4double energyDeposit = 0.0;
873 G4double energyDeposit = 0.0; << 1015 G4double finalGlobalTime = theTrack.GetGlobalTime();
874 G4double finalGlobalTime = theTrack.GetG << 1016 G4double finalLocalTime = theTrack.GetLocalTime();
875 G4double finalLocalTime = theTrack.GetLo << 1017 G4int index;
876 G4int index; << 1018 G4ThreeVector currentPosition;
877 G4ThreeVector currentPosition; << 1019 currentPosition = theTrack.GetPosition();
878 currentPosition = theTrack.GetPosition() << 1020
879 << 1021 G4DecayProducts* products = DoDecay(*theParticleDef);
880 G4IonTable* theIonTable; << 1022
881 G4ParticleDefinition* parentNucleus; << 1023 // If the product is the same as the input kill the track if
882 << 1024 // necessary to prevent infinite loop (11/05/10, F.Lei)
883 // Get decay chains for the given nuclid << 1025 if (products->entries() == 1) {
884 if (!IsRateTableReady(*theParticleDef)) << 1026 fParticleChangeForRadDecay.SetNumberOfSecondaries(0);
885 CalculateChainsFromParent(*theParticleDef); << 1027 fParticleChangeForRadDecay.ProposeTrackStatus(fStopAndKill);
886 GetChainsFromParent(*theParticleDef); << 1028 fParticleChangeForRadDecay.ProposeLocalEnergyDeposit(0.0);
887 << 1029 ClearNumberOfInteractionLengthLeft();
888 // Declare some of the variables require << 1030 return &fParticleChangeForRadDecay;
889 G4int PZ; << 1031 }
890 G4int PA; <<
891 G4double PE; <<
892 G4String keyName; <<
893 std::vector<G4double> PT; <<
894 std::vector<G4double> PR; <<
895 G4double taotime; <<
896 long double decayRate; <<
897 <<
898 std::size_t i; <<
899 G4int numberOfSecondaries; <<
900 G4int totalNumberOfSecondaries = 0; <<
901 G4double currentTime = 0.; <<
902 G4int ndecaych; <<
903 G4DynamicParticle* asecondaryparticle; <<
904 std::vector<G4DynamicParticle*> secondar <<
905 std::vector<G4double> pw; <<
906 std::vector<G4double> ptime; <<
907 pw.clear(); <<
908 ptime.clear(); <<
909 <<
910 // Now apply the nucleus splitting <<
911 for (G4int n = 0; n < NSplit; ++n) { <<
912 // Get the decay time following the de <<
913 // supplied by user <<
914 G4double theDecayTime = GetDecayTime() <<
915 G4int nbin = GetDecayTimeBin(theDecayT <<
916 <<
917 // calculate the first part of the wei <<
918 G4double weight1 = 1.; <<
919 if (nbin == 1) { <<
920 weight1 = 1./DProfile[nbin-1] <<
921 *(DBin[nbin]-DBin[nbin-1]) <<
922 } else if (nbin > 1) { <<
923 // Go from nbin to nbin-2 because fl <<
924 weight1 = 1./(DProfile[nbin]-DProfil <<
925 *(DBin[nbin]-DBin[nbin-1]) <<
926 // weight1 = (probability of choosin <<
927 } <<
928 // it should be calculated in seconds <<
929 weight1 /= s ; <<
930 <<
931 // loop over all the possible secondar <<
932 // the first one is itself. <<
933 for (i = 0; i < theDecayRateVector.siz <<
934 PZ = theDecayRateVector[i].GetZ(); <<
935 PA = theDecayRateVector[i].GetA(); <<
936 PE = theDecayRateVector[i].GetE(); <<
937 PT = theDecayRateVector[i].GetTaos() <<
938 PR = theDecayRateVector[i].GetDecayR <<
939 <<
940 // The array of arrays theDecayRateV <<
941 // chains of a given parent nucleus <<
942 // nuclide (Z,A,E). <<
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 <<
981 // decay products of this isotope <<
982 <<
983 // For each nuclide, calculate all t <<
984 // the parent nuclide <<
985 decayRate = 0.L; <<
986 for (G4int j = 0; j < G4int(PT.size( <<
987 taotime = ConvolveSourceTimeProfil <<
988 decayRate -= PR[j] * (long double) <<
989 // Eq.4.23 of of the TN <<
990 // note the negative here is requi <<
991 // equation is defined to be negat <<
992 // i.e. decay away, but we need po <<
993 1032
994 // G4cout << j << "\t"<< PT[j]/s < << 1033 // Get parent particle information and boost the decay products to the
995 } << 1034 // laboratory frame based on this information.
996 1035
997 // At this point any negative decay << 1036 // The Parent Energy used for the boost should be the total energy of
998 // (order 10**-30) that negative val << 1037 // the nucleus of the parent ion without the energy of the shell electrons
999 // errors. Set them to zero. << 1038 // (correction for bug 1359 by L. Desorgher)
1000 if (decayRate < 0.0) decayRate = 0. << 1039 G4double ParentEnergy = theParticle->GetKineticEnergy()
1001 << 1040 + theParticle->GetParticleDefinition()->GetPDGMass();
1002 // G4cout <<theDecayTime/s <<"\t"< << 1041 G4ThreeVector ParentDirection(theParticle->GetMomentumDirection());
1003 // G4cout << theTrack.GetWeight() << 1042
1004 << 1043 if (theTrack.GetTrackStatus() == fStopButAlive) {
1005 // Add isotope to the radioactivity << 1044 // This condition seems to be always True, further investigation is needed
1006 // One table for each observation t << 1045 // (L.Desorgher)
1007 // SetDecayBias(G4String filename) << 1046 // The particle is decayed at rest.
1008 << 1047 // since the time is still for rest particle in G4 we need to add the
1009 theRadioactivityTables[decayWindows << 1048 // additional time lapsed between the particle come to rest and the
1010 ->AddIsotope(PZ,PA,PE,weight1 << 1049 // actual decay. This time is simply sampled with the mean-life of
1011 << 1050 // the particle. But we need to protect the case PDGTime < 0.
1012 // Now calculate the statistical we << 1051 // (F.Lei 11/05/10)
1013 // One needs to fold the source bia << 1052 G4double temptime = -std::log( G4UniformRand())
1014 // also need to include the track w << 1053 *theParticleDef->GetPDGLifeTime();
1015 G4double weight = weight1*decayRate << 1054 if (temptime < 0.) temptime = 0.;
1016 << 1055 finalGlobalTime += temptime;
1017 // decay the isotope << 1056 finalLocalTime += temptime;
1018 theIonTable = (G4IonTable *)(G4Part << 1057 energyDeposit += theParticle->GetKineticEnergy();
1019 parentNucleus = theIonTable->GetIon << 1058 }
1020 << 1059 products->Boost(ParentEnergy, ParentDirection);
1021 // Create a temprary products buffe <<
1022 // Its contents to be transfered to <<
1023 G4DecayProducts* tempprods = nullpt <<
1024 <<
1025 // Decide whether to apply branchin <<
1026 if (BRBias) { <<
1027 G4DecayTable* decayTable = GetDec <<
1028 G4VDecayChannel* theDecayChannel = null <<
1029 if (nullptr != decayTable) { <<
1030 ndecaych = G4int(decayTable->entries( <<
1031 theDecayChannel = decayTable->G <<
1032 } <<
1033 <<
1034 if (theDecayChannel == nullptr) { <<
1035 // Decay channel not found. <<
1036 <<
1037 if (GetVerboseLevel() > 0) { <<
1038 G4cout << " G4RadioactiveDeca <<
1039 G4cout << " for this nucleus; <<
1040 G4cout << G4endl; <<
1041 if (nullptr != decayTable) { <<
1042 } <<
1043 // DHW 6 Dec 2010 - do decay as if no <<
1044 tempprods = DoDecay(*parentNucl <<
1045 } else { <<
1046 // A decay channel has been ide <<
1047 G4double tempmass = parentNucle <<
1048 tempprods = theDecayChannel->De <<
1049 weight *= (theDecayChannel->Get <<
1050 } <<
1051 } else { <<
1052 tempprods = DoDecay(*parentNucleu <<
1053 } <<
1054 1060
1055 // save the secondaries for buffers << 1061 // Add products in theParticleChangeForRadDecay.
1056 numberOfSecondaries = tempprods->en << 1062 G4int numberOfSecondaries = products->entries();
1057 currentTime = finalGlobalTime + the << 1063 fParticleChangeForRadDecay.SetNumberOfSecondaries(numberOfSecondaries);
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 1064
1074 delete tempprods; << 1065 #ifdef G4VERBOSE
1075 } // end of i loop << 1066 if (GetVerboseLevel()>1) {
1076 } // end of n loop << 1067 G4cout <<"G4RadioactiveDecay::DecayIt : Decay vertex :";
1077 << 1068 G4cout <<" Time: " <<finalGlobalTime/ns <<"[ns]";
1078 // now deal with the secondaries in the << 1069 G4cout <<" X:" <<(theTrack.GetPosition()).x() /cm <<"[cm]";
1079 // and submmit them back to the trackin << 1070 G4cout <<" Y:" <<(theTrack.GetPosition()).y() /cm <<"[cm]";
1080 totalNumberOfSecondaries = (G4int)pw.si << 1071 G4cout <<" Z:" <<(theTrack.GetPosition()).z() /cm <<"[cm]";
1081 fParticleChangeForRadDecay.SetNumberOfS << 1072 G4cout << G4endl;
1082 for (index=0; index < totalNumberOfSeco << 1073 G4cout <<"G4Decay::DecayIt : decay products in Lab. Frame" <<G4endl;
1083 G4Track* secondary = new G4Track(seco << 1074 products->DumpInfo();
1084 ptim << 1075 products->IsChecked();
1085 secondary->SetGoodForTrackingFlag(); << 1076 }
1086 secondary->SetTouchableHandle(theTrac << 1077 #endif
1087 secondary->SetWeight(pw[index]); << 1078 for (index=0; index < numberOfSecondaries; index++) {
1088 fParticleChangeForRadDecay.AddSeconda << 1079 G4Track* secondary = new G4Track(products->PopProducts(),
1089 } << 1080 finalGlobalTime, currentPosition);
>> 1081 secondary->SetGoodForTrackingFlag();
>> 1082 secondary->SetTouchableHandle(theTrack.GetTouchableHandle());
>> 1083 fParticleChangeForRadDecay.AddSecondary(secondary);
>> 1084 }
>> 1085 delete products;
>> 1086
>> 1087 // Kill the parent particle
>> 1088 fParticleChangeForRadDecay.ProposeTrackStatus(fStopAndKill) ;
>> 1089 fParticleChangeForRadDecay.ProposeLocalEnergyDeposit(energyDeposit);
>> 1090 fParticleChangeForRadDecay.ProposeLocalTime(finalLocalTime);
>> 1091 // Reset NumberOfInteractionLengthLeft.
>> 1092 ClearNumberOfInteractionLengthLeft();
>> 1093 */
>> 1094 // Decay without variance reduction
>> 1095 DecayAnalog(theTrack);
>> 1096 return &fParticleChangeForRadDecay ;
>> 1097 }
>> 1098 }
1090 1099
1091 // Kill the parent particle << 1100
>> 1101 void G4RadioactiveDecay::DecayAnalog(const G4Track& theTrack)
>> 1102 {
>> 1103 const G4DynamicParticle* theParticle = theTrack.GetDynamicParticle();
>> 1104 const G4ParticleDefinition* theParticleDef = theParticle->GetDefinition();
>> 1105 G4DecayProducts* products = DoDecay(*theParticleDef);
>> 1106
>> 1107 // Check if the product is the same as input and kill the track if
>> 1108 // necessary to prevent infinite loop (11/05/10, F.Lei)
>> 1109 if (products->entries() == 1) {
>> 1110 fParticleChangeForRadDecay.SetNumberOfSecondaries(0);
>> 1111 fParticleChangeForRadDecay.ProposeTrackStatus(fStopAndKill);
>> 1112 fParticleChangeForRadDecay.ProposeLocalEnergyDeposit(0.0);
>> 1113 ClearNumberOfInteractionLengthLeft();
>> 1114 delete products;
>> 1115 return;
>> 1116 }
>> 1117
>> 1118 G4double energyDeposit = 0.0;
>> 1119 G4double finalGlobalTime = theTrack.GetGlobalTime();
>> 1120 G4double finalLocalTime = theTrack.GetLocalTime();
>> 1121
>> 1122 // Get parent particle information and boost the decay products to the
>> 1123 // laboratory frame
>> 1124
>> 1125 // ParentEnergy used for the boost should be the total energy of the nucleus
>> 1126 // of the parent ion without the energy of the shell electrons
>> 1127 // (correction for bug 1359 by L. Desorgher)
>> 1128 G4double ParentEnergy = theParticle->GetKineticEnergy()
>> 1129 + theParticle->GetParticleDefinition()->GetPDGMass();
>> 1130 G4ThreeVector ParentDirection(theParticle->GetMomentumDirection());
>> 1131
>> 1132 if (theTrack.GetTrackStatus() == fStopButAlive) {
>> 1133 // this condition seems to be always True, further investigation is needed (L.Desorgher)
>> 1134
>> 1135 // The particle is decayed at rest
>> 1136 // Since the time is for the particle at rest, need to add additional time
>> 1137 // lapsed between particle coming to rest and the actual decay. This time
>> 1138 // is sampled with the mean-life of the particle. Need to protect the case
>> 1139 // PDGTime < 0. (F.Lei 11/05/10)
>> 1140 G4double temptime = -std::log(G4UniformRand() ) *
>> 1141 theParticleDef->GetPDGLifeTime();
>> 1142 if (temptime < 0.) temptime = 0.;
>> 1143 finalGlobalTime += temptime;
>> 1144 finalLocalTime += temptime;
>> 1145 energyDeposit += theParticle->GetKineticEnergy();
>> 1146
>> 1147 // Kill the parent particle, and ignore its decay, if it decays later than the
>> 1148 // threshold fThresholdForVeryLongDecayTime (whose default value corresponds
>> 1149 // to more than twice the age of the universe).
>> 1150 // This kind of cut has been introduced (in April 2021) in order to avoid to
>> 1151 // account energy depositions happening after many billions of years in
>> 1152 // ordinary materials used in calorimetry, in particular Tungsten and Lead
>> 1153 // (via their natural unstable, but very long lived, isotopes, such as
>> 1154 // W183, W180 and Pb204).
>> 1155 // Note that the cut is not on the average, mean lifetime, but on the actual
>> 1156 // sampled global decay time.
>> 1157 if ( finalGlobalTime > fThresholdForVeryLongDecayTime ) {
>> 1158 fParticleChangeForRadDecay.SetNumberOfSecondaries(0);
1092 fParticleChangeForRadDecay.ProposeTrack 1159 fParticleChangeForRadDecay.ProposeTrackStatus(fStopAndKill) ;
1093 fParticleChangeForRadDecay.ProposeLocal << 1160 fParticleChangeForRadDecay.ProposeLocalEnergyDeposit(0.0);
1094 fParticleChangeForRadDecay.ProposeLocal <<
1095 // Reset NumberOfInteractionLengthLeft. <<
1096 ClearNumberOfInteractionLengthLeft(); 1161 ClearNumberOfInteractionLengthLeft();
1097 } // end VR decay << 1162 delete products;
>> 1163 return;
>> 1164 }
>> 1165 }
>> 1166 products->Boost(ParentEnergy, ParentDirection);
1098 1167
1099 return &fParticleChangeForRadDecay; << 1168 // Add products in theParticleChangeForRadDecay.
1100 } // end of data found branch << 1169 G4int numberOfSecondaries = products->entries();
1101 } << 1170 fParticleChangeForRadDecay.SetNumberOfSecondaries(numberOfSecondaries);
1102 1171
>> 1172 if (GetVerboseLevel() > 1) {
>> 1173 G4cout << "G4RadioactiveDecay::DecayAnalog: Decay vertex :";
>> 1174 G4cout << " Time: " << finalGlobalTime/ns << "[ns]";
>> 1175 G4cout << " X:" << (theTrack.GetPosition()).x() /cm << "[cm]";
>> 1176 G4cout << " Y:" << (theTrack.GetPosition()).y() /cm << "[cm]";
>> 1177 G4cout << " Z:" << (theTrack.GetPosition()).z() /cm << "[cm]";
>> 1178 G4cout << G4endl;
>> 1179 G4cout << "G4Decay::DecayIt : decay products in Lab. Frame" << G4endl;
>> 1180 products->DumpInfo();
>> 1181 products->IsChecked();
>> 1182 }
1103 1183
1104 // Add gamma, X-ray, conversion and auger ele << 1184 const G4int modelID_forIT = G4PhysicsModelCatalog::GetModelID( "model_RDM_IT" );
1105 void << 1185 G4int modelID = modelID_forIT + 10*theRadDecayMode;
1106 G4RadioactiveDecay::AddDeexcitationSpectrumFo << 1186 const G4int modelID_forAtomicRelaxation =
1107 G4dou << 1187 G4PhysicsModelCatalog::GetModelID( "model_RDM_AtomicRelaxation" );
1108 std:: << 1188 for ( G4int index = 0; index < numberOfSecondaries; ++index ) {
1109 std:: << 1189 G4Track* secondary = new G4Track( products->PopProducts(), finalGlobalTime,
1110 std:: << 1190 theTrack.GetPosition() );
1111 { << 1191 secondary->SetWeight( theTrack.GetWeight() );
1112 G4double elevel=((const G4Ions*)(apartDef)) << 1192 secondary->SetCreatorModelID( modelID );
1113 G4double life_time=apartDef->GetPDGLifeTime << 1193 // Change for atomics relaxation
1114 G4ITDecay* anITChannel = 0; << 1194 if ( theRadDecayMode == IT && index > 0 ) {
1115 << 1195 if ( index == numberOfSecondaries-1 ) {
1116 while (life_time < halflifethreshold && ele << 1196 secondary->SetCreatorModelID( modelID_forIT );
1117 decayIT->SetupDecay(apartDef); <<
1118 G4DecayProducts* pevap_products = decayIT <<
1119 G4int nb_pevapSecondaries = pevap_product <<
1120 <<
1121 G4DynamicParticle* a_pevap_secondary = 0; <<
1122 G4ParticleDefinition* secDef = 0; <<
1123 for (G4int ind = 0; ind < nb_pevapSeconda <<
1124 a_pevap_secondary= pevap_products->PopP <<
1125 secDef = a_pevap_secondary->GetDefiniti <<
1126 <<
1127 if (secDef->GetBaryonNumber() > 4) { <<
1128 elevel = ((const G4Ions*)(secDef))->G <<
1129 life_time = secDef->GetPDGLifeTime(); <<
1130 apartDef = secDef; <<
1131 if (secDef->GetPDGStable() ) { <<
1132 weights_v.push_back(weight); <<
1133 times_v.push_back(currentTime); <<
1134 secondaries_v.push_back(a_pevap_sec <<
1135 } <<
1136 } else { 1197 } else {
1137 weights_v.push_back(weight); << 1198 secondary->SetCreatorModelID( modelID_forAtomicRelaxation) ;
1138 times_v.push_back(currentTime); <<
1139 secondaries_v.push_back(a_pevap_secon <<
1140 } 1199 }
>> 1200 } else if ( theRadDecayMode >= KshellEC && theRadDecayMode <= NshellEC &&
>> 1201 index < numberOfSecondaries-1 ) {
>> 1202 secondary->SetCreatorModelID( modelID_forAtomicRelaxation );
>> 1203 }
>> 1204 secondary->SetGoodForTrackingFlag();
>> 1205 secondary->SetTouchableHandle( theTrack.GetTouchableHandle() );
>> 1206 fParticleChangeForRadDecay.AddSecondary( secondary );
>> 1207 }
>> 1208
>> 1209 delete products;
>> 1210
>> 1211 // Kill the parent particle
>> 1212 fParticleChangeForRadDecay.ProposeTrackStatus(fStopAndKill) ;
>> 1213 fParticleChangeForRadDecay.ProposeLocalEnergyDeposit(energyDeposit);
>> 1214 fParticleChangeForRadDecay.ProposeLocalTime(finalLocalTime);
>> 1215
>> 1216 // Reset NumberOfInteractionLengthLeft.
>> 1217 ClearNumberOfInteractionLengthLeft();
>> 1218 }
>> 1219
>> 1220
>> 1221 G4DecayProducts*
>> 1222 G4RadioactiveDecay::DoDecay(const G4ParticleDefinition& theParticleDef)
>> 1223 {
>> 1224 G4DecayProducts* products = 0;
>> 1225 G4DecayTable* theDecayTable = GetDecayTable(&theParticleDef);
>> 1226
>> 1227 // Choose a decay channel.
>> 1228 // G4DecayTable::SelectADecayChannel checks to see if sum of daughter masses
>> 1229 // exceeds parent mass. Pass it the parent mass + maximum Q value to account
>> 1230 // for difference in mass defect.
>> 1231 G4double parentPlusQ = theParticleDef.GetPDGMass() + 30.*MeV;
>> 1232 G4VDecayChannel* theDecayChannel = theDecayTable->SelectADecayChannel(parentPlusQ);
>> 1233
>> 1234 if (theDecayChannel == 0) {
>> 1235 // Decay channel not found.
>> 1236 G4ExceptionDescription ed;
>> 1237 ed << " Cannot determine decay channel for " << theParticleDef.GetParticleName() << G4endl;
>> 1238 G4Exception("G4RadioactiveDecay::DoDecay", "HAD_RDM_013",
>> 1239 FatalException, ed);
>> 1240 } else {
>> 1241 // A decay channel has been identified, so execute the DecayIt.
>> 1242 #ifdef G4VERBOSE
>> 1243 if (GetVerboseLevel() > 1) {
>> 1244 G4cout << "G4RadioactiveDecay::DoIt : selected decay channel addr: "
>> 1245 << theDecayChannel << G4endl;
1141 } 1246 }
>> 1247 #endif
>> 1248 theRadDecayMode = (static_cast<G4NuclearDecay*>(theDecayChannel))->GetDecayMode();
>> 1249 products = theDecayChannel->DecayIt(theParticleDef.GetPDGMass() );
>> 1250
>> 1251 // Apply directional bias if requested by user
>> 1252 CollimateDecay(products);
>> 1253 }
>> 1254
>> 1255 return products;
>> 1256 }
>> 1257
>> 1258
>> 1259 // Apply directional bias for "visible" daughters (e+-, gamma, n, p, alpha)
>> 1260
>> 1261 void G4RadioactiveDecay::CollimateDecay(G4DecayProducts* products) {
>> 1262
>> 1263 if (origin == forceDecayDirection) return; // No collimation requested
>> 1264 if (180.*deg == forceDecayHalfAngle) return;
>> 1265 if (0 == products || 0 == products->entries()) return;
>> 1266
>> 1267 #ifdef G4VERBOSE
>> 1268 if (GetVerboseLevel() > 1) G4cout << "Begin of CollimateDecay..." << G4endl;
>> 1269 #endif
1142 1270
1143 delete anITChannel; << 1271 // Particles suitable for directional biasing (for if-blocks below)
1144 delete pevap_products; << 1272 static const G4ParticleDefinition* electron = G4Electron::Definition();
>> 1273 static const G4ParticleDefinition* positron = G4Positron::Definition();
>> 1274 static const G4ParticleDefinition* neutron = G4Neutron::Definition();
>> 1275 static const G4ParticleDefinition* gamma = G4Gamma::Definition();
>> 1276 static const G4ParticleDefinition* alpha = G4Alpha::Definition();
>> 1277 static const G4ParticleDefinition* triton = G4Triton::Definition();
>> 1278 static const G4ParticleDefinition* proton = G4Proton::Definition();
>> 1279
>> 1280 G4ThreeVector newDirection; // Re-use to avoid memory churn
>> 1281 for (G4int i=0; i<products->entries(); i++) {
>> 1282 G4DynamicParticle* daughter = (*products)[i];
>> 1283 const G4ParticleDefinition* daughterType =
>> 1284 daughter->GetParticleDefinition();
>> 1285 if (daughterType == electron || daughterType == positron ||
>> 1286 daughterType == neutron || daughterType == gamma ||
>> 1287 daughterType == alpha || daughterType == triton || daughterType == proton) CollimateDecayProduct(daughter);
1145 } 1288 }
>> 1289 }
>> 1290
>> 1291 void G4RadioactiveDecay::CollimateDecayProduct(G4DynamicParticle* daughter) {
>> 1292 #ifdef G4VERBOSE
>> 1293 if (GetVerboseLevel() > 1) {
>> 1294 G4cout << "CollimateDecayProduct for daughter "
>> 1295 << daughter->GetParticleDefinition()->GetParticleName() << G4endl;
>> 1296 }
>> 1297 #endif
>> 1298
>> 1299 G4ThreeVector collimate = ChooseCollimationDirection();
>> 1300 if (origin != collimate) daughter->SetMomentumDirection(collimate);
>> 1301 }
>> 1302
>> 1303
>> 1304 // Choose random direction within collimation cone
>> 1305
>> 1306 G4ThreeVector G4RadioactiveDecay::ChooseCollimationDirection() const {
>> 1307 if (origin == forceDecayDirection) return origin; // Don't do collimation
>> 1308 if (forceDecayHalfAngle == 180.*deg) return origin;
>> 1309
>> 1310 G4ThreeVector dir = forceDecayDirection;
>> 1311
>> 1312 // Return direction offset by random throw
>> 1313 if (forceDecayHalfAngle > 0.) {
>> 1314 // Generate uniform direction around central axis
>> 1315 G4double phi = 2.*pi*G4UniformRand();
>> 1316 G4double cosMin = std::cos(forceDecayHalfAngle);
>> 1317 G4double cosTheta = (1.-cosMin)*G4UniformRand() + cosMin; // [cosMin,1.)
>> 1318
>> 1319 dir.setPhi(dir.phi()+phi);
>> 1320 dir.setTheta(dir.theta()+std::acos(cosTheta));
>> 1321 }
>> 1322
>> 1323 #ifdef G4VERBOSE
>> 1324 if (GetVerboseLevel()>1)
>> 1325 G4cout << " ChooseCollimationDirection returns " << dir << G4endl;
>> 1326 #endif
>> 1327
>> 1328 return dir;
1146 } 1329 }
1147 1330
1148 1331