This commit was generated by cvs2svn to compensate for changes in r2,
[u/mrichter/AliRoot.git] / GEANT321 / gphys / gprobi.F
1 *
2 * $Id$
3 *
4 * $Log$
5 * Revision 1.1.1.1  1995/10/24 10:21:33  cernlib
6 * Geant
7 *
8 *
9 #include "geant321/pilot.h"
10 *CMZ :  3.21/02 29/03/94  15.41.23  by  S.Giani
11 *-- Author :
12       SUBROUTINE GPROBI
13 C.
14 C.    ******************************************************************
15 C.    *                                                                *
16 C.    *  Initialise material constants used in the computation of      *
17 C.    *  the probability for various interactions.                     *
18 C.    *                                                                *
19 C.    *    ==>Called by : GPHYSI                                       *
20 C.    *       Authors    R.Brun, G.Patrick, L.Urban  *********         *
21 C.    *                                                                *
22 C.    ******************************************************************
23 C.
24 #include "geant321/gcbank.inc"
25 #include "geant321/gconsp.inc"
26 #include "geant321/gcmate.inc"
27 #include "geant321/gcjloc.inc"
28       DIMENSION EK(4),EL1(4),EL2(4)
29       DATA EK / 0.66644E-8 , 0.22077E-9 ,-0.32552E-11, 0.18199E-13/
30       DATA EL1/-0.29179E-9 , 0.87983E-10,-0.12589E-11, 0.69602E-14/
31       DATA EL2/-0.68606E-9 , 0.10078E-9 ,-0.14496E-11, 0.78809E-14/
32       DATA ALFA/7.29735E-3/
33       DATA REL/0.2817938/
34 C.
35 C.    ------------------------------------------------------------------
36 C.
37       IF(Z.LT.1.) GOTO 999
38       AEFF=A
39       JPROB = LQ(JMA-4)
40       IF(JMIXT.GT.0)THEN
41          JMI1=LQ(JMIXT-1)
42          AEFF=Q(JMI1+1)
43       ENDIF
44 C
45 C             store constants for PAIR/BREMS routines
46 C
47       X      = (Z*ALFA)**2
48       FC     = (( - 0.002 * X + 0.0083) * X - 0.0369) * X + 0.20206
49       FC     = X * (FC + 1. / (1. + X))
50       C1=Z**0.333333
51       C2=LOG(C1)
52       C3=LOG(183./C1)-FC
53       C4=LOG(1440./(C1*C1))/C3
54       Q(JPROB+1)=Z*(Z+C4)*C3/A
55       Q(JPROB+2)=C1
56       Q(JPROB+3)=C2
57       Q(JPROB+4)=FC
58 C
59 C             constants for PHOTOEFFECT
60 C
61       Z2   = Z*Z
62       EKZ  = Z2*(EK(1) +Z*(EK(2) +Z*(EK(3) +Z*EK(4))))
63       EL1Z = Z2*(EL1(1)+Z*(EL1(2)+Z*(EL1(3)+Z*EL1(4))))
64       EL2Z = Z2*(EL2(1)+Z*(EL2(2)+Z*(EL2(3)+Z*EL2(4))))
65       Q(JPROB+5)=EKZ
66       Q(JPROB+6)=EL1Z
67       Q(JPROB+7)=EL2Z
68 C
69 C             Constants for Hadronic interactions
70 C
71       Q(JPROB+8)= 1000.*AEFF/(AVO*DENS)
72 C
73 C             Constants for electron/positron ionisation losses
74 C             and S5 for one-photon annihilation
75 C
76       IF(JMIXT.LE.0)THEN
77          POTI=16.E-9*Z**0.9
78          S1=Z/A
79          S5=Z**5/A*ALFA**4
80       ELSE
81          NLMAT=Q(JMA+11)
82          NLM=IABS(NLMAT)
83          S1=0.
84          S2=0.
85          S5=0.
86          DO 10 J=1,NLM
87             AJ=Q(JMIXT+J)
88             ZJ=Q(JMIXT+NLM+J)
89             WJ=Q(JMIXT+2*NLM+J)
90             S1=S1+WJ*ZJ/AJ
91             S2=S2+WJ*ZJ*LOG(ZJ)/AJ
92             S5=S5+WJ*ZJ**5/AJ*ALFA**4
93    10    CONTINUE
94          POTI=16.E-9*EXP(0.9*S2/S1)
95       ENDIF
96       Q(JPROB+9) = POTI
97       Q(JPROB+10) = LOG(POTI)
98 C
99       CON1=LOG(POTI/EMASS)
100       CON2=DENS*S1
101       CON3=1.+2.*LOG(POTI/(28.8E-9*SQRT(CON2)))
102 C
103 C             Condensed material ?
104 C             (at present that means: DENS.GT.0.05 g/cm**3)
105 C
106       IF(DENS.GT.0.05)THEN
107          IF(POTI.LT.1.E-7)THEN
108             IF(CON3.LT.3.681)THEN
109                CON4=0.2
110             ELSE
111                CON4=0.326*CON3-1.
112             ENDIF
113             CON5=2.
114          ELSE
115             IF(CON3.LT.5.215)THEN
116                CON4=0.2
117             ELSE
118                CON4=0.326*CON3-1.5
119             ENDIF
120             CON5=3.
121          ENDIF
122       ELSE
123 C
124 C             Gas (T=0 C, P= 1 ATM)
125 C             if T.NE. 0 C and/or P.NE. 1 ATM
126 C             you have to modify the variable X
127 C             X=>X+0.5*LOG((273+T C)/(273*P ATM))
128 C             in the function GDRELE
129 C             ------------------------
130 C
131          IF(CON3.LE.12.25)THEN
132             IP=INT((CON3-10.)/0.5)+1
133             IF(IP.LT.0) IP=0
134             IF(IP.GT.4) IP=4
135             CON4=1.6+0.1*FLOAT(IP)
136             CON5=4.
137          ELSE
138             IF(CON3.LE.13.804)THEN
139                CON4=2.
140                CON5=5.
141             ELSE
142                CON4=0.326*CON3-2.5
143                CON5=5.
144             ENDIF
145          ENDIF
146       ENDIF
147 C
148       XA=CON3/4.606
149       CON6=4.606*(XA-CON4)/(CON5-CON4)**3.
150       Q(JPROB+11)=CON1
151       Q(JPROB+12)=CON2
152       Q(JPROB+13)=-CON3
153       Q(JPROB+14)=CON4
154       Q(JPROB+15)=CON5
155       Q(JPROB+16)=CON6
156 C
157 C            constant for delta rays
158 C            (the same constant is used in the Compton
159 C              and Annihilation subroutines )
160 C            and for one-photon annihilation
161 C
162       Q(JPROB+17)=AVO*TWOPI*REL*REL*DENS*S1
163       Q(JPROB+18)=AVO*TWOPI*REL*REL*DENS*S5
164 C
165 C            Constants for Moliere scattering
166 C
167       IF(JMIXT.LE.0)THEN
168          CALL GMOLI(A,Z,1.,1,DENS,Q(JPROB+21),Q(JPROB+25))
169       ELSE
170          CALL GMOLI(Q(JMIXT+1),Q(JMIXT+NLM+1),Q(JMIXT+2*NLM+1),
171      +              NLM,DENS,Q(JPROB+21),Q(JPROB+25))
172       ENDIF
173 C
174 C                Constants for muon bremsstrahlung
175 C
176       Q(JPROB+31)=LOG(189.*EMMU/(EMASS*C1))
177       IF(Z.GT.10)Q(JPROB+31)=Q(JPROB+31)+LOG(0.666666/C1)
178       SE         =SQRT(2.71828)
179       Q(JPROB+32)=189.*SE*EMMU*EMMU/(2.*EMASS*C1)
180       Q(JPROB+33)=0.75*SE*EMMU*C1
181 C
182   999 END