C++ Interface to Tauola
frho_pi_belle.f
1 Complex*16 FUNCTION fpibel(W,flpar)
2c****************************************************************************
3c pion form factor from belle
4c
5c rho(770)+rho(1450)+rho(1700) model is used as in:
6c m.fujikawa et al., phys.rev.d 78 (2008) 072006
7c -------------------
8c * flpar: = 0 - (all free) fit result for par(1...11)
9c = 1 - (par(1)=f_pi(0)=1-fixed) fit result for par(2...11)
10c
11c * parameter : par(1)=overall norm. factor(=f_pi(0))
12c : par(2)=rho(770) mass
13c : par(3)=rho(770) width
14c : par(4)=rho(1450) mass
15c : par(5)=rho(1450) width
16c : par(6)=rho(1450) admixture abs. value(|beta|)
17c : par(7)=rho(1450) admixture phase(phi=arg(beta))
18c : par(8)=rho(1700) mass
19c : par(9)=rho(1700) width
20c : par(10)=rho(1700) admixture abs. value(|gamma|)
21c : par(11)=rho(1700) admixture phase(phi3=arg(gamma))
22c
23c * x:=s=(m_pipi0)^2=w^2
24c
25c * notice: the following code was extracted(and checked) directly
26c from the belle fit function with minor cosmetic changes
27c(this is internal comment of belle)
28c
29c * notice: the following code was extracted(and checked) directly
30c from tauola_with_belle_fpi.f sent by d. epifanov and
31c it is an identical copy of
32c Complex*16 FUNCTION fpibel(W,flpar)
33c used in tauola_with_belle_fpi.f
34c
35c * called: curr_pipi0 if (ff2pirho = 2 or ff2pirho =3)
36c in ../rchl-currents/value_parameter.f
37c ff2pirho = 2 is for flpar = 0
38c ff2pirho = 3 is for flpar = 1
39c
40c****************************************************************************
41 implicit none
42 integer flpar
43 real W
44 real*8 x,par(11),pi,pimas,pi0mas,taumas
45 real*8 beta,berho1,berho2,berho3
46 real*8 ps,prho1,prho2,prho3
47 real*8 gamma1,gamma2,gamma3
48 real*8 hs,h1,h2,h3,dhds1,dhds2,dhds3
49 real*8 d1,d2,d3,fs1,fs2,fs3
50 real*8 a1,b1,c1,bw_re1,bw_im1
51 real*8 a2,b2,c2,bw_re2,bw_im2
52 real*8 a3,b3,c3,bw_re3,bw_im3
53 real*8 sinphi,cosphi,sinphi3,cosphi3
54 complex*16 mbw1,mbw2,mbw3,mff,mephi,mephi3
55 parameter(pi=3.141592653589)
56
57 x=dble(w**2)
58c------------- parameters --------------------
59c particle mass(unit: gev)
60
61 pimas = 0.1395702
62 pi0mas = 0.13498
63 taumas = 1.77699
64
65c---------------------------------------------
66c opt. par. from table vii of prd78(2008) 072006
67
68 if(flpar.eq.0) then
69 par(1) = 1.02
70 par(2) = 0.7749
71 par(3) = 0.1486
72 par(4) = 1.428
73 par(5) = 0.413
74 par(6) = 0.13
75 par(7) = 197
76 par(8) = 1.694
77 par(9) = 0.135
78 par(10)= 0.028
79 par(11)= -3
80 else
81 par(1) = 1.0
82 par(2) = 0.7746
83 par(3) = 0.1481
84 par(4) = 1.446
85 par(5) = 0.434
86 par(6) = 0.15
87 par(7) = 202
88 par(8) = 1.728
89 par(9) = 0.164
90 par(10)= 0.037
91 par(11)= 24
92 endif
93
94c========= beta(s, rho(770), rho(1450), rho(1700)) ==================
95
96 beta =sqrt( (1.0 - (pimas - pi0mas)**2/x)
97 + *(1.0 - (pimas + pi0mas)**2/x) )
98
99c ----- rho(770) -----
100 berho1=sqrt( ( 1.0 - (pimas - pi0mas)**2/(par(2)*par(2)) )
101 + *( 1.0 - (pimas + pi0mas)**2/(par(2)*par(2)) ) )
102
103c ----- rho(1450) -----
104 berho2=sqrt( ( 1.0 - (pimas - pi0mas)**2/(par(4)*par(4)) )
105 + *( 1.0 - (pimas + pi0mas)**2/(par(4)*par(4)) ) )
106
107c ----- rho(1700) -----
108 berho3=sqrt( ( 1.0 - (pimas - pi0mas)**2/(par(8)*par(8)) )
109 + *( 1.0 - (pimas + pi0mas)**2/(par(8)*par(8)) ) )
110
111c========= momentum(s, rho(770), rho(1450), rho(1700)) ==============
112c
113c ps : momentum of pi at s
114c prho1 : momentum of pi at s=rho(770)**2
115c prho2 : momentum of pi at s=rho(1450)**2
116c prho3 : momentum of pi at s=rho(1700)**2
117
118 ps = 0.5 * sqrt(x) * beta
119c ----- rho(770) -----
120 prho1 = 0.5* par(2) * berho1
121c ----- rho(1450) -----
122 prho2 = 0.5* par(4) * berho2
123c ----- rho(1700) -----
124 prho3 = 0.5* par(8) * berho3
125
126c========= width(rho(770), rho(1450), rho(1700)) ====================
127
128c ----- rho(770) -----
129 gamma1 = par(3) * ( par(2)*par(2)/x ) * (ps/prho1)**3
130c ----- rho(1450) -----
131 gamma2 = par(5) * ( par(4)*par(4)/x ) * (ps/prho2)**3
132c ----- rho(1700) -----
133 gamma3 = par(9) * ( par(8)*par(8)/x ) * (ps/prho3)**3
134
135c============== h(rho(770), rho(1450), rho(1700)) ====================
136c
137c hs : h(s)...s(s=x)
138c h1 : h(s) for s=rho(770)**2
139c h2 : h(s) for s=rho(1450)**2
140c h3 : h(s) for s=rho(1700)**2
141
142 hs = (2.0/pi)*(ps/sqrt(x))*
143 + log( (sqrt(x) + 2.0*ps)/(2.0*pimas) )
144c ----- rho(770) -----
145 h1 = (2.0/pi)*(prho1/par(2))*
146 + log( (par(2) + 2.0*prho1)/(2.0*pimas) )
147c ----- rho(1450) -----
148 h2 = (2.0/pi)*(prho2/par(4))*
149 + log( (par(4) + 2.0*prho2)/(2.0*pimas) )
150c ----- rho(1700) -----
151 h3 = (2.0/pi)*(prho3/par(8))*
152 + log( (par(8) + 2.0*prho3)/(2.0*pimas) )
153
154c============ dhds(rho(770), rho(1450), rho(1700)) ===================
155c
156c dhds = dh/ds = h(m_rho)[1/(8*p(m_rho)^2) - 1/(2m_rho^2)
157c + 1/(2pi*m_rho^2)
158c dhds1 : dh/ds for m_rho=rho(770)
159c dhds2 : dh/ds for m_rho=rho(1450)
160c dhds3 : dh/ds for m_rho=rho(1700)
161
162c ----- rho(770) -----
163 dhds1 = h1*( 1.0/(8.0*prho1**2) - 1.0/(2.0*par(2)**2) )
164 + + 1.0/(2.0*pi*par(2)**2)
165c ----- rho(1450) -----
166 dhds2 = h2*( 1.0/(8.0*prho2**2) - 1.0/(2.0*par(4)**2) )
167 + + 1.0/(2.0*pi*par(4)**2)
168c ----- rho(1700) -----
169 dhds3 = h3*( 1.0/(8.0*prho3**2) - 1.0/(2.0*par(8)**2) )
170 + + 1.0/(2.0*pi*par(8)**2)
171
172c============ d(rho(770), rho(1450), rho(1700)) ======================
173c
174c d = 3/pi * (m_pi^2/p(m_rho)^2) *
175c ln {(m_rho+2*p(m_rho))/2m_pi} + m_rho/(2pi*p(m_rho))
176c - (m_pi^2 * m_rho)/(pi*p(m_rho)^3)
177c d1 : d for m_rho=rho(770)
178c d2 : d for m_rho=rho(1450)
179c d2 : d for m_rho=rho(1700)
180
181c ----- rho(770) -----
182 d1 = (3.0/pi)*(pimas**2/prho1**2)
183 + * log((par(2) + 2.0*prho1)/(2.0*pimas))
184 + + par(2)/(2.0*pi*prho1)
185 + - ((pimas**2)*par(2))/(pi*prho1**3)
186c ----- rho(1450) -----
187 d2 = (3.0/pi)*(pimas**2/prho2**2)
188 + * log((par(4) + 2.0*prho2)/(2.0*pimas))
189 + + par(4)/(2.0*pi*prho2)
190 + - ((pimas**2)*par(4))/(pi*prho2**3)
191c ----- rho(1700) -----
192 d3 = (3.0/pi)*(pimas**2/prho3**2)
193 + * log((par(8) + 2.0*prho3)/(2.0*pimas))
194 + + par(8)/(2.0*pi*prho3)
195 + - ((pimas**2)*par(8))/(pi*prho3**3)
196
197c================ f(s) (rho(770), rho(1450), rho(1700) ) ==========
198c
199c f(s) = gamma-rho * m_rhp**2/p(m_rho)[ p(s)^2(h(s)-h(m_rho)
200c + (m_rho^2 -s)p(m_rho)^2 * dh/ds]
201c fs1 : f(s) for m_rho=rho(770)
202c fs2 : f(s) for m_rho=rho(1450)
203c fs3 : f(s) for m_rho=rho(1700)
204c
205c ----- rho(770) -----
206 fs1 = par(3) * (par(2)**2/prho1**3) *
207 + ( ps**2 *(hs - h1) + (par(2)**2 -x)*prho1**2 * dhds1 )
208c ----- rho(1450) -----
209 fs2 = par(5) * (par(4)**2/prho2**3) *
210 + ( ps**2 *(hs - h2) + (par(4)**2 -x)*prho2**2 * dhds2 )
211c ----- rho(1700) -----
212 fs3 = par(9) * (par(8)**2/prho3**3) *
213 + ( ps**2 *(hs - h3) + (par(8)**2 -x)*prho3**2 * dhds3 )
214
215c====== bw form g&s model(rho(770) + rho(1450) + rho(1700)) ======
216c bw_re - real part ; bw_im - imaginary part
217
218c ----- rho(770) -----
219 a1=par(2)*par(2) - x + fs1
220 b1=sqrt(x)*gamma1
221 c1=par(2)*par(2) + d1*par(2)*par(3)
222
223 bw_re1= (a1*c1)/(a1**2 + b1**2)
224 bw_im1= (b1*c1)/(a1**2 + b1**2)
225 mbw1=dcmplx(bw_re1,bw_im1)
226
227c ----- rho(1450) -----
228 a2=par(4)*par(4) - x + fs2
229 b2=sqrt(x)*gamma2
230 c2=par(4)*par(4)+ d2*par(4)*par(5)
231
232 bw_re2= (a2*c2)/(a2**2 + b2**2)
233 bw_im2= (b2*c2)/(a2**2 + b2**2)
234 mbw2=dcmplx(bw_re2,bw_im2)
235
236c ----- rho(1700) -----
237 a3=par(8)*par(8) - x + fs3
238 b3=sqrt(x)*gamma3
239 c3=par(8)*par(8)+ d3*par(8)*par(9)
240
241 bw_re3= (a3*c3)/(a3**2 + b3**2)
242 bw_im3= (b3*c3)/(a3**2 + b3**2)
243 mbw3=dcmplx(bw_re3,bw_im3)
244
245c====== admixtures of rho(1450) and rho(1700) ======
246
247c --------- exp(i*arg(beta)) ----------
248 sinphi=sin(par(7)*(pi/180.0)) ! use degree
249 cosphi=cos(par(7)*(pi/180.0)) ! use degree
250 mephi=dcmplx(cosphi,sinphi)
251
252c --------- exp(i*arg(gamma)) ---------
253 sinphi3=sin(par(11)*(pi/180.0)) ! use degree
254 cosphi3=cos(par(11)*(pi/180.0)) ! use degree
255 mephi3=dcmplx(cosphi3,sinphi3)
256
257c================= form factor =====================
258
259 mff = 1.0/(1.0+par(6)*mephi+par(10)*mephi3)
260 fpibel=par(1)*mff*(mbw1+par(6)*mephi*mbw2+par(10)*mephi3*mbw3)
261 RETURN
262 END