Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/processes/electromagnetic/highenergy/src/G4GammaConversionToMuons.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 //         ------------ G4GammaConversionToMuons physics process ------
 28 //         by H.Burkhardt, S. Kelner and R. Kokoulin, April 2002
 29 //
 30 //
 31 // 07-08-02: missprint in OR condition in DoIt : f1<0 || f1>f1_max ..etc ...
 32 // 25-10-04: migrade to new interfaces of ParticleChange (vi)
 33 // ---------------------------------------------------------------------------
 34 
 35 #include "G4GammaConversionToMuons.hh"
 36 
 37 #include "G4BetheHeitler5DModel.hh"
 38 #include "G4Electron.hh"
 39 #include "G4EmParameters.hh"
 40 #include "G4EmProcessSubType.hh"
 41 #include "G4Exp.hh"
 42 #include "G4Gamma.hh"
 43 #include "G4Log.hh"
 44 #include "G4LossTableManager.hh"
 45 #include "G4MuonMinus.hh"
 46 #include "G4MuonPlus.hh"
 47 #include "G4NistManager.hh"
 48 #include "G4PhysicalConstants.hh"
 49 #include "G4Positron.hh"
 50 #include "G4ProductionCutsTable.hh"
 51 #include "G4SystemOfUnits.hh"
 52 #include "G4UnitsTable.hh"
 53 
 54 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
 55 
 56 static const G4double sqrte = std::sqrt(std::exp(1.));
 57 static const G4double PowSat = -0.88;
 58 
 59 G4GammaConversionToMuons::G4GammaConversionToMuons(const G4String& processName,
 60                G4ProcessType type)
 61   : G4VDiscreteProcess (processName, type),
 62     Mmuon(G4MuonPlus::MuonPlus()->GetPDGMass()),
 63     Rc(CLHEP::elm_coupling / Mmuon),
 64     LimitEnergy(5. * Mmuon),
 65     LowestEnergyLimit(2. * Mmuon),
 66     HighestEnergyLimit(1e12 * CLHEP::GeV),  // ok to 1e12GeV, then LPM suppression
 67     theGamma(G4Gamma::Gamma()),
 68     theMuonPlus(G4MuonPlus::MuonPlus()),
 69     theMuonMinus(G4MuonMinus::MuonMinus())
 70 {
 71   SetProcessSubType(fGammaConversionToMuMu);
 72   fManager = G4LossTableManager::Instance();
 73   fManager->Register(this);
 74 }
 75 
 76 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
 77 
 78 G4GammaConversionToMuons::~G4GammaConversionToMuons()
 79 {
 80   fManager->DeRegister(this);
 81 }
 82 
 83 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
 84 
 85 G4bool G4GammaConversionToMuons::IsApplicable(const G4ParticleDefinition& part)
 86 {
 87   return (&part == theGamma);
 88 }
 89 
 90 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 91 
 92 void G4GammaConversionToMuons::BuildPhysicsTable(const G4ParticleDefinition& p)
 93 {
 94   Energy5DLimit = G4EmParameters::Instance()->MaxEnergyFor5DMuPair();
 95 
 96   auto table = G4Material::GetMaterialTable();
 97   std::size_t nelm = 0;
 98   for (auto const& mat : *table) {
 99     std::size_t n = mat->GetNumberOfElements();
100     nelm = std::max(nelm, n);
101   }
102   temp.resize(nelm, 0);
103 
104   if (Energy5DLimit > 0.0 && nullptr != f5Dmodel) {
105     f5Dmodel = new G4BetheHeitler5DModel();
106     f5Dmodel->SetLeptonPair(theMuonPlus, theMuonMinus);
107     const std::size_t numElems = G4ProductionCutsTable::GetProductionCutsTable()->GetTableSize();
108     const G4DataVector cuts(numElems);
109     f5Dmodel->Initialise(&p, cuts);
110   }
111   PrintInfoDefinition();
112 }
113 
114 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
115 
116 G4double G4GammaConversionToMuons::GetMeanFreePath(const G4Track& aTrack, G4double,
117                                                    G4ForceCondition*)
118 // returns the photon mean free path in GEANT4 internal units
119 {
120   const G4DynamicParticle* aDynamicGamma = aTrack.GetDynamicParticle();
121   G4double GammaEnergy = aDynamicGamma->GetKineticEnergy();
122   const G4Material* aMaterial = aTrack.GetMaterial();
123   return ComputeMeanFreePath(GammaEnergy, aMaterial);
124 }
125 
126 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
127 
128 G4double 
129 G4GammaConversionToMuons::ComputeMeanFreePath(G4double GammaEnergy,
130                                               const G4Material* aMaterial)
131 
132 // computes and returns the photon mean free path in GEANT4 internal units
133 {
134   if(GammaEnergy <= LowestEnergyLimit) { return DBL_MAX; }
135   const G4ElementVector* theElementVector = aMaterial->GetElementVector();
136   const G4double* NbOfAtomsPerVolume = aMaterial->GetVecNbOfAtomsPerVolume();
137 
138   G4double SIGMA = 0.0;
139   G4double fact  = 1.0;
140   G4double e = GammaEnergy;
141   // low energy approximation as in Bethe-Heitler model
142   if(e < LimitEnergy) {
143     G4double y = (e - LowestEnergyLimit)/(LimitEnergy - LowestEnergyLimit);
144     fact = y*y;
145     e = LimitEnergy;
146   } 
147 
148   for ( std::size_t i=0 ; i < aMaterial->GetNumberOfElements(); ++i)
149   {
150     SIGMA += NbOfAtomsPerVolume[i] * fact *
151       ComputeCrossSectionPerAtom(e, (*theElementVector)[i]->GetZasInt());
152   }
153   return (SIGMA > 0.0) ? 1./SIGMA : DBL_MAX;
154 }
155 
156 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
157 
158 G4double G4GammaConversionToMuons::GetCrossSectionPerAtom(
159                                    const G4DynamicParticle* aDynamicGamma,
160                                    const G4Element* anElement)
161 
162 // gives the total cross section per atom in GEANT4 internal units
163 {
164    return ComputeCrossSectionPerAtom(aDynamicGamma->GetKineticEnergy(),
165                                      anElement->GetZasInt());
166 }
167 
168 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
169 
170 G4double G4GammaConversionToMuons::ComputeCrossSectionPerAtom(
171                          G4double Egam, G4int Z)
172        
173 // Calculates the microscopic cross section in GEANT4 internal units.
174 // Total cross section parametrisation from H.Burkhardt
175 // It gives a good description at any energy (from 0 to 10**21 eV)
176 { 
177   if(Egam <= LowestEnergyLimit) { return 0.0; }
178 
179   G4NistManager* nist = G4NistManager::Instance();
180 
181   G4double PowThres, Ecor, B, Dn, Zthird, Winfty, WMedAppr, Wsatur, sigfac;
182 
183   if (Z == 1) {  // special case of Hydrogen
184     B = 202.4;
185     Dn = 1.49;
186   }
187   else {
188     B = 183.;
189     Dn = 1.54 * nist->GetA27(Z);
190   }
191   Zthird = 1. / nist->GetZ13(Z);  // Z**(-1/3)
192   Winfty = B * Zthird * Mmuon / (Dn * electron_mass_c2);
193   WMedAppr = 1. / (4. * Dn * sqrte * Mmuon);
194   Wsatur = Winfty / WMedAppr;
195   sigfac = 4. * fine_structure_const * Z * Z * Rc * Rc;
196   PowThres = 1.479 + 0.00799 * Dn;
197   Ecor = -18. + 4347. / (B * Zthird);
198 
199   G4double CorFuc = 1. + .04 * G4Log(1. + Ecor / Egam);
200   G4double Eg =
201     G4Exp(G4Log(1. - 4. * Mmuon / Egam) * PowThres)
202     * G4Exp(G4Log(G4Exp(G4Log(Wsatur) * PowSat) + G4Exp(G4Log(Egam) * PowSat)) / PowSat);
203 
204   G4double CrossSection = 7. / 9. * sigfac * G4Log(1. + WMedAppr * CorFuc * Eg);
205   CrossSection *= CrossSecFactor;  // increase the CrossSection by  (by default 1)
206   return CrossSection;
207 }
208 
209 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
210 
211 void G4GammaConversionToMuons::SetCrossSecFactor(G4double fac)
212 // Set the factor to artificially increase the cross section
213 {
214   if (fac < 0.0) return;
215   CrossSecFactor = fac;
216   if (verboseLevel > 1) {
217     G4cout << "The cross section for GammaConversionToMuons is artificially "
218            << "increased by the CrossSecFactor=" << CrossSecFactor << G4endl;
219   }
220 }
221 
222 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
223 
224 G4VParticleChange* G4GammaConversionToMuons::PostStepDoIt(
225                                                         const G4Track& aTrack,
226                                                         const G4Step&  aStep)
227 //
228 // generation of gamma->mu+mu-
229 //
230 {
231   aParticleChange.Initialize(aTrack);
232   const G4Material* aMaterial = aTrack.GetMaterial();
233 
234   // current Gamma energy and direction, return if energy too low
235   const G4DynamicParticle* aDynamicGamma = aTrack.GetDynamicParticle();
236   G4double Egam = aDynamicGamma->GetKineticEnergy();
237   if (Egam <= LowestEnergyLimit) {
238     return G4VDiscreteProcess::PostStepDoIt(aTrack,aStep);
239   }
240   //
241   // Kill the incident photon
242   //
243   aParticleChange.ProposeMomentumDirection( 0., 0., 0. ) ;
244   aParticleChange.ProposeEnergy( 0. ) ;
245   aParticleChange.ProposeTrackStatus( fStopAndKill ) ;
246 
247   if (Egam <= Energy5DLimit) {
248     std::vector<G4DynamicParticle*> fvect;
249     f5Dmodel->SampleSecondaries(&fvect, aTrack.GetMaterialCutsCouple(), 
250         aTrack.GetDynamicParticle(), 0.0, DBL_MAX);
251     for(auto dp : fvect) { aParticleChange.AddSecondary(dp); }
252     return G4VDiscreteProcess::PostStepDoIt(aTrack,aStep);
253   }  
254 
255   G4ParticleMomentum GammaDirection = aDynamicGamma->GetMomentumDirection();
256 
257   // select randomly one element constituting the material
258   const G4Element* anElement = SelectRandomAtom(aDynamicGamma, aMaterial);
259   G4int Z = anElement->GetZasInt();
260   G4NistManager* nist = G4NistManager::Instance();
261 
262   G4double B, Dn;
263   G4double A027 = nist->GetA27(Z);
264 
265   if (Z == 1) {  // special case of Hydrogen
266     B = 202.4;
267     Dn = 1.49;
268   }
269   else {
270     B = 183.;
271     Dn = 1.54 * A027;
272   }
273   G4double Zthird = 1. / nist->GetZ13(Z);  // Z**(-1/3)
274   G4double Winfty = B * Zthird * Mmuon / (Dn * electron_mass_c2);
275 
276   G4double C1Num = 0.138 * A027;
277   G4double C1Num2 = C1Num * C1Num;
278   G4double C2Term2 = electron_mass_c2 / (183. * Zthird * Mmuon);
279 
280   G4double GammaMuonInv = Mmuon / Egam;
281 
282   // generate xPlus according to the differential cross section by rejection
283   G4double xmin = (Egam <= LimitEnergy) ? 0.5 : 0.5 - std::sqrt(0.25 - GammaMuonInv);
284   G4double xmax = 1. - xmin;
285 
286   G4double Ds2 = (Dn * sqrte - 2.);
287   G4double sBZ = sqrte * B * Zthird / electron_mass_c2;
288   G4double LogWmaxInv =
289     1. / G4Log(Winfty * (1. + 2. * Ds2 * GammaMuonInv) / (1. + 2. * sBZ * Mmuon * GammaMuonInv));
290   G4double xPlus = 0.5;
291   G4double xMinus = 0.5;
292   G4double xPM = 0.25;
293 
294   G4int nn = 0;
295   const G4int nmax = 1000;
296 
297   // sampling for Egam > LimitEnergy
298   if (xmin < 0.5) {
299     G4double result, W;
300     do {
301       xPlus = xmin + G4UniformRand() * (xmax - xmin);
302       xMinus = 1. - xPlus;
303       xPM = xPlus * xMinus;
304       G4double del = Mmuon * Mmuon / (2. * Egam * xPM);
305       W = Winfty * (1. + Ds2 * del / Mmuon) / (1. + sBZ * del);
306       G4double xxp = 1. - 4. / 3. * xPM;  // the main xPlus dependence
307       result = (xxp > 0.) ? xxp * G4Log(W) * LogWmaxInv : 0.0;
308       if (result > 1.) {
309         G4cout << "G4GammaConversionToMuons::PostStepDoIt WARNING:"
310                << " in dSigxPlusGen, result=" << result << " > 1" << G4endl;
311       }
312       ++nn;
313       if(nn >= nmax) { break; }
314     }
315     // Loop checking, 07-Aug-2015, Vladimir Ivanchenko
316     while (G4UniformRand() > result);
317   }
318 
319   // now generate the angular variables via the auxilary variables t,psi,rho
320   G4double t;
321   G4double psi;
322   G4double rho;
323 
324   G4double a3 = (GammaMuonInv / (2. * xPM));
325   G4double a33 = a3 * a3;
326   G4double f1;
327   G4double b1  = 1./(4.*C1Num2);
328   G4double b3  = b1*b1*b1;
329   G4double a21 = a33 + b1;
330   
331   G4double f1_max=-(1.-xPM)*(2.*b1+(a21+a33)*G4Log(a33/a21))/(2*b3);  
332 
333   G4double thetaPlus,thetaMinus,phiHalf; // final angular variables
334   nn = 0;
335   // t, psi, rho generation start  (while angle < pi)
336   do {
337     //generate t by the rejection method
338     do { 
339       ++nn;
340       t=G4UniformRand();
341       G4double a34=a33/(t*t);
342       G4double a22 = a34 + b1;
343       if(std::abs(b1)<0.0001*a34) {
344         // special case of a34=a22 because of logarithm accuracy
345         f1=(1.-2.*xPM+4.*xPM*t*(1.-t))/(12.*a34*a34*a34*a34);
346       }
347       else {
348         f1=-(1.-2.*xPM+4.*xPM*t*(1.-t))*(2.*b1+(a22+a34)*G4Log(a34/a22))/(2*b3);
349       }
350       if (f1 < 0.0 || f1 > f1_max) {  // should never happend
351         G4cout << "G4GammaConversionToMuons::PostStepDoIt WARNING:"
352         << "outside allowed range f1=" << f1
353         << " is set to zero, a34 = "<< a34 << " a22 = "<<a22<<"."
354         << G4endl;
355         f1 = 0.0;
356       }
357       if(nn > nmax) { break; }
358       // Loop checking, 07-Aug-2015, Vladimir Ivanchenko  
359     } while ( G4UniformRand()*f1_max > f1);
360     // generate psi by the rejection method
361     G4double f2_max=1.-2.*xPM*(1.-4.*t*(1.-t));
362     // long version
363     G4double f2;
364     do { 
365       ++nn;
366       psi=twopi*G4UniformRand();
367       f2=1.-2.*xPM+4.*xPM*t*(1.-t)*(1.+std::cos(2.*psi));
368       if(f2<0 || f2> f2_max) { // should never happend
369         G4cout << "G4GammaConversionToMuons::PostStepDoIt WARNING:"
370                << "outside allowed range f2=" << f2 << " is set to zero" << G4endl;
371         f2 = 0.0;
372       }
373       if(nn >= nmax) { break; }
374       // Loop checking, 07-Aug-2015, Vladimir Ivanchenko
375     } while ( G4UniformRand()*f2_max > f2);
376 
377     // generate rho by direct transformation
378     G4double C2Term1=GammaMuonInv/(2.*xPM*t);
379     G4double C22 = C2Term1*C2Term1+C2Term2*C2Term2;
380     G4double C2=4.*C22*C22/std::sqrt(xPM);
381     G4double rhomax=(1./t-1.)*1.9/A027;
382     G4double beta=G4Log( (C2+rhomax*rhomax*rhomax*rhomax)/C2 );
383     rho=G4Exp(G4Log(C2 *( G4Exp(beta*G4UniformRand())-1. ))*0.25);
384 
385     //now get from t and psi the kinematical variables
386     G4double u=std::sqrt(1./t-1.);
387     G4double xiHalf=0.5*rho*std::cos(psi);
388     phiHalf=0.5*rho/u*std::sin(psi);
389 
390     thetaPlus =GammaMuonInv*(u+xiHalf)/xPlus;
391     thetaMinus=GammaMuonInv*(u-xiHalf)/xMinus;
392 
393     // protection against infinite loop
394     if(nn > nmax) {
395       if(std::abs(thetaPlus)>pi) { thetaPlus = 0.0; }
396       if(std::abs(thetaMinus)>pi) { thetaMinus = 0.0; }
397     }
398 
399     // Loop checking, 07-Aug-2015, Vladimir Ivanchenko
400   } while ( std::abs(thetaPlus)>pi || std::abs(thetaMinus) >pi);
401 
402   // now construct the vectors
403   // azimuthal symmetry, take phi0 at random between 0 and 2 pi
404   G4double phi0=twopi*G4UniformRand(); 
405   G4double EPlus=xPlus*Egam;
406   G4double EMinus=xMinus*Egam;
407 
408   // mu+ mu- directions for gamma in z-direction
409   G4ThreeVector MuPlusDirection  ( std::sin(thetaPlus) *std::cos(phi0+phiHalf),
410                    std::sin(thetaPlus)  *std::sin(phi0+phiHalf), std::cos(thetaPlus) );
411   G4ThreeVector MuMinusDirection (-std::sin(thetaMinus)*std::cos(phi0-phiHalf),
412                   -std::sin(thetaMinus) *std::sin(phi0-phiHalf), std::cos(thetaMinus) );
413   // rotate to actual gamma direction
414   MuPlusDirection.rotateUz(GammaDirection);
415   MuMinusDirection.rotateUz(GammaDirection);
416 
417   // create G4DynamicParticle object for the particle1
418   auto aParticle1 = new G4DynamicParticle(theMuonPlus, MuPlusDirection, EPlus - Mmuon);
419   aParticleChange.AddSecondary(aParticle1);
420   // create G4DynamicParticle object for the particle2
421   auto aParticle2 = new G4DynamicParticle(theMuonMinus, MuMinusDirection, EMinus - Mmuon);
422   aParticleChange.AddSecondary(aParticle2);
423   //  Reset NbOfInteractionLengthLeft and return aParticleChange
424   return G4VDiscreteProcess::PostStepDoIt( aTrack, aStep );
425 }
426 
427 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
428 
429 const G4Element* G4GammaConversionToMuons::SelectRandomAtom(
430       const G4DynamicParticle* aDynamicGamma,
431       const G4Material* aMaterial)
432 {
433   // select randomly 1 element within the material, invoked by PostStepDoIt
434 
435   const std::size_t NumberOfElements      = aMaterial->GetNumberOfElements();
436   const G4ElementVector* theElementVector = aMaterial->GetElementVector();
437   const G4Element* elm = (*theElementVector)[0];
438 
439   if (NumberOfElements > 1) {
440     G4double e = std::max(aDynamicGamma->GetKineticEnergy(), LimitEnergy);
441     const G4double* natom = aMaterial->GetVecNbOfAtomsPerVolume();
442 
443     G4double sum = 0.;
444     for (std::size_t i=0; i<NumberOfElements; ++i) {
445       elm = (*theElementVector)[i];
446       sum += natom[i]*ComputeCrossSectionPerAtom(e, elm->GetZasInt());
447       temp[i] = sum;
448     }
449     sum *= G4UniformRand();
450     for (std::size_t i=0; i<NumberOfElements; ++i) {
451       if(sum <= temp[i]) {
452         elm = (*theElementVector)[i];
453         break;
454       }
455     }
456   }
457   return elm;
458 }
459 
460 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
461 
462 void G4GammaConversionToMuons::PrintInfoDefinition()
463 {
464   G4String comments = "gamma->mu+mu- Bethe Heitler process, SubType= ";
465   G4cout << G4endl << GetProcessName() << ":  " << comments << GetProcessSubType() << G4endl;
466   G4cout << "        good cross section parametrization from "
467          << G4BestUnit(LowestEnergyLimit, "Energy") << " to " << HighestEnergyLimit / GeV
468          << " GeV for all Z." << G4endl;
469   G4cout << "        cross section factor: " << CrossSecFactor << G4endl;
470 }
471 
472 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
473