Geant4  10.03.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: G4MuPairProductionModel.cc 97392 2016-06-02 10:10:32Z gcosmo $
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 
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"
86 #include "G4Log.hh"
87 #include "G4Exp.hh"
88 
89 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
90 
91 // static members
92 //
93 static const G4double ak1 = 6.9;
94 static const G4double ak2 = 1.0;
95 static const G4int zdat[5] = {1, 4, 13, 29, 92};
96 
98  { 0.0199, 0.1017, 0.2372, 0.4083, 0.5917, 0.7628, 0.8983, 0.9801 };
100  { 0.0506, 0.1112, 0.1569, 0.1813, 0.1813, 0.1569, 0.1112, 0.0506 };
101 
102 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
103 
104 using namespace std;
105 
107  const G4String& nam)
108  : G4VEmModel(nam),
109  particle(nullptr),
112  sqrte(sqrt(G4Exp(1.))),
113  currentZ(0),
114  fParticleChange(nullptr),
115  minPairEnergy(4.*electron_mass_c2),
116  lowestKinEnergy(1.0*GeV),
117  nzdat(5),
118  nYBinPerDecade(4),
119  nbiny(1000),
120  nbine(0),
121  ymin(-5.),
122  dy(0.005)
123 {
125 
126  theElectron = G4Electron::Electron();
127  thePositron = G4Positron::Positron();
128 
129  particleMass = lnZ = z13 = z23 = 0.;
130 
131  // setup lowest limit dependent on particle mass
132  if(p) {
133  SetParticle(p);
134  lowestKinEnergy = std::max(lowestKinEnergy,p->GetPDGMass()*8.0);
135  }
136  emin = lowestKinEnergy;
137  emax = 10.*TeV;
138 }
139 
140 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
141 
143 {}
144 
145 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
146 
148  const G4ParticleDefinition*,
149  G4double cut)
150 {
151  return std::max(lowestKinEnergy,cut);
152 }
153 
154 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
155 
157  const G4DataVector& cuts)
158 {
159  SetParticle(p);
160  if(!fParticleChange) { fParticleChange = GetParticleChangeForLoss(); }
161 
162  // for low-energy application this process should not work
163  if(lowestKinEnergy >= HighEnergyLimit()) { return; }
164 
165  // define scale of internal table for each thread only once
166  if(0 == nbine) {
167  emin = std::max(lowestKinEnergy, LowEnergyLimit());
168  emax = std::max(HighEnergyLimit(), emin*2);
169  nbine = size_t(nYBinPerDecade*std::log10(emax/emin));
170  if(nbine < 3) { nbine = 3; }
171 
172  ymin = G4Log(minPairEnergy/emin);
173  dy = -ymin/G4double(nbiny);
174  }
175 
176  if(IsMaster() && p == particle) {
177 
178  if(!fElementData) {
179  fElementData = new G4ElementData();
180  MakeSamplingTables();
181  }
182  InitialiseElementSelectors(p, cuts);
183  }
184 }
185 
186 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
187 
189  G4VEmModel* masterModel)
190 {
191  if(p == particle && lowestKinEnergy < HighEnergyLimit()) {
193  fElementData = masterModel->GetElementData();
194  }
195 }
196 
197 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
198 
200  const G4Material* material,
201  const G4ParticleDefinition*,
202  G4double kineticEnergy,
203  G4double cutEnergy)
204 {
205  G4double dedx = 0.0;
206  if (cutEnergy <= minPairEnergy || kineticEnergy <= lowestKinEnergy)
207  { return dedx; }
208 
209  const G4ElementVector* theElementVector = material->GetElementVector();
210  const G4double* theAtomicNumDensityVector =
211  material->GetAtomicNumDensityVector();
212 
213  // loop for elements in the material
214  for (size_t i=0; i<material->GetNumberOfElements(); ++i) {
215  G4double Z = (*theElementVector)[i]->GetZ();
216  G4double tmax = MaxSecondaryEnergyForElement(kineticEnergy, Z);
217  G4double loss = ComputMuPairLoss(Z, kineticEnergy, cutEnergy, tmax);
218  dedx += loss*theAtomicNumDensityVector[i];
219  }
220  dedx = std::max(dedx, 0.0);
221  return dedx;
222 }
223 
224 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
225 
227  G4double tkin,
228  G4double cutEnergy,
229  G4double tmax)
230 {
231  G4double loss = 0.0;
232 
233  G4double cut = std::min(cutEnergy,tmax);
234  if(cut <= minPairEnergy) { return loss; }
235 
236  // calculate the rectricted loss
237  // numerical integration in log(PairEnergy)
238  G4double aaa = G4Log(minPairEnergy);
239  G4double bbb = G4Log(cut);
240 
241  G4int kkk = (G4int)((bbb-aaa)/ak1+ak2);
242  if (kkk > 8) { kkk = 8; }
243  else if (kkk < 1) { kkk = 1; }
244 
245  G4double hhh = (bbb-aaa)/(G4double)kkk;
246  G4double x = aaa;
247 
248  for (G4int l=0 ; l<kkk; l++)
249  {
250 
251  for (G4int ll=0; ll<8; ll++)
252  {
253  G4double ep = G4Exp(x+xgi[ll]*hhh);
254  loss += wgi[ll]*ep*ep*ComputeDMicroscopicCrossSection(tkin, Z, ep);
255  }
256  x += hhh;
257  }
258  loss *= hhh;
259  loss = std::max(loss, 0.0);
260  return loss;
261 }
262 
263 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
264 
266  G4double tkin,
267  G4double Z,
268  G4double cutEnergy)
269 {
270  G4double cross = 0.;
271  G4double tmax = MaxSecondaryEnergyForElement(tkin, Z);
272  G4double cut = std::max(cutEnergy, minPairEnergy);
273  if (tmax <= cut) { return cross; }
274 
275  // G4double ak1=6.9 ;
276  // G4double ak2=1.0 ;
277  G4double aaa = G4Log(cut);
278  G4double bbb = G4Log(tmax);
279  G4int kkk = (G4int)((bbb-aaa)/ak1 + ak2);
280  if(kkk > 8) { kkk = 8; }
281  else if (kkk < 1) { kkk = 1; }
282 
283  G4double hhh = (bbb-aaa)/G4double(kkk);
284  G4double x = aaa;
285 
286  for(G4int l=0; l<kkk; ++l)
287  {
288  for(G4int i=0; i<8; ++i)
289  {
290  G4double ep = G4Exp(x + xgi[i]*hhh);
291  cross += ep*wgi[i]*ComputeDMicroscopicCrossSection(tkin, Z, ep);
292  }
293  x += hhh;
294  }
295 
296  cross *= hhh;
297  cross = std::max(cross, 0.0);
298  return cross;
299 }
300 
301 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
302 
304  G4double tkin,
305  G4double Z,
306  G4double pairEnergy)
307 // Calculates the differential (D) microscopic cross section
308 // using the cross section formula of R.P. Kokoulin (18/01/98)
309 // Code modified by R.P. Kokoulin, V.N. Ivanchenko (27/01/04)
310 {
311  static const G4double bbbtf= 183. ;
312  static const G4double bbbh = 202.4 ;
313  static const G4double g1tf = 1.95e-5 ;
314  static const G4double g2tf = 5.3e-5 ;
315  static const G4double g1h = 4.4e-5 ;
316  static const G4double g2h = 4.8e-5 ;
317 
318  G4double totalEnergy = tkin + particleMass;
319  G4double residEnergy = totalEnergy - pairEnergy;
320  G4double massratio = particleMass/electron_mass_c2;
321  G4double massratio2 = massratio*massratio;
322  G4double cross = 0.;
323 
324  G4double c3 = 0.75*sqrte*particleMass;
325  if (residEnergy <= c3*z13) { return cross; }
326 
328  G4double c8 = 6.*particleMass*particleMass;
329  G4double alf = c7/pairEnergy;
330  G4double a3 = 1. - alf;
331  if (a3 <= 0.) { return cross; }
332 
333  // zeta calculation
334  G4double bbb,g1,g2;
335  if( Z < 1.5 ) { bbb = bbbh ; g1 = g1h ; g2 = g2h ; }
336  else { bbb = bbbtf; g1 = g1tf; g2 = g2tf; }
337 
338  G4double zeta = 0;
339  G4double zeta1 =
340  0.073*G4Log(totalEnergy/(particleMass+g1*z23*totalEnergy))-0.26;
341  if ( zeta1 > 0.)
342  {
343  G4double zeta2 =
344  0.058*G4Log(totalEnergy/(particleMass+g2*z13*totalEnergy))-0.14;
345  zeta = zeta1/zeta2 ;
346  }
347 
348  G4double z2 = Z*(Z+zeta);
349  G4double screen0 = 2.*electron_mass_c2*sqrte*bbb/(z13*pairEnergy);
350  G4double a0 = totalEnergy*residEnergy;
351  G4double a1 = pairEnergy*pairEnergy/a0;
352  G4double bet = 0.5*a1;
353  G4double xi0 = 0.25*massratio2*a1;
354  G4double del = c8/a0;
355 
356  G4double rta3 = sqrt(a3);
357  G4double tmnexp = alf/(1. + rta3) + del*rta3;
358  if(tmnexp >= 1.0) { return cross; }
359 
360  G4double tmn = G4Log(tmnexp);
361  G4double sum = 0.;
362 
363  // Gaussian integration in ln(1-ro) ( with 8 points)
364  for (G4int i=0; i<8; ++i)
365  {
366  G4double a4 = G4Exp(tmn*xgi[i]); // a4 = (1.-asymmetry)
367  G4double a5 = a4*(2.-a4) ;
368  G4double a6 = 1.-a5 ;
369  G4double a7 = 1.+a6 ;
370  G4double a9 = 3.+a6 ;
371  G4double xi = xi0*a5 ;
372  G4double xii = 1./xi ;
373  G4double xi1 = 1.+xi ;
374  G4double screen = screen0*xi1/a5 ;
375  G4double yeu = 5.-a6+4.*bet*a7 ;
376  G4double yed = 2.*(1.+3.*bet)*G4Log(3.+xii)-a6-a1*(2.-a6) ;
377  G4double ye1 = 1.+yeu/yed ;
378  G4double ale = G4Log(bbb/z13*sqrt(xi1*ye1)/(1.+screen*ye1)) ;
379  G4double cre = 0.5*G4Log(1.+2.25*z23*xi1*ye1/massratio2) ;
380  G4double be;
381 
382  if (xi <= 1.e3) {
383  be = ((2.+a6)*(1.+bet)+xi*a9)*G4Log(1.+xii)+(a5-bet)/xi1-a9;
384  } else {
385  be = (3.-a6+a1*a7)/(2.*xi);
386  }
387  G4double fe = (ale-cre)*be;
388  if ( fe < 0.) fe = 0. ;
389 
390  G4double ymu = 4.+a6 +3.*bet*a7 ;
391  G4double ymd = a7*(1.5+a1)*G4Log(3.+xi)+1.-1.5*a6 ;
392  G4double ym1 = 1.+ymu/ymd ;
393  G4double alm_crm = G4Log(bbb*massratio/(1.5*z23*(1.+screen*ym1)));
394  G4double a10,bm;
395  if ( xi >= 1.e-3)
396  {
397  a10 = (1.+a1)*a5 ;
398  bm = (a7*(1.+1.5*bet)-a10*xii)*G4Log(xi1)+xi*(a5-bet)/xi1+a10;
399  } else {
400  bm = (5.-a6+bet*a9)*(xi/2.);
401  }
402 
403  G4double fm = alm_crm*bm;
404  if ( fm < 0.) { fm = 0.; }
405 
406  sum += wgi[i]*a4*(fe+fm/massratio2);
407  }
408 
409  cross = -tmn*sum*factorForCross*z2*residEnergy/(totalEnergy*pairEnergy);
410  cross = std::max(cross, 0.0);
411  return cross;
412 }
413 
414 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
415 
417  const G4ParticleDefinition*,
418  G4double kineticEnergy,
420  G4double cutEnergy,
421  G4double maxEnergy)
422 {
423  G4double cross = 0.0;
424  if (kineticEnergy <= lowestKinEnergy) { return cross; }
425 
426  G4double maxPairEnergy = MaxSecondaryEnergyForElement(kineticEnergy, Z);
427  G4double tmax = std::min(maxEnergy, maxPairEnergy);
428  G4double cut = std::max(cutEnergy, minPairEnergy);
429  if (cut >= tmax) { return cross; }
430 
431  cross = ComputeMicroscopicCrossSection(kineticEnergy, Z, cut);
432  if(tmax < kineticEnergy) {
433  cross -= ComputeMicroscopicCrossSection(kineticEnergy, Z, tmax);
434  }
435  return cross;
436 }
437 
438 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
439 
440 void G4MuPairProductionModel::MakeSamplingTables()
441 {
442  G4double factore = G4Exp(G4Log(emax/emin)/G4double(nbine));
443 
444  for (G4int iz=0; iz<nzdat; ++iz) {
445 
446  G4double Z = zdat[iz];
447  G4Physics2DVector* pv = new G4Physics2DVector(nbiny+1,nbine+1);
448  G4double kinEnergy = emin;
449 
450  for (size_t it=0; it<=nbine; ++it) {
451 
452  pv->PutY(it, G4Log(kinEnergy/MeV));
453  G4double maxPairEnergy = MaxSecondaryEnergyForElement(kinEnergy, Z);
454  /*
455  G4cout << "it= " << it << " E= " << kinEnergy
456  << " " << particle->GetParticleName()
457  << " maxE= " << maxPairEnergy << " minE= " << minPairEnergy
458  << " ymin= " << ymin << G4endl;
459  */
460  G4double coef = G4Log(minPairEnergy/kinEnergy)/ymin;
461  G4double ymax = G4Log(maxPairEnergy/kinEnergy)/coef;
462  G4double fac = (ymax - ymin)/dy;
463  size_t imax = (size_t)fac;
464  fac -= (G4double)imax;
465 
466  G4double xSec = 0.0;
467  G4double x = ymin;
468  /*
469  G4cout << "Z= " << currentZ << " z13= " << z13
470  << " mE= " << maxPairEnergy << " ymin= " << ymin
471  << " dy= " << dy << " c= " << coef << G4endl;
472  */
473  // start from zero
474  pv->PutValue(0, it, 0.0);
475  if(0 == it) { pv->PutX(nbiny, 0.0); }
476 
477  for (size_t i=0; i<nbiny; ++i) {
478 
479  if(0 == it) { pv->PutX(i, x); }
480 
481  if(i < imax) {
482  G4double ep = kinEnergy*G4Exp(coef*(x + dy*0.5));
483 
484  // not multiplied by interval, because table
485  // will be used only for sampling
486  //G4cout << "i= " << i << " x= " << x << "E= " << kinEnergy
487  // << " Egamma= " << ep << G4endl;
488  xSec += ep*ComputeDMicroscopicCrossSection(kinEnergy, Z, ep);
489 
490  // last bin before the kinematic limit
491  } else if(i == imax) {
492  G4double ep = kinEnergy*G4Exp(coef*(x + fac*dy*0.5));
493  xSec += ep*fac*ComputeDMicroscopicCrossSection(kinEnergy, Z, ep);
494  }
495  pv->PutValue(i + 1, it, xSec);
496  x += dy;
497  }
498  kinEnergy *= factore;
499 
500  // to avoid precision lost
501  if(it+1 == nbine) { kinEnergy = emax; }
502  }
504  }
505 }
506 
507 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
508 
510  std::vector<G4DynamicParticle*>* vdp,
511  const G4MaterialCutsCouple* couple,
512  const G4DynamicParticle* aDynamicParticle,
513  G4double tmin,
514  G4double tmax)
515 {
516  G4double kineticEnergy = aDynamicParticle->GetKineticEnergy();
517  //G4cout << "------- G4MuPairProductionModel::SampleSecondaries E(MeV)= "
518  // << kineticEnergy << " "
519  // << aDynamicParticle->GetDefinition()->GetParticleName() << G4endl;
520  G4double totalEnergy = kineticEnergy + particleMass;
521  G4double totalMomentum =
522  sqrt(kineticEnergy*(kineticEnergy + 2.0*particleMass));
523 
524  G4ThreeVector partDirection = aDynamicParticle->GetMomentumDirection();
525 
526  // select randomly one element constituing the material
527  const G4Element* anElement = SelectRandomAtom(couple,particle,kineticEnergy);
528 
529  // define interval of energy transfer
530  G4double maxPairEnergy = MaxSecondaryEnergyForElement(kineticEnergy,
531  anElement->GetZ());
532  G4double maxEnergy = std::min(tmax, maxPairEnergy);
533  G4double minEnergy = std::max(tmin, minPairEnergy);
534 
535  if(minEnergy >= maxEnergy) { return; }
536  //G4cout << "emin= " << minEnergy << " emax= " << maxEnergy
537  // << " minPair= " << minPairEnergy << " maxpair= " << maxPairEnergy
538  // << " ymin= " << ymin << " dy= " << dy << G4endl;
539 
540  G4double coeff = G4Log(minPairEnergy/kineticEnergy)/ymin;
541 
542  // compute limits
543  G4double yymin = G4Log(minEnergy/kineticEnergy)/coeff;
544  G4double yymax = G4Log(maxEnergy/kineticEnergy)/coeff;
545 
546  //G4cout << "yymin= " << yymin << " yymax= " << yymax << G4endl;
547 
548  // units should not be used, bacause table was built without
549  G4double logTkin = G4Log(kineticEnergy/MeV);
550 
551  // sample e-e+ energy, pair energy first
552 
553  // select sample table via Z
554  G4int iz1(0), iz2(0);
555  for(G4int iz=0; iz<nzdat; ++iz) {
556  if(currentZ == zdat[iz]) {
557  iz1 = iz2 = currentZ;
558  break;
559  } else if(currentZ < zdat[iz]) {
560  iz2 = zdat[iz];
561  if(iz > 0) { iz1 = zdat[iz-1]; }
562  else { iz1 = iz2; }
563  break;
564  }
565  }
566  if(0 == iz1) { iz1 = iz2 = zdat[nzdat-1]; }
567 
568  G4double PairEnergy = 0.0;
569  G4int count = 0;
570  //G4cout << "start loop Z1= " << iz1 << " Z2= " << iz2 << G4endl;
571  do {
572  ++count;
573  // sampling using only one random number
574  G4double rand = G4UniformRand();
575 
576  G4double x = FindScaledEnergy(iz1, rand, logTkin, yymin, yymax);
577  if(iz1 != iz2) {
578  G4double x2 = FindScaledEnergy(iz2, rand, logTkin, yymin, yymax);
579  G4double lz1= nist->GetLOGZ(iz1);
580  G4double lz2= nist->GetLOGZ(iz2);
581  //G4cout << count << ". x= " << x << " x2= " << x2
582  // << " Z1= " << iz1 << " Z2= " << iz2 << G4endl;
583  x += (x2 - x)*(lnZ - lz1)/(lz2 - lz1);
584  }
585  //G4cout << "x= " << x << " coeff= " << coeff << G4endl;
586  PairEnergy = kineticEnergy*G4Exp(x*coeff);
587 
588  // Loop checking, 03-Aug-2015, Vladimir Ivanchenko
589  } while((PairEnergy < minEnergy || PairEnergy > maxEnergy) && 10 > count);
590 
591  //G4cout << "## PairEnergy(GeV)= " << PairEnergy/GeV
592  // << " Etot(GeV)= " << totalEnergy/GeV << G4endl;
593 
594  // sample r=(E+-E-)/PairEnergy ( uniformly .....)
595  G4double rmax =
596  (1.-6.*particleMass*particleMass/(totalEnergy*(totalEnergy-PairEnergy)))
597  *sqrt(1.-minPairEnergy/PairEnergy);
598  G4double r = rmax * (-1.+2.*G4UniformRand()) ;
599 
600  // compute energies from PairEnergy,r
601  G4double ElectronEnergy = (1.-r)*PairEnergy*0.5;
602  G4double PositronEnergy = PairEnergy - ElectronEnergy;
603 
604  // The angle of the emitted virtual photon is sampled
605  // according to the muon bremsstrahlung model
606 
607  G4double gam = totalEnergy/particleMass;
608  G4double gmax = gam*std::min(1.0, totalEnergy/PairEnergy - 1.0);
609  G4double gmax2= gmax*gmax;
610  G4double x = G4UniformRand()*gmax2/(1.0 + gmax2);
611 
612  G4double theta = sqrt(x/(1.0 - x))/gam;
613  G4double sint = sin(theta);
614  G4double phi = twopi * G4UniformRand() ;
615  G4double dirx = sint*cos(phi), diry = sint*sin(phi), dirz = cos(theta) ;
616 
617  G4ThreeVector gDirection(dirx, diry, dirz);
618  gDirection.rotateUz(partDirection);
619 
620  // the angles of e- and e+ assumed to be the same as virtual gamma
621 
622  // create G4DynamicParticle object for the particle1
623  G4DynamicParticle* aParticle1 =
624  new G4DynamicParticle(theElectron, gDirection,
625  ElectronEnergy - electron_mass_c2);
626 
627  // create G4DynamicParticle object for the particle2
628  G4DynamicParticle* aParticle2 =
629  new G4DynamicParticle(thePositron, gDirection,
630  PositronEnergy - electron_mass_c2);
631 
632  // primary change
633  kineticEnergy -= (ElectronEnergy + PositronEnergy);
634  fParticleChange->SetProposedKineticEnergy(kineticEnergy);
635 
636  partDirection *= totalMomentum;
637  partDirection -= (aParticle1->GetMomentum() + aParticle2->GetMomentum());
638  partDirection = partDirection.unit();
639  fParticleChange->SetProposedMomentumDirection(partDirection);
640 
641  // add secondary
642  vdp->push_back(aParticle1);
643  vdp->push_back(aParticle2);
644  //G4cout << "-- G4MuPairProductionModel::SampleSecondaries done" << G4endl;
645 }
646 
647 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
648 
649 void G4MuPairProductionModel::DataCorrupted(G4int Z, G4double logTkin)
650 {
652  ed << "G4ElementData is not properly initialized Z= " << Z
653  << " Ekin(MeV)= " << G4Exp(logTkin)
654  << " IsMasterThread= " << IsMaster()
655  << " Model " << GetName();
656  G4Exception("G4MuPairProductionModel::()","em0033",FatalException,
657  ed,"");
658 }
659 
660 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
661 
662 
const G4double a0
G4double MaxSecondaryEnergyForElement(G4double kineticEnergy, G4double Z)
G4double LowEnergyLimit() const
Definition: G4VEmModel.hh:643
virtual void Initialise(const G4ParticleDefinition *, const G4DataVector &) override
static const G4double xgi[8]
G4ParticleChangeForLoss * GetParticleChangeForLoss()
Definition: G4VEmModel.cc:118
std::vector< G4Element * > G4ElementVector
std::ostringstream G4ExceptionDescription
Definition: globals.hh:76
G4double GetKineticEnergy() const
void InitialiseElementSelectors(const G4ParticleDefinition *, const G4DataVector &)
Definition: G4VEmModel.cc:146
const G4ParticleDefinition * particle
G4double HighEnergyLimit() const
Definition: G4VEmModel.hh:636
static const G4double ak2
G4MuPairProductionModel(const G4ParticleDefinition *p=nullptr, const G4String &nam="muPairProd")
const char * p
Definition: xmltok.h:285
G4double GetZ() const
Definition: G4Element.hh:131
G4bool IsMaster() const
Definition: G4VEmModel.hh:725
void SetParticle(const G4ParticleDefinition *)
tuple x
Definition: test.py:50
G4ElementData * GetElementData()
Definition: G4VEmModel.hh:826
virtual void SampleSecondaries(std::vector< G4DynamicParticle * > *, const G4MaterialCutsCouple *, const G4DynamicParticle *, G4double tmin, G4double maxEnergy) override
const G4ElementVector * GetElementVector() const
Definition: G4Material.hh:190
int G4int
Definition: G4Types.hh:78
static G4NistManager * Instance()
static constexpr double twopi
Definition: G4SIunits.hh:76
static constexpr double electron_mass_c2
G4double ComputMuPairLoss(G4double Z, G4double tkin, G4double cut, G4double tmax)
static constexpr double TeV
Definition: G4SIunits.hh:218
void PutValue(size_t idx, size_t idy, G4double value)
string material
Definition: eplot.py:19
G4double ComputeMicroscopicCrossSection(G4double tkin, G4double Z, G4double cut)
static const G4int zdat[5]
void InitialiseForElement(G4int Z, G4PhysicsVector *v)
#define G4UniformRand()
Definition: Randomize.hh:97
virtual G4double ComputeDEDXPerVolume(const G4Material *, const G4ParticleDefinition *, G4double kineticEnergy, G4double cutEnergy) override
static const G4double sqrte
const G4ThreeVector & GetMomentumDirection() const
static const G4double ak1
Hep3Vector & rotateUz(const Hep3Vector &)
Definition: ThreeVector.cc:38
void SetProposedKineticEnergy(G4double proposedKinEnergy)
G4double GetLOGZ(G4int Z) const
float electron_mass_c2
Definition: hepunit.py:274
std::vector< G4EmElementSelector * > * GetElementSelectors()
Definition: G4VEmModel.hh:809
void SetProposedMomentumDirection(const G4ThreeVector &dir)
const G4double * GetAtomicNumDensityVector() const
Definition: G4Material.hh:216
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
G4double G4Log(G4double x)
Definition: G4Log.hh:230
G4double G4Exp(G4double initial_x)
Exponential Function double precision.
Definition: G4Exp.hh:183
static G4Positron * Positron()
Definition: G4Positron.cc:94
void SetElementSelectors(std::vector< G4EmElementSelector * > *)
Definition: G4VEmModel.hh:817
G4double GetPDGMass() const
void PutX(size_t idx, G4double value)
T max(const T t1, const T t2)
brief Return the largest of the two arguments
Hep3Vector unit() const
virtual G4double ComputeDMicroscopicCrossSection(G4double tkin, G4double Z, G4double pairEnergy)
T min(const T t1, const T t2)
brief Return the smallest of the two arguments
static const G4double wgi[8]
static constexpr double GeV
Definition: G4SIunits.hh:217
static G4Electron * Electron()
Definition: G4Electron.cc:94
static const G4double fac
static constexpr double MeV
Definition: G4SIunits.hh:214
const G4String & GetName() const
Definition: G4VEmModel.hh:802
static constexpr double pi
Definition: G4SIunits.hh:75
size_t GetNumberOfElements() const
Definition: G4Material.hh:186
static const G4int imax
double G4double
Definition: G4Types.hh:76
void PutY(size_t idy, G4double value)
G4ElementData * fElementData
Definition: G4VEmModel.hh:423
virtual G4double ComputeCrossSectionPerAtom(const G4ParticleDefinition *, G4double kineticEnergy, G4double Z, G4double A, G4double cutEnergy, G4double maxEnergy) override
const G4Element * SelectRandomAtom(const G4MaterialCutsCouple *, const G4ParticleDefinition *, G4double kineticEnergy, G4double cutEnergy=0.0, G4double maxEnergy=DBL_MAX)
Definition: G4VEmModel.hh:544
virtual void InitialiseLocal(const G4ParticleDefinition *, G4VEmModel *masterModel) override
G4ThreeVector GetMomentum() const
G4fissionEvent * fe
virtual G4double MinPrimaryEnergy(const G4Material *, const G4ParticleDefinition *, G4double) override