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

Last change on this file since 924 was 819, checked in by garnier, 16 years ago

import all except CVS
File size: 23.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: G4MuPairProductionModel.cc,v 1.35 2007/10/11 13:52:04 vnivanch Exp $
27// GEANT4 tag $Name: geant4-09-01-patch-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  minPairEnergy(4.*electron_mass_c2),
107  lowestKinEnergy(1.*GeV),
108  factorForCross(4.*fine_structure_const*fine_structure_const
109                   *classic_electr_radius*classic_electr_radius/(3.*pi)),
110    sqrte(sqrt(exp(1.))),
111    currentZ(0),
112    particle(0),
113    nzdat(5),
114    ntdat(8),
115    nbiny(1000),
116    nmaxElements(0),
117    ymin(-5.),
118    ymax(0.),
119    dy((ymax-ymin)/nbiny),
120    ignoreCut(false),
121    samplingTablesAreFilled(false)
122{
123  SetLowEnergyLimit(minPairEnergy);
124
125  theElectron = G4Electron::Electron();
126  thePositron = G4Positron::Positron();
127
128  if(p) SetParticle(p);
129}
130
131//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
132
133G4MuPairProductionModel::~G4MuPairProductionModel()
134{}
135
136//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
137
138G4double G4MuPairProductionModel::MinEnergyCut(const G4ParticleDefinition*,
139                                               const G4MaterialCutsCouple* )
140{
141  return minPairEnergy;
142}
143
144//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
145
146void G4MuPairProductionModel::SetParticle(const G4ParticleDefinition* p)
147{
148  if(!particle) {
149    particle = p;
150    particleMass = particle->GetPDGMass();
151  }
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(pParticleChange) {
164    if(ignoreCut) {
165      gParticleChange = 
166        reinterpret_cast<G4ParticleChangeForGamma*>(pParticleChange);
167      fParticleChange = 0;
168    } else {
169      fParticleChange = 
170        reinterpret_cast<G4ParticleChangeForLoss*>(pParticleChange);
171      gParticleChange = 0;
172    }
173  } else {
174    fParticleChange = new G4ParticleChangeForLoss();
175    gParticleChange = 0;
176  }
177}
178
179//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
180
181G4double G4MuPairProductionModel::ComputeDEDXPerVolume(
182                                              const G4Material* material,
183                                              const G4ParticleDefinition*,
184                                                    G4double kineticEnergy,
185                                                    G4double cutEnergy)
186{
187  G4double dedx = 0.0;
188  if (cutEnergy <= minPairEnergy || kineticEnergy <= lowestKinEnergy
189      || ignoreCut)
190    return dedx;
191
192  const G4ElementVector* theElementVector = material->GetElementVector();
193  const G4double* theAtomicNumDensityVector =
194                                   material->GetAtomicNumDensityVector();
195
196  //  loop for elements in the material
197  for (size_t i=0; i<material->GetNumberOfElements(); i++) {
198     G4double Z = (*theElementVector)[i]->GetZ();
199     SetCurrentElement(Z);
200     G4double tmax = MaxSecondaryEnergy(particle, kineticEnergy);
201     G4double loss = ComputMuPairLoss(Z, kineticEnergy, cutEnergy, tmax);
202     dedx += loss*theAtomicNumDensityVector[i];
203  }
204  if (dedx < 0.) dedx = 0.;
205  return dedx;
206}
207
208//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
209
210G4double G4MuPairProductionModel::ComputMuPairLoss(G4double Z, 
211                                                   G4double tkin,
212                                                   G4double cutEnergy, 
213                                                   G4double tmax)
214{
215  SetCurrentElement(Z);
216  G4double loss = 0.0;
217
218  G4double cut  = min(cutEnergy,tmax);
219  if(cut <= minPairEnergy) return loss;
220
221  // calculate the rectricted loss
222  // numerical integration in log(PairEnergy)
223  G4double ak1=6.9;
224  G4double ak2=1.0;
225  G4double aaa = log(minPairEnergy);
226  G4double bbb = log(cut);
227  G4int    kkk = (G4int)((bbb-aaa)/ak1+ak2);
228  if (kkk > 8) kkk = 8;
229  G4double hhh = (bbb-aaa)/(G4double)kkk;
230  G4double x = aaa;
231
232  for (G4int l=0 ; l<kkk; l++)
233  {
234
235    for (G4int ll=0; ll<8; ll++)
236    {
237      G4double ep = exp(x+xgi[ll]*hhh);
238      loss += wgi[ll]*ep*ep*ComputeDMicroscopicCrossSection(tkin, Z, ep);
239    }
240    x += hhh;
241  }
242  loss *= hhh;
243  if (loss < 0.) loss = 0.;
244  return loss;
245}
246
247//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
248
249G4double G4MuPairProductionModel::ComputeMicroscopicCrossSection(
250                                           G4double tkin,
251                                           G4double Z,
252                                           G4double cut)
253
254{
255  G4double cross = 0. ;
256
257  SetCurrentElement(Z);
258  G4double tmax = MaxSecondaryEnergy(particle, tkin);
259
260  if (tmax <= cut) return cross;
261
262  G4double ak1=6.9 ;
263  G4double ak2=1.0 ;
264  G4double aaa = log(cut);
265  G4double bbb = log(tmax);
266  G4int kkk = (G4int)((bbb-aaa)/ak1 + ak2);
267  if(kkk > 8) kkk = 8;
268  G4double hhh = (bbb-aaa)/float(kkk);
269  G4double x = aaa;
270
271  for(G4int l=0; l<kkk; l++)
272  {
273    for(G4int i=0; i<8; i++)
274    {
275      G4double ep = exp(x + xgi[i]*hhh);
276      cross += ep*wgi[i]*ComputeDMicroscopicCrossSection(tkin, Z, ep);
277    }
278    x += hhh;
279  }
280
281  cross *=hhh;
282  if(cross < 0.0) cross = 0.0;
283  return cross;
284}
285
286G4double G4MuPairProductionModel::ComputeDMicroscopicCrossSection(
287                                           G4double tkin,
288                                           G4double Z,
289                                           G4double pairEnergy)
290 // Calculates the  differential (D) microscopic cross section
291 // using the cross section formula of R.P. Kokoulin (18/01/98)
292 // Code modified by R.P. Kokoulin, V.N. Ivanchenko (27/01/04)
293{
294  G4double bbbtf= 183. ;
295  G4double bbbh = 202.4 ;
296  G4double g1tf = 1.95e-5 ;
297  G4double g2tf = 5.3e-5 ;
298  G4double g1h  = 4.4e-5 ;
299  G4double g2h  = 4.8e-5 ;
300
301  G4double totalEnergy  = tkin + particleMass;
302  G4double residEnergy  = totalEnergy - pairEnergy;
303  G4double massratio    = particleMass/electron_mass_c2 ;
304  G4double massratio2   = massratio*massratio ;
305  G4double cross = 0.;
306
307  SetCurrentElement(Z);
308
309  G4double c3 = 0.75*sqrte*particleMass;
310  if (residEnergy <= c3*z13) return cross;
311
312  G4double c7 = 4.*electron_mass_c2;
313  G4double c8 = 6.*particleMass*particleMass;
314  G4double alf = c7/pairEnergy;
315  G4double a3 = 1. - alf;
316  if (a3 <= 0.) return cross;
317
318  // zeta calculation
319  G4double bbb,g1,g2;
320  if( Z < 1.5 ) { bbb = bbbh ; g1 = g1h ; g2 = g2h ; }
321  else          { bbb = bbbtf; g1 = g1tf; g2 = g2tf; }
322
323  G4double zeta = 0;
324  G4double zeta1 = 0.073*log(totalEnergy/(particleMass+g1*z23*totalEnergy))-0.26;
325  if ( zeta1 > 0.)
326  {
327    G4double zeta2 = 0.058*log(totalEnergy/(particleMass+g2*z13*totalEnergy))-0.14;
328    zeta  = zeta1/zeta2 ;
329  }
330
331  G4double z2 = Z*(Z+zeta);
332  G4double screen0 = 2.*electron_mass_c2*sqrte*bbb/(z13*pairEnergy);
333  G4double a0 = totalEnergy*residEnergy;
334  G4double a1 = pairEnergy*pairEnergy/a0;
335  G4double bet = 0.5*a1;
336  G4double xi0 = 0.25*massratio2*a1;
337  G4double del = c8/a0;
338
339  G4double rta3 = sqrt(a3);
340  G4double tmnexp = alf/(1. + rta3) + del*rta3;
341  if(tmnexp >= 1.0) return cross;
342
343  G4double tmn = log(tmnexp);
344  G4double sum = 0.;
345
346  // Gaussian integration in ln(1-ro) ( with 8 points)
347  for (G4int i=0; i<8; i++)
348  {
349    G4double a4 = exp(tmn*xgi[i]);     // a4 = (1.-asymmetry)
350    G4double a5 = a4*(2.-a4) ;
351    G4double a6 = 1.-a5 ;
352    G4double a7 = 1.+a6 ;
353    G4double a9 = 3.+a6 ;
354    G4double xi = xi0*a5 ;
355    G4double xii = 1./xi ;
356    G4double xi1 = 1.+xi ;
357    G4double screen = screen0*xi1/a5 ;
358    G4double yeu = 5.-a6+4.*bet*a7 ;
359    G4double yed = 2.*(1.+3.*bet)*log(3.+xii)-a6-a1*(2.-a6) ;
360    G4double ye1 = 1.+yeu/yed ;
361    G4double ale=log(bbb/z13*sqrt(xi1*ye1)/(1.+screen*ye1)) ;
362    G4double cre = 0.5*log(1.+2.25*z23*xi1*ye1/massratio2) ;
363    G4double be;
364
365    if (xi <= 1.e3) be = ((2.+a6)*(1.+bet)+xi*a9)*log(1.+xii)+(a5-bet)/xi1-a9;
366    else            be = (3.-a6+a1*a7)/(2.*xi);
367
368    G4double fe = (ale-cre)*be;
369    if ( fe < 0.) fe = 0. ;
370
371    G4double ymu = 4.+a6 +3.*bet*a7 ;
372    G4double ymd = a7*(1.5+a1)*log(3.+xi)+1.-1.5*a6 ;
373    G4double ym1 = 1.+ymu/ymd ;
374    G4double alm_crm = log(bbb*massratio/(1.5*z23*(1.+screen*ym1)));
375    G4double a10,bm;
376    if ( xi >= 1.e-3)
377    {
378      a10 = (1.+a1)*a5 ;
379      bm  = (a7*(1.+1.5*bet)-a10*xii)*log(xi1)+xi*(a5-bet)/xi1+a10;
380    } else {
381      bm = (5.-a6+bet*a9)*(xi/2.);
382    }
383
384    G4double fm = alm_crm*bm;
385    if ( fm < 0.) fm = 0. ;
386
387    sum += wgi[i]*a4*(fe+fm/massratio2);
388  }
389
390  cross = -tmn*sum*factorForCross*z2*residEnergy/(totalEnergy*pairEnergy);
391
392  return cross;
393}
394
395//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
396
397G4double G4MuPairProductionModel::ComputeCrossSectionPerAtom(
398                                           const G4ParticleDefinition*,
399                                                 G4double kineticEnergy,
400                                                 G4double Z, G4double,
401                                                 G4double cutEnergy,
402                                                 G4double)
403{
404  G4double cut  = max(minPairEnergy,cutEnergy);
405  if(ignoreCut) cut = minPairEnergy;
406  G4double cross = ComputeMicroscopicCrossSection (kineticEnergy, Z, cut);
407  return cross;
408}
409
410//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
411
412G4double G4MuPairProductionModel::CrossSectionPerVolume(
413                                               const G4Material* material,
414                                               const G4ParticleDefinition*,
415                                                     G4double kineticEnergy,
416                                                     G4double cutEnergy,
417                                                     G4double maxEnergy)
418{
419  G4double cross = 0.0;
420  if (kineticEnergy <= lowestKinEnergy) return cross;
421
422  maxEnergy += particleMass;
423
424  const G4ElementVector* theElementVector = material->GetElementVector();
425  const G4double* theAtomNumDensityVector = material->
426                                                    GetAtomicNumDensityVector();
427
428  for (size_t i=0; i<material->GetNumberOfElements(); i++) {
429    G4double Z = (*theElementVector)[i]->GetZ();
430    SetCurrentElement(Z);
431    G4double tmax = min(maxEnergy,MaxSecondaryEnergy(particle, kineticEnergy));
432    G4double cut  = max(minPairEnergy,cutEnergy);
433    if(ignoreCut) cut = minPairEnergy;
434    if(cut < tmax) {
435      G4double cr = ComputeMicroscopicCrossSection(kineticEnergy, Z, cut)
436                  - ComputeMicroscopicCrossSection(kineticEnergy, Z, tmax);
437
438      cross += theAtomNumDensityVector[i] * cr;
439    }
440  }
441  return cross;
442}
443
444//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
445
446void G4MuPairProductionModel::MakeSamplingTables()
447{
448  for (G4int iz=0; iz<nzdat; iz++)
449  {
450    G4double Z = zdat[iz];
451    SetCurrentElement(Z);
452
453    for (G4int it=0; it<ntdat; it++)
454    {
455      G4double kineticEnergy = tdat[it];
456      G4double maxPairEnergy = MaxSecondaryEnergy(particle,kineticEnergy);
457
458      G4double CrossSection = 0.0 ;
459
460      G4double y = ymin - 0.5*dy ;
461      G4double yy = ymin - dy ;
462      G4double x = exp(y);
463      G4double fac = exp(dy);
464      G4double dx = exp(yy)*(fac - 1.0);
465
466      G4double c = log(maxPairEnergy/minPairEnergy);
467
468      for (G4int i=0 ; i<nbiny; i++)
469      {
470        y += dy ;
471        if(c > 0.0) {
472          x *= fac;
473          dx*= fac;
474          G4double ep = minPairEnergy*exp(c*x) ;
475          CrossSection += ep*dx*ComputeDMicroscopicCrossSection(
476                                      kineticEnergy, Z, ep);
477        }
478        ya[i] = y;
479        proba[iz][it][i] = CrossSection;
480      }
481
482      ya[nbiny]=ymax;
483
484      proba[iz][it][nbiny] = CrossSection;
485
486    }
487  }
488  samplingTablesAreFilled = true;
489}
490
491//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
492
493void G4MuPairProductionModel::SampleSecondaries(std::vector<G4DynamicParticle*>* vdp, 
494                                                const G4MaterialCutsCouple* couple,
495                                                const G4DynamicParticle* aDynamicParticle,
496                                                G4double tmin,
497                                                G4double tmax)
498{
499  G4double kineticEnergy = aDynamicParticle->GetKineticEnergy();
500  G4double totalEnergy   = kineticEnergy + particleMass ;
501  G4ParticleMomentum ParticleDirection = 
502    aDynamicParticle->GetMomentumDirection();
503
504  G4int it;
505  for(it=1; it<ntdat; it++) {if(kineticEnergy <= tdat[it]) break;}
506  if(it == ntdat) it--;
507  G4double dt = log(kineticEnergy/tdat[it-1])/log(tdat[it]/tdat[it-1]);
508
509  // select randomly one element constituing the material
510  const G4Element* anElement = SelectRandomAtom(kineticEnergy, dt, it, couple, tmin);
511  SetCurrentElement(anElement->GetZ());
512
513  // define interval of enegry transfer
514  G4double maxPairEnergy = MaxSecondaryEnergy(particle,kineticEnergy);
515  G4double maxEnergy     = std::min(tmax, maxPairEnergy);
516  G4double minEnergy     = std::max(tmin, minPairEnergy);
517  if(ignoreCut)minEnergy = minPairEnergy;
518  if(minEnergy >= maxEnergy) return;
519  //G4cout << "emin= " << minEnergy << " emax= " << maxEnergy
520  //     << " minPair= " << minPairEnergy << " maxpair= " << maxPairEnergy
521  //       << " ymin= " << ymin << " dy= " << dy << G4endl;
522
523  // select bins
524  G4int iymin = 0;
525  G4int iymax = nbiny-1;
526  if( minEnergy > minPairEnergy)
527  {
528    G4double xc = log(minEnergy/minPairEnergy)/log(maxPairEnergy/minPairEnergy);
529    iymin = (G4int)((log(xc) - ymin)/dy);
530    if(iymin >= nbiny) iymin = nbiny-1;
531    else if(iymin < 0) iymin = 0;
532    xc = log(maxEnergy/minPairEnergy)/log(maxPairEnergy/minPairEnergy);
533    iymax = (G4int)((log(xc) - ymin)/dy) + 1;
534    if(iymax >= nbiny) iymax = nbiny-1;
535    else if(iymax < 0) iymax = 0;
536  }
537
538  // sample e-e+ energy, pair energy first
539  G4int iz, iy;
540
541  for(iz=1; iz<nzdat; iz++) {if(currentZ <= zdat[iz]) break;}
542  if(iz == nzdat) iz--;
543
544  G4double dz = log(currentZ/zdat[iz-1])/log(zdat[iz]/zdat[iz-1]);
545
546  G4double pmin = InterpolatedIntegralCrossSection(dt,dz,iz,it,iymin,currentZ);
547  G4double pmax = InterpolatedIntegralCrossSection(dt,dz,iz,it,iymax,currentZ);
548
549  G4double p = pmin+G4UniformRand()*(pmax - pmin);
550
551  // interpolate sampling vector;
552  G4double p1 = pmin;
553  G4double p2 = pmin;
554  for(iy=iymin+1; iy<=iymax; iy++) {
555    p1 = p2;
556    p2 = InterpolatedIntegralCrossSection(dt, dz, iz, it, iy, currentZ);
557    if(p <= p2) break;
558  }
559  // G4cout << "iy= " << iy << " iymin= " << iymin << " iymax= "
560  //        << iymax << " Z= " << currentZ << G4endl;
561  G4double y = ya[iy-1] + dy*(p - p1)/(p2 - p1);
562
563  G4double PairEnergy = minPairEnergy*exp(exp(y)
564                       *log(maxPairEnergy/minPairEnergy));
565                       
566  if(PairEnergy < minEnergy) PairEnergy = minEnergy;
567  if(PairEnergy > maxEnergy) PairEnergy = maxEnergy;
568
569  // sample r=(E+-E-)/PairEnergy  ( uniformly .....)
570  G4double rmax =
571    (1.-6.*particleMass*particleMass/(totalEnergy*(totalEnergy-PairEnergy)))
572                                       *sqrt(1.-minPairEnergy/PairEnergy);
573  G4double r = rmax * (-1.+2.*G4UniformRand()) ;
574
575  // compute energies from PairEnergy,r
576  G4double ElectronEnergy = (1.-r)*PairEnergy*0.5;
577  G4double PositronEnergy = PairEnergy - ElectronEnergy;
578
579  //  angles of the emitted particles ( Z - axis along the parent particle)
580  //      (mean theta for the moment)
581
582  //
583  // scattered electron (positron) angles. ( Z - axis along the parent photon)
584  //
585  //  universal distribution suggested by L. Urban
586  // (Geant3 manual (1993) Phys211),
587  //  derived from Tsai distribution (Rev Mod Phys 49,421(1977))
588  //  G4cout << "Ee= " << ElectronEnergy << " Ep= " << PositronEnergy << G4endl;
589  G4double u;
590  const G4double a1 = 0.625 , a2 = 3.*a1 , d = 27. ;
591
592  if (9./(9.+d) >G4UniformRand()) u= - log(G4UniformRand()*G4UniformRand())/a1;
593  else                            u= - log(G4UniformRand()*G4UniformRand())/a2;
594
595  G4double TetEl = u*electron_mass_c2/ElectronEnergy;
596  G4double TetPo = u*electron_mass_c2/PositronEnergy;
597  G4double Phi  = twopi * G4UniformRand();
598  G4double dxEl= sin(TetEl)*cos(Phi),dyEl= sin(TetEl)*sin(Phi),dzEl=cos(TetEl);
599  G4double dxPo=-sin(TetPo)*cos(Phi),dyPo=-sin(TetPo)*sin(Phi),dzPo=cos(TetPo);
600
601  G4ThreeVector ElectDirection (dxEl, dyEl, dzEl);
602  ElectDirection.rotateUz(ParticleDirection);
603
604  // create G4DynamicParticle object for the particle1
605  G4DynamicParticle* aParticle1= new G4DynamicParticle(theElectron,
606                                                       ElectDirection,
607                                             ElectronEnergy - electron_mass_c2);
608
609  G4ThreeVector PositDirection (dxPo, dyPo, dzPo);
610  PositDirection.rotateUz(ParticleDirection);
611
612  // create G4DynamicParticle object for the particle2
613  G4DynamicParticle* aParticle2 = 
614    new G4DynamicParticle(thePositron,
615                          PositDirection,
616                          PositronEnergy - electron_mass_c2);
617
618  // primary change
619  kineticEnergy -= (ElectronEnergy + PositronEnergy);
620  if(fParticleChange)
621    fParticleChange->SetProposedKineticEnergy(kineticEnergy);
622  else 
623    gParticleChange->SetProposedKineticEnergy(kineticEnergy);
624
625  vdp->push_back(aParticle1);
626  vdp->push_back(aParticle2);
627}
628
629//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
630
631const G4Element* G4MuPairProductionModel::SelectRandomAtom(
632                 G4double kinEnergy, G4double dt, G4int it,
633           const G4MaterialCutsCouple* couple, G4double tmin)
634{
635  // select randomly 1 element within the material
636
637  const G4Material* material = couple->GetMaterial();
638  size_t nElements = material->GetNumberOfElements();
639  const G4ElementVector* theElementVector = material->GetElementVector();
640  if (nElements == 1) return (*theElementVector)[0];
641
642  if(nElements > nmaxElements) {
643    nmaxElements = nElements;
644    partialSum.resize(nmaxElements);
645  }
646
647  const G4double* theAtomNumDensityVector=material->GetAtomicNumDensityVector();
648
649  G4double sum = 0.0;
650
651  size_t i;
652  for (i=0; i<nElements; i++) {
653    G4double Z = ((*theElementVector)[i])->GetZ();
654    SetCurrentElement(Z);
655    G4double maxPairEnergy = MaxSecondaryEnergy(particle,kinEnergy);
656    G4double minEnergy     = std::max(tmin, minPairEnergy);
657    if(ignoreCut)minEnergy = minPairEnergy;
658
659    G4int iz;
660    for(iz=1; iz<nzdat; iz++) {if(Z <= zdat[iz]) break;}
661    if(iz == nzdat) iz--;
662    G4double dz = log(Z/zdat[iz-1])/log(zdat[iz]/zdat[iz-1]);
663
664    G4double sigcut;
665    if(minEnergy <= minPairEnergy)
666      sigcut = 0.;
667    else
668    {
669      G4double xc = log(minEnergy/minPairEnergy)/log(maxPairEnergy/minPairEnergy);
670      G4int iy = (G4int)((log(xc) - ymin)/dy);
671      if(iy < 0) iy = 0;
672      if(iy >= nbiny) iy = nbiny-1;
673      sigcut = InterpolatedIntegralCrossSection(dt,dz,iz,it,iy,   Z);
674    }
675
676    G4double sigtot = InterpolatedIntegralCrossSection(dt,dz,iz,it,nbiny,Z);
677    G4double dl = (sigtot - sigcut)*theAtomNumDensityVector[i];
678
679    // protection
680    if(dl < 0.0) dl = 0.0;
681    sum += dl;
682    partialSum[i] = sum;
683  }
684
685  G4double rval = G4UniformRand()*sum;
686  for (i=0; i<nElements; i++) {
687    if(rval<=partialSum[i]) return (*theElementVector)[i];
688  }
689
690  return (*theElementVector)[nElements - 1];
691
692}
693
694//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
695
696
Note: See TracBrowser for help on using the repository browser.