Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/processes/electromagnetic/muons/src/G4MuPairProductionModel.cc

Version: [ ReleaseNotes ] [ 1.0 ] [ 1.1 ] [ 2.0 ] [ 3.0 ] [ 3.1 ] [ 3.2 ] [ 4.0 ] [ 4.0.p1 ] [ 4.0.p2 ] [ 4.1 ] [ 4.1.p1 ] [ 5.0 ] [ 5.0.p1 ] [ 5.1 ] [ 5.1.p1 ] [ 5.2 ] [ 5.2.p1 ] [ 5.2.p2 ] [ 6.0 ] [ 6.0.p1 ] [ 6.1 ] [ 6.2 ] [ 6.2.p1 ] [ 6.2.p2 ] [ 7.0 ] [ 7.0.p1 ] [ 7.1 ] [ 7.1.p1 ] [ 8.0 ] [ 8.0.p1 ] [ 8.1 ] [ 8.1.p1 ] [ 8.1.p2 ] [ 8.2 ] [ 8.2.p1 ] [ 8.3 ] [ 8.3.p1 ] [ 8.3.p2 ] [ 9.0 ] [ 9.0.p1 ] [ 9.0.p2 ] [ 9.1 ] [ 9.1.p1 ] [ 9.1.p2 ] [ 9.1.p3 ] [ 9.2 ] [ 9.2.p1 ] [ 9.2.p2 ] [ 9.2.p3 ] [ 9.2.p4 ] [ 9.3 ] [ 9.3.p1 ] [ 9.3.p2 ] [ 9.4 ] [ 9.4.p1 ] [ 9.4.p2 ] [ 9.4.p3 ] [ 9.4.p4 ] [ 9.5 ] [ 9.5.p1 ] [ 9.5.p2 ] [ 9.6 ] [ 9.6.p1 ] [ 9.6.p2 ] [ 9.6.p3 ] [ 9.6.p4 ] [ 10.0 ] [ 10.0.p1 ] [ 10.0.p2 ] [ 10.0.p3 ] [ 10.0.p4 ] [ 10.1 ] [ 10.1.p1 ] [ 10.1.p2 ] [ 10.1.p3 ] [ 10.2 ] [ 10.2.p1 ] [ 10.2.p2 ] [ 10.2.p3 ] [ 10.3 ] [ 10.3.p1 ] [ 10.3.p2 ] [ 10.3.p3 ] [ 10.4 ] [ 10.4.p1 ] [ 10.4.p2 ] [ 10.4.p3 ] [ 10.5 ] [ 10.5.p1 ] [ 10.6 ] [ 10.6.p1 ] [ 10.6.p2 ] [ 10.6.p3 ] [ 10.7 ] [ 10.7.p1 ] [ 10.7.p2 ] [ 10.7.p3 ] [ 10.7.p4 ] [ 11.0 ] [ 11.0.p1 ] [ 11.0.p2 ] [ 11.0.p3, ] [ 11.0.p4 ] [ 11.1 ] [ 11.1.1 ] [ 11.1.2 ] [ 11.1.3 ] [ 11.2 ] [ 11.2.1 ] [ 11.2.2 ] [ 11.3.0 ]

  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 "G4ElementDataRegistry.hh"
 84 #include "G4ProductionCutsTable.hh"
 85 #include "G4ParticleChangeForLoss.hh"
 86 #include "G4ModifiedMephi.hh"
 87 #include "G4Log.hh"
 88 #include "G4Exp.hh"
 89 #include "G4AutoLock.hh"
 90 
 91 #include <iostream>
 92 #include <fstream>
 93 
 94 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 95 
 96 const G4int G4MuPairProductionModel::ZDATPAIR[] = {1, 4, 13, 29, 92};
 97 
 98 const G4double G4MuPairProductionModel::xgi[] = {
 99     0.0198550717512320, 0.1016667612931865, 0.2372337950418355, 0.4082826787521750,
100     0.5917173212478250, 0.7627662049581645, 0.8983332387068135, 0.9801449282487680
101   };
102 
103 const G4double G4MuPairProductionModel::wgi[] = {
104     0.0506142681451880, 0.1111905172266870, 0.1568533229389435, 0.1813418916891810,
105     0.1813418916891810, 0.1568533229389435, 0.1111905172266870, 0.0506142681451880
106   };
107 
108 namespace
109 {
110   G4Mutex theMuPairMutex = G4MUTEX_INITIALIZER;
111 
112   const G4double ak1 = 6.9;
113   const G4double ak2 = 1.0;
114 }
115 
116 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
117 
118 G4MuPairProductionModel::G4MuPairProductionModel(const G4ParticleDefinition* p,
119                                                  const G4String& nam)
120   : G4VEmModel(nam),
121   factorForCross(CLHEP::fine_structure_const*CLHEP::fine_structure_const*
122      CLHEP::classic_electr_radius*CLHEP::classic_electr_radius*
123      4./(3.*CLHEP::pi)),
124   sqrte(std::sqrt(G4Exp(1.))),
125   minPairEnergy(4.*CLHEP::electron_mass_c2),
126   lowestKinEnergy(0.85*CLHEP::GeV)
127 {
128   nist = G4NistManager::Instance();
129 
130   theElectron = G4Electron::Electron();
131   thePositron = G4Positron::Positron();
132 
133   if(nullptr != p) { 
134     SetParticle(p); 
135     lowestKinEnergy = std::max(lowestKinEnergy, p->GetPDGMass()*8.0);  
136   }
137   emin = lowestKinEnergy;
138   emax = emin*10000.;
139   SetAngularDistribution(new G4ModifiedMephi());
140 }
141 
142 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
143 
144 G4double G4MuPairProductionModel::MinPrimaryEnergy(const G4Material*,
145                                                    const G4ParticleDefinition*,
146                                                    G4double cut)
147 {
148   return std::max(lowestKinEnergy, cut);
149 }
150 
151 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
152 
153 void G4MuPairProductionModel::Initialise(const G4ParticleDefinition* p,
154                                          const G4DataVector& cuts)
155 { 
156   SetParticle(p); 
157 
158   if (nullptr == fParticleChange) { 
159     fParticleChange = GetParticleChangeForLoss();
160 
161     // define scale of internal table for each thread only once
162     if (0 == nbine) {
163       emin = std::max(lowestKinEnergy, LowEnergyLimit());
164       emax = std::max(HighEnergyLimit(), emin*2);
165       nbine = std::size_t(nYBinPerDecade*std::log10(emax/emin));
166       if(nbine < 3) { nbine = 3; }
167 
168       ymin = G4Log(minPairEnergy/emin);
169       dy = -ymin/G4double(nbiny);
170     }
171     if (p == particle) {
172       G4int pdg = std::abs(p->GetPDGEncoding());
173       if (pdg == 2212) {
174   dataName = "pEEPairProd";
175       } else if (pdg == 321) {
176   dataName = "kaonEEPairProd";
177       } else if (pdg == 211) {
178   dataName = "pionEEPairProd";
179       } else if (pdg == 11) {
180   dataName = "eEEPairProd";
181       } else if (pdg == 13) {
182         if (GetName() == "muToMuonPairProd") {
183           dataName = "muMuMuPairProd";
184   } else {
185     dataName = "muEEPairProd";
186   }
187       } 
188     }
189   }
190 
191   // for low-energy application this process should not work
192   if(lowestKinEnergy >= HighEnergyLimit()) { return; }
193 
194   if (p == particle) {
195     fElementData =
196       G4ElementDataRegistry::Instance()->GetElementDataByName(dataName);
197     if (nullptr == fElementData) { 
198       G4AutoLock l(&theMuPairMutex);
199       fElementData =
200   G4ElementDataRegistry::Instance()->GetElementDataByName(dataName);
201       if (nullptr == fElementData) { 
202   fElementData = new G4ElementData(NZDATPAIR);
203   fElementData->SetName(dataName);
204       }
205       G4bool useDataFile = G4EmParameters::Instance()->RetrieveMuDataFromFile();
206       if (useDataFile)  { useDataFile = RetrieveTables(); }
207       if (!useDataFile) { MakeSamplingTables(); }
208       if (fTableToFile) { StoreTables(); }
209       l.unlock();
210     }
211     if (IsMaster()) {
212       InitialiseElementSelectors(p, cuts);
213     }
214   }
215 }
216 
217 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
218 
219 void G4MuPairProductionModel::InitialiseLocal(const G4ParticleDefinition* p,
220                                               G4VEmModel* masterModel)
221 {
222   if(p == particle && lowestKinEnergy < HighEnergyLimit()) {
223     SetElementSelectors(masterModel->GetElementSelectors());
224   }
225 }
226 
227 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
228 
229 G4double G4MuPairProductionModel::ComputeDEDXPerVolume(
230                                               const G4Material* material,
231                                               const G4ParticleDefinition*,
232                                                     G4double kineticEnergy,
233                                                     G4double cutEnergy)
234 {
235   G4double dedx = 0.0;
236   if (cutEnergy <= minPairEnergy || kineticEnergy <= lowestKinEnergy)
237     { return dedx; }
238 
239   const G4ElementVector* theElementVector = material->GetElementVector();
240   const G4double* theAtomicNumDensityVector =
241                                    material->GetAtomicNumDensityVector();
242 
243   //  loop for elements in the material
244   for (std::size_t i=0; i<material->GetNumberOfElements(); ++i) {
245      G4double Z = (*theElementVector)[i]->GetZ();
246      G4double tmax = MaxSecondaryEnergyForElement(kineticEnergy, Z);
247      G4double loss = ComputMuPairLoss(Z, kineticEnergy, cutEnergy, tmax);
248      dedx += loss*theAtomicNumDensityVector[i];
249   }
250   dedx = std::max(dedx, 0.0);
251   return dedx;
252 }
253 
254 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
255 
256 G4double G4MuPairProductionModel::ComputMuPairLoss(G4double Z, 
257                                                    G4double tkin,
258                                                    G4double cutEnergy, 
259                                                    G4double tmax)
260 {
261   G4double loss = 0.0;
262 
263   G4double cut = std::min(cutEnergy, tmax);
264   if(cut <= minPairEnergy) { return loss; }
265 
266   // calculate the rectricted loss
267   // numerical integration in log(PairEnergy)
268   G4double aaa = G4Log(minPairEnergy);
269   G4double bbb = G4Log(cut);
270 
271   G4int kkk = std::min(std::max(G4lrint((bbb-aaa)/ak1 + ak2), 8), 1);
272   G4double hhh = (bbb-aaa)/kkk;
273   G4double x = aaa;
274 
275   for (G4int l=0 ; l<kkk; ++l) {
276     for (G4int ll=0; ll<NINTPAIR; ++ll) {
277       G4double ep = G4Exp(x+xgi[ll]*hhh);
278       loss += wgi[ll]*ep*ep*ComputeDMicroscopicCrossSection(tkin, Z, ep);
279     }
280     x += hhh;
281   }
282   loss *= hhh;
283   loss = std::max(loss, 0.0);
284   return loss;
285 }
286 
287 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
288 
289 G4double G4MuPairProductionModel::ComputeMicroscopicCrossSection(
290                                            G4double tkin,
291                                            G4double Z,
292                                            G4double cutEnergy)
293 {
294   G4double cross = 0.;
295   G4double tmax = MaxSecondaryEnergyForElement(tkin, Z);
296   G4double cut  = std::max(cutEnergy, minPairEnergy);
297   if (tmax <= cut) { return cross; }
298 
299   G4double aaa = G4Log(cut);
300   G4double bbb = G4Log(tmax);
301   G4int kkk = std::min(std::max(G4lrint((bbb-aaa)/ak1 + ak2), 8), 1);
302 
303   G4double hhh = (bbb-aaa)/(kkk);
304   G4double x = aaa;
305 
306   for (G4int l=0; l<kkk; ++l) {
307     for (G4int i=0; i<NINTPAIR; ++i) {
308       G4double ep = G4Exp(x + xgi[i]*hhh);
309       cross += ep*wgi[i]*ComputeDMicroscopicCrossSection(tkin, Z, ep);
310     }
311     x += hhh;
312   }
313 
314   cross *= hhh;
315   cross = std::max(cross, 0.0);
316   return cross;
317 }
318 
319 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
320 
321 G4double G4MuPairProductionModel::ComputeDMicroscopicCrossSection(
322                                            G4double tkin,
323                                            G4double Z,
324                                            G4double pairEnergy)
325 // Calculates the  differential (D) microscopic cross section
326 // using the cross section formula of R.P. Kokoulin (18/01/98)
327 // Code modified by R.P. Kokoulin, V.N. Ivanchenko (27/01/04)
328 {
329   static const G4double bbbtf= 183. ;
330   static const G4double bbbh = 202.4 ;
331   static const G4double g1tf = 1.95e-5 ;
332   static const G4double g2tf = 5.3e-5 ;
333   static const G4double g1h  = 4.4e-5 ;
334   static const G4double g2h  = 4.8e-5 ;
335 
336   if (pairEnergy <= minPairEnergy)
337     return 0.0;
338 
339   G4double totalEnergy  = tkin + particleMass;
340   G4double residEnergy  = totalEnergy - pairEnergy;
341 
342   if (residEnergy <= 0.75*sqrte*z13*particleMass)
343     return 0.0;
344 
345   G4double a0 = 1.0 / (totalEnergy * residEnergy);
346   G4double alf = 4.0 * electron_mass_c2 / pairEnergy;
347   G4double rt = std::sqrt(1.0 - alf);
348   G4double delta = 6.0 * particleMass * particleMass * a0;
349   G4double tmnexp = alf/(1.0 + rt) + delta*rt;
350 
351   if(tmnexp >= 1.0) { return 0.0; }
352 
353   G4double tmn = G4Log(tmnexp);
354 
355   G4double massratio = particleMass/CLHEP::electron_mass_c2;
356   G4double massratio2 = massratio*massratio;
357   G4double inv_massratio2 = 1.0 / massratio2;
358 
359   // zeta calculation
360   G4double bbb,g1,g2;
361   if( Z < 1.5 ) { bbb = bbbh ; g1 = g1h ; g2 = g2h ; }
362   else          { bbb = bbbtf; g1 = g1tf; g2 = g2tf; }
363 
364   G4double zeta = 0.0;
365   G4double z1exp = totalEnergy / (particleMass + g1*z23*totalEnergy);
366 
367   // 35.221047195922 is the root of zeta1(x) = 0.073 * log(x) - 0.26, so the
368   // condition below is the same as zeta1 > 0.0, but without calling log(x)
369   if (z1exp > 35.221047195922)
370   {
371     G4double z2exp = totalEnergy / (particleMass + g2*z13*totalEnergy);
372     zeta = (0.073 * G4Log(z1exp) - 0.26) / (0.058 * G4Log(z2exp) - 0.14);
373   }
374 
375   G4double z2 = Z*(Z+zeta);
376   G4double screen0 = 2.*electron_mass_c2*sqrte*bbb/(z13*pairEnergy);
377   G4double beta = 0.5*pairEnergy*pairEnergy*a0;
378   G4double xi0 = 0.5*massratio2*beta;
379 
380   // Gaussian integration in ln(1-ro) ( with 8 points)
381   G4double rho[NINTPAIR];
382   G4double rho2[NINTPAIR];
383   G4double xi[NINTPAIR];
384   G4double xi1[NINTPAIR];
385   G4double xii[NINTPAIR];
386 
387   for (G4int i = 0; i < NINTPAIR; ++i)
388   {
389     rho[i] = G4Exp(tmn*xgi[i]) - 1.0; // rho = -asymmetry
390     rho2[i] = rho[i] * rho[i];
391     xi[i] = xi0*(1.0-rho2[i]);
392     xi1[i] = 1.0 + xi[i];
393     xii[i] = 1.0 / xi[i];
394   }
395 
396   G4double ye1[NINTPAIR];
397   G4double ym1[NINTPAIR];
398 
399   G4double b40 = 4.0 * beta;
400   G4double b62 = 6.0 * beta + 2.0;
401 
402   for (G4int i = 0; i < NINTPAIR; ++i)
403   {
404     G4double yeu = (b40 + 5.0) + (b40 - 1.0) * rho2[i];
405     G4double yed = b62*G4Log(3.0 + xii[i]) + (2.0 * beta - 1.0)*rho2[i] - b40;
406 
407     G4double ymu = b62 * (1.0 + rho2[i]) + 6.0;
408     G4double ymd = (b40 + 3.0)*(1.0 + rho2[i])*G4Log(3.0 + xi[i])
409       + 2.0 - 3.0 * rho2[i];
410 
411     ye1[i] = 1.0 + yeu / yed;
412     ym1[i] = 1.0 + ymu / ymd;
413   }
414 
415   G4double be[NINTPAIR];
416   G4double bm[NINTPAIR];
417 
418   for(G4int i = 0; i < NINTPAIR; ++i) {
419     if(xi[i] <= 1000.0) {
420       be[i] = ((2.0 + rho2[i])*(1.0 + beta) +
421          xi[i]*(3.0 + rho2[i]))*G4Log(1.0 + xii[i]) +
422   (1.0 - rho2[i] - beta)/xi1[i] - (3.0 + rho2[i]);
423     } else {
424       be[i] = 0.5*(3.0 - rho2[i] + 2.0*beta*(1.0 + rho2[i]))*xii[i];
425     }
426 
427     if(xi[i] >= 0.001) {
428       G4double a10 = (1.0 + 2.0 * beta) * (1.0 - rho2[i]);
429       bm[i] = ((1.0 + rho2[i])*(1.0 + 1.5 * beta) - a10*xii[i])*G4Log(xi1[i]) +
430                 xi[i] * (1.0 - rho2[i] - beta)/xi1[i] + a10;
431     } else {
432       bm[i] = 0.5*(5.0 - rho2[i] + beta * (3.0 + rho2[i]))*xi[i];
433     }
434   }
435 
436   G4double sum = 0.0;
437 
438   for (G4int i = 0; i < NINTPAIR; ++i) {
439     G4double screen = screen0*xi1[i]/(1.0 - rho2[i]);
440     G4double ale = G4Log(bbb/z13*std::sqrt(xi1[i]*ye1[i])/(1. + screen*ye1[i]));
441     G4double cre = 0.5*G4Log(1. + 2.25*z23*xi1[i]*ye1[i]*inv_massratio2);
442 
443     G4double fe = (ale-cre)*be[i];
444     fe = std::max(fe, 0.0);
445 
446     G4double alm_crm = G4Log(bbb*massratio/(1.5*z23*(1. + screen*ym1[i])));
447     G4double fm = std::max(alm_crm*bm[i], 0.0)*inv_massratio2;
448 
449     sum += wgi[i]*(1.0 + rho[i])*(fe + fm);
450   }
451 
452   return -tmn*sum*factorForCross*z2*residEnergy/(totalEnergy*pairEnergy);
453 }
454 
455 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
456 
457 G4double G4MuPairProductionModel::ComputeCrossSectionPerAtom(
458                                            const G4ParticleDefinition*,
459                                                  G4double kineticEnergy,
460                                                  G4double Z, G4double,
461                                                  G4double cutEnergy,
462                                                  G4double maxEnergy)
463 {
464   G4double cross = 0.0;
465   if (kineticEnergy <= lowestKinEnergy) { return cross; }
466 
467   G4double maxPairEnergy = MaxSecondaryEnergyForElement(kineticEnergy, Z);
468   G4double tmax = std::min(maxEnergy, maxPairEnergy);
469   G4double cut  = std::max(cutEnergy, minPairEnergy);
470   if (cut >= tmax) { return cross; }
471 
472   cross = ComputeMicroscopicCrossSection(kineticEnergy, Z, cut);
473   if(tmax < kineticEnergy) {
474     cross -= ComputeMicroscopicCrossSection(kineticEnergy, Z, tmax);
475   }
476   return cross;
477 }
478 
479 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
480 
481 void G4MuPairProductionModel::MakeSamplingTables()
482 {
483   G4double factore = G4Exp(G4Log(emax/emin)/G4double(nbine));
484 
485   for (G4int iz=0; iz<NZDATPAIR; ++iz) {
486 
487     G4double Z = ZDATPAIR[iz];
488     G4Physics2DVector* pv = new G4Physics2DVector(nbiny+1,nbine+1);
489     G4double kinEnergy = emin;
490 
491     for (std::size_t it=0; it<=nbine; ++it) {
492 
493       pv->PutY(it, G4Log(kinEnergy/CLHEP::MeV));
494       G4double maxPairEnergy = MaxSecondaryEnergyForElement(kinEnergy, Z);
495       /*
496       G4cout << "it= " << it << " E= " << kinEnergy 
497              << "  " << particle->GetParticleName()   
498              << " maxE= " << maxPairEnergy << "  minE= " << minPairEnergy 
499              << " ymin= " << ymin << G4endl;
500       */
501       G4double coef = G4Log(minPairEnergy/kinEnergy)/ymin;
502       G4double ymax = G4Log(maxPairEnergy/kinEnergy)/coef;
503       G4double fac  = (ymax - ymin)/dy;
504       std::size_t imax   = (std::size_t)fac;
505       fac -= (G4double)imax;
506    
507       G4double xSec = 0.0;
508       G4double x = ymin;
509       /*
510       G4cout << "Z= " << currentZ << " z13= " << z13 
511              << " mE= " << maxPairEnergy << "  ymin= " << ymin 
512              << " dy= " << dy << "  c= " << coef << G4endl;
513       */
514       // start from zero
515       pv->PutValue(0, it, 0.0);
516       if(0 == it) { pv->PutX(nbiny, 0.0); }
517 
518       for (std::size_t i=0; i<nbiny; ++i) {
519 
520         if(0 == it) { pv->PutX(i, x); }
521 
522         if(i < imax) {
523           G4double ep = kinEnergy*G4Exp(coef*(x + dy*0.5));
524 
525           // not multiplied by interval, because table 
526           // will be used only for sampling
527           //G4cout << "i= " << i << " x= " << x << "E= " << kinEnergy  
528           //         << " Egamma= " << ep << G4endl;
529           xSec += ep*ComputeDMicroscopicCrossSection(kinEnergy, Z, ep);
530 
531           // last bin before the kinematic limit
532         } else if(i == imax) {
533           G4double ep = kinEnergy*G4Exp(coef*(x + fac*dy*0.5));
534           xSec += ep*fac*ComputeDMicroscopicCrossSection(kinEnergy, Z, ep);
535         }
536         pv->PutValue(i + 1, it, xSec);
537         x += dy;
538       } 
539       kinEnergy *= factore;
540 
541       // to avoid precision lost
542       if(it+1 == nbine) { kinEnergy = emax; }
543     }
544     fElementData->InitialiseForElement(iz, pv);
545   }
546 }
547 
548 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
549 
550 void G4MuPairProductionModel::SampleSecondaries(
551                               std::vector<G4DynamicParticle*>* vdp, 
552                               const G4MaterialCutsCouple* couple,
553                               const G4DynamicParticle* aDynamicParticle,
554                               G4double tmin,
555                               G4double tmax)
556 {
557   G4double kinEnergy = aDynamicParticle->GetKineticEnergy();
558   //G4cout << "------- G4MuPairProductionModel::SampleSecondaries E(MeV)= " 
559   //         << kinEnergy << "  " 
560   //         << aDynamicParticle->GetDefinition()->GetParticleName() << G4endl;
561   G4double totalEnergy   = kinEnergy + particleMass;
562   G4double totalMomentum = 
563     std::sqrt(kinEnergy*(kinEnergy + 2.0*particleMass));
564 
565   G4ThreeVector partDirection = aDynamicParticle->GetMomentumDirection();
566 
567   // select randomly one element constituing the material
568   const G4Element* anElement = SelectRandomAtom(couple,particle,kinEnergy);
569 
570   // define interval of energy transfer
571   G4double maxPairEnergy = MaxSecondaryEnergyForElement(kinEnergy, 
572                                                         anElement->GetZ());
573   G4double maxEnergy = std::min(tmax, maxPairEnergy);
574   G4double minEnergy = std::max(tmin, minPairEnergy);
575 
576   if (minEnergy >= maxEnergy) { return; }
577   //G4cout << "emin= " << minEnergy << " emax= " << maxEnergy 
578   // << " minPair= " << minPairEnergy << " maxpair= " << maxPairEnergy 
579   //    << " ymin= " << ymin << " dy= " << dy << G4endl;
580 
581   G4double coeff = G4Log(minPairEnergy/kinEnergy)/ymin;
582 
583   // compute limits 
584   G4double yymin = G4Log(minEnergy/kinEnergy)/coeff;
585   G4double yymax = G4Log(maxEnergy/kinEnergy)/coeff;
586  
587   //G4cout << "yymin= " << yymin << "  yymax= " << yymax << G4endl;
588 
589   // units should not be used, bacause table was built without
590   G4double logTkin = G4Log(kinEnergy/CLHEP::MeV);
591 
592   // sample e-e+ energy, pair energy first
593 
594   // select sample table via Z
595   G4int iz1(0), iz2(0);
596   for (G4int iz=0; iz<NZDATPAIR; ++iz) { 
597     if(currentZ == ZDATPAIR[iz]) {
598       iz1 = iz2 = iz; 
599       break;
600     } else if(currentZ < ZDATPAIR[iz]) {
601       iz2 = iz;
602       if(iz > 0) { iz1 = iz-1; }
603       else { iz1 = iz2; }
604       break;
605     } 
606   }
607   if (0 == iz1) { iz1 = iz2 = NZDATPAIR-1; }
608 
609   G4double pairEnergy = 0.0;
610   G4int count = 0;
611   //G4cout << "start loop Z1= " << iz1 << " Z2= " << iz2 << G4endl;
612   do {
613     ++count;
614     // sampling using only one random number
615     G4double rand = G4UniformRand();
616   
617     G4double x = FindScaledEnergy(iz1, rand, logTkin, yymin, yymax);
618     if(iz1 != iz2) {
619       G4double x2 = FindScaledEnergy(iz2, rand, logTkin, yymin, yymax);
620       G4double lz1= nist->GetLOGZ(ZDATPAIR[iz1]);
621       G4double lz2= nist->GetLOGZ(ZDATPAIR[iz2]);
622       //G4cout << count << ".  x= " << x << "  x2= " << x2 
623       //             << " Z1= " << iz1 << " Z2= " << iz2 << G4endl;
624       x += (x2 - x)*(lnZ - lz1)/(lz2 - lz1);
625     }
626     //G4cout << "x= " << x << "  coeff= " << coeff << G4endl;
627     pairEnergy = kinEnergy*G4Exp(x*coeff);
628     
629     // Loop checking, 03-Aug-2015, Vladimir Ivanchenko
630   } while((pairEnergy < minEnergy || pairEnergy > maxEnergy) && 10 > count);
631 
632   //G4cout << "## pairEnergy(GeV)= " << pairEnergy/GeV 
633   //         << " Etot(GeV)= " << totalEnergy/GeV << G4endl; 
634 
635   // sample r=(E+-E-)/pairEnergy  ( uniformly .....)
636   G4double rmax =
637     (1.-6.*particleMass*particleMass/(totalEnergy*(totalEnergy-pairEnergy)))
638     *std::sqrt(1.-minPairEnergy/pairEnergy);
639   G4double r = rmax * (-1.+2.*G4UniformRand()) ;
640 
641   // compute energies from pairEnergy,r
642   G4double eEnergy = (1.-r)*pairEnergy*0.5;
643   G4double pEnergy = pairEnergy - eEnergy;
644 
645   // Sample angles 
646   G4ThreeVector eDirection, pDirection;
647   //
648   GetAngularDistribution()->SamplePairDirections(aDynamicParticle, 
649                                                  eEnergy, pEnergy,
650                                                  eDirection, pDirection);
651   // create G4DynamicParticle object for e+e-
652   eEnergy = std::max(eEnergy - CLHEP::electron_mass_c2, 0.0);
653   pEnergy = std::max(pEnergy - CLHEP::electron_mass_c2, 0.0);
654   auto aParticle1 = new G4DynamicParticle(theElectron,eDirection,eEnergy);
655   auto aParticle2 = new G4DynamicParticle(thePositron,pDirection,pEnergy);
656   // Fill output vector
657   vdp->push_back(aParticle1);
658   vdp->push_back(aParticle2);
659 
660   // primary change
661   kinEnergy -= pairEnergy;
662   partDirection *= totalMomentum;
663   partDirection -= (aParticle1->GetMomentum() + aParticle2->GetMomentum());
664   partDirection = partDirection.unit();
665 
666   // if energy transfer is higher than threshold (very high by default)
667   // then stop tracking the primary particle and create a new secondary
668   if (pairEnergy > SecondaryThreshold()) {
669     fParticleChange->ProposeTrackStatus(fStopAndKill);
670     fParticleChange->SetProposedKineticEnergy(0.0);
671     auto newdp = new G4DynamicParticle(particle, partDirection, kinEnergy);
672     vdp->push_back(newdp);
673   } else { // continue tracking the primary e-/e+ otherwise
674     fParticleChange->SetProposedMomentumDirection(partDirection);
675     fParticleChange->SetProposedKineticEnergy(kinEnergy);
676   }
677   //G4cout << "-- G4MuPairProductionModel::SampleSecondaries done" << G4endl; 
678 }
679 
680 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
681 
682 G4double 
683 G4MuPairProductionModel::FindScaledEnergy(G4int iz, G4double rand,
684             G4double logTkin,
685             G4double yymin, G4double yymax)
686 {
687   G4double res = yymin;
688   G4Physics2DVector* pv = fElementData->GetElement2DData(iz);
689   if (nullptr != pv) { 
690     G4double pmin = pv->Value(yymin, logTkin);
691     G4double pmax = pv->Value(yymax, logTkin);
692     G4double p0   = pv->Value(0.0, logTkin);
693     if(p0 <= 0.0) { DataCorrupted(ZDATPAIR[iz], logTkin); }
694     else { res = pv->FindLinearX((pmin + rand*(pmax - pmin))/p0, logTkin); }
695   } else {
696     DataCorrupted(ZDATPAIR[iz], logTkin); 
697   }
698   return res;
699 }
700 
701 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
702 
703 void G4MuPairProductionModel::DataCorrupted(G4int Z, G4double logTkin) const
704 {
705   G4ExceptionDescription ed;
706   ed << "G4ElementData is not properly initialized Z= " << Z
707      << " Ekin(MeV)= " << G4Exp(logTkin)
708      << " IsMasterThread= " << IsMaster() 
709      << " Model " << GetName();
710   G4Exception("G4MuPairProductionModel::()", "em0033", FatalException, ed, "");
711 }
712 
713 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
714 
715 void G4MuPairProductionModel::StoreTables() const
716 {
717   for (G4int iz=0; iz<NZDATPAIR; ++iz) {
718     G4int Z = ZDATPAIR[iz];
719     G4Physics2DVector* pv = fElementData->GetElement2DData(Z);
720     if(nullptr == pv) { 
721       DataCorrupted(Z, 1.0);
722       return;
723     }
724     std::ostringstream ss;
725     ss << "mupair/" << particle->GetParticleName() << Z << ".dat";
726     std::ofstream outfile(ss.str());
727     pv->Store(outfile);
728   }
729 }
730 
731 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
732 
733 G4bool G4MuPairProductionModel::RetrieveTables()
734 {
735   for (G4int iz=0; iz<NZDATPAIR; ++iz) {
736     G4double Z = ZDATPAIR[iz];
737     G4Physics2DVector* pv = new G4Physics2DVector(nbiny+1,nbine+1);
738     std::ostringstream ss;
739     ss << G4EmParameters::Instance()->GetDirLEDATA() << "/mupair/"
740        << particle->GetParticleName() << Z << ".dat";
741     std::ifstream infile(ss.str(), std::ios::in);
742     if(!pv->Retrieve(infile)) { 
743       delete pv;
744       return false; 
745     }
746     fElementData->InitialiseForElement(iz, pv);
747   }
748   return true;
749 }
750 
751 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
752