Ignore:
Timestamp:
Nov 5, 2010, 3:45:55 PM (14 years ago)
Author:
garnier
Message:

update ti head

File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/source/processes/hadronic/models/coherent_elastic/include/G4NuclNuclDiffuseElastic.hh

    r1337 r1340  
    2525//
    2626//
    27 // $Id: G4NuclNuclDiffuseElastic.hh,v 1.8 2009/04/10 13:22:25 grichine Exp $
    28 // GEANT4 tag $Name: geant4-09-04-beta-01 $
     27// $Id: G4NuclNuclDiffuseElastic.hh,v 1.13 2010/09/28 16:28:58 gcosmo Exp $
     28// GEANT4 tag $Name: geant4-09-03-ref-09 $
    2929//
    3030//
    … …  
    189189
    190190  G4complex GammaLogarithm(G4complex xx);
     191  G4complex GammaLogB2n(G4complex xx);
    191192
    192193  G4double  GetErf(G4double x);
    … …  
    203204  G4complex GetErfInt(G4complex z); // , G4int nMax);
    204205
     206  G4double GetLegendrePol(G4int n, G4double x);
    205207
    206208  G4complex TestErfcComp(G4complex z, G4int nMax);
    … …  
    210212  G4complex CoulombAmplitude(G4double theta);
    211213  void CalculateCoulombPhaseZero();
     214  G4double CalculateCoulombPhase(G4int n);
    212215  void CalculateRutherfordAnglePar();
    213216
    … …  
    225228  G4complex Amplitude(G4double theta);
    226229  G4double  AmplitudeMod2(G4double theta);
     230
     231  G4complex AmplitudeGla(G4double theta);
     232  G4double  AmplitudeGlaMod2(G4double theta);
     233
     234  G4complex AmplitudeGG(G4double theta);
     235  G4double  AmplitudeGGMod2(G4double theta);
     236
    227237  void      InitParameters(const G4ParticleDefinition* theParticle, 
    228238                              G4double partMom, G4double Z, G4double A);
     239
     240  void      InitParametersGla(const G4DynamicParticle* aParticle, 
     241                              G4double partMom, G4double Z, G4double A);
     242
     243  G4double GetHadronNucleonXscNS( G4ParticleDefinition* pParticle,
     244                                                 G4double pTkin,
     245                                  G4ParticleDefinition* tParticle);
     246
     247  G4double CalcMandelstamS( const G4double mp ,
     248                                                       const G4double mt ,
     249                                                    const G4double Plab );
    229250
    230251  G4double GetProfileLambda(){return fProfileLambda;};
    … …  
    233254  void SetProfileDelta(G4double pd) {fProfileDelta = pd;};
    234255  void SetProfileAlpha(G4double pa){fProfileAlpha = pa;};
     256  void SetCofLambda(G4double pa){fCofLambda = pa;};
     257  void SetCofAlpha(G4double pa){fCofAlpha = pa;};
     258  void SetCofDelta(G4double pa){fCofDelta = pa;};
     259  void SetCofPhase(G4double pa){fCofPhase = pa;};
     260  void SetCofFar(G4double pa){fCofFar = pa;};
     261  void SetEtaRatio(G4double pa){fEtaRatio = pa;};
     262  void SetMaxL(G4int l){fMaxL = l;};
    235263
    236264private:
    … …  
    269297  G4double fNuclearRadius2;
    270298  G4double fNuclearRadius;
     299  G4double fNuclearRadiusSquare;
    271300
    272301  G4double fBeta;
    … …  
    277306  G4double fCoulombPhase0;
    278307  G4double fHalfRutThetaTg;
     308  G4double fHalfRutThetaTg2;
    279309  G4double fRutherfordTheta;
    280310
    … …  
    282312  G4double fProfileDelta;
    283313  G4double fProfileAlpha;
     314
     315  G4double fCofLambda;
     316  G4double fCofAlpha;
     317  G4double fCofDelta;
     318  G4double fCofPhase;
     319  G4double fCofFar;
     320
     321  G4int    fMaxL;
     322  G4double fSumSigma;
     323  G4double fEtaRatio;
    284324
    285325  G4double fReZ;
    … …  
    631671}
    632672
     673///////////////////////////////////////////////////////////////////
     674//
     675// For the calculation of arg Gamma(z) one needs complex extension
     676// of ln(Gamma(z)) here is approximate algorithm
     677
     678inline G4complex G4NuclNuclDiffuseElastic::GammaLogB2n(G4complex z)
     679{
     680  G4complex z1 = 12.*z;
     681  G4complex z2 = z*z;
     682  G4complex z3 = z2*z;
     683  G4complex z5 = z2*z3;
     684  G4complex z7 = z2*z5;
     685
     686  z3 *= 360.;
     687  z5 *= 1260.;
     688  z7 *= 1680.;
     689
     690  G4complex result  = (z-0.5)*std::log(z)-z+0.5*std::log(twopi);
     691            result += 1./z1 - 1./z3 +1./z5 -1./z7;
     692  return result;
     693}
     694
    633695/////////////////////////////////////////////////////////////////
    634696//
    … …  
    681743  return erfcz;
    682744}
     745
     746inline  G4double G4NuclNuclDiffuseElastic::GetLegendrePol(G4int n, G4double theta)
     747{
     748  G4double legPol, epsilon = 1.e-6;
     749  G4double x = std::cos(theta);
     750
     751  if     ( n  < 0 ) legPol = 0.;
     752  else if( n == 0 ) legPol = 1.;
     753  else if( n == 1 ) legPol = x;
     754  else if( n == 2 ) legPol = (3.*x*x-1.)/2.;
     755  else if( n == 3 ) legPol = (5.*x*x*x-3.*x)/2.;
     756  else if( n == 4 ) legPol = (35.*x*x*x*x-30.*x*x+3.)/8.;
     757  else if( n == 5 ) legPol = (63.*x*x*x*x*x-70.*x*x*x+15.*x)/8.;
     758  else if( n == 6 ) legPol = (231.*x*x*x*x*x*x-315.*x*x*x*x+105.*x*x-5.)/16.;
     759  else           
     760  {
     761    // legPol = ( (2*n-1)*x*GetLegendrePol(n-1,x) - (n-1)*GetLegendrePol(n-2,x) )/n;
     762
     763    legPol = std::sqrt( 2./(n*pi*std::sin(theta+epsilon)) )*std::sin( (n+0.5)*theta+0.25*pi );
     764  }
     765  return legPol;
     766}
     767
     768
    683769
    684770/////////////////////////////////////////////////////////////////
    … …  
    864950  G4double sinHalfTheta2 = sinHalfTheta*sinHalfTheta;
    865951  sinHalfTheta2         += fAm;
     952
    866953  G4double order         = 2.*fCoulombPhase0 - fZommerfeld*std::log(sinHalfTheta2);
    867954  G4complex z            = G4complex(0., order);
    868955  ca                     = std::exp(z);
     956
    869957  ca                    *= -fZommerfeld/(2.*fWaveVector*sinHalfTheta2);
    870958
    … …  
    880968{
    881969  G4complex z        = G4complex(1,fZommerfeld);
    882   G4complex gammalog = GammaLogarithm(z);
    883   fCoulombPhase0      = gammalog.imag();
     970  // G4complex gammalog = GammaLogarithm(z);
     971  G4complex gammalog = GammaLogB2n(z);
     972  fCoulombPhase0     = gammalog.imag();
     973}
     974
     975/////////////////////////////////////////////////////////////////
     976//
     977//
     978
     979
     980inline  G4double G4NuclNuclDiffuseElastic::CalculateCoulombPhase(G4int n)
     981{
     982  G4complex z        = G4complex(1. + n, fZommerfeld);
     983  // G4complex gammalog = GammaLogarithm(z);
     984  G4complex gammalog = GammaLogB2n(z);
     985  return gammalog.imag();
    884986}
    885987
    … …  
    892994inline  void G4NuclNuclDiffuseElastic::CalculateRutherfordAnglePar()
    893995{
    894   fHalfRutThetaTg  = fZommerfeld/(fWaveVector*fNuclearRadius);
    895   fRutherfordTheta = 2.*std::atan(fHalfRutThetaTg);
     996  fHalfRutThetaTg   = fZommerfeld/fProfileLambda;  // (fWaveVector*fNuclearRadius);
     997  fRutherfordTheta  = 2.*std::atan(fHalfRutThetaTg);
     998  fHalfRutThetaTg2  = fHalfRutThetaTg*fHalfRutThetaTg;
    896999  G4cout<<"fRutherfordTheta = "<<fRutherfordTheta/degree<<" degree"<<G4endl;
    8971000
    … …  
    9421045{
    9431046  G4double twosigma = 2.*fCoulombPhase0;
    944   twosigma -= fZommerfeld*std::log(fHalfRutThetaTg/(1.+fHalfRutThetaTg*fHalfRutThetaTg));
     1047  twosigma -= fZommerfeld*std::log(fHalfRutThetaTg2/(1.+fHalfRutThetaTg2));
    9451048  twosigma += fRutherfordTheta*fZommerfeld/fHalfRutThetaTg - halfpi;
    9461049  twosigma -= fProfileLambda*theta - 0.25*pi;
    9471050
     1051  twosigma *= fCofPhase;
     1052
    9481053  G4complex z = G4complex(0., twosigma);
    9491054
    … …  
    9581063{
    9591064  G4double twosigma = 2.*fCoulombPhase0;
    960   twosigma -= fZommerfeld*std::log(fHalfRutThetaTg/(1.+fHalfRutThetaTg*fHalfRutThetaTg));
     1065  twosigma -= fZommerfeld*std::log(fHalfRutThetaTg2/(1.+fHalfRutThetaTg2));
    9611066  twosigma += fRutherfordTheta*fZommerfeld/fHalfRutThetaTg - halfpi;
    9621067  twosigma += fProfileLambda*theta - 0.25*pi;
    9631068
     1069  twosigma *= fCofPhase;
     1070
    9641071  G4complex z = G4complex(0., twosigma);
    9651072
    … …  
    9741081inline G4complex G4NuclNuclDiffuseElastic::GammaLess(G4double theta)
    9751082{
    976   G4double sinThetaR      = 2.*fHalfRutThetaTg/(1. + fHalfRutThetaTg*fHalfRutThetaTg);
    977   G4double cosHalfThetaR2 = 1./(1. + fHalfRutThetaTg*fHalfRutThetaTg);
     1083  G4double sinThetaR      = 2.*fHalfRutThetaTg/(1. + fHalfRutThetaTg2);
     1084  G4double cosHalfThetaR2 = 1./(1. + fHalfRutThetaTg2);
    9781085
    9791086  G4double u              = std::sqrt(0.5*fProfileLambda/sinThetaR);
    … …  
    9871094  G4complex order         = G4complex(u,u);
    9881095  order                  /= std::sqrt(2.);
     1096
    9891097  G4complex gamma         = pi*kappa*GetErfcInt(-order)*std::exp(im*(u*u+0.25*pi));
    9901098  G4complex a0            = 0.5*(1. + 4.*(1.+im*u2)*cosHalfThetaR2/3.)/sinThetaR;
    9911099  G4complex a1            = 0.5*(1. + 2.*(1.+im*u2m2p3)*cosHalfThetaR2)/sinThetaR;
    9921100  G4complex out           = gamma*(1. - a1*dTheta) - a0;
     1101
    9931102  return out;
    9941103}
    … …  
    10001109inline G4complex G4NuclNuclDiffuseElastic::GammaMore(G4double theta)
    10011110{
    1002   G4double sinThetaR      = 2.*fHalfRutThetaTg/(1. + fHalfRutThetaTg*fHalfRutThetaTg);
    1003   G4double cosHalfThetaR2 = 1./(1. + fHalfRutThetaTg*fHalfRutThetaTg);
     1111  G4double sinThetaR      = 2.*fHalfRutThetaTg/(1. + fHalfRutThetaTg2);
     1112  G4double cosHalfThetaR2 = 1./(1. + fHalfRutThetaTg2);
    10041113
    10051114  G4double u              = std::sqrt(0.5*fProfileLambda/sinThetaR);
    … …  
    10141123  order                  /= std::sqrt(2.);
    10151124  G4complex gamma         = pi*kappa*GetErfcInt(order)*std::exp(im*(u*u+0.25*pi));
    1016   G4complex a0            = 0.5*(1. + 3.*(1.+im*u2)*cosHalfThetaR2/3.)/sinThetaR;
     1125  G4complex a0            = 0.5*(1. + 4.*(1.+im*u2)*cosHalfThetaR2/3.)/sinThetaR;
    10171126  G4complex a1            = 0.5*(1. + 2.*(1.+im*u2m2p3)*cosHalfThetaR2)/sinThetaR;
    10181127  G4complex out           = -gamma*(1. - a1*dTheta) - a0;
     1128
    10191129  return out;
    10201130}
    … …  
    10281138  G4double kappa = std::sqrt(0.5*fProfileLambda/std::sin(theta)/pi);
    10291139  G4complex out = G4complex(kappa/fWaveVector,0.);
     1140
    10301141  out *= PhaseNear(theta);
    10311142
    1032   if(theta <= fRutherfordTheta)
     1143  if( theta <= fRutherfordTheta )
    10331144  {
    10341145    out *= GammaLess(theta) + ProfileNear(theta);
     1146    // out *= GammaMore(theta) + ProfileNear(theta);
    10351147    out += CoulombAmplitude(theta);
    10361148  }
    … …  
    10381150  {
    10391151    out *= GammaMore(theta) + ProfileNear(theta);
     1152    // out *= GammaLess(theta) + ProfileNear(theta);
    10401153  }
    10411154  return out;
    … …  
    10501163  G4double kappa = std::sqrt(0.5*fProfileLambda/std::sin(theta)/pi);
    10511164  G4complex out = G4complex(kappa/fWaveVector,0.);
    1052   out *= ProfileFar(theta)*PhaseFar(theta);
     1165  out *= ProfileFar(theta);
     1166  out *= PhaseFar(theta);
    10531167  return out;
    10541168}
    … …  
    10621176{
    10631177 
    1064   G4complex out = AmplitudeNear(theta) + AmplitudeFar(theta);
     1178  G4complex out = AmplitudeNear(theta) + fCofFar*AmplitudeFar(theta);
     1179  // G4complex out = AmplitudeNear(theta);
     1180  // G4complex out = AmplitudeFar(theta);
    10651181  return    out;
    10661182}
    … …  
    10771193}
    10781194
     1195/////////////////////////////////////////////////////////////////
     1196//
     1197//
     1198
     1199inline  G4complex G4NuclNuclDiffuseElastic::AmplitudeGla(G4double theta)
     1200{
     1201  G4int n;
     1202  G4double T12b, b, b2; // cosTheta = std::cos(theta);
     1203  G4complex out = G4complex(0.,0.), shiftC, shiftN;
     1204  G4complex im  = G4complex(0.,1.);
     1205
     1206  for( n = 0; n < fMaxL; n++)
     1207  {
     1208    shiftC = std::exp( im*2.*CalculateCoulombPhase(n) );
     1209    // b = ( fZommerfeld + std::sqrt( fZommerfeld*fZommerfeld + n*(n+1) ) )/fWaveVector;
     1210    b = ( std::sqrt( G4double(n*(n+1)) ) )/fWaveVector;
     1211    b2 = b*b;
     1212    T12b = fSumSigma*std::exp(-b2/fNuclearRadiusSquare)/pi/fNuclearRadiusSquare;         
     1213    shiftN = std::exp( -0.5*(1.-im*fEtaRatio)*T12b ) - 1.;
     1214    out +=  (2.*n+1.)*shiftC*shiftN*GetLegendrePol(n, theta);   
     1215  }
     1216  out /= 2.*im*fWaveVector;
     1217  out += CoulombAmplitude(theta);
     1218  return out;
     1219}
     1220
     1221/////////////////////////////////////////////////////////////////
     1222//
     1223//
     1224
     1225inline  G4double  G4NuclNuclDiffuseElastic::AmplitudeGlaMod2(G4double theta)
     1226{
     1227  G4complex out = AmplitudeGla(theta);
     1228  G4double mod2 = out.real()*out.real() + out.imag()*out.imag();
     1229  return   mod2;
     1230}
     1231
     1232
     1233/////////////////////////////////////////////////////////////////
     1234//
     1235//
     1236
     1237inline  G4complex G4NuclNuclDiffuseElastic::AmplitudeGG(G4double theta)
     1238{
     1239  G4int n;
     1240  G4double T12b, a, aTemp, b2, sinThetaH = std::sin(0.5*theta);
     1241  G4double sinThetaH2 = sinThetaH*sinThetaH;
     1242  G4complex out = G4complex(0.,0.);
     1243  G4complex im  = G4complex(0.,1.);
     1244
     1245  a  = -fSumSigma/twopi/fNuclearRadiusSquare;
     1246  b2 = fWaveVector*fWaveVector*fNuclearRadiusSquare*sinThetaH2;
     1247
     1248  aTemp = a;
     1249
     1250  for( n = 1; n < fMaxL; n++)
     1251  {   
     1252    T12b   = aTemp*std::exp(-b2/n)/n;         
     1253    aTemp *= a; 
     1254    out   += T12b;
     1255    G4cout<<"out = "<<out<<G4endl; 
     1256  }
     1257  out *= -4.*im*fWaveVector/pi;
     1258  out += CoulombAmplitude(theta);
     1259  return out;
     1260}
     1261
     1262/////////////////////////////////////////////////////////////////
     1263//
     1264//
     1265
     1266inline  G4double  G4NuclNuclDiffuseElastic::AmplitudeGGMod2(G4double theta)
     1267{
     1268  G4complex out = AmplitudeGG(theta);
     1269  G4double mod2 = out.real()*out.real() + out.imag()*out.imag();
     1270  return   mod2;
     1271}
     1272
    10791273
    10801274///////////////////////////////////////////////////////////////////////////////
    10811275//
    10821276// Test for given particle and element table of momentum, angle probability.
    1083 // For the moment in lab system.
     1277// For the partMom in CMS.
    10841278
    10851279inline void G4NuclNuclDiffuseElastic::InitParameters(const G4ParticleDefinition* theParticle, 
    … …  
    11011295  fWaveVector = partMom/hbarc;
    11021296
    1103   G4double lambda = fWaveVector*fNuclearRadius;
     1297  G4double lambda = fCofLambda*fWaveVector*fNuclearRadius;
     1298  G4cout<<"kR = "<<lambda<<G4endl;
     1299
     1300  if( z )
     1301  {
     1302    a           = partMom/m1; // beta*gamma for m1
     1303    fBeta       = a/std::sqrt(1+a*a);
     1304    fZommerfeld = CalculateZommerfeld( fBeta, z, fAtomicNumber);
     1305    fAm         = CalculateAm( partMom, fZommerfeld, fAtomicNumber);
     1306  }
     1307  fProfileLambda = lambda*std::sqrt(1.-2*fZommerfeld/lambda);
     1308  G4cout<<"fProfileLambda = "<<fProfileLambda<<G4endl;
     1309  fProfileDelta  = fCofDelta*fProfileLambda;
     1310  fProfileAlpha  = fCofAlpha*fProfileLambda;
     1311
     1312  CalculateCoulombPhaseZero();
     1313  CalculateRutherfordAnglePar();
     1314
     1315  return;
     1316}
     1317
     1318
     1319///////////////////////////////////////////////////////////////////////////////
     1320//
     1321// Test for given particle and element table of momentum, angle probability.
     1322// For the partMom in CMS.
     1323
     1324inline void G4NuclNuclDiffuseElastic::InitParametersGla(const G4DynamicParticle* aParticle, 
     1325                                          G4double partMom, G4double Z, G4double A)
     1326{
     1327  fAtomicNumber  = Z;     // target atomic number
     1328  fAtomicWeight  = A;     // target number of nucleons
     1329
     1330  fNuclearRadius2 = CalculateNuclearRad(fAtomicWeight); // target nucleus radius
     1331  G4double A1     = G4double( aParticle->GetDefinition()->GetBaryonNumber() );   
     1332  fNuclearRadius1 = CalculateNuclearRad(A1); // projectile nucleus radius
     1333  fNuclearRadiusSquare = fNuclearRadius1*fNuclearRadius1+fNuclearRadius2*fNuclearRadius2;
     1334 
     1335
     1336  G4double a  = 0., kR12;
     1337  G4double z  = aParticle->GetDefinition()->GetPDGCharge();
     1338  G4double m1 = aParticle->GetDefinition()->GetPDGMass();
     1339
     1340  fWaveVector = partMom/hbarc;
     1341
     1342  G4double pN = A1 - z;
     1343  if( pN < 0. ) pN = 0.;
     1344
     1345  G4double tN = A - Z;
     1346  if( tN < 0. ) tN = 0.;
     1347
     1348  G4double pTkin = aParticle->GetKineticEnergy(); 
     1349  pTkin /= A1;
     1350
     1351
     1352  fSumSigma = (Z*z+pN*tN)*GetHadronNucleonXscNS(theProton, pTkin, theProton) +
     1353              (z*tN+pN*Z)*GetHadronNucleonXscNS(theProton, pTkin, theNeutron);
     1354
     1355  G4cout<<"fSumSigma = "<<fSumSigma/millibarn<<" mb"<<G4endl;
     1356  G4cout<<"pi*R2 = "<<pi*fNuclearRadiusSquare/millibarn<<" mb"<<G4endl;
     1357  kR12 = fWaveVector*std::sqrt(fNuclearRadiusSquare);
     1358  G4cout<<"k*sqrt(R2) = "<<kR12<<" "<<G4endl;
     1359  fMaxL = (G4int(kR12)+1)*4;
     1360  G4cout<<"fMaxL = "<<fMaxL<<" "<<G4endl;
    11041361
    11051362  if( z )
    … …  
    11101367      fAm         = CalculateAm( partMom, fZommerfeld, fAtomicNumber);
    11111368  }
    1112   fProfileLambda = lambda*std::sqrt(1.-2*fZommerfeld/lambda);
    1113   G4cout<<"fProfileLambda = "<<fProfileLambda<<G4endl;
    1114   fProfileDelta  = 0.1*fProfileLambda;
    1115   fProfileAlpha  = 0.05*fProfileLambda;
    11161369
    11171370  CalculateCoulombPhaseZero();
    1118   CalculateRutherfordAnglePar();
     1371 
    11191372
    11201373  return;
    … …  
    11221375
    11231376
     1377/////////////////////////////////////////////////////////////////////////////////////
     1378//
     1379// Returns nucleon-nucleon cross-section based on N. Starkov parametrisation of
     1380// data from mainly http://wwwppds.ihep.su:8001/c5-6A.html database
     1381// projectile nucleon is pParticle with pTkin shooting target nucleon tParticle
     1382
     1383inline G4double
     1384G4NuclNuclDiffuseElastic::GetHadronNucleonXscNS( G4ParticleDefinition* pParticle,
     1385                                                 G4double pTkin,
     1386                                                 G4ParticleDefinition* tParticle)
     1387{
     1388  G4double xsection(0), Delta, A0, B0;
     1389  G4double hpXsc(0);
     1390  G4double hnXsc(0);
     1391
     1392
     1393  G4double targ_mass     = tParticle->GetPDGMass();
     1394  G4double proj_mass     = pParticle->GetPDGMass();
     1395
     1396  G4double proj_energy   = proj_mass + pTkin;
     1397  G4double proj_momentum = std::sqrt(pTkin*(pTkin+2*proj_mass));
     1398
     1399  G4double sMand = CalcMandelstamS ( proj_mass , targ_mass , proj_momentum );
     1400
     1401  sMand         /= GeV*GeV;  // in GeV for parametrisation
     1402  proj_momentum /= GeV;
     1403  proj_energy   /= GeV;
     1404  proj_mass     /= GeV;
     1405  G4double logS = std::log(sMand);
     1406
     1407  // General PDG fit constants
     1408
     1409
     1410  // fEtaRatio=Re[f(0)]/Im[f(0)]
     1411
     1412  if( proj_momentum >= 1.2 )
     1413  {
     1414    fEtaRatio  = 0.13*(logS - 5.8579332)*std::pow(sMand,-0.18);
     1415  }       
     1416  else if( proj_momentum >= 0.6 )
     1417  {
     1418    fEtaRatio = -75.5*(std::pow(proj_momentum,0.25)-0.95)/
     1419          (std::pow(3*proj_momentum,2.2)+1);     
     1420  }
     1421  else
     1422  {
     1423    fEtaRatio = 15.5*proj_momentum/(27*proj_momentum*proj_momentum*proj_momentum+2);
     1424  }
     1425  G4cout<<"fEtaRatio = "<<fEtaRatio<<G4endl;
     1426
     1427  // xsc
     1428 
     1429  if( proj_momentum >= 10. ) // high energy: pp = nn = np
     1430    // if( proj_momentum >= 2.)
     1431  {
     1432    Delta = 1.;
     1433
     1434    if( proj_energy < 40. ) Delta = 0.916+0.0021*proj_energy;
     1435
     1436    if( proj_momentum >= 10.)
     1437    {
     1438        B0 = 7.5;
     1439        A0 = 100. - B0*std::log(3.0e7);
     1440
     1441        xsection = A0 + B0*std::log(proj_energy) - 11
     1442                  + 103*std::pow(2*0.93827*proj_energy + proj_mass*proj_mass+
     1443                     0.93827*0.93827,-0.165);        //  mb
     1444    }
     1445  }
     1446  else // low energy pp = nn != np
     1447  {
     1448      if(pParticle == tParticle) // pp or nn      // nn to be pp
     1449      {
     1450        if( proj_momentum < 0.73 )
     1451        {
     1452          hnXsc = 23 + 50*( std::pow( std::log(0.73/proj_momentum), 3.5 ) );
     1453        }
     1454        else if( proj_momentum < 1.05  )
     1455        {
     1456          hnXsc = 23 + 40*(std::log(proj_momentum/0.73))*
     1457                         (std::log(proj_momentum/0.73));
     1458        }
     1459        else  // if( proj_momentum < 10.  )
     1460        {
     1461          hnXsc = 39.0 +
     1462              75*(proj_momentum - 1.2)/(std::pow(proj_momentum,3.0) + 0.15);
     1463        }
     1464        xsection = hnXsc;
     1465      }
     1466      else  // pn to be np
     1467      {
     1468        if( proj_momentum < 0.8 )
     1469        {
     1470          hpXsc = 33+30*std::pow(std::log(proj_momentum/1.3),4.0);
     1471        }     
     1472        else if( proj_momentum < 1.4 )
     1473        {
     1474          hpXsc = 33+30*std::pow(std::log(proj_momentum/0.95),2.0);
     1475        }
     1476        else    // if( proj_momentum < 10.  )
     1477        {
     1478          hpXsc = 33.3+
     1479              20.8*(std::pow(proj_momentum,2.0)-1.35)/
     1480                 (std::pow(proj_momentum,2.50)+0.95);
     1481        }
     1482        xsection = hpXsc;
     1483      }
     1484  }
     1485  xsection *= millibarn; // parametrised in mb
     1486  G4cout<<"xsection = "<<xsection/millibarn<<" mb"<<G4endl;
     1487  return xsection;
     1488}
     1489
     1490////////////////////////////////////////////////////////////////////////////////////
     1491//
     1492//
     1493
     1494inline G4double G4NuclNuclDiffuseElastic::CalcMandelstamS( const G4double mp ,
     1495                                                       const G4double mt ,
     1496                                                       const G4double Plab )
     1497{
     1498  G4double Elab = std::sqrt ( mp * mp + Plab * Plab );
     1499  G4double sMand  = mp*mp + mt*mt + 2*Elab*mt ;
     1500
     1501  return sMand;
     1502}
    11241503
    11251504
Note: See TracChangeset for help on using the changeset viewer.