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