]>
Commit | Line | Data |
---|---|---|
fe4da5cc | 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 |