]>
Commit | Line | Data |
---|---|---|
fe4da5cc | 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 02/07/94 18.28.12 by S.Giani | |
11 | *-- Author : | |
12 | SUBROUTINE GDRAY | |
13 | C. | |
14 | C. ****************************************************************** | |
15 | C. * * | |
16 | C. * Generates Delta rays * | |
17 | C. * * | |
18 | C. * ==>Called by : GTELEC,GTHADR,GTMUON * | |
19 | C. * Authors D.Ward, L.Urban ******** * | |
20 | C. * * | |
21 | C. ****************************************************************** | |
22 | C. | |
23 | #include "geant321/gcphys.inc" | |
24 | #include "geant321/gctrak.inc" | |
25 | #include "geant321/gconsp.inc" | |
26 | #include "geant321/gckine.inc" | |
27 | #include "geant321/gcking.inc" | |
28 | #include "geant321/gccuts.inc" | |
29 | DIMENSION PELS(3) | |
30 | DIMENSION RNDM(2) | |
31 | LOGICAL ROTATE | |
32 | C. | |
33 | C. ------------------------------------------------------------------ | |
34 | C. | |
35 | P=VECT(7) | |
36 | XE=GETOT | |
37 | TE=GEKIN | |
38 | GAM=XE/AMASS | |
39 | GAM2=GAM*GAM | |
40 | T=GAM-1. | |
41 | X=DCUTE/(T*EMASS) | |
42 | C | |
43 | KCASE = NAMEC(10) | |
44 | IF(IPART.EQ.3) THEN | |
45 | C | |
46 | C======> Moller scattering | |
47 | C | |
48 | IF(X.GE.0.5) GO TO 90 | |
49 | CC=1.-2.*X | |
50 | C | |
51 | 10 CALL GRNDM(RNDM,2) | |
52 | E=X/(1.-CC*RNDM(1)) | |
53 | C | |
54 | B1=4./(9.*GAM2-10.*GAM+5.) | |
55 | B2=T*T*B1 | |
56 | B3=(2.*GAM2+2.*GAM-1.)*B1 | |
57 | E1=1.-E | |
58 | C | |
59 | SCREJ=B2*E*E-B3*E/E1+B1*GAM2/(E1*E1) | |
60 | C | |
61 | IF(RNDM(2).GT.SCREJ) GOTO 10 | |
62 | C | |
63 | ELSEIF(IPART.EQ.2)THEN | |
64 | C | |
65 | C======> Bhabha scattering | |
66 | C | |
67 | IF(X.GE.1.) GO TO 90 | |
68 | X1=1.-X | |
69 | 20 CALL GRNDM(RNDM,2) | |
70 | E=X/(1.-X1*RNDM(1)) | |
71 | C | |
72 | Y=1./(GAM+1.) | |
73 | Y2=Y*Y | |
74 | C=1.-2.*Y | |
75 | B1=2.-Y2 | |
76 | B2=C*(3.+Y2) | |
77 | C2=C*C | |
78 | B4=C2*C | |
79 | B3=C2+B4 | |
80 | B0=GAM2/(GAM2-1.) | |
81 | C | |
82 | SCREJ=(((B4*E-B3)*E+B2)*E-B1)*E+B0 | |
83 | SCREJ=SCREJ/((((B4*X-B3)*X+B2)*X-B1)*X+B0) | |
84 | IF(RNDM(2).GT.SCREJ) GOTO 20 | |
85 | C | |
86 | ELSE | |
87 | C | |
88 | C======> Heavy particle. | |
89 | C | |
90 | TMAX=2.*EMASS*(GAM2-1.)/ | |
91 | + (1.+2.*GAM*EMASS/AMASS+(EMASS/AMASS)**2) | |
92 | IF(TMAX.LE.DCUTM) GOTO 90 | |
93 | 40 CALL GRNDM(RNDM,2) | |
94 | E=1./DCUTM+RNDM(1)*(1./TMAX-1./DCUTM) | |
95 | E=1./E | |
96 | BET2=1.-1./GAM2 | |
97 | SCREJ=1.-BET2*(E/TMAX) | |
98 | C --- extra term for spin 1/2 parent. | |
99 | IF(AMASS.GT.0.9 .OR. AMASS.LT.0.12) | |
100 | + SCREJ=SCREJ+0.5*(E/GETOT)**2 | |
101 | IF(RNDM(2).GT.SCREJ) GO TO 40 | |
102 | E=E/(T*EMASS) | |
103 | C | |
104 | ENDIF | |
105 | C | |
106 | EEL=(T*E+1.)*EMASS | |
107 | TEL=EEL-EMASS | |
108 | PEL=SQRT(ABS((EEL+EMASS)*TEL)) | |
109 | COSTH=(XE*EEL+EMASS*(TEL-XE))/(P*PEL) | |
110 | IF(COSTH.GE.1.) THEN | |
111 | COSTH=1. | |
112 | SINTH=0. | |
113 | ELSEIF(COSTH.LE.-1.) THEN | |
114 | COSTH=-1. | |
115 | SINTH=0. | |
116 | ELSE | |
117 | SINTH=SQRT((1.+COSTH)*(1.-COSTH)) | |
118 | ENDIF | |
119 | CALL GRNDM(RNDM,1) | |
120 | PHI = TWOPI*RNDM(1) | |
121 | COSPHI = COS(PHI) | |
122 | SINPHI = SIN(PHI) | |
123 | C | |
124 | C Polar co-ordinates to momentum components. | |
125 | C | |
126 | NGKINE = 1 | |
127 | GKIN(1,1)=PEL*SINTH*COSPHI | |
128 | GKIN(2,1)=PEL*SINTH*SINPHI | |
129 | GKIN(3,1)=PEL*COSTH | |
130 | GKIN(4,1)=EEL | |
131 | GKIN(5,1)=3 | |
132 | TOFD(NGKINE)=0. | |
133 | GPOS(1,NGKINE) = VECT(1) | |
134 | GPOS(2,NGKINE) = VECT(2) | |
135 | GPOS(3,NGKINE) = VECT(3) | |
136 | C | |
137 | PELS(1)=-GKIN(1,1) | |
138 | PELS(2)=-GKIN(2,1) | |
139 | PELS(3)=P-GKIN(3,1) | |
140 | C | |
141 | CALL GFANG(VECT(4),COSTH,SINTH,COSPH,SINPH,ROTATE) | |
142 | IF(ROTATE) THEN | |
143 | CALL GDROT(PELS(1),COSTH,SINTH,COSPH,SINPH) | |
144 | CALL GDROT(GKIN,COSTH,SINTH,COSPH,SINPH) | |
145 | ENDIF | |
146 | C | |
147 | C Correct track vector for lost energy and scattered angles | |
148 | C | |
149 | TELS=TE-TEL | |
150 | EELS=TELS+AMASS | |
151 | PEELS=SQRT(ABS((EELS+AMASS)*TELS)) | |
152 | IF(PEELS.GT.0.)THEN | |
153 | DO 55 I=1,3 | |
154 | VECT(I+3) = PELS(I)/PEELS | |
155 | 55 CONTINUE | |
156 | ENDIF | |
157 | VECT(7) = PEELS | |
158 | GEKIN=TELS | |
159 | GETOT=EELS | |
160 | CALL GEKBIN | |
161 | IF((IDRAY.NE.1).OR.(TEL.LE.CUTELE)) THEN | |
162 | NGKINE = 0 | |
163 | DESTEP = DESTEP + TEL | |
164 | ENDIF | |
165 | C | |
166 | C Update probabilities | |
167 | C | |
168 | 90 CALL GRNDM(RNDM,1) | |
169 | ZINTDR=-LOG(RNDM(1)) | |
170 | SLDRAY=SLENG | |
171 | STEPDR=BIG | |
172 | C | |
173 | END |