C++ Interface to Tauola
tauola-F/jetset-F/tauface-jetset.f
1/* copyright(c) 1991-2023 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)
44C *************************************
45C general tauola interface, should work in every case until
46C hepevt is OK, does not check if hepevt is 'clean'
47C in particular will decay decayed taus...
48C only longitudinal spin effects are included.
49C in W decay v-a vertex is assumed
50C date: 12 DEC 1998. date: 21 June 1999. date: 24 Jan 2001 date: 24 Aug 2001
51C 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
85C 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
93C store decay vertexes
94 DIMENSION IMOTHER (20)
95 INTEGER KFHIGGS(3)
96
97C store daughter pointers
98 INTEGER ISON
99 COMMON /ISONS_TAU/ISON(2)
100 SAVE /ISONS_TAU/
101
102.EQ. IF(MODE-1) THEN
103C ***********************
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
121C 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)
128C 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
134C 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
142C ***********************
143C
144C..... find tau-s and fill common block /TAUPOS/
145C this is to avoid LUND history fillings. This call is optional
146 CALL PHYFIX(NSTOP,NSTART)
147C 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
155C ... and to find mothers giving taus.
156 NDEC = 0
157C(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
182C(BPK)--<
183 NDEC=NDEC+1
184 IMOTHER(NDEC)= JMOTH
185 ENDIF
186 9999 CONTINUE
187 ENDDO
188
189C ... 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
197C(BPK)-->
198C 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
223C(BPK)--<
224
225
226C ... we correct HEPEVT (fix developped with Catherine BISCARAT)
227.EQ.c IF (JDAHEP(2,IM)0) THEN ! ID of second daughter was missing
228c ISECU=1
229c 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
231c JDAHEP(2,IM)=I
232.EQ..AND..NE.c ELSEIF (JMOHEP(1,I)IMISECU1) THEN ! we have found one after there
233c JDAHEP(2,IM)=0 ! was something else, lets kill game
234c ENDIF
235.NE.c IF (JMOHEP(1,I)IM) ISECU=0 ! other stuff starts
236c ENDDO
237c ENDIF
238
239C ... 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
245C ... if there will be more we will come here again
246 666 CONTINUE
247
248C(BPK)-->
249 DO I=MAX(NP1+1,ISON(1)),ISON(2)
250C(BPK)--<
251.EQ..OR..EQ. IF(IDHEP(I)-KFTAUIDHEP(I)-KFNUE) NP1=I
252 ENDDO
253C(BPK)-->
254 DO I=MAX(NP2+1,ISON(1)),ISON(2)
255C(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
268c.....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
281C there is no angular dependence in gamma/Z polarization
282C there is no s-dependence in gamma/Z polarization at all
283C there is even no Z polarization in any form
284C main reason is that nobody asked ...
285C but it is prepared and longitudinal correlations
286C 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
312c.....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
342C 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
350C 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
361C ***********************
362C
363 CALL DEXAY(100,POL1)
364 CALL DEKAY(100,POL1x)
365 WRITE(IOUT,7002)
366 ENDIF
367C *****
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)
457C IM0 NP1 NP2 are the positions of Z/gamma tau tau in hepevt common block.
458C the purpose of this routine is to calculate polarization of Z along tau direction.
459C this is highly non-trivial due to necessity of reading infromation from hard process
460C 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
467C 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
489C(BPK)--> BROTHERS OF TAU ALREADY FOUND
490 INTEGER ISON
491 COMMON /ISONS_TAU/ISON(2)
492C(BPK)--<
493C >>
494C >> STEP 1: find where are particles in hepevent and pick them up
495C >>
496 HOPE=HOPEin
497C sometimes shade Z of Z is its mother ...
498 IM=IM0
499 IM00=JMOHEP(1,IM0)
500C 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
504C
505C find (host generator-level) incoming beam-bare-particles which form Z and co.
506 IMO1=JMOHEP(1,IM)
507 IMO2=JMOHEP(2,IM)
508
509C(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
515C(BPK)--<
516
517 IFFULL=0
518C 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)
521C(BPK)-->
522.EQ. IF (IDHEP(IMO1)IDHEP(NP1)) IMO1=JMOHEP(1,IMO1) ! PROTECT AGAINST COPIES
523C(BPK)--<
524 IMO2=JMOHEP(2,NP1)
525C(BPK)-->
526.EQ. IF (IDHEP(IMO2)IDHEP(NP2)) IMO2=JMOHEP(1,IMO2) ! PROTECT AGAINST COPIES
527C(BPK)--<
528 IFFULL=1
529C 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)
533C(BPK)-->
534.EQ. IF (IDHEP(IMO1)IDHEP(NP1)) IMO1=JMOHEP(1,IMO1) ! PROTECT AGAINST COPIES
535C(BPK)--<
536 IMO2=JMOHEP(2,NP1)
537C(BPK)-->
538.EQ. IF (IDHEP(IMO2)IDHEP(NP2)) IMO2=JMOHEP(1,IMO2) ! PROTECT AGAINST COPIES
539C(BPK)--<
540 IFFULL=1
541 ENDIF
542
543
544C 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
549C
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
555C corrections due to possible differences in 4-momentum of shadow vs true Z.
556C(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
563C(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
588C options which may be useful in some cases of two heavy boson production
589C need individual considerations. To be developed.
590C PLZAPX=PLZAP0(11,IDFQ1,SVAR,0D0) ! gamma/Z mixture as if produced from e beam
591C PLZAPX=PLZAP0(12,IDFQ1,SVAR,0D0) ! pure Z
592 PLZAPX=0.5 ! pure gamma
593 RETURN
594 ENDIF
595C >>
596C >> STEP 2 look for brothers of Z which have to be included in effective incoming particles
597C >>
598C let us define beginning and end of particles which are produced in parallel to Z
599C in principle following should work
600
601C(BPK)--> ACCOMMODATE FOR HERWIG - IM00 POINTS TO BEAM PARTICLE OR HARD CMS
602 NX1=JDAHEP(1,IM00)
603 NX2=JDAHEP(2,IM00)
604C 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
608C(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
631C case of annihilation of two bosons is hopeles
632.GE..AND..GE. IF (ABS(IDFP1)20ABS(IDFP2)20) HOPE=.FALSE.
633C 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
637C options which may be useful in some cases of two heavy boson production
638C need individual considerations. To be developed.
639C PLZAPX=PLZAP0(11,IDFQ1,SVAR,0D0) ! gamma/Z mixture as if produced from e beam
640C 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
648C >>
649C >> STEP 3 we combine gluons, photons into incoming effective beams
650C >>
651
652C in the following we will ignore the possibility of photon emission from taus
653C 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*--------------------------------------------------------------------------
668c IFLAV=min(ABS(IDFP1),ABS(IDFP2))
669c that means that always quark or lepton i.e. process like
670c f g(gamma) --> f Z0 --> tau tau
671c we glue fermions to effective beams that is f f --> Z0 --> tau tau
672c 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
697C 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
718C 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
740C now spectators ...
741
742C(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.
747C(BPK)--<
748
749 DO k=NX1,NX2
750.NE..AND..NE..AND. IF (ABS(IDHEP(K))IFLAVKIM
751C(BPK)-->
752.NE..AND..NE. $ KNPH1KNPH2) THEN
753C(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
821C >>
822C >> STEP 4 look for brothers of tau (sons of Z!) which have to be included in
823c >> effective outcoming taus
824C >>
825C let us define beginning and end of particles which are produced in
826c parallel to tau
827
828
829
830C find outcoming particles which come from Z
831
832
833
834
835C(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
839C(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
870C 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)
878C take effective beam which is less massive, it should be irrelevant
879C but in case HEPEVT is particulary dirty may help.
880C this routine calculate reduced system transver and cosine of scattering
881C 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
896C 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
913C 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)
922C this function calculates probability for the helicity +1 +1 configuration
923C 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
928C >>>>> 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)
941C ----------------------------------------------------------------------
942C THIS ROUTINE PROVIDES BORN CROSS SECTION. IT HAS THE SAME
943C STRUCTURE AS FUNTIS AND FUNTIH, THUS CAN BE USED AS SIMPLER
944C EXAMPLE OF THE METHOD APPLIED THERE
945C INPUT PARAMETERS ARE: SVAR -- transfer
946C COSTHE -- cosine of angle between tau+ and 1st beam
947C TA,TB -- helicity states of tau+ tau-
948C
949C called by : BORNY, BORAS, BORNV, WAGA, WEIGHT
950C ----------------------------------------------------------------------
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
961C=====================================================================
962 COMMON / T_GSWPRM /SWSQ,AMW,AMZ,AMH,AMTOP,GAMMZ
963 REAL*8 SWSQ,AMW,AMZ,AMH,AMTOP,GAMMZ
964C SWSQ = sin2 (theta Weinberg)
965C AMW,AMZ = W & Z boson masses respectively
966C AMH = the Higgs mass
967C AMTOP = the top mass
968C 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/
982C
983C MEMORIZATION =========================================================
984.NE..OR..NE..OR..NE. IF ( MODEMODE0SVARSVAR0COSTHECOST0
985.OR..NE. $ IDE0IDE)THEN
986C
987 KEYGSW=1
988C ** 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))
995C 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))
1000C 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
1020C
1021C******************
1022C* IN CALCULATING CROSS SECTION ONLY DIAGONAL ELEMENTS
1023C* OF THE SPIN DENSITY MATRICES ENTER (LONGITUD. POL. ONLY.)
1024C* 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
1037C MASS TERM IN BORN
1038.GE. IF (MODE1) THEN
1039 BORN=BORN+CDABS(ABORNM(I,J))**2*FACTOM
1040 ENDIF
1041
1042 150 CONTINUE
1043C************
1044 FUNT=BORN
1045.LT. IF(FUNT0.D0) FUNT=BORN
1046
1047C
1048.GT. IF (SVAR4D0*AMFIN**2) THEN
1049C 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
1056C 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
1074C SWSQ = sin2 (theta Weinberg)
1075C AMW,AMZ = W & Z boson masses respectively
1076C AMH = the Higgs mass
1077C AMTOP = the top mass
1078C GAMMZ = Z0 width
1079C
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
1162C ----------------------------------------------------------------------
1163C
1164C INITIALISATION OF COUPLING CONSTANTS AND FERMION-GAMMA / Z0 VERTEX
1165C
1166C called by : KORALZ
1167C ----------------------------------------------------------------------
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))
1186C
1187 NDIAG0=2
1188 NDIAGA=11
1189 KEYA = 1
1190 KEYZ = 1
1191C
1192C
1193 RETURN
1194 END
1195
1196 SUBROUTINE T_GIVIZO(IDFERM,IHELIC,SIZO3,CHARGE,KOLOR)
1197C ----------------------------------------------------------------------
1198C PROVIDES ELECTRIC CHARGE AND WEAK IZOSPIN OF A FAMILY FERMION
1199C IDFERM=1,2,3,4 DENOTES NEUTRINO, LEPTON, UP AND DOWN QUARK
1200C NEGATIVE IDFERM=-1,-2,-3,-4, DENOTES ANTIPARTICLE
1201C IHELIC=+1,-1 DENOTES RIGHT AND LEFT HANDEDNES ( CHIRALITY)
1202C SIZO3 IS THIRD PROJECTION OF WEAK IZOSPIN (PLUS MINUS HALF)
1203C AND CHARGE IS ELECTRIC CHARGE IN UNITS OF ELECTRON CHARGE
1204C KOLOR IS A QCD COLOUR, 1 FOR LEPTON, 3 FOR QUARKS
1205C
1206C called by : EVENTE, EVENTM, FUNTIH, .....
1207C ----------------------------------------------------------------------
1208 IMPLICIT REAL*8(A-H,O-Z)
1209C
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
1220C** NOTE THAT CONVENTIONALY Z0 COUPLING IS
1221C** 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/
1229C 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)
1242C ----------------------------------------------------------------------
1243C this subroutine fills one entry into the HEPEVT common
1244C and updates the information for affected mother entries
1245C
1246C written by Martin W. Gruenewald (91/01/28)
1247C
1248C called by : ZTOHEP,BTOHEP,DWLUxy
1249C ----------------------------------------------------------------------
1250C
1251C this is the hepevt class in old style. No d_h_ class pre-name
1252C 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
1274C
1275 REAL*4 P4(4)
1276C
1277C check address mode
1278.EQ. IF (N0) THEN
1279C
1280C append mode
1281 IHEP=NHEP+1
1282.GT. ELSE IF (N0) THEN
1283C
1284C absolute position
1285 IHEP=N
1286 ELSE
1287C
1288C relative position
1289 IHEP=NHEP+N
1290 END IF
1291C
1292C check on IHEP
1293.LE..OR..GT. IF ((IHEP0)(IHEPNMXHEP)) RETURN
1294C
1295C 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
1305C
1306 DO I=1,4
1307 PHEP(I,IHEP)=P4(I)
1308C
1309C 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
1313C FLAG FOR PHOTOS...
1314 QEDRAD(IHEP)=PHFLAG
1315C
1316C update process:
1317 DO IP=JMOHEP(1,IHEP),JMOHEP(2,IHEP)
1318.GT. IF(IP0)THEN
1319C
1320C if there is a daughter at IHEP, mother entry at IP has decayed
1321.EQ. IF(ISTHEP(IP)1)ISTHEP(IP)=2
1322C
1323C 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
1332C
1333 RETURN
1334 END
1335
1336
1337 FUNCTION IHEPDIM(DUM)
1338C this is the hepevt class in old style. No d_h_ class pre-name
1339C 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)
1373C no initialization required. Must be called once after every:
1374C 1) CALL DEKAY(1+10,...)
1375C 2) CALL DEKAY(2+10,...)
1376C 3) CALL DEXAY(1,...)
1377C 4) CALL DEXAY(2,...)
1378C subroutine to decay originating from TAUOLA's taus:
1379c 1) etas(with CALL taueta(jak))
1380c 2) later pi0's from taus.
1381C 3) extensions to other applications possible.
1382C this routine belongs to >tauola universal interface<, but uses
1383C routines from >tauola< utilities as well. 25.08.2005
1384C 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
1406C position of taus, must be defined by host program:
1407 COMMON /TAUPOS/ NP1,NP2
1408c
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)
1423C 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
1451c 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
1456c
1457 END
1458 SUBROUTINE taueta(JAK)
1459c subroutine to decay etas's from taus.
1460C this routine belongs to tauola universal interface, but uses
1461C routines from tauola utilities. Just flat phase space, but 4 channels.
1462C it is called at the beginning of SUBR. TAUPI0(JAK)
1463C and as far as hepevt search it is basically the same as TAUPI0. 25.08.2005
1464C 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
1493C position of taus, must be defined by host program:
1494 COMMON /TAUPOS/ NP1,NP2
1495c
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))
1503C 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
1515c 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
1538c 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))
1570c 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)
1586c 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)
1590c 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)
1593c 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)
1598c 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
1610c 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
1618c
1619 END
1620 SUBROUTINE tauk0s(JAK)
1621c subroutine to decay k0s's from taus.
1622C this routine belongs to tauola universal interface, but uses
1623C routines from tauola utilities. Just flat phase space, but 4 channels.
1624C it is called at the beginning of SUBR. TAUPI0(JAK)
1625C and as far as hepevt search it is basically the same as TAUPI0. 25.08.2005
1626C 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
1656C position of taus, must be defined by host program:
1657 COMMON /TAUPOS/ NP1,NP2
1658c
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))
1666C 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
1680c 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
1724c 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
1728c
1729 ENDIF
1730 ENDDO
1731
1732 END