source: trunk/source/processes/hadronic/models/coherent_elastic/include/G4NuclNuclDiffuseElastic.hh @ 1337

Last change on this file since 1337 was 1337, checked in by garnier, 14 years ago

tag geant4.9.4 beta 1 + modifs locales
File size: 31.4 KB
Line 
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: G4NuclNuclDiffuseElastic.hh,v 1.8 2009/04/10 13:22:25 grichine Exp $
28// GEANT4 tag $Name: geant4-09-04-beta-01 $
29//
30//
31// G4 Model: optical elastic scattering with 4-momentum balance
32//
33// Class Description
34// Final state production model for nucleus-nucleus elastic scattering;
35// Coulomb amplitude is not considered as correction
36// (as in G4DiffuseElastic)
37// Class Description - End
38//
39//
40// 17.03.09 V. Grichine implementation for Coulomb elastic scattering
41
42
43#ifndef G4NuclNuclDiffuseElastic_h
44#define G4NuclNuclDiffuseElastic_h 1
45 
46#include "globals.hh"
47#include <complex>
48#include "G4Integrator.hh"
49
50#include "G4HadronicInteraction.hh"
51#include "G4HadProjectile.hh"
52#include "G4Nucleus.hh"
53
54using namespace std;
55
56class G4ParticleDefinition;
57class G4PhysicsTable;
58class G4PhysicsLogVector;
59
60class G4NuclNuclDiffuseElastic : public G4HadronicInteraction
61{
62public:
63
64  G4NuclNuclDiffuseElastic();
65
66  G4NuclNuclDiffuseElastic(const G4ParticleDefinition* aParticle);
67
68
69
70
71
72  virtual ~G4NuclNuclDiffuseElastic();
73
74  void Initialise();
75
76  void InitialiseOnFly(G4double Z, G4double A);
77
78  void BuildAngleTable();
79
80 
81  G4HadFinalState * ApplyYourself(const G4HadProjectile & aTrack, 
82                                  G4Nucleus & targetNucleus);
83
84
85  void SetPlabLowLimit(G4double value);
86
87  void SetHEModelLowLimit(G4double value);
88
89  void SetQModelLowLimit(G4double value);
90
91  void SetLowestEnergyLimit(G4double value);
92
93  void SetRecoilKinEnergyLimit(G4double value);
94
95  G4double SampleT(const G4ParticleDefinition* aParticle, 
96                         G4double p, G4double A);
97
98  G4double SampleTableT(const G4ParticleDefinition* aParticle, 
99                         G4double p, G4double Z, G4double A);
100
101  G4double SampleThetaCMS(const G4ParticleDefinition* aParticle, G4double p, G4double A);
102
103  G4double SampleTableThetaCMS(const G4ParticleDefinition* aParticle, G4double p, 
104                                     G4double Z, G4double A);
105
106  G4double GetScatteringAngle(G4int iMomentum, G4int iAngle, G4double position);
107
108  G4double SampleThetaLab(const G4HadProjectile* aParticle, 
109                                G4double tmass, G4double A);
110
111  G4double GetDiffuseElasticXsc( const G4ParticleDefinition* particle, 
112                                 G4double theta, 
113                                 G4double momentum, 
114                                 G4double A         );
115
116  G4double GetInvElasticXsc( const G4ParticleDefinition* particle, 
117                                 G4double theta, 
118                                 G4double momentum, 
119                                 G4double A, G4double Z );
120
121  G4double GetDiffuseElasticSumXsc( const G4ParticleDefinition* particle, 
122                                 G4double theta, 
123                                 G4double momentum, 
124                                 G4double A, G4double Z );
125
126  G4double GetInvElasticSumXsc( const G4ParticleDefinition* particle, 
127                                 G4double tMand, 
128                                 G4double momentum, 
129                                 G4double A, G4double Z );
130
131  G4double IntegralElasticProb( const G4ParticleDefinition* particle, 
132                                 G4double theta, 
133                                 G4double momentum, 
134                                 G4double A            );
135 
136
137  G4double GetCoulombElasticXsc( const G4ParticleDefinition* particle, 
138                                 G4double theta, 
139                                 G4double momentum, 
140                                 G4double Z         );
141
142  G4double GetInvCoulombElasticXsc( const G4ParticleDefinition* particle, 
143                                 G4double tMand, 
144                                 G4double momentum, 
145                                 G4double A, G4double Z         );
146
147  G4double GetCoulombTotalXsc( const G4ParticleDefinition* particle, 
148                                 G4double momentum, G4double Z       );
149
150  G4double GetCoulombIntegralXsc( const G4ParticleDefinition* particle, 
151                                 G4double momentum, G4double Z, 
152                                 G4double theta1, G4double theta2         );
153
154
155  G4double CalculateParticleBeta( const G4ParticleDefinition* particle, 
156                                        G4double momentum    );
157
158  G4double CalculateZommerfeld( G4double beta, G4double Z1, G4double Z2 );
159
160  G4double CalculateAm( G4double momentum, G4double n, G4double Z);
161
162  G4double CalculateNuclearRad( G4double A);
163
164  G4double ThetaCMStoThetaLab(const G4DynamicParticle* aParticle, 
165                                G4double tmass, G4double thetaCMS);
166
167  G4double ThetaLabToThetaCMS(const G4DynamicParticle* aParticle, 
168                                G4double tmass, G4double thetaLab);
169
170  void TestAngleTable(const G4ParticleDefinition* theParticle, G4double partMom,
171                      G4double Z, G4double A);
172
173
174
175  G4double BesselJzero(G4double z);
176  G4double BesselJone(G4double z);
177  G4double DampFactor(G4double z);
178  G4double BesselOneByArg(G4double z);
179
180  G4double GetDiffElasticProb(G4double theta);
181  G4double GetDiffElasticSumProb(G4double theta);
182  G4double GetDiffElasticSumProbA(G4double alpha);
183  G4double GetIntegrandFunction(G4double theta);
184
185  G4double GetNuclearRadius(){return fNuclearRadius;};
186
187
188  // Technical math functions for strong Coulomb contribution
189
190  G4complex GammaLogarithm(G4complex xx);
191
192  G4double  GetErf(G4double x);
193
194  G4complex GetErfcComp(G4complex z, G4int nMax);
195  G4complex GetErfcSer(G4complex z, G4int nMax);
196  G4complex GetErfcInt(G4complex z); // , G4int nMax);
197
198  G4complex GetErfComp(G4complex z, G4int nMax);  // AandS algorithm != Ser, Int
199  G4complex GetErfSer(G4complex z, G4int nMax);
200
201  G4double GetExpCos(G4double x);
202  G4double GetExpSin(G4double x);
203  G4complex GetErfInt(G4complex z); // , G4int nMax);
204
205
206  G4complex TestErfcComp(G4complex z, G4int nMax);
207  G4complex TestErfcSer(G4complex z, G4int nMax);
208  G4complex TestErfcInt(G4complex z); // , G4int nMax);
209
210  G4complex CoulombAmplitude(G4double theta);
211  void CalculateCoulombPhaseZero();
212  void CalculateRutherfordAnglePar();
213
214  G4double ProfileNear(G4double theta);
215  G4double ProfileFar(G4double theta);
216
217  G4complex PhaseNear(G4double theta);
218  G4complex PhaseFar(G4double theta);
219
220  G4complex GammaLess(G4double theta);
221  G4complex GammaMore(G4double theta);
222
223  G4complex AmplitudeNear(G4double theta);
224  G4complex AmplitudeFar(G4double theta);
225  G4complex Amplitude(G4double theta);
226  G4double  AmplitudeMod2(G4double theta);
227  void      InitParameters(const G4ParticleDefinition* theParticle, 
228                              G4double partMom, G4double Z, G4double A); 
229
230  G4double GetProfileLambda(){return fProfileLambda;};
231
232  void SetProfileLambda(G4double pl) {fProfileLambda = pl;};
233  void SetProfileDelta(G4double pd) {fProfileDelta = pd;};
234  void SetProfileAlpha(G4double pa){fProfileAlpha = pa;};
235
236private:
237
238
239  G4ParticleDefinition* theProton;
240  G4ParticleDefinition* theNeutron;
241  G4ParticleDefinition* theDeuteron;
242  G4ParticleDefinition* theAlpha;
243
244  const G4ParticleDefinition* thePionPlus;
245  const G4ParticleDefinition* thePionMinus;
246
247  G4double lowEnergyRecoilLimit; 
248  G4double lowEnergyLimitHE; 
249  G4double lowEnergyLimitQ; 
250  G4double lowestEnergyLimit; 
251  G4double plabLowLimit;
252
253  G4int fEnergyBin;
254  G4int fAngleBin;
255
256  G4PhysicsLogVector*           fEnergyVector;
257  G4PhysicsTable*               fAngleTable;
258  std::vector<G4PhysicsTable*>  fAngleBank;
259
260  std::vector<G4double> fElementNumberVector;
261  std::vector<G4String> fElementNameVector;
262
263  const G4ParticleDefinition* fParticle;
264  G4double fWaveVector;
265  G4double fAtomicWeight;
266  G4double fAtomicNumber;
267
268  G4double fNuclearRadius1;
269  G4double fNuclearRadius2;
270  G4double fNuclearRadius;
271
272  G4double fBeta;
273  G4double fZommerfeld;
274  G4double fAm;
275  G4bool   fAddCoulomb;
276
277  G4double fCoulombPhase0;
278  G4double fHalfRutThetaTg;
279  G4double fRutherfordTheta;
280
281  G4double fProfileLambda;
282  G4double fProfileDelta;
283  G4double fProfileAlpha;
284
285  G4double fReZ;
286
287};
288
289
290inline void G4NuclNuclDiffuseElastic::SetRecoilKinEnergyLimit(G4double value)
291{
292  lowEnergyRecoilLimit = value;
293}
294
295inline void G4NuclNuclDiffuseElastic::SetPlabLowLimit(G4double value)
296{
297  plabLowLimit = value;
298}
299
300inline void G4NuclNuclDiffuseElastic::SetHEModelLowLimit(G4double value)
301{
302  lowEnergyLimitHE = value;
303}
304
305inline void G4NuclNuclDiffuseElastic::SetQModelLowLimit(G4double value)
306{
307  lowEnergyLimitQ = value;
308}
309
310inline void G4NuclNuclDiffuseElastic::SetLowestEnergyLimit(G4double value)
311{
312  lowestEnergyLimit = value;
313}
314
315
316/////////////////////////////////////////////////////////////
317//
318// Bessel J0 function based on rational approximation from
319// J.F. Hart, Computer Approximations, New York, Willey 1968, p. 141
320
321inline G4double G4NuclNuclDiffuseElastic::BesselJzero(G4double value)
322{
323  G4double modvalue, value2, fact1, fact2, arg, shift, bessel;
324
325  modvalue = fabs(value);
326
327  if ( value < 8.0 && value > -8.0 )
328  {
329    value2 = value*value;
330
331    fact1  = 57568490574.0 + value2*(-13362590354.0 
332                           + value2*( 651619640.7 
333                           + value2*(-11214424.18 
334                           + value2*( 77392.33017 
335                           + value2*(-184.9052456   ) ) ) ) );
336
337    fact2  = 57568490411.0 + value2*( 1029532985.0 
338                           + value2*( 9494680.718
339                           + value2*(59272.64853
340                           + value2*(267.8532712 
341                           + value2*1.0               ) ) ) );
342
343    bessel = fact1/fact2;
344  } 
345  else 
346  {
347    arg    = 8.0/modvalue;
348
349    value2 = arg*arg;
350
351    shift  = modvalue-0.785398164;
352
353    fact1  = 1.0 + value2*(-0.1098628627e-2 
354                 + value2*(0.2734510407e-4
355                 + value2*(-0.2073370639e-5 
356                 + value2*0.2093887211e-6    ) ) );
357
358    fact2  = -0.1562499995e-1 + value2*(0.1430488765e-3
359                              + value2*(-0.6911147651e-5 
360                              + value2*(0.7621095161e-6
361                              - value2*0.934945152e-7    ) ) );
362
363    bessel = sqrt(0.636619772/modvalue)*(cos(shift)*fact1 - arg*sin(shift)*fact2 );
364  }
365  return bessel;
366}
367
368/////////////////////////////////////////////////////////////
369//
370// Bessel J1 function based on rational approximation from
371// J.F. Hart, Computer Approximations, New York, Willey 1968, p. 141
372
373inline G4double G4NuclNuclDiffuseElastic::BesselJone(G4double value)
374{
375  G4double modvalue, value2, fact1, fact2, arg, shift, bessel;
376
377  modvalue = fabs(value);
378
379  if ( modvalue < 8.0 ) 
380  {
381    value2 = value*value;
382
383    fact1  = value*(72362614232.0 + value2*(-7895059235.0 
384                                  + value2*( 242396853.1
385                                  + value2*(-2972611.439 
386                                  + value2*( 15704.48260 
387                                  + value2*(-30.16036606  ) ) ) ) ) );
388
389    fact2  = 144725228442.0 + value2*(2300535178.0 
390                            + value2*(18583304.74
391                            + value2*(99447.43394 
392                            + value2*(376.9991397 
393                            + value2*1.0             ) ) ) );
394    bessel = fact1/fact2;
395  } 
396  else 
397  {
398    arg    = 8.0/modvalue;
399
400    value2 = arg*arg;
401
402    shift  = modvalue - 2.356194491;
403
404    fact1  = 1.0 + value2*( 0.183105e-2 
405                 + value2*(-0.3516396496e-4
406                 + value2*(0.2457520174e-5 
407                 + value2*(-0.240337019e-6          ) ) ) );
408
409    fact2 = 0.04687499995 + value2*(-0.2002690873e-3
410                          + value2*( 0.8449199096e-5
411                          + value2*(-0.88228987e-6
412                          + value2*0.105787412e-6       ) ) );
413
414    bessel = sqrt( 0.636619772/modvalue)*(cos(shift)*fact1 - arg*sin(shift)*fact2);
415
416    if (value < 0.0) bessel = -bessel;
417  }
418  return bessel;
419}
420
421////////////////////////////////////////////////////////////////////
422//
423// damp factor in diffraction x/sh(x), x was already *pi
424
425inline G4double G4NuclNuclDiffuseElastic::DampFactor(G4double x)
426{
427  G4double df;
428  G4double f2 = 2., f3 = 6., f4 = 24.; // first factorials
429
430  // x *= pi;
431
432  if( std::fabs(x) < 0.01 )
433  { 
434    df = 1./(1. + x/f2 + x*x/f3 + x*x*x/f4);
435  }
436  else
437  {
438    df = x/std::sinh(x); 
439  }
440  return df;
441}
442
443
444////////////////////////////////////////////////////////////////////
445//
446// return J1(x)/x with special case for small x
447
448inline G4double G4NuclNuclDiffuseElastic::BesselOneByArg(G4double x)
449{
450  G4double x2, result;
451 
452  if( std::fabs(x) < 0.01 )
453  { 
454   x     *= 0.5;
455   x2     = x*x;
456   result = 2. - x2 + x2*x2/6.;
457  }
458  else
459  {
460    result = BesselJone(x)/x; 
461  }
462  return result;
463}
464
465////////////////////////////////////////////////////////////////////
466//
467// return particle beta
468
469inline  G4double G4NuclNuclDiffuseElastic::CalculateParticleBeta( const G4ParticleDefinition* particle, 
470                                        G4double momentum    )
471{
472  G4double mass = particle->GetPDGMass();
473  G4double a    = momentum/mass;
474  fBeta         = a/std::sqrt(1+a*a);
475
476  return fBeta; 
477}
478
479////////////////////////////////////////////////////////////////////
480//
481// return Zommerfeld parameter for Coulomb scattering
482
483inline  G4double G4NuclNuclDiffuseElastic::CalculateZommerfeld( G4double beta, G4double Z1, G4double Z2 )
484{
485  fZommerfeld = fine_structure_const*Z1*Z2/beta;
486
487  return fZommerfeld; 
488}
489
490////////////////////////////////////////////////////////////////////
491//
492// return Wentzel correction for Coulomb scattering
493
494inline  G4double G4NuclNuclDiffuseElastic::CalculateAm( G4double momentum, G4double n, G4double Z)
495{
496  G4double k   = momentum/hbarc;
497  G4double ch  = 1.13 + 3.76*n*n;
498  G4double zn  = 1.77*k*std::pow(Z,-1./3.)*Bohr_radius;
499  G4double zn2 = zn*zn;
500  fAm          = ch/zn2;
501
502  return fAm;
503}
504
505////////////////////////////////////////////////////////////////////
506//
507// calculate nuclear radius for different atomic weights using different approximations
508
509inline  G4double G4NuclNuclDiffuseElastic::CalculateNuclearRad( G4double A)
510{
511  G4double r0, radius;
512
513  if( A < 50. )
514  {
515    if( A > 10. ) r0  = 1.16*( 1 - std::pow(A, -2./3.) )*fermi;   // 1.08*fermi;
516    else          r0  = 1.1*fermi;
517
518    radius = r0*std::pow(A, 1./3.);
519  }
520  else
521  {
522    r0 = 1.7*fermi;   // 1.7*fermi;
523
524    radius = r0*std::pow(A, 0.27); // 0.27);
525  }
526  return radius;
527}
528
529////////////////////////////////////////////////////////////////////
530//
531// return Coulomb scattering differential xsc with Wentzel correction 
532
533inline  G4double G4NuclNuclDiffuseElastic::GetCoulombElasticXsc( const G4ParticleDefinition* particle, 
534                                 G4double theta, 
535                                 G4double momentum, 
536                                 G4double Z         )
537{
538  G4double sinHalfTheta  = std::sin(0.5*theta);
539  G4double sinHalfTheta2 = sinHalfTheta*sinHalfTheta;
540  G4double beta          = CalculateParticleBeta( particle, momentum);
541  G4double z             = particle->GetPDGCharge();
542  G4double n             = CalculateZommerfeld( beta, z, Z );
543  G4double am            = CalculateAm( momentum, n, Z);
544  G4double k             = momentum/hbarc;
545  G4double ch            = 0.5*n/k;
546  G4double ch2           = ch*ch;
547  G4double xsc           = ch2/(sinHalfTheta2+am)/(sinHalfTheta2+am);
548
549  return xsc;
550}
551
552
553////////////////////////////////////////////////////////////////////
554//
555// return Coulomb scattering total xsc with Wentzel correction 
556
557inline  G4double G4NuclNuclDiffuseElastic::GetCoulombTotalXsc( const G4ParticleDefinition* particle, 
558                                                             G4double momentum, G4double Z  )
559{
560  G4double beta          = CalculateParticleBeta( particle, momentum);
561  G4cout<<"beta = "<<beta<<G4endl;
562  G4double z             = particle->GetPDGCharge();
563  G4double n             = CalculateZommerfeld( beta, z, Z );
564  G4cout<<"fZomerfeld = "<<n<<G4endl;
565  G4double am            = CalculateAm( momentum, n, Z);
566  G4cout<<"cof Am = "<<am<<G4endl;
567  G4double k             = momentum/hbarc;
568  G4cout<<"k = "<<k*fermi<<" 1/fermi"<<G4endl;
569  G4cout<<"k*Bohr_radius = "<<k*Bohr_radius<<G4endl;
570  G4double ch            = n/k;
571  G4double ch2           = ch*ch;
572  G4double xsc           = ch2*pi/(am +am*am);
573
574  return xsc;
575}
576
577////////////////////////////////////////////////////////////////////
578//
579// return Coulomb scattering xsc with Wentzel correction  integrated between
580// theta1 and < theta2
581
582inline  G4double G4NuclNuclDiffuseElastic::GetCoulombIntegralXsc( const G4ParticleDefinition* particle, 
583                                 G4double momentum, G4double Z, 
584                                 G4double theta1, G4double theta2 )
585{
586  G4double c1 = std::cos(theta1);
587  G4cout<<"c1 = "<<c1<<G4endl;
588  G4double c2 = std::cos(theta2);
589  G4cout<<"c2 = "<<c2<<G4endl;
590  G4double beta          = CalculateParticleBeta( particle, momentum);
591  // G4cout<<"beta = "<<beta<<G4endl;
592  G4double z             = particle->GetPDGCharge();
593  G4double n             = CalculateZommerfeld( beta, z, Z );
594  // G4cout<<"fZomerfeld = "<<n<<G4endl;
595  G4double am            = CalculateAm( momentum, n, Z);
596  // G4cout<<"cof Am = "<<am<<G4endl;
597  G4double k             = momentum/hbarc;
598  // G4cout<<"k = "<<k*fermi<<" 1/fermi"<<G4endl;
599  // G4cout<<"k*Bohr_radius = "<<k*Bohr_radius<<G4endl;
600  G4double ch            = n/k;
601  G4double ch2           = ch*ch;
602  am *= 2.;
603  G4double xsc           = ch2*twopi*(c1-c2);
604           xsc          /= (1 - c1 + am)*(1 - c2 + am);
605
606  return xsc;
607}
608
609///////////////////////////////////////////////////////////////////
610//
611// For the calculation of arg Gamma(z) one needs complex extension
612// of ln(Gamma(z))
613
614inline G4complex G4NuclNuclDiffuseElastic::GammaLogarithm(G4complex zz)
615{
616  static G4double cof[6] = { 76.18009172947146,     -86.50532032941677,
617                             24.01409824083091,      -1.231739572450155,
618                              0.1208650973866179e-2, -0.5395239384953e-5  } ;
619  register G4int j;
620  G4complex z = zz - 1.0;
621  G4complex tmp = z + 5.5;
622  tmp -= (z + 0.5) * std::log(tmp);
623  G4complex ser = G4complex(1.000000000190015,0.);
624
625  for ( j = 0; j <= 5; j++ )
626  {
627    z += 1.0;
628    ser += cof[j]/z;
629  }
630  return -tmp + std::log(2.5066282746310005*ser);
631}
632
633/////////////////////////////////////////////////////////////////
634//
635//
636
637inline G4double  G4NuclNuclDiffuseElastic::GetErf(G4double x)
638{
639  G4double t, z, tmp, result;
640
641  z   = std::fabs(x);
642  t   = 1.0/(1.0+0.5*z);
643
644  tmp = t*exp(-z*z-1.26551223+t*(1.00002368+t*(0.37409196+t*(0.09678418+
645                t*(-0.18628806+t*(0.27886807+t*(-1.13520398+t*(1.48851587+
646                t*(-0.82215223+t*0.17087277)))))))));
647
648  if( x >= 0.) result = 1. - tmp;
649  else         result = 1. + tmp; 
650   
651  return result;
652}
653
654/////////////////////////////////////////////////////////////////
655//
656//
657
658inline G4complex G4NuclNuclDiffuseElastic::GetErfcComp(G4complex z, G4int nMax)
659{
660  G4complex erfcz = 1. - GetErfComp( z, nMax);
661  return erfcz;
662}
663
664/////////////////////////////////////////////////////////////////
665//
666//
667
668inline G4complex G4NuclNuclDiffuseElastic::GetErfcSer(G4complex z, G4int nMax)
669{
670  G4complex erfcz = 1. - GetErfSer( z, nMax);
671  return erfcz;
672}
673
674/////////////////////////////////////////////////////////////////
675//
676//
677
678inline G4complex G4NuclNuclDiffuseElastic::GetErfcInt(G4complex z) // , G4int nMax)
679{
680  G4complex erfcz = 1. - GetErfInt( z); // , nMax);
681  return erfcz;
682}
683
684/////////////////////////////////////////////////////////////////
685//
686//
687
688inline G4complex G4NuclNuclDiffuseElastic::TestErfcComp(G4complex z, G4int nMax)
689{
690  G4complex miz = G4complex( z.imag(), -z.real() ); 
691  G4complex erfcz = 1. - GetErfComp( miz, nMax);
692  G4complex w = std::exp(-z*z)*erfcz;
693  return w;
694}
695
696/////////////////////////////////////////////////////////////////
697//
698//
699
700inline G4complex G4NuclNuclDiffuseElastic::TestErfcSer(G4complex z, G4int nMax)
701{
702  G4complex miz = G4complex( z.imag(), -z.real() ); 
703  G4complex erfcz = 1. - GetErfSer( miz, nMax);
704  G4complex w = std::exp(-z*z)*erfcz;
705  return w;
706}
707
708/////////////////////////////////////////////////////////////////
709//
710//
711
712inline G4complex G4NuclNuclDiffuseElastic::TestErfcInt(G4complex z) // , G4int nMax)
713{
714  G4complex miz = G4complex( z.imag(), -z.real() ); 
715  G4complex erfcz = 1. - GetErfInt( miz); // , nMax);
716  G4complex w = std::exp(-z*z)*erfcz;
717  return w;
718}
719
720/////////////////////////////////////////////////////////////////
721//
722//
723
724inline G4complex G4NuclNuclDiffuseElastic::GetErfComp(G4complex z, G4int nMax)
725{
726  G4int n;
727  G4double n2, cofn, shny, chny, fn, gn;
728
729  G4double x = z.real();
730  G4double y = z.imag();
731
732  G4double outRe = 0., outIm = 0.;
733
734  G4double twox  = 2.*x;
735  G4double twoxy = twox*y;
736  G4double twox2 = twox*twox;
737
738  G4double cof1 = std::exp(-x*x)/pi;
739
740  G4double cos2xy = std::cos(twoxy);
741  G4double sin2xy = std::sin(twoxy);
742
743  G4double twoxcos2xy = twox*cos2xy;
744  G4double twoxsin2xy = twox*sin2xy;
745
746  for( n = 1; n <= nMax; n++)
747  {
748    n2   = n*n;
749
750    cofn = std::exp(-0.5*n2)/(n2+twox2);  // /(n2+0.5*twox2);
751
752    chny = std::cosh(n*y);
753    shny = std::sinh(n*y);
754
755    fn  = twox - twoxcos2xy*chny + n*sin2xy*shny;
756    gn  =        twoxsin2xy*chny + n*cos2xy*shny;
757
758    fn *= cofn;
759    gn *= cofn;
760
761    outRe += fn;
762    outIm += gn;
763  }
764  outRe *= 2*cof1;
765  outIm *= 2*cof1;
766
767  if(std::abs(x) < 0.0001)
768  {
769    outRe += GetErf(x);
770    outIm += cof1*y;
771  }
772  else
773  {
774    outRe += GetErf(x) + cof1*(1-cos2xy)/twox;
775    outIm += cof1*sin2xy/twox;
776  }
777  return G4complex(outRe, outIm);
778}
779
780/////////////////////////////////////////////////////////////////
781//
782//
783
784inline G4complex G4NuclNuclDiffuseElastic::GetErfSer(G4complex z, G4int nMax)
785{
786  G4int n;
787  G4double a =1., b = 1., tmp;
788  G4complex sum = z, d = z;
789
790  for( n = 1; n <= nMax; n++)
791  {
792    a *= 2.;
793    b *= 2.*n +1.;
794    d *= z*z;
795
796    tmp = a/b;
797
798    sum += tmp*d;
799  }
800  sum *= 2.*std::exp(-z*z)/std::sqrt(pi);
801
802  return sum;
803}
804
805/////////////////////////////////////////////////////////////////////
806
807inline  G4double  G4NuclNuclDiffuseElastic::GetExpCos(G4double x)
808{
809  G4double result;
810
811  result = std::exp(x*x-fReZ*fReZ);
812  result *= std::cos(2.*x*fReZ);
813  return result;
814}
815
816/////////////////////////////////////////////////////////////////////
817
818inline  G4double  G4NuclNuclDiffuseElastic::GetExpSin(G4double x)
819{
820  G4double result;
821
822  result = std::exp(x*x-fReZ*fReZ);
823  result *= std::sin(2.*x*fReZ);
824  return result;
825}
826
827
828
829/////////////////////////////////////////////////////////////////
830//
831//
832
833inline G4complex G4NuclNuclDiffuseElastic::GetErfInt(G4complex z) // , G4int nMax)
834{
835  G4double outRe, outIm;
836
837  G4double x = z.real();
838  G4double y = z.imag();
839  fReZ       = x;
840
841  G4Integrator<G4NuclNuclDiffuseElastic,G4double(G4NuclNuclDiffuseElastic::*)(G4double)> integral;
842
843  outRe = integral.Legendre96(this,&G4NuclNuclDiffuseElastic::GetExpSin, 0., y );
844  outIm = integral.Legendre96(this,&G4NuclNuclDiffuseElastic::GetExpCos, 0., y );
845
846  outRe *= 2./sqrt(pi);
847  outIm *= 2./sqrt(pi);
848
849  outRe += GetErf(x);
850
851  return G4complex(outRe, outIm);
852}
853
854
855/////////////////////////////////////////////////////////////////
856//
857//
858
859inline  G4complex G4NuclNuclDiffuseElastic::CoulombAmplitude(G4double theta)
860{
861  G4complex ca;
862 
863  G4double sinHalfTheta  = std::sin(0.5*theta);
864  G4double sinHalfTheta2 = sinHalfTheta*sinHalfTheta; 
865  sinHalfTheta2         += fAm;
866  G4double order         = 2.*fCoulombPhase0 - fZommerfeld*std::log(sinHalfTheta2);
867  G4complex z            = G4complex(0., order);
868  ca                     = std::exp(z);
869  ca                    *= -fZommerfeld/(2.*fWaveVector*sinHalfTheta2);
870
871  return ca; 
872}
873
874/////////////////////////////////////////////////////////////////
875//
876//
877
878
879inline  void G4NuclNuclDiffuseElastic::CalculateCoulombPhaseZero()
880{
881  G4complex z        = G4complex(1,fZommerfeld); 
882  G4complex gammalog = GammaLogarithm(z);
883  fCoulombPhase0      = gammalog.imag();
884}
885
886
887/////////////////////////////////////////////////////////////////
888//
889//
890
891
892inline  void G4NuclNuclDiffuseElastic::CalculateRutherfordAnglePar()
893{
894  fHalfRutThetaTg  = fZommerfeld/(fWaveVector*fNuclearRadius);
895  fRutherfordTheta = 2.*std::atan(fHalfRutThetaTg);
896  G4cout<<"fRutherfordTheta = "<<fRutherfordTheta/degree<<" degree"<<G4endl;
897
898}
899
900/////////////////////////////////////////////////////////////////
901//
902//
903
904inline   G4double G4NuclNuclDiffuseElastic::ProfileNear(G4double theta)
905{
906  G4double dTheta = fRutherfordTheta - theta;
907  G4double result = 0., argument = 0.;
908
909  if(std::abs(dTheta) < 0.001) result = fProfileAlpha*fProfileDelta;
910  else
911  {
912    argument = fProfileDelta*dTheta;
913    result   = pi*argument*std::exp(fProfileAlpha*argument);
914    result  /= std::sinh(pi*argument);
915    result  -= 1.;
916    result  /= dTheta;
917  }
918  return result;
919}
920
921/////////////////////////////////////////////////////////////////
922//
923//
924
925inline   G4double G4NuclNuclDiffuseElastic::ProfileFar(G4double theta)
926{
927  G4double dTheta   = fRutherfordTheta + theta;
928  G4double argument = fProfileDelta*dTheta;
929
930  G4double result   = pi*argument*std::exp(fProfileAlpha*argument);
931  result           /= std::sinh(pi*argument);
932  result           /= dTheta;
933
934  return result;
935}
936
937/////////////////////////////////////////////////////////////////
938//
939//
940
941inline   G4complex G4NuclNuclDiffuseElastic::PhaseNear(G4double theta)
942{
943  G4double twosigma = 2.*fCoulombPhase0; 
944  twosigma -= fZommerfeld*std::log(fHalfRutThetaTg/(1.+fHalfRutThetaTg*fHalfRutThetaTg));
945  twosigma += fRutherfordTheta*fZommerfeld/fHalfRutThetaTg - halfpi;
946  twosigma -= fProfileLambda*theta - 0.25*pi;
947
948  G4complex z = G4complex(0., twosigma);
949
950  return std::exp(z);
951}
952
953/////////////////////////////////////////////////////////////////
954//
955//
956
957inline   G4complex G4NuclNuclDiffuseElastic::PhaseFar(G4double theta)
958{
959  G4double twosigma = 2.*fCoulombPhase0; 
960  twosigma -= fZommerfeld*std::log(fHalfRutThetaTg/(1.+fHalfRutThetaTg*fHalfRutThetaTg));
961  twosigma += fRutherfordTheta*fZommerfeld/fHalfRutThetaTg - halfpi;
962  twosigma += fProfileLambda*theta - 0.25*pi;
963
964  G4complex z = G4complex(0., twosigma);
965
966  return std::exp(z);
967}
968
969/////////////////////////////////////////////////////////////////
970//
971//
972
973
974inline G4complex G4NuclNuclDiffuseElastic::GammaLess(G4double theta)
975{
976  G4double sinThetaR      = 2.*fHalfRutThetaTg/(1. + fHalfRutThetaTg*fHalfRutThetaTg);
977  G4double cosHalfThetaR2 = 1./(1. + fHalfRutThetaTg*fHalfRutThetaTg);
978
979  G4double u              = std::sqrt(0.5*fProfileLambda/sinThetaR);
980  G4double kappa          = u/std::sqrt(pi);
981  G4double dTheta         = theta - fRutherfordTheta;
982  u                      *= dTheta;
983  G4double u2             = u*u;
984  G4double u2m2p3         = u2*2./3.;
985
986  G4complex im            = G4complex(0.,1.);
987  G4complex order         = G4complex(u,u);
988  order                  /= std::sqrt(2.);
989  G4complex gamma         = pi*kappa*GetErfcInt(-order)*std::exp(im*(u*u+0.25*pi));
990  G4complex a0            = 0.5*(1. + 4.*(1.+im*u2)*cosHalfThetaR2/3.)/sinThetaR;
991  G4complex a1            = 0.5*(1. + 2.*(1.+im*u2m2p3)*cosHalfThetaR2)/sinThetaR;
992  G4complex out           = gamma*(1. - a1*dTheta) - a0;
993  return out;
994}
995
996/////////////////////////////////////////////////////////////////
997//
998//
999
1000inline G4complex G4NuclNuclDiffuseElastic::GammaMore(G4double theta)
1001{
1002  G4double sinThetaR      = 2.*fHalfRutThetaTg/(1. + fHalfRutThetaTg*fHalfRutThetaTg);
1003  G4double cosHalfThetaR2 = 1./(1. + fHalfRutThetaTg*fHalfRutThetaTg);
1004
1005  G4double u              = std::sqrt(0.5*fProfileLambda/sinThetaR);
1006  G4double kappa          = u/std::sqrt(pi);
1007  G4double dTheta         = theta - fRutherfordTheta;
1008  u                      *= dTheta;
1009  G4double u2             = u*u;
1010  G4double u2m2p3         = u2*2./3.;
1011
1012  G4complex im            = G4complex(0.,1.);
1013  G4complex order         = G4complex(u,u);
1014  order                  /= std::sqrt(2.);
1015  G4complex gamma         = pi*kappa*GetErfcInt(order)*std::exp(im*(u*u+0.25*pi));
1016  G4complex a0            = 0.5*(1. + 3.*(1.+im*u2)*cosHalfThetaR2/3.)/sinThetaR;
1017  G4complex a1            = 0.5*(1. + 2.*(1.+im*u2m2p3)*cosHalfThetaR2)/sinThetaR;
1018  G4complex out           = -gamma*(1. - a1*dTheta) - a0;
1019  return out;
1020}
1021
1022/////////////////////////////////////////////////////////////////
1023//
1024//
1025
1026inline  G4complex G4NuclNuclDiffuseElastic::AmplitudeNear(G4double theta)
1027{
1028  G4double kappa = std::sqrt(0.5*fProfileLambda/std::sin(theta)/pi);
1029  G4complex out = G4complex(kappa/fWaveVector,0.);
1030  out *= PhaseNear(theta);
1031
1032  if(theta <= fRutherfordTheta)
1033  {
1034    out *= GammaLess(theta) + ProfileNear(theta);
1035    out += CoulombAmplitude(theta);
1036  }
1037  else
1038  {
1039    out *= GammaMore(theta) + ProfileNear(theta);
1040  }
1041  return out;
1042}
1043
1044/////////////////////////////////////////////////////////////////
1045//
1046//
1047
1048inline  G4complex G4NuclNuclDiffuseElastic::AmplitudeFar(G4double theta)
1049{
1050  G4double kappa = std::sqrt(0.5*fProfileLambda/std::sin(theta)/pi);
1051  G4complex out = G4complex(kappa/fWaveVector,0.);
1052  out *= ProfileFar(theta)*PhaseFar(theta);
1053  return out;
1054}
1055
1056
1057/////////////////////////////////////////////////////////////////
1058//
1059//
1060
1061inline  G4complex G4NuclNuclDiffuseElastic::Amplitude(G4double theta)
1062{
1063 
1064  G4complex out = AmplitudeNear(theta) + AmplitudeFar(theta);
1065  return    out;
1066}
1067
1068/////////////////////////////////////////////////////////////////
1069//
1070//
1071
1072inline  G4double  G4NuclNuclDiffuseElastic::AmplitudeMod2(G4double theta)
1073{
1074  G4complex out = Amplitude(theta);
1075  G4double mod2 = out.real()*out.real() + out.imag()*out.imag();
1076  return   mod2;
1077}
1078
1079
1080///////////////////////////////////////////////////////////////////////////////
1081//
1082// Test for given particle and element table of momentum, angle probability.
1083// For the moment in lab system.
1084
1085inline void G4NuclNuclDiffuseElastic::InitParameters(const G4ParticleDefinition* theParticle, 
1086                                          G4double partMom, G4double Z, G4double A) 
1087{
1088  fAtomicNumber  = Z;     // atomic number
1089  fAtomicWeight  = A;     // number of nucleons
1090
1091  fNuclearRadius2 = CalculateNuclearRad(fAtomicWeight);
1092  G4double A1     = G4double( theParticle->GetBaryonNumber() );   
1093  fNuclearRadius1 = CalculateNuclearRad(A1);
1094  // fNuclearRadius = std::sqrt(fNuclearRadius1*fNuclearRadius1+fNuclearRadius2*fNuclearRadius2);
1095  fNuclearRadius = fNuclearRadius1 + fNuclearRadius2;
1096
1097  G4double a = 0.;
1098  G4double z = theParticle->GetPDGCharge();
1099  G4double m1 = theParticle->GetPDGMass();
1100
1101  fWaveVector = partMom/hbarc;
1102
1103  G4double lambda = fWaveVector*fNuclearRadius;
1104
1105  if( z )
1106  {
1107      a           = partMom/m1; // beta*gamma for m1
1108      fBeta       = a/std::sqrt(1+a*a);
1109      fZommerfeld = CalculateZommerfeld( fBeta, z, fAtomicNumber);
1110      fAm         = CalculateAm( partMom, fZommerfeld, fAtomicNumber);
1111  }
1112  fProfileLambda = lambda*std::sqrt(1.-2*fZommerfeld/lambda);
1113  G4cout<<"fProfileLambda = "<<fProfileLambda<<G4endl;
1114  fProfileDelta  = 0.1*fProfileLambda;
1115  fProfileAlpha  = 0.05*fProfileLambda;
1116
1117  CalculateCoulombPhaseZero();
1118  CalculateRutherfordAnglePar();
1119
1120  return;
1121}
1122
1123
1124
1125
1126
1127#endif
Note: See TracBrowser for help on using the repository browser.