source: trunk/source/processes/hadronic/models/radioactive_decay/src/G4NuclearDecayChannel.cc@ 1036

Last change on this file since 1036 was 962, checked in by garnier, 17 years ago

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