| 1 | // |
|---|
| 2 | // ******************************************************************** |
|---|
| 3 | // * License and Disclaimer * |
|---|
| 4 | // * * |
|---|
| 5 | // * The Geant4 software is copyright of the Copyright Holders of * |
|---|
| 6 | // * the Geant4 Collaboration. It is provided under the terms and * |
|---|
| 7 | // * conditions of the Geant4 Software License, included in the file * |
|---|
| 8 | // * LICENSE and available at http://cern.ch/geant4/license . These * |
|---|
| 9 | // * include a list of copyright holders. * |
|---|
| 10 | // * * |
|---|
| 11 | // * Neither the authors of this software system, nor their employing * |
|---|
| 12 | // * institutes,nor the agencies providing financial support for this * |
|---|
| 13 | // * work make any representation or warranty, express or implied, * |
|---|
| 14 | // * regarding this software system or assume any liability for its * |
|---|
| 15 | // * use. Please see the license in the file LICENSE and URL above * |
|---|
| 16 | // * for the full disclaimer and the limitation of liability. * |
|---|
| 17 | // * * |
|---|
| 18 | // * This code implementation is the result of the scientific and * |
|---|
| 19 | // * technical work of the GEANT4 collaboration. * |
|---|
| 20 | // * By using, copying, modifying or distributing the software (or * |
|---|
| 21 | // * any work based on the software) you agree to acknowledge its * |
|---|
| 22 | // * use in resulting scientific publications, and indicate your * |
|---|
| 23 | // * acceptance of all terms of the Geant4 Software license. * |
|---|
| 24 | // ******************************************************************** |
|---|
| 25 | // |
|---|
| 26 | #include "G4InuclSpecialFunctions.hh" |
|---|
| 27 | |
|---|
| 28 | G4double G4InuclSpecialFunctions::bindingEnergyKummel(G4double A, |
|---|
| 29 | G4double Z) { |
|---|
| 30 | G4int verboseLevel = 2; |
|---|
| 31 | if (verboseLevel > 3) { |
|---|
| 32 | G4cout << " >>> G4InuclSpecialFunctions::bindingEnergyKummel" << G4endl; |
|---|
| 33 | } |
|---|
| 34 | |
|---|
| 35 | // calculates the nuclei binding energy using Kummel mass formula |
|---|
| 36 | |
|---|
| 37 | // data for Kummel mass formula |
|---|
| 38 | const G4double ASH[6] = {20.0, 28.0, 50.0, 82.0, 126.0, 184.0}; |
|---|
| 39 | const G4double OMN[3] = {2.511464e0, 7.640294e-2, 2.618317e-2}; |
|---|
| 40 | const G4double OMP[3] = {2.511464e0, 7.640294e-2, 6.243892e-3}; |
|---|
| 41 | const G4double TMET[5] = {9.477193e0, 7.166809e-1, 1.2169e-1, 4.07179e-2, 2.927847e-2}; |
|---|
| 42 | const G4double TMET1[4] = {-4.775384e-2, 5.880883e-5, -1.98964e-3, 3.427961e-6}; |
|---|
| 43 | |
|---|
| 44 | const G4double BET[5] = {1.535628e4, 1.502538e4, 1.481831e4, 1.459232e4, 1.40664e4}; |
|---|
| 45 | |
|---|
| 46 | const G4double SHN[5] = {8.0, 22.0, 32.0, 44.0, 58.0}; |
|---|
| 47 | |
|---|
| 48 | const G4double TNUN[5] = {8.728778e2, 8.45789e1, 9.210738e1, 5.98e1, 3.724e1}; |
|---|
| 49 | |
|---|
| 50 | const G4double TNUP[5] = {8.728778e2, 8.45789e1, 3.74e1, 5.42e1, 3.724e1}; |
|---|
| 51 | |
|---|
| 52 | const G4double TKSN[5] = {3.34582e1, 3.858216e0, -1.489218e-1, -9.2e-1, -6.36e-1}; |
|---|
| 53 | |
|---|
| 54 | const G4double TKSP[5] = {3.34582e1, 3.858216e0, 7.2e-1, -8.2e-1, -6.36e-1}; |
|---|
| 55 | |
|---|
| 56 | const G4double RON[5] = {-6.441735e3, -5.422058e3, -3.722645e3, -2.437172e4, -2.645867e4}; |
|---|
| 57 | |
|---|
| 58 | const G4double ROP[5] = {-6.75601e3, -5.877588e3, -3.216382e4, -3.010688e4, -2.645867e4}; |
|---|
| 59 | |
|---|
| 60 | const G4double SINGM[4] = {-3.578531e0, -2.6911e0, -7.487347e-1, 0.0}; |
|---|
| 61 | |
|---|
| 62 | const G4double C[2] = {6.1e3, 8.31e3}; |
|---|
| 63 | |
|---|
| 64 | const G4double AKU = 6.04122e5; |
|---|
| 65 | |
|---|
| 66 | const G4double US = 1.661835e4; |
|---|
| 67 | |
|---|
| 68 | const G4double UC = 6.218614e2; |
|---|
| 69 | |
|---|
| 70 | const G4double UT = 1.983151e4; |
|---|
| 71 | |
|---|
| 72 | const G4double SIPG = -2.067547e0; |
|---|
| 73 | |
|---|
| 74 | const G4double TAU = 2.231241e0; |
|---|
| 75 | |
|---|
| 76 | const G4double AL0 = 0.151; |
|---|
| 77 | |
|---|
| 78 | const G4double USB = 1.95114e4; |
|---|
| 79 | |
|---|
| 80 | const G4double UCB = 753.3; |
|---|
| 81 | |
|---|
| 82 | // const G4double ALD = 15.4941; |
|---|
| 83 | // const G4double ALD1 = 17.9439; |
|---|
| 84 | |
|---|
| 85 | // const G4double C3 = 0.7059; |
|---|
| 86 | |
|---|
| 87 | //const G4double C4 = 1.21129; |
|---|
| 88 | |
|---|
| 89 | // const G4double PKLD = 1.7826; |
|---|
| 90 | |
|---|
| 91 | const G4double DMU = 3.132902e4; |
|---|
| 92 | |
|---|
| 93 | G4double DM; |
|---|
| 94 | |
|---|
| 95 | G4double AN = A - Z; |
|---|
| 96 | |
|---|
| 97 | G4int INS = 0; |
|---|
| 98 | G4double ANE = 0.0; |
|---|
| 99 | G4double HNE = 0.0; |
|---|
| 100 | G4int i(0); |
|---|
| 101 | |
|---|
| 102 | for (i = 1; i < 6; i++) { |
|---|
| 103 | |
|---|
| 104 | if (AN <= ASH[i]) { |
|---|
| 105 | ANE = AN - ASH[i - 1]; |
|---|
| 106 | HNE = ASH[i] - AN; |
|---|
| 107 | INS = i - 1; |
|---|
| 108 | |
|---|
| 109 | break; |
|---|
| 110 | }; |
|---|
| 111 | }; |
|---|
| 112 | |
|---|
| 113 | G4int IPS = 0; |
|---|
| 114 | |
|---|
| 115 | G4double APR = 0.0; |
|---|
| 116 | |
|---|
| 117 | G4double HPR = 0.0; |
|---|
| 118 | |
|---|
| 119 | for (i = 1; i < 6; i++) { |
|---|
| 120 | |
|---|
| 121 | if (Z <= ASH[i]) { |
|---|
| 122 | APR = Z - ASH[i - 1]; |
|---|
| 123 | HPR = ASH[i] - Z; |
|---|
| 124 | IPS = i - 1; |
|---|
| 125 | |
|---|
| 126 | break; |
|---|
| 127 | }; |
|---|
| 128 | }; |
|---|
| 129 | |
|---|
| 130 | G4double PPHN = ANE * HNE; |
|---|
| 131 | G4double PPHZ = APR * HPR; |
|---|
| 132 | |
|---|
| 133 | // omega terms |
|---|
| 134 | G4double OMT = 0.0; |
|---|
| 135 | |
|---|
| 136 | if (INS <= 2) OMT += OMN[INS] * PPHN; |
|---|
| 137 | |
|---|
| 138 | if (IPS <= 2) OMT += OMP[IPS] * PPHZ; |
|---|
| 139 | |
|---|
| 140 | if (verboseLevel > 3) { |
|---|
| 141 | G4cout << " OMT " << OMT << G4endl; |
|---|
| 142 | } |
|---|
| 143 | |
|---|
| 144 | // theta term 1 |
|---|
| 145 | |
|---|
| 146 | G4double TET = 0.0; |
|---|
| 147 | |
|---|
| 148 | if (verboseLevel> 3) { |
|---|
| 149 | G4cout << " PPHN " << PPHN << " PPHZ " << PPHZ << G4endl; |
|---|
| 150 | G4cout << " INS " << INS << " IPS " << IPS << G4endl; |
|---|
| 151 | } |
|---|
| 152 | |
|---|
| 153 | if (PPHN > 0.5 && PPHZ > 0.5) { |
|---|
| 154 | G4int IT = 0; |
|---|
| 155 | |
|---|
| 156 | if ((INS + 1) * (IPS + 1) > 0) { |
|---|
| 157 | |
|---|
| 158 | if (INS - IPS == 1) { |
|---|
| 159 | |
|---|
| 160 | if (INS < 5) { |
|---|
| 161 | IT = INS; |
|---|
| 162 | |
|---|
| 163 | } else { |
|---|
| 164 | IT = -1; |
|---|
| 165 | }; |
|---|
| 166 | |
|---|
| 167 | } else { |
|---|
| 168 | IT = -1; |
|---|
| 169 | }; |
|---|
| 170 | }; |
|---|
| 171 | |
|---|
| 172 | if (verboseLevel > 3){ |
|---|
| 173 | G4cout << " IT " << IT << G4endl; |
|---|
| 174 | } |
|---|
| 175 | |
|---|
| 176 | if (IT >= 0) TET = TMET[IT] * PPHN * PPHZ; |
|---|
| 177 | }; |
|---|
| 178 | |
|---|
| 179 | if (verboseLevel> 3) { |
|---|
| 180 | G4cout << " TET " << TET << G4endl; |
|---|
| 181 | } |
|---|
| 182 | |
|---|
| 183 | // theta term 2 |
|---|
| 184 | G4double TET1 = 0.0; |
|---|
| 185 | |
|---|
| 186 | if (IPS == 2) { |
|---|
| 187 | TET1 += TMET1[0] * PPHZ * PPHZ; |
|---|
| 188 | |
|---|
| 189 | if (INS == 3) TET1 += TMET1[1] * PPHN * PPHZ * HNE * HPR; |
|---|
| 190 | }; |
|---|
| 191 | |
|---|
| 192 | if (INS == 3) { |
|---|
| 193 | |
|---|
| 194 | G4double TVSP = PPHN * PPHN * HNE; |
|---|
| 195 | TET1 += TMET1[2] * TVSP; |
|---|
| 196 | TET1 += TMET1[3] * TVSP * PPHN; |
|---|
| 197 | }; |
|---|
| 198 | |
|---|
| 199 | if (verboseLevel > 3) { |
|---|
| 200 | G4cout << " TET1 " << TET1 << G4endl; |
|---|
| 201 | } |
|---|
| 202 | |
|---|
| 203 | // betta, nu, ksy terms |
|---|
| 204 | G4double TBET = 0.0; |
|---|
| 205 | |
|---|
| 206 | if (INS != 0) { |
|---|
| 207 | |
|---|
| 208 | for (G4int i = 0; i <= INS - 1; i++) TBET += BET[i] * SHN[i]; |
|---|
| 209 | }; |
|---|
| 210 | TBET += BET[INS] * ANE; |
|---|
| 211 | |
|---|
| 212 | if (IPS != 0) { |
|---|
| 213 | |
|---|
| 214 | for (G4int i = 0; i <= IPS - 1; i++) TBET += BET[i] * SHN[i]; |
|---|
| 215 | }; |
|---|
| 216 | TBET += BET[IPS] * APR; |
|---|
| 217 | |
|---|
| 218 | if (verboseLevel> 3) { |
|---|
| 219 | G4cout << " TBET " << TBET << G4endl; |
|---|
| 220 | } |
|---|
| 221 | |
|---|
| 222 | G4double TBET1 = 0.0; |
|---|
| 223 | |
|---|
| 224 | if (PPHN > 0.0 || PPHZ > 0.0) |
|---|
| 225 | TBET1 = 0.5 * ((TNUN[INS] + TKSN[INS] * ANE) * PPHN + |
|---|
| 226 | (TNUP[IPS] + TKSP[IPS] * APR) * PPHZ); |
|---|
| 227 | TBET -= TBET1; |
|---|
| 228 | |
|---|
| 229 | if (verboseLevel > 3) { |
|---|
| 230 | G4cout << " TBET1 " << TBET1 << G4endl; |
|---|
| 231 | } |
|---|
| 232 | |
|---|
| 233 | // deformation |
|---|
| 234 | G4double TDEF = 0.0; |
|---|
| 235 | G4double X = Z * Z / A; |
|---|
| 236 | G4double X1 = std::pow(A, 0.3333333); |
|---|
| 237 | G4double X2 = X1 * X1; |
|---|
| 238 | |
|---|
| 239 | if (IPS != INS && INS >= 3 && IPS >= 2) { |
|---|
| 240 | G4double X3 = 2.0 * USB - UCB * X; |
|---|
| 241 | G4double DNZ = 0.0; |
|---|
| 242 | |
|---|
| 243 | if (!(INS < 3 || (INS - IPS) != 1)) DNZ = C[INS - 3]; |
|---|
| 244 | DNZ = TBET1 - DNZ; |
|---|
| 245 | G4double X4 = 0.2 * X2 * AL0 * AL0 * X3; |
|---|
| 246 | |
|---|
| 247 | if (DNZ > X4) { |
|---|
| 248 | G4double X5 = USB + X * UCB; |
|---|
| 249 | G4double X6 = std::log(DNZ / X4); |
|---|
| 250 | G4double X7 = std::sqrt(X6); |
|---|
| 251 | |
|---|
| 252 | // G4double ALM = AL0 * (X7 + 0.143 * AL0 * X5 / X3); |
|---|
| 253 | TDEF = -X4 * (X6 + 1.0) + DNZ + 0.038 * X2 * std::pow(AL0 * X7, 3.0) * X5; |
|---|
| 254 | }; |
|---|
| 255 | }; |
|---|
| 256 | |
|---|
| 257 | if (verboseLevel > 3) { |
|---|
| 258 | G4cout << " TDEF " << TDEF << G4endl; |
|---|
| 259 | } |
|---|
| 260 | |
|---|
| 261 | // pairing |
|---|
| 262 | G4double TPE = 0.0; |
|---|
| 263 | G4double TPEN = 0.0; |
|---|
| 264 | G4double TPEP = 0.0; |
|---|
| 265 | G4double DN0 = 0.0; |
|---|
| 266 | G4double DZ0 = 0.0; |
|---|
| 267 | G4double AV = 2.0 * int(0.5 * AN + 0.1); |
|---|
| 268 | |
|---|
| 269 | if (AN > AV) DN0 = 1.0; |
|---|
| 270 | |
|---|
| 271 | AV = 2.0 * int(0.5 * Z + 0.1); |
|---|
| 272 | |
|---|
| 273 | if (Z > AV) DZ0 = 1.0; |
|---|
| 274 | |
|---|
| 275 | if (DN0 * DZ0 > 0.0) TPE = DMU / A; |
|---|
| 276 | |
|---|
| 277 | if (DN0 + DZ0 > 0.0) { |
|---|
| 278 | |
|---|
| 279 | if (DN0 > 0.0) { |
|---|
| 280 | TPEN = RON[INS] / X1; |
|---|
| 281 | |
|---|
| 282 | if (INS >= 1) { |
|---|
| 283 | |
|---|
| 284 | if (AN >= 82.0) { |
|---|
| 285 | TPEN /= X1; |
|---|
| 286 | |
|---|
| 287 | if (INS == 3) { |
|---|
| 288 | G4double PN = 0.0; |
|---|
| 289 | G4double HN = 0.0; |
|---|
| 290 | |
|---|
| 291 | if (AN > 90.0) PN = AN - 90.0; |
|---|
| 292 | |
|---|
| 293 | if (AN < 116.0) HN = 116.0 - AN; |
|---|
| 294 | TPEN += TAU * PN * HN; |
|---|
| 295 | }; |
|---|
| 296 | TPEN += SINGM[INS - 1] * PPHN; |
|---|
| 297 | }; |
|---|
| 298 | }; |
|---|
| 299 | }; |
|---|
| 300 | |
|---|
| 301 | if (DZ0 > 0.0) { |
|---|
| 302 | TPEP = ROP[IPS] / X1; |
|---|
| 303 | |
|---|
| 304 | if (IPS == 1) TPEP += SIPG * PPHZ; |
|---|
| 305 | |
|---|
| 306 | if (IPS > 1) TPEP = TPEP / X1; |
|---|
| 307 | |
|---|
| 308 | }; |
|---|
| 309 | }; |
|---|
| 310 | TPE += TPEN + TPEP; |
|---|
| 311 | |
|---|
| 312 | if (verboseLevel > 3) { |
|---|
| 313 | G4cout << " TPE " << TPE << " TPEN " << TPEN << " TPEP " << TPEP << G4endl; |
|---|
| 314 | } |
|---|
| 315 | |
|---|
| 316 | // collect everything |
|---|
| 317 | DM = (AKU - US * X2 + Z * (Z - 1.0) / X1 * (OMT - UC) - UT / A * |
|---|
| 318 | (AN - Z) * (AN - Z) + TET + TET1 + TBET + TDEF + TPE) * 0.001; |
|---|
| 319 | |
|---|
| 320 | if (verboseLevel > 3) { |
|---|
| 321 | G4cout << " kummel " << G4endl; |
|---|
| 322 | } |
|---|
| 323 | |
|---|
| 324 | return DM; |
|---|
| 325 | } |
|---|
