Bugfix in AliPoints2Memory
[u/mrichter/AliRoot.git] / GEANT321 / gphys / gcoeff.F
1 *
2 * $Id$
3 *
4 * $Log$
5 * Revision 1.1.1.1  1995/10/24 10:21:23  cernlib
6 * Geant
7 *
8 *
9 #include "geant321/pilot.h"
10 *CMZ :  3.21/02 29/03/94  15.41.21  by  S.Giani
11 *-- Author :
12       SUBROUTINE GCOEFF
13 C.
14 C.    ******************************************************************
15 C.    *                                                                *
16 C.    *  Calculates the coefficients for the energy loss               *
17 C.    *     interpolation                                              *
18 C.    *     There are 4 tables : electron,positron,muon,proton         *
19 C.    *                                                                *
20 C.    *    ==>Called by : GPHYSI                                       *
21 C.    *       Author      F.Carminati *********                        *
22 C.    *                                                                *
23 C.    ******************************************************************
24 C.
25 #include "geant321/gcbank.inc"
26 #include "geant321/gctrak.inc"
27 #include "geant321/gcjloc.inc"
28 #include "geant321/gcmulo.inc"
29 #include "geant321/gconsp.inc"
30 #include "geant321/gcmate.inc"
31 #if !defined(CERNLIB_SINGLE)
32       DOUBLE PRECISION CX1,CX2,CX3,CY1,CY2,CY3,CDEN1,CDEN2,CDEN3
33       DOUBLE PRECISION ACOEFF,BCOEFF,CCOEFF,XRAT,CCOEF1,CCOEF3
34       DOUBLE PRECISION SQEPSM,CFACT
35       PARAMETER (EPSMAC=1E-6)
36 #endif
37 #if defined(CERNLIB_SINGLE)
38       PARAMETER (EPSMAC=1E-11)
39 #endif
40 *
41       SQEPSM = MAX(1.,91./NEK1)*10.*SQRT(EPSMAC)
42       DO 10 IEKBIN=1,NEK1-2
43 *
44          I1 = IEKBIN
45          I2 = I1 + 1
46          I3 = I2 + 1
47          CY1 = ELOW(I1)
48          CY2 = ELOW(I2)
49          CY3 = ELOW(I3)
50          IECOEF = 3*(IEKBIN-1)
51 *
52 * *** Electrons
53 *
54          JRANG = LQ(JMA-15)
55          JCOEF = LQ(JMA-17)
56 *
57          CX1 = Q(JRANG+I1)
58          CX2 = Q(JRANG+I2)
59          CX3 = Q(JRANG+I3)
60          IF(CX1.NE.CX2.AND.CX1.NE.CX3.AND.CX2.NE.CX3) THEN
61             CDEN1 = 1./((CX1-CX2)*(CX1-CX3))
62             CDEN2 = 1./((CX2-CX1)*(CX2-CX3))
63             CDEN3 = 1./((CX3-CX1)*(CX3-CX2))
64             ACOEFF = CY1*CDEN1+CY2*CDEN2+CY3*CDEN3
65             BCOEFF = -(CY1*(CX2+CX3)*CDEN1+CY2*(CX1+CX3)*CDEN2+
66      +                 CY3*(CX1+CX2)*CDEN3)
67             CCOEFF = CY1*CX2*CX3*CDEN1+CX1*CY2*CX3*CDEN2+
68      +               CX1*CX2*CY3*CDEN3
69             IF(ACOEFF.EQ.0.) THEN
70                XRAT = 0.
71             ELSEIF(BCOEFF.GT.0.) THEN
72                CFACT  = SQRT(ABS(ACOEFF))
73                CCOEF1 = SQRT(ABS(CCOEFF-CY1))*CFACT
74                CCOEF3 = SQRT(ABS(CCOEFF-CY3))*CFACT
75                XRAT   = MAX(CCOEF1,CCOEF3)/BCOEFF
76             ELSE
77                XRAT=1.
78             ENDIF
79             IF(XRAT.LE.SQEPSM) THEN
80                Q(JCOEF+IECOEF+1) = 0.
81                Q(JCOEF+IECOEF+2) = BCOEFF
82                Q(JCOEF+IECOEF+3) = CCOEFF
83             ELSE
84                Q(JCOEF+IECOEF+1) = ACOEFF
85                Q(JCOEF+IECOEF+2) = 0.5*BCOEFF/ACOEFF
86                Q(JCOEF+IECOEF+3) = CCOEFF/ACOEFF
87             ENDIF
88          ENDIF
89 *
90 * *** Positons
91 *
92          JRANG = LQ(JMA-15) + NEK1
93          JCOEF = LQ(JMA-17) +3*NEK1
94 *
95          CX1 = Q(JRANG+I1)
96          CX2 = Q(JRANG+I2)
97          CX3 = Q(JRANG+I3)
98          IF(CX1.NE.CX2.AND.CX1.NE.CX3.AND.CX2.NE.CX3) THEN
99             CDEN1 = 1./((CX1-CX2)*(CX1-CX3))
100             CDEN2 = 1./((CX2-CX1)*(CX2-CX3))
101             CDEN3 = 1./((CX3-CX1)*(CX3-CX2))
102             ACOEFF = CY1*CDEN1+CY2*CDEN2+CY3*CDEN3
103             BCOEFF = -(CY1*(CX2+CX3)*CDEN1+CY2*(CX1+CX3)*CDEN2+
104      +                 CY3*(CX1+CX2)*CDEN3)
105             CCOEFF = CY1*CX2*CX3*CDEN1+CX1*CY2*CX3*CDEN2+
106      +               CX1*CX2*CY3*CDEN3
107             IF(ACOEFF.EQ.0.) THEN
108                XRAT = 0.
109             ELSEIF(BCOEFF.GT.0.) THEN
110                CFACT  = SQRT(ABS(ACOEFF))
111                CCOEF1 = SQRT(ABS(CCOEFF-CY1))*CFACT
112                CCOEF3 = SQRT(ABS(CCOEFF-CY3))*CFACT
113                XRAT   = MAX(CCOEF1,CCOEF3)/BCOEFF
114             ELSE
115                XRAT=1.
116             ENDIF
117             IF(XRAT.LE.SQEPSM) THEN
118                Q(JCOEF+IECOEF+1) = 0.
119                Q(JCOEF+IECOEF+2) = BCOEFF
120                Q(JCOEF+IECOEF+3) = CCOEFF
121             ELSE
122                Q(JCOEF+IECOEF+1) = ACOEFF
123                Q(JCOEF+IECOEF+2) = 0.5*BCOEFF/ACOEFF
124                Q(JCOEF+IECOEF+3) = CCOEFF/ACOEFF
125             ENDIF
126          ENDIF
127 *
128 * *** Muons
129 *
130          JRANG = LQ(JMA-16)
131          JCOEF = LQ(JMA-18)
132 *
133          CX1 = Q(JRANG+I1)
134          CX2 = Q(JRANG+I2)
135          CX3 = Q(JRANG+I3)
136          IF(CX1.NE.CX2.AND.CX1.NE.CX3.AND.CX2.NE.CX3) THEN
137             CDEN1 = 1./((CX1-CX2)*(CX1-CX3))
138             CDEN2 = 1./((CX2-CX1)*(CX2-CX3))
139             CDEN3 = 1./((CX3-CX1)*(CX3-CX2))
140             ACOEFF = CY1*CDEN1+CY2*CDEN2+CY3*CDEN3
141             BCOEFF = -(CY1*(CX2+CX3)*CDEN1+CY2*(CX1+CX3)*CDEN2+
142      +                 CY3*(CX1+CX2)*CDEN3)
143             CCOEFF = CY1*CX2*CX3*CDEN1+CX1*CY2*CX3*CDEN2+
144      +               CX1*CX2*CY3*CDEN3
145             IF(ACOEFF.EQ.0.) THEN
146                XRAT = 0.
147             ELSEIF(BCOEFF.GT.0.) THEN
148                CFACT  = SQRT(ABS(ACOEFF))
149                CCOEF1 = SQRT(ABS(CCOEFF-CY1))*CFACT
150                CCOEF3 = SQRT(ABS(CCOEFF-CY3))*CFACT
151                XRAT   = MAX(CCOEF1,CCOEF3)/BCOEFF
152             ELSE
153                XRAT=1.
154             ENDIF
155             IF(XRAT.LE.SQEPSM) THEN
156                Q(JCOEF+IECOEF+1) = 0.
157                Q(JCOEF+IECOEF+2) = BCOEFF
158                Q(JCOEF+IECOEF+3) = CCOEFF
159             ELSE
160                Q(JCOEF+IECOEF+1) = ACOEFF
161                Q(JCOEF+IECOEF+2) = 0.5*BCOEFF/ACOEFF
162                Q(JCOEF+IECOEF+3) = CCOEFF/ACOEFF
163             ENDIF
164          ENDIF
165 *
166 * *** Protons
167 *
168          JRANG = LQ(JMA-16) + NEK1
169          JCOEF = LQ(JMA-18) +3*NEK1
170 *
171          CX1 = Q(JRANG+I1)
172          CX2 = Q(JRANG+I2)
173          CX3 = Q(JRANG+I3)
174          IF(CX1.NE.CX2.AND.CX1.NE.CX3.AND.CX2.NE.CX3) THEN
175             CDEN1 = 1./((CX1-CX2)*(CX1-CX3))
176             CDEN2 = 1./((CX2-CX1)*(CX2-CX3))
177             CDEN3 = 1./((CX3-CX1)*(CX3-CX2))
178             ACOEFF = CY1*CDEN1+CY2*CDEN2+CY3*CDEN3
179             BCOEFF = -(CY1*(CX2+CX3)*CDEN1+CY2*(CX1+CX3)*CDEN2+
180      +                 CY3*(CX1+CX2)*CDEN3)
181             CCOEFF = CY1*CX2*CX3*CDEN1+CX1*CY2*CX3*CDEN2+
182      +               CX1*CX2*CY3*CDEN3
183             IF(ACOEFF.EQ.0.) THEN
184                XRAT = 0.
185             ELSEIF(BCOEFF.GT.0.) THEN
186                CFACT  = SQRT(ABS(ACOEFF))
187                CCOEF1 = SQRT(ABS(CCOEFF-CY1))*CFACT
188                CCOEF3 = SQRT(ABS(CCOEFF-CY3))*CFACT
189                XRAT   = MAX(CCOEF1,CCOEF3)/BCOEFF
190             ELSE
191                XRAT=1.
192             ENDIF
193             IF(XRAT.LE.SQEPSM) THEN
194                Q(JCOEF+IECOEF+1) = 0.
195                Q(JCOEF+IECOEF+2) = BCOEFF
196                Q(JCOEF+IECOEF+3) = CCOEFF
197             ELSE
198                Q(JCOEF+IECOEF+1) = ACOEFF
199                Q(JCOEF+IECOEF+2) = 0.5*BCOEFF/ACOEFF
200                Q(JCOEF+IECOEF+3) = CCOEFF/ACOEFF
201             ENDIF
202          ENDIF
203 *
204    10 CONTINUE
205 *
206       END