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

Last change on this file since 1337 was 1337, checked in by garnier, 14 years ago

tag geant4.9.4 beta 1 + modifs locales
File size: 8.5 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// $Id: G4MuMinusCaptureCascade.cc,v 1.16 2008/05/05 09:09:06 vnivanch Exp $
27// GEANT4 tag $Name: geant4-09-04-beta-01 $
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.