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

Last change on this file since 843 was 819, checked in by garnier, 16 years ago

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