C++ Interface to Tauola
ffwid3pi.f
1 DOUBLE PRECISION FUNCTION ffwid3pi(QQ,S1,S3)
2 IMPLICIT NONE
3 DOUBLE PRECISION QQ,S1,S3
4C **************************************************************
5C Input: QQ S1 S3 ! mpi-pi-pi+**2 mpi-pi+**2 mpi-pi-**2
6C Calls: functions FORM1, FORM2, FORM4
7C Uses constants: tau mass, pi mass, normalization constant.
8C Load: Initialized tauola library,
9C Output: d\Gamma(tau --> 3pi nu)/(dQQ dS1 dS3)
10C Remark: If QQ S1 S3 are outside of the phase space
11C function FFWID3PI returns zero.
12C **************************************************************
13 COMPLEX F1,F2,F4, FORM1,FORM2, FORM4
14 DOUBLE PRECISION V11,V12,V22,GGF2,VUD2,ABS1,QQMIN,
15 & QQMAX,S3MAX,S3MIN,S1MIN,S1MAX
16 REAL XQQ,XS1,XS3, XS2,RQQ
17 COMMON / decpar / gfermi,gv,ga,ccabib,scabib,gamel
18 real*4 gfermi,gv,ga,ccabib,scabib,gamel
19 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
20 * ,ampiz,ampi,amro,gamro,ama1,gama1
21 * ,amk,amkz,amkst,gamkst
22 DOUBLE PRECISION XLAM,X,Y,Z
23 DOUBLE PRECISION XAMPI2
24 DOUBLE PRECISION GETFPIRPT
25 real*4 amtau,amnuta,amel,amnue,ammu,amnumu
26 * ,ampiz,ampi,amro,gamro,ama1,gama1
27 * ,amk,amkz,amkst,gamkst
28 DOUBLE PRECISION PI
29 DATA pi /3.141592653589793238462643d0/
30 INTEGER IMODE,IDUM,IFRCHL
31 REAL RRQ,RCHLWIDA1PI
32
33
34 xlam(x,y,z)= sqrt(abs((x-y-z)**2 - 4.*y*z))
35
36 abs1 = 1.d-5
37
38
39 ggf2 = gfermi**2
40 vud2 = ccabib**2
41
42
43C TO CHANGE THE VARIABLES TO SINGLE PRECISION
44C INPUT FOR FORM1,FORM2,FORM4 IS SINGLE PRECISION
45
46 xqq = qq
47 xs1 = s1
48 xs2 = qq -s1-s3 + 3.*ampi**2
49 xs3 = s3
50 xampi2 = ampi**2
51C Limits for PHASE SPACE
52C Limits for QQ
53 qqmin = 9.d0*ampi**2
54 qqmax = (amtau-amnuta)**2
55
56C Limits for S1
57 s1max=(dsqrt(qq) - ampi)**2 -abs1
58 s1min=4.d0*ampi**2 +abs1
59
60C LIMIT FOR XS3
61 s3max = (qq - ampi**2)**2 -
62 & ( xlam(qq,s1,xampi2)
63 & - xlam(s1,xampi2,xampi2) )**2
64 s3min = (qq - ampi**2)**2 -
65 & (xlam(qq,s1,xampi2)
66 & + xlam(s1,xampi2,xampi2) )**2
67
68 s3max = s3max/4./s1
69 s3min = s3min/4./s1
70
71C Check on PHASE SPACE
72C
73 IF((xs2.LE.0.) .OR.(s3max.LE.s3min)
74 & .OR.(xs1.LE.s1min).OR.(xs1.GE.s1max)
75 & .OR.(xs3.LE.s3min).OR.(xs3.GE.s3max)
76 & .OR.(qq.LE.qqmin).OR.(qq.GE.qqmax)
77 & ) THEN
78 ffwid3pi = 0.d0
79 RETURN
80 ENDIF
81
82C
83 v11 = -xs1+4.d0*ampi**2 -(xs2-xs3)**2/(4.d0*xqq)
84 v22 = -xs2+4.d0*ampi**2 - (xs3-xs1)**2/(4.d0*xqq)
85 v12 = 0.5d0*(xs3-xs1-xs2+4.d0*ampi**2)-0.25d0*(xs3-xs2)*(xs3-xs1)/xqq
86
87
88 f1 = form1(0,xqq,xs1,xs2)
89 f2 = form2(0,xqq,xs2,xs1)
90 f4 = form4(0,xqq,xs2,xs1,xs3)
91
92C formula 3.46 of [3]
93 ffwid3pi = abs(f1*conjg(f1))*v11+abs(f2*conjg(f2))*v22+
94 $ 2.d0*real(f1*conjg(f2))*v12
95
96 CALL ifgfact(2,imode,idum)
97 CALL inirchlget(ifrchl)
98 IF (imode.EQ.0) THEN
99C VERSION A: The 3 pion contribution to the a1 width
100C factor of a1 phase space and zeroing a1 propagator etc.is
101C done in RCHLWIDA1PI
102 rqq=qq
103 IF (ifrchl.EQ.1) THEN
104C CASE OF RCHL
105 ffwid3pi = rchlwida1pi(rqq,ffwid3pi)
106 ELSE
107C NOT READY YET
108 WRITE(*,*) 'FFWID3PI is not ready for non rchl currents'
109 stop
110 ENDIF
111C to get the total a1 width the contribution from (KKPI)- and K-K0pi0
112C channels have to be added
113 ELSE
114C VERSION B: calculation of 3 pion spectra in tau to 3pi nu channel.
115
116C factor for phase space of tau to XQQ nu decay and contribution from F4
117C (formula 3.21 of [3])
118
119 ffwid3pi = (- ffwid3pi/3.d0*(1.d0+2.d0*xqq/(amtau**2))+
120 $ xqq*abs(f4*conjg(f4)))*(amtau**2/xqq-1.d0)**2
121
122C Flux factor and normalization const.
123 ffwid3pi =
124 & ggf2*vud2/(128.d0*(2.d0*pi)**5*amtau)/2.d0
125 & *ffwid3pi
126 IF (ifrchl.EQ.1) THEN
127C CASE OF RCHL
128C RChL normalization constant
129 ffwid3pi =ffwid3pi/getfpirpt(1)**2
130 ELSE
131C NOT READY YET
132 WRITE(*,*) 'FFWID3PI is not ready for non rchl currents'
133 stop
134 ENDIF
135
136 ENDIF
137
138 RETURN
139 END
140
141 REAL FUNCTION RCHLWIDA1PI(RQQ,FFWID3PI)
142C The 3 pion contribution to the a1 width
143C (in [3] simple pretabulation is used through formula 3.48)
144C for calculation of g(QQ) of 3.45 3.46 of [3] in RChL style
145C a1 propagator has to be taken with the zero width.
146
147 IMPLICIT NONE
148 COMPLEX FA1RCHL
149 REAL RQQ
150 DOUBLE PRECISION FFWID3PI
151 COMMON / decpar / gfermi,gv,ga,ccabib,scabib,gamel
152 real*4 gfermi,gv,ga,ccabib,scabib,gamel
153 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
154 * ,ampiz,ampi,amro,gamro,ama1,gama1
155 * ,amk,amkz,amkst,gamkst
156 DOUBLE PRECISION XLAM,X,Y,Z
157 DOUBLE PRECISION XAMPI2
158 DOUBLE PRECISION GETFPIRPT
159 real*4 amtau,amnuta,amel,amnue,ammu,amnumu
160 * ,ampiz,ampi,amro,gamro,ama1,gama1
161 * ,amk,amkz,amkst,gamkst
162
163 common/rcht_3pi/ fpi_rpt,fv_rpt,gv_rpt,fa_rpt,beta_rho,fk_rpt
164 & ,fv1_rpt,gv1_rpt
165 DOUBLE PRECISION FPI_RPT,FV_RPT,GV_RPT,FA_RPT,BETA_RHO,FK_RPT
166 & ,FV1_RPT,GV1_RPT
167 DOUBLE PRECISION PI
168 DATA pi /3.141592653589793238462643d0/
169
170C AMA1 should be replaced by variable from the rchl namespace.
171 rchlwida1pi=- 1.0/real(fa1rchl(rqq)*conjg(fa1rchl(rqq)))/rqq**2
172 $ /(96.d0*8.d0*pi**3*ama1)/(fa_rpt**2*fpi_rpt**2)
173 $ *ffwid3pi/2.d0
174 END
175
176 DOUBLE PRECISION FUNCTION getfpirpt(I)
177 IMPLICIT NONE
178 common/rcht_3pi/ fpi_rpt,fv_rpt,gv_rpt,fa_rpt,beta_rho,fk_rpt
179 & ,fv1_rpt,gv1_rpt
180 DOUBLE PRECISION FPI_RPT,FV_RPT,GV_RPT,FA_RPT,BETA_RHO,FK_RPT
181 & ,FV1_RPT,GV1_RPT
182 INTEGER I
183 getfpirpt=fpi_rpt
184 END