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

Last change on this file since 1337 was 1315, checked in by garnier, 15 years ago

update geant4-09-04-beta-cand-01 interfaces-V09-03-09 vis-V09-03-08
File size: 25.8 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 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//
112G4NuclearDecayChannel::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//
145G4NuclearDecayChannel::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
191void 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//
232G4DecayProducts *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
479G4DecayProducts *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}
Note: See TracBrowser for help on using the repository browser.