1 | // |
---|
2 | // ******************************************************************** |
---|
3 | // * License and Disclaimer * |
---|
4 | // * * |
---|
5 | // * The Geant4 software is copyright of the Copyright Holders of * |
---|
6 | // * the Geant4 Collaboration. It is provided under the terms and * |
---|
7 | // * conditions of the Geant4 Software License, included in the file * |
---|
8 | // * LICENSE and available at http://cern.ch/geant4/license . These * |
---|
9 | // * include a list of copyright holders. * |
---|
10 | // * * |
---|
11 | // * Neither the authors of this software system, nor their employing * |
---|
12 | // * institutes,nor the agencies providing financial support for this * |
---|
13 | // * work make any representation or warranty, express or implied, * |
---|
14 | // * regarding this software system or assume any liability for its * |
---|
15 | // * use. Please see the license in the file LICENSE and URL above * |
---|
16 | // * for the full disclaimer and the limitation of liability. * |
---|
17 | // * * |
---|
18 | // * This code implementation is the result of the scientific and * |
---|
19 | // * technical work of the GEANT4 collaboration. * |
---|
20 | // * By using, copying, modifying or distributing the software (or * |
---|
21 | // * any work based on the software) you agree to acknowledge its * |
---|
22 | // * use in resulting scientific publications, and indicate your * |
---|
23 | // * acceptance of all terms of the Geant4 Software license. * |
---|
24 | // ******************************************************************** |
---|
25 | // |
---|
26 | // %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
---|
27 | // |
---|
28 | // MODULES: G4NuclearDecayChannel.cc |
---|
29 | // |
---|
30 | // Version: 0.b.4 |
---|
31 | // Date: 14/04/00 |
---|
32 | // Author: F Lei & P R Truscott |
---|
33 | // Organisation: DERA UK |
---|
34 | // Customer: ESA/ESTEC, NOORDWIJK |
---|
35 | // Contract: 12115/96/JG/NL Work Order No. 3 |
---|
36 | // |
---|
37 | // %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
---|
38 | // |
---|
39 | // CHANGE HISTORY |
---|
40 | // -------------- |
---|
41 | // |
---|
42 | // 29 February 2000, P R Truscott, DERA UK |
---|
43 | // 0.b.3 release. |
---|
44 | // |
---|
45 | // 18 October 2002, F Lei |
---|
46 | // modified link metheds in DecayIt() to G4PhotoEvaporation() in order to |
---|
47 | // use the new Internal Coversion feature. |
---|
48 | // 13 April 2000, F Lei, DERA UK |
---|
49 | // Changes made are: |
---|
50 | // 1) Use PhotonEvaporation instead of DiscreteGammaDeexcitation |
---|
51 | // 2) verbose control |
---|
52 | // |
---|
53 | // %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
---|
54 | /////////////////////////////////////////////////////////////////////////////// |
---|
55 | // |
---|
56 | #include "G4NuclearLevelManager.hh" |
---|
57 | #include "G4NuclearLevelStore.hh" |
---|
58 | #include "G4NuclearDecayChannel.hh" |
---|
59 | #include "G4DynamicParticle.hh" |
---|
60 | #include "G4DecayProducts.hh" |
---|
61 | #include "G4DecayTable.hh" |
---|
62 | #include "G4PhysicsLogVector.hh" |
---|
63 | #include "G4ParticleChangeForRadDecay.hh" |
---|
64 | #include "G4IonTable.hh" |
---|
65 | |
---|
66 | #include "G4BetaFermiFunction.hh" |
---|
67 | #include "G4PhotonEvaporation.hh" |
---|
68 | #include "G4AtomicTransitionManager.hh" |
---|
69 | #include "G4AtomicShell.hh" |
---|
70 | #include "G4AtomicDeexcitation.hh" |
---|
71 | |
---|
72 | const G4double G4NuclearDecayChannel:: pTolerance = 0.001; |
---|
73 | const G4double G4NuclearDecayChannel:: levelTolerance = 2.0*keV; |
---|
74 | //const G4bool G4NuclearDecayChannel:: FermiOn = true; |
---|
75 | |
---|
76 | /////////////////////////////////////////////////////////////////////////////// |
---|
77 | // |
---|
78 | // |
---|
79 | // Constructor for one decay product (the nucleus). |
---|
80 | // |
---|
81 | G4NuclearDecayChannel::G4NuclearDecayChannel |
---|
82 | (const G4RadioactiveDecayMode &theMode, |
---|
83 | G4int Verbose, |
---|
84 | const G4ParticleDefinition *theParentNucleus, |
---|
85 | G4double theBR, |
---|
86 | G4double theQtransition, |
---|
87 | G4int A, |
---|
88 | G4int Z, |
---|
89 | G4double theDaughterExcitation) : |
---|
90 | G4GeneralPhaseSpaceDecay(Verbose), decayMode(theMode) |
---|
91 | { |
---|
92 | #ifdef G4VERBOSE |
---|
93 | if (GetVerboseLevel()>1) |
---|
94 | {G4cout <<"G4NuclearDecayChannel constructor for " <<G4int(theMode) <<G4endl;} |
---|
95 | #endif |
---|
96 | SetParent(theParentNucleus); |
---|
97 | FillParent(); |
---|
98 | parent_mass = theParentNucleus->GetPDGMass(); |
---|
99 | SetBR (theBR); |
---|
100 | SetNumberOfDaughters (1); |
---|
101 | FillDaughterNucleus (0, A, Z, theDaughterExcitation); |
---|
102 | Qtransition = theQtransition; |
---|
103 | halflifethreshold = 1e-6*second; |
---|
104 | applyICM = true; |
---|
105 | applyARM = true; |
---|
106 | } |
---|
107 | /////////////////////////////////////////////////////////////////////////////// |
---|
108 | // |
---|
109 | // |
---|
110 | // Constructor for a daughter nucleus and one other particle. |
---|
111 | // |
---|
112 | G4NuclearDecayChannel::G4NuclearDecayChannel |
---|
113 | (const G4RadioactiveDecayMode &theMode, |
---|
114 | G4int Verbose, |
---|
115 | const G4ParticleDefinition *theParentNucleus, |
---|
116 | G4double theBR, |
---|
117 | G4double theQtransition, |
---|
118 | G4int A, |
---|
119 | G4int Z, |
---|
120 | G4double theDaughterExcitation, |
---|
121 | const G4String theDaughterName1) : |
---|
122 | G4GeneralPhaseSpaceDecay(Verbose), decayMode(theMode) |
---|
123 | { |
---|
124 | #ifdef G4VERBOSE |
---|
125 | if (GetVerboseLevel()>1) |
---|
126 | {G4cout <<"G4NuclearDecayChannel constructor for " <<G4int(theMode) <<G4endl;} |
---|
127 | #endif |
---|
128 | SetParent (theParentNucleus); |
---|
129 | FillParent(); |
---|
130 | parent_mass = theParentNucleus->GetPDGMass(); |
---|
131 | SetBR (theBR); |
---|
132 | SetNumberOfDaughters (2); |
---|
133 | SetDaughter(0, theDaughterName1); |
---|
134 | FillDaughterNucleus (1, A, Z, theDaughterExcitation); |
---|
135 | Qtransition = theQtransition; |
---|
136 | halflifethreshold = 1e-6*second; |
---|
137 | applyICM = true; |
---|
138 | applyARM = true; |
---|
139 | } |
---|
140 | /////////////////////////////////////////////////////////////////////////////// |
---|
141 | // |
---|
142 | // |
---|
143 | // Constructor for a daughter nucleus and two other particles. |
---|
144 | // |
---|
145 | G4NuclearDecayChannel::G4NuclearDecayChannel |
---|
146 | (const G4RadioactiveDecayMode &theMode, |
---|
147 | G4int Verbose, |
---|
148 | const G4ParticleDefinition *theParentNucleus, |
---|
149 | G4double theBR, |
---|
150 | G4double theFFN, |
---|
151 | G4bool betaS, |
---|
152 | CLHEP::RandGeneral* randBeta, |
---|
153 | G4double theQtransition, |
---|
154 | G4int A, |
---|
155 | G4int Z, |
---|
156 | G4double theDaughterExcitation, |
---|
157 | const G4String theDaughterName1, |
---|
158 | const G4String theDaughterName2) : |
---|
159 | G4GeneralPhaseSpaceDecay(Verbose), decayMode(theMode) |
---|
160 | //,BetaSimple(betaS), |
---|
161 | // RandomEnergy(randBeta), Qtransition(theQtransition),FermiFN(theFFN) |
---|
162 | { |
---|
163 | #ifdef G4VERBOSE |
---|
164 | if (GetVerboseLevel()>1) |
---|
165 | {G4cout <<"G4NuclearDecayChannel constructor for " <<G4int(theMode) <<G4endl;} |
---|
166 | #endif |
---|
167 | SetParent (theParentNucleus); |
---|
168 | FillParent(); |
---|
169 | parent_mass = theParentNucleus->GetPDGMass(); |
---|
170 | SetBR (theBR); |
---|
171 | SetNumberOfDaughters (3); |
---|
172 | SetDaughter(0, theDaughterName1); |
---|
173 | SetDaughter(2, theDaughterName2); |
---|
174 | FillDaughterNucleus(1, A, Z, theDaughterExcitation); |
---|
175 | BetaSimple = betaS; |
---|
176 | RandomEnergy = randBeta; |
---|
177 | Qtransition = theQtransition; |
---|
178 | FermiFN = theFFN; |
---|
179 | halflifethreshold = 1e-6*second; |
---|
180 | applyICM = true; |
---|
181 | applyARM = true; |
---|
182 | } |
---|
183 | |
---|
184 | //////////////////////////////////////////////////////////////////////////////// |
---|
185 | // |
---|
186 | // |
---|
187 | // |
---|
188 | // |
---|
189 | #include "G4HadTmpUtil.hh" |
---|
190 | |
---|
191 | void G4NuclearDecayChannel::FillDaughterNucleus (G4int index, G4int A, G4int Z, |
---|
192 | G4double theDaughterExcitation) |
---|
193 | { |
---|
194 | // |
---|
195 | // |
---|
196 | // Determine if the proposed daughter nucleus has a sensible A, Z and excitation |
---|
197 | // energy. |
---|
198 | // |
---|
199 | if (A<1 || Z<0 || theDaughterExcitation <0.0) |
---|
200 | { |
---|
201 | G4cerr <<"Error in G4NuclearDecayChannel::FillDaughterNucleus"; |
---|
202 | G4cerr <<"Inappropriate values of daughter A, Z or excitation" <<G4endl; |
---|
203 | G4cerr <<"A = " <<A <<" and Z = " <<Z; |
---|
204 | G4cerr <<" Ex = " <<theDaughterExcitation*MeV <<"MeV" <<G4endl; |
---|
205 | G4Exception(__FILE__, G4inttostring(__LINE__), FatalException, "G4NuclearDecayChannel::FillDaughterNucleus"); |
---|
206 | } |
---|
207 | // |
---|
208 | // |
---|
209 | // Save A and Z to local variables. Find the GROUND STATE of the daughter |
---|
210 | // nucleus and save this, as an ion, in the array of daughters. |
---|
211 | // |
---|
212 | daughterA = A; |
---|
213 | daughterZ = Z; |
---|
214 | if (Z == 1 && A == 1) { |
---|
215 | daughterNucleus = G4Proton::Definition(); |
---|
216 | } else if (Z == 0 && A == 1) { |
---|
217 | daughterNucleus = G4Neutron::Definition(); |
---|
218 | } else { |
---|
219 | G4IonTable *theIonTable = |
---|
220 | (G4IonTable*)(G4ParticleTable::GetParticleTable()->GetIonTable()); |
---|
221 | daughterNucleus = theIonTable->GetIon(daughterZ, daughterA, theDaughterExcitation*MeV); |
---|
222 | } |
---|
223 | daughterExcitation = theDaughterExcitation; |
---|
224 | SetDaughter(index, daughterNucleus); |
---|
225 | } |
---|
226 | |
---|
227 | /////////////////////////////////////////////////////////////////////////////// |
---|
228 | // |
---|
229 | // |
---|
230 | // |
---|
231 | // |
---|
232 | G4DecayProducts *G4NuclearDecayChannel::DecayIt (G4double theParentMass) |
---|
233 | { |
---|
234 | // |
---|
235 | // |
---|
236 | // Load-up the details of the parent and daughter particles if they have not |
---|
237 | // been defined properly. |
---|
238 | // |
---|
239 | if (parent == 0) FillParent(); |
---|
240 | if (daughters == 0) FillDaughters(); |
---|
241 | // |
---|
242 | // |
---|
243 | // We want to ensure that the difference between the total |
---|
244 | // parent and daughter masses equals the energy liberated by the transition. |
---|
245 | // |
---|
246 | theParentMass = 0.0; |
---|
247 | for( G4int index=0; index < numberOfDaughters; index++) |
---|
248 | {theParentMass += daughters[index]->GetPDGMass();} |
---|
249 | theParentMass += Qtransition ; |
---|
250 | // bug fix for beta+ decay (flei 25/09/01) |
---|
251 | if (decayMode == 2) theParentMass -= 2*0.511 * MeV; |
---|
252 | // |
---|
253 | #ifdef G4VERBOSE |
---|
254 | if (GetVerboseLevel()>1) { |
---|
255 | G4cout << "G4NuclearDecayChannel::DecayIt "<< G4endl; |
---|
256 | G4cout << "the decay mass = " << theParentMass << G4endl; |
---|
257 | } |
---|
258 | #endif |
---|
259 | |
---|
260 | SetParentMass (theParentMass); |
---|
261 | |
---|
262 | // |
---|
263 | // |
---|
264 | // Define a product vector. |
---|
265 | // |
---|
266 | G4DecayProducts *products = 0; |
---|
267 | // |
---|
268 | // |
---|
269 | // Depending upon the number of daughters, select the appropriate decay |
---|
270 | // kinematics scheme. |
---|
271 | // |
---|
272 | switch (numberOfDaughters) |
---|
273 | { |
---|
274 | case 0: |
---|
275 | G4cerr << "G4NuclearDecayChannel::DecayIt "; |
---|
276 | G4cerr << " daughters not defined " <<G4endl; |
---|
277 | break; |
---|
278 | case 1: |
---|
279 | products = OneBodyDecayIt(); |
---|
280 | break; |
---|
281 | case 2: |
---|
282 | products = TwoBodyDecayIt(); |
---|
283 | break; |
---|
284 | case 3: |
---|
285 | products = BetaDecayIt(); |
---|
286 | break; |
---|
287 | default: |
---|
288 | G4cerr <<"Error in G4NuclearDecayChannel::DecayIt" <<G4endl; |
---|
289 | G4cerr <<"Number of daughters in decay = " <<numberOfDaughters <<G4endl; |
---|
290 | G4Exception(__FILE__, G4inttostring(__LINE__), FatalException, "G4NuclearDecayChannel::DecayIt"); |
---|
291 | } |
---|
292 | if (products == 0) { |
---|
293 | G4cerr << "G4NuclearDecayChannel::DecayIt "; |
---|
294 | G4cerr << *parent_name << " can not decay " << G4endl; |
---|
295 | DumpInfo(); |
---|
296 | } |
---|
297 | // |
---|
298 | // If the decay is to an excited state of the daughter nuclide, we need |
---|
299 | // to apply the photo-evaporation process. This includes the IT decay mode itself. |
---|
300 | // |
---|
301 | // needed to hold the shell idex after ICM |
---|
302 | G4int shellIndex = -1; |
---|
303 | // |
---|
304 | if (daughterExcitation > 0.0) |
---|
305 | { |
---|
306 | // |
---|
307 | // Pop the daughter nucleus off the product vector - we need to retain |
---|
308 | // the momentum of this particle. |
---|
309 | // |
---|
310 | dynamicDaughter = products->PopProducts(); |
---|
311 | G4LorentzVector daughterMomentum = dynamicDaughter->Get4Momentum(); |
---|
312 | G4ThreeVector const daughterMomentum1(static_cast<const G4LorentzVector> (daughterMomentum)); |
---|
313 | // |
---|
314 | // |
---|
315 | // Now define a G4Fragment with the correct A, Z and excitation, and declare and |
---|
316 | // initialise a G4PhotonEvaporation object. |
---|
317 | // |
---|
318 | G4Fragment nucleus(daughterA, daughterZ, daughterMomentum); |
---|
319 | G4PhotonEvaporation* deexcitation = new G4PhotonEvaporation; |
---|
320 | deexcitation->SetVerboseLevel(GetVerboseLevel()); |
---|
321 | // switch on/off internal electron conversion |
---|
322 | deexcitation->SetICM(applyICM); |
---|
323 | // set the maximum life-time for a level that will be treated. Level with life-time longer than this |
---|
324 | // will be outputed as meta-stable isotope |
---|
325 | // |
---|
326 | deexcitation->SetMaxHalfLife(halflifethreshold); |
---|
327 | // but in IT mode, we need to force the transition |
---|
328 | if (decayMode == 0) { |
---|
329 | deexcitation->RDMForced(true); |
---|
330 | } else { |
---|
331 | deexcitation->RDMForced(false); |
---|
332 | } |
---|
333 | // |
---|
334 | // Now apply the photo-evaporation |
---|
335 | // Use BreakUp() so limit to one transition at a time, if ICM is requested |
---|
336 | // this change is realted to bug#1001 (F.Lei 07/05/2010) |
---|
337 | G4FragmentVector* gammas = 0; |
---|
338 | if (applyICM) { |
---|
339 | gammas = deexcitation->BreakUp(nucleus); |
---|
340 | } else { |
---|
341 | gammas = deexcitation->BreakItUp(nucleus); |
---|
342 | } |
---|
343 | // the returned G4FragmentVector contains the residual nuclide |
---|
344 | // as its last entry. |
---|
345 | G4int nGammas=gammas->size()-1; |
---|
346 | // |
---|
347 | // Go through each gamma/e- and add it to the decay product. The angular distribution |
---|
348 | // of the gammas is isotropic, and the residual nucleus is assumed not to have suffered |
---|
349 | // any recoil as a result of this de-excitation. |
---|
350 | // |
---|
351 | for (G4int ig=0; ig<nGammas; ig++) |
---|
352 | { |
---|
353 | G4DynamicParticle *theGammaRay = new |
---|
354 | G4DynamicParticle (gammas->operator[](ig)->GetParticleDefinition(), |
---|
355 | gammas->operator[](ig)->GetMomentum()); |
---|
356 | theGammaRay -> SetProperTime(gammas->operator[](ig)->GetCreationTime()); |
---|
357 | products->PushProducts (theGammaRay); |
---|
358 | } |
---|
359 | // |
---|
360 | // now the nucleus |
---|
361 | G4double finalDaughterExcitation = gammas->operator[](nGammas)->GetExcitationEnergy(); |
---|
362 | // f.lei (03/01/03) this is needed to fix the crach in test18 |
---|
363 | if (finalDaughterExcitation <= 1.0*keV) finalDaughterExcitation = 0 ; |
---|
364 | |
---|
365 | // f.lei (07/03/05) added the delete to fix bug#711 |
---|
366 | if (dynamicDaughter) delete dynamicDaughter; |
---|
367 | |
---|
368 | G4IonTable *theIonTable = (G4IonTable*)(G4ParticleTable::GetParticleTable()->GetIonTable()); |
---|
369 | dynamicDaughter = new G4DynamicParticle |
---|
370 | (theIonTable->GetIon(daughterZ,daughterA,finalDaughterExcitation), |
---|
371 | daughterMomentum1); |
---|
372 | products->PushProducts (dynamicDaughter); |
---|
373 | |
---|
374 | // retrive the ICM shell index |
---|
375 | shellIndex = deexcitation->GetVacantShellNumber(); |
---|
376 | |
---|
377 | // |
---|
378 | // Delete/reset variables associated with the gammas. |
---|
379 | // |
---|
380 | while (!gammas->empty()) { |
---|
381 | delete *(gammas->end()-1); |
---|
382 | gammas->pop_back(); |
---|
383 | } |
---|
384 | // gammas->clearAndDestroy(); |
---|
385 | delete gammas; |
---|
386 | delete deexcitation; |
---|
387 | } |
---|
388 | // |
---|
389 | // now we have to take care of the EC product which have to go through the ARM |
---|
390 | // |
---|
391 | G4int eShell = -1; |
---|
392 | if (decayMode == 3 || decayMode == 4 || decayMode == 5) { |
---|
393 | switch (decayMode) |
---|
394 | { |
---|
395 | case KshellEC: |
---|
396 | // |
---|
397 | { |
---|
398 | eShell = 0; // --> 0 from 1 (f.lei 30/4/2008) |
---|
399 | } |
---|
400 | break; |
---|
401 | case LshellEC: |
---|
402 | // |
---|
403 | { |
---|
404 | eShell = G4int(G4UniformRand()*3)+1; |
---|
405 | } |
---|
406 | break; |
---|
407 | case MshellEC: |
---|
408 | // |
---|
409 | { |
---|
410 | // limit the shell index to 6 as specified by the ARM (F.Lei 06/05/2010) |
---|
411 | // eShell = G4int(G4UniformRand()*5)+4; |
---|
412 | eShell = G4int(G4UniformRand()*3)+4; |
---|
413 | } |
---|
414 | break; |
---|
415 | case ERROR: |
---|
416 | default: |
---|
417 | G4Exception("G4NuclearDecayChannel::DecayIt()", "601", |
---|
418 | FatalException, "Error in decay mode selection"); |
---|
419 | } |
---|
420 | } |
---|
421 | // Also have to deal with the IT case where ICM may have been applied |
---|
422 | // |
---|
423 | if (decayMode == 0) { |
---|
424 | eShell = shellIndex; |
---|
425 | } |
---|
426 | // |
---|
427 | // now apply ARM if it is requested and there is a vaccancy |
---|
428 | // |
---|
429 | if (applyARM && eShell != -1) { |
---|
430 | G4int aZ = daughterZ; |
---|
431 | if (aZ > 5 && aZ < 100) { // only applies to 5< Z <100 |
---|
432 | // Retrieve the corresponding identifier and binding energy of the selected shell |
---|
433 | const G4AtomicTransitionManager* transitionManager = G4AtomicTransitionManager::Instance(); |
---|
434 | const G4AtomicShell* shell = transitionManager->Shell(aZ, eShell); |
---|
435 | G4double bindingEnergy = shell->BindingEnergy(); |
---|
436 | G4int shellId = shell->ShellId(); |
---|
437 | |
---|
438 | G4AtomicDeexcitation* atomDeex = new G4AtomicDeexcitation(); |
---|
439 | //the default is no Auger electron generation. |
---|
440 | // Switch it on/off here! |
---|
441 | atomDeex->ActivateAugerElectronProduction(true); |
---|
442 | std::vector<G4DynamicParticle*>* armProducts = atomDeex->GenerateParticles(aZ,shellId); |
---|
443 | |
---|
444 | // pop up the daughter before insertion; |
---|
445 | // f.lei (30/04/2008) check if the total kinetic energy is less than |
---|
446 | // the shell binding energy; if true add the difference to the daughter to conserve the energy |
---|
447 | dynamicDaughter = products->PopProducts(); |
---|
448 | G4double tARMEnergy = 0.0; |
---|
449 | for (size_t i = 0; i < armProducts->size(); i++) { |
---|
450 | products->PushProducts ((*armProducts)[i]); |
---|
451 | tARMEnergy += (*armProducts)[i]->GetKineticEnergy(); |
---|
452 | } |
---|
453 | if ((bindingEnergy - tARMEnergy) > 0.1*keV){ |
---|
454 | G4double dEnergy = dynamicDaughter->GetKineticEnergy() + (bindingEnergy - tARMEnergy); |
---|
455 | dynamicDaughter->SetKineticEnergy(dEnergy); |
---|
456 | } |
---|
457 | products->PushProducts(dynamicDaughter); |
---|
458 | |
---|
459 | #ifdef G4VERBOSE |
---|
460 | if (GetVerboseLevel()>0) |
---|
461 | { |
---|
462 | G4cout <<"G4NuclearDecayChannel::Selected shell number for ARM = " <<shellId <<G4endl; |
---|
463 | G4cout <<"G4NuclearDecayChannel::ARM products = " <<armProducts->size()<<G4endl; |
---|
464 | G4cout <<" The binding energy = " << bindingEnergy << G4endl; |
---|
465 | G4cout <<" Total ARM particle kinetic energy = " << tARMEnergy << G4endl; |
---|
466 | } |
---|
467 | #endif |
---|
468 | |
---|
469 | delete armProducts; |
---|
470 | delete atomDeex; |
---|
471 | } |
---|
472 | } |
---|
473 | |
---|
474 | return products; |
---|
475 | } |
---|
476 | //////////////////////////////////////////////////////////////////////////////// |
---|
477 | // |
---|
478 | |
---|
479 | G4DecayProducts *G4NuclearDecayChannel::BetaDecayIt() |
---|
480 | |
---|
481 | { |
---|
482 | if (GetVerboseLevel()>1) G4cout << "G4Decay::BetaDecayIt()"<<G4endl; |
---|
483 | |
---|
484 | //daughters'mass |
---|
485 | G4double daughtermass[3]; |
---|
486 | G4double sumofdaughtermass = 0.0; |
---|
487 | G4double pmass = GetParentMass(); |
---|
488 | for (G4int index=0; index<3; index++) |
---|
489 | { |
---|
490 | daughtermass[index] = daughters[index]->GetPDGMass(); |
---|
491 | sumofdaughtermass += daughtermass[index]; |
---|
492 | } |
---|
493 | |
---|
494 | //create parent G4DynamicParticle at rest |
---|
495 | G4ParticleMomentum dummy; |
---|
496 | G4DynamicParticle * parentparticle = new G4DynamicParticle( parent, dummy, 0.0); |
---|
497 | |
---|
498 | //create G4Decayproducts |
---|
499 | G4DecayProducts *products = new G4DecayProducts(*parentparticle); |
---|
500 | delete parentparticle; |
---|
501 | |
---|
502 | G4double Q = pmass - sumofdaughtermass; |
---|
503 | |
---|
504 | // 09/11/2004 flei |
---|
505 | // All Beta decays are now treated with the improved 3 body decay algorithm. No more slow/fast modes |
---|
506 | /* |
---|
507 | if (BetaSimple == true) { |
---|
508 | |
---|
509 | // Use the histogramed distribution to generate the beta energy |
---|
510 | G4double daughtermomentum[2]; |
---|
511 | G4double daughterenergy[2]; |
---|
512 | daughterenergy[0] = RandomEnergy->shoot() * Q; |
---|
513 | daughtermomentum[0] = std::sqrt(daughterenergy[0]*daughterenergy[0] + |
---|
514 | 2.0*daughterenergy[0] * daughtermass[0]); |
---|
515 | // the recoil neuleus is asummed to have a maximum energy of Q/daughterA/1000. |
---|
516 | daughterenergy[1] = G4UniformRand() * Q/(1000.*daughterA); |
---|
517 | G4double recoilmomentumsquared = daughterenergy[1]*daughterenergy[1] + |
---|
518 | 2.0*daughterenergy[1] * daughtermass[1]; |
---|
519 | if (recoilmomentumsquared < 0.0) recoilmomentumsquared = 0.0; |
---|
520 | daughtermomentum[1] = std::sqrt(recoilmomentumsquared); |
---|
521 | // |
---|
522 | //create daughter G4DynamicParticle |
---|
523 | G4double costheta, sintheta, phi, sinphi, cosphi; |
---|
524 | // G4double costhetan, sinthetan, phin, sinphin, cosphin; |
---|
525 | costheta = 2.*G4UniformRand()-1.0; |
---|
526 | sintheta = std::sqrt((1.0-costheta)*(1.0+costheta)); |
---|
527 | phi = twopi*G4UniformRand()*rad; |
---|
528 | sinphi = std::sin(phi); |
---|
529 | cosphi = std::cos(phi); |
---|
530 | G4ParticleMomentum direction0(sintheta*cosphi,sintheta*sinphi,costheta); |
---|
531 | G4DynamicParticle * daughterparticle |
---|
532 | = new G4DynamicParticle( daughters[0], direction0*daughtermomentum[0]); |
---|
533 | products->PushProducts(daughterparticle); |
---|
534 | // The two products are independent in directions |
---|
535 | costheta = 2.*G4UniformRand()-1.0; |
---|
536 | sintheta = std::sqrt((1.0-costheta)*(1.0+costheta)); |
---|
537 | phi = twopi*G4UniformRand()*rad; |
---|
538 | sinphi = std::sin(phi); |
---|
539 | cosphi = std::cos(phi); |
---|
540 | G4ParticleMomentum direction1(sintheta*cosphi,sintheta*sinphi,costheta); |
---|
541 | daughterparticle |
---|
542 | = new G4DynamicParticle( daughters[1], direction1*daughtermomentum[1]); |
---|
543 | products->PushProducts(daughterparticle); |
---|
544 | |
---|
545 | // the neutrino is igored in this case |
---|
546 | |
---|
547 | } else { |
---|
548 | */ |
---|
549 | /* original slow method |
---|
550 | //calculate daughter momentum |
---|
551 | // Generate two |
---|
552 | G4double rd1, rd2; |
---|
553 | G4double daughtermomentum[3]; |
---|
554 | G4double daughterenergy[3]; |
---|
555 | G4double momentummax=0.0, momentumsum = 0.0; |
---|
556 | G4double fermif; |
---|
557 | G4BetaFermiFunction* aBetaFermiFunction; |
---|
558 | if (decayMode == 1) { |
---|
559 | // beta-decay |
---|
560 | aBetaFermiFunction = new G4BetaFermiFunction (daughterA, daughterZ); |
---|
561 | } else { |
---|
562 | // beta+decay |
---|
563 | aBetaFermiFunction = new G4BetaFermiFunction (daughterA, -daughterZ); |
---|
564 | } |
---|
565 | if (GetVerboseLevel()>1) { |
---|
566 | G4cout<< " Q = " <<Q<<G4endl; |
---|
567 | G4cout<< " daughterA = " <<daughterA<<G4endl; |
---|
568 | G4cout<< " daughterZ = " <<daughterZ<<G4endl; |
---|
569 | G4cout<< " decayMode = " <<static_cast<G4int>(decayMode) << G4endl; |
---|
570 | G4cout<< " FermiFN = " <<FermiFN<<G4endl; |
---|
571 | } |
---|
572 | do |
---|
573 | { |
---|
574 | rd1 = G4UniformRand(); |
---|
575 | rd2 = G4UniformRand(); |
---|
576 | |
---|
577 | momentummax = 0.0; |
---|
578 | momentumsum = 0.0; |
---|
579 | |
---|
580 | // daughter 0 |
---|
581 | |
---|
582 | // energy = rd2*(pmass - sumofdaughtermass); |
---|
583 | daughtermomentum[0] = std::sqrt(rd2) * std::sqrt((Q + 2.0*daughtermass[0])*Q); |
---|
584 | daughterenergy[0] = std::sqrt(daughtermomentum[0]*daughtermomentum[0] + |
---|
585 | daughtermass[0] * daughtermass[0]) - daughtermass[0]; |
---|
586 | if ( daughtermomentum[0] >momentummax )momentummax = daughtermomentum[0]; |
---|
587 | momentumsum += daughtermomentum[0]; |
---|
588 | |
---|
589 | // daughter 2 |
---|
590 | // energy = (1.-rd1)*(pmass - sumofdaughtermass); |
---|
591 | daughtermomentum[2] = std::sqrt(rd1)*std::sqrt((Q + 2.0*daughtermass[2])*Q); |
---|
592 | daughterenergy[2] = std::sqrt(daughtermomentum[2]*daughtermomentum[2] + |
---|
593 | daughtermass[2] * daughtermass[2]) - daughtermass[2]; |
---|
594 | if ( daughtermomentum[2] >momentummax )momentummax = daughtermomentum[2]; |
---|
595 | momentumsum += daughtermomentum[2]; |
---|
596 | |
---|
597 | // daughter 1 |
---|
598 | |
---|
599 | daughterenergy[1] = Q - daughterenergy[0] - daughterenergy[2]; |
---|
600 | if (daughterenergy[1] > 0.0) { |
---|
601 | daughtermomentum[1] = std::sqrt(daughterenergy[1]*daughterenergy[1] + |
---|
602 | 2.0*daughterenergy[1] * daughtermass[1]); |
---|
603 | if ( daughtermomentum[1] >momentummax ) momentummax = |
---|
604 | daughtermomentum[1]; |
---|
605 | momentumsum += daughtermomentum[1]; |
---|
606 | } else { |
---|
607 | momentummax = momentumsum = Q; |
---|
608 | } |
---|
609 | // beta particles is sampled with no coulomb effects applied above. Now |
---|
610 | // apply the Fermi function using rejection method. |
---|
611 | daughterenergy[0] = daughterenergy[0]*MeV/0.511; |
---|
612 | fermif = aBetaFermiFunction->GetFF(daughterenergy[0])/FermiFN; |
---|
613 | // fermif: normalised Fermi factor |
---|
614 | if (G4UniformRand() > fermif) momentummax = momentumsum = Q; |
---|
615 | // rejection method |
---|
616 | } while (momentummax > momentumsum - momentummax ); |
---|
617 | delete aBetaFermiFunction; |
---|
618 | |
---|
619 | // output message |
---|
620 | if (GetVerboseLevel()>1) { |
---|
621 | G4cout <<" daughter 0:" <<daughtermomentum[0]/GeV <<"[GeV/c]" <<G4endl; |
---|
622 | G4cout <<" daughter 1:" <<daughtermomentum[1]/GeV <<"[GeV/c]" <<G4endl; |
---|
623 | G4cout <<" daughter 2:" <<daughtermomentum[2]/GeV <<"[GeV/c]" <<G4endl; |
---|
624 | G4cout <<" momentum sum:" <<momentumsum/GeV <<"[GeV/c]" <<G4endl; |
---|
625 | } |
---|
626 | */ |
---|
627 | // faster method as suggested by Dirk Kruecker of FZ-Julich |
---|
628 | G4double daughtermomentum[3]; |
---|
629 | G4double daughterenergy[3]; |
---|
630 | // Use the histogramed distribution to generate the beta energy |
---|
631 | daughterenergy[0] = RandomEnergy->shoot() * Q; |
---|
632 | daughtermomentum[0] = std::sqrt(daughterenergy[0]*daughterenergy[0] + |
---|
633 | 2.0*daughterenergy[0] * daughtermass[0]); |
---|
634 | |
---|
635 | //neutrino energy distribution is flat within the kinematical limits |
---|
636 | G4double rd = 2*G4UniformRand()-1; |
---|
637 | // limits |
---|
638 | G4double Mme=pmass-daughtermass[0]; |
---|
639 | G4double K=0.5-daughtermass[1]*daughtermass[1]/(2*Mme*Mme-4*pmass*daughterenergy[0]); |
---|
640 | |
---|
641 | daughterenergy[2]=K*(Mme-daughterenergy[0]+rd*daughtermomentum[0]); |
---|
642 | daughtermomentum[2] = daughterenergy[2] ; |
---|
643 | |
---|
644 | // the recoil neuleus |
---|
645 | daughterenergy[1] = Q-daughterenergy[0]-daughterenergy[2]; |
---|
646 | G4double recoilmomentumsquared = daughterenergy[1]*daughterenergy[1] + |
---|
647 | 2.0*daughterenergy[1] * daughtermass[1]; |
---|
648 | if (recoilmomentumsquared < 0.0) recoilmomentumsquared = 0.0; |
---|
649 | daughtermomentum[1] = std::sqrt(recoilmomentumsquared); |
---|
650 | |
---|
651 | // output message |
---|
652 | if (GetVerboseLevel()>1) { |
---|
653 | G4cout <<" daughter 0:" <<daughtermomentum[0]/GeV <<"[GeV/c]" <<G4endl; |
---|
654 | G4cout <<" daughter 1:" <<daughtermomentum[1]/GeV <<"[GeV/c]" <<G4endl; |
---|
655 | G4cout <<" daughter 2:" <<daughtermomentum[2]/GeV <<"[GeV/c]" <<G4endl; |
---|
656 | } |
---|
657 | //create daughter G4DynamicParticle |
---|
658 | G4double costheta, sintheta, phi, sinphi, cosphi; |
---|
659 | G4double costhetan, sinthetan, phin, sinphin, cosphin; |
---|
660 | costheta = 2.*G4UniformRand()-1.0; |
---|
661 | sintheta = std::sqrt((1.0-costheta)*(1.0+costheta)); |
---|
662 | phi = twopi*G4UniformRand()*rad; |
---|
663 | sinphi = std::sin(phi); |
---|
664 | cosphi = std::cos(phi); |
---|
665 | G4ParticleMomentum direction0(sintheta*cosphi,sintheta*sinphi,costheta); |
---|
666 | G4DynamicParticle * daughterparticle |
---|
667 | = new G4DynamicParticle( daughters[0], direction0*daughtermomentum[0]); |
---|
668 | products->PushProducts(daughterparticle); |
---|
669 | |
---|
670 | costhetan = (daughtermomentum[1]*daughtermomentum[1]- |
---|
671 | daughtermomentum[2]*daughtermomentum[2]- |
---|
672 | daughtermomentum[0]*daughtermomentum[0])/ |
---|
673 | (2.0*daughtermomentum[2]*daughtermomentum[0]); |
---|
674 | // added the following test to avoid rounding erros. A problem |
---|
675 | // reported bye Ben Morgan of Uni.Warwick |
---|
676 | if (costhetan > 1.) costhetan = 1.; |
---|
677 | if (costhetan < -1.) costhetan = -1.; |
---|
678 | sinthetan = std::sqrt((1.0-costhetan)*(1.0+costhetan)); |
---|
679 | phin = twopi*G4UniformRand()*rad; |
---|
680 | sinphin = std::sin(phin); |
---|
681 | cosphin = std::cos(phin); |
---|
682 | G4ParticleMomentum direction2; |
---|
683 | direction2.setX( sinthetan*cosphin*costheta*cosphi - |
---|
684 | sinthetan*sinphin*sinphi + costhetan*sintheta*cosphi); |
---|
685 | direction2.setY( sinthetan*cosphin*costheta*sinphi + |
---|
686 | sinthetan*sinphin*cosphi + costhetan*sintheta*sinphi); |
---|
687 | direction2.setZ( -sinthetan*cosphin*sintheta + |
---|
688 | costhetan*costheta); |
---|
689 | daughterparticle = new G4DynamicParticle |
---|
690 | ( daughters[2], direction2*(daughtermomentum[2]/direction2.mag())); |
---|
691 | products->PushProducts(daughterparticle); |
---|
692 | |
---|
693 | daughterparticle = |
---|
694 | new G4DynamicParticle (daughters[1], |
---|
695 | (direction0*daughtermomentum[0] + |
---|
696 | direction2*(daughtermomentum[2]/direction2.mag()))*(-1.0)); |
---|
697 | products->PushProducts(daughterparticle); |
---|
698 | // } |
---|
699 | // delete daughterparticle; |
---|
700 | |
---|
701 | if (GetVerboseLevel()>1) { |
---|
702 | G4cout << "G4NuclearDecayChannel::BetaDecayIt "; |
---|
703 | G4cout << " create decay products in rest frame " <<G4endl; |
---|
704 | products->DumpInfo(); |
---|
705 | } |
---|
706 | return products; |
---|
707 | } |
---|