ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4BetaPlusDecay.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4BetaPlusDecay.cc
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 //
27 // //
28 // File: G4BetaPlusDecay.cc //
29 // Author: D.H. Wright (SLAC) //
30 // Date: 14 November 2014 //
31 // //
33 
34 #include "G4BetaPlusDecay.hh"
36 #include "G4IonTable.hh"
37 #include "G4ThreeVector.hh"
38 #include "G4DynamicParticle.hh"
39 #include "G4DecayProducts.hh"
40 #include "G4PhysicalConstants.hh"
41 #include "G4SystemOfUnits.hh"
42 #include <iostream>
43 #include <iomanip>
44 
46  const G4double& branch, const G4double& e0,
47  const G4double& excitationE,
48  const G4Ions::G4FloatLevelBase& flb,
49  const G4BetaDecayType& betaType)
50  : G4NuclearDecay("beta+ decay", BetaPlus, excitationE, flb),
51  endpointEnergy(e0 - 2.*CLHEP::electron_mass_c2)
52 {
53  SetParent(theParentNucleus); // Store name of parent nucleus, delete G4MT_parent
54  SetBR(branch);
55 
57  G4IonTable* theIonTable =
59  G4int daughterZ = theParentNucleus->GetAtomicNumber() - 1;
60  G4int daughterA = theParentNucleus->GetAtomicMass();
61  SetDaughter(0, theIonTable->GetIon(daughterZ, daughterA, excitationE, flb) );
62  SetUpBetaSpectrumSampler(daughterZ, daughterA, betaType);
63  SetDaughter(1, "e+");
64  SetDaughter(2, "nu_e");
65 }
66 
67 
69 {
70  delete spectrumSampler;
71 }
72 
73 
75 {
76  // Fill G4MT_parent with theParentNucleus (stored by SetParent in ctor)
78 
79  // Fill G4MT_daughters with e-, nu and residual nucleus (stored by SetDaughter)
81 
82  G4double parentMass = G4MT_parent->GetPDGMass();
83  G4double eMass = G4MT_daughters[1]->GetPDGMass();
84  G4double nucleusMass = G4MT_daughters[0]->GetPDGMass();
85  // Set up final state
86  // parentParticle is set at rest here because boost with correct momentum
87  // is done later
88  G4DynamicParticle parentParticle(G4MT_parent, G4ThreeVector(0,0,0), 0.0);
89  G4DecayProducts* products = new G4DecayProducts(parentParticle);
90 
91  if (spectrumSampler) {
92  // Generate positron isotropic in angle, with energy from stored spectrum
93  G4double eKE = endpointEnergy*spectrumSampler->shoot(G4Random::getTheEngine() );
94  G4double eMomentum = std::sqrt(eKE*(eKE + 2.*eMass) );
95 
96  G4double cosTheta = 2.*G4UniformRand() - 1.0;
97  G4double sinTheta = std::sqrt(1.0 - cosTheta*cosTheta);
99  G4double sinPhi = std::sin(phi);
100  G4double cosPhi = std::cos(phi);
101 
102  G4ParticleMomentum eDirection(sinTheta*cosPhi, sinTheta*sinPhi, cosTheta);
103  G4DynamicParticle* dynamicPositron
104  = new G4DynamicParticle(G4MT_daughters[1], eDirection*eMomentum);
105  products->PushProducts(dynamicPositron);
106 
107  // Generate neutrino with angle relative to positron, and energy from
108  // energy-momentum conservation using endpoint energy of reaction
109  G4double cosThetaENu = 2.*G4UniformRand() - 1.;
110  G4double eTE = eMass + eKE;
111  G4double nuEnergy = ((endpointEnergy - eKE)*(parentMass + nucleusMass - eTE)
112  - eMomentum*eMomentum)/(parentMass - eTE + eMomentum*cosThetaENu)/2.;
113 
114  G4double sinThetaENu = std::sqrt(1.0 - cosThetaENu*cosThetaENu);
115  phi = twopi*G4UniformRand()*rad;
116  G4double sinPhiNu = std::sin(phi);
117  G4double cosPhiNu = std::cos(phi);
118 
119  G4ParticleMomentum nuDirection;
120  nuDirection.setX(sinThetaENu*cosPhiNu*cosTheta*cosPhi -
121  sinThetaENu*sinPhiNu*sinPhi + cosThetaENu*sinTheta*cosPhi);
122  nuDirection.setY(sinThetaENu*cosPhiNu*cosTheta*sinPhi +
123  sinThetaENu*sinPhiNu*cosPhi + cosThetaENu*sinTheta*sinPhi);
124  nuDirection.setZ(-sinThetaENu*cosPhiNu*sinTheta + cosThetaENu*cosTheta);
125 
126  G4DynamicParticle* dynamicNeutrino
127  = new G4DynamicParticle(G4MT_daughters[2], nuDirection*nuEnergy);
128  products->PushProducts(dynamicNeutrino);
129 
130  // Generate daughter nucleus from sum of positron and neutrino 4-vectors:
131  // p_D = - p_e - p_nu
132  G4DynamicParticle* dynamicDaughter =
134  -eDirection*eMomentum - nuDirection*nuEnergy);
135  products->PushProducts(dynamicDaughter);
136 
137  } else {
138  // positron energy below threshold -> no decay
139  G4DynamicParticle* noDecay =
140  new G4DynamicParticle(G4MT_parent, G4ThreeVector(0,0,0), 0.0);
141  products->PushProducts(noDecay);
142  }
143 
144  // Check energy conservation against endpoint value, not nuclear masses
145  /*
146  G4int nProd = products->entries();
147  G4DynamicParticle* temp = 0;
148  G4double Esum = 0.0;
149  for (G4int i = 0; i < nProd; i++) {
150  temp = products->operator[](i);
151  Esum += temp->GetKineticEnergy();
152  }
153  G4double eCons = (endpointEnergy - Esum)/keV;
154  if (eCons > 0.001) G4cout << " Beta+ check: eCons (keV) = " << eCons << G4endl;
155  */
156  return products;
157 }
158 
159 
160 void
162  const G4int& daughterA,
163  const G4BetaDecayType& betaType)
164 {
166  G4BetaDecayCorrections corrections(-daughterZ, daughterA);
167  spectrumSampler = 0;
168 
169  // Check for cases in which Q < 2Me (e.g. z67.a162)
170  if (e0 > 0.) {
171  // Array to store spectrum pdf
172  G4int npti = 100;
173  G4double* pdf = new G4double[npti];
174 
175  G4double e; // Total positron energy in units of electron mass
176  G4double p; // Positron momentum in units of electron mass
177  G4double f; // Spectral shap function
178  for (G4int ptn = 0; ptn < npti; ptn++) {
179  // Calculate simple phase space
180  e = 1. + e0*(ptn + 0.5)/G4double(npti);
181  p = std::sqrt(e*e - 1.);
182  f = p*e*(e0 - e + 1.)*(e0 - e + 1.);
183 
184  // Apply Fermi factor to get allowed shape
185  f *= corrections.FermiFunction(e);
186 
187  // Apply shape factor for forbidden transitions
188  f *= corrections.ShapeFactor(betaType, p, e0-e+1.);
189  pdf[ptn] = f;
190  }
191  spectrumSampler = new G4RandGeneral(pdf, npti);
192  delete[] pdf;
193  }
194 }
195 
196 
198 {
199  G4cout << " G4BetaPlusDecay for parent nucleus " << GetParentName() << G4endl;
200  G4cout << " decays to " << GetDaughterName(0) << " , " << GetDaughterName(1)
201  << " and " << GetDaughterName(2) << " with branching ratio " << GetBR()
202  << "% and endpoint energy " << endpointEnergy/keV << " keV " << G4endl;
203 }
204