source: trunk/source/processes/electromagnetic/standard/src/G4eCoulombScatteringModel.cc@ 1036

Last change on this file since 1036 was 1007, checked in by garnier, 17 years ago

update to geant4.9.2
File size: 11.6 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: G4eCoulombScatteringModel.cc,v 1.59 2008/10/22 18:39:29 vnivanch Exp $
27// GEANT4 tag $Name: geant4-09-02 $
28//
29// -------------------------------------------------------------------
30//
31// GEANT4 Class file
32//
33//
34// File name: G4eCoulombScatteringModel
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.08.06 V.Ivanchenko add inline function ScreeningParameter
45// 09.10.07 V.Ivanchenko reorganized methods, add cut dependence in scattering off e-
46// 09.06.08 V.Ivanchenko add SelectIsotope and sampling of the recoil ion
47//
48// Class Description:
49//
50// -------------------------------------------------------------------
51//
52//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
53//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
54
55#include "G4eCoulombScatteringModel.hh"
56#include "Randomize.hh"
57#include "G4DataVector.hh"
58#include "G4ElementTable.hh"
59#include "G4PhysicsLogVector.hh"
60#include "G4ParticleChangeForGamma.hh"
61#include "G4Electron.hh"
62#include "G4Positron.hh"
63#include "G4Proton.hh"
64#include "G4ParticleTable.hh"
65
66//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
67
68using namespace std;
69
70G4eCoulombScatteringModel::G4eCoulombScatteringModel(const G4String& nam)
71 : G4VEmModel(nam),
72 cosThetaMin(1.0),
73 cosThetaMax(-1.0),
74 q2Limit(TeV*TeV),
75 alpha2(fine_structure_const*fine_structure_const),
76 faclim(100.0),
77 isInitialised(false)
78{
79 fNistManager = G4NistManager::Instance();
80 theParticleTable = G4ParticleTable::GetParticleTable();
81 theElectron = G4Electron::Electron();
82 thePositron = G4Positron::Positron();
83 theProton = G4Proton::Proton();
84 currentMaterial = 0;
85 currentElement = 0;
86 a0 = alpha2*electron_mass_c2*electron_mass_c2/(0.885*0.885);
87 G4double p0 = electron_mass_c2*classic_electr_radius;
88 coeff = twopi*p0*p0;
89 constn = 6.937e-6/(MeV*MeV);
90 tkin = targetZ = mom2 = DBL_MIN;
91 elecXSection = nucXSection = 0.0;
92 recoilThreshold = DBL_MAX;
93 ecut = DBL_MAX;
94 particle = 0;
95 currentCouple = 0;
96 for(size_t j=0; j<100; j++) {
97 FF[j] = 0.0;
98 }
99}
100
101//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
102
103G4eCoulombScatteringModel::~G4eCoulombScatteringModel()
104{}
105
106//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
107
108void G4eCoulombScatteringModel::Initialise(const G4ParticleDefinition* p,
109 const G4DataVector& cuts)
110{
111 SetupParticle(p);
112 currentCouple = 0;
113 elecXSection = nucXSection = 0.0;
114 tkin = targetZ = mom2 = DBL_MIN;
115 ecut = etag = DBL_MAX;
116 cosThetaMin = cos(PolarAngleLimit());
117 currentCuts = &cuts;
118 //G4cout << "!!! G4eCoulombScatteringModel::Initialise for "
119 // << p->GetParticleName() << " cos(TetMin)= " << cosThetaMin
120 // << " cos(TetMax)= " << cosThetaMax <<G4endl;
121 if(!isInitialised) {
122 isInitialised = true;
123
124 if(pParticleChange)
125 fParticleChange =
126 reinterpret_cast<G4ParticleChangeForGamma*>(pParticleChange);
127 else
128 fParticleChange = new G4ParticleChangeForGamma();
129 }
130 if(mass < GeV && particle->GetParticleType() != "nucleus") {
131 InitialiseElementSelectors(p,cuts);
132 }
133}
134
135//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
136
137void G4eCoulombScatteringModel::ComputeMaxElectronScattering(G4double cutEnergy)
138{
139 ecut = cutEnergy;
140 G4double tmax = tkin;
141 cosTetMaxElec = 1.0;
142 if(mass > MeV) {
143 G4double ratio = electron_mass_c2/mass;
144 G4double tau = tkin/mass;
145 tmax = 2.0*electron_mass_c2*tau*(tau + 2.)/
146 (1.0 + 2.0*ratio*(tau + 1.0) + ratio*ratio);
147 cosTetMaxElec = 1.0 - std::min(cutEnergy, tmax)*electron_mass_c2/mom2;
148 } else {
149
150 if(particle == theElectron) tmax *= 0.5;
151 G4double t = std::min(cutEnergy, tmax);
152 G4double mom21 = t*(t + 2.0*electron_mass_c2);
153 G4double t1 = tkin - t;
154 //G4cout << "tkin= " << tkin << " t= " << t << " t1= " << t1 << G4endl;
155 if(t1 > 0.0) {
156 G4double mom22 = t1*(t1 + 2.0*mass);
157 G4double ctm = (mom2 + mom22 - mom21)*0.5/sqrt(mom2*mom22);
158 //G4cout << "ctm= " << ctm << G4endl;
159 if(ctm < 1.0) cosTetMaxElec = ctm;
160 }
161 }
162}
163
164//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
165
166G4double G4eCoulombScatteringModel::ComputeCrossSectionPerAtom(
167 const G4ParticleDefinition* p,
168 G4double kinEnergy,
169 G4double Z, G4double,
170 G4double cutEnergy, G4double)
171{
172 //G4cout << "### G4eCoulombScatteringModel::ComputeCrossSectionPerAtom for "
173 // << p->GetParticleName()<<" Z= "<<Z<<" e(MeV)= "<< kinEnergy/MeV << G4endl;
174 G4double xsec = 0.0;
175 SetupParticle(p);
176 G4double ekin = std::max(lowEnergyLimit, kinEnergy);
177 SetupKinematic(ekin, cutEnergy);
178 if(cosTetMaxNuc < cosTetMinNuc) {
179 SetupTarget(Z, ekin);
180 xsec = CrossSectionPerAtom();
181 }
182 /*
183 G4cout << "e(MeV)= " << ekin/MeV << "cosTetMinNuc= " << cosTetMinNuc
184 << " cosTetMaxNuc= " << cosTetMaxNuc
185 << " cosTetMaxElec= " << cosTetMaxElec
186 << " screenZ= " << screenZ
187 << " formfactA= " << formfactA
188 << " cosTetMaxHad= " << cosTetMaxHad << G4endl;
189 */
190 return xsec;
191}
192
193//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
194
195G4double G4eCoulombScatteringModel::CrossSectionPerAtom()
196{
197 // This method needs initialisation before be called
198
199 G4double fac = coeff*targetZ*chargeSquare*invbeta2/mom2;
200 elecXSection = 0.0;
201 nucXSection = 0.0;
202
203 G4double x = 1.0 - cosTetMinNuc;
204 G4double x1 = x + screenZ;
205
206 if(cosTetMaxElec2 < cosTetMinNuc) {
207 elecXSection = fac*(cosTetMinNuc - cosTetMaxElec2)/
208 (x1*(1.0 - cosTetMaxElec2 + screenZ));
209 nucXSection = elecXSection;
210 }
211
212 //G4cout << "XS tkin(MeV)= " << tkin<<" xs= " <<nucXSection
213 // << " costmax= " << cosTetMaxNuc2
214 // << " costmin= " << cosTetMinNuc << " Z= " << targetZ <<G4endl;
215 if(cosTetMaxNuc2 < cosTetMinNuc) {
216 G4double s = screenZ*formfactA;
217 G4double z1 = 1.0 - cosTetMaxNuc2 + screenZ;
218 G4double d = (1.0 - s)/formfactA;
219 //G4cout <<"x1= "<<x1<<" z1= " <<z1<<" s= "<<s << " d= " <<d <<G4endl;
220 if(d < 0.2*x1) {
221 G4double x2 = x1*x1;
222 G4double z2 = z1*z1;
223 x = (1.0/(x1*x2) - 1.0/(z1*z2) - d*1.5*(1.0/(x2*x2) - 1.0/(z2*z2)))/
224 (3.0*formfactA*formfactA);
225 } else {
226 G4double x2 = x1 + d;
227 G4double z2 = z1 + d;
228 x = (1.0 + 2.0*s)*((cosTetMinNuc - cosTetMaxNuc2)*(1.0/(x1*z1) + 1.0/(x2*z2)) -
229 2.0*log(z1*x2/(z2*x1))/d);
230 }
231 nucXSection += fac*targetZ*x;
232 }
233
234 //G4cout<<" cross(bn)= "<<nucXSection/barn<<" xsElec(bn)= "<<elecXSection/barn
235 // << " Asc= " << screenZ << G4endl;
236
237 return nucXSection;
238}
239
240//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
241
242void G4eCoulombScatteringModel::SampleSecondaries(
243 std::vector<G4DynamicParticle*>* fvect,
244 const G4MaterialCutsCouple* couple,
245 const G4DynamicParticle* dp,
246 G4double cutEnergy,
247 G4double)
248{
249 G4double kinEnergy = dp->GetKineticEnergy();
250 if(kinEnergy <= DBL_MIN) return;
251 DefineMaterial(couple);
252 SetupParticle(dp->GetDefinition());
253 G4double ekin = std::max(lowEnergyLimit, kinEnergy);
254 SetupKinematic(ekin, cutEnergy);
255 //G4cout << "G4eCoulombScatteringModel::SampleSecondaries e(MeV)= "
256 // << kinEnergy << " " << particle->GetParticleName() << G4endl;
257
258 // Choose nucleus
259 currentElement = SelectRandomAtom(couple,particle,ekin,cutEnergy,ekin);
260
261 SetupTarget(currentElement->GetZ(),ekin);
262
263 G4double cost = SampleCosineTheta();
264 G4double z1 = 1.0 - cost;
265 if(z1 < 0.0) return;
266
267 G4double sint = sqrt(z1*(1.0 + cost));
268
269 //G4cout<<"## Sampled sint= " << sint << " Z= " << targetZ
270 // << " screenZ= " << screenZ << " cn= " << formfactA << G4endl;
271
272 G4double phi = twopi * G4UniformRand();
273
274 G4ThreeVector direction = dp->GetMomentumDirection();
275 G4ThreeVector newDirection(cos(phi)*sint,sin(phi)*sint,cost);
276 newDirection.rotateUz(direction);
277
278 fParticleChange->ProposeMomentumDirection(newDirection);
279
280 // recoil sampling assuming a small recoil
281 // and first order correction to primary 4-momentum
282 if(lowEnergyLimit < kinEnergy) {
283 G4int ia = SelectIsotopeNumber(currentElement);
284 G4double Trec = z1*mom2/(amu_c2*G4double(ia));
285 G4double th =
286 std::min(recoilThreshold,
287 targetZ*currentElement->GetIonisation()->GetMeanExcitationEnergy());
288
289 if(Trec > th) {
290 G4int iz = G4int(targetZ);
291 G4ParticleDefinition* ion = theParticleTable->FindIon(iz, ia, 0, iz);
292 Trec = z1*mom2/ion->GetPDGMass();
293 if(Trec < kinEnergy) {
294 G4ThreeVector dir = (direction - newDirection).unit();
295 G4DynamicParticle* newdp = new G4DynamicParticle(ion, dir, Trec);
296 fvect->push_back(newdp);
297 fParticleChange->SetProposedKineticEnergy(kinEnergy - Trec);
298 }
299 }
300 }
301
302 return;
303}
304
305//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
306
307G4double G4eCoulombScatteringModel::SampleCosineTheta()
308{
309 G4double costm = cosTetMaxNuc2;
310 G4double formf = formfactA;
311 G4double prob = 0.0;
312 G4double xs = CrossSectionPerAtom();
313 if(xs > 0.0) prob = elecXSection/xs;
314
315 // scattering off e or A?
316 if(G4UniformRand() < prob) {
317 costm = cosTetMaxElec2;
318 formf = 0.0;
319 }
320
321 /*
322 G4cout << "SampleCost: e(MeV)= " << tkin
323 << " ctmin= " << cosThetaMin
324 << " ctmaxN= " << cosTetMaxNuc
325 << " ctmax= " << costm
326 << " Z= " << targetZ << " A= " << targetA
327 << G4endl;
328 */
329 if(costm >= cosTetMinNuc) return 2.0;
330
331 G4double x1 = 1. - cosTetMinNuc + screenZ;
332 G4double x2 = 1. - costm + screenZ;
333 G4double x3 = cosTetMinNuc - costm;
334 G4double grej, z1;
335 do {
336 z1 = x1*x2/(x1 + G4UniformRand()*x3) - screenZ;
337 grej = 1.0/(1.0 + formf*z1);
338 } while ( G4UniformRand() > grej*grej );
339
340 //G4cout << "z= " << z1 << " cross= " << nucXSection/barn
341 // << " crossE= " << elecXSection/barn << G4endl;
342
343 return 1.0 - z1;
344}
345
346//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
347
348
Note: See TracBrowser for help on using the repository browser.