source: trunk/source/processes/hadronic/stopping/src/G4MuMinusCaptureCascade.cc@ 1196

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

update CVS release candidate geant4.9.3.01
File size: 8.5 KB
RevLine 
[819]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//
[962]26// $Id: G4MuMinusCaptureCascade.cc,v 1.16 2008/05/05 09:09:06 vnivanch Exp $
[1196]27// GEANT4 tag $Name: geant4-09-03-cand-01 $
[819]28//
29// G4MuonMinusCaptureAtRest physics process
30//
31// E-mail: Vladimir.Ivantchenko@cern.ch
32//
33// Created: 02.04.00 V.Ivanchenko
34//
35// Modified:
36// 06.04.01 V.Ivanchenko Bug in theta distribution fixed
37// 13.02.07 V.Ivanchenko Fixes in decay - add random distribution of e-
38// direction; factor 2 in potential energy
39//
40//----------------------------------------------------------------------
41
42#include "G4MuMinusCaptureCascade.hh"
43#include "G4LorentzVector.hh"
44#include "G4ParticleMomentum.hh"
45#include "G4MuonMinus.hh"
46#include "G4Electron.hh"
47#include "G4Gamma.hh"
48#include "G4NeutrinoMu.hh"
49#include "G4AntiNeutrinoE.hh"
50#include "G4GHEKinematicsVector.hh"
51
52//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
53
54G4MuMinusCaptureCascade::G4MuMinusCaptureCascade()
55{
56 theElectron = G4Electron::Electron();
57 theGamma = G4Gamma::Gamma();
58 Emass = theElectron->GetPDGMass();
59 MuMass = G4MuonMinus::MuonMinus()->GetPDGMass();
60}
61
62//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
63
64G4MuMinusCaptureCascade::~G4MuMinusCaptureCascade()
65{ }
66
67//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
68
69G4double G4MuMinusCaptureCascade::GetKShellEnergy(G4double Z)
70{
71 // Calculate the Energy of K Mesoatom Level for this Element using
72 // the Energy of Hydrogen Atom taken into account finite size of the
73 // nucleus (V.Ivanchenko)
74 const G4int ListK = 28;
75 static G4double ListZK[ListK] = {
76 1., 2., 4., 6., 8., 11., 14., 17., 18., 21., 24.,
77 26., 29., 32., 38., 40., 41., 44., 49., 53., 55.,
78 60., 65., 70., 75., 81., 85., 92.};
79 static G4double ListKEnergy[ListK] = {
80 0.00275, 0.011, 0.043, 0.098, 0.173, 0.326,
81 0.524, 0.765, 0.853, 1.146, 1.472,
82 1.708, 2.081, 2.475, 3.323, 3.627,
83 3.779, 4.237, 5.016, 5.647, 5.966,
84 6.793, 7.602, 8.421, 9.249, 10.222,
85 10.923,11.984};
86
87 // Energy with finit size corrections
88 G4double KEnergy = GetLinApprox(ListK,ListZK,ListKEnergy,Z);
89
90 return KEnergy;
91}
92
93//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
94
95void G4MuMinusCaptureCascade::AddNewParticle(G4ParticleDefinition* aParticle,
96 G4ThreeVector& Momentum,
97 G4double mass,
98 G4int* nParticle,
99 G4GHEKinematicsVector* Cascade)
100{
101 // Store particle in the HEK vector and increment counter
102 Cascade[*nParticle].SetZero();
103 Cascade[*nParticle].SetMass( mass );
104 Cascade[*nParticle].SetMomentumAndUpdate(Momentum.x(), Momentum.y(), Momentum.z());
105 Cascade[*nParticle].SetParticleDef( aParticle );
106 (*nParticle)++;
107
108 return;
109}
110
111//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
112
113G4int G4MuMinusCaptureCascade::DoCascade(const G4double Z, const G4double massA,
114 G4GHEKinematicsVector* Cascade)
115{
116 // Inicialization - cascade start from 14th level
117 // N.C.Mukhopadhyay Phy. Rep. 30 (1977) 1.
118 G4int nPart = 0;
119 G4double EnergyLevel[14];
120
121 G4double mass = MuMass * massA / (MuMass + massA) ;
122
123 const G4double KEnergy = 13.6 * eV * Z * Z * mass/ electron_mass_c2;
124
125 EnergyLevel[0] = GetKShellEnergy(Z);
126 for( G4int i = 2; i < 15; i++ ) {
127 EnergyLevel[i-1] = KEnergy / (i*i) ;
128 }
129
130 G4int nElec = G4int(Z);
131 G4int nAuger = 1;
132 G4int nLevel = 13;
133 G4double DeltaE;
134 G4double pGamma = Z*Z*Z*Z;
135
136 // Capture on 14-th level
137 G4double ptot = std::sqrt(EnergyLevel[13]*(EnergyLevel[13] + 2.0*Emass));
138 G4ThreeVector moment = ptot * GetRandomVec();
139
140 AddNewParticle(theElectron,moment,Emass,&nPart,Cascade);
141
142 // Emit new photon or electron
143 // Simplified model for probabilities
144 // N.C.Mukhopadhyay Phy. Rep. 30 (1977) 1.
145 do {
146
147 // case of Auger electrons
148 if((nAuger < nElec) && ((pGamma + 10000.0) * G4UniformRand() < 10000.0) ) {
149 nAuger++;
150 DeltaE = EnergyLevel[nLevel-1] - EnergyLevel[nLevel];
151 nLevel--;
152
153 ptot = std::sqrt(DeltaE * (DeltaE + 2.0*Emass));
154 moment = ptot * GetRandomVec();
155
156 AddNewParticle(theElectron, moment, Emass, &nPart, Cascade);
157
158 } else {
159
160 // Case of photon cascade, probabilities from
161 // C.S.Wu and L.Wilets, Ann. Rev. Nuclear Sci. 19 (1969) 527.
162
163 G4double var = (10.0 + G4double(nLevel - 1) ) * G4UniformRand();
164 G4int iLevel = nLevel - 1 ;
165 if(var > 10.0) iLevel -= G4int(var-10.0) + 1;
166 if( iLevel < 0 ) iLevel = 0;
167 DeltaE = EnergyLevel[iLevel] - EnergyLevel[nLevel];
168 nLevel = iLevel;
169 moment = DeltaE * GetRandomVec();
170 AddNewParticle(theGamma, moment, 0.0, &nPart, Cascade);
171 }
172
173 } while( nLevel > 0 );
174
175 return nPart;
176}
177
178//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
179
180void G4MuMinusCaptureCascade::DoBoundMuonMinusDecay(G4double Z,
181 G4int* nCascade,
182 G4GHEKinematicsVector* Cascade)
183{
184 // Simulation on Decay of mu- on a K-shell of the muonic atom
185 G4double xmax = ( 1.0 + Emass*Emass/ (MuMass*MuMass) );
186 G4double xmin = 2.0*Emass/MuMass;
187 G4double KEnergy = GetKShellEnergy(Z);
188 /*
189 G4cout << "G4MuMinusCaptureCascade::DoBoundMuonMinusDecay"
190 << " XMAX= " << xmax
191 << " Ebound= " << KEnergy
192 << G4endl;
193 */
194 G4double pmu = std::sqrt(KEnergy*(KEnergy + 2.0*MuMass));
195 G4double emu = KEnergy + MuMass;
196 G4ThreeVector moment = GetRandomVec();
197 G4LorentzVector MU(pmu*moment,emu);
198 G4ThreeVector bst = MU.boostVector();
199
200 G4double Eelect, Pelect, x, ecm;
201 G4LorentzVector EL, NN;
202 // Calculate electron energy
203 do {
204 do {
205 x = xmin + (xmax-xmin)*G4UniformRand();
206 } while (G4UniformRand() > (3.0 - 2.0*x)*x*x );
207 Eelect = x*MuMass*0.5;
208 Pelect = 0.0;
209 if(Eelect > Emass) {
210 Pelect = std::sqrt( Eelect*Eelect - Emass*Emass );
211 } else {
212 Pelect = 0.0;
213 Eelect = Emass;
214 }
215 G4ThreeVector e_mom = GetRandomVec();
216 EL = G4LorentzVector(Pelect*e_mom,Eelect);
217 EL.boost(bst);
218 Eelect = EL.e() - Emass - 2.0*KEnergy;
219 //
220 // Calculate rest frame parameters of 2 neutrinos
221 //
222 NN = MU - EL;
223 ecm = NN.mag2();
224 } while (Eelect < 0.0 || ecm < 0.0);
225
226 //
227 // Create electron
228 //
229 moment = std::sqrt(Eelect * (Eelect + 2.0*Emass))*(EL.vect().unit());
230 AddNewParticle(theElectron, moment, Emass, nCascade, Cascade);
231 //
232 // Create Neutrinos
233 //
234 ecm = 0.5*std::sqrt(ecm);
235 bst = NN.boostVector();
236 G4ThreeVector p1 = ecm * GetRandomVec();
237 G4LorentzVector N1 = G4LorentzVector(p1,ecm);
238 N1.boost(bst);
239 G4ThreeVector p1lab = N1.vect();
240 AddNewParticle(G4AntiNeutrinoE::AntiNeutrinoE(),p1lab,0.0,nCascade,Cascade);
241 NN -= N1;
242 G4ThreeVector p2lab = NN.vect();
243 AddNewParticle(G4NeutrinoMu::NeutrinoMu(),p2lab,0.0,nCascade,Cascade);
244
245 return;
246}
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
Note: See TracBrowser for help on using the repository browser.