G4DiffuseElastic Class Reference

#include <G4DiffuseElastic.hh>

Inheritance diagram for G4DiffuseElastic:

G4HadronElastic G4HadronicInteraction

Public Member Functions

 G4DiffuseElastic ()
virtual ~G4DiffuseElastic ()
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 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 ()

Detailed Description

Definition at line 58 of file G4DiffuseElastic.hh.

Constructor & Destructor Documentation

G4DiffuseElastic::G4DiffuseElastic (  ) 

Definition at line 67 of file G4DiffuseElastic.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("DiffuseElastic"), fParticle(0)
00069 {
00070   SetMinEnergy( 0.01*GeV );
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 }

G4DiffuseElastic::~G4DiffuseElastic (  )  [virtual]

Definition at line 107 of file G4DiffuseElastic.cc.

References G4PhysicsTable::clearAndDestroy().

00108 {
00109   if(fEnergyVector) delete fEnergyVector;
00110 
00111   if( fAngleTable )
00112   {
00113       fAngleTable->clearAndDestroy();
00114       delete fAngleTable ;
00115   }
00116 }

Member Function Documentation

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

Definition at line 311 of file G4DiffuseElastic.hh.

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

00312 {
00313   G4double modvalue, value2, fact1, fact2, arg, shift, bessel;
00314 
00315   modvalue = fabs(value);
00316 
00317   if ( modvalue < 8.0 ) 
00318   {
00319     value2 = value*value;
00320 
00321     fact1  = value*(72362614232.0 + value2*(-7895059235.0 
00322                                   + value2*( 242396853.1
00323                                   + value2*(-2972611.439 
00324                                   + value2*( 15704.48260 
00325                                   + value2*(-30.16036606  ) ) ) ) ) );
00326 
00327     fact2  = 144725228442.0 + value2*(2300535178.0 
00328                             + value2*(18583304.74
00329                             + value2*(99447.43394 
00330                             + value2*(376.9991397 
00331                             + value2*1.0             ) ) ) );
00332     bessel = fact1/fact2;
00333   } 
00334   else 
00335   {
00336     arg    = 8.0/modvalue;
00337 
00338     value2 = arg*arg;
00339 
00340     shift  = modvalue - 2.356194491;
00341 
00342     fact1  = 1.0 + value2*( 0.183105e-2 
00343                  + value2*(-0.3516396496e-4
00344                  + value2*(0.2457520174e-5 
00345                  + value2*(-0.240337019e-6          ) ) ) );
00346 
00347     fact2 = 0.04687499995 + value2*(-0.2002690873e-3
00348                           + value2*( 0.8449199096e-5
00349                           + value2*(-0.88228987e-6
00350                           + value2*0.105787412e-6       ) ) );
00351 
00352     bessel = sqrt( 0.636619772/modvalue)*(cos(shift)*fact1 - arg*sin(shift)*fact2);
00353 
00354     if (value < 0.0) bessel = -bessel;
00355   }
00356   return bessel;
00357 }

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

Definition at line 259 of file G4DiffuseElastic.hh.

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

00260 {
00261   G4double modvalue, value2, fact1, fact2, arg, shift, bessel;
00262 
00263   modvalue = fabs(value);
00264 
00265   if ( value < 8.0 && value > -8.0 )
00266   {
00267     value2 = value*value;
00268 
00269     fact1  = 57568490574.0 + value2*(-13362590354.0 
00270                            + value2*( 651619640.7 
00271                            + value2*(-11214424.18 
00272                            + value2*( 77392.33017 
00273                            + value2*(-184.9052456   ) ) ) ) );
00274 
00275     fact2  = 57568490411.0 + value2*( 1029532985.0 
00276                            + value2*( 9494680.718
00277                            + value2*(59272.64853
00278                            + value2*(267.8532712 
00279                            + value2*1.0               ) ) ) );
00280 
00281     bessel = fact1/fact2;
00282   } 
00283   else 
00284   {
00285     arg    = 8.0/modvalue;
00286 
00287     value2 = arg*arg;
00288 
00289     shift  = modvalue-0.785398164;
00290 
00291     fact1  = 1.0 + value2*(-0.1098628627e-2 
00292                  + value2*(0.2734510407e-4
00293                  + value2*(-0.2073370639e-5 
00294                  + value2*0.2093887211e-6    ) ) );
00295 
00296     fact2  = -0.1562499995e-1 + value2*(0.1430488765e-3
00297                               + value2*(-0.6911147651e-5 
00298                               + value2*(0.7621095161e-6
00299                               - value2*0.934945152e-7    ) ) );
00300 
00301     bessel = sqrt(0.636619772/modvalue)*(cos(shift)*fact1 - arg*sin(shift)*fact2 );
00302   }
00303   return bessel;
00304 }

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

Definition at line 386 of file G4DiffuseElastic.hh.

References BesselJone().

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

00387 {
00388   G4double x2, result;
00389   
00390   if( std::fabs(x) < 0.01 )
00391   { 
00392    x     *= 0.5;
00393    x2     = x*x;
00394    result = 2. - x2 + x2*x2/6.;
00395   }
00396   else
00397   {
00398     result = BesselJone(x)/x; 
00399   }
00400   return result;
00401 }

void G4DiffuseElastic::BuildAngleTable (  ) 

Definition at line 922 of file G4DiffuseElastic.cc.

References CalculateAm(), CalculateZommerfeld(), GetIntegrandFunction(), G4PhysicsVector::GetLowEdgeEnergy(), G4ParticleDefinition::GetPDGCharge(), G4ParticleDefinition::GetPDGMass(), G4PhysicsTable::insertAt(), G4Integrator< T, F >::Legendre10(), and G4PhysicsFreeVector::PutValue().

Referenced by Initialise(), and InitialiseOnFly().

00923 {
00924   G4int i, j;
00925   G4double partMom, kinE, a = 0., z = fParticle->GetPDGCharge(), m1 = fParticle->GetPDGMass();
00926   G4double alpha1, alpha2, alphaMax, alphaCoulomb, delta = 0., sum = 0.;
00927 
00928   G4Integrator<G4DiffuseElastic,G4double(G4DiffuseElastic::*)(G4double)> integral;
00929   
00930   fAngleTable = new G4PhysicsTable(fEnergyBin);
00931 
00932   for( i = 0; i < fEnergyBin; i++)
00933   {
00934     kinE        = fEnergyVector->GetLowEdgeEnergy(i);
00935     partMom     = std::sqrt( kinE*(kinE + 2*m1) );
00936 
00937     fWaveVector = partMom/hbarc;
00938 
00939     G4double kR     = fWaveVector*fNuclearRadius;
00940     G4double kR2    = kR*kR;
00941     G4double kRmax  = 18.6; // 10.6; 10.6, 18, 10.174; ~ 3 maxima of J1 or 15., 25.
00942     G4double kRcoul = 1.9; // 1.2; 1.4, 2.5; // on the first slope of J1
00943     // G4double kRlim  = 1.2;
00944     // G4double kRlim2 = kRlim*kRlim/kR2;
00945 
00946     alphaMax = kRmax*kRmax/kR2;
00947 
00948     if (alphaMax > 4.) alphaMax = 4.;  // vmg05-02-09: was pi2 
00949 
00950     alphaCoulomb = kRcoul*kRcoul/kR2;
00951 
00952     if( z )
00953     {
00954       a           = partMom/m1;         // beta*gamma for m1
00955       fBeta       = a/std::sqrt(1+a*a);
00956       fZommerfeld = CalculateZommerfeld( fBeta, z, fAtomicNumber);
00957       fAm         = CalculateAm( partMom, fZommerfeld, fAtomicNumber);
00958     }
00959     G4PhysicsFreeVector* angleVector = new G4PhysicsFreeVector(fAngleBin-1);
00960 
00961     // G4PhysicsLogVector*  angleBins = new G4PhysicsLogVector( 0.001*alphaMax, alphaMax, fAngleBin );
00962 
00963     G4double delth = alphaMax/fAngleBin;
00964         
00965     sum = 0.;
00966 
00967     // fAddCoulomb = false;
00968     fAddCoulomb = true;
00969 
00970     // for(j = 1; j < fAngleBin; j++)
00971     for(j = fAngleBin-1; j >= 1; j--)
00972     {
00973       // alpha1 = angleBins->GetLowEdgeEnergy(j-1);
00974       // alpha2 = angleBins->GetLowEdgeEnergy(j);
00975 
00976       // alpha1 = alphaMax*(j-1)/fAngleBin;
00977       // alpha2 = alphaMax*( j )/fAngleBin;
00978 
00979       alpha1 = delth*(j-1);
00980       // if(alpha1 < kRlim2) alpha1 = kRlim2;
00981       alpha2 = alpha1 + delth;
00982 
00983       // if( ( alpha2 > alphaCoulomb ) && z ) fAddCoulomb = true;
00984       if( ( alpha1 < alphaCoulomb ) && z ) fAddCoulomb = false;
00985 
00986       delta = integral.Legendre10(this, &G4DiffuseElastic::GetIntegrandFunction, alpha1, alpha2);
00987       // delta = integral.Legendre96(this, &G4DiffuseElastic::GetIntegrandFunction, alpha1, alpha2);
00988 
00989       sum += delta;
00990       
00991       angleVector->PutValue( j-1 , alpha1, sum ); // alpha2
00992       // G4cout<<"j-1 = "<<j-1<<"; alpha2 = "<<alpha2<<"; sum = "<<sum<<G4endl;
00993     }
00994     fAngleTable->insertAt(i,angleVector);
00995 
00996     // delete[] angleVector; 
00997     // delete[] angleBins; 
00998   }
00999   return;
01000 }

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

Definition at line 432 of file G4DiffuseElastic.hh.

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

00433 {
00434   G4double k   = momentum/CLHEP::hbarc;
00435   G4double ch  = 1.13 + 3.76*n*n;
00436   G4double zn  = 1.77*k*std::pow(Z,-1./3.)*CLHEP::Bohr_radius;
00437   G4double zn2 = zn*zn;
00438   fAm          = ch/zn2;
00439 
00440   return fAm;
00441 }

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

Definition at line 447 of file G4DiffuseElastic.hh.

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

00448 {
00449   G4double r0;
00450 
00451   if( A < 50. )
00452   {
00453     if( A > 10. ) r0  = 1.16*( 1 - std::pow(A, -2./3.) )*CLHEP::fermi;   // 1.08*fermi;
00454     else          r0  = 1.1*CLHEP::fermi;
00455 
00456     fNuclearRadius = r0*std::pow(A, 1./3.);
00457   }
00458   else
00459   {
00460     r0 = 1.7*CLHEP::fermi;   // 1.7*fermi;
00461 
00462     fNuclearRadius = r0*std::pow(A, 0.27); // 0.27);
00463   }
00464   return fNuclearRadius;
00465 }

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

Definition at line 407 of file G4DiffuseElastic.hh.

References G4ParticleDefinition::GetPDGMass().

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

00409 {
00410   G4double mass = particle->GetPDGMass();
00411   G4double a    = momentum/mass;
00412   fBeta         = a/std::sqrt(1+a*a);
00413 
00414   return fBeta; 
00415 }

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

Definition at line 421 of file G4DiffuseElastic.hh.

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

00422 {
00423   fZommerfeld = CLHEP::fine_structure_const*Z1*Z2/beta;
00424 
00425   return fZommerfeld; 
00426 }

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

Definition at line 363 of file G4DiffuseElastic.hh.

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

00364 {
00365   G4double df;
00366   G4double f2 = 2., f3 = 6., f4 = 24.; // first factorials
00367 
00368   // x *= pi;
00369 
00370   if( std::fabs(x) < 0.01 )
00371   { 
00372     df = 1./(1. + x/f2 + x*x/f3 + x*x*x/f4);
00373   }
00374   else
00375   {
00376     df = x/std::sinh(x); 
00377   }
00378   return df;
00379 }

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

Definition at line 471 of file G4DiffuseElastic.hh.

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

Referenced by GetInvCoulombElasticXsc().

00475 {
00476   G4double sinHalfTheta  = std::sin(0.5*theta);
00477   G4double sinHalfTheta2 = sinHalfTheta*sinHalfTheta;
00478   G4double beta          = CalculateParticleBeta( particle, momentum);
00479   G4double z             = particle->GetPDGCharge();
00480   G4double n             = CalculateZommerfeld( beta, z, Z );
00481   G4double am            = CalculateAm( momentum, n, Z);
00482   G4double k             = momentum/CLHEP::hbarc;
00483   G4double ch            = 0.5*n/k;
00484   G4double ch2           = ch*ch;
00485   G4double xsc           = ch2/(sinHalfTheta2+am)/(sinHalfTheta2+am);
00486 
00487   return xsc;
00488 }

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

Definition at line 520 of file G4DiffuseElastic.hh.

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

00523 {
00524   G4double c1 = std::cos(theta1);
00525   G4cout<<"c1 = "<<c1<<G4endl;
00526   G4double c2 = std::cos(theta2);
00527   G4cout<<"c2 = "<<c2<<G4endl;
00528   G4double beta          = CalculateParticleBeta( particle, momentum);
00529   // G4cout<<"beta = "<<beta<<G4endl;
00530   G4double z             = particle->GetPDGCharge();
00531   G4double n             = CalculateZommerfeld( beta, z, Z );
00532   // G4cout<<"fZomerfeld = "<<n<<G4endl;
00533   G4double am            = CalculateAm( momentum, n, Z);
00534   // G4cout<<"cof Am = "<<am<<G4endl;
00535   G4double k             = momentum/CLHEP::hbarc;
00536   // G4cout<<"k = "<<k*CLHEP::fermi<<" 1/fermi"<<G4endl;
00537   // G4cout<<"k*Bohr_radius = "<<k*CLHEP::Bohr_radius<<G4endl;
00538   G4double ch            = n/k;
00539   G4double ch2           = ch*ch;
00540   am *= 2.;
00541   G4double xsc           = ch2*CLHEP::twopi*(c1-c2);
00542            xsc          /= (1 - c1 + am)*(1 - c2 + am);
00543 
00544   return xsc;
00545 }

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

Definition at line 495 of file G4DiffuseElastic.hh.

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

00497 {
00498   G4double beta          = CalculateParticleBeta( particle, momentum);
00499   G4cout<<"beta = "<<beta<<G4endl;
00500   G4double z             = particle->GetPDGCharge();
00501   G4double n             = CalculateZommerfeld( beta, z, Z );
00502   G4cout<<"fZomerfeld = "<<n<<G4endl;
00503   G4double am            = CalculateAm( momentum, n, Z);
00504   G4cout<<"cof Am = "<<am<<G4endl;
00505   G4double k             = momentum/CLHEP::hbarc;
00506   G4cout<<"k = "<<k*CLHEP::fermi<<" 1/fermi"<<G4endl;
00507   G4cout<<"k*Bohr_radius = "<<k*CLHEP::Bohr_radius<<G4endl;
00508   G4double ch            = n/k;
00509   G4double ch2           = ch*ch;
00510   G4double xsc           = ch2*CLHEP::pi/(am +am*am);
00511 
00512   return xsc;
00513 }

G4double G4DiffuseElastic::GetDiffElasticProb ( G4double  theta  ) 

Definition at line 363 of file G4DiffuseElastic.cc.

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

Referenced by GetDiffuseElasticXsc().

00368 {
00369   G4double sigma, bzero, bzero2, bonebyarg, bonebyarg2, damp, damp2;
00370   G4double delta, diffuse, gamma;
00371   G4double e1, e2, bone, bone2;
00372 
00373   // G4double wavek = momentum/hbarc;  // wave vector
00374   // G4double r0    = 1.08*fermi;
00375   // G4double rad   = r0*std::pow(A, 1./3.);
00376 
00377   G4double kr    = fWaveVector*fNuclearRadius; // wavek*rad;
00378   G4double kr2   = kr*kr;
00379   G4double krt   = kr*theta;
00380 
00381   bzero      = BesselJzero(krt);
00382   bzero2     = bzero*bzero;    
00383   bone       = BesselJone(krt);
00384   bone2      = bone*bone;
00385   bonebyarg  = BesselOneByArg(krt);
00386   bonebyarg2 = bonebyarg*bonebyarg;  
00387 
00388   if (fParticle == theProton)
00389   {
00390     diffuse = 0.63*fermi;
00391     gamma   = 0.3*fermi;
00392     delta   = 0.1*fermi*fermi;
00393     e1      = 0.3*fermi;
00394     e2      = 0.35*fermi;
00395   }
00396   else // as proton, if were not defined 
00397   {
00398     diffuse = 0.63*fermi;
00399     gamma   = 0.3*fermi;
00400     delta   = 0.1*fermi*fermi;
00401     e1      = 0.3*fermi;
00402     e2      = 0.35*fermi;
00403   }
00404   G4double lambda = 15.; // 15 ok
00405 
00406   //  G4double kgamma    = fWaveVector*gamma;   // wavek*delta;
00407 
00408   G4double kgamma    = lambda*(1.-std::exp(-fWaveVector*gamma/lambda));   // wavek*delta;
00409   G4double kgamma2   = kgamma*kgamma;
00410 
00411   // G4double dk2t  = delta*fWaveVector*fWaveVector*theta; // delta*wavek*wavek*theta;
00412   // G4double dk2t2 = dk2t*dk2t;
00413   // G4double pikdt = pi*fWaveVector*diffuse*theta;// pi*wavek*diffuse*theta;
00414 
00415   G4double pikdt    = lambda*(1.-std::exp(-pi*fWaveVector*diffuse*theta/lambda));   // wavek*delta;
00416 
00417   damp           = DampFactor(pikdt);
00418   damp2          = damp*damp;
00419 
00420   G4double mode2k2 = (e1*e1+e2*e2)*fWaveVector*fWaveVector;  
00421   G4double e2dk3t  = -2.*e2*delta*fWaveVector*fWaveVector*fWaveVector*theta;
00422 
00423 
00424   sigma  = kgamma2;
00425   // sigma  += dk2t2;
00426   sigma *= bzero2;
00427   sigma += mode2k2*bone2 + e2dk3t*bzero*bone;
00428   sigma += kr2*bonebyarg2;
00429   sigma *= damp2;          // *rad*rad;
00430 
00431   return sigma;
00432 }

G4double G4DiffuseElastic::GetDiffElasticSumProb ( G4double  theta  ) 

Definition at line 440 of file G4DiffuseElastic.cc.

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

Referenced by GetDiffuseElasticSumXsc().

00445 {
00446   G4double sigma, bzero, bzero2, bonebyarg, bonebyarg2, damp, damp2;
00447   G4double delta, diffuse, gamma;
00448   G4double e1, e2, bone, bone2;
00449 
00450   // G4double wavek = momentum/hbarc;  // wave vector
00451   // G4double r0    = 1.08*fermi;
00452   // G4double rad   = r0*std::pow(A, 1./3.);
00453 
00454   G4double kr    = fWaveVector*fNuclearRadius; // wavek*rad;
00455   G4double kr2   = kr*kr;
00456   G4double krt   = kr*theta;
00457 
00458   bzero      = BesselJzero(krt);
00459   bzero2     = bzero*bzero;    
00460   bone       = BesselJone(krt);
00461   bone2      = bone*bone;
00462   bonebyarg  = BesselOneByArg(krt);
00463   bonebyarg2 = bonebyarg*bonebyarg;  
00464 
00465   if (fParticle == theProton)
00466   {
00467     diffuse = 0.63*fermi;
00468     // diffuse = 0.6*fermi;
00469     gamma   = 0.3*fermi;
00470     delta   = 0.1*fermi*fermi;
00471     e1      = 0.3*fermi;
00472     e2      = 0.35*fermi;
00473   }
00474   else // as proton, if were not defined 
00475   {
00476     diffuse = 0.63*fermi;
00477     gamma   = 0.3*fermi;
00478     delta   = 0.1*fermi*fermi;
00479     e1      = 0.3*fermi;
00480     e2      = 0.35*fermi;
00481   }
00482   G4double lambda = 15.; // 15 ok
00483   // G4double kgamma    = fWaveVector*gamma;   // wavek*delta;
00484   G4double kgamma    = lambda*(1.-std::exp(-fWaveVector*gamma/lambda));   // wavek*delta;
00485 
00486   // G4cout<<"kgamma = "<<kgamma<<G4endl;
00487 
00488   if(fAddCoulomb)  // add Coulomb correction
00489   {
00490     G4double sinHalfTheta  = std::sin(0.5*theta);
00491     G4double sinHalfTheta2 = sinHalfTheta*sinHalfTheta;
00492 
00493     kgamma += 0.5*fZommerfeld/kr/(sinHalfTheta2+fAm); // correction at J0()
00494   // kgamma += 0.65*fZommerfeld/kr/(sinHalfTheta2+fAm); // correction at J0()
00495   }
00496 
00497   G4double kgamma2   = kgamma*kgamma;
00498 
00499   // G4double dk2t  = delta*fWaveVector*fWaveVector*theta; // delta*wavek*wavek*theta;
00500   //   G4cout<<"dk2t = "<<dk2t<<G4endl;
00501   // G4double dk2t2 = dk2t*dk2t;
00502   // G4double pikdt = pi*fWaveVector*diffuse*theta;// pi*wavek*diffuse*theta;
00503 
00504   G4double pikdt    = lambda*(1.-std::exp(-pi*fWaveVector*diffuse*theta/lambda));   // wavek*delta;
00505 
00506   // G4cout<<"pikdt = "<<pikdt<<G4endl;
00507 
00508   damp           = DampFactor(pikdt);
00509   damp2          = damp*damp;
00510 
00511   G4double mode2k2 = (e1*e1+e2*e2)*fWaveVector*fWaveVector;  
00512   G4double e2dk3t  = -2.*e2*delta*fWaveVector*fWaveVector*fWaveVector*theta;
00513 
00514   sigma  = kgamma2;
00515   // sigma += dk2t2;
00516   sigma *= bzero2;
00517   sigma += mode2k2*bone2; 
00518   sigma += e2dk3t*bzero*bone;
00519 
00520   // sigma += kr2*(1 + 8.*fZommerfeld*fZommerfeld/kr2)*bonebyarg2;  // correction at J1()/()
00521   sigma += kr2*bonebyarg2;  // correction at J1()/()
00522 
00523   sigma *= damp2;          // *rad*rad;
00524 
00525   return sigma;
00526 }

G4double G4DiffuseElastic::GetDiffElasticSumProbA ( G4double  alpha  ) 

Definition at line 535 of file G4DiffuseElastic.cc.

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

Referenced by GetIntegrandFunction().

00536 {
00537   G4double theta; 
00538 
00539   theta = std::sqrt(alpha);
00540 
00541   // theta = std::acos( 1 - alpha/2. );
00542 
00543   G4double sigma, bzero, bzero2, bonebyarg, bonebyarg2, damp, damp2;
00544   G4double delta, diffuse, gamma;
00545   G4double e1, e2, bone, bone2;
00546 
00547   // G4double wavek = momentum/hbarc;  // wave vector
00548   // G4double r0    = 1.08*fermi;
00549   // G4double rad   = r0*std::pow(A, 1./3.);
00550 
00551   G4double kr    = fWaveVector*fNuclearRadius; // wavek*rad;
00552   G4double kr2   = kr*kr;
00553   G4double krt   = kr*theta;
00554 
00555   bzero      = BesselJzero(krt);
00556   bzero2     = bzero*bzero;    
00557   bone       = BesselJone(krt);
00558   bone2      = bone*bone;
00559   bonebyarg  = BesselOneByArg(krt);
00560   bonebyarg2 = bonebyarg*bonebyarg;  
00561 
00562   if (fParticle == theProton)
00563   {
00564     diffuse = 0.63*fermi;
00565     // diffuse = 0.6*fermi;
00566     gamma   = 0.3*fermi;
00567     delta   = 0.1*fermi*fermi;
00568     e1      = 0.3*fermi;
00569     e2      = 0.35*fermi;
00570   }
00571   else // as proton, if were not defined 
00572   {
00573     diffuse = 0.63*fermi;
00574     gamma   = 0.3*fermi;
00575     delta   = 0.1*fermi*fermi;
00576     e1      = 0.3*fermi;
00577     e2      = 0.35*fermi;
00578   }
00579   G4double lambda = 15.; // 15 ok
00580   // G4double kgamma    = fWaveVector*gamma;   // wavek*delta;
00581   G4double kgamma    = lambda*(1.-std::exp(-fWaveVector*gamma/lambda));   // wavek*delta;
00582 
00583   // G4cout<<"kgamma = "<<kgamma<<G4endl;
00584 
00585   if(fAddCoulomb)  // add Coulomb correction
00586   {
00587     G4double sinHalfTheta  = theta*0.5; // std::sin(0.5*theta);
00588     G4double sinHalfTheta2 = sinHalfTheta*sinHalfTheta;
00589 
00590     kgamma += 0.5*fZommerfeld/kr/(sinHalfTheta2+fAm); // correction at J0()
00591   // kgamma += 0.65*fZommerfeld/kr/(sinHalfTheta2+fAm); // correction at J0()
00592   }
00593 
00594   G4double kgamma2   = kgamma*kgamma;
00595 
00596   // G4double dk2t  = delta*fWaveVector*fWaveVector*theta; // delta*wavek*wavek*theta;
00597   //   G4cout<<"dk2t = "<<dk2t<<G4endl;
00598   // G4double dk2t2 = dk2t*dk2t;
00599   // G4double pikdt = pi*fWaveVector*diffuse*theta;// pi*wavek*diffuse*theta;
00600 
00601   G4double pikdt    = lambda*(1.-std::exp(-pi*fWaveVector*diffuse*theta/lambda));   // wavek*delta;
00602 
00603   // G4cout<<"pikdt = "<<pikdt<<G4endl;
00604 
00605   damp           = DampFactor(pikdt);
00606   damp2          = damp*damp;
00607 
00608   G4double mode2k2 = (e1*e1+e2*e2)*fWaveVector*fWaveVector;  
00609   G4double e2dk3t  = -2.*e2*delta*fWaveVector*fWaveVector*fWaveVector*theta;
00610 
00611   sigma  = kgamma2;
00612   // sigma += dk2t2;
00613   sigma *= bzero2;
00614   sigma += mode2k2*bone2; 
00615   sigma += e2dk3t*bzero*bone;
00616 
00617   // sigma += kr2*(1 + 8.*fZommerfeld*fZommerfeld/kr2)*bonebyarg2;  // correction at J1()/()
00618   sigma += kr2*bonebyarg2;  // correction at J1()/()
00619 
00620   sigma *= damp2;          // *rad*rad;
00621 
00622   return sigma;
00623 }

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

Definition at line 227 of file G4DiffuseElastic.cc.

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

Referenced by GetInvElasticSumXsc().

00231 {
00232   fParticle      = particle;
00233   fWaveVector    = momentum/hbarc;
00234   fAtomicWeight  = A;
00235   fAtomicNumber  = Z;
00236   fNuclearRadius = CalculateNuclearRad(A);
00237   fAddCoulomb    = false;
00238 
00239   G4double z     = particle->GetPDGCharge();
00240 
00241   G4double kRt   = fWaveVector*fNuclearRadius*theta;
00242   G4double kRtC  = 1.9;
00243 
00244   if( z && (kRt > kRtC) )
00245   {
00246     fAddCoulomb = true;
00247     fBeta       = CalculateParticleBeta( particle, momentum);
00248     fZommerfeld = CalculateZommerfeld( fBeta, z, fAtomicNumber);
00249     fAm         = CalculateAm( momentum, fZommerfeld, fAtomicNumber);
00250   }
00251   G4double sigma = fNuclearRadius*fNuclearRadius*GetDiffElasticSumProb(theta);
00252 
00253   return sigma;
00254 }

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

Definition at line 156 of file G4DiffuseElastic.cc.

References CalculateNuclearRad(), and GetDiffElasticProb().

Referenced by GetInvElasticXsc().

00160 {
00161   fParticle      = particle;
00162   fWaveVector    = momentum/hbarc;
00163   fAtomicWeight  = A;
00164   fAddCoulomb    = false;
00165   fNuclearRadius = CalculateNuclearRad(A);
00166 
00167   G4double sigma = fNuclearRadius*fNuclearRadius*GetDiffElasticProb(theta);
00168 
00169   return sigma;
00170 }

G4double G4DiffuseElastic::GetIntegrandFunction ( G4double  theta  ) 

Definition at line 631 of file G4DiffuseElastic.cc.

References GetDiffElasticSumProbA().

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

00632 {
00633   G4double result;
00634 
00635   result  = GetDiffElasticSumProbA(alpha);
00636 
00637   // result *= 2*pi*std::sin(theta);
00638 
00639   return result;
00640 }

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

Definition at line 314 of file G4DiffuseElastic.cc.

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

00318 {
00319   G4double m1 = particle->GetPDGMass();
00320   G4LorentzVector lv1(0.,0.,plab,std::sqrt(plab*plab+m1*m1));
00321 
00322   G4int iZ = static_cast<G4int>(Z+0.5);
00323   G4int iA = static_cast<G4int>(A+0.5);
00324   G4ParticleDefinition * theDef = 0;
00325 
00326   if      (iZ == 1 && iA == 1) theDef = theProton;
00327   else if (iZ == 1 && iA == 2) theDef = theDeuteron;
00328   else if (iZ == 1 && iA == 3) theDef = G4Triton::Triton();
00329   else if (iZ == 2 && iA == 3) theDef = G4He3::He3();
00330   else if (iZ == 2 && iA == 4) theDef = theAlpha;
00331   else theDef = G4ParticleTable::GetParticleTable()->FindIon(iZ,iA,0,iZ);
00332  
00333   G4double tmass = theDef->GetPDGMass();
00334 
00335   G4LorentzVector lv(0.0,0.0,0.0,tmass);   
00336   lv += lv1;
00337 
00338   G4ThreeVector bst = lv.boostVector();
00339   lv1.boost(-bst);
00340 
00341   G4ThreeVector p1 = lv1.vect();
00342   G4double ptot    = p1.mag();
00343   G4double ptot2 = ptot*ptot;
00344   G4double cost = 1 - 0.5*std::fabs(tMand)/ptot2;
00345 
00346   if( cost >= 1.0 )      cost = 1.0;  
00347   else if( cost <= -1.0) cost = -1.0;
00348   
00349   G4double thetaCMS = std::acos(cost);
00350 
00351   G4double sigma = GetCoulombElasticXsc( particle, thetaCMS, ptot, Z );
00352 
00353   sigma *= pi/ptot2;
00354 
00355   return sigma;
00356 }

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

Definition at line 262 of file G4DiffuseElastic.cc.

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

00266 {
00267   G4double m1 = particle->GetPDGMass();
00268 
00269   G4LorentzVector lv1(0.,0.,plab,std::sqrt(plab*plab+m1*m1));
00270 
00271   G4int iZ = static_cast<G4int>(Z+0.5);
00272   G4int iA = static_cast<G4int>(A+0.5);
00273 
00274   G4ParticleDefinition* theDef = 0;
00275 
00276   if      (iZ == 1 && iA == 1) theDef = theProton;
00277   else if (iZ == 1 && iA == 2) theDef = theDeuteron;
00278   else if (iZ == 1 && iA == 3) theDef = G4Triton::Triton();
00279   else if (iZ == 2 && iA == 3) theDef = G4He3::He3();
00280   else if (iZ == 2 && iA == 4) theDef = theAlpha;
00281   else theDef = G4ParticleTable::GetParticleTable()->FindIon(iZ,iA,0,iZ);
00282  
00283   G4double tmass = theDef->GetPDGMass();
00284 
00285   G4LorentzVector lv(0.0,0.0,0.0,tmass);   
00286   lv += lv1;
00287 
00288   G4ThreeVector bst = lv.boostVector();
00289   lv1.boost(-bst);
00290 
00291   G4ThreeVector p1 = lv1.vect();
00292   G4double ptot    = p1.mag();
00293   G4double ptot2   = ptot*ptot;
00294   G4double cost    = 1 - 0.5*std::fabs(tMand)/ptot2;
00295 
00296   if( cost >= 1.0 )      cost = 1.0;  
00297   else if( cost <= -1.0) cost = -1.0;
00298   
00299   G4double thetaCMS = std::acos(cost);
00300 
00301   G4double sigma = GetDiffuseElasticSumXsc( particle, thetaCMS, ptot, A, Z );
00302 
00303   sigma *= pi/ptot2;
00304 
00305   return sigma;
00306 }

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

Definition at line 177 of file G4DiffuseElastic.cc.

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

00181 {
00182   G4double m1 = particle->GetPDGMass();
00183   G4LorentzVector lv1(0.,0.,plab,std::sqrt(plab*plab+m1*m1));
00184 
00185   G4int iZ = static_cast<G4int>(Z+0.5);
00186   G4int iA = static_cast<G4int>(A+0.5);
00187   G4ParticleDefinition * theDef = 0;
00188 
00189   if      (iZ == 1 && iA == 1) theDef = theProton;
00190   else if (iZ == 1 && iA == 2) theDef = theDeuteron;
00191   else if (iZ == 1 && iA == 3) theDef = G4Triton::Triton();
00192   else if (iZ == 2 && iA == 3) theDef = G4He3::He3();
00193   else if (iZ == 2 && iA == 4) theDef = theAlpha;
00194   else theDef = G4ParticleTable::GetParticleTable()->FindIon(iZ,iA,0,iZ);
00195  
00196   G4double tmass = theDef->GetPDGMass();
00197 
00198   G4LorentzVector lv(0.0,0.0,0.0,tmass);   
00199   lv += lv1;
00200 
00201   G4ThreeVector bst = lv.boostVector();
00202   lv1.boost(-bst);
00203 
00204   G4ThreeVector p1 = lv1.vect();
00205   G4double ptot    = p1.mag();
00206   G4double ptot2 = ptot*ptot;
00207   G4double cost = 1 - 0.5*std::fabs(tMand)/ptot2;
00208 
00209   if( cost >= 1.0 )      cost = 1.0;  
00210   else if( cost <= -1.0) cost = -1.0;
00211   
00212   G4double thetaCMS = std::acos(cost);
00213 
00214   G4double sigma = GetDiffuseElasticXsc( particle, thetaCMS, ptot, A);
00215 
00216   sigma *= pi/ptot2;
00217 
00218   return sigma;
00219 }

G4double G4DiffuseElastic::GetNuclearRadius (  )  [inline]

Definition at line 186 of file G4DiffuseElastic.hh.

00186 {return fNuclearRadius;};

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

Definition at line 1007 of file G4DiffuseElastic.cc.

References G4UniformRand.

Referenced by SampleTableThetaCMS().

01008 {
01009  G4double x1, x2, y1, y2, randAngle;
01010 
01011   if( iAngle == 0 )
01012   {
01013     randAngle = (*fAngleTable)(iMomentum)->GetLowEdgeEnergy(iAngle);
01014     // iAngle++;
01015   }
01016   else
01017   {
01018     if ( iAngle >= G4int((*fAngleTable)(iMomentum)->GetVectorLength()) )
01019     {
01020       iAngle = (*fAngleTable)(iMomentum)->GetVectorLength() - 1;
01021     }
01022     y1 = (*(*fAngleTable)(iMomentum))(iAngle-1);
01023     y2 = (*(*fAngleTable)(iMomentum))(iAngle);
01024 
01025     x1 = (*fAngleTable)(iMomentum)->GetLowEdgeEnergy(iAngle-1);
01026     x2 = (*fAngleTable)(iMomentum)->GetLowEdgeEnergy(iAngle);
01027 
01028     if ( x1 == x2 )   randAngle = x2;
01029     else
01030     {
01031       if ( y1 == y2 ) randAngle = x1 + ( x2 - x1 )*G4UniformRand();
01032       else
01033       {
01034         randAngle = x1 + ( position - y1 )*( x2 - x1 )/( y2 - y1 );
01035       }
01036     }
01037   }
01038   return randAngle;
01039 }

void G4DiffuseElastic::Initialise (  ) 

Definition at line 122 of file G4DiffuseElastic.cc.

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

00123 {
00124 
00125   // fEnergyVector = new G4PhysicsLogVector( theMinEnergy, theMaxEnergy, fEnergyBin );
00126 
00127   const G4ElementTable* theElementTable = G4Element::GetElementTable();
00128 
00129   size_t jEl, numOfEl = G4Element::GetNumberOfElements();
00130 
00131   for(jEl = 0 ; jEl < numOfEl; ++jEl) // application element loop
00132   {
00133     fAtomicNumber = (*theElementTable)[jEl]->GetZ();     // atomic number
00134     fAtomicWeight = (*theElementTable)[jEl]->GetN();     // number of nucleons
00135     fNuclearRadius = CalculateNuclearRad(fAtomicWeight);
00136 
00137     if(verboseLevel > 0) 
00138     {   
00139       G4cout<<"G4DiffuseElastic::Initialise() the element: "
00140             <<(*theElementTable)[jEl]->GetName()<<G4endl;
00141     }
00142     fElementNumberVector.push_back(fAtomicNumber);
00143     fElementNameVector.push_back((*theElementTable)[jEl]->GetName());
00144 
00145     BuildAngleTable();
00146     fAngleBank.push_back(fAngleTable);
00147   }  
00148   return;
00149 }

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

Definition at line 896 of file G4DiffuseElastic.cc.

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

Referenced by SampleTableThetaCMS().

00897 {
00898   fAtomicNumber  = Z;     // atomic number
00899   fAtomicWeight  = A;     // number of nucleons
00900 
00901   fNuclearRadius = CalculateNuclearRad(fAtomicWeight);
00902   
00903   if( verboseLevel > 0 )    
00904   {
00905     G4cout<<"G4DiffuseElastic::Initialise() the element with Z = "
00906           <<Z<<"; and A = "<<A<<G4endl;
00907   }
00908   fElementNumberVector.push_back(fAtomicNumber);
00909 
00910   BuildAngleTable();
00911 
00912   fAngleBank.push_back(fAngleTable);
00913 
00914   return;
00915 }

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

Definition at line 647 of file G4DiffuseElastic.cc.

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

00651 {
00652   G4double result;
00653   fParticle      = particle;
00654   fWaveVector    = momentum/hbarc;
00655   fAtomicWeight  = A;
00656 
00657   fNuclearRadius = CalculateNuclearRad(A);
00658 
00659 
00660   G4Integrator<G4DiffuseElastic,G4double(G4DiffuseElastic::*)(G4double)> integral;
00661 
00662   // result = integral.Legendre10(this,&G4DiffuseElastic::GetIntegrandFunction, 0., theta ); 
00663   result = integral.Legendre96(this,&G4DiffuseElastic::GetIntegrandFunction, 0., theta ); 
00664 
00665   return result;
00666 }

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

Reimplemented from G4HadronElastic.

Definition at line 738 of file G4DiffuseElastic.cc.

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

00740 {
00741   fParticle = aParticle;
00742   G4double m1 = fParticle->GetPDGMass();
00743   G4double totElab = std::sqrt(m1*m1+p*p);
00744   G4double mass2 = G4NucleiProperties::GetNuclearMass(A, Z);
00745   G4LorentzVector lv1(p,0.0,0.0,totElab);
00746   G4LorentzVector  lv(0.0,0.0,0.0,mass2);   
00747   lv += lv1;
00748 
00749   G4ThreeVector bst = lv.boostVector();
00750   lv1.boost(-bst);
00751 
00752   G4ThreeVector p1 = lv1.vect();
00753   G4double momentumCMS = p1.mag();
00754 
00755   G4double t = SampleTableT( aParticle,  momentumCMS, G4double(Z), G4double(A) ); // sample theta2 in cms
00756   return t;
00757 }

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

Definition at line 672 of file G4DiffuseElastic.cc.

References SampleThetaCMS().

Referenced by SampleThetaLab().

00673 {
00674   G4double theta = SampleThetaCMS( aParticle,  p, A); // sample theta in cms
00675   G4double t     = 2*p*p*( 1 - std::cos(theta) ); // -t !!!
00676   return t;
00677 }

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

Definition at line 763 of file G4DiffuseElastic.cc.

References G4InuclParticleNames::alpha, and SampleTableThetaCMS().

Referenced by SampleInvariantT().

00765 {
00766   G4double alpha = SampleTableThetaCMS( aParticle,  p, Z, A); // sample theta2 in cms
00767   // G4double t     = 2*p*p*( 1 - std::cos(std::sqrt(alpha)) );             // -t !!!
00768   G4double t     = p*p*alpha;             // -t !!!
00769   return t;
00770 }

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

Definition at line 778 of file G4DiffuseElastic.cc.

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

Referenced by SampleTableT().

00780 {
00781   size_t iElement;
00782   G4int iMomentum, iAngle;  
00783   G4double randAngle, position, theta1, theta2, E1, E2, W1, W2, W;  
00784   G4double m1 = particle->GetPDGMass();
00785 
00786   for(iElement = 0; iElement < fElementNumberVector.size(); iElement++)
00787   {
00788     if( std::fabs(Z - fElementNumberVector[iElement]) < 0.5) break;
00789   }
00790   if ( iElement == fElementNumberVector.size() ) 
00791   {
00792     InitialiseOnFly(Z,A); // table preparation, if needed
00793 
00794     // iElement--;
00795 
00796     // G4cout << "G4DiffuseElastic: Element with atomic number " << Z
00797     // << " is not found, return zero angle" << G4endl;
00798     // return 0.; // no table for this element
00799   }
00800   // G4cout<<"iElement = "<<iElement<<G4endl;
00801 
00802   fAngleTable = fAngleBank[iElement];
00803 
00804   G4double kinE = std::sqrt(momentum*momentum + m1*m1) - m1;
00805 
00806   for( iMomentum = 0; iMomentum < fEnergyBin; iMomentum++)
00807   {
00808     if( kinE < fEnergyVector->GetLowEdgeEnergy(iMomentum) ) break;
00809   }
00810   if ( iMomentum >= fEnergyBin ) iMomentum = fEnergyBin-1;   // kinE is more then theMaxEnergy
00811   if ( iMomentum < 0 )           iMomentum = 0; // against negative index, kinE < theMinEnergy
00812 
00813   // G4cout<<"iMomentum = "<<iMomentum<<G4endl;
00814 
00815   if (iMomentum == fEnergyBin -1 || iMomentum == 0 )   // the table edges
00816   {
00817     position = (*(*fAngleTable)(iMomentum))(fAngleBin-2)*G4UniformRand();
00818 
00819     // G4cout<<"position = "<<position<<G4endl;
00820 
00821     for(iAngle = 0; iAngle < fAngleBin-1; iAngle++)
00822     {
00823       if( position < (*(*fAngleTable)(iMomentum))(iAngle) ) break;
00824     }
00825     if (iAngle >= fAngleBin-1) iAngle = fAngleBin-2;
00826 
00827     // G4cout<<"iAngle = "<<iAngle<<G4endl;
00828 
00829     randAngle = GetScatteringAngle(iMomentum, iAngle, position);
00830 
00831     // G4cout<<"randAngle = "<<randAngle<<G4endl;
00832   }
00833   else  // kinE inside between energy table edges
00834   {
00835     // position = (*(*fAngleTable)(iMomentum))(fAngleBin-2)*G4UniformRand();
00836     position = (*(*fAngleTable)(iMomentum))(0)*G4UniformRand();
00837 
00838     // G4cout<<"position = "<<position<<G4endl;
00839 
00840     for(iAngle = 0; iAngle < fAngleBin-1; iAngle++)
00841     {
00842       // if( position < (*(*fAngleTable)(iMomentum))(iAngle) ) break;
00843       if( position > (*(*fAngleTable)(iMomentum))(iAngle) ) break;
00844     }
00845     if (iAngle >= fAngleBin-1) iAngle = fAngleBin-2;
00846 
00847     // G4cout<<"iAngle = "<<iAngle<<G4endl;
00848 
00849     theta2  = GetScatteringAngle(iMomentum, iAngle, position);
00850 
00851     // G4cout<<"theta2 = "<<theta2<<G4endl;
00852     E2 = fEnergyVector->GetLowEdgeEnergy(iMomentum);
00853 
00854     // G4cout<<"E2 = "<<E2<<G4endl;
00855     
00856     iMomentum--;
00857     
00858     // position = (*(*fAngleTable)(iMomentum))(fAngleBin-2)*G4UniformRand();
00859 
00860     // G4cout<<"position = "<<position<<G4endl;
00861 
00862     for(iAngle = 0; iAngle < fAngleBin-1; iAngle++)
00863     {
00864       // if( position < (*(*fAngleTable)(iMomentum))(iAngle) ) break;
00865       if( position > (*(*fAngleTable)(iMomentum))(iAngle) ) break;
00866     }
00867     if (iAngle >= fAngleBin-1) iAngle = fAngleBin-2;
00868     
00869     theta1  = GetScatteringAngle(iMomentum, iAngle, position);
00870 
00871     // G4cout<<"theta1 = "<<theta1<<G4endl;
00872 
00873     E1 = fEnergyVector->GetLowEdgeEnergy(iMomentum);
00874 
00875     // G4cout<<"E1 = "<<E1<<G4endl;
00876 
00877     W  = 1.0/(E2 - E1);
00878     W1 = (E2 - kinE)*W;
00879     W2 = (kinE - E1)*W;
00880 
00881     randAngle = W1*theta1 + W2*theta2;
00882     
00883     // randAngle = theta2;
00884     // G4cout<<"randAngle = "<<randAngle<<G4endl;
00885   }
00886   // G4double angle = randAngle;
00887   // if (randAngle > 0.) randAngle /= 2*pi*std::sin(angle);
00888 
00889   return randAngle;
00890 }

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

Definition at line 685 of file G4DiffuseElastic.cc.

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

Referenced by SampleT().

00687 {
00688   G4int i, iMax = 100;  
00689   G4double norm, result, theta1, theta2, thetaMax, sum = 0.;
00690 
00691   fParticle      = particle;
00692   fWaveVector    = momentum/hbarc;
00693   fAtomicWeight  = A;
00694 
00695   fNuclearRadius = CalculateNuclearRad(A);
00696   
00697   thetaMax = 10.174/fWaveVector/fNuclearRadius;
00698 
00699   if (thetaMax > pi) thetaMax = pi;
00700 
00701   G4Integrator<G4DiffuseElastic,G4double(G4DiffuseElastic::*)(G4double)> integral;
00702 
00703   // result = integral.Legendre10(this,&G4DiffuseElastic::GetIntegrandFunction, 0., theta ); 
00704   norm = integral.Legendre96(this,&G4DiffuseElastic::GetIntegrandFunction, 0., thetaMax );
00705 
00706   norm *= G4UniformRand();
00707 
00708   for(i = 1; i <= iMax; i++)
00709   {
00710     theta1 = (i-1)*thetaMax/iMax; 
00711     theta2 = i*thetaMax/iMax;
00712     sum   += integral.Legendre10(this,&G4DiffuseElastic::GetIntegrandFunction, theta1, theta2);
00713 
00714     if ( sum >= norm ) 
00715     {
00716       result = 0.5*(theta1 + theta2);
00717       break;
00718     }
00719   }
00720   if (i > iMax ) result = 0.5*(theta1 + theta2);
00721 
00722   G4double sigma = pi*thetaMax/iMax;
00723 
00724   result += G4RandGauss::shoot(0.,sigma);
00725 
00726   if(result < 0.) result = 0.;
00727   if(result > thetaMax) result = thetaMax;
00728 
00729   return result;
00730 }

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

Definition at line 1050 of file G4DiffuseElastic.cc.

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

01052 {
01053   const G4ParticleDefinition* theParticle = aParticle->GetDefinition();
01054   G4double m1 = theParticle->GetPDGMass();
01055   G4double plab = aParticle->GetTotalMomentum();
01056   G4LorentzVector lv1 = aParticle->Get4Momentum();
01057   G4LorentzVector lv(0.0,0.0,0.0,tmass);   
01058   lv += lv1;
01059 
01060   G4ThreeVector bst = lv.boostVector();
01061   lv1.boost(-bst);
01062 
01063   G4ThreeVector p1 = lv1.vect();
01064   G4double ptot    = p1.mag();
01065   G4double tmax    = 4.0*ptot*ptot;
01066   G4double t       = 0.0;
01067 
01068 
01069   //
01070   // Sample t
01071   //
01072   
01073   t = SampleT( theParticle, ptot, A);
01074 
01075   // NaN finder
01076   if(!(t < 0.0 || t >= 0.0)) 
01077   {
01078     if (verboseLevel > 0) 
01079     {
01080       G4cout << "G4DiffuseElastic:WARNING: A = " << A 
01081              << " mom(GeV)= " << plab/GeV 
01082              << " S-wave will be sampled" 
01083              << G4endl; 
01084     }
01085     t = G4UniformRand()*tmax; 
01086   }
01087   if(verboseLevel>1)
01088   {
01089     G4cout <<" t= " << t << " tmax= " << tmax 
01090            << " ptot= " << ptot << G4endl;
01091   }
01092   // Sampling of angles in CM system
01093 
01094   G4double phi  = G4UniformRand()*twopi;
01095   G4double cost = 1. - 2.0*t/tmax;
01096   G4double sint;
01097 
01098   if( cost >= 1.0 ) 
01099   {
01100     cost = 1.0;
01101     sint = 0.0;
01102   }
01103   else if( cost <= -1.0) 
01104   {
01105     cost = -1.0;
01106     sint =  0.0;
01107   }
01108   else  
01109   {
01110     sint = std::sqrt((1.0-cost)*(1.0+cost));
01111   }    
01112   if (verboseLevel>1) 
01113   {
01114     G4cout << "cos(t)=" << cost << " std::sin(t)=" << sint << G4endl;
01115   }
01116   G4ThreeVector v1(sint*std::cos(phi),sint*std::sin(phi),cost);
01117   v1 *= ptot;
01118   G4LorentzVector nlv1(v1.x(),v1.y(),v1.z(),std::sqrt(ptot*ptot + m1*m1));
01119 
01120   nlv1.boost(bst); 
01121 
01122   G4ThreeVector np1 = nlv1.vect();
01123 
01124     // G4double theta = std::acos( np1.z()/np1.mag() );  // degree;
01125 
01126   G4double theta = np1.theta();
01127 
01128   return theta;
01129 }

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

Definition at line 238 of file G4DiffuseElastic.hh.

00239 {
00240   lowEnergyLimitHE = value;
00241 }

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

Reimplemented from G4HadronElastic.

Definition at line 248 of file G4DiffuseElastic.hh.

00249 {
00250   lowestEnergyLimit = value;
00251 }

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

Definition at line 233 of file G4DiffuseElastic.hh.

00234 {
00235   plabLowLimit = value;
00236 }

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

Definition at line 243 of file G4DiffuseElastic.hh.

00244 {
00245   lowEnergyLimitQ = value;
00246 }

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

Definition at line 228 of file G4DiffuseElastic.hh.

00229 {
00230   lowEnergyRecoilLimit = value;
00231 }

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

Definition at line 1257 of file G4DiffuseElastic.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().

01259 {
01260   fAtomicNumber  = Z;     // atomic number
01261   fAtomicWeight  = A;     // number of nucleons
01262   fNuclearRadius = CalculateNuclearRad(fAtomicWeight);
01263   
01264      
01265   
01266   G4cout<<"G4DiffuseElastic::TestAngleTable() init the element with Z = "
01267           <<Z<<"; and A = "<<A<<G4endl;
01268  
01269   fElementNumberVector.push_back(fAtomicNumber);
01270 
01271  
01272 
01273 
01274   G4int i=0, j;
01275   G4double a = 0., z = theParticle->GetPDGCharge(),  m1 = fParticle->GetPDGMass();
01276   G4double alpha1=0., alpha2=0., alphaMax=0., alphaCoulomb=0.;
01277   G4double deltaL10 = 0., deltaL96 = 0., deltaAG = 0.;
01278   G4double sumL10 = 0.,sumL96 = 0.,sumAG = 0.;
01279   G4double epsilon = 0.001;
01280 
01281   G4Integrator<G4DiffuseElastic,G4double(G4DiffuseElastic::*)(G4double)> integral;
01282   
01283   fAngleTable = new G4PhysicsTable(fEnergyBin);
01284 
01285   fWaveVector = partMom/hbarc;
01286 
01287   G4double kR     = fWaveVector*fNuclearRadius;
01288   G4double kR2    = kR*kR;
01289   G4double kRmax  = 10.6; // 10.6, 18, 10.174; ~ 3 maxima of J1 or 15., 25.
01290   G4double kRcoul = 1.2; // 1.4, 2.5; // on the first slope of J1
01291 
01292   alphaMax = kRmax*kRmax/kR2;
01293 
01294   if (alphaMax > 4.) alphaMax = 4.;  // vmg05-02-09: was pi2 
01295 
01296   alphaCoulomb = kRcoul*kRcoul/kR2;
01297 
01298   if( z )
01299   {
01300       a           = partMom/m1; // beta*gamma for m1
01301       fBeta       = a/std::sqrt(1+a*a);
01302       fZommerfeld = CalculateZommerfeld( fBeta, z, fAtomicNumber);
01303       fAm         = CalculateAm( partMom, fZommerfeld, fAtomicNumber);
01304   }
01305   G4PhysicsFreeVector* angleVector = new G4PhysicsFreeVector(fAngleBin-1);
01306 
01307   // G4PhysicsLogVector*  angleBins = new G4PhysicsLogVector( 0.001*alphaMax, alphaMax, fAngleBin );
01308     
01309   
01310   fAddCoulomb = false;
01311 
01312   for(j = 1; j < fAngleBin; j++)
01313   {
01314       // alpha1 = angleBins->GetLowEdgeEnergy(j-1);
01315       // alpha2 = angleBins->GetLowEdgeEnergy(j);
01316 
01317     alpha1 = alphaMax*(j-1)/fAngleBin;
01318     alpha2 = alphaMax*( j )/fAngleBin;
01319 
01320     if( ( alpha2 > alphaCoulomb ) && z ) fAddCoulomb = true;
01321 
01322     deltaL10 = integral.Legendre10(this, &G4DiffuseElastic::GetIntegrandFunction, alpha1, alpha2);
01323     deltaL96 = integral.Legendre96(this, &G4DiffuseElastic::GetIntegrandFunction, alpha1, alpha2);
01324     deltaAG  = integral.AdaptiveGauss(this, &G4DiffuseElastic::GetIntegrandFunction, 
01325                                        alpha1, alpha2,epsilon);
01326 
01327       // G4cout<<alpha1<<"\t"<<std::sqrt(alpha1)/degree<<"\t"
01328       //     <<deltaL10<<"\t"<<deltaL96<<"\t"<<deltaAG<<G4endl;
01329 
01330     sumL10 += deltaL10;
01331     sumL96 += deltaL96;
01332     sumAG  += deltaAG;
01333 
01334     G4cout<<alpha1<<"\t"<<std::sqrt(alpha1)/degree<<"\t"
01335             <<sumL10<<"\t"<<sumL96<<"\t"<<sumAG<<G4endl;
01336       
01337     angleVector->PutValue( j-1 , alpha1, sumL10 ); // alpha2
01338   }
01339   fAngleTable->insertAt(i,angleVector);
01340   fAngleBank.push_back(fAngleTable);
01341 
01342   /*
01343   // Integral over all angle range - Bad accuracy !!!
01344 
01345   sumL10 = integral.Legendre10(this, &G4DiffuseElastic::GetIntegrandFunction, 0., alpha2);
01346   sumL96 = integral.Legendre96(this, &G4DiffuseElastic::GetIntegrandFunction, 0., alpha2);
01347   sumAG  = integral.AdaptiveGauss(this, &G4DiffuseElastic::GetIntegrandFunction, 
01348                                        0., alpha2,epsilon);
01349   G4cout<<G4endl;
01350   G4cout<<alpha2<<"\t"<<std::sqrt(alpha2)/degree<<"\t"
01351             <<sumL10<<"\t"<<sumL96<<"\t"<<sumAG<<G4endl;
01352   */
01353   return;
01354 }

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

Definition at line 1138 of file G4DiffuseElastic.cc.

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

01140 {
01141   const G4ParticleDefinition* theParticle = aParticle->GetDefinition();
01142   G4double m1 = theParticle->GetPDGMass();
01143   // G4double plab = aParticle->GetTotalMomentum();
01144   G4LorentzVector lv1 = aParticle->Get4Momentum();
01145   G4LorentzVector lv(0.0,0.0,0.0,tmass);   
01146 
01147   lv += lv1;
01148 
01149   G4ThreeVector bst = lv.boostVector();
01150 
01151   lv1.boost(-bst);
01152 
01153   G4ThreeVector p1 = lv1.vect();
01154   G4double ptot    = p1.mag();
01155 
01156   G4double phi  = G4UniformRand()*twopi;
01157   G4double cost = std::cos(thetaCMS);
01158   G4double sint;
01159 
01160   if( cost >= 1.0 ) 
01161   {
01162     cost = 1.0;
01163     sint = 0.0;
01164   }
01165   else if( cost <= -1.0) 
01166   {
01167     cost = -1.0;
01168     sint =  0.0;
01169   }
01170   else  
01171   {
01172     sint = std::sqrt((1.0-cost)*(1.0+cost));
01173   }    
01174   if (verboseLevel>1) 
01175   {
01176     G4cout << "cos(tcms)=" << cost << " std::sin(tcms)=" << sint << G4endl;
01177   }
01178   G4ThreeVector v1(sint*std::cos(phi),sint*std::sin(phi),cost);
01179   v1 *= ptot;
01180   G4LorentzVector nlv1(v1.x(),v1.y(),v1.z(),std::sqrt(ptot*ptot + m1*m1));
01181 
01182   nlv1.boost(bst); 
01183 
01184   G4ThreeVector np1 = nlv1.vect();
01185 
01186 
01187   G4double thetaLab = np1.theta();
01188 
01189   return thetaLab;
01190 }

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

Definition at line 1198 of file G4DiffuseElastic.cc.

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

01200 {
01201   const G4ParticleDefinition* theParticle = aParticle->GetDefinition();
01202   G4double m1 = theParticle->GetPDGMass();
01203   G4double plab = aParticle->GetTotalMomentum();
01204   G4LorentzVector lv1 = aParticle->Get4Momentum();
01205   G4LorentzVector lv(0.0,0.0,0.0,tmass);   
01206 
01207   lv += lv1;
01208 
01209   G4ThreeVector bst = lv.boostVector();
01210 
01211   // lv1.boost(-bst);
01212 
01213   // G4ThreeVector p1 = lv1.vect();
01214   // G4double ptot    = p1.mag();
01215 
01216   G4double phi  = G4UniformRand()*twopi;
01217   G4double cost = std::cos(thetaLab);
01218   G4double sint;
01219 
01220   if( cost >= 1.0 ) 
01221   {
01222     cost = 1.0;
01223     sint = 0.0;
01224   }
01225   else if( cost <= -1.0) 
01226   {
01227     cost = -1.0;
01228     sint =  0.0;
01229   }
01230   else  
01231   {
01232     sint = std::sqrt((1.0-cost)*(1.0+cost));
01233   }    
01234   if (verboseLevel>1) 
01235   {
01236     G4cout << "cos(tlab)=" << cost << " std::sin(tlab)=" << sint << G4endl;
01237   }
01238   G4ThreeVector v1(sint*std::cos(phi),sint*std::sin(phi),cost);
01239   v1 *= plab;
01240   G4LorentzVector nlv1(v1.x(),v1.y(),v1.z(),std::sqrt(plab*plab + m1*m1));
01241 
01242   nlv1.boost(-bst); 
01243 
01244   G4ThreeVector np1 = nlv1.vect();
01245 
01246 
01247   G4double thetaCMS = np1.theta();
01248 
01249   return thetaCMS;
01250 }

The documentation for this class was generated from the following files:
Generated on Mon May 27 17:51:46 2013 for Geant4 by  doxygen 1.4.7