- Timestamp:
- May 28, 2009, 4:26:57 PM (15 years ago)
- Location:
- trunk/source/processes/hadronic/models/coherent_elastic/src
- Files:
-
- 2 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/source/processes/hadronic/models/coherent_elastic/src/G4ChargeExchange.cc
r1007 r1055 25 25 // 26 26 // 27 // $Id: G4ChargeExchange.cc,v 1.1 4 2008/12/18 13:01:48 gunterExp $28 // GEANT4 tag $Name: geant4-09-0 2$27 // $Id: G4ChargeExchange.cc,v 1.13 2008/11/20 12:35:19 vnivanch Exp $ 28 // GEANT4 tag $Name: geant4-09-03-beta-cand-01 $ 29 29 // 30 30 // … … 312 312 dd = 10.; 313 313 } 314 G4double x1 = (1.0 - std::exp(-tmax*bb))*aa/bb;315 G4double x2 = (1.0 - std::exp(-tmax*dd))*cc/dd;314 G4double x1 = (1.0 - exp(-tmax*bb))*aa/bb; 315 G4double x2 = (1.0 - exp(-tmax*dd))*cc/dd; 316 316 317 317 G4double t; -
trunk/source/processes/hadronic/models/coherent_elastic/src/G4DiffuseElastic.cc
r1007 r1055 24 24 // ******************************************************************** 25 25 // 26 // $Id: G4DiffuseElastic.cc,v 1.2 0 2008/01/14 10:39:13 vnivanchExp $27 // GEANT4 tag $Name: geant4-09-0 2$26 // $Id: G4DiffuseElastic.cc,v 1.23 2009/03/02 09:17:43 grichine Exp $ 27 // GEANT4 tag $Name: geant4-09-03-beta-cand-01 $ 28 28 // 29 29 // … … 67 67 { 68 68 SetMinEnergy( 0.01*GeV ); 69 SetMaxEnergy( 1 00.*TeV );69 SetMaxEnergy( 1.*TeV ); 70 70 verboseLevel = 0; 71 71 lowEnergyRecoilLimit = 100.*keV; … … 83 83 84 84 fEnergyBin = 200; 85 fAngleBin = 100;86 87 fEnergyVector = 0;85 fAngleBin = 200; 86 87 fEnergyVector = new G4PhysicsLogVector( theMinEnergy, theMaxEnergy, fEnergyBin ); 88 88 fAngleTable = 0; 89 89 … … 107 107 { 108 108 SetMinEnergy( 0.01*GeV ); 109 SetMaxEnergy( 1 00.*TeV );109 SetMaxEnergy( 1.*TeV ); 110 110 verboseLevel = 0; 111 111 lowEnergyRecoilLimit = 100.*keV; … … 122 122 thePionMinus= G4PionMinus::PionMinus(); 123 123 124 fEnergyBin = 200; 125 fAngleBin =100;124 fEnergyBin = 200; // 200; // 100; 125 fAngleBin = 400; // 200; // 100; 126 126 127 127 // fEnergyVector = 0; … … 174 174 fAtomicWeight = (*theElementTable)[jEl]->GetN(); // number of nucleons 175 175 fNuclearRadius = CalculateNuclearRad(fAtomicWeight); 176 if(verboseLevel > 0) 176 177 if(verboseLevel > 0) 178 { 177 179 G4cout<<"G4DiffuseElastic::Initialise() the element: " 178 180 <<(*theElementTable)[jEl]->GetName()<<G4endl; 181 } 179 182 fElementNumberVector.push_back(fAtomicNumber); 180 183 fElementNameVector.push_back((*theElementTable)[jEl]->GetName()); … … 183 186 fAngleBank.push_back(fAngleTable); 184 187 } 185 return;186 }187 188 //////////////////////////////////////////////////////////////////////////////189 //190 // Initialisation for given particle on fly using new element number191 192 void G4DiffuseElastic::InitialiseOnFly(G4double Z, G4double A)193 {194 fAtomicNumber = Z; // atomic number195 fAtomicWeight = A; // number of nucleons196 fNuclearRadius = CalculateNuclearRad(fAtomicWeight);197 if(verboseLevel > 0)198 G4cout<<"G4DiffuseElastic::Initialise() the element with Z = "199 <<Z<<"; and A = "<<A<<G4endl;200 fElementNumberVector.push_back(fAtomicNumber);201 202 BuildAngleTable();203 fAngleBank.push_back(fAngleTable);204 return;205 }206 207 ///////////////////////////////////////////////////////////////////////////////208 //209 // Build for given particle and element table of momentum, angle probability.210 // For the moment in lab system.211 212 void G4DiffuseElastic::BuildAngleTable()213 {214 G4int i, j;215 G4double partMom, kinE, a=0., z = fParticle->GetPDGCharge(), m1 = fParticle->GetPDGMass();216 G4double theta1, theta2, thetaMax, thetaCoulomb, sum = 0.;217 218 G4Integrator<G4DiffuseElastic,G4double(G4DiffuseElastic::*)(G4double)> integral;219 220 fAngleTable = new G4PhysicsTable(fEnergyBin);221 222 for(i = 0; i < fEnergyBin; i++)223 {224 kinE = fEnergyVector->GetLowEdgeEnergy(i);225 partMom = std::sqrt( kinE*(kinE + 2*m1) );226 fWaveVector = partMom/hbarc;227 228 thetaMax = 10.174/fWaveVector/fNuclearRadius;229 230 if (thetaMax > pi) thetaMax = pi;231 232 thetaCoulomb = 0.2*thetaMax;233 234 if(z)235 {236 a = partMom/m1;237 fBeta = a/std::sqrt(1+a*a);238 fZommerfeld = CalculateZommerfeld( fBeta, z, fAtomicNumber);239 fAm = CalculateAm( partMom, fZommerfeld, fAtomicNumber);240 }241 G4PhysicsFreeVector* angleVector = new G4PhysicsFreeVector(fAngleBin);242 243 G4PhysicsLogVector* angleBins = new G4PhysicsLogVector( 0.01*thetaMax, thetaMax, fAngleBin );244 245 for(j = 1; j < fAngleBin; j++)246 {247 theta1 = angleBins->GetLowEdgeEnergy(j-1);248 theta2 = angleBins->GetLowEdgeEnergy(j);249 250 if(theta2 > thetaCoulomb && z) fAddCoulomb = true;251 252 sum += integral.Legendre10(this,&G4DiffuseElastic::GetIntegrandFunction, theta1,theta2);253 254 angleVector->PutValue( j-1 , theta2, sum );255 // G4cout<<"j-1 = "<<j-1<<"; theta2 = "<<theta2<<"; sum = "<<sum<<G4endl;256 }257 fAddCoulomb = false;258 259 fAngleTable->insertAt(i,angleVector);260 261 // delete[] angleVector;262 // delete[] angleBins;263 }264 188 return; 265 189 } … … 525 449 fAtomicWeight = A; 526 450 fAtomicNumber = Z; 527 G4double z = particle->GetPDGCharge(); 528 if(z) 451 fNuclearRadius = CalculateNuclearRad(A); 452 fAddCoulomb = false; 453 454 G4double z = particle->GetPDGCharge(); 455 456 G4double kRt = fWaveVector*fNuclearRadius*theta; 457 G4double kRtC = 1.9; 458 459 if( z && (kRt > kRtC) ) 529 460 { 530 461 fAddCoulomb = true; 531 fBeta = CalculateParticleBeta( particle, momentum); 532 fZommerfeld = CalculateZommerfeld( fBeta, z, fAtomicNumber); 533 fAm = CalculateAm( momentum, fZommerfeld, fAtomicNumber); 534 } 535 fNuclearRadius = CalculateNuclearRad(A); 536 462 fBeta = CalculateParticleBeta( particle, momentum); 463 fZommerfeld = CalculateZommerfeld( fBeta, z, fAtomicNumber); 464 fAm = CalculateAm( momentum, fZommerfeld, fAtomicNumber); 465 } 537 466 G4double sigma = fNuclearRadius*fNuclearRadius*GetDiffElasticSumProb(theta); 538 467 … … 552 481 { 553 482 G4double m1 = particle->GetPDGMass(); 483 554 484 G4LorentzVector lv1(0.,0.,plab,std::sqrt(plab*plab+m1*m1)); 555 485 556 486 G4int iZ = static_cast<G4int>(Z+0.5); 557 487 G4int iA = static_cast<G4int>(A+0.5); 558 G4ParticleDefinition * theDef = 0; 488 489 G4ParticleDefinition* theDef = 0; 559 490 560 491 if (iZ == 1 && iA == 1) theDef = theProton; … … 575 506 G4ThreeVector p1 = lv1.vect(); 576 507 G4double ptot = p1.mag(); 577 G4double ptot2 = ptot*ptot;578 G4double cost = 1 - 0.5*std::fabs(tMand)/ptot2;508 G4double ptot2 = ptot*ptot; 509 G4double cost = 1 - 0.5*std::fabs(tMand)/ptot2; 579 510 580 511 if( cost >= 1.0 ) cost = 1.0; … … 658 589 // G4double r0 = 1.08*fermi; 659 590 // G4double rad = r0*std::pow(A, 1./3.); 591 660 592 G4double kr = fWaveVector*fNuclearRadius; // wavek*rad; 661 593 G4double kr2 = kr*kr; … … 686 618 } 687 619 G4double lambda = 15.; // 15 ok 620 688 621 // G4double kg = fWaveVector*gamma; // wavek*delta; 622 689 623 G4double kg = lambda*(1.-std::exp(-fWaveVector*gamma/lambda)); // wavek*delta; 690 624 G4double kg2 = kg*kg; 625 691 626 // G4double dk2t = delta*fWaveVector*fWaveVector*theta; // delta*wavek*wavek*theta; 692 627 // G4double dk2t2 = dk2t*dk2t; 693 628 // G4double pikdt = pi*fWaveVector*diffuse*theta;// pi*wavek*diffuse*theta; 629 694 630 G4double pikdt = lambda*(1.-std::exp(-pi*fWaveVector*diffuse*theta/lambda)); // wavek*delta; 695 631 … … 730 666 // G4double r0 = 1.08*fermi; 731 667 // G4double rad = r0*std::pow(A, 1./3.); 668 732 669 G4double kr = fWaveVector*fNuclearRadius; // wavek*rad; 733 670 G4double kr2 = kr*kr; … … 774 711 775 712 G4double kg2 = kg*kg; 713 776 714 // G4double dk2t = delta*fWaveVector*fWaveVector*theta; // delta*wavek*wavek*theta; 777 778 715 // G4cout<<"dk2t = "<<dk2t<<G4endl; 779 780 716 // G4double dk2t2 = dk2t*dk2t; 781 782 717 // G4double pikdt = pi*fWaveVector*diffuse*theta;// pi*wavek*diffuse*theta; 718 783 719 G4double pikdt = lambda*(1.-std::exp(-pi*fWaveVector*diffuse*theta/lambda)); // wavek*delta; 784 720 … … 808 744 //////////////////////////////////////////////////////////////////////////// 809 745 // 746 // return differential elastic probability d(probability)/d(t) with 747 // Coulomb correction 748 749 G4double 750 G4DiffuseElastic::GetDiffElasticSumProbA( G4double alpha ) 751 { 752 G4double theta; 753 754 theta = std::sqrt(alpha); 755 756 // theta = std::acos( 1 - alpha/2. ); 757 758 G4double sigma, bzero, bzero2, bonebyarg, bonebyarg2, damp, damp2; 759 G4double delta, diffuse, gamma; 760 G4double e1, e2, bone, bone2; 761 762 // G4double wavek = momentum/hbarc; // wave vector 763 // G4double r0 = 1.08*fermi; 764 // G4double rad = r0*std::pow(A, 1./3.); 765 766 G4double kr = fWaveVector*fNuclearRadius; // wavek*rad; 767 G4double kr2 = kr*kr; 768 G4double krt = kr*theta; 769 770 bzero = BesselJzero(krt); 771 bzero2 = bzero*bzero; 772 bone = BesselJone(krt); 773 bone2 = bone*bone; 774 bonebyarg = BesselOneByArg(krt); 775 bonebyarg2 = bonebyarg*bonebyarg; 776 777 if (fParticle == theProton) 778 { 779 diffuse = 0.63*fermi; 780 // diffuse = 0.6*fermi; 781 gamma = 0.3*fermi; 782 delta = 0.1*fermi*fermi; 783 e1 = 0.3*fermi; 784 e2 = 0.35*fermi; 785 } 786 else // as proton, if were not defined 787 { 788 diffuse = 0.63*fermi; 789 gamma = 0.3*fermi; 790 delta = 0.1*fermi*fermi; 791 e1 = 0.3*fermi; 792 e2 = 0.35*fermi; 793 } 794 G4double lambda = 15.; // 15 ok 795 // G4double kg = fWaveVector*gamma; // wavek*delta; 796 G4double kg = lambda*(1.-std::exp(-fWaveVector*gamma/lambda)); // wavek*delta; 797 798 // G4cout<<"kg = "<<kg<<G4endl; 799 800 if(fAddCoulomb) // add Coulomb correction 801 { 802 G4double sinHalfTheta = theta*0.5; // std::sin(0.5*theta); 803 G4double sinHalfTheta2 = sinHalfTheta*sinHalfTheta; 804 805 kg += 0.5*fZommerfeld/kr/(sinHalfTheta2+fAm); // correction at J0() 806 // kg += 0.65*fZommerfeld/kr/(sinHalfTheta2+fAm); // correction at J0() 807 } 808 809 G4double kg2 = kg*kg; 810 811 // G4double dk2t = delta*fWaveVector*fWaveVector*theta; // delta*wavek*wavek*theta; 812 // G4cout<<"dk2t = "<<dk2t<<G4endl; 813 // G4double dk2t2 = dk2t*dk2t; 814 // G4double pikdt = pi*fWaveVector*diffuse*theta;// pi*wavek*diffuse*theta; 815 816 G4double pikdt = lambda*(1.-std::exp(-pi*fWaveVector*diffuse*theta/lambda)); // wavek*delta; 817 818 // G4cout<<"pikdt = "<<pikdt<<G4endl; 819 820 damp = DampFactor(pikdt); 821 damp2 = damp*damp; 822 823 G4double mode2k2 = (e1*e1+e2*e2)*fWaveVector*fWaveVector; 824 G4double e2dk3t = -2.*e2*delta*fWaveVector*fWaveVector*fWaveVector*theta; 825 826 sigma = kg2; 827 // sigma += dk2t2; 828 sigma *= bzero2; 829 sigma += mode2k2*bone2; 830 sigma += e2dk3t*bzero*bone; 831 832 // sigma += kr2*(1 + 8.*fZommerfeld*fZommerfeld/kr2)*bonebyarg2; // correction at J1()/() 833 sigma += kr2*bonebyarg2; // correction at J1()/() 834 835 sigma *= damp2; // *rad*rad; 836 837 return sigma; 838 } 839 840 841 //////////////////////////////////////////////////////////////////////////// 842 // 810 843 // return differential elastic probability 2*pi*sin(theta)*d(probability)/d(omega) 811 844 812 845 G4double 813 G4DiffuseElastic::GetIntegrandFunction( G4double theta )846 G4DiffuseElastic::GetIntegrandFunction( G4double alpha ) 814 847 { 815 848 G4double result; 816 849 817 result = 2*pi*std::sin(theta); 818 result *= GetDiffElasticSumProb(theta); 850 result = GetDiffElasticSumProbA(alpha); 851 852 // result *= 2*pi*std::sin(theta); 853 819 854 return result; 820 855 } … … 859 894 //////////////////////////////////////////////////////////////////////////// 860 895 // 861 // Return inv momentum transfer -t > 0 from initialisation table862 863 G4double G4DiffuseElastic::SampleTableT( const G4ParticleDefinition* aParticle, G4double p,864 G4double Z, G4double A)865 {866 G4double theta = SampleTableThetaCMS( aParticle, p, Z, A); // sample theta in cms867 G4double t = 2*p*p*( 1 - std::cos(theta) ); // -t !!!868 return t;869 }870 871 ////////////////////////////////////////////////////////////////////////////872 //873 896 // Return scattering angle sampled in cms 874 897 … … 922 945 } 923 946 924 //////////////////////////////////////////////////////////////////////////// 925 // 926 // Return scattering angle sampled in cms according to precalculated table. 947 ///////////////////////////////////////////////////////////////////////////// 948 ///////////////////// Table preparation and reading //////////////////////// 949 //////////////////////////////////////////////////////////////////////////// 950 // 951 // Return inv momentum transfer -t > 0 from initialisation table 952 953 G4double G4DiffuseElastic::SampleTableT( const G4ParticleDefinition* aParticle, G4double p, 954 G4double Z, G4double A) 955 { 956 G4double alpha = SampleTableThetaCMS( aParticle, p, Z, A); // sample theta2 in cms 957 // G4double t = 2*p*p*( 1 - std::cos(std::sqrt(alpha)) ); // -t !!! 958 G4double t = p*p*alpha; // -t !!! 959 return t; 960 } 961 962 //////////////////////////////////////////////////////////////////////////// 963 // 964 // Return scattering angle2 sampled in cms according to precalculated table. 927 965 928 966 … … 942 980 if ( iElement == fElementNumberVector.size() ) 943 981 { 944 InitialiseOnFly(Z,A); 982 InitialiseOnFly(Z,A); // table preparation, if needed 983 945 984 // iElement--; 946 985 … … 955 994 G4double kinE = std::sqrt(momentum*momentum + m1*m1) - m1; 956 995 957 for( iMomentum = 0; iMomentum < fEnergyBin; iMomentum++)996 for( iMomentum = 0; iMomentum < fEnergyBin; iMomentum++) 958 997 { 959 998 if( kinE < fEnergyVector->GetLowEdgeEnergy(iMomentum) ) break; 960 999 } 961 if ( iMomentum == fEnergyBin ) iMomentum--; // kinE is more then theMaxEnergy1000 if ( iMomentum >= fEnergyBin ) iMomentum = fEnergyBin-1; // kinE is more then theMaxEnergy 962 1001 if ( iMomentum < 0 ) iMomentum = 0; // against negative index, kinE < theMinEnergy 1002 963 1003 // G4cout<<"iMomentum = "<<iMomentum<<G4endl; 964 1004 965 if (iMomentum == fEnergyBin -1 || iMomentum == 0 ) // the table edges1005 if (iMomentum == fEnergyBin -1 || iMomentum == 0 ) // the table edges 966 1006 { 967 1007 position = (*(*fAngleTable)(iMomentum))(fAngleBin-2)*G4UniformRand(); 1008 968 1009 // G4cout<<"position = "<<position<<G4endl; 969 1010 970 for(iAngle = 0; iAngle < fAngleBin ; iAngle++)1011 for(iAngle = 0; iAngle < fAngleBin-1; iAngle++) 971 1012 { 972 1013 if( position < (*(*fAngleTable)(iMomentum))(iAngle) ) break; 973 1014 } 974 if (iAngle == fAngleBin) iAngle--; 1015 if (iAngle >= fAngleBin-1) iAngle = fAngleBin-2; 1016 975 1017 // G4cout<<"iAngle = "<<iAngle<<G4endl; 976 1018 977 1019 randAngle = GetScatteringAngle(iMomentum, iAngle, position); 1020 978 1021 // G4cout<<"randAngle = "<<randAngle<<G4endl; 979 1022 } 980 else 981 { 982 position = (*(*fAngleTable)(iMomentum))(fAngleBin-2)*G4UniformRand(); 1023 else // kinE inside between energy table edges 1024 { 1025 // position = (*(*fAngleTable)(iMomentum))(fAngleBin-2)*G4UniformRand(); 1026 position = (*(*fAngleTable)(iMomentum))(0)*G4UniformRand(); 1027 983 1028 // G4cout<<"position = "<<position<<G4endl; 984 1029 985 for(iAngle = 0; iAngle < fAngleBin ; iAngle++)1030 for(iAngle = 0; iAngle < fAngleBin-1; iAngle++) 986 1031 { 987 if( position < (*(*fAngleTable)(iMomentum))(iAngle) ) break; 1032 // if( position < (*(*fAngleTable)(iMomentum))(iAngle) ) break; 1033 if( position > (*(*fAngleTable)(iMomentum))(iAngle) ) break; 988 1034 } 989 if (iAngle == fAngleBin) iAngle--; 1035 if (iAngle >= fAngleBin-1) iAngle = fAngleBin-2; 1036 990 1037 // G4cout<<"iAngle = "<<iAngle<<G4endl; 991 1038 992 1039 theta2 = GetScatteringAngle(iMomentum, iAngle, position); 1040 993 1041 // G4cout<<"theta2 = "<<theta2<<G4endl; 994 1042 E2 = fEnergyVector->GetLowEdgeEnergy(iMomentum); 1043 995 1044 // G4cout<<"E2 = "<<E2<<G4endl; 996 1045 997 1046 iMomentum--; 998 999 position = (*(*fAngleTable)(iMomentum))(fAngleBin-2)*G4UniformRand(); 1047 1048 // position = (*(*fAngleTable)(iMomentum))(fAngleBin-2)*G4UniformRand(); 1049 1000 1050 // G4cout<<"position = "<<position<<G4endl; 1001 1051 1002 for(iAngle = 0; iAngle < fAngleBin ; iAngle++)1052 for(iAngle = 0; iAngle < fAngleBin-1; iAngle++) 1003 1053 { 1004 if( position < (*(*fAngleTable)(iMomentum))(iAngle) ) break; 1054 // if( position < (*(*fAngleTable)(iMomentum))(iAngle) ) break; 1055 if( position > (*(*fAngleTable)(iMomentum))(iAngle) ) break; 1005 1056 } 1006 if (iAngle == fAngleBin) iAngle--;1007 1057 if (iAngle >= fAngleBin-1) iAngle = fAngleBin-2; 1058 1008 1059 theta1 = GetScatteringAngle(iMomentum, iAngle, position); 1060 1009 1061 // G4cout<<"theta1 = "<<theta1<<G4endl; 1062 1010 1063 E1 = fEnergyVector->GetLowEdgeEnergy(iMomentum); 1064 1011 1065 // G4cout<<"E1 = "<<E1<<G4endl; 1012 1066 … … 1016 1070 1017 1071 randAngle = W1*theta1 + W2*theta2; 1072 1073 // randAngle = theta2; 1018 1074 // G4cout<<"randAngle = "<<randAngle<<G4endl; 1019 1075 } 1076 // G4double angle = randAngle; 1077 // if (randAngle > 0.) randAngle /= 2*pi*std::sin(angle); 1078 1020 1079 return randAngle; 1021 1080 } 1022 1081 1082 ////////////////////////////////////////////////////////////////////////////// 1083 // 1084 // Initialisation for given particle on fly using new element number 1085 1086 void G4DiffuseElastic::InitialiseOnFly(G4double Z, G4double A) 1087 { 1088 fAtomicNumber = Z; // atomic number 1089 fAtomicWeight = A; // number of nucleons 1090 1091 fNuclearRadius = CalculateNuclearRad(fAtomicWeight); 1092 1093 if( verboseLevel > 0 ) 1094 { 1095 G4cout<<"G4DiffuseElastic::Initialise() the element with Z = " 1096 <<Z<<"; and A = "<<A<<G4endl; 1097 } 1098 fElementNumberVector.push_back(fAtomicNumber); 1099 1100 BuildAngleTable(); 1101 1102 fAngleBank.push_back(fAngleTable); 1103 1104 return; 1105 } 1106 1107 /////////////////////////////////////////////////////////////////////////////// 1108 // 1109 // Build for given particle and element table of momentum, angle probability. 1110 // For the moment in lab system. 1111 1112 void G4DiffuseElastic::BuildAngleTable() 1113 { 1114 G4int i, j; 1115 G4double partMom, kinE, a = 0., z = fParticle->GetPDGCharge(), m1 = fParticle->GetPDGMass(); 1116 G4double alpha1, alpha2, alphaMax, alphaCoulomb, delta = 0., sum = 0.; 1117 1118 G4Integrator<G4DiffuseElastic,G4double(G4DiffuseElastic::*)(G4double)> integral; 1119 1120 fAngleTable = new G4PhysicsTable(fEnergyBin); 1121 1122 for( i = 0; i < fEnergyBin; i++) 1123 { 1124 kinE = fEnergyVector->GetLowEdgeEnergy(i); 1125 partMom = std::sqrt( kinE*(kinE + 2*m1) ); 1126 1127 fWaveVector = partMom/hbarc; 1128 1129 G4double kR = fWaveVector*fNuclearRadius; 1130 G4double kR2 = kR*kR; 1131 G4double kRmax = 18.6; // 10.6; 10.6, 18, 10.174; ~ 3 maxima of J1 or 15., 25. 1132 G4double kRcoul = 1.9; // 1.2; 1.4, 2.5; // on the first slope of J1 1133 // G4double kRlim = 1.2; 1134 // G4double kRlim2 = kRlim*kRlim/kR2; 1135 1136 alphaMax = kRmax*kRmax/kR2; 1137 1138 if (alphaMax > 4.) alphaMax = 4.; // vmg05-02-09: was pi2 1139 1140 alphaCoulomb = kRcoul*kRcoul/kR2; 1141 1142 if( z ) 1143 { 1144 a = partMom/m1; // beta*gamma for m1 1145 fBeta = a/std::sqrt(1+a*a); 1146 fZommerfeld = CalculateZommerfeld( fBeta, z, fAtomicNumber); 1147 fAm = CalculateAm( partMom, fZommerfeld, fAtomicNumber); 1148 } 1149 G4PhysicsFreeVector* angleVector = new G4PhysicsFreeVector(fAngleBin-1); 1150 1151 // G4PhysicsLogVector* angleBins = new G4PhysicsLogVector( 0.001*alphaMax, alphaMax, fAngleBin ); 1152 1153 G4double delth = alphaMax/fAngleBin; 1154 1155 sum = 0.; 1156 1157 // fAddCoulomb = false; 1158 fAddCoulomb = true; 1159 1160 // for(j = 1; j < fAngleBin; j++) 1161 for(j = fAngleBin-1; j >= 1; j--) 1162 { 1163 // alpha1 = angleBins->GetLowEdgeEnergy(j-1); 1164 // alpha2 = angleBins->GetLowEdgeEnergy(j); 1165 1166 // alpha1 = alphaMax*(j-1)/fAngleBin; 1167 // alpha2 = alphaMax*( j )/fAngleBin; 1168 1169 alpha1 = delth*(j-1); 1170 // if(alpha1 < kRlim2) alpha1 = kRlim2; 1171 alpha2 = alpha1 + delth; 1172 1173 // if( ( alpha2 > alphaCoulomb ) && z ) fAddCoulomb = true; 1174 if( ( alpha1 < alphaCoulomb ) && z ) fAddCoulomb = false; 1175 1176 delta = integral.Legendre10(this, &G4DiffuseElastic::GetIntegrandFunction, alpha1, alpha2); 1177 // delta = integral.Legendre96(this, &G4DiffuseElastic::GetIntegrandFunction, alpha1, alpha2); 1178 1179 sum += delta; 1180 1181 angleVector->PutValue( j-1 , alpha1, sum ); // alpha2 1182 // G4cout<<"j-1 = "<<j-1<<"; alpha2 = "<<alpha2<<"; sum = "<<sum<<G4endl; 1183 } 1184 fAngleTable->insertAt(i,angleVector); 1185 1186 // delete[] angleVector; 1187 // delete[] angleBins; 1188 } 1189 return; 1190 } 1191 1023 1192 ///////////////////////////////////////////////////////////////////////////////// 1024 1193 // … … 1026 1195 1027 1196 G4double 1028 G4DiffuseElastic:: GetScatteringAngle( G4int iMomentum, G4int iAngle, G4double position)1197 G4DiffuseElastic:: GetScatteringAngle( G4int iMomentum, G4int iAngle, G4double position ) 1029 1198 { 1030 1199 G4double x1, x2, y1, y2, randAngle; … … 1033 1202 { 1034 1203 randAngle = (*fAngleTable)(iMomentum)->GetLowEdgeEnergy(iAngle); 1204 // iAngle++; 1035 1205 } 1036 1206 else … … 1046 1216 x2 = (*fAngleTable)(iMomentum)->GetLowEdgeEnergy(iAngle); 1047 1217 1048 if ( x1 == x2 ) 1218 if ( x1 == x2 ) randAngle = x2; 1049 1219 else 1050 1220 { 1051 if ( y1 == y2 ) randAngle = x1 + (x2 - x1)*G4UniformRand();1221 if ( y1 == y2 ) randAngle = x1 + ( x2 - x1 )*G4UniformRand(); 1052 1222 else 1053 1223 { 1054 randAngle = x1 + ( position - y1)*(x2 - x1)/(y2 - y1);1224 randAngle = x1 + ( position - y1 )*( x2 - x1 )/( y2 - y1 ); 1055 1225 } 1056 1226 } … … 1270 1440 } 1271 1441 1442 /////////////////////////////////////////////////////////////////////////////// 1443 // 1444 // Test for given particle and element table of momentum, angle probability. 1445 // For the moment in lab system. 1446 1447 void G4DiffuseElastic::TestAngleTable(const G4ParticleDefinition* theParticle, G4double partMom, 1448 G4double Z, G4double A) 1449 { 1450 fAtomicNumber = Z; // atomic number 1451 fAtomicWeight = A; // number of nucleons 1452 fNuclearRadius = CalculateNuclearRad(fAtomicWeight); 1453 1454 1455 1456 G4cout<<"G4DiffuseElastic::TestAngleTable() init the element with Z = " 1457 <<Z<<"; and A = "<<A<<G4endl; 1458 1459 fElementNumberVector.push_back(fAtomicNumber); 1460 1461 1462 1463 1464 G4int i=0, j; 1465 G4double a = 0., z = theParticle->GetPDGCharge(), m1 = fParticle->GetPDGMass(); 1466 G4double alpha1=0., alpha2=0., alphaMax=0., alphaCoulomb=0.; 1467 G4double deltaL10 = 0., deltaL96 = 0., deltaAG = 0.; 1468 G4double sumL10 = 0.,sumL96 = 0.,sumAG = 0.; 1469 G4double epsilon = 0.001; 1470 1471 G4Integrator<G4DiffuseElastic,G4double(G4DiffuseElastic::*)(G4double)> integral; 1472 1473 fAngleTable = new G4PhysicsTable(fEnergyBin); 1474 1475 fWaveVector = partMom/hbarc; 1476 1477 G4double kR = fWaveVector*fNuclearRadius; 1478 G4double kR2 = kR*kR; 1479 G4double kRmax = 10.6; // 10.6, 18, 10.174; ~ 3 maxima of J1 or 15., 25. 1480 G4double kRcoul = 1.2; // 1.4, 2.5; // on the first slope of J1 1481 1482 alphaMax = kRmax*kRmax/kR2; 1483 1484 if (alphaMax > 4.) alphaMax = 4.; // vmg05-02-09: was pi2 1485 1486 alphaCoulomb = kRcoul*kRcoul/kR2; 1487 1488 if( z ) 1489 { 1490 a = partMom/m1; // beta*gamma for m1 1491 fBeta = a/std::sqrt(1+a*a); 1492 fZommerfeld = CalculateZommerfeld( fBeta, z, fAtomicNumber); 1493 fAm = CalculateAm( partMom, fZommerfeld, fAtomicNumber); 1494 } 1495 G4PhysicsFreeVector* angleVector = new G4PhysicsFreeVector(fAngleBin-1); 1496 1497 // G4PhysicsLogVector* angleBins = new G4PhysicsLogVector( 0.001*alphaMax, alphaMax, fAngleBin ); 1498 1499 1500 fAddCoulomb = false; 1501 1502 for(j = 1; j < fAngleBin; j++) 1503 { 1504 // alpha1 = angleBins->GetLowEdgeEnergy(j-1); 1505 // alpha2 = angleBins->GetLowEdgeEnergy(j); 1506 1507 alpha1 = alphaMax*(j-1)/fAngleBin; 1508 alpha2 = alphaMax*( j )/fAngleBin; 1509 1510 if( ( alpha2 > alphaCoulomb ) && z ) fAddCoulomb = true; 1511 1512 deltaL10 = integral.Legendre10(this, &G4DiffuseElastic::GetIntegrandFunction, alpha1, alpha2); 1513 deltaL96 = integral.Legendre96(this, &G4DiffuseElastic::GetIntegrandFunction, alpha1, alpha2); 1514 deltaAG = integral.AdaptiveGauss(this, &G4DiffuseElastic::GetIntegrandFunction, 1515 alpha1, alpha2,epsilon); 1516 1517 // G4cout<<alpha1<<"\t"<<std::sqrt(alpha1)/degree<<"\t" 1518 // <<deltaL10<<"\t"<<deltaL96<<"\t"<<deltaAG<<G4endl; 1519 1520 sumL10 += deltaL10; 1521 sumL96 += deltaL96; 1522 sumAG += deltaAG; 1523 1524 G4cout<<alpha1<<"\t"<<std::sqrt(alpha1)/degree<<"\t" 1525 <<sumL10<<"\t"<<sumL96<<"\t"<<sumAG<<G4endl; 1526 1527 angleVector->PutValue( j-1 , alpha1, sumL10 ); // alpha2 1528 } 1529 fAngleTable->insertAt(i,angleVector); 1530 fAngleBank.push_back(fAngleTable); 1531 1532 /* 1533 // Integral over all angle range - Bad accuracy !!! 1534 1535 sumL10 = integral.Legendre10(this, &G4DiffuseElastic::GetIntegrandFunction, 0., alpha2); 1536 sumL96 = integral.Legendre96(this, &G4DiffuseElastic::GetIntegrandFunction, 0., alpha2); 1537 sumAG = integral.AdaptiveGauss(this, &G4DiffuseElastic::GetIntegrandFunction, 1538 0., alpha2,epsilon); 1539 G4cout<<G4endl; 1540 G4cout<<alpha2<<"\t"<<std::sqrt(alpha2)/degree<<"\t" 1541 <<sumL10<<"\t"<<sumL96<<"\t"<<sumAG<<G4endl; 1542 */ 1543 return; 1544 } 1545 1272 1546 // 1273 1547 //
Note: See TracChangeset
for help on using the changeset viewer.