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