C++InterfacetoTauola
F/jetset-F/tauface-jetset.f
1 /* copyright(c) 1991-2012 free software foundation, inc.
2  this file is part of the gnu c library.
3 
4  the gnu c library is free software; you can redistribute it and/or
5  modify it under the terms of the gnu lesser general Public
6  license as published by the free software foundation; either
7  version 2.1 of the license, or(at your option) any later version.
8 
9  the gnu c library is distributed in the hope that it will be useful,
10  but without any warranty; without even the implied warranty of
11  merchantability or fitness for a particular purpose. see the gnu
12  lesser general Public license for more details.
13 
14  you should have received a copy of the gnu lesser general Public
15  license along with the gnu c library; if not, see
16  <http://www.gnu.org/licenses/>. */
17 
18 
19 /* this header is separate from features.h so that the compiler can
20  include it implicitly at the start of every compilation. it must
21  not itself include <features.h> or any other header that includes
22  <features.h> because the implicit include comes before any feature
23  test macros that may be defined in a source file before it first
24  explicitly includes a system header. gcc knows the name of this
25  header in order to preinclude it. */
26 
27 /* we do support the iec 559 math functionality, real and complex. */
28 
29 /* wchar_t uses iso/iec 10646 (2nd ed., published 2011-03-15) /
30  unicode 6.0. */
31 
32 /* we do not support c11 <threads.h>. */
33 
34  SUBROUTINE tauola(MODE,KEYPOL)
35 c *************************************
36 c general tauola interface, should work in every case until
37 c hepevt is ok, does not check if hepevt is 'clean'
38 c in particular will decay decayed taus...
39 c only longitudinal spin effects are included.
40 c in w decay v-a vertex is assumed
41 c date: 12 dec 1998. date: 21 june 1999. date: 24 jan 2001 date: 24 aug 2001
42 c this is the hepevt class in old style. no d_h_ class pre-name
43  INTEGER nmxhep
44  parameter(nmxhep=10000)
45  REAL*8 phep, vhep ! to be real*4/ *8 depending on host
46  INTEGER nevhep,nhep,isthep,idhep,jmohep,
47  $ jdahep
48  COMMON /hepevt/
49  $ nevhep, ! serial number
50  $ nhep, ! number of particles
51  $ isthep(nmxhep), ! status code
52  $ idhep(nmxhep), ! particle ident KF
53  $ jmohep(2,nmxhep), ! parent particles
54  $ jdahep(2,nmxhep), ! childreen particles
55  $ phep(5,nmxhep), ! four-momentum, mass [GeV]
56  $ vhep(4,nmxhep) ! vertex [mm]
57 * ----------------------------------------------------------------------
58  LOGICAL qedrad
59  COMMON /phoqed/
60  $ qedrad(nmxhep) ! Photos flag
61 * ----------------------------------------------------------------------
62  SAVE hepevt,phoqed
63  COMMON /taupos/ np1, np2
64  REAL*4 phoi(4),phof(4)
65  double precision q1(4),q2(4),p1(4),p2(4),p3(4),p4(4)
66  COMMON / momdec / q1,q2,p1,p2,p3,p4
67 * tauola, photos and jetset overall switches
68  COMMON /libra/ jak1,jak2,itdkrc,ifphot,ifhadm,ifhadp
69  REAL*4 rrr(1)
70  LOGICAL ifpseudo
71  common /pseudocoup/ csc,ssc
72  REAL*4 csc,ssc
73  save pseudocoup
74  COMMON / inout / inut,iout
75 
76 c to switch tau polarization off in taus
77  dimension pol1(4), pol2(4)
78  double precision pol1x(4), pol2x(4)
79  INTEGER ion(3)
80  DATA pol1 /0.0,0.0,0.0,0.0/
81  DATA pol2 /0.0,0.0,0.0,0.0/
82  DATA pi /3.141592653589793238462643d0/
83 
84 c store decay vertexes
85  dimension imother(20)
86  INTEGER kfhiggs(3)
87 
88 c store daughter pointers
89  INTEGER ison
90  COMMON /isons_tau/ison(2)
91  SAVE /isons_tau/
92 
93  IF(mode.EQ.-1) THEN
94 c ***********************
95 
96  jak1 = 0 ! decay mode first tau
97  jak2 = 0 ! decay mode second tau
98  itdkrc=1.0 ! switch of radiative corrections in decay
99  ifphot=1.0 ! PHOTOS switch
100  ifhadm=1.0
101  ifhadp=1.0
102  pol=1.0 ! tau polarization dipswitch must be 1 or 0
103 
104  kfhiggs(1) = 25
105  kfhiggs(2) = 35
106  kfhiggs(3) = 36
107  kfhigch = 37
108  kfz0 = 23
109  kfgam = 22
110  kftau = 15
111  kfnue = 16
112 c couplings of the 'pseudoscalar higgs' as in cern-th/2003-166
113  psi=0.5*pi ! 0.15*PI
114  xmtau=1.777 ! tau mass
115  xmh=120 ! higgs boson mass
116  betah=sqrt(1d0-4*xmtau**2/xmh**2)
117  csc=cos(psi)*betah
118  ssc=sin(psi)
119 c write(*,*) ' scalar component=',csc,' pseudo-scalar component=',ssc
120  IF (ifphot.EQ.1) CALL phoini ! this if PHOTOS was not initialized earlier
121  CALL inietc(jak1,jak2,itdkrc,ifphot)
122  CALL inimas
123  CALL iniphx(0.01d0)
124  CALL initdk
125 c activation of pi0 and eta decays: (1) means on, (0) off
126  ion(1)=0
127  ion(2)=0
128  ion(3)=0
129  CALL taupi0(-1,1,ion)
130  CALL dekay(-1,pol1x)
131  WRITE(iout,7001) pol,psi,ion(1),ion(2),ion(3)
132  ELSEIF(mode.EQ.0) THEN
133 c ***********************
134 c
135 c..... find tau-s and fill common block /taupos/
136 c this is to avoid lund history fillings. this call is optional
137  CALL phyfix(nstop,nstart)
138 c clear mothers of the previous event
139  DO ii=1,20
140  imother(ii)=0
141  ENDDO
142 
143  DO ii=1,2
144  ison(ii)=0
145  ENDDO
146 c ... and to find mothers giving taus.
147  ndec = 0
148 c(bpk)--> look for mother, check that it is not the history entry(e.g. mstp(128)=0)
149  DO i=nstart,nhep
150  IF(abs(idhep(i)).EQ.kftau.AND.isthep(i).EQ.1.AND.
151  $ (isthep(i).GE.125.OR.isthep(i).LT.120)) THEN
152  imoth=jmohep(1,i)
153  DO WHILE (abs(idhep(imoth)).EQ.kftau) ! KEEP WALKING UP
154  imoth=jmohep(1,imoth)
155  ENDDO
156  IF (isthep(imoth).EQ.3.OR.
157  $ (isthep(imoth).GE.120.AND.isthep(imoth).LE.125)) THEN
158  DO j=nstart,nhep ! WE HAVE WALKED INTO HARD RECORD
159  IF (idhep(j).EQ.idhep(imoth).AND.
160  $ jmohep(1,j).EQ.imoth.AND.
161  $ isthep(j).EQ.2) THEN
162  jmoth=j
163  goto 66
164  ENDIF
165  ENDDO
166  ELSE
167  jmoth=imoth
168  ENDIF
169  66 CONTINUE
170  DO ii=1,ndec
171  IF(jmoth.EQ.imother(ii)) goto 9999
172  ENDDO
173 c(bpk)--<
174  ndec=ndec+1
175  imother(ndec)= jmoth
176  ENDIF
177  9999 CONTINUE
178  ENDDO
179 
180 c ... taus of every mother are treated in this main loop
181  DO ii=1,ndec
182  im=imother(ii)
183  ncount=0
184  np1=0
185  np2=0
186 
187 
188 c(bpk)-->
189 c correcting hepevt is out of question at this point..
190  im0=im
191  IF (idhep(jmohep(1,im0)).EQ.idhep(im0)) im0=jmohep(1,im0)
192  isel=-1
193  DO i=nstart,nhep
194  IF (isthep(i).EQ.3.OR.
195  $ (isthep(i).GE.120.AND.isthep(i).LE.125)) THEN ! HARD RECORD
196  goto 76
197  ENDIF
198  imoth=jmohep(1,i)
199  DO WHILE (idhep(imoth).EQ.idhep(i).OR.abs(idhep(imoth)).EQ.kftau) ! KEEP WALKING UP
200  imoth=jmohep(1,imoth)
201  ENDDO
202  IF ((imoth.EQ.im0.OR.imoth.EQ.im).AND.isel.EQ.-1) THEN
203  ison(1)=i
204  isel=0
205  ELSEIF ((imoth.EQ.im0.OR.imoth.EQ.im).AND.isel.EQ.0) THEN
206  ison(2)=i
207  ELSEIF ((imoth.NE.im0.AND.imoth.NE.im).AND.isel.EQ.0) THEN
208  isel=1
209  goto 77
210  ENDIF
211  76 CONTINUE
212  ENDDO
213  77 CONTINUE
214 c(bpk)--<
215 
216 
217 c ... we correct hepevt(fix developped with catherine biscarat)
218 c IF (jdahep(2,im).EQ.0) THEN ! ID of second daughter was missing
219 c isecu=1
220 c DO i=jdahep(1,im)+1,nhep ! OK lets look for it
221 c IF (jmohep(1,i).EQ.im.AND.isecu.EQ.1) THEN ! we have found one
222 c jdahep(2,im)=i
223 c ELSEIF (jmohep(1,i).EQ.im.AND.isecu.NE.1) THEN ! we have found one after there
224 c jdahep(2,im)=0 ! was something else, lets kill game
225 c ENDIF
226 c IF (jmohep(1,i).NE.im) isecu=0 ! other stuff starts
227 c ENDDO
228 c ENDIF
229 
230 c ... we check whether there are just two or more tau-likes
231  DO i=ison(1),ison(2)
232  IF(idhep(i).EQ.-kftau.OR.idhep(i).EQ.-kfnue) ncount=ncount+1
233  IF(idhep(i).EQ. kftau.OR.idhep(i).EQ. kfnue) ncount=ncount+1
234  ENDDO
235 
236 c ... if there will be more we will come here again
237  666 CONTINUE
238 
239 c(bpk)-->
240  DO i=max(np1+1,ison(1)),ison(2)
241 c(bpk)--<
242  IF(idhep(i).EQ.-kftau.OR.idhep(i).EQ.-kfnue) np1=i
243  ENDDO
244 c(bpk)-->
245  DO i=max(np2+1,ison(1)),ison(2)
246 c(bpk)--<
247  IF(idhep(i).EQ. kftau.OR.idhep(i).EQ. kfnue) np2=i
248  ENDDO
249  DO i=1,4
250  p1(i)= phep(i,np1) !momentum of tau+
251  p2(i)= phep(i,np2) !momentum of tau-
252  q1(i)= p1(i)+p2(i)
253  ENDDO
254 
255  pol1(3)= 0d0
256  pol2(3)= 0d0
257 
258  IF(keypol.EQ.1) THEN
259 c.....include polarisation effect
260  CALL ranmar(rrr,1)
261 
262  IF(idhep(im).EQ.kfhiggs(1).OR.idhep(im).EQ.kfhiggs(2).OR.
263  $ idhep(im).EQ.kfhiggs(3)) THEN ! case of Higgs
264  IF(rrr(1).LT.0.5) THEN
265  pol1(3)= pol
266  pol2(3)=-pol
267  ELSE
268  pol1(3)=-pol
269  pol2(3)= pol
270  ENDIF
271  ELSEIF((idhep(im).EQ.kfz0).OR.(idhep(im).EQ.kfgam)) THEN ! case of gamma/Z
272 c there is no angular dependence in gamma/z polarization
273 c there is no s-dependence in gamma/z polarization at all
274 c there is even no z polarization in any form
275 c main reason is that nobody asked ...
276 c but it is prepared and longitudinal correlations
277 c can be included up to koralz standards
278 
279  polz0=plzapx(.true.,im,np1,np2)
280  IF(rrr(1).LT.polz0) THEN
281  pol1(3)= pol
282  pol2(3)= pol
283  ELSE
284  pol1(3)=-pol
285  pol2(3)=-pol
286  ENDIF
287  ELSEIF(idhep(np1).EQ.-idhep(np2))THEN ! undef orig. only s-dep poss.
288  polz0=plzapx(.true.,im,np1,np2)
289  IF(rrr(1).LT.polz0) THEN
290  pol1(3)= pol
291  pol2(3)= pol
292  ELSE
293  pol1(3)=-pol
294  pol2(3)=-pol
295  ENDIF
296  ELSEIF(abs(idhep(im)).EQ.kfhigch) THEN ! case of charged Higgs
297  pol1(3)= pol
298  pol2(3)= pol
299  ELSE ! case of W+ or W-
300  pol1(3)= -pol
301  pol2(3)= -pol
302  ENDIF
303 c.....include polarisation effect
304  ENDIF
305 
306  IF(idhep(im).EQ.kfhiggs(1).OR.idhep(im).EQ.kfhiggs(2).OR.
307  $ idhep(im).EQ.kfhiggs(3)) THEN
308  IF(idhep(np1).EQ.-kftau .AND.
309  $ (jdahep(1,np1).LE.np1.OR.jdahep(1,np1).GT.nhep) .AND.
310  $ idhep(np2).EQ. kftau .AND.
311  $ (jdahep(1,np2).LE.np2.OR.jdahep(1,np2).GT.nhep)
312  $ ) THEN
313  IF (idhep(im).EQ.kfhiggs(1)) THEN
314  ifpseudo= .false.
315  ELSEIF (idhep(im).EQ.kfhiggs(2)) THEN
316  ifpseudo= .false.
317  ELSEIF (idhep(im).EQ.kfhiggs(3)) THEN
318  ifpseudo= .true.
319  ELSE
320  WRITE(*,*) 'Warning from TAUOLA:'
321  WRITE(*,*) 'I stop this run, wrong IDHEP(IM)=',idhep(im)
322  stop
323  ENDIF
324  CALL spinhiggs(im,np1,np2,ifpseudo,pol1,pol2)
325  IF (ifphot.EQ.1) CALL photos(im) ! Bremsstrahlung in Higgs decay
326  ! AFTER adding taus !!
327 
328 
329  ENDIF
330  ELSE
331  IF(idhep(np1).EQ.-kftau.AND.
332  $ (jdahep(1,np1).LE.np1.OR.jdahep(1,np1).GT.nhep)) THEN
333 c here check on if np1 was not decayed should be verified
334  CALL dexay(1,pol1)
335  IF (ifphot.EQ.1) CALL photos(np1)
336  CALL taupi0(0,1,ion)
337  ENDIF
338 
339  IF(idhep(np2).EQ. kftau.AND.
340  $ (jdahep(1,np2).LE.np2.OR.jdahep(1,np2).GT.nhep)) THEN
341 c here check on if np2 was not decayed should be added
342  CALL dexay(2,pol2)
343  IF (ifphot.EQ.1) CALL photos(np2)
344  CALL taupi0(0,2,ion)
345  ENDIF
346  ENDIF
347  ncount=ncount-2
348  IF (ncount.GT.0) goto 666
349  ENDDO
350 
351  ELSEIF(mode.EQ.1) THEN
352 c ***********************
353 c
354  CALL dexay(100,pol1)
355  CALL dekay(100,pol1x)
356  WRITE(iout,7002)
357  ENDIF
358 c *****
359  7001 FORMAT(///1x,15(5h*****)
360  $ /,' *', 25x,'*****TAUOLA UNIVERSAL INTERFACE: ******',9x,1h*,
361  $ /,' *', 25x,'*****VERSION 1.22, April 2009 (gfort)**',9x,1h*,
362  $ /,' *', 25x,'**AUTHORS: P. Golonka, B. Kersevan, ***',9x,1h*,
363  $ /,' *', 25x,'**T. Pierzchala, E. Richter-Was, ******',9x,1h*,
364  $ /,' *', 25x,'****** Z. Was, M. Worek ***************',9x,1h*,
365  $ /,' *', 25x,'**USEFUL DISCUSSIONS, IN PARTICULAR ***',9x,1h*,
366  $ /,' *', 25x,'*WITH C. Biscarat and S. Slabospitzky**',9x,1h*,
367  $ /,' *', 25x,'****** are warmly acknowledged ********',9x,1h*,
368  $ /,' *', 25x,' ',9x,1h*,
369  $ /,' *', 25x,'********** INITIALIZATION ************',9x,1h*,
370  $ /,' *',f20.5,5x,'tau polarization switch must be 1 or 0 ',9x,1h*,
371  $ /,' *',f20.5,5x,'Higs scalar/pseudo mix CERN-TH/2003-166',9x,1h*,
372  $ /,' *',i10, 15x,'PI0 decay switch must be 1 or 0 ',9x,1h*,
373  $ /,' *',i10, 15x,'ETA decay switch must be 1 or 0 ',9x,1h*,
374  $ /,' *',i10, 15x,'K0S decay switch must be 1 or 0 ',9x,1h*,
375  $ /,1x,15(5h*****)/)
376 
377  7002 FORMAT(///1x,15(5h*****)
378  $ /,' *', 25x,'*****TAUOLA UNIVERSAL INTERFACE: ******',9x,1h*,
379  $ /,' *', 25x,'*****VERSION 1.22, April 2009 (gfort)**',9x,1h*,
380  $ /,' *', 25x,'**AUTHORS: P. Golonka, B. Kersevan, ***',9x,1h*,
381  $ /,' *', 25x,'**T. Pierzchala, E. Richter-Was, ******',9x,1h*,
382  $ /,' *', 25x,'****** Z. Was, M. Worek ***************',9x,1h*,
383  $ /,' *', 25x,'**USEFUL DISCUSSIONS, IN PARTICULAR ***',9x,1h*,
384  $ /,' *', 25x,'*WITH C. Biscarat and S. Slabospitzky**',9x,1h*,
385  $ /,' *', 25x,'****** are warmly acknowledged ********',9x,1h*,
386  $ /,' *', 25x,'****** END OF MODULE OPERATION ********',9x,1h*,
387  $ /,1x,15(5h*****)/)
388 
389  END
390 
391  SUBROUTINE spinhiggs(IM,NP1,NP2,IFPSEUDO,Pol1,Pol2)
392  LOGICAL ifpseudo
393  REAL*8 hh1,hh2,wthiggs
394  dimension pol1(4), pol2(4),hh1(4),hh2(4), rrr(1)
395 ! CALL DEXAY(1,POL1) ! Kept for tests
396 ! CALL DEXAY(2,POL2) ! Kept for tests
397  INTEGER ion(3)
398  10 CONTINUE
399  CALL ranmar(rrr,1)
400  CALL dekay(1,hh1)
401  CALL dekay(2,hh2)
402  wt=wthiggs(ifpseudo,hh1,hh2)
403  IF (rrr(1).GT.wt) goto 10
404  CALL dekay(1+10,hh1)
405  CALL taupi0(0,1,ion)
406  CALL dekay(2+10,hh2)
407  CALL taupi0(0,2,ion)
408  END
409  FUNCTION wthiggs(IFPSEUDO,HH1,HH2)
410  LOGICAL ifpseudo
411  common /pseudocoup/ csc,ssc
412  REAL*4 csc,ssc
413  save pseudocoup
414  REAL*8 hh1(4),hh2(4),r(4,4),wthiggs
415  DO k=1,4
416  DO l=1,4
417  r(k,l)=0
418  ENDDO
419  ENDDO
420  wthiggs=0d0
421 
422  r(4,4)= 1d0 ! unpolarized part
423  r(3,3)=-1d0 ! longitudinal
424  ! other missing
425  IF (ifpseudo) THEN
426  r(1,1)=-1
427  r(2,2)= -1
428  r(1,1)=(csc**2-ssc**2)/(csc**2+ssc**2)
429  r(2,2)=(csc**2-ssc**2)/(csc**2+ssc**2)
430  r(1,2)=2*csc*ssc/(csc**2+ssc**2)
431  r(2,1)=-2*csc*ssc/(csc**2+ssc**2)
432  ELSE
433  r(1,1)=1
434  r(2,2)=1
435  ENDIF
436 
437 
438 
439  DO k=1,4
440  DO l=1,4
441  wthiggs=wthiggs+r(k,l)*hh1(k)*hh2(l)
442  ENDDO
443  ENDDO
444  wthiggs=wthiggs/4d0
445  END
446 
447  FUNCTION plzapx(HOPEin,IM0,NP1,NP2)
448 c im0 np1 np2 are the positions of z/gamma tau tau in hepevt common block.
449 c the purpose of this routine is to calculate polarization of z along tau direction.
450 c this is highly non-trivial due to necessity of reading infromation from hard process
451 c history in hepevt, which is often written not up to the gramatic rules.
452  REAL*8 plzap0,svar,costhe,sini,sfin,zprop2,
453  $ p1(4),p2(4),q1(4),q2(4),qq(4),ph(4),pd1(4),pd2(4),
454  $ pq1(4),pq2(4),pb(4),pa(4)
455  REAL*4 plzapx
456  INTEGER im
457  LOGICAL hope,hopein
458 c this is the hepevt class in old style. no d_h_ class pre-name
459  INTEGER nmxhep
460  parameter(nmxhep=10000)
461  REAL*8 phep, vhep ! to be real*4/ *8 depending on host
462  INTEGER nevhep,nhep,isthep,idhep,jmohep,
463  $ jdahep
464  COMMON /hepevt/
465  $ nevhep, ! serial number
466  $ nhep, ! number of particles
467  $ isthep(nmxhep), ! status code
468  $ idhep(nmxhep), ! particle ident KF
469  $ jmohep(2,nmxhep), ! parent particles
470  $ jdahep(2,nmxhep), ! childreen particles
471  $ phep(5,nmxhep), ! four-momentum, mass [GeV]
472  $ vhep(4,nmxhep) ! vertex [mm]
473 * ----------------------------------------------------------------------
474  LOGICAL qedrad
475  COMMON /phoqed/
476  $ qedrad(nmxhep) ! Photos flag
477 * ----------------------------------------------------------------------
478  SAVE hepevt,phoqed
479 
480 c(bpk)--> brothers of tau already found
481  INTEGER ison
482  COMMON /isons_tau/ison(2)
483 c(bpk)--<
484 c >>
485 c >> step 1: find where are particles in hepevent and pick them up
486 c >>
487  hope=hopein
488 c sometimes shade z of z is its mother ...
489  im=im0
490  im00=jmohep(1,im0)
491 c to protect against check on mother of beam particles.
492  IF (im00.GT.0) THEN
493  IF (idhep(im0).EQ.idhep(im00)) im=jmohep(1,im0)
494  ENDIF
495 c
496 c find(host generator-level) incoming beam-bare-particles which form z and co.
497  imo1=jmohep(1,im)
498  imo2=jmohep(2,im)
499 
500 c(bpk)--> in herwig the POINTER might be to hard cms
501  im00=imo1
502  IF (isthep(im00).EQ.120) THEN
503  imo1=jmohep(1,im00)
504  imo2=jmohep(2,im00)
505  ENDIF
506 c(bpk)--<
507 
508  iffull=0
509 c case when it was like e+e- --> tau+tau- gammas and e+e- were 'first' in hepevt.
510  IF (imo1.EQ.0.AND.imo2.EQ.0) THEN
511  imo1=jmohep(1,np1)
512 c(bpk)-->
513  IF (idhep(imo1).EQ.idhep(np1)) imo1=jmohep(1,imo1) ! PROTECT AGAINST COPIES
514 c(bpk)--<
515  imo2=jmohep(2,np1)
516 c(bpk)-->
517  IF (idhep(imo2).EQ.idhep(np2)) imo2=jmohep(1,imo2) ! PROTECT AGAINST COPIES
518 c(bpk)--<
519  iffull=1
520 c case when it was like qq --> tau+tau- gammas and qq were not 'first' in hepevt.
521 
522  ELSEIF (idhep(im).NE.22.AND.idhep(im).NE.23) THEN
523  imo1=jmohep(1,np1)
524 c(bpk)-->
525  IF (idhep(imo1).EQ.idhep(np1)) imo1=jmohep(1,imo1) ! PROTECT AGAINST COPIES
526 c(bpk)--<
527  imo2=jmohep(2,np1)
528 c(bpk)-->
529  IF (idhep(imo2).EQ.idhep(np2)) imo2=jmohep(1,imo2) ! PROTECT AGAINST COPIES
530 c(bpk)--<
531  iffull=1
532  ENDIF
533 
534 
535 c and check if it really happened
536  IF (imo1.EQ.0) hope=.false.
537  IF (imo2.EQ.0) hope=.false.
538  IF (imo1.EQ.imo2) hope=.false.
539 
540 c
541  DO i=1,4
542  q1(i)= phep(i,np1) !momentum of tau+
543  q2(i)= phep(i,np2) !momentum of tau-
544  ENDDO
545 
546 c corrections due to possible differences in 4-momentum of shadow vs true z.
547 c(bpk)-->
548  IF (im.EQ.jmohep(1,im0).AND.
549  $ (idhep(im).EQ.22.OR.idhep(im).EQ.23)) THEN
550  DO k=1,4
551  pb(k)=phep(k,im)
552  pa(k)=phep(k,im0)
553  ENDDO
554 c(bpk)--<
555 
556  CALL bostdq( 1,pa, q1, q1)
557  CALL bostdq( 1,pa, q2, q2)
558  CALL bostdq(-1,pb, q1, q1)
559  CALL bostdq(-1,pb, q2, q2)
560 
561  ENDIF
562 
563  DO i=1,4
564  qq(i)= q1(i)+q2(i) !momentum of Z
565  IF (hope) p1(i)=phep(i,imo1) !momentum of beam1
566  IF (hope) p2(i)=phep(i,imo2) !momentum of beam2
567  ph(i)=0d0
568  pd1(i)=0d0
569  pd2(i)=0d0
570  ENDDO
571 ! These momenta correspond to quarks, gluons photons or taus
572  idfq1=idhep(np1)
573  idfq2=idhep(np2)
574  IF (hope) idfp1=idhep(imo1)
575  IF (hope) idfp2=idhep(imo2)
576 
577  svar=qq(4)**2-qq(3)**2-qq(2)**2-qq(1)**2
578  IF (.NOT.hope) THEN
579 c options which may be useful in some cases of two heavy boson production
580 c need individual considerations. to be developed.
581 c plzapx=plzap0(11,idfq1,svar,0d0) ! gamma/Z mixture as if produced from e beam
582 c plzapx=plzap0(12,idfq1,svar,0d0) ! pure Z
583  plzapx=0.5 ! pure gamma
584  RETURN
585  ENDIF
586 c >>
587 c >> step 2 look for brothers of z which have to be included in effective incoming particles
588 c >>
589 c let us define beginning and end of particles which are produced in parallel to z
590 c in principle following should work
591 
592 c(bpk)--> accommodate for herwig - im00 points to beam particle or hard cms
593  nx1=jdahep(1,im00)
594  nx2=jdahep(2,im00)
595 c but ...
596  inbr=im ! OK, HARD RECORD Z/GAMMA
597  IF (iffull.EQ.1) inbr=np1 ! OK, NO Z/GAMMA
598  IF (idhep(jmohep(1,inbr)).EQ.idhep(inbr)) inbr=jmohep(1,inbr) ! FORCE HARD RECORD
599 c(bpk)--<
600  IF(nx1.EQ.0.OR.nx2.EQ.0) THEN
601  nx1=inbr
602  nx2=inbr
603  DO k=1,inbr-1
604  IF(jmohep(1,inbr-k).EQ.jmohep(1,inbr)) THEN
605  nx1=inbr-k
606  ELSE
607  goto 7
608  ENDIF
609  ENDDO
610  7 CONTINUE
611 
612  DO k=inbr+1,nhep
613  IF(jmohep(1,k).EQ.jmohep(1,inbr)) THEN
614  nx2=k
615  ELSE
616  goto 8
617  ENDIF
618  ENDDO
619  8 CONTINUE
620  ENDIF
621 
622 c case of annihilation of two bosons is hopeles
623  IF (abs(idfp1).GE.20.AND.abs(idfp2).GE.20) hope=.false.
624 c case of annihilation of non-matching flavors is hopeless
625  IF (abs(idfp1).LE.20.AND.abs(idfp2).LE.20.AND.idfp1+idfp2.NE.0)
626  $ hope=.false.
627  IF (.NOT.hope) THEN
628 c options which may be useful in some cases of two heavy boson production
629 c need individual considerations. to be developed.
630 c plzapx=plzap0(11,idfq1,svar,0d0) ! gamma/Z mixture as if produced from e beam
631 c plzapx=plzap0(12,idfq1,svar,0d0) ! pure Z
632  plzapx=0.5 ! pure gamma
633  RETURN
634  ENDIF
635  IF (abs(idfp1).LT.20) ide= idfp1
636  IF (abs(idfp2).LT.20) ide=-idfp2
637 
638 
639 c >>
640 c >> step 3 we combine gluons, photons into incoming effective beams
641 c >>
642 
643 c in the following we will ignore the possibility of photon emission from taus
644 c however at certain moment it will be necessary to take care of
645 
646  DO l=1,4
647  pd1(l)=p1(l)
648  pd2(l)=p2(l)
649  ENDDO
650 
651  DO l=1,4
652  pq1(l)=q1(l)
653  pq2(l)=q2(l)
654  ENDDO
655 
656  iflav=min(abs(idfp1),abs(idfp2))
657 
658 *--------------------------------------------------------------------------
659 c iflav=min(abs(idfp1),abs(idfp2))
660 c that means that always quark or lepton i.e. process like
661 c f g(gamma) --> f z0 --> tau tau
662 c we glue fermions to effective beams that is f f --> z0 --> tau tau
663 c with gamma/g emission from initial fermion.
664 *---------------------------------------------------------------------------
665 
666  IF (abs(idfp1).GE.20) THEN
667  DO k=nx1,nx2
668  idp=idhep(k)
669  IF (abs(idp).EQ.iflav) THEN
670  DO l=1,4
671  pd1(l)=-phep(l,k)
672  ENDDO
673  ENDIF
674  ENDDO
675  ENDIF
676 
677  IF (abs(idfp2).GE.20) THEN
678  DO k=nx1,nx2
679  idp=idhep(k)
680  IF (abs(idp).EQ.iflav) THEN
681  DO l=1,4
682  pd2(l)=-phep(l,k)
683  ENDDO
684  ENDIF
685  ENDDO
686  ENDIF
687 
688 c if first beam was boson: gluing
689 
690  IF (abs(idfp1).GE.20) THEN
691  DO l=1,4
692  ph(l)=p1(l)
693  ENDDO
694  xm1=abs((pd1(4)+ph(4))**2-(pd1(3)+ph(3))**2
695  $ -(pd1(2)+ph(2))**2-(pd1(1)+ph(1))**2)
696  xm2=abs((pd2(4)+ph(4))**2-(pd2(3)+ph(3))**2
697  $ -(pd2(2)+ph(2))**2-(pd2(1)+ph(1))**2)
698  IF (xm1.LT.xm2) THEN
699  DO l=1,4
700  pd1(l)=pd1(l)+p1(l)
701  ENDDO
702  ELSE
703  DO l=1,4
704  pd2(l)=pd2(l)+p1(l)
705  ENDDO
706  ENDIF
707  ENDIF
708 
709 c if second beam was boson: gluing
710 
711 
712  IF (abs(idfp2).GE.20) THEN
713  DO l=1,4
714  ph(l)=p2(l)
715  ENDDO
716  xm1=abs((pd1(4)+ph(4))**2-(pd1(3)+ph(3))**2
717  $ -(pd1(2)+ph(2))**2-(pd1(1)+ph(1))**2)
718  xm2=abs((pd2(4)+ph(4))**2-(pd2(3)+ph(3))**2
719  $ -(pd2(2)+ph(2))**2-(pd2(1)+ph(1))**2)
720  IF (xm1.LT.xm2) THEN
721  DO l=1,4
722  pd1(l)=pd1(l)+p2(l)
723  ENDDO
724  ELSE
725  DO l=1,4
726  pd2(l)=pd2(l)+p2(l)
727  ENDDO
728  ENDIF
729  ENDIF
730 
731 c now spectators ...
732 
733 c(bpk)-->
734  nph1=np1
735  nph2=np2
736  IF (idhep(jmohep(1,np1)).EQ.idhep(np1)) nph1=jmohep(1,np1) ! SHOULD PUT US IN HARD REC.
737  IF (idhep(jmohep(1,np2)).EQ.idhep(np2)) nph2=jmohep(1,np2) ! SHOULD PUT US IN HARD REC.
738 c(bpk)--<
739 
740  DO k=nx1,nx2
741  IF (abs(idhep(k)).NE.iflav.AND.k.NE.im.AND.
742 c(bpk)-->
743  $ k.NE.nph1.AND.k.NE.nph2) THEN
744 c(bpk)--<
745  IF(idhep(k).EQ.22.AND.iffull.EQ.1) THEN
746  DO l=1,4
747  ph(l)=phep(l,k)
748  ENDDO
749  xm1=abs((pd1(4)-ph(4))**2-(pd1(3)-ph(3))**2
750  $ -(pd1(2)-ph(2))**2-(pd1(1)-ph(1))**2)
751  xm2=abs((pd2(4)-ph(4))**2-(pd2(3)-ph(3))**2
752  $ -(pd2(2)-ph(2))**2-(pd2(1)-ph(1))**2)
753  xm3=abs((pq1(4)+ph(4))**2-(pq1(3)+ph(3))**2
754  $ -(pq1(2)+ph(2))**2-(pq1(1)+ph(1))**2)
755  xm4=abs((pq2(4)+ph(4))**2-(pq2(3)+ph(3))**2
756  $ -(pq2(2)+ph(2))**2-(pq2(1)+ph(1))**2)
757 
758 
759  sini=abs((pd1(4)+pd2(4)-ph(4))**2-(pd1(3)+pd2(3)-ph(3))**2
760  $ -(pd1(2)+pd2(2)-ph(2))**2-(pd1(1)+pd2(1)-ph(1))**2)
761  sfin=abs((pd1(4)+pd2(4) )**2-(pd1(3)+pd2(3) )**2
762  $ -(pd1(2)+pd2(2) )**2-(pd1(1)+pd2(1) )**2)
763 
764  facini=zprop2(sini)
765  facfin=zprop2(sfin)
766 
767  xm1=xm1/facini
768  xm2=xm2/facini
769  xm3=xm3/facfin
770  xm4=xm4/facfin
771 
772  xm=min(xm1,xm2,xm3,xm4)
773  IF (xm1.EQ.xm) THEN
774  DO l=1,4
775  pd1(l)=pd1(l)-ph(l)
776  ENDDO
777  ELSEIF (xm2.EQ.xm) THEN
778  DO l=1,4
779  pd2(l)=pd2(l)-ph(l)
780  ENDDO
781  ELSEIF (xm3.EQ.xm) THEN
782  DO l=1,4
783  q1(l)=pq1(l)+ph(l)
784  ENDDO
785  ELSE
786  DO l=1,4
787  q2(l)=pq2(l)+ph(l)
788  ENDDO
789  ENDIF
790  ELSE
791  DO l=1,4
792  ph(l)=phep(l,k)
793  ENDDO
794  xm1=abs((pd1(4)-ph(4))**2-(pd1(3)-ph(3))**2
795  $ -(pd1(2)-ph(2))**2-(pd1(1)-ph(1))**2)
796  xm2=abs((pd2(4)-ph(4))**2-(pd2(3)-ph(3))**2
797  $ -(pd2(2)-ph(2))**2-(pd2(1)-ph(1))**2)
798  IF (xm1.LT.xm2) THEN
799  DO l=1,4
800  pd1(l)=pd1(l)-ph(l)
801  ENDDO
802  ELSE
803  DO l=1,4
804  pd2(l)=pd2(l)-ph(l)
805  ENDDO
806  ENDIF
807  ENDIF
808  ENDIF
809  ENDDO
810 
811 
812 c >>
813 c >> step 4 look for brothers of tau(sons of z!) which have to be included in
814 c >> effective outcoming taus
815 c >>
816 c let us define beginning and end of particles which are produced in
817 c parallel to tau
818 
819 
820 
821 c find outcoming particles which come from z
822 
823 
824 
825 
826 c(bpk)--> ok, it would have to be along taus in hard record with the same mother
827  IF (abs(idhep(im0)).EQ.22.OR.abs(idhep(im0)).EQ.23) THEN
828  DO k=ison(1),ison(2)
829  IF(abs(idhep(k)).EQ.22) THEN
830 c(bpk)--<
831 
832  do l=1,4
833  ph(l)=phep(l,k)
834  enddo
835 
836  xm3=abs((pq1(4)+ph(4))**2-(pq1(3)+ph(3))**2
837  $ -(pq1(2)+ph(2))**2-(pq1(1)+ph(1))**2)
838  xm4=abs((pq2(4)+ph(4))**2-(pq2(3)+ph(3))**2
839  $ -(pq2(2)+ph(2))**2-(pq2(1)+ph(1))**2)
840 
841  xm=min(xm3,xm4)
842 
843  IF (xm3.EQ.xm) THEN
844  DO l=1,4
845  q1(l)=pq1(l)+ph(l)
846  ENDDO
847  ELSE
848  DO l=1,4
849  q2(l)=pq2(l)+ph(l)
850  ENDDO
851  ENDIF
852  endif
853  enddo
854  ENDIF
855 
856 
857 
858 *------------------------------------------------------------------------
859 
860 
861 c out of effective momenta we calculate costhe and later polarization
862  CALL angulu(pd1,pd2,q1,q2,costhe)
863 
864  plzapx=plzap0(ide,idfq1,svar,costhe)
865  END
866 
867  SUBROUTINE angulu(PD1,PD2,Q1,Q2,COSTHE)
868  REAL*8 pd1(4),pd2(4),q1(4),q2(4),costhe,p(4),qq(4),qt(4)
869 c take effective beam which is less massive, it should be irrelevant
870 c but in case hepevt is particulary dirty may help.
871 c this routine calculate reduced system transver and cosine of scattering
872 c angle.
873 
874  xm1=abs(pd1(4)**2-pd1(3)**2-pd1(2)**2-pd1(1)**2)
875  xm2=abs(pd2(4)**2-pd2(3)**2-pd2(2)**2-pd2(1)**2)
876  IF (xm1.LT.xm2) THEN
877  sign=1d0
878  DO k=1,4
879  p(k)=pd1(k)
880  ENDDO
881  ELSE
882  sign=-1d0
883  DO k=1,4
884  p(k)=pd2(k)
885  ENDDO
886  ENDIF
887 c calculate space like part of p(in z restframe)
888  DO k=1,4
889  qq(k)=q1(k)+q2(k)
890  qt(k)=q1(k)-q2(k)
891  ENDDO
892 
893  xmqq=sqrt(qq(4)**2-qq(3)**2-qq(2)**2-qq(1)**2)
894 
895  qtxqq=qt(4)*qq(4)-qt(3)*qq(3)-qt(2)*qq(2)-qt(1)*qq(1)
896  DO k=1,4
897  qt(k)=qt(k)-qq(k)*qtxqq/xmqq**2
898  ENDDO
899 
900  pxqq=p(4)*qq(4)-p(3)*qq(3)-p(2)*qq(2)-p(1)*qq(1)
901  DO k=1,4
902  p(k)=p(k)-qq(k)*pxqq/xmqq**2
903  ENDDO
904 c calculate costhe
905  pxp =sqrt(p(1)**2+p(2)**2+p(3)**2-p(4)**2)
906  qtxqt=sqrt(qt(3)**2+qt(2)**2+qt(1)**2-qt(4)**2)
907  pxqt =p(3)*qt(3)+p(2)*qt(2)+p(1)*qt(1)-p(4)*qt(4)
908  costhe=pxqt/pxp/qtxqt
909  costhe=costhe*sign
910  END
911 
912  FUNCTION plzap0(IDE,IDF,SVAR,COSTH0)
913 c this function calculates probability for the helicity +1 +1 configuration
914 c of taus for given z/gamma transfer and costh0 cosine of scattering angle
915  REAL*8 plzap0,svar,costhe,costh0,t_born
916 
917  costhe=costh0
918 c >>>>> IF (ide*idf.LT.0) costhe=-costh0 ! this is probably not needed ID
919 c >>>>> of first beam is used by t_giviz0 including sign
920 
921  IF (idf.GT.0) THEN
922  CALL initwk(ide,idf,svar)
923  ELSE
924  CALL initwk(-ide,-idf,svar)
925  ENDIF
926  plzap0=t_born(0,svar,costhe,1d0,1d0)
927  $ /(t_born(0,svar,costhe,1d0,1d0)+t_born(0,svar,costhe,-1d0,-1d0))
928 
929 ! PLZAP0=0.5
930  END
931  FUNCTION t_born(MODE,SVAR,COSTHE,TA,TB)
932 c ----------------------------------------------------------------------
933 c this routine provides born cross section. it has the same
934 c structure as funtis and funtih, thus can be used as simpler
935 c example of the method applied there
936 c input parameters are: svar -- transfer
937 c costhe -- cosine of angle between tau+ and 1st beam
938 c ta,tb -- helicity states of tau+ tau-
939 c
940 c called by : borny, boras, bornv, waga, weight
941 c ----------------------------------------------------------------------
942  IMPLICIT REAL*8(a-h,o-z)
943  COMMON / t_beampm / ene ,amin,amfin,ide,idf
944  REAL*8 ene ,amin,amfin
945  COMMON / t_gauspm /ss,poln,t3e,qe,t3f,qf
946  & ,xupgi ,xupzi ,xupgf ,xupzf
947  & ,ndiag0,ndiaga,keya,keyz
948  & ,itce,jtce,itcf,jtcf,kolor
949  REAL*8 ss,poln,t3e,qe,t3f,qf
950  & ,xupgi(2),xupzi(2),xupgf(2),xupzf(2)
951  REAL*8 seps1,seps2
952 c=====================================================================
953  COMMON / t_gswprm /swsq,amw,amz,amh,amtop,gammz
954  REAL*8 swsq,amw,amz,amh,amtop,gammz
955 c swsq = sin2(theta weinberg)
956 c amw,amz = w & z boson masses respectively
957 c amh = the higgs mass
958 c amtop = the top mass
959 c gammz = z0 width
960  COMPLEX*16 aborn(2,2),aphot(2,2),azett(2,2)
961  COMPLEX*16 xupzfp(2),xupzip(2)
962  COMPLEX*16 abornm(2,2),aphotm(2,2),azettm(2,2)
963  COMPLEX*16 propa,propz
964  COMPLEX*16 xr,xi
965  COMPLEX*16 xupf,xupi,xff(4),xfem,xfota,xrho,xke,xkf,xkef
966  COMPLEX*16 xthing,xve,xvf,xvef
967  DATA xi/(0.d0,1.d0)/,xr/(1.d0,0.d0)/
968  DATA mode0 /-5/
969  DATA ide0 /-55/
970  DATA svar0,cost0 /-5.d0,-6.d0/
971  DATA pi /3.141592653589793238462643d0/
972  DATA seps1,seps2 /0d0,0d0/
973 c
974 c memorization =========================================================
975  IF ( mode.NE.mode0.OR.svar.NE.svar0.OR.costhe.NE.cost0
976  $ .OR.ide0.NE.ide)THEN
977 c
978  keygsw=1
979 c ** propagators
980  ide0=ide
981  mode0=mode
982  svar0=svar
983  cost0=costhe
984  sinthe=sqrt(1.d0-costhe**2)
985  beta=sqrt(max(0d0,1d0-4d0*amfin**2/svar))
986 c i multiply axial coupling by beta factor.
987  xupzfp(1)=0.5d0*(xupzf(1)+xupzf(2))+0.5*beta*(xupzf(1)-xupzf(2))
988  xupzfp(2)=0.5d0*(xupzf(1)+xupzf(2))-0.5*beta*(xupzf(1)-xupzf(2))
989  xupzip(1)=0.5d0*(xupzi(1)+xupzi(2))+0.5*(xupzi(1)-xupzi(2))
990  xupzip(2)=0.5d0*(xupzi(1)+xupzi(2))-0.5*(xupzi(1)-xupzi(2))
991 c final state vector coupling
992  xupf =0.5d0*(xupzf(1)+xupzf(2))
993  xupi =0.5d0*(xupzi(1)+xupzi(2))
994  xthing =0d0
995 
996  propa =1d0/svar
997  propz =1d0/dcmplx(svar-amz**2,svar/amz*gammz)
998  IF (keygsw.EQ.0) propz=0.d0
999  DO 50 i=1,2
1000  DO 50 j=1,2
1001  regula= (3-2*i)*(3-2*j) + costhe
1002  regulm=-(3-2*i)*(3-2*j) * sinthe *2.d0*amfin/sqrt(svar)
1003  aphot(i,j)=propa*(xupgi(i)*xupgf(j)*regula)
1004  azett(i,j)=propz*(xupzip(i)*xupzfp(j)+xthing)*regula
1005  aborn(i,j)=aphot(i,j)+azett(i,j)
1006  aphotm(i,j)=propa*dcmplx(0d0,1d0)*xupgi(i)*xupgf(j)*regulm
1007  azettm(i,j)=propz*dcmplx(0d0,1d0)*(xupzip(i)*xupf+xthing)*regulm
1008  abornm(i,j)=aphotm(i,j)+azettm(i,j)
1009  50 CONTINUE
1010  ENDIF
1011 c
1012 c******************
1013 c* in calculating cross section only diagonal elements
1014 c* of the spin density matrices enter(longitud. pol. only.)
1015 c* helicity conservation explicitly obeyed
1016  polar1= (seps1)
1017  polar2= (-seps2)
1018  born=0d0
1019  DO 150 i=1,2
1020  helic= 3-2*i
1021  DO 150 j=1,2
1022  helit=3-2*j
1023  factor=kolor*(1d0+helic*polar1)*(1d0-helic*polar2)/4d0
1024  factom=factor*(1+helit*ta)*(1-helit*tb)
1025  factor=factor*(1+helit*ta)*(1+helit*tb)
1026 
1027  born=born+cdabs(aborn(i,j))**2*factor
1028 c mass term in born
1029  IF (mode.GE.1) THEN
1030  born=born+cdabs(abornm(i,j))**2*factom
1031  ENDIF
1032 
1033  150 CONTINUE
1034 c************
1035  funt=born
1036  IF(funt.LT.0.d0) funt=born
1037 
1038 c
1039  IF (svar.GT.4d0*amfin**2) THEN
1040 c phase space threshold factor
1041  thresh=sqrt(1-4d0*amfin**2/svar)
1042  t_born= funt*svar**2*thresh
1043  ELSE
1044  thresh=0.d0
1045  t_born=0.d0
1046  ENDIF
1047 c zw here was an error 19. 05. 1989
1048 ! write(*,*) 'KKKK ',PROPA,PROPZ,XUPGI,XUPGF,XUPZI,XUPZF
1049 ! write(*,*) 'KKKK X',svar,costhe,TA,TB,T_BORN
1050  END
1051 
1052  SUBROUTINE initwk(IDEX,IDFX,SVAR)
1053 ! initialization routine coupling masses etc.
1054  IMPLICIT REAL*8 (a-h,o-z)
1055  COMMON / t_beampm / ene ,amin,amfin,ide,idf
1056  REAL*8 ene ,amin,amfin
1057  COMMON / t_gauspm /ss,poln,t3e,qe,t3f,qf
1058  & ,xupgi ,xupzi ,xupgf ,xupzf
1059  & ,ndiag0,ndiaga,keya,keyz
1060  & ,itce,jtce,itcf,jtcf,kolor
1061  REAL*8 ss,poln,t3e,qe,t3f,qf
1062  & ,xupgi(2),xupzi(2),xupgf(2),xupzf(2)
1063  COMMON / t_gswprm /swsq,amw,amz,amh,amtop,gammz
1064  REAL*8 swsq,amw,amz,amh,amtop,gammz
1065 c swsq = sin2(theta weinberg)
1066 c amw,amz = w & z boson masses respectively
1067 c amh = the higgs mass
1068 c amtop = the top mass
1069 c gammz = z0 width
1070 c
1071  ene=sqrt(svar)/2
1072  amin=0.511d-3
1073  swsq=0.23147
1074  amz=91.1882
1075  gammz=2.4952
1076  IF (idfx.EQ. 15) then
1077  idf=2 ! denotes tau +2 tau-
1078  amfin=1.77703 !this mass is irrelevant if small, used in ME only
1079  ELSEIF (idfx.EQ.-15) then
1080  idf=-2 ! denotes tau -2 tau-
1081  amfin=1.77703 !this mass is irrelevant if small, used in ME only
1082  ELSE
1083  WRITE(*,*) 'INITWK: WRONG IDFX'
1084  stop
1085  ENDIF
1086 
1087  IF (idex.EQ. 11) then !electron
1088  ide= 2
1089  amin=0.511d-3
1090  ELSEIF (idex.EQ.-11) then !positron
1091  ide=-2
1092  amin=0.511d-3
1093  ELSEIF (idex.EQ. 13) then !mu+
1094  ide= 2
1095  amin=0.105659
1096  ELSEIF (idex.EQ.-13) then !mu-
1097  ide=-2
1098  amin=0.105659
1099  ELSEIF (idex.EQ. 1) then !d
1100  ide= 4
1101  amin=0.05
1102  ELSEIF (idex.EQ.- 1) then !d~
1103  ide=-4
1104  amin=0.05
1105  ELSEIF (idex.EQ. 2) then !u
1106  ide= 3
1107  amin=0.02
1108  ELSEIF (idex.EQ.- 2) then !u~
1109  ide=-3
1110  amin=0.02
1111  ELSEIF (idex.EQ. 3) then !s
1112  ide= 4
1113  amin=0.3
1114  ELSEIF (idex.EQ.- 3) then !s~
1115  ide=-4
1116  amin=0.3
1117  ELSEIF (idex.EQ. 4) then !c
1118  ide= 3
1119  amin=1.3
1120  ELSEIF (idex.EQ.- 4) then !c~
1121  ide=-3
1122  amin=1.3
1123  ELSEIF (idex.EQ. 5) then !b
1124  ide= 4
1125  amin=4.5
1126  ELSEIF (idex.EQ.- 5) then !b~
1127  ide=-4
1128  amin=4.5
1129  ELSEIF (idex.EQ. 12) then !nu_e
1130  ide= 1
1131  amin=0.1d-3
1132  ELSEIF (idex.EQ.- 12) then !nu_e~
1133  ide=-1
1134  amin=0.1d-3
1135  ELSEIF (idex.EQ. 14) then !nu_mu
1136  ide= 1
1137  amin=0.1d-3
1138  ELSEIF (idex.EQ.- 14) then !nu_mu~
1139  ide=-1
1140  amin=0.1d-3
1141  ELSEIF (idex.EQ. 16) then !nu_tau
1142  ide= 1
1143  amin=0.1d-3
1144  ELSEIF (idex.EQ.- 16) then !nu_tau~
1145  ide=-1
1146  amin=0.1d-3
1147 
1148  ELSE
1149  WRITE(*,*) 'INITWK: WRONG IDEX'
1150  stop
1151  ENDIF
1152 
1153 c ----------------------------------------------------------------------
1154 c
1155 c initialisation of coupling constants and fermion-gamma / z0 vertex
1156 c
1157 c called by : koralz
1158 c ----------------------------------------------------------------------
1159  itce=ide/iabs(ide)
1160  jtce=(1-itce)/2
1161  itcf=idf/iabs(idf)
1162  jtcf=(1-itcf)/2
1163  CALL t_givizo( ide, 1,aizor,qe,kdumm)
1164  CALL t_givizo( ide,-1,aizol,qe,kdumm)
1165  xupgi(1)=qe
1166  xupgi(2)=qe
1167  t3e = aizol+aizor
1168  xupzi(1)=(aizor-qe*swsq)/sqrt(swsq*(1-swsq))
1169  xupzi(2)=(aizol-qe*swsq)/sqrt(swsq*(1-swsq))
1170  CALL t_givizo( idf, 1,aizor,qf,kolor)
1171  CALL t_givizo( idf,-1,aizol,qf,kolor)
1172  xupgf(1)=qf
1173  xupgf(2)=qf
1174  t3f = aizol+aizor
1175  xupzf(1)=(aizor-qf*swsq)/sqrt(swsq*(1-swsq))
1176  xupzf(2)=(aizol-qf*swsq)/sqrt(swsq*(1-swsq))
1177 c
1178  ndiag0=2
1179  ndiaga=11
1180  keya = 1
1181  keyz = 1
1182 c
1183 c
1184  RETURN
1185  END
1186 
1187  SUBROUTINE t_givizo(IDFERM,IHELIC,SIZO3,CHARGE,KOLOR)
1188 c ----------------------------------------------------------------------
1189 c provides electric charge and weak izospin of a family fermion
1190 c idferm=1,2,3,4 denotes neutrino, lepton, up and down quark
1191 c negative idferm=-1,-2,-3,-4, denotes antiparticle
1192 c ihelic=+1,-1 denotes right and left handednes( chirality)
1193 c sizo3 is third projection of weak izospin(plus minus half)
1194 c and charge is electric charge in units of electron charge
1195 c kolor is a qcd colour, 1 for lepton, 3 for quarks
1196 c
1197 c called by : evente, eventm, funtih, .....
1198 c ----------------------------------------------------------------------
1199  IMPLICIT REAL*8(a-h,o-z)
1200 c
1201  IF(idferm.EQ.0.OR.iabs(idferm).GT.4) goto 901
1202  IF(iabs(ihelic).NE.1) goto 901
1203  ih =ihelic
1204  idtype =iabs(idferm)
1205  ic =idferm/idtype
1206  lepqua=int(idtype*0.4999999d0)
1207  iupdow=idtype-2*lepqua-1
1208  charge =(-iupdow+2d0/3d0*lepqua)*ic
1209  sizo3 =0.25d0*(ic-ih)*(1-2*iupdow)
1210  kolor=1+2*lepqua
1211 c** note that conventionaly z0 coupling is
1212 c** xoupz=(sizo3-charge*swsq)/sqrt(swsq*(1-swsq))
1213  RETURN
1214  901 print *,' STOP IN GIVIZO: WRONG PARAMS.'
1215  stop
1216  END
1217  SUBROUTINE phyfix(NSTOP,NSTART)
1218  common/lujets/n,k(4000,5),p(4000,5),v(4000,5)
1219  SAVE /lujets/
1220 c nstop nstart : when phytia history ends and event starts.
1221  nstop=0
1222  nstart=1
1223  DO i=1, n
1224  IF(k(i,1).NE.21) THEN
1225  nstop = i-1
1226  nstart= i
1227  goto 500
1228  ENDIF
1229  ENDDO
1230  500 CONTINUE
1231  END
1232  SUBROUTINE filhep(N,IST,ID,JMO1,JMO2,JDA1,JDA2,P4,PINV,PHFLAG)
1233 c ----------------------------------------------------------------------
1234 c this subroutine fills one entry into the hepevt common
1235 c and updates the information for affected mother entries
1236 c
1237 c written by martin w. gruenewald(91/01/28)
1238 c
1239 c called by : ztohep,btohep,dwluxy
1240 c ----------------------------------------------------------------------
1241 c
1242 c this is the hepevt class in old style. no d_h_ class pre-name
1243 c this is the hepevt class in old style. no d_h_ class pre-name
1244  INTEGER nmxhep
1245  parameter(nmxhep=10000)
1246  REAL*8 phep, vhep ! to be real*4/ *8 depending on host
1247  INTEGER nevhep,nhep,isthep,idhep,jmohep,
1248  $ jdahep
1249  COMMON /hepevt/
1250  $ nevhep, ! serial number
1251  $ nhep, ! number of particles
1252  $ isthep(nmxhep), ! status code
1253  $ idhep(nmxhep), ! particle ident KF
1254  $ jmohep(2,nmxhep), ! parent particles
1255  $ jdahep(2,nmxhep), ! childreen particles
1256  $ phep(5,nmxhep), ! four-momentum, mass [GeV]
1257  $ vhep(4,nmxhep) ! vertex [mm]
1258 * ----------------------------------------------------------------------
1259  LOGICAL qedrad
1260  COMMON /phoqed/
1261  $ qedrad(nmxhep) ! Photos flag
1262 * ----------------------------------------------------------------------
1263  SAVE hepevt,phoqed
1264  LOGICAL phflag
1265 c
1266  REAL*4 p4(4)
1267 c
1268 c check address mode
1269  IF (n.EQ.0) THEN
1270 c
1271 c append mode
1272  ihep=nhep+1
1273  ELSE IF (n.GT.0) THEN
1274 c
1275 c absolute position
1276  ihep=n
1277  ELSE
1278 c
1279 c relative position
1280  ihep=nhep+n
1281  END IF
1282 c
1283 c check on ihep
1284  IF ((ihep.LE.0).OR.(ihep.GT.nmxhep)) RETURN
1285 c
1286 c add entry
1287  nhep=ihep
1288  isthep(ihep)=ist
1289  idhep(ihep)=id
1290  jmohep(1,ihep)=jmo1
1291  IF(jmo1.LT.0)jmohep(1,ihep)=jmohep(1,ihep)+ihep
1292  jmohep(2,ihep)=jmo2
1293  IF(jmo2.LT.0)jmohep(2,ihep)=jmohep(2,ihep)+ihep
1294  jdahep(1,ihep)=jda1
1295  jdahep(2,ihep)=jda2
1296 c
1297  DO i=1,4
1298  phep(i,ihep)=p4(i)
1299 c
1300 c koral-b and koral-z do not provide vertex and/or lifetime informations
1301  vhep(i,ihep)=0.0
1302  END DO
1303  phep(5,ihep)=pinv
1304 c flag for photos...
1305  qedrad(ihep)=phflag
1306 c
1307 c update process:
1308  DO ip=jmohep(1,ihep),jmohep(2,ihep)
1309  IF(ip.GT.0)THEN
1310 c
1311 c if there is a daughter at ihep, mother entry at ip has decayed
1312  IF(isthep(ip).EQ.1)isthep(ip)=2
1313 c
1314 c and daughter pointers of mother entry must be updated
1315  IF(jdahep(1,ip).EQ.0)THEN
1316  jdahep(1,ip)=ihep
1317  jdahep(2,ip)=ihep
1318  ELSE
1319  jdahep(2,ip)=max(ihep,jdahep(2,ip))
1320  END IF
1321  END IF
1322  END DO
1323 c
1324  RETURN
1325  END
1326 
1327 
1328  FUNCTION ihepdim(DUM)
1329 c this is the hepevt class in old style. no d_h_ class pre-name
1330 c this is the hepevt class in old style. no d_h_ class pre-name
1331  INTEGER nmxhep
1332  parameter(nmxhep=10000)
1333  REAL*8 phep, vhep ! to be real*4/ *8 depending on host
1334  INTEGER nevhep,nhep,isthep,idhep,jmohep,
1335  $ jdahep
1336  COMMON /hepevt/
1337  $ nevhep, ! serial number
1338  $ nhep, ! number of particles
1339  $ isthep(nmxhep), ! status code
1340  $ idhep(nmxhep), ! particle ident KF
1341  $ jmohep(2,nmxhep), ! parent particles
1342  $ jdahep(2,nmxhep), ! childreen particles
1343  $ phep(5,nmxhep), ! four-momentum, mass [GeV]
1344  $ vhep(4,nmxhep) ! vertex [mm]
1345 * ----------------------------------------------------------------------
1346  LOGICAL qedrad
1347  COMMON /phoqed/
1348  $ qedrad(nmxhep) ! Photos flag
1349 * ----------------------------------------------------------------------
1350  SAVE hepevt,phoqed
1351  ihepdim=nhep
1352  END
1353  FUNCTION zprop2(S)
1354  IMPLICIT REAL*8(a-h,o-z)
1355  COMPLEX*16 cprz0,cprz0m
1356  amz=91.1882
1357  gammz=2.49
1358  cprz0=dcmplx((s-amz**2),s/amz*gammz)
1359  cprz0m=1/cprz0
1360  zprop2=(abs(cprz0m))**2
1361  END
1362 
1363  SUBROUTINE taupi0(MODE,JAK,ION)
1364 c no initialization required. must be called once after every:
1365 c 1) CALL dekay(1+10,...)
1366 c 2) CALL dekay(2+10,...)
1367 c 3) CALL dexay(1,...)
1368 c 4) CALL dexay(2,...)
1369 c subroutine to decay originating from tauola's taus:
1370 C 1) etas (with CALL TAUETA(JAK))
1371 C 2) later pi0's from taus.
1372 c 3) extensions to other applications possible.
1373 c this routine belongs to >tauola universal interface<, but uses
1374 c routines from >tauola< utilities as well. 25.08.2005
1375 c this is the hepevt class in old style. no d_h_ class pre-name
1376  INTEGER nmxhep
1377  parameter(nmxhep=10000)
1378  REAL*8 phep, vhep ! to be real*4/ *8 depending on host
1379  INTEGER nevhep,nhep,isthep,idhep,jmohep,
1380  $ jdahep
1381  COMMON /hepevt/
1382  $ nevhep, ! serial number
1383  $ nhep, ! number of particles
1384  $ isthep(nmxhep), ! status code
1385  $ idhep(nmxhep), ! particle ident KF
1386  $ jmohep(2,nmxhep), ! parent particles
1387  $ jdahep(2,nmxhep), ! childreen particles
1388  $ phep(5,nmxhep), ! four-momentum, mass [GeV]
1389  $ vhep(4,nmxhep) ! vertex [mm]
1390 * ----------------------------------------------------------------------
1391  LOGICAL qedrad
1392  COMMON /phoqed/
1393  $ qedrad(nmxhep) ! Photos flag
1394 * ----------------------------------------------------------------------
1395  SAVE hepevt,phoqed
1396 
1397 c position of taus, must be defined by host program:
1398  COMMON /taupos/ np1,np2
1399 c
1400  REAL phot1(4),phot2(4)
1401  REAL*8 r,x(4),y(4),pi0(4)
1402  INTEGER jezeli(3),ion(3)
1403  DATA jezeli /0,0,0/
1404  SAVE jezeli
1405  IF (mode.EQ.-1) THEN
1406  jezeli(1)=ion(1)
1407  jezeli(2)=ion(2)
1408  jezeli(3)=ion(3)
1409  RETURN
1410  ENDIF
1411  IF (jezeli(1).EQ.0) RETURN
1412  IF (jezeli(2).EQ.1) CALL taueta(jak)
1413  IF (jezeli(3).EQ.1) CALL tauk0s(jak)
1414 c position of decaying particle:
1415  IF((kto.EQ. 1).OR.(kto.EQ.11)) THEN
1416  nps=np1
1417  ELSE
1418  nps=np2
1419  ENDIF
1420  nhepm=nhep ! to avoid infinite loop
1421  DO k=jdahep(1,nps),nhepm ! we search for pi0's from tau till eor.
1422  IF (idhep(k).EQ.111.AND.jdahep(1,k).LE.k) THEN ! IF we found pi0
1423  DO l=1,4
1424  pi0(l)= phep(l,k)
1425  ENDDO
1426 ! random 3 vector on the sphere, masless
1427  r=sqrt(pi0(4)**2-pi0(3)**2-pi0(2)**2-pi0(1)**2)/2d0
1428  CALL spherd(r,x)
1429  x(4)=r
1430  y(4)=r
1431 
1432  y(1)=-x(1)
1433  y(2)=-x(2)
1434  y(3)=-x(3)
1435 ! boost to lab and to real*4
1436  CALL bostdq(-1,pi0,x,x)
1437  CALL bostdq(-1,pi0,y,y)
1438  DO l=1,4
1439  phot1(l)=x(l)
1440  phot2(l)=y(l)
1441  ENDDO
1442 c to hepevt
1443  CALL filhep(0,1,22,k,k,0,0,phot1,0.0,.true.)
1444  CALL filhep(0,1,22,k,k,0,0,phot2,0.0,.true.)
1445  ENDIF
1446  ENDDO
1447 c
1448  END
1449  SUBROUTINE taueta(JAK)
1450 c subroutine to decay etas's from taus.
1451 C this routine belongs to tauola universal interface, but uses
1452 C routines from tauola utilities. Just flat phase space, but 4 channels.
1453 C it is called at the beginning of SUBR. TAUPI0(JAK)
1454 C and as far as hepevt search it is basically the same as TAUPI0. 25.08.2005
1455 C this is the hepevt class in old style. No d_h_ class pre-name
1456  INTEGER NMXHEP
1457  PARAMETER (NMXHEP=10000)
1458  REAL*8 phep, vhep ! to be real*4/ *8 depending on host
1459  INTEGER nevhep,nhep,isthep,idhep,jmohep,
1460  $ jdahep
1461  COMMON /hepevt/
1462  $ nevhep, ! serial number
1463  $ nhep, ! number of particles
1464  $ isthep(nmxhep), ! status code
1465  $ idhep(nmxhep), ! particle ident KF
1466  $ jmohep(2,nmxhep), ! parent particles
1467  $ jdahep(2,nmxhep), ! childreen particles
1468  $ phep(5,nmxhep), ! four-momentum, mass [GeV]
1469  $ vhep(4,nmxhep) ! vertex [mm]
1470 * ----------------------------------------------------------------------
1471  LOGICAL qedrad
1472  COMMON /phoqed/
1473  $ qedrad(nmxhep) ! Photos flag
1474 * ----------------------------------------------------------------------
1475  SAVE hepevt,phoqed
1476  COMMON / PARMAS / AMTAU,AMNUTA,AMEL,AMNUE,AMMU,AMNUMU
1477  * ,AMPIZ,AMPI,AMRO,GAMRO,AMA1,GAMA1
1478  * ,AMK,AMKZ,AMKST,GAMKST
1479 *
1480  REAL*4 AMTAU,AMNUTA,AMEL,AMNUE,AMMU,AMNUMU
1481  * ,AMPIZ,AMPI,AMRO,GAMRO,AMA1,GAMA1
1482  * ,AMK,AMKZ,AMKST,GAMKST
1483 
1484 C position of taus, must be defined by host program:
1485  COMMON /TAUPOS/ NP1,NP2
1486 c
1487  REAL RRR(1),BRSUM(3), RR(2)
1488  REAL PHOT1(4),PHOT2(4),PHOT3(4)
1489  REAL*8 X(4), Y(4), Z(4)
1490  REAL YM1,YM2,YM3
1491  REAL*8 R,RU,PETA(4),XM1,XM2,XM3,XM,XLAM
1492  REAL*8 a,b,c
1493  XLAM(a,b,c)=SQRT(ABS((a-b-c)**2-4.0*b*c))
1494 C position of decaying particle:
1495 .EQ..OR..EQ. IF((KTO 1)(KTO11)) THEN
1496  NPS=NP1
1497  ELSE
1498  NPS=NP2
1499  ENDIF
1500  nhepM=nhep ! to avoid infinite loop
1501  DO K=JDAHEP(1,NPS),nhepM ! we search for etas's from tau till eor.
1502  IF (idhep(k).EQ.221.AND.jdahep(1,k).LE.k) THEN ! IF we found eta
1503  DO l=1,4
1504  peta(l)= phep(l,k) ! eta 4 momentum
1505  ENDDO
1506 c eta cumulated branching ratios:
1507  brsum(1)=0.389 ! gamma gamma
1508  brsum(2)=brsum(1)+0.319 ! 3 pi0
1509  brsum(3)=brsum(2)+0.237 ! pi+ pi- pi0 rest is thus pi+pi-gamma
1510  CALL ranmar(rrr,1)
1511 
1512  IF (rrr(1).LT.brsum(1)) THEN ! gamma gamma channel exactly like pi0
1513 ! random 3 vector on the sphere, masless
1514  r=sqrt(peta(4)**2-peta(3)**2-peta(2)**2-peta(1)**2)/2d0
1515  CALL spherd(r,x)
1516  x(4)=r
1517  y(4)=r
1518 
1519  y(1)=-x(1)
1520  y(2)=-x(2)
1521  y(3)=-x(3)
1522 ! boost to lab and to real*4
1523  CALL bostdq(-1,peta,x,x)
1524  CALL bostdq(-1,peta,y,y)
1525  DO l=1,4
1526  phot1(l)=x(l)
1527  phot2(l)=y(l)
1528  ENDDO
1529 c to hepevt
1530  CALL filhep(0,1,22,k,k,0,0,phot1,0.0,.true.)
1531  CALL filhep(0,1,22,k,k,0,0,phot2,0.0,.true.)
1532  ELSE ! 3 body channels
1533  IF(rrr(1).LT.brsum(2)) THEN ! 3 pi0
1534  id1= 111
1535  id2= 111
1536  id3= 111
1537  xm1=ampiz ! masses
1538  xm2=ampiz
1539  xm3=ampiz
1540  ELSEIF(rrr(1).LT.brsum(3)) THEN ! pi+ pi- pi0
1541  id1= 211
1542  id2=-211
1543  id3= 111
1544  xm1=ampi ! masses
1545  xm2=ampi
1546  xm3=ampiz
1547  ELSE ! pi+ pi- gamma
1548  id1= 211
1549  id2=-211
1550  id3= 22
1551  xm1=ampi ! masses
1552  xm2=ampi
1553  xm3=0.0
1554  ENDIF
1555  7 CONTINUE ! we generate mass of the first pair:
1556  CALL ranmar(rr,2)
1557  r=sqrt(peta(4)**2-peta(3)**2-peta(2)**2-peta(1)**2)
1558  amin=xm1+xm2
1559  amax=r-xm3
1560  am2=sqrt(amin**2+rr(1)*(amax**2-amin**2))
1561 c weight for flat phase space
1562  wt=xlam(1d0*r**2,1d0*am2**2,1d0*xm3**2)
1563  & *xlam(1d0*am2**2,1d0*xm1**2,1d0*xm2**2)
1564  & /r**2 /am2**2
1565  IF (rr(2).GT.wt) goto 7
1566 
1567  ru=xlam(1d0*am2**2,1d0*xm1**2,1d0*xm2**2)/am2/2 ! momenta of the
1568  ! first two products
1569  ! in the rest frame of that pair
1570  CALL spherd(ru,x)
1571  x(4)=sqrt(ru**2+xm1**2)
1572  y(4)=sqrt(ru**2+xm2**2)
1573 
1574  y(1)=-x(1)
1575  y(2)=-x(2)
1576  y(3)=-x(3)
1577 c generate momentum of that pair in rest frame of eta:
1578  ru=xlam(1d0*r**2,1d0*am2**2,1d0*xm3**2)/r/2
1579  CALL spherd(ru,z)
1580  z(4)=sqrt(ru**2+am2**2)
1581 c and boost first two decay products to rest frame of eta.
1582  CALL bostdq(-1,z,x,x)
1583  CALL bostdq(-1,z,y,y)
1584 c redefine z(4) to 4-momentum of the last decay product:
1585  z(1)=-z(1)
1586  z(2)=-z(2)
1587  z(3)=-z(3)
1588  z(4)=sqrt(ru**2+xm3**2)
1589 c boost all to lab and move to real*4; also masses
1590  CALL bostdq(-1,peta,x,x)
1591  CALL bostdq(-1,peta,y,y)
1592  CALL bostdq(-1,peta,z,z)
1593  DO l=1,4
1594  phot1(l)=x(l)
1595  phot2(l)=y(l)
1596  phot3(l)=z(l)
1597  ENDDO
1598  ym1=xm1
1599  ym2=xm2
1600  ym3=xm3
1601 c to hepevt
1602  CALL filhep(0,1,id1,k,k,0,0,phot1,ym1,.true.)
1603  CALL filhep(0,1,id2,k,k,0,0,phot2,ym2,.true.)
1604  CALL filhep(0,1,id3,k,k,0,0,phot3,ym3,.true.)
1605  ENDIF
1606 
1607  ENDIF
1608  ENDDO
1609 c
1610  END
1611  SUBROUTINE tauk0s(JAK)
1612 c subroutine to decay k0s's from taus.
1613 C this routine belongs to tauola universal interface, but uses
1614 C routines from tauola utilities. Just flat phase space, but 4 channels.
1615 C it is called at the beginning of SUBR. TAUPI0(JAK)
1616 C and as far as hepevt search it is basically the same as TAUPI0. 25.08.2005
1617 C this is the hepevt class in old style. No d_h_ class pre-name
1618  INTEGER NMXHEP
1619  PARAMETER (NMXHEP=10000)
1620  REAL*8 phep, vhep ! to be real*4/ *8 depending on host
1621  INTEGER nevhep,nhep,isthep,idhep,jmohep,
1622  $ jdahep
1623  COMMON /hepevt/
1624  $ nevhep, ! serial number
1625  $ nhep, ! number of particles
1626  $ isthep(nmxhep), ! status code
1627  $ idhep(nmxhep), ! particle ident KF
1628  $ jmohep(2,nmxhep), ! parent particles
1629  $ jdahep(2,nmxhep), ! childreen particles
1630  $ phep(5,nmxhep), ! four-momentum, mass [GeV]
1631  $ vhep(4,nmxhep) ! vertex [mm]
1632 * ----------------------------------------------------------------------
1633  LOGICAL qedrad
1634  COMMON /phoqed/
1635  $ qedrad(nmxhep) ! Photos flag
1636 * ----------------------------------------------------------------------
1637  SAVE hepevt,phoqed
1638 
1639  COMMON / PARMAS / AMTAU,AMNUTA,AMEL,AMNUE,AMMU,AMNUMU
1640  * ,AMPIZ,AMPI,AMRO,GAMRO,AMA1,GAMA1
1641  * ,AMK,AMKZ,AMKST,GAMKST
1642 *
1643  REAL*4 AMTAU,AMNUTA,AMEL,AMNUE,AMMU,AMNUMU
1644  * ,AMPIZ,AMPI,AMRO,GAMRO,AMA1,GAMA1
1645  * ,AMK,AMKZ,AMKST,GAMKST
1646 
1647 C position of taus, must be defined by host program:
1648  COMMON /TAUPOS/ NP1,NP2
1649 c
1650  REAL RRR(1),BRSUM(3), RR(2)
1651  REAL PHOT1(4),PHOT2(4),PHOT3(4)
1652  REAL*8 X(4), Y(4), Z(4)
1653  REAL YM1,YM2,YM3
1654  REAL*8 R,RU,PETA(4),XM1,XM2,XM3,XM,XLAM
1655  REAL*8 a,b,c
1656  XLAM(a,b,c)=SQRT(ABS((a-b-c)**2-4.0*b*c))
1657 C position of decaying particle:
1658 .EQ..OR..EQ. IF((KTO 1)(KTO11)) THEN
1659  NPS=NP1
1660  ELSE
1661  NPS=NP2
1662  ENDIF
1663  nhepM=nhep ! to avoid infinite loop
1664  DO K=JDAHEP(1,NPS),nhepM ! we search for K0S's from tau till eor.
1665  IF (idhep(k).EQ.310.AND.jdahep(1,k).LE.k) THEN ! IF we found K0S
1666 
1667 
1668  DO l=1,4
1669  peta(l)= phep(l,k) ! K0S 4 momentum (this is cloned from eta decay)
1670  ENDDO
1671 c k0s cumulated branching ratios:
1672  brsum(1)=0.313 ! 2 PI0
1673  brsum(2)=1.0 ! BRSUM(1)+0.319 ! Pi+ PI-
1674  brsum(3)=brsum(2)+0.237 ! pi+ pi- pi0 rest is thus pi+pi-gamma
1675  CALL ranmar(rrr,1)
1676 
1677  IF(rrr(1).LT.brsum(1)) THEN ! 2 pi0
1678  id1= 111
1679  id2= 111
1680  xm1=ampiz ! masses
1681  xm2=ampiz
1682  ELSEIF(rrr(1).LT.brsum(2)) THEN ! pi+ pi-
1683  id1= 211
1684  id2=-211
1685  xm1=ampi ! masses
1686  xm2=ampi
1687  ELSE ! gamma gamma unused !!!
1688  id1= 22
1689  id2= 22
1690  xm1= 0.0 ! masses
1691  xm2= 0.0
1692  ENDIF
1693 
1694 ! random 3 vector on the sphere, of equal mass !!
1695  r=sqrt(peta(4)**2-peta(3)**2-peta(2)**2-peta(1)**2)/2d0
1696  r4=r
1697  r=sqrt(abs(r**2-xm1**2))
1698  CALL spherd(r,x)
1699  x(4)=r4
1700  y(4)=r4
1701 
1702  y(1)=-x(1)
1703  y(2)=-x(2)
1704  y(3)=-x(3)
1705 ! boost to lab and to real*4
1706  CALL bostdq(-1,peta,x,x)
1707  CALL bostdq(-1,peta,y,y)
1708  DO l=1,4
1709  phot1(l)=x(l)
1710  phot2(l)=y(l)
1711  ENDDO
1712 
1713  ym1=xm1
1714  ym2=xm2
1715 c to hepevt
1716  CALL filhep(0,1,id1,k,k,0,0,phot1,ym1,.true.)
1717  CALL filhep(0,1,id2,k,k,0,0,phot2,ym2,.true.)
1718 
1719 c
1720  ENDIF
1721  ENDDO
1722 
1723  END