source: trunk/source/processes/electromagnetic/standard/src/G4CoulombScatteringModel.cc@ 1005

Last change on this file since 1005 was 991, checked in by garnier, 17 years ago

update
File size: 6.8 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: G4CoulombScatteringModel.cc,v 1.37 2008/07/31 13:11:34 vnivanch Exp $
27// GEANT4 tag $Name: geant4-09-02 $
28//
29// -------------------------------------------------------------------
30//
31// GEANT4 Class file
32//
33//
34// File name: G4CoulombScatteringModel
35//
36// Author: Vladimir Ivanchenko
37//
38// Creation date: 22.08.2005
39//
40// Modifications:
41// 01.08.06 V.Ivanchenko extend upper limit of table to TeV and review the
42// logic of building - only elements from G4ElementTable
43// 08.08.06 V.Ivanchenko build internal table in ekin scale, introduce faclim
44// 19.10.06 V.Ivanchenko use inheritance from G4eCoulombScatteringModel
45// 09.10.07 V.Ivanchenko reorganized methods, add cut dependence in scattering off e-
46// 09.06.08 V.Ivanchenko SelectIsotope is moved to the base class
47//
48// Class Description:
49//
50// -------------------------------------------------------------------
51//
52//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
53//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
54
55#include "G4CoulombScatteringModel.hh"
56#include "Randomize.hh"
57#include "G4ParticleChangeForGamma.hh"
58#include "G4ParticleTable.hh"
59#include "G4IonTable.hh"
60#include "G4Proton.hh"
61
62//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
63
64using namespace std;
65
66G4CoulombScatteringModel::G4CoulombScatteringModel(const G4String& nam)
67 : G4eCoulombScatteringModel(nam)
68{}
69
70//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
71
72G4CoulombScatteringModel::~G4CoulombScatteringModel()
73{}
74
75//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
76
77G4double G4CoulombScatteringModel::ComputeCrossSectionPerAtom(
78 const G4ParticleDefinition* p,
79 G4double kinEnergy,
80 G4double Z,
81 G4double,
82 G4double cutEnergy,
83 G4double)
84{
85 SetupParticle(p);
86 G4double ekin = std::max(lowEnergyLimit, kinEnergy);
87 SetupKinematic(ekin, cutEnergy);
88
89 // save lab system kinematics
90 G4double xtkin = tkin;
91 G4double xmom2 = mom2;
92 G4double xinvb = invbeta2;
93
94 // CM system
95 G4int iz = G4int(Z);
96 G4double m2 = fNistManager->GetAtomicMassAmu(iz)*amu_c2;
97 G4double etot = tkin + mass;
98 G4double ptot = sqrt(mom2);
99
100 G4double m12 = mass*mass;
101 G4double momCM= ptot*m2/sqrt(m12 + m2*m2 + 2.0*etot*m2);
102
103 mom2 = momCM*momCM;
104 tkin = sqrt(mom2 + m12) - mass;
105 invbeta2 = 1.0 + m12/mom2;
106
107 SetupTarget(Z, tkin);
108
109 G4double xsec = CrossSectionPerAtom();
110
111 // restore Lab system kinematics
112 tkin = xtkin;
113 mom2 = xmom2;
114 invbeta2 = xinvb;
115
116 return xsec;
117}
118
119//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
120
121void G4CoulombScatteringModel::SampleSecondaries(
122 std::vector<G4DynamicParticle*>* fvect,
123 const G4MaterialCutsCouple* couple,
124 const G4DynamicParticle* dp,
125 G4double cutEnergy,
126 G4double)
127{
128 G4double kinEnergy = dp->GetKineticEnergy();
129 if(kinEnergy <= DBL_MIN) return;
130 DefineMaterial(couple);
131 SetupParticle(dp->GetDefinition());
132 G4double ekin = std::max(lowEnergyLimit, kinEnergy);
133 SetupKinematic(ekin, cutEnergy);
134
135 // Choose nucleus
136 currentElement = SelectRandomAtom(couple,particle,ekin,ecut,tkin);
137
138 G4double Z = currentElement->GetZ();
139 G4int iz = G4int(Z);
140 G4int ia = SelectIsotopeNumber(currentElement);
141 G4double m2 = theParticleTable->GetIonTable()->GetNucleusMass(iz, ia);
142
143 // CM system
144 G4double etot = tkin + mass;
145 G4double ptot = sqrt(mom2);
146
147 G4double momCM= ptot*m2/sqrt(mass*mass + m2*m2 + 2.0*etot*m2);
148 mom2 = momCM*momCM;
149 G4double m12 = mass*mass;
150 G4double eCM = sqrt(mom2 + m12);
151
152 // a correction for heavy projectile
153 G4double fm = m2/(mass + m2);
154 invbeta2 = 1.0 + m12*fm*fm/mom2;
155
156 // sample scattering angle in CM system
157 SetupTarget(Z, eCM - mass);
158
159 G4double cost = SampleCosineTheta();
160 G4double z1 = 1.0 - cost;
161 if(z1 < 0.0) return;
162
163 G4double sint = sqrt(z1*(1.0 + cost));
164 G4double phi = twopi * G4UniformRand();
165
166 // kinematics in the Lab system
167 G4double bet = ptot/(etot + m2);
168 G4double gam = 1.0/sqrt((1.0 - bet)*(1.0 + bet));
169 G4double pzCM = momCM*cost;
170
171 G4ThreeVector v1(momCM*cos(phi)*sint,momCM*sin(phi)*sint,gam*(pzCM + bet*eCM));
172 G4ThreeVector dir = dp->GetMomentumDirection();
173 G4ThreeVector newDirection = v1.unit();
174 newDirection.rotateUz(dir);
175 fParticleChange->ProposeMomentumDirection(newDirection);
176
177 G4double elab = gam*(eCM + bet*pzCM);
178 ekin = elab - mass;
179 if(ekin < 0.0) ekin = 0.0;
180 fParticleChange->SetProposedKineticEnergy(ekin);
181
182 // recoil
183 G4double erec = kinEnergy - ekin;
184 G4double th =
185 std::min(recoilThreshold,
186 Z*currentElement->GetIonisation()->GetMeanExcitationEnergy());
187
188 if(erec > th) {
189 G4ParticleDefinition* ion = theParticleTable->FindIon(iz, ia, 0, iz);
190 G4double plab = sqrt(ekin*(ekin + 2.0*mass));
191 G4ThreeVector p2 = (ptot*dir - plab*newDirection).unit();
192 G4DynamicParticle* newdp = new G4DynamicParticle(ion, p2, erec);
193 fvect->push_back(newdp);
194 } else if(erec > 0.0) {
195 fParticleChange->ProposeLocalEnergyDeposit(erec);
196 fParticleChange->ProposeNonIonizingEnergyDeposit(erec);
197 }
198}
199
200//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
201
202
Note: See TracBrowser for help on using the repository browser.