Changeset 1337 for trunk/source/processes/hadronic/cross_sections/src/G4GlauberGribovCrossSection.cc
- Timestamp:
- Sep 30, 2010, 2:47:17 PM (14 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/source/processes/hadronic/cross_sections/src/G4GlauberGribovCrossSection.cc
r1228 r1337 24 24 // ******************************************************************** 25 25 // 26 // 26 // author: V. Grichine 27 // 27 28 // 17.07.06 V. Grichine - first implementation 28 29 // 22.01.07 V.Ivanchenko - add interface with Z and A 29 30 // 05.03.07 V.Ivanchenko - add IfZAApplicable 30 // 31 // 11.06.10 V. Grichine - update for antiprotons 31 32 32 33 #include "G4GlauberGribovCrossSection.hh" … … 361 362 // cofInelastic = 2.0; 362 363 363 364 nucleusSquare = cofTotal*pi*R*R; // basically 2piRR 365 ratio = sigma/nucleusSquare; 366 367 xsection = nucleusSquare*std::log( 1. + ratio ); 368 369 xsection *= GetParticleBarCorTot(theParticle, Z); 370 371 fTotalXsc = xsection; 364 if( A > 1.5 ) 365 { 366 nucleusSquare = cofTotal*pi*R*R; // basically 2piRR 367 ratio = sigma/nucleusSquare; 368 369 xsection = nucleusSquare*std::log( 1. + ratio ); 370 371 xsection *= GetParticleBarCorTot(theParticle, Z); 372 373 fTotalXsc = xsection; 372 374 373 375 374 376 375 fInelasticXsc = nucleusSquare*std::log( 1. + cofInelastic*ratio )/cofInelastic;376 377 fInelasticXsc *= GetParticleBarCorIn(theParticle, Z);378 379 fElasticXsc = fTotalXsc - fInelasticXsc;377 fInelasticXsc = nucleusSquare*std::log( 1. + cofInelastic*ratio )/cofInelastic; 378 379 fInelasticXsc *= GetParticleBarCorIn(theParticle, Z); 380 381 fElasticXsc = fTotalXsc - fInelasticXsc; 380 382 381 383 382 G4double difratio = ratio/(1.+ratio); 383 384 fDiffractionXsc = 0.5*nucleusSquare*( difratio - std::log( 1. + difratio ) ); 385 386 387 sigma = GetHNinelasticXsc(aParticle, A, Z); 388 ratio = sigma/nucleusSquare; 389 390 fProductionXsc = nucleusSquare*std::log( 1. + cofInelastic*ratio )/cofInelastic; 391 392 if (fElasticXsc < 0.) fElasticXsc = 0.; 393 384 G4double difratio = ratio/(1.+ratio); 385 386 fDiffractionXsc = 0.5*nucleusSquare*( difratio - std::log( 1. + difratio ) ); 387 388 389 sigma = GetHNinelasticXsc(aParticle, A, Z); 390 ratio = sigma/nucleusSquare; 391 392 fProductionXsc = nucleusSquare*std::log( 1. + cofInelastic*ratio )/cofInelastic; 393 394 if (fElasticXsc < 0.) fElasticXsc = 0.; 395 } 396 else // H 397 { 398 fTotalXsc = sigma; 399 xsection = sigma; 400 401 if ( theParticle != theAProton ) 402 { 403 sigma = GetHNinelasticXsc(aParticle, A, Z); 404 fInelasticXsc = sigma; 405 fElasticXsc = fTotalXsc - fInelasticXsc; 406 } 407 else 408 { 409 fElasticXsc = fTotalXsc - fInelasticXsc; 410 } 411 if (fElasticXsc < 0.) fElasticXsc = 0.; 412 413 } 394 414 return xsection; 395 415 } … … 607 627 608 628 G4double Nt = At-Zt; // number of neutrons 629 609 630 if (Nt < 0.) Nt = 0.; 610 631 … … 740 761 741 762 G4double Nt = At-Zt; // number of neutrons 763 742 764 if (Nt < 0.) Nt = 0.; 743 765 … … 887 909 // xsection *= 0.95; 888 910 } 889 else if(theParticle == theAProton) 890 { 891 xsection = Zt*( 35.45 + B*std::pow(std::log(sMand/s0),2.) 892 + 42.53*std::pow(sMand,-eta1) + 33.34*std::pow(sMand,-eta2)); 893 894 xsection += Nt*( 35.80 + B*std::pow(std::log(sMand/s0),2.) 895 + 40.15*std::pow(sMand,-eta1) + 30.*std::pow(sMand,-eta2)); 896 } 897 else if(theParticle == thePiPlus) 911 else if( theParticle == theAProton ) 912 { 913 // xsection = Zt*( 35.45 + B*std::pow(std::log(sMand/s0),2.) 914 // + 42.53*std::pow(sMand,-eta1) + 33.34*std::pow(sMand,-eta2)); 915 916 // xsection += Nt*( 35.80 + B*std::pow(std::log(sMand/s0),2.) 917 // + 40.15*std::pow(sMand,-eta1) + 30.*std::pow(sMand,-eta2)); 918 919 G4double logP = std::log(proj_momentum); 920 921 if( proj_momentum <= 1.0 ) 922 { 923 xsection = Zt*(65.55 + 53.84/(proj_momentum+1.e-6) ); 924 } 925 else 926 { 927 xsection = Zt*( 41.1 + 77.2*std::pow( proj_momentum, -0.68) 928 + 0.293*logP*logP - 1.82*logP ); 929 } 930 if ( Nt > 0.) 931 { 932 xsection += Nt*( 41.9 + 96.2*std::pow( proj_momentum, -0.99) - 0.154*logP); 933 } 934 else // H 935 { 936 fInelasticXsc = 38.0 + 38.0*std::pow( proj_momentum, -0.96) 937 - 0.169*logP*logP; 938 fInelasticXsc *= millibarn; 939 } 940 } 941 else if( theParticle == thePiPlus ) 898 942 { 899 943 if(proj_momentum < 0.4) … … 902 946 hpXsc = Ex3+20.0; 903 947 } 904 else if( proj_momentum < 1.15)948 else if( proj_momentum < 1.15 ) 905 949 { 906 950 G4double Ex4 = 88*(std::log(proj_momentum/0.75))*(std::log(proj_momentum/0.75)); … … 1107 1151 G4double NumberOfTargetNeutrons = NumberOfTargetNucleons - NumberOfTargetProtons; 1108 1152 1109 if( NumberOfTargetNeutrons < 0.) NumberOfTargetNeutrons = 0.;1153 if( NumberOfTargetNeutrons < 0.) NumberOfTargetNeutrons = 0.; 1110 1154 1111 1155 G4double Xtotal, Xelastic, Xinelastic; … … 1286 1330 Xinelastic = Xtotal - Xelastic; 1287 1331 1288 if( Xinelastic < 0.) Xinelastic = 0.;1332 if( Xinelastic < 0.) Xinelastic = 0.; 1289 1333 1290 1334 return Xinelastic*= millibarn;
Note: See TracChangeset
for help on using the changeset viewer.