Bugfix in AliPoints2Memory
[u/mrichter/AliRoot.git] / GEANT321 / gphys / gcomp.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 GCOMP
13 C.
14 C.    ******************************************************************
15 C.    *                                                                *
16 C.    *  Simulates photon-electron COMPTON scattering.                 *
17 C.    *                                                                *
18 C.    *  The scattered photon energy is sampled using the quantum-     *
19 C.    *  mechanical KLEIN-NISHINA formula. For this, the random       *
20 C.    *  number techniques of BUTCHER and MESSEL(NUC. PHYS.20(1960),   *
21 C.    *  15) are employed.                                             *
22 C.    *  NOTE :                                                        *
23 C.    *  (1) Effects due to binding of atomic electrons are            *
24 C.    *      ignored(recoil electron energy assumed large compared     *
25 C.    *      with binding energy).                                     *
26 C.    *                                                                *
27 C.    *    ==>Called by : GTGAMA                                       *
28 C.    *       Authors    G.Patrick, L.Urban  *********                 *
29 C.    *                                                                *
30 C.    ******************************************************************
31 C.
32 #include "geant321/gcphys.inc"
33 #include "geant321/gctrak.inc"
34 #include "geant321/gcking.inc"
35 #include "geant321/gconsp.inc"
36 #include "geant321/gccuts.inc"
37       DIMENSION PGAM(3)
38       REAL RNDM(5)
39       LOGICAL ROTATE
40 C.
41 C.    ------------------------------------------------------------------
42 C.
43       KCASE  = NAMEC(7)
44       EGAM1  = VECT(7)
45       EZERO  = EGAM1/EMASS
46       EMINI  = 1.+2.*EZERO
47       EMIN   = 1./EMINI
48       DSIG1  = LOG(EMINI)
49       DSIG2  = 0.5*(1.-EMIN*EMIN)
50       DSIGT  = DSIG1+DSIG2
51 C
52 C             Decide which part of F(E)=1/E+E to sample from.
53 C
54    10 CALL GRNDM(RNDM,3)
55       IF (DSIG1.LT.DSIGT*RNDM(1))THEN
56 C
57 C             Sample from F2(E) distribution.
58 C
59          BRD    = RNDM(2)
60          CALL GRNDM(RNDM(4),1)
61          IF (EZERO.GE.(EZERO+1.)*RNDM(4))THEN
62             CALL GRNDM(RNDM(5),1)
63             BRD    = MAX(BRD,RNDM(5))
64          ENDIF
65 C
66          BR     = EMIN+(1.-EMIN)*BRD
67       ELSE
68          BR     = EMIN*EXP(DSIG1*RNDM(2))
69       ENDIF
70 C
71 C             Scattered photon energy.
72 C
73       EGAM2  = BR*EGAM1
74 C
75 C             Calculate rejection function G(E).
76 C
77       T      = EMASS*(1.-BR)/EGAM2
78       SINTH  = MAX(0.,T*(2.-T))
79       REJ    = 1.0-(BR*SINTH)/(1.+BR*BR)
80       IF (RNDM(3).GT.REJ)                        GO TO 10
81 C
82 C             Successful sampling of scattered photon.
83 C
84 C             CUTS ON ENERGY THRESHOLDS ?
85 C
86       TEL     = EGAM1-EGAM2
87       IF((EGAM2.LE.CUTGAM).AND.(TEL.LE.CUTELE)) THEN
88          ISTOP = 2
89          VECT(7) = 0.
90          GEKIN = 0.
91          GETOT = 0.
92          NGKINE = 0.
93          DESTEP = DESTEP + EGAM1
94          RETURN
95       ENDIF
96 C
97 C             Generate photon angles with respect to a Z-axis
98 C             defined along the parent photon.
99 C             PHI is generated isotropically
100 C
101       SINTH  = SQRT(SINTH)
102       COSTH  = 1.-T
103       CALL GRNDM(RNDM,1)
104       PHI    = TWOPI*RNDM(1)
105       COSPHI = COS(PHI)
106       SINPHI = SIN(PHI)
107 C
108 C             Polar co-ordinates to momentum components.
109 C
110       PGAM(1) = EGAM2*SINTH*COSPHI
111       PGAM(2) = EGAM2*SINTH*SINPHI
112       PGAM(3) = EGAM2*COSTH
113 C
114 C             Momentum vector of recoil electron.
115 C
116       NGKINE = 1
117       EEL    = TEL + EMASS
118       GKIN(1,1) = -PGAM(1)
119       GKIN(2,1) = -PGAM(2)
120       GKIN(3,1) = EGAM1-PGAM(3)
121       GKIN(4,1)=EEL
122       GKIN(5,1)=3
123       TOFD(NGKINE)=0.
124       GPOS(1,1) = VECT(1)
125       GPOS(2,1) = VECT(2)
126       GPOS(3,1) = VECT(3)
127 C
128 C             Rotate electron and scattered photon into GEANT system
129 C
130       CALL GFANG(VECT(4),COSTH,SINTH,COSPH,SINPH,ROTATE)
131       IF(ROTATE) THEN
132          CALL GDROT(PGAM(1),COSTH,SINTH,COSPH,SINPH)
133          CALL GDROT(GKIN,COSTH,SINTH,COSPH,SINPH)
134       ENDIF
135 C
136 C             Correct photon for energy lost and scattered angle
137 C
138       DO 60 I=1,3
139    60 VECT(I+3) = PGAM(I)/EGAM2
140       VECT(7) = EGAM2
141       GETOT = EGAM2
142       GEKIN = EGAM2
143       CALL GEKBIN
144 C
145 C              Stop electron ?
146 C
147       IF((ICOMP.NE.1).OR.(TEL.LE.CUTELE)) THEN
148          NGKINE = 0
149          DESTEP = DESTEP + TEL
150       ENDIF
151 C
152 C             Update probabilities
153 C
154       CALL GRNDM(RNDM,1)
155       ZINTCO=-LOG(RNDM(1))
156 C
157       END