source: trunk/source/processes/electromagnetic/muons/src/G4MuPairProductionModel.cc@ 1006

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

fichier oublies
File size: 22.0 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: G4MuPairProductionModel.cc,v 1.40 2009/02/20 14:48:16 vnivanch Exp $
27// GEANT4 tag $Name: geant4-09-02-ref-02 $
28//
29// -------------------------------------------------------------------
30//
31// GEANT4 Class file
32//
33//
34// File name: G4MuPairProductionModel
35//
36// Author: Vladimir Ivanchenko on base of Laszlo Urban code
37//
38// Creation date: 24.06.2002
39//
40// Modifications:
41//
42// 04-12-02 Change G4DynamicParticle constructor in PostStep (V.Ivanchenko)
43// 23-12-02 Change interface in order to move to cut per region (V.Ivanchenko)
44// 24-01-03 Fix for compounds (V.Ivanchenko)
45// 27-01-03 Make models region aware (V.Ivanchenko)
46// 13-02-03 Add model (V.Ivanchenko)
47// 06-06-03 Fix in cross section calculation for high energy (V.Ivanchenko)
48// 20-10-03 2*xi in ComputeDDMicroscopicCrossSection (R.Kokoulin)
49// 8 integration points in ComputeDMicroscopicCrossSection
50// 12-01-04 Take min cut of e- and e+ not its sum (V.Ivanchenko)
51// 10-02-04 Update parameterisation using R.Kokoulin model (V.Ivanchenko)
52// 28-04-04 For complex materials repeat calculation of max energy for each
53// material (V.Ivanchenko)
54// 01-11-04 Fix bug inside ComputeDMicroscopicCrossSection (R.Kokoulin)
55// 08-04-05 Major optimisation of internal interfaces (V.Ivantchenko)
56// 03-08-05 Add SetParticle method (V.Ivantchenko)
57// 23-10-05 Add protection in sampling of e+e- pair energy needed for
58// low cuts (V.Ivantchenko)
59// 13-02-06 Add ComputeCrossSectionPerAtom (mma)
60// 24-04-07 Add protection in SelectRandomAtom method (V.Ivantchenko)
61// 12-05-06 Updated sampling (use cut) in SelectRandomAtom (A.Bogdanov)
62// 11-10-07 Add ignoreCut flag (V.Ivanchenko)
63
64//
65// Class Description:
66//
67//
68// -------------------------------------------------------------------
69//
70//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
71//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
72
73#include "G4MuPairProductionModel.hh"
74#include "G4Electron.hh"
75#include "G4Positron.hh"
76#include "G4MuonMinus.hh"
77#include "G4MuonPlus.hh"
78#include "Randomize.hh"
79#include "G4Material.hh"
80#include "G4Element.hh"
81#include "G4ElementVector.hh"
82#include "G4ProductionCutsTable.hh"
83#include "G4ParticleChangeForLoss.hh"
84#include "G4ParticleChangeForGamma.hh"
85
86//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
87
88// static members
89//
90G4double G4MuPairProductionModel::zdat[]={1., 4., 13., 29., 92.};
91G4double G4MuPairProductionModel::adat[]={1.01, 9.01, 26.98, 63.55, 238.03};
92G4double G4MuPairProductionModel::tdat[]={1.e3, 1.e4, 1.e5, 1.e6, 1.e7, 1.e8,
93 1.e9, 1.e10};
94G4double G4MuPairProductionModel::xgi[]={ 0.0199, 0.1017, 0.2372, 0.4083,
95 0.5917, 0.7628, 0.8983, 0.9801 };
96G4double G4MuPairProductionModel::wgi[]={ 0.0506, 0.1112, 0.1569, 0.1813,
97 0.1813, 0.1569, 0.1112, 0.0506 };
98
99//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
100
101using namespace std;
102
103G4MuPairProductionModel::G4MuPairProductionModel(const G4ParticleDefinition* p,
104 const G4String& nam)
105 : G4VEmModel(nam),
106 particle(0),
107 factorForCross(4.*fine_structure_const*fine_structure_const
108 *classic_electr_radius*classic_electr_radius/(3.*pi)),
109 sqrte(sqrt(exp(1.))),
110 currentZ(0),
111 fParticleChange(0),
112 minPairEnergy(4.*electron_mass_c2),
113 lowestKinEnergy(1.*GeV),
114 nzdat(5),
115 ntdat(8),
116 nbiny(1000),
117 nmaxElements(0),
118 ymin(-5.),
119 ymax(0.),
120 dy((ymax-ymin)/nbiny),
121 samplingTablesAreFilled(false)
122{
123 SetLowEnergyLimit(minPairEnergy);
124 nist = G4NistManager::Instance();
125
126 theElectron = G4Electron::Electron();
127 thePositron = G4Positron::Positron();
128
129 if(p) SetParticle(p);
130}
131
132//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
133
134G4MuPairProductionModel::~G4MuPairProductionModel()
135{}
136
137//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
138
139G4double G4MuPairProductionModel::MinEnergyCut(const G4ParticleDefinition*,
140 const G4MaterialCutsCouple* )
141{
142 return minPairEnergy;
143}
144
145//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
146
147G4double G4MuPairProductionModel::MaxSecondaryEnergy(const G4ParticleDefinition*,
148 G4double kineticEnergy)
149{
150 G4double maxPairEnergy = kineticEnergy + particleMass*(1.0 - 0.75*sqrte*z13);
151 return maxPairEnergy;
152}
153
154//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
155
156void G4MuPairProductionModel::Initialise(const G4ParticleDefinition* p,
157 const G4DataVector&)
158{
159 if (!samplingTablesAreFilled) {
160 if(p) SetParticle(p);
161 MakeSamplingTables();
162 }
163 if(!fParticleChange) {
164 if(pParticleChange)
165 fParticleChange =
166 reinterpret_cast<G4ParticleChangeForLoss*>(pParticleChange);
167 else
168 fParticleChange = new G4ParticleChangeForLoss();
169 }
170}
171
172//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
173
174G4double G4MuPairProductionModel::ComputeDEDXPerVolume(
175 const G4Material* material,
176 const G4ParticleDefinition*,
177 G4double kineticEnergy,
178 G4double cutEnergy)
179{
180 G4double dedx = 0.0;
181 if (cutEnergy <= minPairEnergy || kineticEnergy <= lowestKinEnergy)
182 return dedx;
183
184 const G4ElementVector* theElementVector = material->GetElementVector();
185 const G4double* theAtomicNumDensityVector =
186 material->GetAtomicNumDensityVector();
187
188 // loop for elements in the material
189 for (size_t i=0; i<material->GetNumberOfElements(); i++) {
190 G4double Z = (*theElementVector)[i]->GetZ();
191 SetCurrentElement(Z);
192 G4double tmax = MaxSecondaryEnergy(particle, kineticEnergy);
193 G4double loss = ComputMuPairLoss(Z, kineticEnergy, cutEnergy, tmax);
194 dedx += loss*theAtomicNumDensityVector[i];
195 }
196 if (dedx < 0.) dedx = 0.;
197 return dedx;
198}
199
200//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
201
202G4double G4MuPairProductionModel::ComputMuPairLoss(G4double Z,
203 G4double tkin,
204 G4double cutEnergy,
205 G4double tmax)
206{
207 SetCurrentElement(Z);
208 G4double loss = 0.0;
209
210 G4double cut = std::min(cutEnergy,tmax);
211 if(cut <= minPairEnergy) return loss;
212
213 // calculate the rectricted loss
214 // numerical integration in log(PairEnergy)
215 G4double ak1=6.9;
216 G4double ak2=1.0;
217 G4double aaa = log(minPairEnergy);
218 G4double bbb = log(cut);
219 G4int kkk = (G4int)((bbb-aaa)/ak1+ak2);
220 if (kkk > 8) kkk = 8;
221 G4double hhh = (bbb-aaa)/(G4double)kkk;
222 G4double x = aaa;
223
224 for (G4int l=0 ; l<kkk; l++)
225 {
226
227 for (G4int ll=0; ll<8; ll++)
228 {
229 G4double ep = exp(x+xgi[ll]*hhh);
230 loss += wgi[ll]*ep*ep*ComputeDMicroscopicCrossSection(tkin, Z, ep);
231 }
232 x += hhh;
233 }
234 loss *= hhh;
235 if (loss < 0.) loss = 0.;
236 return loss;
237}
238
239//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
240
241G4double G4MuPairProductionModel::ComputeMicroscopicCrossSection(
242 G4double tkin,
243 G4double Z,
244 G4double cut)
245{
246 G4double cross = 0.;
247 SetCurrentElement(Z);
248 G4double tmax = MaxSecondaryEnergy(particle, tkin);
249 if (tmax <= cut) return cross;
250
251 G4double ak1=6.9 ;
252 G4double ak2=1.0 ;
253 G4double aaa = log(cut);
254 G4double bbb = log(tmax);
255 G4int kkk = (G4int)((bbb-aaa)/ak1 + ak2);
256 if(kkk > 8) kkk = 8;
257 G4double hhh = (bbb-aaa)/float(kkk);
258 G4double x = aaa;
259
260 for(G4int l=0; l<kkk; l++)
261 {
262 for(G4int i=0; i<8; i++)
263 {
264 G4double ep = exp(x + xgi[i]*hhh);
265 cross += ep*wgi[i]*ComputeDMicroscopicCrossSection(tkin, Z, ep);
266 }
267 x += hhh;
268 }
269
270 cross *=hhh;
271 if(cross < 0.0) cross = 0.0;
272 return cross;
273}
274
275G4double G4MuPairProductionModel::ComputeDMicroscopicCrossSection(
276 G4double tkin,
277 G4double Z,
278 G4double pairEnergy)
279 // Calculates the differential (D) microscopic cross section
280 // using the cross section formula of R.P. Kokoulin (18/01/98)
281 // Code modified by R.P. Kokoulin, V.N. Ivanchenko (27/01/04)
282{
283 G4double bbbtf= 183. ;
284 G4double bbbh = 202.4 ;
285 G4double g1tf = 1.95e-5 ;
286 G4double g2tf = 5.3e-5 ;
287 G4double g1h = 4.4e-5 ;
288 G4double g2h = 4.8e-5 ;
289
290 G4double totalEnergy = tkin + particleMass;
291 G4double residEnergy = totalEnergy - pairEnergy;
292 G4double massratio = particleMass/electron_mass_c2 ;
293 G4double massratio2 = massratio*massratio ;
294 G4double cross = 0.;
295
296 SetCurrentElement(Z);
297
298 G4double c3 = 0.75*sqrte*particleMass;
299 if (residEnergy <= c3*z13) return cross;
300
301 G4double c7 = 4.*electron_mass_c2;
302 G4double c8 = 6.*particleMass*particleMass;
303 G4double alf = c7/pairEnergy;
304 G4double a3 = 1. - alf;
305 if (a3 <= 0.) return cross;
306
307 // zeta calculation
308 G4double bbb,g1,g2;
309 if( Z < 1.5 ) { bbb = bbbh ; g1 = g1h ; g2 = g2h ; }
310 else { bbb = bbbtf; g1 = g1tf; g2 = g2tf; }
311
312 G4double zeta = 0;
313 G4double zeta1 = 0.073*log(totalEnergy/(particleMass+g1*z23*totalEnergy))-0.26;
314 if ( zeta1 > 0.)
315 {
316 G4double zeta2 = 0.058*log(totalEnergy/(particleMass+g2*z13*totalEnergy))-0.14;
317 zeta = zeta1/zeta2 ;
318 }
319
320 G4double z2 = Z*(Z+zeta);
321 G4double screen0 = 2.*electron_mass_c2*sqrte*bbb/(z13*pairEnergy);
322 G4double a0 = totalEnergy*residEnergy;
323 G4double a1 = pairEnergy*pairEnergy/a0;
324 G4double bet = 0.5*a1;
325 G4double xi0 = 0.25*massratio2*a1;
326 G4double del = c8/a0;
327
328 G4double rta3 = sqrt(a3);
329 G4double tmnexp = alf/(1. + rta3) + del*rta3;
330 if(tmnexp >= 1.0) return cross;
331
332 G4double tmn = log(tmnexp);
333 G4double sum = 0.;
334
335 // Gaussian integration in ln(1-ro) ( with 8 points)
336 for (G4int i=0; i<8; i++)
337 {
338 G4double a4 = exp(tmn*xgi[i]); // a4 = (1.-asymmetry)
339 G4double a5 = a4*(2.-a4) ;
340 G4double a6 = 1.-a5 ;
341 G4double a7 = 1.+a6 ;
342 G4double a9 = 3.+a6 ;
343 G4double xi = xi0*a5 ;
344 G4double xii = 1./xi ;
345 G4double xi1 = 1.+xi ;
346 G4double screen = screen0*xi1/a5 ;
347 G4double yeu = 5.-a6+4.*bet*a7 ;
348 G4double yed = 2.*(1.+3.*bet)*log(3.+xii)-a6-a1*(2.-a6) ;
349 G4double ye1 = 1.+yeu/yed ;
350 G4double ale=log(bbb/z13*sqrt(xi1*ye1)/(1.+screen*ye1)) ;
351 G4double cre = 0.5*log(1.+2.25*z23*xi1*ye1/massratio2) ;
352 G4double be;
353
354 if (xi <= 1.e3) be = ((2.+a6)*(1.+bet)+xi*a9)*log(1.+xii)+(a5-bet)/xi1-a9;
355 else be = (3.-a6+a1*a7)/(2.*xi);
356
357 G4double fe = (ale-cre)*be;
358 if ( fe < 0.) fe = 0. ;
359
360 G4double ymu = 4.+a6 +3.*bet*a7 ;
361 G4double ymd = a7*(1.5+a1)*log(3.+xi)+1.-1.5*a6 ;
362 G4double ym1 = 1.+ymu/ymd ;
363 G4double alm_crm = log(bbb*massratio/(1.5*z23*(1.+screen*ym1)));
364 G4double a10,bm;
365 if ( xi >= 1.e-3)
366 {
367 a10 = (1.+a1)*a5 ;
368 bm = (a7*(1.+1.5*bet)-a10*xii)*log(xi1)+xi*(a5-bet)/xi1+a10;
369 } else {
370 bm = (5.-a6+bet*a9)*(xi/2.);
371 }
372
373 G4double fm = alm_crm*bm;
374 if ( fm < 0.) fm = 0. ;
375
376 sum += wgi[i]*a4*(fe+fm/massratio2);
377 }
378
379 cross = -tmn*sum*factorForCross*z2*residEnergy/(totalEnergy*pairEnergy);
380
381 return cross;
382}
383
384//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
385
386G4double G4MuPairProductionModel::ComputeCrossSectionPerAtom(
387 const G4ParticleDefinition*,
388 G4double kineticEnergy,
389 G4double Z, G4double,
390 G4double cutEnergy,
391 G4double maxEnergy)
392{
393 G4double cross = 0.0;
394 if (kineticEnergy <= lowestKinEnergy) return cross;
395
396 SetCurrentElement(Z);
397 G4double tmax = std::min(maxEnergy, kineticEnergy);
398 G4double cut = std::min(cutEnergy, kineticEnergy);
399 if(cut < minPairEnergy) cut = minPairEnergy;
400 if (cut >= tmax) return cross;
401
402 cross = ComputeMicroscopicCrossSection (kineticEnergy, Z, cut);
403 if(tmax < kineticEnergy) {
404 cross -= ComputeMicroscopicCrossSection(kineticEnergy, Z, tmax);
405 }
406 return cross;
407}
408
409//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
410
411void G4MuPairProductionModel::MakeSamplingTables()
412{
413 for (G4int iz=0; iz<nzdat; iz++)
414 {
415 G4double Z = zdat[iz];
416 SetCurrentElement(Z);
417
418 for (G4int it=0; it<ntdat; it++) {
419
420 G4double kineticEnergy = tdat[it];
421 G4double maxPairEnergy = MaxSecondaryEnergy(particle,kineticEnergy);
422 // G4cout << "Z= " << currentZ << " z13= " << z13
423 //<< " mE= " << maxPairEnergy << G4endl;
424 G4double CrossSection = 0.0 ;
425
426 if(maxPairEnergy > minPairEnergy) {
427
428 G4double y = ymin - 0.5*dy ;
429 G4double yy = ymin - dy ;
430 G4double x = exp(y);
431 G4double fac = exp(dy);
432 G4double dx = exp(yy)*(fac - 1.0);
433
434 G4double c = log(maxPairEnergy/minPairEnergy);
435
436 for (G4int i=0 ; i<nbiny; i++) {
437 y += dy ;
438 if(c > 0.0) {
439 x *= fac;
440 dx*= fac;
441 G4double ep = minPairEnergy*exp(c*x) ;
442 CrossSection +=
443 ep*dx*ComputeDMicroscopicCrossSection(kineticEnergy, Z, ep);
444 }
445 ya[i] = y;
446 proba[iz][it][i] = CrossSection;
447 }
448
449 } else {
450 for (G4int i=0 ; i<nbiny; i++) {
451 proba[iz][it][i] = CrossSection;
452 }
453 }
454
455 ya[nbiny]=ymax;
456 proba[iz][it][nbiny] = CrossSection;
457
458 }
459 }
460 samplingTablesAreFilled = true;
461}
462
463//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
464
465void G4MuPairProductionModel::SampleSecondaries(std::vector<G4DynamicParticle*>* vdp,
466 const G4MaterialCutsCouple* couple,
467 const G4DynamicParticle* aDynamicParticle,
468 G4double tmin,
469 G4double tmax)
470{
471 G4double kineticEnergy = aDynamicParticle->GetKineticEnergy();
472 G4double totalEnergy = kineticEnergy + particleMass ;
473 G4ParticleMomentum ParticleDirection =
474 aDynamicParticle->GetMomentumDirection();
475
476 G4int it;
477 for(it=1; it<ntdat; it++) {if(kineticEnergy <= tdat[it]) break;}
478 if(it == ntdat) it--;
479 G4double dt = log(kineticEnergy/tdat[it-1])/log(tdat[it]/tdat[it-1]);
480
481 // select randomly one element constituing the material
482 const G4Element* anElement = SelectRandomAtom(kineticEnergy, dt, it, couple, tmin);
483 SetCurrentElement(anElement->GetZ());
484
485 // define interval of enegry transfer
486 G4double maxPairEnergy = MaxSecondaryEnergy(particle,kineticEnergy);
487 G4double maxEnergy = std::min(tmax, maxPairEnergy);
488 G4double minEnergy = std::max(tmin, minPairEnergy);
489
490 if(minEnergy >= maxEnergy) return;
491 //G4cout << "emin= " << minEnergy << " emax= " << maxEnergy
492 // << " minPair= " << minPairEnergy << " maxpair= " << maxPairEnergy
493 // << " ymin= " << ymin << " dy= " << dy << G4endl;
494
495 // select bins
496 G4int iymin = 0;
497 G4int iymax = nbiny-1;
498 if( minEnergy > minPairEnergy)
499 {
500 G4double xc = log(minEnergy/minPairEnergy)/log(maxPairEnergy/minPairEnergy);
501 iymin = (G4int)((log(xc) - ymin)/dy);
502 if(iymin >= nbiny) iymin = nbiny-1;
503 else if(iymin < 0) iymin = 0;
504 xc = log(maxEnergy/minPairEnergy)/log(maxPairEnergy/minPairEnergy);
505 iymax = (G4int)((log(xc) - ymin)/dy) + 1;
506 if(iymax >= nbiny) iymax = nbiny-1;
507 else if(iymax < 0) iymax = 0;
508 }
509
510 // sample e-e+ energy, pair energy first
511 G4int iz, iy;
512
513 for(iz=1; iz<nzdat; iz++) {if(currentZ <= zdat[iz]) break;}
514 if(iz == nzdat) iz--;
515
516 G4double dz = log(currentZ/zdat[iz-1])/log(zdat[iz]/zdat[iz-1]);
517
518 G4double pmin = InterpolatedIntegralCrossSection(dt,dz,iz,it,iymin,currentZ);
519 G4double pmax = InterpolatedIntegralCrossSection(dt,dz,iz,it,iymax,currentZ);
520
521 G4double p = pmin+G4UniformRand()*(pmax - pmin);
522
523 // interpolate sampling vector;
524 G4double p1 = pmin;
525 G4double p2 = pmin;
526 for(iy=iymin+1; iy<=iymax; iy++) {
527 p1 = p2;
528 p2 = InterpolatedIntegralCrossSection(dt, dz, iz, it, iy, currentZ);
529 if(p <= p2) break;
530 }
531 // G4cout << "iy= " << iy << " iymin= " << iymin << " iymax= "
532 // << iymax << " Z= " << currentZ << G4endl;
533 G4double y = ya[iy-1] + dy*(p - p1)/(p2 - p1);
534
535 G4double PairEnergy = minPairEnergy*exp(exp(y)
536 *log(maxPairEnergy/minPairEnergy));
537
538 if(PairEnergy < minEnergy) PairEnergy = minEnergy;
539 if(PairEnergy > maxEnergy) PairEnergy = maxEnergy;
540
541 // sample r=(E+-E-)/PairEnergy ( uniformly .....)
542 G4double rmax =
543 (1.-6.*particleMass*particleMass/(totalEnergy*(totalEnergy-PairEnergy)))
544 *sqrt(1.-minPairEnergy/PairEnergy);
545 G4double r = rmax * (-1.+2.*G4UniformRand()) ;
546
547 // compute energies from PairEnergy,r
548 G4double ElectronEnergy = (1.-r)*PairEnergy*0.5;
549 G4double PositronEnergy = PairEnergy - ElectronEnergy;
550
551 // angles of the emitted particles ( Z - axis along the parent particle)
552 // (mean theta for the moment)
553
554 //
555 // scattered electron (positron) angles. ( Z - axis along the parent photon)
556 //
557 // universal distribution suggested by L. Urban
558 // (Geant3 manual (1993) Phys211),
559 // derived from Tsai distribution (Rev Mod Phys 49,421(1977))
560 // G4cout << "Ee= " << ElectronEnergy << " Ep= " << PositronEnergy << G4endl;
561 G4double u;
562 const G4double a1 = 0.625 , a2 = 3.*a1 , d = 27. ;
563
564 if (9./(9.+d) >G4UniformRand()) u= - log(G4UniformRand()*G4UniformRand())/a1;
565 else u= - log(G4UniformRand()*G4UniformRand())/a2;
566
567 G4double TetEl = u*electron_mass_c2/ElectronEnergy;
568 G4double TetPo = u*electron_mass_c2/PositronEnergy;
569 G4double Phi = twopi * G4UniformRand();
570 G4double dxEl= sin(TetEl)*cos(Phi),dyEl= sin(TetEl)*sin(Phi),dzEl=cos(TetEl);
571 G4double dxPo=-sin(TetPo)*cos(Phi),dyPo=-sin(TetPo)*sin(Phi),dzPo=cos(TetPo);
572
573 G4ThreeVector ElectDirection (dxEl, dyEl, dzEl);
574 ElectDirection.rotateUz(ParticleDirection);
575
576 // create G4DynamicParticle object for the particle1
577 G4DynamicParticle* aParticle1= new G4DynamicParticle(theElectron,
578 ElectDirection,
579 ElectronEnergy - electron_mass_c2);
580
581 G4ThreeVector PositDirection (dxPo, dyPo, dzPo);
582 PositDirection.rotateUz(ParticleDirection);
583
584 // create G4DynamicParticle object for the particle2
585 G4DynamicParticle* aParticle2 =
586 new G4DynamicParticle(thePositron,
587 PositDirection,
588 PositronEnergy - electron_mass_c2);
589
590 // primary change
591 kineticEnergy -= (ElectronEnergy + PositronEnergy);
592 fParticleChange->SetProposedKineticEnergy(kineticEnergy);
593
594 vdp->push_back(aParticle1);
595 vdp->push_back(aParticle2);
596}
597
598//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
599
600const G4Element* G4MuPairProductionModel::SelectRandomAtom(
601 G4double kinEnergy, G4double dt, G4int it,
602 const G4MaterialCutsCouple* couple, G4double tmin)
603{
604 // select randomly 1 element within the material
605
606 const G4Material* material = couple->GetMaterial();
607 size_t nElements = material->GetNumberOfElements();
608 const G4ElementVector* theElementVector = material->GetElementVector();
609 if (nElements == 1) return (*theElementVector)[0];
610
611 if(nElements > nmaxElements) {
612 nmaxElements = nElements;
613 partialSum.resize(nmaxElements);
614 }
615
616 const G4double* theAtomNumDensityVector=material->GetAtomicNumDensityVector();
617
618 G4double sum = 0.0;
619
620 size_t i;
621 for (i=0; i<nElements; i++) {
622 G4double Z = ((*theElementVector)[i])->GetZ();
623 SetCurrentElement(Z);
624 G4double maxPairEnergy = MaxSecondaryEnergy(particle,kinEnergy);
625 G4double minEnergy = std::max(tmin, minPairEnergy);
626
627 G4int iz;
628 for(iz=1; iz<nzdat; iz++) {if(Z <= zdat[iz]) break;}
629 if(iz == nzdat) iz--;
630 G4double dz = log(Z/zdat[iz-1])/log(zdat[iz]/zdat[iz-1]);
631
632 G4double sigcut;
633 if(minEnergy <= minPairEnergy)
634 sigcut = 0.;
635 else
636 {
637 G4double xc = log(minEnergy/minPairEnergy)/log(maxPairEnergy/minPairEnergy);
638 G4int iy = (G4int)((log(xc) - ymin)/dy);
639 if(iy < 0) iy = 0;
640 if(iy >= nbiny) iy = nbiny-1;
641 sigcut = InterpolatedIntegralCrossSection(dt,dz,iz,it,iy, Z);
642 }
643
644 G4double sigtot = InterpolatedIntegralCrossSection(dt,dz,iz,it,nbiny,Z);
645 G4double dl = (sigtot - sigcut)*theAtomNumDensityVector[i];
646
647 // protection
648 if(dl < 0.0) dl = 0.0;
649 sum += dl;
650 partialSum[i] = sum;
651 }
652
653 G4double rval = G4UniformRand()*sum;
654 for (i=0; i<nElements; i++) {
655 if(rval<=partialSum[i]) return (*theElementVector)[i];
656 }
657
658 return (*theElementVector)[nElements - 1];
659
660}
661
662//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
663
664
Note: See TracBrowser for help on using the repository browser.