1 /* copyright(c) 1991-2012 free software foundation, inc.
2 this file is part of the gnu c library.
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.
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.
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 <http://www.gnu.org/licenses/>. */
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. */
27 /* we
do support the iec 559 math functionality,
real and complex. */
29 /* wchar_t uses iso/iec 10646 (2nd ed., published 2011-03-15) /
32 /* we
do not support c11 <threads.h>. */
34 SUBROUTINE tauola(MODE,KEYPOL)
35 c *************************************
36 c general tauola interface, should work in every
case until
37 c hepevt is ok, does not check
if hepevt is
'clean'
38 c in particular will decay decayed taus...
39 c only longitudinal spin effects are included.
40 c in w decay v-a vertex is assumed
41 c date: 12 dec 1998. date: 21 june 1999. date: 24 jan 2001 date: 24 aug 2001
42 c this is the hepevt class in old style. no d_h_ class pre-name
44 parameter(nmxhep=10000)
46 INTEGER nevhep,nhep,isthep,idhep,jmohep,
57 * ----------------------------------------------------------------------
61 * ----------------------------------------------------------------------
63 COMMON /taupos/ np1, np2
64 REAL*4 phoi(4),phof(4)
65 double precision q1(4),q2(4),p1(4),p2(4),p3(4),p4(4)
66 COMMON / momdec / q1,q2,p1,p2,p3,p4
67 * tauola, photos and jetset overall switches
68 COMMON /libra/ jak1,jak2,itdkrc,ifphot,ifhadm,ifhadp
71 common /pseudocoup/ csc,ssc
74 COMMON / inout / inut,iout
76 c to switch tau polarization off in taus
77 dimension pol1(4), pol2(4)
78 double precision pol1x(4), pol2x(4)
80 DATA pol1 /0.0,0.0,0.0,0.0/
81 DATA pol2 /0.0,0.0,0.0,0.0/
82 DATA pi /3.141592653589793238462643d0/
84 c store decay vertexes
88 c store daughter pointers
90 COMMON /isons_tau/ison(2)
94 c ***********************
112 c couplings of the
'pseudoscalar higgs' as in cern-th/2003-166
116 betah=sqrt(1d0-4*xmtau**2/xmh**2)
119 c
write(*,*)
' scalar component=',csc,
' pseudo-scalar component=',ssc
120 IF (ifphot.EQ.1) CALL phoini
121 CALL inietc(jak1,jak2,itdkrc,ifphot)
125 c activation of pi0 and eta decays: (1) means on, (0) off
129 CALL taupi0(-1,1,ion)
131 WRITE(iout,7001) pol,psi,ion(1),ion(2),ion(3)
132 ELSEIF(mode.EQ.0)
THEN
133 c ***********************
135 c..... find tau-s and fill
common block /taupos/
136 c this is to avoid lund history fillings. this call is
optional
137 CALL phyfix(nstop,nstart)
138 c clear mothers of the previous event
146 c ... and to find mothers giving taus.
148 c(bpk)--> look for mother, check that it is not the history entry(e.g. mstp(128)=0)
150 IF(abs(idhep(i)).EQ.kftau.AND.isthep(i).EQ.1.AND.
151 $ (isthep(i).GE.125.OR.isthep(i).LT.120))
THEN
153 DO WHILE (abs(idhep(imoth)).EQ.kftau)
154 imoth=jmohep(1,imoth)
156 IF (isthep(imoth).EQ.3.OR.
157 $ (isthep(imoth).GE.120.AND.isthep(imoth).LE.125))
THEN
159 IF (idhep(j).EQ.idhep(imoth).AND.
160 $ jmohep(1,j).EQ.imoth.AND.
161 $ isthep(j).EQ.2)
THEN
171 IF(jmoth.EQ.imother(ii)) goto 9999
180 c ... taus of every mother are treated in this main loop
189 c correcting hepevt is out of question at this point..
191 IF (idhep(jmohep(1,im0)).EQ.idhep(im0)) im0=jmohep(1,im0)
194 IF (isthep(i).EQ.3.OR.
195 $ (isthep(i).GE.120.AND.isthep(i).LE.125))
THEN
199 DO WHILE (idhep(imoth).EQ.idhep(i).OR.abs(idhep(imoth)).EQ.kftau)
200 imoth=jmohep(1,imoth)
202 IF ((imoth.EQ.im0.OR.imoth.EQ.im).AND.isel.EQ.-1)
THEN
205 ELSEIF ((imoth.EQ.im0.OR.imoth.EQ.im).AND.isel.EQ.0)
THEN
207 ELSEIF ((imoth.NE.im0.AND.imoth.NE.im).AND.isel.EQ.0)
THEN
217 c ... we correct hepevt(fix developped with catherine biscarat)
218 c
IF (jdahep(2,im).EQ.0)
THEN
220 c
DO i=jdahep(1,im)+1,nhep
221 c
IF (jmohep(1,i).EQ.im.AND.isecu.EQ.1)
THEN
223 c
ELSEIF (jmohep(1,i).EQ.im.AND.isecu.NE.1)
THEN
226 c
IF (jmohep(1,i).NE.im) isecu=0
230 c ... we check whether there are just two or more tau-likes
232 IF(idhep(i).EQ.-kftau.OR.idhep(i).EQ.-kfnue) ncount=ncount+1
233 IF(idhep(i).EQ. kftau.OR.idhep(i).EQ. kfnue) ncount=ncount+1
236 c ...
if there will be more we will come here again
240 DO i=max(np1+1,ison(1)),ison(2)
242 IF(idhep(i).EQ.-kftau.OR.idhep(i).EQ.-kfnue) np1=i
245 DO i=max(np2+1,ison(1)),ison(2)
247 IF(idhep(i).EQ. kftau.OR.idhep(i).EQ. kfnue) np2=i
259 c.....include polarisation effect
262 IF(idhep(im).EQ.kfhiggs(1).OR.idhep(im).EQ.kfhiggs(2).OR.
263 $ idhep(im).EQ.kfhiggs(3))
THEN
264 IF(rrr(1).LT.0.5)
THEN
271 ELSEIF((idhep(im).EQ.kfz0).OR.(idhep(im).EQ.kfgam))
THEN
272 c there is no angular dependence in gamma/z polarization
273 c there is no s-dependence in gamma/z polarization at all
274 c there is even no z polarization in any form
275 c main reason is that nobody asked ...
276 c but it is prepared and longitudinal correlations
277 c can be included up to koralz standards
279 polz0=plzapx(.true.,im,np1,np2)
280 IF(rrr(1).LT.polz0)
THEN
287 ELSEIF(idhep(np1).EQ.-idhep(np2))
THEN
288 polz0=plzapx(.true.,im,np1,np2)
289 IF(rrr(1).LT.polz0)
THEN
296 ELSEIF(abs(idhep(im)).EQ.kfhigch)
THEN
303 c.....include polarisation effect
306 IF(idhep(im).EQ.kfhiggs(1).OR.idhep(im).EQ.kfhiggs(2).OR.
307 $ idhep(im).EQ.kfhiggs(3))
THEN
308 IF(idhep(np1).EQ.-kftau .AND.
309 $ (jdahep(1,np1).LE.np1.OR.jdahep(1,np1).GT.nhep) .AND.
310 $ idhep(np2).EQ. kftau .AND.
311 $ (jdahep(1,np2).LE.np2.OR.jdahep(1,np2).GT.nhep)
313 IF (idhep(im).EQ.kfhiggs(1))
THEN
315 ELSEIF (idhep(im).EQ.kfhiggs(2))
THEN
317 ELSEIF (idhep(im).EQ.kfhiggs(3))
THEN
320 WRITE(*,*)
'Warning from TAUOLA:'
321 WRITE(*,*)
'I stop this run, wrong IDHEP(IM)=',idhep(im)
324 CALL spinhiggs(im,np1,np2,ifpseudo,pol1,pol2)
325 IF (ifphot.EQ.1) CALL photos(im)
331 IF(idhep(np1).EQ.-kftau.AND.
332 $ (jdahep(1,np1).LE.np1.OR.jdahep(1,np1).GT.nhep))
THEN
333 c here check on
if np1 was not decayed should be verified
335 IF (ifphot.EQ.1) CALL photos(np1)
339 IF(idhep(np2).EQ. kftau.AND.
340 $ (jdahep(1,np2).LE.np2.OR.jdahep(1,np2).GT.nhep))
THEN
341 c here check on
if np2 was not decayed should be added
343 IF (ifphot.EQ.1) CALL photos(np2)
348 IF (ncount.GT.0) goto 666
351 ELSEIF(mode.EQ.1)
THEN
352 c ***********************
355 CALL dekay(100,pol1x)
359 7001
FORMAT(///1x,15(5h*****)
360 $ /,
' *', 25x,
'*****TAUOLA UNIVERSAL INTERFACE: ******',9x,1h*,
361 $ /,
' *', 25x,
'*****VERSION 1.22, April 2009 (gfort)**',9x,1h*,
362 $ /,
' *', 25x,
'**AUTHORS: P. Golonka, B. Kersevan, ***',9x,1h*,
363 $ /,
' *', 25x,
'**T. Pierzchala, E. Richter-Was, ******',9x,1h*,
364 $ /,
' *', 25x,
'****** Z. Was, M. Worek ***************',9x,1h*,
365 $ /,
' *', 25x,
'**USEFUL DISCUSSIONS, IN PARTICULAR ***',9x,1h*,
366 $ /,
' *', 25x,
'*WITH C. Biscarat and S. Slabospitzky**',9x,1h*,
367 $ /,
' *', 25x,
'****** are warmly acknowledged ********',9x,1h*,
368 $ /,
' *', 25x,
' ',9x,1h*,
369 $ /,
' *', 25x,
'********** INITIALIZATION ************',9x,1h*,
370 $ /,
' *',f20.5,5x,
'tau polarization switch must be 1 or 0 ',9x,1h*,
371 $ /,
' *',f20.5,5x,
'Higs scalar/pseudo mix CERN-TH/2003-166',9x,1h*,
372 $ /,
' *',i10, 15x,
'PI0 decay switch must be 1 or 0 ',9x,1h*,
373 $ /,
' *',i10, 15x,
'ETA decay switch must be 1 or 0 ',9x,1h*,
374 $ /,
' *',i10, 15x,
'K0S decay switch must be 1 or 0 ',9x,1h*,
377 7002
FORMAT(///1x,15(5h*****)
378 $ /,
' *', 25x,
'*****TAUOLA UNIVERSAL INTERFACE: ******',9x,1h*,
379 $ /,
' *', 25x,
'*****VERSION 1.22, April 2009 (gfort)**',9x,1h*,
380 $ /,
' *', 25x,
'**AUTHORS: P. Golonka, B. Kersevan, ***',9x,1h*,
381 $ /,
' *', 25x,
'**T. Pierzchala, E. Richter-Was, ******',9x,1h*,
382 $ /,
' *', 25x,
'****** Z. Was, M. Worek ***************',9x,1h*,
383 $ /,
' *', 25x,
'**USEFUL DISCUSSIONS, IN PARTICULAR ***',9x,1h*,
384 $ /,
' *', 25x,
'*WITH C. Biscarat and S. Slabospitzky**',9x,1h*,
385 $ /,
' *', 25x,
'****** are warmly acknowledged ********',9x,1h*,
386 $ /,
' *', 25x,
'****** END OF MODULE OPERATION ********',9x,1h*,
391 SUBROUTINE spinhiggs(IM,NP1,NP2,IFPSEUDO,Pol1,Pol2)
393 REAL*8 hh1,hh2,wthiggs
394 dimension pol1(4), pol2(4),hh1(4),hh2(4), rrr(1)
402 wt=wthiggs(ifpseudo,hh1,hh2)
403 IF (rrr(1).GT.wt) goto 10
409 FUNCTION wthiggs(IFPSEUDO,HH1,HH2)
411 common /pseudocoup/ csc,ssc
414 REAL*8 hh1(4),hh2(4),r(4,4),wthiggs
428 r(1,1)=(csc**2-ssc**2)/(csc**2+ssc**2)
429 r(2,2)=(csc**2-ssc**2)/(csc**2+ssc**2)
430 r(1,2)=2*csc*ssc/(csc**2+ssc**2)
431 r(2,1)=-2*csc*ssc/(csc**2+ssc**2)
441 wthiggs=wthiggs+r(k,l)*hh1(k)*hh2(l)
447 FUNCTION plzapx(HOPEin,IM0,NP1,NP2)
448 c im0 np1 np2 are the positions of z/gamma tau tau in hepevt
common block.
449 c the purpose of this routine is to calculate polarization of z along tau direction.
450 c this is highly non-trivial due to necessity of reading infromation from hard process
451 c history in hepevt, which is often written not up to the gramatic rules.
452 REAL*8 plzap0,svar,costhe,sini,sfin,zprop2,
453 $ p1(4),p2(4),q1(4),q2(4),qq(4),ph(4),pd1(4),pd2(4),
454 $ pq1(4),pq2(4),pb(4),pa(4)
458 c this is the hepevt class in old style. no d_h_ class pre-name
460 parameter(nmxhep=10000)
462 INTEGER nevhep,nhep,isthep,idhep,jmohep,
473 * ----------------------------------------------------------------------
477 * ----------------------------------------------------------------------
480 c(bpk)--> brothers of tau already found
482 COMMON /isons_tau/ison(2)
485 c >> step 1: find
where are particles in hepevent and pick them up
488 c sometimes shade z of z is its mother ...
491 c to protect against check on mother of beam particles.
493 IF (idhep(im0).EQ.idhep(im00)) im=jmohep(1,im0)
496 c find(host generator-level) incoming beam-bare-particles which form z and co.
500 c(bpk)--> in herwig the
POINTER might be to hard cms
502 IF (isthep(im00).EQ.120)
THEN
509 c
case when it was like e+e- --> tau+tau- gammas and e+e- were
'first' in hepevt.
510 IF (imo1.EQ.0.AND.imo2.EQ.0)
THEN
513 IF (idhep(imo1).EQ.idhep(np1)) imo1=jmohep(1,imo1)
517 IF (idhep(imo2).EQ.idhep(np2)) imo2=jmohep(1,imo2)
520 c
case when it was like qq --> tau+tau- gammas and qq were not
'first' in hepevt.
522 ELSEIF (idhep(im).NE.22.AND.idhep(im).NE.23)
THEN
525 IF (idhep(imo1).EQ.idhep(np1)) imo1=jmohep(1,imo1)
529 IF (idhep(imo2).EQ.idhep(np2)) imo2=jmohep(1,imo2)
535 c and check
if it really happened
536 IF (imo1.EQ.0) hope=.false.
537 IF (imo2.EQ.0) hope=.false.
538 IF (imo1.EQ.imo2) hope=.false.
546 c corrections due to possible differences in 4-momentum of shadow vs true z.
548 IF (im.EQ.jmohep(1,im0).AND.
549 $ (idhep(im).EQ.22.OR.idhep(im).EQ.23))
THEN
556 CALL bostdq( 1,pa, q1, q1)
557 CALL bostdq( 1,pa, q2, q2)
558 CALL bostdq(-1,pb, q1, q1)
559 CALL bostdq(-1,pb, q2, q2)
565 IF (hope) p1(i)=phep(i,imo1)
566 IF (hope) p2(i)=phep(i,imo2)
574 IF (hope) idfp1=idhep(imo1)
575 IF (hope) idfp2=idhep(imo2)
577 svar=qq(4)**2-qq(3)**2-qq(2)**2-qq(1)**2
579 c options which may be useful in some cases of two heavy boson production
580 c need individual considerations. to be developed.
581 c plzapx=plzap0(11,idfq1,svar,0d0)
582 c plzapx=plzap0(12,idfq1,svar,0d0)
587 c >> step 2 look for brothers of z which have to be included in effective incoming particles
589 c let us define beginning and
end of particles which are produced in parallel to z
590 c in principle following should work
592 c(bpk)--> accommodate for herwig - im00 points to beam particle or hard cms
597 IF (iffull.EQ.1) inbr=np1
598 IF (idhep(jmohep(1,inbr)).EQ.idhep(inbr)) inbr=jmohep(1,inbr)
600 IF(nx1.EQ.0.OR.nx2.EQ.0)
THEN
604 IF(jmohep(1,inbr-k).EQ.jmohep(1,inbr))
THEN
613 IF(jmohep(1,k).EQ.jmohep(1,inbr))
THEN
622 c
case of annihilation of two bosons is hopeles
623 IF (abs(idfp1).GE.20.AND.abs(idfp2).GE.20) hope=.false.
624 c
case of annihilation of non-matching flavors is hopeless
625 IF (abs(idfp1).LE.20.AND.abs(idfp2).LE.20.AND.idfp1+idfp2.NE.0)
628 c options which may be useful in some cases of two heavy boson production
629 c need individual considerations. to be developed.
630 c plzapx=plzap0(11,idfq1,svar,0d0)
631 c plzapx=plzap0(12,idfq1,svar,0d0)
635 IF (abs(idfp1).LT.20) ide= idfp1
636 IF (abs(idfp2).LT.20) ide=-idfp2
640 c >> step 3 we combine gluons, photons into incoming effective beams
643 c in the following we will ignore the possibility of photon emission from taus
644 c however at certain moment it will be necessary to take care of
656 iflav=min(abs(idfp1),abs(idfp2))
658 *--------------------------------------------------------------------------
659 c iflav=min(abs(idfp1),abs(idfp2))
660 c that means that always quark or lepton i.e. process like
661 c f g(gamma) --> f z0 --> tau tau
662 c we glue fermions to effective beams that is f f --> z0 --> tau tau
663 c with gamma/g emission from initial fermion.
664 *---------------------------------------------------------------------------
666 IF (abs(idfp1).GE.20)
THEN
669 IF (abs(idp).EQ.iflav)
THEN
677 IF (abs(idfp2).GE.20)
THEN
680 IF (abs(idp).EQ.iflav)
THEN
688 c
if first beam was boson: gluing
690 IF (abs(idfp1).GE.20)
THEN
694 xm1=abs((pd1(4)+ph(4))**2-(pd1(3)+ph(3))**2
695 $ -(pd1(2)+ph(2))**2-(pd1(1)+ph(1))**2)
696 xm2=abs((pd2(4)+ph(4))**2-(pd2(3)+ph(3))**2
697 $ -(pd2(2)+ph(2))**2-(pd2(1)+ph(1))**2)
709 c
if second beam was boson: gluing
712 IF (abs(idfp2).GE.20)
THEN
716 xm1=abs((pd1(4)+ph(4))**2-(pd1(3)+ph(3))**2
717 $ -(pd1(2)+ph(2))**2-(pd1(1)+ph(1))**2)
718 xm2=abs((pd2(4)+ph(4))**2-(pd2(3)+ph(3))**2
719 $ -(pd2(2)+ph(2))**2-(pd2(1)+ph(1))**2)
736 IF (idhep(jmohep(1,np1)).EQ.idhep(np1)) nph1=jmohep(1,np1)
737 IF (idhep(jmohep(1,np2)).EQ.idhep(np2)) nph2=jmohep(1,np2)
741 IF (abs(idhep(k)).NE.iflav.AND.k.NE.im.AND.
743 $ k.NE.nph1.AND.k.NE.nph2)
THEN
745 IF(idhep(k).EQ.22.AND.iffull.EQ.1)
THEN
749 xm1=abs((pd1(4)-ph(4))**2-(pd1(3)-ph(3))**2
750 $ -(pd1(2)-ph(2))**2-(pd1(1)-ph(1))**2)
751 xm2=abs((pd2(4)-ph(4))**2-(pd2(3)-ph(3))**2
752 $ -(pd2(2)-ph(2))**2-(pd2(1)-ph(1))**2)
753 xm3=abs((pq1(4)+ph(4))**2-(pq1(3)+ph(3))**2
754 $ -(pq1(2)+ph(2))**2-(pq1(1)+ph(1))**2)
755 xm4=abs((pq2(4)+ph(4))**2-(pq2(3)+ph(3))**2
756 $ -(pq2(2)+ph(2))**2-(pq2(1)+ph(1))**2)
759 sini=abs((pd1(4)+pd2(4)-ph(4))**2-(pd1(3)+pd2(3)-ph(3))**2
760 $ -(pd1(2)+pd2(2)-ph(2))**2-(pd1(1)+pd2(1)-ph(1))**2)
761 sfin=abs((pd1(4)+pd2(4) )**2-(pd1(3)+pd2(3) )**2
762 $ -(pd1(2)+pd2(2) )**2-(pd1(1)+pd2(1) )**2)
772 xm=min(xm1,xm2,xm3,xm4)
777 ELSEIF (xm2.EQ.xm)
THEN
781 ELSEIF (xm3.EQ.xm)
THEN
794 xm1=abs((pd1(4)-ph(4))**2-(pd1(3)-ph(3))**2
795 $ -(pd1(2)-ph(2))**2-(pd1(1)-ph(1))**2)
796 xm2=abs((pd2(4)-ph(4))**2-(pd2(3)-ph(3))**2
797 $ -(pd2(2)-ph(2))**2-(pd2(1)-ph(1))**2)
813 c >> step 4 look for brothers of tau(sons of z
814 c >> effective outcoming taus
816 c let us define beginning and
end of particles which are produced in
821 c find outcoming particles which come from z
826 c(bpk)--> ok, it would have to be along taus in hard record with the same mother
827 IF (abs(idhep(im0)).EQ.22.OR.abs(idhep(im0)).EQ.23)
THEN
829 IF(abs(idhep(k)).EQ.22)
THEN
836 xm3=abs((pq1(4)+ph(4))**2-(pq1(3)+ph(3))**2
837 $ -(pq1(2)+ph(2))**2-(pq1(1)+ph(1))**2)
838 xm4=abs((pq2(4)+ph(4))**2-(pq2(3)+ph(3))**2
839 $ -(pq2(2)+ph(2))**2-(pq2(1)+ph(1))**2)
858 *------------------------------------------------------------------------
861 c out of effective momenta we calculate costhe and later polarization
862 CALL angulu(pd1,pd2,q1,q2,costhe)
864 plzapx=plzap0(ide,idfq1,svar,costhe)
867 SUBROUTINE angulu(PD1,PD2,Q1,Q2,COSTHE)
868 REAL*8 pd1(4),pd2(4),q1(4),q2(4),costhe,p(4),qq(4),qt(4)
869 c take effective beam which is less massive, it should be irrelevant
870 c but in
case hepevt is particulary dirty may help.
871 c this routine calculate reduced system transver and cosine of scattering
874 xm1=abs(pd1(4)**2-pd1(3)**2-pd1(2)**2-pd1(1)**2)
875 xm2=abs(pd2(4)**2-pd2(3)**2-pd2(2)**2-pd2(1)**2)
887 c calculate space like part of p(in z restframe)
893 xmqq=sqrt(qq(4)**2-qq(3)**2-qq(2)**2-qq(1)**2)
895 qtxqq=qt(4)*qq(4)-qt(3)*qq(3)-qt(2)*qq(2)-qt(1)*qq(1)
897 qt(k)=qt(k)-qq(k)*qtxqq/xmqq**2
900 pxqq=p(4)*qq(4)-p(3)*qq(3)-p(2)*qq(2)-p(1)*qq(1)
902 p(k)=p(k)-qq(k)*pxqq/xmqq**2
905 pxp =sqrt(p(1)**2+p(2)**2+p(3)**2-p(4)**2)
906 qtxqt=sqrt(qt(3)**2+qt(2)**2+qt(1)**2-qt(4)**2)
907 pxqt =p(3)*qt(3)+p(2)*qt(2)+p(1)*qt(1)-p(4)*qt(4)
908 costhe=pxqt/pxp/qtxqt
912 FUNCTION plzap0(IDE,IDF,SVAR,COSTH0)
913 c this
function calculates probability for the helicity +1 +1 configuration
914 c of taus for given z/gamma transfer and costh0 cosine of scattering angle
915 REAL*8 plzap0,svar,costhe,costh0,t_born
918 c >>>>>
IF (ide*idf.LT.0) costhe=-costh0
919 c >>>>> of first beam is used by t_giviz0 including sign
922 CALL initwk(ide,idf,svar)
924 CALL initwk(-ide,-idf,svar)
926 plzap0=t_born(0,svar,costhe,1d0,1d0)
927 $ /(t_born(0,svar,costhe,1d0,1d0)+t_born(0,svar,costhe,-1d0,-1d0))
931 FUNCTION t_born(MODE,SVAR,COSTHE,TA,TB)
932 c ----------------------------------------------------------------------
933 c this routine provides born cross section. it has the same
934 c structure as funtis and funtih, thus can be used as simpler
935 c example of the method applied there
936 c input parameters are: svar -- transfer
937 c costhe -- cosine of angle between tau+ and 1st beam
938 c ta,tb -- helicity states of tau+ tau-
940 c called by : borny, boras, bornv, waga, weight
941 c ----------------------------------------------------------------------
942 IMPLICIT REAL*8(a-h,o-z)
943 COMMON / t_beampm / ene ,amin,amfin,ide,idf
944 REAL*8 ene ,amin,amfin
945 COMMON / t_gauspm /ss,poln,t3e,qe,t3f,qf
946 & ,xupgi ,xupzi ,xupgf ,xupzf
947 & ,ndiag0,ndiaga,keya,keyz
948 & ,itce,jtce,itcf,jtcf,kolor
949 REAL*8 ss,poln,t3e,qe,t3f,qf
950 & ,xupgi(2),xupzi(2),xupgf(2),xupzf(2)
952 c=====================================================================
953 COMMON / t_gswprm /swsq,amw,amz,amh,amtop,gammz
954 REAL*8 swsq,amw,amz,amh,amtop,gammz
955 c swsq = sin2(theta weinberg)
956 c amw,amz = w & z boson masses respectively
957 c amh = the higgs mass
958 c amtop = the top mass
960 COMPLEX*16 aborn(2,2),aphot(2,2),azett(2,2)
961 COMPLEX*16 xupzfp(2),xupzip(2)
962 COMPLEX*16 abornm(2,2),aphotm(2,2),azettm(2,2)
963 COMPLEX*16 propa,propz
965 COMPLEX*16 xupf,xupi,xff(4),xfem,xfota,xrho,xke,xkf,xkef
966 COMPLEX*16 xthing,xve,xvf,xvef
967 DATA xi/(0.d0,1.d0)/,xr/(1.d0,0.d0)/
970 DATA svar0,cost0 /-5.d0,-6.d0/
971 DATA pi /3.141592653589793238462643d0/
972 DATA seps1,seps2 /0d0,0d0/
974 c memorization =========================================================
975 IF ( mode.NE.mode0.OR.svar.NE.svar0.OR.costhe.NE.cost0
976 $ .OR.ide0.NE.ide)
THEN
984 sinthe=sqrt(1.d0-costhe**2)
985 beta=sqrt(max(0d0,1d0-4d0*amfin**2/svar))
986 c i multiply axial coupling by beta factor.
987 xupzfp(1)=0.5d0*(xupzf(1)+xupzf(2))+0.5*beta*(xupzf(1)-xupzf(2))
988 xupzfp(2)=0.5d0*(xupzf(1)+xupzf(2))-0.5*beta*(xupzf(1)-xupzf(2))
989 xupzip(1)=0.5d0*(xupzi(1)+xupzi(2))+0.5*(xupzi(1)-xupzi(2))
990 xupzip(2)=0.5d0*(xupzi(1)+xupzi(2))-0.5*(xupzi(1)-xupzi(2))
991 c final state vector coupling
992 xupf =0.5d0*(xupzf(1)+xupzf(2))
993 xupi =0.5d0*(xupzi(1)+xupzi(2))
997 propz =1d0/dcmplx(svar-amz**2,svar/amz*gammz)
998 IF (keygsw.EQ.0) propz=0.d0
1001 regula= (3-2*i)*(3-2*j) + costhe
1002 regulm=-(3-2*i)*(3-2*j) * sinthe *2.d0*amfin/sqrt(svar)
1003 aphot(i,j)=propa*(xupgi(i)*xupgf(j)*regula)
1004 azett(i,j)=propz*(xupzip(i)*xupzfp(j)+xthing)*regula
1005 aborn(i,j)=aphot(i,j)+azett(i,j)
1006 aphotm(i,j)=propa*dcmplx(0d0,1d0)*xupgi(i)*xupgf(j)*regulm
1007 azettm(i,j)=propz*dcmplx(0d0,1d0)*(xupzip(i)*xupf+xthing)*regulm
1008 abornm(i,j)=aphotm(i,j)+azettm(i,j)
1013 c* in calculating cross section only diagonal elements
1014 c* of the spin density matrices enter(longitud. pol. only.)
1015 c* helicity conservation explicitly obeyed
1023 factor=kolor*(1d0+helic*polar1)*(1d0-helic*polar2)/4d0
1024 factom=factor*(1+helit*ta)*(1-helit*tb)
1025 factor=factor*(1+helit*ta)*(1+helit*tb)
1027 born=born+cdabs(aborn(i,j))**2*factor
1030 born=born+cdabs(abornm(i,j))**2*factom
1036 IF(funt.LT.0.d0) funt=born
1039 IF (svar.GT.4d0*amfin**2)
THEN
1040 c phase space threshold factor
1041 thresh=sqrt(1-4d0*amfin**2/svar)
1042 t_born= funt*svar**2*thresh
1047 c zw here was an error 19. 05. 1989
1052 SUBROUTINE initwk(IDEX,IDFX,SVAR)
1054 IMPLICIT REAL*8 (a-h,o-z)
1055 COMMON / t_beampm / ene ,amin,amfin,ide,idf
1056 REAL*8 ene ,amin,amfin
1057 COMMON / t_gauspm /ss,poln,t3e,qe,t3f,qf
1058 & ,xupgi ,xupzi ,xupgf ,xupzf
1059 & ,ndiag0,ndiaga,keya,keyz
1060 & ,itce,jtce,itcf,jtcf,kolor
1061 REAL*8 ss,poln,t3e,qe,t3f,qf
1062 & ,xupgi(2),xupzi(2),xupgf(2),xupzf(2)
1063 COMMON / t_gswprm /swsq,amw,amz,amh,amtop,gammz
1064 REAL*8 swsq,amw,amz,amh,amtop,gammz
1065 c swsq = sin2(theta weinberg)
1066 c amw,amz = w & z boson masses respectively
1067 c amh = the higgs mass
1068 c amtop = the top mass
1076 IF (idfx.EQ. 15)
then
1079 ELSEIF (idfx.EQ.-15)
then
1083 WRITE(*,*)
'INITWK: WRONG IDFX'
1087 IF (idex.EQ. 11)
then
1090 ELSEIF (idex.EQ.-11)
then
1093 ELSEIF (idex.EQ. 13)
then
1096 ELSEIF (idex.EQ.-13)
then
1099 ELSEIF (idex.EQ. 1)
then
1102 ELSEIF (idex.EQ.- 1)
then
1105 ELSEIF (idex.EQ. 2)
then
1108 ELSEIF (idex.EQ.- 2)
then
1111 ELSEIF (idex.EQ. 3)
then
1114 ELSEIF (idex.EQ.- 3)
then
1117 ELSEIF (idex.EQ. 4)
then
1120 ELSEIF (idex.EQ.- 4)
then
1123 ELSEIF (idex.EQ. 5)
then
1126 ELSEIF (idex.EQ.- 5)
then
1129 ELSEIF (idex.EQ. 12)
then
1132 ELSEIF (idex.EQ.- 12)
then
1135 ELSEIF (idex.EQ. 14)
then
1138 ELSEIF (idex.EQ.- 14)
then
1141 ELSEIF (idex.EQ. 16)
then
1144 ELSEIF (idex.EQ.- 16)
then
1149 WRITE(*,*)
'INITWK: WRONG IDEX'
1153 c ----------------------------------------------------------------------
1155 c initialisation of coupling constants and fermion-gamma / z0 vertex
1157 c called by : koralz
1158 c ----------------------------------------------------------------------
1163 CALL t_givizo( ide, 1,aizor,qe,kdumm)
1164 CALL t_givizo( ide,-1,aizol,qe,kdumm)
1168 xupzi(1)=(aizor-qe*swsq)/sqrt(swsq*(1-swsq))
1169 xupzi(2)=(aizol-qe*swsq)/sqrt(swsq*(1-swsq))
1170 CALL t_givizo( idf, 1,aizor,qf,kolor)
1171 CALL t_givizo( idf,-1,aizol,qf,kolor)
1175 xupzf(1)=(aizor-qf*swsq)/sqrt(swsq*(1-swsq))
1176 xupzf(2)=(aizol-qf*swsq)/sqrt(swsq*(1-swsq))
1187 SUBROUTINE t_givizo(IDFERM,IHELIC,SIZO3,CHARGE,KOLOR)
1188 c ----------------------------------------------------------------------
1189 c provides electric charge and weak izospin of a family fermion
1190 c idferm=1,2,3,4 denotes neutrino, lepton, up and down quark
1191 c negative idferm=-1,-2,-3,-4, denotes antiparticle
1192 c ihelic=+1,-1 denotes right and left handednes( chirality)
1193 c sizo3 is third projection of weak izospin(plus minus half)
1194 c and charge is electric charge in units of electron charge
1195 c kolor is a qcd colour, 1 for lepton, 3 for quarks
1197 c called by : evente, eventm, funtih, .....
1198 c ----------------------------------------------------------------------
1199 IMPLICIT REAL*8(a-h,o-z)
1201 IF(idferm.EQ.0.OR.iabs(idferm).GT.4) goto 901
1202 IF(iabs(ihelic).NE.1) goto 901
1204 idtype =iabs(idferm)
1206 lepqua=int(idtype*0.4999999d0)
1207 iupdow=idtype-2*lepqua-1
1208 charge =(-iupdow+2d0/3d0*lepqua)*ic
1209 sizo3 =0.25d0*(ic-ih)*(1-2*iupdow)
1211 c** note that conventionaly z0 coupling is
1212 c** xoupz=(sizo3-charge*swsq)/sqrt(swsq*(1-swsq))
1214 901 print *,
' STOP IN GIVIZO: WRONG PARAMS.'
1217 SUBROUTINE phyfix(NSTOP,NSTART)
1218 common/lujets/n,k(4000,5),p(4000,5),v(4000,5)
1220 c nstop nstart : when phytia history ends and event starts.
1224 IF(k(i,1).NE.21)
THEN
1232 SUBROUTINE filhep(N,IST,ID,JMO1,JMO2,JDA1,JDA2,P4,PINV,PHFLAG)
1233 c ----------------------------------------------------------------------
1234 c this
subroutine fills one entry into the hepevt common
1235 c and updates the information for affected mother entries
1237 c written by martin w. gruenewald(91/01/28)
1239 c called by : ztohep,btohep,dwluxy
1240 c ----------------------------------------------------------------------
1242 c this is the hepevt class in old style. no d_h_ class pre-name
1243 c this is the hepevt class in old style. no d_h_ class pre-name
1245 parameter(nmxhep=10000)
1247 INTEGER nevhep,nhep,isthep,idhep,jmohep,
1258 * ----------------------------------------------------------------------
1262 * ----------------------------------------------------------------------
1268 c check address mode
1273 ELSE IF (n.GT.0)
THEN
1284 IF ((ihep.LE.0).OR.(ihep.GT.nmxhep))
RETURN
1291 IF(jmo1.LT.0)jmohep(1,ihep)=jmohep(1,ihep)+ihep
1293 IF(jmo2.LT.0)jmohep(2,ihep)=jmohep(2,ihep)+ihep
1300 c koral-b and koral-z
do not provide vertex and/or lifetime informations
1304 c flag for photos...
1308 DO ip=jmohep(1,ihep),jmohep(2,ihep)
1311 c
if there is a daughter at ihep, mother entry at ip has decayed
1312 IF(isthep(ip).EQ.1)isthep(ip)=2
1314 c and daughter pointers of mother entry must be updated
1315 IF(jdahep(1,ip).EQ.0)
THEN
1319 jdahep(2,ip)=max(ihep,jdahep(2,ip))
1328 FUNCTION ihepdim(DUM)
1329 c this is the hepevt class in old style. no d_h_ class pre-name
1330 c this is the hepevt class in old style. no d_h_ class pre-name
1332 parameter(nmxhep=10000)
1334 INTEGER nevhep,nhep,isthep,idhep,jmohep,
1345 * ----------------------------------------------------------------------
1349 * ----------------------------------------------------------------------
1354 IMPLICIT REAL*8(a-h,o-z)
1355 COMPLEX*16 cprz0,cprz0m
1358 cprz0=dcmplx((s-amz**2),s/amz*gammz)
1360 zprop2=(abs(cprz0m))**2
1363 SUBROUTINE taupi0(MODE,JAK,ION)
1364 c no initialization required. must be called once after every:
1365 c 1) CALL dekay(1+10,...)
1366 c 2) CALL dekay(2+10,...)
1367 c 3) CALL dexay(1,...)
1368 c 4) CALL dexay(2,...)
1369 c
subroutine to decay originating from tauola
's taus:
1370 C 1) etas (with CALL TAUETA(JAK))
1371 C 2) later pi0's from taus.
1372 c 3) extensions to other applications possible.
1373 c this routine belongs to >tauola universal interface<, but uses
1374 c routines from >tauola< utilities as well. 25.08.2005
1375 c this is the hepevt class in old style. no d_h_ class pre-name
1377 parameter(nmxhep=10000)
1379 INTEGER nevhep,nhep,isthep,idhep,jmohep,
1390 * ----------------------------------------------------------------------
1394 * ----------------------------------------------------------------------
1397 c position of taus, must be defined by host program:
1398 COMMON /taupos/ np1,np2
1400 REAL phot1(4),phot2(4)
1401 REAL*8 r,x(4),y(4),pi0(4)
1402 INTEGER jezeli(3),ion(3)
1405 IF (mode.EQ.-1)
THEN
1411 IF (jezeli(1).EQ.0)
RETURN
1412 IF (jezeli(2).EQ.1) CALL taueta(jak)
1413 IF (jezeli(3).EQ.1) CALL tauk0s(jak)
1414 c position of decaying particle:
1415 IF((kto.EQ. 1).OR.(kto.EQ.11))
THEN
1421 DO k=jdahep(1,nps),nhepm
1422 IF (idhep(k).EQ.111.AND.jdahep(1,k).LE.k)
THEN
1427 r=sqrt(pi0(4)**2-pi0(3)**2-pi0(2)**2-pi0(1)**2)/2d0
1436 CALL bostdq(-1,pi0,x,x)
1437 CALL bostdq(-1,pi0,y,y)
1443 CALL filhep(0,1,22,k,k,0,0,phot1,0.0,.true.)
1444 CALL filhep(0,1,22,k,k,0,0,phot2,0.0,.true.)
1449 SUBROUTINE taueta(JAK)
1450 c
subroutine to decay etas
's from taus.
1451 C this routine belongs to tauola universal interface, but uses
1452 C routines from tauola utilities. Just flat phase space, but 4 channels.
1453 C it is called at the beginning of SUBR. TAUPI0(JAK)
1454 C and as far as hepevt search it is basically the same as TAUPI0. 25.08.2005
1455 C this is the hepevt class in old style. No d_h_ class pre-name
1457 PARAMETER (NMXHEP=10000)
1458 REAL*8 phep, vhep ! to be real*4/ *8 depending on host
1459 INTEGER nevhep,nhep,isthep,idhep,jmohep,
1462 $ nevhep, ! serial number
1463 $ nhep, ! number of particles
1464 $ isthep(nmxhep), ! status code
1465 $ idhep(nmxhep), ! particle ident KF
1466 $ jmohep(2,nmxhep), ! parent particles
1467 $ jdahep(2,nmxhep), ! childreen particles
1468 $ phep(5,nmxhep), ! four-momentum, mass [GeV]
1469 $ vhep(4,nmxhep) ! vertex [mm]
1470 * ----------------------------------------------------------------------
1473 $ qedrad(nmxhep) ! Photos flag
1474 * ----------------------------------------------------------------------
1476 COMMON / PARMAS / AMTAU,AMNUTA,AMEL,AMNUE,AMMU,AMNUMU
1477 * ,AMPIZ,AMPI,AMRO,GAMRO,AMA1,GAMA1
1478 * ,AMK,AMKZ,AMKST,GAMKST
1480 REAL*4 AMTAU,AMNUTA,AMEL,AMNUE,AMMU,AMNUMU
1481 * ,AMPIZ,AMPI,AMRO,GAMRO,AMA1,GAMA1
1482 * ,AMK,AMKZ,AMKST,GAMKST
1484 C position of taus, must be defined by host program:
1485 COMMON /TAUPOS/ NP1,NP2
1487 REAL RRR(1),BRSUM(3), RR(2)
1488 REAL PHOT1(4),PHOT2(4),PHOT3(4)
1489 REAL*8 X(4), Y(4), Z(4)
1491 REAL*8 R,RU,PETA(4),XM1,XM2,XM3,XM,XLAM
1493 XLAM(a,b,c)=SQRT(ABS((a-b-c)**2-4.0*b*c))
1494 C position of decaying particle:
1495 .EQ..OR..EQ.
IF((KTO 1)(KTO11)) THEN
1500 nhepM=nhep ! to avoid infinite loop
1501 DO K=JDAHEP(1,NPS),nhepM ! we search for etas's from tau till eor.
1502 IF (idhep(k).EQ.221.AND.jdahep(1,k).LE.k)
THEN
1506 c eta cumulated branching ratios:
1508 brsum(2)=brsum(1)+0.319
1509 brsum(3)=brsum(2)+0.237
1512 IF (rrr(1).LT.brsum(1))
THEN
1514 r=sqrt(peta(4)**2-peta(3)**2-peta(2)**2-peta(1)**2)/2d0
1523 CALL bostdq(-1,peta,x,x)
1524 CALL bostdq(-1,peta,y,y)
1530 CALL filhep(0,1,22,k,k,0,0,phot1,0.0,.true.)
1531 CALL filhep(0,1,22,k,k,0,0,phot2,0.0,.true.)
1533 IF(rrr(1).LT.brsum(2))
THEN
1540 ELSEIF(rrr(1).LT.brsum(3))
THEN
1557 r=sqrt(peta(4)**2-peta(3)**2-peta(2)**2-peta(1)**2)
1560 am2=sqrt(amin**2+rr(1)*(amax**2-amin**2))
1561 c weight for flat phase space
1562 wt=xlam(1d0*r**2,1d0*am2**2,1d0*xm3**2)
1563 & *xlam(1d0*am2**2,1d0*xm1**2,1d0*xm2**2)
1565 IF (rr(2).GT.wt) goto 7
1567 ru=xlam(1d0*am2**2,1d0*xm1**2,1d0*xm2**2)/am2/2
1571 x(4)=sqrt(ru**2+xm1**2)
1572 y(4)=sqrt(ru**2+xm2**2)
1577 c generate momentum of that pair in rest frame of eta:
1578 ru=xlam(1d0*r**2,1d0*am2**2,1d0*xm3**2)/r/2
1580 z(4)=sqrt(ru**2+am2**2)
1581 c and boost first two decay products to rest frame of eta.
1582 CALL bostdq(-1,z,x,x)
1583 CALL bostdq(-1,z,y,y)
1584 c redefine z(4) to 4-momentum of the last decay product:
1588 z(4)=sqrt(ru**2+xm3**2)
1589 c boost all to lab and move to real*4; also masses
1590 CALL bostdq(-1,peta,x,x)
1591 CALL bostdq(-1,peta,y,y)
1592 CALL bostdq(-1,peta,z,z)
1602 CALL filhep(0,1,id1,k,k,0,0,phot1,ym1,.true.)
1603 CALL filhep(0,1,id2,k,k,0,0,phot2,ym2,.true.)
1604 CALL filhep(0,1,id3,k,k,0,0,phot3,ym3,.true.)
1611 SUBROUTINE tauk0s(JAK)
1612 c
subroutine to decay k0s
's from taus.
1613 C this routine belongs to tauola universal interface, but uses
1614 C routines from tauola utilities. Just flat phase space, but 4 channels.
1615 C it is called at the beginning of SUBR. TAUPI0(JAK)
1616 C and as far as hepevt search it is basically the same as TAUPI0. 25.08.2005
1617 C this is the hepevt class in old style. No d_h_ class pre-name
1619 PARAMETER (NMXHEP=10000)
1620 REAL*8 phep, vhep ! to be real*4/ *8 depending on host
1621 INTEGER nevhep,nhep,isthep,idhep,jmohep,
1624 $ nevhep, ! serial number
1625 $ nhep, ! number of particles
1626 $ isthep(nmxhep), ! status code
1627 $ idhep(nmxhep), ! particle ident KF
1628 $ jmohep(2,nmxhep), ! parent particles
1629 $ jdahep(2,nmxhep), ! childreen particles
1630 $ phep(5,nmxhep), ! four-momentum, mass [GeV]
1631 $ vhep(4,nmxhep) ! vertex [mm]
1632 * ----------------------------------------------------------------------
1635 $ qedrad(nmxhep) ! Photos flag
1636 * ----------------------------------------------------------------------
1639 COMMON / PARMAS / AMTAU,AMNUTA,AMEL,AMNUE,AMMU,AMNUMU
1640 * ,AMPIZ,AMPI,AMRO,GAMRO,AMA1,GAMA1
1641 * ,AMK,AMKZ,AMKST,GAMKST
1643 REAL*4 AMTAU,AMNUTA,AMEL,AMNUE,AMMU,AMNUMU
1644 * ,AMPIZ,AMPI,AMRO,GAMRO,AMA1,GAMA1
1645 * ,AMK,AMKZ,AMKST,GAMKST
1647 C position of taus, must be defined by host program:
1648 COMMON /TAUPOS/ NP1,NP2
1650 REAL RRR(1),BRSUM(3), RR(2)
1651 REAL PHOT1(4),PHOT2(4),PHOT3(4)
1652 REAL*8 X(4), Y(4), Z(4)
1654 REAL*8 R,RU,PETA(4),XM1,XM2,XM3,XM,XLAM
1656 XLAM(a,b,c)=SQRT(ABS((a-b-c)**2-4.0*b*c))
1657 C position of decaying particle:
1658 .EQ..OR..EQ.
IF((KTO 1)(KTO11)) THEN
1663 nhepM=nhep ! to avoid infinite loop
1664 DO K=JDAHEP(1,NPS),nhepM ! we search for K0S's from tau till eor.
1665 IF (idhep(k).EQ.310.AND.jdahep(1,k).LE.k)
THEN
1671 c k0s cumulated branching ratios:
1674 brsum(3)=brsum(2)+0.237
1677 IF(rrr(1).LT.brsum(1))
THEN
1682 ELSEIF(rrr(1).LT.brsum(2))
THEN
1695 r=sqrt(peta(4)**2-peta(3)**2-peta(2)**2-peta(1)**2)/2d0
1697 r=sqrt(abs(r**2-xm1**2))
1706 CALL bostdq(-1,peta,x,x)
1707 CALL bostdq(-1,peta,y,y)
1716 CALL filhep(0,1,id1,k,k,0,0,phot1,ym1,.true.)
1717 CALL filhep(0,1,id2,k,k,0,0,phot2,ym2,.true.)