Loading...
Searching...
No Matches
G4NuclNuclDiffuseElastic Class Reference
#include <G4NuclNuclDiffuseElastic.hh>
Inheritance diagram for G4NuclNuclDiffuseElastic:
Additional Inherited Members | |
| Protected Member Functions inherited from G4HadronicInteraction | |
| void | SetModelName (const G4String &nam) |
| G4bool | IsBlocked () const |
| void | Block () |
| Protected Attributes inherited from G4HadronicInteraction | |
| G4HadFinalState | theParticleChange |
| G4int | verboseLevel |
| G4double | theMinEnergy |
| G4double | theMaxEnergy |
| G4bool | isBlocked |
Detailed Description
Definition at line 59 of file G4NuclNuclDiffuseElastic.hh.
Constructor & Destructor Documentation
◆ G4NuclNuclDiffuseElastic()
| G4NuclNuclDiffuseElastic::G4NuclNuclDiffuseElastic | ( | ) |
Definition at line 67 of file G4NuclNuclDiffuseElastic.cc.
69{
70 SetMinEnergy( 50*MeV );
71 SetMaxEnergy( 1.*TeV );
72 verboseLevel = 0;
73 lowEnergyRecoilLimit = 100.*keV;
74 lowEnergyLimitQ = 0.0*GeV;
75 lowEnergyLimitHE = 0.0*GeV;
76 lowestEnergyLimit= 0.0*keV;
77 plabLowLimit = 20.0*MeV;
78
79 theProton = G4Proton::Proton();
80 theNeutron = G4Neutron::Neutron();
81 theDeuteron = G4Deuteron::Deuteron();
82 theAlpha = G4Alpha::Alpha();
83 thePionPlus = G4PionPlus::PionPlus();
84 thePionMinus= G4PionMinus::PionMinus();
85
86 fEnergyBin = 200;
87 fAngleBin = 200;
88
90 fAngleTable = 0;
91
92 fParticle = 0;
93 fWaveVector = 0.;
94 fAtomicWeight = 0.;
95 fAtomicNumber = 0.;
96 fNuclearRadius = 0.;
97 fBeta = 0.;
98 fZommerfeld = 0.;
99 fAm = 0.;
100 fAddCoulomb = false;
101 // Ranges of angle table relative to current Rutherford (Coulomb grazing) angle
102
103 // Empirical parameters
104
105 fCofAlphaMax = 1.5;
106 fCofAlphaCoulomb = 0.5;
107
108 fProfileDelta = 1.;
109 fProfileAlpha = 0.5;
110
111 fCofLambda = 1.0;
112 fCofDelta = 0.04;
113 fCofAlpha = 0.095;
114
115 fNuclearRadius1 = fNuclearRadius2 = fNuclearRadiusSquare = fNuclearRadiusCof
116 = fRutherfordRatio = fCoulombPhase0 = fHalfRutThetaTg = fHalfRutThetaTg2
117 = fRutherfordTheta = fProfileLambda = fCofPhase = fCofFar = fCofAlphaMax
118 = fCofAlphaCoulomb = fSumSigma = fEtaRatio = fReZ = 0.0;
119 fMaxL = 0;
120
121}
Definition: G4HadronElastic.hh:51
void SetMinEnergy(G4double anEnergy)
Definition: G4HadronicInteraction.hh:90
void SetMaxEnergy(const G4double anEnergy)
Definition: G4HadronicInteraction.hh:103
Definition: G4PhysicsLogVector.hh:60
Referenced by BuildAngleTable(), GetCint(), GetErfInt(), GetSint(), IntegralElasticProb(), SampleThetaCMS(), and TestAngleTable().
◆ ~G4NuclNuclDiffuseElastic()
|
virtual |
Definition at line 127 of file G4NuclNuclDiffuseElastic.cc.
128{
129 if(fEnergyVector) delete fEnergyVector;
130
131 if( fAngleTable )
132 {
133 fAngleTable->clearAndDestroy();
134 delete fAngleTable ;
135 }
136}
void clearAndDestroy()
Member Function Documentation
◆ Amplitude()
Definition at line 1304 of file G4NuclNuclDiffuseElastic.hh.
1305{
1306
1308 // G4complex out = AmplitudeNear(theta);
1309 // G4complex out = AmplitudeFar(theta);
1310 return out;
1311}
G4complex AmplitudeFar(G4double theta)
Definition: G4NuclNuclDiffuseElastic.hh:1290
G4complex AmplitudeNear(G4double theta)
Definition: G4NuclNuclDiffuseElastic.hh:1265
Referenced by AmplitudeMod2().
◆ AmplitudeFar()
Definition at line 1290 of file G4NuclNuclDiffuseElastic.hh.
1291{
1292 G4double kappa = std::sqrt(0.5*fProfileLambda/std::sin(theta)/CLHEP::pi);
1294 out *= ProfileFar(theta);
1295 out *= PhaseFar(theta);
1296 return out;
1297}
G4complex PhaseFar(G4double theta)
Definition: G4NuclNuclDiffuseElastic.hh:1191
G4double ProfileFar(G4double theta)
Definition: G4NuclNuclDiffuseElastic.hh:1138
Referenced by Amplitude().
◆ AmplitudeGG()
Definition at line 1489 of file G4NuclNuclDiffuseElastic.hh.
1490{
1492 G4double T12b, a, aTemp, b2, sinThetaH = std::sin(0.5*theta);
1493 G4double sinThetaH2 = sinThetaH*sinThetaH;
1496
1497 a = -fSumSigma/CLHEP::twopi/fNuclearRadiusSquare;
1498 b2 = fWaveVector*fWaveVector*fNuclearRadiusSquare*sinThetaH2;
1499
1500 aTemp = a;
1501
1503 {
1504 T12b = aTemp*std::exp(-b2/n)/n;
1505 aTemp *= a;
1506 out += T12b;
1508 }
1509 out *= -4.*im*fWaveVector/CLHEP::pi;
1510 out += CoulombAmplitude(theta);
1511 return out;
1512}
G4DLLIMPORT std::ostream G4cout
G4complex CoulombAmplitude(G4double theta)
Definition: G4NuclNuclDiffuseElastic.hh:1043
Referenced by AmplitudeGGMod2().
◆ AmplitudeGGMod2()
Definition at line 1518 of file G4NuclNuclDiffuseElastic.hh.
1519{
1521 G4double mod2 = out.real()*out.real() + out.imag()*out.imag();
1522 return mod2;
1523}
G4complex AmplitudeGG(G4double theta)
Definition: G4NuclNuclDiffuseElastic.hh:1489
◆ AmplitudeGla()
Definition at line 1451 of file G4NuclNuclDiffuseElastic.hh.
1452{
1457
1459 {
1460 shiftC = std::exp( im*2.*CalculateCoulombPhase(n) );
1461 // b = ( fZommerfeld + std::sqrt( fZommerfeld*fZommerfeld + n*(n+1) ) )/fWaveVector;
1462 b = ( std::sqrt( G4double(n*(n+1)) ) )/fWaveVector;
1463 b2 = b*b;
1464 T12b = fSumSigma*std::exp(-b2/fNuclearRadiusSquare)/CLHEP::pi/fNuclearRadiusSquare;
1465 shiftN = std::exp( -0.5*(1.-im*fEtaRatio)*T12b ) - 1.;
1467 }
1468 out /= 2.*im*fWaveVector;
1469 out += CoulombAmplitude(theta);
1470 return out;
1471}
G4double CalculateCoulombPhase(G4int n)
Definition: G4NuclNuclDiffuseElastic.hh:1090
G4double GetLegendrePol(G4int n, G4double x)
Definition: G4NuclNuclDiffuseElastic.hh:814
Referenced by AmplitudeGlaMod2().
◆ AmplitudeGlaMod2()
Definition at line 1477 of file G4NuclNuclDiffuseElastic.hh.
1478{
1480 G4double mod2 = out.real()*out.real() + out.imag()*out.imag();
1481 return mod2;
1482}
G4complex AmplitudeGla(G4double theta)
Definition: G4NuclNuclDiffuseElastic.hh:1451
◆ AmplitudeMod2()
Definition at line 1317 of file G4NuclNuclDiffuseElastic.hh.
1318{
1320 G4double mod2 = out.real()*out.real() + out.imag()*out.imag();
1321 return mod2;
1322}
G4complex Amplitude(G4double theta)
Definition: G4NuclNuclDiffuseElastic.hh:1304
◆ AmplitudeNear()
Definition at line 1265 of file G4NuclNuclDiffuseElastic.hh.
1266{
1267 G4double kappa = std::sqrt(0.5*fProfileLambda/std::sin(theta)/CLHEP::pi);
1269
1270 out *= PhaseNear(theta);
1271
1272 if( theta <= fRutherfordTheta )
1273 {
1275 // out *= GammaMore(theta) + ProfileNear(theta);
1276 out += CoulombAmplitude(theta);
1277 }
1278 else
1279 {
1281 // out *= GammaLess(theta) + ProfileNear(theta);
1282 }
1283 return out;
1284}
G4double ProfileNear(G4double theta)
Definition: G4NuclNuclDiffuseElastic.hh:1117
G4complex PhaseNear(G4double theta)
Definition: G4NuclNuclDiffuseElastic.hh:1173
G4complex GammaMore(G4double theta)
Definition: G4NuclNuclDiffuseElastic.hh:1238
G4complex GammaLess(G4double theta)
Definition: G4NuclNuclDiffuseElastic.hh:1210
Referenced by Amplitude().
◆ AmplitudeSim()
Definition at line 1328 of file G4NuclNuclDiffuseElastic.hh.
1329{
1330 G4double sinThetaR = 2.*fHalfRutThetaTg/(1. + fHalfRutThetaTg2);
1331 G4double dTheta = 0.5*(theta - fRutherfordTheta);
1332 G4double sindTheta = std::sin(dTheta);
1333 G4double persqrt2 = std::sqrt(0.5);
1334
1336 order *= std::sqrt(0.5*fProfileLambda/sinThetaR)*2.*sindTheta;
1337 // order *= std::sqrt(0.5*fProfileLambda/sinThetaR)*2.*dTheta;
1338
1339 G4complex out;
1340
1341 if ( theta <= fRutherfordTheta )
1342 {
1344 }
1345 else
1346 {
1348 }
1349
1350 out *= CoulombAmplitude(theta);
1351
1352 return out;
1353}
G4complex GetErfcInt(G4complex z)
Definition: G4NuclNuclDiffuseElastic.hh:808
Referenced by AmplitudeSimMod2().
◆ AmplitudeSimMod2()
Definition at line 1440 of file G4NuclNuclDiffuseElastic.hh.
1441{
1443 G4double mod2 = out.real()*out.real() + out.imag()*out.imag();
1444 return mod2;
1445}
G4complex AmplitudeSim(G4double theta)
Definition: G4NuclNuclDiffuseElastic.hh:1328
◆ BesselJone()
Definition at line 458 of file G4NuclNuclDiffuseElastic.hh.
459{
460 G4double modvalue, value2, fact1, fact2, arg, shift, bessel;
461
462 modvalue = fabs(value);
463
464 if ( modvalue < 8.0 )
465 {
466 value2 = value*value;
467
468 fact1 = value*(72362614232.0 + value2*(-7895059235.0
469 + value2*( 242396853.1
470 + value2*(-2972611.439
471 + value2*( 15704.48260
472 + value2*(-30.16036606 ) ) ) ) ) );
473
474 fact2 = 144725228442.0 + value2*(2300535178.0
475 + value2*(18583304.74
476 + value2*(99447.43394
477 + value2*(376.9991397
478 + value2*1.0 ) ) ) );
479 bessel = fact1/fact2;
480 }
481 else
482 {
483 arg = 8.0/modvalue;
484
485 value2 = arg*arg;
486
487 shift = modvalue - 2.356194491;
488
489 fact1 = 1.0 + value2*( 0.183105e-2
490 + value2*(-0.3516396496e-4
491 + value2*(0.2457520174e-5
492 + value2*(-0.240337019e-6 ) ) ) );
493
494 fact2 = 0.04687499995 + value2*(-0.2002690873e-3
495 + value2*( 0.8449199096e-5
496 + value2*(-0.88228987e-6
497 + value2*0.105787412e-6 ) ) );
498
499 bessel = sqrt( 0.636619772/modvalue)*(cos(shift)*fact1 - arg*sin(shift)*fact2);
500
501 if (value < 0.0) bessel = -bessel;
502 }
503 return bessel;
504}
Referenced by BesselOneByArg(), GetDiffElasticProb(), GetDiffElasticSumProb(), and GetDiffElasticSumProbA().
◆ BesselJzero()
Definition at line 406 of file G4NuclNuclDiffuseElastic.hh.
407{
408 G4double modvalue, value2, fact1, fact2, arg, shift, bessel;
409
410 modvalue = fabs(value);
411
412 if ( value < 8.0 && value > -8.0 )
413 {
414 value2 = value*value;
415
416 fact1 = 57568490574.0 + value2*(-13362590354.0
417 + value2*( 651619640.7
418 + value2*(-11214424.18
419 + value2*( 77392.33017
420 + value2*(-184.9052456 ) ) ) ) );
421
422 fact2 = 57568490411.0 + value2*( 1029532985.0
423 + value2*( 9494680.718
424 + value2*(59272.64853
425 + value2*(267.8532712
426 + value2*1.0 ) ) ) );
427
428 bessel = fact1/fact2;
429 }
430 else
431 {
432 arg = 8.0/modvalue;
433
434 value2 = arg*arg;
435
436 shift = modvalue-0.785398164;
437
438 fact1 = 1.0 + value2*(-0.1098628627e-2
439 + value2*(0.2734510407e-4
440 + value2*(-0.2073370639e-5
441 + value2*0.2093887211e-6 ) ) );
442
443 fact2 = -0.1562499995e-1 + value2*(0.1430488765e-3
444 + value2*(-0.6911147651e-5
445 + value2*(0.7621095161e-6
446 - value2*0.934945152e-7 ) ) );
447
448 bessel = sqrt(0.636619772/modvalue)*(cos(shift)*fact1 - arg*sin(shift)*fact2 );
449 }
450 return bessel;
451}
Referenced by GetDiffElasticProb(), GetDiffElasticSumProb(), and GetDiffElasticSumProbA().
◆ BesselOneByArg()
Definition at line 533 of file G4NuclNuclDiffuseElastic.hh.
534{
535 G4double x2, result;
536
537 if( std::fabs(x) < 0.01 )
538 {
539 x *= 0.5;
540 x2 = x*x;
541 result = 2. - x2 + x2*x2/6.;
542 }
543 else
544 {
545 result = BesselJone(x)/x;
546 }
547 return result;
548}
G4double BesselJone(G4double z)
Definition: G4NuclNuclDiffuseElastic.hh:458
Referenced by GetDiffElasticProb(), GetDiffElasticSumProb(), and GetDiffElasticSumProbA().
◆ BuildAngleTable()
| void G4NuclNuclDiffuseElastic::BuildAngleTable | ( | ) |
Definition at line 958 of file G4NuclNuclDiffuseElastic.cc.
959{
960 G4int i, j;
962 G4double alpha1, alpha2, alphaMax, alphaCoulomb, delta = 0., sum = 0.;
963
964 // G4cout<<"particle z = "<<z<<"; particle m1 = "<<m1/GeV<<" GeV"<<G4endl;
965
966 G4Integrator<G4NuclNuclDiffuseElastic,G4double(G4NuclNuclDiffuseElastic::*)(G4double)> integral;
967
969
970 for( i = 0; i < fEnergyBin; i++)
971 {
972 kinE = fEnergyVector->GetLowEdgeEnergy(i);
973
974 // G4cout<<G4endl;
975 // G4cout<<"kinE = "<<kinE/MeV<<" MeV"<<G4endl;
976
977 partMom = std::sqrt( kinE*(kinE + 2*m1) );
978
979 InitDynParameters(fParticle, partMom);
980
981 alphaMax = fRutherfordTheta*fCofAlphaMax;
982
984
985
986 alphaCoulomb = fRutherfordTheta*fCofAlphaCoulomb;
987
988 // G4cout<<"alphaCoulomb = "<<alphaCoulomb/degree<<"; alphaMax = "<<alphaMax/degree<<G4endl;
989
991
992 // G4PhysicsLogVector* angleBins = new G4PhysicsLogVector( 0.001*alphaMax, alphaMax, fAngleBin );
993
994 G4double delth = (alphaMax-alphaCoulomb)/fAngleBin;
995
996 sum = 0.;
997
998 // fAddCoulomb = false;
999 fAddCoulomb = true;
1000
1001 // for(j = 1; j < fAngleBin; j++)
1002 for(j = fAngleBin-1; j >= 1; j--)
1003 {
1004 // alpha1 = angleBins->GetLowEdgeEnergy(j-1);
1005 // alpha2 = angleBins->GetLowEdgeEnergy(j);
1006
1007 // alpha1 = alphaMax*(j-1)/fAngleBin;
1008 // alpha2 = alphaMax*( j )/fAngleBin;
1009
1010 alpha1 = alphaCoulomb + delth*(j-1);
1011 // if(alpha1 < kRlim2) alpha1 = kRlim2;
1012 alpha2 = alpha1 + delth;
1013
1014 delta = integral.Legendre10(this, &G4NuclNuclDiffuseElastic::GetFresnelIntegrandXsc, alpha1, alpha2);
1015 // delta = integral.Legendre96(this, &G4NuclNuclDiffuseElastic::GetIntegrandFunction, alpha1, alpha2);
1016
1017 sum += delta;
1018
1020 // G4cout<<"j-1 = "<<j-1<<"; alpha2 = "<<alpha2/degree<<"; sum = "<<sum<<G4endl;
1021 }
1022 fAngleTable->insertAt(i,angleVector);
1023
1024 // delete[] angleVector;
1025 // delete[] angleBins;
1026 }
1027 return;
1028}
Definition: G4Integrator.hh:50
G4NuclNuclDiffuseElastic()
Definition: G4NuclNuclDiffuseElastic.cc:67
void InitDynParameters(const G4ParticleDefinition *theParticle, G4double partMom)
Definition: G4NuclNuclDiffuseElastic.hh:1576
G4double GetFresnelIntegrandXsc(G4double alpha)
Definition: G4NuclNuclDiffuseElastic.hh:1427
Definition: G4PhysicsFreeVector.hh:63
void PutValue(size_t binNumber, G4double binValue, G4double dataValue)
Definition: G4PhysicsFreeVector.cc:94
Definition: G4PhysicsTable.hh:67
virtual G4double GetLowEdgeEnergy(size_t binNumber) const
Definition: G4PhysicsVector.cc:158
Referenced by Initialise(), and InitialiseOnFly().
◆ CalcMandelstamS()
|
inline |
Definition at line 1781 of file G4NuclNuclDiffuseElastic.hh.
1784{
1785 G4double Elab = std::sqrt ( mp * mp + Plab * Plab );
1786 G4double sMand = mp*mp + mt*mt + 2*Elab*mt ;
1787
1788 return sMand;
1789}
Referenced by GetHadronNucleonXscNS().
◆ CalculateAm()
Definition at line 579 of file G4NuclNuclDiffuseElastic.hh.
580{
581 G4double k = momentum/CLHEP::hbarc;
583 G4double zn = 1.77*k*std::pow(Z,-1./3.)*CLHEP::Bohr_radius;
584 G4double zn2 = zn*zn;
585 fAm = ch/zn2;
586
587 return fAm;
588}
Referenced by GetCoulombElasticXsc(), GetCoulombIntegralXsc(), GetCoulombTotalXsc(), GetDiffuseElasticSumXsc(), InitDynParameters(), InitParameters(), InitParametersGla(), and TestAngleTable().
◆ CalculateCoulombPhase()
Definition at line 1090 of file G4NuclNuclDiffuseElastic.hh.
1091{
1093 // G4complex gammalog = GammaLogarithm(z);
1095 return gammalog.imag();
1096}
G4complex GammaLogB2n(G4complex xx)
Definition: G4NuclNuclDiffuseElastic.hh:746
Referenced by AmplitudeGla().
◆ CalculateCoulombPhaseZero()
|
inline |
Definition at line 1077 of file G4NuclNuclDiffuseElastic.hh.
1078{
1080 // G4complex gammalog = GammaLogarithm(z);
1082 fCoulombPhase0 = gammalog.imag();
1083}
Referenced by InitDynParameters(), InitParameters(), and InitParametersGla().
◆ CalculateNuclearRad()
Definition at line 594 of file G4NuclNuclDiffuseElastic.hh.
595{
596 G4double r0 = 1.*CLHEP::fermi, radius;
597 // r0 *= 1.12;
598 // r0 *= 1.44;
599 r0 *= fNuclearRadiusCof;
600
601 /*
602 if( A < 50. )
603 {
604 if( A > 10. ) r0 = 1.16*( 1 - std::pow(A, -2./3.) )*CLHEP::fermi; // 1.08*fermi;
605 else r0 = 1.1*CLHEP::fermi;
606
607 radius = r0*std::pow(A, 1./3.);
608 }
609 else
610 {
611 r0 = 1.7*CLHEP::fermi; // 1.7*fermi;
612
613 radius = r0*std::pow(A, 0.27); // 0.27);
614 }
615 */
616 radius = r0*std::pow(A, 1./3.);
617
618 return radius;
619}
Referenced by GetDiffuseElasticSumXsc(), GetDiffuseElasticXsc(), Initialise(), InitialiseOnFly(), InitParameters(), InitParametersGla(), IntegralElasticProb(), SampleThetaCMS(), and TestAngleTable().
◆ CalculateParticleBeta()
|
inline |
Definition at line 554 of file G4NuclNuclDiffuseElastic.hh.
Referenced by GetCoulombElasticXsc(), GetCoulombIntegralXsc(), GetCoulombTotalXsc(), and GetDiffuseElasticSumXsc().
◆ CalculateRutherfordAnglePar()
|
inline |
Definition at line 1104 of file G4NuclNuclDiffuseElastic.hh.
1105{
1106 fHalfRutThetaTg = fZommerfeld/fProfileLambda; // (fWaveVector*fNuclearRadius);
1107 fRutherfordTheta = 2.*std::atan(fHalfRutThetaTg);
1108 fHalfRutThetaTg2 = fHalfRutThetaTg*fHalfRutThetaTg;
1109 // G4cout<<"fRutherfordTheta = "<<fRutherfordTheta/degree<<" degree"<<G4endl;
1110
1111}
Referenced by InitDynParameters(), and InitParameters().
◆ CalculateZommerfeld()
|
inline |
Definition at line 568 of file G4NuclNuclDiffuseElastic.hh.
569{
570 fZommerfeld = CLHEP::fine_structure_const*Z1*Z2/beta;
571
572 return fZommerfeld;
573}
Referenced by GetCoulombElasticXsc(), GetCoulombIntegralXsc(), GetCoulombTotalXsc(), GetDiffuseElasticSumXsc(), InitDynParameters(), InitParameters(), InitParametersGla(), and TestAngleTable().
◆ CoulombAmplitude()
Definition at line 1043 of file G4NuclNuclDiffuseElastic.hh.
1044{
1045 G4complex ca;
1046
1047 G4double sinHalfTheta = std::sin(0.5*theta);
1048 G4double sinHalfTheta2 = sinHalfTheta*sinHalfTheta;
1049 sinHalfTheta2 += fAm;
1050
1051 G4double order = 2.*fCoulombPhase0 - fZommerfeld*std::log(sinHalfTheta2);
1053 ca = std::exp(z);
1054
1055 ca *= -fZommerfeld/(2.*fWaveVector*sinHalfTheta2);
1056
1057 return ca;
1058}
Referenced by AmplitudeGG(), AmplitudeGla(), AmplitudeNear(), AmplitudeSim(), and CoulombAmplitudeMod2().
◆ CoulombAmplitudeMod2()
Definition at line 1064 of file G4NuclNuclDiffuseElastic.hh.
◆ DampFactor()
Definition at line 510 of file G4NuclNuclDiffuseElastic.hh.
511{
512 G4double df;
514
515 // x *= pi;
516
517 if( std::fabs(x) < 0.01 )
518 {
519 df = 1./(1. + x/f2 + x*x/f3 + x*x*x/f4);
520 }
521 else
522 {
523 df = x/std::sinh(x);
524 }
525 return df;
526}
Referenced by GetDiffElasticProb(), GetDiffElasticSumProb(), and GetDiffElasticSumProbA().
◆ GammaLess()
Definition at line 1210 of file G4NuclNuclDiffuseElastic.hh.
1211{
1212 G4double sinThetaR = 2.*fHalfRutThetaTg/(1. + fHalfRutThetaTg2);
1213 G4double cosHalfThetaR2 = 1./(1. + fHalfRutThetaTg2);
1214
1215 G4double u = std::sqrt(0.5*fProfileLambda/sinThetaR);
1216 G4double kappa = u/std::sqrt(CLHEP::pi);
1217 G4double dTheta = theta - fRutherfordTheta;
1218 u *= dTheta;
1219 G4double u2 = u*u;
1220 G4double u2m2p3 = u2*2./3.;
1221
1224 order /= std::sqrt(2.);
1225
1227 G4complex a0 = 0.5*(1. + 4.*(1.+im*u2)*cosHalfThetaR2/3.)/sinThetaR;
1228 G4complex a1 = 0.5*(1. + 2.*(1.+im*u2m2p3)*cosHalfThetaR2)/sinThetaR;
1229 G4complex out = gamma*(1. - a1*dTheta) - a0;
1230
1231 return out;
1232}
Referenced by AmplitudeNear().
◆ GammaLogarithm()
Definition at line 722 of file G4NuclNuclDiffuseElastic.hh.
723{
725 24.01409824083091, -1.231739572450155,
726 0.1208650973866179e-2, -0.5395239384953e-5 } ;
728 G4complex z = zz - 1.0;
729 G4complex tmp = z + 5.5;
730 tmp -= (z + 0.5) * std::log(tmp);
732
733 for ( j = 0; j <= 5; j++ )
734 {
735 z += 1.0;
736 ser += cof[j]/z;
737 }
738 return -tmp + std::log(2.5066282746310005*ser);
739}
◆ GammaLogB2n()
Definition at line 746 of file G4NuclNuclDiffuseElastic.hh.
747{
748 G4complex z1 = 12.*z;
749 G4complex z2 = z*z;
750 G4complex z3 = z2*z;
751 G4complex z5 = z2*z3;
752 G4complex z7 = z2*z5;
753
754 z3 *= 360.;
755 z5 *= 1260.;
756 z7 *= 1680.;
757
758 G4complex result = (z-0.5)*std::log(z)-z+0.5*std::log(CLHEP::twopi);
759 result += 1./z1 - 1./z3 +1./z5 -1./z7;
760 return result;
761}
Referenced by CalculateCoulombPhase(), and CalculateCoulombPhaseZero().
◆ GammaMore()
Definition at line 1238 of file G4NuclNuclDiffuseElastic.hh.
1239{
1240 G4double sinThetaR = 2.*fHalfRutThetaTg/(1. + fHalfRutThetaTg2);
1241 G4double cosHalfThetaR2 = 1./(1. + fHalfRutThetaTg2);
1242
1243 G4double u = std::sqrt(0.5*fProfileLambda/sinThetaR);
1244 G4double kappa = u/std::sqrt(CLHEP::pi);
1245 G4double dTheta = theta - fRutherfordTheta;
1246 u *= dTheta;
1247 G4double u2 = u*u;
1248 G4double u2m2p3 = u2*2./3.;
1249
1252 order /= std::sqrt(2.);
1254 G4complex a0 = 0.5*(1. + 4.*(1.+im*u2)*cosHalfThetaR2/3.)/sinThetaR;
1255 G4complex a1 = 0.5*(1. + 2.*(1.+im*u2m2p3)*cosHalfThetaR2)/sinThetaR;
1256 G4complex out = -gamma*(1. - a1*dTheta) - a0;
1257
1258 return out;
1259}
Referenced by AmplitudeNear().
◆ GetCint()
Definition at line 1012 of file G4NuclNuclDiffuseElastic.hh.
1013{
1014 G4double out;
1015
1016 G4Integrator<G4NuclNuclDiffuseElastic,G4double(G4NuclNuclDiffuseElastic::*)(G4double)> integral;
1017
1019
1020 return out;
1021}
G4double GetCosHaPit2(G4double t)
Definition: G4NuclNuclDiffuseElastic.hh:198
Referenced by GetRatioGen(), and GetRatioSim().
◆ GetCofAlphaCoulomb()
|
inline |
Definition at line 301 of file G4NuclNuclDiffuseElastic.hh.
301{return fCofAlphaCoulomb;};
◆ GetCofAlphaMax()
|
inline |
Definition at line 300 of file G4NuclNuclDiffuseElastic.hh.
300{return fCofAlphaMax;};
◆ GetCosHaPit2()
Definition at line 198 of file G4NuclNuclDiffuseElastic.hh.
198{return std::cos(CLHEP::halfpi*t*t);};
Referenced by GetCint().
◆ GetCoulombElasticXsc()
|
inline |
Definition at line 625 of file G4NuclNuclDiffuseElastic.hh.
629{
630 G4double sinHalfTheta = std::sin(0.5*theta);
631 G4double sinHalfTheta2 = sinHalfTheta*sinHalfTheta;
636 G4double k = momentum/CLHEP::hbarc;
638 G4double ch2 = ch*ch;
639 G4double xsc = ch2/(sinHalfTheta2+am)/(sinHalfTheta2+am);
640
641 return xsc;
642}
G4double CalculateZommerfeld(G4double beta, G4double Z1, G4double Z2)
Definition: G4NuclNuclDiffuseElastic.hh:568
G4double CalculateParticleBeta(const G4ParticleDefinition *particle, G4double momentum)
Definition: G4NuclNuclDiffuseElastic.hh:554
G4double CalculateAm(G4double momentum, G4double n, G4double Z)
Definition: G4NuclNuclDiffuseElastic.hh:579
G4double GetPDGCharge() const
Definition: G4ParticleDefinition.hh:119
Referenced by GetInvCoulombElasticXsc().
◆ GetCoulombIntegralXsc()
|
inline |
Definition at line 690 of file G4NuclNuclDiffuseElastic.hh.
693{
694 G4double c1 = std::cos(theta1);
696 G4double c2 = std::cos(theta2);
699 // G4cout<<"beta = "<<beta<<G4endl;
702 // G4cout<<"fZomerfeld = "<<n<<G4endl;
704 // G4cout<<"cof Am = "<<am<<G4endl;
705 G4double k = momentum/CLHEP::hbarc;
706 // G4cout<<"k = "<<k*CLHEP::fermi<<" 1/fermi"<<G4endl;
707 // G4cout<<"k*Bohr_radius = "<<k*CLHEP::Bohr_radius<<G4endl;
709 G4double ch2 = ch*ch;
710 am *= 2.;
711 G4double xsc = ch2*CLHEP::twopi*(c1-c2);
712 xsc /= (1 - c1 + am)*(1 - c2 + am);
713
714 return xsc;
715}
◆ GetCoulombTotalXsc()
|
inline |
Definition at line 665 of file G4NuclNuclDiffuseElastic.hh.
◆ GetDiffElasticProb()
Definition at line 389 of file G4NuclNuclDiffuseElastic.cc.
394{
395 G4double sigma, bzero, bzero2, bonebyarg, bonebyarg2, damp, damp2;
396 G4double delta, diffuse, gamma;
397 G4double e1, e2, bone, bone2;
398
399 // G4double wavek = momentum/hbarc; // wave vector
400 // G4double r0 = 1.08*fermi;
401 // G4double rad = r0*std::pow(A, 1./3.);
402
404 G4double kr2 = kr*kr;
405 G4double krt = kr*theta;
406
407 bzero = BesselJzero(krt);
408 bzero2 = bzero*bzero;
409 bone = BesselJone(krt);
410 bone2 = bone*bone;
411 bonebyarg = BesselOneByArg(krt);
412 bonebyarg2 = bonebyarg*bonebyarg;
413
414 if (fParticle == theProton)
415 {
416 diffuse = 0.63*fermi;
417 gamma = 0.3*fermi;
418 delta = 0.1*fermi*fermi;
419 e1 = 0.3*fermi;
420 e2 = 0.35*fermi;
421 }
422 else // as proton, if were not defined
423 {
424 diffuse = 0.63*fermi;
425 gamma = 0.3*fermi;
426 delta = 0.1*fermi*fermi;
427 e1 = 0.3*fermi;
428 e2 = 0.35*fermi;
429 }
431
432 // G4double kgamma = fWaveVector*gamma; // wavek*delta;
433
435 G4double kgamma2 = kgamma*kgamma;
436
437 // G4double dk2t = delta*fWaveVector*fWaveVector*theta; // delta*wavek*wavek*theta;
438 // G4double dk2t2 = dk2t*dk2t;
439 // G4double pikdt = pi*fWaveVector*diffuse*theta;// pi*wavek*diffuse*theta;
440
442
443 damp = DampFactor(pikdt);
444 damp2 = damp*damp;
445
446 G4double mode2k2 = (e1*e1+e2*e2)*fWaveVector*fWaveVector;
447 G4double e2dk3t = -2.*e2*delta*fWaveVector*fWaveVector*fWaveVector*theta;
448
449
450 sigma = kgamma2;
451 // sigma += dk2t2;
452 sigma *= bzero2;
453 sigma += mode2k2*bone2 + e2dk3t*bzero*bone;
454 sigma += kr2*bonebyarg2;
455 sigma *= damp2; // *rad*rad;
456
457 return sigma;
458}
G4double BesselJzero(G4double z)
Definition: G4NuclNuclDiffuseElastic.hh:406
G4double BesselOneByArg(G4double z)
Definition: G4NuclNuclDiffuseElastic.hh:533
G4double DampFactor(G4double z)
Definition: G4NuclNuclDiffuseElastic.hh:510
Referenced by GetDiffuseElasticXsc().
◆ GetDiffElasticSumProb()
Definition at line 466 of file G4NuclNuclDiffuseElastic.cc.
471{
472 G4double sigma, bzero, bzero2, bonebyarg, bonebyarg2, damp, damp2;
473 G4double delta, diffuse, gamma;
474 G4double e1, e2, bone, bone2;
475
476 // G4double wavek = momentum/hbarc; // wave vector
477 // G4double r0 = 1.08*fermi;
478 // G4double rad = r0*std::pow(A, 1./3.);
479
481 G4double kr2 = kr*kr;
482 G4double krt = kr*theta;
483
484 bzero = BesselJzero(krt);
485 bzero2 = bzero*bzero;
486 bone = BesselJone(krt);
487 bone2 = bone*bone;
488 bonebyarg = BesselOneByArg(krt);
489 bonebyarg2 = bonebyarg*bonebyarg;
490
491 if (fParticle == theProton)
492 {
493 diffuse = 0.63*fermi;
494 // diffuse = 0.6*fermi;
495 gamma = 0.3*fermi;
496 delta = 0.1*fermi*fermi;
497 e1 = 0.3*fermi;
498 e2 = 0.35*fermi;
499 }
500 else // as proton, if were not defined
501 {
502 diffuse = 0.63*fermi;
503 gamma = 0.3*fermi;
504 delta = 0.1*fermi*fermi;
505 e1 = 0.3*fermi;
506 e2 = 0.35*fermi;
507 }
509 // G4double kgamma = fWaveVector*gamma; // wavek*delta;
511
512 // G4cout<<"kgamma = "<<kgamma<<G4endl;
513
514 if(fAddCoulomb) // add Coulomb correction
515 {
516 G4double sinHalfTheta = std::sin(0.5*theta);
517 G4double sinHalfTheta2 = sinHalfTheta*sinHalfTheta;
518
519 kgamma += 0.5*fZommerfeld/kr/(sinHalfTheta2+fAm); // correction at J0()
520 // kgamma += 0.65*fZommerfeld/kr/(sinHalfTheta2+fAm); // correction at J0()
521 }
522
523 G4double kgamma2 = kgamma*kgamma;
524
525 // G4double dk2t = delta*fWaveVector*fWaveVector*theta; // delta*wavek*wavek*theta;
526 // G4cout<<"dk2t = "<<dk2t<<G4endl;
527 // G4double dk2t2 = dk2t*dk2t;
528 // G4double pikdt = pi*fWaveVector*diffuse*theta;// pi*wavek*diffuse*theta;
529
531
532 // G4cout<<"pikdt = "<<pikdt<<G4endl;
533
534 damp = DampFactor(pikdt);
535 damp2 = damp*damp;
536
537 G4double mode2k2 = (e1*e1+e2*e2)*fWaveVector*fWaveVector;
538 G4double e2dk3t = -2.*e2*delta*fWaveVector*fWaveVector*fWaveVector*theta;
539
540 sigma = kgamma2;
541 // sigma += dk2t2;
542 sigma *= bzero2;
543 sigma += mode2k2*bone2;
544 sigma += e2dk3t*bzero*bone;
545
546 // sigma += kr2*(1 + 8.*fZommerfeld*fZommerfeld/kr2)*bonebyarg2; // correction at J1()/()
547 sigma += kr2*bonebyarg2; // correction at J1()/()
548
549 sigma *= damp2; // *rad*rad;
550
551 return sigma;
552}
Referenced by GetDiffuseElasticSumXsc().
◆ GetDiffElasticSumProbA()
Definition at line 561 of file G4NuclNuclDiffuseElastic.cc.
562{
563 G4double theta;
564
565 theta = std::sqrt(alpha);
566
567 // theta = std::acos( 1 - alpha/2. );
568
569 G4double sigma, bzero, bzero2, bonebyarg, bonebyarg2, damp, damp2;
570 G4double delta, diffuse, gamma;
571 G4double e1, e2, bone, bone2;
572
573 // G4double wavek = momentum/hbarc; // wave vector
574 // G4double r0 = 1.08*fermi;
575 // G4double rad = r0*std::pow(A, 1./3.);
576
578 G4double kr2 = kr*kr;
579 G4double krt = kr*theta;
580
581 bzero = BesselJzero(krt);
582 bzero2 = bzero*bzero;
583 bone = BesselJone(krt);
584 bone2 = bone*bone;
585 bonebyarg = BesselOneByArg(krt);
586 bonebyarg2 = bonebyarg*bonebyarg;
587
588 if (fParticle == theProton)
589 {
590 diffuse = 0.63*fermi;
591 // diffuse = 0.6*fermi;
592 gamma = 0.3*fermi;
593 delta = 0.1*fermi*fermi;
594 e1 = 0.3*fermi;
595 e2 = 0.35*fermi;
596 }
597 else // as proton, if were not defined
598 {
599 diffuse = 0.63*fermi;
600 gamma = 0.3*fermi;
601 delta = 0.1*fermi*fermi;
602 e1 = 0.3*fermi;
603 e2 = 0.35*fermi;
604 }
606 // G4double kgamma = fWaveVector*gamma; // wavek*delta;
608
609 // G4cout<<"kgamma = "<<kgamma<<G4endl;
610
611 if(fAddCoulomb) // add Coulomb correction
612 {
614 G4double sinHalfTheta2 = sinHalfTheta*sinHalfTheta;
615
616 kgamma += 0.5*fZommerfeld/kr/(sinHalfTheta2+fAm); // correction at J0()
617 // kgamma += 0.65*fZommerfeld/kr/(sinHalfTheta2+fAm); // correction at J0()
618 }
619
620 G4double kgamma2 = kgamma*kgamma;
621
622 // G4double dk2t = delta*fWaveVector*fWaveVector*theta; // delta*wavek*wavek*theta;
623 // G4cout<<"dk2t = "<<dk2t<<G4endl;
624 // G4double dk2t2 = dk2t*dk2t;
625 // G4double pikdt = pi*fWaveVector*diffuse*theta;// pi*wavek*diffuse*theta;
626
628
629 // G4cout<<"pikdt = "<<pikdt<<G4endl;
630
631 damp = DampFactor(pikdt);
632 damp2 = damp*damp;
633
634 G4double mode2k2 = (e1*e1+e2*e2)*fWaveVector*fWaveVector;
635 G4double e2dk3t = -2.*e2*delta*fWaveVector*fWaveVector*fWaveVector*theta;
636
637 sigma = kgamma2;
638 // sigma += dk2t2;
639 sigma *= bzero2;
640 sigma += mode2k2*bone2;
641 sigma += e2dk3t*bzero*bone;
642
643 // sigma += kr2*(1 + 8.*fZommerfeld*fZommerfeld/kr2)*bonebyarg2; // correction at J1()/()
644 sigma += kr2*bonebyarg2; // correction at J1()/()
645
646 sigma *= damp2; // *rad*rad;
647
648 return sigma;
649}
Referenced by GetIntegrandFunction().
◆ GetDiffuseElasticSumXsc()
| G4double G4NuclNuclDiffuseElastic::GetDiffuseElasticSumXsc | ( | const G4ParticleDefinition * | particle, |
| G4double | theta, | ||
| G4double | momentum, | ||
| G4double | A, | ||
| G4double | Z | ||
| ) |
Definition at line 253 of file G4NuclNuclDiffuseElastic.cc.
257{
258 fParticle = particle;
259 fWaveVector = momentum/hbarc;
260 fAtomicWeight = A;
261 fAtomicNumber = Z;
262 fNuclearRadius = CalculateNuclearRad(A);
263 fAddCoulomb = false;
264
266
267 G4double kRt = fWaveVector*fNuclearRadius*theta;
268 G4double kRtC = 1.9;
269
270 if( z && (kRt > kRtC) )
271 {
272 fAddCoulomb = true;
273 fBeta = CalculateParticleBeta( particle, momentum);
274 fZommerfeld = CalculateZommerfeld( fBeta, z, fAtomicNumber);
275 fAm = CalculateAm( momentum, fZommerfeld, fAtomicNumber);
276 }
278
279 return sigma;
280}
G4double GetDiffElasticSumProb(G4double theta)
Definition: G4NuclNuclDiffuseElastic.cc:466
G4double CalculateNuclearRad(G4double A)
Definition: G4NuclNuclDiffuseElastic.hh:594
Referenced by GetInvElasticSumXsc().
◆ GetDiffuseElasticXsc()
| G4double G4NuclNuclDiffuseElastic::GetDiffuseElasticXsc | ( | const G4ParticleDefinition * | particle, |
| G4double | theta, | ||
| G4double | momentum, | ||
| G4double | A | ||
| ) |
Definition at line 182 of file G4NuclNuclDiffuseElastic.cc.
186{
187 fParticle = particle;
188 fWaveVector = momentum/hbarc;
189 fAtomicWeight = A;
190 fAddCoulomb = false;
191 fNuclearRadius = CalculateNuclearRad(A);
192
194
195 return sigma;
196}
G4double GetDiffElasticProb(G4double theta)
Definition: G4NuclNuclDiffuseElastic.cc:389
Referenced by GetInvElasticXsc().
◆ GetErf()
Definition at line 767 of file G4NuclNuclDiffuseElastic.hh.
768{
769 G4double t, z, tmp, result;
770
771 z = std::fabs(x);
772 t = 1.0/(1.0+0.5*z);
773
774 tmp = t*exp(-z*z-1.26551223+t*(1.00002368+t*(0.37409196+t*(0.09678418+
775 t*(-0.18628806+t*(0.27886807+t*(-1.13520398+t*(1.48851587+
776 t*(-0.82215223+t*0.17087277)))))))));
777
778 if( x >= 0.) result = 1. - tmp;
779 else result = 1. + tmp;
780
781 return result;
782}
Referenced by GetErfComp(), and GetErfInt().
◆ GetErfcComp()
Definition at line 788 of file G4NuclNuclDiffuseElastic.hh.
789{
791 return erfcz;
792}
G4complex GetErfComp(G4complex z, G4int nMax)
Definition: G4NuclNuclDiffuseElastic.hh:878
◆ GetErfcInt()
Definition at line 808 of file G4NuclNuclDiffuseElastic.hh.
809{
811 return erfcz;
812}
G4complex GetErfInt(G4complex z)
Definition: G4NuclNuclDiffuseElastic.hh:987
Referenced by AmplitudeSim(), GammaLess(), and GammaMore().
◆ GetErfComp()
Definition at line 878 of file G4NuclNuclDiffuseElastic.hh.
879{
881 G4double n2, cofn, shny, chny, fn, gn;
882
883 G4double x = z.real();
884 G4double y = z.imag();
885
886 G4double outRe = 0., outIm = 0.;
887
888 G4double twox = 2.*x;
889 G4double twoxy = twox*y;
890 G4double twox2 = twox*twox;
891
892 G4double cof1 = std::exp(-x*x)/CLHEP::pi;
893
894 G4double cos2xy = std::cos(twoxy);
895 G4double sin2xy = std::sin(twoxy);
896
897 G4double twoxcos2xy = twox*cos2xy;
898 G4double twoxsin2xy = twox*sin2xy;
899
901 {
903
904 cofn = std::exp(-0.5*n2)/(n2+twox2); // /(n2+0.5*twox2);
905
906 chny = std::cosh(n*y);
907 shny = std::sinh(n*y);
908
909 fn = twox - twoxcos2xy*chny + n*sin2xy*shny;
910 gn = twoxsin2xy*chny + n*cos2xy*shny;
911
912 fn *= cofn;
913 gn *= cofn;
914
915 outRe += fn;
916 outIm += gn;
917 }
918 outRe *= 2*cof1;
919 outIm *= 2*cof1;
920
921 if(std::abs(x) < 0.0001)
922 {
923 outRe += GetErf(x);
924 outIm += cof1*y;
925 }
926 else
927 {
928 outRe += GetErf(x) + cof1*(1-cos2xy)/twox;
929 outIm += cof1*sin2xy/twox;
930 }
932}
G4double GetErf(G4double x)
Definition: G4NuclNuclDiffuseElastic.hh:767
Referenced by GetErfcComp(), and TestErfcComp().
◆ GetErfcSer()
Definition at line 798 of file G4NuclNuclDiffuseElastic.hh.
799{
801 return erfcz;
802}
G4complex GetErfSer(G4complex z, G4int nMax)
Definition: G4NuclNuclDiffuseElastic.hh:938
◆ GetErfInt()
Definition at line 987 of file G4NuclNuclDiffuseElastic.hh.
988{
989 G4double outRe, outIm;
990
991 G4double x = z.real();
992 G4double y = z.imag();
993 fReZ = x;
994
995 G4Integrator<G4NuclNuclDiffuseElastic,G4double(G4NuclNuclDiffuseElastic::*)(G4double)> integral;
996
999
1000 outRe *= 2./sqrt(CLHEP::pi);
1001 outIm *= 2./sqrt(CLHEP::pi);
1002
1003 outRe += GetErf(x);
1004
1006}
G4double GetExpSin(G4double x)
Definition: G4NuclNuclDiffuseElastic.hh:972
G4double GetExpCos(G4double x)
Definition: G4NuclNuclDiffuseElastic.hh:961
Referenced by GetErfcInt(), and TestErfcInt().
◆ GetErfSer()
Definition at line 938 of file G4NuclNuclDiffuseElastic.hh.
939{
941 G4double a =1., b = 1., tmp;
942 G4complex sum = z, d = z;
943
945 {
946 a *= 2.;
947 b *= 2.*n +1.;
948 d *= z*z;
949
950 tmp = a/b;
951
952 sum += tmp*d;
953 }
954 sum *= 2.*std::exp(-z*z)/std::sqrt(CLHEP::pi);
955
956 return sum;
957}
Referenced by GetErfcSer(), and TestErfcSer().
◆ GetExpCos()
Definition at line 961 of file G4NuclNuclDiffuseElastic.hh.
962{
963 G4double result;
964
965 result = std::exp(x*x-fReZ*fReZ);
966 result *= std::cos(2.*x*fReZ);
967 return result;
968}
Referenced by GetErfInt().
◆ GetExpSin()
Definition at line 972 of file G4NuclNuclDiffuseElastic.hh.
973{
974 G4double result;
975
976 result = std::exp(x*x-fReZ*fReZ);
977 result *= std::sin(2.*x*fReZ);
978 return result;
979}
Referenced by GetErfInt().
◆ GetFresnelDiffuseXsc()
Definition at line 1415 of file G4NuclNuclDiffuseElastic.hh.
1416{
1419 G4double xsc = ratio*ruthXsc;
1420 return xsc;
1421}
G4double GetRatioGen(G4double theta)
Definition: G4NuclNuclDiffuseElastic.hh:1379
G4double GetRutherfordXsc(G4double theta)
Definition: G4NuclNuclDiffuseElastic.hh:648
Referenced by GetFresnelIntegrandXsc().
◆ GetFresnelIntegrandXsc()
Definition at line 1427 of file G4NuclNuclDiffuseElastic.hh.
1428{
1429 G4double theta = std::sqrt(alpha);
1431 return xsc;
1432}
G4double GetFresnelDiffuseXsc(G4double theta)
Definition: G4NuclNuclDiffuseElastic.hh:1415
Referenced by BuildAngleTable().
◆ GetHadronNucleonXscNS()
|
inline |
Definition at line 1671 of file G4NuclNuclDiffuseElastic.hh.
1674{
1676 G4double hpXsc(0);
1677 G4double hnXsc(0);
1678
1679
1682
1683 G4double proj_energy = proj_mass + pTkin;
1684 G4double proj_momentum = std::sqrt(pTkin*(pTkin+2*proj_mass));
1685
1687
1688 sMand /= CLHEP::GeV*CLHEP::GeV; // in GeV for parametrisation
1689 proj_momentum /= CLHEP::GeV;
1690 proj_energy /= CLHEP::GeV;
1691 proj_mass /= CLHEP::GeV;
1692 G4double logS = std::log(sMand);
1693
1694 // General PDG fit constants
1695
1696
1697 // fEtaRatio=Re[f(0)]/Im[f(0)]
1698
1699 if( proj_momentum >= 1.2 )
1700 {
1701 fEtaRatio = 0.13*(logS - 5.8579332)*std::pow(sMand,-0.18);
1702 }
1703 else if( proj_momentum >= 0.6 )
1704 {
1705 fEtaRatio = -75.5*(std::pow(proj_momentum,0.25)-0.95)/
1706 (std::pow(3*proj_momentum,2.2)+1);
1707 }
1708 else
1709 {
1710 fEtaRatio = 15.5*proj_momentum/(27*proj_momentum*proj_momentum*proj_momentum+2);
1711 }
1713
1714 // xsc
1715
1716 if( proj_momentum >= 10. ) // high energy: pp = nn = np
1717 // if( proj_momentum >= 2.)
1718 {
1719 //Delta = 1.;
1720
1721 //if( proj_energy < 40. ) Delta = 0.916+0.0021*proj_energy;
1722
1723 if( proj_momentum >= 10.)
1724 {
1725 B0 = 7.5;
1726 A0 = 100. - B0*std::log(3.0e7);
1727
1728 xsection = A0 + B0*std::log(proj_energy) - 11
1729 + 103*std::pow(2*0.93827*proj_energy + proj_mass*proj_mass+
1730 0.93827*0.93827,-0.165); // mb
1731 }
1732 }
1733 else // low energy pp = nn != np
1734 {
1735 if(pParticle == tParticle) // pp or nn // nn to be pp
1736 {
1737 if( proj_momentum < 0.73 )
1738 {
1739 hnXsc = 23 + 50*( std::pow( std::log(0.73/proj_momentum), 3.5 ) );
1740 }
1741 else if( proj_momentum < 1.05 )
1742 {
1743 hnXsc = 23 + 40*(std::log(proj_momentum/0.73))*
1744 (std::log(proj_momentum/0.73));
1745 }
1746 else // if( proj_momentum < 10. )
1747 {
1748 hnXsc = 39.0 +
1749 75*(proj_momentum - 1.2)/(std::pow(proj_momentum,3.0) + 0.15);
1750 }
1751 xsection = hnXsc;
1752 }
1753 else // pn to be np
1754 {
1755 if( proj_momentum < 0.8 )
1756 {
1757 hpXsc = 33+30*std::pow(std::log(proj_momentum/1.3),4.0);
1758 }
1759 else if( proj_momentum < 1.4 )
1760 {
1761 hpXsc = 33+30*std::pow(std::log(proj_momentum/0.95),2.0);
1762 }
1763 else // if( proj_momentum < 10. )
1764 {
1765 hpXsc = 33.3+
1766 20.8*(std::pow(proj_momentum,2.0)-1.35)/
1767 (std::pow(proj_momentum,2.50)+0.95);
1768 }
1769 xsection = hpXsc;
1770 }
1771 }
1772 xsection *= CLHEP::millibarn; // parametrised in mb
1774 return xsection;
1775}
G4double CalcMandelstamS(const G4double mp, const G4double mt, const G4double Plab)
Definition: G4NuclNuclDiffuseElastic.hh:1781
Referenced by InitParametersGla().
◆ GetIntegrandFunction()
Definition at line 657 of file G4NuclNuclDiffuseElastic.cc.
658{
659 G4double result;
660
661 result = GetDiffElasticSumProbA(alpha);
662
663 // result *= 2*pi*std::sin(theta);
664
665 return result;
666}
G4double GetDiffElasticSumProbA(G4double alpha)
Definition: G4NuclNuclDiffuseElastic.cc:561
Referenced by IntegralElasticProb(), SampleThetaCMS(), and TestAngleTable().
◆ GetInvCoulombElasticXsc()
| G4double G4NuclNuclDiffuseElastic::GetInvCoulombElasticXsc | ( | const G4ParticleDefinition * | particle, |
| G4double | tMand, | ||
| G4double | momentum, | ||
| G4double | A, | ||
| G4double | Z | ||
| ) |
Definition at line 340 of file G4NuclNuclDiffuseElastic.cc.
344{
346 G4LorentzVector lv1(0.,0.,plab,std::sqrt(plab*plab+m1*m1));
347
350 G4ParticleDefinition * theDef = 0;
351
352 if (iZ == 1 && iA == 1) theDef = theProton;
353 else if (iZ == 1 && iA == 2) theDef = theDeuteron;
356 else if (iZ == 2 && iA == 4) theDef = theAlpha;
358
360
361 G4LorentzVector lv(0.0,0.0,0.0,tmass);
362 lv += lv1;
363
364 G4ThreeVector bst = lv.boostVector();
365 lv1.boost(-bst);
366
367 G4ThreeVector p1 = lv1.vect();
369 G4double ptot2 = ptot*ptot;
370 G4double cost = 1 - 0.5*std::fabs(tMand)/ptot2;
371
372 if( cost >= 1.0 ) cost = 1.0;
373 else if( cost <= -1.0) cost = -1.0;
374
375 G4double thetaCMS = std::acos(cost);
376
378
379 sigma *= pi/ptot2;
380
381 return sigma;
382}
Definition: ThreeVector.h:41
double mag() const
Definition: LorentzVector.h:72
G4double GetCoulombElasticXsc(const G4ParticleDefinition *particle, G4double theta, G4double momentum, G4double Z)
Definition: G4NuclNuclDiffuseElastic.hh:625
Definition: G4ParticleDefinition.hh:68
G4ParticleDefinition * FindIon(G4int atomicNumber, G4int atomicMass, G4double excitationEnergy)
Definition: G4ParticleTable.cc:399
static G4ParticleTable * GetParticleTable()
Definition: G4ParticleTable.cc:63
◆ GetInvElasticSumXsc()
| G4double G4NuclNuclDiffuseElastic::GetInvElasticSumXsc | ( | const G4ParticleDefinition * | particle, |
| G4double | tMand, | ||
| G4double | momentum, | ||
| G4double | A, | ||
| G4double | Z | ||
| ) |
Definition at line 288 of file G4NuclNuclDiffuseElastic.cc.
292{
294
295 G4LorentzVector lv1(0.,0.,plab,std::sqrt(plab*plab+m1*m1));
296
299
300 G4ParticleDefinition* theDef = 0;
301
302 if (iZ == 1 && iA == 1) theDef = theProton;
303 else if (iZ == 1 && iA == 2) theDef = theDeuteron;
306 else if (iZ == 2 && iA == 4) theDef = theAlpha;
308
310
311 G4LorentzVector lv(0.0,0.0,0.0,tmass);
312 lv += lv1;
313
314 G4ThreeVector bst = lv.boostVector();
315 lv1.boost(-bst);
316
317 G4ThreeVector p1 = lv1.vect();
319 G4double ptot2 = ptot*ptot;
320 G4double cost = 1 - 0.5*std::fabs(tMand)/ptot2;
321
322 if( cost >= 1.0 ) cost = 1.0;
323 else if( cost <= -1.0) cost = -1.0;
324
325 G4double thetaCMS = std::acos(cost);
326
328
329 sigma *= pi/ptot2;
330
331 return sigma;
332}
G4double GetDiffuseElasticSumXsc(const G4ParticleDefinition *particle, G4double theta, G4double momentum, G4double A, G4double Z)
Definition: G4NuclNuclDiffuseElastic.cc:253
◆ GetInvElasticXsc()
| G4double G4NuclNuclDiffuseElastic::GetInvElasticXsc | ( | const G4ParticleDefinition * | particle, |
| G4double | theta, | ||
| G4double | momentum, | ||
| G4double | A, | ||
| G4double | Z | ||
| ) |
Definition at line 203 of file G4NuclNuclDiffuseElastic.cc.
207{
209 G4LorentzVector lv1(0.,0.,plab,std::sqrt(plab*plab+m1*m1));
210
213 G4ParticleDefinition * theDef = 0;
214
215 if (iZ == 1 && iA == 1) theDef = theProton;
216 else if (iZ == 1 && iA == 2) theDef = theDeuteron;
219 else if (iZ == 2 && iA == 4) theDef = theAlpha;
221
223
224 G4LorentzVector lv(0.0,0.0,0.0,tmass);
225 lv += lv1;
226
227 G4ThreeVector bst = lv.boostVector();
228 lv1.boost(-bst);
229
230 G4ThreeVector p1 = lv1.vect();
232 G4double ptot2 = ptot*ptot;
233 G4double cost = 1 - 0.5*std::fabs(tMand)/ptot2;
234
235 if( cost >= 1.0 ) cost = 1.0;
236 else if( cost <= -1.0) cost = -1.0;
237
238 G4double thetaCMS = std::acos(cost);
239
241
242 sigma *= pi/ptot2;
243
244 return sigma;
245}
G4double GetDiffuseElasticXsc(const G4ParticleDefinition *particle, G4double theta, G4double momentum, G4double A)
Definition: G4NuclNuclDiffuseElastic.cc:182
◆ GetLegendrePol()
Definition at line 814 of file G4NuclNuclDiffuseElastic.hh.
815{
816 G4double legPol, epsilon = 1.e-6;
817 G4double x = std::cos(theta);
818
819 if ( n < 0 ) legPol = 0.;
820 else if( n == 0 ) legPol = 1.;
821 else if( n == 1 ) legPol = x;
822 else if( n == 2 ) legPol = (3.*x*x-1.)/2.;
823 else if( n == 3 ) legPol = (5.*x*x*x-3.*x)/2.;
824 else if( n == 4 ) legPol = (35.*x*x*x*x-30.*x*x+3.)/8.;
825 else if( n == 5 ) legPol = (63.*x*x*x*x*x-70.*x*x*x+15.*x)/8.;
826 else if( n == 6 ) legPol = (231.*x*x*x*x*x*x-315.*x*x*x*x+105.*x*x-5.)/16.;
827 else
828 {
829 // legPol = ( (2*n-1)*x*GetLegendrePol(n-1,x) - (n-1)*GetLegendrePol(n-2,x) )/n;
830
831 legPol = std::sqrt( 2./(n*CLHEP::pi*std::sin(theta+epsilon)) )*std::sin( (n+0.5)*theta+0.25*CLHEP::pi );
832 }
833 return legPol;
834}
Referenced by AmplitudeGla().
◆ GetNuclearRadius()
|
inline |
Definition at line 188 of file G4NuclNuclDiffuseElastic.hh.
188{return fNuclearRadius;};
◆ GetProfileLambda()
|
inline |
Definition at line 282 of file G4NuclNuclDiffuseElastic.hh.
282{return fProfileLambda;};
◆ GetRatioGen()
Definition at line 1379 of file G4NuclNuclDiffuseElastic.hh.
1380{
1381 G4double sinThetaR = 2.*fHalfRutThetaTg/(1. + fHalfRutThetaTg2);
1382 G4double dTheta = 0.5*(theta - fRutherfordTheta);
1383 G4double sindTheta = std::sin(dTheta);
1384
1386 G4double prof2 = prof*prof;
1387 // G4double profmod = std::abs(prof);
1388 G4double order = std::sqrt(fProfileLambda/sinThetaR/CLHEP::pi)*2.*sindTheta;
1389
1390 order = std::abs(order); // since sin changes sign!
1391 // G4cout<<"order = "<<order<<G4endl;
1392
1395
1396 G4double out;
1397
1398 if ( theta <= fRutherfordTheta )
1399 {
1400 out = 1. + 0.5*( (0.5-cosFresnel)*(0.5-cosFresnel)+(0.5-sinFresnel)*(0.5-sinFresnel) )*prof2;
1401 out += ( cosFresnel + sinFresnel - 1. )*prof;
1402 }
1403 else
1404 {
1405 out = 0.5*( (0.5-cosFresnel)*(0.5-cosFresnel)+(0.5-sinFresnel)*(0.5-sinFresnel) )*prof2;
1406 }
1407
1408 return out;
1409}
G4double Profile(G4double theta)
Definition: G4NuclNuclDiffuseElastic.hh:1154
G4double GetSint(G4double x)
Definition: G4NuclNuclDiffuseElastic.hh:1027
G4double GetCint(G4double x)
Definition: G4NuclNuclDiffuseElastic.hh:1012
Referenced by GetFresnelDiffuseXsc().
◆ GetRatioSim()
Definition at line 1359 of file G4NuclNuclDiffuseElastic.hh.
1360{
1361 G4double sinThetaR = 2.*fHalfRutThetaTg/(1. + fHalfRutThetaTg2);
1362 G4double dTheta = 0.5*(theta - fRutherfordTheta);
1363 G4double sindTheta = std::sin(dTheta);
1364
1365 G4double order = std::sqrt(fProfileLambda/sinThetaR/CLHEP::pi)*2.*sindTheta;
1366 // G4cout<<"order = "<<order<<G4endl;
1369
1370 G4double out = 0.5*( cosFresnel*cosFresnel + sinFresnel*sinFresnel );
1371
1372 return out;
1373}
◆ GetRutherfordXsc()
Definition at line 648 of file G4NuclNuclDiffuseElastic.hh.
649{
650 G4double sinHalfTheta = std::sin(0.5*theta);
651 G4double sinHalfTheta2 = sinHalfTheta*sinHalfTheta;
652
653 G4double ch2 = fRutherfordRatio*fRutherfordRatio;
654
655 G4double xsc = ch2/(sinHalfTheta2+fAm)/(sinHalfTheta2+fAm);
656
657 return xsc;
658}
Referenced by GetFresnelDiffuseXsc().
◆ GetScatteringAngle()
| G4double G4NuclNuclDiffuseElastic::GetScatteringAngle | ( | G4int | iMomentum, |
| G4int | iAngle, | ||
| G4double | position | ||
| ) |
Definition at line 1035 of file G4NuclNuclDiffuseElastic.cc.
1036{
1037 G4double x1, x2, y1, y2, randAngle;
1038
1039 if( iAngle == 0 )
1040 {
1041 randAngle = (*fAngleTable)(iMomentum)->GetLowEdgeEnergy(iAngle);
1042 // iAngle++;
1043 }
1044 else
1045 {
1047 {
1048 iAngle = (*fAngleTable)(iMomentum)->GetVectorLength() - 1;
1049 }
1050 y1 = (*(*fAngleTable)(iMomentum))(iAngle-1);
1051 y2 = (*(*fAngleTable)(iMomentum))(iAngle);
1052
1053 x1 = (*fAngleTable)(iMomentum)->GetLowEdgeEnergy(iAngle-1);
1054 x2 = (*fAngleTable)(iMomentum)->GetLowEdgeEnergy(iAngle);
1055
1056 if ( x1 == x2 ) randAngle = x2;
1057 else
1058 {
1060 else
1061 {
1062 randAngle = x1 + ( position - y1 )*( x2 - x1 )/( y2 - y1 );
1063 }
1064 }
1065 }
1066 return randAngle;
1067}
Definition: xmltok.h:110
Referenced by SampleTableThetaCMS().
◆ GetSinHaPit2()
Definition at line 199 of file G4NuclNuclDiffuseElastic.hh.
199{return std::sin(CLHEP::halfpi*t*t);};
Referenced by GetSint().
◆ GetSint()
Definition at line 1027 of file G4NuclNuclDiffuseElastic.hh.
1028{
1029 G4double out;
1030
1031 G4Integrator<G4NuclNuclDiffuseElastic,G4double(G4NuclNuclDiffuseElastic::*)(G4double)> integral;
1032
1034
1035 return out;
1036}
G4double GetSinHaPit2(G4double t)
Definition: G4NuclNuclDiffuseElastic.hh:199
Referenced by GetRatioGen(), and GetRatioSim().
◆ InitDynParameters()
|
inline |
Definition at line 1576 of file G4NuclNuclDiffuseElastic.hh.
1578{
1579 G4double a = 0.;
1582
1583 fWaveVector = partMom/CLHEP::hbarc;
1584
1586
1587 if( z )
1588 {
1589 a = partMom/m1; // beta*gamma for m1
1590 fBeta = a/std::sqrt(1+a*a);
1591 fZommerfeld = CalculateZommerfeld( fBeta, z, fAtomicNumber);
1592 fRutherfordRatio = fZommerfeld/fWaveVector;
1593 fAm = CalculateAm( partMom, fZommerfeld, fAtomicNumber);
1594 }
1596 fProfileDelta = fCofDelta*fProfileLambda;
1597 fProfileAlpha = fCofAlpha*fProfileLambda;
1598
1599 CalculateCoulombPhaseZero();
1600 CalculateRutherfordAnglePar();
1601
1602 return;
1603}
void CalculateRutherfordAnglePar()
Definition: G4NuclNuclDiffuseElastic.hh:1104
void CalculateCoulombPhaseZero()
Definition: G4NuclNuclDiffuseElastic.hh:1077
Referenced by BuildAngleTable().
◆ Initialise()
| void G4NuclNuclDiffuseElastic::Initialise | ( | ) |
Definition at line 142 of file G4NuclNuclDiffuseElastic.cc.
143{
144
145 // fEnergyVector = new G4PhysicsLogVector( theMinEnergy, theMaxEnergy, fEnergyBin );
146
149
150 // projectile radius
151
154
155 for(jEl = 0 ; jEl < numOfEl; ++jEl) // application element loop
156 {
157 fAtomicNumber = (*theElementTable)[jEl]->GetZ(); // atomic number
158 fAtomicWeight = (*theElementTable)[jEl]->GetN(); // number of nucleons
159
160 fNuclearRadius = CalculateNuclearRad(fAtomicWeight);
161 fNuclearRadius += R1;
162
164 {
166 <<(*theElementTable)[jEl]->GetName()<<G4endl;
167 }
168 fElementNumberVector.push_back(fAtomicNumber);
169 fElementNameVector.push_back((*theElementTable)[jEl]->GetName());
170
171 BuildAngleTable();
172 fAngleBank.push_back(fAngleTable);
173 }
174}
static const G4ElementTable * GetElementTable()
Definition: G4Element.cc:399
void BuildAngleTable()
Definition: G4NuclNuclDiffuseElastic.cc:958
G4int GetBaryonNumber() const
Definition: G4ParticleDefinition.hh:139
◆ InitialiseOnFly()
Definition at line 929 of file G4NuclNuclDiffuseElastic.cc.
930{
931 fAtomicNumber = Z; // atomic number
932 fAtomicWeight = A; // number of nucleons
933
936
937 fNuclearRadius = CalculateNuclearRad(fAtomicWeight) + R1;
938
940 {
943 }
944 fElementNumberVector.push_back(fAtomicNumber);
945
946 BuildAngleTable();
947
948 fAngleBank.push_back(fAngleTable);
949
950 return;
951}
Referenced by SampleTableThetaCMS().
◆ InitParameters()
|
inline |
Definition at line 1531 of file G4NuclNuclDiffuseElastic.hh.
1533{
1534 fAtomicNumber = Z; // atomic number
1535 fAtomicWeight = A; // number of nucleons
1536
1537 fNuclearRadius2 = CalculateNuclearRad(fAtomicWeight);
1539 fNuclearRadius1 = CalculateNuclearRad(A1);
1540 // fNuclearRadius = std::sqrt(fNuclearRadius1*fNuclearRadius1+fNuclearRadius2*fNuclearRadius2);
1541 fNuclearRadius = fNuclearRadius1 + fNuclearRadius2;
1542
1543 G4double a = 0.;
1546
1547 fWaveVector = partMom/CLHEP::hbarc;
1548
1551
1552 if( z )
1553 {
1554 a = partMom/m1; // beta*gamma for m1
1555 fBeta = a/std::sqrt(1+a*a);
1556 fZommerfeld = CalculateZommerfeld( fBeta, z, fAtomicNumber);
1557 fRutherfordRatio = fZommerfeld/fWaveVector;
1558 fAm = CalculateAm( partMom, fZommerfeld, fAtomicNumber);
1559 }
1563 fProfileDelta = fCofDelta*fProfileLambda;
1564 fProfileAlpha = fCofAlpha*fProfileLambda;
1565
1566 CalculateCoulombPhaseZero();
1567 CalculateRutherfordAnglePar();
1568
1569 return;
1570}
◆ InitParametersGla()
|
inline |
Definition at line 1611 of file G4NuclNuclDiffuseElastic.hh.
1613{
1614 fAtomicNumber = Z; // target atomic number
1615 fAtomicWeight = A; // target number of nucleons
1616
1620 fNuclearRadiusSquare = fNuclearRadius1*fNuclearRadius1+fNuclearRadius2*fNuclearRadius2;
1621
1622
1623 G4double a = 0., kR12;
1626
1627 fWaveVector = partMom/CLHEP::hbarc;
1628
1629 G4double pN = A1 - z;
1630 if( pN < 0. ) pN = 0.;
1631
1632 G4double tN = A - Z;
1633 if( tN < 0. ) tN = 0.;
1634
1636 pTkin /= A1;
1637
1638
1639 fSumSigma = (Z*z+pN*tN)*GetHadronNucleonXscNS(theProton, pTkin, theProton) +
1640 (z*tN+pN*Z)*GetHadronNucleonXscNS(theProton, pTkin, theNeutron);
1641
1644 kR12 = fWaveVector*std::sqrt(fNuclearRadiusSquare);
1646 fMaxL = (G4int(kR12)+1)*4;
1648
1649 if( z )
1650 {
1651 a = partMom/m1; // beta*gamma for m1
1652 fBeta = a/std::sqrt(1+a*a);
1653 fZommerfeld = CalculateZommerfeld( fBeta, z, fAtomicNumber);
1654 fAm = CalculateAm( partMom, fZommerfeld, fAtomicNumber);
1655 }
1656
1657 CalculateCoulombPhaseZero();
1658
1659
1660 return;
1661}
G4ParticleDefinition * GetDefinition() const
G4double GetKineticEnergy() const
G4double GetHadronNucleonXscNS(G4ParticleDefinition *pParticle, G4double pTkin, G4ParticleDefinition *tParticle)
Definition: G4NuclNuclDiffuseElastic.hh:1671
◆ IntegralElasticProb()
| G4double G4NuclNuclDiffuseElastic::IntegralElasticProb | ( | const G4ParticleDefinition * | particle, |
| G4double | theta, | ||
| G4double | momentum, | ||
| G4double | A | ||
| ) |
Definition at line 673 of file G4NuclNuclDiffuseElastic.cc.
677{
678 G4double result;
679 fParticle = particle;
680 fWaveVector = momentum/hbarc;
681 fAtomicWeight = A;
682
683 fNuclearRadius = CalculateNuclearRad(A);
684
685
686 G4Integrator<G4NuclNuclDiffuseElastic,G4double(G4NuclNuclDiffuseElastic::*)(G4double)> integral;
687
688 // result = integral.Legendre10(this,&G4NuclNuclDiffuseElastic::GetIntegrandFunction, 0., theta );
690
691 return result;
692}
G4double GetIntegrandFunction(G4double theta)
Definition: G4NuclNuclDiffuseElastic.cc:657
◆ PhaseFar()
Definition at line 1191 of file G4NuclNuclDiffuseElastic.hh.
1192{
1193 G4double twosigma = 2.*fCoulombPhase0;
1194 twosigma -= fZommerfeld*std::log(fHalfRutThetaTg2/(1.+fHalfRutThetaTg2));
1195 twosigma += fRutherfordTheta*fZommerfeld/fHalfRutThetaTg - CLHEP::halfpi;
1196 twosigma += fProfileLambda*theta - 0.25*CLHEP::pi;
1197
1198 twosigma *= fCofPhase;
1199
1201
1202 return std::exp(z);
1203}
Referenced by AmplitudeFar().
◆ PhaseNear()
Definition at line 1173 of file G4NuclNuclDiffuseElastic.hh.
1174{
1175 G4double twosigma = 2.*fCoulombPhase0;
1176 twosigma -= fZommerfeld*std::log(fHalfRutThetaTg2/(1.+fHalfRutThetaTg2));
1177 twosigma += fRutherfordTheta*fZommerfeld/fHalfRutThetaTg - CLHEP::halfpi;
1178 twosigma -= fProfileLambda*theta - 0.25*CLHEP::pi;
1179
1180 twosigma *= fCofPhase;
1181
1183
1184 return std::exp(z);
1185}
Referenced by AmplitudeNear().
◆ Profile()
Definition at line 1154 of file G4NuclNuclDiffuseElastic.hh.
1155{
1156 G4double dTheta = fRutherfordTheta - theta;
1157 G4double result = 0., argument = 0.;
1158
1159 if(std::abs(dTheta) < 0.001) result = 1.;
1160 else
1161 {
1162 argument = fProfileDelta*dTheta;
1163 result = CLHEP::pi*argument;
1164 result /= std::sinh(CLHEP::pi*argument);
1165 }
1166 return result;
1167}
Referenced by GetRatioGen().
◆ ProfileFar()
Definition at line 1138 of file G4NuclNuclDiffuseElastic.hh.
1139{
1140 G4double dTheta = fRutherfordTheta + theta;
1141 G4double argument = fProfileDelta*dTheta;
1142
1143 G4double result = CLHEP::pi*argument*std::exp(fProfileAlpha*argument);
1144 result /= std::sinh(CLHEP::pi*argument);
1145 result /= dTheta;
1146
1147 return result;
1148}
Referenced by AmplitudeFar().
◆ ProfileNear()
Definition at line 1117 of file G4NuclNuclDiffuseElastic.hh.
1118{
1119 G4double dTheta = fRutherfordTheta - theta;
1120 G4double result = 0., argument = 0.;
1121
1122 if(std::abs(dTheta) < 0.001) result = fProfileAlpha*fProfileDelta;
1123 else
1124 {
1125 argument = fProfileDelta*dTheta;
1126 result = CLHEP::pi*argument*std::exp(fProfileAlpha*argument);
1127 result /= std::sinh(CLHEP::pi*argument);
1128 result -= 1.;
1129 result /= dTheta;
1130 }
1131 return result;
1132}
Referenced by AmplitudeNear(), and AmplitudeSim().
◆ SampleInvariantT()
|
virtual |
Reimplemented from G4HadronElastic.
Definition at line 766 of file G4NuclNuclDiffuseElastic.cc.
768{
769 fParticle = aParticle;
771 G4double totElab = std::sqrt(m1*m1+p*p);
773 G4LorentzVector lv1(p,0.0,0.0,totElab);
774 G4LorentzVector lv(0.0,0.0,0.0,mass2);
775 lv += lv1;
776
777 G4ThreeVector bst = lv.boostVector();
778 lv1.boost(-bst);
779
780 G4ThreeVector p1 = lv1.vect();
782
783 G4double t = SampleTableT( aParticle, momentumCMS, G4double(Z), G4double(A) ); // sample theta2 in cms
784 return t;
785}
G4double SampleTableT(const G4ParticleDefinition *aParticle, G4double p, G4double Z, G4double A)
Definition: G4NuclNuclDiffuseElastic.cc:791
static G4double GetNuclearMass(const G4double A, const G4double Z)
Definition: G4NucleiProperties.cc:53
◆ SampleT()
| G4double G4NuclNuclDiffuseElastic::SampleT | ( | const G4ParticleDefinition * | aParticle, |
| G4double | p, | ||
| G4double | A | ||
| ) |
Definition at line 698 of file G4NuclNuclDiffuseElastic.cc.
700{
703 return t;
704}
G4double SampleThetaCMS(const G4ParticleDefinition *aParticle, G4double p, G4double A)
Definition: G4NuclNuclDiffuseElastic.cc:712
Referenced by SampleThetaLab().
◆ SampleTableT()
| G4double G4NuclNuclDiffuseElastic::SampleTableT | ( | const G4ParticleDefinition * | aParticle, |
| G4double | p, | ||
| G4double | Z, | ||
| G4double | A | ||
| ) |
Definition at line 791 of file G4NuclNuclDiffuseElastic.cc.
793{
795 // G4double t = 2*p*p*( 1 - std::cos(std::sqrt(alpha)) ); // -t !!!
797 return t;
798}
G4double SampleTableThetaCMS(const G4ParticleDefinition *aParticle, G4double p, G4double Z, G4double A)
Definition: G4NuclNuclDiffuseElastic.cc:806
Referenced by SampleInvariantT().
◆ SampleTableThetaCMS()
| G4double G4NuclNuclDiffuseElastic::SampleTableThetaCMS | ( | const G4ParticleDefinition * | aParticle, |
| G4double | p, | ||
| G4double | Z, | ||
| G4double | A | ||
| ) |
Definition at line 806 of file G4NuclNuclDiffuseElastic.cc.
808{
809 size_t iElement;
810 G4int iMomentum, iAngle;
812 G4double m1 = particle->GetPDGMass();
813
814 for(iElement = 0; iElement < fElementNumberVector.size(); iElement++)
815 {
816 if( std::fabs(Z - fElementNumberVector[iElement]) < 0.5) break;
817 }
818 if ( iElement == fElementNumberVector.size() )
819 {
821
822 // iElement--;
823
824 // G4cout << "G4NuclNuclDiffuseElastic: Element with atomic number " << Z
825 // << " is not found, return zero angle" << G4endl;
826 // return 0.; // no table for this element
827 }
828 // G4cout<<"iElement = "<<iElement<<G4endl;
829
830 fAngleTable = fAngleBank[iElement];
831
832 G4double kinE = std::sqrt(momentum*momentum + m1*m1) - m1;
833
834 for( iMomentum = 0; iMomentum < fEnergyBin; iMomentum++)
835 {
836 // G4cout<<iMomentum<<", kinE = "<<kinE/MeV<<", vectorE = "<<fEnergyVector->GetLowEdgeEnergy(iMomentum)/MeV<<G4endl;
837
838 if( kinE < fEnergyVector->GetLowEdgeEnergy(iMomentum) ) break;
839 }
840 // G4cout<<"iMomentum = "<<iMomentum<<", fEnergyBin -1 = "<<fEnergyBin -1<<G4endl;
841
842
843 if ( iMomentum >= fEnergyBin ) iMomentum = fEnergyBin-1; // kinE is more then theMaxEnergy
844 if ( iMomentum < 0 ) iMomentum = 0; // against negative index, kinE < theMinEnergy
845
846
847 if (iMomentum == fEnergyBin -1 || iMomentum == 0 ) // the table edges
848 {
850
851 // G4cout<<"position = "<<position<<G4endl;
852
853 for(iAngle = 0; iAngle < fAngleBin-1; iAngle++)
854 {
856 }
857 if (iAngle >= fAngleBin-1) iAngle = fAngleBin-2;
858
859 // G4cout<<"iAngle = "<<iAngle<<G4endl;
860
862
863 // G4cout<<"randAngle = "<<randAngle<<G4endl;
864 }
865 else // kinE inside between energy table edges
866 {
867 // position = (*(*fAngleTable)(iMomentum))(fAngleBin-2)*G4UniformRand();
869
870 // G4cout<<"position = "<<position<<G4endl;
871
872 for(iAngle = 0; iAngle < fAngleBin-1; iAngle++)
873 {
874 // if( position < (*(*fAngleTable)(iMomentum))(iAngle) ) break;
876 }
877 if (iAngle >= fAngleBin-1) iAngle = fAngleBin-2;
878
879 // G4cout<<"iAngle = "<<iAngle<<G4endl;
880
882
883 // G4cout<<"theta2 = "<<theta2<<G4endl;
884
885 E2 = fEnergyVector->GetLowEdgeEnergy(iMomentum);
886
887 // G4cout<<"E2 = "<<E2<<G4endl;
888
889 iMomentum--;
890
891 // position = (*(*fAngleTable)(iMomentum))(fAngleBin-2)*G4UniformRand();
892
893 // G4cout<<"position = "<<position<<G4endl;
894
895 for(iAngle = 0; iAngle < fAngleBin-1; iAngle++)
896 {
897 // if( position < (*(*fAngleTable)(iMomentum))(iAngle) ) break;
899 }
900 if (iAngle >= fAngleBin-1) iAngle = fAngleBin-2;
901
903
904 // G4cout<<"theta1 = "<<theta1<<G4endl;
905
906 E1 = fEnergyVector->GetLowEdgeEnergy(iMomentum);
907
908 // G4cout<<"E1 = "<<E1<<G4endl;
909
910 W = 1.0/(E2 - E1);
911 W1 = (E2 - kinE)*W;
912 W2 = (kinE - E1)*W;
913
914 randAngle = W1*theta1 + W2*theta2;
915
916 // randAngle = theta2;
917 }
918 // G4double angle = randAngle;
919 // if (randAngle > 0.) randAngle /= 2*pi*std::sin(angle);
920 // G4cout<<"randAngle = "<<randAngle/degree<<G4endl;
921
922 return randAngle;
923}
void InitialiseOnFly(G4double Z, G4double A)
Definition: G4NuclNuclDiffuseElastic.cc:929
G4double GetScatteringAngle(G4int iMomentum, G4int iAngle, G4double position)
Definition: G4NuclNuclDiffuseElastic.cc:1035
Referenced by SampleTableT().
◆ SampleThetaCMS()
| G4double G4NuclNuclDiffuseElastic::SampleThetaCMS | ( | const G4ParticleDefinition * | aParticle, |
| G4double | p, | ||
| G4double | A | ||
| ) |
Definition at line 712 of file G4NuclNuclDiffuseElastic.cc.
714{
715 G4int i, iMax = 100;
716 G4double norm, result, theta1, theta2, thetaMax, sum = 0.;
717
718 fParticle = particle;
719 fWaveVector = momentum/hbarc;
720 fAtomicWeight = A;
721
722 fNuclearRadius = CalculateNuclearRad(A);
723
724 thetaMax = 10.174/fWaveVector/fNuclearRadius;
725
727
728 G4Integrator<G4NuclNuclDiffuseElastic,G4double(G4NuclNuclDiffuseElastic::*)(G4double)> integral;
729
730 // result = integral.Legendre10(this,&G4NuclNuclDiffuseElastic::GetIntegrandFunction, 0., theta );
731 norm = integral.Legendre96(this,&G4NuclNuclDiffuseElastic::GetIntegrandFunction, 0., thetaMax );
732
733 norm *= G4UniformRand();
734
735 for(i = 1; i <= iMax; i++)
736 {
737 theta1 = (i-1)*thetaMax/iMax;
738 theta2 = i*thetaMax/iMax;
739 sum += integral.Legendre10(this,&G4NuclNuclDiffuseElastic::GetIntegrandFunction, theta1, theta2);
740
741 if ( sum >= norm )
742 {
743 result = 0.5*(theta1 + theta2);
744 break;
745 }
746 }
747 if (i > iMax ) result = 0.5*(theta1 + theta2);
748
750
751 result += G4RandGauss::shoot(0.,sigma);
752
753 if(result < 0.) result = 0.;
754 if(result > thetaMax) result = thetaMax;
755
756 return result;
757}
Referenced by SampleT().
◆ SampleThetaLab()
| G4double G4NuclNuclDiffuseElastic::SampleThetaLab | ( | const G4HadProjectile * | aParticle, |
| G4double | tmass, | ||
| G4double | A | ||
| ) |
Definition at line 1078 of file G4NuclNuclDiffuseElastic.cc.
1080{
1085 G4LorentzVector lv(0.0,0.0,0.0,tmass);
1086 lv += lv1;
1087
1088 G4ThreeVector bst = lv.boostVector();
1089 lv1.boost(-bst);
1090
1093 G4double tmax = 4.0*ptot*ptot;
1094 G4double t = 0.0;
1095
1096
1097 //
1098 // Sample t
1099 //
1100
1101 t = SampleT( theParticle, ptot, A);
1102
1103 // NaN finder
1104 if(!(t < 0.0 || t >= 0.0))
1105 {
1107 {
1109 << " mom(GeV)= " << plab/GeV
1110 << " S-wave will be sampled"
1111 << G4endl;
1112 }
1113 t = G4UniformRand()*tmax;
1114 }
1116 {
1119 }
1120 // Sampling of angles in CM system
1121
1123 G4double cost = 1. - 2.0*t/tmax;
1124 G4double sint;
1125
1126 if( cost >= 1.0 )
1127 {
1128 cost = 1.0;
1129 sint = 0.0;
1130 }
1131 else if( cost <= -1.0)
1132 {
1133 cost = -1.0;
1134 sint = 0.0;
1135 }
1136 else
1137 {
1138 sint = std::sqrt((1.0-cost)*(1.0+cost));
1139 }
1141 {
1143 }
1144 G4ThreeVector v1(sint*std::cos(phi),sint*std::sin(phi),cost);
1145 v1 *= ptot;
1146 G4LorentzVector nlv1(v1.x(),v1.y(),v1.z(),std::sqrt(ptot*ptot + m1*m1));
1147
1148 nlv1.boost(bst);
1149
1150 G4ThreeVector np1 = nlv1.vect();
1151
1152 // G4double theta = std::acos( np1.z()/np1.mag() ); // degree;
1153
1155
1156 return theta;
1157}
double theta() const
HepLorentzVector & boost(double, double, double)
Definition: LorentzVector.cc:59
Hep3Vector vect() const
G4double GetTotalMomentum() const
Definition: G4HadProjectile.hh:101
const G4ParticleDefinition * GetDefinition() const
Definition: G4HadProjectile.hh:81
const G4LorentzVector & Get4Momentum() const
Definition: G4HadProjectile.hh:86
G4double SampleT(const G4ParticleDefinition *aParticle, G4double p, G4double A)
Definition: G4NuclNuclDiffuseElastic.cc:698
◆ SetCofAlpha()
|
inline |
Definition at line 289 of file G4NuclNuclDiffuseElastic.hh.
289{fCofAlpha = pa;};
◆ SetCofAlphaCoulomb()
|
inline |
Definition at line 291 of file G4NuclNuclDiffuseElastic.hh.
291{fCofAlphaCoulomb = pa;};
◆ SetCofAlphaMax()
|
inline |
Definition at line 290 of file G4NuclNuclDiffuseElastic.hh.
290{fCofAlphaMax = pa;};
◆ SetCofDelta()
|
inline |
Definition at line 293 of file G4NuclNuclDiffuseElastic.hh.
293{fCofDelta = pa;};
◆ SetCofFar()
|
inline |
Definition at line 295 of file G4NuclNuclDiffuseElastic.hh.
295{fCofFar = pa;};
◆ SetCofLambda()
|
inline |
Definition at line 287 of file G4NuclNuclDiffuseElastic.hh.
287{fCofLambda = pa;};
◆ SetCofPhase()
|
inline |
Definition at line 294 of file G4NuclNuclDiffuseElastic.hh.
294{fCofPhase = pa;};
◆ SetEtaRatio()
|
inline |
Definition at line 296 of file G4NuclNuclDiffuseElastic.hh.
296{fEtaRatio = pa;};
◆ SetHEModelLowLimit()
|
inline |
Definition at line 385 of file G4NuclNuclDiffuseElastic.hh.
386{
387 lowEnergyLimitHE = value;
388}
◆ SetLowestEnergyLimit()
|
inline |
Definition at line 395 of file G4NuclNuclDiffuseElastic.hh.
396{
397 lowestEnergyLimit = value;
398}
◆ SetMaxL()
|
inline |
Definition at line 297 of file G4NuclNuclDiffuseElastic.hh.
297{fMaxL = l;};
◆ SetNuclearRadiusCof()
|
inline |
Definition at line 298 of file G4NuclNuclDiffuseElastic.hh.
298{fNuclearRadiusCof = r;};
◆ SetPlabLowLimit()
|
inline |
Definition at line 380 of file G4NuclNuclDiffuseElastic.hh.
381{
382 plabLowLimit = value;
383}
◆ SetProfileAlpha()
|
inline |
Definition at line 286 of file G4NuclNuclDiffuseElastic.hh.
286{fProfileAlpha = pa;};
◆ SetProfileDelta()
|
inline |
Definition at line 285 of file G4NuclNuclDiffuseElastic.hh.
285{fProfileDelta = pd;};
◆ SetProfileLambda()
|
inline |
Definition at line 284 of file G4NuclNuclDiffuseElastic.hh.
284{fProfileLambda = pl;};
◆ SetQModelLowLimit()
|
inline |
Definition at line 390 of file G4NuclNuclDiffuseElastic.hh.
391{
392 lowEnergyLimitQ = value;
393}
◆ SetRecoilKinEnergyLimit()
|
inline |
Definition at line 375 of file G4NuclNuclDiffuseElastic.hh.
376{
377 lowEnergyRecoilLimit = value;
378}
◆ TestAngleTable()
| void G4NuclNuclDiffuseElastic::TestAngleTable | ( | const G4ParticleDefinition * | theParticle, |
| G4double | partMom, | ||
| G4double | Z, | ||
| G4double | A | ||
| ) |
Definition at line 1286 of file G4NuclNuclDiffuseElastic.cc.
1288{
1289 fAtomicNumber = Z; // atomic number
1290 fAtomicWeight = A; // number of nucleons
1291 fNuclearRadius = CalculateNuclearRad(fAtomicWeight);
1292
1293
1294
1297
1298 fElementNumberVector.push_back(fAtomicNumber);
1299
1300
1301
1302
1303 G4int i=0, j;
1305 G4double alpha1=0., alpha2=0., alphaMax=0., alphaCoulomb=0.;
1306 G4double deltaL10 = 0., deltaL96 = 0., deltaAG = 0.;
1307 G4double sumL10 = 0.,sumL96 = 0.,sumAG = 0.;
1308 G4double epsilon = 0.001;
1309
1310 G4Integrator<G4NuclNuclDiffuseElastic,G4double(G4NuclNuclDiffuseElastic::*)(G4double)> integral;
1311
1313
1314 fWaveVector = partMom/hbarc;
1315
1316 G4double kR = fWaveVector*fNuclearRadius;
1317 G4double kR2 = kR*kR;
1320
1321 alphaMax = kRmax*kRmax/kR2;
1322
1323 if (alphaMax > 4.) alphaMax = 4.; // vmg05-02-09: was pi2
1324
1325 alphaCoulomb = kRcoul*kRcoul/kR2;
1326
1327 if( z )
1328 {
1329 a = partMom/m1; // beta*gamma for m1
1330 fBeta = a/std::sqrt(1+a*a);
1331 fZommerfeld = CalculateZommerfeld( fBeta, z, fAtomicNumber);
1332 fAm = CalculateAm( partMom, fZommerfeld, fAtomicNumber);
1333 }
1335
1336 // G4PhysicsLogVector* angleBins = new G4PhysicsLogVector( 0.001*alphaMax, alphaMax, fAngleBin );
1337
1338
1339 fAddCoulomb = false;
1340
1341 for(j = 1; j < fAngleBin; j++)
1342 {
1343 // alpha1 = angleBins->GetLowEdgeEnergy(j-1);
1344 // alpha2 = angleBins->GetLowEdgeEnergy(j);
1345
1346 alpha1 = alphaMax*(j-1)/fAngleBin;
1347 alpha2 = alphaMax*( j )/fAngleBin;
1348
1349 if( ( alpha2 > alphaCoulomb ) && z ) fAddCoulomb = true;
1350
1351 deltaL10 = integral.Legendre10(this, &G4NuclNuclDiffuseElastic::GetIntegrandFunction, alpha1, alpha2);
1352 deltaL96 = integral.Legendre96(this, &G4NuclNuclDiffuseElastic::GetIntegrandFunction, alpha1, alpha2);
1354 alpha1, alpha2,epsilon);
1355
1356 // G4cout<<alpha1<<"\t"<<std::sqrt(alpha1)/degree<<"\t"
1357 // <<deltaL10<<"\t"<<deltaL96<<"\t"<<deltaAG<<G4endl;
1358
1359 sumL10 += deltaL10;
1360 sumL96 += deltaL96;
1361 sumAG += deltaAG;
1362
1365
1367 }
1368 fAngleTable->insertAt(i,angleVector);
1369 fAngleBank.push_back(fAngleTable);
1370
1371 /*
1372 // Integral over all angle range - Bad accuracy !!!
1373
1374 sumL10 = integral.Legendre10(this, &G4NuclNuclDiffuseElastic::GetIntegrandFunction, 0., alpha2);
1375 sumL96 = integral.Legendre96(this, &G4NuclNuclDiffuseElastic::GetIntegrandFunction, 0., alpha2);
1376 sumAG = integral.AdaptiveGauss(this, &G4NuclNuclDiffuseElastic::GetIntegrandFunction,
1377 0., alpha2,epsilon);
1378 G4cout<<G4endl;
1379 G4cout<<alpha2<<"\t"<<std::sqrt(alpha2)/degree<<"\t"
1380 <<sumL10<<"\t"<<sumL96<<"\t"<<sumAG<<G4endl;
1381 */
1382 return;
1383}
◆ TestErfcComp()
Definition at line 842 of file G4NuclNuclDiffuseElastic.hh.
◆ TestErfcInt()
◆ TestErfcSer()
◆ ThetaCMStoThetaLab()
| G4double G4NuclNuclDiffuseElastic::ThetaCMStoThetaLab | ( | const G4DynamicParticle * | aParticle, |
| G4double | tmass, | ||
| G4double | thetaCMS | ||
| ) |
Definition at line 1166 of file G4NuclNuclDiffuseElastic.cc.
1168{
1171 // G4double plab = aParticle->GetTotalMomentum();
1173 G4LorentzVector lv(0.0,0.0,0.0,tmass);
1174
1175 lv += lv1;
1176
1177 G4ThreeVector bst = lv.boostVector();
1178
1179 lv1.boost(-bst);
1180
1183
1185 G4double cost = std::cos(thetaCMS);
1186 G4double sint;
1187
1188 if( cost >= 1.0 )
1189 {
1190 cost = 1.0;
1191 sint = 0.0;
1192 }
1193 else if( cost <= -1.0)
1194 {
1195 cost = -1.0;
1196 sint = 0.0;
1197 }
1198 else
1199 {
1200 sint = std::sqrt((1.0-cost)*(1.0+cost));
1201 }
1203 {
1205 }
1206 G4ThreeVector v1(sint*std::cos(phi),sint*std::sin(phi),cost);
1207 v1 *= ptot;
1208 G4LorentzVector nlv1(v1.x(),v1.y(),v1.z(),std::sqrt(ptot*ptot + m1*m1));
1209
1210 nlv1.boost(bst);
1211
1212 G4ThreeVector np1 = nlv1.vect();
1213
1214
1216
1217 return thetaLab;
1218}
G4LorentzVector Get4Momentum() const
◆ ThetaLabToThetaCMS()
| G4double G4NuclNuclDiffuseElastic::ThetaLabToThetaCMS | ( | const G4DynamicParticle * | aParticle, |
| G4double | tmass, | ||
| G4double | thetaLab | ||
| ) |
Definition at line 1227 of file G4NuclNuclDiffuseElastic.cc.
1229{
1234 G4LorentzVector lv(0.0,0.0,0.0,tmass);
1235
1236 lv += lv1;
1237
1238 G4ThreeVector bst = lv.boostVector();
1239
1240 // lv1.boost(-bst);
1241
1242 // G4ThreeVector p1 = lv1.vect();
1243 // G4double ptot = p1.mag();
1244
1246 G4double cost = std::cos(thetaLab);
1247 G4double sint;
1248
1249 if( cost >= 1.0 )
1250 {
1251 cost = 1.0;
1252 sint = 0.0;
1253 }
1254 else if( cost <= -1.0)
1255 {
1256 cost = -1.0;
1257 sint = 0.0;
1258 }
1259 else
1260 {
1261 sint = std::sqrt((1.0-cost)*(1.0+cost));
1262 }
1264 {
1266 }
1267 G4ThreeVector v1(sint*std::cos(phi),sint*std::sin(phi),cost);
1268 v1 *= plab;
1269 G4LorentzVector nlv1(v1.x(),v1.y(),v1.z(),std::sqrt(plab*plab + m1*m1));
1270
1271 nlv1.boost(-bst);
1272
1273 G4ThreeVector np1 = nlv1.vect();
1274
1275
1277
1278 return thetaCMS;
1279}
G4double GetTotalMomentum() const
The documentation for this class was generated from the following files:
