This commit was generated by cvs2svn to compensate for changes in r2,
[u/mrichter/AliRoot.git] / GEANT321 / gphys / gphot.F
1 *
2 * $Id$
3 *
4 * $Log$
5 * Revision 1.1.1.1  1995/10/24 10:21:29  cernlib
6 * Geant
7 *
8 *
9 #include "geant321/pilot.h"
10 *CMZ :  3.21/02 29/03/94  15.41.22  by  S.Giani
11 *-- Author :
12       SUBROUTINE GPHOT
13 C.
14 C.    ******************************************************************
15 C.    *                                                                *
16 C.    *  GENERATES PHOTO ELECTRIC MECHANISM                            *
17 C.    *  Corrected version of L. Urban's routine.                      *
18 C.    *  Improvements:                                                 *
19 C.    *    1. Angular distributions of photoelectrons from K-L3 shells *
20 C.    *    2. Generation of shell decay mode                           *
21 C.    *    3. Probability of interactioon with a shell = function      *
22 C.    *       of photon energy                                         *
23 C.    *                                                                *
24 C.    *    ==>CALLED BY : GTGAMA                                       *
25 C.    *       AUTHOR    : J. Chwastowski                               *
26 C.    *                                                                *
27 C.    ******************************************************************
28 C.
29 #include "geant321/gcbank.inc"
30 #include "geant321/gctrak.inc"
31 #include "geant321/gcphys.inc"
32 #include "geant321/gconsp.inc"
33 #include "geant321/gcking.inc"
34 #include "geant321/gccuts.inc"
35 #include "geant321/gcjloc.inc"
36 #include "geant321/gcunit.inc"
37       DIMENSION POT(4),PROB(4),RNA(9)
38       EQUIVALENCE (RNA(1),RN01),(RNA(2),RN02),(RNA(3),RN03)
39       EQUIVALENCE (RNA(4),RN04),(RNA(5),RN05),(RNA(6),RN06)
40       EQUIVALENCE (RNA(7),RN07),(RNA(8),RN08),(RNA(9),RN09)
41       EQUIVALENCE (POT(1),POTK),(POT(2),POTL1)
42       EQUIVALENCE (POT(3),POTL2),(POT(4),POTL3)
43       EQUIVALENCE (PROB(1),PROBK),(PROB(2),PROBL1)
44       EQUIVALENCE (PROB(3),PROBL2),(PROB(4),PROBL3)
45       SAVE ZINOLD,POT,NSHELL
46       DATA ZINOLD / 0.0 /
47 C.
48 C.    ------------------------------------------------------------------
49 C.
50       KCASE = NAMEC(8)
51 C
52 C             STOP ELECTRON ?
53 C
54 C Check if the photoelectric effect was activated. If not deposit
55 C gamma & return
56       IF(IPHOT.NE.1) THEN
57          ISTOP = 2
58          NGKINE= 0
59          DESTEP = DESTEP + VECT(7)
60          VECT(7) = 0.
61          GEKIN = 0.
62          GETOT = 0.
63       ELSE
64          E=VECT(7)
65          CALL GRNDM(RNA,9)
66          JPHXS = LQ(JPHOT-1)
67          NZ = Q(JPHXS+1)
68          IF(NZ.GT.1) THEN
69             QS = 0.0
70             QS2 = GPHSG1(E)*RN01
71             DO 10 I = 1,NZ-1
72                QS1 = GPHSGP(I,E)
73                QS = QS+QS1
74                IF(QS2.LE.QS) THEN
75                   K = I
76                   GO TO 20
77                ENDIF
78    10       CONTINUE
79             K = NZ
80    20       CONTINUE
81             JPHFN = LQ(JPHXS-K)
82             NUSED = Q(JPHFN+1)*5+1
83             JFN = JPHFN+NUSED
84             ZINT = Q(JPHXS+1+K)
85          ELSE
86             JPHFN = LQ(JPHXS-1)
87             NUSED = Q(JPHFN+1)*5+1
88             JFN = JPHFN+NUSED
89             ZINT = Q(JPHXS+1+1)
90          ENDIF
91 C COPY SHELLS POTENTIALS FROM THE ZEBRA STUCTURE
92 C Check if this atom was used in last entry
93          IF(ZINT.NE.ZINOLD) THEN
94             NSHELL = Q(JFN+1)
95             DO 30 I = 1,NSHELL
96                POT(I) = Q(JFN+1+I)
97    30       CONTINUE
98             ZINOLD = ZINT
99          ENDIF
100 C Check if E-gamma is bigger than the L3 ionization potential.
101 C This will make GPHOT a little faster.
102          ISHELL = 0
103          PROB(1) = 0.
104          PROB(2) = 0.
105          PROB(3) = 0.
106          PROB(4) = 0.
107          IF(E.GE.POTL3) THEN
108 C If ZINT < 5 we can have K shell only, so
109             IF(ZINT.LT.5) THEN
110                IF(E.GT.POTK) THEN
111                   PROBK = 1.
112                   TK = E-POTK
113                   ISHELL = 1
114                ENDIF
115             ELSE
116 C The probabilities given below come from crude approximation
117 C It uses the jump ratios and assumes that they are valid for the whole energy
118 C range.
119                IF(E.LT.POTL2) THEN
120                   PROBL3 = 1.0
121                   TK = E-POTL3
122                   ISHELL = 4
123                ELSE
124                   E3 = E-POTL3
125                   GAMAL3 = E3/EMASS+1.
126                   BETAL3 = SQRT(E3*(E3+2.0*EMASS))/(E+EMASS)
127                   E2 = E-POTL2
128                   GAMAL2 = E2/EMASS+1.
129                   BETAL2 = SQRT(E2*(E2+2.0*EMASS))/(E+EMASS)
130                   EFRAC = EMASS/E
131                   PROBL3 = GAVRL3(GAMAL3,BETAL3,EFRAC)
132                   PROBL2 = GAVRL2(GAMAL2,BETAL2,EFRAC)
133                   ANOR = 1./(PROBL3+PROBL2)
134                   PROBL3 = PROBL3*ANOR
135                   PROBL2 = PROBL2*ANOR
136                   IF(E.LT.POTL1) THEN
137                      IF(RN02.LT.PROBL3) THEN
138                         ISHELL = 4
139                         TK = E-POTL3
140                      ELSE
141                         ISHELL = 3
142                         TK = E-POTL2
143                      ENDIF
144                   ELSE
145 C Parametrization of L1 jump ratio gives constant 1.2
146                      PROBL1 = 1.-1./1.2
147                      IF(E.LT.POTK) THEN
148                         PROBL2 = (1.-PROBL1)*PROBL2
149                         PROBL3 = (1.-PROBL1)*PROBL3
150                      ELSE
151                         PROBK = 125./ZINT+3.5
152                         PROBK = 1.-1/PROBK
153                         PROBL1 = (1.-PROBK)*PROBL1
154                         PROBL2 = (1.-PROBK-PROBL1)*PROBL2
155                         PROBL3 = (1.-PROBK-PROBL1)*PROBL3
156                      ENDIF
157                      IF(POTL3.LE.0.0) PROBL3 = 0.0
158                      IF(POTL2.LE.0.0) PROBL2 = 0.0
159                      IF(POTL1.LE.0.0) PROBL1 = 0.0
160                      ANOR = PROBK+PROBL1+PROBL2+PROBL3
161                      IF(ANOR.GT.0.0) THEN
162                         ANOR = 1./ANOR
163                         PROBK = PROBK*ANOR
164                         PROBL1 = PROBL1*ANOR+PROBK
165                         PROBL2 = PROBL2*ANOR+PROBL1
166                         PROBL3 = PROBL3*ANOR+PROBL2
167                         ISHELL = 4
168                         TK = E-POTL3
169                         IF(RN02.LE.PROBK) THEN
170                            ISHELL = 1
171                            TK = E-POTK
172                         ELSEIF(RN02.LE.PROBL1) THEN
173                            ISHELL = 2
174                            TK = E-POTL1
175                         ELSEIF(RN02.LE.PROBL2) THEN
176                            ISHELL = 3
177                            TK = E-POTL2
178                         ENDIF
179                      ENDIF
180                   ENDIF
181                ENDIF
182             ENDIF
183             IF(TK.LE.CUTELE) ISHELL = -ISHELL
184          ENDIF
185          IF(ISHELL.LT.1) THEN
186 C None of the shells was chosen because of the CUTELE
187             ISTOP = 2
188             IF(ISHELL.LT.0) THEN
189                DESTEP = DESTEP+TK
190             ELSEIF(ISHELL.EQ.0) THEN
191                DESTEP = DESTEP+VECT(7)
192             ENDIF
193             NGKINE= 0
194             VECT(7) = 0.
195             GEKIN = 0.
196             GETOT = 0.
197          ELSE
198 C
199 C             ENERGY AND MOMENTUM OF PHOTOELECTRON
200 C
201             EEL=TK + EMASS
202             PEL=SQRT((TK+2.*EMASS)*TK)
203             BETA = PEL/EEL
204             ISTOP = 1
205             NGKINE = 1
206             IF(ISHELL.EQ.1) THEN
207                COST = GPHAK(BETA)
208             ELSEIF(ISHELL.EQ.2) THEN
209                COST = GPHAL1(BETA)
210             ELSEIF(ISHELL.EQ.3) THEN
211                COST = GPHAL2(BETA)
212             ELSEIF(ISHELL.EQ.4) THEN
213                COST = GPHAL3(BETA)
214             ENDIF
215             PHI = TWOPI*RN03
216             COSPHI = COS(PHI)
217             SINPHI = SIN(PHI)
218             SINT = SQRT((1.-COST)*(1.+COST))
219             GKIN(1,NGKINE) = PEL*SINT*COSPHI
220             GKIN(2,NGKINE) = PEL*SINT*SINPHI
221             GKIN(3,NGKINE) = PEL*COST
222             GKIN(4,NGKINE) = EEL
223             GKIN(5,NGKINE) = 3.
224             TOFD(NGKINE) = 0
225             GPOS(1,NGKINE) = VECT(1)
226             GPOS(2,NGKINE) = VECT(2)
227             GPOS(3,NGKINE) = VECT(3)
228 C
229 C             ROTATE ELECTRON AND SCATTERED PHOTON INTO GEANT SYSTEM
230 C
231             CALL GVROT(VECT(4),GKIN)
232          ENDIF
233          IF(ISHELL.NE.0) THEN
234             ISHELL = ABS(ISHELL)
235             IF(ZINT.GE.5.AND.POT(ISHELL).GT.MIN(CUTGAM,CUTELE)) THEN
236 C Generate shell decay mode
237                IF(RN04.LE.Q(JFN+1+NSHELL+ISHELL)) THEN
238                   IF(POT(ISHELL).LE.CUTGAM) THEN
239                      DESTEP = DESTEP+POT(ISHELL)
240                   ELSE
241 C Radiative shell decay
242                      JS = JFN+1+2*NSHELL+ISHELL
243                      JS = JPHFN+Q(JS)
244                      NPOINT = Q(JS)
245                      DO 40 I = 1,NPOINT
246                         IF(RN05.LT.Q(JS+I)) THEN
247                            TSEC = Q(JS+NPOINT+I)
248                            IF(TSEC.GT.CUTGAM) THEN
249                               NGKINE = NGKINE+1
250                               PHI = TWOPI*RN06
251                               COST = 2.*RN07-1.
252                               COSPHI = COS(PHI)
253                               SINPHI = SIN(PHI)
254                               SINT = SQRT((1.-COST)*(1.+COST))
255                               GKIN(1,NGKINE) = TSEC*SINT*COSPHI
256                               GKIN(2,NGKINE) = TSEC*SINT*SINPHI
257                               GKIN(3,NGKINE) = TSEC*COST
258                               GKIN(4,NGKINE) = TSEC
259                               GKIN(5,NGKINE) = 1.
260                               TOFD(NGKINE) = 0.
261                               GPOS(1,NGKINE) = VECT(1)
262                               GPOS(2,NGKINE) = VECT(2)
263                               GPOS(3,NGKINE) = VECT(3)
264                            ELSE
265                               DESTEP = DESTEP+ABS(TSEC)
266                            ENDIF
267 C The following particle forces the energy conservation
268                            TSEC = POT(ISHELL)-ABS(TSEC)
269                            IF(TSEC.GT.CUTGAM) THEN
270                               NGKINE = NGKINE+1
271                               PHI = TWOPI*RN08
272                               COST = 2.*RN09-1.
273                               COSPHI = COS(PHI)
274                               SINPHI = SIN(PHI)
275                               SINT = SQRT((1.-COST)*(1.+COST))
276                               GKIN(1,NGKINE) = TSEC*SINT*COSPHI
277                               GKIN(2,NGKINE) = TSEC*SINT*SINPHI
278                               GKIN(3,NGKINE) = TSEC*COST
279                               GKIN(4,NGKINE) = TSEC
280                               GKIN(5,NGKINE) = 1.
281                               TOFD(NGKINE) = 0.
282                               GPOS(1,NGKINE) = VECT(1)
283                               GPOS(2,NGKINE) = VECT(2)
284                               GPOS(3,NGKINE) = VECT(3)
285                            ELSE
286                               DESTEP = DESTEP+TSEC
287                            ENDIF
288                            GO TO 50
289                         ENDIF
290    40                CONTINUE
291    50                CONTINUE
292                   ENDIF
293                ELSE
294                   IF(POT(ISHELL).LE.CUTELE) THEN
295                      DESTEP = DESTEP+POT(ISHELL)
296                   ELSE
297 c Nonradiative decay
298                      JS = JFN+1+3*NSHELL+ISHELL
299                      JS = JPHFN+Q(JS)
300                      NPOINT = Q(JS)
301                      DO 60 I = 1,NPOINT
302                         IF(RN05.LT.Q(JS+I)) THEN
303                            TSEC = Q(JS+NPOINT+I)
304                            IF(TSEC.GT.CUTELE) THEN
305                               EEL=TSEC + EMASS
306                               PEL=SQRT((TSEC+2.*EMASS)*TSEC)
307                               NGKINE = NGKINE+1
308                               PHI = TWOPI*RN06
309                               COST = 2.*RN07-1.
310                               COSPHI = COS(PHI)
311                               SINPHI = SIN(PHI)
312                               SINT = SQRT((1.-COST)*(1.+COST))
313                               GKIN(1,NGKINE) = PEL*SINT*COSPHI
314                               GKIN(2,NGKINE) = PEL*SINT*SINPHI
315                               GKIN(3,NGKINE) = PEL*COST
316                               GKIN(4,NGKINE) = EEL
317                               GKIN(5,NGKINE) = 3.
318                               TOFD(NGKINE) = 0.
319                               GPOS(1,NGKINE) = VECT(1)
320                               GPOS(2,NGKINE) = VECT(2)
321                               GPOS(3,NGKINE) = VECT(3)
322                            ELSE
323                               DESTEP = DESTEP+ABS(TSEC)
324                            ENDIF
325 C The following particle forces the energy conservation
326                            TSEC = POT(ISHELL)-ABS(TSEC)
327                            IF(TSEC.GT.CUTGAM) THEN
328                               NGKINE = NGKINE+1
329                               PHI = TWOPI*RN08
330                               COST = 2.*RN09-1.
331                               COSPHI = COS(PHI)
332                               SINPHI = SIN(PHI)
333                               SINT = SQRT((1.-COST)*(1.+COST))
334                               GKIN(1,NGKINE) = TSEC*SINT*COSPHI
335                               GKIN(2,NGKINE) = TSEC*SINT*SINPHI
336                               GKIN(3,NGKINE) = TSEC*COST
337                               GKIN(4,NGKINE) = TSEC
338                               GKIN(5,NGKINE) = 1.
339                               TOFD(NGKINE) = 0.
340                               GPOS(1,NGKINE) = VECT(1)
341                               GPOS(2,NGKINE) = VECT(2)
342                               GPOS(3,NGKINE) = VECT(3)
343                            ELSE
344                               DESTEP = DESTEP+TSEC
345                            ENDIF
346                            GO TO 70
347                         ENDIF
348    60                CONTINUE
349    70                CONTINUE
350                   ENDIF
351                ENDIF
352             ELSE
353                DESTEP = DESTEP+POT(ISHELL)
354             ENDIF
355          ENDIF
356       ENDIF
357       END