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 ]

Diff markup

Differences between /processes/electromagnetic/muons/src/G4MuPairProductionModel.cc (Version 11.3.0) and /processes/electromagnetic/muons/src/G4MuPairProductionModel.cc (Version 2.0)


  1 //                                                  1 
  2 // *******************************************    
  3 // * License and Disclaimer                       
  4 // *                                              
  5 // * The  Geant4 software  is  copyright of th    
  6 // * the Geant4 Collaboration.  It is provided    
  7 // * conditions of the Geant4 Software License    
  8 // * LICENSE and available at  http://cern.ch/    
  9 // * include a list of copyright holders.         
 10 // *                                              
 11 // * Neither the authors of this software syst    
 12 // * institutes,nor the agencies providing fin    
 13 // * work  make  any representation or  warran    
 14 // * regarding  this  software system or assum    
 15 // * use.  Please see the license in the file     
 16 // * for the full disclaimer and the limitatio    
 17 // *                                              
 18 // * This  code  implementation is the result     
 19 // * technical work of the GEANT4 collaboratio    
 20 // * By using,  copying,  modifying or  distri    
 21 // * any work based  on the software)  you  ag    
 22 // * use  in  resulting  scientific  publicati    
 23 // * acceptance of all terms of the Geant4 Sof    
 24 // *******************************************    
 25 //                                                
 26 //                                                
 27 // -------------------------------------------    
 28 //                                                
 29 // GEANT4 Class file                              
 30 //                                                
 31 //                                                
 32 // File name:     G4MuPairProductionModel         
 33 //                                                
 34 // Author:        Vladimir Ivanchenko on base     
 35 //                                                
 36 // Creation date: 24.06.2002                      
 37 //                                                
 38 // Modifications:                                 
 39 //                                                
 40 // 04-12-02 Change G4DynamicParticle construct    
 41 // 23-12-02 Change interface in order to move     
 42 // 24-01-03 Fix for compounds (V.Ivanchenko)      
 43 // 27-01-03 Make models region aware (V.Ivanch    
 44 // 13-02-03 Add model (V.Ivanchenko)              
 45 // 06-06-03 Fix in cross section calculation f    
 46 // 20-10-03 2*xi in ComputeDDMicroscopicCrossS    
 47 //          8 integration points in ComputeDMi    
 48 // 12-01-04 Take min cut of e- and e+ not its     
 49 // 10-02-04 Update parameterisation using R.Ko    
 50 // 28-04-04 For complex materials repeat calcu    
 51 //          material (V.Ivanchenko)               
 52 // 01-11-04 Fix bug inside ComputeDMicroscopic    
 53 // 08-04-05 Major optimisation of internal int    
 54 // 03-08-05 Add SetParticle method (V.Ivantche    
 55 // 23-10-05 Add protection in sampling of e+e-    
 56 //          low cuts (V.Ivantchenko)              
 57 // 13-02-06 Add ComputeCrossSectionPerAtom (mm    
 58 // 24-04-07 Add protection in SelectRandomAtom    
 59 // 12-05-06 Updated sampling (use cut) in Sele    
 60 // 11-10-07 Add ignoreCut flag (V.Ivanchenko)     
 61                                                   
 62 //                                                
 63 // Class Description:                             
 64 //                                                
 65 //                                                
 66 // -------------------------------------------    
 67 //                                                
 68 //....oooOO0OOooo........oooOO0OOooo........oo    
 69 //....oooOO0OOooo........oooOO0OOooo........oo    
 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........oo    
 95                                                   
 96 const G4int G4MuPairProductionModel::ZDATPAIR[    
 97                                                   
 98 const G4double G4MuPairProductionModel::xgi[]     
 99     0.0198550717512320, 0.1016667612931865, 0.    
100     0.5917173212478250, 0.7627662049581645, 0.    
101   };                                              
102                                                   
103 const G4double G4MuPairProductionModel::wgi[]     
104     0.0506142681451880, 0.1111905172266870, 0.    
105     0.1813418916891810, 0.1568533229389435, 0.    
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........oo    
117                                                   
118 G4MuPairProductionModel::G4MuPairProductionMod    
119                                                   
120   : G4VEmModel(nam),                              
121   factorForCross(CLHEP::fine_structure_const*C    
122      CLHEP::classic_electr_radius*CLHEP::class    
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    
136   }                                               
137   emin = lowestKinEnergy;                         
138   emax = emin*10000.;                             
139   SetAngularDistribution(new G4ModifiedMephi()    
140 }                                                 
141                                                   
142 //....oooOO0OOooo........oooOO0OOooo........oo    
143                                                   
144 G4double G4MuPairProductionModel::MinPrimaryEn    
145                                                   
146                                                   
147 {                                                 
148   return std::max(lowestKinEnergy, cut);          
149 }                                                 
150                                                   
151 //....oooOO0OOooo........oooOO0OOooo........oo    
152                                                   
153 void G4MuPairProductionModel::Initialise(const    
154                                          const    
155 {                                                 
156   SetParticle(p);                                 
157                                                   
158   if (nullptr == fParticleChange) {               
159     fParticleChange = GetParticleChangeForLoss    
160                                                   
161     // define scale of internal table for each    
162     if (0 == nbine) {                             
163       emin = std::max(lowestKinEnergy, LowEner    
164       emax = std::max(HighEnergyLimit(), emin*    
165       nbine = std::size_t(nYBinPerDecade*std::    
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 s    
192   if(lowestKinEnergy >= HighEnergyLimit()) { r    
193                                                   
194   if (p == particle) {                            
195     fElementData =                                
196       G4ElementDataRegistry::Instance()->GetEl    
197     if (nullptr == fElementData) {                
198       G4AutoLock l(&theMuPairMutex);              
199       fElementData =                              
200   G4ElementDataRegistry::Instance()->GetElemen    
201       if (nullptr == fElementData) {              
202   fElementData = new G4ElementData(NZDATPAIR);    
203   fElementData->SetName(dataName);                
204       }                                           
205       G4bool useDataFile = G4EmParameters::Ins    
206       if (useDataFile)  { useDataFile = Retrie    
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........oo    
218                                                   
219 void G4MuPairProductionModel::InitialiseLocal(    
220                                                   
221 {                                                 
222   if(p == particle && lowestKinEnergy < HighEn    
223     SetElementSelectors(masterModel->GetElemen    
224   }                                               
225 }                                                 
226                                                   
227 //....oooOO0OOooo........oooOO0OOooo........oo    
228                                                   
229 G4double G4MuPairProductionModel::ComputeDEDXP    
230                                                   
231                                                   
232                                                   
233                                                   
234 {                                                 
235   G4double dedx = 0.0;                            
236   if (cutEnergy <= minPairEnergy || kineticEne    
237     { return dedx; }                              
238                                                   
239   const G4ElementVector* theElementVector = ma    
240   const G4double* theAtomicNumDensityVector =     
241                                    material->G    
242                                                   
243   //  loop for elements in the material           
244   for (std::size_t i=0; i<material->GetNumberO    
245      G4double Z = (*theElementVector)[i]->GetZ    
246      G4double tmax = MaxSecondaryEnergyForElem    
247      G4double loss = ComputMuPairLoss(Z, kinet    
248      dedx += loss*theAtomicNumDensityVector[i]    
249   }                                               
250   dedx = std::max(dedx, 0.0);                     
251   return dedx;                                    
252 }                                                 
253                                                   
254 //....oooOO0OOooo........oooOO0OOooo........oo    
255                                                   
256 G4double G4MuPairProductionModel::ComputMuPair    
257                                                   
258                                                   
259                                                   
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-a    
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*ComputeDMicroscopi    
279     }                                             
280     x += hhh;                                     
281   }                                               
282   loss *= hhh;                                    
283   loss = std::max(loss, 0.0);                     
284   return loss;                                    
285 }                                                 
286                                                   
287 //....oooOO0OOooo........oooOO0OOooo........oo    
288                                                   
289 G4double G4MuPairProductionModel::ComputeMicro    
290                                            G4d    
291                                            G4d    
292                                            G4d    
293 {                                                 
294   G4double cross = 0.;                            
295   G4double tmax = MaxSecondaryEnergyForElement    
296   G4double cut  = std::max(cutEnergy, minPairE    
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-a    
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]*ComputeDMicroscopicCr    
310     }                                             
311     x += hhh;                                     
312   }                                               
313                                                   
314   cross *= hhh;                                   
315   cross = std::max(cross, 0.0);                   
316   return cross;                                   
317 }                                                 
318                                                   
319 //....oooOO0OOooo........oooOO0OOooo........oo    
320                                                   
321 G4double G4MuPairProductionModel::ComputeDMicr    
322                                            G4d    
323                                            G4d    
324                                            G4d    
325 // Calculates the  differential (D) microscopi    
326 // using the cross section formula of R.P. Kok    
327 // Code modified by R.P. Kokoulin, V.N. Ivanch    
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 - pairEn    
341                                                   
342   if (residEnergy <= 0.75*sqrte*z13*particleMa    
343     return 0.0;                                   
344                                                   
345   G4double a0 = 1.0 / (totalEnergy * residEner    
346   G4double alf = 4.0 * electron_mass_c2 / pair    
347   G4double rt = std::sqrt(1.0 - alf);             
348   G4double delta = 6.0 * particleMass * partic    
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::ele    
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 =    
362   else          { bbb = bbbtf; g1 = g1tf; g2 =    
363                                                   
364   G4double zeta = 0.0;                            
365   G4double z1exp = totalEnergy / (particleMass    
366                                                   
367   // 35.221047195922 is the root of zeta1(x) =    
368   // condition below is the same as zeta1 > 0.    
369   if (z1exp > 35.221047195922)                    
370   {                                               
371     G4double z2exp = totalEnergy / (particleMa    
372     zeta = (0.073 * G4Log(z1exp) - 0.26) / (0.    
373   }                                               
374                                                   
375   G4double z2 = Z*(Z+zeta);                       
376   G4double screen0 = 2.*electron_mass_c2*sqrte    
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    
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 =    
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) *    
405     G4double yed = b62*G4Log(3.0 + xii[i]) + (    
406                                                   
407     G4double ymu = b62 * (1.0 + rho2[i]) + 6.0    
408     G4double ymd = (b40 + 3.0)*(1.0 + rho2[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 + xi    
422   (1.0 - rho2[i] - beta)/xi1[i] - (3.0 + rho2[    
423     } else {                                      
424       be[i] = 0.5*(3.0 - rho2[i] + 2.0*beta*(1    
425     }                                             
426                                                   
427     if(xi[i] >= 0.001) {                          
428       G4double a10 = (1.0 + 2.0 * beta) * (1.0    
429       bm[i] = ((1.0 + rho2[i])*(1.0 + 1.5 * be    
430                 xi[i] * (1.0 - rho2[i] - beta)    
431     } else {                                      
432       bm[i] = 0.5*(5.0 - rho2[i] + beta * (3.0    
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 - rh    
440     G4double ale = G4Log(bbb/z13*std::sqrt(xi1    
441     G4double cre = 0.5*G4Log(1. + 2.25*z23*xi1    
442                                                   
443     G4double fe = (ale-cre)*be[i];                
444     fe = std::max(fe, 0.0);                       
445                                                   
446     G4double alm_crm = G4Log(bbb*massratio/(1.    
447     G4double fm = std::max(alm_crm*bm[i], 0.0)    
448                                                   
449     sum += wgi[i]*(1.0 + rho[i])*(fe + fm);       
450   }                                               
451                                                   
452   return -tmn*sum*factorForCross*z2*residEnerg    
453 }                                                 
454                                                   
455 //....oooOO0OOooo........oooOO0OOooo........oo    
456                                                   
457 G4double G4MuPairProductionModel::ComputeCross    
458                                            con    
459                                                   
460                                                   
461                                                   
462                                                   
463 {                                                 
464   G4double cross = 0.0;                           
465   if (kineticEnergy <= lowestKinEnergy) { retu    
466                                                   
467   G4double maxPairEnergy = MaxSecondaryEnergyF    
468   G4double tmax = std::min(maxEnergy, maxPairE    
469   G4double cut  = std::max(cutEnergy, minPairE    
470   if (cut >= tmax) { return cross; }              
471                                                   
472   cross = ComputeMicroscopicCrossSection(kinet    
473   if(tmax < kineticEnergy) {                      
474     cross -= ComputeMicroscopicCrossSection(ki    
475   }                                               
476   return cross;                                   
477 }                                                 
478                                                   
479 //....oooOO0OOooo........oooOO0OOooo........oo    
480                                                   
481 void G4MuPairProductionModel::MakeSamplingTabl    
482 {                                                 
483   G4double factore = G4Exp(G4Log(emax/emin)/G4    
484                                                   
485   for (G4int iz=0; iz<NZDATPAIR; ++iz) {          
486                                                   
487     G4double Z = ZDATPAIR[iz];                    
488     G4Physics2DVector* pv = new G4Physics2DVec    
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 = MaxSecondaryEne    
495       /*                                          
496       G4cout << "it= " << it << " E= " << kinE    
497              << "  " << particle->GetParticleN    
498              << " maxE= " << maxPairEnergy <<     
499              << " ymin= " << ymin << G4endl;      
500       */                                          
501       G4double coef = G4Log(minPairEnergy/kinE    
502       G4double ymax = G4Log(maxPairEnergy/kinE    
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= "     
511              << " mE= " << maxPairEnergy << "     
512              << " dy= " << dy << "  c= " << co    
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*(    
524                                                   
525           // not multiplied by interval, becau    
526           // will be used only for sampling       
527           //G4cout << "i= " << i << " x= " <<     
528           //         << " Egamma= " << ep << G    
529           xSec += ep*ComputeDMicroscopicCrossS    
530                                                   
531           // last bin before the kinematic lim    
532         } else if(i == imax) {                    
533           G4double ep = kinEnergy*G4Exp(coef*(    
534           xSec += ep*fac*ComputeDMicroscopicCr    
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........oo    
549                                                   
550 void G4MuPairProductionModel::SampleSecondarie    
551                               std::vector<G4Dy    
552                               const G4Material    
553                               const G4DynamicP    
554                               G4double tmin,      
555                               G4double tmax)      
556 {                                                 
557   G4double kinEnergy = aDynamicParticle->GetKi    
558   //G4cout << "------- G4MuPairProductionModel    
559   //         << kinEnergy << "  "                 
560   //         << aDynamicParticle->GetDefinitio    
561   G4double totalEnergy   = kinEnergy + particl    
562   G4double totalMomentum =                        
563     std::sqrt(kinEnergy*(kinEnergy + 2.0*parti    
564                                                   
565   G4ThreeVector partDirection = aDynamicPartic    
566                                                   
567   // select randomly one element constituing t    
568   const G4Element* anElement = SelectRandomAto    
569                                                   
570   // define interval of energy transfer           
571   G4double maxPairEnergy = MaxSecondaryEnergyF    
572                                                   
573   G4double maxEnergy = std::min(tmax, maxPairE    
574   G4double minEnergy = std::max(tmin, minPairE    
575                                                   
576   if (minEnergy >= maxEnergy) { return; }         
577   //G4cout << "emin= " << minEnergy << " emax=    
578   // << " minPair= " << minPairEnergy << " max    
579   //    << " ymin= " << ymin << " dy= " << dy     
580                                                   
581   G4double coeff = G4Log(minPairEnergy/kinEner    
582                                                   
583   // compute limits                               
584   G4double yymin = G4Log(minEnergy/kinEnergy)/    
585   G4double yymax = G4Log(maxEnergy/kinEnergy)/    
586                                                   
587   //G4cout << "yymin= " << yymin << "  yymax=     
588                                                   
589   // units should not be used, bacause table w    
590   G4double logTkin = G4Log(kinEnergy/CLHEP::Me    
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    
612   do {                                            
613     ++count;                                      
614     // sampling using only one random number      
615     G4double rand = G4UniformRand();              
616                                                   
617     G4double x = FindScaledEnergy(iz1, rand, l    
618     if(iz1 != iz2) {                              
619       G4double x2 = FindScaledEnergy(iz2, rand    
620       G4double lz1= nist->GetLOGZ(ZDATPAIR[iz1    
621       G4double lz2= nist->GetLOGZ(ZDATPAIR[iz2    
622       //G4cout << count << ".  x= " << x << "     
623       //             << " Z1= " << iz1 << " Z2    
624       x += (x2 - x)*(lnZ - lz1)/(lz2 - lz1);      
625     }                                             
626     //G4cout << "x= " << x << "  coeff= " << c    
627     pairEnergy = kinEnergy*G4Exp(x*coeff);        
628                                                   
629     // Loop checking, 03-Aug-2015, Vladimir Iv    
630   } while((pairEnergy < minEnergy || pairEnerg    
631                                                   
632   //G4cout << "## pairEnergy(GeV)= " << pairEn    
633   //         << " Etot(GeV)= " << totalEnergy/    
634                                                   
635   // sample r=(E+-E-)/pairEnergy  ( uniformly     
636   G4double rmax =                                 
637     (1.-6.*particleMass*particleMass/(totalEne    
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()->SamplePairDirectio    
649                                                   
650                                                   
651   // create G4DynamicParticle object for e+e-     
652   eEnergy = std::max(eEnergy - CLHEP::electron    
653   pEnergy = std::max(pEnergy - CLHEP::electron    
654   auto aParticle1 = new G4DynamicParticle(theE    
655   auto aParticle2 = new G4DynamicParticle(theP    
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()     
664   partDirection = partDirection.unit();           
665                                                   
666   // if energy transfer is higher than thresho    
667   // then stop tracking the primary particle a    
668   if (pairEnergy > SecondaryThreshold()) {        
669     fParticleChange->ProposeTrackStatus(fStopA    
670     fParticleChange->SetProposedKineticEnergy(    
671     auto newdp = new G4DynamicParticle(particl    
672     vdp->push_back(newdp);                        
673   } else { // continue tracking the primary e-    
674     fParticleChange->SetProposedMomentumDirect    
675     fParticleChange->SetProposedKineticEnergy(    
676   }                                               
677   //G4cout << "-- G4MuPairProductionModel::Sam    
678 }                                                 
679                                                   
680 //....oooOO0OOooo........oooOO0OOooo........oo    
681                                                   
682 G4double                                          
683 G4MuPairProductionModel::FindScaledEnergy(G4in    
684             G4double logTkin,                     
685             G4double yymin, G4double yymax)       
686 {                                                 
687   G4double res = yymin;                           
688   G4Physics2DVector* pv = fElementData->GetEle    
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]    
694     else { res = pv->FindLinearX((pmin + rand*    
695   } else {                                        
696     DataCorrupted(ZDATPAIR[iz], logTkin);         
697   }                                               
698   return res;                                     
699 }                                                 
700                                                   
701 //....oooOO0OOooo........oooOO0OOooo........oo    
702                                                   
703 void G4MuPairProductionModel::DataCorrupted(G4    
704 {                                                 
705   G4ExceptionDescription ed;                      
706   ed << "G4ElementData is not properly initial    
707      << " Ekin(MeV)= " << G4Exp(logTkin)          
708      << " IsMasterThread= " << IsMaster()         
709      << " Model " << GetName();                   
710   G4Exception("G4MuPairProductionModel::()", "    
711 }                                                 
712                                                   
713 //....oooOO0OOooo........oooOO0OOooo........oo    
714                                                   
715 void G4MuPairProductionModel::StoreTables() co    
716 {                                                 
717   for (G4int iz=0; iz<NZDATPAIR; ++iz) {          
718     G4int Z = ZDATPAIR[iz];                       
719     G4Physics2DVector* pv = fElementData->GetE    
720     if(nullptr == pv) {                           
721       DataCorrupted(Z, 1.0);                      
722       return;                                     
723     }                                             
724     std::ostringstream ss;                        
725     ss << "mupair/" << particle->GetParticleNa    
726     std::ofstream outfile(ss.str());              
727     pv->Store(outfile);                           
728   }                                               
729 }                                                 
730                                                   
731 //....oooOO0OOooo........oooOO0OOooo........oo    
732                                                   
733 G4bool G4MuPairProductionModel::RetrieveTables    
734 {                                                 
735   for (G4int iz=0; iz<NZDATPAIR; ++iz) {          
736     G4double Z = ZDATPAIR[iz];                    
737     G4Physics2DVector* pv = new G4Physics2DVec    
738     std::ostringstream ss;                        
739     ss << G4EmParameters::Instance()->GetDirLE    
740        << particle->GetParticleName() << Z <<     
741     std::ifstream infile(ss.str(), std::ios::i    
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........oo    
752