00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032 #include "G4NeutronHPKallbachMannSyst.hh"
00033 #include "G4SystemOfUnits.hh"
00034 #include "Randomize.hh"
00035 #include "G4HadronicException.hh"
00036
00037 G4double G4NeutronHPKallbachMannSyst::Sample(G4double anEnergy)
00038 {
00039 G4double result;
00040
00041 G4double zero = GetKallbachZero(anEnergy);
00042 if(zero>1) zero=1.;
00043 if(zero<-1)zero=-1.;
00044 G4double max = Kallbach(zero, anEnergy);
00045 double upper = Kallbach(1., anEnergy);
00046 double lower = Kallbach(-1., anEnergy);
00047 if(upper>max) max=upper;
00048 if(lower>max) max=lower;
00049 G4double value, random;
00050 do
00051 {
00052 result = 2.*G4UniformRand()-1;
00053 value = Kallbach(result, anEnergy)/max;
00054 random = G4UniformRand();
00055 }
00056 while(random>value);
00057
00058 return result;
00059 }
00060
00061 G4double G4NeutronHPKallbachMannSyst::Kallbach(G4double cosTh, G4double anEnergy)
00062 {
00063
00064 G4double result;
00065 G4double theX = A(anEnergy)*cosTh;
00066 result = 0.5*(std::exp( theX)*(1+theCompoundFraction)
00067 +std::exp(-theX)*(1-theCompoundFraction));
00068 return result;
00069 }
00070
00071 G4double G4NeutronHPKallbachMannSyst::GetKallbachZero(G4double anEnergy)
00072 {
00073 G4double result;
00074 if ( theCompoundFraction == 1 )
00075 {
00076
00077 theCompoundFraction *= (1-1.0e-15);
00078 }
00079 result = 0.5 * (1./A(anEnergy)) * std::log((1-theCompoundFraction)/(1+theCompoundFraction));
00080 return result;
00081 }
00082
00083 G4double G4NeutronHPKallbachMannSyst::A(G4double anEnergy)
00084 {
00085 G4double result;
00086 G4double C1 = 0.04/MeV;
00087 G4double C2 = 1.8E-6/(MeV*MeV*MeV);
00088 G4double C3 = 6.7E-7/(MeV*MeV*MeV*MeV);
00089
00090 G4double epsa = anEnergy*theTargetMass/(theTargetMass+theIncidentMass);
00091 G4int Ac = theTargetA+1;
00092 G4int Nc = Ac - theTargetZ;
00093 G4int AA = theTargetA;
00094 G4int ZA = theTargetZ;
00095 G4double ea = epsa+SeparationEnergy(Ac, Nc, AA, ZA);
00096 G4double Et1 = 130*MeV;
00097 G4double R1 = std::min(ea, Et1);
00098
00099 G4double epsb = theProductEnergy*(theProductMass+theResidualMass)/theResidualMass;
00100 G4int AB = theResidualA;
00101 G4int ZB = theResidualZ;
00102 G4double eb = epsb+SeparationEnergy(Ac, Nc, AB, ZB );
00103 G4double X1 = R1*eb/ea;
00104 G4double Et3 = 41*MeV;
00105 G4double R3 = std::min(ea, Et3);
00106 G4double X3 = R3*eb/ea;
00107 G4double Ma = 1;
00108 G4double mb(0);
00109 G4int productA = theTargetA+1-theResidualA;
00110 G4int productZ = theTargetZ-theResidualZ;
00111 if(productZ==0)
00112 {
00113 mb = 0.5;
00114 }
00115 else if(productZ==1)
00116 {
00117 mb = 1;
00118 }
00119 else if(productZ==2)
00120 {
00121 mb = 2;
00122 if(productA==3) mb=1;
00123 }
00124 else
00125 {
00126 throw G4HadronicException(__FILE__, __LINE__, "Severe error in the sampling of Kallbach-Mann Systematics");
00127 }
00128
00129 result = C1*X1 + C2*std::pow(X1, 3.) + C3*Ma*mb*std::pow(X3, 4.);
00130 return result;
00131 }
00132
00133 G4double G4NeutronHPKallbachMannSyst::SeparationEnergy(G4int Ac, G4int Nc, G4int AA, G4int ZA)
00134 {
00135 G4double result;
00136 G4int NA = AA-ZA;
00137 G4int Zc = Ac-Nc;
00138 result = 15.68*(Ac-AA);
00139 result += -28.07*((Nc-Zc)*(Nc-Zc)/Ac - (NA-ZA)*(NA-ZA)/AA);
00140 result += -18.56*(std::pow(G4double(Ac), 2./3.) - std::pow(G4double(AA), 2./3.));
00141 result += 33.22*((Nc-Zc)*(Nc-Zc)/std::pow(G4double(Ac), 4./3.) - (NA-ZA)*(NA-ZA)/std::pow(G4double(AA), 4./3.));
00142 result += -0.717*(Zc*Zc/std::pow(G4double(Ac),1./3.)-ZA*ZA/std::pow(G4double(AA),1./3.));
00143 result += 1.211*(Zc*Zc/Ac-ZA*ZA/AA);
00144 G4double totalBinding(0);
00145 G4int productA = theTargetA+1-theResidualA;
00146 G4int productZ = theTargetZ-theResidualZ;
00147 if(productZ==0&&productA==1) totalBinding=0;
00148 if(productZ==1&&productA==1) totalBinding=0;
00149 if(productZ==1&&productA==2) totalBinding=2.22;
00150 if(productZ==1&&productA==3) totalBinding=8.48;
00151 if(productZ==2&&productA==3) totalBinding=7.72;
00152 if(productZ==2&&productA==4) totalBinding=28.3;
00153 result += -totalBinding;
00154 result *= MeV;
00155 return result;
00156 }