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

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

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