00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039
00040
00041
00042 #ifndef G4NuclNuclDiffuseElastic_h
00043 #define G4NuclNuclDiffuseElastic_h 1
00044
00045 #include <complex>
00046 #include <CLHEP/Units/PhysicalConstants.h>
00047 #include "globals.hh"
00048 #include "G4Integrator.hh"
00049 #include "G4HadronElastic.hh"
00050 #include "G4HadProjectile.hh"
00051 #include "G4Nucleus.hh"
00052
00053 using namespace std;
00054
00055 class G4ParticleDefinition;
00056 class G4PhysicsTable;
00057 class G4PhysicsLogVector;
00058
00059 class G4NuclNuclDiffuseElastic : public G4HadronElastic
00060 {
00061 public:
00062
00063 G4NuclNuclDiffuseElastic();
00064
00065
00066
00067
00068
00069
00070
00071 virtual ~G4NuclNuclDiffuseElastic();
00072
00073 void Initialise();
00074
00075 void InitialiseOnFly(G4double Z, G4double A);
00076
00077 void BuildAngleTable();
00078
00079
00080
00081
00082 virtual G4double SampleInvariantT(const G4ParticleDefinition* p,
00083 G4double plab,
00084 G4int Z, G4int A);
00085
00086 void SetPlabLowLimit(G4double value);
00087
00088 void SetHEModelLowLimit(G4double value);
00089
00090 void SetQModelLowLimit(G4double value);
00091
00092 void SetLowestEnergyLimit(G4double value);
00093
00094 void SetRecoilKinEnergyLimit(G4double value);
00095
00096 G4double SampleT(const G4ParticleDefinition* aParticle,
00097 G4double p, G4double A);
00098
00099 G4double SampleTableT(const G4ParticleDefinition* aParticle,
00100 G4double p, G4double Z, G4double A);
00101
00102 G4double SampleThetaCMS(const G4ParticleDefinition* aParticle, G4double p, G4double A);
00103
00104 G4double SampleTableThetaCMS(const G4ParticleDefinition* aParticle, G4double p,
00105 G4double Z, G4double A);
00106
00107 G4double GetScatteringAngle(G4int iMomentum, G4int iAngle, G4double position);
00108
00109 G4double SampleThetaLab(const G4HadProjectile* aParticle,
00110 G4double tmass, G4double A);
00111
00112 G4double GetDiffuseElasticXsc( const G4ParticleDefinition* particle,
00113 G4double theta,
00114 G4double momentum,
00115 G4double A );
00116
00117 G4double GetInvElasticXsc( const G4ParticleDefinition* particle,
00118 G4double theta,
00119 G4double momentum,
00120 G4double A, G4double Z );
00121
00122 G4double GetDiffuseElasticSumXsc( const G4ParticleDefinition* particle,
00123 G4double theta,
00124 G4double momentum,
00125 G4double A, G4double Z );
00126
00127 G4double GetInvElasticSumXsc( const G4ParticleDefinition* particle,
00128 G4double tMand,
00129 G4double momentum,
00130 G4double A, G4double Z );
00131
00132 G4double IntegralElasticProb( const G4ParticleDefinition* particle,
00133 G4double theta,
00134 G4double momentum,
00135 G4double A );
00136
00137
00138 G4double GetCoulombElasticXsc( const G4ParticleDefinition* particle,
00139 G4double theta,
00140 G4double momentum,
00141 G4double Z );
00142
00143 G4double GetRutherfordXsc( G4double theta );
00144
00145 G4double GetInvCoulombElasticXsc( const G4ParticleDefinition* particle,
00146 G4double tMand,
00147 G4double momentum,
00148 G4double A, G4double Z );
00149
00150 G4double GetCoulombTotalXsc( const G4ParticleDefinition* particle,
00151 G4double momentum, G4double Z );
00152
00153 G4double GetCoulombIntegralXsc( const G4ParticleDefinition* particle,
00154 G4double momentum, G4double Z,
00155 G4double theta1, G4double theta2 );
00156
00157
00158 G4double CalculateParticleBeta( const G4ParticleDefinition* particle,
00159 G4double momentum );
00160
00161 G4double CalculateZommerfeld( G4double beta, G4double Z1, G4double Z2 );
00162
00163 G4double CalculateAm( G4double momentum, G4double n, G4double Z);
00164
00165 G4double CalculateNuclearRad( G4double A);
00166
00167 G4double ThetaCMStoThetaLab(const G4DynamicParticle* aParticle,
00168 G4double tmass, G4double thetaCMS);
00169
00170 G4double ThetaLabToThetaCMS(const G4DynamicParticle* aParticle,
00171 G4double tmass, G4double thetaLab);
00172
00173 void TestAngleTable(const G4ParticleDefinition* theParticle, G4double partMom,
00174 G4double Z, G4double A);
00175
00176
00177
00178 G4double BesselJzero(G4double z);
00179 G4double BesselJone(G4double z);
00180 G4double DampFactor(G4double z);
00181 G4double BesselOneByArg(G4double z);
00182
00183 G4double GetDiffElasticProb(G4double theta);
00184 G4double GetDiffElasticSumProb(G4double theta);
00185 G4double GetDiffElasticSumProbA(G4double alpha);
00186 G4double GetIntegrandFunction(G4double theta);
00187
00188 G4double GetNuclearRadius(){return fNuclearRadius;};
00189
00190
00191
00192
00193 G4complex GammaLogarithm(G4complex xx);
00194 G4complex GammaLogB2n(G4complex xx);
00195
00196 G4double GetErf(G4double x);
00197
00198 G4double GetCosHaPit2(G4double t){return std::cos(CLHEP::halfpi*t*t);};
00199 G4double GetSinHaPit2(G4double t){return std::sin(CLHEP::halfpi*t*t);};
00200
00201 G4double GetCint(G4double x);
00202 G4double GetSint(G4double x);
00203
00204
00205 G4complex GetErfcComp(G4complex z, G4int nMax);
00206 G4complex GetErfcSer(G4complex z, G4int nMax);
00207 G4complex GetErfcInt(G4complex z);
00208
00209 G4complex GetErfComp(G4complex z, G4int nMax);
00210 G4complex GetErfSer(G4complex z, G4int nMax);
00211
00212 G4double GetExpCos(G4double x);
00213 G4double GetExpSin(G4double x);
00214 G4complex GetErfInt(G4complex z);
00215
00216
00217
00218
00219 G4double GetLegendrePol(G4int n, G4double x);
00220
00221 G4complex TestErfcComp(G4complex z, G4int nMax);
00222 G4complex TestErfcSer(G4complex z, G4int nMax);
00223 G4complex TestErfcInt(G4complex z);
00224
00225 G4complex CoulombAmplitude(G4double theta);
00226 G4double CoulombAmplitudeMod2(G4double theta);
00227
00228
00229 void CalculateCoulombPhaseZero();
00230 G4double CalculateCoulombPhase(G4int n);
00231 void CalculateRutherfordAnglePar();
00232
00233 G4double ProfileNear(G4double theta);
00234 G4double ProfileFar(G4double theta);
00235 G4double Profile(G4double theta);
00236
00237 G4complex PhaseNear(G4double theta);
00238 G4complex PhaseFar(G4double theta);
00239
00240 G4complex GammaLess(G4double theta);
00241 G4complex GammaMore(G4double theta);
00242
00243 G4complex AmplitudeNear(G4double theta);
00244 G4complex AmplitudeFar(G4double theta);
00245
00246 G4complex Amplitude(G4double theta);
00247 G4double AmplitudeMod2(G4double theta);
00248
00249 G4complex AmplitudeSim(G4double theta);
00250 G4double AmplitudeSimMod2(G4double theta);
00251
00252 G4double GetRatioSim(G4double theta);
00253 G4double GetRatioGen(G4double theta);
00254
00255 G4double GetFresnelDiffuseXsc(G4double theta);
00256 G4double GetFresnelIntegrandXsc(G4double alpha);
00257
00258
00259 G4complex AmplitudeGla(G4double theta);
00260 G4double AmplitudeGlaMod2(G4double theta);
00261
00262 G4complex AmplitudeGG(G4double theta);
00263 G4double AmplitudeGGMod2(G4double theta);
00264
00265 void InitParameters(const G4ParticleDefinition* theParticle,
00266 G4double partMom, G4double Z, G4double A);
00267
00268 void InitDynParameters(const G4ParticleDefinition* theParticle,
00269 G4double partMom);
00270
00271 void InitParametersGla(const G4DynamicParticle* aParticle,
00272 G4double partMom, G4double Z, G4double A);
00273
00274 G4double GetHadronNucleonXscNS( G4ParticleDefinition* pParticle,
00275 G4double pTkin,
00276 G4ParticleDefinition* tParticle);
00277
00278 G4double CalcMandelstamS( const G4double mp ,
00279 const G4double mt ,
00280 const G4double Plab );
00281
00282 G4double GetProfileLambda(){return fProfileLambda;};
00283
00284 void SetProfileLambda(G4double pl) {fProfileLambda = pl;};
00285 void SetProfileDelta(G4double pd) {fProfileDelta = pd;};
00286 void SetProfileAlpha(G4double pa){fProfileAlpha = pa;};
00287 void SetCofLambda(G4double pa){fCofLambda = pa;};
00288
00289 void SetCofAlpha(G4double pa){fCofAlpha = pa;};
00290 void SetCofAlphaMax(G4double pa){fCofAlphaMax = pa;};
00291 void SetCofAlphaCoulomb(G4double pa){fCofAlphaCoulomb = pa;};
00292
00293 void SetCofDelta(G4double pa){fCofDelta = pa;};
00294 void SetCofPhase(G4double pa){fCofPhase = pa;};
00295 void SetCofFar(G4double pa){fCofFar = pa;};
00296 void SetEtaRatio(G4double pa){fEtaRatio = pa;};
00297 void SetMaxL(G4int l){fMaxL = l;};
00298 void SetNuclearRadiusCof(G4double r){fNuclearRadiusCof = r;};
00299
00300 G4double GetCofAlphaMax(){return fCofAlphaMax;};
00301 G4double GetCofAlphaCoulomb(){return fCofAlphaCoulomb;};
00302
00303 private:
00304
00305
00306 G4ParticleDefinition* theProton;
00307 G4ParticleDefinition* theNeutron;
00308 G4ParticleDefinition* theDeuteron;
00309 G4ParticleDefinition* theAlpha;
00310
00311 const G4ParticleDefinition* thePionPlus;
00312 const G4ParticleDefinition* thePionMinus;
00313
00314 G4double lowEnergyRecoilLimit;
00315 G4double lowEnergyLimitHE;
00316 G4double lowEnergyLimitQ;
00317 G4double lowestEnergyLimit;
00318 G4double plabLowLimit;
00319
00320 G4int fEnergyBin;
00321 G4int fAngleBin;
00322
00323 G4PhysicsLogVector* fEnergyVector;
00324 G4PhysicsTable* fAngleTable;
00325 std::vector<G4PhysicsTable*> fAngleBank;
00326
00327 std::vector<G4double> fElementNumberVector;
00328 std::vector<G4String> fElementNameVector;
00329
00330 const G4ParticleDefinition* fParticle;
00331
00332 G4double fWaveVector;
00333 G4double fAtomicWeight;
00334 G4double fAtomicNumber;
00335
00336 G4double fNuclearRadius1;
00337 G4double fNuclearRadius2;
00338 G4double fNuclearRadius;
00339 G4double fNuclearRadiusSquare;
00340 G4double fNuclearRadiusCof;
00341
00342 G4double fBeta;
00343 G4double fZommerfeld;
00344 G4double fRutherfordRatio;
00345 G4double fAm;
00346 G4bool fAddCoulomb;
00347
00348 G4double fCoulombPhase0;
00349 G4double fHalfRutThetaTg;
00350 G4double fHalfRutThetaTg2;
00351 G4double fRutherfordTheta;
00352
00353 G4double fProfileLambda;
00354 G4double fProfileDelta;
00355 G4double fProfileAlpha;
00356
00357 G4double fCofLambda;
00358 G4double fCofAlpha;
00359 G4double fCofDelta;
00360 G4double fCofPhase;
00361 G4double fCofFar;
00362
00363 G4double fCofAlphaMax;
00364 G4double fCofAlphaCoulomb;
00365
00366 G4int fMaxL;
00367 G4double fSumSigma;
00368 G4double fEtaRatio;
00369
00370 G4double fReZ;
00371
00372 };
00373
00374
00375 inline void G4NuclNuclDiffuseElastic::SetRecoilKinEnergyLimit(G4double value)
00376 {
00377 lowEnergyRecoilLimit = value;
00378 }
00379
00380 inline void G4NuclNuclDiffuseElastic::SetPlabLowLimit(G4double value)
00381 {
00382 plabLowLimit = value;
00383 }
00384
00385 inline void G4NuclNuclDiffuseElastic::SetHEModelLowLimit(G4double value)
00386 {
00387 lowEnergyLimitHE = value;
00388 }
00389
00390 inline void G4NuclNuclDiffuseElastic::SetQModelLowLimit(G4double value)
00391 {
00392 lowEnergyLimitQ = value;
00393 }
00394
00395 inline void G4NuclNuclDiffuseElastic::SetLowestEnergyLimit(G4double value)
00396 {
00397 lowestEnergyLimit = value;
00398 }
00399
00400
00402
00403
00404
00405
00406 inline G4double G4NuclNuclDiffuseElastic::BesselJzero(G4double value)
00407 {
00408 G4double modvalue, value2, fact1, fact2, arg, shift, bessel;
00409
00410 modvalue = fabs(value);
00411
00412 if ( value < 8.0 && value > -8.0 )
00413 {
00414 value2 = value*value;
00415
00416 fact1 = 57568490574.0 + value2*(-13362590354.0
00417 + value2*( 651619640.7
00418 + value2*(-11214424.18
00419 + value2*( 77392.33017
00420 + value2*(-184.9052456 ) ) ) ) );
00421
00422 fact2 = 57568490411.0 + value2*( 1029532985.0
00423 + value2*( 9494680.718
00424 + value2*(59272.64853
00425 + value2*(267.8532712
00426 + value2*1.0 ) ) ) );
00427
00428 bessel = fact1/fact2;
00429 }
00430 else
00431 {
00432 arg = 8.0/modvalue;
00433
00434 value2 = arg*arg;
00435
00436 shift = modvalue-0.785398164;
00437
00438 fact1 = 1.0 + value2*(-0.1098628627e-2
00439 + value2*(0.2734510407e-4
00440 + value2*(-0.2073370639e-5
00441 + value2*0.2093887211e-6 ) ) );
00442
00443 fact2 = -0.1562499995e-1 + value2*(0.1430488765e-3
00444 + value2*(-0.6911147651e-5
00445 + value2*(0.7621095161e-6
00446 - value2*0.934945152e-7 ) ) );
00447
00448 bessel = sqrt(0.636619772/modvalue)*(cos(shift)*fact1 - arg*sin(shift)*fact2 );
00449 }
00450 return bessel;
00451 }
00452
00454
00455
00456
00457
00458 inline G4double G4NuclNuclDiffuseElastic::BesselJone(G4double value)
00459 {
00460 G4double modvalue, value2, fact1, fact2, arg, shift, bessel;
00461
00462 modvalue = fabs(value);
00463
00464 if ( modvalue < 8.0 )
00465 {
00466 value2 = value*value;
00467
00468 fact1 = value*(72362614232.0 + value2*(-7895059235.0
00469 + value2*( 242396853.1
00470 + value2*(-2972611.439
00471 + value2*( 15704.48260
00472 + value2*(-30.16036606 ) ) ) ) ) );
00473
00474 fact2 = 144725228442.0 + value2*(2300535178.0
00475 + value2*(18583304.74
00476 + value2*(99447.43394
00477 + value2*(376.9991397
00478 + value2*1.0 ) ) ) );
00479 bessel = fact1/fact2;
00480 }
00481 else
00482 {
00483 arg = 8.0/modvalue;
00484
00485 value2 = arg*arg;
00486
00487 shift = modvalue - 2.356194491;
00488
00489 fact1 = 1.0 + value2*( 0.183105e-2
00490 + value2*(-0.3516396496e-4
00491 + value2*(0.2457520174e-5
00492 + value2*(-0.240337019e-6 ) ) ) );
00493
00494 fact2 = 0.04687499995 + value2*(-0.2002690873e-3
00495 + value2*( 0.8449199096e-5
00496 + value2*(-0.88228987e-6
00497 + value2*0.105787412e-6 ) ) );
00498
00499 bessel = sqrt( 0.636619772/modvalue)*(cos(shift)*fact1 - arg*sin(shift)*fact2);
00500
00501 if (value < 0.0) bessel = -bessel;
00502 }
00503 return bessel;
00504 }
00505
00507
00508
00509
00510 inline G4double G4NuclNuclDiffuseElastic::DampFactor(G4double x)
00511 {
00512 G4double df;
00513 G4double f2 = 2., f3 = 6., f4 = 24.;
00514
00515
00516
00517 if( std::fabs(x) < 0.01 )
00518 {
00519 df = 1./(1. + x/f2 + x*x/f3 + x*x*x/f4);
00520 }
00521 else
00522 {
00523 df = x/std::sinh(x);
00524 }
00525 return df;
00526 }
00527
00528
00530
00531
00532
00533 inline G4double G4NuclNuclDiffuseElastic::BesselOneByArg(G4double x)
00534 {
00535 G4double x2, result;
00536
00537 if( std::fabs(x) < 0.01 )
00538 {
00539 x *= 0.5;
00540 x2 = x*x;
00541 result = 2. - x2 + x2*x2/6.;
00542 }
00543 else
00544 {
00545 result = BesselJone(x)/x;
00546 }
00547 return result;
00548 }
00549
00551
00552
00553
00554 inline G4double G4NuclNuclDiffuseElastic::CalculateParticleBeta( const G4ParticleDefinition* particle,
00555 G4double momentum )
00556 {
00557 G4double mass = particle->GetPDGMass();
00558 G4double a = momentum/mass;
00559 fBeta = a/std::sqrt(1+a*a);
00560
00561 return fBeta;
00562 }
00563
00565
00566
00567
00568 inline G4double G4NuclNuclDiffuseElastic::CalculateZommerfeld( G4double beta, G4double Z1, G4double Z2 )
00569 {
00570 fZommerfeld = CLHEP::fine_structure_const*Z1*Z2/beta;
00571
00572 return fZommerfeld;
00573 }
00574
00576
00577
00578
00579 inline G4double G4NuclNuclDiffuseElastic::CalculateAm( G4double momentum, G4double n, G4double Z)
00580 {
00581 G4double k = momentum/CLHEP::hbarc;
00582 G4double ch = 1.13 + 3.76*n*n;
00583 G4double zn = 1.77*k*std::pow(Z,-1./3.)*CLHEP::Bohr_radius;
00584 G4double zn2 = zn*zn;
00585 fAm = ch/zn2;
00586
00587 return fAm;
00588 }
00589
00591
00592
00593
00594 inline G4double G4NuclNuclDiffuseElastic::CalculateNuclearRad( G4double A)
00595 {
00596 G4double r0 = 1.*CLHEP::fermi, radius;
00597
00598
00599 r0 *= fNuclearRadiusCof;
00600
00601
00602
00603
00604
00605
00606
00607
00608
00609
00610
00611
00612
00613
00614
00615
00616 radius = r0*std::pow(A, 1./3.);
00617
00618 return radius;
00619 }
00620
00622
00623
00624
00625 inline G4double G4NuclNuclDiffuseElastic::GetCoulombElasticXsc( const G4ParticleDefinition* particle,
00626 G4double theta,
00627 G4double momentum,
00628 G4double Z )
00629 {
00630 G4double sinHalfTheta = std::sin(0.5*theta);
00631 G4double sinHalfTheta2 = sinHalfTheta*sinHalfTheta;
00632 G4double beta = CalculateParticleBeta( particle, momentum);
00633 G4double z = particle->GetPDGCharge();
00634 G4double n = CalculateZommerfeld( beta, z, Z );
00635 G4double am = CalculateAm( momentum, n, Z);
00636 G4double k = momentum/CLHEP::hbarc;
00637 G4double ch = 0.5*n/k;
00638 G4double ch2 = ch*ch;
00639 G4double xsc = ch2/(sinHalfTheta2+am)/(sinHalfTheta2+am);
00640
00641 return xsc;
00642 }
00643
00645
00646
00647
00648 inline G4double G4NuclNuclDiffuseElastic::GetRutherfordXsc( G4double theta )
00649 {
00650 G4double sinHalfTheta = std::sin(0.5*theta);
00651 G4double sinHalfTheta2 = sinHalfTheta*sinHalfTheta;
00652
00653 G4double ch2 = fRutherfordRatio*fRutherfordRatio;
00654
00655 G4double xsc = ch2/(sinHalfTheta2+fAm)/(sinHalfTheta2+fAm);
00656
00657 return xsc;
00658 }
00659
00660
00662
00663
00664
00665 inline G4double G4NuclNuclDiffuseElastic::GetCoulombTotalXsc( const G4ParticleDefinition* particle,
00666 G4double momentum, G4double Z )
00667 {
00668 G4double beta = CalculateParticleBeta( particle, momentum);
00669 G4cout<<"beta = "<<beta<<G4endl;
00670 G4double z = particle->GetPDGCharge();
00671 G4double n = CalculateZommerfeld( beta, z, Z );
00672 G4cout<<"fZomerfeld = "<<n<<G4endl;
00673 G4double am = CalculateAm( momentum, n, Z);
00674 G4cout<<"cof Am = "<<am<<G4endl;
00675 G4double k = momentum/CLHEP::hbarc;
00676 G4cout<<"k = "<<k*CLHEP::fermi<<" 1/fermi"<<G4endl;
00677 G4cout<<"k*Bohr_radius = "<<k*CLHEP::Bohr_radius<<G4endl;
00678 G4double ch = n/k;
00679 G4double ch2 = ch*ch;
00680 G4double xsc = ch2*CLHEP::pi/(am +am*am);
00681
00682 return xsc;
00683 }
00684
00686
00687
00688
00689
00690 inline G4double G4NuclNuclDiffuseElastic::GetCoulombIntegralXsc( const G4ParticleDefinition* particle,
00691 G4double momentum, G4double Z,
00692 G4double theta1, G4double theta2 )
00693 {
00694 G4double c1 = std::cos(theta1);
00695 G4cout<<"c1 = "<<c1<<G4endl;
00696 G4double c2 = std::cos(theta2);
00697 G4cout<<"c2 = "<<c2<<G4endl;
00698 G4double beta = CalculateParticleBeta( particle, momentum);
00699
00700 G4double z = particle->GetPDGCharge();
00701 G4double n = CalculateZommerfeld( beta, z, Z );
00702
00703 G4double am = CalculateAm( momentum, n, Z);
00704
00705 G4double k = momentum/CLHEP::hbarc;
00706
00707
00708 G4double ch = n/k;
00709 G4double ch2 = ch*ch;
00710 am *= 2.;
00711 G4double xsc = ch2*CLHEP::twopi*(c1-c2);
00712 xsc /= (1 - c1 + am)*(1 - c2 + am);
00713
00714 return xsc;
00715 }
00716
00718
00719
00720
00721
00722 inline G4complex G4NuclNuclDiffuseElastic::GammaLogarithm(G4complex zz)
00723 {
00724 static G4double cof[6] = { 76.18009172947146, -86.50532032941677,
00725 24.01409824083091, -1.231739572450155,
00726 0.1208650973866179e-2, -0.5395239384953e-5 } ;
00727 register G4int j;
00728 G4complex z = zz - 1.0;
00729 G4complex tmp = z + 5.5;
00730 tmp -= (z + 0.5) * std::log(tmp);
00731 G4complex ser = G4complex(1.000000000190015,0.);
00732
00733 for ( j = 0; j <= 5; j++ )
00734 {
00735 z += 1.0;
00736 ser += cof[j]/z;
00737 }
00738 return -tmp + std::log(2.5066282746310005*ser);
00739 }
00740
00742
00743
00744
00745
00746 inline G4complex G4NuclNuclDiffuseElastic::GammaLogB2n(G4complex z)
00747 {
00748 G4complex z1 = 12.*z;
00749 G4complex z2 = z*z;
00750 G4complex z3 = z2*z;
00751 G4complex z5 = z2*z3;
00752 G4complex z7 = z2*z5;
00753
00754 z3 *= 360.;
00755 z5 *= 1260.;
00756 z7 *= 1680.;
00757
00758 G4complex result = (z-0.5)*std::log(z)-z+0.5*std::log(CLHEP::twopi);
00759 result += 1./z1 - 1./z3 +1./z5 -1./z7;
00760 return result;
00761 }
00762
00764
00765
00766
00767 inline G4double G4NuclNuclDiffuseElastic::GetErf(G4double x)
00768 {
00769 G4double t, z, tmp, result;
00770
00771 z = std::fabs(x);
00772 t = 1.0/(1.0+0.5*z);
00773
00774 tmp = t*exp(-z*z-1.26551223+t*(1.00002368+t*(0.37409196+t*(0.09678418+
00775 t*(-0.18628806+t*(0.27886807+t*(-1.13520398+t*(1.48851587+
00776 t*(-0.82215223+t*0.17087277)))))))));
00777
00778 if( x >= 0.) result = 1. - tmp;
00779 else result = 1. + tmp;
00780
00781 return result;
00782 }
00783
00785
00786
00787
00788 inline G4complex G4NuclNuclDiffuseElastic::GetErfcComp(G4complex z, G4int nMax)
00789 {
00790 G4complex erfcz = 1. - GetErfComp( z, nMax);
00791 return erfcz;
00792 }
00793
00795
00796
00797
00798 inline G4complex G4NuclNuclDiffuseElastic::GetErfcSer(G4complex z, G4int nMax)
00799 {
00800 G4complex erfcz = 1. - GetErfSer( z, nMax);
00801 return erfcz;
00802 }
00803
00805
00806
00807
00808 inline G4complex G4NuclNuclDiffuseElastic::GetErfcInt(G4complex z)
00809 {
00810 G4complex erfcz = 1. - GetErfInt( z);
00811 return erfcz;
00812 }
00813
00814 inline G4double G4NuclNuclDiffuseElastic::GetLegendrePol(G4int n, G4double theta)
00815 {
00816 G4double legPol, epsilon = 1.e-6;
00817 G4double x = std::cos(theta);
00818
00819 if ( n < 0 ) legPol = 0.;
00820 else if( n == 0 ) legPol = 1.;
00821 else if( n == 1 ) legPol = x;
00822 else if( n == 2 ) legPol = (3.*x*x-1.)/2.;
00823 else if( n == 3 ) legPol = (5.*x*x*x-3.*x)/2.;
00824 else if( n == 4 ) legPol = (35.*x*x*x*x-30.*x*x+3.)/8.;
00825 else if( n == 5 ) legPol = (63.*x*x*x*x*x-70.*x*x*x+15.*x)/8.;
00826 else if( n == 6 ) legPol = (231.*x*x*x*x*x*x-315.*x*x*x*x+105.*x*x-5.)/16.;
00827 else
00828 {
00829
00830
00831 legPol = std::sqrt( 2./(n*CLHEP::pi*std::sin(theta+epsilon)) )*std::sin( (n+0.5)*theta+0.25*CLHEP::pi );
00832 }
00833 return legPol;
00834 }
00835
00836
00837
00839
00840
00841
00842 inline G4complex G4NuclNuclDiffuseElastic::TestErfcComp(G4complex z, G4int nMax)
00843 {
00844 G4complex miz = G4complex( z.imag(), -z.real() );
00845 G4complex erfcz = 1. - GetErfComp( miz, nMax);
00846 G4complex w = std::exp(-z*z)*erfcz;
00847 return w;
00848 }
00849
00851
00852
00853
00854 inline G4complex G4NuclNuclDiffuseElastic::TestErfcSer(G4complex z, G4int nMax)
00855 {
00856 G4complex miz = G4complex( z.imag(), -z.real() );
00857 G4complex erfcz = 1. - GetErfSer( miz, nMax);
00858 G4complex w = std::exp(-z*z)*erfcz;
00859 return w;
00860 }
00861
00863
00864
00865
00866 inline G4complex G4NuclNuclDiffuseElastic::TestErfcInt(G4complex z)
00867 {
00868 G4complex miz = G4complex( z.imag(), -z.real() );
00869 G4complex erfcz = 1. - GetErfInt( miz);
00870 G4complex w = std::exp(-z*z)*erfcz;
00871 return w;
00872 }
00873
00875
00876
00877
00878 inline G4complex G4NuclNuclDiffuseElastic::GetErfComp(G4complex z, G4int nMax)
00879 {
00880 G4int n;
00881 G4double n2, cofn, shny, chny, fn, gn;
00882
00883 G4double x = z.real();
00884 G4double y = z.imag();
00885
00886 G4double outRe = 0., outIm = 0.;
00887
00888 G4double twox = 2.*x;
00889 G4double twoxy = twox*y;
00890 G4double twox2 = twox*twox;
00891
00892 G4double cof1 = std::exp(-x*x)/CLHEP::pi;
00893
00894 G4double cos2xy = std::cos(twoxy);
00895 G4double sin2xy = std::sin(twoxy);
00896
00897 G4double twoxcos2xy = twox*cos2xy;
00898 G4double twoxsin2xy = twox*sin2xy;
00899
00900 for( n = 1; n <= nMax; n++)
00901 {
00902 n2 = n*n;
00903
00904 cofn = std::exp(-0.5*n2)/(n2+twox2);
00905
00906 chny = std::cosh(n*y);
00907 shny = std::sinh(n*y);
00908
00909 fn = twox - twoxcos2xy*chny + n*sin2xy*shny;
00910 gn = twoxsin2xy*chny + n*cos2xy*shny;
00911
00912 fn *= cofn;
00913 gn *= cofn;
00914
00915 outRe += fn;
00916 outIm += gn;
00917 }
00918 outRe *= 2*cof1;
00919 outIm *= 2*cof1;
00920
00921 if(std::abs(x) < 0.0001)
00922 {
00923 outRe += GetErf(x);
00924 outIm += cof1*y;
00925 }
00926 else
00927 {
00928 outRe += GetErf(x) + cof1*(1-cos2xy)/twox;
00929 outIm += cof1*sin2xy/twox;
00930 }
00931 return G4complex(outRe, outIm);
00932 }
00933
00935
00936
00937
00938 inline G4complex G4NuclNuclDiffuseElastic::GetErfSer(G4complex z, G4int nMax)
00939 {
00940 G4int n;
00941 G4double a =1., b = 1., tmp;
00942 G4complex sum = z, d = z;
00943
00944 for( n = 1; n <= nMax; n++)
00945 {
00946 a *= 2.;
00947 b *= 2.*n +1.;
00948 d *= z*z;
00949
00950 tmp = a/b;
00951
00952 sum += tmp*d;
00953 }
00954 sum *= 2.*std::exp(-z*z)/std::sqrt(CLHEP::pi);
00955
00956 return sum;
00957 }
00958
00960
00961 inline G4double G4NuclNuclDiffuseElastic::GetExpCos(G4double x)
00962 {
00963 G4double result;
00964
00965 result = std::exp(x*x-fReZ*fReZ);
00966 result *= std::cos(2.*x*fReZ);
00967 return result;
00968 }
00969
00971
00972 inline G4double G4NuclNuclDiffuseElastic::GetExpSin(G4double x)
00973 {
00974 G4double result;
00975
00976 result = std::exp(x*x-fReZ*fReZ);
00977 result *= std::sin(2.*x*fReZ);
00978 return result;
00979 }
00980
00981
00982
00984
00985
00986
00987 inline G4complex G4NuclNuclDiffuseElastic::GetErfInt(G4complex z)
00988 {
00989 G4double outRe, outIm;
00990
00991 G4double x = z.real();
00992 G4double y = z.imag();
00993 fReZ = x;
00994
00995 G4Integrator<G4NuclNuclDiffuseElastic,G4double(G4NuclNuclDiffuseElastic::*)(G4double)> integral;
00996
00997 outRe = integral.Legendre96(this,&G4NuclNuclDiffuseElastic::GetExpSin, 0., y );
00998 outIm = integral.Legendre96(this,&G4NuclNuclDiffuseElastic::GetExpCos, 0., y );
00999
01000 outRe *= 2./sqrt(CLHEP::pi);
01001 outIm *= 2./sqrt(CLHEP::pi);
01002
01003 outRe += GetErf(x);
01004
01005 return G4complex(outRe, outIm);
01006 }
01007
01009
01010
01011
01012 inline G4double G4NuclNuclDiffuseElastic::GetCint(G4double x)
01013 {
01014 G4double out;
01015
01016 G4Integrator<G4NuclNuclDiffuseElastic,G4double(G4NuclNuclDiffuseElastic::*)(G4double)> integral;
01017
01018 out= integral.Legendre96(this,&G4NuclNuclDiffuseElastic::GetCosHaPit2, 0., x );
01019
01020 return out;
01021 }
01022
01024
01025
01026
01027 inline G4double G4NuclNuclDiffuseElastic::GetSint(G4double x)
01028 {
01029 G4double out;
01030
01031 G4Integrator<G4NuclNuclDiffuseElastic,G4double(G4NuclNuclDiffuseElastic::*)(G4double)> integral;
01032
01033 out= integral.Legendre96(this,&G4NuclNuclDiffuseElastic::GetSinHaPit2, 0., x );
01034
01035 return out;
01036 }
01037
01038
01040
01041
01042
01043 inline G4complex G4NuclNuclDiffuseElastic::CoulombAmplitude(G4double theta)
01044 {
01045 G4complex ca;
01046
01047 G4double sinHalfTheta = std::sin(0.5*theta);
01048 G4double sinHalfTheta2 = sinHalfTheta*sinHalfTheta;
01049 sinHalfTheta2 += fAm;
01050
01051 G4double order = 2.*fCoulombPhase0 - fZommerfeld*std::log(sinHalfTheta2);
01052 G4complex z = G4complex(0., order);
01053 ca = std::exp(z);
01054
01055 ca *= -fZommerfeld/(2.*fWaveVector*sinHalfTheta2);
01056
01057 return ca;
01058 }
01059
01061
01062
01063
01064 inline G4double G4NuclNuclDiffuseElastic::CoulombAmplitudeMod2(G4double theta)
01065 {
01066 G4complex ca = CoulombAmplitude(theta);
01067 G4double out = ca.real()*ca.real() + ca.imag()*ca.imag();
01068
01069 return out;
01070 }
01071
01073
01074
01075
01076
01077 inline void G4NuclNuclDiffuseElastic::CalculateCoulombPhaseZero()
01078 {
01079 G4complex z = G4complex(1,fZommerfeld);
01080
01081 G4complex gammalog = GammaLogB2n(z);
01082 fCoulombPhase0 = gammalog.imag();
01083 }
01084
01086
01087
01088
01089
01090 inline G4double G4NuclNuclDiffuseElastic::CalculateCoulombPhase(G4int n)
01091 {
01092 G4complex z = G4complex(1. + n, fZommerfeld);
01093
01094 G4complex gammalog = GammaLogB2n(z);
01095 return gammalog.imag();
01096 }
01097
01098
01100
01101
01102
01103
01104 inline void G4NuclNuclDiffuseElastic::CalculateRutherfordAnglePar()
01105 {
01106 fHalfRutThetaTg = fZommerfeld/fProfileLambda;
01107 fRutherfordTheta = 2.*std::atan(fHalfRutThetaTg);
01108 fHalfRutThetaTg2 = fHalfRutThetaTg*fHalfRutThetaTg;
01109
01110
01111 }
01112
01114
01115
01116
01117 inline G4double G4NuclNuclDiffuseElastic::ProfileNear(G4double theta)
01118 {
01119 G4double dTheta = fRutherfordTheta - theta;
01120 G4double result = 0., argument = 0.;
01121
01122 if(std::abs(dTheta) < 0.001) result = fProfileAlpha*fProfileDelta;
01123 else
01124 {
01125 argument = fProfileDelta*dTheta;
01126 result = CLHEP::pi*argument*std::exp(fProfileAlpha*argument);
01127 result /= std::sinh(CLHEP::pi*argument);
01128 result -= 1.;
01129 result /= dTheta;
01130 }
01131 return result;
01132 }
01133
01135
01136
01137
01138 inline G4double G4NuclNuclDiffuseElastic::ProfileFar(G4double theta)
01139 {
01140 G4double dTheta = fRutherfordTheta + theta;
01141 G4double argument = fProfileDelta*dTheta;
01142
01143 G4double result = CLHEP::pi*argument*std::exp(fProfileAlpha*argument);
01144 result /= std::sinh(CLHEP::pi*argument);
01145 result /= dTheta;
01146
01147 return result;
01148 }
01149
01151
01152
01153
01154 inline G4double G4NuclNuclDiffuseElastic::Profile(G4double theta)
01155 {
01156 G4double dTheta = fRutherfordTheta - theta;
01157 G4double result = 0., argument = 0.;
01158
01159 if(std::abs(dTheta) < 0.001) result = 1.;
01160 else
01161 {
01162 argument = fProfileDelta*dTheta;
01163 result = CLHEP::pi*argument;
01164 result /= std::sinh(CLHEP::pi*argument);
01165 }
01166 return result;
01167 }
01168
01170
01171
01172
01173 inline G4complex G4NuclNuclDiffuseElastic::PhaseNear(G4double theta)
01174 {
01175 G4double twosigma = 2.*fCoulombPhase0;
01176 twosigma -= fZommerfeld*std::log(fHalfRutThetaTg2/(1.+fHalfRutThetaTg2));
01177 twosigma += fRutherfordTheta*fZommerfeld/fHalfRutThetaTg - CLHEP::halfpi;
01178 twosigma -= fProfileLambda*theta - 0.25*CLHEP::pi;
01179
01180 twosigma *= fCofPhase;
01181
01182 G4complex z = G4complex(0., twosigma);
01183
01184 return std::exp(z);
01185 }
01186
01188
01189
01190
01191 inline G4complex G4NuclNuclDiffuseElastic::PhaseFar(G4double theta)
01192 {
01193 G4double twosigma = 2.*fCoulombPhase0;
01194 twosigma -= fZommerfeld*std::log(fHalfRutThetaTg2/(1.+fHalfRutThetaTg2));
01195 twosigma += fRutherfordTheta*fZommerfeld/fHalfRutThetaTg - CLHEP::halfpi;
01196 twosigma += fProfileLambda*theta - 0.25*CLHEP::pi;
01197
01198 twosigma *= fCofPhase;
01199
01200 G4complex z = G4complex(0., twosigma);
01201
01202 return std::exp(z);
01203 }
01204
01206
01207
01208
01209
01210 inline G4complex G4NuclNuclDiffuseElastic::GammaLess(G4double theta)
01211 {
01212 G4double sinThetaR = 2.*fHalfRutThetaTg/(1. + fHalfRutThetaTg2);
01213 G4double cosHalfThetaR2 = 1./(1. + fHalfRutThetaTg2);
01214
01215 G4double u = std::sqrt(0.5*fProfileLambda/sinThetaR);
01216 G4double kappa = u/std::sqrt(CLHEP::pi);
01217 G4double dTheta = theta - fRutherfordTheta;
01218 u *= dTheta;
01219 G4double u2 = u*u;
01220 G4double u2m2p3 = u2*2./3.;
01221
01222 G4complex im = G4complex(0.,1.);
01223 G4complex order = G4complex(u,u);
01224 order /= std::sqrt(2.);
01225
01226 G4complex gamma = CLHEP::pi*kappa*GetErfcInt(-order)*std::exp(im*(u*u+0.25*CLHEP::pi));
01227 G4complex a0 = 0.5*(1. + 4.*(1.+im*u2)*cosHalfThetaR2/3.)/sinThetaR;
01228 G4complex a1 = 0.5*(1. + 2.*(1.+im*u2m2p3)*cosHalfThetaR2)/sinThetaR;
01229 G4complex out = gamma*(1. - a1*dTheta) - a0;
01230
01231 return out;
01232 }
01233
01235
01236
01237
01238 inline G4complex G4NuclNuclDiffuseElastic::GammaMore(G4double theta)
01239 {
01240 G4double sinThetaR = 2.*fHalfRutThetaTg/(1. + fHalfRutThetaTg2);
01241 G4double cosHalfThetaR2 = 1./(1. + fHalfRutThetaTg2);
01242
01243 G4double u = std::sqrt(0.5*fProfileLambda/sinThetaR);
01244 G4double kappa = u/std::sqrt(CLHEP::pi);
01245 G4double dTheta = theta - fRutherfordTheta;
01246 u *= dTheta;
01247 G4double u2 = u*u;
01248 G4double u2m2p3 = u2*2./3.;
01249
01250 G4complex im = G4complex(0.,1.);
01251 G4complex order = G4complex(u,u);
01252 order /= std::sqrt(2.);
01253 G4complex gamma = CLHEP::pi*kappa*GetErfcInt(order)*std::exp(im*(u*u+0.25*CLHEP::pi));
01254 G4complex a0 = 0.5*(1. + 4.*(1.+im*u2)*cosHalfThetaR2/3.)/sinThetaR;
01255 G4complex a1 = 0.5*(1. + 2.*(1.+im*u2m2p3)*cosHalfThetaR2)/sinThetaR;
01256 G4complex out = -gamma*(1. - a1*dTheta) - a0;
01257
01258 return out;
01259 }
01260
01262
01263
01264
01265 inline G4complex G4NuclNuclDiffuseElastic::AmplitudeNear(G4double theta)
01266 {
01267 G4double kappa = std::sqrt(0.5*fProfileLambda/std::sin(theta)/CLHEP::pi);
01268 G4complex out = G4complex(kappa/fWaveVector,0.);
01269
01270 out *= PhaseNear(theta);
01271
01272 if( theta <= fRutherfordTheta )
01273 {
01274 out *= GammaLess(theta) + ProfileNear(theta);
01275
01276 out += CoulombAmplitude(theta);
01277 }
01278 else
01279 {
01280 out *= GammaMore(theta) + ProfileNear(theta);
01281
01282 }
01283 return out;
01284 }
01285
01287
01288
01289
01290 inline G4complex G4NuclNuclDiffuseElastic::AmplitudeFar(G4double theta)
01291 {
01292 G4double kappa = std::sqrt(0.5*fProfileLambda/std::sin(theta)/CLHEP::pi);
01293 G4complex out = G4complex(kappa/fWaveVector,0.);
01294 out *= ProfileFar(theta);
01295 out *= PhaseFar(theta);
01296 return out;
01297 }
01298
01299
01301
01302
01303
01304 inline G4complex G4NuclNuclDiffuseElastic::Amplitude(G4double theta)
01305 {
01306
01307 G4complex out = AmplitudeNear(theta) + fCofFar*AmplitudeFar(theta);
01308
01309
01310 return out;
01311 }
01312
01314
01315
01316
01317 inline G4double G4NuclNuclDiffuseElastic::AmplitudeMod2(G4double theta)
01318 {
01319 G4complex out = Amplitude(theta);
01320 G4double mod2 = out.real()*out.real() + out.imag()*out.imag();
01321 return mod2;
01322 }
01323
01325
01326
01327
01328 inline G4complex G4NuclNuclDiffuseElastic::AmplitudeSim(G4double theta)
01329 {
01330 G4double sinThetaR = 2.*fHalfRutThetaTg/(1. + fHalfRutThetaTg2);
01331 G4double dTheta = 0.5*(theta - fRutherfordTheta);
01332 G4double sindTheta = std::sin(dTheta);
01333 G4double persqrt2 = std::sqrt(0.5);
01334
01335 G4complex order = G4complex(persqrt2,persqrt2);
01336 order *= std::sqrt(0.5*fProfileLambda/sinThetaR)*2.*sindTheta;
01337
01338
01339 G4complex out;
01340
01341 if ( theta <= fRutherfordTheta )
01342 {
01343 out = 1. - 0.5*GetErfcInt(-order)*ProfileNear(theta);
01344 }
01345 else
01346 {
01347 out = 0.5*GetErfcInt(order)*ProfileNear(theta);
01348 }
01349
01350 out *= CoulombAmplitude(theta);
01351
01352 return out;
01353 }
01354
01356
01357
01358
01359 inline G4double G4NuclNuclDiffuseElastic::GetRatioSim(G4double theta)
01360 {
01361 G4double sinThetaR = 2.*fHalfRutThetaTg/(1. + fHalfRutThetaTg2);
01362 G4double dTheta = 0.5*(theta - fRutherfordTheta);
01363 G4double sindTheta = std::sin(dTheta);
01364
01365 G4double order = std::sqrt(fProfileLambda/sinThetaR/CLHEP::pi)*2.*sindTheta;
01366
01367 G4double cosFresnel = 0.5 - GetCint(order);
01368 G4double sinFresnel = 0.5 - GetSint(order);
01369
01370 G4double out = 0.5*( cosFresnel*cosFresnel + sinFresnel*sinFresnel );
01371
01372 return out;
01373 }
01374
01376
01377
01378
01379 inline G4double G4NuclNuclDiffuseElastic::GetRatioGen(G4double theta)
01380 {
01381 G4double sinThetaR = 2.*fHalfRutThetaTg/(1. + fHalfRutThetaTg2);
01382 G4double dTheta = 0.5*(theta - fRutherfordTheta);
01383 G4double sindTheta = std::sin(dTheta);
01384
01385 G4double prof = Profile(theta);
01386 G4double prof2 = prof*prof;
01387
01388 G4double order = std::sqrt(fProfileLambda/sinThetaR/CLHEP::pi)*2.*sindTheta;
01389
01390 order = std::abs(order);
01391
01392
01393 G4double cosFresnel = GetCint(order);
01394 G4double sinFresnel = GetSint(order);
01395
01396 G4double out;
01397
01398 if ( theta <= fRutherfordTheta )
01399 {
01400 out = 1. + 0.5*( (0.5-cosFresnel)*(0.5-cosFresnel)+(0.5-sinFresnel)*(0.5-sinFresnel) )*prof2;
01401 out += ( cosFresnel + sinFresnel - 1. )*prof;
01402 }
01403 else
01404 {
01405 out = 0.5*( (0.5-cosFresnel)*(0.5-cosFresnel)+(0.5-sinFresnel)*(0.5-sinFresnel) )*prof2;
01406 }
01407
01408 return out;
01409 }
01410
01412
01413
01414
01415 inline G4double G4NuclNuclDiffuseElastic::GetFresnelDiffuseXsc(G4double theta)
01416 {
01417 G4double ratio = GetRatioGen(theta);
01418 G4double ruthXsc = GetRutherfordXsc(theta);
01419 G4double xsc = ratio*ruthXsc;
01420 return xsc;
01421 }
01422
01424
01425
01426
01427 inline G4double G4NuclNuclDiffuseElastic::GetFresnelIntegrandXsc(G4double alpha)
01428 {
01429 G4double theta = std::sqrt(alpha);
01430 G4double xsc = GetFresnelDiffuseXsc(theta);
01431 return xsc;
01432 }
01433
01434
01435
01437
01438
01439
01440 inline G4double G4NuclNuclDiffuseElastic::AmplitudeSimMod2(G4double theta)
01441 {
01442 G4complex out = AmplitudeSim(theta);
01443 G4double mod2 = out.real()*out.real() + out.imag()*out.imag();
01444 return mod2;
01445 }
01446
01448
01449
01450
01451 inline G4complex G4NuclNuclDiffuseElastic::AmplitudeGla(G4double theta)
01452 {
01453 G4int n;
01454 G4double T12b, b, b2;
01455 G4complex out = G4complex(0.,0.), shiftC, shiftN;
01456 G4complex im = G4complex(0.,1.);
01457
01458 for( n = 0; n < fMaxL; n++)
01459 {
01460 shiftC = std::exp( im*2.*CalculateCoulombPhase(n) );
01461
01462 b = ( std::sqrt( G4double(n*(n+1)) ) )/fWaveVector;
01463 b2 = b*b;
01464 T12b = fSumSigma*std::exp(-b2/fNuclearRadiusSquare)/CLHEP::pi/fNuclearRadiusSquare;
01465 shiftN = std::exp( -0.5*(1.-im*fEtaRatio)*T12b ) - 1.;
01466 out += (2.*n+1.)*shiftC*shiftN*GetLegendrePol(n, theta);
01467 }
01468 out /= 2.*im*fWaveVector;
01469 out += CoulombAmplitude(theta);
01470 return out;
01471 }
01472
01474
01475
01476
01477 inline G4double G4NuclNuclDiffuseElastic::AmplitudeGlaMod2(G4double theta)
01478 {
01479 G4complex out = AmplitudeGla(theta);
01480 G4double mod2 = out.real()*out.real() + out.imag()*out.imag();
01481 return mod2;
01482 }
01483
01484
01486
01487
01488
01489 inline G4complex G4NuclNuclDiffuseElastic::AmplitudeGG(G4double theta)
01490 {
01491 G4int n;
01492 G4double T12b, a, aTemp, b2, sinThetaH = std::sin(0.5*theta);
01493 G4double sinThetaH2 = sinThetaH*sinThetaH;
01494 G4complex out = G4complex(0.,0.);
01495 G4complex im = G4complex(0.,1.);
01496
01497 a = -fSumSigma/CLHEP::twopi/fNuclearRadiusSquare;
01498 b2 = fWaveVector*fWaveVector*fNuclearRadiusSquare*sinThetaH2;
01499
01500 aTemp = a;
01501
01502 for( n = 1; n < fMaxL; n++)
01503 {
01504 T12b = aTemp*std::exp(-b2/n)/n;
01505 aTemp *= a;
01506 out += T12b;
01507 G4cout<<"out = "<<out<<G4endl;
01508 }
01509 out *= -4.*im*fWaveVector/CLHEP::pi;
01510 out += CoulombAmplitude(theta);
01511 return out;
01512 }
01513
01515
01516
01517
01518 inline G4double G4NuclNuclDiffuseElastic::AmplitudeGGMod2(G4double theta)
01519 {
01520 G4complex out = AmplitudeGG(theta);
01521 G4double mod2 = out.real()*out.real() + out.imag()*out.imag();
01522 return mod2;
01523 }
01524
01525
01527
01528
01529
01530
01531 inline void G4NuclNuclDiffuseElastic::InitParameters(const G4ParticleDefinition* theParticle,
01532 G4double partMom, G4double Z, G4double A)
01533 {
01534 fAtomicNumber = Z;
01535 fAtomicWeight = A;
01536
01537 fNuclearRadius2 = CalculateNuclearRad(fAtomicWeight);
01538 G4double A1 = G4double( theParticle->GetBaryonNumber() );
01539 fNuclearRadius1 = CalculateNuclearRad(A1);
01540
01541 fNuclearRadius = fNuclearRadius1 + fNuclearRadius2;
01542
01543 G4double a = 0.;
01544 G4double z = theParticle->GetPDGCharge();
01545 G4double m1 = theParticle->GetPDGMass();
01546
01547 fWaveVector = partMom/CLHEP::hbarc;
01548
01549 G4double lambda = fCofLambda*fWaveVector*fNuclearRadius;
01550 G4cout<<"kR = "<<lambda<<G4endl;
01551
01552 if( z )
01553 {
01554 a = partMom/m1;
01555 fBeta = a/std::sqrt(1+a*a);
01556 fZommerfeld = CalculateZommerfeld( fBeta, z, fAtomicNumber);
01557 fRutherfordRatio = fZommerfeld/fWaveVector;
01558 fAm = CalculateAm( partMom, fZommerfeld, fAtomicNumber);
01559 }
01560 G4cout<<"fZommerfeld = "<<fZommerfeld<<G4endl;
01561 fProfileLambda = lambda;
01562 G4cout<<"fProfileLambda = "<<fProfileLambda<<G4endl;
01563 fProfileDelta = fCofDelta*fProfileLambda;
01564 fProfileAlpha = fCofAlpha*fProfileLambda;
01565
01566 CalculateCoulombPhaseZero();
01567 CalculateRutherfordAnglePar();
01568
01569 return;
01570 }
01572
01573
01574
01575
01576 inline void G4NuclNuclDiffuseElastic::InitDynParameters(const G4ParticleDefinition* theParticle,
01577 G4double partMom)
01578 {
01579 G4double a = 0.;
01580 G4double z = theParticle->GetPDGCharge();
01581 G4double m1 = theParticle->GetPDGMass();
01582
01583 fWaveVector = partMom/CLHEP::hbarc;
01584
01585 G4double lambda = fCofLambda*fWaveVector*fNuclearRadius;
01586
01587 if( z )
01588 {
01589 a = partMom/m1;
01590 fBeta = a/std::sqrt(1+a*a);
01591 fZommerfeld = CalculateZommerfeld( fBeta, z, fAtomicNumber);
01592 fRutherfordRatio = fZommerfeld/fWaveVector;
01593 fAm = CalculateAm( partMom, fZommerfeld, fAtomicNumber);
01594 }
01595 fProfileLambda = lambda;
01596 fProfileDelta = fCofDelta*fProfileLambda;
01597 fProfileAlpha = fCofAlpha*fProfileLambda;
01598
01599 CalculateCoulombPhaseZero();
01600 CalculateRutherfordAnglePar();
01601
01602 return;
01603 }
01604
01605
01607
01608
01609
01610
01611 inline void G4NuclNuclDiffuseElastic::InitParametersGla(const G4DynamicParticle* aParticle,
01612 G4double partMom, G4double Z, G4double A)
01613 {
01614 fAtomicNumber = Z;
01615 fAtomicWeight = A;
01616
01617 fNuclearRadius2 = CalculateNuclearRad(fAtomicWeight);
01618 G4double A1 = G4double( aParticle->GetDefinition()->GetBaryonNumber() );
01619 fNuclearRadius1 = CalculateNuclearRad(A1);
01620 fNuclearRadiusSquare = fNuclearRadius1*fNuclearRadius1+fNuclearRadius2*fNuclearRadius2;
01621
01622
01623 G4double a = 0., kR12;
01624 G4double z = aParticle->GetDefinition()->GetPDGCharge();
01625 G4double m1 = aParticle->GetDefinition()->GetPDGMass();
01626
01627 fWaveVector = partMom/CLHEP::hbarc;
01628
01629 G4double pN = A1 - z;
01630 if( pN < 0. ) pN = 0.;
01631
01632 G4double tN = A - Z;
01633 if( tN < 0. ) tN = 0.;
01634
01635 G4double pTkin = aParticle->GetKineticEnergy();
01636 pTkin /= A1;
01637
01638
01639 fSumSigma = (Z*z+pN*tN)*GetHadronNucleonXscNS(theProton, pTkin, theProton) +
01640 (z*tN+pN*Z)*GetHadronNucleonXscNS(theProton, pTkin, theNeutron);
01641
01642 G4cout<<"fSumSigma = "<<fSumSigma/CLHEP::millibarn<<" mb"<<G4endl;
01643 G4cout<<"pi*R2 = "<<CLHEP::pi*fNuclearRadiusSquare/CLHEP::millibarn<<" mb"<<G4endl;
01644 kR12 = fWaveVector*std::sqrt(fNuclearRadiusSquare);
01645 G4cout<<"k*sqrt(R2) = "<<kR12<<" "<<G4endl;
01646 fMaxL = (G4int(kR12)+1)*4;
01647 G4cout<<"fMaxL = "<<fMaxL<<" "<<G4endl;
01648
01649 if( z )
01650 {
01651 a = partMom/m1;
01652 fBeta = a/std::sqrt(1+a*a);
01653 fZommerfeld = CalculateZommerfeld( fBeta, z, fAtomicNumber);
01654 fAm = CalculateAm( partMom, fZommerfeld, fAtomicNumber);
01655 }
01656
01657 CalculateCoulombPhaseZero();
01658
01659
01660 return;
01661 }
01662
01663
01665
01666
01667
01668
01669
01670 inline G4double
01671 G4NuclNuclDiffuseElastic::GetHadronNucleonXscNS( G4ParticleDefinition* pParticle,
01672 G4double pTkin,
01673 G4ParticleDefinition* tParticle)
01674 {
01675 G4double xsection(0), A0, B0;
01676 G4double hpXsc(0);
01677 G4double hnXsc(0);
01678
01679
01680 G4double targ_mass = tParticle->GetPDGMass();
01681 G4double proj_mass = pParticle->GetPDGMass();
01682
01683 G4double proj_energy = proj_mass + pTkin;
01684 G4double proj_momentum = std::sqrt(pTkin*(pTkin+2*proj_mass));
01685
01686 G4double sMand = CalcMandelstamS ( proj_mass , targ_mass , proj_momentum );
01687
01688 sMand /= CLHEP::GeV*CLHEP::GeV;
01689 proj_momentum /= CLHEP::GeV;
01690 proj_energy /= CLHEP::GeV;
01691 proj_mass /= CLHEP::GeV;
01692 G4double logS = std::log(sMand);
01693
01694
01695
01696
01697
01698
01699 if( proj_momentum >= 1.2 )
01700 {
01701 fEtaRatio = 0.13*(logS - 5.8579332)*std::pow(sMand,-0.18);
01702 }
01703 else if( proj_momentum >= 0.6 )
01704 {
01705 fEtaRatio = -75.5*(std::pow(proj_momentum,0.25)-0.95)/
01706 (std::pow(3*proj_momentum,2.2)+1);
01707 }
01708 else
01709 {
01710 fEtaRatio = 15.5*proj_momentum/(27*proj_momentum*proj_momentum*proj_momentum+2);
01711 }
01712 G4cout<<"fEtaRatio = "<<fEtaRatio<<G4endl;
01713
01714
01715
01716 if( proj_momentum >= 10. )
01717
01718 {
01719
01720
01721
01722
01723 if( proj_momentum >= 10.)
01724 {
01725 B0 = 7.5;
01726 A0 = 100. - B0*std::log(3.0e7);
01727
01728 xsection = A0 + B0*std::log(proj_energy) - 11
01729 + 103*std::pow(2*0.93827*proj_energy + proj_mass*proj_mass+
01730 0.93827*0.93827,-0.165);
01731 }
01732 }
01733 else
01734 {
01735 if(pParticle == tParticle)
01736 {
01737 if( proj_momentum < 0.73 )
01738 {
01739 hnXsc = 23 + 50*( std::pow( std::log(0.73/proj_momentum), 3.5 ) );
01740 }
01741 else if( proj_momentum < 1.05 )
01742 {
01743 hnXsc = 23 + 40*(std::log(proj_momentum/0.73))*
01744 (std::log(proj_momentum/0.73));
01745 }
01746 else
01747 {
01748 hnXsc = 39.0 +
01749 75*(proj_momentum - 1.2)/(std::pow(proj_momentum,3.0) + 0.15);
01750 }
01751 xsection = hnXsc;
01752 }
01753 else
01754 {
01755 if( proj_momentum < 0.8 )
01756 {
01757 hpXsc = 33+30*std::pow(std::log(proj_momentum/1.3),4.0);
01758 }
01759 else if( proj_momentum < 1.4 )
01760 {
01761 hpXsc = 33+30*std::pow(std::log(proj_momentum/0.95),2.0);
01762 }
01763 else
01764 {
01765 hpXsc = 33.3+
01766 20.8*(std::pow(proj_momentum,2.0)-1.35)/
01767 (std::pow(proj_momentum,2.50)+0.95);
01768 }
01769 xsection = hpXsc;
01770 }
01771 }
01772 xsection *= CLHEP::millibarn;
01773 G4cout<<"xsection = "<<xsection/CLHEP::millibarn<<" mb"<<G4endl;
01774 return xsection;
01775 }
01776
01778
01779
01780
01781 inline G4double G4NuclNuclDiffuseElastic::CalcMandelstamS( const G4double mp ,
01782 const G4double mt ,
01783 const G4double Plab )
01784 {
01785 G4double Elab = std::sqrt ( mp * mp + Plab * Plab );
01786 G4double sMand = mp*mp + mt*mt + 2*Elab*mt ;
01787
01788 return sMand;
01789 }
01790
01791
01792
01793 #endif