source: trunk/source/processes/hadronic/models/cascade/cascade/src/bindingEnergyKummel.cc @ 1315

Last change on this file since 1315 was 1315, checked in by garnier, 14 years ago

update geant4-09-04-beta-cand-01 interfaces-V09-03-09 vis-V09-03-08
File size: 8.0 KB
Line 
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// $Id: bindingEnergyKummel.cc,v 1.12 2010/03/19 05:03:23 mkelsey Exp $
26// Geant4 tag: $Name: geant4-09-04-beta-cand-01 $
27//
28// 20100202  M. Kelsey -- Eliminate unnecessary use of std::pow()
29
30#include "G4InuclSpecialFunctions.hh"
31
32using namespace G4InuclSpecialFunctions;
33
34
35G4double G4InuclSpecialFunctions::bindingEnergyKummel(G4double A, 
36                                                      G4double Z) {
37  G4int verboseLevel = 2;
38  if (verboseLevel > 3) {
39    G4cout << " >>> G4InuclSpecialFunctions::bindingEnergyKummel" << G4endl;
40  }
41
42  // calculates the nuclei binding energy using Kummel mass formula
43
44// data for Kummel mass formula
45  const G4double ASH[6] = {20.0, 28.0, 50.0, 82.0, 126.0, 184.0};
46  const G4double OMN[3] = {2.511464e0, 7.640294e-2, 2.618317e-2};
47  const G4double OMP[3] = {2.511464e0, 7.640294e-2, 6.243892e-3};
48  const G4double TMET[5] = {9.477193e0, 7.166809e-1, 1.2169e-1, 4.07179e-2, 2.927847e-2}; 
49  const G4double TMET1[4] = {-4.775384e-2, 5.880883e-5, -1.98964e-3, 3.427961e-6}; 
50
51  const G4double BET[5] = {1.535628e4, 1.502538e4, 1.481831e4, 1.459232e4, 1.40664e4};
52
53  const G4double SHN[5] = {8.0, 22.0, 32.0, 44.0, 58.0};
54
55  const G4double TNUN[5] = {8.728778e2, 8.45789e1, 9.210738e1, 5.98e1, 3.724e1}; 
56
57  const G4double TNUP[5] = {8.728778e2, 8.45789e1, 3.74e1, 5.42e1, 3.724e1};
58
59  const G4double TKSN[5] = {3.34582e1, 3.858216e0, -1.489218e-1, -9.2e-1, -6.36e-1};
60
61  const G4double TKSP[5] = {3.34582e1, 3.858216e0, 7.2e-1, -8.2e-1, -6.36e-1};
62
63  const G4double RON[5] = {-6.441735e3, -5.422058e3, -3.722645e3, -2.437172e4, -2.645867e4};
64
65  const G4double ROP[5] = {-6.75601e3, -5.877588e3, -3.216382e4, -3.010688e4, -2.645867e4}; 
66
67  const G4double SINGM[4] = {-3.578531e0, -2.6911e0, -7.487347e-1, 0.0};
68
69  const G4double C[2] = {6.1e3, 8.31e3};
70
71  const G4double AKU = 6.04122e5;
72
73  const G4double US = 1.661835e4;
74
75  const G4double UC = 6.218614e2;
76
77  const G4double UT = 1.983151e4;
78
79  const G4double SIPG = -2.067547e0;
80
81  const G4double TAU = 2.231241e0;
82
83  const G4double AL0 = 0.151;
84
85  const G4double USB = 1.95114e4;
86
87  const G4double UCB = 753.3;
88
89  //  const G4double ALD = 15.4941;
90  // const G4double ALD1 = 17.9439;
91
92  // const G4double C3 = 0.7059;
93
94  //const G4double C4 = 1.21129;
95
96  //  const G4double PKLD = 1.7826;
97
98  const G4double DMU = 3.132902e4;
99       
100  G4double DM;
101
102  G4double AN = A - Z;
103
104  G4int INS = 0;
105  G4double ANE = 0.0;
106  G4double HNE = 0.0;
107  G4int i(0);
108
109  for (i = 1; i < 6; i++) {
110
111    if (AN <= ASH[i]) {
112      ANE = AN - ASH[i - 1];
113      HNE = ASH[i] - AN;
114      INS = i - 1;
115
116      break;
117    };
118  };
119 
120  G4int IPS = 0;
121
122  G4double APR = 0.0;
123
124  G4double HPR = 0.0;
125
126  for (i = 1; i < 6; i++) {
127
128    if (Z <= ASH[i]) {
129      APR = Z - ASH[i - 1];
130      HPR = ASH[i] - Z;
131      IPS = i - 1;
132
133      break;
134    };
135  };
136 
137  G4double PPHN = ANE * HNE;
138  G4double PPHZ = APR * HPR;
139
140  // omega terms
141  G4double OMT = 0.0;
142
143  if (INS <= 2) OMT += OMN[INS] * PPHN;
144
145  if (IPS <= 2) OMT += OMP[IPS] * PPHZ;
146
147  if (verboseLevel > 3) {
148    G4cout << " OMT " << OMT << G4endl;
149  }
150
151  // theta term 1
152
153  G4double TET = 0.0;
154
155  if (verboseLevel> 3) {
156    G4cout << " PPHN " << PPHN << " PPHZ " << PPHZ << G4endl;
157    G4cout << " INS " << INS << " IPS " << IPS << G4endl;
158  }
159
160  if (PPHN > 0.5 && PPHZ > 0.5) {
161    G4int IT = 0;
162
163    if ((INS + 1) * (IPS + 1) > 0) {
164
165      if (INS - IPS == 1) {
166
167        if (INS < 5) {
168          IT = INS;
169
170        } else {
171          IT = -1; 
172        };               
173
174      } else {
175        IT = -1;
176      }; 
177    };
178
179    if (verboseLevel > 3){
180      G4cout << " IT " << IT << G4endl; 
181    }
182
183    if (IT >= 0) TET = TMET[IT] * PPHN * PPHZ;
184  };
185
186  if (verboseLevel> 3) {
187    G4cout << " TET " << TET << G4endl;
188  }
189
190  // theta term 2
191  G4double TET1 = 0.0;
192
193  if (IPS == 2) {
194    TET1 += TMET1[0] * PPHZ * PPHZ;
195
196    if (INS == 3) TET1 += TMET1[1] * PPHN * PPHZ * HNE * HPR;
197  };
198
199  if (INS == 3) {
200
201    G4double TVSP = PPHN * PPHN * HNE;
202    TET1 += TMET1[2] * TVSP;
203    TET1 += TMET1[3] * TVSP * PPHN;
204  };
205
206  if (verboseLevel > 3) {
207    G4cout << " TET1 " << TET1 << G4endl;
208  }
209
210  // betta, nu, ksy terms
211  G4double TBET = 0.0;
212
213  if (INS != 0) {
214
215    for (G4int i = 0; i <= INS - 1; i++) TBET += BET[i] * SHN[i];
216  };
217  TBET += BET[INS] * ANE;
218       
219  if (IPS != 0) {
220
221    for (G4int i = 0; i <= IPS - 1; i++) TBET += BET[i] * SHN[i];
222  };
223  TBET += BET[IPS] * APR;
224
225  if (verboseLevel> 3) {
226    G4cout << " TBET " << TBET << G4endl;
227  }
228
229  G4double TBET1 = 0.0;
230
231  if (PPHN > 0.0 || PPHZ > 0.0) 
232    TBET1 = 0.5 * ((TNUN[INS] + TKSN[INS] * ANE) * PPHN + 
233                   (TNUP[IPS] + TKSP[IPS] * APR) * PPHZ);
234  TBET -= TBET1;
235
236  if (verboseLevel > 3) {
237    G4cout << " TBET1 " << TBET1 << G4endl;
238  }
239
240  // deformation
241  G4double TDEF = 0.0;
242  G4double X = Z * Z / A;
243  G4double X1 = G4cbrt(A);
244  G4double X2 = X1 * X1;
245
246  if (IPS != INS && INS >= 3 && IPS >= 2) {
247    G4double X3 = 2.0 * USB - UCB * X;
248    G4double DNZ = 0.0;
249
250    if (!(INS < 3 || (INS - IPS) != 1)) DNZ = C[INS - 3];
251    DNZ = TBET1 - DNZ;
252    G4double X4 = 0.2 * X2 * AL0 * AL0 * X3;
253 
254    if (DNZ > X4) {
255      G4double X5 = USB + X * UCB;
256      G4double X6 = std::log(DNZ / X4);
257      G4double X7 = std::sqrt(X6);
258
259      // G4double ALM = AL0 * (X7 + 0.143 * AL0 * X5 / X3);
260      TDEF = -X4 * (X6 + 1.0) + DNZ + 0.038 * X2 * (AL0*X7 * AL0*X7 * AL0*X7) * X5;
261    };
262  };
263
264  if (verboseLevel > 3) {
265    G4cout << " TDEF " << TDEF << G4endl;
266  }
267
268  // pairing
269  G4double TPE = 0.0; 
270  G4double TPEN = 0.0;
271  G4double TPEP = 0.0;
272  G4double DN0 = 0.0;
273  G4double DZ0 = 0.0;
274  G4double AV = 2.0 * int(0.5 * AN + 0.1);     
275
276  if (AN > AV) DN0 = 1.0;
277
278  AV = 2.0 * int(0.5 * Z + 0.1); 
279
280  if (Z > AV) DZ0 = 1.0;
281
282  if (DN0 * DZ0 > 0.0) TPE = DMU / A;
283
284  if (DN0 + DZ0 > 0.0) {
285
286    if (DN0 > 0.0) {
287      TPEN = RON[INS] / X1;
288
289      if (INS >= 1) {
290
291        if (AN >= 82.0) {
292          TPEN /= X1;
293
294          if (INS == 3) {
295            G4double PN = 0.0;
296            G4double HN = 0.0;
297
298            if (AN > 90.0) PN = AN - 90.0;
299
300            if (AN < 116.0) HN = 116.0 - AN;
301            TPEN += TAU * PN * HN;
302          };
303          TPEN += SINGM[INS - 1] * PPHN;
304        }; 
305      };
306    };
307
308    if (DZ0 > 0.0) {
309      TPEP = ROP[IPS] / X1;
310
311      if (IPS == 1) TPEP += SIPG * PPHZ;
312
313      if (IPS > 1) TPEP = TPEP / X1;
314
315    }; 
316  };
317  TPE += TPEN + TPEP;
318
319  if (verboseLevel > 3) {
320    G4cout << " TPE " << TPE << " TPEN " << TPEN << " TPEP " << TPEP << G4endl;
321  }
322
323  // collect everything
324  DM = (AKU - US * X2 + Z * (Z - 1.0) / X1 * (OMT - UC) - UT / A *
325        (AN - Z) * (AN - Z) + TET + TET1 + TBET + TDEF + TPE) * 0.001;   
326
327  if (verboseLevel > 3) {
328    G4cout << " kummel " << G4endl; 
329  }
330
331  return DM;
332}
Note: See TracBrowser for help on using the repository browser.