C++InterfacetoTauola
modified/curr_cleo.f
1 
2 
3 *AJW 1 version of CURR from KORALB.
4  SUBROUTINE curr_cleo(MNUM,PIM1,PIM2,PIM3,PIM4,HADCUR)
5 C ==================================================================
6 C AJW, 11/97 - based on original CURR from TAUOLA:
7 C hadronic current for 4 pi final state
8 C R. Fisher, J. Wess and F. Wagner Z. Phys C3 (1980) 313
9 C R. Decker Z. Phys C36 (1987) 487.
10 C M. Gell-Mann, D. Sharp, W. Wagner Phys. Rev. Lett 8 (1962) 261.
11 C BUT, rewritten to be more general and less "theoretical",
12 C using parameters tuned by Vasia and DSC.
13 C ==================================================================
14 
15  COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
16  * ,ampiz,ampi,amro,gamro,ama1,gama1
17  * ,amk,amkz,amkst,gamkst
18 C
19  REAL*4 amtau,amnuta,amel,amnue,ammu,amnumu
20  * ,ampiz,ampi,amro,gamro,ama1,gama1
21  * ,amk,amkz,amkst,gamkst
22 C
23  REAL pim1(4),pim2(4),pim3(4),pim4(4)
24  COMPLEX hadcur(4)
25 
26  INTEGER k,l,mnum,k1,k2,iro,i,j,kk
27  REAL pa(4),pb(4),paa(4)
28  REAL aa(4,4),pp(4,4)
29  REAL a,xm,xg,g1,g2,g,amro2,gamro2,amro3,gamro3,amom,gamom
30  REAL fro,coef1,fpi,coef2,qq,sk,denom,sig,qqa,ss23,ss24,ss34,qp1p2
31  REAL qp1p3,qp1p4,p1p2,p1p3,p1p4,sign
32  REAL pkorb,ampa
33  COMPLEX alf0,alf1,alf2,alf3
34  COMPLEX lam0,lam1,lam2,lam3
35  COMPLEX bet1,bet2,bet3
36  COMPLEX form1,form2,form3,form4,form2pi
37  COMPLEX bwigm,wigfor,fpikm,fpikmd
38  COMPLEX ampl(7),ampr
39  COMPLEX bwign
40 C
41  bwign(a,xm,xg)=1.0/cmplx(a-xm**2,xm*xg)
42 C*******************************************************************************
43 C
44 C --- masses and constants
45  IF (g1.NE.12.924) THEN
46  g1=12.924
47  g2=1475.98
48  fpi=93.3e-3
49  g =g1*g2
50  fro=0.266*amro**2
51  coef1=2.0*sqrt(3.0)/fpi**2
52  coef2=fro*g ! overall constant for the omega current
53  coef2= coef2*0.56 ! factor 0.56 reduces contribution of omega from 68% to 40 %
54 
55 C masses and widths for for rho-prim and rho-bis:
56  amro2 = 1.465
57  gamro2= 0.310
58  amro3=1.700
59  gamro3=0.235
60 C
61  amom = pkorb(1,14)
62  gamom = pkorb(2,14)
63  amro2 = pkorb(1,21)
64  gamro2= pkorb(2,21)
65  amro3 = pkorb(1,22)
66  gamro3= pkorb(2,22)
67 C
68 C Amplitudes for (pi-pi-pi0pi+) -> PS, rho0, rho-, rho+, omega.
69  ampl(1) = cmplx(pkorb(3,31)*coef1,0.)
70  ampl(2) = cmplx(pkorb(3,32)*coef1,0.)*cexp(cmplx(0.,pkorb(3,42)))
71  ampl(3) = cmplx(pkorb(3,33)*coef1,0.)*cexp(cmplx(0.,pkorb(3,43)))
72  ampl(4) = cmplx(pkorb(3,34)*coef1,0.)*cexp(cmplx(0.,pkorb(3,44)))
73  ampl(5) = cmplx(pkorb(3,35)*coef2,0.)*cexp(cmplx(0.,pkorb(3,45)))
74 C Amplitudes for (pi0pi0pi0pi-) -> PS, rho-.
75  ampl(6) = cmplx(pkorb(3,36)*coef1)
76  ampl(7) = cmplx(pkorb(3,37)*coef1)
77 C
78 C rho' contributions to rho' -> pi-omega:
79  alf0 = cmplx(pkorb(3,51),0.0)
80  alf1 = cmplx(pkorb(3,52)*amro**2,0.0)
81  alf2 = cmplx(pkorb(3,53)*amro2**2,0.0)
82  alf3 = cmplx(pkorb(3,54)*amro3**2,0.0)
83 C rho' contribtions to rho' -> rhopipi:
84  lam0 = cmplx(pkorb(3,55),0.0)
85  lam1 = cmplx(pkorb(3,56)*amro**2,0.0)
86  lam2 = cmplx(pkorb(3,57)*amro2**2,0.0)
87  lam3 = cmplx(pkorb(3,58)*amro3**2,0.0)
88 C rho contributions to rhopipi, rho -> 2pi:
89  bet1 = cmplx(pkorb(3,59)*amro**2,0.0)
90  bet2 = cmplx(pkorb(3,60)*amro2**2,0.0)
91  bet3 = cmplx(pkorb(3,61)*amro3**2,0.0)
92 C
93  END IF
94 C**************************************************
95 C
96 C --- initialization of four vectors
97  DO 7 k=1,4
98  DO 8 l=1,4
99  8 aa(k,l)=0.0
100  hadcur(k)=cmplx(0.0)
101  paa(k)=pim1(k)+pim2(k)+pim3(k)+pim4(k)
102  pp(1,k)=pim1(k)
103  pp(2,k)=pim2(k)
104  pp(3,k)=pim3(k)
105  7 pp(4,k)=pim4(k)
106 C
107  IF (mnum.EQ.1) THEN
108 C ===================================================================
109 C pi- pi- p0 pi+ case ====
110 C ===================================================================
111  qq=paa(4)**2-paa(3)**2-paa(2)**2-paa(1)**2
112 
113 C Add M(4pi)-dependence to rhopipi channels:
114  form4= lam0+lam1*bwign(qq,amro,gamro)
115  * +lam2*bwign(qq,amro2,gamro2)
116  * +lam3*bwign(qq,amro3,gamro3)
117 
118 C --- loop over five contributions of the rho-pi-pi
119  DO 201 k1=1,3
120  DO 201 k2=3,4
121 C
122  IF (k2.EQ.k1) THEN
123  goto 201
124  ELSEIF (k2.EQ.3) THEN
125 C rho-
126  ampr = ampl(3)
127  ampa = ampiz
128  ELSEIF (k1.EQ.3) THEN
129 C rho+
130  ampr = ampl(4)
131  ampa = ampiz
132  ELSE
133 C rho0
134  ampr = ampl(2)
135  ampa = ampi
136  END IF
137 C
138  sk=(pp(k1,4)+pp(k2,4))**2-(pp(k1,3)+pp(k2,3))**2
139  $ -(pp(k1,2)+pp(k2,2))**2-(pp(k1,1)+pp(k2,1))**2
140 
141 C -- definition of AA matrix
142 C -- cronecker delta
143  DO 202 i=1,4
144  DO 203 j=1,4
145  203 aa(i,j)=0.0
146  202 aa(i,i)=1.0
147 C ... and the rest ...
148  DO 204 l=1,4
149  IF (l.NE.k1.AND.l.NE.k2) THEN
150  denom=(paa(4)-pp(l,4))**2-(paa(3)-pp(l,3))**2
151  $ -(paa(2)-pp(l,2))**2-(paa(1)-pp(l,1))**2
152  DO 205 i=1,4
153  DO 205 j=1,4
154  sig= 1.0
155  IF(j.NE.4) sig=-sig
156  aa(i,j)=aa(i,j)
157  $ -sig*(paa(i)-2.0*pp(l,i))*(paa(j)-pp(l,j))/denom
158  205 CONTINUE
159  ENDIF
160  204 CONTINUE
161 C
162 C --- lets add something to HADCURR
163 C FORM1= FPIKM(SQRT(SK),AMPI,AMPI) *FPIKM(SQRT(QQ),AMPI,AMPI)
164 C FORM1= AMPL(1)+AMPR*FPIKM(SQRT(SK),AMPI,AMPI)
165 
166  form2pi= bet1*bwigm(sk,amro,gamro,ampa,ampi)
167  1 +bet2*bwigm(sk,amro2,gamro2,ampa,ampi)
168  2 +bet3*bwigm(sk,amro3,gamro3,ampa,ampi)
169  form1= ampl(1)+ampr*form2pi
170 C
171  DO 206 i=1,4
172  DO 206 j=1,4
173  hadcur(i)=hadcur(i)+form1*form4*aa(i,j)*(pp(k1,j)-pp(k2,j))
174  206 CONTINUE
175 C --- end of the rho-pi-pi current (5 possibilities)
176  201 CONTINUE
177 C
178 C ===================================================================
179 C Now modify the coefficient for the omega-pi current: =
180 C ===================================================================
181  IF (ampl(5).EQ.cmplx(0.,0.)) goto 311
182 
183 C Overall rho+rhoprime for the 4pi system:
184 C FORM2=AMPL(5)*(BWIGN(QQ,AMRO,GAMRO)+ELPHA*BWIGN(QQ,AMROP,GAMROP))
185 C Modified M(4pi)-dependence:
186  form2=ampl(5)*(alf0+alf1*bwign(qq,amro,gamro)
187  * +alf2*bwign(qq,amro2,gamro2)
188  * +alf3*bwign(qq,amro3,gamro3))
189 C
190 C --- there are two possibilities for omega current
191 C --- PA PB are corresponding first and second pi-s
192  DO 301 kk=1,2
193  DO 302 i=1,4
194  pa(i)=pp(kk,i)
195  pb(i)=pp(3-kk,i)
196  302 CONTINUE
197 C --- lorentz invariants
198  qqa=0.0
199  ss23=0.0
200  ss24=0.0
201  ss34=0.0
202  qp1p2=0.0
203  qp1p3=0.0
204  qp1p4=0.0
205  p1p2 =0.0
206  p1p3 =0.0
207  p1p4 =0.0
208  DO 303 k=1,4
209  sign=-1.0
210  IF (k.EQ.4) sign= 1.0
211  qqa=qqa+sign*(paa(k)-pa(k))**2
212  ss23=ss23+sign*(pb(k) +pim3(k))**2
213  ss24=ss24+sign*(pb(k) +pim4(k))**2
214  ss34=ss34+sign*(pim3(k)+pim4(k))**2
215  qp1p2=qp1p2+sign*(paa(k)-pa(k))*pb(k)
216  qp1p3=qp1p3+sign*(paa(k)-pa(k))*pim3(k)
217  qp1p4=qp1p4+sign*(paa(k)-pa(k))*pim4(k)
218  p1p2=p1p2+sign*pa(k)*pb(k)
219  p1p3=p1p3+sign*pa(k)*pim3(k)
220  p1p4=p1p4+sign*pa(k)*pim4(k)
221  303 CONTINUE
222 C
223 C omega -> rho pi for the 3pi system:
224 C FORM3=BWIGN(QQA,AMOM,GAMOM)*(BWIGN(SS23,AMRO,GAMRO)+
225 C $ BWIGN(SS24,AMRO,GAMRO)+BWIGN(SS34,AMRO,GAMRO))
226 C No omega -> rho pi; just straight omega:
227  form3=bwign(qqa,amom,gamom)
228 C
229  DO 304 k=1,4
230  hadcur(k)=hadcur(k)+form2*form3*(
231  $ pb(k)*(qp1p3*p1p4-qp1p4*p1p3)
232  $ +pim3(k)*(qp1p4*p1p2-qp1p2*p1p4)
233  $ +pim4(k)*(qp1p2*p1p3-qp1p3*p1p2) )
234  304 CONTINUE
235  301 CONTINUE
236  311 CONTINUE
237 C
238  ELSE
239 C ===================================================================
240 C pi0 pi0 p0 pi- case ====
241 C ===================================================================
242  qq=paa(4)**2-paa(3)**2-paa(2)**2-paa(1)**2
243 
244 C --- loop over three contribution of the non-omega current
245  DO 101 k=1,3
246  sk=(pp(k,4)+pim4(4))**2-(pp(k,3)+pim4(3))**2
247  $ -(pp(k,2)+pim4(2))**2-(pp(k,1)+pim4(1))**2
248 
249 C -- definition of AA matrix
250 C -- cronecker delta
251  DO 102 i=1,4
252  DO 103 j=1,4
253  103 aa(i,j)=0.0
254  102 aa(i,i)=1.0
255 C
256 C ... and the rest ...
257  DO 104 l=1,3
258  IF (l.NE.k) THEN
259  denom=(paa(4)-pp(l,4))**2-(paa(3)-pp(l,3))**2
260  $ -(paa(2)-pp(l,2))**2-(paa(1)-pp(l,1))**2
261  DO 105 i=1,4
262  DO 105 j=1,4
263  sig=1.0
264  IF(j.NE.4) sig=-sig
265  aa(i,j)=aa(i,j)
266  $ -sig*(paa(i)-2.0*pp(l,i))*(paa(j)-pp(l,j))/denom
267  105 CONTINUE
268  ENDIF
269  104 CONTINUE
270 
271 C --- lets add something to HADCURR
272 C FORM1= FPIKM(SQRT(SK),AMPI,AMPI) *FPIKMD(SQRT(QQ),AMPI,AMPI)
273 CCCCCCCCCCCCC FORM1=WIGFOR(SK,AMRO,GAMRO) (tests)
274 C FORM1= FPIKM(SQRT(SK),AMPI,AMPI) *FPIKM(SQRT(QQ),AMPI,AMPI)
275  form1 = ampl(6)+ampl(7)*fpikm(sqrt(sk),ampi,ampi)
276 
277  DO 106 i=1,4
278  DO 106 j=1,4
279  hadcur(i)=hadcur(i)+form1*aa(i,j)*(pp(k,j)-pp(4,j))
280  106 CONTINUE
281 C --- end of the non omega current (3 possibilities)
282  101 CONTINUE
283 
284  ENDIF
285  END
286 
287 
288