Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/processes/hadronic/models/coherent_elastic/include/G4DiffuseElastic.hh

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 // Author: V. Grichine (Vladimir,Grichine@cern.ch)
 29 //
 30 //
 31 // G4 Model: diffuse optical elastic scattering with 4-momentum balance 
 32 //
 33 // Class Description
 34 // Final state production model for hadron nuclear elastic scattering; 
 35 // Class Description - End
 36 //
 37 //
 38 // 24.05.07 V. Grichine, first implementation for hadron (no Coulomb) elastic scattering
 39 // 04.09.07 V. Grichine, implementation for Coulomb elastic scattering
 40 // 12.06.11 V. Grichine, new interface to G4hadronElastic
 41 
 42 #ifndef G4DiffuseElastic_h
 43 #define G4DiffuseElastic_h 1
 44 
 45 #include <CLHEP/Units/PhysicalConstants.h>
 46 #include "globals.hh"
 47 #include "G4HadronElastic.hh"
 48 #include "G4HadProjectile.hh"
 49 #include "G4Nucleus.hh"
 50 
 51 #include "G4Pow.hh"
 52 
 53 
 54 class G4ParticleDefinition;
 55 class G4PhysicsTable;
 56 class G4PhysicsLogVector;
 57 
 58 class G4DiffuseElastic : public G4HadronElastic // G4HadronicInteraction
 59 {
 60 public:
 61 
 62   G4DiffuseElastic();
 63 
 64   // G4DiffuseElastic(const G4ParticleDefinition* aParticle);
 65 
 66 
 67 
 68 
 69 
 70   virtual ~G4DiffuseElastic();
 71 
 72   virtual G4bool IsApplicable(const G4HadProjectile &/*aTrack*/, 
 73             G4Nucleus & /*targetNucleus*/);
 74 
 75   void Initialise();
 76 
 77   void InitialiseOnFly(G4double Z, G4double A);
 78 
 79   void BuildAngleTable();
 80 
 81  
 82   // G4HadFinalState* ApplyYourself(const G4HadProjectile & aTrack, G4Nucleus & targetNucleus);
 83 
 84   virtual G4double SampleInvariantT(const G4ParticleDefinition* p, 
 85             G4double plab,
 86             G4int Z, G4int A);
 87 
 88   G4double NeutronTuniform(G4int Z);
 89 
 90   void SetPlabLowLimit(G4double value);
 91 
 92   void SetHEModelLowLimit(G4double value);
 93 
 94   void SetQModelLowLimit(G4double value);
 95 
 96   void SetLowestEnergyLimit(G4double value);
 97 
 98   void SetRecoilKinEnergyLimit(G4double value);
 99 
100   G4double SampleT(const G4ParticleDefinition* aParticle, 
101                          G4double p, G4double A);
102 
103   G4double SampleTableT(const G4ParticleDefinition* aParticle, 
104                          G4double p, G4double Z, G4double A);
105 
106   G4double SampleThetaCMS(const G4ParticleDefinition* aParticle, G4double p, G4double A);
107 
108   G4double SampleTableThetaCMS(const G4ParticleDefinition* aParticle, G4double p, 
109                                      G4double Z, G4double A);
110 
111   G4double GetScatteringAngle(G4int iMomentum, G4int iAngle, G4double position);
112 
113   G4double SampleThetaLab(const G4HadProjectile* aParticle, 
114                                 G4double tmass, G4double A);
115 
116   G4double GetDiffuseElasticXsc( const G4ParticleDefinition* particle, 
117                                  G4double theta, 
118                G4double momentum, 
119          G4double A         );
120 
121   G4double GetInvElasticXsc( const G4ParticleDefinition* particle, 
122                                  G4double theta, 
123                G4double momentum, 
124          G4double A, G4double Z );
125 
126   G4double GetDiffuseElasticSumXsc( const G4ParticleDefinition* particle, 
127                                  G4double theta, 
128                G4double momentum, 
129          G4double A, G4double Z );
130 
131   G4double GetInvElasticSumXsc( const G4ParticleDefinition* particle, 
132                                  G4double tMand, 
133                G4double momentum, 
134          G4double A, G4double Z );
135 
136   G4double IntegralElasticProb( const G4ParticleDefinition* particle, 
137                                  G4double theta, 
138                G4double momentum, 
139          G4double A            );
140   
141 
142   G4double GetCoulombElasticXsc( const G4ParticleDefinition* particle, 
143                                  G4double theta, 
144                G4double momentum, 
145          G4double Z         );
146 
147   G4double GetInvCoulombElasticXsc( const G4ParticleDefinition* particle, 
148                                  G4double tMand, 
149                G4double momentum, 
150          G4double A, G4double Z         );
151 
152   G4double GetCoulombTotalXsc( const G4ParticleDefinition* particle,  
153                G4double momentum, G4double Z       );
154 
155   G4double GetCoulombIntegralXsc( const G4ParticleDefinition* particle,  
156                G4double momentum, G4double Z, 
157                                  G4double theta1, G4double theta2         );
158 
159 
160   G4double CalculateParticleBeta( const G4ParticleDefinition* particle, 
161                                   G4double momentum    );
162 
163   G4double CalculateZommerfeld( G4double beta, G4double Z1, G4double Z2 );
164 
165   G4double CalculateAm( G4double momentum, G4double n, G4double Z);
166 
167   G4double CalculateNuclearRad( G4double A);
168 
169   G4double ThetaCMStoThetaLab(const G4DynamicParticle* aParticle, 
170                                 G4double tmass, G4double thetaCMS);
171 
172   G4double ThetaLabToThetaCMS(const G4DynamicParticle* aParticle, 
173                                 G4double tmass, G4double thetaLab);
174 
175   void TestAngleTable(const G4ParticleDefinition* theParticle, G4double partMom,
176           G4double Z, G4double A);
177 
178 
179 
180   G4double BesselJzero(G4double z);
181   G4double BesselJone(G4double z);
182   G4double DampFactor(G4double z);
183   G4double BesselOneByArg(G4double z);
184 
185   G4double GetDiffElasticProb(G4double theta);
186   G4double GetDiffElasticSumProb(G4double theta);
187   G4double GetDiffElasticSumProbA(G4double alpha);
188   G4double GetIntegrandFunction(G4double theta);
189 
190 
191   G4double GetNuclearRadius(){return fNuclearRadius;};
192 
193 private:
194 
195 
196   G4ParticleDefinition* theProton;
197   G4ParticleDefinition* theNeutron;
198   G4ParticleDefinition* theDeuteron;
199   G4ParticleDefinition* theAlpha;
200 
201   const G4ParticleDefinition* thePionPlus;
202   const G4ParticleDefinition* thePionMinus;
203 
204   G4double lowEnergyRecoilLimit;  
205   G4double lowEnergyLimitHE;  
206   G4double lowEnergyLimitQ;  
207   G4double lowestEnergyLimit;  
208   G4double plabLowLimit;
209 
210   G4int fEnergyBin;
211   G4int fAngleBin;
212 
213   G4PhysicsLogVector*           fEnergyVector;
214   G4PhysicsTable*               fAngleTable;
215   std::vector<G4PhysicsTable*>  fAngleBank;
216 
217   std::vector<G4double> fElementNumberVector;
218   std::vector<G4String> fElementNameVector;
219 
220   const G4ParticleDefinition* fParticle;
221   G4double fWaveVector;
222   G4double fAtomicWeight;
223   G4double fAtomicNumber;
224   G4double fNuclearRadius;
225   G4double fBeta;
226   G4double fZommerfeld;
227   G4double fAm;
228   G4bool fAddCoulomb;
229 
230 };
231 
232 inline G4bool G4DiffuseElastic::IsApplicable(const G4HadProjectile & projectile, 
233             G4Nucleus & nucleus)
234 {
235   if( ( projectile.GetDefinition() == G4Proton::Proton() ||
236         projectile.GetDefinition() == G4Neutron::Neutron() ||
237         projectile.GetDefinition() == G4PionPlus::PionPlus() ||
238         projectile.GetDefinition() == G4PionMinus::PionMinus() ||
239         projectile.GetDefinition() == G4KaonPlus::KaonPlus() ||
240         projectile.GetDefinition() == G4KaonMinus::KaonMinus() ) &&
241 
242         nucleus.GetZ_asInt() >= 2 ) return true;
243   else                              return false;
244 }
245 
246 inline void G4DiffuseElastic::SetRecoilKinEnergyLimit(G4double value)
247 {
248   lowEnergyRecoilLimit = value;
249 }
250 
251 inline void G4DiffuseElastic::SetPlabLowLimit(G4double value)
252 {
253   plabLowLimit = value;
254 }
255 
256 inline void G4DiffuseElastic::SetHEModelLowLimit(G4double value)
257 {
258   lowEnergyLimitHE = value;
259 }
260 
261 inline void G4DiffuseElastic::SetQModelLowLimit(G4double value)
262 {
263   lowEnergyLimitQ = value;
264 }
265 
266 inline void G4DiffuseElastic::SetLowestEnergyLimit(G4double value)
267 {
268   lowestEnergyLimit = value;
269 }
270 
271 
272 /////////////////////////////////////////////////////////////
273 //
274 // Bessel J0 function based on rational approximation from 
275 // J.F. Hart, Computer Approximations, New York, Willey 1968, p. 141 
276 
277 inline G4double G4DiffuseElastic::BesselJzero(G4double value)
278 {
279   G4double modvalue, value2, fact1, fact2, arg, shift, bessel;
280 
281   modvalue = std::fabs(value);
282 
283   if ( value < 8.0 && value > -8.0 )
284   {
285     value2 = value*value;
286 
287     fact1  = 57568490574.0 + value2*(-13362590354.0 
288                            + value2*( 651619640.7 
289                            + value2*(-11214424.18 
290                            + value2*( 77392.33017 
291                            + value2*(-184.9052456   ) ) ) ) );
292 
293     fact2  = 57568490411.0 + value2*( 1029532985.0 
294                            + value2*( 9494680.718
295                            + value2*(59272.64853
296                            + value2*(267.8532712 
297                            + value2*1.0               ) ) ) );
298 
299     bessel = fact1/fact2;
300   } 
301   else 
302   {
303     arg    = 8.0/modvalue;
304 
305     value2 = arg*arg;
306 
307     shift  = modvalue-0.785398164;
308 
309     fact1  = 1.0 + value2*(-0.1098628627e-2 
310                  + value2*(0.2734510407e-4
311                  + value2*(-0.2073370639e-5 
312                  + value2*0.2093887211e-6    ) ) );
313 
314     fact2  = -0.1562499995e-1 + value2*(0.1430488765e-3
315                               + value2*(-0.6911147651e-5 
316                               + value2*(0.7621095161e-6
317                               - value2*0.934945152e-7    ) ) );
318 
319     bessel = std::sqrt(0.636619772/modvalue)*(std::cos(shift)*fact1 - arg*std::sin(shift)*fact2 );
320   }
321   return bessel;
322 }
323 
324 /////////////////////////////////////////////////////////////
325 //
326 // Bessel J1 function based on rational approximation from 
327 // J.F. Hart, Computer Approximations, New York, Willey 1968, p. 141 
328 
329 inline G4double G4DiffuseElastic::BesselJone(G4double value)
330 {
331   G4double modvalue, value2, fact1, fact2, arg, shift, bessel;
332 
333   modvalue = std::fabs(value);
334 
335   if ( modvalue < 8.0 ) 
336   {
337     value2 = value*value;
338 
339     fact1  = value*(72362614232.0 + value2*(-7895059235.0 
340                                   + value2*( 242396853.1
341                                   + value2*(-2972611.439 
342                                   + value2*( 15704.48260 
343                                   + value2*(-30.16036606  ) ) ) ) ) );
344 
345     fact2  = 144725228442.0 + value2*(2300535178.0 
346                             + value2*(18583304.74
347                             + value2*(99447.43394 
348                             + value2*(376.9991397 
349                             + value2*1.0             ) ) ) );
350     bessel = fact1/fact2;
351   } 
352   else 
353   {
354     arg    = 8.0/modvalue;
355 
356     value2 = arg*arg;
357 
358     shift  = modvalue - 2.356194491;
359 
360     fact1  = 1.0 + value2*( 0.183105e-2 
361                  + value2*(-0.3516396496e-4
362                  + value2*(0.2457520174e-5 
363                  + value2*(-0.240337019e-6          ) ) ) );
364 
365     fact2 = 0.04687499995 + value2*(-0.2002690873e-3
366                           + value2*( 0.8449199096e-5
367                           + value2*(-0.88228987e-6
368                           + value2*0.105787412e-6       ) ) );
369 
370     bessel = std::sqrt( 0.636619772/modvalue)*(std::cos(shift)*fact1 - arg*std::sin(shift)*fact2);
371 
372     if (value < 0.0) bessel = -bessel;
373   }
374   return bessel;
375 }
376 
377 ////////////////////////////////////////////////////////////////////
378 //
379 // damp factor in diffraction x/sh(x), x was already *pi
380 
381 inline G4double G4DiffuseElastic::DampFactor(G4double x)
382 {
383   G4double df;
384   G4double f2 = 2., f3 = 6., f4 = 24.; // first factorials
385 
386   // x *= pi;
387 
388   if( std::fabs(x) < 0.01 )
389   { 
390     df = 1./(1. + x/f2 + x*x/f3 + x*x*x/f4);
391   }
392   else
393   {
394     df = x/std::sinh(x); 
395   }
396   return df;
397 }
398 
399 
400 ////////////////////////////////////////////////////////////////////
401 //
402 // return J1(x)/x with special case for small x
403 
404 inline G4double G4DiffuseElastic::BesselOneByArg(G4double x)
405 {
406   G4double x2, result;
407   
408   if( std::fabs(x) < 0.01 )
409   { 
410    x     *= 0.5;
411    x2     = x*x;
412    result = 2. - x2 + x2*x2/6.;
413   }
414   else
415   {
416     result = BesselJone(x)/x; 
417   }
418   return result;
419 }
420 
421 ////////////////////////////////////////////////////////////////////
422 //
423 // return particle beta
424 
425 inline  G4double G4DiffuseElastic::CalculateParticleBeta( const G4ParticleDefinition* particle, 
426                                   G4double momentum    )
427 {
428   G4double mass = particle->GetPDGMass();
429   G4double a    = momentum/mass;
430   fBeta         = a/std::sqrt(1+a*a);
431 
432   return fBeta; 
433 }
434 
435 ////////////////////////////////////////////////////////////////////
436 //
437 // return Zommerfeld parameter for Coulomb scattering
438 
439 inline  G4double G4DiffuseElastic::CalculateZommerfeld( G4double beta, G4double Z1, G4double Z2 )
440 {
441   fZommerfeld = CLHEP::fine_structure_const*Z1*Z2/beta;
442 
443   return fZommerfeld; 
444 }
445 
446 ////////////////////////////////////////////////////////////////////
447 //
448 // return Wentzel correction for Coulomb scattering
449 
450 inline  G4double G4DiffuseElastic::CalculateAm( G4double momentum, G4double n, G4double Z)
451 {
452   G4double k   = momentum/CLHEP::hbarc;
453   G4double ch  = 1.13 + 3.76*n*n;
454   G4double zn  = 1.77*k*(1.0/G4Pow::GetInstance()->A13(Z))*CLHEP::Bohr_radius;
455   G4double zn2 = zn*zn;
456   fAm          = ch/zn2;
457 
458   return fAm;
459 }
460 
461 ////////////////////////////////////////////////////////////////////
462 //
463 // calculate nuclear radius for different atomic weights using different approximations
464 
465 inline  G4double G4DiffuseElastic::CalculateNuclearRad( G4double A)
466 {
467   G4double R, r0, a11, a12, a13, a2, a3;
468 
469   a11 = 1.26;  // 1.08, 1.16
470   a12 = 1.;  // 1.08, 1.16
471   a13 = 1.12;  // 1.08, 1.16
472   a2 = 1.1;
473   a3 = 1.;
474 
475   // Special rms radii for light nucleii
476 
477   if (A < 50.)
478   {
479     if     (std::abs(A-1.) < 0.5)                         return 0.89*CLHEP::fermi; // p
480     else if(std::abs(A-2.) < 0.5)                         return 2.13*CLHEP::fermi; // d
481     else if(  // std::abs(Z-1.) < 0.5 && 
482 std::abs(A-3.) < 0.5) return 1.80*CLHEP::fermi; // t
483 
484     // else if(std::abs(Z-2.) < 0.5 && std::abs(A-3.) < 0.5) return 1.96CLHEP::fermi; // He3
485     else if( // std::abs(Z-2.) < 0.5 && 
486 std::abs(A-4.) < 0.5) return 1.68*CLHEP::fermi; // He4
487 
488     else if(  // std::abs(Z-3.) < 0.5
489         std::abs(A-7.) < 0.5   )                         return 2.40*CLHEP::fermi; // Li7
490     else if(  // std::abs(Z-4.) < 0.5 
491 std::abs(A-9.) < 0.5)                         return 2.51*CLHEP::fermi; // Be9
492 
493     else if( 10.  < A && A <= 16. ) r0  = a11*( 1 - (1.0/G4Pow::GetInstance()->A23(A)) )*CLHEP::fermi;   // 1.08CLHEP::fermi;
494     else if( 15.  < A && A <= 20. ) r0  = a12*( 1 - (1.0/G4Pow::GetInstance()->A23(A)) )*CLHEP::fermi;
495     else if( 20.  < A && A <= 30. ) r0  = a13*( 1 - (1.0/G4Pow::GetInstance()->A23(A)) )*CLHEP::fermi;
496     else                            r0  = a2*CLHEP::fermi;
497 
498     R = r0*G4Pow::GetInstance()->A13(A);
499   }
500   else
501   {
502     r0 = a3*CLHEP::fermi;
503 
504     R  = r0*G4Pow::GetInstance()->powA(A, 0.27);
505   }
506   fNuclearRadius = R;
507   return R;
508   /*
509   G4double r0;
510   if( A < 50. )
511   {
512     if( A > 10. ) r0  = 1.16*( 1 - (1.0/G4Pow::GetInstance()->A23(A)) )*CLHEP::fermi;   // 1.08CLHEP::fermi;
513     else          r0  = 1.1*CLHEP::fermi;
514     fNuclearRadius = r0*G4Pow::GetInstance()->A13(A);
515   }
516   else
517   {
518     r0 = 1.7*CLHEP::fermi;   // 1.7*CLHEP::fermi;
519     fNuclearRadius = r0*G4Pow::GetInstance()->powA(A, 0.27); // 0.27);
520   }
521   return fNuclearRadius;
522   */
523 }
524 
525 ////////////////////////////////////////////////////////////////////
526 //
527 // return Coulomb scattering differential xsc with Wentzel correction  
528 
529 inline  G4double G4DiffuseElastic::GetCoulombElasticXsc( const G4ParticleDefinition* particle, 
530                                  G4double theta, 
531                G4double momentum, 
532          G4double Z         )
533 {
534   G4double sinHalfTheta  = std::sin(0.5*theta);
535   G4double sinHalfTheta2 = sinHalfTheta*sinHalfTheta;
536   G4double beta          = CalculateParticleBeta( particle, momentum);
537   G4double z             = particle->GetPDGCharge();
538   G4double n             = CalculateZommerfeld( beta, z, Z );
539   G4double am            = CalculateAm( momentum, n, Z);
540   G4double k             = momentum/CLHEP::hbarc;
541   G4double ch            = 0.5*n/k;
542   G4double ch2           = ch*ch;
543   G4double xsc           = ch2/(sinHalfTheta2+am)/(sinHalfTheta2+am);
544 
545   return xsc;
546 }
547 
548 
549 ////////////////////////////////////////////////////////////////////
550 //
551 // return Coulomb scattering total xsc with Wentzel correction  
552 
553 inline  G4double G4DiffuseElastic::GetCoulombTotalXsc( const G4ParticleDefinition* particle,  
554                                            G4double momentum, G4double Z  )
555 {
556   G4double beta          = CalculateParticleBeta( particle, momentum);
557   G4cout<<"beta = "<<beta<<G4endl;
558   G4double z             = particle->GetPDGCharge();
559   G4double n             = CalculateZommerfeld( beta, z, Z );
560   G4cout<<"fZomerfeld = "<<n<<G4endl;
561   G4double am            = CalculateAm( momentum, n, Z);
562   G4cout<<"cof Am = "<<am<<G4endl;
563   G4double k             = momentum/CLHEP::hbarc;
564   G4cout<<"k = "<<k*CLHEP::fermi<<" 1/fermi"<<G4endl;
565   G4cout<<"k*Bohr_radius = "<<k*CLHEP::Bohr_radius<<G4endl;
566   G4double ch            = n/k;
567   G4double ch2           = ch*ch;
568   G4double xsc           = ch2*CLHEP::pi/(am +am*am);
569 
570   return xsc;
571 }
572 
573 ////////////////////////////////////////////////////////////////////
574 //
575 // return Coulomb scattering xsc with Wentzel correction  integrated between
576 // theta1 and < theta2
577 
578 inline  G4double G4DiffuseElastic::GetCoulombIntegralXsc( const G4ParticleDefinition* particle,  
579                G4double momentum, G4double Z, 
580                                  G4double theta1, G4double theta2 )
581 {
582   G4double c1 = std::cos(theta1);
583   G4cout<<"c1 = "<<c1<<G4endl;
584   G4double c2 = std::cos(theta2);
585   G4cout<<"c2 = "<<c2<<G4endl;
586   G4double beta          = CalculateParticleBeta( particle, momentum);
587   // G4cout<<"beta = "<<beta<<G4endl;
588   G4double z             = particle->GetPDGCharge();
589   G4double n             = CalculateZommerfeld( beta, z, Z );
590   // G4cout<<"fZomerfeld = "<<n<<G4endl;
591   G4double am            = CalculateAm( momentum, n, Z);
592   // G4cout<<"cof Am = "<<am<<G4endl;
593   G4double k             = momentum/CLHEP::hbarc;
594   // G4cout<<"k = "<<k*CLHEP::fermi<<" 1/fermi"<<G4endl;
595   // G4cout<<"k*Bohr_radius = "<<k*CLHEP::Bohr_radius<<G4endl;
596   G4double ch            = n/k;
597   G4double ch2           = ch*ch;
598   am *= 2.;
599   G4double xsc           = ch2*CLHEP::twopi*(c1-c2);
600            xsc          /= (1 - c1 + am)*(1 - c2 + am);
601 
602   return xsc;
603 }
604 
605 #endif
606