]>
Commit | Line | Data |
---|---|---|
fe4da5cc | 1 | * |
2 | * $Id$ | |
3 | * | |
4 | * $Log$ | |
5 | * Revision 1.1.1.1 1995/10/24 10:20:01 cernlib | |
6 | * Geant | |
7 | * | |
8 | * | |
9 | #include "geant321/pilot.h" | |
10 | *CMZ : 3.21/02 29/03/94 15.41.44 by S.Giani | |
11 | *-- Author : | |
12 | *$ CREATE THREPD.FOR | |
13 | *COPY THREPD | |
14 | * | |
15 | *=== threpd ===========================================================* | |
16 | * | |
17 | SUBROUTINE THREPD(UMO,ECM1,ECM2,ECM3,PCM1,PCM2,PCM3,COD1,COF1, | |
18 | *SIF1,COD2,COF2,SIF2,COD3,COF3,SIF3,AM1,AM2,AM3) | |
19 | ||
20 | #include "geant321/dblprc.inc" | |
21 | #include "geant321/dimpar.inc" | |
22 | #include "geant321/iounit.inc" | |
23 | * | |
24 | *----------------------------------------------------------------------* | |
25 | * Threpd89: slight revision by A. Ferrari * | |
26 | *----------------------------------------------------------------------* | |
27 | * | |
28 | DIMENSION F(5),XX(5) | |
29 | C***THREE PARTICLE DECAY IN THE CM - SYSTEM | |
30 | COMMON /FKGAMR/ REDU,AMO,AMM(15 ) | |
31 | COMMON/FKDREI/UUMO,AAM1,AAM2,AAM3,S22,UMO2,AM11,AM22,AM33,S2SUP | |
32 | *,S2SAP(2) | |
33 | COMMON/FKPRUN/ISYS | |
34 | REAL RNDM(2) | |
35 | SAVE EPS | |
36 | DATA EPS/1.D-16/ | |
37 | UMOO=UMO+UMO | |
38 | C***S1, S2, S3 ARE THE INVARIANT MASSES OF THE PARTICLES 1, 2, 3 | |
39 | C***J. VON NEUMANN - RANDOM - SELECTION OF S2 | |
40 | C***CALCULATION OF THE MAXIMUM OF THE S2 - DISTRIBUTION | |
41 | UUMO=UMO | |
42 | AAM1=AM1 | |
43 | AAM2=AM2 | |
44 | AAM3=AM3 | |
45 | ||
46 | GU=(AM2+AM3)**2 | |
47 | GO=(UMO-AM1)**2 | |
48 | UFAK=1.0000000000001D0 | |
49 | IF (GU.GT.GO) UFAK=0.9999999999999D0 | |
50 | OFAK=2.D0-UFAK | |
51 | GU=GU*UFAK | |
52 | GO=GO*OFAK | |
53 | DS2=(GO-GU)/99.D0 | |
54 | AM11=AM1*AM1 | |
55 | AM22=AM2*AM2 | |
56 | AM33=AM3*AM3 | |
57 | UMO2=UMO*UMO | |
58 | RHO2=0.D0 | |
59 | S22=GU | |
60 | DO 124 I=1,100 | |
61 | S21=S22 | |
62 | S22=GU+(I-1.D0)*DS2 | |
63 | RHO1=RHO2 | |
64 | RHO2=XLAMB(S22,UMO2,AM11)*XLAMB(S22,AM22,AM33)/(S22+EPS) | |
65 | IF(RHO2.LT.RHO1) GO TO 125 | |
66 | 124 CONTINUE | |
67 | 125 S2SUP=(S22-S21)*.5D0+S21 | |
68 | SUPRHO=XLAMB(S2SUP,UMO2,AM11)*XLAMB(S2SUP,AM22,AM33)/(S2SUP+EPS) | |
69 | SUPRHO=SUPRHO*1.05D0 | |
70 | XO=S21-DS2 | |
71 | IF (GU.LT.GO.AND.XO.LT.GU) XO=GU | |
72 | IF (GU.GT.GO.AND.XO.GT.GU) XO=GU | |
73 | XX(1)=XO | |
74 | XX(3)=S22 | |
75 | X1=(XO+S22)*0.5D0 | |
76 | XX(2)=X1 | |
77 | F(3)=RHO2 | |
78 | F(1)=XLAMB(XO,UMO2,AM11)*XLAMB(XO,AM22,AM33)/(XO+EPS) | |
79 | F(2)=XLAMB(X1,UMO2,AM11)*XLAMB(X1,AM22,AM33)/(X1+EPS) | |
80 | DO 126 I=1,16 | |
81 | X4=(XX(1)+XX(2))*0.5D0 | |
82 | X5=(XX(2)+XX(3))*0.5D0 | |
83 | F(4)=XLAMB(X4,UMO2,AM11)*XLAMB(X4,AM22,AM33)/(X4+EPS) | |
84 | F(5)=XLAMB(X5,UMO2,AM11)*XLAMB(X5,AM22,AM33)/(X5+EPS) | |
85 | XX(4)=X4 | |
86 | XX(5)=X5 | |
87 | DO 128 II=1,5 | |
88 | IA=II | |
89 | DO 128 III=IA,5 | |
90 | IF (F (II).GE.F (III)) GO TO 128 | |
91 | FH=F(II) | |
92 | F(II)=F(III) | |
93 | F(III)=FH | |
94 | FH=XX(II) | |
95 | XX(II)=XX(III) | |
96 | XX(III)=FH | |
97 | 128 CONTINUE | |
98 | SUPRHO=F(1) | |
99 | S2SUP=XX(1) | |
100 | DO 129 II=1,3 | |
101 | IA=II | |
102 | DO 129 III=IA,3 | |
103 | IF (XX(II).GE.XX(III)) GO TO 129 | |
104 | FH=F(II) | |
105 | F(II)=F(III) | |
106 | F(III)=FH | |
107 | FH=XX(II) | |
108 | XX(II)=XX(III) | |
109 | XX(III)=FH | |
110 | 129 CONTINUE | |
111 | 126 CONTINUE | |
112 | AM23=(AM2+AM3)**2 | |
113 | ITH=0 | |
114 | REDU=2.D0 | |
115 | 1 CONTINUE | |
116 | ITH=ITH+1 | |
117 | IF (ITH.GT.200) REDU=-9.D0 | |
118 | IF (ITH.GT.200) GO TO 400 | |
119 | CALL GRNDM(RNDM,2) | |
120 | C=RNDM(1) | |
121 | S2=AM23+C*((UMO-AM1)**2-AM23) | |
122 | Y=RNDM(2) | |
123 | Y=Y*SUPRHO | |
124 | RHO=XLAMB(S2,UMO2,AM11)*XLAMB(S2,AM22,AM33)/S2 | |
125 | IF(Y.GT.RHO) GO TO 1 | |
126 | C***RANDOM SELECTION OF S3 AND CALCULATION OF S1 | |
127 | CALL GRNDM(RNDM,1) | |
128 | S1=RNDM(1) | |
129 | S1=S1*RHO+AM11+AM22-(S2-UMO2+AM11)*(S2+AM22-AM33)/(2.D0*S2)- | |
130 | &RHO*.5D0 | |
131 | S3=UMO2+AM11+AM22+AM33-S1-S2 | |
132 | ECM1=(UMO2+AM11-S2)/UMOO | |
133 | ECM2=(UMO2+AM22-S3)/UMOO | |
134 | ECM3=(UMO2+AM33-S1)/UMOO | |
135 | PCM1=SQRT((ECM1+AM1)*(ECM1-AM1)) | |
136 | PCM2=SQRT((ECM2+AM2)*(ECM2-AM2)) | |
137 | PCM3=SQRT((ECM3+AM3)*(ECM3-AM3)) | |
138 | CALL SFECFE(SFE,CFE) | |
139 | C***TH IS THE ANGLE BETWEEN PARTICLES 1 AND 2 | |
140 | C***TH1, TH2 ARE THE ANGLES BETWEEN PARTICLES 1, 2 AND THE DIRECTION OF | |
141 | IF ((PCM1.LE.1.D-3).OR.(PCM2.LE.1.D-3)) GO TO 200 | |
142 | COSTH=(ECM1*ECM2+0.5D0*(AM11+AM22-S1))/(PCM1*PCM2) | |
143 | GO TO 300 | |
144 | 200 CALL GRNDM(RNDM,1) | |
145 | UW=RNDM(1) | |
146 | COSTH=(UW-.5D0)*2.D0 | |
147 | 300 CONTINUE | |
148 | TMPONE = 0.9999999999999999D0 | |
149 | IF(ABS(COSTH).GT.0.9999999999999999D0) | |
150 | &COSTH=SIGN(TMPONE,COSTH) | |
151 | IF (REDU.LT.1.D0) RETURN | |
152 | COSTH2=(PCM3*PCM3+PCM2*PCM2-PCM1*PCM1)/(2.D0*PCM2*PCM3) | |
153 | IF(ABS(COSTH2).GT.0.9999999999999999D0) | |
154 | &COSTH2=SIGN(TMPONE,COSTH2) | |
155 | SINTH2=SQRT(1.D0-COSTH2*COSTH2) | |
156 | SINTH1=COSTH2*SQRT(1.D0-COSTH*COSTH)-COSTH*SINTH2 | |
157 | COSTH1=COSTH*COSTH2+SINTH2*SQRT(1.D0-COSTH*COSTH) | |
158 | C***RANDOM SELECTION OF THE SPHERICAL COORDINATES OF THE DIRECTION OF PA | |
159 | C***CFE, SFE ARE COS AND SIN OF THE ROTATION ANGLE OF THE SYSTEM 1, 2 AR | |
160 | C***THE DIRECTION OF PARTICLE 3 | |
161 | C***CALCULATION OF THE SPHERICAL COORDINATES OF PARTICLES 1, 2 | |
162 | CX11=-COSTH1 | |
163 | CY11=SINTH1*CFE | |
164 | CZ11=SINTH1*SFE | |
165 | CX22=-COSTH2 | |
166 | CY22=-SINTH2*CFE | |
167 | CZ22=-SINTH2*SFE | |
168 | CALL SFECFE(SIF3,COF3) | |
169 | CALL GRNDM(RNDM,1) | |
170 | COD3=RNDM(1) | |
171 | COD3=2.D0*COD3-1.D0 | |
172 | SID3=SQRT(1.D0-COD3*COD3) | |
173 | 2 FORMAT(5F20.15) | |
174 | COD1=CX11*COD3+CZ11*SID3 | |
175 | IF((1.D0-COD1*COD1).LT.1.D-14)WRITE(ISYS,2)COD1,COF3,SID3, | |
176 | &CX11,CZ11 | |
177 | SID1=SQRT(1.D0-COD1*COD1) | |
178 | COF1=(CX11*SID3*COF3-CY11*SIF3-CZ11*COD3*COF3)/SID1 | |
179 | SIF1=(CX11*SID3*SIF3+CY11*COF3-CZ11*COD3*SIF3)/SID1 | |
180 | COD2=CX22*COD3+CZ22*SID3 | |
181 | SID2=SQRT(1.D0-COD2*COD2) | |
182 | COF2=(CX22*SID3*COF3-CY22*SIF3-CZ22*COD3*COF3)/SID2 | |
183 | SIF2=(CX22*SID3*SIF3+CY22*COF3-CZ22*COD3*SIF3)/SID2 | |
184 | 400 RETURN | |
185 | END |