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
00033
00034
00035
00036
00037
00038
00039
00040
00041
00042
00043
00044
00045
00046
00047
00048
00049
00050
00051
00052
00053
00054
00055
00056
00057
00058
00059
00060
00061
00062
00063
00064
00066
00067 #include "G4NuclearAbrasionGeometry.hh"
00068 #include "G4WilsonRadius.hh"
00069 #include "G4PhysicalConstants.hh"
00070 #include "G4SystemOfUnits.hh"
00072
00073 G4NuclearAbrasionGeometry::G4NuclearAbrasionGeometry (G4double AP1,
00074 G4double AT1, G4double r1)
00075 {
00076
00077
00078
00079
00080 G4WilsonRadius aR;
00081 AP = AP1;
00082 AT = AT1;
00083 rP = aR.GetWilsonRadius(AP);
00084 rT = aR.GetWilsonRadius(AT);
00085 r = r1;
00086 n = rP / (rP + rT);
00087 b = r / (rP + rT);
00088 m = rT / rP;
00089 Q = (1.0 - b)/n;
00090 S = Q * Q;
00091 T = S * Q;
00092 R = std::sqrt(m*n);
00093 U = 1.0/m - 2.0;
00094
00095
00096
00097
00098
00099 rth = 2.0/3.0;
00100 B = 10.0 * MeV;
00101 }
00103
00104 G4NuclearAbrasionGeometry::~G4NuclearAbrasionGeometry ()
00105 {;}
00107
00108 void G4NuclearAbrasionGeometry::SetPeripheralThreshold (G4double rth1)
00109 {if (rth1 > 0.0 && rth1 <= 1.0) rth = rth1;}
00111
00112 G4double G4NuclearAbrasionGeometry::GetPeripheralThreshold ()
00113 {return rth;}
00115
00116 G4double G4NuclearAbrasionGeometry::P ()
00117 {
00118
00119
00120
00121
00122
00123
00124 G4double valueP = 0.0;
00125
00126 if (rT > rP)
00127 {
00128 if (rT-rP<=r && r<=rT+rP) valueP = 0.125*R*U*S - 0.125*(0.5*R*U+1.0)*T;
00129 else valueP = -1.0;
00130 }
00131 else
00132 {
00133 if (rP-rT<=r && r<=rP+rT) valueP = 0.125*R*U*S - 0.125*(0.5*std::sqrt(n/m)*U-
00134 (std::sqrt(1.0-m*m)/n - 1.0)*std::sqrt((2.0-m)/std::pow(m,5.0)))*T;
00135 else valueP = (std::sqrt(1.0-m*m)/n-1.0)*std::sqrt(1.0-b*b/n/n);
00136 }
00137
00138 if (!(valueP <= 1.0 && valueP>= -1.0))
00139 {
00140 if (valueP > 1.0) valueP = 1.0;
00141 else valueP = -1.0;
00142 }
00143 return valueP;
00144 }
00146
00147 G4double G4NuclearAbrasionGeometry::F ()
00148 {
00149
00150
00151
00152
00153
00154
00155 G4double valueF = 0.0;
00156
00157 if (rT > rP)
00158 {
00159 if (rT-rP<=r && r<=rT+rP) valueF = 0.75*R*S - 0.125*(3.0*R-1.0)*T;
00160 else valueF = 1.0;
00161 }
00162 else
00163 {
00164 if (rP-rT<=r && r<=rP+rT) valueF = 0.75*R*S - 0.125*(3.0*std::sqrt(n/m)-
00165 (1.0-std::pow(1.0-m*m,3.0/2.0))*std::sqrt(1.0-std::pow(1.0-m,2.0))/std::pow(m,3.0))*T;
00166 else valueF = (1.0-std::pow(1.0-m*m,3.0/2.0))*std::sqrt(1.0-b*b/n/n);
00167 }
00168
00169 if (!(valueF <= 1.0 && valueF>= 0.0))
00170 {
00171 if (valueF > 1.0) valueF = 1.0;
00172 else valueF = 0.0;
00173 }
00174 return valueF;
00175 }
00177
00178 G4double G4NuclearAbrasionGeometry::GetExcitationEnergyOfProjectile ()
00179 {
00180 G4double F1 = F();
00181 G4double P1 = P();
00182 G4double Es = 0.0;
00183
00184 Es = 0.95 * MeV * 4.0 * pi * rP*rP/fermi/fermi *
00185 (1.0+P1-std::pow(1.0-F1,2.0/3.0));
00186
00187 if ((r-rP)/rT < rth)
00188 {
00189 G4double omega = 0.0;
00190 if (AP < 12.0) omega = 1500.0;
00191 else if (AP <= 16.0) omega = 1500.0 - 320.0*(AP-12.0);
00192 Es *= 1.0 + F1*(5.0+omega*F1*F1);
00193 }
00194
00195 if (Es < 0.0)
00196 Es = 0.0;
00197 else if (Es > B * AP)
00198 Es = B * AP;
00199 return Es;
00200 }
00201
00202
00203 G4double G4NuclearAbrasionGeometry::GetExcitationEnergyOfTarget ()
00204 {
00205
00206
00207
00208
00209
00210 G4NuclearAbrasionGeometry* revAbrasionGeometry =
00211 new G4NuclearAbrasionGeometry(AT, AP, r);
00212 G4double F1 = revAbrasionGeometry->F();
00213 G4double P1 = revAbrasionGeometry->P();
00214 G4double Es = 0.0;
00215
00216 Es = 0.95 * MeV * 4.0 * pi * rT*rT/fermi/fermi *
00217 (1.0+P1-std::pow(1.0-F1,2.0/3.0));
00218
00219
00220 if ((r-rT)/rP < rth) {
00221 G4double omega = 0.0;
00222 if (AT < 12.0) omega = 1500.0;
00223 else if (AT <= 16.0) omega = 1500.0 - 320.0*(AT-12.0);
00224 Es *= 1.0 + F1*(5.0+omega*F1*F1);
00225 }
00226
00227 if (Es < 0.0)
00228 Es = 0.0;
00229 else if (Es > B * AT)
00230 Es = B * AT;
00231
00232 delete revAbrasionGeometry;
00233
00234 return Es;
00235 }