Geant4_10
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: G4NuclNuclDiffuseElastic.hh 68039 2013-03-13 14:17:40Z gcosmo $
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  virtual ~G4NuclNuclDiffuseElastic();
68 
69  void Initialise();
70 
71  void InitialiseOnFly(G4double Z, G4double A);
72 
73  void BuildAngleTable();
74 
75  virtual G4double SampleInvariantT(const G4ParticleDefinition* p,
76  G4double plab,
77  G4int Z, G4int A);
78 
79  void SetPlabLowLimit(G4double value);
80 
81  void SetHEModelLowLimit(G4double value);
82 
83  void SetQModelLowLimit(G4double value);
84 
85  void SetLowestEnergyLimit(G4double value);
86 
87  void SetRecoilKinEnergyLimit(G4double value);
88 
89  G4double SampleT(const G4ParticleDefinition* aParticle,
90  G4double p, G4double A);
91 
92  G4double SampleTableT(const G4ParticleDefinition* aParticle,
93  G4double p, G4double Z, G4double A);
94 
95  G4double SampleThetaCMS(const G4ParticleDefinition* aParticle, G4double p, G4double A);
96 
97  G4double SampleTableThetaCMS(const G4ParticleDefinition* aParticle, G4double p,
98  G4double Z, G4double A);
99 
100  G4double GetScatteringAngle(G4int iMomentum, G4int iAngle, G4double position);
101 
102  G4double SampleThetaLab(const G4HadProjectile* aParticle,
103  G4double tmass, G4double A);
104 
105  G4double GetDiffuseElasticXsc( const G4ParticleDefinition* particle,
106  G4double theta,
107  G4double momentum,
108  G4double A );
109 
110  G4double GetInvElasticXsc( const G4ParticleDefinition* particle,
111  G4double theta,
112  G4double momentum,
113  G4double A, G4double Z );
114 
115  G4double GetDiffuseElasticSumXsc( const G4ParticleDefinition* particle,
116  G4double theta,
117  G4double momentum,
118  G4double A, G4double Z );
119 
120  G4double GetInvElasticSumXsc( const G4ParticleDefinition* particle,
121  G4double tMand,
122  G4double momentum,
123  G4double A, G4double Z );
124 
125  G4double IntegralElasticProb( const G4ParticleDefinition* particle,
126  G4double theta,
127  G4double momentum,
128  G4double A );
129 
130 
131  G4double GetCoulombElasticXsc( const G4ParticleDefinition* particle,
132  G4double theta,
133  G4double momentum,
134  G4double Z );
135 
136  G4double GetRutherfordXsc( G4double theta );
137 
138  G4double GetInvCoulombElasticXsc( const G4ParticleDefinition* particle,
139  G4double tMand,
140  G4double momentum,
141  G4double A, G4double Z );
142 
143  G4double GetCoulombTotalXsc( const G4ParticleDefinition* particle,
144  G4double momentum, G4double Z );
145 
146  G4double GetCoulombIntegralXsc( const G4ParticleDefinition* particle,
147  G4double momentum, G4double Z,
148  G4double theta1, G4double theta2 );
149 
150 
151  G4double CalculateParticleBeta( const G4ParticleDefinition* particle,
152  G4double momentum );
153 
154  G4double CalculateZommerfeld( G4double beta, G4double Z1, G4double Z2 );
155 
156  G4double CalculateAm( G4double momentum, G4double n, G4double Z);
157 
158  G4double CalculateNuclearRad( G4double A);
159 
160  G4double ThetaCMStoThetaLab(const G4DynamicParticle* aParticle,
161  G4double tmass, G4double thetaCMS);
162 
163  G4double ThetaLabToThetaCMS(const G4DynamicParticle* aParticle,
164  G4double tmass, G4double thetaLab);
165 
166  void TestAngleTable(const G4ParticleDefinition* theParticle, G4double partMom,
167  G4double Z, G4double A);
168 
169 
170 
171  G4double BesselJzero(G4double z);
172  G4double BesselJone(G4double z);
173  G4double DampFactor(G4double z);
174  G4double BesselOneByArg(G4double z);
175 
176  G4double GetDiffElasticProb(G4double theta);
177  G4double GetDiffElasticSumProb(G4double theta);
178  G4double GetDiffElasticSumProbA(G4double alpha);
179  G4double GetIntegrandFunction(G4double theta);
180 
181  G4double GetNuclearRadius(){return fNuclearRadius;};
182 
183 
184  // Technical math functions for strong Coulomb contribution
185 
186  G4complex GammaLogarithm(G4complex xx);
187  G4complex GammaLogB2n(G4complex xx);
188 
189  G4double GetErf(G4double x);
190 
191  G4double GetCosHaPit2(G4double t){return std::cos(CLHEP::halfpi*t*t);};
192  G4double GetSinHaPit2(G4double t){return std::sin(CLHEP::halfpi*t*t);};
193 
194  G4double GetCint(G4double x);
195  G4double GetSint(G4double x);
196 
197 
198  G4complex GetErfcComp(G4complex z, G4int nMax);
199  G4complex GetErfcSer(G4complex z, G4int nMax);
200  G4complex GetErfcInt(G4complex z); // , G4int nMax);
201 
202  G4complex GetErfComp(G4complex z, G4int nMax); // AandS algorithm != Ser, Int
203  G4complex GetErfSer(G4complex z, G4int nMax);
204 
205  G4double GetExpCos(G4double x);
206  G4double GetExpSin(G4double x);
207  G4complex GetErfInt(G4complex z); // , G4int nMax);
208 
209  G4double GetLegendrePol(G4int n, G4double x);
210 
211  G4complex TestErfcComp(G4complex z, G4int nMax);
212  G4complex TestErfcSer(G4complex z, G4int nMax);
213  G4complex TestErfcInt(G4complex z); // , G4int nMax);
214 
215  G4complex CoulombAmplitude(G4double theta);
216  G4double CoulombAmplitudeMod2(G4double theta);
217 
218  void CalculateCoulombPhaseZero();
219  G4double CalculateCoulombPhase(G4int n);
220  void CalculateRutherfordAnglePar();
221 
222  G4double ProfileNear(G4double theta);
223  G4double ProfileFar(G4double theta);
224  G4double Profile(G4double theta);
225 
226  G4complex PhaseNear(G4double theta);
227  G4complex PhaseFar(G4double theta);
228 
229  G4complex GammaLess(G4double theta);
230  G4complex GammaMore(G4double theta);
231 
232  G4complex AmplitudeNear(G4double theta);
233  G4complex AmplitudeFar(G4double theta);
234 
235  G4complex Amplitude(G4double theta);
236  G4double AmplitudeMod2(G4double theta);
237 
238  G4complex AmplitudeSim(G4double theta);
239  G4double AmplitudeSimMod2(G4double theta);
240 
241  G4double GetRatioSim(G4double theta);
242  G4double GetRatioGen(G4double theta);
243 
244  G4double GetFresnelDiffuseXsc(G4double theta);
245  G4double GetFresnelIntegrandXsc(G4double alpha);
246 
247 
248  G4complex AmplitudeGla(G4double theta);
249  G4double AmplitudeGlaMod2(G4double theta);
250 
251  G4complex AmplitudeGG(G4double theta);
252  G4double AmplitudeGGMod2(G4double theta);
253 
254  void InitParameters(const G4ParticleDefinition* theParticle,
255  G4double partMom, G4double Z, G4double A);
256 
257  void InitDynParameters(const G4ParticleDefinition* theParticle,
258  G4double partMom);
259 
260  void InitParametersGla(const G4DynamicParticle* aParticle,
261  G4double partMom, G4double Z, G4double A);
262 
263  G4double GetHadronNucleonXscNS( G4ParticleDefinition* pParticle,
264  G4double pTkin,
265  G4ParticleDefinition* tParticle);
266 
267  G4double CalcMandelstamS( const G4double mp , const G4double mt ,
268  const G4double Plab );
269 
270  G4double GetProfileLambda(){return fProfileLambda;};
271 
272  inline void SetProfileLambda(G4double pl) {fProfileLambda = pl;};
273  inline void SetProfileDelta(G4double pd) {fProfileDelta = pd;};
274  inline void SetProfileAlpha(G4double pa){fProfileAlpha = pa;};
275  inline void SetCofLambda(G4double pa){fCofLambda = pa;};
276 
277  inline void SetCofAlpha(G4double pa){fCofAlpha = pa;};
278  inline void SetCofAlphaMax(G4double pa){fCofAlphaMax = pa;};
279  inline void SetCofAlphaCoulomb(G4double pa){fCofAlphaCoulomb = pa;};
280 
281  inline void SetCofDelta(G4double pa){fCofDelta = pa;};
282  inline void SetCofPhase(G4double pa){fCofPhase = pa;};
283  inline void SetCofFar(G4double pa){fCofFar = pa;};
284  inline void SetEtaRatio(G4double pa){fEtaRatio = pa;};
285  inline void SetMaxL(G4int l){fMaxL = l;};
286  inline void SetNuclearRadiusCof(G4double r){fNuclearRadiusCof = r;};
287 
288  inline G4double GetCofAlphaMax(){return fCofAlphaMax;};
289  inline G4double GetCofAlphaCoulomb(){return fCofAlphaCoulomb;};
290 
291 private:
292 
293  G4ParticleDefinition* theProton;
294  G4ParticleDefinition* theNeutron;
295  G4ParticleDefinition* theDeuteron;
296  G4ParticleDefinition* theAlpha;
297 
298  const G4ParticleDefinition* thePionPlus;
299  const G4ParticleDefinition* thePionMinus;
300 
301  G4double lowEnergyRecoilLimit;
302  G4double lowEnergyLimitHE;
303  G4double lowEnergyLimitQ;
304  G4double lowestEnergyLimit;
305  G4double plabLowLimit;
306 
307  G4int fEnergyBin;
308  G4int fAngleBin;
309 
310  G4PhysicsLogVector* fEnergyVector;
311  G4PhysicsTable* fAngleTable;
312  std::vector<G4PhysicsTable*> fAngleBank;
313 
314  std::vector<G4double> fElementNumberVector;
315  std::vector<G4String> fElementNameVector;
316 
317  const G4ParticleDefinition* fParticle;
318 
319  G4double fWaveVector;
320  G4double fAtomicWeight;
321  G4double fAtomicNumber;
322 
323  G4double fNuclearRadius1;
324  G4double fNuclearRadius2;
325  G4double fNuclearRadius;
326  G4double fNuclearRadiusSquare;
327  G4double fNuclearRadiusCof;
328 
329  G4double fBeta;
330  G4double fZommerfeld;
331  G4double fRutherfordRatio;
332  G4double fAm;
333  G4bool fAddCoulomb;
334 
335  G4double fCoulombPhase0;
336  G4double fHalfRutThetaTg;
337  G4double fHalfRutThetaTg2;
338  G4double fRutherfordTheta;
339 
340  G4double fProfileLambda;
341  G4double fProfileDelta;
342  G4double fProfileAlpha;
343 
344  G4double fCofLambda;
345  G4double fCofAlpha;
346  G4double fCofDelta;
347  G4double fCofPhase;
348  G4double fCofFar;
349 
350  G4double fCofAlphaMax;
351  G4double fCofAlphaCoulomb;
352 
353  G4int fMaxL;
354  G4double fSumSigma;
355  G4double fEtaRatio;
356 
357  G4double fReZ;
358 
359 };
360 
361 
362 inline void G4NuclNuclDiffuseElastic::SetRecoilKinEnergyLimit(G4double value)
363 {
364  lowEnergyRecoilLimit = value;
365 }
366 
367 inline void G4NuclNuclDiffuseElastic::SetPlabLowLimit(G4double value)
368 {
369  plabLowLimit = value;
370 }
371 
372 inline void G4NuclNuclDiffuseElastic::SetHEModelLowLimit(G4double value)
373 {
374  lowEnergyLimitHE = value;
375 }
376 
377 inline void G4NuclNuclDiffuseElastic::SetQModelLowLimit(G4double value)
378 {
379  lowEnergyLimitQ = value;
380 }
381 
382 inline void G4NuclNuclDiffuseElastic::SetLowestEnergyLimit(G4double value)
383 {
384  lowestEnergyLimit = value;
385 }
386 
388 //
389 // damp factor in diffraction x/sh(x), x was already *pi
390 
391 inline G4double G4NuclNuclDiffuseElastic::DampFactor(G4double x)
392 {
393  G4double df;
394  G4double f2 = 2., f3 = 6., f4 = 24.; // first factorials
395 
396  // x *= pi;
397 
398  if( std::fabs(x) < 0.01 )
399  {
400  df = 1./(1. + x/f2 + x*x/f3 + x*x*x/f4);
401  }
402  else
403  {
404  df = x/std::sinh(x);
405  }
406  return df;
407 }
408 
409 
411 //
412 // return J1(x)/x with special case for small x
413 
414 inline G4double G4NuclNuclDiffuseElastic::BesselOneByArg(G4double x)
415 {
416  G4double x2, result;
417 
418  if( std::fabs(x) < 0.01 )
419  {
420  x *= 0.5;
421  x2 = x*x;
422  result = 2. - x2 + x2*x2/6.;
423  }
424  else
425  {
426  result = BesselJone(x)/x;
427  }
428  return result;
429 }
430 
432 //
433 // return particle beta
434 
435 inline G4double
436 G4NuclNuclDiffuseElastic::CalculateParticleBeta(const G4ParticleDefinition* particle,
437  G4double momentum)
438 {
439  G4double mass = particle->GetPDGMass();
440  G4double a = momentum/mass;
441  fBeta = a/std::sqrt(1+a*a);
442 
443  return fBeta;
444 }
445 
447 //
448 // return Zommerfeld parameter for Coulomb scattering
449 
450 inline G4double
451 G4NuclNuclDiffuseElastic::CalculateZommerfeld(G4double beta, G4double Z1, G4double Z2 )
452 {
453  fZommerfeld = CLHEP::fine_structure_const*Z1*Z2/beta;
454 
455  return fZommerfeld;
456 }
457 
459 //
460 // return Wentzel correction for Coulomb scattering
461 
462 inline G4double
463 G4NuclNuclDiffuseElastic::CalculateAm( G4double momentum, G4double n, G4double Z)
464 {
465  G4double k = momentum/CLHEP::hbarc;
466  G4double ch = 1.13 + 3.76*n*n;
467  G4double zn = 1.77*k*std::pow(Z,-1./3.)*CLHEP::Bohr_radius;
468  G4double zn2 = zn*zn;
469  fAm = ch/zn2;
470 
471  return fAm;
472 }
473 
475 //
476 // calculate nuclear radius for different atomic weights using different approximations
477 
478 inline G4double G4NuclNuclDiffuseElastic::CalculateNuclearRad( G4double A)
479 {
480  G4double r0 = 1.*CLHEP::fermi, radius;
481  // r0 *= 1.12;
482  // r0 *= 1.44;
483  r0 *= fNuclearRadiusCof;
484 
485  /*
486  if( A < 50. )
487  {
488  if( A > 10. ) r0 = 1.16*( 1 - std::pow(A, -2./3.) )*CLHEP::fermi; // 1.08*fermi;
489  else r0 = 1.1*CLHEP::fermi;
490 
491  radius = r0*std::pow(A, 1./3.);
492  }
493  else
494  {
495  r0 = 1.7*CLHEP::fermi; // 1.7*fermi;
496 
497  radius = r0*std::pow(A, 0.27); // 0.27);
498  }
499  */
500  radius = r0*std::pow(A, 1./3.);
501 
502  return radius;
503 }
504 
506 //
507 // return Coulomb scattering differential xsc with Wentzel correction. Test function
508 
509 inline G4double G4NuclNuclDiffuseElastic::GetCoulombElasticXsc( const G4ParticleDefinition* particle,
510  G4double theta,
511  G4double momentum,
512  G4double Z )
513 {
514  G4double sinHalfTheta = std::sin(0.5*theta);
515  G4double sinHalfTheta2 = sinHalfTheta*sinHalfTheta;
516  G4double beta = CalculateParticleBeta( particle, momentum);
517  G4double z = particle->GetPDGCharge();
518  G4double n = CalculateZommerfeld( beta, z, Z );
519  G4double am = CalculateAm( momentum, n, Z);
520  G4double k = momentum/CLHEP::hbarc;
521  G4double ch = 0.5*n/k;
522  G4double ch2 = ch*ch;
523  G4double xsc = ch2/((sinHalfTheta2+am)*(sinHalfTheta2+am));
524 
525  return xsc;
526 }
527 
529 //
530 // return Rutherford scattering differential xsc with Wentzel correction. For Sampling.
531 
532 inline G4double G4NuclNuclDiffuseElastic::GetRutherfordXsc(G4double theta)
533 {
534  G4double sinHalfTheta = std::sin(0.5*theta);
535  G4double sinHalfTheta2 = sinHalfTheta*sinHalfTheta;
536 
537  G4double ch2 = fRutherfordRatio*fRutherfordRatio;
538 
539  G4double xsc = ch2/(sinHalfTheta2+fAm)/(sinHalfTheta2+fAm);
540 
541  return xsc;
542 }
543 
545 //
546 // return Coulomb scattering total xsc with Wentzel correction
547 
548 inline G4double G4NuclNuclDiffuseElastic::GetCoulombTotalXsc( const G4ParticleDefinition* particle,
549  G4double momentum, G4double Z )
550 {
551  G4double beta = CalculateParticleBeta( particle, momentum);
552  G4cout<<"beta = "<<beta<<G4endl;
553  G4double z = particle->GetPDGCharge();
554  G4double n = CalculateZommerfeld( beta, z, Z );
555  G4cout<<"fZomerfeld = "<<n<<G4endl;
556  G4double am = CalculateAm( momentum, n, Z);
557  G4cout<<"cof Am = "<<am<<G4endl;
558  G4double k = momentum/CLHEP::hbarc;
559  G4cout<<"k = "<<k*CLHEP::fermi<<" 1/fermi"<<G4endl;
560  G4cout<<"k*Bohr_radius = "<<k*CLHEP::Bohr_radius<<G4endl;
561  G4double ch = n/k;
562  G4double ch2 = ch*ch;
563  G4double xsc = ch2*CLHEP::pi/(am +am*am);
564 
565  return xsc;
566 }
567 
569 //
570 // return Coulomb scattering xsc with Wentzel correction integrated between
571 // theta1 and < theta2
572 
573 inline G4double
574 G4NuclNuclDiffuseElastic::GetCoulombIntegralXsc(const G4ParticleDefinition* particle,
575  G4double momentum, G4double Z,
576  G4double theta1, G4double theta2 )
577 {
578  G4double c1 = std::cos(theta1);
579  //G4cout<<"c1 = "<<c1<<G4endl;
580  G4double c2 = std::cos(theta2);
581  // G4cout<<"c2 = "<<c2<<G4endl;
582  G4double beta = CalculateParticleBeta( particle, momentum);
583  // G4cout<<"beta = "<<beta<<G4endl;
584  G4double z = particle->GetPDGCharge();
585  G4double n = CalculateZommerfeld( beta, z, Z );
586  // G4cout<<"fZomerfeld = "<<n<<G4endl;
587  G4double am = CalculateAm( momentum, n, Z);
588  // G4cout<<"cof Am = "<<am<<G4endl;
589  G4double k = momentum/CLHEP::hbarc;
590  // G4cout<<"k = "<<k*CLHEP::fermi<<" 1/fermi"<<G4endl;
591  // G4cout<<"k*Bohr_radius = "<<k*CLHEP::Bohr_radius<<G4endl;
592  G4double ch = n/k;
593  G4double ch2 = ch*ch;
594  am *= 2.;
595  G4double xsc = ch2*CLHEP::twopi*(c1-c2)/((1 - c1 + am)*(1 - c2 + am));
596 
597  return xsc;
598 }
599 
601 //
602 // For the calculation of arg Gamma(z) one needs complex extension
603 // of ln(Gamma(z)) here is approximate algorithm
604 
605 inline G4complex G4NuclNuclDiffuseElastic::GammaLogB2n(G4complex z)
606 {
607  G4complex z1 = 12.*z;
608  G4complex z2 = z*z;
609  G4complex z3 = z2*z;
610  G4complex z5 = z2*z3;
611  G4complex z7 = z2*z5;
612 
613  z3 *= 360.;
614  z5 *= 1260.;
615  z7 *= 1680.;
616 
617  G4complex result = (z-0.5)*std::log(z)-z+0.5*std::log(CLHEP::twopi);
618  result += 1./z1 - 1./z3 +1./z5 -1./z7;
619  return result;
620 }
621 
623 //
624 //
625 
626 inline G4double G4NuclNuclDiffuseElastic::GetErf(G4double x)
627 {
628  G4double t, z, tmp, result;
629 
630  z = std::fabs(x);
631  t = 1.0/(1.0+0.5*z);
632 
633  tmp = t*exp(-z*z-1.26551223+t*(1.00002368+t*(0.37409196+t*(0.09678418+
634  t*(-0.18628806+t*(0.27886807+t*(-1.13520398+t*(1.48851587+
635  t*(-0.82215223+t*0.17087277)))))))));
636 
637  if( x >= 0.) result = 1. - tmp;
638  else result = 1. + tmp;
639 
640  return result;
641 }
642 
644 //
645 //
646 
647 inline G4complex G4NuclNuclDiffuseElastic::GetErfcComp(G4complex z, G4int nMax)
648 {
649  G4complex erfcz = 1. - GetErfComp( z, nMax);
650  return erfcz;
651 }
652 
654 //
655 //
656 
657 inline G4complex G4NuclNuclDiffuseElastic::GetErfcSer(G4complex z, G4int nMax)
658 {
659  G4complex erfcz = 1. - GetErfSer( z, nMax);
660  return erfcz;
661 }
662 
664 //
665 //
666 
667 inline G4complex G4NuclNuclDiffuseElastic::GetErfcInt(G4complex z) // , G4int nMax)
668 {
669  G4complex erfcz = 1. - GetErfInt( z); // , nMax);
670  return erfcz;
671 }
672 
673 
675 //
676 //
677 
678 inline G4complex G4NuclNuclDiffuseElastic::TestErfcComp(G4complex z, G4int nMax)
679 {
680  G4complex miz = G4complex( z.imag(), -z.real() );
681  G4complex erfcz = 1. - GetErfComp( miz, nMax);
682  G4complex w = std::exp(-z*z)*erfcz;
683  return w;
684 }
685 
687 //
688 //
689 
690 inline G4complex G4NuclNuclDiffuseElastic::TestErfcSer(G4complex z, G4int nMax)
691 {
692  G4complex miz = G4complex( z.imag(), -z.real() );
693  G4complex erfcz = 1. - GetErfSer( miz, nMax);
694  G4complex w = std::exp(-z*z)*erfcz;
695  return w;
696 }
697 
699 //
700 //
701 
702 inline G4complex G4NuclNuclDiffuseElastic::TestErfcInt(G4complex z) // , G4int nMax)
703 {
704  G4complex miz = G4complex( z.imag(), -z.real() );
705  G4complex erfcz = 1. - GetErfInt( miz); // , nMax);
706  G4complex w = std::exp(-z*z)*erfcz;
707  return w;
708 }
709 
711 //
712 //
713 
714 inline G4complex G4NuclNuclDiffuseElastic::GetErfSer(G4complex z, G4int nMax)
715 {
716  G4int n;
717  G4double a =1., b = 1., tmp;
718  G4complex sum = z, d = z;
719 
720  for( n = 1; n <= nMax; n++)
721  {
722  a *= 2.;
723  b *= 2.*n +1.;
724  d *= z*z;
725 
726  tmp = a/b;
727 
728  sum += tmp*d;
729  }
730  sum *= 2.*std::exp(-z*z)/std::sqrt(CLHEP::pi);
731 
732  return sum;
733 }
734 
736 
737 inline G4double G4NuclNuclDiffuseElastic::GetExpCos(G4double x)
738 {
739  G4double result;
740 
741  result = std::exp(x*x-fReZ*fReZ);
742  result *= std::cos(2.*x*fReZ);
743  return result;
744 }
745 
747 
748 inline G4double G4NuclNuclDiffuseElastic::GetExpSin(G4double x)
749 {
750  G4double result;
751 
752  result = std::exp(x*x-fReZ*fReZ);
753  result *= std::sin(2.*x*fReZ);
754  return result;
755 }
756 
757 
759 //
760 //
761 
762 inline G4double G4NuclNuclDiffuseElastic::GetCint(G4double x)
763 {
764  G4double out;
765 
766  G4Integrator<G4NuclNuclDiffuseElastic,G4double(G4NuclNuclDiffuseElastic::*)(G4double)> integral;
767 
768  out= integral.Legendre96(this,&G4NuclNuclDiffuseElastic::GetCosHaPit2, 0., x );
769 
770  return out;
771 }
772 
774 //
775 //
776 
777 inline G4double G4NuclNuclDiffuseElastic::GetSint(G4double x)
778 {
779  G4Integrator<G4NuclNuclDiffuseElastic,G4double(G4NuclNuclDiffuseElastic::*)(G4double)> integral;
780 
781  G4double out =
782  integral.Legendre96(this,&G4NuclNuclDiffuseElastic::GetSinHaPit2, 0., x );
783 
784  return out;
785 }
786 
787 
789 //
790 //
791 
792 inline G4complex G4NuclNuclDiffuseElastic::CoulombAmplitude(G4double theta)
793 {
794  G4complex ca;
795 
796  G4double sinHalfTheta = std::sin(0.5*theta);
797  G4double sinHalfTheta2 = sinHalfTheta*sinHalfTheta;
798  sinHalfTheta2 += fAm;
799 
800  G4double order = 2.*fCoulombPhase0 - fZommerfeld*std::log(sinHalfTheta2);
801  G4complex z = G4complex(0., order);
802  ca = std::exp(z);
803 
804  ca *= -fZommerfeld/(2.*fWaveVector*sinHalfTheta2);
805 
806  return ca;
807 }
808 
810 //
811 //
812 
813 inline G4double G4NuclNuclDiffuseElastic::CoulombAmplitudeMod2(G4double theta)
814 {
815  G4complex ca = CoulombAmplitude(theta);
816  G4double out = ca.real()*ca.real() + ca.imag()*ca.imag();
817 
818  return out;
819 }
820 
822 //
823 //
824 
825 
826 inline void G4NuclNuclDiffuseElastic::CalculateCoulombPhaseZero()
827 {
828  G4complex z = G4complex(1,fZommerfeld);
829  // G4complex gammalog = GammaLogarithm(z);
830  G4complex gammalog = GammaLogB2n(z);
831  fCoulombPhase0 = gammalog.imag();
832 }
833 
835 //
836 //
837 
838 
839 inline G4double G4NuclNuclDiffuseElastic::CalculateCoulombPhase(G4int n)
840 {
841  G4complex z = G4complex(1. + n, fZommerfeld);
842  // G4complex gammalog = GammaLogarithm(z);
843  G4complex gammalog = GammaLogB2n(z);
844  return gammalog.imag();
845 }
846 
847 
849 //
850 //
851 
852 
853 inline void G4NuclNuclDiffuseElastic::CalculateRutherfordAnglePar()
854 {
855  fHalfRutThetaTg = fZommerfeld/fProfileLambda; // (fWaveVector*fNuclearRadius);
856  fRutherfordTheta = 2.*std::atan(fHalfRutThetaTg);
857  fHalfRutThetaTg2 = fHalfRutThetaTg*fHalfRutThetaTg;
858  // G4cout<<"fRutherfordTheta = "<<fRutherfordTheta/degree<<" degree"<<G4endl;
859 
860 }
861 
863 //
864 //
865 
866 inline G4double G4NuclNuclDiffuseElastic::ProfileNear(G4double theta)
867 {
868  G4double dTheta = fRutherfordTheta - theta;
869  G4double result = 0., argument = 0.;
870 
871  if(std::abs(dTheta) < 0.001) result = fProfileAlpha*fProfileDelta;
872  else
873  {
874  argument = fProfileDelta*dTheta;
875  result = CLHEP::pi*argument*std::exp(fProfileAlpha*argument);
876  result /= std::sinh(CLHEP::pi*argument);
877  result -= 1.;
878  result /= dTheta;
879  }
880  return result;
881 }
882 
884 //
885 //
886 
887 inline G4double G4NuclNuclDiffuseElastic::ProfileFar(G4double theta)
888 {
889  G4double dTheta = fRutherfordTheta + theta;
890  G4double argument = fProfileDelta*dTheta;
891 
892  G4double result = CLHEP::pi*argument*std::exp(fProfileAlpha*argument);
893  result /= std::sinh(CLHEP::pi*argument);
894  result /= dTheta;
895 
896  return result;
897 }
898 
900 //
901 //
902 
903 inline G4double G4NuclNuclDiffuseElastic::Profile(G4double theta)
904 {
905  G4double dTheta = fRutherfordTheta - theta;
906  G4double result = 0., argument = 0.;
907 
908  if(std::abs(dTheta) < 0.001) result = 1.;
909  else
910  {
911  argument = fProfileDelta*dTheta;
912  result = CLHEP::pi*argument;
913  result /= std::sinh(result);
914  }
915  return result;
916 }
917 
919 //
920 //
921 
922 inline G4complex G4NuclNuclDiffuseElastic::PhaseNear(G4double theta)
923 {
924  G4double twosigma = 2.*fCoulombPhase0;
925  twosigma -= fZommerfeld*std::log(fHalfRutThetaTg2/(1.+fHalfRutThetaTg2));
926  twosigma += fRutherfordTheta*fZommerfeld/fHalfRutThetaTg - CLHEP::halfpi;
927  twosigma -= fProfileLambda*theta - 0.25*CLHEP::pi;
928 
929  twosigma *= fCofPhase;
930 
931  G4complex z = G4complex(0., twosigma);
932 
933  return std::exp(z);
934 }
935 
937 //
938 //
939 
940 inline G4complex G4NuclNuclDiffuseElastic::PhaseFar(G4double theta)
941 {
942  G4double twosigma = 2.*fCoulombPhase0;
943  twosigma -= fZommerfeld*std::log(fHalfRutThetaTg2/(1.+fHalfRutThetaTg2));
944  twosigma += fRutherfordTheta*fZommerfeld/fHalfRutThetaTg - CLHEP::halfpi;
945  twosigma += fProfileLambda*theta - 0.25*CLHEP::pi;
946 
947  twosigma *= fCofPhase;
948 
949  G4complex z = G4complex(0., twosigma);
950 
951  return std::exp(z);
952 }
953 
954 
956 //
957 //
958 
959 inline G4complex G4NuclNuclDiffuseElastic::AmplitudeFar(G4double theta)
960 {
961  G4double kappa = std::sqrt(0.5*fProfileLambda/std::sin(theta)/CLHEP::pi);
962  G4complex out = G4complex(kappa/fWaveVector,0.);
963  out *= ProfileFar(theta);
964  out *= PhaseFar(theta);
965  return out;
966 }
967 
968 
970 //
971 //
972 
973 inline G4complex G4NuclNuclDiffuseElastic::Amplitude(G4double theta)
974 {
975 
976  G4complex out = AmplitudeNear(theta) + fCofFar*AmplitudeFar(theta);
977  // G4complex out = AmplitudeNear(theta);
978  // G4complex out = AmplitudeFar(theta);
979  return out;
980 }
981 
983 //
984 //
985 
986 inline G4double G4NuclNuclDiffuseElastic::AmplitudeMod2(G4double theta)
987 {
988  G4complex out = Amplitude(theta);
989  G4double mod2 = out.real()*out.real() + out.imag()*out.imag();
990  return mod2;
991 }
992 
993 
995 //
996 //
997 
998 inline G4double G4NuclNuclDiffuseElastic::GetRatioSim(G4double theta)
999 {
1000  G4double sinThetaR = 2.*fHalfRutThetaTg/(1. + fHalfRutThetaTg2);
1001  G4double dTheta = 0.5*(theta - fRutherfordTheta);
1002  G4double sindTheta = std::sin(dTheta);
1003 
1004  G4double order = std::sqrt(fProfileLambda/sinThetaR/CLHEP::pi)*2.*sindTheta;
1005  // G4cout<<"order = "<<order<<G4endl;
1006  G4double cosFresnel = 0.5 - GetCint(order);
1007  G4double sinFresnel = 0.5 - GetSint(order);
1008 
1009  G4double out = 0.5*( cosFresnel*cosFresnel + sinFresnel*sinFresnel );
1010 
1011  return out;
1012 }
1013 
1015 //
1016 // The xsc for Fresnel smooth nucleus profile
1017 
1018 inline G4double G4NuclNuclDiffuseElastic::GetFresnelDiffuseXsc(G4double theta)
1019 {
1020  G4double ratio = GetRatioGen(theta);
1021  G4double ruthXsc = GetRutherfordXsc(theta);
1022  G4double xsc = ratio*ruthXsc;
1023  return xsc;
1024 }
1025 
1027 //
1028 // The xsc for Fresnel smooth nucleus profile for integration
1029 
1030 inline G4double G4NuclNuclDiffuseElastic::GetFresnelIntegrandXsc(G4double alpha)
1031 {
1032  G4double theta = std::sqrt(alpha);
1033  G4double xsc = GetFresnelDiffuseXsc(theta);
1034  return xsc;
1035 }
1036 
1038 //
1039 //
1040 
1041 inline G4double G4NuclNuclDiffuseElastic::AmplitudeSimMod2(G4double theta)
1042 {
1043  G4complex out = AmplitudeSim(theta);
1044  G4double mod2 = out.real()*out.real() + out.imag()*out.imag();
1045  return mod2;
1046 }
1047 
1048 
1050 //
1051 //
1052 
1053 inline G4double G4NuclNuclDiffuseElastic::AmplitudeGlaMod2(G4double theta)
1054 {
1055  G4complex out = AmplitudeGla(theta);
1056  G4double mod2 = out.real()*out.real() + out.imag()*out.imag();
1057  return mod2;
1058 }
1059 
1061 //
1062 //
1063 
1064 inline G4double G4NuclNuclDiffuseElastic::AmplitudeGGMod2(G4double theta)
1065 {
1066  G4complex out = AmplitudeGG(theta);
1067  G4double mod2 = out.real()*out.real() + out.imag()*out.imag();
1068  return mod2;
1069 }
1070 
1072 //
1073 //
1074 
1075 inline G4double G4NuclNuclDiffuseElastic::CalcMandelstamS( const G4double mp ,
1076  const G4double mt ,
1077  const G4double Plab )
1078 {
1079  G4double Elab = std::sqrt ( mp * mp + Plab * Plab );
1080  G4double sMand = mp*mp + mt*mt + 2*Elab*mt ;
1081 
1082  return sMand;
1083 }
1084 
1085 
1086 
1087 #endif
Z2
Double_t Z2
Definition: plot.C:266
G4NuclNuclDiffuseElastic::GetCoulombIntegralXsc
G4double GetCoulombIntegralXsc(const G4ParticleDefinition *particle, G4double momentum, G4double Z, G4double theta1, G4double theta2)
Definition: G4NuclNuclDiffuseElastic.hh:574
tmp
Float_t tmp
Definition: plot.C:37
G4NuclNuclDiffuseElastic::CalculateNuclearRad
G4double CalculateNuclearRad(G4double A)
Definition: G4NuclNuclDiffuseElastic.hh:478
G4Integrator::Legendre96
G4double Legendre96(T &typeT, F f, G4double a, G4double b)
G4NuclNuclDiffuseElastic::GetErfcComp
G4complex GetErfcComp(G4complex z, G4int nMax)
Definition: G4NuclNuclDiffuseElastic.hh:647
x2
Double_t x2[nxs]
Definition: Style.C:19
test.a
tuple a
Definition: test.py:11
d
Float_t d
Definition: plot.C:237
G4NuclNuclDiffuseElastic::AmplitudeSimMod2
G4double AmplitudeSimMod2(G4double theta)
Definition: G4NuclNuclDiffuseElastic.hh:1041
G4NuclNuclDiffuseElastic::CalculateAm
G4double CalculateAm(G4double momentum, G4double n, G4double Z)
Definition: G4NuclNuclDiffuseElastic.hh:463
G4HadProjectile.hh
G4HadronElastic.hh
p
const char * p
Definition: xmltok.h:285
xx
Double_t xx
Definition: macro.C:10
result
G4double G4NeutronHPJENDLHEData::G4double result
Definition: G4NeutronHPJENDLHEData.cc:239
G4NuclNuclDiffuseElastic::PhaseNear
G4complex PhaseNear(G4double theta)
Definition: G4NuclNuclDiffuseElastic.hh:922
c1
TCanvas * c1
Definition: plotHisto.C:7
f4
Float_t f4
Definition: comparison_ascii.C:56
G4NuclNuclDiffuseElastic::CoulombAmplitude
G4complex CoulombAmplitude(G4double theta)
Definition: G4NuclNuclDiffuseElastic.hh:792
test.x
tuple x
Definition: test.py:50
G4NuclNuclDiffuseElastic::GetRatioSim
G4double GetRatioSim(G4double theta)
Definition: G4NuclNuclDiffuseElastic.hh:998
G4NuclNuclDiffuseElastic::GetCoulombElasticXsc
G4double GetCoulombElasticXsc(const G4ParticleDefinition *particle, G4double theta, G4double momentum, G4double Z)
Definition: G4NuclNuclDiffuseElastic.hh:509
G4NuclNuclDiffuseElastic::CalculateParticleBeta
G4double CalculateParticleBeta(const G4ParticleDefinition *particle, G4double momentum)
Definition: G4NuclNuclDiffuseElastic.hh:436
G4int
int G4int
Definition: G4Types.hh:78
G4NuclNuclDiffuseElastic::GetCoulombTotalXsc
G4double GetCoulombTotalXsc(const G4ParticleDefinition *particle, G4double momentum, G4double Z)
Definition: G4NuclNuclDiffuseElastic.hh:548
G4NuclNuclDiffuseElastic::Profile
G4double Profile(G4double theta)
Definition: G4NuclNuclDiffuseElastic.hh:903
G4NuclNuclDiffuseElastic::GammaLogB2n
G4complex GammaLogB2n(G4complex xx)
Definition: G4NuclNuclDiffuseElastic.hh:605
G4NuclNuclDiffuseElastic::SetPlabLowLimit
void SetPlabLowLimit(G4double value)
Definition: G4NuclNuclDiffuseElastic.hh:367
test.b
tuple b
Definition: test.py:12
n
Char_t n[5]
Definition: comparison_ascii.C:55
readPY.pl
tuple pl
Definition: readPY.py:5
G4complex
std::complex< G4double > G4complex
Definition: G4Types.hh:81
G4cout
G4GLOB_DLL std::ostream G4cout
f2
Float_t f2
Definition: comparison_ascii.C:56
Z
Float_t Z
Definition: plot.C:39
G4NuclNuclDiffuseElastic::SetRecoilKinEnergyLimit
void SetRecoilKinEnergyLimit(G4double value)
Definition: G4NuclNuclDiffuseElastic.hh:362
G4NuclNuclDiffuseElastic::SetLowestEnergyLimit
void SetLowestEnergyLimit(G4double value)
Definition: G4NuclNuclDiffuseElastic.hh:382
G4bool
bool G4bool
Definition: G4Types.hh:79
G4Nucleus.hh
G4NuclNuclDiffuseElastic::CoulombAmplitudeMod2
G4double CoulombAmplitudeMod2(G4double theta)
Definition: G4NuclNuclDiffuseElastic.hh:813
G4NuclNuclDiffuseElastic::GetRutherfordXsc
G4double GetRutherfordXsc(G4double theta)
Definition: G4NuclNuclDiffuseElastic.hh:532
PhysicalConstants.h
G4NuclNuclDiffuseElastic::TestErfcInt
G4complex TestErfcInt(G4complex z)
Definition: G4NuclNuclDiffuseElastic.hh:702
G4NuclNuclDiffuseElastic::SetHEModelLowLimit
void SetHEModelLowLimit(G4double value)
Definition: G4NuclNuclDiffuseElastic.hh:372
G4NuclNuclDiffuseElastic::GetFresnelDiffuseXsc
G4double GetFresnelDiffuseXsc(G4double theta)
Definition: G4NuclNuclDiffuseElastic.hh:1018
f3
Float_t f3
Definition: comparison_ascii.C:56
r
jump r
Definition: plot.C:36
G4NuclNuclDiffuseElastic::CalculateZommerfeld
G4double CalculateZommerfeld(G4double beta, G4double Z1, G4double Z2)
Definition: G4NuclNuclDiffuseElastic.hh:451
G4ParticleDefinition::GetPDGMass
G4double GetPDGMass() const
Definition: G4ParticleDefinition.hh:161
G4NuclNuclDiffuseElastic::SetQModelLowLimit
void SetQModelLowLimit(G4double value)
Definition: G4NuclNuclDiffuseElastic.hh:377
G4NuclNuclDiffuseElastic::ProfileNear
G4double ProfileNear(G4double theta)
Definition: G4NuclNuclDiffuseElastic.hh:866
G4NuclNuclDiffuseElastic::GetFresnelIntegrandXsc
G4double GetFresnelIntegrandXsc(G4double alpha)
Definition: G4NuclNuclDiffuseElastic.hh:1030
G4NuclNuclDiffuseElastic::GetErfcSer
G4complex GetErfcSer(G4complex z, G4int nMax)
Definition: G4NuclNuclDiffuseElastic.hh:657
G4NuclNuclDiffuseElastic::AmplitudeMod2
G4double AmplitudeMod2(G4double theta)
Definition: G4NuclNuclDiffuseElastic.hh:986
G4NuclNuclDiffuseElastic::GetErfcInt
G4complex GetErfcInt(G4complex z)
Definition: G4NuclNuclDiffuseElastic.hh:667
G4Integrator.hh
test.z
tuple z
Definition: test.py:28
G4NuclNuclDiffuseElastic::Amplitude
G4complex Amplitude(G4double theta)
Definition: G4NuclNuclDiffuseElastic.hh:973
value
const XML_Char int const XML_Char * value
Definition: expat.h:331
G4endl
#define G4endl
Definition: G4ios.hh:61
G4NuclNuclDiffuseElastic::AmplitudeFar
G4complex AmplitudeFar(G4double theta)
Definition: G4NuclNuclDiffuseElastic.hh:959
G4NuclNuclDiffuseElastic::CalcMandelstamS
G4double CalcMandelstamS(const G4double mp, const G4double mt, const G4double Plab)
Definition: G4NuclNuclDiffuseElastic.hh:1075
G4NuclNuclDiffuseElastic::AmplitudeGGMod2
G4double AmplitudeGGMod2(G4double theta)
Definition: G4NuclNuclDiffuseElastic.hh:1064
G4NuclNuclDiffuseElastic::TestErfcComp
G4complex TestErfcComp(G4complex z, G4int nMax)
Definition: G4NuclNuclDiffuseElastic.hh:678
G4double
double G4double
Definition: G4Types.hh:76
G4ParticleDefinition::GetPDGCharge
G4double GetPDGCharge() const
Definition: G4ParticleDefinition.hh:163
G4NuclNuclDiffuseElastic::AmplitudeGlaMod2
G4double AmplitudeGlaMod2(G4double theta)
Definition: G4NuclNuclDiffuseElastic.hh:1053
G4NuclNuclDiffuseElastic::PhaseFar
G4complex PhaseFar(G4double theta)
Definition: G4NuclNuclDiffuseElastic.hh:940
G4NuclNuclDiffuseElastic::ProfileFar
G4double ProfileFar(G4double theta)
Definition: G4NuclNuclDiffuseElastic.hh:887
G4NuclNuclDiffuseElastic::GetErfSer
G4complex GetErfSer(G4complex z, G4int nMax)
Definition: G4NuclNuclDiffuseElastic.hh:714
Z1
Double_t Z1
Definition: plot.C:266
G4NuclNuclDiffuseElastic::TestErfcSer
G4complex TestErfcSer(G4complex z, G4int nMax)
Definition: G4NuclNuclDiffuseElastic.hh:690
Generated on Sat Dec 14 2013 14:34:46 for Geant4_10 by   doxygen 1.8.5