This commit was generated by cvs2svn to compensate for changes in r2,
[u/mrichter/AliRoot.git] / GEANT321 / gphys / gdrelx.F
1 *
2 * $Id$
3 *
4 * $Log$
5 * Revision 1.1.1.1  1995/10/24 10:21:24  cernlib
6 * Geant
7 *
8 *
9 #include "geant321/pilot.h"
10 *CMZ :  3.21/03 06/10/94  18.22.43  by  S.Giani
11 *-- Author :
12       SUBROUTINE GDRELX(A,Z,DENS,T,HMASS,DEDX)
13 C.
14 C.    ******************************************************************
15 C.    *                                                                *
16 C.    *  Calculates the mean 1/DENS*dE/dx of a particle with kinetic   *
17 C.    *  energy T in an element of atomic number Z, atomic weight A    *
18 C.    *  and density DENS ( the density is just used for the           *
19 C.    *  calculation of the density effect in the case of high T).     *
20 C.    *  The routine reproduces the experimental and/or tabulated      *
21 C.    *  energy losses rather well down to T -> 0.                     *
22 C.    *  Simple parametrization is used for  T .le. T2L=2 MeV (see     *
23 C.    *  H.H.Andersen,J.F.Ziegler:Hydrogen stopping powers and         *
24 C.    *  ranges in all element,Pergamon Press,1977.).                  *
25 C.    *  For T .gt. T2L=2 MeV the corrected Bethe-Bloch stopping       *
26 C.    *  power / restricted energy loss formula is used.               *
27 C.    *                                                                *
28 C.    *                                                                *
29 C.    *    ==>Called by : GDRELA                                       *
30 C.    *       Author    L.Urban    *********                           *
31 C.    *                                                                *
32 C.    ******************************************************************
33 C.
34 #include "geant321/gconsp.inc"
35 #include "geant321/gccuts.inc"
36 #include "geant321/gcunit.inc"
37       PARAMETER (AMUKEV=931494.32,D=0.000153537,T1L=0.00001,T2L=0.002)
38       PARAMETER (EMPROT=0.9382723)
39       DIMENSION B(6,92),C(6,92),CECOF(6)
40 *
41       DATA ((B(I,J),I=1,6),J=1,19) /
42      +     1.262,1.44,242.6,12000.,0.1159,18.8,
43      +     1.229,1.397,484.5,5873.,0.05225,41.7,
44      +     1.411,1.6,725.6,3013.,0.04578,47.6,
45      +     2.248,2.59,966.,153.8,0.03475,62.7,
46      +     2.474,2.815,1206.,1060.,0.02855,76.0,
47      +     2.631,2.989,1445.,957.2,0.02819,77.3,
48      +     2.954,3.35,1683.,1900.,0.02513,86.7,
49      +     2.652,3.,1920.,2000.,0.0223,97.7,
50      +     2.085,2.352,2157.,2634.,0.01816,120.,
51      +     1.951,2.199,2393.,2699.,0.01568,139.,
52      +     2.542,2.869,2628.,1854.,0.01472,148.,
53      +     3.792,4.293,2862.,1009.,0.01397,156.,
54      +     4.154,4.739,2766.,164.5,0.02023,162.,
55      +     4.15,4.7,3329.,550.,0.01321,165.,
56      +     3.232,3.647,3561.,1560.,0.01267,172.,
57      +     3.447,3.891,3792.,1219.,0.01211,180.,
58      +     5.047,5.714,4023.,878.6,0.01178,185.,
59      +     5.731,6.5,4253.,530.,0.01123,194.,
60      +     5.151,5.833,4482.,545.7,0.01129,193./
61       DATA ((B(I,J),I=1,6),J=20,38) /
62      +     5.521,6.252,4710.,553.3,0.01112,196.,
63      +     5.201,5.884,4938.,560.9,0.009995,218.,
64      +     4.862,5.496,5165.,568.5,0.009474,230.,
65      +     4.48,5.055,5391.,952.3,0.009117,239.,
66      +     3.983,4.489,5616.,1336.,0.008413,259.,
67      +     3.469,3.907,5725.,1461.,0.008829,270.,
68      +     3.519,3.963,6065.,1243.,0.007782,280.,
69      +     3.14,3.535,6288.,1372.,0.007361,296.,
70      +     3.553,4.004,6205.,555.1,0.008763,310.,
71      +     3.696,4.175,4673.,387.8,0.02188,322.,
72      +     4.21,4.75,6953.,295.2,0.006809,320.,
73      +     5.041,5.697,7173.,202.6,0.006725,324.,
74      +     5.554,6.3,6496.,110.,0.009689,330.,
75      +     5.323,6.012,7611.,292.5,0.006447,338.,
76      +     5.874,6.656,7395.,117.5,0.007684,340.,
77      +     5.611,6.335,8046.,365.2,0.006244,349.,
78      +     6.411,7.25,8262.,220.,0.006087,358.,
79      +     5.694,6.429,8478.,292.9,0.006087,358.,
80      +     6.339,7.159,8693.,330.3,0.006003,363./
81       DATA ((B(I,J),I=1,6),J=39,57)  /
82      +     6.407,7.234,8907.,367.8,0.005889,370.,
83      +     6.734,7.603,9120.,405.2,0.005765,378.,
84      +     6.902,7.791,9333.,442.7,0.005587,390.,
85      +     6.425,7.248,9545.,480.2,0.005367,406.,
86      +     6.799,7.671,9756.,517.6,0.005315,410.,
87      +     6.108,6.887,9966.,555.1,0.005151,423.,
88      +     5.924,6.677,10180.,592.5,0.004919,443.,
89      +     5.238,5.9,10380.,630.,0.004758,458.,
90      +     5.623,6.354,7160.,337.6,0.01394,466.,
91      +     5.814,6.554,10800.,355.5,0.004626,471.,
92      +     6.23,7.024,11010.,370.9,0.00454,480.,
93      +     6.41,7.227,11210.,386.4,0.004474,487.,
94      +     7.5,8.48,8608.,348.,0.009074,494.,
95      +     6.979,7.871,11620.,392.4,0.004402,495.,
96      +     7.725,8.716,11830.,394.8,0.004376,498.,
97      +     8.231,9.289,12030.,397.3,0.004384,497.,
98      +     7.287,8.218,12230.,399.7,0.004447,490.,
99      +     7.899,8.911,12430.,402.1,0.004511,483.,
100      +     8.041,9.071,12630.,404.5,0.00454,480./
101       DATA ((B(I,J),I=1,6),J=58,76)  /
102      +     7.489,8.444,12830.,406.9,0.00442,493.,
103      +     7.291,8.219,13030.,409.3,0.004298,507.,
104      +     7.098,8.,13230.,411.8,0.004182,521.,
105      +     6.91,7.786,13430.,414.2,0.004058,537.,
106      +     6.728,7.58,13620.,416.6,0.003976,548.,
107      +     6.551,7.38,13820.,419.,0.003877,562.,
108      +     6.739,7.592,14020.,421.4,0.003863,564.,
109      +     6.212,6.996,14210.,423.9,0.003725,585.,
110      +     5.517,6.21,14400.,426.3,0.003632,600.,
111      +     5.219,5.874,14600.,428.7,0.003498,623.,
112      +     5.071,5.706,14790.,433.,0.003405,640.,
113      +     4.926,5.542,14980.,433.5,0.003342,652.,
114      +     4.787,5.386,15170.,435.9,0.003292,662.,
115      +     4.893,5.505,15360.,438.4,0.003243,672.,
116      +     5.028,5.657,15550.,440.8,0.003195,682.,
117      +     4.738,5.329,15740.,443.2,0.003186,684.,
118      +     4.574,5.144,15930.,442.4,0.003144,693.,
119      +     5.2,5.851,16120.,441.6,0.003122,698.,
120      +     5.07,5.704,16300.,440.9,0.003082,707./
121       DATA ((B(I,J),I=1,6),J=77,92)  /
122      +     4.945,5.563,16490.,440.1,0.002965,735.,
123      +     4.476,5.034,16670.,439.3,0.002871,759.,
124      +     4.856,5.46,18320.,438.5,0.002542,755.,
125      +     4.308,4.843,17040.,487.8,0.002882,756.,
126      +     4.723,5.311,17220.,537.,0.002913,748.,
127      +     5.319,5.982,17400.,586.3,0.002871,759.,
128      +     5.956,6.7,17800.,677.,0.00266,765.,
129      +     6.158,6.928,17770.,586.3,0.002812,775.,
130      +     6.204,6.979,17950.,586.3,0.002776,785.,
131      +     6.181,6.954,18120.,586.3,0.002748,793.,
132      +     6.949,7.82,18300.,586.3,0.002737,796.,
133      +     7.506,8.448,18480.,586.3,0.002727,799.,
134      +     7.649,8.609,18660.,586.3,0.002697,808.,
135      +     7.71,8.679,18830.,586.3,0.002641,825.,
136      +     7.407,8.336,19010.,586.3,0.002603,837.,
137      +     7.29,8.204,19180.,586.3,0.002573,847./
138       DATA C/92*0.,92*0.,92*0.,92*0.,92*0.,92*0./
139       DATA CECOF/0.42237,0.0304,-0.00038,3.858,-0.1668,0.00158/
140       DATA EPS/0.000001/
141 *     ------------------------------------------------------------------
142 *      in the case of non-integer Z the low energy parameters
143 *      and the ionization potential are taken at INT(Z) !
144 *
145       IZ=INT(Z+EPS)
146       IF((IZ.LT.1).OR.(IZ.GT.92)) GOTO 10
147 *
148 *     Calculate coefficients C(I,J) if it has not been done already
149 *
150       IF(C(1,IZ).LT.EPS) THEN
151          FAC=AVO/A
152          C(1,IZ)=FAC*AMUKEV**0.5*B(1,IZ)
153          C(2,IZ)=FAC*AMUKEV**0.45*B(2,IZ)
154          C(3,IZ)=FAC*B(3,IZ)/AMUKEV
155          C(4,IZ)=B(4,IZ)/AMUKEV
156          C(5,IZ)=AMUKEV*B(5,IZ)
157 *
158 *                     POTI=16.E-9*Z**0.9
159 *
160          POTI=B(6,IZ)*1.E-9
161 *
162 *
163          C(6,IZ)=POTI
164 *
165       ENDIF
166 *     ----------------------------------------------------------------
167       T1LIM=HMASS*T1L/EMPROT
168       T2LIM=HMASS*T2L/EMPROT
169 *
170 *     Calculate dE/dx
171 *     ---> for T .le. T1LIM (very low energy)
172 *
173       IF(T.LE.T1LIM) THEN
174          TAU=T/HMASS
175          DEDX=C(1,IZ)*TAU**0.5
176       ELSE
177 *
178 *     ---> for T1LIM .lt. T   and  T .le. T2LIM (low energy)
179 *
180          IF(T.LE.T2LIM) THEN
181             TAU=T/HMASS
182             SL=C(2,IZ)*TAU**0.45
183             SH=C(3,IZ)*LOG(1.+C(4,IZ)/TAU+C(5,IZ)*TAU)/TAU
184             DEDX=SL*SH/(SL+SH)
185 *
186 *     ---> for T .gt. T2LIM ( "high " energy , Bethe-Bloch formula)
187 *
188          ELSE
189             P=SQRT(T*(T+2.*HMASS))
190             E=T+HMASS
191             BET2=(P/E)**2
192             ETA=P/HMASS
193             ETA2=ETA*ETA
194 *+++ new line follows.....
195             B2G2=ETA*ETA
196 *+++ end of correction
197             TMAX=2.*EMASS*T*(T+2.*HMASS)
198 *+++  correction of the next line
199 *           TMAX=TMAX/(HMASS**2+EMASS**2+EMASS*(T+HMASS))
200             TMAX=TMAX/(HMASS**2+EMASS**2+2.*EMASS*E)
201 *+++ end of correction
202 *
203 *         density correction
204 *
205             POTI=C(6,IZ)
206             CC=1.+2.*LOG(POTI/(28.8E-9*SQRT(DENS*Z/A)))
207 *         condensed material ? ( dens .gt. 0.05 ? )
208             IF(DENS.GT.0.05) THEN
209                IF(POTI.LT.1.E-7) THEN
210                   IF(CC.LT.3.681) THEN
211                      X0=0.2
212                   ELSE
213                      X0=0.326*CC-1.
214                   ENDIF
215                   X1=2.
216                ELSE
217                   IF(CC.LT.5.215) THEN
218                      X0=0.2
219                   ELSE
220                      X0=0.326*CC-1.5
221                   ENDIF
222                   X1=3.
223                ENDIF
224 *         gas ?   ( dens . le . 0.05 ? )
225             ELSE
226                IF(CC.LE.12.25) THEN
227                   IP=INT((CC-10.)/0.5)+1
228                   IF(IP.LT.0) IP=0
229                   IF(IP.GT.4) IP=4
230                   X0=1.6+0.1*FLOAT(IP)
231                   X1=4.
232                ELSE
233                   IF(CC.LE.13.804) THEN
234                      X0=2.
235                      X1=5.
236                   ELSE
237                      X0=0.326*CC-2.5
238                      X1=5.
239                   ENDIF
240                ENDIF
241             ENDIF
242 *
243             XA=CC/4.606
244             XM=3.
245             AA=4.606*(XA-X0)/(X1-X0)**XM
246 *
247             X=LOG10(ETA)
248             DELTA=0.
249             IF(X.GT.X0) THEN
250                DELTA=4.606*X-CC
251                IF(X.LT.X1) DELTA=DELTA+AA*(X1-X)**XM
252             ENDIF
253 *
254 *         calculate shell correction
255 *
256             POTSQ=POTI*POTI
257             IF(ETA.GT.0.13) THEN
258                F1=1./ETA2
259                F2=F1*F1
260                F3=F1*F2
261                F4=(F1*CECOF(1)+F2*CECOF(2)+F3*CECOF(3))*1.E+12
262                F5=(F1*CECOF(4)+F2*CECOF(5)+F3*CECOF(6))*1.E+18
263                CE=F4*POTSQ+F5*POTSQ*POTI
264             ELSE
265                ETA2=0.0169
266                F1=1./ETA2
267                F2=F1*F1
268                F3=F1*F2
269                F4=(F1*CECOF(1)+F2*CECOF(2)+F3*CECOF(3))*1.E+12
270                F5=(F1*CECOF(4)+F2*CECOF(5)+F3*CECOF(6))*1.E+18
271                CE=F4*POTSQ+F5*POTSQ*POTI
272                CE=CE*LOG(T/T2LIM)/LOG(0.0079/T2LIM)
273             ENDIF
274 *
275             F1=D*Z/(A*BET2)
276 *
277 *         stopping power or restricted dE/dx ?
278 *
279 *+++  correction of the next few lines
280 *           IF(DCUTM.GE.TMAX) THEN
281 *              F2=2.*(LOG(TMAX/POTI)-BET2)
282 *           ELSE
283 *              F2=LOG(TMAX*DCUTM/POTSQ)-BET2*(1.+DCUTM/TMAX)
284 *           ENDIF
285             TUPP=DCUTM
286             IF(TMAX.LT.DCUTM) TUPP=TMAX
287             F2=LOG(2.*EMASS*B2G2/POTI)+LOG(TUPP/POTI)
288      +         -BET2*(1.+TUPP/TMAX)
289 *+++ end of correction
290             DEDX=F1*(F2-DELTA-2.*CE/Z)
291 *
292 *
293             TAU=T2LIM/HMASS
294             SL=C(2,IZ)*TAU**0.45
295             SH=C(3,IZ)*LOG(1.+C(4,IZ)/TAU+C(5,IZ)*TAU)/TAU
296             ST=SL*SH/(SL+SH)
297 *
298             TMAX=2.*EMASS*T2LIM*(T2LIM+2.*HMASS)
299 *+++  correction of the next line
300 *           TMAX=TMAX/(HMASS**2+EMASS**2+EMASS*(T2LIM+HMASS))
301             TMAX=TMAX/(HMASS**2+EMASS**2+2.*EMASS*(T2LIM+HMASS))
302 *+++  end of correction
303             BET2=T2LIM*(T2LIM+2.*HMASS)/(T2LIM+HMASS)**2
304             SBB=2.*(LOG(TMAX/POTI)-BET2)
305             SBB=D*Z*SBB/(A*BET2)
306             CORBB=(ST/SBB-1.)*T2LIM
307 *
308             DEDX=DEDX*(1.+CORBB/T)
309 *
310          ENDIF
311       ENDIF
312 *
313       RETURN
314    10 WRITE (CHMAIL,10000) Z
315       CALL GMAIL(1,1)
316 10000 FORMAT(5X,'***GDRELP  (Z.LT.1.).OR.(Z.GT.92.) ==> Z=',F10.1)
317       END