Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/processes/hadronic/models/radioactive_decay/src/G4BetaMinusDecay.cc

Version: [ ReleaseNotes ] [ 1.0 ] [ 1.1 ] [ 2.0 ] [ 3.0 ] [ 3.1 ] [ 3.2 ] [ 4.0 ] [ 4.0.p1 ] [ 4.0.p2 ] [ 4.1 ] [ 4.1.p1 ] [ 5.0 ] [ 5.0.p1 ] [ 5.1 ] [ 5.1.p1 ] [ 5.2 ] [ 5.2.p1 ] [ 5.2.p2 ] [ 6.0 ] [ 6.0.p1 ] [ 6.1 ] [ 6.2 ] [ 6.2.p1 ] [ 6.2.p2 ] [ 7.0 ] [ 7.0.p1 ] [ 7.1 ] [ 7.1.p1 ] [ 8.0 ] [ 8.0.p1 ] [ 8.1 ] [ 8.1.p1 ] [ 8.1.p2 ] [ 8.2 ] [ 8.2.p1 ] [ 8.3 ] [ 8.3.p1 ] [ 8.3.p2 ] [ 9.0 ] [ 9.0.p1 ] [ 9.0.p2 ] [ 9.1 ] [ 9.1.p1 ] [ 9.1.p2 ] [ 9.1.p3 ] [ 9.2 ] [ 9.2.p1 ] [ 9.2.p2 ] [ 9.2.p3 ] [ 9.2.p4 ] [ 9.3 ] [ 9.3.p1 ] [ 9.3.p2 ] [ 9.4 ] [ 9.4.p1 ] [ 9.4.p2 ] [ 9.4.p3 ] [ 9.4.p4 ] [ 9.5 ] [ 9.5.p1 ] [ 9.5.p2 ] [ 9.6 ] [ 9.6.p1 ] [ 9.6.p2 ] [ 9.6.p3 ] [ 9.6.p4 ] [ 10.0 ] [ 10.0.p1 ] [ 10.0.p2 ] [ 10.0.p3 ] [ 10.0.p4 ] [ 10.1 ] [ 10.1.p1 ] [ 10.1.p2 ] [ 10.1.p3 ] [ 10.2 ] [ 10.2.p1 ] [ 10.2.p2 ] [ 10.2.p3 ] [ 10.3 ] [ 10.3.p1 ] [ 10.3.p2 ] [ 10.3.p3 ] [ 10.4 ] [ 10.4.p1 ] [ 10.4.p2 ] [ 10.4.p3 ] [ 10.5 ] [ 10.5.p1 ] [ 10.6 ] [ 10.6.p1 ] [ 10.6.p2 ] [ 10.6.p3 ] [ 10.7 ] [ 10.7.p1 ] [ 10.7.p2 ] [ 10.7.p3 ] [ 10.7.p4 ] [ 11.0 ] [ 11.0.p1 ] [ 11.0.p2 ] [ 11.0.p3, ] [ 11.0.p4 ] [ 11.1 ] [ 11.1.1 ] [ 11.1.2 ] [ 11.1.3 ] [ 11.2 ] [ 11.2.1 ] [ 11.2.2 ] [ 11.3.0 ]

Diff markup

Differences between /processes/hadronic/models/radioactive_decay/src/G4BetaMinusDecay.cc (Version 11.3.0) and /processes/hadronic/models/radioactive_decay/src/G4BetaMinusDecay.cc (Version 10.7)


  1 //                                                  1 //
  2 // *******************************************      2 // ********************************************************************
  3 // * License and Disclaimer                         3 // * License and Disclaimer                                           *
  4 // *                                                4 // *                                                                  *
  5 // * The  Geant4 software  is  copyright of th      5 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
  6 // * the Geant4 Collaboration.  It is provided      6 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
  7 // * conditions of the Geant4 Software License      7 // * conditions of the Geant4 Software License,  included in the file *
  8 // * LICENSE and available at  http://cern.ch/      8 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
  9 // * include a list of copyright holders.           9 // * include a list of copyright holders.                             *
 10 // *                                               10 // *                                                                  *
 11 // * Neither the authors of this software syst     11 // * Neither the authors of this software system, nor their employing *
 12 // * institutes,nor the agencies providing fin     12 // * institutes,nor the agencies providing financial support for this *
 13 // * work  make  any representation or  warran     13 // * work  make  any representation or  warranty, express or implied, *
 14 // * regarding  this  software system or assum     14 // * regarding  this  software system or assume any liability for its *
 15 // * use.  Please see the license in the file      15 // * use.  Please see the license in the file  LICENSE  and URL above *
 16 // * for the full disclaimer and the limitatio     16 // * for the full disclaimer and the limitation of liability.         *
 17 // *                                               17 // *                                                                  *
 18 // * This  code  implementation is the result      18 // * This  code  implementation is the result of  the  scientific and *
 19 // * technical work of the GEANT4 collaboratio     19 // * technical work of the GEANT4 collaboration.                      *
 20 // * By using,  copying,  modifying or  distri     20 // * By using,  copying,  modifying or  distributing the software (or *
 21 // * any work based  on the software)  you  ag     21 // * any work based  on the software)  you  agree  to acknowledge its *
 22 // * use  in  resulting  scientific  publicati     22 // * use  in  resulting  scientific  publications,  and indicate your *
 23 // * acceptance of all terms of the Geant4 Sof     23 // * acceptance of all terms of the Geant4 Software license.          *
 24 // *******************************************     24 // ********************************************************************
 25 //                                                 25 //
 26 //////////////////////////////////////////////     26 ////////////////////////////////////////////////////////////////////////////////
 27 //                                                 27 //                                                                            //
 28 //  File:   G4BetaMinusDecay.cc                    28 //  File:   G4BetaMinusDecay.cc                                               //
 29 //  Author: D.H. Wright (SLAC)                     29 //  Author: D.H. Wright (SLAC)                                                //
 30 //  Date:   25 October 2014                        30 //  Date:   25 October 2014                                                   //
 31 //  Modifications:                             << 
 32 //    23.08.2023 V.Ivanchenko make it thread s << 
 33 //                                                 31 //                                                                            //
 34 //////////////////////////////////////////////     32 ////////////////////////////////////////////////////////////////////////////////
 35                                                    33 
 36 #include "G4BetaMinusDecay.hh"                     34 #include "G4BetaMinusDecay.hh"
 37 #include "G4BetaDecayCorrections.hh"               35 #include "G4BetaDecayCorrections.hh"
 38 #include "G4ThreeVector.hh"                        36 #include "G4ThreeVector.hh"
 39 #include "G4LorentzVector.hh"                  << 
 40 #include "G4DynamicParticle.hh"                    37 #include "G4DynamicParticle.hh"
 41 #include "G4DecayProducts.hh"                      38 #include "G4DecayProducts.hh"
 42 #include "G4PhysicalConstants.hh"                  39 #include "G4PhysicalConstants.hh"
 43 #include "G4SystemOfUnits.hh"                      40 #include "G4SystemOfUnits.hh"
 44 #include "G4Electron.hh"                       << 
 45 #include "G4AntiNeutrinoE.hh"                  << 
 46 #include "G4RandomDirection.hh"                << 
 47 #include "G4BetaSpectrumSampler.hh"            << 
 48 #include <iostream>                                41 #include <iostream>
 49 #include <iomanip>                                 42 #include <iomanip>
 50                                                    43 
 51 namespace {                                    << 
 52   const G4double eMass = CLHEP::electron_mass_ << 
 53 }                                              << 
 54                                                << 
 55 G4BetaMinusDecay::G4BetaMinusDecay(const G4Par     44 G4BetaMinusDecay::G4BetaMinusDecay(const G4ParticleDefinition* theParentNucleus,
 56                                    const G4dou     45                                    const G4double& branch, const G4double& e0,
 57                                    const G4dou     46                                    const G4double& excitationE,
 58                                    const G4Ion     47                                    const G4Ions::G4FloatLevelBase& flb,
 59                                    const G4Bet     48                                    const G4BetaDecayType& betaType)
 60  : G4NuclearDecay("beta- decay", BetaMinus, ex <<  49  : G4NuclearDecay("beta- decay", BetaMinus, excitationE, flb), endpointEnergy(e0)
 61    maxEnergy(e0/eMass),                        << 
 62    estep(maxEnergy/(G4double)(npti - 1))       << 
 63 {                                                  50 {
 64   SetParent(theParentNucleus);  // Store name      51   SetParent(theParentNucleus);  // Store name of parent nucleus, delete G4MT_parent 
 65   SetBR(branch);                                   52   SetBR(branch);
 66   SetNumberOfDaughters(3);                     << 
 67                                                    53 
 68   fPrimaryIon = theParentNucleus;              <<  54   SetNumberOfDaughters(3);
 69   fLepton = G4Electron::Electron();            <<  55   G4IonTable* theIonTable =
 70   fNeutrino = G4AntiNeutrinoE::AntiNeutrinoE() <<  56     (G4IonTable*)(G4ParticleTable::GetParticleTable()->GetIonTable());
 71                                                << 
 72   G4IonTable* theIonTable = G4ParticleTable::G << 
 73   G4int daughterZ = theParentNucleus->GetAtomi     57   G4int daughterZ = theParentNucleus->GetAtomicNumber() + 1;
 74   G4int daughterA = theParentNucleus->GetAtomi     58   G4int daughterA = theParentNucleus->GetAtomicMass();
 75   fResIon = const_cast<const G4ParticleDefinit <<  59   SetDaughter(0, theIonTable->GetIon(daughterZ, daughterA, excitationE, flb) );
 76                                                <<  60   SetDaughter(1, "e-");
 77   parentMass = theParentNucleus->GetPDGMass(); <<  61   SetDaughter(2, "anti_nu_e");
 78   resMass = fResIon->GetPDGMass();             << 
 79                                                    62 
 80   SetUpBetaSpectrumSampler(daughterZ, daughter     63   SetUpBetaSpectrumSampler(daughterZ, daughterA, betaType);
                                                   >>  64 }
 81                                                    65 
 82   SetDaughter(0, fResIon);                     << 
 83   SetDaughter(1, fLepton);                     << 
 84   SetDaughter(2, fNeutrino);                   << 
 85                                                    66 
                                                   >>  67 G4BetaMinusDecay::~G4BetaMinusDecay()
                                                   >>  68 {
                                                   >>  69   delete spectrumSampler;
                                                   >>  70 }
                                                   >>  71 
                                                   >>  72 
                                                   >>  73 G4DecayProducts* G4BetaMinusDecay::DecayIt(G4double)
                                                   >>  74 {
 86   // Fill G4MT_parent with theParentNucleus (s     75   // Fill G4MT_parent with theParentNucleus (stored by SetParent in ctor)  
 87   CheckAndFillParent();                            76   CheckAndFillParent();
 88                                                    77 
 89   // Fill G4MT_daughters with e-, nu and resid     78   // Fill G4MT_daughters with e-, nu and residual nucleus (stored by SetDaughter)  
 90   CheckAndFillDaughters();                         79   CheckAndFillDaughters();
 91 }                                              << 
 92                                                    80 
 93 G4DecayProducts* G4BetaMinusDecay::DecayIt(G4d <<  81   G4double parentMass = G4MT_parent->GetPDGMass();
 94 {                                              <<  82   G4double eMass = G4MT_daughters[1]->GetPDGMass();
                                                   >>  83   G4double nucleusMass = G4MT_daughters[0]->GetPDGMass();
 95   // Set up final state                            84   // Set up final state
 96   // parentParticle is set at rest here becaus     85   // parentParticle is set at rest here because boost with correct momentum 
 97   // is done later                                 86   // is done later
 98   G4DynamicParticle prim(fPrimaryIon, G4ThreeV <<  87   G4DynamicParticle parentParticle(G4MT_parent, G4ThreeVector(0,0,0), 0.0);
 99   G4DecayProducts* products = new G4DecayProdu <<  88   G4DecayProducts* products = new G4DecayProducts(parentParticle);
100                                                << 
101   // Generate electron isotropic in angle, wit << 
102   const G4double eKE = eMass*G4BetaSpectrumSam << 
103                                                << 
104   G4double eMomentum = std::sqrt(eKE*(eKE + 2. << 
105   G4ThreeVector dir = G4RandomDirection();     << 
106   G4DynamicParticle* dp = new G4DynamicParticl << 
107   products->PushProducts(dp);                  << 
108   /*                                           << 
109   G4cout << "G4BetaPlusDecay::DecayIt: " << fP << 
110    << " -> " << fResIon->GetParticleName() <<  << 
111    << " + " << fNeutrino->GetParticleName() << << 
112          << G4endl;                            << 
113   */                                           << 
114                                                    89 
115   // 4-momentum of residual ion and neutrino   <<  90   if (spectrumSampler) {
116   G4LorentzVector lv(-eMomentum*dir.x(), -eMom <<  91     // Electron, neutrino and daughter nucleus energies
117                      parentMass - eKE - eMass) <<  92     G4double eKE = endpointEnergy*spectrumSampler->shoot(G4Random::getTheEngine() );
118                                                <<  93     G4double eMomentum = std::sqrt(eKE*(eKE + 2.*eMass) );
119   // centrum of mass system                    <<  94 
120   G4double M = lv.mag();                       <<  95     G4double cosThetaENu = 2.*G4UniformRand() - 1.;
121   const G4double elim = CLHEP::eV;             <<  96     G4double eTE = eMass + eKE;
122   G4double edel = M - resMass;                 <<  97     G4double nuEnergy = ((endpointEnergy - eKE)*(parentMass + nucleusMass - eTE)
123   // Free energy should be above limit         <<  98             - eMomentum*eMomentum)/(parentMass - eTE + eMomentum*cosThetaENu)/2.;
124   if (edel >= elim) {                          <<  99 
125     // neutrino                                << 100     // Electron 4-vector, isotropic angular distribution
126     G4double eNu = 0.5*(M - resMass*resMass/M) << 101     G4double cosTheta = 2.*G4UniformRand() - 1.0;
127     G4LorentzVector lvnu(eNu*G4RandomDirection << 102     G4double sinTheta = std::sqrt(1.0 - cosTheta*cosTheta);
128     lvnu.boost(lv.boostVector());              << 103 
129     dir = lvnu.vect().unit();                  << 104     G4double phi = twopi*G4UniformRand()*rad;
130     dp = new G4DynamicParticle(fNeutrino, dir, << 105     G4double sinPhi = std::sin(phi);
131     products->PushProducts(dp);                << 106     G4double cosPhi = std::cos(phi);
132                                                << 107 
133     // residual                                << 108     G4ParticleMomentum eDirection(sinTheta*cosPhi, sinTheta*sinPhi, cosTheta);
134     lv -= lvnu;                                << 109     G4DynamicParticle* dynamicElectron
135     dir = lv.vect().unit();                    << 110       = new G4DynamicParticle(G4MT_daughters[1], eDirection*eMomentum);
136     G4double ekin = std::max(lv.e() - resMass, << 111     products->PushProducts(dynamicElectron);
137     dp = new G4DynamicParticle(fResIon, dir, e << 112 
138     products->PushProducts(dp);                << 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 =
                                                   >> 133       new G4DynamicParticle(G4MT_daughters[0],
                                                   >> 134                             -eDirection*eMomentum - nuDirection*nuEnergy);
                                                   >> 135     products->PushProducts(dynamicDaughter);
139                                                   136 
140   } else {                                        137   } else {
141     // neglecting relativistic kinematic and g << 138     // electron energy below threshold -> no decay
142     dp = new G4DynamicParticle(fNeutrino, G4Ra << 139     G4DynamicParticle* noDecay =
143     products->PushProducts(dp);                << 140       new G4DynamicParticle(G4MT_parent, G4ThreeVector(0,0,0), 0.0);
144     dp = new G4DynamicParticle(fResIon, G4Thre << 141     products->PushProducts(noDecay);
145     products->PushProducts(dp);                << 
146   }                                               142   }
147                                                   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   */
148   return products;                                158   return products;
149 }                                                 159 }
150                                                   160 
151                                                   161 
152 void                                              162 void
153 G4BetaMinusDecay::SetUpBetaSpectrumSampler(con    163 G4BetaMinusDecay::SetUpBetaSpectrumSampler(const G4int& daughterZ,
154                                            con    164                                            const G4int& daughterA,
155                                            con    165                                            const G4BetaDecayType& betaType)
156 {                                                 166 {
157   cdf[0] = 0.0;                                << 167   G4double e0 = endpointEnergy/CLHEP::electron_mass_c2;
158                                                << 168   G4BetaDecayCorrections corrections(daughterZ, daughterA);
159   // Check for cases in which Q < 0            << 169   spectrumSampler = 0;
160   if (maxEnergy > 0.) {                        << 170 
161     G4BetaDecayCorrections corrections(daughte << 171   if (e0 > 0) {
162                                                << 172     // Array to store spectrum pdf
163     // Fill array to store cumulative spectrum << 173     G4int npti = 100;
164     G4double ex;  // Kinetic energy normalized << 174     G4double* pdf = new G4double[npti];
165     G4double p;   // Momentum in units of elec << 175 
166     G4double f;   // Spectral shape function   << 176     G4double e;  // Total electron energy in units of electron mass
167     G4double f0 = 0.0;                         << 177     G4double p;  // Electron momentum in units of electron mass
168     G4double sum = 0.0;                        << 178     G4double f;  // Spectral shape function
169     for (G4int i = 1; i < npti-1; ++i) {       << 179     for (G4int ptn = 0; ptn < npti; ptn++) {
170       ex = estep*i;                            << 180       // Calculate simple phase space
171       p = std::sqrt(ex*(ex + 2.));             << 181       e = 1. + e0*(G4double(ptn) + 0.5)/G4double(npti);
172       f = p*(1. + ex)*(maxEnergy - ex)*(maxEne << 182       p = std::sqrt(e*e - 1.);
                                                   >> 183       f = p*e*(e0 - e + 1.)*(e0 - e + 1.);
173                                                   184 
174       // Apply Fermi factor to get allowed sha    185       // Apply Fermi factor to get allowed shape
175       f *= corrections.FermiFunction(1. + ex); << 186       f *= corrections.FermiFunction(e);
176                                                   187 
177       // Apply shape factor for forbidden tran    188       // Apply shape factor for forbidden transitions
178       f *= corrections.ShapeFactor(betaType, p << 189       f *= corrections.ShapeFactor(betaType, p, e0-e+1.);
179       sum += f + f0;                           << 190       pdf[ptn] = f;
180       cdf[i] = sum;                            << 
181       f0 = f;                                  << 
182     }                                             191     }
183     cdf[npti-1] = sum + f0;                    << 192     spectrumSampler = new G4RandGeneral(pdf, npti);
184   } else {                                     << 193     delete[] pdf;
185     for (G4int i = 1; i < npti; ++i) { cdf[i]  << 
186   }                                               194   }
187 }                                                 195 }
188                                                   196 
189                                                   197 
190 void G4BetaMinusDecay::DumpNuclearInfo()          198 void G4BetaMinusDecay::DumpNuclearInfo()
191 {                                                 199 {
192   G4cout << " G4BetaMinusDecay  " << fPrimaryI << 200   G4cout << " G4BetaMinusDecay for parent nucleus " << GetParentName() << G4endl;
193    << " -> " << fResIon->GetParticleName() <<  << 201   G4cout << " decays to " << GetDaughterName(0) << " , " << GetDaughterName(1) 
194    << " + " << fNeutrino->GetParticleName() << << 202          << " and " << GetDaughterName(2) << " with branching ratio " << GetBR()
195    << maxEnergy*eMass << " BR=" << GetBR() <<  << 203          << "% and endpoint energy " << endpointEnergy/keV << " keV " << G4endl;
196 }                                                 204 }
197                                                   205 
198                                                   206