ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4BetaMinusDecay.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4BetaMinusDecay.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: G4BetaMinusDecay.cc //
29 // Author: D.H. Wright (SLAC) //
30 // Date: 25 October 2014 //
31 // //
33 
34 #include "G4BetaMinusDecay.hh"
36 #include "G4ThreeVector.hh"
37 #include "G4DynamicParticle.hh"
38 #include "G4DecayProducts.hh"
39 #include "G4PhysicalConstants.hh"
40 #include "G4SystemOfUnits.hh"
41 #include <iostream>
42 #include <iomanip>
43 
45  const G4double& branch, const G4double& e0,
46  const G4double& excitationE,
47  const G4Ions::G4FloatLevelBase& flb,
48  const G4BetaDecayType& betaType)
49  : G4NuclearDecay("beta- decay", BetaMinus, excitationE, flb), endpointEnergy(e0)
50 {
51  SetParent(theParentNucleus); // Store name of parent nucleus, delete G4MT_parent
52  SetBR(branch);
53 
55  G4IonTable* theIonTable =
57  G4int daughterZ = theParentNucleus->GetAtomicNumber() + 1;
58  G4int daughterA = theParentNucleus->GetAtomicMass();
59  SetDaughter(0, theIonTable->GetIon(daughterZ, daughterA, excitationE, flb) );
60  SetDaughter(1, "e-");
61  SetDaughter(2, "anti_nu_e");
62 
63  SetUpBetaSpectrumSampler(daughterZ, daughterA, betaType);
64 }
65 
66 
68 {
69  delete spectrumSampler;
70 }
71 
72 
74 {
75  // Fill G4MT_parent with theParentNucleus (stored by SetParent in ctor)
77 
78  // Fill G4MT_daughters with e-, nu and residual nucleus (stored by SetDaughter)
80 
81  G4double parentMass = G4MT_parent->GetPDGMass();
82  G4double eMass = G4MT_daughters[1]->GetPDGMass();
83  G4double nucleusMass = G4MT_daughters[0]->GetPDGMass();
84  // Set up final state
85  // parentParticle is set at rest here because boost with correct momentum
86  // is done later
87  G4DynamicParticle parentParticle(G4MT_parent, G4ThreeVector(0,0,0), 0.0);
88  G4DecayProducts* products = new G4DecayProducts(parentParticle);
89 
90  if (spectrumSampler) {
91  // Electron, neutrino and daughter nucleus energies
92  G4double eKE = endpointEnergy*spectrumSampler->shoot(G4Random::getTheEngine() );
93  G4double eMomentum = std::sqrt(eKE*(eKE + 2.*eMass) );
94 
95  G4double cosThetaENu = 2.*G4UniformRand() - 1.;
96  G4double eTE = eMass + eKE;
97  G4double nuEnergy = ((endpointEnergy - eKE)*(parentMass + nucleusMass - eTE)
98  - eMomentum*eMomentum)/(parentMass - eTE + eMomentum*cosThetaENu)/2.;
99 
100  // Electron 4-vector, isotropic angular distribution
101  G4double cosTheta = 2.*G4UniformRand() - 1.0;
102  G4double sinTheta = std::sqrt(1.0 - cosTheta*cosTheta);
103 
105  G4double sinPhi = std::sin(phi);
106  G4double cosPhi = std::cos(phi);
107 
108  G4ParticleMomentum eDirection(sinTheta*cosPhi, sinTheta*sinPhi, cosTheta);
109  G4DynamicParticle* dynamicElectron
110  = new G4DynamicParticle(G4MT_daughters[1], eDirection*eMomentum);
111  products->PushProducts(dynamicElectron);
112 
113  // Neutrino 4-vector
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  // Daughter nucleus 4-vector
131  // p_D = - p_e - p_nu
132  G4DynamicParticle* dynamicDaughter =
134  -eDirection*eMomentum - nuDirection*nuEnergy);
135  products->PushProducts(dynamicDaughter);
136 
137  } else {
138  // electron 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 Q 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  // G4cout << temp->GetParticleDefinition()->GetParticleName() << " has "
152  // << temp->GetTotalEnergy()/keV << " keV " << G4endl;
153  Esum += temp->GetKineticEnergy();
154  }
155  G4double eCons = (endpointEnergy - Esum)/keV;
156  if (std::abs(eCons) > 0.001) G4cout << " Beta- check: eCons = " << eCons << G4endl;
157  */
158  return products;
159 }
160 
161 
162 void
164  const G4int& daughterA,
165  const G4BetaDecayType& betaType)
166 {
168  G4BetaDecayCorrections corrections(daughterZ, daughterA);
169  spectrumSampler = 0;
170 
171  if (e0 > 0) {
172  // Array to store spectrum pdf
173  G4int npti = 100;
174  G4double* pdf = new G4double[npti];
175 
176  G4double e; // Total electron energy in units of electron mass
177  G4double p; // Electron momentum in units of electron mass
178  G4double f; // Spectral shape function
179  for (G4int ptn = 0; ptn < npti; ptn++) {
180  // Calculate simple phase space
181  e = 1. + e0*(G4double(ptn) + 0.5)/G4double(npti);
182  p = std::sqrt(e*e - 1.);
183  f = p*e*(e0 - e + 1.)*(e0 - e + 1.);
184 
185  // Apply Fermi factor to get allowed shape
186  f *= corrections.FermiFunction(e);
187 
188  // Apply shape factor for forbidden transitions
189  f *= corrections.ShapeFactor(betaType, p, e0-e+1.);
190  pdf[ptn] = f;
191  }
192  spectrumSampler = new G4RandGeneral(pdf, npti);
193  delete[] pdf;
194  }
195 }
196 
197 
199 {
200  G4cout << " G4BetaMinusDecay for parent nucleus " << GetParentName() << G4endl;
201  G4cout << " decays to " << GetDaughterName(0) << " , " << GetDaughterName(1)
202  << " and " << GetDaughterName(2) << " with branching ratio " << GetBR()
203  << "% and endpoint energy " << endpointEnergy/keV << " keV " << G4endl;
204 }
205