C++ Interface to Tauola
tauola-BBB/jetset-F/tauface-jetset.F
1  SUBROUTINE tauola(MODE,KEYPOL)
2 C *************************************
3 C general tauola interface, should work in every case until
4 C hepevt is OK, does not check if hepevt is 'clean'
5 C in particular will decay decayed taus...
6 C only longitudinal spin effects are included.
7 C in W decay v-a vertex is assumed
8 C date: 12 DEC 1998. date: 21 June 1999. date: 24 Jan 2001 date: 24 Aug 2001
9 #include "../../include/HEPEVT.h"
10  COMMON /taupos/ np1, np2
11  real*4 phoi(4),phof(4)
12  double precision Q1(4),Q2(4),P1(4),P2(4),P3(4),P4(4)
13  COMMON / momdec / q1,q2,p1,p2,p3,p4
14 * tauola, photos and jetset overall switches
15  COMMON /libra/ jak1,jak2,itdkrc,ifphot,ifhadm,ifhadp
16  real*4 rrr(1)
17  LOGICAL IFPSEUDO
18  common /pseudocoup/ csc,ssc
19  real*4 csc,ssc
20  save pseudocoup
21  COMMON / inout / inut,iout
22 
23 C to switch tau polarization OFF in taus
24  dimension pol1(4), pol2(4)
25  double precision POL1x(4), POL2x(4)
26  INTEGER ION(3)
27  DATA pol1 /0.0,0.0,0.0,0.0/
28  DATA pol2 /0.0,0.0,0.0,0.0/
29  DATA pi /3.141592653589793238462643d0/
30 
31 C store decay vertexes
32  dimension imother(20)
33  INTEGER KFHIGGS(3)
34 
35 C store daughter pointers
36  INTEGER ISON
37  COMMON /isons_tau/ison(2)
38  SAVE /isons_tau/
39 
40  IF(mode.EQ.-1) THEN
41 C ***********************
42 
43  jak1 = 0 ! decay mode first tau
44  jak2 = 0 ! decay mode second tau
45  itdkrc=1.0 ! switch of radiative corrections in decay
46  ifphot=1.0 ! PHOTOS switch
47  ifhadm=1.0
48  ifhadp=1.0
49  pol=1.0 ! tau polarization dipswitch must be 1 or 0
50 
51  kfhiggs(1) = 25
52  kfhiggs(2) = 35
53  kfhiggs(3) = 36
54  kfhigch = 37
55  kfz0 = 23
56  kfgam = 22
57  kftau = 15
58  kfnue = 16
59 C couplings of the 'pseudoscalar higgs' as in CERN-TH/2003-166
60  psi=0.5*pi ! 0.15*PI
61  xmtau=1.777 ! tau mass
62  xmh=120 ! higgs boson mass
63  betah=sqrt(1d0-4*xmtau**2/xmh**2)
64  csc=cos(psi)*betah
65  ssc=sin(psi)
66 C write(*,*) ' scalar component=',csc,' pseudo-scalar component=',ssc
67  IF (ifphot.EQ.1) CALL phoini ! this if PHOTOS was not initialized earlier
68  CALL inietc(jak1,jak2,itdkrc,ifphot)
69  CALL inimas
70  CALL iniphx(0.01d0)
71  CALL initdk
72 C activation of pi0 and eta decays: (-1,1) means on, (-1,0) off
73  ion(1)=1
74  ion(2)=1
75  ion(3)=1
76  CALL taupi0(-1,1,ion)
77  CALL dekay(-1,pol1x)
78  WRITE(iout,7001) pol,psi,ion(1),ion(2),ion(3)
79  ELSEIF(mode.EQ.0) THEN
80 C ***********************
81 C
82 C..... find tau-s and fill common block /TAUPOS/
83 C this is to avoid LUND history fillings. This call is optional
84  CALL phyfix(nstop,nstart)
85 C clear mothers of the previous event
86  DO ii=1,20
87  imother(ii)=0
88  ENDDO
89 
90  DO ii=1,2
91  ison(ii)=0
92  ENDDO
93 C ... and to find mothers giving taus.
94  ndec = 0
95 C(BPK)--> LOOK FOR MOTHER, CHECK THAT IT IS NOT THE HISTORY ENTRY (E.G. MSTP(128)=0)
96  DO i=nstart,nhep
97  IF(abs(idhep(i)).EQ.kftau.AND.isthep(i).EQ.1.AND.
98  $ (isthep(i).GE.125.OR.isthep(i).LT.120)) THEN
99  imoth=jmohep(1,i)
100  DO WHILE (abs(idhep(imoth)).EQ.kftau) ! KEEP WALKING UP
101  imoth=jmohep(1,imoth)
102  ENDDO
103  IF (isthep(imoth).EQ.3.OR.
104  $ (isthep(imoth).GE.120.AND.isthep(imoth).LE.125)) THEN
105  DO j=nstart,nhep ! WE HAVE WALKED INTO HARD RECORD
106  IF (idhep(j).EQ.idhep(imoth).AND.
107  $ jmohep(1,j).EQ.imoth.AND.
108  $ isthep(j).EQ.2) THEN
109  jmoth=j
110  GOTO 66
111  ENDIF
112  ENDDO
113  ELSE
114  jmoth=imoth
115  ENDIF
116  66 CONTINUE
117  DO ii=1,ndec
118  IF(jmoth.EQ.imother(ii)) GOTO 9999
119  ENDDO
120 C(BPK)--<
121  ndec=ndec+1
122  imother(ndec)= jmoth
123  ENDIF
124  9999 CONTINUE
125  ENDDO
126 
127 C ... taus of every mother are treated in this main loop
128  DO ii=1,ndec
129  im=imother(ii)
130  ncount=0
131  np1=0
132  np2=0
133 
134 
135 C(BPK)-->
136 C CORRECTING HEPEVT IS OUT OF QUESTION AT THIS POINT..
137  im0=im
138  IF (idhep(jmohep(1,im0)).EQ.idhep(im0)) im0=jmohep(1,im0)
139  isel=-1
140  DO i=nstart,nhep
141  IF (isthep(i).EQ.3.OR.
142  $ (isthep(i).GE.120.AND.isthep(i).LE.125)) THEN ! HARD RECORD
143  GOTO 76
144  ENDIF
145  imoth=jmohep(1,i)
146  DO WHILE (idhep(imoth).EQ.idhep(i).OR.abs(idhep(imoth)).EQ.kftau) ! KEEP WALKING UP
147  imoth=jmohep(1,imoth)
148  ENDDO
149  IF ((imoth.EQ.im0.OR.imoth.EQ.im).AND.isel.EQ.-1) THEN
150  ison(1)=i
151  isel=0
152  ELSEIF ((imoth.EQ.im0.OR.imoth.EQ.im).AND.isel.EQ.0) THEN
153  ison(2)=i
154  ELSEIF ((imoth.NE.im0.AND.imoth.NE.im).AND.isel.EQ.0) THEN
155  isel=1
156  GOTO 77
157  ENDIF
158  76 CONTINUE
159  ENDDO
160  77 CONTINUE
161 C(BPK)--<
162 
163 
164 C ... we correct HEPEVT (fix developped with Catherine BISCARAT)
165 c IF (JDAHEP(2,IM).EQ.0) THEN ! ID of second daughter was missing
166 c ISECU=1
167 c DO I=JDAHEP(1,IM)+1,NHEP ! OK lets look for it
168 c IF (JMOHEP(1,I).EQ.IM.AND.ISECU.EQ.1) THEN ! we have found one
169 c JDAHEP(2,IM)=I
170 c ELSEIF (JMOHEP(1,I).EQ.IM.AND.ISECU.NE.1) THEN ! we have found one after there
171 c JDAHEP(2,IM)=0 ! was something else, lets kill game
172 c ENDIF
173 c IF (JMOHEP(1,I).NE.IM) ISECU=0 ! other stuff starts
174 c ENDDO
175 c ENDIF
176 
177 C ... we check whether there are just two or more tau-likes
178  DO i=ison(1),ison(2)
179  IF(idhep(i).EQ.-kftau.OR.idhep(i).EQ.-kfnue) ncount=ncount+1
180  IF(idhep(i).EQ. kftau.OR.idhep(i).EQ. kfnue) ncount=ncount+1
181  ENDDO
182 
183 C ... if there will be more we will come here again
184  666 CONTINUE
185 
186 C(BPK)-->
187  DO i=max(np1+1,ison(1)),ison(2)
188 C(BPK)--<
189  IF(idhep(i).EQ.-kftau.OR.idhep(i).EQ.-kfnue) np1=i
190  ENDDO
191 C(BPK)-->
192  DO i=max(np2+1,ison(1)),ison(2)
193 C(BPK)--<
194  IF(idhep(i).EQ. kftau.OR.idhep(i).EQ. kfnue) np2=i
195  ENDDO
196  DO i=1,4
197  p1(i)= phep(i,np1) !momentum of tau+
198  p2(i)= phep(i,np2) !momentum of tau-
199  q1(i)= p1(i)+p2(i)
200  ENDDO
201 
202  pol1(3)= 0d0
203  pol2(3)= 0d0
204 
205  IF(keypol.EQ.1) THEN
206 c.....include polarisation effect
207  CALL ranmar(rrr,1)
208 
209  IF(idhep(im).EQ.kfhiggs(1).OR.idhep(im).EQ.kfhiggs(2).OR.
210  $ idhep(im).EQ.kfhiggs(3)) THEN ! case of Higgs
211  IF(rrr(1).LT.0.5) THEN
212  pol1(3)= pol
213  pol2(3)=-pol
214  ELSE
215  pol1(3)=-pol
216  pol2(3)= pol
217  ENDIF
218  ELSEIF((idhep(im).EQ.kfz0).OR.(idhep(im).EQ.kfgam)) THEN ! case of gamma/Z
219 C there is no angular dependence in gamma/Z polarization
220 C there is no s-dependence in gamma/Z polarization at all
221 C there is even no Z polarization in any form
222 C main reason is that nobody asked ...
223 C but it is prepared and longitudinal correlations
224 C can be included up to KORALZ standards
225 
226  polz0=plzapx(.true.,im,np1,np2)
227  IF(rrr(1).LT.polz0) THEN
228  pol1(3)= pol
229  pol2(3)= pol
230  ELSE
231  pol1(3)=-pol
232  pol2(3)=-pol
233  ENDIF
234  ELSEIF(idhep(np1).EQ.-idhep(np2))THEN ! undef orig. only s-dep poss.
235  polz0=plzapx(.true.,im,np1,np2)
236  IF(rrr(1).LT.polz0) THEN
237  pol1(3)= pol
238  pol2(3)= pol
239  ELSE
240  pol1(3)=-pol
241  pol2(3)=-pol
242  ENDIF
243  ELSEIF(abs(idhep(im)).EQ.kfhigch) THEN ! case of charged Higgs
244  pol1(3)= pol
245  pol2(3)= pol
246  ELSE ! case of W+ or W-
247  pol1(3)= -pol
248  pol2(3)= -pol
249  ENDIF
250 c.....include polarisation effect
251  ENDIF
252 
253  IF(idhep(im).EQ.kfhiggs(1).OR.idhep(im).EQ.kfhiggs(2).OR.
254  $ idhep(im).EQ.kfhiggs(3)) THEN
255  IF(idhep(np1).EQ.-kftau .AND.
256  $ (jdahep(1,np1).LE.np1.OR.jdahep(1,np1).GT.nhep) .AND.
257  $ idhep(np2).EQ. kftau .AND.
258  $ (jdahep(1,np2).LE.np2.OR.jdahep(1,np2).GT.nhep)
259  $ ) THEN
260  IF (idhep(im).EQ.kfhiggs(1)) THEN
261  ifpseudo= .false.
262  ELSEIF (idhep(im).EQ.kfhiggs(2)) THEN
263  ifpseudo= .false.
264  ELSEIF (idhep(im).EQ.kfhiggs(3)) THEN
265  ifpseudo= .true.
266  ELSE
267  WRITE(*,*) 'Warning from TAUOLA:'
268  WRITE(*,*) 'I stop this run, wrong IDHEP(IM)=',idhep(im)
269  stop
270  ENDIF
271  CALL spinhiggs(im,np1,np2,ifpseudo,pol1,pol2)
272  IF (ifphot.EQ.1) CALL photos(im) ! Bremsstrahlung in Higgs decay
273  ! AFTER adding taus !!
274 
275 
276  ENDIF
277  ELSE
278  IF(idhep(np1).EQ.-kftau.AND.
279  $ (jdahep(1,np1).LE.np1.OR.jdahep(1,np1).GT.nhep)) THEN
280 C here check on if NP1 was not decayed should be verified
281  CALL dexay(1,pol1)
282  IF (ifphot.EQ.1) CALL photos(np1)
283  CALL taupi0(0,1,ion)
284  ENDIF
285 
286  IF(idhep(np2).EQ. kftau.AND.
287  $ (jdahep(1,np2).LE.np2.OR.jdahep(1,np2).GT.nhep)) THEN
288 C here check on if NP2 was not decayed should be added
289  CALL dexay(2,pol2)
290  IF (ifphot.EQ.1) CALL photos(np2)
291  CALL taupi0(0,2,ion)
292  ENDIF
293  ENDIF
294  ncount=ncount-2
295  IF (ncount.GT.0) GOTO 666
296  ENDDO
297 
298  ELSEIF(mode.EQ.1) THEN
299 C ***********************
300 C
301  CALL dexay(100,pol1)
302  CALL dekay(100,pol1x)
303  WRITE(iout,7002)
304  ENDIF
305 C *****
306  7001 FORMAT(///1x,15(5h*****)
307  $ /,' *', 25x,'*****TAUOLA UNIVERSAL INTERFACE: ******',9x,1h*,
308  $ /,' *', 25x,'*****VERSION 1.22, April 2009 (gfort)**',9x,1h*,
309  $ /,' *', 25x,'**AUTHORS: P. Golonka, B. Kersevan, ***',9x,1h*,
310  $ /,' *', 25x,'**T. Pierzchala, E. Richter-Was, ******',9x,1h*,
311  $ /,' *', 25x,'****** Z. Was, M. Worek ***************',9x,1h*,
312  $ /,' *', 25x,'**USEFUL DISCUSSIONS, IN PARTICULAR ***',9x,1h*,
313  $ /,' *', 25x,'*WITH C. Biscarat and S. Slabospitzky**',9x,1h*,
314  $ /,' *', 25x,'****** are warmly acknowledged ********',9x,1h*,
315  $ /,' *', 25x,' ',9x,1h*,
316  $ /,' *', 25x,'********** INITIALIZATION ************',9x,1h*,
317  $ /,' *',f20.5,5x,'tau polarization switch must be 1 or 0 ',9x,1h*,
318  $ /,' *',f20.5,5x,'Higs scalar/pseudo mix CERN-TH/2003-166',9x,1h*,
319  $ /,' *',i10, 15x,'PI0 decay switch must be 1 or 0 ',9x,1h*,
320  $ /,' *',i10, 15x,'ETA decay switch must be 1 or 0 ',9x,1h*,
321  $ /,' *',i10, 15x,'K0S decay switch must be 1 or 0 ',9x,1h*,
322  $ /,1x,15(5h*****)/)
323 
324  7002 FORMAT(///1x,15(5h*****)
325  $ /,' *', 25x,'*****TAUOLA UNIVERSAL INTERFACE: ******',9x,1h*,
326  $ /,' *', 25x,'*****VERSION 1.22, April 2009 (gfort)**',9x,1h*,
327  $ /,' *', 25x,'**AUTHORS: P. Golonka, B. Kersevan, ***',9x,1h*,
328  $ /,' *', 25x,'**T. Pierzchala, E. Richter-Was, ******',9x,1h*,
329  $ /,' *', 25x,'****** Z. Was, M. Worek ***************',9x,1h*,
330  $ /,' *', 25x,'**USEFUL DISCUSSIONS, IN PARTICULAR ***',9x,1h*,
331  $ /,' *', 25x,'*WITH C. Biscarat and S. Slabospitzky**',9x,1h*,
332  $ /,' *', 25x,'****** are warmly acknowledged ********',9x,1h*,
333  $ /,' *', 25x,'****** END OF MODULE OPERATION ********',9x,1h*,
334  $ /,1x,15(5h*****)/)
335 
336  END
337 
338  SUBROUTINE spinhiggs(IM,NP1,NP2,IFPSEUDO,Pol1,Pol2)
339  LOGICAL IFPSEUDO
340  real*8 hh1,hh2,wthiggs
341  dimension pol1(4), pol2(4),hh1(4),hh2(4), rrr(1)
342 ! CALL DEXAY(1,POL1) ! Kept for tests
343 ! CALL DEXAY(2,POL2) ! Kept for tests
344  INTEGER ION(3)
345  10 CONTINUE
346  CALL ranmar(rrr,1)
347  CALL dekay(1,hh1)
348  CALL dekay(2,hh2)
349  wt=wthiggs(ifpseudo,hh1,hh2)
350  IF (rrr(1).GT.wt) GOTO 10
351  CALL dekay(1+10,hh1)
352  CALL taupi0(0,1,ion)
353  CALL dekay(2+10,hh2)
354  CALL taupi0(0,2,ion)
355  END
356  FUNCTION wthiggs(IFPSEUDO,HH1,HH2)
357  LOGICAL IFPSEUDO
358  common /pseudocoup/ csc,ssc
359  real*4 csc,ssc
360  save pseudocoup
361  real*8 hh1(4),hh2(4),r(4,4),wthiggs
362  DO k=1,4
363  DO l=1,4
364  r(k,l)=0
365  ENDDO
366  ENDDO
367  wthiggs=0d0
368 
369  r(4,4)= 1d0 ! unpolarized part
370  r(3,3)=-1d0 ! longitudinal
371  ! other missing
372  IF (ifpseudo) THEN
373  r(1,1)=-1
374  r(2,2)= -1
375  r(1,1)=(csc**2-ssc**2)/(csc**2+ssc**2)
376  r(2,2)=(csc**2-ssc**2)/(csc**2+ssc**2)
377  r(1,2)=2*csc*ssc/(csc**2+ssc**2)
378  r(2,1)=-2*csc*ssc/(csc**2+ssc**2)
379  ELSE
380  r(1,1)=1
381  r(2,2)=1
382  ENDIF
383 
384 
385 
386  DO k=1,4
387  DO l=1,4
388  wthiggs=wthiggs+r(k,l)*hh1(k)*hh2(l)
389  ENDDO
390  ENDDO
391  wthiggs=wthiggs/4d0
392  END
393 
394  FUNCTION plzapx(HOPEin,IM0,NP1,NP2)
395 C IM0 NP1 NP2 are the positions of Z/gamma tau tau in hepevt common block.
396 C the purpose of this routine is to calculate polarization of Z along tau direction.
397 C this is highly non-trivial due to necessity of reading infromation from hard process
398 C history in HEPEVT, which is often written not up to the gramatic rules.
399  real*8 plzap0,svar,costhe,sini,sfin,zprop2,
400  $ p1(4),p2(4),q1(4),q2(4),qq(4),ph(4),pd1(4),pd2(4),
401  $ pq1(4),pq2(4),pb(4),pa(4)
402  real*4 plzapx
403  INTEGER IM
404  LOGICAL HOPE,HOPEin
405 #include "../../include/HEPEVT.h"
406 
407 C(BPK)--> BROTHERS OF TAU ALREADY FOUND
408  INTEGER ISON
409  COMMON /isons_tau/ison(2)
410 C(BPK)--<
411 C >>
412 C >> STEP 1: find where are particles in hepevent and pick them up
413 C >>
414  hope=hopein
415 C sometimes shade Z of Z is its mother ...
416  im=im0
417  im00=jmohep(1,im0)
418 C to protect against check on mother of beam particles.
419  IF (im00.GT.0) THEN
420  IF (idhep(im0).EQ.idhep(im00)) im=jmohep(1,im0)
421  ENDIF
422 C
423 C find (host generator-level) incoming beam-bare-particles which form Z and co.
424  imo1=jmohep(1,im)
425  imo2=jmohep(2,im)
426 
427 C(BPK)--> IN HERWIG THE POINTER MIGHT BE TO HARD CMS
428  im00=imo1
429  IF (isthep(im00).EQ.120) THEN
430  imo1=jmohep(1,im00)
431  imo2=jmohep(2,im00)
432  ENDIF
433 C(BPK)--<
434 
435  iffull=0
436 C case when it was like e+e- --> tau+tau- gammas and e+e- were 'first' in hepevt.
437  IF (imo1.EQ.0.AND.imo2.EQ.0) THEN
438  imo1=jmohep(1,np1)
439 C(BPK)-->
440  IF (idhep(imo1).EQ.idhep(np1)) imo1=jmohep(1,imo1) ! PROTECT AGAINST COPIES
441 C(BPK)--<
442  imo2=jmohep(2,np1)
443 C(BPK)-->
444  IF (idhep(imo2).EQ.idhep(np2)) imo2=jmohep(1,imo2) ! PROTECT AGAINST COPIES
445 C(BPK)--<
446  iffull=1
447 C case when it was like qq --> tau+tau- gammas and qq were NOT 'first' in hepevt.
448 
449  ELSEIF (idhep(im).NE.22.AND.idhep(im).NE.23) THEN
450  imo1=jmohep(1,np1)
451 C(BPK)-->
452  IF (idhep(imo1).EQ.idhep(np1)) imo1=jmohep(1,imo1) ! PROTECT AGAINST COPIES
453 C(BPK)--<
454  imo2=jmohep(2,np1)
455 C(BPK)-->
456  IF (idhep(imo2).EQ.idhep(np2)) imo2=jmohep(1,imo2) ! PROTECT AGAINST COPIES
457 C(BPK)--<
458  iffull=1
459  ENDIF
460 
461 
462 C and check if it really happened
463  IF (imo1.EQ.0) hope=.false.
464  IF (imo2.EQ.0) hope=.false.
465  IF (imo1.EQ.imo2) hope=.false.
466 
467 C
468  DO i=1,4
469  q1(i)= phep(i,np1) !momentum of tau+
470  q2(i)= phep(i,np2) !momentum of tau-
471  ENDDO
472 
473 C corrections due to possible differences in 4-momentum of shadow vs true Z.
474 C(BPK)-->
475  IF (im.EQ.jmohep(1,im0).AND.
476  $ (idhep(im).EQ.22.OR.idhep(im).EQ.23)) THEN
477  DO k=1,4
478  pb(k)=phep(k,im)
479  pa(k)=phep(k,im0)
480  ENDDO
481 C(BPK)--<
482 
483  CALL bostdq( 1,pa, q1, q1)
484  CALL bostdq( 1,pa, q2, q2)
485  CALL bostdq(-1,pb, q1, q1)
486  CALL bostdq(-1,pb, q2, q2)
487 
488  ENDIF
489 
490  DO i=1,4
491  qq(i)= q1(i)+q2(i) !momentum of Z
492  IF (hope) p1(i)=phep(i,imo1) !momentum of beam1
493  IF (hope) p2(i)=phep(i,imo2) !momentum of beam2
494  ph(i)=0d0
495  pd1(i)=0d0
496  pd2(i)=0d0
497  ENDDO
498 ! These momenta correspond to quarks, gluons photons or taus
499  idfq1=idhep(np1)
500  idfq2=idhep(np2)
501  IF (hope) idfp1=idhep(imo1)
502  IF (hope) idfp2=idhep(imo2)
503 
504  svar=qq(4)**2-qq(3)**2-qq(2)**2-qq(1)**2
505  IF (.NOT.hope) THEN
506 C options which may be useful in some cases of two heavy boson production
507 C need individual considerations. To be developed.
508 C PLZAPX=PLZAP0(11,IDFQ1,SVAR,0D0) ! gamma/Z mixture as if produced from e beam
509 C PLZAPX=PLZAP0(12,IDFQ1,SVAR,0D0) ! pure Z
510  plzapx=0.5 ! pure gamma
511  RETURN
512  ENDIF
513 C >>
514 C >> STEP 2 look for brothers of Z which have to be included in effective incoming particles
515 C >>
516 C let us define beginning and end of particles which are produced in parallel to Z
517 C in principle following should work
518 
519 C(BPK)--> ACCOMMODATE FOR HERWIG - IM00 POINTS TO BEAM PARTICLE OR HARD CMS
520  nx1=jdahep(1,im00)
521  nx2=jdahep(2,im00)
522 C but ...
523  inbr=im ! OK, HARD RECORD Z/GAMMA
524  IF (iffull.EQ.1) inbr=np1 ! OK, NO Z/GAMMA
525  IF (idhep(jmohep(1,inbr)).EQ.idhep(inbr)) inbr=jmohep(1,inbr) ! FORCE HARD RECORD
526 C(BPK)--<
527  IF(nx1.EQ.0.OR.nx2.EQ.0) THEN
528  nx1=inbr
529  nx2=inbr
530  DO k=1,inbr-1
531  IF(jmohep(1,inbr-k).EQ.jmohep(1,inbr)) THEN
532  nx1=inbr-k
533  ELSE
534  GOTO 7
535  ENDIF
536  ENDDO
537  7 CONTINUE
538 
539  DO k=inbr+1,nhep
540  IF(jmohep(1,k).EQ.jmohep(1,inbr)) THEN
541  nx2=k
542  ELSE
543  GOTO 8
544  ENDIF
545  ENDDO
546  8 CONTINUE
547  ENDIF
548 
549 C case of annihilation of two bosons is hopeles
550  IF (abs(idfp1).GE.20.AND.abs(idfp2).GE.20) hope=.false.
551 C case of annihilation of non-matching flavors is hopeless
552  IF (abs(idfp1).LE.20.AND.abs(idfp2).LE.20.AND.idfp1+idfp2.NE.0)
553  $ hope=.false.
554  IF (.NOT.hope) THEN
555 C options which may be useful in some cases of two heavy boson production
556 C need individual considerations. To be developed.
557 C PLZAPX=PLZAP0(11,IDFQ1,SVAR,0D0) ! gamma/Z mixture as if produced from e beam
558 C PLZAPX=PLZAP0(12,IDFQ1,SVAR,0D0) ! pure Z
559  plzapx=0.5 ! pure gamma
560  RETURN
561  ENDIF
562  IF (abs(idfp1).LT.20) ide= idfp1
563  IF (abs(idfp2).LT.20) ide=-idfp2
564 
565 
566 C >>
567 C >> STEP 3 we combine gluons, photons into incoming effective beams
568 C >>
569 
570 C in the following we will ignore the possibility of photon emission from taus
571 C however at certain moment it will be necessary to take care of
572 
573  DO l=1,4
574  pd1(l)=p1(l)
575  pd2(l)=p2(l)
576  ENDDO
577 
578  DO l=1,4
579  pq1(l)=q1(l)
580  pq2(l)=q2(l)
581  ENDDO
582 
583  iflav=min(abs(idfp1),abs(idfp2))
584 
585 *--------------------------------------------------------------------------
586 c IFLAV=min(ABS(IDFP1),ABS(IDFP2))
587 c that means that always quark or lepton i.e. process like
588 c f g(gamma) --> f Z0 --> tau tau
589 c we glue fermions to effective beams that is f f --> Z0 --> tau tau
590 c with gamma/g emission from initial fermion.
591 *---------------------------------------------------------------------------
592 
593  IF (abs(idfp1).GE.20) THEN
594  DO k=nx1,nx2
595  idp=idhep(k)
596  IF (abs(idp).EQ.iflav) THEN
597  DO l=1,4
598  pd1(l)=-phep(l,k)
599  ENDDO
600  ENDIF
601  ENDDO
602  ENDIF
603 
604  IF (abs(idfp2).GE.20) THEN
605  DO k=nx1,nx2
606  idp=idhep(k)
607  IF (abs(idp).EQ.iflav) THEN
608  DO l=1,4
609  pd2(l)=-phep(l,k)
610  ENDDO
611  ENDIF
612  ENDDO
613  ENDIF
614 
615 C if first beam was boson: gluing
616 
617  IF (abs(idfp1).GE.20) THEN
618  DO l=1,4
619  ph(l)=p1(l)
620  ENDDO
621  xm1=abs((pd1(4)+ph(4))**2-(pd1(3)+ph(3))**2
622  $ -(pd1(2)+ph(2))**2-(pd1(1)+ph(1))**2)
623  xm2=abs((pd2(4)+ph(4))**2-(pd2(3)+ph(3))**2
624  $ -(pd2(2)+ph(2))**2-(pd2(1)+ph(1))**2)
625  IF (xm1.LT.xm2) THEN
626  DO l=1,4
627  pd1(l)=pd1(l)+p1(l)
628  ENDDO
629  ELSE
630  DO l=1,4
631  pd2(l)=pd2(l)+p1(l)
632  ENDDO
633  ENDIF
634  ENDIF
635 
636 C if second beam was boson: gluing
637 
638 
639  IF (abs(idfp2).GE.20) THEN
640  DO l=1,4
641  ph(l)=p2(l)
642  ENDDO
643  xm1=abs((pd1(4)+ph(4))**2-(pd1(3)+ph(3))**2
644  $ -(pd1(2)+ph(2))**2-(pd1(1)+ph(1))**2)
645  xm2=abs((pd2(4)+ph(4))**2-(pd2(3)+ph(3))**2
646  $ -(pd2(2)+ph(2))**2-(pd2(1)+ph(1))**2)
647  IF (xm1.LT.xm2) THEN
648  DO l=1,4
649  pd1(l)=pd1(l)+p2(l)
650  ENDDO
651  ELSE
652  DO l=1,4
653  pd2(l)=pd2(l)+p2(l)
654  ENDDO
655  ENDIF
656  ENDIF
657 
658 C now spectators ...
659 
660 C(BPK)-->
661  nph1=np1
662  nph2=np2
663  IF (idhep(jmohep(1,np1)).EQ.idhep(np1)) nph1=jmohep(1,np1) ! SHOULD PUT US IN HARD REC.
664  IF (idhep(jmohep(1,np2)).EQ.idhep(np2)) nph2=jmohep(1,np2) ! SHOULD PUT US IN HARD REC.
665 C(BPK)--<
666 
667  DO k=nx1,nx2
668  IF (abs(idhep(k)).NE.iflav.AND.k.NE.im.AND.
669 C(BPK)-->
670  $ k.NE.nph1.AND.k.NE.nph2) THEN
671 C(BPK)--<
672  IF(idhep(k).EQ.22.AND.iffull.EQ.1) THEN
673  DO l=1,4
674  ph(l)=phep(l,k)
675  ENDDO
676  xm1=abs((pd1(4)-ph(4))**2-(pd1(3)-ph(3))**2
677  $ -(pd1(2)-ph(2))**2-(pd1(1)-ph(1))**2)
678  xm2=abs((pd2(4)-ph(4))**2-(pd2(3)-ph(3))**2
679  $ -(pd2(2)-ph(2))**2-(pd2(1)-ph(1))**2)
680  xm3=abs((pq1(4)+ph(4))**2-(pq1(3)+ph(3))**2
681  $ -(pq1(2)+ph(2))**2-(pq1(1)+ph(1))**2)
682  xm4=abs((pq2(4)+ph(4))**2-(pq2(3)+ph(3))**2
683  $ -(pq2(2)+ph(2))**2-(pq2(1)+ph(1))**2)
684 
685 
686  sini=abs((pd1(4)+pd2(4)-ph(4))**2-(pd1(3)+pd2(3)-ph(3))**2
687  $ -(pd1(2)+pd2(2)-ph(2))**2-(pd1(1)+pd2(1)-ph(1))**2)
688  sfin=abs((pd1(4)+pd2(4) )**2-(pd1(3)+pd2(3) )**2
689  $ -(pd1(2)+pd2(2) )**2-(pd1(1)+pd2(1) )**2)
690 
691  facini=zprop2(sini)
692  facfin=zprop2(sfin)
693 
694  xm1=xm1/facini
695  xm2=xm2/facini
696  xm3=xm3/facfin
697  xm4=xm4/facfin
698 
699  xm=min(xm1,xm2,xm3,xm4)
700  IF (xm1.EQ.xm) THEN
701  DO l=1,4
702  pd1(l)=pd1(l)-ph(l)
703  ENDDO
704  ELSEIF (xm2.EQ.xm) THEN
705  DO l=1,4
706  pd2(l)=pd2(l)-ph(l)
707  ENDDO
708  ELSEIF (xm3.EQ.xm) THEN
709  DO l=1,4
710  q1(l)=pq1(l)+ph(l)
711  ENDDO
712  ELSE
713  DO l=1,4
714  q2(l)=pq2(l)+ph(l)
715  ENDDO
716  ENDIF
717  ELSE
718  DO l=1,4
719  ph(l)=phep(l,k)
720  ENDDO
721  xm1=abs((pd1(4)-ph(4))**2-(pd1(3)-ph(3))**2
722  $ -(pd1(2)-ph(2))**2-(pd1(1)-ph(1))**2)
723  xm2=abs((pd2(4)-ph(4))**2-(pd2(3)-ph(3))**2
724  $ -(pd2(2)-ph(2))**2-(pd2(1)-ph(1))**2)
725  IF (xm1.LT.xm2) THEN
726  DO l=1,4
727  pd1(l)=pd1(l)-ph(l)
728  ENDDO
729  ELSE
730  DO l=1,4
731  pd2(l)=pd2(l)-ph(l)
732  ENDDO
733  ENDIF
734  ENDIF
735  ENDIF
736  ENDDO
737 
738 
739 C >>
740 C >> STEP 4 look for brothers of tau (sons of Z!) which have to be included in
741 c >> effective outcoming taus
742 C >>
743 C let us define beginning and end of particles which are produced in
744 c parallel to tau
745 
746 
747 
748 C find outcoming particles which come from Z
749 
750 
751 
752 
753 C(BPK)--> OK, IT WOULD HAVE TO BE ALONG TAUS IN HARD RECORD WITH THE SAME MOTHER
754  IF (abs(idhep(im0)).EQ.22.OR.abs(idhep(im0)).EQ.23) THEN
755  DO k=ison(1),ison(2)
756  IF(abs(idhep(k)).EQ.22) THEN
757 C(BPK)--<
758 
759  do l=1,4
760  ph(l)=phep(l,k)
761  enddo
762 
763  xm3=abs((pq1(4)+ph(4))**2-(pq1(3)+ph(3))**2
764  $ -(pq1(2)+ph(2))**2-(pq1(1)+ph(1))**2)
765  xm4=abs((pq2(4)+ph(4))**2-(pq2(3)+ph(3))**2
766  $ -(pq2(2)+ph(2))**2-(pq2(1)+ph(1))**2)
767 
768  xm=min(xm3,xm4)
769 
770  IF (xm3.EQ.xm) THEN
771  DO l=1,4
772  q1(l)=pq1(l)+ph(l)
773  ENDDO
774  ELSE
775  DO l=1,4
776  q2(l)=pq2(l)+ph(l)
777  ENDDO
778  ENDIF
779  endif
780  enddo
781  ENDIF
782 
783 
784 
785 *------------------------------------------------------------------------
786 
787 
788 C out of effective momenta we calculate COSTHE and later polarization
789  CALL angulu(pd1,pd2,q1,q2,costhe)
790 
791  plzapx=plzap0(ide,idfq1,svar,costhe)
792  END
793 
794  SUBROUTINE angulu(PD1,PD2,Q1,Q2,COSTHE)
795  real*8 pd1(4),pd2(4),q1(4),q2(4),costhe,p(4),qq(4),qt(4)
796 C take effective beam which is less massive, it should be irrelevant
797 C but in case HEPEVT is particulary dirty may help.
798 C this routine calculate reduced system transver and cosine of scattering
799 C angle.
800 
801  xm1=abs(pd1(4)**2-pd1(3)**2-pd1(2)**2-pd1(1)**2)
802  xm2=abs(pd2(4)**2-pd2(3)**2-pd2(2)**2-pd2(1)**2)
803  IF (xm1.LT.xm2) THEN
804  sign=1d0
805  DO k=1,4
806  p(k)=pd1(k)
807  ENDDO
808  ELSE
809  sign=-1d0
810  DO k=1,4
811  p(k)=pd2(k)
812  ENDDO
813  ENDIF
814 C calculate space like part of P (in Z restframe)
815  DO k=1,4
816  qq(k)=q1(k)+q2(k)
817  qt(k)=q1(k)-q2(k)
818  ENDDO
819 
820  xmqq=sqrt(qq(4)**2-qq(3)**2-qq(2)**2-qq(1)**2)
821 
822  qtxqq=qt(4)*qq(4)-qt(3)*qq(3)-qt(2)*qq(2)-qt(1)*qq(1)
823  DO k=1,4
824  qt(k)=qt(k)-qq(k)*qtxqq/xmqq**2
825  ENDDO
826 
827  pxqq=p(4)*qq(4)-p(3)*qq(3)-p(2)*qq(2)-p(1)*qq(1)
828  DO k=1,4
829  p(k)=p(k)-qq(k)*pxqq/xmqq**2
830  ENDDO
831 C calculate costhe
832  pxp =sqrt(p(1)**2+p(2)**2+p(3)**2-p(4)**2)
833  qtxqt=sqrt(qt(3)**2+qt(2)**2+qt(1)**2-qt(4)**2)
834  pxqt =p(3)*qt(3)+p(2)*qt(2)+p(1)*qt(1)-p(4)*qt(4)
835  costhe=pxqt/pxp/qtxqt
836  costhe=costhe*sign
837  END
838 
839  FUNCTION plzap0(IDE,IDF,SVAR,COSTH0)
840 C this function calculates probability for the helicity +1 +1 configuration
841 C of taus for given Z/gamma transfer and COSTH0 cosine of scattering angle
842  real*8 plzap0,svar,costhe,costh0,t_born
843 
844  costhe=costh0
845 C >>>>> IF (IDE*IDF.LT.0) COSTHE=-COSTH0 ! this is probably not needed ID
846 C >>>>> of first beam is used by T_GIVIZ0 including sign
847 
848  IF (idf.GT.0) THEN
849  CALL initwk(ide,idf,svar)
850  ELSE
851  CALL initwk(-ide,-idf,svar)
852  ENDIF
853  plzap0=t_born(0,svar,costhe,1d0,1d0)
854  $ /(t_born(0,svar,costhe,1d0,1d0)+t_born(0,svar,costhe,-1d0,-1d0))
855 
856 ! PLZAP0=0.5
857  END
858  FUNCTION t_born(MODE,SVAR,COSTHE,TA,TB)
859 C ----------------------------------------------------------------------
860 C THIS ROUTINE PROVIDES BORN CROSS SECTION. IT HAS THE SAME
861 C STRUCTURE AS FUNTIS AND FUNTIH, THUS CAN BE USED AS SIMPLER
862 C EXAMPLE OF THE METHOD APPLIED THERE
863 C INPUT PARAMETERS ARE: SVAR -- transfer
864 C COSTHE -- cosine of angle between tau+ and 1st beam
865 C TA,TB -- helicity states of tau+ tau-
866 C
867 C called by : BORNY, BORAS, BORNV, WAGA, WEIGHT
868 C ----------------------------------------------------------------------
869  IMPLICIT REAL*8(a-h,o-z)
870  COMMON / t_beampm / ene ,amin,amfin,ide,idf
871  real*8 ene ,amin,amfin
872  COMMON / t_gauspm /ss,poln,t3e,qe,t3f,qf
873  & ,xupgi ,xupzi ,xupgf ,xupzf
874  & ,ndiag0,ndiaga,keya,keyz
875  & ,itce,jtce,itcf,jtcf,kolor
876  real*8 ss,poln,t3e,qe,t3f,qf
877  & ,xupgi(2),xupzi(2),xupgf(2),xupzf(2)
878  real*8 seps1,seps2
879 C=====================================================================
880  COMMON / t_gswprm /swsq,amw,amz,amh,amtop,gammz
881  real*8 swsq,amw,amz,amh,amtop,gammz
882 C SWSQ = sin2 (theta Weinberg)
883 C AMW,AMZ = W & Z boson masses respectively
884 C AMH = the Higgs mass
885 C AMTOP = the top mass
886 C GAMMZ = Z0 width
887  COMPLEX*16 ABORN(2,2),APHOT(2,2),AZETT(2,2)
888  COMPLEX*16 XUPZFP(2),XUPZIP(2)
889  COMPLEX*16 ABORNM(2,2),APHOTM(2,2),AZETTM(2,2)
890  COMPLEX*16 PROPA,PROPZ
891  COMPLEX*16 XR,XI
892  COMPLEX*16 XUPF,XUPI,XFF(4),XFEM,XFOTA,XRHO,XKE,XKF,XKEF
893  COMPLEX*16 XTHING,XVE,XVF,XVEF
894  DATA xi/(0.d0,1.d0)/,xr/(1.d0,0.d0)/
895  DATA mode0 /-5/
896  DATA ide0 /-55/
897  DATA svar0,cost0 /-5.d0,-6.d0/
898  DATA pi /3.141592653589793238462643d0/
899  DATA seps1,seps2 /0d0,0d0/
900 C
901 C MEMORIZATION =========================================================
902  IF ( mode.NE.mode0.OR.svar.NE.svar0.OR.costhe.NE.cost0
903  $ .OR.ide0.NE.ide)THEN
904 C
905  keygsw=1
906 C ** PROPAGATORS
907  ide0=ide
908  mode0=mode
909  svar0=svar
910  cost0=costhe
911  sinthe=sqrt(1.d0-costhe**2)
912  beta=sqrt(max(0d0,1d0-4d0*amfin**2/svar))
913 C I MULTIPLY AXIAL COUPLING BY BETA FACTOR.
914  xupzfp(1)=0.5d0*(xupzf(1)+xupzf(2))+0.5*beta*(xupzf(1)-xupzf(2))
915  xupzfp(2)=0.5d0*(xupzf(1)+xupzf(2))-0.5*beta*(xupzf(1)-xupzf(2))
916  xupzip(1)=0.5d0*(xupzi(1)+xupzi(2))+0.5*(xupzi(1)-xupzi(2))
917  xupzip(2)=0.5d0*(xupzi(1)+xupzi(2))-0.5*(xupzi(1)-xupzi(2))
918 C FINAL STATE VECTOR COUPLING
919  xupf =0.5d0*(xupzf(1)+xupzf(2))
920  xupi =0.5d0*(xupzi(1)+xupzi(2))
921  xthing =0d0
922 
923  propa =1d0/svar
924  propz =1d0/dcmplx(svar-amz**2,svar/amz*gammz)
925  IF (keygsw.EQ.0) propz=0.d0
926  DO 50 i=1,2
927  DO 50 j=1,2
928  regula= (3-2*i)*(3-2*j) + costhe
929  regulm=-(3-2*i)*(3-2*j) * sinthe *2.d0*amfin/sqrt(svar)
930  aphot(i,j)=propa*(xupgi(i)*xupgf(j)*regula)
931  azett(i,j)=propz*(xupzip(i)*xupzfp(j)+xthing)*regula
932  aborn(i,j)=aphot(i,j)+azett(i,j)
933  aphotm(i,j)=propa*dcmplx(0d0,1d0)*xupgi(i)*xupgf(j)*regulm
934  azettm(i,j)=propz*dcmplx(0d0,1d0)*(xupzip(i)*xupf+xthing)*regulm
935  abornm(i,j)=aphotm(i,j)+azettm(i,j)
936  50 CONTINUE
937  ENDIF
938 C
939 C******************
940 C* IN CALCULATING CROSS SECTION ONLY DIAGONAL ELEMENTS
941 C* OF THE SPIN DENSITY MATRICES ENTER (LONGITUD. POL. ONLY.)
942 C* HELICITY CONSERVATION EXPLICITLY OBEYED
943  polar1= (seps1)
944  polar2= (-seps2)
945  born=0d0
946  DO 150 i=1,2
947  helic= 3-2*i
948  DO 150 j=1,2
949  helit=3-2*j
950  factor=kolor*(1d0+helic*polar1)*(1d0-helic*polar2)/4d0
951  factom=factor*(1+helit*ta)*(1-helit*tb)
952  factor=factor*(1+helit*ta)*(1+helit*tb)
953 
954  born=born+cdabs(aborn(i,j))**2*factor
955 C MASS TERM IN BORN
956  IF (mode.GE.1) THEN
957  born=born+cdabs(abornm(i,j))**2*factom
958  ENDIF
959 
960  150 CONTINUE
961 C************
962  funt=born
963  IF(funt.LT.0.d0) funt=born
964 
965 C
966  IF (svar.GT.4d0*amfin**2) THEN
967 C PHASE SPACE THRESHOLD FACTOR
968  thresh=sqrt(1-4d0*amfin**2/svar)
969  t_born= funt*svar**2*thresh
970  ELSE
971  thresh=0.d0
972  t_born=0.d0
973  ENDIF
974 C ZW HERE WAS AN ERROR 19. 05. 1989
975 ! write(*,*) 'KKKK ',PROPA,PROPZ,XUPGI,XUPGF,XUPZI,XUPZF
976 ! write(*,*) 'KKKK X',svar,costhe,TA,TB,T_BORN
977  END
978 
979  SUBROUTINE initwk(IDEX,IDFX,SVAR)
980 ! initialization routine coupling masses etc.
981  IMPLICIT REAL*8 (a-h,o-z)
982  COMMON / t_beampm / ene ,amin,amfin,ide,idf
983  real*8 ene ,amin,amfin
984  COMMON / t_gauspm /ss,poln,t3e,qe,t3f,qf
985  & ,xupgi ,xupzi ,xupgf ,xupzf
986  & ,ndiag0,ndiaga,keya,keyz
987  & ,itce,jtce,itcf,jtcf,kolor
988  real*8 ss,poln,t3e,qe,t3f,qf
989  & ,xupgi(2),xupzi(2),xupgf(2),xupzf(2)
990  COMMON / t_gswprm /swsq,amw,amz,amh,amtop,gammz
991  real*8 swsq,amw,amz,amh,amtop,gammz
992 C SWSQ = sin2 (theta Weinberg)
993 C AMW,AMZ = W & Z boson masses respectively
994 C AMH = the Higgs mass
995 C AMTOP = the top mass
996 C GAMMZ = Z0 width
997 C
998  ene=sqrt(svar)/2
999  amin=0.511d-3
1000  swsq=0.23147
1001  amz=91.1882
1002  gammz=2.4952
1003  IF (idfx.EQ. 15) then
1004  idf=2 ! denotes tau +2 tau-
1005  amfin=1.77703 !this mass is irrelevant if small, used in ME only
1006  ELSEIF (idfx.EQ.-15) then
1007  idf=-2 ! denotes tau -2 tau-
1008  amfin=1.77703 !this mass is irrelevant if small, used in ME only
1009  ELSE
1010  WRITE(*,*) 'INITWK: WRONG IDFX'
1011  stop
1012  ENDIF
1013 
1014  IF (idex.EQ. 11) then !electron
1015  ide= 2
1016  amin=0.511d-3
1017  ELSEIF (idex.EQ.-11) then !positron
1018  ide=-2
1019  amin=0.511d-3
1020  ELSEIF (idex.EQ. 13) then !mu+
1021  ide= 2
1022  amin=0.105659
1023  ELSEIF (idex.EQ.-13) then !mu-
1024  ide=-2
1025  amin=0.105659
1026  ELSEIF (idex.EQ. 1) then !d
1027  ide= 4
1028  amin=0.05
1029  ELSEIF (idex.EQ.- 1) then !d~
1030  ide=-4
1031  amin=0.05
1032  ELSEIF (idex.EQ. 2) then !u
1033  ide= 3
1034  amin=0.02
1035  ELSEIF (idex.EQ.- 2) then !u~
1036  ide=-3
1037  amin=0.02
1038  ELSEIF (idex.EQ. 3) then !s
1039  ide= 4
1040  amin=0.3
1041  ELSEIF (idex.EQ.- 3) then !s~
1042  ide=-4
1043  amin=0.3
1044  ELSEIF (idex.EQ. 4) then !c
1045  ide= 3
1046  amin=1.3
1047  ELSEIF (idex.EQ.- 4) then !c~
1048  ide=-3
1049  amin=1.3
1050  ELSEIF (idex.EQ. 5) then !b
1051  ide= 4
1052  amin=4.5
1053  ELSEIF (idex.EQ.- 5) then !b~
1054  ide=-4
1055  amin=4.5
1056  ELSEIF (idex.EQ. 12) then !nu_e
1057  ide= 1
1058  amin=0.1d-3
1059  ELSEIF (idex.EQ.- 12) then !nu_e~
1060  ide=-1
1061  amin=0.1d-3
1062  ELSEIF (idex.EQ. 14) then !nu_mu
1063  ide= 1
1064  amin=0.1d-3
1065  ELSEIF (idex.EQ.- 14) then !nu_mu~
1066  ide=-1
1067  amin=0.1d-3
1068  ELSEIF (idex.EQ. 16) then !nu_tau
1069  ide= 1
1070  amin=0.1d-3
1071  ELSEIF (idex.EQ.- 16) then !nu_tau~
1072  ide=-1
1073  amin=0.1d-3
1074 
1075  ELSE
1076  WRITE(*,*) 'INITWK: WRONG IDEX'
1077  stop
1078  ENDIF
1079 
1080 C ----------------------------------------------------------------------
1081 C
1082 C INITIALISATION OF COUPLING CONSTANTS AND FERMION-GAMMA / Z0 VERTEX
1083 C
1084 C called by : KORALZ
1085 C ----------------------------------------------------------------------
1086  itce=ide/iabs(ide)
1087  jtce=(1-itce)/2
1088  itcf=idf/iabs(idf)
1089  jtcf=(1-itcf)/2
1090  CALL t_givizo( ide, 1,aizor,qe,kdumm)
1091  CALL t_givizo( ide,-1,aizol,qe,kdumm)
1092  xupgi(1)=qe
1093  xupgi(2)=qe
1094  t3e = aizol+aizor
1095  xupzi(1)=(aizor-qe*swsq)/sqrt(swsq*(1-swsq))
1096  xupzi(2)=(aizol-qe*swsq)/sqrt(swsq*(1-swsq))
1097  CALL t_givizo( idf, 1,aizor,qf,kolor)
1098  CALL t_givizo( idf,-1,aizol,qf,kolor)
1099  xupgf(1)=qf
1100  xupgf(2)=qf
1101  t3f = aizol+aizor
1102  xupzf(1)=(aizor-qf*swsq)/sqrt(swsq*(1-swsq))
1103  xupzf(2)=(aizol-qf*swsq)/sqrt(swsq*(1-swsq))
1104 C
1105  ndiag0=2
1106  ndiaga=11
1107  keya = 1
1108  keyz = 1
1109 C
1110 C
1111  RETURN
1112  END
1113 
1114  SUBROUTINE t_givizo(IDFERM,IHELIC,SIZO3,CHARGE,KOLOR)
1115 C ----------------------------------------------------------------------
1116 C PROVIDES ELECTRIC CHARGE AND WEAK IZOSPIN OF A FAMILY FERMION
1117 C IDFERM=1,2,3,4 DENOTES NEUTRINO, LEPTON, UP AND DOWN QUARK
1118 C NEGATIVE IDFERM=-1,-2,-3,-4, DENOTES ANTIPARTICLE
1119 C IHELIC=+1,-1 DENOTES RIGHT AND LEFT HANDEDNES ( CHIRALITY)
1120 C SIZO3 IS THIRD PROJECTION OF WEAK IZOSPIN (PLUS MINUS HALF)
1121 C AND CHARGE IS ELECTRIC CHARGE IN UNITS OF ELECTRON CHARGE
1122 C KOLOR IS A QCD COLOUR, 1 FOR LEPTON, 3 FOR QUARKS
1123 C
1124 C called by : EVENTE, EVENTM, FUNTIH, .....
1125 C ----------------------------------------------------------------------
1126  IMPLICIT REAL*8(a-h,o-z)
1127 C
1128  IF(idferm.EQ.0.OR.iabs(idferm).GT.4) GOTO 901
1129  IF(iabs(ihelic).NE.1) GOTO 901
1130  ih =ihelic
1131  idtype =iabs(idferm)
1132  ic =idferm/idtype
1133  lepqua=int(idtype*0.4999999d0)
1134  iupdow=idtype-2*lepqua-1
1135  charge =(-iupdow+2d0/3d0*lepqua)*ic
1136  sizo3 =0.25d0*(ic-ih)*(1-2*iupdow)
1137  kolor=1+2*lepqua
1138 C** NOTE THAT CONVENTIONALY Z0 COUPLING IS
1139 C** XOUPZ=(SIZO3-CHARGE*SWSQ)/SQRT(SWSQ*(1-SWSQ))
1140  RETURN
1141  901 print *,' STOP IN GIVIZO: WRONG PARAMS.'
1142  stop
1143  END
1144 #include "../../include/phyfix.h"
1145  SUBROUTINE filhep(N,IST,ID,JMO1,JMO2,JDA1,JDA2,P4,PINV,PHFLAG)
1146 C ----------------------------------------------------------------------
1147 C this subroutine fills one entry into the HEPEVT common
1148 C and updates the information for affected mother entries
1149 C
1150 C written by Martin W. Gruenewald (91/01/28)
1151 C
1152 C called by : ZTOHEP,BTOHEP,DWLUxy
1153 C ----------------------------------------------------------------------
1154 C
1155 C this is the hepevt class in old style. No d_h_ class pre-name
1156 #include "../../include/HEPEVT.h"
1157  LOGICAL PHFLAG
1158 C
1159  real*4 p4(4)
1160 C
1161 C check address mode
1162  IF (n.EQ.0) THEN
1163 C
1164 C append mode
1165  ihep=nhep+1
1166  ELSE IF (n.GT.0) THEN
1167 C
1168 C absolute position
1169  ihep=n
1170  ELSE
1171 C
1172 C relative position
1173  ihep=nhep+n
1174  END IF
1175 C
1176 C check on IHEP
1177  IF ((ihep.LE.0).OR.(ihep.GT.nmxhep)) RETURN
1178 C
1179 C add entry
1180  nhep=ihep
1181  isthep(ihep)=ist
1182  idhep(ihep)=id
1183  jmohep(1,ihep)=jmo1
1184  IF(jmo1.LT.0)jmohep(1,ihep)=jmohep(1,ihep)+ihep
1185  jmohep(2,ihep)=jmo2
1186  IF(jmo2.LT.0)jmohep(2,ihep)=jmohep(2,ihep)+ihep
1187  jdahep(1,ihep)=jda1
1188  jdahep(2,ihep)=jda2
1189 C
1190  DO i=1,4
1191  phep(i,ihep)=p4(i)
1192 C
1193 C KORAL-B and KORAL-Z do not provide vertex and/or lifetime informations
1194  vhep(i,ihep)=0.0
1195  END DO
1196  phep(5,ihep)=pinv
1197 C FLAG FOR PHOTOS...
1198  qedrad(ihep)=phflag
1199 C
1200 C update process:
1201  DO ip=jmohep(1,ihep),jmohep(2,ihep)
1202  IF(ip.GT.0)THEN
1203 C
1204 C if there is a daughter at IHEP, mother entry at IP has decayed
1205  IF(isthep(ip).EQ.1)isthep(ip)=2
1206 C
1207 C and daughter pointers of mother entry must be updated
1208  IF(jdahep(1,ip).EQ.0)THEN
1209  jdahep(1,ip)=ihep
1210  jdahep(2,ip)=ihep
1211  ELSE
1212  jdahep(2,ip)=max(ihep,jdahep(2,ip))
1213  END IF
1214  END IF
1215  END DO
1216 C
1217  RETURN
1218  END
1219 
1220 
1221  FUNCTION ihepdim(DUM)
1222 C this is the hepevt class in old style. No d_h_ class pre-name
1223 #include "../../include/HEPEVT.h"
1224  ihepdim=nhep
1225  END
1226  FUNCTION zprop2(S)
1227  IMPLICIT REAL*8(a-h,o-z)
1228  COMPLEX*16 CPRZ0,CPRZ0M
1229  amz=91.1882
1230  gammz=2.49
1231  cprz0=dcmplx((s-amz**2),s/amz*gammz)
1232  cprz0m=1/cprz0
1233  zprop2=(abs(cprz0m))**2
1234  END
1235 
1236  SUBROUTINE taupi0(MODE,JAK,ION)
1237 C no initialization required. Must be called once after every:
1238 C 1) CALL DEKAY(1+10,...)
1239 C 2) CALL DEKAY(2+10,...)
1240 C 3) CALL DEXAY(1,...)
1241 C 4) CALL DEXAY(2,...)
1242 C subroutine to decay originating from TAUOLA's taus:
1243 C 1) etas (with CALL TAUETA(JAK))
1244 C 2) later pi0's from taus.
1245 C 3) extensions to other applications possible.
1246 C this routine belongs to >tauola universal interface<, but uses
1247 C routines from >tauola< utilities as well. 25.08.2005
1248 #include "../../include/HEPEVT.h"
1249 
1250 C position of taus, must be defined by host program:
1251  COMMON /taupos/ np1,np2
1252 c
1253  REAL PHOT1(4),PHOT2(4)
1254  real*8 r,x(4),y(4),pi0(4)
1255  INTEGER JEZELI(3),ION(3)
1256  DATA jezeli /0,0,0/
1257  SAVE jezeli
1258  IF (mode.EQ.-1) THEN
1259  jezeli(1)=ion(1)
1260  jezeli(2)=ion(2)
1261  jezeli(3)=ion(3)
1262  RETURN
1263  ENDIF
1264  IF (jezeli(1).EQ.0) RETURN
1265  IF (jezeli(2).EQ.1) CALL taueta(jak)
1266  IF (jezeli(3).EQ.1) CALL tauk0s(jak)
1267 C position of decaying particle:
1268  IF((kto.EQ. 1).OR.(kto.EQ.11)) THEN
1269  nps=np1
1270  ELSE
1271  nps=np2
1272  ENDIF
1273  nhepm=nhep ! to avoid infinite loop
1274  DO k=jdahep(1,nps),nhepm ! we search for pi0's from tau till eor.
1275  IF (idhep(k).EQ.111.AND.jdahep(1,k).LE.k) THEN ! IF we found pi0
1276  DO l=1,4
1277  pi0(l)= phep(l,k)
1278  ENDDO
1279 ! random 3 vector on the sphere, masless
1280  r=sqrt(pi0(4)**2-pi0(3)**2-pi0(2)**2-pi0(1)**2)/2d0
1281  CALL spherd(r,x)
1282  x(4)=r
1283  y(4)=r
1284 
1285  y(1)=-x(1)
1286  y(2)=-x(2)
1287  y(3)=-x(3)
1288 ! boost to lab and to real*4
1289  CALL bostdq(-1,pi0,x,x)
1290  CALL bostdq(-1,pi0,y,y)
1291  DO l=1,4
1292  phot1(l)=x(l)
1293  phot2(l)=y(l)
1294  ENDDO
1295 C to hepevt
1296  CALL filhep(0,1,22,k,k,0,0,phot1,0.0,.true.)
1297  CALL filhep(0,1,22,k,k,0,0,phot2,0.0,.true.)
1298  ENDIF
1299  ENDDO
1300 C
1301  END
1302  SUBROUTINE taueta(JAK)
1303 C subroutine to decay etas's from taus.
1304 C this routine belongs to tauola universal interface, but uses
1305 C routines from tauola utilities. Just flat phase space, but 4 channels.
1306 C it is called at the beginning of SUBR. TAUPI0(JAK)
1307 C and as far as hepevt search it is basically the same as TAUPI0. 25.08.2005
1308 #include "../../include/HEPEVT.h"
1309  COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
1310  * ,ampiz,ampi,amro,gamro,ama1,gama1
1311  * ,amk,amkz,amkst,gamkst
1312 *
1313  real*4 amtau,amnuta,amel,amnue,ammu,amnumu
1314  * ,ampiz,ampi,amro,gamro,ama1,gama1
1315  * ,amk,amkz,amkst,gamkst
1316 
1317 C position of taus, must be defined by host program:
1318  COMMON /taupos/ np1,np2
1319 c
1320  REAL RRR(1),BRSUM(3), RR(2)
1321  REAL PHOT1(4),PHOT2(4),PHOT3(4)
1322  real*8 x(4), y(4), z(4)
1323  REAL YM1,YM2,YM3
1324  real*8 r,ru,peta(4),xm1,xm2,xm3,xm,xlam
1325  real*8 a,b,c
1326  xlam(a,b,c)=sqrt(abs((a-b-c)**2-4.0*b*c))
1327 C position of decaying particle:
1328  IF((kto.EQ. 1).OR.(kto.EQ.11)) THEN
1329  nps=np1
1330  ELSE
1331  nps=np2
1332  ENDIF
1333  nhepm=nhep ! to avoid infinite loop
1334  DO k=jdahep(1,nps),nhepm ! we search for etas's from tau till eor.
1335  IF (idhep(k).EQ.221.AND.jdahep(1,k).LE.k) THEN ! IF we found eta
1336  DO l=1,4
1337  peta(l)= phep(l,k) ! eta 4 momentum
1338  ENDDO
1339 C eta cumulated branching ratios:
1340  brsum(1)=0.389 ! gamma gamma
1341  brsum(2)=brsum(1)+0.319 ! 3 pi0
1342  brsum(3)=brsum(2)+0.237 ! pi+ pi- pi0 rest is thus pi+pi-gamma
1343  CALL ranmar(rrr,1)
1344 
1345  IF (rrr(1).LT.brsum(1)) THEN ! gamma gamma channel exactly like pi0
1346 ! random 3 vector on the sphere, masless
1347  r=sqrt(peta(4)**2-peta(3)**2-peta(2)**2-peta(1)**2)/2d0
1348  CALL spherd(r,x)
1349  x(4)=r
1350  y(4)=r
1351 
1352  y(1)=-x(1)
1353  y(2)=-x(2)
1354  y(3)=-x(3)
1355 ! boost to lab and to real*4
1356  CALL bostdq(-1,peta,x,x)
1357  CALL bostdq(-1,peta,y,y)
1358  DO l=1,4
1359  phot1(l)=x(l)
1360  phot2(l)=y(l)
1361  ENDDO
1362 C to hepevt
1363  CALL filhep(0,1,22,k,k,0,0,phot1,0.0,.true.)
1364  CALL filhep(0,1,22,k,k,0,0,phot2,0.0,.true.)
1365  ELSE ! 3 body channels
1366  IF(rrr(1).LT.brsum(2)) THEN ! 3 pi0
1367  id1= 111
1368  id2= 111
1369  id3= 111
1370  xm1=ampiz ! masses
1371  xm2=ampiz
1372  xm3=ampiz
1373  ELSEIF(rrr(1).LT.brsum(3)) THEN ! pi+ pi- pi0
1374  id1= 211
1375  id2=-211
1376  id3= 111
1377  xm1=ampi ! masses
1378  xm2=ampi
1379  xm3=ampiz
1380  ELSE ! pi+ pi- gamma
1381  id1= 211
1382  id2=-211
1383  id3= 22
1384  xm1=ampi ! masses
1385  xm2=ampi
1386  xm3=0.0
1387  ENDIF
1388  7 CONTINUE ! we generate mass of the first pair:
1389  CALL ranmar(rr,2)
1390  r=sqrt(peta(4)**2-peta(3)**2-peta(2)**2-peta(1)**2)
1391  amin=xm1+xm2
1392  amax=r-xm3
1393  am2=sqrt(amin**2+rr(1)*(amax**2-amin**2))
1394 C weight for flat phase space
1395  wt=xlam(1d0*r**2,1d0*am2**2,1d0*xm3**2)
1396  & *xlam(1d0*am2**2,1d0*xm1**2,1d0*xm2**2)
1397  & /r**2 /am2**2
1398  IF (rr(2).GT.wt) GOTO 7
1399 
1400  ru=xlam(1d0*am2**2,1d0*xm1**2,1d0*xm2**2)/am2/2 ! momenta of the
1401  ! first two products
1402  ! in the rest frame of that pair
1403  CALL spherd(ru,x)
1404  x(4)=sqrt(ru**2+xm1**2)
1405  y(4)=sqrt(ru**2+xm2**2)
1406 
1407  y(1)=-x(1)
1408  y(2)=-x(2)
1409  y(3)=-x(3)
1410 C generate momentum of that pair in rest frame of eta:
1411  ru=xlam(1d0*r**2,1d0*am2**2,1d0*xm3**2)/r/2
1412  CALL spherd(ru,z)
1413  z(4)=sqrt(ru**2+am2**2)
1414 C and boost first two decay products to rest frame of eta.
1415  CALL bostdq(-1,z,x,x)
1416  CALL bostdq(-1,z,y,y)
1417 C redefine Z(4) to 4-momentum of the last decay product:
1418  z(1)=-z(1)
1419  z(2)=-z(2)
1420  z(3)=-z(3)
1421  z(4)=sqrt(ru**2+xm3**2)
1422 C boost all to lab and move to real*4; also masses
1423  CALL bostdq(-1,peta,x,x)
1424  CALL bostdq(-1,peta,y,y)
1425  CALL bostdq(-1,peta,z,z)
1426  DO l=1,4
1427  phot1(l)=x(l)
1428  phot2(l)=y(l)
1429  phot3(l)=z(l)
1430  ENDDO
1431  ym1=xm1
1432  ym2=xm2
1433  ym3=xm3
1434 C to hepevt
1435  CALL filhep(0,1,id1,k,k,0,0,phot1,ym1,.true.)
1436  CALL filhep(0,1,id2,k,k,0,0,phot2,ym2,.true.)
1437  CALL filhep(0,1,id3,k,k,0,0,phot3,ym3,.true.)
1438  ENDIF
1439 
1440  ENDIF
1441  ENDDO
1442 C
1443  END
1444  SUBROUTINE tauk0s(JAK)
1445 C subroutine to decay K0S's from taus.
1446 C this routine belongs to tauola universal interface, but uses
1447 C routines from tauola utilities. Just flat phase space, but 4 channels.
1448 C it is called at the beginning of SUBR. TAUPI0(JAK)
1449 C and as far as hepevt search it is basically the same as TAUPI0. 25.08.2005
1450 #include "../../include/HEPEVT.h"
1451 
1452  COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
1453  * ,ampiz,ampi,amro,gamro,ama1,gama1
1454  * ,amk,amkz,amkst,gamkst
1455 *
1456  real*4 amtau,amnuta,amel,amnue,ammu,amnumu
1457  * ,ampiz,ampi,amro,gamro,ama1,gama1
1458  * ,amk,amkz,amkst,gamkst
1459 
1460 C position of taus, must be defined by host program:
1461  COMMON /taupos/ np1,np2
1462 c
1463  REAL RRR(1),BRSUM(3), RR(2)
1464  REAL PHOT1(4),PHOT2(4),PHOT3(4)
1465  real*8 x(4), y(4), z(4)
1466  REAL YM1,YM2,YM3
1467  real*8 r,ru,peta(4),xm1,xm2,xm3,xm,xlam
1468  real*8 a,b,c
1469  xlam(a,b,c)=sqrt(abs((a-b-c)**2-4.0*b*c))
1470 C position of decaying particle:
1471  IF((kto.EQ. 1).OR.(kto.EQ.11)) THEN
1472  nps=np1
1473  ELSE
1474  nps=np2
1475  ENDIF
1476  nhepm=nhep ! to avoid infinite loop
1477  DO k=jdahep(1,nps),nhepm ! we search for K0S's from tau till eor.
1478  IF (idhep(k).EQ.310.AND.jdahep(1,k).LE.k) THEN ! IF we found K0S
1479 
1480 
1481  DO l=1,4
1482  peta(l)= phep(l,k) ! K0S 4 momentum (this is cloned from eta decay)
1483  ENDDO
1484 C K0S cumulated branching ratios:
1485  brsum(1)=0.313 ! 2 PI0
1486  brsum(2)=1.0 ! BRSUM(1)+0.319 ! Pi+ PI-
1487  brsum(3)=brsum(2)+0.237 ! pi+ pi- pi0 rest is thus pi+pi-gamma
1488  CALL ranmar(rrr,1)
1489 
1490  IF(rrr(1).LT.brsum(1)) THEN ! 2 pi0
1491  id1= 111
1492  id2= 111
1493  xm1=ampiz ! masses
1494  xm2=ampiz
1495  ELSEIF(rrr(1).LT.brsum(2)) THEN ! pi+ pi-
1496  id1= 211
1497  id2=-211
1498  xm1=ampi ! masses
1499  xm2=ampi
1500  ELSE ! gamma gamma unused !!!
1501  id1= 22
1502  id2= 22
1503  xm1= 0.0 ! masses
1504  xm2= 0.0
1505  ENDIF
1506 
1507 ! random 3 vector on the sphere, of equal mass !!
1508  r=sqrt(peta(4)**2-peta(3)**2-peta(2)**2-peta(1)**2)/2d0
1509  r4=r
1510  r=sqrt(abs(r**2-xm1**2))
1511  CALL spherd(r,x)
1512  x(4)=r4
1513  y(4)=r4
1514 
1515  y(1)=-x(1)
1516  y(2)=-x(2)
1517  y(3)=-x(3)
1518 ! boost to lab and to real*4
1519  CALL bostdq(-1,peta,x,x)
1520  CALL bostdq(-1,peta,y,y)
1521  DO l=1,4
1522  phot1(l)=x(l)
1523  phot2(l)=y(l)
1524  ENDDO
1525 
1526  ym1=xm1
1527  ym2=xm2
1528 C to hepevt
1529  CALL filhep(0,1,id1,k,k,0,0,phot1,ym1,.true.)
1530  CALL filhep(0,1,id2,k,k,0,0,phot2,ym2,.true.)
1531 
1532 C
1533  ENDIF
1534  ENDDO
1535 
1536  END