Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
pytbdy.f
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file pytbdy.f
1 
2 C*********************************************************************
3 
4 C...PYTBDY
5 C...Generates 3-body decays of gauginos.
6 
7  SUBROUTINE pytbdy(IDIN)
8 
9 C...Double precision and integer declarations.
10  IMPLICIT DOUBLE PRECISION(a-h, o-z)
11  IMPLICIT INTEGER(i-n)
12  INTEGER pyk,pychge,pycomp
13 C...Parameter statement to help give large particle numbers.
14  parameter(ksusy1=1000000,ksusy2=2000000,ktechn=3000000,
15  &kexcit=4000000,kdimen=5000000)
16 C...Commonblocks.
17  common/pyjets/n,npad,k(4000,5),p(4000,5),v(4000,5)
18  common/pydat1/mstu(200),paru(200),mstj(200),parj(200)
19  common/pydat2/kchg(500,4),pmas(500,4),parf(2000),vckm(4,4)
20 C COMMON/PYDAT3/MDCY(500,3),MDME(8000,2),BRAT(8000),KFDP(8000,5)
21 C COMMON/PYPARS/MSTP(200),PARP(200),MSTI(200),PARI(200)
22  common/pyssmt/zmix(4,4),umix(2,2),vmix(2,2),smz(4),smw(2),
23  &sfmix(16,4),zmixi(4,4),umixi(2,2),vmixi(2,2)
24 C SAVE /PYJETS/,/PYDAT1/,/PYDAT2/,/PYDAT3/,/PYPARS/,/PYSSMT/
25  SAVE /pyjets/,/pydat1/,/pydat2/,/pyssmt/
26 
27 C...Local variables.
28  DOUBLE PRECISION xm(5)
29  COMPLEX*16 olpp,orpp,qll,qlr,qrr,qrl,glij,grij,propz
30  COMPLEX*16 qlls,qrrs,qlrs,qrls,qllu,qrru,qlrt,qrlt
31  COMPLEX*16 zmixc(4,4),umixc(2,2),vmixc(2,2)
32  DOUBLE PRECISION s12min,s12max,yjaco1,s23ave,s23df1,s23df2
33  DOUBLE PRECISION d1,d2,d3,p1,p2,p3,cthe1,sthe1,cthe3,sthe3
34  DOUBLE PRECISION cphi1,sphi1
35  DOUBLE PRECISION s23del,eps
36  DOUBLE PRECISION golden,ax,bx,cx,tol,xmin,r,c
37  parameter(r=0.61803399d0,c=1d0-r,tol=1d-3)
38  DOUBLE PRECISION f1,f2,x0,x1,x2,x3
39  INTEGER inoid(4)
40  DATA inoid/22,23,25,35/
41  DATA eps/1d-6/
42 
43  id=idin
44  iskip=1
45  xm(1)=p(n+1,5)
46  xm(2)=p(n+2,5)
47  xm(3)=p(n+3,5)
48  xm(5)=p(id,5)
49 
50 C...GENERATE S12
51  s12min=(xm(1)+xm(2))**2
52  s12max=(xm(5)-xm(3))**2
53  yjaco1=s12max-s12min
54 
55 C...Initialize some parameters
56  xw=paru(102)
57  xw1=1d0-xw
58  tanw=sqrt(xw/xw1)
59  izid1=0
60  iwid1=0
61  izid2=0
62  iwid2=0
63 
64  ia=k(n+2,2)
65  ja=k(n+3,2)
66 
67 C...Mrenna: check that we are indeed decaying a SUSY particle
68  IF(iabs(k(id,2)).LT.ksusy1.OR.iabs(k(id,2)).GE.3000000) THEN
69 
70  ELSE
71  DO 100 i1=1,4
72  IF(mod(k(n+1,2),ksusy1).EQ.inoid(i1)) izid1=i1
73  IF(mod(k(id,2),ksusy1).EQ.inoid(i1)) izid2=i1
74  100 CONTINUE
75  IF(mod(k(n+1,2),ksusy1).EQ.24) iwid1=1
76  IF(mod(k(n+1,2),ksusy1).EQ.37) iwid1=2
77  IF(mod(k(id,2),ksusy1).EQ.24) iwid2=1
78  IF(mod(k(id,2),ksusy1).EQ.37) iwid2=2
79  zm12=xm(5)**2
80  zm22=xm(1)**2
81  ei=kchg(pycomp(iabs(ia)),1)/3d0
82  t3i=sign(1d0,ei+1d-6)/2d0
83  ENDIF
84 
85  IF(max(abs(ia),abs(ja)).EQ.6) THEN
86  iskip=0
87  ELSEIF(izid1*izid2.NE.0) THEN
88  sqmz=pmas(23,1)**2
89  gmmz=pmas(23,1)*pmas(23,2)
90  DO 110 i=1,4
91  zmixc(izid1,i)=dcmplx(zmix(izid1,i),zmixi(izid1,i))
92  zmixc(izid2,i)=dcmplx(zmix(izid2,i),zmixi(izid2,i))
93  110 CONTINUE
94  olpp=(zmixc(izid1,3)*dconjg(zmixc(izid2,3))-
95  & zmixc(izid1,4)*dconjg(zmixc(izid2,4)))/2d0
96  orpp=dconjg(olpp)
97  xll2=pmas(pycomp(ksusy1+iabs(ia)),1)**2
98  xlr2=xll2
99  xrr2=pmas(pycomp(ksusy2+iabs(ia)),1)**2
100  xrl2=xrr2
101  glij=(t3i*zmixc(izid1,2)-tanw*(t3i-ei)*zmixc(izid1,1))*
102  & dconjg(t3i*zmixc(izid2,2)-tanw*(t3i-ei)*zmixc(izid2,1))
103  grij=zmixc(izid1,1)*dconjg(zmixc(izid2,1))*(ei*tanw)**2
104  xm1m2=smz(izid1)*smz(izid2)
105  qlls=dcmplx((t3i-ei*xw)/xw1)*olpp
106  qllu=-glij
107  qlrs=-dcmplx((t3i-ei*xw)/xw1)*orpp
108  qlrt=dconjg(glij)
109  qrls=-dcmplx((ei*xw)/xw1)*olpp
110  qrlt=grij
111  qrrs=dcmplx((ei*xw)/xw1)*orpp
112  qrru=-dconjg(grij)
113  ELSEIF(izid1*iwid2.NE.0.OR.izid2*iwid1.NE.0) THEN
114  IF(izid1.NE.0) THEN
115  xm1m2=smz(izid1)*smw(iwid2)
116  izid1=iwid2
117  izid2=izid1
118  ELSE
119  xm1m2=smz(izid2)*smw(iwid1)
120  izid1=iwid1
121  ENDIF
122  rt2i = 1d0/sqrt(2d0)
123  sqmz=pmas(24,1)**2
124  gmmz=pmas(24,1)*pmas(24,2)
125  DO 120 i=1,2
126  vmixc(izid1,i)=dcmplx(vmix(izid1,i),vmixi(izid1,i))
127  umixc(izid1,i)=dcmplx(umix(izid1,i),umixi(izid1,i))
128  120 CONTINUE
129  DO 130 i=1,4
130  zmixc(izid2,i)=dcmplx(zmix(izid2,i),zmixi(izid2,i))
131  130 CONTINUE
132  qlls=(dconjg(zmixc(izid2,2))*vmixc(izid1,1)-
133  & dconjg(zmixc(izid2,4))*vmixc(izid1,2)*rt2i)
134  qlrs=(zmixc(izid2,2)*dconjg(umixc(izid1,1))+
135  & zmixc(izid2,3)*dconjg(umixc(izid1,2))*rt2i)
136  ej=kchg(iabs(ja),1)/3d0
137  t3j=sign(1d0,ej+1d-6)/2d0
138  qrls=dcmplx(0d0,0d0)
139  qrlt=qrls
140  qrrs=qrls
141  qrru=qrls
142  xrr2=1d6**2
143  xrl2=xrr2
144  xlr2 = pmas(pycomp(ksusy1+iabs(ja)),1)**2
145  xll2 = pmas(pycomp(ksusy1+iabs(ia)),1)**2
146  IF(mod(ia,2).EQ.0) THEN
147  qllu=vmixc(izid1,1)*dconjg(zmixc(izid2,1)*(ei-t3i)*
148  & tanw+zmixc(izid2,2)*t3i)
149  qlrt=-dconjg(umixc(izid1,1))*(
150  & zmixc(izid2,1)*(ej-t3j)*tanw+zmixc(izid2,2)*t3j)
151  ELSE
152  qllu=vmixc(izid1,1)*dconjg(zmixc(izid2,1)*(ej-t3j)*
153  & tanw+zmixc(izid2,2)*t3j)
154  qlrt=-dconjg(umixc(izid1,1))*(
155  & zmixc(izid2,1)*(ei-t3i)*tanw+zmixc(izid2,2)*t3i)
156  ENDIF
157  ELSEIF(iwid1*iwid2.NE.0) THEN
158  izid1=iwid1
159  izid2=iwid2
160  xm1m2=smw(iwid1)*smw(iwid2)
161  sqmz=pmas(23,1)**2
162  gmmz=pmas(23,1)*pmas(23,2)
163  DO 140 i=1,2
164  vmixc(izid1,i)=dcmplx(vmix(izid1,i),vmixi(izid1,i))
165  umixc(izid1,i)=dcmplx(umix(izid1,i),umixi(izid1,i))
166  vmixc(izid2,i)=dcmplx(vmix(izid2,i),vmixi(izid2,i))
167  umixc(izid2,i)=dcmplx(umix(izid2,i),umixi(izid2,i))
168  140 CONTINUE
169  olpp=-vmixc(izid2,1)*dconjg(vmixc(izid1,1))-
170  & vmixc(izid2,2)*dconjg(vmixc(izid1,2))/2d0
171  orpp=-umixc(izid1,1)*dconjg(umixc(izid2,1))-
172  & umixc(izid1,2)*dconjg(umixc(izid2,2))/2d0
173  qrls=-dcmplx(ei/xw1)*orpp
174  qlls=dcmplx((t3i-xw*ei)/xw/xw1)*orpp
175  qrrs=-dcmplx(ei/xw1)*olpp
176  qlrs=dcmplx((t3i-xw*ei)/xw/xw1)*olpp
177  IF(mod(ia,2).EQ.0) THEN
178  xlr2=pmas(pycomp(ksusy1+iabs(ia)-1),1)**2
179  qlrt=-umixc(izid2,1)*dconjg(umixc(izid1,1))*dcmplx(t3i/xw)
180  ELSE
181  xlr2=pmas(pycomp(ksusy1+iabs(ia)+1),1)**2
182  qlrt=-vmixc(izid2,1)*dconjg(vmixc(izid1,1))*dcmplx(t3i/xw)
183  ENDIF
184  ELSEIF(mod(k(n+1,2),ksusy1).EQ.21.OR.mod(k(id,2),ksusy1).EQ.21)
185  &THEN
186  iskip=0
187  ELSE
188  iskip=0
189  ENDIF
190 
191  IF(iskip.NE.0) THEN
192  wtmax=0d0
193  DO 160 kt=1,100
194  s12=s12min+yjaco1*(kt-1)/99
195  s23ave=xm(2)**2+xm(3)**2-(s12+xm(2)**2-xm(1)**2)
196  & *(s12+xm(3)**2-xm(5)**2)/(2d0*s12)
197  s23df1=(s12-xm(2)**2-xm(1)**2)**2
198  & -(2d0*xm(1)*xm(2))**2
199  s23df2=(s12-xm(3)**2-xm(5)**2)**2
200  & -(2d0*xm(3)*xm(5))**2
201  s23df1=s23df1*eps
202  s23df2=s23df2*eps
203  s23del=sqrt(max(0d0,s23df1*s23df2))/(2d0*s12)
204  s23del=s23del/eps
205  s23min=s23ave-s23del
206  s23max=s23ave+s23del
207  yjaco2=s23max-s23min
208  th=s12
209  DO 150 ks=1,100
210  s23=s23min+yjaco2*(ks-1)/99
211  sh=s23
212  uh=zm12+zm22-sh-th
213  wu2 = (uh-zm12)*(uh-zm22)
214  wt2 = (th-zm12)*(th-zm22)
215  ws2 = xm1m2*sh
216  propz2 = (sh-sqmz)**2 + gmmz**2
217  propz=dcmplx(sh-sqmz,-gmmz)/dcmplx(propz2)
218  qll=qlls*propz+qllu/dcmplx(uh-xll2)
219  qlr=qlrs*propz+qlrt/dcmplx(th-xlr2)
220  qrl=qrls*propz+qrlt/dcmplx(th-xrl2)
221  qrr=qrrs*propz+qrru/dcmplx(uh-xrr2)
222  wt0=-((abs(qll)**2+abs(qrr)**2)*wu2+
223  & (abs(qrl)**2+abs(qlr)**2)*wt2+
224  & 2d0*dble(qlr*dconjg(qll)+qrl*dconjg(qrr))*ws2)
225  IF(wt0.GT.wtmax) wtmax=wt0
226  150 CONTINUE
227  160 CONTINUE
228 
229  wtmax=wtmax*1.05d0
230  ENDIF
231 
232 C...FIND S12*
233  ax=s12min
234  cx=s12max
235  bx=s12min+0.5d0*yjaco1
236  x0=ax
237  x3=cx
238  IF(abs(cx-bx).GT.abs(bx-ax))THEN
239  x1=bx
240  x2=bx+c*(cx-bx)
241  ELSE
242  x2=bx
243  x1=bx-c*(bx-ax)
244  ENDIF
245 
246 C...SOLVE FOR F1 AND F2
247  s23df1=(x1-xm(2)**2-xm(1)**2)**2
248  &-(2d0*xm(1)*xm(2))**2
249  s23df2=(x1-xm(3)**2-xm(5)**2)**2
250  &-(2d0*xm(3)*xm(5))**2
251  s23df1=s23df1*eps
252  s23df2=s23df2*eps
253  s23del=sqrt(max(0d0,s23df1*s23df2))/(2d0*x1)
254  f1=-2d0*s23del/eps
255  s23df1=(x2-xm(2)**2-xm(1)**2)**2
256  &-(2d0*xm(1)*xm(2))**2
257  s23df2=(x2-xm(3)**2-xm(5)**2)**2
258  &-(2d0*xm(3)*xm(5))**2
259  s23df1=s23df1*eps
260  s23df2=s23df2*eps
261  s23del=sqrt(max(0d0,s23df1*s23df2))/(2d0*x2)
262  f2=-2d0*s23del/eps
263 
264  170 IF(abs(x3-x0).GT.tol*(abs(x1)+abs(x2)))THEN
265 C...Possibility of infinite loop with .LT.; changed to .LE. (SKANDS)
266  IF(f2.LE.f1)THEN
267  x0=x1
268  x1=x2
269  x2=r*x1+c*x3
270  f1=f2
271  s23df1=(x2-xm(2)**2-xm(1)**2)**2
272  & -(2d0*xm(1)*xm(2))**2
273  s23df2=(x2-xm(3)**2-xm(5)**2)**2
274  & -(2d0*xm(3)*xm(5))**2
275  s23df1=s23df1*eps
276  s23df2=s23df2*eps
277  s23del=sqrt(max(0d0,s23df1*s23df2))/(2d0*x2)
278  f2=-2d0*s23del/eps
279  ELSE
280  x3=x2
281  x2=x1
282  x1=r*x2+c*x0
283  f2=f1
284  s23df1=(x1-xm(2)**2-xm(1)**2)**2
285  & -(2d0*xm(1)*xm(2))**2
286  s23df2=(x1-xm(3)**2-xm(5)**2)**2
287  & -(2d0*xm(3)*xm(5))**2
288  s23df1=s23df1*eps
289  s23df2=s23df2*eps
290  s23del=sqrt(max(0d0,s23df1*s23df2))/(2d0*x1)
291  f1=-2d0*s23del/eps
292  ENDIF
293  goto 170
294  ENDIF
295 C...WE WANT THE MAXIMUM, NOT THE MINIMUM
296  IF(f1.LT.f2)THEN
297  golden=-f1
298  xmin=x1
299  ELSE
300  golden=-f2
301  xmin=x2
302  ENDIF
303 
304  iknt=0
305  180 s12=s12min+pyr(0)*yjaco1
306  iknt=iknt+1
307 C...GENERATE S23
308  s23ave=xm(2)**2+xm(3)**2-(s12+xm(2)**2-xm(1)**2)
309  &*(s12+xm(3)**2-xm(5)**2)/(2d0*s12)
310  s23df1=(s12-xm(2)**2-xm(1)**2)**2
311  &-(2d0*xm(1)*xm(2))**2
312  s23df2=(s12-xm(3)**2-xm(5)**2)**2
313  &-(2d0*xm(3)*xm(5))**2
314  s23df1=s23df1*eps
315  s23df2=s23df2*eps
316  s23del=sqrt(max(0d0,s23df1*s23df2))/(2d0*s12)
317  s23del=s23del/eps
318  s23min=s23ave-s23del
319  s23max=s23ave+s23del
320  yjaco2=s23max-s23min
321  s23=s23min+pyr(0)*yjaco2
322 
323 C...CHECK THE SAMPLING
324  IF(iknt.GT.100) THEN
325  WRITE(mstu(11),*) ' IKNT > 100 IN PYTBDY '
326  goto 190
327  ENDIF
328  IF(yjaco2.LT.pyr(0)*golden) goto 180
329 
330  IF(iskip.EQ.0) goto 190
331 
332  sh=s23
333  th=s12
334  uh=zm12+zm22-sh-th
335 
336  wu2 = (uh-zm12)*(uh-zm22)
337  wt2 = (th-zm12)*(th-zm22)
338  ws2 = xm1m2*sh
339  propz2 = (sh-sqmz)**2 + gmmz**2
340  propz=dcmplx(sh-sqmz,-gmmz)/dcmplx(propz2)
341 
342  qll=qlls*propz+qllu/dcmplx(uh-xll2)
343  qlr=qlrs*propz+qlrt/dcmplx(th-xlr2)
344  qrl=qrls*propz+qrlt/dcmplx(th-xrl2)
345  qrr=qrrs*propz+qrru/dcmplx(uh-xrr2)
346 c QLL=DCMPLX((T3I-EI*XW)/XW1)*OLPP*PROPZ-GLIJ/DCMPLX(UH-XML2)
347 c QLR=-DCMPLX((T3I-EI*XW)/XW1)*ORPP*PROPZ+DCONJG(GLIJ)
348 c &/DCMPLX(TH-XML2)
349 c QRL=-DCMPLX((EI*XW)/XW1)*OLPP*PROPZ+GRIJ/DCMPLX(TH-XMR2)
350 c QRR=DCMPLX((EI*XW)/XW1)*ORPP*PROPZ
351 c &-DCONJG(GRIJ)/DCMPLX(UH-XMR2)
352  wt=-((abs(qll)**2+abs(qrr)**2)*wu2+
353  &(abs(qrl)**2+abs(qlr)**2)*wt2+
354  &2d0*dble(qlr*dconjg(qll)+qrl*dconjg(qrr))*ws2)
355 
356  IF(wt.LT.pyr(0)*wtmax) goto 180
357  IF(wt.GT.wtmax) print*,' WT > WTMAX ',wt,wtmax
358 
359  190 d3=(xm(5)**2+xm(3)**2-s12)/(2d0*xm(5))
360  d1=(xm(5)**2+xm(1)**2-s23)/(2d0*xm(5))
361  d2=xm(5)-d1-d3
362  p1=sqrt(d1*d1-xm(1)**2)
363  p2=sqrt(d2*d2-xm(2)**2)
364  p3=sqrt(d3*d3-xm(3)**2)
365  cthe1=2d0*pyr(0)-1d0
366  ang1=2d0*pyr(0)*paru(1)
367  cphi1=cos(ang1)
368  sphi1=sin(ang1)
369  arg=1d0-cthe1**2
370  IF(arg.LT.0d0.AND.arg.GT.-1d-3) arg=0d0
371  sthe1=sqrt(arg)
372  p(n+1,1)=p1*sthe1*cphi1
373  p(n+1,2)=p1*sthe1*sphi1
374  p(n+1,3)=p1*cthe1
375  p(n+1,4)=d1
376 
377 C...GET CPHI3
378  ang3=2d0*pyr(0)*paru(1)
379  cphi3=cos(ang3)
380  sphi3=sin(ang3)
381  cthe3=(p2**2-p1**2-p3**2)/2d0/p1/p3
382  arg=1d0-cthe3**2
383  IF(arg.LT.0d0.AND.arg.GT.-1d-3) arg=0d0
384  sthe3=sqrt(arg)
385  p(n+3,1)=-p3*sthe3*cphi3*cthe1*cphi1
386  &+p3*sthe3*sphi3*sphi1
387  &+p3*cthe3*sthe1*cphi1
388  p(n+3,2)=-p3*sthe3*cphi3*cthe1*sphi1
389  &-p3*sthe3*sphi3*cphi1
390  &+p3*cthe3*sthe1*sphi1
391  p(n+3,3)=p3*sthe3*cphi3*sthe1
392  &+p3*cthe3*cthe1
393  p(n+3,4)=d3
394 
395  DO 200 i=1,3
396  p(n+2,i)=-p(n+1,i)-p(n+3,i)
397  200 CONTINUE
398  p(n+2,4)=d2
399 
400  RETURN
401  END