Some function moved to AliZDC
[u/mrichter/AliRoot.git] / GEANT321 / gstrag / gkokri.F
1 *
2 * $Id$
3 *
4 * $Log$
5 * Revision 1.1.1.1  1995/10/24 10:21:37  cernlib
6 * Geant
7 *
8 *
9 #include "geant321/pilot.h"
10 *CMZ :  3.21/02 29/03/94  15.41.25  by  S.Giani
11 *-- Author :
12       FUNCTION GKOKRI(E,EMINEV,EMAXEV)
13 C.
14 C.    ******************************************************************
15 C.    *                                                                *
16 C.    *     Input energy in eV, Sandia tables in keV                   *
17 C.    *                                                                *
18 C.    *                                                                *
19 C.    ******************************************************************
20 C.
21 #include "geant321/gcbank.inc"
22 #include "geant321/gcjloc.inc"
23 #include "geant321/gconsp.inc"
24 #include "geant321/gcmate.inc"
25 #include "geant321/gc10ev.inc"
26 #if !defined(CERNLIB_SINGLE)
27       DOUBLE PRECISION EINV,ECUR,ONE,RES,ZERO,EBEG,EEND
28       DOUBLE PRECISION ALPH,BETA,GAMM,E1,E2,GPSCIN,POLE
29       DOUBLE PRECISION EMAX,EMIN,REST
30 #endif
31       PARAMETER (ONE=1,ZERO=0)
32 C.
33 C.    ------------------------------------------------------------------
34 C.
35       REST   = ZERO
36       EMAX   = EMAXEV*1E-3
37       EMIN   = EMINEV*1E-3
38       ECUR   = E*1E-3
39       EINV   = ONE/ECUR
40 C Use Sandia data
41       JPHXS  = LQ(JPHOT-1)
42       NZ     = Q(JPHXS+1)
43       JWEIGH = JPHXS+1+2*NZ
44       DO 30 JZ=1,NZ
45          RES    = ZERO
46          EBEG   = EMIN
47          WEIGHT = Q(JWEIGH+JZ)
48          JPHFN  = LQ(JPHXS-JZ)
49          IPOINT = JPHFN+1
50          IMAX   = Q(IPOINT)
51          IPOINT = IPOINT-4
52          DO 10 I = 1,IMAX
53             IPOINT = IPOINT+5
54             EEND   = Q(IPOINT)
55             IF(EEND.GT.EMIN) THEN
56                E1     = MAX(EBEG,EMIN)
57                E2     = MIN(EEND,EMAX)
58                J      = IPOINT+1
59                IF(ECUR.LE.E2.AND.ECUR.GE.E1) THEN
60 *
61 * *** The pole of the integration is in this interval
62                   EPS1 = (ECUR-EBEG)/ECUR
63                   EPS2 = (EEND-ECUR)/EEND
64                   IF(EPS1.LT.EPS2) THEN
65 *
66 * *** First the pole and then a simple integration
67                      ALPH = ONE-EPS1
68                      E1 = ECUR/ALPH
69                      E2 = EEND
70                   ELSE
71 *
72 * *** First a simple integration and then the pole
73                      ALPH = ONE-EPS2
74                      E1 = EBEG
75                      E2 = ECUR*ALPH
76                   ENDIF
77                   BETA = -LOG(ALPH)
78                   GAMM = ONE/ALPH
79                   POLE = EINV*(Q(J)*BETA+EINV*( Q(J+1)*(GAMM-ALPH)+
80      +            EINV*( Q(J+2)*(0.5*GAMM**2+BETA-0.5*ALPH**2)+EINV*
81      +            Q(J+3)*(GAMM**3/3+GAMM-ALPH-ALPH**3/3))))
82                   RES = RES - POLE
83                ENDIF
84 *
85 * *** This is a normal integration
86                RES = RES + GPSCIN(E1,E2,ECUR,Q(J))
87             ENDIF
88             EBEG = EEND
89             IF(EBEG.GE.EMAX) GO TO 20
90    10    CONTINUE
91    20    REST = REST+WEIGHT*RES
92    30 CONTINUE
93 C RES value is in cm**2/(g keV)
94       GKOKRI = REST*1E-6*A/AVO
95 C Now in Megabarns/eV
96 C
97       END