Geant4  9.6.p02
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4NuclNuclDiffuseElastic.hh
Go to the documentation of this file.
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 // $Id$
28 //
29 //
30 // G4 Model: optical elastic scattering with 4-momentum balance
31 //
32 // Class Description
33 // Final state production model for nucleus-nucleus elastic scattering;
34 // Coulomb amplitude is not considered as correction
35 // (as in G4DiffuseElastic)
36 // Class Description - End
37 //
38 //
39 // 17.03.09 V. Grichine implementation for Coulomb elastic scattering
40 
41 
42 #ifndef G4NuclNuclDiffuseElastic_h
43 #define G4NuclNuclDiffuseElastic_h 1
44 
45 #include <complex>
46 #include <CLHEP/Units/PhysicalConstants.h>
47 #include "globals.hh"
48 #include "G4Integrator.hh"
49 #include "G4HadronElastic.hh"
50 #include "G4HadProjectile.hh"
51 #include "G4Nucleus.hh"
52 
53 using namespace std;
54 
55 class G4ParticleDefinition;
56 class G4PhysicsTable;
57 class G4PhysicsLogVector;
58 
59 class G4NuclNuclDiffuseElastic : public G4HadronElastic // G4HadronicInteraction
60 {
61 public:
62 
63  G4NuclNuclDiffuseElastic();
64 
65  // G4NuclNuclDiffuseElastic(const G4ParticleDefinition* aParticle);
66 
67 
68 
69 
70 
71  virtual ~G4NuclNuclDiffuseElastic();
72 
73  void Initialise();
74 
75  void InitialiseOnFly(G4double Z, G4double A);
76 
77  void BuildAngleTable();
78 
79 
80  // G4HadFinalState * ApplyYourself(const G4HadProjectile & aTrack, G4Nucleus & targetNucleus);
81 
82  virtual G4double SampleInvariantT(const G4ParticleDefinition* p,
83  G4double plab,
84  G4int Z, G4int A);
85 
86  void SetPlabLowLimit(G4double value);
87 
88  void SetHEModelLowLimit(G4double value);
89 
90  void SetQModelLowLimit(G4double value);
91 
92  void SetLowestEnergyLimit(G4double value);
93 
94  void SetRecoilKinEnergyLimit(G4double value);
95 
96  G4double SampleT(const G4ParticleDefinition* aParticle,
97  G4double p, G4double A);
98 
99  G4double SampleTableT(const G4ParticleDefinition* aParticle,
100  G4double p, G4double Z, G4double A);
101 
102  G4double SampleThetaCMS(const G4ParticleDefinition* aParticle, G4double p, G4double A);
103 
104  G4double SampleTableThetaCMS(const G4ParticleDefinition* aParticle, G4double p,
105  G4double Z, G4double A);
106 
107  G4double GetScatteringAngle(G4int iMomentum, G4int iAngle, G4double position);
108 
109  G4double SampleThetaLab(const G4HadProjectile* aParticle,
110  G4double tmass, G4double A);
111 
112  G4double GetDiffuseElasticXsc( const G4ParticleDefinition* particle,
113  G4double theta,
114  G4double momentum,
115  G4double A );
116 
117  G4double GetInvElasticXsc( const G4ParticleDefinition* particle,
118  G4double theta,
119  G4double momentum,
120  G4double A, G4double Z );
121 
122  G4double GetDiffuseElasticSumXsc( const G4ParticleDefinition* particle,
123  G4double theta,
124  G4double momentum,
125  G4double A, G4double Z );
126 
127  G4double GetInvElasticSumXsc( const G4ParticleDefinition* particle,
128  G4double tMand,
129  G4double momentum,
130  G4double A, G4double Z );
131 
132  G4double IntegralElasticProb( const G4ParticleDefinition* particle,
133  G4double theta,
134  G4double momentum,
135  G4double A );
136 
137 
138  G4double GetCoulombElasticXsc( const G4ParticleDefinition* particle,
139  G4double theta,
140  G4double momentum,
141  G4double Z );
142 
143  G4double GetRutherfordXsc( G4double theta );
144 
145  G4double GetInvCoulombElasticXsc( const G4ParticleDefinition* particle,
146  G4double tMand,
147  G4double momentum,
148  G4double A, G4double Z );
149 
150  G4double GetCoulombTotalXsc( const G4ParticleDefinition* particle,
151  G4double momentum, G4double Z );
152 
153  G4double GetCoulombIntegralXsc( const G4ParticleDefinition* particle,
154  G4double momentum, G4double Z,
155  G4double theta1, G4double theta2 );
156 
157 
158  G4double CalculateParticleBeta( const G4ParticleDefinition* particle,
159  G4double momentum );
160 
161  G4double CalculateZommerfeld( G4double beta, G4double Z1, G4double Z2 );
162 
163  G4double CalculateAm( G4double momentum, G4double n, G4double Z);
164 
165  G4double CalculateNuclearRad( G4double A);
166 
167  G4double ThetaCMStoThetaLab(const G4DynamicParticle* aParticle,
168  G4double tmass, G4double thetaCMS);
169 
170  G4double ThetaLabToThetaCMS(const G4DynamicParticle* aParticle,
171  G4double tmass, G4double thetaLab);
172 
173  void TestAngleTable(const G4ParticleDefinition* theParticle, G4double partMom,
174  G4double Z, G4double A);
175 
176 
177 
178  G4double BesselJzero(G4double z);
179  G4double BesselJone(G4double z);
180  G4double DampFactor(G4double z);
181  G4double BesselOneByArg(G4double z);
182 
183  G4double GetDiffElasticProb(G4double theta);
184  G4double GetDiffElasticSumProb(G4double theta);
185  G4double GetDiffElasticSumProbA(G4double alpha);
186  G4double GetIntegrandFunction(G4double theta);
187 
188  G4double GetNuclearRadius(){return fNuclearRadius;};
189 
190 
191  // Technical math functions for strong Coulomb contribution
192 
193  G4complex GammaLogarithm(G4complex xx);
194  G4complex GammaLogB2n(G4complex xx);
195 
196  G4double GetErf(G4double x);
197 
198  G4double GetCosHaPit2(G4double t){return std::cos(CLHEP::halfpi*t*t);};
199  G4double GetSinHaPit2(G4double t){return std::sin(CLHEP::halfpi*t*t);};
200 
201  G4double GetCint(G4double x);
202  G4double GetSint(G4double x);
203 
204 
205  G4complex GetErfcComp(G4complex z, G4int nMax);
206  G4complex GetErfcSer(G4complex z, G4int nMax);
207  G4complex GetErfcInt(G4complex z); // , G4int nMax);
208 
209  G4complex GetErfComp(G4complex z, G4int nMax); // AandS algorithm != Ser, Int
210  G4complex GetErfSer(G4complex z, G4int nMax);
211 
212  G4double GetExpCos(G4double x);
213  G4double GetExpSin(G4double x);
214  G4complex GetErfInt(G4complex z); // , G4int nMax);
215 
216 
217 
218 
219  G4double GetLegendrePol(G4int n, G4double x);
220 
221  G4complex TestErfcComp(G4complex z, G4int nMax);
222  G4complex TestErfcSer(G4complex z, G4int nMax);
223  G4complex TestErfcInt(G4complex z); // , G4int nMax);
224 
225  G4complex CoulombAmplitude(G4double theta);
226  G4double CoulombAmplitudeMod2(G4double theta);
227 
228 
229  void CalculateCoulombPhaseZero();
230  G4double CalculateCoulombPhase(G4int n);
231  void CalculateRutherfordAnglePar();
232 
233  G4double ProfileNear(G4double theta);
234  G4double ProfileFar(G4double theta);
235  G4double Profile(G4double theta);
236 
237  G4complex PhaseNear(G4double theta);
238  G4complex PhaseFar(G4double theta);
239 
240  G4complex GammaLess(G4double theta);
241  G4complex GammaMore(G4double theta);
242 
243  G4complex AmplitudeNear(G4double theta);
244  G4complex AmplitudeFar(G4double theta);
245 
246  G4complex Amplitude(G4double theta);
247  G4double AmplitudeMod2(G4double theta);
248 
249  G4complex AmplitudeSim(G4double theta);
250  G4double AmplitudeSimMod2(G4double theta);
251 
252  G4double GetRatioSim(G4double theta);
253  G4double GetRatioGen(G4double theta);
254 
255  G4double GetFresnelDiffuseXsc(G4double theta);
256  G4double GetFresnelIntegrandXsc(G4double alpha);
257 
258 
259  G4complex AmplitudeGla(G4double theta);
260  G4double AmplitudeGlaMod2(G4double theta);
261 
262  G4complex AmplitudeGG(G4double theta);
263  G4double AmplitudeGGMod2(G4double theta);
264 
265  void InitParameters(const G4ParticleDefinition* theParticle,
266  G4double partMom, G4double Z, G4double A);
267 
268  void InitDynParameters(const G4ParticleDefinition* theParticle,
269  G4double partMom);
270 
271  void InitParametersGla(const G4DynamicParticle* aParticle,
272  G4double partMom, G4double Z, G4double A);
273 
274  G4double GetHadronNucleonXscNS( G4ParticleDefinition* pParticle,
275  G4double pTkin,
276  G4ParticleDefinition* tParticle);
277 
278  G4double CalcMandelstamS( const G4double mp ,
279  const G4double mt ,
280  const G4double Plab );
281 
282  G4double GetProfileLambda(){return fProfileLambda;};
283 
284  void SetProfileLambda(G4double pl) {fProfileLambda = pl;};
285  void SetProfileDelta(G4double pd) {fProfileDelta = pd;};
286  void SetProfileAlpha(G4double pa){fProfileAlpha = pa;};
287  void SetCofLambda(G4double pa){fCofLambda = pa;};
288 
289  void SetCofAlpha(G4double pa){fCofAlpha = pa;};
290  void SetCofAlphaMax(G4double pa){fCofAlphaMax = pa;};
291  void SetCofAlphaCoulomb(G4double pa){fCofAlphaCoulomb = pa;};
292 
293  void SetCofDelta(G4double pa){fCofDelta = pa;};
294  void SetCofPhase(G4double pa){fCofPhase = pa;};
295  void SetCofFar(G4double pa){fCofFar = pa;};
296  void SetEtaRatio(G4double pa){fEtaRatio = pa;};
297  void SetMaxL(G4int l){fMaxL = l;};
298  void SetNuclearRadiusCof(G4double r){fNuclearRadiusCof = r;};
299 
300  G4double GetCofAlphaMax(){return fCofAlphaMax;};
301  G4double GetCofAlphaCoulomb(){return fCofAlphaCoulomb;};
302 
303 private:
304 
305 
306  G4ParticleDefinition* theProton;
307  G4ParticleDefinition* theNeutron;
308  G4ParticleDefinition* theDeuteron;
309  G4ParticleDefinition* theAlpha;
310 
311  const G4ParticleDefinition* thePionPlus;
312  const G4ParticleDefinition* thePionMinus;
313 
314  G4double lowEnergyRecoilLimit;
315  G4double lowEnergyLimitHE;
316  G4double lowEnergyLimitQ;
317  G4double lowestEnergyLimit;
318  G4double plabLowLimit;
319 
320  G4int fEnergyBin;
321  G4int fAngleBin;
322 
323  G4PhysicsLogVector* fEnergyVector;
324  G4PhysicsTable* fAngleTable;
325  std::vector<G4PhysicsTable*> fAngleBank;
326 
327  std::vector<G4double> fElementNumberVector;
328  std::vector<G4String> fElementNameVector;
329 
330  const G4ParticleDefinition* fParticle;
331 
332  G4double fWaveVector;
333  G4double fAtomicWeight;
334  G4double fAtomicNumber;
335 
336  G4double fNuclearRadius1;
337  G4double fNuclearRadius2;
338  G4double fNuclearRadius;
339  G4double fNuclearRadiusSquare;
340  G4double fNuclearRadiusCof;
341 
342  G4double fBeta;
343  G4double fZommerfeld;
344  G4double fRutherfordRatio;
345  G4double fAm;
346  G4bool fAddCoulomb;
347 
348  G4double fCoulombPhase0;
349  G4double fHalfRutThetaTg;
350  G4double fHalfRutThetaTg2;
351  G4double fRutherfordTheta;
352 
353  G4double fProfileLambda;
354  G4double fProfileDelta;
355  G4double fProfileAlpha;
356 
357  G4double fCofLambda;
358  G4double fCofAlpha;
359  G4double fCofDelta;
360  G4double fCofPhase;
361  G4double fCofFar;
362 
363  G4double fCofAlphaMax;
364  G4double fCofAlphaCoulomb;
365 
366  G4int fMaxL;
367  G4double fSumSigma;
368  G4double fEtaRatio;
369 
370  G4double fReZ;
371 
372 };
373 
374 
375 inline void G4NuclNuclDiffuseElastic::SetRecoilKinEnergyLimit(G4double value)
376 {
377  lowEnergyRecoilLimit = value;
378 }
379 
380 inline void G4NuclNuclDiffuseElastic::SetPlabLowLimit(G4double value)
381 {
382  plabLowLimit = value;
383 }
384 
385 inline void G4NuclNuclDiffuseElastic::SetHEModelLowLimit(G4double value)
386 {
387  lowEnergyLimitHE = value;
388 }
389 
390 inline void G4NuclNuclDiffuseElastic::SetQModelLowLimit(G4double value)
391 {
392  lowEnergyLimitQ = value;
393 }
394 
395 inline void G4NuclNuclDiffuseElastic::SetLowestEnergyLimit(G4double value)
396 {
397  lowestEnergyLimit = value;
398 }
399 
400 
402 //
403 // Bessel J0 function based on rational approximation from
404 // J.F. Hart, Computer Approximations, New York, Willey 1968, p. 141
405 
406 inline G4double G4NuclNuclDiffuseElastic::BesselJzero(G4double value)
407 {
408  G4double modvalue, value2, fact1, fact2, arg, shift, bessel;
409 
410  modvalue = fabs(value);
411 
412  if ( value < 8.0 && value > -8.0 )
413  {
414  value2 = value*value;
415 
416  fact1 = 57568490574.0 + value2*(-13362590354.0
417  + value2*( 651619640.7
418  + value2*(-11214424.18
419  + value2*( 77392.33017
420  + value2*(-184.9052456 ) ) ) ) );
421 
422  fact2 = 57568490411.0 + value2*( 1029532985.0
423  + value2*( 9494680.718
424  + value2*(59272.64853
425  + value2*(267.8532712
426  + value2*1.0 ) ) ) );
427 
428  bessel = fact1/fact2;
429  }
430  else
431  {
432  arg = 8.0/modvalue;
433 
434  value2 = arg*arg;
435 
436  shift = modvalue-0.785398164;
437 
438  fact1 = 1.0 + value2*(-0.1098628627e-2
439  + value2*(0.2734510407e-4
440  + value2*(-0.2073370639e-5
441  + value2*0.2093887211e-6 ) ) );
442 
443  fact2 = -0.1562499995e-1 + value2*(0.1430488765e-3
444  + value2*(-0.6911147651e-5
445  + value2*(0.7621095161e-6
446  - value2*0.934945152e-7 ) ) );
447 
448  bessel = sqrt(0.636619772/modvalue)*(cos(shift)*fact1 - arg*sin(shift)*fact2 );
449  }
450  return bessel;
451 }
452 
454 //
455 // Bessel J1 function based on rational approximation from
456 // J.F. Hart, Computer Approximations, New York, Willey 1968, p. 141
457 
458 inline G4double G4NuclNuclDiffuseElastic::BesselJone(G4double value)
459 {
460  G4double modvalue, value2, fact1, fact2, arg, shift, bessel;
461 
462  modvalue = fabs(value);
463 
464  if ( modvalue < 8.0 )
465  {
466  value2 = value*value;
467 
468  fact1 = value*(72362614232.0 + value2*(-7895059235.0
469  + value2*( 242396853.1
470  + value2*(-2972611.439
471  + value2*( 15704.48260
472  + value2*(-30.16036606 ) ) ) ) ) );
473 
474  fact2 = 144725228442.0 + value2*(2300535178.0
475  + value2*(18583304.74
476  + value2*(99447.43394
477  + value2*(376.9991397
478  + value2*1.0 ) ) ) );
479  bessel = fact1/fact2;
480  }
481  else
482  {
483  arg = 8.0/modvalue;
484 
485  value2 = arg*arg;
486 
487  shift = modvalue - 2.356194491;
488 
489  fact1 = 1.0 + value2*( 0.183105e-2
490  + value2*(-0.3516396496e-4
491  + value2*(0.2457520174e-5
492  + value2*(-0.240337019e-6 ) ) ) );
493 
494  fact2 = 0.04687499995 + value2*(-0.2002690873e-3
495  + value2*( 0.8449199096e-5
496  + value2*(-0.88228987e-6
497  + value2*0.105787412e-6 ) ) );
498 
499  bessel = sqrt( 0.636619772/modvalue)*(cos(shift)*fact1 - arg*sin(shift)*fact2);
500 
501  if (value < 0.0) bessel = -bessel;
502  }
503  return bessel;
504 }
505 
507 //
508 // damp factor in diffraction x/sh(x), x was already *pi
509 
510 inline G4double G4NuclNuclDiffuseElastic::DampFactor(G4double x)
511 {
512  G4double df;
513  G4double f2 = 2., f3 = 6., f4 = 24.; // first factorials
514 
515  // x *= pi;
516 
517  if( std::fabs(x) < 0.01 )
518  {
519  df = 1./(1. + x/f2 + x*x/f3 + x*x*x/f4);
520  }
521  else
522  {
523  df = x/std::sinh(x);
524  }
525  return df;
526 }
527 
528 
530 //
531 // return J1(x)/x with special case for small x
532 
533 inline G4double G4NuclNuclDiffuseElastic::BesselOneByArg(G4double x)
534 {
535  G4double x2, result;
536 
537  if( std::fabs(x) < 0.01 )
538  {
539  x *= 0.5;
540  x2 = x*x;
541  result = 2. - x2 + x2*x2/6.;
542  }
543  else
544  {
545  result = BesselJone(x)/x;
546  }
547  return result;
548 }
549 
551 //
552 // return particle beta
553 
554 inline G4double G4NuclNuclDiffuseElastic::CalculateParticleBeta( const G4ParticleDefinition* particle,
555  G4double momentum )
556 {
557  G4double mass = particle->GetPDGMass();
558  G4double a = momentum/mass;
559  fBeta = a/std::sqrt(1+a*a);
560 
561  return fBeta;
562 }
563 
565 //
566 // return Zommerfeld parameter for Coulomb scattering
567 
568 inline G4double G4NuclNuclDiffuseElastic::CalculateZommerfeld( G4double beta, G4double Z1, G4double Z2 )
569 {
570  fZommerfeld = CLHEP::fine_structure_const*Z1*Z2/beta;
571 
572  return fZommerfeld;
573 }
574 
576 //
577 // return Wentzel correction for Coulomb scattering
578 
579 inline G4double G4NuclNuclDiffuseElastic::CalculateAm( G4double momentum, G4double n, G4double Z)
580 {
581  G4double k = momentum/CLHEP::hbarc;
582  G4double ch = 1.13 + 3.76*n*n;
583  G4double zn = 1.77*k*std::pow(Z,-1./3.)*CLHEP::Bohr_radius;
584  G4double zn2 = zn*zn;
585  fAm = ch/zn2;
586 
587  return fAm;
588 }
589 
591 //
592 // calculate nuclear radius for different atomic weights using different approximations
593 
594 inline G4double G4NuclNuclDiffuseElastic::CalculateNuclearRad( G4double A)
595 {
596  G4double r0 = 1.*CLHEP::fermi, radius;
597  // r0 *= 1.12;
598  // r0 *= 1.44;
599  r0 *= fNuclearRadiusCof;
600 
601  /*
602  if( A < 50. )
603  {
604  if( A > 10. ) r0 = 1.16*( 1 - std::pow(A, -2./3.) )*CLHEP::fermi; // 1.08*fermi;
605  else r0 = 1.1*CLHEP::fermi;
606 
607  radius = r0*std::pow(A, 1./3.);
608  }
609  else
610  {
611  r0 = 1.7*CLHEP::fermi; // 1.7*fermi;
612 
613  radius = r0*std::pow(A, 0.27); // 0.27);
614  }
615  */
616  radius = r0*std::pow(A, 1./3.);
617 
618  return radius;
619 }
620 
622 //
623 // return Coulomb scattering differential xsc with Wentzel correction. Test function
624 
625 inline G4double G4NuclNuclDiffuseElastic::GetCoulombElasticXsc( const G4ParticleDefinition* particle,
626  G4double theta,
627  G4double momentum,
628  G4double Z )
629 {
630  G4double sinHalfTheta = std::sin(0.5*theta);
631  G4double sinHalfTheta2 = sinHalfTheta*sinHalfTheta;
632  G4double beta = CalculateParticleBeta( particle, momentum);
633  G4double z = particle->GetPDGCharge();
634  G4double n = CalculateZommerfeld( beta, z, Z );
635  G4double am = CalculateAm( momentum, n, Z);
636  G4double k = momentum/CLHEP::hbarc;
637  G4double ch = 0.5*n/k;
638  G4double ch2 = ch*ch;
639  G4double xsc = ch2/(sinHalfTheta2+am)/(sinHalfTheta2+am);
640 
641  return xsc;
642 }
643 
645 //
646 // return Rutherford scattering differential xsc with Wentzel correction. For Sampling.
647 
648 inline G4double G4NuclNuclDiffuseElastic::GetRutherfordXsc( G4double theta )
649 {
650  G4double sinHalfTheta = std::sin(0.5*theta);
651  G4double sinHalfTheta2 = sinHalfTheta*sinHalfTheta;
652 
653  G4double ch2 = fRutherfordRatio*fRutherfordRatio;
654 
655  G4double xsc = ch2/(sinHalfTheta2+fAm)/(sinHalfTheta2+fAm);
656 
657  return xsc;
658 }
659 
660 
662 //
663 // return Coulomb scattering total xsc with Wentzel correction
664 
665 inline G4double G4NuclNuclDiffuseElastic::GetCoulombTotalXsc( const G4ParticleDefinition* particle,
666  G4double momentum, G4double Z )
667 {
668  G4double beta = CalculateParticleBeta( particle, momentum);
669  G4cout<<"beta = "<<beta<<G4endl;
670  G4double z = particle->GetPDGCharge();
671  G4double n = CalculateZommerfeld( beta, z, Z );
672  G4cout<<"fZomerfeld = "<<n<<G4endl;
673  G4double am = CalculateAm( momentum, n, Z);
674  G4cout<<"cof Am = "<<am<<G4endl;
675  G4double k = momentum/CLHEP::hbarc;
676  G4cout<<"k = "<<k*CLHEP::fermi<<" 1/fermi"<<G4endl;
677  G4cout<<"k*Bohr_radius = "<<k*CLHEP::Bohr_radius<<G4endl;
678  G4double ch = n/k;
679  G4double ch2 = ch*ch;
680  G4double xsc = ch2*CLHEP::pi/(am +am*am);
681 
682  return xsc;
683 }
684 
686 //
687 // return Coulomb scattering xsc with Wentzel correction integrated between
688 // theta1 and < theta2
689 
690 inline G4double G4NuclNuclDiffuseElastic::GetCoulombIntegralXsc( const G4ParticleDefinition* particle,
691  G4double momentum, G4double Z,
692  G4double theta1, G4double theta2 )
693 {
694  G4double c1 = std::cos(theta1);
695  G4cout<<"c1 = "<<c1<<G4endl;
696  G4double c2 = std::cos(theta2);
697  G4cout<<"c2 = "<<c2<<G4endl;
698  G4double beta = CalculateParticleBeta( particle, momentum);
699  // G4cout<<"beta = "<<beta<<G4endl;
700  G4double z = particle->GetPDGCharge();
701  G4double n = CalculateZommerfeld( beta, z, Z );
702  // G4cout<<"fZomerfeld = "<<n<<G4endl;
703  G4double am = CalculateAm( momentum, n, Z);
704  // G4cout<<"cof Am = "<<am<<G4endl;
705  G4double k = momentum/CLHEP::hbarc;
706  // G4cout<<"k = "<<k*CLHEP::fermi<<" 1/fermi"<<G4endl;
707  // G4cout<<"k*Bohr_radius = "<<k*CLHEP::Bohr_radius<<G4endl;
708  G4double ch = n/k;
709  G4double ch2 = ch*ch;
710  am *= 2.;
711  G4double xsc = ch2*CLHEP::twopi*(c1-c2);
712  xsc /= (1 - c1 + am)*(1 - c2 + am);
713 
714  return xsc;
715 }
716 
718 //
719 // For the calculation of arg Gamma(z) one needs complex extension
720 // of ln(Gamma(z))
721 
722 inline G4complex G4NuclNuclDiffuseElastic::GammaLogarithm(G4complex zz)
723 {
724  static G4double cof[6] = { 76.18009172947146, -86.50532032941677,
725  24.01409824083091, -1.231739572450155,
726  0.1208650973866179e-2, -0.5395239384953e-5 } ;
727  register G4int j;
728  G4complex z = zz - 1.0;
729  G4complex tmp = z + 5.5;
730  tmp -= (z + 0.5) * std::log(tmp);
731  G4complex ser = G4complex(1.000000000190015,0.);
732 
733  for ( j = 0; j <= 5; j++ )
734  {
735  z += 1.0;
736  ser += cof[j]/z;
737  }
738  return -tmp + std::log(2.5066282746310005*ser);
739 }
740 
742 //
743 // For the calculation of arg Gamma(z) one needs complex extension
744 // of ln(Gamma(z)) here is approximate algorithm
745 
746 inline G4complex G4NuclNuclDiffuseElastic::GammaLogB2n(G4complex z)
747 {
748  G4complex z1 = 12.*z;
749  G4complex z2 = z*z;
750  G4complex z3 = z2*z;
751  G4complex z5 = z2*z3;
752  G4complex z7 = z2*z5;
753 
754  z3 *= 360.;
755  z5 *= 1260.;
756  z7 *= 1680.;
757 
758  G4complex result = (z-0.5)*std::log(z)-z+0.5*std::log(CLHEP::twopi);
759  result += 1./z1 - 1./z3 +1./z5 -1./z7;
760  return result;
761 }
762 
764 //
765 //
766 
767 inline G4double G4NuclNuclDiffuseElastic::GetErf(G4double x)
768 {
769  G4double t, z, tmp, result;
770 
771  z = std::fabs(x);
772  t = 1.0/(1.0+0.5*z);
773 
774  tmp = t*exp(-z*z-1.26551223+t*(1.00002368+t*(0.37409196+t*(0.09678418+
775  t*(-0.18628806+t*(0.27886807+t*(-1.13520398+t*(1.48851587+
776  t*(-0.82215223+t*0.17087277)))))))));
777 
778  if( x >= 0.) result = 1. - tmp;
779  else result = 1. + tmp;
780 
781  return result;
782 }
783 
785 //
786 //
787 
788 inline G4complex G4NuclNuclDiffuseElastic::GetErfcComp(G4complex z, G4int nMax)
789 {
790  G4complex erfcz = 1. - GetErfComp( z, nMax);
791  return erfcz;
792 }
793 
795 //
796 //
797 
798 inline G4complex G4NuclNuclDiffuseElastic::GetErfcSer(G4complex z, G4int nMax)
799 {
800  G4complex erfcz = 1. - GetErfSer( z, nMax);
801  return erfcz;
802 }
803 
805 //
806 //
807 
808 inline G4complex G4NuclNuclDiffuseElastic::GetErfcInt(G4complex z) // , G4int nMax)
809 {
810  G4complex erfcz = 1. - GetErfInt( z); // , nMax);
811  return erfcz;
812 }
813 
814 inline G4double G4NuclNuclDiffuseElastic::GetLegendrePol(G4int n, G4double theta)
815 {
816  G4double legPol, epsilon = 1.e-6;
817  G4double x = std::cos(theta);
818 
819  if ( n < 0 ) legPol = 0.;
820  else if( n == 0 ) legPol = 1.;
821  else if( n == 1 ) legPol = x;
822  else if( n == 2 ) legPol = (3.*x*x-1.)/2.;
823  else if( n == 3 ) legPol = (5.*x*x*x-3.*x)/2.;
824  else if( n == 4 ) legPol = (35.*x*x*x*x-30.*x*x+3.)/8.;
825  else if( n == 5 ) legPol = (63.*x*x*x*x*x-70.*x*x*x+15.*x)/8.;
826  else if( n == 6 ) legPol = (231.*x*x*x*x*x*x-315.*x*x*x*x+105.*x*x-5.)/16.;
827  else
828  {
829  // legPol = ( (2*n-1)*x*GetLegendrePol(n-1,x) - (n-1)*GetLegendrePol(n-2,x) )/n;
830 
831  legPol = std::sqrt( 2./(n*CLHEP::pi*std::sin(theta+epsilon)) )*std::sin( (n+0.5)*theta+0.25*CLHEP::pi );
832  }
833  return legPol;
834 }
835 
836 
837 
839 //
840 //
841 
842 inline G4complex G4NuclNuclDiffuseElastic::TestErfcComp(G4complex z, G4int nMax)
843 {
844  G4complex miz = G4complex( z.imag(), -z.real() );
845  G4complex erfcz = 1. - GetErfComp( miz, nMax);
846  G4complex w = std::exp(-z*z)*erfcz;
847  return w;
848 }
849 
851 //
852 //
853 
854 inline G4complex G4NuclNuclDiffuseElastic::TestErfcSer(G4complex z, G4int nMax)
855 {
856  G4complex miz = G4complex( z.imag(), -z.real() );
857  G4complex erfcz = 1. - GetErfSer( miz, nMax);
858  G4complex w = std::exp(-z*z)*erfcz;
859  return w;
860 }
861 
863 //
864 //
865 
866 inline G4complex G4NuclNuclDiffuseElastic::TestErfcInt(G4complex z) // , G4int nMax)
867 {
868  G4complex miz = G4complex( z.imag(), -z.real() );
869  G4complex erfcz = 1. - GetErfInt( miz); // , nMax);
870  G4complex w = std::exp(-z*z)*erfcz;
871  return w;
872 }
873 
875 //
876 //
877 
878 inline G4complex G4NuclNuclDiffuseElastic::GetErfComp(G4complex z, G4int nMax)
879 {
880  G4int n;
881  G4double n2, cofn, shny, chny, fn, gn;
882 
883  G4double x = z.real();
884  G4double y = z.imag();
885 
886  G4double outRe = 0., outIm = 0.;
887 
888  G4double twox = 2.*x;
889  G4double twoxy = twox*y;
890  G4double twox2 = twox*twox;
891 
892  G4double cof1 = std::exp(-x*x)/CLHEP::pi;
893 
894  G4double cos2xy = std::cos(twoxy);
895  G4double sin2xy = std::sin(twoxy);
896 
897  G4double twoxcos2xy = twox*cos2xy;
898  G4double twoxsin2xy = twox*sin2xy;
899 
900  for( n = 1; n <= nMax; n++)
901  {
902  n2 = n*n;
903 
904  cofn = std::exp(-0.5*n2)/(n2+twox2); // /(n2+0.5*twox2);
905 
906  chny = std::cosh(n*y);
907  shny = std::sinh(n*y);
908 
909  fn = twox - twoxcos2xy*chny + n*sin2xy*shny;
910  gn = twoxsin2xy*chny + n*cos2xy*shny;
911 
912  fn *= cofn;
913  gn *= cofn;
914 
915  outRe += fn;
916  outIm += gn;
917  }
918  outRe *= 2*cof1;
919  outIm *= 2*cof1;
920 
921  if(std::abs(x) < 0.0001)
922  {
923  outRe += GetErf(x);
924  outIm += cof1*y;
925  }
926  else
927  {
928  outRe += GetErf(x) + cof1*(1-cos2xy)/twox;
929  outIm += cof1*sin2xy/twox;
930  }
931  return G4complex(outRe, outIm);
932 }
933 
935 //
936 //
937 
938 inline G4complex G4NuclNuclDiffuseElastic::GetErfSer(G4complex z, G4int nMax)
939 {
940  G4int n;
941  G4double a =1., b = 1., tmp;
942  G4complex sum = z, d = z;
943 
944  for( n = 1; n <= nMax; n++)
945  {
946  a *= 2.;
947  b *= 2.*n +1.;
948  d *= z*z;
949 
950  tmp = a/b;
951 
952  sum += tmp*d;
953  }
954  sum *= 2.*std::exp(-z*z)/std::sqrt(CLHEP::pi);
955 
956  return sum;
957 }
958 
960 
961 inline G4double G4NuclNuclDiffuseElastic::GetExpCos(G4double x)
962 {
963  G4double result;
964 
965  result = std::exp(x*x-fReZ*fReZ);
966  result *= std::cos(2.*x*fReZ);
967  return result;
968 }
969 
971 
972 inline G4double G4NuclNuclDiffuseElastic::GetExpSin(G4double x)
973 {
974  G4double result;
975 
976  result = std::exp(x*x-fReZ*fReZ);
977  result *= std::sin(2.*x*fReZ);
978  return result;
979 }
980 
981 
982 
984 //
985 //
986 
987 inline G4complex G4NuclNuclDiffuseElastic::GetErfInt(G4complex z) // , G4int nMax)
988 {
989  G4double outRe, outIm;
990 
991  G4double x = z.real();
992  G4double y = z.imag();
993  fReZ = x;
994 
995  G4Integrator<G4NuclNuclDiffuseElastic,G4double(G4NuclNuclDiffuseElastic::*)(G4double)> integral;
996 
997  outRe = integral.Legendre96(this,&G4NuclNuclDiffuseElastic::GetExpSin, 0., y );
998  outIm = integral.Legendre96(this,&G4NuclNuclDiffuseElastic::GetExpCos, 0., y );
999 
1000  outRe *= 2./sqrt(CLHEP::pi);
1001  outIm *= 2./sqrt(CLHEP::pi);
1002 
1003  outRe += GetErf(x);
1004 
1005  return G4complex(outRe, outIm);
1006 }
1007 
1009 //
1010 //
1011 
1012 inline G4double G4NuclNuclDiffuseElastic::GetCint(G4double x)
1013 {
1014  G4double out;
1015 
1016  G4Integrator<G4NuclNuclDiffuseElastic,G4double(G4NuclNuclDiffuseElastic::*)(G4double)> integral;
1017 
1018  out= integral.Legendre96(this,&G4NuclNuclDiffuseElastic::GetCosHaPit2, 0., x );
1019 
1020  return out;
1021 }
1022 
1024 //
1025 //
1026 
1027 inline G4double G4NuclNuclDiffuseElastic::GetSint(G4double x)
1028 {
1029  G4double out;
1030 
1031  G4Integrator<G4NuclNuclDiffuseElastic,G4double(G4NuclNuclDiffuseElastic::*)(G4double)> integral;
1032 
1033  out= integral.Legendre96(this,&G4NuclNuclDiffuseElastic::GetSinHaPit2, 0., x );
1034 
1035  return out;
1036 }
1037 
1038 
1040 //
1041 //
1042 
1043 inline G4complex G4NuclNuclDiffuseElastic::CoulombAmplitude(G4double theta)
1044 {
1045  G4complex ca;
1046 
1047  G4double sinHalfTheta = std::sin(0.5*theta);
1048  G4double sinHalfTheta2 = sinHalfTheta*sinHalfTheta;
1049  sinHalfTheta2 += fAm;
1050 
1051  G4double order = 2.*fCoulombPhase0 - fZommerfeld*std::log(sinHalfTheta2);
1052  G4complex z = G4complex(0., order);
1053  ca = std::exp(z);
1054 
1055  ca *= -fZommerfeld/(2.*fWaveVector*sinHalfTheta2);
1056 
1057  return ca;
1058 }
1059 
1061 //
1062 //
1063 
1064 inline G4double G4NuclNuclDiffuseElastic::CoulombAmplitudeMod2(G4double theta)
1065 {
1066  G4complex ca = CoulombAmplitude(theta);
1067  G4double out = ca.real()*ca.real() + ca.imag()*ca.imag();
1068 
1069  return out;
1070 }
1071 
1073 //
1074 //
1075 
1076 
1077 inline void G4NuclNuclDiffuseElastic::CalculateCoulombPhaseZero()
1078 {
1079  G4complex z = G4complex(1,fZommerfeld);
1080  // G4complex gammalog = GammaLogarithm(z);
1081  G4complex gammalog = GammaLogB2n(z);
1082  fCoulombPhase0 = gammalog.imag();
1083 }
1084 
1086 //
1087 //
1088 
1089 
1090 inline G4double G4NuclNuclDiffuseElastic::CalculateCoulombPhase(G4int n)
1091 {
1092  G4complex z = G4complex(1. + n, fZommerfeld);
1093  // G4complex gammalog = GammaLogarithm(z);
1094  G4complex gammalog = GammaLogB2n(z);
1095  return gammalog.imag();
1096 }
1097 
1098 
1100 //
1101 //
1102 
1103 
1104 inline void G4NuclNuclDiffuseElastic::CalculateRutherfordAnglePar()
1105 {
1106  fHalfRutThetaTg = fZommerfeld/fProfileLambda; // (fWaveVector*fNuclearRadius);
1107  fRutherfordTheta = 2.*std::atan(fHalfRutThetaTg);
1108  fHalfRutThetaTg2 = fHalfRutThetaTg*fHalfRutThetaTg;
1109  // G4cout<<"fRutherfordTheta = "<<fRutherfordTheta/degree<<" degree"<<G4endl;
1110 
1111 }
1112 
1114 //
1115 //
1116 
1117 inline G4double G4NuclNuclDiffuseElastic::ProfileNear(G4double theta)
1118 {
1119  G4double dTheta = fRutherfordTheta - theta;
1120  G4double result = 0., argument = 0.;
1121 
1122  if(std::abs(dTheta) < 0.001) result = fProfileAlpha*fProfileDelta;
1123  else
1124  {
1125  argument = fProfileDelta*dTheta;
1126  result = CLHEP::pi*argument*std::exp(fProfileAlpha*argument);
1127  result /= std::sinh(CLHEP::pi*argument);
1128  result -= 1.;
1129  result /= dTheta;
1130  }
1131  return result;
1132 }
1133 
1135 //
1136 //
1137 
1138 inline G4double G4NuclNuclDiffuseElastic::ProfileFar(G4double theta)
1139 {
1140  G4double dTheta = fRutherfordTheta + theta;
1141  G4double argument = fProfileDelta*dTheta;
1142 
1143  G4double result = CLHEP::pi*argument*std::exp(fProfileAlpha*argument);
1144  result /= std::sinh(CLHEP::pi*argument);
1145  result /= dTheta;
1146 
1147  return result;
1148 }
1149 
1151 //
1152 //
1153 
1154 inline G4double G4NuclNuclDiffuseElastic::Profile(G4double theta)
1155 {
1156  G4double dTheta = fRutherfordTheta - theta;
1157  G4double result = 0., argument = 0.;
1158 
1159  if(std::abs(dTheta) < 0.001) result = 1.;
1160  else
1161  {
1162  argument = fProfileDelta*dTheta;
1163  result = CLHEP::pi*argument;
1164  result /= std::sinh(CLHEP::pi*argument);
1165  }
1166  return result;
1167 }
1168 
1170 //
1171 //
1172 
1173 inline G4complex G4NuclNuclDiffuseElastic::PhaseNear(G4double theta)
1174 {
1175  G4double twosigma = 2.*fCoulombPhase0;
1176  twosigma -= fZommerfeld*std::log(fHalfRutThetaTg2/(1.+fHalfRutThetaTg2));
1177  twosigma += fRutherfordTheta*fZommerfeld/fHalfRutThetaTg - CLHEP::halfpi;
1178  twosigma -= fProfileLambda*theta - 0.25*CLHEP::pi;
1179 
1180  twosigma *= fCofPhase;
1181 
1182  G4complex z = G4complex(0., twosigma);
1183 
1184  return std::exp(z);
1185 }
1186 
1188 //
1189 //
1190 
1191 inline G4complex G4NuclNuclDiffuseElastic::PhaseFar(G4double theta)
1192 {
1193  G4double twosigma = 2.*fCoulombPhase0;
1194  twosigma -= fZommerfeld*std::log(fHalfRutThetaTg2/(1.+fHalfRutThetaTg2));
1195  twosigma += fRutherfordTheta*fZommerfeld/fHalfRutThetaTg - CLHEP::halfpi;
1196  twosigma += fProfileLambda*theta - 0.25*CLHEP::pi;
1197 
1198  twosigma *= fCofPhase;
1199 
1200  G4complex z = G4complex(0., twosigma);
1201 
1202  return std::exp(z);
1203 }
1204 
1206 //
1207 //
1208 
1209 
1210 inline G4complex G4NuclNuclDiffuseElastic::GammaLess(G4double theta)
1211 {
1212  G4double sinThetaR = 2.*fHalfRutThetaTg/(1. + fHalfRutThetaTg2);
1213  G4double cosHalfThetaR2 = 1./(1. + fHalfRutThetaTg2);
1214 
1215  G4double u = std::sqrt(0.5*fProfileLambda/sinThetaR);
1216  G4double kappa = u/std::sqrt(CLHEP::pi);
1217  G4double dTheta = theta - fRutherfordTheta;
1218  u *= dTheta;
1219  G4double u2 = u*u;
1220  G4double u2m2p3 = u2*2./3.;
1221 
1222  G4complex im = G4complex(0.,1.);
1223  G4complex order = G4complex(u,u);
1224  order /= std::sqrt(2.);
1225 
1226  G4complex gamma = CLHEP::pi*kappa*GetErfcInt(-order)*std::exp(im*(u*u+0.25*CLHEP::pi));
1227  G4complex a0 = 0.5*(1. + 4.*(1.+im*u2)*cosHalfThetaR2/3.)/sinThetaR;
1228  G4complex a1 = 0.5*(1. + 2.*(1.+im*u2m2p3)*cosHalfThetaR2)/sinThetaR;
1229  G4complex out = gamma*(1. - a1*dTheta) - a0;
1230 
1231  return out;
1232 }
1233 
1235 //
1236 //
1237 
1238 inline G4complex G4NuclNuclDiffuseElastic::GammaMore(G4double theta)
1239 {
1240  G4double sinThetaR = 2.*fHalfRutThetaTg/(1. + fHalfRutThetaTg2);
1241  G4double cosHalfThetaR2 = 1./(1. + fHalfRutThetaTg2);
1242 
1243  G4double u = std::sqrt(0.5*fProfileLambda/sinThetaR);
1244  G4double kappa = u/std::sqrt(CLHEP::pi);
1245  G4double dTheta = theta - fRutherfordTheta;
1246  u *= dTheta;
1247  G4double u2 = u*u;
1248  G4double u2m2p3 = u2*2./3.;
1249 
1250  G4complex im = G4complex(0.,1.);
1251  G4complex order = G4complex(u,u);
1252  order /= std::sqrt(2.);
1253  G4complex gamma = CLHEP::pi*kappa*GetErfcInt(order)*std::exp(im*(u*u+0.25*CLHEP::pi));
1254  G4complex a0 = 0.5*(1. + 4.*(1.+im*u2)*cosHalfThetaR2/3.)/sinThetaR;
1255  G4complex a1 = 0.5*(1. + 2.*(1.+im*u2m2p3)*cosHalfThetaR2)/sinThetaR;
1256  G4complex out = -gamma*(1. - a1*dTheta) - a0;
1257 
1258  return out;
1259 }
1260 
1262 //
1263 //
1264 
1265 inline G4complex G4NuclNuclDiffuseElastic::AmplitudeNear(G4double theta)
1266 {
1267  G4double kappa = std::sqrt(0.5*fProfileLambda/std::sin(theta)/CLHEP::pi);
1268  G4complex out = G4complex(kappa/fWaveVector,0.);
1269 
1270  out *= PhaseNear(theta);
1271 
1272  if( theta <= fRutherfordTheta )
1273  {
1274  out *= GammaLess(theta) + ProfileNear(theta);
1275  // out *= GammaMore(theta) + ProfileNear(theta);
1276  out += CoulombAmplitude(theta);
1277  }
1278  else
1279  {
1280  out *= GammaMore(theta) + ProfileNear(theta);
1281  // out *= GammaLess(theta) + ProfileNear(theta);
1282  }
1283  return out;
1284 }
1285 
1287 //
1288 //
1289 
1290 inline G4complex G4NuclNuclDiffuseElastic::AmplitudeFar(G4double theta)
1291 {
1292  G4double kappa = std::sqrt(0.5*fProfileLambda/std::sin(theta)/CLHEP::pi);
1293  G4complex out = G4complex(kappa/fWaveVector,0.);
1294  out *= ProfileFar(theta);
1295  out *= PhaseFar(theta);
1296  return out;
1297 }
1298 
1299 
1301 //
1302 //
1303 
1304 inline G4complex G4NuclNuclDiffuseElastic::Amplitude(G4double theta)
1305 {
1306 
1307  G4complex out = AmplitudeNear(theta) + fCofFar*AmplitudeFar(theta);
1308  // G4complex out = AmplitudeNear(theta);
1309  // G4complex out = AmplitudeFar(theta);
1310  return out;
1311 }
1312 
1314 //
1315 //
1316 
1317 inline G4double G4NuclNuclDiffuseElastic::AmplitudeMod2(G4double theta)
1318 {
1319  G4complex out = Amplitude(theta);
1320  G4double mod2 = out.real()*out.real() + out.imag()*out.imag();
1321  return mod2;
1322 }
1323 
1325 //
1326 //
1327 
1328 inline G4complex G4NuclNuclDiffuseElastic::AmplitudeSim(G4double theta)
1329 {
1330  G4double sinThetaR = 2.*fHalfRutThetaTg/(1. + fHalfRutThetaTg2);
1331  G4double dTheta = 0.5*(theta - fRutherfordTheta);
1332  G4double sindTheta = std::sin(dTheta);
1333  G4double persqrt2 = std::sqrt(0.5);
1334 
1335  G4complex order = G4complex(persqrt2,persqrt2);
1336  order *= std::sqrt(0.5*fProfileLambda/sinThetaR)*2.*sindTheta;
1337  // order *= std::sqrt(0.5*fProfileLambda/sinThetaR)*2.*dTheta;
1338 
1339  G4complex out;
1340 
1341  if ( theta <= fRutherfordTheta )
1342  {
1343  out = 1. - 0.5*GetErfcInt(-order)*ProfileNear(theta);
1344  }
1345  else
1346  {
1347  out = 0.5*GetErfcInt(order)*ProfileNear(theta);
1348  }
1349 
1350  out *= CoulombAmplitude(theta);
1351 
1352  return out;
1353 }
1354 
1356 //
1357 //
1358 
1359 inline G4double G4NuclNuclDiffuseElastic::GetRatioSim(G4double theta)
1360 {
1361  G4double sinThetaR = 2.*fHalfRutThetaTg/(1. + fHalfRutThetaTg2);
1362  G4double dTheta = 0.5*(theta - fRutherfordTheta);
1363  G4double sindTheta = std::sin(dTheta);
1364 
1365  G4double order = std::sqrt(fProfileLambda/sinThetaR/CLHEP::pi)*2.*sindTheta;
1366  // G4cout<<"order = "<<order<<G4endl;
1367  G4double cosFresnel = 0.5 - GetCint(order);
1368  G4double sinFresnel = 0.5 - GetSint(order);
1369 
1370  G4double out = 0.5*( cosFresnel*cosFresnel + sinFresnel*sinFresnel );
1371 
1372  return out;
1373 }
1374 
1376 //
1377 // The ratio el/ruth for Fresnel smooth nucleus profile
1378 
1379 inline G4double G4NuclNuclDiffuseElastic::GetRatioGen(G4double theta)
1380 {
1381  G4double sinThetaR = 2.*fHalfRutThetaTg/(1. + fHalfRutThetaTg2);
1382  G4double dTheta = 0.5*(theta - fRutherfordTheta);
1383  G4double sindTheta = std::sin(dTheta);
1384 
1385  G4double prof = Profile(theta);
1386  G4double prof2 = prof*prof;
1387  // G4double profmod = std::abs(prof);
1388  G4double order = std::sqrt(fProfileLambda/sinThetaR/CLHEP::pi)*2.*sindTheta;
1389 
1390  order = std::abs(order); // since sin changes sign!
1391  // G4cout<<"order = "<<order<<G4endl;
1392 
1393  G4double cosFresnel = GetCint(order);
1394  G4double sinFresnel = GetSint(order);
1395 
1396  G4double out;
1397 
1398  if ( theta <= fRutherfordTheta )
1399  {
1400  out = 1. + 0.5*( (0.5-cosFresnel)*(0.5-cosFresnel)+(0.5-sinFresnel)*(0.5-sinFresnel) )*prof2;
1401  out += ( cosFresnel + sinFresnel - 1. )*prof;
1402  }
1403  else
1404  {
1405  out = 0.5*( (0.5-cosFresnel)*(0.5-cosFresnel)+(0.5-sinFresnel)*(0.5-sinFresnel) )*prof2;
1406  }
1407 
1408  return out;
1409 }
1410 
1412 //
1413 // The xsc for Fresnel smooth nucleus profile
1414 
1415 inline G4double G4NuclNuclDiffuseElastic::GetFresnelDiffuseXsc(G4double theta)
1416 {
1417  G4double ratio = GetRatioGen(theta);
1418  G4double ruthXsc = GetRutherfordXsc(theta);
1419  G4double xsc = ratio*ruthXsc;
1420  return xsc;
1421 }
1422 
1424 //
1425 // The xsc for Fresnel smooth nucleus profile for integration
1426 
1427 inline G4double G4NuclNuclDiffuseElastic::GetFresnelIntegrandXsc(G4double alpha)
1428 {
1429  G4double theta = std::sqrt(alpha);
1430  G4double xsc = GetFresnelDiffuseXsc(theta);
1431  return xsc;
1432 }
1433 
1434 
1435 
1437 //
1438 //
1439 
1440 inline G4double G4NuclNuclDiffuseElastic::AmplitudeSimMod2(G4double theta)
1441 {
1442  G4complex out = AmplitudeSim(theta);
1443  G4double mod2 = out.real()*out.real() + out.imag()*out.imag();
1444  return mod2;
1445 }
1446 
1448 //
1449 //
1450 
1451 inline G4complex G4NuclNuclDiffuseElastic::AmplitudeGla(G4double theta)
1452 {
1453  G4int n;
1454  G4double T12b, b, b2; // cosTheta = std::cos(theta);
1455  G4complex out = G4complex(0.,0.), shiftC, shiftN;
1456  G4complex im = G4complex(0.,1.);
1457 
1458  for( n = 0; n < fMaxL; n++)
1459  {
1460  shiftC = std::exp( im*2.*CalculateCoulombPhase(n) );
1461  // b = ( fZommerfeld + std::sqrt( fZommerfeld*fZommerfeld + n*(n+1) ) )/fWaveVector;
1462  b = ( std::sqrt( G4double(n*(n+1)) ) )/fWaveVector;
1463  b2 = b*b;
1464  T12b = fSumSigma*std::exp(-b2/fNuclearRadiusSquare)/CLHEP::pi/fNuclearRadiusSquare;
1465  shiftN = std::exp( -0.5*(1.-im*fEtaRatio)*T12b ) - 1.;
1466  out += (2.*n+1.)*shiftC*shiftN*GetLegendrePol(n, theta);
1467  }
1468  out /= 2.*im*fWaveVector;
1469  out += CoulombAmplitude(theta);
1470  return out;
1471 }
1472 
1474 //
1475 //
1476 
1477 inline G4double G4NuclNuclDiffuseElastic::AmplitudeGlaMod2(G4double theta)
1478 {
1479  G4complex out = AmplitudeGla(theta);
1480  G4double mod2 = out.real()*out.real() + out.imag()*out.imag();
1481  return mod2;
1482 }
1483 
1484 
1486 //
1487 //
1488 
1489 inline G4complex G4NuclNuclDiffuseElastic::AmplitudeGG(G4double theta)
1490 {
1491  G4int n;
1492  G4double T12b, a, aTemp, b2, sinThetaH = std::sin(0.5*theta);
1493  G4double sinThetaH2 = sinThetaH*sinThetaH;
1494  G4complex out = G4complex(0.,0.);
1495  G4complex im = G4complex(0.,1.);
1496 
1497  a = -fSumSigma/CLHEP::twopi/fNuclearRadiusSquare;
1498  b2 = fWaveVector*fWaveVector*fNuclearRadiusSquare*sinThetaH2;
1499 
1500  aTemp = a;
1501 
1502  for( n = 1; n < fMaxL; n++)
1503  {
1504  T12b = aTemp*std::exp(-b2/n)/n;
1505  aTemp *= a;
1506  out += T12b;
1507  G4cout<<"out = "<<out<<G4endl;
1508  }
1509  out *= -4.*im*fWaveVector/CLHEP::pi;
1510  out += CoulombAmplitude(theta);
1511  return out;
1512 }
1513 
1515 //
1516 //
1517 
1518 inline G4double G4NuclNuclDiffuseElastic::AmplitudeGGMod2(G4double theta)
1519 {
1520  G4complex out = AmplitudeGG(theta);
1521  G4double mod2 = out.real()*out.real() + out.imag()*out.imag();
1522  return mod2;
1523 }
1524 
1525 
1527 //
1528 // Test for given particle and element table of momentum, angle probability.
1529 // For the partMom in CMS.
1530 
1531 inline void G4NuclNuclDiffuseElastic::InitParameters(const G4ParticleDefinition* theParticle,
1532  G4double partMom, G4double Z, G4double A)
1533 {
1534  fAtomicNumber = Z; // atomic number
1535  fAtomicWeight = A; // number of nucleons
1536 
1537  fNuclearRadius2 = CalculateNuclearRad(fAtomicWeight);
1538  G4double A1 = G4double( theParticle->GetBaryonNumber() );
1539  fNuclearRadius1 = CalculateNuclearRad(A1);
1540  // fNuclearRadius = std::sqrt(fNuclearRadius1*fNuclearRadius1+fNuclearRadius2*fNuclearRadius2);
1541  fNuclearRadius = fNuclearRadius1 + fNuclearRadius2;
1542 
1543  G4double a = 0.;
1544  G4double z = theParticle->GetPDGCharge();
1545  G4double m1 = theParticle->GetPDGMass();
1546 
1547  fWaveVector = partMom/CLHEP::hbarc;
1548 
1549  G4double lambda = fCofLambda*fWaveVector*fNuclearRadius;
1550  G4cout<<"kR = "<<lambda<<G4endl;
1551 
1552  if( z )
1553  {
1554  a = partMom/m1; // beta*gamma for m1
1555  fBeta = a/std::sqrt(1+a*a);
1556  fZommerfeld = CalculateZommerfeld( fBeta, z, fAtomicNumber);
1557  fRutherfordRatio = fZommerfeld/fWaveVector;
1558  fAm = CalculateAm( partMom, fZommerfeld, fAtomicNumber);
1559  }
1560  G4cout<<"fZommerfeld = "<<fZommerfeld<<G4endl;
1561  fProfileLambda = lambda; // *std::sqrt(1.-2*fZommerfeld/lambda);
1562  G4cout<<"fProfileLambda = "<<fProfileLambda<<G4endl;
1563  fProfileDelta = fCofDelta*fProfileLambda;
1564  fProfileAlpha = fCofAlpha*fProfileLambda;
1565 
1566  CalculateCoulombPhaseZero();
1567  CalculateRutherfordAnglePar();
1568 
1569  return;
1570 }
1572 //
1573 // Test for given particle and element table of momentum, angle probability.
1574 // For the partMom in CMS.
1575 
1576 inline void G4NuclNuclDiffuseElastic::InitDynParameters(const G4ParticleDefinition* theParticle,
1577  G4double partMom)
1578 {
1579  G4double a = 0.;
1580  G4double z = theParticle->GetPDGCharge();
1581  G4double m1 = theParticle->GetPDGMass();
1582 
1583  fWaveVector = partMom/CLHEP::hbarc;
1584 
1585  G4double lambda = fCofLambda*fWaveVector*fNuclearRadius;
1586 
1587  if( z )
1588  {
1589  a = partMom/m1; // beta*gamma for m1
1590  fBeta = a/std::sqrt(1+a*a);
1591  fZommerfeld = CalculateZommerfeld( fBeta, z, fAtomicNumber);
1592  fRutherfordRatio = fZommerfeld/fWaveVector;
1593  fAm = CalculateAm( partMom, fZommerfeld, fAtomicNumber);
1594  }
1595  fProfileLambda = lambda; // *std::sqrt(1.-2*fZommerfeld/lambda);
1596  fProfileDelta = fCofDelta*fProfileLambda;
1597  fProfileAlpha = fCofAlpha*fProfileLambda;
1598 
1599  CalculateCoulombPhaseZero();
1600  CalculateRutherfordAnglePar();
1601 
1602  return;
1603 }
1604 
1605 
1607 //
1608 // Test for given particle and element table of momentum, angle probability.
1609 // For the partMom in CMS.
1610 
1611 inline void G4NuclNuclDiffuseElastic::InitParametersGla(const G4DynamicParticle* aParticle,
1612  G4double partMom, G4double Z, G4double A)
1613 {
1614  fAtomicNumber = Z; // target atomic number
1615  fAtomicWeight = A; // target number of nucleons
1616 
1617  fNuclearRadius2 = CalculateNuclearRad(fAtomicWeight); // target nucleus radius
1618  G4double A1 = G4double( aParticle->GetDefinition()->GetBaryonNumber() );
1619  fNuclearRadius1 = CalculateNuclearRad(A1); // projectile nucleus radius
1620  fNuclearRadiusSquare = fNuclearRadius1*fNuclearRadius1+fNuclearRadius2*fNuclearRadius2;
1621 
1622 
1623  G4double a = 0., kR12;
1624  G4double z = aParticle->GetDefinition()->GetPDGCharge();
1625  G4double m1 = aParticle->GetDefinition()->GetPDGMass();
1626 
1627  fWaveVector = partMom/CLHEP::hbarc;
1628 
1629  G4double pN = A1 - z;
1630  if( pN < 0. ) pN = 0.;
1631 
1632  G4double tN = A - Z;
1633  if( tN < 0. ) tN = 0.;
1634 
1635  G4double pTkin = aParticle->GetKineticEnergy();
1636  pTkin /= A1;
1637 
1638 
1639  fSumSigma = (Z*z+pN*tN)*GetHadronNucleonXscNS(theProton, pTkin, theProton) +
1640  (z*tN+pN*Z)*GetHadronNucleonXscNS(theProton, pTkin, theNeutron);
1641 
1642  G4cout<<"fSumSigma = "<<fSumSigma/CLHEP::millibarn<<" mb"<<G4endl;
1643  G4cout<<"pi*R2 = "<<CLHEP::pi*fNuclearRadiusSquare/CLHEP::millibarn<<" mb"<<G4endl;
1644  kR12 = fWaveVector*std::sqrt(fNuclearRadiusSquare);
1645  G4cout<<"k*sqrt(R2) = "<<kR12<<" "<<G4endl;
1646  fMaxL = (G4int(kR12)+1)*4;
1647  G4cout<<"fMaxL = "<<fMaxL<<" "<<G4endl;
1648 
1649  if( z )
1650  {
1651  a = partMom/m1; // beta*gamma for m1
1652  fBeta = a/std::sqrt(1+a*a);
1653  fZommerfeld = CalculateZommerfeld( fBeta, z, fAtomicNumber);
1654  fAm = CalculateAm( partMom, fZommerfeld, fAtomicNumber);
1655  }
1656 
1657  CalculateCoulombPhaseZero();
1658 
1659 
1660  return;
1661 }
1662 
1663 
1665 //
1666 // Returns nucleon-nucleon cross-section based on N. Starkov parametrisation of
1667 // data from mainly http://wwwppds.ihep.su:8001/c5-6A.html database
1668 // projectile nucleon is pParticle with pTkin shooting target nucleon tParticle
1669 
1670 inline G4double
1671 G4NuclNuclDiffuseElastic::GetHadronNucleonXscNS( G4ParticleDefinition* pParticle,
1672  G4double pTkin,
1673  G4ParticleDefinition* tParticle)
1674 {
1675  G4double xsection(0), /*Delta,*/ A0, B0;
1676  G4double hpXsc(0);
1677  G4double hnXsc(0);
1678 
1679 
1680  G4double targ_mass = tParticle->GetPDGMass();
1681  G4double proj_mass = pParticle->GetPDGMass();
1682 
1683  G4double proj_energy = proj_mass + pTkin;
1684  G4double proj_momentum = std::sqrt(pTkin*(pTkin+2*proj_mass));
1685 
1686  G4double sMand = CalcMandelstamS ( proj_mass , targ_mass , proj_momentum );
1687 
1688  sMand /= CLHEP::GeV*CLHEP::GeV; // in GeV for parametrisation
1689  proj_momentum /= CLHEP::GeV;
1690  proj_energy /= CLHEP::GeV;
1691  proj_mass /= CLHEP::GeV;
1692  G4double logS = std::log(sMand);
1693 
1694  // General PDG fit constants
1695 
1696 
1697  // fEtaRatio=Re[f(0)]/Im[f(0)]
1698 
1699  if( proj_momentum >= 1.2 )
1700  {
1701  fEtaRatio = 0.13*(logS - 5.8579332)*std::pow(sMand,-0.18);
1702  }
1703  else if( proj_momentum >= 0.6 )
1704  {
1705  fEtaRatio = -75.5*(std::pow(proj_momentum,0.25)-0.95)/
1706  (std::pow(3*proj_momentum,2.2)+1);
1707  }
1708  else
1709  {
1710  fEtaRatio = 15.5*proj_momentum/(27*proj_momentum*proj_momentum*proj_momentum+2);
1711  }
1712  G4cout<<"fEtaRatio = "<<fEtaRatio<<G4endl;
1713 
1714  // xsc
1715 
1716  if( proj_momentum >= 10. ) // high energy: pp = nn = np
1717  // if( proj_momentum >= 2.)
1718  {
1719  //Delta = 1.;
1720 
1721  //if( proj_energy < 40. ) Delta = 0.916+0.0021*proj_energy;
1722 
1723  if( proj_momentum >= 10.)
1724  {
1725  B0 = 7.5;
1726  A0 = 100. - B0*std::log(3.0e7);
1727 
1728  xsection = A0 + B0*std::log(proj_energy) - 11
1729  + 103*std::pow(2*0.93827*proj_energy + proj_mass*proj_mass+
1730  0.93827*0.93827,-0.165); // mb
1731  }
1732  }
1733  else // low energy pp = nn != np
1734  {
1735  if(pParticle == tParticle) // pp or nn // nn to be pp
1736  {
1737  if( proj_momentum < 0.73 )
1738  {
1739  hnXsc = 23 + 50*( std::pow( std::log(0.73/proj_momentum), 3.5 ) );
1740  }
1741  else if( proj_momentum < 1.05 )
1742  {
1743  hnXsc = 23 + 40*(std::log(proj_momentum/0.73))*
1744  (std::log(proj_momentum/0.73));
1745  }
1746  else // if( proj_momentum < 10. )
1747  {
1748  hnXsc = 39.0 +
1749  75*(proj_momentum - 1.2)/(std::pow(proj_momentum,3.0) + 0.15);
1750  }
1751  xsection = hnXsc;
1752  }
1753  else // pn to be np
1754  {
1755  if( proj_momentum < 0.8 )
1756  {
1757  hpXsc = 33+30*std::pow(std::log(proj_momentum/1.3),4.0);
1758  }
1759  else if( proj_momentum < 1.4 )
1760  {
1761  hpXsc = 33+30*std::pow(std::log(proj_momentum/0.95),2.0);
1762  }
1763  else // if( proj_momentum < 10. )
1764  {
1765  hpXsc = 33.3+
1766  20.8*(std::pow(proj_momentum,2.0)-1.35)/
1767  (std::pow(proj_momentum,2.50)+0.95);
1768  }
1769  xsection = hpXsc;
1770  }
1771  }
1772  xsection *= CLHEP::millibarn; // parametrised in mb
1773  G4cout<<"xsection = "<<xsection/CLHEP::millibarn<<" mb"<<G4endl;
1774  return xsection;
1775 }
1776 
1778 //
1779 //
1780 
1781 inline G4double G4NuclNuclDiffuseElastic::CalcMandelstamS( const G4double mp ,
1782  const G4double mt ,
1783  const G4double Plab )
1784 {
1785  G4double Elab = std::sqrt ( mp * mp + Plab * Plab );
1786  G4double sMand = mp*mp + mt*mt + 2*Elab*mt ;
1787 
1788  return sMand;
1789 }
1790 
1791 
1792 
1793 #endif
Generated on Sat May 25 2013 14:33:51 for Geant4 by   doxygen 1.8.4