// // ******************************************************************** // * License and Disclaimer * // * * // * The Geant4 software is copyright of the Copyright Holders of * // * the Geant4 Collaboration. It is provided under the terms and * // * conditions of the Geant4 Software License, included in the file * // * LICENSE and available at http://cern.ch/geant4/license . These * // * include a list of copyright holders. * // * * // * Neither the authors of this software system, nor their employing * // * institutes,nor the agencies providing financial support for this * // * work make any representation or warranty, express or implied, * // * regarding this software system or assume any liability for its * // * use. Please see the license in the file LICENSE and URL above * // * for the full disclaimer and the limitation of liability. * // * * // * This code implementation is the result of the scientific and * // * technical work of the GEANT4 collaboration. * // * By using, copying, modifying or distributing the software (or * // * any work based on the software) you agree to acknowledge its * // * use in resulting scientific publications, and indicate your * // * acceptance of all terms of the Geant4 Software license. * // ******************************************************************** // // $Id: bindingEnergyKummel.cc,v 1.13 2010/06/25 09:45:14 gunter Exp $ // Geant4 tag: $Name: geant4-09-04-beta-01 $ // // 20100202 M. Kelsey -- Eliminate unnecessary use of std::pow() #include "G4InuclSpecialFunctions.hh" using namespace G4InuclSpecialFunctions; G4double G4InuclSpecialFunctions::bindingEnergyKummel(G4double A, G4double Z) { G4int verboseLevel = 2; if (verboseLevel > 3) { G4cout << " >>> G4InuclSpecialFunctions::bindingEnergyKummel" << G4endl; } // calculates the nuclei binding energy using Kummel mass formula // data for Kummel mass formula const G4double ASH[6] = {20.0, 28.0, 50.0, 82.0, 126.0, 184.0}; const G4double OMN[3] = {2.511464e0, 7.640294e-2, 2.618317e-2}; const G4double OMP[3] = {2.511464e0, 7.640294e-2, 6.243892e-3}; const G4double TMET[5] = {9.477193e0, 7.166809e-1, 1.2169e-1, 4.07179e-2, 2.927847e-2}; const G4double TMET1[4] = {-4.775384e-2, 5.880883e-5, -1.98964e-3, 3.427961e-6}; const G4double BET[5] = {1.535628e4, 1.502538e4, 1.481831e4, 1.459232e4, 1.40664e4}; const G4double SHN[5] = {8.0, 22.0, 32.0, 44.0, 58.0}; const G4double TNUN[5] = {8.728778e2, 8.45789e1, 9.210738e1, 5.98e1, 3.724e1}; const G4double TNUP[5] = {8.728778e2, 8.45789e1, 3.74e1, 5.42e1, 3.724e1}; const G4double TKSN[5] = {3.34582e1, 3.858216e0, -1.489218e-1, -9.2e-1, -6.36e-1}; const G4double TKSP[5] = {3.34582e1, 3.858216e0, 7.2e-1, -8.2e-1, -6.36e-1}; const G4double RON[5] = {-6.441735e3, -5.422058e3, -3.722645e3, -2.437172e4, -2.645867e4}; const G4double ROP[5] = {-6.75601e3, -5.877588e3, -3.216382e4, -3.010688e4, -2.645867e4}; const G4double SINGM[4] = {-3.578531e0, -2.6911e0, -7.487347e-1, 0.0}; const G4double C[2] = {6.1e3, 8.31e3}; const G4double AKU = 6.04122e5; const G4double US = 1.661835e4; const G4double UC = 6.218614e2; const G4double UT = 1.983151e4; const G4double SIPG = -2.067547e0; const G4double TAU = 2.231241e0; const G4double AL0 = 0.151; const G4double USB = 1.95114e4; const G4double UCB = 753.3; // const G4double ALD = 15.4941; // const G4double ALD1 = 17.9439; // const G4double C3 = 0.7059; //const G4double C4 = 1.21129; // const G4double PKLD = 1.7826; const G4double DMU = 3.132902e4; G4double DM; G4double AN = A - Z; G4int INS = 0; G4double ANE = 0.0; G4double HNE = 0.0; G4int i(0); for (i = 1; i < 6; i++) { if (AN <= ASH[i]) { ANE = AN - ASH[i - 1]; HNE = ASH[i] - AN; INS = i - 1; break; }; }; G4int IPS = 0; G4double APR = 0.0; G4double HPR = 0.0; for (i = 1; i < 6; i++) { if (Z <= ASH[i]) { APR = Z - ASH[i - 1]; HPR = ASH[i] - Z; IPS = i - 1; break; }; }; G4double PPHN = ANE * HNE; G4double PPHZ = APR * HPR; // omega terms G4double OMT = 0.0; if (INS <= 2) OMT += OMN[INS] * PPHN; if (IPS <= 2) OMT += OMP[IPS] * PPHZ; if (verboseLevel > 3) { G4cout << " OMT " << OMT << G4endl; } // theta term 1 G4double TET = 0.0; if (verboseLevel> 3) { G4cout << " PPHN " << PPHN << " PPHZ " << PPHZ << G4endl; G4cout << " INS " << INS << " IPS " << IPS << G4endl; } if (PPHN > 0.5 && PPHZ > 0.5) { G4int IT = 0; if ((INS + 1) * (IPS + 1) > 0) { if (INS - IPS == 1) { if (INS < 5) { IT = INS; } else { IT = -1; }; } else { IT = -1; }; }; if (verboseLevel > 3){ G4cout << " IT " << IT << G4endl; } if (IT >= 0) TET = TMET[IT] * PPHN * PPHZ; }; if (verboseLevel> 3) { G4cout << " TET " << TET << G4endl; } // theta term 2 G4double TET1 = 0.0; if (IPS == 2) { TET1 += TMET1[0] * PPHZ * PPHZ; if (INS == 3) TET1 += TMET1[1] * PPHN * PPHZ * HNE * HPR; }; if (INS == 3) { G4double TVSP = PPHN * PPHN * HNE; TET1 += TMET1[2] * TVSP; TET1 += TMET1[3] * TVSP * PPHN; }; if (verboseLevel > 3) { G4cout << " TET1 " << TET1 << G4endl; } // betta, nu, ksy terms G4double TBET = 0.0; if (INS != 0) { for (G4int i = 0; i <= INS - 1; i++) TBET += BET[i] * SHN[i]; }; TBET += BET[INS] * ANE; if (IPS != 0) { for (G4int i = 0; i <= IPS - 1; i++) TBET += BET[i] * SHN[i]; }; TBET += BET[IPS] * APR; if (verboseLevel> 3) { G4cout << " TBET " << TBET << G4endl; } G4double TBET1 = 0.0; if (PPHN > 0.0 || PPHZ > 0.0) TBET1 = 0.5 * ((TNUN[INS] + TKSN[INS] * ANE) * PPHN + (TNUP[IPS] + TKSP[IPS] * APR) * PPHZ); TBET -= TBET1; if (verboseLevel > 3) { G4cout << " TBET1 " << TBET1 << G4endl; } // deformation G4double TDEF = 0.0; G4double X = Z * Z / A; G4double X1 = G4cbrt(A); G4double X2 = X1 * X1; if (IPS != INS && INS >= 3 && IPS >= 2) { G4double X3 = 2.0 * USB - UCB * X; G4double DNZ = 0.0; if (!(INS < 3 || (INS - IPS) != 1)) DNZ = C[INS - 3]; DNZ = TBET1 - DNZ; G4double X4 = 0.2 * X2 * AL0 * AL0 * X3; if (DNZ > X4) { G4double X5 = USB + X * UCB; G4double X6 = std::log(DNZ / X4); G4double X7 = std::sqrt(X6); // G4double ALM = AL0 * (X7 + 0.143 * AL0 * X5 / X3); TDEF = -X4 * (X6 + 1.0) + DNZ + 0.038 * X2 * (AL0*X7 * AL0*X7 * AL0*X7) * X5; }; }; if (verboseLevel > 3) { G4cout << " TDEF " << TDEF << G4endl; } // pairing G4double TPE = 0.0; G4double TPEN = 0.0; G4double TPEP = 0.0; G4double DN0 = 0.0; G4double DZ0 = 0.0; G4double AV = 2.0 * int(0.5 * AN + 0.1); if (AN > AV) DN0 = 1.0; AV = 2.0 * int(0.5 * Z + 0.1); if (Z > AV) DZ0 = 1.0; if (DN0 * DZ0 > 0.0) TPE = DMU / A; if (DN0 + DZ0 > 0.0) { if (DN0 > 0.0) { TPEN = RON[INS] / X1; if (INS >= 1) { if (AN >= 82.0) { TPEN /= X1; if (INS == 3) { G4double PN = 0.0; G4double HN = 0.0; if (AN > 90.0) PN = AN - 90.0; if (AN < 116.0) HN = 116.0 - AN; TPEN += TAU * PN * HN; }; TPEN += SINGM[INS - 1] * PPHN; }; }; }; if (DZ0 > 0.0) { TPEP = ROP[IPS] / X1; if (IPS == 1) TPEP += SIPG * PPHZ; if (IPS > 1) TPEP = TPEP / X1; }; }; TPE += TPEN + TPEP; if (verboseLevel > 3) { G4cout << " TPE " << TPE << " TPEN " << TPEN << " TPEP " << TPEP << G4endl; } // collect everything DM = (AKU - US * X2 + Z * (Z - 1.0) / X1 * (OMT - UC) - UT / A * (AN - Z) * (AN - Z) + TET + TET1 + TBET + TDEF + TPE) * 0.001; if (verboseLevel > 3) { G4cout << " kummel " << G4endl; } return DM; }