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

Last change on this file was 1315, checked in by garnier, 14 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.