This commit was generated by cvs2svn to compensate for changes in r2,
[u/mrichter/AliRoot.git] / GEANT321 / gphys / gpairg.F
1 *
2 * $Id$
3 *
4 * $Log$
5 * Revision 1.5  1998/02/09 15:59:47  japost
6 *   Fixed a problem on AIX 4 xlf, caused by max(double,float).
7 *
8 * Revision 1.4  1998/02/06 16:46:57  japost
9 * Fix a wrong parenthesis.
10 *
11 * Revision 1.3  1998/02/06 16:22:24  japost
12 *   Protected a square root from a negative argument.
13 *   This root was added there in previous changes, and not deleted from its
14 * old position. In its old position it was protected from being negative, but in
15 * its new position it was not.
16 *
17 *   Deleted the same square root from its old position, as it was redundant.
18 *
19 * Revision 1.2  1996/03/13 12:03:24  ravndal
20 * Tranverse momentum conservation
21 *
22 * Revision 1.1.1.1  1995/10/24 10:21:28  cernlib
23 * Geant
24 *
25 *
26 #include "geant321/pilot.h"
27 *CMZ :  3.21/04 21/02/95  11.53.59  by  S.Giani
28 *-- Author :
29 #if defined(CERNLIB_HPUX)
30 $OPTIMIZE OFF
31 #endif
32       SUBROUTINE GPAIRG
33 C.
34 C.    ******************************************************************
35 C.    *                                                                *
36 C.    *  Simulates e+e- pair production by photons.                    *
37 C.    *                                                                *
38 C.    *  The secondary electron energies are sampled using the         *
39 C.    *  Coulomb corrected BETHE-HEITLER cross-sections.For this the   *
40 C.    *   modified version of the random number techniques of          *
41 C.    *   BUTCHER and MESSEL (NUCL.PHYS,20(1960),15) are employed.     *
42 C.    *                                                                *
43 C.    *  NOTE :                                                        *
44 C.    *  (1) Effects due to the breakdown of the BORN approximation at *
45 C.    *      low energies are ignored.                                 *
46 C.    *  (2) The differential cross-section implicitly takes account   *
47 C.    *      of pair production in both nuclear and atomic electron    *
48 C.    *      fields. However, triplet production is not generated.     *
49 C.    *                                                                *
50 C.    *    ==>Called by : GTGAMA                                       *
51 C.    *       Authors    G.Patrick, L.Urban  *********                 *
52 C.    *                                                                *
53 C.    ******************************************************************
54 C.
55 #include "geant321/gcbank.inc"
56 #include "geant321/gcjloc.inc"
57 #include "geant321/gconsp.inc"
58 #include "geant321/gctrak.inc"
59 #include "geant321/gcking.inc"
60 #include "geant321/gcphys.inc"
61 #include "geant321/gccuts.inc"
62       DIMENSION NTYPEL(2)
63       DIMENSION RNDM(2)
64       LOGICAL ROTATE
65       PARAMETER (ONE=1,ONETHR=ONE/3,EMAS2=2*EMASS)
66 C.
67 C.    ------------------------------------------------------------------
68 C.
69 C             If not enough energy : no pair production
70 C
71       EGAM   = VECT(7)
72       IF (EGAM.LT.EMAS2) GO TO 999
73 C
74       KCASE  = NAMEC(6)
75       IF(IPAIR.NE.1) THEN
76          ISTOP  = 2
77          NGKINE = 0
78          DESTEP = DESTEP + EGAM
79          VECT(7)= 0.
80          GEKIN  = 0.
81          GETOT  = 0.
82          GO TO 999
83       ENDIF
84 C
85 C             For low energy photons approximate the electron energy by
86 C             sampling from a uniform distribution in the interval
87 C             EMASS -> EGAM/2.
88 C
89       IF (EGAM.LE.2.1E - 03)THEN
90          CALL GRNDM(RNDM,1)
91          EEL1   = EMASS + (RNDM(1)*(0.5*EGAM - EMASS))
92          X=EEL1/EGAM
93          GO TO 20
94       ENDIF
95 C
96       Z3=Q(JPROB+2)
97       F=8.*Q(JPROB+3)
98       IF(EGAM.GT.0.05) F=F+8.*Q(JPROB+4)
99       X0=EMASS/EGAM
100       DX=0.5-X0
101       DMIN=544.*X0/Z3
102       DMIN2=DMIN*DMIN
103       IF(DMIN.LE.1.)THEN
104          F10=42.392-7.796*DMIN+1.961*DMIN2-F
105          F20=41.405-5.828*DMIN+0.8945*DMIN2-F
106       ELSE
107          F10=42.24-8.368*LOG(DMIN+0.952)-F
108          F20=F10
109       ENDIF
110 C
111 C             Calculate limit for screening variable,DELTA, to ensure
112 C             that screening rejection functions always remain
113 C             positive.
114 C
115       DMAX=EXP((42.24-F)/8.368)-0.952
116 C
117 C             Differential cross-section factors which form
118 C             the coefficients of the screening functions.
119 C
120       DSIG1=DX*DX*F10/3.
121       DSIG2=0.5*F20
122       BPAR   = DSIG1 / (DSIG1 + DSIG2)
123 C
124 C             Decide which screening rejection function to use and
125 C             sample the electron/photon fractional energy BR.
126 C
127    10 CALL GRNDM(RNDM,2)
128       IF(RNDM(1).LT.BPAR)THEN
129          X=0.5-DX*RNDM(2)**ONETHR
130          IREJ=1
131       ELSE
132          X=X0+DX*RNDM(2)
133          IREJ   = 2
134       ENDIF
135 C
136 C             Calculate DELTA ensuring positivity.
137 C
138       D=0.25*DMIN/(X*(1.-X))
139       IF(D.GE.DMAX) GOTO 10
140       D2=D*D
141 C
142 C             Calculate F1 and F2 functions using approximations.
143 C             F10 and F20 are the F1 and F2 functions calculated for the
144 C             DELTA=DELTA minimum.
145 C
146       IF(D.LE.1.)THEN
147          F1=42.392-7.796*D+1.961*D2-F
148          F2=41.405-5.828*D+0.8945*D2-F
149       ELSE
150          F1=42.24-8.368*LOG(D+0.952)-F
151          F2=F1
152       ENDIF
153       IF(IREJ.NE.2)THEN
154          SCREJ=F1/F10
155       ELSE
156          SCREJ=F2/F20
157       ENDIF
158 C
159 C             Accept or reject on basis of random variate.
160 C
161       CALL GRNDM(RNDM,1)
162       IF(RNDM(1).GT.SCREJ) GOTO 10
163       EEL1=X*EGAM
164 C
165 C             Successful sampling of first electron energy.
166 C
167 C             Select charges randomly.
168 C
169    20 NTYPEL(1) = 2
170       CALL GRNDM(RNDM,2)
171       IF (RNDM(1).GT.0.5) NTYPEL(1) = 3
172       NTYPEL(2) = 5 - NTYPEL(1)
173 C
174 C             Generate electron decay angles with respect to a Z-axis
175 C             defined along the parent photon.
176 C             PHI is generated isotropically and THETA is assigned
177 C             a universal angular distribution
178 C
179       EMASS1 = EMASS
180       THETA  = GBTETH(EEL1, EMASS1, X)*EMASS/EEL1
181       SINTH  = SIN(THETA)
182       COSTH  = COS(THETA)
183       PHI    = TWOPI*RNDM(2)
184       COSPHI = COS(PHI)
185       SINPHI = SIN(PHI)
186 C
187 C             Rotate tracks into GEANT system
188 C
189       CALL GFANG(VECT(4),COSAL,SINAL,COSBT,SINBT,ROTATE)
190 C
191 C            Polar co-ordinates to momentum components.
192 C
193       NGKINE = 0
194       TEL1 = EEL1 - EMASS
195       PEL1 = SQRT(MAX((EEL1+REAL(EMASS))*TEL1,0.))
196       IF(TEL1.GT.CUTELE) THEN
197          NGKINE = NGKINE + 1
198          GKIN(1,NGKINE) = PEL1 * SINTH * COSPHI
199          GKIN(2,NGKINE) = PEL1 * SINTH * SINPHI
200          GKIN(3,NGKINE) = PEL1 * COSTH
201          GKIN(4,NGKINE) = EEL1
202          GKIN(5,NGKINE) = NTYPEL(1)
203          TOFD(NGKINE)=0.
204          GPOS(1,NGKINE) = VECT(1)
205          GPOS(2,NGKINE) = VECT(2)
206          GPOS(3,NGKINE) = VECT(3)
207          IF(ROTATE)
208      +   CALL GDROT(GKIN(1,NGKINE),COSAL,SINAL,COSBT,SINBT)
209       ELSE
210          DESTEP = DESTEP + TEL1
211          IF(NTYPEL(1).EQ.2) CALL GANNI2
212       ENDIF
213 C
214 C             Momentum vector of second electron. Recoil momentum of
215 C             target nucleus/electron ignored.
216 C
217       EEL2=EGAM-EEL1
218       TEL2=EEL2-EMASS
219       IF(TEL2.GT.CUTELE) THEN
220          PEL2 = SQRT((EEL2+EMASS)*TEL2)
221          NGKINE = NGKINE + 1
222          SINTH=SINTH*PEL1/PEL2
223          COSTH=SQRT(MAX(0.,1.-SINTH**2))
224          GKIN(1,NGKINE)=-PEL2*SINTH*COSPHI
225          GKIN(2,NGKINE)=-PEL2*SINTH*SINPHI
226          GKIN(3,NGKINE)=PEL2*COSTH
227          GKIN(4,NGKINE)=EEL2
228          GKIN(5,NGKINE) = NTYPEL(2)
229          TOFD(NGKINE)=0.
230          GPOS(1,NGKINE) = VECT(1)
231          GPOS(2,NGKINE) = VECT(2)
232          GPOS(3,NGKINE) = VECT(3)
233          IF(ROTATE)
234      +   CALL GDROT(GKIN(1,NGKINE),COSAL,SINAL,COSBT,SINBT)
235       ELSE
236          DESTEP = DESTEP + TEL2
237          IF(NTYPEL(2).EQ.2) CALL GANNI2
238       ENDIF
239       ISTOP = 1
240       IF(NGKINE.EQ.0) ISTOP = 2
241  999  END
242 #if defined(CERNLIB_HPUX)
243 $OPTIMIZE ON
244 #endif