ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4MuPairProductionModel.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4MuPairProductionModel.cc
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 //
27 // -------------------------------------------------------------------
28 //
29 // GEANT4 Class file
30 //
31 //
32 // File name: G4MuPairProductionModel
33 //
34 // Author: Vladimir Ivanchenko on base of Laszlo Urban code
35 //
36 // Creation date: 24.06.2002
37 //
38 // Modifications:
39 //
40 // 04-12-02 Change G4DynamicParticle constructor in PostStep (V.Ivanchenko)
41 // 23-12-02 Change interface in order to move to cut per region (V.Ivanchenko)
42 // 24-01-03 Fix for compounds (V.Ivanchenko)
43 // 27-01-03 Make models region aware (V.Ivanchenko)
44 // 13-02-03 Add model (V.Ivanchenko)
45 // 06-06-03 Fix in cross section calculation for high energy (V.Ivanchenko)
46 // 20-10-03 2*xi in ComputeDDMicroscopicCrossSection (R.Kokoulin)
47 // 8 integration points in ComputeDMicroscopicCrossSection
48 // 12-01-04 Take min cut of e- and e+ not its sum (V.Ivanchenko)
49 // 10-02-04 Update parameterisation using R.Kokoulin model (V.Ivanchenko)
50 // 28-04-04 For complex materials repeat calculation of max energy for each
51 // material (V.Ivanchenko)
52 // 01-11-04 Fix bug inside ComputeDMicroscopicCrossSection (R.Kokoulin)
53 // 08-04-05 Major optimisation of internal interfaces (V.Ivantchenko)
54 // 03-08-05 Add SetParticle method (V.Ivantchenko)
55 // 23-10-05 Add protection in sampling of e+e- pair energy needed for
56 // low cuts (V.Ivantchenko)
57 // 13-02-06 Add ComputeCrossSectionPerAtom (mma)
58 // 24-04-07 Add protection in SelectRandomAtom method (V.Ivantchenko)
59 // 12-05-06 Updated sampling (use cut) in SelectRandomAtom (A.Bogdanov)
60 // 11-10-07 Add ignoreCut flag (V.Ivanchenko)
61 
62 //
63 // Class Description:
64 //
65 //
66 // -------------------------------------------------------------------
67 //
68 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
69 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
70 
71 #include "G4MuPairProductionModel.hh"
72 #include "G4PhysicalConstants.hh"
73 #include "G4SystemOfUnits.hh"
74 #include "G4EmParameters.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 #include "G4Log.hh"
87 #include "G4Exp.hh"
88 #include <iostream>
89 #include <fstream>
90 
91 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
92 
93 // static members
94 //
95 static const G4double ak1 = 6.9;
96 static const G4double ak2 = 1.0;
97 static const G4int zdat[5] = {1, 4, 13, 29, 92};
98 
99 const G4double G4MuPairProductionModel::xgi[] =
100  { 0.0199, 0.1017, 0.2372, 0.4083, 0.5917, 0.7628, 0.8983, 0.9801 };
101 const G4double G4MuPairProductionModel::wgi[8] =
102  { 0.0506, 0.1112, 0.1569, 0.1813, 0.1813, 0.1569, 0.1112, 0.0506 };
103 
104 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
105 
106 using namespace std;
107 
108 G4MuPairProductionModel::G4MuPairProductionModel(const G4ParticleDefinition* p,
109  const G4String& nam)
110  : G4VEmModel(nam),
111  particle(nullptr),
112  factorForCross(4.*fine_structure_const*fine_structure_const
113  *classic_electr_radius*classic_electr_radius/(3.*pi)),
114  sqrte(sqrt(G4Exp(1.))),
115  currentZ(0),
116  fParticleChange(nullptr),
117  minPairEnergy(4.*electron_mass_c2),
118  lowestKinEnergy(1.0*GeV),
119  nzdat(5),
120  nYBinPerDecade(4),
121  nbiny(1000),
122  nbine(0),
123  ymin(-5.),
124  dy(0.005),
125  fTableToFile(false)
126 {
127  nist = G4NistManager::Instance();
128 
129  theElectron = G4Electron::Electron();
130  thePositron = G4Positron::Positron();
131 
132  particleMass = lnZ = z13 = z23 = 0.;
133 
134  // setup lowest limit dependent on particle mass
135  if(p) {
136  SetParticle(p);
137  lowestKinEnergy = std::max(lowestKinEnergy,p->GetPDGMass()*8.0);
138  }
139  emin = lowestKinEnergy;
140  emax = 10.*TeV;
141 }
142 
143 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
144 
145 G4MuPairProductionModel::~G4MuPairProductionModel()
146 {}
147 
148 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
149 
150 G4double G4MuPairProductionModel::MinPrimaryEnergy(const G4Material*,
151  const G4ParticleDefinition*,
152  G4double cut)
153 {
154  return std::max(lowestKinEnergy,cut);
155 }
156 
157 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
158 
159 void G4MuPairProductionModel::Initialise(const G4ParticleDefinition* p,
160  const G4DataVector& cuts)
161 {
162  SetParticle(p);
163  if(!fParticleChange) { fParticleChange = GetParticleChangeForLoss(); }
164 
165  // for low-energy application this process should not work
166  if(lowestKinEnergy >= HighEnergyLimit()) { return; }
167 
168  // define scale of internal table for each thread only once
169  if(0 == nbine) {
170  emin = std::max(lowestKinEnergy, LowEnergyLimit());
171  emax = std::max(HighEnergyLimit(), emin*2);
172  nbine = size_t(nYBinPerDecade*std::log10(emax/emin));
173  if(nbine < 3) { nbine = 3; }
174 
175  ymin = G4Log(minPairEnergy/emin);
176  dy = -ymin/G4double(nbiny);
177  }
178 
179  if(IsMaster() && p == particle) {
180  if(!fElementData) {
181  fElementData = new G4ElementData();
182  G4bool dataFile = G4EmParameters::Instance()->RetrieveMuDataFromFile();
183  if(dataFile) { dataFile = RetrieveTables(); }
184  if(!dataFile) { MakeSamplingTables(); }
185  if(fTableToFile) { StoreTables(); }
186  }
187  InitialiseElementSelectors(p, cuts);
188  }
189 }
190 
191 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
192 
193 void G4MuPairProductionModel::InitialiseLocal(const G4ParticleDefinition* p,
194  G4VEmModel* masterModel)
195 {
196  if(p == particle && lowestKinEnergy < HighEnergyLimit()) {
197  SetElementSelectors(masterModel->GetElementSelectors());
198  fElementData = masterModel->GetElementData();
199  }
200 }
201 
202 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
203 
204 G4double G4MuPairProductionModel::ComputeDEDXPerVolume(
205  const G4Material* material,
206  const G4ParticleDefinition*,
207  G4double kineticEnergy,
208  G4double cutEnergy)
209 {
210  G4double dedx = 0.0;
211  if (cutEnergy <= minPairEnergy || kineticEnergy <= lowestKinEnergy)
212  { return dedx; }
213 
214  const G4ElementVector* theElementVector = material->GetElementVector();
215  const G4double* theAtomicNumDensityVector =
216  material->GetAtomicNumDensityVector();
217 
218  // loop for elements in the material
219  for (size_t i=0; i<material->GetNumberOfElements(); ++i) {
220  G4double Z = (*theElementVector)[i]->GetZ();
221  G4double tmax = MaxSecondaryEnergyForElement(kineticEnergy, Z);
222  G4double loss = ComputMuPairLoss(Z, kineticEnergy, cutEnergy, tmax);
223  dedx += loss*theAtomicNumDensityVector[i];
224  }
225  dedx = std::max(dedx, 0.0);
226  return dedx;
227 }
228 
229 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
230 
231 G4double G4MuPairProductionModel::ComputMuPairLoss(G4double Z,
232  G4double tkin,
233  G4double cutEnergy,
234  G4double tmax)
235 {
236  G4double loss = 0.0;
237 
238  G4double cut = std::min(cutEnergy,tmax);
239  if(cut <= minPairEnergy) { return loss; }
240 
241  // calculate the rectricted loss
242  // numerical integration in log(PairEnergy)
243  G4double aaa = G4Log(minPairEnergy);
244  G4double bbb = G4Log(cut);
245 
246  G4int kkk = (G4int)((bbb-aaa)/ak1+ak2);
247  if (kkk > 8) { kkk = 8; }
248  else if (kkk < 1) { kkk = 1; }
249 
250  G4double hhh = (bbb-aaa)/(G4double)kkk;
251  G4double x = aaa;
252 
253  for (G4int l=0 ; l<kkk; ++l) {
254  for (G4int ll=0; ll<8; ++ll)
255  {
256  G4double ep = G4Exp(x+xgi[ll]*hhh);
257  loss += wgi[ll]*ep*ep*ComputeDMicroscopicCrossSection(tkin, Z, ep);
258  }
259  x += hhh;
260  }
261  loss *= hhh;
262  loss = std::max(loss, 0.0);
263  return loss;
264 }
265 
266 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
267 
268 G4double G4MuPairProductionModel::ComputeMicroscopicCrossSection(
269  G4double tkin,
270  G4double Z,
271  G4double cutEnergy)
272 {
273  G4double cross = 0.;
274  G4double tmax = MaxSecondaryEnergyForElement(tkin, Z);
275  G4double cut = std::max(cutEnergy, minPairEnergy);
276  if (tmax <= cut) { return cross; }
277 
278  // G4double ak1=6.9 ;
279  // G4double ak2=1.0 ;
280  G4double aaa = G4Log(cut);
281  G4double bbb = G4Log(tmax);
282  G4int kkk = (G4int)((bbb-aaa)/ak1 + ak2);
283  if(kkk > 8) { kkk = 8; }
284  else if (kkk < 1) { kkk = 1; }
285 
286  G4double hhh = (bbb-aaa)/G4double(kkk);
287  G4double x = aaa;
288 
289  for(G4int l=0; l<kkk; ++l)
290  {
291  for(G4int i=0; i<8; ++i)
292  {
293  G4double ep = G4Exp(x + xgi[i]*hhh);
294  cross += ep*wgi[i]*ComputeDMicroscopicCrossSection(tkin, Z, ep);
295  }
296  x += hhh;
297  }
298 
299  cross *= hhh;
300  cross = std::max(cross, 0.0);
301  return cross;
302 }
303 
304 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
305 
306 G4double G4MuPairProductionModel::ComputeDMicroscopicCrossSection(
307  G4double tkin,
308  G4double Z,
309  G4double pairEnergy)
310 // Calculates the differential (D) microscopic cross section
311 // using the cross section formula of R.P. Kokoulin (18/01/98)
312 // Code modified by R.P. Kokoulin, V.N. Ivanchenko (27/01/04)
313 {
314  static const G4double bbbtf= 183. ;
315  static const G4double bbbh = 202.4 ;
316  static const G4double g1tf = 1.95e-5 ;
317  static const G4double g2tf = 5.3e-5 ;
318  static const G4double g1h = 4.4e-5 ;
319  static const G4double g2h = 4.8e-5 ;
320 
321  G4double totalEnergy = tkin + particleMass;
322  G4double residEnergy = totalEnergy - pairEnergy;
323  G4double massratio = particleMass/electron_mass_c2;
324  G4double massratio2 = massratio*massratio;
325  G4double cross = 0.;
326 
327  G4double c3 = 0.75*sqrte*particleMass;
328  if (residEnergy <= c3*z13) { return cross; }
329 
330  static const G4double c7 = 4.*CLHEP::electron_mass_c2;
331  G4double c8 = 6.*particleMass*particleMass;
332  G4double alf = c7/pairEnergy;
333  G4double a3 = 1. - alf;
334  if (a3 <= 0.) { return cross; }
335 
336  // zeta calculation
337  G4double bbb,g1,g2;
338  if( Z < 1.5 ) { bbb = bbbh ; g1 = g1h ; g2 = g2h ; }
339  else { bbb = bbbtf; g1 = g1tf; g2 = g2tf; }
340 
341  G4double zeta = 0;
342  G4double zeta1 =
343  0.073*G4Log(totalEnergy/(particleMass+g1*z23*totalEnergy))-0.26;
344  if ( zeta1 > 0.)
345  {
346  G4double zeta2 =
347  0.058*G4Log(totalEnergy/(particleMass+g2*z13*totalEnergy))-0.14;
348  zeta = zeta1/zeta2 ;
349  }
350 
351  G4double z2 = Z*(Z+zeta);
352  G4double screen0 = 2.*electron_mass_c2*sqrte*bbb/(z13*pairEnergy);
353  G4double a0 = totalEnergy*residEnergy;
354  G4double a1 = pairEnergy*pairEnergy/a0;
355  G4double bet = 0.5*a1;
356  G4double xi0 = 0.25*massratio2*a1;
357  G4double del = c8/a0;
358 
359  G4double rta3 = sqrt(a3);
360  G4double tmnexp = alf/(1. + rta3) + del*rta3;
361  if(tmnexp >= 1.0) { return cross; }
362 
363  G4double tmn = G4Log(tmnexp);
364  G4double sum = 0.;
365 
366  // Gaussian integration in ln(1-ro) ( with 8 points)
367  for (G4int i=0; i<8; ++i)
368  {
369  G4double a4 = G4Exp(tmn*xgi[i]); // a4 = (1.-asymmetry)
370  G4double a5 = a4*(2.-a4) ;
371  G4double a6 = 1.-a5 ;
372  G4double a7 = 1.+a6 ;
373  G4double a9 = 3.+a6 ;
374  G4double xi = xi0*a5 ;
375  G4double xii = 1./xi ;
376  G4double xi1 = 1.+xi ;
377  G4double screen = screen0*xi1/a5 ;
378  G4double yeu = 5.-a6+4.*bet*a7 ;
379  G4double yed = 2.*(1.+3.*bet)*G4Log(3.+xii)-a6-a1*(2.-a6) ;
380  G4double ye1 = 1.+yeu/yed ;
381  G4double ale = G4Log(bbb/z13*sqrt(xi1*ye1)/(1.+screen*ye1)) ;
382  G4double cre = 0.5*G4Log(1.+2.25*z23*xi1*ye1/massratio2) ;
383  G4double be;
384 
385  if (xi <= 1.e3) {
386  be = ((2.+a6)*(1.+bet)+xi*a9)*G4Log(1.+xii)+(a5-bet)/xi1-a9;
387  } else {
388  be = (3.-a6+a1*a7)/(2.*xi);
389  }
390  G4double fe = (ale-cre)*be;
391  if ( fe < 0.) fe = 0. ;
392 
393  G4double ymu = 4.+a6 +3.*bet*a7 ;
394  G4double ymd = a7*(1.5+a1)*G4Log(3.+xi)+1.-1.5*a6 ;
395  G4double ym1 = 1.+ymu/ymd ;
396  G4double alm_crm = G4Log(bbb*massratio/(1.5*z23*(1.+screen*ym1)));
397  G4double a10,bm;
398  if ( xi >= 1.e-3)
399  {
400  a10 = (1.+a1)*a5 ;
401  bm = (a7*(1.+1.5*bet)-a10*xii)*G4Log(xi1)+xi*(a5-bet)/xi1+a10;
402  } else {
403  bm = (5.-a6+bet*a9)*(xi/2.);
404  }
405 
406  G4double fm = alm_crm*bm;
407  if ( fm < 0.) { fm = 0.; }
408 
409  sum += wgi[i]*a4*(fe+fm/massratio2);
410  }
411 
412  cross = -tmn*sum*factorForCross*z2*residEnergy/(totalEnergy*pairEnergy);
413  cross = std::max(cross, 0.0);
414  return cross;
415 }
416 
417 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
418 
419 G4double G4MuPairProductionModel::ComputeCrossSectionPerAtom(
420  const G4ParticleDefinition*,
421  G4double kineticEnergy,
422  G4double Z, G4double,
423  G4double cutEnergy,
424  G4double maxEnergy)
425 {
426  G4double cross = 0.0;
427  if (kineticEnergy <= lowestKinEnergy) { return cross; }
428 
429  G4double maxPairEnergy = MaxSecondaryEnergyForElement(kineticEnergy, Z);
430  G4double tmax = std::min(maxEnergy, maxPairEnergy);
431  G4double cut = std::max(cutEnergy, minPairEnergy);
432  if (cut >= tmax) { return cross; }
433 
434  cross = ComputeMicroscopicCrossSection(kineticEnergy, Z, cut);
435  if(tmax < kineticEnergy) {
436  cross -= ComputeMicroscopicCrossSection(kineticEnergy, Z, tmax);
437  }
438  return cross;
439 }
440 
441 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
442 
443 void G4MuPairProductionModel::MakeSamplingTables()
444 {
445  G4double factore = G4Exp(G4Log(emax/emin)/G4double(nbine));
446 
447  for (G4int iz=0; iz<nzdat; ++iz) {
448 
449  G4double Z = zdat[iz];
450  G4Physics2DVector* pv = new G4Physics2DVector(nbiny+1,nbine+1);
451  G4double kinEnergy = emin;
452 
453  for (size_t it=0; it<=nbine; ++it) {
454 
455  pv->PutY(it, G4Log(kinEnergy/MeV));
456  G4double maxPairEnergy = MaxSecondaryEnergyForElement(kinEnergy, Z);
457  /*
458  G4cout << "it= " << it << " E= " << kinEnergy
459  << " " << particle->GetParticleName()
460  << " maxE= " << maxPairEnergy << " minE= " << minPairEnergy
461  << " ymin= " << ymin << G4endl;
462  */
463  G4double coef = G4Log(minPairEnergy/kinEnergy)/ymin;
464  G4double ymax = G4Log(maxPairEnergy/kinEnergy)/coef;
465  G4double fac = (ymax - ymin)/dy;
466  size_t imax = (size_t)fac;
467  fac -= (G4double)imax;
468 
469  G4double xSec = 0.0;
470  G4double x = ymin;
471  /*
472  G4cout << "Z= " << currentZ << " z13= " << z13
473  << " mE= " << maxPairEnergy << " ymin= " << ymin
474  << " dy= " << dy << " c= " << coef << G4endl;
475  */
476  // start from zero
477  pv->PutValue(0, it, 0.0);
478  if(0 == it) { pv->PutX(nbiny, 0.0); }
479 
480  for (size_t i=0; i<nbiny; ++i) {
481 
482  if(0 == it) { pv->PutX(i, x); }
483 
484  if(i < imax) {
485  G4double ep = kinEnergy*G4Exp(coef*(x + dy*0.5));
486 
487  // not multiplied by interval, because table
488  // will be used only for sampling
489  //G4cout << "i= " << i << " x= " << x << "E= " << kinEnergy
490  // << " Egamma= " << ep << G4endl;
491  xSec += ep*ComputeDMicroscopicCrossSection(kinEnergy, Z, ep);
492 
493  // last bin before the kinematic limit
494  } else if(i == imax) {
495  G4double ep = kinEnergy*G4Exp(coef*(x + fac*dy*0.5));
496  xSec += ep*fac*ComputeDMicroscopicCrossSection(kinEnergy, Z, ep);
497  }
498  pv->PutValue(i + 1, it, xSec);
499  x += dy;
500  }
501  kinEnergy *= factore;
502 
503  // to avoid precision lost
504  if(it+1 == nbine) { kinEnergy = emax; }
505  }
506  fElementData->InitialiseForElement(zdat[iz], pv);
507  }
508 }
509 
510 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
511 
512 void G4MuPairProductionModel::SampleSecondaries(
513  std::vector<G4DynamicParticle*>* vdp,
514  const G4MaterialCutsCouple* couple,
515  const G4DynamicParticle* aDynamicParticle,
516  G4double tmin,
517  G4double tmax)
518 {
519  G4double kineticEnergy = aDynamicParticle->GetKineticEnergy();
520  //G4cout << "------- G4MuPairProductionModel::SampleSecondaries E(MeV)= "
521  // << kineticEnergy << " "
522  // << aDynamicParticle->GetDefinition()->GetParticleName() << G4endl;
523  G4double totalEnergy = kineticEnergy + particleMass;
524  G4double totalMomentum =
525  sqrt(kineticEnergy*(kineticEnergy + 2.0*particleMass));
526 
527  G4ThreeVector partDirection = aDynamicParticle->GetMomentumDirection();
528 
529  // select randomly one element constituing the material
530  const G4Element* anElement = SelectRandomAtom(couple,particle,kineticEnergy);
531 
532  // define interval of energy transfer
533  G4double maxPairEnergy = MaxSecondaryEnergyForElement(kineticEnergy,
534  anElement->GetZ());
535  G4double maxEnergy = std::min(tmax, maxPairEnergy);
536  G4double minEnergy = std::max(tmin, minPairEnergy);
537 
538  if(minEnergy >= maxEnergy) { return; }
539  //G4cout << "emin= " << minEnergy << " emax= " << maxEnergy
540  // << " minPair= " << minPairEnergy << " maxpair= " << maxPairEnergy
541  // << " ymin= " << ymin << " dy= " << dy << G4endl;
542 
543  G4double coeff = G4Log(minPairEnergy/kineticEnergy)/ymin;
544 
545  // compute limits
546  G4double yymin = G4Log(minEnergy/kineticEnergy)/coeff;
547  G4double yymax = G4Log(maxEnergy/kineticEnergy)/coeff;
548 
549  //G4cout << "yymin= " << yymin << " yymax= " << yymax << G4endl;
550 
551  // units should not be used, bacause table was built without
552  G4double logTkin = G4Log(kineticEnergy/MeV);
553 
554  // sample e-e+ energy, pair energy first
555 
556  // select sample table via Z
557  G4int iz1(0), iz2(0);
558  for(G4int iz=0; iz<nzdat; ++iz) {
559  if(currentZ == zdat[iz]) {
560  iz1 = iz2 = currentZ;
561  break;
562  } else if(currentZ < zdat[iz]) {
563  iz2 = zdat[iz];
564  if(iz > 0) { iz1 = zdat[iz-1]; }
565  else { iz1 = iz2; }
566  break;
567  }
568  }
569  if(0 == iz1) { iz1 = iz2 = zdat[nzdat-1]; }
570 
571  G4double PairEnergy = 0.0;
572  G4int count = 0;
573  //G4cout << "start loop Z1= " << iz1 << " Z2= " << iz2 << G4endl;
574  do {
575  ++count;
576  // sampling using only one random number
577  G4double rand = G4UniformRand();
578 
579  G4double x = FindScaledEnergy(iz1, rand, logTkin, yymin, yymax);
580  if(iz1 != iz2) {
581  G4double x2 = FindScaledEnergy(iz2, rand, logTkin, yymin, yymax);
582  G4double lz1= nist->GetLOGZ(iz1);
583  G4double lz2= nist->GetLOGZ(iz2);
584  //G4cout << count << ". x= " << x << " x2= " << x2
585  // << " Z1= " << iz1 << " Z2= " << iz2 << G4endl;
586  x += (x2 - x)*(lnZ - lz1)/(lz2 - lz1);
587  }
588  //G4cout << "x= " << x << " coeff= " << coeff << G4endl;
589  PairEnergy = kineticEnergy*G4Exp(x*coeff);
590 
591  // Loop checking, 03-Aug-2015, Vladimir Ivanchenko
592  } while((PairEnergy < minEnergy || PairEnergy > maxEnergy) && 10 > count);
593 
594  //G4cout << "## PairEnergy(GeV)= " << PairEnergy/GeV
595  // << " Etot(GeV)= " << totalEnergy/GeV << G4endl;
596 
597  // sample r=(E+-E-)/PairEnergy ( uniformly .....)
598  G4double rmax =
599  (1.-6.*particleMass*particleMass/(totalEnergy*(totalEnergy-PairEnergy)))
600  *sqrt(1.-minPairEnergy/PairEnergy);
601  G4double r = rmax * (-1.+2.*G4UniformRand()) ;
602 
603  // compute energies from PairEnergy,r
604  G4double ElectronEnergy = (1.-r)*PairEnergy*0.5;
605  G4double PositronEnergy = PairEnergy - ElectronEnergy;
606 
607  // The angle of the emitted virtual photon is sampled
608  // according to the muon bremsstrahlung model
609 
610  G4double gam = totalEnergy/particleMass;
611  G4double gmax = gam*std::min(1.0, totalEnergy/PairEnergy - 1.0);
612  G4double gmax2= gmax*gmax;
613  G4double x = G4UniformRand()*gmax2/(1.0 + gmax2);
614 
615  G4double theta = sqrt(x/(1.0 - x))/gam;
616  G4double sint = sin(theta);
617  G4double phi = twopi * G4UniformRand() ;
618  G4double dirx = sint*cos(phi), diry = sint*sin(phi), dirz = cos(theta) ;
619 
620  G4ThreeVector gDirection(dirx, diry, dirz);
621  gDirection.rotateUz(partDirection);
622 
623  // the angles of e- and e+ assumed to be the same as virtual gamma
624 
625  // create G4DynamicParticle object for the particle1
626  G4double ekin = std::max(ElectronEnergy - electron_mass_c2,0.0);
627  G4DynamicParticle* aParticle1 =
628  new G4DynamicParticle(theElectron, gDirection, ekin);
629 
630  // create G4DynamicParticle object for the particle2
631  ekin = std::max(PositronEnergy - electron_mass_c2,0.0);
632  G4DynamicParticle* aParticle2 =
633  new G4DynamicParticle(thePositron, gDirection, ekin);
634 
635  // primary change
636  kineticEnergy -= (ElectronEnergy + PositronEnergy);
637  fParticleChange->SetProposedKineticEnergy(kineticEnergy);
638 
639  partDirection *= totalMomentum;
640  partDirection -= (aParticle1->GetMomentum() + aParticle2->GetMomentum());
641  partDirection = partDirection.unit();
642  fParticleChange->SetProposedMomentumDirection(partDirection);
643 
644  // add secondary
645  vdp->push_back(aParticle1);
646  vdp->push_back(aParticle2);
647  //G4cout << "-- G4MuPairProductionModel::SampleSecondaries done" << G4endl;
648 }
649 
650 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
651 
652 void G4MuPairProductionModel::DataCorrupted(G4int Z, G4double logTkin) const
653 {
654  G4ExceptionDescription ed;
655  ed << "G4ElementData is not properly initialized Z= " << Z
656  << " Ekin(MeV)= " << G4Exp(logTkin)
657  << " IsMasterThread= " << IsMaster()
658  << " Model " << GetName();
659  G4Exception("G4MuPairProductionModel::()","em0033",FatalException,
660  ed,"");
661 }
662 
663 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
664 
665 void G4MuPairProductionModel::StoreTables() const
666 {
667  for (G4int iz=0; iz<nzdat; ++iz) {
668  G4int Z = zdat[iz];
669  G4Physics2DVector* pv = fElementData->GetElement2DData(Z);
670  if(!pv) { DataCorrupted(Z, 1.0); }
671  std::ostringstream ss;
672  ss << "mupair/" << particle->GetParticleName() << Z << ".dat";
673  std::ofstream outfile(ss.str());
674  pv->Store(outfile);
675  }
676 }
677 
678 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
679 
680 G4bool G4MuPairProductionModel::RetrieveTables()
681 {
682  char* path = std::getenv("G4LEDATA");
683  G4String dir("");
684  if (path) {
685  std::ostringstream ost;
686  ost << path << "/mupair/";
687  dir = ost.str();
688  } else {
689  dir = "./mupair/";
690  }
691 
692  for (G4int iz=0; iz<nzdat; ++iz) {
693  G4double Z = zdat[iz];
694  G4Physics2DVector* pv = new G4Physics2DVector(nbiny+1,nbine+1);
695  if(!pv) {
696  DataCorrupted(Z, 2.0);
697  return false;
698  }
699  std::ostringstream ss;
700  ss << dir << particle->GetParticleName() << Z << ".dat";
701  std::ifstream infile(ss.str(), std::ios::in);
702  if(!pv->Retrieve(infile)) { return false; }
703  fElementData->InitialiseForElement(Z, pv);
704  }
705  return true;
706 }
707 
708 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
709 
710