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

Last change on this file since 900 was 819, checked in by garnier, 17 years ago

import all except CVS
File size: 9.1 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.29 2007/11/09 11:45:45 vnivanch Exp $
27// GEANT4 tag $Name: geant4-09-01-patch-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//
47// Class Description:
48//
49// -------------------------------------------------------------------
50//
51//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
52//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
53
54#include "G4CoulombScatteringModel.hh"
55#include "Randomize.hh"
56#include "G4ParticleChangeForGamma.hh"
57#include "G4NistManager.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(
67 G4double thetaMin, G4double thetaMax, G4bool build,
68 G4double tlim, const G4String& nam)
69 : G4eCoulombScatteringModel(thetaMin,thetaMax,build,tlim,nam)
70{
71 theMatManager = G4NistManager::Instance();
72 theParticleTable = G4ParticleTable::GetParticleTable();
73}
74
75//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
76
77G4CoulombScatteringModel::~G4CoulombScatteringModel()
78{}
79
80//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
81
82G4double G4CoulombScatteringModel::ComputeCrossSectionPerAtom(
83 const G4ParticleDefinition* p,
84 G4double kinEnergy,
85 G4double Z,
86 G4double A,
87 G4double cutEnergy,
88 G4double)
89{
90 if(p == particle && kinEnergy == tkin && Z == targetZ &&
91 A == targetA && cutEnergy == ecut) return nucXSection;
92
93 // Lab system
94 G4double ekin = std::max(keV, kinEnergy);
95 nucXSection = ComputeElectronXSectionPerAtom(p,ekin,Z,A,cutEnergy);
96
97 // CM system
98 G4int iz = G4int(Z);
99 G4double m1 = theMatManager->GetAtomicMassAmu(iz)*amu_c2;
100 G4double etot = tkin + mass;
101 G4double ptot = sqrt(mom2);
102 G4double bet = ptot/(etot + m1);
103 G4double gam = 1.0/sqrt((1.0 - bet)*(1.0 + bet));
104 G4double momCM= gam*(ptot - bet*etot);
105
106 // G4cout << "ptot= " << ptot << " etot= " << etot << " beta= "
107 // << bet << " gam= " << gam << " Z= " << Z << " A= " << A << G4endl;
108 // G4cout << " CM. mom= " << momCM << " m= " << m
109 // << " m1= " << m1 << " iz= " << iz <<G4endl;
110
111 G4double momCM2 = momCM*momCM;
112 cosTetMaxNuc = std::max(cosThetaMax, 1.0 - 0.5*q2Limit/momCM2);
113 if(1.5 > targetA && p == theProton && cosTetMaxNuc < 0.0) cosTetMaxNuc = 0.0;
114 //G4cout << " ctmax= " << cosTetMaxNuc
115 //<< " ctmin= " << cosThetaMin << G4endl;
116
117 // Cross section in CM system
118 if(cosTetMaxNuc < cosThetaMin) {
119 G4double effmass = mass*m1/(mass + m1);
120 G4double x1 = 1.0 - cosThetaMin;
121 G4double x2 = 1.0 - cosTetMaxNuc;
122 G4double z1 = x1 + screenZ;
123 G4double z2 = x2 + screenZ;
124 G4double d = 1.0/formfactA;
125 G4double zn1= x1 + d;
126 G4double zn2= x2 + d;
127 nucXSection += coeff*Z*Z*chargeSquare*(1.0 + effmass*effmass/momCM2)
128 *(1./z1 - 1./z2 + 1./zn1 - 1./zn2 +
129 2.0*formfactA*std::log(z1*zn2/(z2*zn1)))/momCM2;
130 //G4cout << "XS: x1= " << x1 << " x2= " << x2
131 //<< " cross= " << cross << G4endl;
132 //G4cout << "momCM2= " << momCM2 << " invbeta2= " << invbeta2
133 // << " coeff= " << coeff << G4endl;
134 }
135 if(nucXSection < 0.0) nucXSection = 0.0;
136 return nucXSection;
137}
138
139//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
140
141G4double G4CoulombScatteringModel::SelectIsotope(const G4Element* elm)
142{
143 G4double N = elm->GetN();
144 G4int ni = elm->GetNumberOfIsotopes();
145 if(ni > 0) {
146 G4double* ab = elm->GetRelativeAbundanceVector();
147 G4double x = G4UniformRand();
148 G4int idx;
149 for(idx=0; idx<ni; idx++) {
150 x -= ab[idx];
151 if (x <= 0.0) break;
152 }
153 if(idx >= ni) {
154 G4cout << "G4CoulombScatteringModel::SelectIsotope WARNING: "
155 << "abandance vector for"
156 << elm->GetName() << " is not normalised to unit" << G4endl;
157 } else {
158 N = G4double(elm->GetIsotope(idx)->GetN());
159 }
160 }
161 return N;
162}
163
164//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
165
166void G4CoulombScatteringModel::SampleSecondaries(
167 std::vector<G4DynamicParticle*>* fvect,
168 const G4MaterialCutsCouple* couple,
169 const G4DynamicParticle* dp,
170 G4double cutEnergy,
171 G4double maxEnergy)
172{
173 const G4Material* aMaterial = couple->GetMaterial();
174 const G4ParticleDefinition* p = dp->GetDefinition();
175 G4double kinEnergy = dp->GetKineticEnergy();
176
177 // Select isotope and setup
178 SetupParticle(p);
179 const G4Element* elm =
180 SelectRandomAtom(aMaterial,p,kinEnergy,cutEnergy,maxEnergy);
181 G4double Z = elm->GetZ();
182 G4double A = SelectIsotope(elm);
183 G4int iz = G4int(Z);
184 G4int ia = G4int(A + 0.5);
185
186 G4double cross =
187 ComputeCrossSectionPerAtom(p,kinEnergy,Z,A,cutEnergy,maxEnergy);
188
189 G4double costm = cosTetMaxNuc;
190 G4double formf = formfactA;
191 if(G4UniformRand()*cross < elecXSection) {
192 costm = cosTetMaxElec;
193 formf = 0.0;
194 }
195
196 // G4cout << "SampleSec: Ekin= " << kinEnergy << " m1= " << m1
197 // << " Z= "<< Z << " A= " <<A<< G4endl;
198
199 if(costm >= cosThetaMin) return;
200
201 // kinematics in CM system
202 G4double m1 = theParticleTable->GetIonTable()->GetNucleusMass(iz, ia);
203 G4double etot = kinEnergy + mass;
204 G4double ptot = sqrt(mom2);
205 G4double bet = ptot/(etot + m1);
206 G4double gam = 1.0/sqrt((1.0 - bet)*(1.0 + bet));
207 G4double pCM = gam*(ptot - bet*etot);
208 G4double eCM = gam*(etot - bet*ptot);
209
210 G4double x1 = 1. - cosThetaMin + screenZ;
211 G4double x2 = 1. - costm;
212 G4double x3 = cosThetaMin - costm;
213
214 G4double grej, z, z1;
215 do {
216 z = G4UniformRand()*x3;
217 z1 = (x1*x2 - screenZ*z)/(x1 + z);
218 if(z1 < 0.0) z1 = 0.0;
219 else if(z1 > 2.0) z1 = 2.0;
220 grej = 1.0/(1.0 + formf*z1);
221 } while ( G4UniformRand() > grej*grej );
222
223 G4double cost = 1.0 - z1;
224 G4double sint= sqrt(z1*(2.0 - z1));
225
226 G4double phi = twopi * G4UniformRand();
227
228 // projectile after scattering
229 G4double pzCM = pCM*cost;
230 G4ThreeVector v1(pCM*cos(phi)*sint,pCM*sin(phi)*sint,gam*(pzCM + bet*eCM));
231 G4ThreeVector dir = dp->GetMomentumDirection();
232 G4ThreeVector newDirection = v1.unit();
233 newDirection.rotateUz(dir);
234 fParticleChange->ProposeMomentumDirection(newDirection);
235 G4double elab = gam*(eCM + bet*pzCM);
236 G4double ekin = elab - mass;
237 if(ekin < 0.0) ekin = 0.0;
238 G4double plab = sqrt(ekin*(ekin + 2.0*mass));
239 fParticleChange->SetProposedKineticEnergy(ekin);
240
241 // recoil
242 G4double erec = kinEnergy - ekin;
243 if(erec > Z*aMaterial->GetIonisation()->GetMeanExcitationEnergy()) {
244 G4ParticleDefinition* ion = theParticleTable->FindIon(iz, ia, 0, iz);
245 G4ThreeVector p2 = (ptot*dir - plab*newDirection).unit();
246 G4DynamicParticle* newdp = new G4DynamicParticle(ion, p2, erec);
247 fvect->push_back(newdp);
248 } else if(erec > 0.0) {
249 fParticleChange->ProposeLocalEnergyDeposit(erec);
250 fParticleChange->ProposeNonIonizingEnergyDeposit(erec);
251 }
252}
253
254//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
255
256
Note: See TracBrowser for help on using the repository browser.