Geant4  9.6.p02
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4MuPairProductionModel.cc
Go to the documentation of this file.
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$
27 //
28 // -------------------------------------------------------------------
29 //
30 // GEANT4 Class file
31 //
32 //
33 // File name: G4MuPairProductionModel
34 //
35 // Author: Vladimir Ivanchenko on base of Laszlo Urban code
36 //
37 // Creation date: 24.06.2002
38 //
39 // Modifications:
40 //
41 // 04-12-02 Change G4DynamicParticle constructor in PostStep (V.Ivanchenko)
42 // 23-12-02 Change interface in order to move to cut per region (V.Ivanchenko)
43 // 24-01-03 Fix for compounds (V.Ivanchenko)
44 // 27-01-03 Make models region aware (V.Ivanchenko)
45 // 13-02-03 Add model (V.Ivanchenko)
46 // 06-06-03 Fix in cross section calculation for high energy (V.Ivanchenko)
47 // 20-10-03 2*xi in ComputeDDMicroscopicCrossSection (R.Kokoulin)
48 // 8 integration points in ComputeDMicroscopicCrossSection
49 // 12-01-04 Take min cut of e- and e+ not its sum (V.Ivanchenko)
50 // 10-02-04 Update parameterisation using R.Kokoulin model (V.Ivanchenko)
51 // 28-04-04 For complex materials repeat calculation of max energy for each
52 // material (V.Ivanchenko)
53 // 01-11-04 Fix bug inside ComputeDMicroscopicCrossSection (R.Kokoulin)
54 // 08-04-05 Major optimisation of internal interfaces (V.Ivantchenko)
55 // 03-08-05 Add SetParticle method (V.Ivantchenko)
56 // 23-10-05 Add protection in sampling of e+e- pair energy needed for
57 // low cuts (V.Ivantchenko)
58 // 13-02-06 Add ComputeCrossSectionPerAtom (mma)
59 // 24-04-07 Add protection in SelectRandomAtom method (V.Ivantchenko)
60 // 12-05-06 Updated sampling (use cut) in SelectRandomAtom (A.Bogdanov)
61 // 11-10-07 Add ignoreCut flag (V.Ivanchenko)
62 
63 //
64 // Class Description:
65 //
66 //
67 // -------------------------------------------------------------------
68 //
69 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
70 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
71 
72 #include "G4MuPairProductionModel.hh"
73 #include "G4PhysicalConstants.hh"
74 #include "G4SystemOfUnits.hh"
75 #include "G4Electron.hh"
76 #include "G4Positron.hh"
77 #include "G4MuonMinus.hh"
78 #include "G4MuonPlus.hh"
79 #include "Randomize.hh"
80 #include "G4Material.hh"
81 #include "G4Element.hh"
82 #include "G4ElementVector.hh"
83 #include "G4ProductionCutsTable.hh"
84 #include "G4ParticleChangeForLoss.hh"
85 #include "G4ParticleChangeForGamma.hh"
86 
87 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
88 
89 // static members
90 //
91 G4double G4MuPairProductionModel::zdat[]={1., 4., 13., 29., 92.};
92 G4double G4MuPairProductionModel::adat[]={1.01, 9.01, 26.98, 63.55, 238.03};
93 G4double G4MuPairProductionModel::tdat[]={1.e3, 1.e4, 1.e5, 1.e6, 1.e7, 1.e8,
94  1.e9, 1.e10};
95 G4double G4MuPairProductionModel::xgi[]={ 0.0199, 0.1017, 0.2372, 0.4083,
96  0.5917, 0.7628, 0.8983, 0.9801 };
97 G4double G4MuPairProductionModel::wgi[]={ 0.0506, 0.1112, 0.1569, 0.1813,
98  0.1813, 0.1569, 0.1112, 0.0506 };
99 
100 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
101 
102 using namespace std;
103 
104 G4MuPairProductionModel::G4MuPairProductionModel(const G4ParticleDefinition* p,
105  const G4String& nam)
106  : G4VEmModel(nam),
107  particle(0),
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  fParticleChange(0),
113  minPairEnergy(4.*electron_mass_c2),
114  lowestKinEnergy(GeV),
115  nzdat(5),
116  ntdat(8),
117  nbiny(1000),
118  nmaxElements(0),
119  ymin(-5.),
120  ymax(0.),
121  dy((ymax-ymin)/nbiny),
122  samplingTablesAreFilled(false)
123 {
124  SetLowEnergyLimit(minPairEnergy);
125  nist = G4NistManager::Instance();
126 
127  theElectron = G4Electron::Electron();
128  thePositron = G4Positron::Positron();
129 
130  particleMass = lnZ = z13 = z23 = 0;
131 
132  for(size_t i=0; i<1001; ++i) { ya[i] = 0.0; }
133 
134  if(p) { SetParticle(p); }
135 }
136 
137 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
138 
139 G4MuPairProductionModel::~G4MuPairProductionModel()
140 {}
141 
142 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
143 
144 G4double G4MuPairProductionModel::MinEnergyCut(const G4ParticleDefinition*,
145  const G4MaterialCutsCouple* )
146 {
147  return minPairEnergy;
148 }
149 
150 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
151 
152 G4double G4MuPairProductionModel::MaxSecondaryEnergy(const G4ParticleDefinition*,
153  G4double kineticEnergy)
154 {
155  G4double maxPairEnergy = kineticEnergy + particleMass*(1.0 - 0.75*sqrte*z13);
156  return maxPairEnergy;
157 }
158 
159 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
160 
161 void G4MuPairProductionModel::Initialise(const G4ParticleDefinition* p,
162  const G4DataVector&)
163 {
164  if (!samplingTablesAreFilled) {
165  if(p) { SetParticle(p); }
166  MakeSamplingTables();
167  }
168  if(!fParticleChange) { fParticleChange = GetParticleChangeForLoss(); }
169 }
170 
171 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
172 
173 G4double G4MuPairProductionModel::ComputeDEDXPerVolume(
174  const G4Material* material,
175  const G4ParticleDefinition*,
176  G4double kineticEnergy,
177  G4double cutEnergy)
178 {
179  G4double dedx = 0.0;
180  if (cutEnergy <= minPairEnergy || kineticEnergy <= lowestKinEnergy)
181  { return dedx; }
182 
183  const G4ElementVector* theElementVector = material->GetElementVector();
184  const G4double* theAtomicNumDensityVector =
185  material->GetAtomicNumDensityVector();
186 
187  // loop for elements in the material
188  for (size_t i=0; i<material->GetNumberOfElements(); ++i) {
189  G4double Z = (*theElementVector)[i]->GetZ();
190  SetCurrentElement(Z);
191  G4double tmax = MaxSecondaryEnergy(particle, kineticEnergy);
192  G4double loss = ComputMuPairLoss(Z, kineticEnergy, cutEnergy, tmax);
193  dedx += loss*theAtomicNumDensityVector[i];
194  }
195  if (dedx < 0.) { dedx = 0.; }
196  return dedx;
197 }
198 
199 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
200 
201 G4double G4MuPairProductionModel::ComputMuPairLoss(G4double Z,
202  G4double tkin,
203  G4double cutEnergy,
204  G4double tmax)
205 {
206  SetCurrentElement(Z);
207  G4double loss = 0.0;
208 
209  G4double cut = std::min(cutEnergy,tmax);
210  if(cut <= minPairEnergy) { return loss; }
211 
212  // calculate the rectricted loss
213  // numerical integration in log(PairEnergy)
214  G4double ak1=6.9;
215  G4double ak2=1.0;
216  G4double aaa = log(minPairEnergy);
217  G4double bbb = log(cut);
218  G4int kkk = (G4int)((bbb-aaa)/ak1+ak2);
219  if (kkk > 8) kkk = 8;
220  else if (kkk < 1) { kkk = 1; }
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 
241 G4double 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)/G4double(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 
275 G4double 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 
386 G4double 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 
398  G4double maxPairEnergy = MaxSecondaryEnergy(particle,kineticEnergy);
399  G4double tmax = std::min(maxEnergy, maxPairEnergy);
400  G4double cut = std::max(cutEnergy, minPairEnergy);
401  if (cut >= tmax) return cross;
402 
403  cross = ComputeMicroscopicCrossSection (kineticEnergy, Z, cut);
404  if(tmax < kineticEnergy) {
405  cross -= ComputeMicroscopicCrossSection(kineticEnergy, Z, tmax);
406  }
407  return cross;
408 }
409 
410 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
411 
412 void G4MuPairProductionModel::MakeSamplingTables()
413 {
414  for (G4int iz=0; iz<nzdat; ++iz)
415  {
416  G4double Z = zdat[iz];
417  SetCurrentElement(Z);
418 
419  for (G4int it=0; it<ntdat; ++it) {
420 
421  G4double kineticEnergy = tdat[it];
422  G4double maxPairEnergy = MaxSecondaryEnergy(particle,kineticEnergy);
423  // G4cout << "Z= " << currentZ << " z13= " << z13
424  //<< " mE= " << maxPairEnergy << G4endl;
425  G4double xSec = 0.0 ;
426 
427  if(maxPairEnergy > minPairEnergy) {
428 
429  G4double y = ymin - 0.5*dy ;
430  G4double yy = ymin - dy ;
431  G4double x = exp(y);
432  G4double fac = exp(dy);
433  G4double dx = exp(yy)*(fac - 1.0);
434 
435  G4double c = log(maxPairEnergy/minPairEnergy);
436 
437  for (G4int i=0 ; i<nbiny; ++i) {
438  y += dy ;
439  if(c > 0.0) {
440  x *= fac;
441  dx*= fac;
442  G4double ep = minPairEnergy*exp(c*x) ;
443  xSec +=
444  ep*dx*ComputeDMicroscopicCrossSection(kineticEnergy, Z, ep);
445  }
446  ya[i] = y;
447  proba[iz][it][i] = xSec;
448  }
449 
450  } else {
451  for (G4int i=0 ; i<nbiny; ++i) {
452  proba[iz][it][i] = xSec;
453  }
454  }
455 
456  ya[nbiny]=ymax;
457  proba[iz][it][nbiny] = xSec;
458 
459  }
460  }
461  samplingTablesAreFilled = true;
462 }
463 
464 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
465 
466 void
467 G4MuPairProductionModel::SampleSecondaries(std::vector<G4DynamicParticle*>* vdp,
468  const G4MaterialCutsCouple* couple,
469  const G4DynamicParticle* aDynamicParticle,
470  G4double tmin,
471  G4double tmax)
472 {
473  G4double kineticEnergy = aDynamicParticle->GetKineticEnergy();
474  G4double totalEnergy = kineticEnergy + particleMass;
475  G4double totalMomentum =
476  sqrt(kineticEnergy*(kineticEnergy + 2.0*particleMass));
477 
478  G4ThreeVector partDirection = aDynamicParticle->GetMomentumDirection();
479 
480  G4int it;
481  for(it=1; it<ntdat; ++it) { if(kineticEnergy <= tdat[it]) { break; } }
482  if(it == ntdat) { --it; }
483  G4double dt = log(kineticEnergy/tdat[it-1])/log(tdat[it]/tdat[it-1]);
484 
485  // select randomly one element constituing the material
486  const G4Element* anElement =
487  SelectRandomAtom(kineticEnergy, dt, it, couple, tmin);
488  SetCurrentElement(anElement->GetZ());
489 
490  // define interval of enegry transfer
491  G4double maxPairEnergy = MaxSecondaryEnergy(particle,kineticEnergy);
492  G4double maxEnergy = std::min(tmax, maxPairEnergy);
493  G4double minEnergy = std::max(tmin, minPairEnergy);
494 
495  if(minEnergy >= maxEnergy) { return; }
496  //G4cout << "emin= " << minEnergy << " emax= " << maxEnergy
497  // << " minPair= " << minPairEnergy << " maxpair= " << maxPairEnergy
498  // << " ymin= " << ymin << " dy= " << dy << G4endl;
499 
500  G4double logmaxmin = log(maxPairEnergy/minPairEnergy);
501 
502  // select bins
503  G4int iymin = 0;
504  G4int iymax = nbiny-1;
505  if( minEnergy > minPairEnergy)
506  {
507  G4double xc = log(minEnergy/minPairEnergy)/logmaxmin;
508  iymin = (G4int)((log(xc) - ymin)/dy);
509  if(iymin >= nbiny) iymin = nbiny-1;
510  else if(iymin < 0) iymin = 0;
511  xc = log(maxEnergy/minPairEnergy)/logmaxmin;
512  iymax = (G4int)((log(xc) - ymin)/dy) + 1;
513  if(iymax >= nbiny) iymax = nbiny-1;
514  else if(iymax < 0) iymax = 0;
515  }
516 
517  // sample e-e+ energy, pair energy first
518  G4int iz, iy;
519 
520  for(iz=1; iz<nzdat; ++iz) { if(currentZ <= zdat[iz]) { break; } }
521  if(iz == nzdat) { --iz; }
522 
523  G4double dz = log(currentZ/zdat[iz-1])/log(zdat[iz]/zdat[iz-1]);
524 
525  G4double pmin = InterpolatedIntegralCrossSection(dt,dz,iz,it,iymin,currentZ);
526  G4double pmax = InterpolatedIntegralCrossSection(dt,dz,iz,it,iymax,currentZ);
527 
528  G4double p = pmin+G4UniformRand()*(pmax - pmin);
529 
530  // interpolate sampling vector;
531  G4double p1 = pmin;
532  G4double p2 = pmin;
533  for(iy=iymin+1; iy<=iymax; ++iy) {
534  p1 = p2;
535  p2 = InterpolatedIntegralCrossSection(dt, dz, iz, it, iy, currentZ);
536  if(p <= p2) { break; }
537  }
538  // G4cout << "iy= " << iy << " iymin= " << iymin << " iymax= "
539  // << iymax << " Z= " << currentZ << G4endl;
540  G4double y = ya[iy-1] + dy*(p - p1)/(p2 - p1);
541 
542  G4double PairEnergy = minPairEnergy*exp( exp(y)*logmaxmin );
543 
544  if(PairEnergy < minEnergy) { PairEnergy = minEnergy; }
545  if(PairEnergy > maxEnergy) { PairEnergy = maxEnergy; }
546 
547  // sample r=(E+-E-)/PairEnergy ( uniformly .....)
548  G4double rmax =
549  (1.-6.*particleMass*particleMass/(totalEnergy*(totalEnergy-PairEnergy)))
550  *sqrt(1.-minPairEnergy/PairEnergy);
551  G4double r = rmax * (-1.+2.*G4UniformRand()) ;
552 
553  // compute energies from PairEnergy,r
554  G4double ElectronEnergy = (1.-r)*PairEnergy*0.5;
555  G4double PositronEnergy = PairEnergy - ElectronEnergy;
556 
557  // The angle of the emitted virtual photon is sampled
558  // according to the muon bremsstrahlung model
559 
560  G4double gam = totalEnergy/particleMass;
561  G4double gmax = gam*std::min(1.0, totalEnergy/PairEnergy - 1.0);
562  G4double gmax2= gmax*gmax;
563  G4double x = G4UniformRand()*gmax2/(1.0 + gmax2);
564 
565  G4double theta = sqrt(x/(1.0 - x))/gam;
566  G4double sint = sin(theta);
567  G4double phi = twopi * G4UniformRand() ;
568  G4double dirx = sint*cos(phi), diry = sint*sin(phi), dirz = cos(theta) ;
569 
570  G4ThreeVector gDirection(dirx, diry, dirz);
571  gDirection.rotateUz(partDirection);
572 
573  // the angles of e- and e+ assumed to be the same as virtual gamma
574 
575  // create G4DynamicParticle object for the particle1
576  G4DynamicParticle* aParticle1 =
577  new G4DynamicParticle(theElectron, gDirection,
578  ElectronEnergy - electron_mass_c2);
579 
580  // create G4DynamicParticle object for the particle2
581  G4DynamicParticle* aParticle2 =
582  new G4DynamicParticle(thePositron, gDirection,
583  PositronEnergy - electron_mass_c2);
584 
585  // primary change
586  kineticEnergy -= (ElectronEnergy + PositronEnergy);
587  fParticleChange->SetProposedKineticEnergy(kineticEnergy);
588 
589  partDirection *= totalMomentum;
590  partDirection -= (aParticle1->GetMomentum() + aParticle2->GetMomentum());
591  partDirection = partDirection.unit();
592  fParticleChange->SetProposedMomentumDirection(partDirection);
593 
594  // add secondary
595  vdp->push_back(aParticle1);
596  vdp->push_back(aParticle2);
597 }
598 
599 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
600 
601 const G4Element* G4MuPairProductionModel::SelectRandomAtom(
602  G4double kinEnergy, G4double dt, G4int it,
603  const G4MaterialCutsCouple* couple, G4double tmin)
604 {
605  // select randomly 1 element within the material
606 
607  const G4Material* material = couple->GetMaterial();
608  size_t nElements = material->GetNumberOfElements();
609  const G4ElementVector* theElementVector = material->GetElementVector();
610  if (nElements == 1) { return (*theElementVector)[0]; }
611 
612  if(nElements > nmaxElements) {
613  nmaxElements = nElements;
614  partialSum.resize(nmaxElements);
615  }
616 
617  const G4double* theAtomNumDensityVector=material->GetAtomicNumDensityVector();
618 
619  G4double sum = 0.0;
620  G4double dl;
621 
622  size_t i;
623  for (i=0; i<nElements; ++i) {
624  G4double Z = ((*theElementVector)[i])->GetZ();
625  SetCurrentElement(Z);
626  G4double maxPairEnergy = MaxSecondaryEnergy(particle,kinEnergy);
627  G4double minEnergy = std::max(tmin, minPairEnergy);
628  dl = 0.0;
629  if(minEnergy < maxPairEnergy) {
630 
631  G4int iz;
632  for(iz=1; iz<nzdat; ++iz) {if(Z <= zdat[iz]) { break; } }
633  if(iz == nzdat) { --iz; }
634  G4double dz = log(Z/zdat[iz-1])/log(zdat[iz]/zdat[iz-1]);
635 
636  G4double sigcut;
637  if(minEnergy <= minPairEnergy)
638  sigcut = 0.;
639  else
640  {
641  G4double xc = log(minEnergy/minPairEnergy)/log(maxPairEnergy/minPairEnergy);
642  G4int iy = (G4int)((log(xc) - ymin)/dy);
643  if(iy < 0) { iy = 0; }
644  if(iy >= nbiny) { iy = nbiny-1; }
645  sigcut = InterpolatedIntegralCrossSection(dt,dz,iz,it,iy, Z);
646  }
647 
648  G4double sigtot = InterpolatedIntegralCrossSection(dt,dz,iz,it,nbiny,Z);
649  dl = (sigtot - sigcut)*theAtomNumDensityVector[i];
650  }
651  // protection
652  if(dl < 0.0) { dl = 0.0; }
653  sum += dl;
654  partialSum[i] = sum;
655  }
656 
657  G4double rval = G4UniformRand()*sum;
658  for (i=0; i<nElements; ++i) {
659  if(rval<=partialSum[i]) { return (*theElementVector)[i]; }
660  }
661 
662  return (*theElementVector)[nElements - 1];
663 
664 }
665 
666 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
667 
668 
Generated on Sat May 25 2013 14:33:35 for Geant4 by   doxygen 1.8.4