C++InterfacetoTauola
F/prod/curr_cleo.f
1 /* copyright(c) 1991-2012 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  <http://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 /* we do support the iec 559 math functionality, real and complex. */
28 
29 /* wchar_t uses iso/iec 10646 (2nd ed., published 2011-03-15) /
30  unicode 6.0. */
31 
32 /* we do not support c11 <threads.h>. */
33 
34 *ajw 1 version of curr from koralb.
35  SUBROUTINE curr_cleo(MNUM,PIM1,PIM2,PIM3,PIM4,HADCUR)
36 c ==================================================================
37 c ajw, 11/97 - based on original curr from tauola:
38 c hadronic current for 4 pi final state
39 c r. fisher, j. wess and f. wagner z. phys c3(1980) 313
40 c r. decker z. phys c36(1987) 487.
41 c m. gell-mann, d. sharp, w. wagner phys. rev. lett 8 (1962) 261.
42 c but, rewritten to be more general and less "theoretical",
43 c using parameters tuned by vasia and dsc.
44 c ==================================================================
45 
46  COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
47  * ,ampiz,ampi,amro,gamro,ama1,gama1
48  * ,amk,amkz,amkst,gamkst
49 c
50  REAL*4 amtau,amnuta,amel,amnue,ammu,amnumu
51  * ,ampiz,ampi,amro,gamro,ama1,gama1
52  * ,amk,amkz,amkst,gamkst
53 c
54  REAL pim1(4),pim2(4),pim3(4),pim4(4)
55  COMPLEX hadcur(4)
56 
57  INTEGER k,l,mnum,k1,k2,iro,i,j,kk
58  REAL pa(4),pb(4),paa(4)
59  REAL aa(4,4),pp(4,4)
60  REAL a,xm,xg,g1,g2,g,amro2,gamro2,amro3,gamro3,amom,gamom
61  REAL fro,coef1,fpi,coef2,qq,sk,denom,sig,qqa,ss23,ss24,ss34,qp1p2
62  REAL qp1p3,qp1p4,p1p2,p1p3,p1p4,sign
63  REAL pkorb,ampa
64  COMPLEX alf0,alf1,alf2,alf3
65  COMPLEX lam0,lam1,lam2,lam3
66  COMPLEX bet1,bet2,bet3
67  COMPLEX form1,form2,form3,form4,form2pi
68  COMPLEX bwigm,wigfor,fpikm,fpikmd
69  COMPLEX ampl(7),ampr
70  COMPLEX bwign
71 c
72  bwign(a,xm,xg)=1.0/cmplx(a-xm**2,xm*xg)
73 c*******************************************************************************
74 c
75 c --- masses and constants
76  IF (g1.NE.12.924) THEN
77  g1=12.924
78  g2=1475.98
79  fpi=93.3e-3
80  g =g1*g2
81  fro=0.266*amro**2
82  coef1=2.0*sqrt(3.0)/fpi**2
83  coef2=fro*g ! overall constant for the omega current
84  coef2= coef2*0.56 ! factor 0.56 reduces contribution of omega from 68% to 40 %
85 
86 c masses and widths for for rho-prim and rho-bis:
87  amro2 = 1.465
88  gamro2= 0.310
89  amro3=1.700
90  gamro3=0.235
91 c
92  amom = pkorb(1,14)
93  gamom = pkorb(2,14)
94  amro2 = pkorb(1,21)
95  gamro2= pkorb(2,21)
96  amro3 = pkorb(1,22)
97  gamro3= pkorb(2,22)
98 c
99 c amplitudes for(pi-pi-pi0pi+) -> ps, rho0, rho-, rho+, omega.
100  ampl(1) = cmplx(pkorb(3,31)*coef1,0.)
101  ampl(2) = cmplx(pkorb(3,32)*coef1,0.)*cexp(cmplx(0.,pkorb(3,42)))
102  ampl(3) = cmplx(pkorb(3,33)*coef1,0.)*cexp(cmplx(0.,pkorb(3,43)))
103  ampl(4) = cmplx(pkorb(3,34)*coef1,0.)*cexp(cmplx(0.,pkorb(3,44)))
104  ampl(5) = cmplx(pkorb(3,35)*coef2,0.)*cexp(cmplx(0.,pkorb(3,45)))
105 c amplitudes for(pi0pi0pi0pi-) -> ps, rho-.
106  ampl(6) = cmplx(pkorb(3,36)*coef1)
107  ampl(7) = cmplx(pkorb(3,37)*coef1)
108 c
109 c rho' contributions to rho' -> pi-omega:
110  alf0 = cmplx(pkorb(3,51),0.0)
111  alf1 = cmplx(pkorb(3,52)*amro**2,0.0)
112  alf2 = cmplx(pkorb(3,53)*amro2**2,0.0)
113  alf3 = cmplx(pkorb(3,54)*amro3**2,0.0)
114 c rho' contribtions to rho' -> rhopipi:
115  lam0 = cmplx(pkorb(3,55),0.0)
116  lam1 = cmplx(pkorb(3,56)*amro**2,0.0)
117  lam2 = cmplx(pkorb(3,57)*amro2**2,0.0)
118  lam3 = cmplx(pkorb(3,58)*amro3**2,0.0)
119 c rho contributions to rhopipi, rho -> 2pi:
120  bet1 = cmplx(pkorb(3,59)*amro**2,0.0)
121  bet2 = cmplx(pkorb(3,60)*amro2**2,0.0)
122  bet3 = cmplx(pkorb(3,61)*amro3**2,0.0)
123 c
124  END IF
125 c**************************************************
126 c
127 c --- initialization of four vectors
128  DO 7 k=1,4
129  DO 8 l=1,4
130  8 aa(k,l)=0.0
131  hadcur(k)=cmplx(0.0)
132  paa(k)=pim1(k)+pim2(k)+pim3(k)+pim4(k)
133  pp(1,k)=pim1(k)
134  pp(2,k)=pim2(k)
135  pp(3,k)=pim3(k)
136  7 pp(4,k)=pim4(k)
137 c
138  IF (mnum.EQ.1) THEN
139 c ===================================================================
140 c pi- pi- p0 pi+ case ====
141 c ===================================================================
142  qq=paa(4)**2-paa(3)**2-paa(2)**2-paa(1)**2
143 
144 c add m(4pi)-dependence to rhopipi channels:
145  form4= lam0+lam1*bwign(qq,amro,gamro)
146  * +lam2*bwign(qq,amro2,gamro2)
147  * +lam3*bwign(qq,amro3,gamro3)
148 
149 c --- loop over five contributions of the rho-pi-pi
150  DO 201 k1=1,3
151  DO 201 k2=3,4
152 c
153  IF (k2.EQ.k1) THEN
154  goto 201
155  ELSEIF (k2.EQ.3) THEN
156 c rho-
157  ampr = ampl(3)
158  ampa = ampiz
159  ELSEIF (k1.EQ.3) THEN
160 c rho+
161  ampr = ampl(4)
162  ampa = ampiz
163  ELSE
164 c rho0
165  ampr = ampl(2)
166  ampa = ampi
167  END IF
168 c
169  sk=(pp(k1,4)+pp(k2,4))**2-(pp(k1,3)+pp(k2,3))**2
170  $ -(pp(k1,2)+pp(k2,2))**2-(pp(k1,1)+pp(k2,1))**2
171 
172 c -- definition of aa matrix
173 c -- cronecker delta
174  DO 202 i=1,4
175  DO 203 j=1,4
176  203 aa(i,j)=0.0
177  202 aa(i,i)=1.0
178 c ... and the rest ...
179  DO 204 l=1,4
180  IF (l.NE.k1.AND.l.NE.k2) THEN
181  denom=(paa(4)-pp(l,4))**2-(paa(3)-pp(l,3))**2
182  $ -(paa(2)-pp(l,2))**2-(paa(1)-pp(l,1))**2
183  DO 205 i=1,4
184  DO 205 j=1,4
185  sig= 1.0
186  IF(j.NE.4) sig=-sig
187  aa(i,j)=aa(i,j)
188  $ -sig*(paa(i)-2.0*pp(l,i))*(paa(j)-pp(l,j))/denom
189  205 CONTINUE
190  ENDIF
191  204 CONTINUE
192 c
193 c --- lets add something to hadcurr
194 c form1= fpikm(sqrt(sk),ampi,ampi) *fpikm(sqrt(qq),ampi,ampi)
195 c form1= ampl(1)+ampr*fpikm(sqrt(sk),ampi,ampi)
196 
197  form2pi= bet1*bwigm(sk,amro,gamro,ampa,ampi)
198  1 +bet2*bwigm(sk,amro2,gamro2,ampa,ampi)
199  2 +bet3*bwigm(sk,amro3,gamro3,ampa,ampi)
200  form1= ampl(1)+ampr*form2pi
201 c
202  DO 206 i=1,4
203  DO 206 j=1,4
204  hadcur(i)=hadcur(i)+form1*form4*aa(i,j)*(pp(k1,j)-pp(k2,j))
205  206 CONTINUE
206 c --- end of the rho-pi-pi current(5 possibilities)
207  201 CONTINUE
208 c
209 c ===================================================================
210 c now modify the coefficient for the omega-pi current: =
211 c ===================================================================
212  IF (ampl(5).EQ.cmplx(0.,0.)) goto 311
213 
214 c overall rho+rhoprime for the 4pi system:
215 c form2=ampl(5)*(bwign(qq,amro,gamro)+elpha*bwign(qq,amrop,gamrop))
216 c modified m(4pi)-dependence:
217  form2=ampl(5)*(alf0+alf1*bwign(qq,amro,gamro)
218  * +alf2*bwign(qq,amro2,gamro2)
219  * +alf3*bwign(qq,amro3,gamro3))
220 c
221 c --- there are two possibilities for omega current
222 c --- pa pb are corresponding first and second pi-s
223  DO 301 kk=1,2
224  DO 302 i=1,4
225  pa(i)=pp(kk,i)
226  pb(i)=pp(3-kk,i)
227  302 CONTINUE
228 c --- lorentz invariants
229  qqa=0.0
230  ss23=0.0
231  ss24=0.0
232  ss34=0.0
233  qp1p2=0.0
234  qp1p3=0.0
235  qp1p4=0.0
236  p1p2 =0.0
237  p1p3 =0.0
238  p1p4 =0.0
239  DO 303 k=1,4
240  sign=-1.0
241  IF (k.EQ.4) sign= 1.0
242  qqa=qqa+sign*(paa(k)-pa(k))**2
243  ss23=ss23+sign*(pb(k) +pim3(k))**2
244  ss24=ss24+sign*(pb(k) +pim4(k))**2
245  ss34=ss34+sign*(pim3(k)+pim4(k))**2
246  qp1p2=qp1p2+sign*(paa(k)-pa(k))*pb(k)
247  qp1p3=qp1p3+sign*(paa(k)-pa(k))*pim3(k)
248  qp1p4=qp1p4+sign*(paa(k)-pa(k))*pim4(k)
249  p1p2=p1p2+sign*pa(k)*pb(k)
250  p1p3=p1p3+sign*pa(k)*pim3(k)
251  p1p4=p1p4+sign*pa(k)*pim4(k)
252  303 CONTINUE
253 c
254 c omega -> rho pi for the 3pi system:
255 c form3=bwign(qqa,amom,gamom)*(bwign(ss23,amro,gamro)+
256 c $ bwign(ss24,amro,gamro)+bwign(ss34,amro,gamro))
257 c no omega -> rho pi; just straight omega:
258  form3=bwign(qqa,amom,gamom)
259 c
260  DO 304 k=1,4
261  hadcur(k)=hadcur(k)+form2*form3*(
262  $ pb(k)*(qp1p3*p1p4-qp1p4*p1p3)
263  $ +pim3(k)*(qp1p4*p1p2-qp1p2*p1p4)
264  $ +pim4(k)*(qp1p2*p1p3-qp1p3*p1p2) )
265  304 CONTINUE
266  301 CONTINUE
267  311 CONTINUE
268 c
269  ELSE
270 c ===================================================================
271 c pi0 pi0 p0 pi- case ====
272 c ===================================================================
273  qq=paa(4)**2-paa(3)**2-paa(2)**2-paa(1)**2
274 
275 c --- loop over three contribution of the non-omega current
276  DO 101 k=1,3
277  sk=(pp(k,4)+pim4(4))**2-(pp(k,3)+pim4(3))**2
278  $ -(pp(k,2)+pim4(2))**2-(pp(k,1)+pim4(1))**2
279 
280 c -- definition of aa matrix
281 c -- cronecker delta
282  DO 102 i=1,4
283  DO 103 j=1,4
284  103 aa(i,j)=0.0
285  102 aa(i,i)=1.0
286 c
287 c ... and the rest ...
288  DO 104 l=1,3
289  IF (l.NE.k) THEN
290  denom=(paa(4)-pp(l,4))**2-(paa(3)-pp(l,3))**2
291  $ -(paa(2)-pp(l,2))**2-(paa(1)-pp(l,1))**2
292  DO 105 i=1,4
293  DO 105 j=1,4
294  sig=1.0
295  IF(j.NE.4) sig=-sig
296  aa(i,j)=aa(i,j)
297  $ -sig*(paa(i)-2.0*pp(l,i))*(paa(j)-pp(l,j))/denom
298  105 CONTINUE
299  ENDIF
300  104 CONTINUE
301 
302 c --- lets add something to hadcurr
303 c form1= fpikm(sqrt(sk),ampi,ampi) *fpikmd(sqrt(qq),ampi,ampi)
304 ccccccccccccc form1=wigfor(sk,amro,gamro) (tests)
305 c form1= fpikm(sqrt(sk),ampi,ampi) *fpikm(sqrt(qq),ampi,ampi)
306  form1 = ampl(6)+ampl(7)*fpikm(sqrt(sk),ampi,ampi)
307 
308  DO 106 i=1,4
309  DO 106 j=1,4
310  hadcur(i)=hadcur(i)+form1*aa(i,j)*(pp(k,j)-pp(4,j))
311  106 CONTINUE
312 c --- end of the non omega current(3 possibilities)
313  101 CONTINUE
314 
315  ENDIF
316  END
317 
318 
319