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

*
* $Id$
*
* $Log$
* Revision 1.5  1998/02/09 15:59:47  japost
*   Fixed a problem on AIX 4 xlf, caused by max(double,float).
*
* Revision 1.4  1998/02/06 16:46:57  japost
* Fix a wrong parenthesis.
*
* Revision 1.3  1998/02/06 16:22:24  japost
*   Protected a square root from a negative argument.
*   This root was added there in previous changes, and not deleted from its
* old position. In its old position it was protected from being negative, but in
* its new position it was not.
*
*   Deleted the same square root from its old position, as it was redundant.
*
* Revision 1.2  1996/03/13 12:03:24  ravndal
* Tranverse momentum conservation
*
* Revision 1.1.1.1  1995/10/24 10:21:28  cernlib
* Geant
*
*
#include "geant321/pilot.h"
*CMZ :  3.21/04 21/02/95  11.53.59  by  S.Giani
*-- Author :
#if defined(CERNLIB_HPUX)
$OPTIMIZE OFF
#endif
      SUBROUTINE GPAIRG
C.
C.    ******************************************************************
C.    *                                                                *
C.    *  Simulates e+e- pair production by photons.                    *
C.    *                                                                *
C.    *  The secondary electron energies are sampled using the         *
C.    *  Coulomb corrected BETHE-HEITLER cross-sections.For this the   *
C.    *   modified version of the random number techniques of          *
C.    *   BUTCHER and MESSEL (NUCL.PHYS,20(1960),15) are employed.     *
C.    *                                                                *
C.    *  NOTE :                                                        *
C.    *  (1) Effects due to the breakdown of the BORN approximation at *
C.    *      low energies are ignored.                                 *
C.    *  (2) The differential cross-section implicitly takes account   *
C.    *      of pair production in both nuclear and atomic electron    *
C.    *      fields. However, triplet production is not generated.     *
C.    *                                                                *
C.    *    ==>Called by : GTGAMA                                       *
C.    *       Authors    G.Patrick, L.Urban  *********                 *
C.    *                                                                *
C.    ******************************************************************
C.
#include "geant321/gcbank.inc"
#include "geant321/gcjloc.inc"
#include "geant321/gconsp.inc"
#include "geant321/gctrak.inc"
#include "geant321/gcking.inc"
#include "geant321/gcphys.inc"
#include "geant321/gccuts.inc"
      DIMENSION NTYPEL(2)
      DIMENSION RNDM(2)
      LOGICAL ROTATE
      PARAMETER (ONE=1,ONETHR=ONE/3,EMAS2=2*EMASS)
C.
C.    ------------------------------------------------------------------
C.
C             If not enough energy : no pair production
C
      EGAM   = VECT(7)
      IF (EGAM.LT.EMAS2) GO TO 999
C
      KCASE  = NAMEC(6)
      IF(IPAIR.NE.1) THEN
         ISTOP  = 2
         NGKINE = 0
         DESTEP = DESTEP + EGAM
         VECT(7)= 0.
         GEKIN  = 0.
         GETOT  = 0.
         GO TO 999
      ENDIF
C
C             For low energy photons approximate the electron energy by
C             sampling from a uniform distribution in the interval
C             EMASS -> EGAM/2.
C
      IF (EGAM.LE.2.1E - 03)THEN
         CALL GRNDM(RNDM,1)
         EEL1   = EMASS + (RNDM(1)*(0.5*EGAM - EMASS))
         X=EEL1/EGAM
         GO TO 20
      ENDIF
C
      Z3=Q(JPROB+2)
      F=8.*Q(JPROB+3)
      IF(EGAM.GT.0.05) F=F+8.*Q(JPROB+4)
      X0=EMASS/EGAM
      DX=0.5-X0
      DMIN=544.*X0/Z3
      DMIN2=DMIN*DMIN
      IF(DMIN.LE.1.)THEN
         F10=42.392-7.796*DMIN+1.961*DMIN2-F
         F20=41.405-5.828*DMIN+0.8945*DMIN2-F
      ELSE
         F10=42.24-8.368*LOG(DMIN+0.952)-F
         F20=F10
      ENDIF
C
C             Calculate limit for screening variable,DELTA, to ensure
C             that screening rejection functions always remain
C             positive.
C
      DMAX=EXP((42.24-F)/8.368)-0.952
C
C             Differential cross-section factors which form
C             the coefficients of the screening functions.
C
      DSIG1=DX*DX*F10/3.
      DSIG2=0.5*F20
      BPAR   = DSIG1 / (DSIG1 + DSIG2)
C
C             Decide which screening rejection function to use and
C             sample the electron/photon fractional energy BR.
C
   10 CALL GRNDM(RNDM,2)
      IF(RNDM(1).LT.BPAR)THEN
         X=0.5-DX*RNDM(2)**ONETHR
         IREJ=1
      ELSE
         X=X0+DX*RNDM(2)
         IREJ   = 2
      ENDIF
C
C             Calculate DELTA ensuring positivity.
C
      D=0.25*DMIN/(X*(1.-X))
      IF(D.GE.DMAX) GOTO 10
      D2=D*D
C
C             Calculate F1 and F2 functions using approximations.
C             F10 and F20 are the F1 and F2 functions calculated for the
C             DELTA=DELTA minimum.
C
      IF(D.LE.1.)THEN
         F1=42.392-7.796*D+1.961*D2-F
         F2=41.405-5.828*D+0.8945*D2-F
      ELSE
         F1=42.24-8.368*LOG(D+0.952)-F
         F2=F1
      ENDIF
      IF(IREJ.NE.2)THEN
         SCREJ=F1/F10
      ELSE
         SCREJ=F2/F20
      ENDIF
C
C             Accept or reject on basis of random variate.
C
      CALL GRNDM(RNDM,1)
      IF(RNDM(1).GT.SCREJ) GOTO 10
      EEL1=X*EGAM
C
C             Successful sampling of first electron energy.
C
C             Select charges randomly.
C
   20 NTYPEL(1) = 2
      CALL GRNDM(RNDM,2)
      IF (RNDM(1).GT.0.5) NTYPEL(1) = 3
      NTYPEL(2) = 5 - NTYPEL(1)
C
C             Generate electron decay angles with respect to a Z-axis
C             defined along the parent photon.
C             PHI is generated isotropically and THETA is assigned
C             a universal angular distribution
C
      EMASS1 = EMASS
      THETA  = GBTETH(EEL1, EMASS1, X)*EMASS/EEL1
      SINTH  = SIN(THETA)
      COSTH  = COS(THETA)
      PHI    = TWOPI*RNDM(2)
      COSPHI = COS(PHI)
      SINPHI = SIN(PHI)
C
C             Rotate tracks into GEANT system
C
      CALL GFANG(VECT(4),COSAL,SINAL,COSBT,SINBT,ROTATE)
C
C            Polar co-ordinates to momentum components.
C
      NGKINE = 0
      TEL1 = EEL1 - EMASS
      PEL1 = SQRT(MAX((EEL1+REAL(EMASS))*TEL1,0.))
      IF(TEL1.GT.CUTELE) THEN
         NGKINE = NGKINE + 1
         GKIN(1,NGKINE) = PEL1 * SINTH * COSPHI
         GKIN(2,NGKINE) = PEL1 * SINTH * SINPHI
         GKIN(3,NGKINE) = PEL1 * COSTH
         GKIN(4,NGKINE) = EEL1
         GKIN(5,NGKINE) = NTYPEL(1)
         TOFD(NGKINE)=0.
         GPOS(1,NGKINE) = VECT(1)
         GPOS(2,NGKINE) = VECT(2)
         GPOS(3,NGKINE) = VECT(3)
         IF(ROTATE)
     +   CALL GDROT(GKIN(1,NGKINE),COSAL,SINAL,COSBT,SINBT)
      ELSE
         DESTEP = DESTEP + TEL1
         IF(NTYPEL(1).EQ.2) CALL GANNI2
      ENDIF
C
C             Momentum vector of second electron. Recoil momentum of
C             target nucleus/electron ignored.
C
      EEL2=EGAM-EEL1
      TEL2=EEL2-EMASS
      IF(TEL2.GT.CUTELE) THEN
         PEL2 = SQRT((EEL2+EMASS)*TEL2)
         NGKINE = NGKINE + 1
         SINTH=SINTH*PEL1/PEL2
         COSTH=SQRT(MAX(0.,1.-SINTH**2))
         GKIN(1,NGKINE)=-PEL2*SINTH*COSPHI
         GKIN(2,NGKINE)=-PEL2*SINTH*SINPHI
         GKIN(3,NGKINE)=PEL2*COSTH
         GKIN(4,NGKINE)=EEL2
         GKIN(5,NGKINE) = NTYPEL(2)
         TOFD(NGKINE)=0.
         GPOS(1,NGKINE) = VECT(1)
         GPOS(2,NGKINE) = VECT(2)
         GPOS(3,NGKINE) = VECT(3)
         IF(ROTATE)
     +   CALL GDROT(GKIN(1,NGKINE),COSAL,SINAL,COSBT,SINBT)
      ELSE
         DESTEP = DESTEP + TEL2
         IF(NTYPEL(2).EQ.2) CALL GANNI2
      ENDIF
      ISTOP = 1
      IF(NGKINE.EQ.0) ISTOP = 2
 999  END
#if defined(CERNLIB_HPUX)
$OPTIMIZE ON
#endif
Commit	Line	Data
fe4da5cc	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.5EGAM - 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.796DMIN+1.961DMIN2-F
105	F20=41.405-5.828DMIN+0.8945DMIN2-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=DXDXF10/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-DXRNDM(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.25DMIN/(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.796D+1.961D2-F
148	F2=41.405-5.828D+0.8945D2-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)=-PEL2SINTHCOSPHI
225	GKIN(2,NGKINE)=-PEL2SINTHSINPHI
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