[u/mrichter/AliRoot.git] / GEANT321 / gphys / gphot.F

*
* $Id$
*
* $Log$
* Revision 1.1.1.1  1995/10/24 10:21:29  cernlib
* Geant
*
*
#include "geant321/pilot.h"
*CMZ :  3.21/02 29/03/94  15.41.22  by  S.Giani
*-- Author :
      SUBROUTINE GPHOT
C.
C.    ******************************************************************
C.    *                                                                *
C.    *  GENERATES PHOTO ELECTRIC MECHANISM                            *
C.    *  Corrected version of L. Urban's routine.                      *
C.    *  Improvements:                                                 *
C.    *    1. Angular distributions of photoelectrons from K-L3 shells *
C.    *    2. Generation of shell decay mode                           *
C.    *    3. Probability of interactioon with a shell = function      *
C.    *       of photon energy                                         *
C.    *                                                                *
C.    *    ==>CALLED BY : GTGAMA                                       *
C.    *       AUTHOR    : J. Chwastowski                               *
C.    *                                                                *
C.    ******************************************************************
C.
#include "geant321/gcbank.inc"
#include "geant321/gctrak.inc"
#include "geant321/gcphys.inc"
#include "geant321/gconsp.inc"
#include "geant321/gcking.inc"
#include "geant321/gccuts.inc"
#include "geant321/gcjloc.inc"
#include "geant321/gcunit.inc"
      DIMENSION POT(4),PROB(4),RNA(9)
      EQUIVALENCE (RNA(1),RN01),(RNA(2),RN02),(RNA(3),RN03)
      EQUIVALENCE (RNA(4),RN04),(RNA(5),RN05),(RNA(6),RN06)
      EQUIVALENCE (RNA(7),RN07),(RNA(8),RN08),(RNA(9),RN09)
      EQUIVALENCE (POT(1),POTK),(POT(2),POTL1)
      EQUIVALENCE (POT(3),POTL2),(POT(4),POTL3)
      EQUIVALENCE (PROB(1),PROBK),(PROB(2),PROBL1)
      EQUIVALENCE (PROB(3),PROBL2),(PROB(4),PROBL3)
      SAVE ZINOLD,POT,NSHELL
      DATA ZINOLD / 0.0 /
C.
C.    ------------------------------------------------------------------
C.
      KCASE = NAMEC(8)
C
C             STOP ELECTRON ?
C
C Check if the photoelectric effect was activated. If not deposit
C gamma & return
      IF(IPHOT.NE.1) THEN
         ISTOP = 2
         NGKINE= 0
         DESTEP = DESTEP + VECT(7)
         VECT(7) = 0.
         GEKIN = 0.
         GETOT = 0.
      ELSE
         E=VECT(7)
         CALL GRNDM(RNA,9)
         JPHXS = LQ(JPHOT-1)
         NZ = Q(JPHXS+1)
         IF(NZ.GT.1) THEN
            QS = 0.0
            QS2 = GPHSG1(E)*RN01
            DO 10 I = 1,NZ-1
               QS1 = GPHSGP(I,E)
               QS = QS+QS1
               IF(QS2.LE.QS) THEN
                  K = I
                  GO TO 20
               ENDIF
   10       CONTINUE
            K = NZ
   20       CONTINUE
            JPHFN = LQ(JPHXS-K)
            NUSED = Q(JPHFN+1)*5+1
            JFN = JPHFN+NUSED
            ZINT = Q(JPHXS+1+K)
         ELSE
            JPHFN = LQ(JPHXS-1)
            NUSED = Q(JPHFN+1)*5+1
            JFN = JPHFN+NUSED
            ZINT = Q(JPHXS+1+1)
         ENDIF
C COPY SHELLS POTENTIALS FROM THE ZEBRA STUCTURE
C Check if this atom was used in last entry
         IF(ZINT.NE.ZINOLD) THEN
            NSHELL = Q(JFN+1)
            DO 30 I = 1,NSHELL
               POT(I) = Q(JFN+1+I)
   30       CONTINUE
            ZINOLD = ZINT
         ENDIF
C Check if E-gamma is bigger than the L3 ionization potential.
C This will make GPHOT a little faster.
         ISHELL = 0
         PROB(1) = 0.
         PROB(2) = 0.
         PROB(3) = 0.
         PROB(4) = 0.
         IF(E.GE.POTL3) THEN
C If ZINT < 5 we can have K shell only, so
            IF(ZINT.LT.5) THEN
               IF(E.GT.POTK) THEN
                  PROBK = 1.
                  TK = E-POTK
                  ISHELL = 1
               ENDIF
            ELSE
C The probabilities given below come from crude approximation
C It uses the jump ratios and assumes that they are valid for the whole energy
C range.
               IF(E.LT.POTL2) THEN
                  PROBL3 = 1.0
                  TK = E-POTL3
                  ISHELL = 4
               ELSE
                  E3 = E-POTL3
                  GAMAL3 = E3/EMASS+1.
                  BETAL3 = SQRT(E3*(E3+2.0*EMASS))/(E+EMASS)
                  E2 = E-POTL2
                  GAMAL2 = E2/EMASS+1.
                  BETAL2 = SQRT(E2*(E2+2.0*EMASS))/(E+EMASS)
                  EFRAC = EMASS/E
                  PROBL3 = GAVRL3(GAMAL3,BETAL3,EFRAC)
                  PROBL2 = GAVRL2(GAMAL2,BETAL2,EFRAC)
                  ANOR = 1./(PROBL3+PROBL2)
                  PROBL3 = PROBL3*ANOR
                  PROBL2 = PROBL2*ANOR
                  IF(E.LT.POTL1) THEN
                     IF(RN02.LT.PROBL3) THEN
                        ISHELL = 4
                        TK = E-POTL3
                     ELSE
                        ISHELL = 3
                        TK = E-POTL2
                     ENDIF
                  ELSE
C Parametrization of L1 jump ratio gives constant 1.2
                     PROBL1 = 1.-1./1.2
                     IF(E.LT.POTK) THEN
                        PROBL2 = (1.-PROBL1)*PROBL2
                        PROBL3 = (1.-PROBL1)*PROBL3
                     ELSE
                        PROBK = 125./ZINT+3.5
                        PROBK = 1.-1/PROBK
                        PROBL1 = (1.-PROBK)*PROBL1
                        PROBL2 = (1.-PROBK-PROBL1)*PROBL2
                        PROBL3 = (1.-PROBK-PROBL1)*PROBL3
                     ENDIF
                     IF(POTL3.LE.0.0) PROBL3 = 0.0
                     IF(POTL2.LE.0.0) PROBL2 = 0.0
                     IF(POTL1.LE.0.0) PROBL1 = 0.0
                     ANOR = PROBK+PROBL1+PROBL2+PROBL3
                     IF(ANOR.GT.0.0) THEN
                        ANOR = 1./ANOR
                        PROBK = PROBK*ANOR
                        PROBL1 = PROBL1*ANOR+PROBK
                        PROBL2 = PROBL2*ANOR+PROBL1
                        PROBL3 = PROBL3*ANOR+PROBL2
                        ISHELL = 4
                        TK = E-POTL3
                        IF(RN02.LE.PROBK) THEN
                           ISHELL = 1
                           TK = E-POTK
                        ELSEIF(RN02.LE.PROBL1) THEN
                           ISHELL = 2
                           TK = E-POTL1
                        ELSEIF(RN02.LE.PROBL2) THEN
                           ISHELL = 3
                           TK = E-POTL2
                        ENDIF
                     ENDIF
                  ENDIF
               ENDIF
            ENDIF
            IF(TK.LE.CUTELE) ISHELL = -ISHELL
         ENDIF
         IF(ISHELL.LT.1) THEN
C None of the shells was chosen because of the CUTELE
            ISTOP = 2
            IF(ISHELL.LT.0) THEN
               DESTEP = DESTEP+TK
            ELSEIF(ISHELL.EQ.0) THEN
               DESTEP = DESTEP+VECT(7)
            ENDIF
            NGKINE= 0
            VECT(7) = 0.
            GEKIN = 0.
            GETOT = 0.
         ELSE
C
C             ENERGY AND MOMENTUM OF PHOTOELECTRON
C
            EEL=TK + EMASS
            PEL=SQRT((TK+2.*EMASS)*TK)
            BETA = PEL/EEL
            ISTOP = 1
            NGKINE = 1
            IF(ISHELL.EQ.1) THEN
               COST = GPHAK(BETA)
            ELSEIF(ISHELL.EQ.2) THEN
               COST = GPHAL1(BETA)
            ELSEIF(ISHELL.EQ.3) THEN
               COST = GPHAL2(BETA)
            ELSEIF(ISHELL.EQ.4) THEN
               COST = GPHAL3(BETA)
            ENDIF
            PHI = TWOPI*RN03
            COSPHI = COS(PHI)
            SINPHI = SIN(PHI)
            SINT = SQRT((1.-COST)*(1.+COST))
            GKIN(1,NGKINE) = PEL*SINT*COSPHI
            GKIN(2,NGKINE) = PEL*SINT*SINPHI
            GKIN(3,NGKINE) = PEL*COST
            GKIN(4,NGKINE) = EEL
            GKIN(5,NGKINE) = 3.
            TOFD(NGKINE) = 0
            GPOS(1,NGKINE) = VECT(1)
            GPOS(2,NGKINE) = VECT(2)
            GPOS(3,NGKINE) = VECT(3)
C
C             ROTATE ELECTRON AND SCATTERED PHOTON INTO GEANT SYSTEM
C
            CALL GVROT(VECT(4),GKIN)
         ENDIF
         IF(ISHELL.NE.0) THEN
            ISHELL = ABS(ISHELL)
            IF(ZINT.GE.5.AND.POT(ISHELL).GT.MIN(CUTGAM,CUTELE)) THEN
C Generate shell decay mode
               IF(RN04.LE.Q(JFN+1+NSHELL+ISHELL)) THEN
                  IF(POT(ISHELL).LE.CUTGAM) THEN
                     DESTEP = DESTEP+POT(ISHELL)
                  ELSE
C Radiative shell decay
                     JS = JFN+1+2*NSHELL+ISHELL
                     JS = JPHFN+Q(JS)
                     NPOINT = Q(JS)
                     DO 40 I = 1,NPOINT
                        IF(RN05.LT.Q(JS+I)) THEN
                           TSEC = Q(JS+NPOINT+I)
                           IF(TSEC.GT.CUTGAM) THEN
                              NGKINE = NGKINE+1
                              PHI = TWOPI*RN06
                              COST = 2.*RN07-1.
                              COSPHI = COS(PHI)
                              SINPHI = SIN(PHI)
                              SINT = SQRT((1.-COST)*(1.+COST))
                              GKIN(1,NGKINE) = TSEC*SINT*COSPHI
                              GKIN(2,NGKINE) = TSEC*SINT*SINPHI
                              GKIN(3,NGKINE) = TSEC*COST
                              GKIN(4,NGKINE) = TSEC
                              GKIN(5,NGKINE) = 1.
                              TOFD(NGKINE) = 0.
                              GPOS(1,NGKINE) = VECT(1)
                              GPOS(2,NGKINE) = VECT(2)
                              GPOS(3,NGKINE) = VECT(3)
                           ELSE
                              DESTEP = DESTEP+ABS(TSEC)
                           ENDIF
C The following particle forces the energy conservation
                           TSEC = POT(ISHELL)-ABS(TSEC)
                           IF(TSEC.GT.CUTGAM) THEN
                              NGKINE = NGKINE+1
                              PHI = TWOPI*RN08
                              COST = 2.*RN09-1.
                              COSPHI = COS(PHI)
                              SINPHI = SIN(PHI)
                              SINT = SQRT((1.-COST)*(1.+COST))
                              GKIN(1,NGKINE) = TSEC*SINT*COSPHI
                              GKIN(2,NGKINE) = TSEC*SINT*SINPHI
                              GKIN(3,NGKINE) = TSEC*COST
                              GKIN(4,NGKINE) = TSEC
                              GKIN(5,NGKINE) = 1.
                              TOFD(NGKINE) = 0.
                              GPOS(1,NGKINE) = VECT(1)
                              GPOS(2,NGKINE) = VECT(2)
                              GPOS(3,NGKINE) = VECT(3)
                           ELSE
                              DESTEP = DESTEP+TSEC
                           ENDIF
                           GO TO 50
                        ENDIF
   40                CONTINUE
   50                CONTINUE
                  ENDIF
               ELSE
                  IF(POT(ISHELL).LE.CUTELE) THEN
                     DESTEP = DESTEP+POT(ISHELL)
                  ELSE
c Nonradiative decay
                     JS = JFN+1+3*NSHELL+ISHELL
                     JS = JPHFN+Q(JS)
                     NPOINT = Q(JS)
                     DO 60 I = 1,NPOINT
                        IF(RN05.LT.Q(JS+I)) THEN
                           TSEC = Q(JS+NPOINT+I)
                           IF(TSEC.GT.CUTELE) THEN
                              EEL=TSEC + EMASS
                              PEL=SQRT((TSEC+2.*EMASS)*TSEC)
                              NGKINE = NGKINE+1
                              PHI = TWOPI*RN06
                              COST = 2.*RN07-1.
                              COSPHI = COS(PHI)
                              SINPHI = SIN(PHI)
                              SINT = SQRT((1.-COST)*(1.+COST))
                              GKIN(1,NGKINE) = PEL*SINT*COSPHI
                              GKIN(2,NGKINE) = PEL*SINT*SINPHI
                              GKIN(3,NGKINE) = PEL*COST
                              GKIN(4,NGKINE) = EEL
                              GKIN(5,NGKINE) = 3.
                              TOFD(NGKINE) = 0.
                              GPOS(1,NGKINE) = VECT(1)
                              GPOS(2,NGKINE) = VECT(2)
                              GPOS(3,NGKINE) = VECT(3)
                           ELSE
                              DESTEP = DESTEP+ABS(TSEC)
                           ENDIF
C The following particle forces the energy conservation
                           TSEC = POT(ISHELL)-ABS(TSEC)
                           IF(TSEC.GT.CUTGAM) THEN
                              NGKINE = NGKINE+1
                              PHI = TWOPI*RN08
                              COST = 2.*RN09-1.
                              COSPHI = COS(PHI)
                              SINPHI = SIN(PHI)
                              SINT = SQRT((1.-COST)*(1.+COST))
                              GKIN(1,NGKINE) = TSEC*SINT*COSPHI
                              GKIN(2,NGKINE) = TSEC*SINT*SINPHI
                              GKIN(3,NGKINE) = TSEC*COST
                              GKIN(4,NGKINE) = TSEC
                              GKIN(5,NGKINE) = 1.
                              TOFD(NGKINE) = 0.
                              GPOS(1,NGKINE) = VECT(1)
                              GPOS(2,NGKINE) = VECT(2)
                              GPOS(3,NGKINE) = VECT(3)
                           ELSE
                              DESTEP = DESTEP+TSEC
                           ENDIF
                           GO TO 70
                        ENDIF
   60                CONTINUE
   70                CONTINUE
                  ENDIF
               ENDIF
            ELSE
               DESTEP = DESTEP+POT(ISHELL)
            ENDIF
         ENDIF
      ENDIF
      END
Commit	Line	Data
fe4da5cc	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.0EMASS))/(E+EMASS)
127	E2 = E-POTL2
128	GAMAL2 = E2/EMASS+1.
129	BETAL2 = SQRT(E2(E2+2.0EMASS))/(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) = PELSINTCOSPHI
220	GKIN(2,NGKINE) = PELSINTSINPHI
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) = TSECSINTCOSPHI
256	GKIN(2,NGKINE) = TSECSINTSINPHI
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) = TSECSINTCOSPHI
277	GKIN(2,NGKINE) = TSECSINTSINPHI
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) = PELSINTCOSPHI
314	GKIN(2,NGKINE) = PELSINTSINPHI
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) = TSECSINTCOSPHI
335	GKIN(2,NGKINE) = TSECSINTSINPHI
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