C++ Interface to Tauola
tauola-F/jetset-F/tauface-jetset.f
1 /* copyright(c) 1991-2021 free software foundation, inc.
2  this file is part of the gnu c library.
3 
4  the gnu c library is free software; you can redistribute it and/or
5  modify it under the terms of the gnu lesser general Public
6  license as published by the free software foundation; either
7  version 2.1 of the license, or(at your option) any later version.
8 
9  the gnu c library is distributed in the hope that it will be useful,
10  but without any warranty; without even the implied warranty of
11  merchantability or fitness for a particular purpose. see the gnu
12  lesser general Public license for more details.
13 
14  you should have received a copy of the gnu lesser general Public
15  license along with the gnu c library; if not, see
16  <https://www.gnu.org/licenses/>. */
17 
18 
19 /* this header is separate from features.h so that the compiler can
20  include it implicitly at the start of every compilation. it must
21  not itself include <features.h> or any other header that includes
22  <features.h> because the implicit include comes before any feature
23  test macros that may be defined in a source file before it first
24  explicitly includes a system header. gcc knows the name of this
25  header in order to preinclude it. */
26 
27 /* glibc's intent is to support the IEC 559 math functionality, real
28  and complex. If the GCC (4.9 and later) predefined macros
29  specifying compiler intent are available, use them to determine
30  whether the overall intent is to support these features; otherwise,
31  presume an older compiler has intent to support these features and
32  define these macros by default. */
33 
34 
35 
36 /* wchar_t uses Unicode 10.0.0. Version 10.0 of the Unicode Standard is
37  synchronized with ISO/IEC 10646:2017, fifth edition, plus
38  the following additions from Amendment 1 to the fifth edition:
39  - 56 emoji characters
40  - 285 hentaigana
41  - 3 additional Zanabazar Square characters */
42 
43  SUBROUTINE TAUOLA(MODE,KEYPOL)
44 C *************************************
45 C general tauola interface, should work in every case until
46 C hepevt is OK, does not check if hepevt is 'clean'
47 C in particular will decay decayed taus...
48 C only longitudinal spin effects are included.
49 C in W decay v-a vertex is assumed
50 C date: 12 DEC 1998. date: 21 June 1999. date: 24 Jan 2001 date: 24 Aug 2001
51 C this is the hepevt class in old style. No d_h_ class pre-name
52  INTEGER NMXHEP
53  PARAMETER (NMXHEP=10000)
54  REAL*8 phep, vhep ! to be real*4/ *8 depending on host
55  INTEGER nevhep,nhep,isthep,idhep,jmohep,
56  $ jdahep
57  COMMON /hepevt/
58  $ nevhep, ! serial number
59  $ nhep, ! number of particles
60  $ isthep(nmxhep), ! status code
61  $ idhep(nmxhep), ! particle ident KF
62  $ jmohep(2,nmxhep), ! parent particles
63  $ jdahep(2,nmxhep), ! childreen particles
64  $ phep(5,nmxhep), ! four-momentum, mass [GeV]
65  $ vhep(4,nmxhep) ! vertex [mm]
66 * ----------------------------------------------------------------------
67  LOGICAL qedrad
68  COMMON /phoqed/
69  $ qedrad(nmxhep) ! Photos flag
70 * ----------------------------------------------------------------------
71  SAVE hepevt,phoqed
72  COMMON /TAUPOS/ NP1, NP2
73  REAL*4 PHOI(4),PHOF(4)
74  double precision Q1(4),Q2(4),P1(4),P2(4),P3(4),P4(4)
75  COMMON / MOMDEC / Q1,Q2,P1,P2,P3,P4
76 * tauola, photos and jetset overall switches
77  COMMON /LIBRA/ JAK1,JAK2,ITDKRC,IFPHOT,IFHADM,IFHADP
78  REAL*4 RRR(1)
79  LOGICAL IFPSEUDO
80  common /pseudocoup/ csc,ssc
81  REAL*4 csc,ssc
82  save pseudocoup
83  COMMON / INOUT / INUT,IOUT
84 
85 C to switch tau polarization OFF in taus
86  DIMENSION POL1(4), POL2(4)
87  double precision POL1x(4), POL2x(4)
88  INTEGER ION(3)
89  DATA POL1 /0.0,0.0,0.0,0.0/
90  DATA POL2 /0.0,0.0,0.0,0.0/
91  DATA PI /3.141592653589793238462643D0/
92 
93 C store decay vertexes
94  DIMENSION IMOTHER (20)
95  INTEGER KFHIGGS(3)
96 
97 C store daughter pointers
98  INTEGER ISON
99  COMMON /ISONS_TAU/ISON(2)
100  SAVE /ISONS_TAU/
101 
102 .EQ. IF(MODE-1) THEN
103 C ***********************
104 
105  JAK1 = 0 ! decay mode first tau
106  JAK2 = 0 ! decay mode second tau
107  ITDKRC=1.0 ! switch of radiative corrections in decay
108  IFPHOT=1.0 ! PHOTOS switch
109  IFHADM=1.0
110  IFHADP=1.0
111  POL=1.0 ! tau polarization dipswitch must be 1 or 0
112 
113  KFHIGGS(1) = 25
114  KFHIGGS(2) = 35
115  KFHIGGS(3) = 36
116  KFHIGCH = 37
117  KFZ0 = 23
118  KFGAM = 22
119  KFTAU = 15
120  KFNUE = 16
121 C couplings of the 'pseudoscalar higgs' as in CERN-TH/2003-166
122  psi=0.5*PI ! 0.15*PI
123  xmtau=1.777 ! tau mass
124  xmh=120 ! higgs boson mass
125  betah=sqrt(1d0-4*xmtau**2/xmh**2)
126  csc=cos(psi)*betah
127  ssc=sin(psi)
128 C write(*,*) ' scalar component=',csc,' pseudo-scalar component=',ssc
129 .EQ. IF (IFPHOT1) CALL PHOINI ! this if PHOTOS was not initialized earlier
130  CALL INIETC(JAK1,JAK2,ITDKRC,IFPHOT)
131  CALL INIMAS
132  CALL INIPHX(0.01d0)
133  CALL INITDK
134 C activation of pi0 and eta decays: (1) means on, (0) off
135  ION(1)=0
136  ION(2)=0
137  ION(3)=0
138  CALL TAUPI0(-1,1,ION)
139  CALL DEKAY(-1,POL1x)
140  WRITE(IOUT,7001) pol,psi,ION(1),ION(2),ION(3)
141 .EQ. ELSEIF(MODE0) THEN
142 C ***********************
143 C
144 C..... find tau-s and fill common block /TAUPOS/
145 C this is to avoid LUND history fillings. This call is optional
146  CALL PHYFIX(NSTOP,NSTART)
147 C clear mothers of the previous event
148  DO II=1,20
149  IMOTHER(II)=0
150  ENDDO
151 
152  DO II=1,2
153  ISON(II)=0
154  ENDDO
155 C ... and to find mothers giving taus.
156  NDEC = 0
157 C(BPK)--> LOOK FOR MOTHER, CHECK THAT IT IS NOT THE HISTORY ENTRY (E.G. MSTP(128)=0)
158  DO I=NSTART,NHEP
159 .EQ..AND..EQ..AND. IF(ABS(IDHEP(I))KFTAUISTHEP(I)1
160 .GE..OR..LT. $ (ISTHEP(I)125ISTHEP(I)120)) THEN
161  IMOTH=JMOHEP(1,I)
162 .EQ. DO WHILE (ABS(IDHEP(IMOTH))KFTAU) ! KEEP WALKING UP
163  IMOTH=JMOHEP(1,IMOTH)
164  ENDDO
165 .EQ..OR. IF (ISTHEP(IMOTH)3
166 .GE..AND..LE. $ (ISTHEP(IMOTH)120ISTHEP(IMOTH)125)) THEN
167  DO J=NSTART,NHEP ! WE HAVE WALKED INTO HARD RECORD
168 .EQ..AND. IF (IDHEP(J)IDHEP(IMOTH)
169 .EQ..AND. $ JMOHEP(1,J)IMOTH
170 .EQ. $ ISTHEP(J)2) THEN
171  JMOTH=J
172  GOTO 66
173  ENDIF
174  ENDDO
175  ELSE
176  JMOTH=IMOTH
177  ENDIF
178  66 CONTINUE
179  DO II=1,NDEC
180 .EQ. IF(JMOTHIMOTHER(II)) GOTO 9999
181  ENDDO
182 C(BPK)--<
183  NDEC=NDEC+1
184  IMOTHER(NDEC)= JMOTH
185  ENDIF
186  9999 CONTINUE
187  ENDDO
188 
189 C ... taus of every mother are treated in this main loop
190  DO II=1,NDEC
191  IM=IMOTHER(II)
192  NCOUNT=0
193  NP1=0
194  NP2=0
195 
196 
197 C(BPK)-->
198 C CORRECTING HEPEVT IS OUT OF QUESTION AT THIS POINT..
199  IM0=IM
200 .EQ. IF (IDHEP(JMOHEP(1,IM0))IDHEP(IM0)) IM0=JMOHEP(1,IM0)
201  ISEL=-1
202  DO I=NSTART,NHEP
203 .EQ..OR. IF (ISTHEP(I)3
204 .GE..AND..LE. $ (ISTHEP(I)120ISTHEP(I)125)) THEN ! HARD RECORD
205  GOTO 76
206  ENDIF
207  IMOTH=JMOHEP(1,I)
208 .EQ..OR..EQ. DO WHILE (IDHEP(IMOTH)IDHEP(I)ABS(IDHEP(IMOTH))KFTAU) ! KEEP WALKING UP
209  IMOTH=JMOHEP(1,IMOTH)
210  ENDDO
211 .EQ..OR..EQ..AND..EQ. IF ((IMOTHIM0IMOTHIM)ISEL-1) THEN
212  ISON(1)=I
213  ISEL=0
214 .EQ..OR..EQ..AND..EQ. ELSEIF ((IMOTHIM0IMOTHIM)ISEL0) THEN
215  ISON(2)=I
216 .NE..AND..NE..AND..EQ. ELSEIF ((IMOTHIM0IMOTHIM)ISEL0) THEN
217  ISEL=1
218  GOTO 77
219  ENDIF
220  76 CONTINUE
221  ENDDO
222  77 CONTINUE
223 C(BPK)--<
224 
225 
226 C ... we correct HEPEVT (fix developped with Catherine BISCARAT)
227 .EQ.c IF (JDAHEP(2,IM)0) THEN ! ID of second daughter was missing
228 c ISECU=1
229 c DO I=JDAHEP(1,IM)+1,NHEP ! OK lets look for it
230 .EQ..AND..EQ.c IF (JMOHEP(1,I)IMISECU1) THEN ! we have found one
231 c JDAHEP(2,IM)=I
232 .EQ..AND..NE.c ELSEIF (JMOHEP(1,I)IMISECU1) THEN ! we have found one after there
233 c JDAHEP(2,IM)=0 ! was something else, lets kill game
234 c ENDIF
235 .NE.c IF (JMOHEP(1,I)IM) ISECU=0 ! other stuff starts
236 c ENDDO
237 c ENDIF
238 
239 C ... we check whether there are just two or more tau-likes
240  DO I=ISON(1),ISON(2)
241 .EQ..OR..EQ. IF(IDHEP(I)-KFTAUIDHEP(I)-KFNUE) NCOUNT=NCOUNT+1
242 .EQ..OR..EQ. IF(IDHEP(I) KFTAUIDHEP(I) KFNUE) NCOUNT=NCOUNT+1
243  ENDDO
244 
245 C ... if there will be more we will come here again
246  666 CONTINUE
247 
248 C(BPK)-->
249  DO I=MAX(NP1+1,ISON(1)),ISON(2)
250 C(BPK)--<
251 .EQ..OR..EQ. IF(IDHEP(I)-KFTAUIDHEP(I)-KFNUE) NP1=I
252  ENDDO
253 C(BPK)-->
254  DO I=MAX(NP2+1,ISON(1)),ISON(2)
255 C(BPK)--<
256 .EQ..OR..EQ. IF(IDHEP(I) KFTAUIDHEP(I) KFNUE) NP2=I
257  ENDDO
258  DO I=1,4
259  P1(I)= PHEP(I,NP1) !momentum of tau+
260  P2(I)= PHEP(I,NP2) !momentum of tau-
261  Q1(I)= P1(I)+P2(I)
262  ENDDO
263 
264  POL1(3)= 0D0
265  POL2(3)= 0D0
266 
267 .EQ. IF(KEYPOL1) THEN
268 c.....include polarisation effect
269  CALL RANMAR(RRR,1)
270 
271 .EQ..OR..EQ..OR. IF(IDHEP(IM)KFHIGGS(1)IDHEP(IM)KFHIGGS(2)
272 .EQ. $ IDHEP(IM)KFHIGGS(3)) THEN ! case of Higgs
273 .LT. IF(RRR(1)0.5) THEN
274  POL1(3)= POL
275  POL2(3)=-POL
276  ELSE
277  POL1(3)=-POL
278  POL2(3)= POL
279  ENDIF
280 .EQ..OR..EQ. ELSEIF((IDHEP(IM)KFZ0)(IDHEP(IM)KFGAM)) THEN ! case of gamma/Z
281 C there is no angular dependence in gamma/Z polarization
282 C there is no s-dependence in gamma/Z polarization at all
283 C there is even no Z polarization in any form
284 C main reason is that nobody asked ...
285 C but it is prepared and longitudinal correlations
286 C can be included up to KORALZ standards
287 
288  POLZ0=PLZAPX(.true.,IM,NP1,NP2)
289 .LT. IF(RRR(1)POLZ0) THEN
290  POL1(3)= POL
291  POL2(3)= POL
292  ELSE
293  POL1(3)=-POL
294  POL2(3)=-POL
295  ENDIF
296 .EQ. ELSEIF(IDHEP(NP1)-IDHEP(NP2))THEN ! undef orig. only s-dep poss.
297  POLZ0=PLZAPX(.true.,IM,NP1,NP2)
298 .LT. IF(RRR(1)POLZ0) THEN
299  POL1(3)= POL
300  POL2(3)= POL
301  ELSE
302  POL1(3)=-POL
303  POL2(3)=-POL
304  ENDIF
305 .EQ. ELSEIF(ABS(IDHEP(IM))KFHIGCH) THEN ! case of charged Higgs
306  POL1(3)= POL
307  POL2(3)= POL
308  ELSE ! case of W+ or W-
309  POL1(3)= -POL
310  POL2(3)= -POL
311  ENDIF
312 c.....include polarisation effect
313  ENDIF
314 
315 .EQ..OR..EQ..OR. IF(IDHEP(IM)KFHIGGS(1)IDHEP(IM)KFHIGGS(2)
316 .EQ. $ IDHEP(IM)KFHIGGS(3)) THEN
317 .EQ..AND. IF(IDHEP(NP1)-KFTAU
318 .LE..OR..GT..AND. $ (JDAHEP(1,NP1)NP1JDAHEP(1,NP1)NHEP)
319 .EQ..AND. $ IDHEP(NP2) KFTAU
320 .LE..OR..GT. $ (JDAHEP(1,NP2)NP2JDAHEP(1,NP2)NHEP)
321  $ ) THEN
322 .EQ. IF (IDHEP(IM)KFHIGGS(1)) THEN
323  IFPSEUDO= .FALSE.
324 .EQ. ELSEIF (IDHEP(IM)KFHIGGS(2)) THEN
325  IFPSEUDO= .FALSE.
326 .EQ. ELSEIF (IDHEP(IM)KFHIGGS(3)) THEN
327  IFPSEUDO= .TRUE.
328  ELSE
329  WRITE(*,*) 'warning from tauola:'
330  WRITE(*,*) 'i stop this run, wrong idhep(im)=',IDHEP(IM)
331  STOP
332  ENDIF
333  CALL SPINHIGGS(IM,NP1,NP2,IFPSEUDO,Pol1,Pol2)
334 .EQ. IF (IFPHOT1) CALL PHOTOS(IM) ! Bremsstrahlung in Higgs decay
335  ! AFTER adding taus !!
336 
337 
338  ENDIF
339  ELSE
340 .EQ..AND. IF(IDHEP(NP1)-KFTAU
341 .LE..OR..GT. $ (JDAHEP(1,NP1)NP1JDAHEP(1,NP1)NHEP)) THEN
342 C here check on if NP1 was not decayed should be verified
343  CALL DEXAY(1,POL1)
344 .EQ. IF (IFPHOT1) CALL PHOTOS(NP1)
345  CALL TAUPI0(0,1,ION)
346  ENDIF
347 
348 .EQ..AND. IF(IDHEP(NP2) KFTAU
349 .LE..OR..GT. $ (JDAHEP(1,NP2)NP2JDAHEP(1,NP2)NHEP)) THEN
350 C here check on if NP2 was not decayed should be added
351  CALL DEXAY(2,POL2)
352 .EQ. IF (IFPHOT1) CALL PHOTOS(NP2)
353  CALL TAUPI0(0,2,ION)
354  ENDIF
355  ENDIF
356  NCOUNT=NCOUNT-2
357 .GT. IF (NCOUNT0) GOTO 666
358  ENDDO
359 
360 .EQ. ELSEIF(MODE1) THEN
361 C ***********************
362 C
363  CALL DEXAY(100,POL1)
364  CALL DEKAY(100,POL1x)
365  WRITE(IOUT,7002)
366  ENDIF
367 C *****
368  7001 FORMAT(///1X,15(5H*****)
369  $ /,' *', 25X,'*****tauola universal interface: ******',9X,1H*,
370  $ /,' *', 25X,'*****version 1.22, april 2009 (gfort)**',9X,1H*,
371  $ /,' *', 25X,'**authors: p. golonka, b. kersevan, ***',9X,1H*,
372  $ /,' *', 25X,'**t. pierzchala, e. richter-was, ******',9X,1H*,
373  $ /,' *', 25X,'****** z. was, m. worek ***************',9X,1H*,
374  $ /,' *', 25X,'**useful discussions, in particular ***',9X,1H*,
375  $ /,' *', 25X,'*with c. biscarat and s. slabospitzky**',9X,1H*,
376  $ /,' *', 25X,'****** are warmly acknowledged ********',9X,1H*,
377  $ /,' *', 25X,' ',9X,1H*,
378  $ /,' *', 25X,'********** initialization ************',9X,1H*,
379  $ /,' *',F20.5,5X,'tau polarization switch must be 1 or 0 ',9X,1H*,
380  $ /,' *',F20.5,5X,'higs scalar/pseudo mix cern-th/2003-166',9X,1H*,
381  $ /,' *',I10, 15X,'pi0 decay switch must be 1 or 0 ',9X,1H*,
382  $ /,' *',I10, 15X,'eta decay switch must be 1 or 0 ',9X,1H*,
383  $ /,' *',I10, 15X,'k0s decay switch must be 1 or 0 ',9X,1H*,
384  $ /,1X,15(5H*****)/)
385 
386  7002 FORMAT(///1X,15(5H*****)
387  $ /,' *', 25X,'*****tauola universal interface: ******',9X,1H*,
388  $ /,' *', 25X,'*****version 1.22, april 2009 (gfort)**',9X,1H*,
389  $ /,' *', 25X,'**authors: p. golonka, b. kersevan, ***',9X,1H*,
390  $ /,' *', 25X,'**t. pierzchala, e. richter-was, ******',9X,1H*,
391  $ /,' *', 25X,'****** z. was, m. worek ***************',9X,1H*,
392  $ /,' *', 25X,'**useful discussions, in particular ***',9X,1H*,
393  $ /,' *', 25X,'*with c. biscarat and s. slabospitzky**',9X,1H*,
394  $ /,' *', 25X,'****** are warmly acknowledged ********',9X,1H*,
395  $ /,' *', 25X,'****** END OF MODULE OPERATION ********',9X,1H*,
396  $ /,1X,15(5H*****)/)
397 
398  END
399 
400  SUBROUTINE SPINHIGGS(IM,NP1,NP2,IFPSEUDO,Pol1,Pol2)
401  LOGICAL IFPSEUDO
402  REAL*8 HH1,HH2,wthiggs
403  DIMENSION POL1(4), POL2(4),HH1(4),HH2(4), RRR(1)
404 ! CALL DEXAY(1,POL1) ! Kept for tests
405 ! CALL DEXAY(2,POL2) ! Kept for tests
406  INTEGER ION(3)
407  10 CONTINUE
408  CALL RANMAR(RRR,1)
409  CALL DEKAY(1,HH1)
410  CALL DEKAY(2,HH2)
411  wt=wthiggs(IFPSEUDO,HH1,HH2)
412 .GT. IF (RRR(1)WT) GOTO 10
413  CALL DEKAY(1+10,HH1)
414  CALL TAUPI0(0,1,ION)
415  CALL DEKAY(2+10,HH2)
416  CALL TAUPI0(0,2,ION)
417  END
418  FUNCTION wthiggs(IFPSEUDO,HH1,HH2)
419  LOGICAL IFPSEUDO
420  common /pseudocoup/ csc,ssc
421  REAL*4 csc,ssc
422  save pseudocoup
423  REAL*8 HH1(4),HH2(4),R(4,4),wthiggs
424  DO K=1,4
425  DO L=1,4
426  R(K,L)=0
427  ENDDO
428  ENDDO
429  WTHIGGS=0D0
430 
431  R(4,4)= 1D0 ! unpolarized part
432  R(3,3)=-1D0 ! longitudinal
433  ! other missing
434  IF (IFPSEUDO) THEN
435  R(1,1)=-1
436  R(2,2)= -1
437  R(1,1)=(csc**2-ssc**2)/(csc**2+ssc**2)
438  R(2,2)=(csc**2-ssc**2)/(csc**2+ssc**2)
439  R(1,2)=2*csc*ssc/(csc**2+ssc**2)
440  R(2,1)=-2*csc*ssc/(csc**2+ssc**2)
441  ELSE
442  R(1,1)=1
443  R(2,2)=1
444  ENDIF
445 
446 
447 
448  DO K=1,4
449  DO L=1,4
450  WTHIGGS=WTHIGGS+R(K,L)*HH1(K)*HH2(L)
451  ENDDO
452  ENDDO
453  WTHIGGS=WTHIGGS/4D0
454  END
455 
456  FUNCTION PLZAPX(HOPEin,IM0,NP1,NP2)
457 C IM0 NP1 NP2 are the positions of Z/gamma tau tau in hepevt common block.
458 C the purpose of this routine is to calculate polarization of Z along tau direction.
459 C this is highly non-trivial due to necessity of reading infromation from hard process
460 C history in HEPEVT, which is often written not up to the gramatic rules.
461  REAL*8 PLZAP0,SVAR,COSTHE,sini,sfin,ZPROP2,
462  $ P1(4),P2(4),Q1(4),Q2(4),QQ(4),PH(4),PD1(4),PD2(4),
463  $ PQ1(4),PQ2(4),PB(4),PA(4)
464  REAL*4 PLZAPX
465  INTEGER IM
466  LOGICAL HOPE,HOPEin
467 C this is the hepevt class in old style. No d_h_ class pre-name
468  INTEGER NMXHEP
469  PARAMETER (NMXHEP=10000)
470  REAL*8 phep, vhep ! to be real*4/ *8 depending on host
471  INTEGER nevhep,nhep,isthep,idhep,jmohep,
472  $ jdahep
473  COMMON /hepevt/
474  $ nevhep, ! serial number
475  $ nhep, ! number of particles
476  $ isthep(nmxhep), ! status code
477  $ idhep(nmxhep), ! particle ident KF
478  $ jmohep(2,nmxhep), ! parent particles
479  $ jdahep(2,nmxhep), ! childreen particles
480  $ phep(5,nmxhep), ! four-momentum, mass [GeV]
481  $ vhep(4,nmxhep) ! vertex [mm]
482 * ----------------------------------------------------------------------
483  LOGICAL qedrad
484  COMMON /phoqed/
485  $ qedrad(nmxhep) ! Photos flag
486 * ----------------------------------------------------------------------
487  SAVE hepevt,phoqed
488 
489 C(BPK)--> BROTHERS OF TAU ALREADY FOUND
490  INTEGER ISON
491  COMMON /ISONS_TAU/ISON(2)
492 C(BPK)--<
493 C >>
494 C >> STEP 1: find where are particles in hepevent and pick them up
495 C >>
496  HOPE=HOPEin
497 C sometimes shade Z of Z is its mother ...
498  IM=IM0
499  IM00=JMOHEP(1,IM0)
500 C to protect against check on mother of beam particles.
501 .GT. IF (IM000) THEN
502 .EQ. IF (IDHEP(IM0)IDHEP(IM00)) IM=JMOHEP(1,IM0)
503  ENDIF
504 C
505 C find (host generator-level) incoming beam-bare-particles which form Z and co.
506  IMO1=JMOHEP(1,IM)
507  IMO2=JMOHEP(2,IM)
508 
509 C(BPK)--> IN HERWIG THE POINTER MIGHT BE TO HARD CMS
510  IM00=IMO1
511 .EQ. IF (ISTHEP(IM00)120) THEN
512  IMO1=JMOHEP(1,IM00)
513  IMO2=JMOHEP(2,IM00)
514  ENDIF
515 C(BPK)--<
516 
517  IFFULL=0
518 C case when it was like e+e- --> tau+tau- gammas and e+e- were 'first' in hepevt.
519 .EQ..AND..EQ. IF (IMO10IMO20) THEN
520  IMO1=JMOHEP(1,NP1)
521 C(BPK)-->
522 .EQ. IF (IDHEP(IMO1)IDHEP(NP1)) IMO1=JMOHEP(1,IMO1) ! PROTECT AGAINST COPIES
523 C(BPK)--<
524  IMO2=JMOHEP(2,NP1)
525 C(BPK)-->
526 .EQ. IF (IDHEP(IMO2)IDHEP(NP2)) IMO2=JMOHEP(1,IMO2) ! PROTECT AGAINST COPIES
527 C(BPK)--<
528  IFFULL=1
529 C case when it was like qq --> tau+tau- gammas and qq were NOT 'first' in hepevt.
530 
531 .NE..AND..NE. ELSEIF (IDHEP(IM)22IDHEP(IM)23) THEN
532  IMO1=JMOHEP(1,NP1)
533 C(BPK)-->
534 .EQ. IF (IDHEP(IMO1)IDHEP(NP1)) IMO1=JMOHEP(1,IMO1) ! PROTECT AGAINST COPIES
535 C(BPK)--<
536  IMO2=JMOHEP(2,NP1)
537 C(BPK)-->
538 .EQ. IF (IDHEP(IMO2)IDHEP(NP2)) IMO2=JMOHEP(1,IMO2) ! PROTECT AGAINST COPIES
539 C(BPK)--<
540  IFFULL=1
541  ENDIF
542 
543 
544 C and check if it really happened
545 .EQ. IF (IMO10) HOPE=.FALSE.
546 .EQ. IF (IMO20) HOPE=.FALSE.
547 .EQ. IF (IMO1IMO2) HOPE=.FALSE.
548 
549 C
550  DO I=1,4
551  Q1(I)= PHEP(I,NP1) !momentum of tau+
552  Q2(I)= PHEP(I,NP2) !momentum of tau-
553  ENDDO
554 
555 C corrections due to possible differences in 4-momentum of shadow vs true Z.
556 C(BPK)-->
557 .EQ..AND. IF (IMJMOHEP(1,IM0)
558 .EQ..OR..EQ. $ (IDHEP(IM)22IDHEP(IM)23)) THEN
559  DO K=1,4
560  PB(K)=PHEP(K,IM)
561  PA(K)=PHEP(K,IM0)
562  ENDDO
563 C(BPK)--<
564 
565  CALL BOSTDQ( 1,PA, Q1, Q1)
566  CALL BOSTDQ( 1,PA, Q2, Q2)
567  CALL BOSTDQ(-1,PB, Q1, Q1)
568  CALL BOSTDQ(-1,PB, Q2, Q2)
569 
570  ENDIF
571 
572  DO I=1,4
573  QQ(I)= Q1(I)+Q2(I) !momentum of Z
574  IF (HOPE) P1(I)=PHEP(I,IMO1) !momentum of beam1
575  IF (HOPE) P2(I)=PHEP(I,IMO2) !momentum of beam2
576  PH(I)=0D0
577  PD1(I)=0D0
578  PD2(I)=0D0
579  ENDDO
580 ! These momenta correspond to quarks, gluons photons or taus
581  IDFQ1=IDHEP(NP1)
582  IDFQ2=IDHEP(NP2)
583  IF (HOPE) IDFP1=IDHEP(IMO1)
584  IF (HOPE) IDFP2=IDHEP(IMO2)
585 
586  SVAR=QQ(4)**2-QQ(3)**2-QQ(2)**2-QQ(1)**2
587 .NOT. IF (HOPE) THEN
588 C options which may be useful in some cases of two heavy boson production
589 C need individual considerations. To be developed.
590 C PLZAPX=PLZAP0(11,IDFQ1,SVAR,0D0) ! gamma/Z mixture as if produced from e beam
591 C PLZAPX=PLZAP0(12,IDFQ1,SVAR,0D0) ! pure Z
592  PLZAPX=0.5 ! pure gamma
593  RETURN
594  ENDIF
595 C >>
596 C >> STEP 2 look for brothers of Z which have to be included in effective incoming particles
597 C >>
598 C let us define beginning and end of particles which are produced in parallel to Z
599 C in principle following should work
600 
601 C(BPK)--> ACCOMMODATE FOR HERWIG - IM00 POINTS TO BEAM PARTICLE OR HARD CMS
602  NX1=JDAHEP(1,IM00)
603  NX2=JDAHEP(2,IM00)
604 C but ...
605  INBR=IM ! OK, HARD RECORD Z/GAMMA
606 .EQ. IF (IFFULL1) INBR=NP1 ! OK, NO Z/GAMMA
607 .EQ. IF (IDHEP(JMOHEP(1,INBR))IDHEP(INBR)) INBR=JMOHEP(1,INBR) ! FORCE HARD RECORD
608 C(BPK)--<
609 .EQ..OR..EQ. IF(NX10NX20) THEN
610  NX1=INBR
611  NX2=INBR
612  DO K=1,INBR-1
613 .EQ. IF(JMOHEP(1,INBR-K)JMOHEP(1,INBR)) THEN
614  NX1=INBR-K
615  ELSE
616  GOTO 7
617  ENDIF
618  ENDDO
619  7 CONTINUE
620 
621  DO K=INBR+1,NHEP
622 .EQ. IF(JMOHEP(1,K)JMOHEP(1,INBR)) THEN
623  NX2=K
624  ELSE
625  GOTO 8
626  ENDIF
627  ENDDO
628  8 CONTINUE
629  ENDIF
630 
631 C case of annihilation of two bosons is hopeles
632 .GE..AND..GE. IF (ABS(IDFP1)20ABS(IDFP2)20) HOPE=.FALSE.
633 C case of annihilation of non-matching flavors is hopeless
634 .LE..AND..LE..AND..NE. IF (ABS(IDFP1)20ABS(IDFP2)20IDFP1+IDFP20)
635  $ HOPE=.FALSE.
636 .NOT. IF (HOPE) THEN
637 C options which may be useful in some cases of two heavy boson production
638 C need individual considerations. To be developed.
639 C PLZAPX=PLZAP0(11,IDFQ1,SVAR,0D0) ! gamma/Z mixture as if produced from e beam
640 C PLZAPX=PLZAP0(12,IDFQ1,SVAR,0D0) ! pure Z
641  PLZAPX=0.5 ! pure gamma
642  RETURN
643  ENDIF
644 .LT. IF (ABS(IDFP1)20) IDE= IDFP1
645 .LT. IF (ABS(IDFP2)20) IDE=-IDFP2
646 
647 
648 C >>
649 C >> STEP 3 we combine gluons, photons into incoming effective beams
650 C >>
651 
652 C in the following we will ignore the possibility of photon emission from taus
653 C however at certain moment it will be necessary to take care of
654 
655  DO L=1,4
656  PD1(L)=P1(L)
657  PD2(L)=P2(L)
658  ENDDO
659 
660  DO L=1,4
661  PQ1(L)=Q1(L)
662  PQ2(L)=Q2(L)
663  ENDDO
664 
665  IFLAV=min(ABS(IDFP1),ABS(IDFP2))
666 
667 *--------------------------------------------------------------------------
668 c IFLAV=min(ABS(IDFP1),ABS(IDFP2))
669 c that means that always quark or lepton i.e. process like
670 c f g(gamma) --> f Z0 --> tau tau
671 c we glue fermions to effective beams that is f f --> Z0 --> tau tau
672 c with gamma/g emission from initial fermion.
673 *---------------------------------------------------------------------------
674 
675 .GE. IF (ABS(IDFP1)20) THEN
676  DO k=NX1,NX2
677  IDP=IDHEP(k)
678 .EQ. IF (ABS(IDP)IFLAV) THEN
679  DO L=1,4
680  PD1(L)=-PHEP(L,K)
681  ENDDO
682  ENDIF
683  ENDDO
684  ENDIF
685 
686 .GE. IF (ABS(IDFP2)20) THEN
687  DO k=NX1,NX2
688  IDP=IDHEP(k)
689 .EQ. IF (ABS(IDP)IFLAV) THEN
690  DO L=1,4
691  PD2(L)=-PHEP(L,K)
692  ENDDO
693  ENDIF
694  ENDDO
695  ENDIF
696 
697 C if first beam was boson: gluing
698 
699 .GE. IF (ABS(IDFP1)20) THEN
700  DO L=1,4
701  PH(L)=P1(L)
702  ENDDO
703  xm1=abs((PD1(4)+PH(4))**2-(PD1(3)+PH(3))**2
704  $ -(PD1(2)+PH(2))**2-(PD1(1)+PH(1))**2)
705  xm2=abs((PD2(4)+PH(4))**2-(PD2(3)+PH(3))**2
706  $ -(PD2(2)+PH(2))**2-(PD2(1)+PH(1))**2)
707 .LT. IF (XM1XM2) THEN
708  DO L=1,4
709  PD1(L)=PD1(L)+P1(L)
710  ENDDO
711  ELSE
712  DO L=1,4
713  PD2(L)=PD2(L)+P1(L)
714  ENDDO
715  ENDIF
716  ENDIF
717 
718 C if second beam was boson: gluing
719 
720 
721 .GE. IF (ABS(IDFP2)20) THEN
722  DO L=1,4
723  PH(L)=P2(L)
724  ENDDO
725  xm1=abs((PD1(4)+PH(4))**2-(PD1(3)+PH(3))**2
726  $ -(PD1(2)+PH(2))**2-(PD1(1)+PH(1))**2)
727  xm2=abs((PD2(4)+PH(4))**2-(PD2(3)+PH(3))**2
728  $ -(PD2(2)+PH(2))**2-(PD2(1)+PH(1))**2)
729 .LT. IF (XM1XM2) THEN
730  DO L=1,4
731  PD1(L)=PD1(L)+P2(L)
732  ENDDO
733  ELSE
734  DO L=1,4
735  PD2(L)=PD2(L)+P2(L)
736  ENDDO
737  ENDIF
738  ENDIF
739 
740 C now spectators ...
741 
742 C(BPK)-->
743  NPH1=NP1
744  NPH2=NP2
745 .EQ. IF (IDHEP(JMOHEP(1,NP1))IDHEP(NP1)) NPH1=JMOHEP(1,NP1) ! SHOULD PUT US IN HARD REC.
746 .EQ. IF (IDHEP(JMOHEP(1,NP2))IDHEP(NP2)) NPH2=JMOHEP(1,NP2) ! SHOULD PUT US IN HARD REC.
747 C(BPK)--<
748 
749  DO k=NX1,NX2
750 .NE..AND..NE..AND. IF (ABS(IDHEP(K))IFLAVKIM
751 C(BPK)-->
752 .NE..AND..NE. $ KNPH1KNPH2) THEN
753 C(BPK)--<
754 .EQ..AND..EQ. IF(IDHEP(K)22IFFULL1) THEN
755  DO L=1,4
756  PH(L)=PHEP(L,K)
757  ENDDO
758  xm1=abs((PD1(4)-PH(4))**2-(PD1(3)-PH(3))**2
759  $ -(PD1(2)-PH(2))**2-(PD1(1)-PH(1))**2)
760  xm2=abs((PD2(4)-PH(4))**2-(PD2(3)-PH(3))**2
761  $ -(PD2(2)-PH(2))**2-(PD2(1)-PH(1))**2)
762  xm3=abs((PQ1(4)+PH(4))**2-(PQ1(3)+PH(3))**2
763  $ -(PQ1(2)+PH(2))**2-(PQ1(1)+PH(1))**2)
764  xm4=abs((PQ2(4)+PH(4))**2-(PQ2(3)+PH(3))**2
765  $ -(PQ2(2)+PH(2))**2-(PQ2(1)+PH(1))**2)
766 
767 
768  sini=abs((PD1(4)+PD2(4)-PH(4))**2-(PD1(3)+PD2(3)-PH(3))**2
769  $ -(PD1(2)+PD2(2)-PH(2))**2-(PD1(1)+PD2(1)-PH(1))**2)
770  sfin=abs((PD1(4)+PD2(4) )**2-(PD1(3)+PD2(3) )**2
771  $ -(PD1(2)+PD2(2) )**2-(PD1(1)+PD2(1) )**2)
772 
773  FACINI=ZPROP2(sini)
774  FACFIN=ZPROP2(sfin)
775 
776  XM1=XM1/FACINI
777  XM2=XM2/FACINI
778  XM3=XM3/FACFIN
779  XM4=XM4/FACFIN
780 
781  XM=MIN(XM1,XM2,XM3,XM4)
782 .EQ. IF (XM1XM) THEN
783  DO L=1,4
784  PD1(L)=PD1(L)-PH(L)
785  ENDDO
786 .EQ. ELSEIF (XM2XM) THEN
787  DO L=1,4
788  PD2(L)=PD2(L)-PH(L)
789  ENDDO
790 .EQ. ELSEIF (XM3XM) THEN
791  DO L=1,4
792  Q1(L)=PQ1(L)+PH(L)
793  ENDDO
794  ELSE
795  DO L=1,4
796  Q2(L)=PQ2(L)+PH(L)
797  ENDDO
798  ENDIF
799  ELSE
800  DO L=1,4
801  PH(L)=PHEP(L,K)
802  ENDDO
803  xm1=abs((PD1(4)-PH(4))**2-(PD1(3)-PH(3))**2
804  $ -(PD1(2)-PH(2))**2-(PD1(1)-PH(1))**2)
805  xm2=abs((PD2(4)-PH(4))**2-(PD2(3)-PH(3))**2
806  $ -(PD2(2)-PH(2))**2-(PD2(1)-PH(1))**2)
807 .LT. IF (XM1XM2) THEN
808  DO L=1,4
809  PD1(L)=PD1(L)-PH(L)
810  ENDDO
811  ELSE
812  DO L=1,4
813  PD2(L)=PD2(L)-PH(L)
814  ENDDO
815  ENDIF
816  ENDIF
817  ENDIF
818  ENDDO
819 
820 
821 C >>
822 C >> STEP 4 look for brothers of tau (sons of Z!) which have to be included in
823 c >> effective outcoming taus
824 C >>
825 C let us define beginning and end of particles which are produced in
826 c parallel to tau
827 
828 
829 
830 C find outcoming particles which come from Z
831 
832 
833 
834 
835 C(BPK)--> OK, IT WOULD HAVE TO BE ALONG TAUS IN HARD RECORD WITH THE SAME MOTHER
836 .EQ..OR..EQ. IF (ABS(IDHEP(IM0))22abs(IDHEP(IM0))23) THEN
837  DO K=ISON(1),ISON(2)
838 .EQ. IF(ABS(IDHEP(K))22) THEN
839 C(BPK)--<
840 
841  do l=1,4
842  ph(l)=phep(l,k)
843  enddo
844 
845  xm3=abs((PQ1(4)+PH(4))**2-(PQ1(3)+PH(3))**2
846  $ -(PQ1(2)+PH(2))**2-(PQ1(1)+PH(1))**2)
847  xm4=abs((PQ2(4)+PH(4))**2-(PQ2(3)+PH(3))**2
848  $ -(PQ2(2)+PH(2))**2-(PQ2(1)+PH(1))**2)
849 
850  XM=MIN(XM3,XM4)
851 
852 .EQ. IF (XM3XM) THEN
853  DO L=1,4
854  Q1(L)=PQ1(L)+PH(L)
855  ENDDO
856  ELSE
857  DO L=1,4
858  Q2(L)=PQ2(L)+PH(L)
859  ENDDO
860  ENDIF
861  endif
862  enddo
863  ENDIF
864 
865 
866 
867 *------------------------------------------------------------------------
868 
869 
870 C out of effective momenta we calculate COSTHE and later polarization
871  CALL ANGULU(PD1,PD2,Q1,Q2,COSTHE)
872 
873  PLZAPX=PLZAP0(IDE,IDFQ1,SVAR,COSTHE)
874  END
875 
876  SUBROUTINE ANGULU(PD1,PD2,Q1,Q2,COSTHE)
877  REAL*8 PD1(4),PD2(4),Q1(4),Q2(4),COSTHE,P(4),QQ(4),QT(4)
878 C take effective beam which is less massive, it should be irrelevant
879 C but in case HEPEVT is particulary dirty may help.
880 C this routine calculate reduced system transver and cosine of scattering
881 C angle.
882 
883  XM1=ABS(PD1(4)**2-PD1(3)**2-PD1(2)**2-PD1(1)**2)
884  XM2=ABS(PD2(4)**2-PD2(3)**2-PD2(2)**2-PD2(1)**2)
885 .LT. IF (XM1XM2) THEN
886  SIGN=1D0
887  DO K=1,4
888  P(K)=PD1(K)
889  ENDDO
890  ELSE
891  SIGN=-1D0
892  DO K=1,4
893  P(K)=PD2(K)
894  ENDDO
895  ENDIF
896 C calculate space like part of P (in Z restframe)
897  DO K=1,4
898  QQ(K)=Q1(k)+Q2(K)
899  QT(K)=Q1(K)-Q2(K)
900  ENDDO
901 
902  XMQQ=SQRT(QQ(4)**2-QQ(3)**2-QQ(2)**2-QQ(1)**2)
903 
904  QTXQQ=QT(4)*QQ(4)-QT(3)*QQ(3)-QT(2)*QQ(2)-QT(1)*QQ(1)
905  DO K=1,4
906  QT(K)=QT(K)-QQ(K)*QTXQQ/XMQQ**2
907  ENDDO
908 
909  PXQQ=P(4)*QQ(4)-P(3)*QQ(3)-P(2)*QQ(2)-P(1)*QQ(1)
910  DO K=1,4
911  P(K)=P(K)-QQ(K)*PXQQ/XMQQ**2
912  ENDDO
913 C calculate costhe
914  PXP =SQRT(p(1)**2+p(2)**2+p(3)**2-p(4)**2)
915  QTXQT=SQRT(QT(3)**2+QT(2)**2+QT(1)**2-QT(4)**2)
916  PXQT =P(3)*QT(3)+P(2)*QT(2)+P(1)*QT(1)-P(4)*QT(4)
917  COSTHE=PXQT/PXP/QTXQT
918  COSTHE=COSTHE*SIGN
919  END
920 
921  FUNCTION PLZAP0(IDE,IDF,SVAR,COSTH0)
922 C this function calculates probability for the helicity +1 +1 configuration
923 C of taus for given Z/gamma transfer and COSTH0 cosine of scattering angle
924  REAL*8 PLZAP0,SVAR,COSTHE,COSTH0,T_BORN
925 
926  COSTHE=COSTH0
927 .LT.C >>>>> IF (IDE*IDF0) COSTHE=-COSTH0 ! this is probably not needed ID
928 C >>>>> of first beam is used by T_GIVIZ0 including sign
929 
930 .GT. IF (IDF0) THEN
931  CALL INITWK(IDE,IDF,SVAR)
932  ELSE
933  CALL INITWK(-IDE,-IDF,SVAR)
934  ENDIF
935  PLZAP0=T_BORN(0,SVAR,COSTHE,1D0,1D0)
936  $ /(T_BORN(0,SVAR,COSTHE,1D0,1D0)+T_BORN(0,SVAR,COSTHE,-1D0,-1D0))
937 
938 ! PLZAP0=0.5
939  END
940  FUNCTION T_BORN(MODE,SVAR,COSTHE,TA,TB)
941 C ----------------------------------------------------------------------
942 C THIS ROUTINE PROVIDES BORN CROSS SECTION. IT HAS THE SAME
943 C STRUCTURE AS FUNTIS AND FUNTIH, THUS CAN BE USED AS SIMPLER
944 C EXAMPLE OF THE METHOD APPLIED THERE
945 C INPUT PARAMETERS ARE: SVAR -- transfer
946 C COSTHE -- cosine of angle between tau+ and 1st beam
947 C TA,TB -- helicity states of tau+ tau-
948 C
949 C called by : BORNY, BORAS, BORNV, WAGA, WEIGHT
950 C ----------------------------------------------------------------------
951  IMPLICIT REAL*8(A-H,O-Z)
952  COMMON / T_BEAMPM / ENE ,AMIN,AMFIN,IDE,IDF
953  REAL*8 ENE ,AMIN,AMFIN
954  COMMON / T_GAUSPM /SS,POLN,T3E,QE,T3F,QF
955  & ,XUPGI ,XUPZI ,XUPGF ,XUPZF
956  & ,NDIAG0,NDIAGA,KEYA,KEYZ
957  & ,ITCE,JTCE,ITCF,JTCF,KOLOR
958  REAL*8 SS,POLN,T3E,QE,T3F,QF
959  & ,XUPGI(2),XUPZI(2),XUPGF(2),XUPZF(2)
960  REAL*8 SEPS1,SEPS2
961 C=====================================================================
962  COMMON / T_GSWPRM /SWSQ,AMW,AMZ,AMH,AMTOP,GAMMZ
963  REAL*8 SWSQ,AMW,AMZ,AMH,AMTOP,GAMMZ
964 C SWSQ = sin2 (theta Weinberg)
965 C AMW,AMZ = W & Z boson masses respectively
966 C AMH = the Higgs mass
967 C AMTOP = the top mass
968 C GAMMZ = Z0 width
969  COMPLEX*16 ABORN(2,2),APHOT(2,2),AZETT(2,2)
970  COMPLEX*16 XUPZFP(2),XUPZIP(2)
971  COMPLEX*16 ABORNM(2,2),APHOTM(2,2),AZETTM(2,2)
972  COMPLEX*16 PROPA,PROPZ
973  COMPLEX*16 XR,XI
974  COMPLEX*16 XUPF,XUPI,XFF(4),XFEM,XFOTA,XRHO,XKE,XKF,XKEF
975  COMPLEX*16 XTHING,XVE,XVF,XVEF
976  DATA XI/(0.D0,1.D0)/,XR/(1.D0,0.D0)/
977  DATA MODE0 /-5/
978  DATA IDE0 /-55/
979  DATA SVAR0,COST0 /-5.D0,-6.D0/
980  DATA PI /3.141592653589793238462643D0/
981  DATA SEPS1,SEPS2 /0D0,0D0/
982 C
983 C MEMORIZATION =========================================================
984 .NE..OR..NE..OR..NE. IF ( MODEMODE0SVARSVAR0COSTHECOST0
985 .OR..NE. $ IDE0IDE)THEN
986 C
987  KEYGSW=1
988 C ** PROPAGATORS
989  IDE0=IDE
990  MODE0=MODE
991  SVAR0=SVAR
992  COST0=COSTHE
993  SINTHE=SQRT(1.D0-COSTHE**2)
994  BETA=SQRT(MAX(0D0,1D0-4D0*AMFIN**2/SVAR))
995 C I MULTIPLY AXIAL COUPLING BY BETA FACTOR.
996  XUPZFP(1)=0.5D0*(XUPZF(1)+XUPZF(2))+0.5*BETA*(XUPZF(1)-XUPZF(2))
997  XUPZFP(2)=0.5D0*(XUPZF(1)+XUPZF(2))-0.5*BETA*(XUPZF(1)-XUPZF(2))
998  XUPZIP(1)=0.5D0*(XUPZI(1)+XUPZI(2))+0.5*(XUPZI(1)-XUPZI(2))
999  XUPZIP(2)=0.5D0*(XUPZI(1)+XUPZI(2))-0.5*(XUPZI(1)-XUPZI(2))
1000 C FINAL STATE VECTOR COUPLING
1001  XUPF =0.5D0*(XUPZF(1)+XUPZF(2))
1002  XUPI =0.5D0*(XUPZI(1)+XUPZI(2))
1003  XTHING =0D0
1004 
1005  PROPA =1D0/SVAR
1006  PROPZ =1D0/DCMPLX(SVAR-AMZ**2,SVAR/AMZ*GAMMZ)
1007 .EQ. IF (KEYGSW0) PROPZ=0.D0
1008  DO 50 I=1,2
1009  DO 50 J=1,2
1010  REGULA= (3-2*I)*(3-2*J) + COSTHE
1011  REGULM=-(3-2*I)*(3-2*J) * SINTHE *2.D0*AMFIN/SQRT(SVAR)
1012  APHOT(I,J)=PROPA*(XUPGI(I)*XUPGF(J)*REGULA)
1013  AZETT(I,J)=PROPZ*(XUPZIP(I)*XUPZFP(J)+XTHING)*REGULA
1014  ABORN(I,J)=APHOT(I,J)+AZETT(I,J)
1015  APHOTM(I,J)=PROPA*DCMPLX(0D0,1D0)*XUPGI(I)*XUPGF(J)*REGULM
1016  AZETTM(I,J)=PROPZ*DCMPLX(0D0,1D0)*(XUPZIP(I)*XUPF+XTHING)*REGULM
1017  ABORNM(I,J)=APHOTM(I,J)+AZETTM(I,J)
1018  50 CONTINUE
1019  ENDIF
1020 C
1021 C******************
1022 C* IN CALCULATING CROSS SECTION ONLY DIAGONAL ELEMENTS
1023 C* OF THE SPIN DENSITY MATRICES ENTER (LONGITUD. POL. ONLY.)
1024 C* HELICITY CONSERVATION EXPLICITLY OBEYED
1025  POLAR1= (SEPS1)
1026  POLAR2= (-SEPS2)
1027  BORN=0D0
1028  DO 150 I=1,2
1029  HELIC= 3-2*I
1030  DO 150 J=1,2
1031  HELIT=3-2*J
1032  FACTOR=KOLOR*(1D0+HELIC*POLAR1)*(1D0-HELIC*POLAR2)/4D0
1033  FACTOM=FACTOR*(1+HELIT*TA)*(1-HELIT*TB)
1034  FACTOR=FACTOR*(1+HELIT*TA)*(1+HELIT*TB)
1035 
1036  BORN=BORN+CDABS(ABORN(I,J))**2*FACTOR
1037 C MASS TERM IN BORN
1038 .GE. IF (MODE1) THEN
1039  BORN=BORN+CDABS(ABORNM(I,J))**2*FACTOM
1040  ENDIF
1041 
1042  150 CONTINUE
1043 C************
1044  FUNT=BORN
1045 .LT. IF(FUNT0.D0) FUNT=BORN
1046 
1047 C
1048 .GT. IF (SVAR4D0*AMFIN**2) THEN
1049 C PHASE SPACE THRESHOLD FACTOR
1050  THRESH=SQRT(1-4D0*AMFIN**2/SVAR)
1051  T_BORN= FUNT*SVAR**2*THRESH
1052  ELSE
1053  THRESH=0.D0
1054  T_BORN=0.D0
1055  ENDIF
1056 C ZW HERE WAS AN ERROR 19. 05. 1989
1057 ! write(*,*) 'KKKK ',PROPA,PROPZ,XUPGI,XUPGF,XUPZI,XUPZF
1058 ! write(*,*) 'KKKK X',svar,costhe,TA,TB,T_BORN
1059  END
1060 
1061  SUBROUTINE INITWK(IDEX,IDFX,SVAR)
1062 ! initialization routine coupling masses etc.
1063  IMPLICIT REAL*8 (A-H,O-Z)
1064  COMMON / T_BEAMPM / ENE ,AMIN,AMFIN,IDE,IDF
1065  REAL*8 ENE ,AMIN,AMFIN
1066  COMMON / T_GAUSPM /SS,POLN,T3E,QE,T3F,QF
1067  & ,XUPGI ,XUPZI ,XUPGF ,XUPZF
1068  & ,NDIAG0,NDIAGA,KEYA,KEYZ
1069  & ,ITCE,JTCE,ITCF,JTCF,KOLOR
1070  REAL*8 SS,POLN,T3E,QE,T3F,QF
1071  & ,XUPGI(2),XUPZI(2),XUPGF(2),XUPZF(2)
1072  COMMON / T_GSWPRM /SWSQ,AMW,AMZ,AMH,AMTOP,GAMMZ
1073  REAL*8 SWSQ,AMW,AMZ,AMH,AMTOP,GAMMZ
1074 C SWSQ = sin2 (theta Weinberg)
1075 C AMW,AMZ = W & Z boson masses respectively
1076 C AMH = the Higgs mass
1077 C AMTOP = the top mass
1078 C GAMMZ = Z0 width
1079 C
1080  ENE=SQRT(SVAR)/2
1081  AMIN=0.511D-3
1082  SWSQ=0.23147
1083  AMZ=91.1882
1084  GAMMZ=2.4952
1085 .EQ. IF (IDFX 15) then
1086  IDF=2 ! denotes tau +2 tau-
1087  AMFIN=1.77703 !this mass is irrelevant if small, used in ME only
1088 .EQ. ELSEIF (IDFX-15) then
1089  IDF=-2 ! denotes tau -2 tau-
1090  AMFIN=1.77703 !this mass is irrelevant if small, used in ME only
1091  ELSE
1092  WRITE(*,*) 'INITWK: WRONG IDFX'
1093  STOP
1094  ENDIF
1095 
1096 .EQ. IF (IDEX 11) then !electron
1097  IDE= 2
1098  AMIN=0.511D-3
1099 .EQ. ELSEIF (IDEX-11) then !positron
1100  IDE=-2
1101  AMIN=0.511D-3
1102 .EQ. ELSEIF (IDEX 13) then !mu+
1103  IDE= 2
1104  AMIN=0.105659
1105 .EQ. ELSEIF (IDEX-13) then !mu-
1106  IDE=-2
1107  AMIN=0.105659
1108 .EQ. ELSEIF (IDEX 1) then !d
1109  IDE= 4
1110  AMIN=0.05
1111 .EQ. ELSEIF (IDEX- 1) then !d~
1112  IDE=-4
1113  AMIN=0.05
1114 .EQ. ELSEIF (IDEX 2) then !u
1115  IDE= 3
1116  AMIN=0.02
1117 .EQ. ELSEIF (IDEX- 2) then !u~
1118  IDE=-3
1119  AMIN=0.02
1120 .EQ. ELSEIF (IDEX 3) then !s
1121  IDE= 4
1122  AMIN=0.3
1123 .EQ. ELSEIF (IDEX- 3) then !s~
1124  IDE=-4
1125  AMIN=0.3
1126 .EQ. ELSEIF (IDEX 4) then !c
1127  IDE= 3
1128  AMIN=1.3
1129 .EQ. ELSEIF (IDEX- 4) then !c~
1130  IDE=-3
1131  AMIN=1.3
1132 .EQ. ELSEIF (IDEX 5) then !b
1133  IDE= 4
1134  AMIN=4.5
1135 .EQ. ELSEIF (IDEX- 5) then !b~
1136  IDE=-4
1137  AMIN=4.5
1138 .EQ. ELSEIF (IDEX 12) then !nu_e
1139  IDE= 1
1140  AMIN=0.1D-3
1141 .EQ. ELSEIF (IDEX- 12) then !nu_e~
1142  IDE=-1
1143  AMIN=0.1D-3
1144 .EQ. ELSEIF (IDEX 14) then !nu_mu
1145  IDE= 1
1146  AMIN=0.1D-3
1147 .EQ. ELSEIF (IDEX- 14) then !nu_mu~
1148  IDE=-1
1149  AMIN=0.1D-3
1150 .EQ. ELSEIF (IDEX 16) then !nu_tau
1151  IDE= 1
1152  AMIN=0.1D-3
1153 .EQ. ELSEIF (IDEX- 16) then !nu_tau~
1154  IDE=-1
1155  AMIN=0.1D-3
1156 
1157  ELSE
1158  WRITE(*,*) 'INITWK: WRONG IDEX'
1159  STOP
1160  ENDIF
1161 
1162 C ----------------------------------------------------------------------
1163 C
1164 C INITIALISATION OF COUPLING CONSTANTS AND FERMION-GAMMA / Z0 VERTEX
1165 C
1166 C called by : KORALZ
1167 C ----------------------------------------------------------------------
1168  ITCE=IDE/IABS(IDE)
1169  JTCE=(1-ITCE)/2
1170  ITCF=IDF/IABS(IDF)
1171  JTCF=(1-ITCF)/2
1172  CALL T_GIVIZO( IDE, 1,AIZOR,QE,KDUMM)
1173  CALL T_GIVIZO( IDE,-1,AIZOL,QE,KDUMM)
1174  XUPGI(1)=QE
1175  XUPGI(2)=QE
1176  T3E = AIZOL+AIZOR
1177  XUPZI(1)=(AIZOR-QE*SWSQ)/SQRT(SWSQ*(1-SWSQ))
1178  XUPZI(2)=(AIZOL-QE*SWSQ)/SQRT(SWSQ*(1-SWSQ))
1179  CALL T_GIVIZO( IDF, 1,AIZOR,QF,KOLOR)
1180  CALL T_GIVIZO( IDF,-1,AIZOL,QF,KOLOR)
1181  XUPGF(1)=QF
1182  XUPGF(2)=QF
1183  T3F = AIZOL+AIZOR
1184  XUPZF(1)=(AIZOR-QF*SWSQ)/SQRT(SWSQ*(1-SWSQ))
1185  XUPZF(2)=(AIZOL-QF*SWSQ)/SQRT(SWSQ*(1-SWSQ))
1186 C
1187  NDIAG0=2
1188  NDIAGA=11
1189  KEYA = 1
1190  KEYZ = 1
1191 C
1192 C
1193  RETURN
1194  END
1195 
1196  SUBROUTINE T_GIVIZO(IDFERM,IHELIC,SIZO3,CHARGE,KOLOR)
1197 C ----------------------------------------------------------------------
1198 C PROVIDES ELECTRIC CHARGE AND WEAK IZOSPIN OF A FAMILY FERMION
1199 C IDFERM=1,2,3,4 DENOTES NEUTRINO, LEPTON, UP AND DOWN QUARK
1200 C NEGATIVE IDFERM=-1,-2,-3,-4, DENOTES ANTIPARTICLE
1201 C IHELIC=+1,-1 DENOTES RIGHT AND LEFT HANDEDNES ( CHIRALITY)
1202 C SIZO3 IS THIRD PROJECTION OF WEAK IZOSPIN (PLUS MINUS HALF)
1203 C AND CHARGE IS ELECTRIC CHARGE IN UNITS OF ELECTRON CHARGE
1204 C KOLOR IS A QCD COLOUR, 1 FOR LEPTON, 3 FOR QUARKS
1205 C
1206 C called by : EVENTE, EVENTM, FUNTIH, .....
1207 C ----------------------------------------------------------------------
1208  IMPLICIT REAL*8(A-H,O-Z)
1209 C
1210 .EQ..OR..GT. IF(IDFERM0IABS(IDFERM)4) GOTO 901
1211 .NE. IF(IABS(IHELIC)1) GOTO 901
1212  IH =IHELIC
1213  IDTYPE =IABS(IDFERM)
1214  IC =IDFERM/IDTYPE
1215  LEPQUA=INT(IDTYPE*0.4999999D0)
1216  IUPDOW=IDTYPE-2*LEPQUA-1
1217  CHARGE =(-IUPDOW+2D0/3D0*LEPQUA)*IC
1218  SIZO3 =0.25D0*(IC-IH)*(1-2*IUPDOW)
1219  KOLOR=1+2*LEPQUA
1220 C** NOTE THAT CONVENTIONALY Z0 COUPLING IS
1221 C** XOUPZ=(SIZO3-CHARGE*SWSQ)/SQRT(SWSQ*(1-SWSQ))
1222  RETURN
1223  901 PRINT *,' STOP IN GIVIZO: WRONG PARAMS.'
1224  STOP
1225  END
1226  SUBROUTINE PHYFIX(NSTOP,NSTART)
1227  COMMON/LUJETS/N,K(4000,5),P(4000,5),V(4000,5)
1228  SAVE /LUJETS/
1229 C NSTOP NSTART : when PHYTIA history ends and event starts.
1230  NSTOP=0
1231  NSTART=1
1232  DO I=1, N
1233 .NE. IF(K(I,1)21) THEN
1234  NSTOP = I-1
1235  NSTART= I
1236  GOTO 500
1237  ENDIF
1238  ENDDO
1239  500 CONTINUE
1240  END
1241  SUBROUTINE FILHEP(N,IST,ID,JMO1,JMO2,JDA1,JDA2,P4,PINV,PHFLAG)
1242 C ----------------------------------------------------------------------
1243 C this subroutine fills one entry into the HEPEVT common
1244 C and updates the information for affected mother entries
1245 C
1246 C written by Martin W. Gruenewald (91/01/28)
1247 C
1248 C called by : ZTOHEP,BTOHEP,DWLUxy
1249 C ----------------------------------------------------------------------
1250 C
1251 C this is the hepevt class in old style. No d_h_ class pre-name
1252 C this is the hepevt class in old style. No d_h_ class pre-name
1253  INTEGER NMXHEP
1254  PARAMETER (NMXHEP=10000)
1255  REAL*8 phep, vhep ! to be real*4/ *8 depending on host
1256  INTEGER nevhep,nhep,isthep,idhep,jmohep,
1257  $ jdahep
1258  COMMON /hepevt/
1259  $ nevhep, ! serial number
1260  $ nhep, ! number of particles
1261  $ isthep(nmxhep), ! status code
1262  $ idhep(nmxhep), ! particle ident KF
1263  $ jmohep(2,nmxhep), ! parent particles
1264  $ jdahep(2,nmxhep), ! childreen particles
1265  $ phep(5,nmxhep), ! four-momentum, mass [GeV]
1266  $ vhep(4,nmxhep) ! vertex [mm]
1267 * ----------------------------------------------------------------------
1268  LOGICAL qedrad
1269  COMMON /phoqed/
1270  $ qedrad(nmxhep) ! Photos flag
1271 * ----------------------------------------------------------------------
1272  SAVE hepevt,phoqed
1273  LOGICAL PHFLAG
1274 C
1275  REAL*4 P4(4)
1276 C
1277 C check address mode
1278 .EQ. IF (N0) THEN
1279 C
1280 C append mode
1281  IHEP=NHEP+1
1282 .GT. ELSE IF (N0) THEN
1283 C
1284 C absolute position
1285  IHEP=N
1286  ELSE
1287 C
1288 C relative position
1289  IHEP=NHEP+N
1290  END IF
1291 C
1292 C check on IHEP
1293 .LE..OR..GT. IF ((IHEP0)(IHEPNMXHEP)) RETURN
1294 C
1295 C add entry
1296  NHEP=IHEP
1297  ISTHEP(IHEP)=IST
1298  IDHEP(IHEP)=ID
1299  JMOHEP(1,IHEP)=JMO1
1300 .LT. IF(JMO10)JMOHEP(1,IHEP)=JMOHEP(1,IHEP)+IHEP
1301  JMOHEP(2,IHEP)=JMO2
1302 .LT. IF(JMO20)JMOHEP(2,IHEP)=JMOHEP(2,IHEP)+IHEP
1303  JDAHEP(1,IHEP)=JDA1
1304  JDAHEP(2,IHEP)=JDA2
1305 C
1306  DO I=1,4
1307  PHEP(I,IHEP)=P4(I)
1308 C
1309 C KORAL-B and KORAL-Z do not provide vertex and/or lifetime informations
1310  VHEP(I,IHEP)=0.0
1311  END DO
1312  PHEP(5,IHEP)=PINV
1313 C FLAG FOR PHOTOS...
1314  QEDRAD(IHEP)=PHFLAG
1315 C
1316 C update process:
1317  DO IP=JMOHEP(1,IHEP),JMOHEP(2,IHEP)
1318 .GT. IF(IP0)THEN
1319 C
1320 C if there is a daughter at IHEP, mother entry at IP has decayed
1321 .EQ. IF(ISTHEP(IP)1)ISTHEP(IP)=2
1322 C
1323 C and daughter pointers of mother entry must be updated
1324 .EQ. IF(JDAHEP(1,IP)0)THEN
1325  JDAHEP(1,IP)=IHEP
1326  JDAHEP(2,IP)=IHEP
1327  ELSE
1328  JDAHEP(2,IP)=MAX(IHEP,JDAHEP(2,IP))
1329  END IF
1330  END IF
1331  END DO
1332 C
1333  RETURN
1334  END
1335 
1336 
1337  FUNCTION IHEPDIM(DUM)
1338 C this is the hepevt class in old style. No d_h_ class pre-name
1339 C this is the hepevt class in old style. No d_h_ class pre-name
1340  INTEGER NMXHEP
1341  PARAMETER (NMXHEP=10000)
1342  REAL*8 phep, vhep ! to be real*4/ *8 depending on host
1343  INTEGER nevhep,nhep,isthep,idhep,jmohep,
1344  $ jdahep
1345  COMMON /hepevt/
1346  $ nevhep, ! serial number
1347  $ nhep, ! number of particles
1348  $ isthep(nmxhep), ! status code
1349  $ idhep(nmxhep), ! particle ident KF
1350  $ jmohep(2,nmxhep), ! parent particles
1351  $ jdahep(2,nmxhep), ! childreen particles
1352  $ phep(5,nmxhep), ! four-momentum, mass [GeV]
1353  $ vhep(4,nmxhep) ! vertex [mm]
1354 * ----------------------------------------------------------------------
1355  LOGICAL qedrad
1356  COMMON /phoqed/
1357  $ qedrad(nmxhep) ! Photos flag
1358 * ----------------------------------------------------------------------
1359  SAVE hepevt,phoqed
1360  IHEPDIM=NHEP
1361  END
1362  FUNCTION ZPROP2(S)
1363  IMPLICIT REAL*8(A-H,O-Z)
1364  COMPLEX*16 CPRZ0,CPRZ0M
1365  AMZ=91.1882
1366  GAMMZ=2.49
1367  CPRZ0=DCMPLX((S-AMZ**2),S/AMZ*GAMMZ)
1368  CPRZ0M=1/CPRZ0
1369  ZPROP2=(ABS(CPRZ0M))**2
1370  END
1371 
1372  SUBROUTINE TAUPI0(MODE,JAK,ION)
1373 C no initialization required. Must be called once after every:
1374 C 1) CALL DEKAY(1+10,...)
1375 C 2) CALL DEKAY(2+10,...)
1376 C 3) CALL DEXAY(1,...)
1377 C 4) CALL DEXAY(2,...)
1378 C subroutine to decay originating from TAUOLA's taus:
1379 c 1) etas(with CALL taueta(jak))
1380 c 2) later pi0's from taus.
1381 C 3) extensions to other applications possible.
1382 C this routine belongs to >tauola universal interface<, but uses
1383 C routines from >tauola< utilities as well. 25.08.2005
1384 C this is the hepevt class in old style. No d_h_ class pre-name
1385  INTEGER NMXHEP
1386  PARAMETER (NMXHEP=10000)
1387  REAL*8 phep, vhep ! to be real*4/ *8 depending on host
1388  INTEGER nevhep,nhep,isthep,idhep,jmohep,
1389  $ jdahep
1390  COMMON /hepevt/
1391  $ nevhep, ! serial number
1392  $ nhep, ! number of particles
1393  $ isthep(nmxhep), ! status code
1394  $ idhep(nmxhep), ! particle ident KF
1395  $ jmohep(2,nmxhep), ! parent particles
1396  $ jdahep(2,nmxhep), ! childreen particles
1397  $ phep(5,nmxhep), ! four-momentum, mass [GeV]
1398  $ vhep(4,nmxhep) ! vertex [mm]
1399 * ----------------------------------------------------------------------
1400  LOGICAL qedrad
1401  COMMON /phoqed/
1402  $ qedrad(nmxhep) ! Photos flag
1403 * ----------------------------------------------------------------------
1404  SAVE hepevt,phoqed
1405 
1406 C position of taus, must be defined by host program:
1407  COMMON /TAUPOS/ NP1,NP2
1408 c
1409  REAL PHOT1(4),PHOT2(4)
1410  REAL*8 R,X(4),Y(4),PI0(4)
1411  INTEGER JEZELI(3),ION(3)
1412  DATA JEZELI /0,0,0/
1413  SAVE JEZELI
1414 .EQ. IF (MODE-1) THEN
1415  JEZELI(1)=ION(1)
1416  JEZELI(2)=ION(2)
1417  JEZELI(3)=ION(3)
1418  RETURN
1419  ENDIF
1420 .EQ. IF (JEZELI(1)0) RETURN
1421 .EQ. IF (JEZELI(2)1) CALL TAUETA(JAK)
1422 .EQ. IF (JEZELI(3)1) CALL TAUK0S(JAK)
1423 C position of decaying particle:
1424 .EQ..OR..EQ. IF((KTO 1)(KTO11)) THEN
1425  NPS=NP1
1426  ELSE
1427  NPS=NP2
1428  ENDIF
1429  nhepM=nhep ! to avoid infinite loop
1430  DO K=JDAHEP(1,NPS),nhepM ! we search for pi0's from tau till eor.
1431  IF (idhep(k).EQ.111.AND.jdahep(1,k).LE.k) THEN ! IF we found pi0
1432  DO l=1,4
1433  pi0(l)= phep(l,k)
1434  ENDDO
1435 ! random 3 vector on the sphere, masless
1436  r=sqrt(pi0(4)**2-pi0(3)**2-pi0(2)**2-pi0(1)**2)/2d0
1437  CALL spherd(r,x)
1438  x(4)=r
1439  y(4)=r
1440 
1441  y(1)=-x(1)
1442  y(2)=-x(2)
1443  y(3)=-x(3)
1444 ! boost to lab and to real*4
1445  CALL bostdq(-1,pi0,x,x)
1446  CALL bostdq(-1,pi0,y,y)
1447  DO l=1,4
1448  phot1(l)=x(l)
1449  phot2(l)=y(l)
1450  ENDDO
1451 c to hepevt
1452  CALL filhep(0,1,22,k,k,0,0,phot1,0.0,.true.)
1453  CALL filhep(0,1,22,k,k,0,0,phot2,0.0,.true.)
1454  ENDIF
1455  ENDDO
1456 c
1457  END
1458  SUBROUTINE taueta(JAK)
1459 c subroutine to decay etas's from taus.
1460 C this routine belongs to tauola universal interface, but uses
1461 C routines from tauola utilities. Just flat phase space, but 4 channels.
1462 C it is called at the beginning of SUBR. TAUPI0(JAK)
1463 C and as far as hepevt search it is basically the same as TAUPI0. 25.08.2005
1464 C this is the hepevt class in old style. No d_h_ class pre-name
1465  INTEGER NMXHEP
1466  PARAMETER (NMXHEP=10000)
1467  REAL*8 phep, vhep ! to be real*4/ *8 depending on host
1468  INTEGER nevhep,nhep,isthep,idhep,jmohep,
1469  $ jdahep
1470  COMMON /hepevt/
1471  $ nevhep, ! serial number
1472  $ nhep, ! number of particles
1473  $ isthep(nmxhep), ! status code
1474  $ idhep(nmxhep), ! particle ident KF
1475  $ jmohep(2,nmxhep), ! parent particles
1476  $ jdahep(2,nmxhep), ! childreen particles
1477  $ phep(5,nmxhep), ! four-momentum, mass [GeV]
1478  $ vhep(4,nmxhep) ! vertex [mm]
1479 * ----------------------------------------------------------------------
1480  LOGICAL qedrad
1481  COMMON /phoqed/
1482  $ qedrad(nmxhep) ! Photos flag
1483 * ----------------------------------------------------------------------
1484  SAVE hepevt,phoqed
1485  COMMON / PARMAS / AMTAU,AMNUTA,AMEL,AMNUE,AMMU,AMNUMU
1486  * ,AMPIZ,AMPI,AMRO,GAMRO,AMA1,GAMA1
1487  * ,AMK,AMKZ,AMKST,GAMKST
1488 *
1489  REAL*4 AMTAU,AMNUTA,AMEL,AMNUE,AMMU,AMNUMU
1490  * ,AMPIZ,AMPI,AMRO,GAMRO,AMA1,GAMA1
1491  * ,AMK,AMKZ,AMKST,GAMKST
1492 
1493 C position of taus, must be defined by host program:
1494  COMMON /TAUPOS/ NP1,NP2
1495 c
1496  REAL RRR(1),BRSUM(3), RR(2)
1497  REAL PHOT1(4),PHOT2(4),PHOT3(4)
1498  REAL*8 X(4), Y(4), Z(4)
1499  REAL YM1,YM2,YM3
1500  REAL*8 R,RU,PETA(4),XM1,XM2,XM3,XM,XLAM
1501  REAL*8 a,b,c
1502  XLAM(a,b,c)=SQRT(ABS((a-b-c)**2-4.0*b*c))
1503 C position of decaying particle:
1504 .EQ..OR..EQ. IF((KTO 1)(KTO11)) THEN
1505  NPS=NP1
1506  ELSE
1507  NPS=NP2
1508  ENDIF
1509  nhepM=nhep ! to avoid infinite loop
1510  DO K=JDAHEP(1,NPS),nhepM ! we search for etas's from tau till eor.
1511  IF (idhep(k).EQ.221.AND.jdahep(1,k).LE.k) THEN ! IF we found eta
1512  DO l=1,4
1513  peta(l)= phep(l,k) ! eta 4 momentum
1514  ENDDO
1515 c eta cumulated branching ratios:
1516  brsum(1)=0.389 ! gamma gamma
1517  brsum(2)=brsum(1)+0.319 ! 3 pi0
1518  brsum(3)=brsum(2)+0.237 ! pi+ pi- pi0 rest is thus pi+pi-gamma
1519  CALL ranmar(rrr,1)
1520 
1521  IF (rrr(1).LT.brsum(1)) THEN ! gamma gamma channel exactly like pi0
1522 ! random 3 vector on the sphere, masless
1523  r=sqrt(peta(4)**2-peta(3)**2-peta(2)**2-peta(1)**2)/2d0
1524  CALL spherd(r,x)
1525  x(4)=r
1526  y(4)=r
1527 
1528  y(1)=-x(1)
1529  y(2)=-x(2)
1530  y(3)=-x(3)
1531 ! boost to lab and to real*4
1532  CALL bostdq(-1,peta,x,x)
1533  CALL bostdq(-1,peta,y,y)
1534  DO l=1,4
1535  phot1(l)=x(l)
1536  phot2(l)=y(l)
1537  ENDDO
1538 c to hepevt
1539  CALL filhep(0,1,22,k,k,0,0,phot1,0.0,.true.)
1540  CALL filhep(0,1,22,k,k,0,0,phot2,0.0,.true.)
1541  ELSE ! 3 body channels
1542  IF(rrr(1).LT.brsum(2)) THEN ! 3 pi0
1543  id1= 111
1544  id2= 111
1545  id3= 111
1546  xm1=ampiz ! masses
1547  xm2=ampiz
1548  xm3=ampiz
1549  ELSEIF(rrr(1).LT.brsum(3)) THEN ! pi+ pi- pi0
1550  id1= 211
1551  id2=-211
1552  id3= 111
1553  xm1=ampi ! masses
1554  xm2=ampi
1555  xm3=ampiz
1556  ELSE ! pi+ pi- gamma
1557  id1= 211
1558  id2=-211
1559  id3= 22
1560  xm1=ampi ! masses
1561  xm2=ampi
1562  xm3=0.0
1563  ENDIF
1564  7 CONTINUE ! we generate mass of the first pair:
1565  CALL ranmar(rr,2)
1566  r=sqrt(peta(4)**2-peta(3)**2-peta(2)**2-peta(1)**2)
1567  amin=xm1+xm2
1568  amax=r-xm3
1569  am2=sqrt(amin**2+rr(1)*(amax**2-amin**2))
1570 c weight for flat phase space
1571  wt=xlam(1d0*r**2,1d0*am2**2,1d0*xm3**2)
1572  & *xlam(1d0*am2**2,1d0*xm1**2,1d0*xm2**2)
1573  & /r**2 /am2**2
1574  IF (rr(2).GT.wt) GOTO 7
1575 
1576  ru=xlam(1d0*am2**2,1d0*xm1**2,1d0*xm2**2)/am2/2 ! momenta of the
1577  ! first two products
1578  ! in the rest frame of that pair
1579  CALL spherd(ru,x)
1580  x(4)=sqrt(ru**2+xm1**2)
1581  y(4)=sqrt(ru**2+xm2**2)
1582 
1583  y(1)=-x(1)
1584  y(2)=-x(2)
1585  y(3)=-x(3)
1586 c generate momentum of that pair in rest frame of eta:
1587  ru=xlam(1d0*r**2,1d0*am2**2,1d0*xm3**2)/r/2
1588  CALL spherd(ru,z)
1589  z(4)=sqrt(ru**2+am2**2)
1590 c and boost first two decay products to rest frame of eta.
1591  CALL bostdq(-1,z,x,x)
1592  CALL bostdq(-1,z,y,y)
1593 c redefine z(4) to 4-momentum of the last decay product:
1594  z(1)=-z(1)
1595  z(2)=-z(2)
1596  z(3)=-z(3)
1597  z(4)=sqrt(ru**2+xm3**2)
1598 c boost all to lab and move to real*4; also masses
1599  CALL bostdq(-1,peta,x,x)
1600  CALL bostdq(-1,peta,y,y)
1601  CALL bostdq(-1,peta,z,z)
1602  DO l=1,4
1603  phot1(l)=x(l)
1604  phot2(l)=y(l)
1605  phot3(l)=z(l)
1606  ENDDO
1607  ym1=xm1
1608  ym2=xm2
1609  ym3=xm3
1610 c to hepevt
1611  CALL filhep(0,1,id1,k,k,0,0,phot1,ym1,.true.)
1612  CALL filhep(0,1,id2,k,k,0,0,phot2,ym2,.true.)
1613  CALL filhep(0,1,id3,k,k,0,0,phot3,ym3,.true.)
1614  ENDIF
1615 
1616  ENDIF
1617  ENDDO
1618 c
1619  END
1620  SUBROUTINE tauk0s(JAK)
1621 c subroutine to decay k0s's from taus.
1622 C this routine belongs to tauola universal interface, but uses
1623 C routines from tauola utilities. Just flat phase space, but 4 channels.
1624 C it is called at the beginning of SUBR. TAUPI0(JAK)
1625 C and as far as hepevt search it is basically the same as TAUPI0. 25.08.2005
1626 C this is the hepevt class in old style. No d_h_ class pre-name
1627  INTEGER NMXHEP
1628  PARAMETER (NMXHEP=10000)
1629  REAL*8 phep, vhep ! to be real*4/ *8 depending on host
1630  INTEGER nevhep,nhep,isthep,idhep,jmohep,
1631  $ jdahep
1632  COMMON /hepevt/
1633  $ nevhep, ! serial number
1634  $ nhep, ! number of particles
1635  $ isthep(nmxhep), ! status code
1636  $ idhep(nmxhep), ! particle ident KF
1637  $ jmohep(2,nmxhep), ! parent particles
1638  $ jdahep(2,nmxhep), ! childreen particles
1639  $ phep(5,nmxhep), ! four-momentum, mass [GeV]
1640  $ vhep(4,nmxhep) ! vertex [mm]
1641 * ----------------------------------------------------------------------
1642  LOGICAL qedrad
1643  COMMON /phoqed/
1644  $ qedrad(nmxhep) ! Photos flag
1645 * ----------------------------------------------------------------------
1646  SAVE hepevt,phoqed
1647 
1648  COMMON / PARMAS / AMTAU,AMNUTA,AMEL,AMNUE,AMMU,AMNUMU
1649  * ,AMPIZ,AMPI,AMRO,GAMRO,AMA1,GAMA1
1650  * ,AMK,AMKZ,AMKST,GAMKST
1651 *
1652  REAL*4 AMTAU,AMNUTA,AMEL,AMNUE,AMMU,AMNUMU
1653  * ,AMPIZ,AMPI,AMRO,GAMRO,AMA1,GAMA1
1654  * ,AMK,AMKZ,AMKST,GAMKST
1655 
1656 C position of taus, must be defined by host program:
1657  COMMON /TAUPOS/ NP1,NP2
1658 c
1659  REAL RRR(1),BRSUM(3), RR(2)
1660  REAL PHOT1(4),PHOT2(4),PHOT3(4)
1661  REAL*8 X(4), Y(4), Z(4)
1662  REAL YM1,YM2,YM3
1663  REAL*8 R,RU,PETA(4),XM1,XM2,XM3,XM,XLAM
1664  REAL*8 a,b,c
1665  XLAM(a,b,c)=SQRT(ABS((a-b-c)**2-4.0*b*c))
1666 C position of decaying particle:
1667 .EQ..OR..EQ. IF((KTO 1)(KTO11)) THEN
1668  NPS=NP1
1669  ELSE
1670  NPS=NP2
1671  ENDIF
1672  nhepM=nhep ! to avoid infinite loop
1673  DO K=JDAHEP(1,NPS),nhepM ! we search for K0S's from tau till eor.
1674  IF (idhep(k).EQ.310.AND.jdahep(1,k).LE.k) THEN ! IF we found K0S
1675 
1676 
1677  DO l=1,4
1678  peta(l)= phep(l,k) ! K0S 4 momentum (this is cloned from eta decay)
1679  ENDDO
1680 c k0s cumulated branching ratios:
1681  brsum(1)=0.313 ! 2 PI0
1682  brsum(2)=1.0 ! BRSUM(1)+0.319 ! Pi+ PI-
1683  brsum(3)=brsum(2)+0.237 ! pi+ pi- pi0 rest is thus pi+pi-gamma
1684  CALL ranmar(rrr,1)
1685 
1686  IF(rrr(1).LT.brsum(1)) THEN ! 2 pi0
1687  id1= 111
1688  id2= 111
1689  xm1=ampiz ! masses
1690  xm2=ampiz
1691  ELSEIF(rrr(1).LT.brsum(2)) THEN ! pi+ pi-
1692  id1= 211
1693  id2=-211
1694  xm1=ampi ! masses
1695  xm2=ampi
1696  ELSE ! gamma gamma unused !!!
1697  id1= 22
1698  id2= 22
1699  xm1= 0.0 ! masses
1700  xm2= 0.0
1701  ENDIF
1702 
1703 ! random 3 vector on the sphere, of equal mass !!
1704  r=sqrt(peta(4)**2-peta(3)**2-peta(2)**2-peta(1)**2)/2d0
1705  r4=r
1706  r=sqrt(abs(r**2-xm1**2))
1707  CALL spherd(r,x)
1708  x(4)=r4
1709  y(4)=r4
1710 
1711  y(1)=-x(1)
1712  y(2)=-x(2)
1713  y(3)=-x(3)
1714 ! boost to lab and to real*4
1715  CALL bostdq(-1,peta,x,x)
1716  CALL bostdq(-1,peta,y,y)
1717  DO l=1,4
1718  phot1(l)=x(l)
1719  phot2(l)=y(l)
1720  ENDDO
1721 
1722  ym1=xm1
1723  ym2=xm2
1724 c to hepevt
1725  CALL filhep(0,1,id1,k,k,0,0,phot1,ym1,.true.)
1726  CALL filhep(0,1,id2,k,k,0,0,phot2,ym2,.true.)
1727 
1728 c
1729  ENDIF
1730  ENDDO
1731 
1732  END