Geant4: G4NuclNuclDiffuseElastic Class Reference

G4NuclNuclDiffuseElastic::G4NuclNuclDiffuseElastic ( )

Definition at line 67 of file G4NuclNuclDiffuseElastic.cc.

References G4Alpha::Alpha(), G4Deuteron::Deuteron(), G4Neutron::Neutron(), G4PionMinus::PionMinus(), G4PionPlus::PionPlus(), G4Proton::Proton(), G4HadronicInteraction::SetMaxEnergy(), G4HadronicInteraction::SetMinEnergy(), G4HadronicInteraction::theMaxEnergy, G4HadronicInteraction::theMinEnergy, and G4HadronicInteraction::verboseLevel.

00068   : G4HadronElastic("NNDiffuseElastic"), fParticle(0)
00069 {
00070   SetMinEnergy( 50*MeV );
00071   SetMaxEnergy( 1.*TeV );
00072   verboseLevel = 0;
00073   lowEnergyRecoilLimit = 100.*keV;  
00074   lowEnergyLimitQ  = 0.0*GeV;  
00075   lowEnergyLimitHE = 0.0*GeV;  
00076   lowestEnergyLimit= 0.0*keV;  
00077   plabLowLimit     = 20.0*MeV;
00078 
00079   theProton   = G4Proton::Proton();
00080   theNeutron  = G4Neutron::Neutron();
00081   theDeuteron = G4Deuteron::Deuteron();
00082   theAlpha    = G4Alpha::Alpha();
00083   thePionPlus = G4PionPlus::PionPlus();
00084   thePionMinus= G4PionMinus::PionMinus();
00085 
00086   fEnergyBin = 200;
00087   fAngleBin  = 200;
00088 
00089   fEnergyVector =  new G4PhysicsLogVector( theMinEnergy, theMaxEnergy, fEnergyBin );
00090   fAngleTable = 0;
00091 
00092   fParticle = 0;
00093   fWaveVector = 0.;
00094   fAtomicWeight = 0.;
00095   fAtomicNumber = 0.;
00096   fNuclearRadius = 0.;
00097   fBeta = 0.;
00098   fZommerfeld = 0.;
00099   fAm = 0.;
00100   fAddCoulomb = false;
00101   // Ranges of angle table relative to current Rutherford (Coulomb grazing) angle
00102 
00103   // Empirical parameters
00104 
00105   fCofAlphaMax     = 1.5;
00106   fCofAlphaCoulomb = 0.5;
00107 
00108   fProfileDelta  = 1.;
00109   fProfileAlpha  = 0.5;
00110 
00111   fCofLambda = 1.0;
00112   fCofDelta  = 0.04;
00113   fCofAlpha  = 0.095;
00114 
00115   fNuclearRadius1 = fNuclearRadius2 = fNuclearRadiusSquare = fNuclearRadiusCof 
00116     = fRutherfordRatio = fCoulombPhase0 = fHalfRutThetaTg = fHalfRutThetaTg2 
00117     = fRutherfordTheta = fProfileLambda = fCofPhase = fCofFar = fCofAlphaMax 
00118     = fCofAlphaCoulomb = fSumSigma = fEtaRatio = fReZ = 0.0;
00119   fMaxL = 0;
00120 
00121 }

G4NuclNuclDiffuseElastic::~G4NuclNuclDiffuseElastic ( ) [virtual]

Definition at line 127 of file G4NuclNuclDiffuseElastic.cc.

References G4PhysicsTable::clearAndDestroy().

00128 {
00129   if(fEnergyVector) delete fEnergyVector;
00130 
00131   if( fAngleTable )
00132   {
00133       fAngleTable->clearAndDestroy();
00134       delete fAngleTable ;
00135   }
00136 }

G4complex G4NuclNuclDiffuseElastic::Amplitude ( G4double theta ) [inline]

Definition at line 1304 of file G4NuclNuclDiffuseElastic.hh.

References AmplitudeFar(), and AmplitudeNear().

Referenced by AmplitudeMod2().

01305 {
01306  
01307   G4complex out = AmplitudeNear(theta) + fCofFar*AmplitudeFar(theta);
01308   // G4complex out = AmplitudeNear(theta);
01309   // G4complex out = AmplitudeFar(theta);
01310   return    out;
01311 }

G4complex G4NuclNuclDiffuseElastic::AmplitudeFar ( G4double theta ) [inline]

Definition at line 1290 of file G4NuclNuclDiffuseElastic.hh.

References PhaseFar(), G4INCL::Math::pi, and ProfileFar().

Referenced by Amplitude().

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 }

G4complex G4NuclNuclDiffuseElastic::AmplitudeGG ( G4double theta ) [inline]

Definition at line 1489 of file G4NuclNuclDiffuseElastic.hh.

References CoulombAmplitude(), G4cout, G4endl, CLHEP::detail::n, and G4INCL::Math::pi.

Referenced by AmplitudeGGMod2().

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 }

G4double G4NuclNuclDiffuseElastic::AmplitudeGGMod2 ( G4double theta ) [inline]

Definition at line 1518 of file G4NuclNuclDiffuseElastic.hh.

References AmplitudeGG().

01519 {
01520   G4complex out = AmplitudeGG(theta);
01521   G4double mod2 = out.real()*out.real() + out.imag()*out.imag();
01522   return   mod2;
01523 }

G4complex G4NuclNuclDiffuseElastic::AmplitudeGla ( G4double theta ) [inline]

Definition at line 1451 of file G4NuclNuclDiffuseElastic.hh.

References CalculateCoulombPhase(), CoulombAmplitude(), GetLegendrePol(), CLHEP::detail::n, and G4INCL::Math::pi.

Referenced by AmplitudeGlaMod2().

01452 {
01453   G4int n;
01454   G4double T12b, b, b2; // cosTheta = std::cos(theta);
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     // b = ( fZommerfeld + std::sqrt( fZommerfeld*fZommerfeld + n*(n+1) ) )/fWaveVector;
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 }

G4double G4NuclNuclDiffuseElastic::AmplitudeGlaMod2 ( G4double theta ) [inline]

Definition at line 1477 of file G4NuclNuclDiffuseElastic.hh.

References AmplitudeGla().

01478 {
01479   G4complex out = AmplitudeGla(theta);
01480   G4double mod2 = out.real()*out.real() + out.imag()*out.imag();
01481   return   mod2;
01482 }

G4double G4NuclNuclDiffuseElastic::AmplitudeMod2 ( G4double theta ) [inline]

Definition at line 1317 of file G4NuclNuclDiffuseElastic.hh.

References Amplitude().

01318 {
01319   G4complex out = Amplitude(theta);
01320   G4double mod2 = out.real()*out.real() + out.imag()*out.imag();
01321   return   mod2;
01322 }

G4complex G4NuclNuclDiffuseElastic::AmplitudeNear ( G4double theta ) [inline]

Definition at line 1265 of file G4NuclNuclDiffuseElastic.hh.

References CoulombAmplitude(), GammaLess(), GammaMore(), PhaseNear(), G4INCL::Math::pi, and ProfileNear().

Referenced by Amplitude().

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     // out *= GammaMore(theta) + ProfileNear(theta);
01276     out += CoulombAmplitude(theta);
01277   }
01278   else
01279   {
01280     out *= GammaMore(theta) + ProfileNear(theta);
01281     // out *= GammaLess(theta) + ProfileNear(theta);
01282   }
01283   return out;
01284 }

G4complex G4NuclNuclDiffuseElastic::AmplitudeSim ( G4double theta ) [inline]

Definition at line 1328 of file G4NuclNuclDiffuseElastic.hh.

References CoulombAmplitude(), GetErfcInt(), and ProfileNear().

Referenced by AmplitudeSimMod2().

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   // order              *= std::sqrt(0.5*fProfileLambda/sinThetaR)*2.*dTheta;
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 }

G4double G4NuclNuclDiffuseElastic::AmplitudeSimMod2 ( G4double theta ) [inline]

Definition at line 1440 of file G4NuclNuclDiffuseElastic.hh.

References AmplitudeSim().

01441 {
01442   G4complex out = AmplitudeSim(theta);
01443   G4double mod2 = out.real()*out.real() + out.imag()*out.imag();
01444   return   mod2;
01445 }

G4double G4NuclNuclDiffuseElastic::BesselJone ( G4double z ) [inline]

Definition at line 458 of file G4NuclNuclDiffuseElastic.hh.

Referenced by BesselOneByArg(), GetDiffElasticProb(), GetDiffElasticSumProb(), and GetDiffElasticSumProbA().

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 }

G4double G4NuclNuclDiffuseElastic::BesselJzero ( G4double z ) [inline]

Definition at line 406 of file G4NuclNuclDiffuseElastic.hh.

Referenced by GetDiffElasticProb(), GetDiffElasticSumProb(), and GetDiffElasticSumProbA().

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 }

G4double G4NuclNuclDiffuseElastic::BesselOneByArg ( G4double z ) [inline]

Definition at line 533 of file G4NuclNuclDiffuseElastic.hh.

References BesselJone().

Referenced by GetDiffElasticProb(), GetDiffElasticSumProb(), and GetDiffElasticSumProbA().

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 }

void G4NuclNuclDiffuseElastic::BuildAngleTable ( )

Definition at line 958 of file G4NuclNuclDiffuseElastic.cc.

References GetFresnelIntegrandXsc(), G4PhysicsVector::GetLowEdgeEnergy(), G4ParticleDefinition::GetPDGMass(), InitDynParameters(), G4PhysicsTable::insertAt(), G4Integrator< T, F >::Legendre10(), G4INCL::Math::pi, and G4PhysicsFreeVector::PutValue().

Referenced by Initialise(), and InitialiseOnFly().

00959 {
00960   G4int i, j;
00961   G4double partMom, kinE, m1 = fParticle->GetPDGMass();
00962   G4double alpha1, alpha2, alphaMax, alphaCoulomb, delta = 0., sum = 0.;
00963 
00964   // G4cout<<"particle z = "<<z<<"; particle m1 = "<<m1/GeV<<" GeV"<<G4endl;
00965 
00966   G4Integrator<G4NuclNuclDiffuseElastic,G4double(G4NuclNuclDiffuseElastic::*)(G4double)> integral;
00967   
00968   fAngleTable = new G4PhysicsTable(fEnergyBin);
00969 
00970   for( i = 0; i < fEnergyBin; i++)
00971   {
00972     kinE        = fEnergyVector->GetLowEdgeEnergy(i);
00973 
00974     // G4cout<<G4endl;
00975     // G4cout<<"kinE = "<<kinE/MeV<<" MeV"<<G4endl;
00976 
00977     partMom     = std::sqrt( kinE*(kinE + 2*m1) );
00978 
00979     InitDynParameters(fParticle, partMom);
00980 
00981     alphaMax = fRutherfordTheta*fCofAlphaMax;
00982 
00983     if(alphaMax > pi) alphaMax = pi;
00984 
00985 
00986     alphaCoulomb = fRutherfordTheta*fCofAlphaCoulomb;
00987 
00988     // G4cout<<"alphaCoulomb = "<<alphaCoulomb/degree<<"; alphaMax = "<<alphaMax/degree<<G4endl;
00989 
00990     G4PhysicsFreeVector* angleVector = new G4PhysicsFreeVector(fAngleBin-1);
00991 
00992     // G4PhysicsLogVector*  angleBins = new G4PhysicsLogVector( 0.001*alphaMax, alphaMax, fAngleBin );
00993 
00994     G4double delth = (alphaMax-alphaCoulomb)/fAngleBin;
00995         
00996     sum = 0.;
00997 
00998     // fAddCoulomb = false;
00999     fAddCoulomb = true;
01000 
01001     // for(j = 1; j < fAngleBin; j++)
01002     for(j = fAngleBin-1; j >= 1; j--)
01003     {
01004       // alpha1 = angleBins->GetLowEdgeEnergy(j-1);
01005       // alpha2 = angleBins->GetLowEdgeEnergy(j);
01006 
01007       // alpha1 = alphaMax*(j-1)/fAngleBin;
01008       // alpha2 = alphaMax*( j )/fAngleBin;
01009 
01010       alpha1 = alphaCoulomb + delth*(j-1);
01011       // if(alpha1 < kRlim2) alpha1 = kRlim2;
01012       alpha2 = alpha1 + delth;
01013 
01014       delta = integral.Legendre10(this, &G4NuclNuclDiffuseElastic::GetFresnelIntegrandXsc, alpha1, alpha2);
01015       // delta = integral.Legendre96(this, &G4NuclNuclDiffuseElastic::GetIntegrandFunction, alpha1, alpha2);
01016 
01017       sum += delta;
01018       
01019       angleVector->PutValue( j-1 , alpha1, sum ); // alpha2
01020       // G4cout<<"j-1 = "<<j-1<<"; alpha2 = "<<alpha2/degree<<"; sum = "<<sum<<G4endl;
01021     }
01022     fAngleTable->insertAt(i,angleVector);
01023 
01024     // delete[] angleVector; 
01025     // delete[] angleBins; 
01026   }
01027   return;
01028 }

G4double G4NuclNuclDiffuseElastic::CalcMandelstamS	(	const G4double	mp,
		const G4double	mt,
		const G4double	Plab
	)			`[inline]`

Definition at line 1781 of file G4NuclNuclDiffuseElastic.hh.

Referenced by GetHadronNucleonXscNS().

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 }

G4double G4NuclNuclDiffuseElastic::CalculateAm	(	G4double	momentum,
		G4double	n,
		G4double	Z
	)			`[inline]`

Definition at line 579 of file G4NuclNuclDiffuseElastic.hh.

Referenced by GetCoulombElasticXsc(), GetCoulombIntegralXsc(), GetCoulombTotalXsc(), GetDiffuseElasticSumXsc(), InitDynParameters(), InitParameters(), InitParametersGla(), and TestAngleTable().

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 }

G4double G4NuclNuclDiffuseElastic::CalculateCoulombPhase ( G4int n ) [inline]

Definition at line 1090 of file G4NuclNuclDiffuseElastic.hh.

References GammaLogB2n().

Referenced by AmplitudeGla().

01091 {
01092   G4complex z        = G4complex(1. + n, fZommerfeld); 
01093   // G4complex gammalog = GammaLogarithm(z);
01094   G4complex gammalog = GammaLogB2n(z);
01095   return gammalog.imag();
01096 }

void G4NuclNuclDiffuseElastic::CalculateCoulombPhaseZero ( ) [inline]

Definition at line 1077 of file G4NuclNuclDiffuseElastic.hh.

References GammaLogB2n().

Referenced by InitDynParameters(), InitParameters(), and InitParametersGla().

01078 {
01079   G4complex z        = G4complex(1,fZommerfeld); 
01080   // G4complex gammalog = GammaLogarithm(z);
01081   G4complex gammalog = GammaLogB2n(z);
01082   fCoulombPhase0     = gammalog.imag();
01083 }

G4double G4NuclNuclDiffuseElastic::CalculateNuclearRad ( G4double A ) [inline]

Definition at line 594 of file G4NuclNuclDiffuseElastic.hh.

Referenced by GetDiffuseElasticSumXsc(), GetDiffuseElasticXsc(), Initialise(), InitialiseOnFly(), InitParameters(), InitParametersGla(), IntegralElasticProb(), SampleThetaCMS(), and TestAngleTable().

00595 {
00596   G4double r0 = 1.*CLHEP::fermi, radius;
00597   // r0 *= 1.12;
00598   // r0 *= 1.44;
00599   r0 *= fNuclearRadiusCof;
00600 
00601   /*
00602   if( A < 50. )
00603   {
00604     if( A > 10. ) r0  = 1.16*( 1 - std::pow(A, -2./3.) )*CLHEP::fermi;   // 1.08*fermi;
00605     else          r0  = 1.1*CLHEP::fermi;
00606 
00607     radius = r0*std::pow(A, 1./3.);
00608   }
00609   else
00610   {
00611     r0 = 1.7*CLHEP::fermi;   // 1.7*fermi;
00612 
00613     radius = r0*std::pow(A, 0.27); // 0.27);
00614   }
00615   */
00616   radius = r0*std::pow(A, 1./3.);
00617 
00618   return radius;
00619 }

G4double G4NuclNuclDiffuseElastic::CalculateParticleBeta	(	const G4ParticleDefinition *	particle,
		G4double	momentum
	)			`[inline]`

Definition at line 554 of file G4NuclNuclDiffuseElastic.hh.

References G4ParticleDefinition::GetPDGMass().

Referenced by GetCoulombElasticXsc(), GetCoulombIntegralXsc(), GetCoulombTotalXsc(), and GetDiffuseElasticSumXsc().

00556 {
00557   G4double mass = particle->GetPDGMass();
00558   G4double a    = momentum/mass;
00559   fBeta         = a/std::sqrt(1+a*a);
00560 
00561   return fBeta; 
00562 }

void G4NuclNuclDiffuseElastic::CalculateRutherfordAnglePar ( ) [inline]

Definition at line 1104 of file G4NuclNuclDiffuseElastic.hh.

Referenced by InitDynParameters(), and InitParameters().

01105 {
01106   fHalfRutThetaTg   = fZommerfeld/fProfileLambda;  // (fWaveVector*fNuclearRadius);
01107   fRutherfordTheta  = 2.*std::atan(fHalfRutThetaTg);
01108   fHalfRutThetaTg2  = fHalfRutThetaTg*fHalfRutThetaTg;
01109   // G4cout<<"fRutherfordTheta = "<<fRutherfordTheta/degree<<" degree"<<G4endl;
01110 
01111 }

G4double G4NuclNuclDiffuseElastic::CalculateZommerfeld	(	G4double	beta,
		G4double	Z1,
		G4double	Z2
	)			`[inline]`

Definition at line 568 of file G4NuclNuclDiffuseElastic.hh.

Referenced by GetCoulombElasticXsc(), GetCoulombIntegralXsc(), GetCoulombTotalXsc(), GetDiffuseElasticSumXsc(), InitDynParameters(), InitParameters(), InitParametersGla(), and TestAngleTable().

00569 {
00570   fZommerfeld = CLHEP::fine_structure_const*Z1*Z2/beta;
00571 
00572   return fZommerfeld; 
00573 }

G4complex G4NuclNuclDiffuseElastic::CoulombAmplitude ( G4double theta ) [inline]

Definition at line 1043 of file G4NuclNuclDiffuseElastic.hh.

Referenced by AmplitudeGG(), AmplitudeGla(), AmplitudeNear(), AmplitudeSim(), and CoulombAmplitudeMod2().

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 }

G4double G4NuclNuclDiffuseElastic::CoulombAmplitudeMod2 ( G4double theta ) [inline]

Definition at line 1064 of file G4NuclNuclDiffuseElastic.hh.

References CoulombAmplitude().

01065 {
01066   G4complex ca = CoulombAmplitude(theta);
01067   G4double out = ca.real()*ca.real() + ca.imag()*ca.imag();
01068 
01069   return  out;
01070 }

G4double G4NuclNuclDiffuseElastic::DampFactor ( G4double z ) [inline]

Definition at line 510 of file G4NuclNuclDiffuseElastic.hh.

Referenced by GetDiffElasticProb(), GetDiffElasticSumProb(), and GetDiffElasticSumProbA().

00511 {
00512   G4double df;
00513   G4double f2 = 2., f3 = 6., f4 = 24.; // first factorials
00514 
00515   // x *= pi;
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 }

G4complex G4NuclNuclDiffuseElastic::GammaLess ( G4double theta ) [inline]

Definition at line 1210 of file G4NuclNuclDiffuseElastic.hh.

References GetErfcInt(), and G4INCL::Math::pi.

Referenced by AmplitudeNear().

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 }

G4complex G4NuclNuclDiffuseElastic::GammaLogarithm ( G4complex xx ) [inline]

Definition at line 722 of file G4NuclNuclDiffuseElastic.hh.

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 }

G4complex G4NuclNuclDiffuseElastic::GammaLogB2n ( G4complex xx ) [inline]

Definition at line 746 of file G4NuclNuclDiffuseElastic.hh.

Referenced by CalculateCoulombPhase(), and CalculateCoulombPhaseZero().

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 }

G4complex G4NuclNuclDiffuseElastic::GammaMore ( G4double theta ) [inline]

Definition at line 1238 of file G4NuclNuclDiffuseElastic.hh.

References GetErfcInt(), and G4INCL::Math::pi.

Referenced by AmplitudeNear().

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 }

G4double G4NuclNuclDiffuseElastic::GetCint ( G4double x ) [inline]

Definition at line 1012 of file G4NuclNuclDiffuseElastic.hh.

References GetCosHaPit2(), and G4Integrator< T, F >::Legendre96().

Referenced by GetRatioGen(), and GetRatioSim().

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 }

G4double G4NuclNuclDiffuseElastic::GetCofAlphaCoulomb ( ) [inline]

Definition at line 301 of file G4NuclNuclDiffuseElastic.hh.

00301 {return fCofAlphaCoulomb;};

G4double G4NuclNuclDiffuseElastic::GetCofAlphaMax ( ) [inline]

Definition at line 300 of file G4NuclNuclDiffuseElastic.hh.

00300 {return fCofAlphaMax;};

G4double G4NuclNuclDiffuseElastic::GetCosHaPit2 ( G4double t ) [inline]

Definition at line 198 of file G4NuclNuclDiffuseElastic.hh.

Referenced by GetCint().

00198 {return std::cos(CLHEP::halfpi*t*t);};

G4double G4NuclNuclDiffuseElastic::GetCoulombElasticXsc	(	const G4ParticleDefinition *	particle,
		G4double	theta,
		G4double	momentum,
		G4double	Z
	)			`[inline]`

Definition at line 625 of file G4NuclNuclDiffuseElastic.hh.

References CalculateAm(), CalculateParticleBeta(), CalculateZommerfeld(), G4ParticleDefinition::GetPDGCharge(), and CLHEP::detail::n.

Referenced by GetInvCoulombElasticXsc().

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 }

G4double G4NuclNuclDiffuseElastic::GetCoulombIntegralXsc	(	const G4ParticleDefinition *	particle,
		G4double	momentum,
		G4double	Z,
		G4double	theta1,
		G4double	theta2
	)			`[inline]`

Definition at line 690 of file G4NuclNuclDiffuseElastic.hh.

References CalculateAm(), CalculateParticleBeta(), CalculateZommerfeld(), G4cout, G4endl, G4ParticleDefinition::GetPDGCharge(), and CLHEP::detail::n.

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   // G4cout<<"beta = "<<beta<<G4endl;
00700   G4double z             = particle->GetPDGCharge();
00701   G4double n             = CalculateZommerfeld( beta, z, Z );
00702   // G4cout<<"fZomerfeld = "<<n<<G4endl;
00703   G4double am            = CalculateAm( momentum, n, Z);
00704   // G4cout<<"cof Am = "<<am<<G4endl;
00705   G4double k             = momentum/CLHEP::hbarc;
00706   // G4cout<<"k = "<<k*CLHEP::fermi<<" 1/fermi"<<G4endl;
00707   // G4cout<<"k*Bohr_radius = "<<k*CLHEP::Bohr_radius<<G4endl;
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 }

G4double G4NuclNuclDiffuseElastic::GetCoulombTotalXsc	(	const G4ParticleDefinition *	particle,
		G4double	momentum,
		G4double	Z
	)			`[inline]`

Definition at line 665 of file G4NuclNuclDiffuseElastic.hh.

References CalculateAm(), CalculateParticleBeta(), CalculateZommerfeld(), G4cout, G4endl, G4ParticleDefinition::GetPDGCharge(), CLHEP::detail::n, and G4INCL::Math::pi.

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 }

G4double G4NuclNuclDiffuseElastic::GetDiffElasticProb ( G4double theta )

Definition at line 389 of file G4NuclNuclDiffuseElastic.cc.

References BesselJone(), BesselJzero(), BesselOneByArg(), DampFactor(), G4InuclParticleNames::lambda, and G4INCL::Math::pi.

Referenced by GetDiffuseElasticXsc().

00394 {
00395   G4double sigma, bzero, bzero2, bonebyarg, bonebyarg2, damp, damp2;
00396   G4double delta, diffuse, gamma;
00397   G4double e1, e2, bone, bone2;
00398 
00399   // G4double wavek = momentum/hbarc;  // wave vector
00400   // G4double r0    = 1.08*fermi;
00401   // G4double rad   = r0*std::pow(A, 1./3.);
00402 
00403   G4double kr    = fWaveVector*fNuclearRadius; // wavek*rad;
00404   G4double kr2   = kr*kr;
00405   G4double krt   = kr*theta;
00406 
00407   bzero      = BesselJzero(krt);
00408   bzero2     = bzero*bzero;    
00409   bone       = BesselJone(krt);
00410   bone2      = bone*bone;
00411   bonebyarg  = BesselOneByArg(krt);
00412   bonebyarg2 = bonebyarg*bonebyarg;  
00413 
00414   if (fParticle == theProton)
00415   {
00416     diffuse = 0.63*fermi;
00417     gamma   = 0.3*fermi;
00418     delta   = 0.1*fermi*fermi;
00419     e1      = 0.3*fermi;
00420     e2      = 0.35*fermi;
00421   }
00422   else // as proton, if were not defined 
00423   {
00424     diffuse = 0.63*fermi;
00425     gamma   = 0.3*fermi;
00426     delta   = 0.1*fermi*fermi;
00427     e1      = 0.3*fermi;
00428     e2      = 0.35*fermi;
00429   }
00430   G4double lambda = 15.; // 15 ok
00431 
00432   //  G4double kgamma    = fWaveVector*gamma;   // wavek*delta;
00433 
00434   G4double kgamma    = lambda*(1.-std::exp(-fWaveVector*gamma/lambda));   // wavek*delta;
00435   G4double kgamma2   = kgamma*kgamma;
00436 
00437   // G4double dk2t  = delta*fWaveVector*fWaveVector*theta; // delta*wavek*wavek*theta;
00438   // G4double dk2t2 = dk2t*dk2t;
00439   // G4double pikdt = pi*fWaveVector*diffuse*theta;// pi*wavek*diffuse*theta;
00440 
00441   G4double pikdt    = lambda*(1.-std::exp(-pi*fWaveVector*diffuse*theta/lambda));   // wavek*delta;
00442 
00443   damp           = DampFactor(pikdt);
00444   damp2          = damp*damp;
00445 
00446   G4double mode2k2 = (e1*e1+e2*e2)*fWaveVector*fWaveVector;  
00447   G4double e2dk3t  = -2.*e2*delta*fWaveVector*fWaveVector*fWaveVector*theta;
00448 
00449 
00450   sigma  = kgamma2;
00451   // sigma  += dk2t2;
00452   sigma *= bzero2;
00453   sigma += mode2k2*bone2 + e2dk3t*bzero*bone;
00454   sigma += kr2*bonebyarg2;
00455   sigma *= damp2;          // *rad*rad;
00456 
00457   return sigma;
00458 }

G4double G4NuclNuclDiffuseElastic::GetDiffElasticSumProb ( G4double theta )

Definition at line 466 of file G4NuclNuclDiffuseElastic.cc.

References BesselJone(), BesselJzero(), BesselOneByArg(), DampFactor(), G4InuclParticleNames::lambda, and G4INCL::Math::pi.

Referenced by GetDiffuseElasticSumXsc().

00471 {
00472   G4double sigma, bzero, bzero2, bonebyarg, bonebyarg2, damp, damp2;
00473   G4double delta, diffuse, gamma;
00474   G4double e1, e2, bone, bone2;
00475 
00476   // G4double wavek = momentum/hbarc;  // wave vector
00477   // G4double r0    = 1.08*fermi;
00478   // G4double rad   = r0*std::pow(A, 1./3.);
00479 
00480   G4double kr    = fWaveVector*fNuclearRadius; // wavek*rad;
00481   G4double kr2   = kr*kr;
00482   G4double krt   = kr*theta;
00483 
00484   bzero      = BesselJzero(krt);
00485   bzero2     = bzero*bzero;    
00486   bone       = BesselJone(krt);
00487   bone2      = bone*bone;
00488   bonebyarg  = BesselOneByArg(krt);
00489   bonebyarg2 = bonebyarg*bonebyarg;  
00490 
00491   if (fParticle == theProton)
00492   {
00493     diffuse = 0.63*fermi;
00494     // diffuse = 0.6*fermi;
00495     gamma   = 0.3*fermi;
00496     delta   = 0.1*fermi*fermi;
00497     e1      = 0.3*fermi;
00498     e2      = 0.35*fermi;
00499   }
00500   else // as proton, if were not defined 
00501   {
00502     diffuse = 0.63*fermi;
00503     gamma   = 0.3*fermi;
00504     delta   = 0.1*fermi*fermi;
00505     e1      = 0.3*fermi;
00506     e2      = 0.35*fermi;
00507   }
00508   G4double lambda = 15.; // 15 ok
00509   // G4double kgamma    = fWaveVector*gamma;   // wavek*delta;
00510   G4double kgamma    = lambda*(1.-std::exp(-fWaveVector*gamma/lambda));   // wavek*delta;
00511 
00512   // G4cout<<"kgamma = "<<kgamma<<G4endl;
00513 
00514   if(fAddCoulomb)  // add Coulomb correction
00515   {
00516     G4double sinHalfTheta  = std::sin(0.5*theta);
00517     G4double sinHalfTheta2 = sinHalfTheta*sinHalfTheta;
00518 
00519     kgamma += 0.5*fZommerfeld/kr/(sinHalfTheta2+fAm); // correction at J0()
00520   // kgamma += 0.65*fZommerfeld/kr/(sinHalfTheta2+fAm); // correction at J0()
00521   }
00522 
00523   G4double kgamma2   = kgamma*kgamma;
00524 
00525   // G4double dk2t  = delta*fWaveVector*fWaveVector*theta; // delta*wavek*wavek*theta;
00526   //   G4cout<<"dk2t = "<<dk2t<<G4endl;
00527   // G4double dk2t2 = dk2t*dk2t;
00528   // G4double pikdt = pi*fWaveVector*diffuse*theta;// pi*wavek*diffuse*theta;
00529 
00530   G4double pikdt    = lambda*(1.-std::exp(-pi*fWaveVector*diffuse*theta/lambda));   // wavek*delta;
00531 
00532   // G4cout<<"pikdt = "<<pikdt<<G4endl;
00533 
00534   damp           = DampFactor(pikdt);
00535   damp2          = damp*damp;
00536 
00537   G4double mode2k2 = (e1*e1+e2*e2)*fWaveVector*fWaveVector;  
00538   G4double e2dk3t  = -2.*e2*delta*fWaveVector*fWaveVector*fWaveVector*theta;
00539 
00540   sigma  = kgamma2;
00541   // sigma += dk2t2;
00542   sigma *= bzero2;
00543   sigma += mode2k2*bone2; 
00544   sigma += e2dk3t*bzero*bone;
00545 
00546   // sigma += kr2*(1 + 8.*fZommerfeld*fZommerfeld/kr2)*bonebyarg2;  // correction at J1()/()
00547   sigma += kr2*bonebyarg2;  // correction at J1()/()
00548 
00549   sigma *= damp2;          // *rad*rad;
00550 
00551   return sigma;
00552 }

G4double G4NuclNuclDiffuseElastic::GetDiffElasticSumProbA ( G4double alpha )

Definition at line 561 of file G4NuclNuclDiffuseElastic.cc.

References BesselJone(), BesselJzero(), BesselOneByArg(), DampFactor(), G4InuclParticleNames::lambda, and G4INCL::Math::pi.

Referenced by GetIntegrandFunction().

00562 {
00563   G4double theta; 
00564 
00565   theta = std::sqrt(alpha);
00566 
00567   // theta = std::acos( 1 - alpha/2. );
00568 
00569   G4double sigma, bzero, bzero2, bonebyarg, bonebyarg2, damp, damp2;
00570   G4double delta, diffuse, gamma;
00571   G4double e1, e2, bone, bone2;
00572 
00573   // G4double wavek = momentum/hbarc;  // wave vector
00574   // G4double r0    = 1.08*fermi;
00575   // G4double rad   = r0*std::pow(A, 1./3.);
00576 
00577   G4double kr    = fWaveVector*fNuclearRadius; // wavek*rad;
00578   G4double kr2   = kr*kr;
00579   G4double krt   = kr*theta;
00580 
00581   bzero      = BesselJzero(krt);
00582   bzero2     = bzero*bzero;    
00583   bone       = BesselJone(krt);
00584   bone2      = bone*bone;
00585   bonebyarg  = BesselOneByArg(krt);
00586   bonebyarg2 = bonebyarg*bonebyarg;  
00587 
00588   if (fParticle == theProton)
00589   {
00590     diffuse = 0.63*fermi;
00591     // diffuse = 0.6*fermi;
00592     gamma   = 0.3*fermi;
00593     delta   = 0.1*fermi*fermi;
00594     e1      = 0.3*fermi;
00595     e2      = 0.35*fermi;
00596   }
00597   else // as proton, if were not defined 
00598   {
00599     diffuse = 0.63*fermi;
00600     gamma   = 0.3*fermi;
00601     delta   = 0.1*fermi*fermi;
00602     e1      = 0.3*fermi;
00603     e2      = 0.35*fermi;
00604   }
00605   G4double lambda = 15.; // 15 ok
00606   // G4double kgamma    = fWaveVector*gamma;   // wavek*delta;
00607   G4double kgamma    = lambda*(1.-std::exp(-fWaveVector*gamma/lambda));   // wavek*delta;
00608 
00609   // G4cout<<"kgamma = "<<kgamma<<G4endl;
00610 
00611   if(fAddCoulomb)  // add Coulomb correction
00612   {
00613     G4double sinHalfTheta  = theta*0.5; // std::sin(0.5*theta);
00614     G4double sinHalfTheta2 = sinHalfTheta*sinHalfTheta;
00615 
00616     kgamma += 0.5*fZommerfeld/kr/(sinHalfTheta2+fAm); // correction at J0()
00617   // kgamma += 0.65*fZommerfeld/kr/(sinHalfTheta2+fAm); // correction at J0()
00618   }
00619 
00620   G4double kgamma2   = kgamma*kgamma;
00621 
00622   // G4double dk2t  = delta*fWaveVector*fWaveVector*theta; // delta*wavek*wavek*theta;
00623   //   G4cout<<"dk2t = "<<dk2t<<G4endl;
00624   // G4double dk2t2 = dk2t*dk2t;
00625   // G4double pikdt = pi*fWaveVector*diffuse*theta;// pi*wavek*diffuse*theta;
00626 
00627   G4double pikdt    = lambda*(1.-std::exp(-pi*fWaveVector*diffuse*theta/lambda));   // wavek*delta;
00628 
00629   // G4cout<<"pikdt = "<<pikdt<<G4endl;
00630 
00631   damp           = DampFactor(pikdt);
00632   damp2          = damp*damp;
00633 
00634   G4double mode2k2 = (e1*e1+e2*e2)*fWaveVector*fWaveVector;  
00635   G4double e2dk3t  = -2.*e2*delta*fWaveVector*fWaveVector*fWaveVector*theta;
00636 
00637   sigma  = kgamma2;
00638   // sigma += dk2t2;
00639   sigma *= bzero2;
00640   sigma += mode2k2*bone2; 
00641   sigma += e2dk3t*bzero*bone;
00642 
00643   // sigma += kr2*(1 + 8.*fZommerfeld*fZommerfeld/kr2)*bonebyarg2;  // correction at J1()/()
00644   sigma += kr2*bonebyarg2;  // correction at J1()/()
00645 
00646   sigma *= damp2;          // *rad*rad;
00647 
00648   return sigma;
00649 }

G4double G4NuclNuclDiffuseElastic::GetDiffuseElasticSumXsc	(	const G4ParticleDefinition *	particle,
		G4double	theta,
		G4double	momentum,
		G4double	A,
		G4double	Z
	)

Definition at line 253 of file G4NuclNuclDiffuseElastic.cc.

References CalculateAm(), CalculateNuclearRad(), CalculateParticleBeta(), CalculateZommerfeld(), GetDiffElasticSumProb(), and G4ParticleDefinition::GetPDGCharge().

Referenced by GetInvElasticSumXsc().

00257 {
00258   fParticle      = particle;
00259   fWaveVector    = momentum/hbarc;
00260   fAtomicWeight  = A;
00261   fAtomicNumber  = Z;
00262   fNuclearRadius = CalculateNuclearRad(A);
00263   fAddCoulomb    = false;
00264 
00265   G4double z     = particle->GetPDGCharge();
00266 
00267   G4double kRt   = fWaveVector*fNuclearRadius*theta;
00268   G4double kRtC  = 1.9;
00269 
00270   if( z && (kRt > kRtC) )
00271   {
00272     fAddCoulomb = true;
00273     fBeta       = CalculateParticleBeta( particle, momentum);
00274     fZommerfeld = CalculateZommerfeld( fBeta, z, fAtomicNumber);
00275     fAm         = CalculateAm( momentum, fZommerfeld, fAtomicNumber);
00276   }
00277   G4double sigma = fNuclearRadius*fNuclearRadius*GetDiffElasticSumProb(theta);
00278 
00279   return sigma;
00280 }

G4double G4NuclNuclDiffuseElastic::GetDiffuseElasticXsc	(	const G4ParticleDefinition *	particle,
		G4double	theta,
		G4double	momentum,
		G4double	A
	)

Definition at line 182 of file G4NuclNuclDiffuseElastic.cc.

References CalculateNuclearRad(), and GetDiffElasticProb().

Referenced by GetInvElasticXsc().

00186 {
00187   fParticle      = particle;
00188   fWaveVector    = momentum/hbarc;
00189   fAtomicWeight  = A;
00190   fAddCoulomb    = false;
00191   fNuclearRadius = CalculateNuclearRad(A);
00192 
00193   G4double sigma = fNuclearRadius*fNuclearRadius*GetDiffElasticProb(theta);
00194 
00195   return sigma;
00196 }

G4double G4NuclNuclDiffuseElastic::GetErf ( G4double x ) [inline]

Definition at line 767 of file G4NuclNuclDiffuseElastic.hh.

Referenced by GetErfComp(), and GetErfInt().

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 }

G4complex G4NuclNuclDiffuseElastic::GetErfcComp	(	G4complex	z,
		G4int	nMax
	)			`[inline]`

Definition at line 788 of file G4NuclNuclDiffuseElastic.hh.

References GetErfComp().

00789 {
00790   G4complex erfcz = 1. - GetErfComp( z, nMax);
00791   return erfcz;
00792 }

G4complex G4NuclNuclDiffuseElastic::GetErfcInt ( G4complex z ) [inline]

Definition at line 808 of file G4NuclNuclDiffuseElastic.hh.

References GetErfInt().

Referenced by AmplitudeSim(), GammaLess(), and GammaMore().

00809 {
00810   G4complex erfcz = 1. - GetErfInt( z); // , nMax);
00811   return erfcz;
00812 }

G4complex G4NuclNuclDiffuseElastic::GetErfComp	(	G4complex	z,
		G4int	nMax
	)			`[inline]`

Definition at line 878 of file G4NuclNuclDiffuseElastic.hh.

References GetErf(), CLHEP::detail::n, and G4INCL::Math::pi.

Referenced by GetErfcComp(), and TestErfcComp().

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);  // /(n2+0.5*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 }

G4complex G4NuclNuclDiffuseElastic::GetErfcSer	(	G4complex	z,
		G4int	nMax
	)			`[inline]`

Definition at line 798 of file G4NuclNuclDiffuseElastic.hh.

References GetErfSer().

00799 {
00800   G4complex erfcz = 1. - GetErfSer( z, nMax);
00801   return erfcz;
00802 }

G4complex G4NuclNuclDiffuseElastic::GetErfInt ( G4complex z ) [inline]

Definition at line 987 of file G4NuclNuclDiffuseElastic.hh.

References GetErf(), GetExpCos(), GetExpSin(), G4Integrator< T, F >::Legendre96(), and G4INCL::Math::pi.

Referenced by GetErfcInt(), and TestErfcInt().

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 }

G4complex G4NuclNuclDiffuseElastic::GetErfSer	(	G4complex	z,
		G4int	nMax
	)			`[inline]`

Definition at line 938 of file G4NuclNuclDiffuseElastic.hh.

References CLHEP::detail::n, and G4INCL::Math::pi.

Referenced by GetErfcSer(), and TestErfcSer().

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 }

G4double G4NuclNuclDiffuseElastic::GetExpCos ( G4double x ) [inline]

Definition at line 961 of file G4NuclNuclDiffuseElastic.hh.

Referenced by GetErfInt().

00962 {
00963   G4double result;
00964 
00965   result = std::exp(x*x-fReZ*fReZ);
00966   result *= std::cos(2.*x*fReZ);
00967   return result;
00968 }

G4double G4NuclNuclDiffuseElastic::GetExpSin ( G4double x ) [inline]

Definition at line 972 of file G4NuclNuclDiffuseElastic.hh.

Referenced by GetErfInt().

00973 {
00974   G4double result;
00975 
00976   result = std::exp(x*x-fReZ*fReZ);
00977   result *= std::sin(2.*x*fReZ);
00978   return result;
00979 }

G4double G4NuclNuclDiffuseElastic::GetFresnelDiffuseXsc ( G4double theta ) [inline]

Definition at line 1415 of file G4NuclNuclDiffuseElastic.hh.

References GetRatioGen(), and GetRutherfordXsc().

Referenced by GetFresnelIntegrandXsc().

01416 {
01417   G4double ratio   = GetRatioGen(theta);
01418   G4double ruthXsc = GetRutherfordXsc(theta);
01419   G4double xsc     = ratio*ruthXsc;
01420   return xsc;
01421 }

G4double G4NuclNuclDiffuseElastic::GetFresnelIntegrandXsc ( G4double alpha ) [inline]

Definition at line 1427 of file G4NuclNuclDiffuseElastic.hh.

References GetFresnelDiffuseXsc().

Referenced by BuildAngleTable().

01428 {
01429   G4double theta = std::sqrt(alpha);
01430   G4double xsc     = GetFresnelDiffuseXsc(theta);
01431   return xsc;
01432 }

G4double G4NuclNuclDiffuseElastic::GetHadronNucleonXscNS	(	G4ParticleDefinition *	pParticle,
		G4double	pTkin,
		G4ParticleDefinition *	tParticle
	)			`[inline]`

Definition at line 1671 of file G4NuclNuclDiffuseElastic.hh.

References CalcMandelstamS(), G4cout, G4endl, and G4ParticleDefinition::GetPDGMass().

Referenced by InitParametersGla().

01674 {
01675   G4double xsection(0), /*Delta,*/ 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;  // in GeV for parametrisation
01689   proj_momentum /= CLHEP::GeV;
01690   proj_energy   /= CLHEP::GeV;
01691   proj_mass     /= CLHEP::GeV;
01692   G4double logS = std::log(sMand);
01693 
01694   // General PDG fit constants
01695 
01696 
01697   // fEtaRatio=Re[f(0)]/Im[f(0)]
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   // xsc
01715   
01716   if( proj_momentum >= 10. ) // high energy: pp = nn = np
01717     // if( proj_momentum >= 2.)
01718   {
01719     //Delta = 1.;
01720 
01721     //if( proj_energy < 40. ) Delta = 0.916+0.0021*proj_energy;
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);        //  mb
01731     }
01732   }
01733   else // low energy pp = nn != np
01734   {
01735       if(pParticle == tParticle) // pp or nn      // nn to be pp
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  // if( proj_momentum < 10.  )
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  // pn to be np
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    // if( proj_momentum < 10.  )
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; // parametrised in mb
01773   G4cout<<"xsection = "<<xsection/CLHEP::millibarn<<" mb"<<G4endl;
01774   return xsection;
01775 }

G4double G4NuclNuclDiffuseElastic::GetIntegrandFunction ( G4double theta )

Definition at line 657 of file G4NuclNuclDiffuseElastic.cc.

References GetDiffElasticSumProbA().

Referenced by IntegralElasticProb(), SampleThetaCMS(), and TestAngleTable().

00658 {
00659   G4double result;
00660 
00661   result  = GetDiffElasticSumProbA(alpha);
00662 
00663   // result *= 2*pi*std::sin(theta);
00664 
00665   return result;
00666 }

G4double G4NuclNuclDiffuseElastic::GetInvCoulombElasticXsc	(	const G4ParticleDefinition *	particle,
		G4double	tMand,
		G4double	momentum,
		G4double	A,
		G4double	Z
	)

Definition at line 340 of file G4NuclNuclDiffuseElastic.cc.

References G4ParticleTable::FindIon(), GetCoulombElasticXsc(), G4ParticleTable::GetParticleTable(), G4ParticleDefinition::GetPDGMass(), G4He3::He3(), G4INCL::Math::pi, and G4Triton::Triton().

00344 {
00345   G4double m1 = particle->GetPDGMass();
00346   G4LorentzVector lv1(0.,0.,plab,std::sqrt(plab*plab+m1*m1));
00347 
00348   G4int iZ = static_cast<G4int>(Z+0.5);
00349   G4int iA = static_cast<G4int>(A+0.5);
00350   G4ParticleDefinition * theDef = 0;
00351 
00352   if      (iZ == 1 && iA == 1) theDef = theProton;
00353   else if (iZ == 1 && iA == 2) theDef = theDeuteron;
00354   else if (iZ == 1 && iA == 3) theDef = G4Triton::Triton();
00355   else if (iZ == 2 && iA == 3) theDef = G4He3::He3();
00356   else if (iZ == 2 && iA == 4) theDef = theAlpha;
00357   else theDef = G4ParticleTable::GetParticleTable()->FindIon(iZ,iA,0,iZ);
00358  
00359   G4double tmass = theDef->GetPDGMass();
00360 
00361   G4LorentzVector lv(0.0,0.0,0.0,tmass);   
00362   lv += lv1;
00363 
00364   G4ThreeVector bst = lv.boostVector();
00365   lv1.boost(-bst);
00366 
00367   G4ThreeVector p1 = lv1.vect();
00368   G4double ptot    = p1.mag();
00369   G4double ptot2 = ptot*ptot;
00370   G4double cost = 1 - 0.5*std::fabs(tMand)/ptot2;
00371 
00372   if( cost >= 1.0 )      cost = 1.0;  
00373   else if( cost <= -1.0) cost = -1.0;
00374   
00375   G4double thetaCMS = std::acos(cost);
00376 
00377   G4double sigma = GetCoulombElasticXsc( particle, thetaCMS, ptot, Z );
00378 
00379   sigma *= pi/ptot2;
00380 
00381   return sigma;
00382 }

G4double G4NuclNuclDiffuseElastic::GetInvElasticSumXsc	(	const G4ParticleDefinition *	particle,
		G4double	tMand,
		G4double	momentum,
		G4double	A,
		G4double	Z
	)

Definition at line 288 of file G4NuclNuclDiffuseElastic.cc.

References G4ParticleTable::FindIon(), GetDiffuseElasticSumXsc(), G4ParticleTable::GetParticleTable(), G4ParticleDefinition::GetPDGMass(), G4He3::He3(), G4INCL::Math::pi, and G4Triton::Triton().

00292 {
00293   G4double m1 = particle->GetPDGMass();
00294 
00295   G4LorentzVector lv1(0.,0.,plab,std::sqrt(plab*plab+m1*m1));
00296 
00297   G4int iZ = static_cast<G4int>(Z+0.5);
00298   G4int iA = static_cast<G4int>(A+0.5);
00299 
00300   G4ParticleDefinition* theDef = 0;
00301 
00302   if      (iZ == 1 && iA == 1) theDef = theProton;
00303   else if (iZ == 1 && iA == 2) theDef = theDeuteron;
00304   else if (iZ == 1 && iA == 3) theDef = G4Triton::Triton();
00305   else if (iZ == 2 && iA == 3) theDef = G4He3::He3();
00306   else if (iZ == 2 && iA == 4) theDef = theAlpha;
00307   else theDef = G4ParticleTable::GetParticleTable()->FindIon(iZ,iA,0,iZ);
00308  
00309   G4double tmass = theDef->GetPDGMass();
00310 
00311   G4LorentzVector lv(0.0,0.0,0.0,tmass);   
00312   lv += lv1;
00313 
00314   G4ThreeVector bst = lv.boostVector();
00315   lv1.boost(-bst);
00316 
00317   G4ThreeVector p1 = lv1.vect();
00318   G4double ptot    = p1.mag();
00319   G4double ptot2   = ptot*ptot;
00320   G4double cost    = 1 - 0.5*std::fabs(tMand)/ptot2;
00321 
00322   if( cost >= 1.0 )      cost = 1.0;  
00323   else if( cost <= -1.0) cost = -1.0;
00324   
00325   G4double thetaCMS = std::acos(cost);
00326 
00327   G4double sigma = GetDiffuseElasticSumXsc( particle, thetaCMS, ptot, A, Z );
00328 
00329   sigma *= pi/ptot2;
00330 
00331   return sigma;
00332 }

G4double G4NuclNuclDiffuseElastic::GetInvElasticXsc	(	const G4ParticleDefinition *	particle,
		G4double	theta,
		G4double	momentum,
		G4double	A,
		G4double	Z
	)

Definition at line 203 of file G4NuclNuclDiffuseElastic.cc.

References G4ParticleTable::FindIon(), GetDiffuseElasticXsc(), G4ParticleTable::GetParticleTable(), G4ParticleDefinition::GetPDGMass(), G4He3::He3(), G4INCL::Math::pi, and G4Triton::Triton().

00207 {
00208   G4double m1 = particle->GetPDGMass();
00209   G4LorentzVector lv1(0.,0.,plab,std::sqrt(plab*plab+m1*m1));
00210 
00211   G4int iZ = static_cast<G4int>(Z+0.5);
00212   G4int iA = static_cast<G4int>(A+0.5);
00213   G4ParticleDefinition * theDef = 0;
00214 
00215   if      (iZ == 1 && iA == 1) theDef = theProton;
00216   else if (iZ == 1 && iA == 2) theDef = theDeuteron;
00217   else if (iZ == 1 && iA == 3) theDef = G4Triton::Triton();
00218   else if (iZ == 2 && iA == 3) theDef = G4He3::He3();
00219   else if (iZ == 2 && iA == 4) theDef = theAlpha;
00220   else theDef = G4ParticleTable::GetParticleTable()->FindIon(iZ,iA,0,iZ);
00221  
00222   G4double tmass = theDef->GetPDGMass();
00223 
00224   G4LorentzVector lv(0.0,0.0,0.0,tmass);   
00225   lv += lv1;
00226 
00227   G4ThreeVector bst = lv.boostVector();
00228   lv1.boost(-bst);
00229 
00230   G4ThreeVector p1 = lv1.vect();
00231   G4double ptot    = p1.mag();
00232   G4double ptot2 = ptot*ptot;
00233   G4double cost = 1 - 0.5*std::fabs(tMand)/ptot2;
00234 
00235   if( cost >= 1.0 )      cost = 1.0;  
00236   else if( cost <= -1.0) cost = -1.0;
00237   
00238   G4double thetaCMS = std::acos(cost);
00239 
00240   G4double sigma = GetDiffuseElasticXsc( particle, thetaCMS, ptot, A);
00241 
00242   sigma *= pi/ptot2;
00243 
00244   return sigma;
00245 }

G4double G4NuclNuclDiffuseElastic::GetLegendrePol	(	G4int	n,
		G4double	x
	)			`[inline]`

Definition at line 814 of file G4NuclNuclDiffuseElastic.hh.

References G4INCL::Math::pi.

Referenced by AmplitudeGla().

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     // legPol = ( (2*n-1)*x*GetLegendrePol(n-1,x) - (n-1)*GetLegendrePol(n-2,x) )/n;
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 }

G4double G4NuclNuclDiffuseElastic::GetNuclearRadius ( ) [inline]

Definition at line 188 of file G4NuclNuclDiffuseElastic.hh.

00188 {return fNuclearRadius;};

G4double G4NuclNuclDiffuseElastic::GetProfileLambda ( ) [inline]

Definition at line 282 of file G4NuclNuclDiffuseElastic.hh.

00282 {return fProfileLambda;};

G4double G4NuclNuclDiffuseElastic::GetRatioGen ( G4double theta ) [inline]

Definition at line 1379 of file G4NuclNuclDiffuseElastic.hh.

References GetCint(), GetSint(), G4INCL::Math::pi, and Profile().

Referenced by GetFresnelDiffuseXsc().

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   // G4double profmod    = std::abs(prof);
01388   G4double order      = std::sqrt(fProfileLambda/sinThetaR/CLHEP::pi)*2.*sindTheta;
01389 
01390   order = std::abs(order); // since sin changes sign!
01391   // G4cout<<"order = "<<order<<G4endl; 
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 }

G4double G4NuclNuclDiffuseElastic::GetRatioSim ( G4double theta ) [inline]

Definition at line 1359 of file G4NuclNuclDiffuseElastic.hh.

References GetCint(), GetSint(), and G4INCL::Math::pi.

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   // G4cout<<"order = "<<order<<G4endl;  
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 }

G4double G4NuclNuclDiffuseElastic::GetRutherfordXsc ( G4double theta ) [inline]

Definition at line 648 of file G4NuclNuclDiffuseElastic.hh.

Referenced by GetFresnelDiffuseXsc().

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 }

G4double G4NuclNuclDiffuseElastic::GetScatteringAngle	(	G4int	iMomentum,
		G4int	iAngle,
		G4double	position
	)

Definition at line 1035 of file G4NuclNuclDiffuseElastic.cc.

References G4UniformRand.

Referenced by SampleTableThetaCMS().

01036 {
01037  G4double x1, x2, y1, y2, randAngle;
01038 
01039   if( iAngle == 0 )
01040   {
01041     randAngle = (*fAngleTable)(iMomentum)->GetLowEdgeEnergy(iAngle);
01042     // iAngle++;
01043   }
01044   else
01045   {
01046     if ( iAngle >= G4int((*fAngleTable)(iMomentum)->GetVectorLength()) )
01047     {
01048       iAngle = (*fAngleTable)(iMomentum)->GetVectorLength() - 1;
01049     }
01050     y1 = (*(*fAngleTable)(iMomentum))(iAngle-1);
01051     y2 = (*(*fAngleTable)(iMomentum))(iAngle);
01052 
01053     x1 = (*fAngleTable)(iMomentum)->GetLowEdgeEnergy(iAngle-1);
01054     x2 = (*fAngleTable)(iMomentum)->GetLowEdgeEnergy(iAngle);
01055 
01056     if ( x1 == x2 )   randAngle = x2;
01057     else
01058     {
01059       if ( y1 == y2 ) randAngle = x1 + ( x2 - x1 )*G4UniformRand();
01060       else
01061       {
01062         randAngle = x1 + ( position - y1 )*( x2 - x1 )/( y2 - y1 );
01063       }
01064     }
01065   }
01066   return randAngle;
01067 }

G4double G4NuclNuclDiffuseElastic::GetSinHaPit2 ( G4double t ) [inline]

Definition at line 199 of file G4NuclNuclDiffuseElastic.hh.

Referenced by GetSint().

00199 {return std::sin(CLHEP::halfpi*t*t);};

G4double G4NuclNuclDiffuseElastic::GetSint ( G4double x ) [inline]

Definition at line 1027 of file G4NuclNuclDiffuseElastic.hh.

References GetSinHaPit2(), and G4Integrator< T, F >::Legendre96().

Referenced by GetRatioGen(), and GetRatioSim().

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 }

void G4NuclNuclDiffuseElastic::InitDynParameters	(	const G4ParticleDefinition *	theParticle,
		G4double	partMom
	)			`[inline]`

Definition at line 1576 of file G4NuclNuclDiffuseElastic.hh.

References CalculateAm(), CalculateCoulombPhaseZero(), CalculateRutherfordAnglePar(), CalculateZommerfeld(), G4ParticleDefinition::GetPDGCharge(), G4ParticleDefinition::GetPDGMass(), and G4InuclParticleNames::lambda.

Referenced by BuildAngleTable().

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; // beta*gamma for 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; // *std::sqrt(1.-2*fZommerfeld/lambda);
01596   fProfileDelta  = fCofDelta*fProfileLambda;
01597   fProfileAlpha  = fCofAlpha*fProfileLambda;
01598 
01599   CalculateCoulombPhaseZero();
01600   CalculateRutherfordAnglePar();
01601 
01602   return;
01603 }

void G4NuclNuclDiffuseElastic::Initialise ( )

Definition at line 142 of file G4NuclNuclDiffuseElastic.cc.

References BuildAngleTable(), CalculateNuclearRad(), G4cout, G4endl, G4ParticleDefinition::GetBaryonNumber(), G4Element::GetElementTable(), G4Element::GetNumberOfElements(), and G4HadronicInteraction::verboseLevel.

00143 {
00144 
00145   // fEnergyVector = new G4PhysicsLogVector( theMinEnergy, theMaxEnergy, fEnergyBin );
00146 
00147   const G4ElementTable* theElementTable = G4Element::GetElementTable();
00148   size_t jEl, numOfEl = G4Element::GetNumberOfElements();
00149 
00150   // projectile radius
00151 
00152   G4double A1 = G4double( fParticle->GetBaryonNumber() );
00153   G4double R1 = CalculateNuclearRad(A1);
00154 
00155   for(jEl = 0 ; jEl < numOfEl; ++jEl) // application element loop
00156   {
00157     fAtomicNumber = (*theElementTable)[jEl]->GetZ();     // atomic number
00158     fAtomicWeight = (*theElementTable)[jEl]->GetN();     // number of nucleons
00159 
00160     fNuclearRadius = CalculateNuclearRad(fAtomicWeight);
00161     fNuclearRadius += R1;
00162 
00163     if(verboseLevel > 0) 
00164     {   
00165       G4cout<<"G4NuclNuclDiffuseElastic::Initialise() the element: "
00166             <<(*theElementTable)[jEl]->GetName()<<G4endl;
00167     }
00168     fElementNumberVector.push_back(fAtomicNumber);
00169     fElementNameVector.push_back((*theElementTable)[jEl]->GetName());
00170 
00171     BuildAngleTable();
00172     fAngleBank.push_back(fAngleTable);
00173   }  
00174 }

void G4NuclNuclDiffuseElastic::InitialiseOnFly	(	G4double	Z,
		G4double	A
	)

Definition at line 929 of file G4NuclNuclDiffuseElastic.cc.

References BuildAngleTable(), CalculateNuclearRad(), G4cout, G4endl, G4ParticleDefinition::GetBaryonNumber(), and G4HadronicInteraction::verboseLevel.

Referenced by SampleTableThetaCMS().

00930 {
00931   fAtomicNumber  = Z;     // atomic number
00932   fAtomicWeight  = A;     // number of nucleons
00933 
00934   G4double A1 = G4double( fParticle->GetBaryonNumber() );
00935   G4double R1 = CalculateNuclearRad(A1);
00936 
00937   fNuclearRadius = CalculateNuclearRad(fAtomicWeight) + R1;
00938   
00939   if( verboseLevel > 0 )    
00940   {
00941     G4cout<<"G4NuclNuclDiffuseElastic::Initialise() the element with Z = "
00942           <<Z<<"; and A = "<<A<<G4endl;
00943   }
00944   fElementNumberVector.push_back(fAtomicNumber);
00945 
00946   BuildAngleTable();
00947 
00948   fAngleBank.push_back(fAngleTable);
00949 
00950   return;
00951 }

void G4NuclNuclDiffuseElastic::InitParameters	(	const G4ParticleDefinition *	theParticle,
		G4double	partMom,
		G4double	Z,
		G4double	A
	)			`[inline]`

Definition at line 1531 of file G4NuclNuclDiffuseElastic.hh.

References CalculateAm(), CalculateCoulombPhaseZero(), CalculateNuclearRad(), CalculateRutherfordAnglePar(), CalculateZommerfeld(), G4cout, G4endl, G4ParticleDefinition::GetBaryonNumber(), G4ParticleDefinition::GetPDGCharge(), G4ParticleDefinition::GetPDGMass(), and G4InuclParticleNames::lambda.

01533 {
01534   fAtomicNumber  = Z;     // atomic number
01535   fAtomicWeight  = A;     // number of nucleons
01536 
01537   fNuclearRadius2 = CalculateNuclearRad(fAtomicWeight);
01538   G4double A1     = G4double( theParticle->GetBaryonNumber() );   
01539   fNuclearRadius1 = CalculateNuclearRad(A1);
01540   // fNuclearRadius = std::sqrt(fNuclearRadius1*fNuclearRadius1+fNuclearRadius2*fNuclearRadius2);
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; // beta*gamma for 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; // *std::sqrt(1.-2*fZommerfeld/lambda);
01562   G4cout<<"fProfileLambda = "<<fProfileLambda<<G4endl;
01563   fProfileDelta  = fCofDelta*fProfileLambda;
01564   fProfileAlpha  = fCofAlpha*fProfileLambda;
01565 
01566   CalculateCoulombPhaseZero();
01567   CalculateRutherfordAnglePar();
01568 
01569   return;
01570 }

void G4NuclNuclDiffuseElastic::InitParametersGla	(	const G4DynamicParticle *	aParticle,
		G4double	partMom,
		G4double	Z,
		G4double	A
	)			`[inline]`

Definition at line 1611 of file G4NuclNuclDiffuseElastic.hh.

References CalculateAm(), CalculateCoulombPhaseZero(), CalculateNuclearRad(), CalculateZommerfeld(), G4cout, G4endl, G4ParticleDefinition::GetBaryonNumber(), G4DynamicParticle::GetDefinition(), GetHadronNucleonXscNS(), G4DynamicParticle::GetKineticEnergy(), G4ParticleDefinition::GetPDGCharge(), G4ParticleDefinition::GetPDGMass(), and G4INCL::Math::pi.

01613 {
01614   fAtomicNumber  = Z;     // target atomic number
01615   fAtomicWeight  = A;     // target number of nucleons
01616 
01617   fNuclearRadius2 = CalculateNuclearRad(fAtomicWeight); // target nucleus radius
01618   G4double A1     = G4double( aParticle->GetDefinition()->GetBaryonNumber() );   
01619   fNuclearRadius1 = CalculateNuclearRad(A1); // projectile nucleus radius
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; // beta*gamma for 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 }

G4double G4NuclNuclDiffuseElastic::IntegralElasticProb	(	const G4ParticleDefinition *	particle,
		G4double	theta,
		G4double	momentum,
		G4double	A
	)

Definition at line 673 of file G4NuclNuclDiffuseElastic.cc.

References CalculateNuclearRad(), GetIntegrandFunction(), and G4Integrator< T, F >::Legendre96().

00677 {
00678   G4double result;
00679   fParticle      = particle;
00680   fWaveVector    = momentum/hbarc;
00681   fAtomicWeight  = A;
00682 
00683   fNuclearRadius = CalculateNuclearRad(A);
00684 
00685 
00686   G4Integrator<G4NuclNuclDiffuseElastic,G4double(G4NuclNuclDiffuseElastic::*)(G4double)> integral;
00687 
00688   // result = integral.Legendre10(this,&G4NuclNuclDiffuseElastic::GetIntegrandFunction, 0., theta ); 
00689   result = integral.Legendre96(this,&G4NuclNuclDiffuseElastic::GetIntegrandFunction, 0., theta ); 
00690 
00691   return result;
00692 }

G4complex G4NuclNuclDiffuseElastic::PhaseFar ( G4double theta ) [inline]

Definition at line 1191 of file G4NuclNuclDiffuseElastic.hh.

References G4INCL::Math::pi.

Referenced by AmplitudeFar().

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 }

G4complex G4NuclNuclDiffuseElastic::PhaseNear ( G4double theta ) [inline]

Definition at line 1173 of file G4NuclNuclDiffuseElastic.hh.

References G4INCL::Math::pi.

Referenced by AmplitudeNear().

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 }

G4double G4NuclNuclDiffuseElastic::Profile ( G4double theta ) [inline]

Definition at line 1154 of file G4NuclNuclDiffuseElastic.hh.

References G4INCL::Math::pi.

Referenced by GetRatioGen().

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 }

G4double G4NuclNuclDiffuseElastic::ProfileFar ( G4double theta ) [inline]

Definition at line 1138 of file G4NuclNuclDiffuseElastic.hh.

References G4INCL::Math::pi.

Referenced by AmplitudeFar().

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 }

G4double G4NuclNuclDiffuseElastic::ProfileNear ( G4double theta ) [inline]

Definition at line 1117 of file G4NuclNuclDiffuseElastic.hh.

References G4INCL::Math::pi.

Referenced by AmplitudeNear(), and AmplitudeSim().

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 }

G4double G4NuclNuclDiffuseElastic::SampleInvariantT	(	const G4ParticleDefinition *	p,
		G4double	plab,
		G4int	Z,
		G4int	A
	)			`[virtual]`

Reimplemented from G4HadronElastic.

Definition at line 766 of file G4NuclNuclDiffuseElastic.cc.

References G4NucleiProperties::GetNuclearMass(), G4ParticleDefinition::GetPDGMass(), and SampleTableT().

00768 {
00769   fParticle = aParticle;
00770   G4double m1 = fParticle->GetPDGMass();
00771   G4double totElab = std::sqrt(m1*m1+p*p);
00772   G4double mass2 = G4NucleiProperties::GetNuclearMass(A, Z);
00773   G4LorentzVector lv1(p,0.0,0.0,totElab);
00774   G4LorentzVector  lv(0.0,0.0,0.0,mass2);   
00775   lv += lv1;
00776 
00777   G4ThreeVector bst = lv.boostVector();
00778   lv1.boost(-bst);
00779 
00780   G4ThreeVector p1 = lv1.vect();
00781   G4double momentumCMS = p1.mag();
00782 
00783   G4double t = SampleTableT( aParticle,  momentumCMS, G4double(Z), G4double(A) ); // sample theta2 in cms
00784   return t;
00785 }

G4double G4NuclNuclDiffuseElastic::SampleT	(	const G4ParticleDefinition *	aParticle,
		G4double	p,
		G4double	A
	)

Definition at line 698 of file G4NuclNuclDiffuseElastic.cc.

References SampleThetaCMS().

Referenced by SampleThetaLab().

00700 {
00701   G4double theta = SampleThetaCMS( aParticle,  p, A); // sample theta in cms
00702   G4double t     = 2*p*p*( 1 - std::cos(theta) ); // -t !!!
00703   return t;
00704 }

G4double G4NuclNuclDiffuseElastic::SampleTableT	(	const G4ParticleDefinition *	aParticle,
		G4double	p,
		G4double	Z,
		G4double	A
	)

Definition at line 791 of file G4NuclNuclDiffuseElastic.cc.

References G4InuclParticleNames::alpha, and SampleTableThetaCMS().

Referenced by SampleInvariantT().

00793 {
00794   G4double alpha = SampleTableThetaCMS( aParticle,  p, Z, A); // sample theta2 in cms
00795   // G4double t     = 2*p*p*( 1 - std::cos(std::sqrt(alpha)) );             // -t !!!
00796   G4double t     = p*p*alpha;             // -t !!!
00797   return t;
00798 }

G4double G4NuclNuclDiffuseElastic::SampleTableThetaCMS	(	const G4ParticleDefinition *	aParticle,
		G4double	p,
		G4double	Z,
		G4double	A
	)

Definition at line 806 of file G4NuclNuclDiffuseElastic.cc.

References G4UniformRand, G4PhysicsVector::GetLowEdgeEnergy(), G4ParticleDefinition::GetPDGMass(), GetScatteringAngle(), InitialiseOnFly(), and position.

Referenced by SampleTableT().

00808 {
00809   size_t iElement;
00810   G4int iMomentum, iAngle;  
00811   G4double randAngle, position, theta1, theta2, E1, E2, W1, W2, W;  
00812   G4double m1 = particle->GetPDGMass();
00813 
00814   for(iElement = 0; iElement < fElementNumberVector.size(); iElement++)
00815   {
00816     if( std::fabs(Z - fElementNumberVector[iElement]) < 0.5) break;
00817   }
00818   if ( iElement == fElementNumberVector.size() ) 
00819   {
00820     InitialiseOnFly(Z,A); // table preparation, if needed
00821 
00822     // iElement--;
00823 
00824     // G4cout << "G4NuclNuclDiffuseElastic: Element with atomic number " << Z
00825     // << " is not found, return zero angle" << G4endl;
00826     // return 0.; // no table for this element
00827   }
00828   // G4cout<<"iElement = "<<iElement<<G4endl;
00829 
00830   fAngleTable = fAngleBank[iElement];
00831 
00832   G4double kinE = std::sqrt(momentum*momentum + m1*m1) - m1;
00833 
00834   for( iMomentum = 0; iMomentum < fEnergyBin; iMomentum++)
00835   {
00836     // G4cout<<iMomentum<<", kinE = "<<kinE/MeV<<", vectorE = "<<fEnergyVector->GetLowEdgeEnergy(iMomentum)/MeV<<G4endl;
00837     
00838     if( kinE < fEnergyVector->GetLowEdgeEnergy(iMomentum) ) break;
00839   }
00840   // G4cout<<"iMomentum = "<<iMomentum<<", fEnergyBin -1 = "<<fEnergyBin -1<<G4endl;
00841 
00842 
00843   if ( iMomentum >= fEnergyBin ) iMomentum = fEnergyBin-1;   // kinE is more then theMaxEnergy
00844   if ( iMomentum < 0 )           iMomentum = 0; // against negative index, kinE < theMinEnergy
00845 
00846 
00847   if (iMomentum == fEnergyBin -1 || iMomentum == 0 )   // the table edges
00848   {
00849     position = (*(*fAngleTable)(iMomentum))(fAngleBin-2)*G4UniformRand();
00850 
00851     // G4cout<<"position = "<<position<<G4endl;
00852 
00853     for(iAngle = 0; iAngle < fAngleBin-1; iAngle++)
00854     {
00855       if( position < (*(*fAngleTable)(iMomentum))(iAngle) ) break;
00856     }
00857     if (iAngle >= fAngleBin-1) iAngle = fAngleBin-2;
00858 
00859     // G4cout<<"iAngle = "<<iAngle<<G4endl;
00860 
00861     randAngle = GetScatteringAngle(iMomentum, iAngle, position);
00862 
00863     // G4cout<<"randAngle = "<<randAngle<<G4endl;
00864   }
00865   else  // kinE inside between energy table edges
00866   {
00867     // position = (*(*fAngleTable)(iMomentum))(fAngleBin-2)*G4UniformRand();
00868     position = (*(*fAngleTable)(iMomentum))(0)*G4UniformRand();
00869 
00870     // G4cout<<"position = "<<position<<G4endl;
00871 
00872     for(iAngle = 0; iAngle < fAngleBin-1; iAngle++)
00873     {
00874       // if( position < (*(*fAngleTable)(iMomentum))(iAngle) ) break;
00875       if( position > (*(*fAngleTable)(iMomentum))(iAngle) ) break;
00876     }
00877     if (iAngle >= fAngleBin-1) iAngle = fAngleBin-2;
00878 
00879     // G4cout<<"iAngle = "<<iAngle<<G4endl;
00880 
00881     theta2  = GetScatteringAngle(iMomentum, iAngle, position);
00882 
00883     // G4cout<<"theta2 = "<<theta2<<G4endl;
00884 
00885     E2 = fEnergyVector->GetLowEdgeEnergy(iMomentum);
00886 
00887     // G4cout<<"E2 = "<<E2<<G4endl;
00888     
00889     iMomentum--;
00890     
00891     // position = (*(*fAngleTable)(iMomentum))(fAngleBin-2)*G4UniformRand();
00892 
00893     // G4cout<<"position = "<<position<<G4endl;
00894 
00895     for(iAngle = 0; iAngle < fAngleBin-1; iAngle++)
00896     {
00897       // if( position < (*(*fAngleTable)(iMomentum))(iAngle) ) break;
00898       if( position > (*(*fAngleTable)(iMomentum))(iAngle) ) break;
00899     }
00900     if (iAngle >= fAngleBin-1) iAngle = fAngleBin-2;
00901     
00902     theta1  = GetScatteringAngle(iMomentum, iAngle, position);
00903 
00904     // G4cout<<"theta1 = "<<theta1<<G4endl;
00905 
00906     E1 = fEnergyVector->GetLowEdgeEnergy(iMomentum);
00907 
00908     // G4cout<<"E1 = "<<E1<<G4endl;
00909 
00910     W  = 1.0/(E2 - E1);
00911     W1 = (E2 - kinE)*W;
00912     W2 = (kinE - E1)*W;
00913 
00914     randAngle = W1*theta1 + W2*theta2;
00915     
00916     // randAngle = theta2;
00917   }
00918   // G4double angle = randAngle;
00919   // if (randAngle > 0.) randAngle /= 2*pi*std::sin(angle);
00920   // G4cout<<"randAngle = "<<randAngle/degree<<G4endl;
00921 
00922   return randAngle;
00923 }

G4double G4NuclNuclDiffuseElastic::SampleThetaCMS	(	const G4ParticleDefinition *	aParticle,
		G4double	p,
		G4double	A
	)

Definition at line 712 of file G4NuclNuclDiffuseElastic.cc.

References CalculateNuclearRad(), G4UniformRand, GetIntegrandFunction(), G4Integrator< T, F >::Legendre10(), G4Integrator< T, F >::Legendre96(), and G4INCL::Math::pi.

Referenced by SampleT().

00714 {
00715   G4int i, iMax = 100;  
00716   G4double norm, result, theta1, theta2, thetaMax, sum = 0.;
00717 
00718   fParticle      = particle;
00719   fWaveVector    = momentum/hbarc;
00720   fAtomicWeight  = A;
00721 
00722   fNuclearRadius = CalculateNuclearRad(A);
00723   
00724   thetaMax = 10.174/fWaveVector/fNuclearRadius;
00725 
00726   if (thetaMax > pi) thetaMax = pi;
00727 
00728   G4Integrator<G4NuclNuclDiffuseElastic,G4double(G4NuclNuclDiffuseElastic::*)(G4double)> integral;
00729 
00730   // result = integral.Legendre10(this,&G4NuclNuclDiffuseElastic::GetIntegrandFunction, 0., theta ); 
00731   norm = integral.Legendre96(this,&G4NuclNuclDiffuseElastic::GetIntegrandFunction, 0., thetaMax );
00732 
00733   norm *= G4UniformRand();
00734 
00735   for(i = 1; i <= iMax; i++)
00736   {
00737     theta1 = (i-1)*thetaMax/iMax; 
00738     theta2 = i*thetaMax/iMax;
00739     sum   += integral.Legendre10(this,&G4NuclNuclDiffuseElastic::GetIntegrandFunction, theta1, theta2);
00740 
00741     if ( sum >= norm ) 
00742     {
00743       result = 0.5*(theta1 + theta2);
00744       break;
00745     }
00746   }
00747   if (i > iMax ) result = 0.5*(theta1 + theta2);
00748 
00749   G4double sigma = pi*thetaMax/iMax;
00750 
00751   result += G4RandGauss::shoot(0.,sigma);
00752 
00753   if(result < 0.) result = 0.;
00754   if(result > thetaMax) result = thetaMax;
00755 
00756   return result;
00757 }

G4double G4NuclNuclDiffuseElastic::SampleThetaLab	(	const G4HadProjectile *	aParticle,
		G4double	tmass,
		G4double	A
	)

Definition at line 1078 of file G4NuclNuclDiffuseElastic.cc.

References G4cout, G4endl, G4UniformRand, G4HadProjectile::Get4Momentum(), G4HadProjectile::GetDefinition(), G4ParticleDefinition::GetPDGMass(), G4HadProjectile::GetTotalMomentum(), SampleT(), and G4HadronicInteraction::verboseLevel.

01080 {
01081   const G4ParticleDefinition* theParticle = aParticle->GetDefinition();
01082   G4double m1 = theParticle->GetPDGMass();
01083   G4double plab = aParticle->GetTotalMomentum();
01084   G4LorentzVector lv1 = aParticle->Get4Momentum();
01085   G4LorentzVector lv(0.0,0.0,0.0,tmass);   
01086   lv += lv1;
01087 
01088   G4ThreeVector bst = lv.boostVector();
01089   lv1.boost(-bst);
01090 
01091   G4ThreeVector p1 = lv1.vect();
01092   G4double ptot    = p1.mag();
01093   G4double tmax    = 4.0*ptot*ptot;
01094   G4double t       = 0.0;
01095 
01096 
01097   //
01098   // Sample t
01099   //
01100   
01101   t = SampleT( theParticle, ptot, A);
01102 
01103   // NaN finder
01104   if(!(t < 0.0 || t >= 0.0)) 
01105   {
01106     if (verboseLevel > 0) 
01107     {
01108       G4cout << "G4NuclNuclDiffuseElastic:WARNING: A = " << A 
01109              << " mom(GeV)= " << plab/GeV 
01110              << " S-wave will be sampled" 
01111              << G4endl; 
01112     }
01113     t = G4UniformRand()*tmax; 
01114   }
01115   if(verboseLevel>1)
01116   {
01117     G4cout <<" t= " << t << " tmax= " << tmax 
01118            << " ptot= " << ptot << G4endl;
01119   }
01120   // Sampling of angles in CM system
01121 
01122   G4double phi  = G4UniformRand()*twopi;
01123   G4double cost = 1. - 2.0*t/tmax;
01124   G4double sint;
01125 
01126   if( cost >= 1.0 ) 
01127   {
01128     cost = 1.0;
01129     sint = 0.0;
01130   }
01131   else if( cost <= -1.0) 
01132   {
01133     cost = -1.0;
01134     sint =  0.0;
01135   }
01136   else  
01137   {
01138     sint = std::sqrt((1.0-cost)*(1.0+cost));
01139   }    
01140   if (verboseLevel>1) 
01141   {
01142     G4cout << "cos(t)=" << cost << " std::sin(t)=" << sint << G4endl;
01143   }
01144   G4ThreeVector v1(sint*std::cos(phi),sint*std::sin(phi),cost);
01145   v1 *= ptot;
01146   G4LorentzVector nlv1(v1.x(),v1.y(),v1.z(),std::sqrt(ptot*ptot + m1*m1));
01147 
01148   nlv1.boost(bst); 
01149 
01150   G4ThreeVector np1 = nlv1.vect();
01151 
01152     // G4double theta = std::acos( np1.z()/np1.mag() );  // degree;
01153 
01154   G4double theta = np1.theta();
01155 
01156   return theta;
01157 }

void G4NuclNuclDiffuseElastic::SetCofAlpha ( G4double pa ) [inline]

Definition at line 289 of file G4NuclNuclDiffuseElastic.hh.

00289 {fCofAlpha = pa;};

void G4NuclNuclDiffuseElastic::SetCofAlphaCoulomb ( G4double pa ) [inline]

Definition at line 291 of file G4NuclNuclDiffuseElastic.hh.

00291 {fCofAlphaCoulomb = pa;};

void G4NuclNuclDiffuseElastic::SetCofAlphaMax ( G4double pa ) [inline]

Definition at line 290 of file G4NuclNuclDiffuseElastic.hh.

00290 {fCofAlphaMax = pa;};

void G4NuclNuclDiffuseElastic::SetCofDelta ( G4double pa ) [inline]

Definition at line 293 of file G4NuclNuclDiffuseElastic.hh.

00293 {fCofDelta = pa;};

void G4NuclNuclDiffuseElastic::SetCofFar ( G4double pa ) [inline]

Definition at line 295 of file G4NuclNuclDiffuseElastic.hh.

00295 {fCofFar = pa;};

void G4NuclNuclDiffuseElastic::SetCofLambda ( G4double pa ) [inline]

Definition at line 287 of file G4NuclNuclDiffuseElastic.hh.

00287 {fCofLambda = pa;};

void G4NuclNuclDiffuseElastic::SetCofPhase ( G4double pa ) [inline]

Definition at line 294 of file G4NuclNuclDiffuseElastic.hh.

00294 {fCofPhase = pa;};

void G4NuclNuclDiffuseElastic::SetEtaRatio ( G4double pa ) [inline]

Definition at line 296 of file G4NuclNuclDiffuseElastic.hh.

00296 {fEtaRatio = pa;};

void G4NuclNuclDiffuseElastic::SetHEModelLowLimit ( G4double value ) [inline]

Definition at line 385 of file G4NuclNuclDiffuseElastic.hh.

00386 {
00387   lowEnergyLimitHE = value;
00388 }

void G4NuclNuclDiffuseElastic::SetLowestEnergyLimit ( G4double value ) [inline]

Reimplemented from G4HadronElastic.

Definition at line 395 of file G4NuclNuclDiffuseElastic.hh.

00396 {
00397   lowestEnergyLimit = value;
00398 }

void G4NuclNuclDiffuseElastic::SetMaxL ( G4int l ) [inline]

Definition at line 297 of file G4NuclNuclDiffuseElastic.hh.

00297 {fMaxL = l;};

void G4NuclNuclDiffuseElastic::SetNuclearRadiusCof ( G4double r ) [inline]

Definition at line 298 of file G4NuclNuclDiffuseElastic.hh.

00298 {fNuclearRadiusCof = r;};

void G4NuclNuclDiffuseElastic::SetPlabLowLimit ( G4double value ) [inline]

Definition at line 380 of file G4NuclNuclDiffuseElastic.hh.

00381 {
00382   plabLowLimit = value;
00383 }

void G4NuclNuclDiffuseElastic::SetProfileAlpha ( G4double pa ) [inline]

Definition at line 286 of file G4NuclNuclDiffuseElastic.hh.

00286 {fProfileAlpha = pa;};

void G4NuclNuclDiffuseElastic::SetProfileDelta ( G4double pd ) [inline]

Definition at line 285 of file G4NuclNuclDiffuseElastic.hh.

00285 {fProfileDelta = pd;};

void G4NuclNuclDiffuseElastic::SetProfileLambda ( G4double pl ) [inline]

Definition at line 284 of file G4NuclNuclDiffuseElastic.hh.

00284 {fProfileLambda = pl;};

void G4NuclNuclDiffuseElastic::SetQModelLowLimit ( G4double value ) [inline]

Definition at line 390 of file G4NuclNuclDiffuseElastic.hh.

00391 {
00392   lowEnergyLimitQ = value;
00393 }

void G4NuclNuclDiffuseElastic::SetRecoilKinEnergyLimit ( G4double value ) [inline]

Definition at line 375 of file G4NuclNuclDiffuseElastic.hh.

00376 {
00377   lowEnergyRecoilLimit = value;
00378 }

void G4NuclNuclDiffuseElastic::TestAngleTable	(	const G4ParticleDefinition *	theParticle,
		G4double	partMom,
		G4double	Z,
		G4double	A
	)

Definition at line 1286 of file G4NuclNuclDiffuseElastic.cc.

References G4Integrator< T, F >::AdaptiveGauss(), CalculateAm(), CalculateNuclearRad(), CalculateZommerfeld(), G4cout, G4endl, GetIntegrandFunction(), G4ParticleDefinition::GetPDGCharge(), G4ParticleDefinition::GetPDGMass(), G4PhysicsTable::insertAt(), G4Integrator< T, F >::Legendre10(), G4Integrator< T, F >::Legendre96(), and G4PhysicsFreeVector::PutValue().

01288 {
01289   fAtomicNumber  = Z;     // atomic number
01290   fAtomicWeight  = A;     // number of nucleons
01291   fNuclearRadius = CalculateNuclearRad(fAtomicWeight);
01292   
01293      
01294   
01295   G4cout<<"G4NuclNuclDiffuseElastic::TestAngleTable() init the element with Z = "
01296           <<Z<<"; and A = "<<A<<G4endl;
01297  
01298   fElementNumberVector.push_back(fAtomicNumber);
01299 
01300  
01301 
01302 
01303   G4int i=0, j;
01304   G4double a = 0., z = theParticle->GetPDGCharge(),  m1 = fParticle->GetPDGMass();
01305   G4double alpha1=0., alpha2=0., alphaMax=0., alphaCoulomb=0.;
01306   G4double deltaL10 = 0., deltaL96 = 0., deltaAG = 0.;
01307   G4double sumL10 = 0.,sumL96 = 0.,sumAG = 0.;
01308   G4double epsilon = 0.001;
01309 
01310   G4Integrator<G4NuclNuclDiffuseElastic,G4double(G4NuclNuclDiffuseElastic::*)(G4double)> integral;
01311   
01312   fAngleTable = new G4PhysicsTable(fEnergyBin);
01313 
01314   fWaveVector = partMom/hbarc;
01315 
01316   G4double kR     = fWaveVector*fNuclearRadius;
01317   G4double kR2    = kR*kR;
01318   G4double kRmax  = 10.6; // 10.6, 18, 10.174; ~ 3 maxima of J1 or 15., 25.
01319   G4double kRcoul = 1.2; // 1.4, 2.5; // on the first slope of J1
01320 
01321   alphaMax = kRmax*kRmax/kR2;
01322 
01323   if (alphaMax > 4.) alphaMax = 4.;  // vmg05-02-09: was pi2 
01324 
01325   alphaCoulomb = kRcoul*kRcoul/kR2;
01326 
01327   if( z )
01328   {
01329       a           = partMom/m1; // beta*gamma for m1
01330       fBeta       = a/std::sqrt(1+a*a);
01331       fZommerfeld = CalculateZommerfeld( fBeta, z, fAtomicNumber);
01332       fAm         = CalculateAm( partMom, fZommerfeld, fAtomicNumber);
01333   }
01334   G4PhysicsFreeVector* angleVector = new G4PhysicsFreeVector(fAngleBin-1);
01335 
01336   // G4PhysicsLogVector*  angleBins = new G4PhysicsLogVector( 0.001*alphaMax, alphaMax, fAngleBin );
01337     
01338   
01339   fAddCoulomb = false;
01340 
01341   for(j = 1; j < fAngleBin; j++)
01342   {
01343       // alpha1 = angleBins->GetLowEdgeEnergy(j-1);
01344       // alpha2 = angleBins->GetLowEdgeEnergy(j);
01345 
01346     alpha1 = alphaMax*(j-1)/fAngleBin;
01347     alpha2 = alphaMax*( j )/fAngleBin;
01348 
01349     if( ( alpha2 > alphaCoulomb ) && z ) fAddCoulomb = true;
01350 
01351     deltaL10 = integral.Legendre10(this, &G4NuclNuclDiffuseElastic::GetIntegrandFunction, alpha1, alpha2);
01352     deltaL96 = integral.Legendre96(this, &G4NuclNuclDiffuseElastic::GetIntegrandFunction, alpha1, alpha2);
01353     deltaAG  = integral.AdaptiveGauss(this, &G4NuclNuclDiffuseElastic::GetIntegrandFunction, 
01354                                        alpha1, alpha2,epsilon);
01355 
01356       // G4cout<<alpha1<<"\t"<<std::sqrt(alpha1)/degree<<"\t"
01357       //     <<deltaL10<<"\t"<<deltaL96<<"\t"<<deltaAG<<G4endl;
01358 
01359     sumL10 += deltaL10;
01360     sumL96 += deltaL96;
01361     sumAG  += deltaAG;
01362 
01363     G4cout<<alpha1<<"\t"<<std::sqrt(alpha1)/degree<<"\t"
01364             <<sumL10<<"\t"<<sumL96<<"\t"<<sumAG<<G4endl;
01365       
01366     angleVector->PutValue( j-1 , alpha1, sumL10 ); // alpha2
01367   }
01368   fAngleTable->insertAt(i,angleVector);
01369   fAngleBank.push_back(fAngleTable);
01370 
01371   /*
01372   // Integral over all angle range - Bad accuracy !!!
01373 
01374   sumL10 = integral.Legendre10(this, &G4NuclNuclDiffuseElastic::GetIntegrandFunction, 0., alpha2);
01375   sumL96 = integral.Legendre96(this, &G4NuclNuclDiffuseElastic::GetIntegrandFunction, 0., alpha2);
01376   sumAG  = integral.AdaptiveGauss(this, &G4NuclNuclDiffuseElastic::GetIntegrandFunction, 
01377                                        0., alpha2,epsilon);
01378   G4cout<<G4endl;
01379   G4cout<<alpha2<<"\t"<<std::sqrt(alpha2)/degree<<"\t"
01380             <<sumL10<<"\t"<<sumL96<<"\t"<<sumAG<<G4endl;
01381   */
01382   return;
01383 }

G4complex G4NuclNuclDiffuseElastic::TestErfcComp	(	G4complex	z,
		G4int	nMax
	)			`[inline]`

Definition at line 842 of file G4NuclNuclDiffuseElastic.hh.

References GetErfComp().

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 }

G4complex G4NuclNuclDiffuseElastic::TestErfcInt ( G4complex z ) [inline]

Definition at line 866 of file G4NuclNuclDiffuseElastic.hh.

References GetErfInt().

00867 {
00868   G4complex miz = G4complex( z.imag(), -z.real() ); 
00869   G4complex erfcz = 1. - GetErfInt( miz); // , nMax);
00870   G4complex w = std::exp(-z*z)*erfcz;
00871   return w;
00872 }

G4complex G4NuclNuclDiffuseElastic::TestErfcSer	(	G4complex	z,
		G4int	nMax
	)			`[inline]`

Definition at line 854 of file G4NuclNuclDiffuseElastic.hh.

References GetErfSer().

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 }

G4double G4NuclNuclDiffuseElastic::ThetaCMStoThetaLab	(	const G4DynamicParticle *	aParticle,
		G4double	tmass,
		G4double	thetaCMS
	)

Definition at line 1166 of file G4NuclNuclDiffuseElastic.cc.

References G4cout, G4endl, G4UniformRand, G4DynamicParticle::Get4Momentum(), G4DynamicParticle::GetDefinition(), G4ParticleDefinition::GetPDGMass(), and G4HadronicInteraction::verboseLevel.

01168 {
01169   const G4ParticleDefinition* theParticle = aParticle->GetDefinition();
01170   G4double m1 = theParticle->GetPDGMass();
01171   // G4double plab = aParticle->GetTotalMomentum();
01172   G4LorentzVector lv1 = aParticle->Get4Momentum();
01173   G4LorentzVector lv(0.0,0.0,0.0,tmass);   
01174 
01175   lv += lv1;
01176 
01177   G4ThreeVector bst = lv.boostVector();
01178 
01179   lv1.boost(-bst);
01180 
01181   G4ThreeVector p1 = lv1.vect();
01182   G4double ptot    = p1.mag();
01183 
01184   G4double phi  = G4UniformRand()*twopi;
01185   G4double cost = std::cos(thetaCMS);
01186   G4double sint;
01187 
01188   if( cost >= 1.0 ) 
01189   {
01190     cost = 1.0;
01191     sint = 0.0;
01192   }
01193   else if( cost <= -1.0) 
01194   {
01195     cost = -1.0;
01196     sint =  0.0;
01197   }
01198   else  
01199   {
01200     sint = std::sqrt((1.0-cost)*(1.0+cost));
01201   }    
01202   if (verboseLevel>1) 
01203   {
01204     G4cout << "cos(tcms)=" << cost << " std::sin(tcms)=" << sint << G4endl;
01205   }
01206   G4ThreeVector v1(sint*std::cos(phi),sint*std::sin(phi),cost);
01207   v1 *= ptot;
01208   G4LorentzVector nlv1(v1.x(),v1.y(),v1.z(),std::sqrt(ptot*ptot + m1*m1));
01209 
01210   nlv1.boost(bst); 
01211 
01212   G4ThreeVector np1 = nlv1.vect();
01213 
01214 
01215   G4double thetaLab = np1.theta();
01216 
01217   return thetaLab;
01218 }

G4double G4NuclNuclDiffuseElastic::ThetaLabToThetaCMS	(	const G4DynamicParticle *	aParticle,
		G4double	tmass,
		G4double	thetaLab
	)

Definition at line 1227 of file G4NuclNuclDiffuseElastic.cc.

References G4cout, G4endl, G4UniformRand, G4DynamicParticle::Get4Momentum(), G4DynamicParticle::GetDefinition(), G4ParticleDefinition::GetPDGMass(), G4DynamicParticle::GetTotalMomentum(), and G4HadronicInteraction::verboseLevel.

01229 {
01230   const G4ParticleDefinition* theParticle = aParticle->GetDefinition();
01231   G4double m1 = theParticle->GetPDGMass();
01232   G4double plab = aParticle->GetTotalMomentum();
01233   G4LorentzVector lv1 = aParticle->Get4Momentum();
01234   G4LorentzVector lv(0.0,0.0,0.0,tmass);   
01235 
01236   lv += lv1;
01237 
01238   G4ThreeVector bst = lv.boostVector();
01239 
01240   // lv1.boost(-bst);
01241 
01242   // G4ThreeVector p1 = lv1.vect();
01243   // G4double ptot    = p1.mag();
01244 
01245   G4double phi  = G4UniformRand()*twopi;
01246   G4double cost = std::cos(thetaLab);
01247   G4double sint;
01248 
01249   if( cost >= 1.0 ) 
01250   {
01251     cost = 1.0;
01252     sint = 0.0;
01253   }
01254   else if( cost <= -1.0) 
01255   {
01256     cost = -1.0;
01257     sint =  0.0;
01258   }
01259   else  
01260   {
01261     sint = std::sqrt((1.0-cost)*(1.0+cost));
01262   }    
01263   if (verboseLevel>1) 
01264   {
01265     G4cout << "cos(tlab)=" << cost << " std::sin(tlab)=" << sint << G4endl;
01266   }
01267   G4ThreeVector v1(sint*std::cos(phi),sint*std::sin(phi),cost);
01268   v1 *= plab;
01269   G4LorentzVector nlv1(v1.x(),v1.y(),v1.z(),std::sqrt(plab*plab + m1*m1));
01270 
01271   nlv1.boost(-bst); 
01272 
01273   G4ThreeVector np1 = nlv1.vect();
01274 
01275 
01276   G4double thetaCMS = np1.theta();
01277 
01278   return thetaCMS;
01279 }


Public Member Functions
	G4NuclNuclDiffuseElastic ()
virtual	~G4NuclNuclDiffuseElastic ()
void	Initialise ()
void	InitialiseOnFly (G4double Z, G4double A)
void	BuildAngleTable ()
virtual G4double	SampleInvariantT (const G4ParticleDefinition *p, G4double plab, G4int Z, G4int A)
void	SetPlabLowLimit (G4double value)
void	SetHEModelLowLimit (G4double value)
void	SetQModelLowLimit (G4double value)
void	SetLowestEnergyLimit (G4double value)
void	SetRecoilKinEnergyLimit (G4double value)
G4double	SampleT (const G4ParticleDefinition *aParticle, G4double p, G4double A)
G4double	SampleTableT (const G4ParticleDefinition *aParticle, G4double p, G4double Z, G4double A)
G4double	SampleThetaCMS (const G4ParticleDefinition *aParticle, G4double p, G4double A)
G4double	SampleTableThetaCMS (const G4ParticleDefinition *aParticle, G4double p, G4double Z, G4double A)
G4double	GetScatteringAngle (G4int iMomentum, G4int iAngle, G4double position)
G4double	SampleThetaLab (const G4HadProjectile *aParticle, G4double tmass, G4double A)
G4double	GetDiffuseElasticXsc (const G4ParticleDefinition *particle, G4double theta, G4double momentum, G4double A)
G4double	GetInvElasticXsc (const G4ParticleDefinition *particle, G4double theta, G4double momentum, G4double A, G4double Z)
G4double	GetDiffuseElasticSumXsc (const G4ParticleDefinition *particle, G4double theta, G4double momentum, G4double A, G4double Z)
G4double	GetInvElasticSumXsc (const G4ParticleDefinition *particle, G4double tMand, G4double momentum, G4double A, G4double Z)
G4double	IntegralElasticProb (const G4ParticleDefinition *particle, G4double theta, G4double momentum, G4double A)
G4double	GetCoulombElasticXsc (const G4ParticleDefinition *particle, G4double theta, G4double momentum, G4double Z)
G4double	GetRutherfordXsc (G4double theta)
G4double	GetInvCoulombElasticXsc (const G4ParticleDefinition *particle, G4double tMand, G4double momentum, G4double A, G4double Z)
G4double	GetCoulombTotalXsc (const G4ParticleDefinition *particle, G4double momentum, G4double Z)
G4double	GetCoulombIntegralXsc (const G4ParticleDefinition *particle, G4double momentum, G4double Z, G4double theta1, G4double theta2)
G4double	CalculateParticleBeta (const G4ParticleDefinition *particle, G4double momentum)
G4double	CalculateZommerfeld (G4double beta, G4double Z1, G4double Z2)
G4double	CalculateAm (G4double momentum, G4double n, G4double Z)
G4double	CalculateNuclearRad (G4double A)
G4double	ThetaCMStoThetaLab (const G4DynamicParticle *aParticle, G4double tmass, G4double thetaCMS)
G4double	ThetaLabToThetaCMS (const G4DynamicParticle *aParticle, G4double tmass, G4double thetaLab)
void	TestAngleTable (const G4ParticleDefinition *theParticle, G4double partMom, G4double Z, G4double A)
G4double	BesselJzero (G4double z)
G4double	BesselJone (G4double z)
G4double	DampFactor (G4double z)
G4double	BesselOneByArg (G4double z)
G4double	GetDiffElasticProb (G4double theta)
G4double	GetDiffElasticSumProb (G4double theta)
G4double	GetDiffElasticSumProbA (G4double alpha)
G4double	GetIntegrandFunction (G4double theta)
G4double	GetNuclearRadius ()
G4complex	GammaLogarithm (G4complex xx)
G4complex	GammaLogB2n (G4complex xx)
G4double	GetErf (G4double x)
G4double	GetCosHaPit2 (G4double t)
G4double	GetSinHaPit2 (G4double t)
G4double	GetCint (G4double x)
G4double	GetSint (G4double x)
G4complex	GetErfcComp (G4complex z, G4int nMax)
G4complex	GetErfcSer (G4complex z, G4int nMax)
G4complex	GetErfcInt (G4complex z)
G4complex	GetErfComp (G4complex z, G4int nMax)
G4complex	GetErfSer (G4complex z, G4int nMax)
G4double	GetExpCos (G4double x)
G4double	GetExpSin (G4double x)
G4complex	GetErfInt (G4complex z)
G4double	GetLegendrePol (G4int n, G4double x)
G4complex	TestErfcComp (G4complex z, G4int nMax)
G4complex	TestErfcSer (G4complex z, G4int nMax)
G4complex	TestErfcInt (G4complex z)
G4complex	CoulombAmplitude (G4double theta)
G4double	CoulombAmplitudeMod2 (G4double theta)
void	CalculateCoulombPhaseZero ()
G4double	CalculateCoulombPhase (G4int n)
void	CalculateRutherfordAnglePar ()
G4double	ProfileNear (G4double theta)
G4double	ProfileFar (G4double theta)
G4double	Profile (G4double theta)
G4complex	PhaseNear (G4double theta)
G4complex	PhaseFar (G4double theta)
G4complex	GammaLess (G4double theta)
G4complex	GammaMore (G4double theta)
G4complex	AmplitudeNear (G4double theta)
G4complex	AmplitudeFar (G4double theta)
G4complex	Amplitude (G4double theta)
G4double	AmplitudeMod2 (G4double theta)
G4complex	AmplitudeSim (G4double theta)
G4double	AmplitudeSimMod2 (G4double theta)
G4double	GetRatioSim (G4double theta)
G4double	GetRatioGen (G4double theta)
G4double	GetFresnelDiffuseXsc (G4double theta)
G4double	GetFresnelIntegrandXsc (G4double alpha)
G4complex	AmplitudeGla (G4double theta)
G4double	AmplitudeGlaMod2 (G4double theta)
G4complex	AmplitudeGG (G4double theta)
G4double	AmplitudeGGMod2 (G4double theta)
void	InitParameters (const G4ParticleDefinition *theParticle, G4double partMom, G4double Z, G4double A)
void	InitDynParameters (const G4ParticleDefinition *theParticle, G4double partMom)
void	InitParametersGla (const G4DynamicParticle *aParticle, G4double partMom, G4double Z, G4double A)
G4double	GetHadronNucleonXscNS (G4ParticleDefinition pParticle, G4double pTkin, G4ParticleDefinition tParticle)
G4double	CalcMandelstamS (const G4double mp, const G4double mt, const G4double Plab)
G4double	GetProfileLambda ()
void	SetProfileLambda (G4double pl)
void	SetProfileDelta (G4double pd)
void	SetProfileAlpha (G4double pa)
void	SetCofLambda (G4double pa)
void	SetCofAlpha (G4double pa)
void	SetCofAlphaMax (G4double pa)
void	SetCofAlphaCoulomb (G4double pa)
void	SetCofDelta (G4double pa)
void	SetCofPhase (G4double pa)
void	SetCofFar (G4double pa)
void	SetEtaRatio (G4double pa)
void	SetMaxL (G4int l)
void	SetNuclearRadiusCof (G4double r)
G4double	GetCofAlphaMax ()
G4double	GetCofAlphaCoulomb ()

G4NuclNuclDiffuseElastic Class Reference

Public Member Functions

Detailed Description

Constructor & Destructor Documentation

Member Function Documentation