C++ Interface to Tauola
tauola/formf.f
1/* copyright(c) 1991-2024 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 <https://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/* glibc's intent is to support the IEC 559 math functionality, real
28 and complex. If the GCC (4.9 and later) predefined macros
29 specifying compiler intent are available, use them to determine
30 whether the overall intent is to support these features; otherwise,
31 presume an older compiler has intent to support these features and
32 define these macros by default. */
33
34
35
36/* wchar_t uses Unicode 10.0.0. Version 10.0 of the Unicode Standard is
37 synchronized with ISO/IEC 10646:2017, fifth edition, plus
38 the following additions from Amendment 1 to the fifth edition:
39 - 56 emoji characters
40 - 285 hentaigana
41 - 3 additional Zanabazar Square characters */
42
43 FUNCTION FORMOM(XMAA,XMOM)
44C ==================================================================
45C formfactorfor pi-pi0 gamma final state
46C R. Decker, Z. Phys C36 (1987) 487.
47C ==================================================================
48 COMMON / PARMAS / AMTAU,AMNUTA,AMEL,AMNUE,AMMU,AMNUMU
49 * ,AMPIZ,AMPI,AMRO,GAMRO,AMA1,GAMA1
50 * ,AMK,AMKZ,AMKST,GAMKST
51C
52 REAL*4 AMTAU,AMNUTA,AMEL,AMNUE,AMMU,AMNUMU
53 * ,AMPIZ,AMPI,AMRO,GAMRO,AMA1,GAMA1
54 * ,AMK,AMKZ,AMKST,GAMKST
55 COMMON / DECPAR / GFERMI,GV,GA,CCABIB,SCABIB,GAMEL
56 REAL*4 GFERMI,GV,GA,CCABIB,SCABIB,GAMEL
57 COMMON /TESTA1/ KEYA1
58 COMPLEX BWIGN,FORMOM
59 DATA ICONT /1/
60* THIS INLINE FUNCT. CALCULATES THE SCALAR PART OF THE PROPAGATOR
61 BWIGN(XM,AM,GAMMA)=1./CMPLX(XM**2-AM**2,GAMMA*AM)
62* HADRON CURRENT
63 FRO =0.266*AMRO**2
64 ELPHA=- 0.1
65 AMROP = 1.7
66 GAMROP= 0.26
67 AMOM =0.782
68 GAMOM =0.0085
69 AROMEG= 1.0
70 GCOUP=12.924
71 GCOUP=GCOUP*AROMEG
72 FQED =SQRT(4.0*3.1415926535/137.03604)
73 FORMOM=FQED*FRO**2/SQRT(2.0)*GCOUP**2*BWIGN(XMOM,AMOM,GAMOM)
74 $ *(BWIGN(XMAA,AMRO,GAMRO)+ELPHA*BWIGN(XMAA,AMROP,GAMROP))
75 $ *(BWIGN( 0.0,AMRO,GAMRO)+ELPHA*BWIGN( 0.0,AMROP,GAMROP))
76 END
77C=======================================================================
78 COMPLEX FUNCTION FK1AB(XMSQ,INDX)
79C ==================================================================
80C complex form-factor for a1+a1prime. AJW 1/98
81C ==================================================================
82
83 COMPLEX F1,F2,AMPA,AMPB
84 INTEGER IFIRST,INDX
85 DATA IFIRST/0/
86
87.EQ. IF (IFIRST0) THEN
88 IFIRST = 1
89 XM1 = PKORB(1,19)
90 XG1 = PKORB(2,19)
91 XM2 = PKORB(1,20)
92 XG2 = PKORB(2,20)
93
94 XM1SQ = XM1*XM1
95 GF1 = GFUN(XM1SQ)
96 GG1 = XM1*XG1/GF1
97 XM2SQ = XM2*XM2
98 GF2 = GFUN(XM2SQ)
99 GG2 = XM2*XG2/GF2
100 END IF
101
102.EQ. IF (INDX1) THEN
103 AMPA = CMPLX(PKORB(3,81),0.)
104 AMPB = CMPLX(PKORB(3,82),0.)
105.EQ. ELSE IF (INDX2) THEN
106 AMPA = CMPLX(PKORB(3,83),0.)
107 AMPB = CMPLX(PKORB(3,84),0.)
108.EQ. ELSEIF (INDX3) THEN
109 AMPA = CMPLX(PKORB(3,85),0.)
110 AMPB = CMPLX(PKORB(3,86),0.)
111.EQ. ELSEIF (INDX4) THEN
112 AMPA = CMPLX(PKORB(3,87),0.)
113 AMPB = CMPLX(PKORB(3,88),0.)
114 END IF
115
116 GF = GFUN(XMSQ)
117 FG1 = GG1*GF
118 FG2 = GG2*GF
119 F1 = CMPLX(-XM1SQ,0.0)/CMPLX(XMSQ-XM1SQ,FG1)
120 F2 = CMPLX(-XM2SQ,0.0)/CMPLX(XMSQ-XM2SQ,FG2)
121 FK1AB = AMPA*F1+AMPB*F2
122
123 RETURN
124 END
125 FUNCTION FORM1(MNUM,QQ,S1,SDWA)
126C ==================================================================
127C formfactorfor F1 for 3 scalar final state
128C R. Fisher, J. Wess and F. Wagner Z. Phys C3 (1980) 313
129C H. Georgi, Weak interactions and modern particle theory,
130C The Benjamin/Cummings Pub. Co., Inc. 1984.
131C R. Decker, E. Mirkes, R. Sauer, Z. Was Karlsruhe preprint TTP92-25
132C and erratum !!!!!!
133C ==================================================================
134C
135 COMPLEX FORM1,WIGNER,WIGFOR,FPIKM,BWIGM
136 COMMON / PARMAS / AMTAU,AMNUTA,AMEL,AMNUE,AMMU,AMNUMU
137 * ,AMPIZ,AMPI,AMRO,GAMRO,AMA1,GAMA1
138 * ,AMK,AMKZ,AMKST,GAMKST
139C
140 REAL*4 AMTAU,AMNUTA,AMEL,AMNUE,AMMU,AMNUMU
141 * ,AMPIZ,AMPI,AMRO,GAMRO,AMA1,GAMA1
142 * ,AMK,AMKZ,AMKST,GAMKST
143 COMPLEX FORMA1,FORMK1,FORMRO,FORMKS
144 COMPLEX FA1A1P,FK1AB,F3PI
145C
146.EQ. IF (MNUM0) THEN
147C ------------ 3 pi hadronic state (a1)
148C FORMRO = FPIKM(SQRT(S1),AMPI,AMPI)
149C FORMRO = F3PI(1,QQ,S1,SDWA)
150C FORMA1 = FA1A1P(QQ)
151C FORM1 = FORMA1*FORMRO
152 FORM1 = F3PI(1,QQ,S1,SDWA)
153
154.EQ. ELSEIF (MNUM1) THEN
155C ------------ K- pi- K+ (K*0 K-)
156 FORMKS = BWIGM(S1,AMKST,GAMKST,AMPI,AMK)
157 FORMA1 = FA1A1P(QQ)
158 FORM1 = FORMA1*FORMKS
159
160.EQ. ELSEIF (MNUM2) THEN
161C ------------ K0 pi- K0B (K*- K0)
162 FORMKS = BWIGM(S1,AMKST,GAMKST,AMPI,AMK)
163 FORMA1 = FA1A1P(QQ)
164 FORM1 = FORMA1*FORMKS
165
166.EQ. ELSEIF (MNUM3) THEN
167C ------------ K- pi0 K0 (K*0 K-)
168 FORMKS = BWIGM(S1,AMKST,GAMKST,AMPI,AMK)
169 FORMA1 = FA1A1P(QQ)
170 FORM1 = FORMA1*FORMKS
171
172.EQ. ELSEIF (MNUM4) THEN
173C ------------ pi0 pi0 K- (K*-pi0)
174 FORMKS = BWIGM(S1,AMKST,GAMKST,AMPI,AMK)
175 FORMK1 = FK1AB(QQ,3)
176 FORM1 = FORMK1*FORMKS
177
178.EQ. ELSEIF (MNUM5) THEN
179C ------------ K- pi- pi+ (rho0 K-)
180 FORMK1 = FK1AB(QQ,4)
181 FORMRO = FPIKM(SQRT(S1),AMPI,AMPI)
182 FORM1 = FORMK1*FORMRO
183
184.EQ. ELSEIF (MNUM6) THEN
185C ------------ pi- K0B pi0 (pi- K*0B)
186 FORMK1 = FK1AB(QQ,1)
187 FORMKS = BWIGM(S1,AMKST,GAMKST,AMK,AMPI)
188 FORM1 = FORMK1*FORMKS
189
190.EQ. ELSEIF (MNUM7) THEN
191C -------------- eta pi- pi0 final state
192 FORM1=0.0
193 ENDIF
194 END
195 FUNCTION FORM2(MNUM,QQ,S1,SDWA)
196C ==================================================================
197C formfactorfor F2 for 3 scalar final state
198C R. Fisher, J. Wess and F. Wagner Z. Phys C3 (1980) 313
199C H. Georgi, Weak interactions and modern particle theory,
200C The Benjamin/Cummings Pub. Co., Inc. 1984.
201C R. Decker, E. Mirkes, R. Sauer, Z. Was Karlsruhe preprint TTP92-25
202C and erratum !!!!!!
203C ==================================================================
204C
205 COMPLEX FORM2,WIGNER,WIGFOR,FPIKM,BWIGM
206 COMMON / PARMAS / AMTAU,AMNUTA,AMEL,AMNUE,AMMU,AMNUMU
207 * ,AMPIZ,AMPI,AMRO,GAMRO,AMA1,GAMA1
208 * ,AMK,AMKZ,AMKST,GAMKST
209C
210 REAL*4 AMTAU,AMNUTA,AMEL,AMNUE,AMMU,AMNUMU
211 * ,AMPIZ,AMPI,AMRO,GAMRO,AMA1,GAMA1
212 * ,AMK,AMKZ,AMKST,GAMKST
213 COMPLEX FORMA1,FORMK1,FORMRO,FORMKS
214 COMPLEX FA1A1P,FK1AB,F3PI
215
216.EQ. IF (MNUM0) THEN
217C ------------ 3 pi hadronic state (a1)
218C FORMRO = FPIKM(SQRT(S1),AMPI,AMPI)
219C FORMRO = F3PI(2,QQ,S1,SDWA)
220C FORMA1 = FA1A1P(QQ)
221C FORM2 = FORMA1*FORMRO
222 FORM2 = F3PI(2,QQ,S1,SDWA)
223
224.EQ. ELSEIF (MNUM1) THEN
225C ------------ K- pi- K+ (rho0 pi-)
226 FORMRO = FPIKM(SQRT(S1),AMK,AMK)
227 FORMA1 = FA1A1P(QQ)
228 FORM2 = FORMA1*FORMRO
229
230.EQ. ELSEIF (MNUM2) THEN
231C ------------ K0 pi- K0B (rho0 pi-)
232 FORMRO = FPIKM(SQRT(S1),AMK,AMK)
233 FORMA1 = FA1A1P(QQ)
234 FORM2 = FORMA1*FORMRO
235
236.EQ. ELSEIF (MNUM3) THEN
237C ------------ K- pi0 K0 (rho- pi0)
238 FORMRO = FPIKM(SQRT(S1),AMK,AMK)
239 FORMA1 = FA1A1P(QQ)
240 FORM2 = FORMA1*FORMRO
241
242.EQ. ELSEIF (MNUM4) THEN
243C ------------ pi0 pi0 K- (K*-pi0)
244 FORMKS = BWIGM(S1,AMKST,GAMKST,AMPI,AMK)
245 FORMK1 = FK1AB(QQ,3)
246 FORM2 = FORMK1*FORMKS
247
248.EQ. ELSEIF (MNUM5) THEN
249C ------------ K- pi- pi+ (K*0B pi-)
250 FORMKS = BWIGM(S1,AMKST,GAMKST,AMPI,AMK)
251 FORMK1 = FK1AB(QQ,1)
252 FORM2 = FORMK1*FORMKS
253C
254.EQ. ELSEIF (MNUM6) THEN
255C ------------ pi- K0B pi0 (rho- K0B)
256 FORMRO = FPIKM(SQRT(S1),AMPI,AMPI)
257 FORMK1 = FK1AB(QQ,2)
258 FORM2 = FORMK1*FORMRO
259C
260.EQ. ELSEIF (MNUM7) THEN
261C -------------- eta pi- pi0 final state
262 FORM2=0.0
263 ENDIF
264C
265 END
266 COMPLEX FUNCTION BWIGM(S,M,G,XM1,XM2)
267C **********************************************************
268C P-WAVE BREIT-WIGNER FOR RHO
269C **********************************************************
270 REAL S,M,G,XM1,XM2
271 REAL PI,QS,QM,W,GS
272 DATA INIT /0/
273C ------------ PARAMETERS --------------------
274.EQ. IF (INIT0) THEN
275 INIT=1
276 PI=3.141592654
277C ------- BREIT-WIGNER -----------------------
278 ENDIF
279.GT. IF (S(XM1+XM2)**2) THEN
280 QS=SQRT(ABS((S -(XM1+XM2)**2)*(S -(XM1-XM2)**2)))/SQRT(S)
281 QM=SQRT(ABS((M**2-(XM1+XM2)**2)*(M**2-(XM1-XM2)**2)))/M
282 W=SQRT(S)
283 GS=G*(M/W)**2*(QS/QM)**3
284 ELSE
285 GS=0.0
286 ENDIF
287 BWIGM=M**2/CMPLX(M**2-S,-SQRT(S)*GS)
288 RETURN
289 END
290 COMPLEX FUNCTION FPIKM(W,XM1,XM2)
291C **********************************************************
292C PION FORM FACTOR
293C **********************************************************
294 COMPLEX BWIGM
295 REAL ROM,ROG,ROM1,ROG1,BETA1,PI,PIM,S,W
296 EXTERNAL BWIG
297 DATA INIT /0/
298C
299C ------------ PARAMETERS --------------------
300.EQ. IF (INIT0 ) THEN
301 INIT=1
302 PI=3.141592654
303 PIM=.140
304 ROM=0.773
305 ROG=0.145
306 ROM1=1.370
307 ROG1=0.510
308 BETA1=-0.145
309 ENDIF
310C -----------------------------------------------
311 S=W**2
312 FPIKM=(BWIGM(S,ROM,ROG,XM1,XM2)+BETA1*BWIGM(S,ROM1,ROG1,XM1,XM2))
313 & /(1+BETA1)
314 RETURN
315 END
316 COMPLEX FUNCTION FPIKMD(W,XM1,XM2)
317C **********************************************************
318C PION FORM FACTOR
319C **********************************************************
320 COMPLEX BWIGM
321 REAL ROM,ROG,ROM1,ROG1,PI,PIM,S,W
322 EXTERNAL BWIG
323 DATA INIT /0/
324C
325C ------------ PARAMETERS --------------------
326.EQ. IF (INIT0 ) THEN
327 INIT=1
328 PI=3.141592654
329 PIM=.140
330 ROM=0.773
331 ROG=0.145
332 ROM1=1.500
333 ROG1=0.220
334 ROM2=1.750
335 ROG2=0.120
336 BETA=6.5
337 DELTA=-26.0
338 ENDIF
339C -----------------------------------------------
340 S=W**2
341 FPIKMD=(DELTA*BWIGM(S,ROM,ROG,XM1,XM2)
342 $ +BETA*BWIGM(S,ROM1,ROG1,XM1,XM2)
343 $ + BWIGM(S,ROM2,ROG2,XM1,XM2))
344 & /(1+BETA+DELTA)
345 RETURN
346 END
347
348 FUNCTION FORM3(MNUM,QQ,S1,SDWA)
349C ==================================================================
350C formfactorfor F3 for 3 scalar final state
351C R. Fisher, J. Wess and F. Wagner Z. Phys C3 (1980) 313
352C H. Georgi, Weak interactions and modern particle theory,
353C The Benjamin/Cummings Pub. Co., Inc. 1984.
354C R. Decker, E. Mirkes, R. Sauer, Z. Was Karlsruhe preprint TTP92-25
355C and erratum !!!!!!
356C ==================================================================
357C
358 COMMON / PARMAS / AMTAU,AMNUTA,AMEL,AMNUE,AMMU,AMNUMU
359 * ,AMPIZ,AMPI,AMRO,GAMRO,AMA1,GAMA1
360 * ,AMK,AMKZ,AMKST,GAMKST
361C
362 REAL*4 AMTAU,AMNUTA,AMEL,AMNUE,AMMU,AMNUMU
363 * ,AMPIZ,AMPI,AMRO,GAMRO,AMA1,GAMA1
364 * ,AMK,AMKZ,AMKST,GAMKST
365 COMPLEX FORM3,BWIGM
366 COMPLEX FORMA1,FORMK1,FORMRO,FORMKS
367 COMPLEX FA1A1P,FK1AB,F3PI
368C
369.EQ. IF (MNUM0) THEN
370C ------------ 3 pi hadronic state (a1)
371C FORMRO = FPIKM(SQRT(S1),AMPI,AMPI)
372C FORMRO = F3PI(3,QQ,S1,SDWA)
373C FORMA1 = FA1A1P(QQ)
374C FORM3 = FORMA1*FORMRO
375 FORM3 = F3PI(3,QQ,S1,SDWA)
376
377.EQ. ELSEIF (MNUM3) THEN
378C ------------ K- pi0 K0 (K*- K0)
379 FORMKS = BWIGM(S1,AMKST,GAMKST,AMPIZ,AMK)
380 FORMA1 = FA1A1P(QQ)
381 FORM3 = FORMA1*FORMKS
382
383.EQ. ELSEIF (MNUM6) THEN
384C ------------ pi- K0B pi0 (K*- pi0)
385 FORMKS = BWIGM(S1,AMKST,GAMKST,AMK,AMPI)
386 FORMK1 = FK1AB(QQ,3)
387 FORM3 = FORMK1*FORMKS
388
389 ELSE
390 FORM3=CMPLX(0.,0.)
391 ENDIF
392 END
393 FUNCTION FORM4(MNUM,QQ,S1,S2,S3)
394C ==================================================================
395C formfactorfor F4 for 3 scalar final state
396C R. Decker, in preparation
397C R. Decker, E. Mirkes, R. Sauer, Z. Was Karlsruhe preprint TTP92-25
398C and erratum !!!!!!
399C ==================================================================
400C
401 COMMON / PARMAS / AMTAU,AMNUTA,AMEL,AMNUE,AMMU,AMNUMU
402 * ,AMPIZ,AMPI,AMRO,GAMRO,AMA1,GAMA1
403 * ,AMK,AMKZ,AMKST,GAMKST
404C
405 REAL*4 AMTAU,AMNUTA,AMEL,AMNUE,AMMU,AMNUMU
406 * ,AMPIZ,AMPI,AMRO,GAMRO,AMA1,GAMA1
407 * ,AMK,AMKZ,AMKST,GAMKST
408 COMPLEX FORM4,WIGNER,FPIKM
409 REAL*4 M
410C ---- this formfactor is switched off .. .
411 FORM4=CMPLX(0.0,0.0)
412 END
413 FUNCTION FORM5(MNUM,QQ,S1,S2)
414C ==================================================================
415C formfactorfor F5 for 3 scalar final state
416C G. Kramer, W. Palmer, S. Pinsky, Phys. Rev. D30 (1984) 89.
417C G. Kramer, W. Palmer Z. Phys. C25 (1984) 195.
418C R. Decker, E. Mirkes, R. Sauer, Z. Was Karlsruhe preprint TTP92-25
419C and erratum !!!!!!
420C ==================================================================
421C
422 COMMON / PARMAS / AMTAU,AMNUTA,AMEL,AMNUE,AMMU,AMNUMU
423 * ,AMPIZ,AMPI,AMRO,GAMRO,AMA1,GAMA1
424 * ,AMK,AMKZ,AMKST,GAMKST
425C
426 REAL*4 AMTAU,AMNUTA,AMEL,AMNUE,AMMU,AMNUMU
427 * ,AMPIZ,AMPI,AMRO,GAMRO,AMA1,GAMA1
428 * ,AMK,AMKZ,AMKST,GAMKST
429 COMPLEX FORM5,WIGNER,FPIKM,FPIKMD,BWIGM
430.EQ. IF (MNUM0) THEN
431C ------------ 3 pi hadronic state (a1)
432 FORM5=0.0
433.EQ. ELSEIF (MNUM1) THEN
434C ------------ K- pi- K+
435 ELPHA=-0.2
436 FORM5=FPIKMD(SQRT(QQ),AMPI,AMPI)/(1+ELPHA)
437 $ *( FPIKM(SQRT(S2),AMPI,AMPI)
438 $ +ELPHA*BWIGM(S1,AMKST,GAMKST,AMPI,AMK))
439.EQ. ELSEIF (MNUM2) THEN
440C ------------ K0 pi- K0B
441 ELPHA=-0.2
442 FORM5=FPIKMD(SQRT(QQ),AMPI,AMPI)/(1+ELPHA)
443 $ *( FPIKM(SQRT(S2),AMPI,AMPI)
444 $ +ELPHA*BWIGM(S1,AMKST,GAMKST,AMPI,AMK))
445.EQ. ELSEIF (MNUM3) THEN
446C ------------ K- K0 pi0
447 FORM5=0.0
448.EQ. ELSEIF (MNUM4) THEN
449C ------------ pi0 pi0 K-
450 FORM5=0.0
451.EQ. ELSEIF (MNUM5) THEN
452C ------------ K- pi- pi+
453 ELPHA=-0.2
454 FORM5=BWIGM(QQ,AMKST,GAMKST,AMPI,AMK)/(1+ELPHA)
455 $ *( FPIKM(SQRT(S1),AMPI,AMPI)
456 $ +ELPHA*BWIGM(S2,AMKST,GAMKST,AMPI,AMK))
457.EQ. ELSEIF (MNUM6) THEN
458C ------------ pi- K0B pi0
459 ELPHA=-0.2
460 FORM5=BWIGM(QQ,AMKST,GAMKST,AMPI,AMKZ)/(1+ELPHA)
461 $ *( FPIKM(SQRT(S2),AMPI,AMPI)
462 $ +ELPHA*BWIGM(S1,AMKST,GAMKST,AMPI,AMK))
463.EQ. ELSEIF (MNUM7) THEN
464C -------------- eta pi- pi0 final state
465 FORM5=FPIKMD(SQRT(QQ),AMPI,AMPI)*FPIKM(SQRT(S1),AMPI,AMPI)
466 ENDIF
467C
468 END
469 SUBROUTINE CURRX(MNUM,PIM1,PIM2,PIM3,PIM4,HADCUR)
470C ==================================================================
471C hadronic current for 4 pi final state
472C R. Fisher, J. Wess and F. Wagner Z. Phys C3 (1980) 313
473C R. Decker Z. Phys C36 (1987) 487.
474C M. Gell-Mann, D. Sharp, W. Wagner Phys. Rev. Lett 8 (1962) 261.
475C ==================================================================
476
477 COMMON / PARMAS / AMTAU,AMNUTA,AMEL,AMNUE,AMMU,AMNUMU
478 * ,AMPIZ,AMPI,AMRO,GAMRO,AMA1,GAMA1
479 * ,AMK,AMKZ,AMKST,GAMKST
480C
481 REAL*4 AMTAU,AMNUTA,AMEL,AMNUE,AMMU,AMNUMU
482 * ,AMPIZ,AMPI,AMRO,GAMRO,AMA1,GAMA1
483 * ,AMK,AMKZ,AMKST,GAMKST
484 COMMON / DECPAR / GFERMI,GV,GA,CCABIB,SCABIB,GAMEL
485 REAL*4 GFERMI,GV,GA,CCABIB,SCABIB,GAMEL
486C ARBITRARY FIXING OF THE FOUR PI X-SECTION NORMALIZATION
487 COMMON /ARBIT/ ARFLAT,AROMEG
488 REAL PIM1(4),PIM2(4),PIM3(4),PIM4(4),PAA(4)
489 COMPLEX HADCUR(4),FORM1,FORM2,FORM3,FPIKM
490 COMPLEX BWIGN
491 REAL PA(4),PB(4)
492 REAL AA(4,4),PP(4,4)
493 DATA PI /3.141592653589793238462643/
494 DATA FPI /93.3E-3/
495 BWIGN(A,XM,XG)=1.0/CMPLX(A-XM**2,XM*XG)
496C
497C --- masses and constants
498 G1=12.924
499 G2=1475.98
500 G =G1*G2
501 ELPHA=-.1
502 AMROP=1.7
503 GAMROP=0.26
504 AMOM=.782
505 GAMOM=0.0085
506 ARFLAT=1.0
507 AROMEG=1.0
508C
509 FRO=0.266*AMRO**2
510 COEF1=2.0*SQRT(3.0)/FPI**2*ARFLAT
511 COEF2=FRO*G*AROMEG
512C --- initialization of four vectors
513 DO 7 K=1,4
514 DO 8 L=1,4
515 8 AA(K,L)=0.0
516 HADCUR(K)=CMPLX(0.0)
517 PAA(K)=PIM1(K)+PIM2(K)+PIM3(K)+PIM4(K)
518 PP(1,K)=PIM1(K)
519 PP(2,K)=PIM2(K)
520 PP(3,K)=PIM3(K)
521 7 PP(4,K)=PIM4(K)
522C
523.EQ. IF (MNUM1) THEN
524C ===================================================================
525C pi- pi- p0 pi+ case ====
526C ===================================================================
527 QQ=PAA(4)**2-PAA(3)**2-PAA(2)**2-PAA(1)**2
528C --- loop over thre contribution of the non-omega current
529 DO 201 K=1,3
530 SK=(PP(K,4)+PIM4(4))**2-(PP(K,3)+PIM4(3))**2
531 $ -(PP(K,2)+PIM4(2))**2-(PP(K,1)+PIM4(1))**2
532C -- definition of AA matrix
533C -- cronecker delta
534 DO 202 I=1,4
535 DO 203 J=1,4
536 203 AA(I,J)=0.0
537 202 AA(I,I)=1.0
538C ... and the rest ...
539 DO 204 L=1,3
540.NE. IF (LK) THEN
541 DENOM=(PAA(4)-PP(L,4))**2-(PAA(3)-PP(L,3))**2
542 $ -(PAA(2)-PP(L,2))**2-(PAA(1)-PP(L,1))**2
543 DO 205 I=1,4
544 DO 205 J=1,4
545 SIG= 1.0
546.NE. IF(J4) SIG=-SIG
547 AA(I,J)=AA(I,J)
548 $ -SIG*(PAA(I)-2.0*PP(L,I))*(PAA(J)-PP(L,J))/DENOM
549 205 CONTINUE
550 ENDIF
551 204 CONTINUE
552C --- lets add something to HADCURR
553 FORM1= FPIKM(SQRT(SK),AMPI,AMPI) *FPIKM(SQRT(QQ),AMPI,AMPI)
554C FORM1= FPIKM(SQRT(SK),AMPI,AMPI) *FPIKMD(SQRT(QQ),AMPI,AMPI)
555CCCCCCCCCCCCCCCCC FORM1=WIGFOR(SK,AMRO,GAMRO) (tests)
556C
557 FIX=1.0
558.EQ. IF (K3) FIX=-2.0
559 DO 206 I=1,4
560 DO 206 J=1,4
561 HADCUR(I)=
562 $ HADCUR(I)+CMPLX(FIX*COEF1)*FORM1*AA(I,J)*(PP(K,J)-PP(4,J))
563 206 CONTINUE
564C --- end of the non omega current (3 possibilities)
565 201 CONTINUE
566C
567C
568C --- there are two possibilities for omega current
569C --- PA PB are corresponding first and second pi-s
570 DO 301 KK=1,2
571 DO 302 I=1,4
572 PA(I)=PP(KK,I)
573 PB(I)=PP(3-KK,I)
574 302 CONTINUE
575C --- lorentz invariants
576 QQA=0.0
577 SS23=0.0
578 SS24=0.0
579 SS34=0.0
580 QP1P2=0.0
581 QP1P3=0.0
582 QP1P4=0.0
583 P1P2 =0.0
584 P1P3 =0.0
585 P1P4 =0.0
586 DO 303 K=1,4
587 SIGN=-1.0
588.EQ. IF (K4) SIGN= 1.0
589 QQA=QQA+SIGN*(PAA(K)-PA(K))**2
590 SS23=SS23+SIGN*(PB(K) +PIM3(K))**2
591 SS24=SS24+SIGN*(PB(K) +PIM4(K))**2
592 SS34=SS34+SIGN*(PIM3(K)+PIM4(K))**2
593 QP1P2=QP1P2+SIGN*(PAA(K)-PA(K))*PB(K)
594 QP1P3=QP1P3+SIGN*(PAA(K)-PA(K))*PIM3(K)
595 QP1P4=QP1P4+SIGN*(PAA(K)-PA(K))*PIM4(K)
596 P1P2=P1P2+SIGN*PA(K)*PB(K)
597 P1P3=P1P3+SIGN*PA(K)*PIM3(K)
598 P1P4=P1P4+SIGN*PA(K)*PIM4(K)
599 303 CONTINUE
600C
601 FORM2=COEF2*(BWIGN(QQ,AMRO,GAMRO)+ELPHA*BWIGN(QQ,AMROP,GAMROP))
602C FORM3=BWIGN(QQA,AMOM,GAMOM)*(BWIGN(SS23,AMRO,GAMRO)+
603C $ BWIGN(SS24,AMRO,GAMRO)+BWIGN(SS34,AMRO,GAMRO))
604 FORM3=BWIGN(QQA,AMOM,GAMOM)
605C
606 DO 304 K=1,4
607 HADCUR(K)=HADCUR(K)+FORM2*FORM3*(
608 $ PB (K)*(QP1P3*P1P4-QP1P4*P1P3)
609 $ +PIM3(K)*(QP1P4*P1P2-QP1P2*P1P4)
610 $ +PIM4(K)*(QP1P2*P1P3-QP1P3*P1P2) )
611 304 CONTINUE
612 301 CONTINUE
613C
614 ELSE
615C ===================================================================
616C pi0 pi0 p0 pi- case ====
617C ===================================================================
618 QQ=PAA(4)**2-PAA(3)**2-PAA(2)**2-PAA(1)**2
619 DO 101 K=1,3
620C --- loop over thre contribution of the non-omega current
621 SK=(PP(K,4)+PIM4(4))**2-(PP(K,3)+PIM4(3))**2
622 $ -(PP(K,2)+PIM4(2))**2-(PP(K,1)+PIM4(1))**2
623C -- definition of AA matrix
624C -- cronecker delta
625 DO 102 I=1,4
626 DO 103 J=1,4
627 103 AA(I,J)=0.0
628 102 AA(I,I)=1.0
629C
630C ... and the rest ...
631 DO 104 L=1,3
632.NE. IF (LK) THEN
633 DENOM=(PAA(4)-PP(L,4))**2-(PAA(3)-PP(L,3))**2
634 $ -(PAA(2)-PP(L,2))**2-(PAA(1)-PP(L,1))**2
635 DO 105 I=1,4
636 DO 105 J=1,4
637 SIG=1.0
638.NE. IF(J4) SIG=-SIG
639 AA(I,J)=AA(I,J)
640 $ -SIG*(PAA(I)-2.0*PP(L,I))*(PAA(J)-PP(L,J))/DENOM
641 105 CONTINUE
642 ENDIF
643 104 CONTINUE
644C --- lets add something to HADCURR
645 FORM1= FPIKM(SQRT(SK),AMPI,AMPI) *FPIKM(SQRT(QQ),AMPI,AMPI)
646C FORM1= FPIKM(SQRT(SK),AMPI,AMPI) *FPIKMD(SQRT(QQ),AMPI,AMPI)
647CCCCCCCCCCCCC FORM1=WIGFOR(SK,AMRO,GAMRO) (tests)
648 DO 106 I=1,4
649 DO 106 J=1,4
650 HADCUR(I)=
651 $ HADCUR(I)+CMPLX(COEF1)*FORM1*AA(I,J)*(PP(K,J)-PP(4,J))
652 106 CONTINUE
653C --- end of the non omega current (3 possibilities)
654 101 CONTINUE
655 ENDIF
656 END
657 FUNCTION WIGFOR(S,XM,XGAM)
658 COMPLEX WIGFOR,WIGNOR
659 WIGNOR=CMPLX(-XM**2,XM*XGAM)
660 WIGFOR=WIGNOR/CMPLX(S-XM**2,XM*XGAM)
661 END
662 SUBROUTINE CURINF
663C HERE the form factors of M. Finkemeier et al. start
664C it ends with the string: M. Finkemeier et al. END
665 COMMON /INOUT/ INUT, IOUT
666 WRITE (UNIT = IOUT,FMT = 99)
667 WRITE (UNIT = IOUT,FMT = 98)
668c print *, 'here is curinf'
669 99 FORMAT(
670 . /, ' *************************************************** ',
671 . /, ' you are using the 4 pion decay mode form factors ',
672 . /, ' which have been described in:',
673 . /, ' r. decker, m. finkemeier, p. heiliger and h.h. jonsson',
674 . /, ' "TAU DECAYS INTO FOUR PIONS" ',
675 . /, ' universitaet karlsruhe preprint ttp 94-13 (1994);',
676 . /, ' lnf-94/066(ir); hep-ph/9410260 ',
677 . /, ' ',
678 . /, ' please note that this routine is using parameters',
679 . /, ' related to the 3 pion decay mode(a1 mode), such as',
680 . /, ' the a1 mass and width(taken from the COMMON /parmas/)',
681 . /, ' and the 2 pion vector resonance form factor(by using',
682 . /, ' the routine fpikm)' ,
683 . /, ' thus IF you decide to change any of these, you will' ,
684 . /, ' have to refit the 4 pion parameters in the common' )
685 98 FORMAT(
686 . ' block /tau4pi/, or you might get a bad discription' ,
687 . /, ' of tau -> 4 pions' ,
688 . /, ' for these formfactors set in routine choice for',
689 . /, ' mnum.eq.102 -- amrx=1.42 and gamrx=.21',
690 . /, ' mnum.eq.101 -- amrx=1.3 and gamrx=.46 prob1,prob2=0.2',
691 . /, ' to optimize phase space parametrization',
692 . /, ' *************************************************** ',
693 . /, ' coded by m. finkemeier and p. heiliger, 29. sept. 1994',
694 . /, ' incorporated to tauola by z. was 17. jan. 1995',
695c . /, ' fitted on(day/month/year) by ... ',
696c . /, ' to .... data ',
697 . /, ' changed by: z. was on 17.01.95',
698 . /, ' changes by: m. finkemeier on 30.01.95' )
699 END
700C
701 SUBROUTINE CURINI
702 COMMON /TAU4PI/ GOMEGA,GAMMA1,GAMMA2,ROM1,ROG1,BETA1,
703 . ROM2,ROG2,BETA2
704 REAL*4 GOMEGA,GAMMA1,GAMMA2,ROM1,ROG1,BETA1,
705 . ROM2,ROG2,BETA2
706 GOMEGA = 1.4
707 GAMMA1 = 0.38
708 GAMMA2 = 0.38
709 ROM1 = 1.35
710 ROG1 = 0.3
711 BETA1 = 0.08
712 ROM2 = 1.70
713 ROG2 = 0.235
714 BETA2 = -0.0075
715 END
716 COMPLEX FUNCTION BWIGA1(QA)
717C ================================================================
718C breit-wigner enhancement of a1
719C ================================================================
720 COMPLEX WIGNER
721 COMMON / PARMAS/ AMTAU,AMNUTA,AMEL,AMNUE,AMMU,AMNUMU,
722 % AMPIZ,AMPI,AMRO,GAMRO,AMA1,GAMA1,
723 % AMK,AMKZ,AMKST,GAMKST
724
725C
726 REAL*4 AMTAU,AMNUTA,AMEL,AMNUE,AMMU,AMNUMU,
727 % AMPIZ,AMPI,AMRO,GAMRO,AMA1,GAMA1,
728 % AMK,AMKZ,AMKST,GAMKST
729
730 WIGNER(A,B,C)=CMPLX(1.0,0.0)/CMPLX(A-B**2,B*C)
731 GAMAX=GAMA1*GFUN(QA)/GFUN(AMA1**2)
732 BWIGA1=-AMA1**2*WIGNER(QA,AMA1,GAMAX)
733 RETURN
734 END
735 COMPLEX FUNCTION BWIGEPS(QEPS)
736C =============================================================
737C breit-wigner enhancement of epsilon
738C =============================================================
739 REAL QEPS,MEPS,GEPS
740 MEPS=1.300
741 GEPS=.600
742 BWIGEPS=CMPLX(MEPS**2,-MEPS*GEPS)/
743 % CMPLX(MEPS**2-QEPS,-MEPS*GEPS)
744 RETURN
745 END
746 COMPLEX FUNCTION FRHO4(W,XM1,XM2)
747C ===========================================================
748C rho-type resonance factor with higher radials, to be used
749C by CURR for the four pion mode
750C ===========================================================
751 COMPLEX BWIGM
752 COMMON /TAU4PI/ GOMEGA,GAMMA1,GAMMA2,ROM1,ROG1,BETA1,
753 . ROM2,ROG2,BETA2
754 REAL*4 GOMEGA,GAMMA1,GAMMA2,ROM1,ROG1,BETA1,
755 . ROM2,ROG2,BETA2
756 REAL ROM,ROG,PI,PIM,S,W
757 EXTERNAL BWIG
758 DATA INIT /0/
759C
760C ------------ PARAMETERS --------------------
761.EQ. IF (INIT0 ) THEN
762 INIT=1
763 PI=3.141592654
764 PIM=.140
765 ROM=0.773
766 ROG=0.145
767 ENDIF
768C -----------------------------------------------
769 S=W**2
770c print *,'rom2,rog2 =',rom2,rog2
771 FRHO4=(BWIGM(S,ROM,ROG,XM1,XM2)+BETA1*BWIGM(S,ROM1,ROG1,XM1,XM2)
772 & +BETA2*BWIGM(S,ROM2,ROG2,XM1,XM2))
773 & /(1+BETA1+BETA2)
774 RETURN
775 END
776 SUBROUTINE CURR(MNUM,PIM1,PIM2,PIM3,PIM4,HADCUR)
777C ==================================================================
778C Hadronic current for 4 pi final state, according to:
779C R. Decker, M. Finkemeier, P. Heiliger, H.H.Jonsson, TTP94-13
780C
781C See also:
782C R. Fisher, J. Wess and F. Wagner Z. Phys C3 (1980) 313
783C R. Decker Z. Phys C36 (1987) 487.
784C M. Gell-Mann, D. Sharp, W. Wagner Phys. Rev. Lett 8 (1962) 261.
785C ==================================================================
786
787 COMMON /TAU4PI/ GOMEGA,GAMMA1,GAMMA2,ROM1,ROG1,BETA1,
788 . ROM2,ROG2,BETA2
789 REAL*4 GOMEGA,GAMMA1,GAMMA2,ROM1,ROG1,BETA1,
790 . ROM2,ROG2,BETA2
791 COMMON / PARMAS / AMTAU,AMNUTA,AMEL,AMNUE,AMMU,AMNUMU
792 * ,AMPIZ,AMPI,AMRO,GAMRO,AMA1,GAMA1
793 * ,AMK,AMKZ,AMKST,GAMKST
794C
795 REAL*4 AMTAU,AMNUTA,AMEL,AMNUE,AMMU,AMNUMU
796 * ,AMPIZ,AMPI,AMRO,GAMRO,AMA1,GAMA1
797 * ,AMK,AMKZ,AMKST,GAMKST
798 COMMON / DECPAR / GFERMI,GV,GA,CCABIB,SCABIB,GAMEL
799 REAL*4 GFERMI,GV,GA,CCABIB,SCABIB,GAMEL
800 REAL PIM1(4),PIM2(4),PIM3(4),PIM4(4),PAA(4)
801 COMPLEX HADCUR(4),FORM1,FORM2,FORM3,FPIKM
802 COMPLEX BWIGN,FRHO4
803 COMPLEX BWIGEPS,BWIGA1
804 COMPLEX HCOMP1(4),HCOMP2(4),HCOMP3(4),HCOMP4(4)
805
806 COMPLEX T243,T213,T143,T123,T341,T342
807 COMPLEX T124,T134,T214,T234,T314,T324
808 COMPLEX S2413,S2314,S1423,S1324,S34
809 COMPLEX S2431,S3421
810 COMPLEX BRACK1,BRACK2,BRACK3,BRACK4A,BRACK4B,BRACK4
811
812 REAL QMP1,QMP2,QMP3,QMP4
813 REAL PS43,PS41,PS42,PS34,PS14,PS13,PS24,PS23
814 REAL PS21,PS31
815
816 REAL PD243,PD241,PD213,PD143,PD142
817 REAL PD123,PD341,PD342,PD413,PD423
818 REAL PD124,PD134,PD214,PD234,PD314,PD324
819 REAL QP1,QP2,QP3,QP4
820
821 REAL PA(4),PB(4)
822 REAL AA(4,4),PP(4,4)
823 DATA PI /3.141592653589793238462643/
824 DATA FPI /93.3E-3/
825 DATA INIT /0/
826 BWIGN(A,XM,XG)=1.0/CMPLX(A-XM**2,XM*XG)
827C
828.EQ. IF (INIT0) THEN
829 CALL CURINI
830 CALL CURINF
831 INIT = 1
832 ENDIF
833C
834C --- MASSES AND CONSTANTS
835 G1=12.924
836 G2=1475.98 * GOMEGA
837 G =G1*G2
838 ELPHA=-.1
839 AMROP=1.7
840 GAMROP=0.26
841 AMOM=.782
842 GAMOM=0.0085
843 ARFLAT=1.0
844 AROMEG=1.0
845C
846 FRO=0.266*AMRO**2
847 COEF1=2.0*SQRT(3.0)/FPI**2*ARFLAT
848 COEF2=FRO*G*AROMEG
849C --- INITIALIZATION OF FOUR VECTORS
850 DO 7 K=1,4
851 DO 8 L=1,4
852 8 AA(K,L)=0.0
853 HADCUR(K)=CMPLX(0.0)
854 PAA(K)=PIM1(K)+PIM2(K)+PIM3(K)+PIM4(K)
855 PP(1,K)=PIM1(K)
856 PP(2,K)=PIM2(K)
857 PP(3,K)=PIM3(K)
858 7 PP(4,K)=PIM4(K)
859C
860.EQ. IF (MNUM1) THEN
861C ===================================================================
862C PI- PI- P0 PI+ CASE ====
863C ===================================================================
864 QQ=PAA(4)**2-PAA(3)**2-PAA(2)**2-PAA(1)**2
865
866C FIRST DEFINITION OF SCALAR PRODUCTS OF MOMENTUM VECTORS
867
868C DEFINE (Q-PI)**2 AS QPI:
869
870 QMP1=(PIM2(4)+PIM3(4)+PIM4(4))**2-(PIM2(3)+PIM3(3)+PIM4(3))**2
871 % -(PIM2(2)+PIM3(2)+PIM4(2))**2-(PIM2(1)+PIM3(1)+PIM4(1))**2
872
873 QMP2=(PIM1(4)+PIM3(4)+PIM4(4))**2-(PIM1(3)+PIM3(3)+PIM4(3))**2
874 % -(PIM1(2)+PIM3(2)+PIM4(2))**2-(PIM1(1)+PIM3(1)+PIM4(1))**2
875
876 QMP3=(PIM1(4)+PIM2(4)+PIM4(4))**2-(PIM1(3)+PIM2(3)+PIM4(3))**2
877 % -(PIM1(2)+PIM2(2)+PIM4(2))**2-(PIM1(1)+PIM2(1)+PIM4(1))**2
878
879 QMP4=(PIM1(4)+PIM2(4)+PIM3(4))**2-(PIM1(3)+PIM2(3)+PIM3(3))**2
880 % -(PIM1(2)+PIM2(2)+PIM3(2))**2-(PIM1(1)+PIM2(1)+PIM3(1))**2
881
882
883C DEFINE (PI+PK)**2 AS PSIK:
884
885 PS43=(PIM4(4)+PIM3(4))**2-(PIM4(3)+PIM3(3))**2
886 % -(PIM4(2)+PIM3(2))**2-(PIM4(1)+PIM3(1))**2
887
888 PS41=(PIM4(4)+PIM1(4))**2-(PIM4(3)+PIM1(3))**2
889 % -(PIM4(2)+PIM1(2))**2-(PIM4(1)+PIM1(1))**2
890
891 PS42=(PIM4(4)+PIM2(4))**2-(PIM4(3)+PIM2(3))**2
892 % -(PIM4(2)+PIM2(2))**2-(PIM4(1)+PIM2(1))**2
893
894 PS34=PS43
895
896 PS14=PS41
897
898 PS13=(PIM1(4)+PIM3(4))**2-(PIM1(3)+PIM3(3))**2
899 % -(PIM1(2)+PIM3(2))**2-(PIM1(1)+PIM3(1))**2
900
901 PS24=PS42
902
903 PS23=(PIM2(4)+PIM3(4))**2-(PIM2(3)+PIM3(3))**2
904 % -(PIM2(2)+PIM3(2))**2-(PIM2(1)+PIM3(1))**2
905
906 PD243=PIM2(4)*(PIM4(4)-PIM3(4))-PIM2(3)*(PIM4(3)-PIM3(3))
907 % -PIM2(2)*(PIM4(2)-PIM3(2))-PIM2(1)*(PIM4(1)-PIM3(1))
908
909 PD241=PIM2(4)*(PIM4(4)-PIM1(4))-PIM2(3)*(PIM4(3)-PIM1(3))
910 % -PIM2(2)*(PIM4(2)-PIM1(2))-PIM2(1)*(PIM4(1)-PIM1(1))
911
912 PD213=PIM2(4)*(PIM1(4)-PIM3(4))-PIM2(3)*(PIM1(3)-PIM3(3))
913 % -PIM2(2)*(PIM1(2)-PIM3(2))-PIM2(1)*(PIM1(1)-PIM3(1))
914
915 PD143=PIM1(4)*(PIM4(4)-PIM3(4))-PIM1(3)*(PIM4(3)-PIM3(3))
916 % -PIM1(2)*(PIM4(2)-PIM3(2))-PIM1(1)*(PIM4(1)-PIM3(1))
917
918 PD142=PIM1(4)*(PIM4(4)-PIM2(4))-PIM1(3)*(PIM4(3)-PIM2(3))
919 % -PIM1(2)*(PIM4(2)-PIM2(2))-PIM1(1)*(PIM4(1)-PIM2(1))
920
921 PD123=PIM1(4)*(PIM2(4)-PIM3(4))-PIM1(3)*(PIM2(3)-PIM3(3))
922 % -PIM1(2)*(PIM2(2)-PIM3(2))-PIM1(1)*(PIM2(1)-PIM3(1))
923
924 PD341=PIM3(4)*(PIM4(4)-PIM1(4))-PIM3(3)*(PIM4(3)-PIM1(3))
925 % -PIM3(2)*(PIM4(2)-PIM1(2))-PIM3(1)*(PIM4(1)-PIM1(1))
926
927 PD342=PIM3(4)*(PIM4(4)-PIM2(4))-PIM3(3)*(PIM4(3)-PIM2(3))
928 % -PIM3(2)*(PIM4(2)-PIM2(2))-PIM3(1)*(PIM4(1)-PIM2(1))
929
930 PD413=PIM4(4)*(PIM1(4)-PIM3(4))-PIM4(3)*(PIM1(3)-PIM3(3))
931 % -PIM4(2)*(PIM1(2)-PIM3(2))-PIM4(1)*(PIM1(1)-PIM3(1))
932
933 PD423=PIM4(4)*(PIM2(4)-PIM3(4))-PIM4(3)*(PIM2(3)-PIM3(3))
934 % -PIM4(2)*(PIM2(2)-PIM3(2))-PIM4(1)*(PIM2(1)-PIM3(1))
935
936C DEFINE Q*PI = QPI:
937
938 QP1=PIM1(4)*(PIM1(4)+PIM2(4)+PIM3(4)+PIM4(4))
939 % -PIM1(3)*(PIM1(3)+PIM2(3)+PIM3(3)+PIM4(3))
940 % -PIM1(2)*(PIM1(2)+PIM2(2)+PIM3(2)+PIM4(2))
941 % -PIM1(1)*(PIM1(1)+PIM2(1)+PIM3(1)+PIM4(1))
942
943 QP2=PIM2(4)*(PIM1(4)+PIM2(4)+PIM3(4)+PIM4(4))
944 % -PIM2(3)*(PIM1(3)+PIM2(3)+PIM3(3)+PIM4(3))
945 % -PIM2(2)*(PIM1(2)+PIM2(2)+PIM3(2)+PIM4(2))
946 % -PIM2(1)*(PIM1(1)+PIM2(1)+PIM3(1)+PIM4(1))
947
948 QP3=PIM3(4)*(PIM1(4)+PIM2(4)+PIM3(4)+PIM4(4))
949 % -PIM3(3)*(PIM1(3)+PIM2(3)+PIM3(3)+PIM4(3))
950 % -PIM3(2)*(PIM1(2)+PIM2(2)+PIM3(2)+PIM4(2))
951 % -PIM3(1)*(PIM1(1)+PIM2(1)+PIM3(1)+PIM4(1))
952
953 QP4=PIM4(4)*(PIM1(4)+PIM2(4)+PIM3(4)+PIM4(4))
954 % -PIM4(3)*(PIM1(3)+PIM2(3)+PIM3(3)+PIM4(3))
955 % -PIM4(2)*(PIM1(2)+PIM2(2)+PIM3(2)+PIM4(2))
956 % -PIM4(1)*(PIM1(1)+PIM2(1)+PIM3(1)+PIM4(1))
957
958
959
960C DEFINE T(PI;PJ,PK)= TIJK:
961
962 T243=BWIGA1(QMP2)*FPIKM(SQRT(PS43),AMPI,AMPI)*GAMMA1
963 T213=BWIGA1(QMP2)*FPIKM(SQRT(PS13),AMPI,AMPI)*GAMMA1
964 T143=BWIGA1(QMP1)*FPIKM(SQRT(PS43),AMPI,AMPI)*GAMMA1
965 T123=BWIGA1(QMP1)*FPIKM(SQRT(PS23),AMPI,AMPI)*GAMMA1
966 T341=BWIGA1(QMP3)*FPIKM(SQRT(PS41),AMPI,AMPI)*GAMMA1
967 T342=BWIGA1(QMP3)*FPIKM(SQRT(PS42),AMPI,AMPI)*GAMMA1
968
969C DEFINE S(I,J;K,L)= SIJKL:
970
971 S2413=FRHO4(SQRT(PS24),AMPI,AMPI)*GAMMA2
972 S2314=FRHO4(SQRT(PS23),AMPI,AMPI)*BWIGEPS(PS14)*GAMMA2
973 S1423=FRHO4(SQRT(PS14),AMPI,AMPI)*GAMMA2
974 S1324=FRHO4(SQRT(PS13),AMPI,AMPI)*BWIGEPS(PS24)*GAMMA2
975 S34=FRHO4(SQRT(PS34),AMPI,AMPI)*GAMMA2
976
977C DEFINITION OF AMPLITUDE, FIRST THE [] BRACKETS:
978
979 BRACK1=2.*T143+2.*T243+T123+T213
980 % +T341*(PD241/QMP3-1.)+T342*(PD142/QMP3-1.)
981 % +3./4.*(S1423+S2413-S2314-S1324)-3.*S34
982
983 BRACK2=2.*T143*PD243/QMP1+3.*T213
984 % +T123*(2.*PD423/QMP1+1.)+T341*(PD241/QMP3+3.)
985 % +T342*(PD142/QMP3+1.)
986 % -3./4.*(S2314+3.*S1324+3.*S1423+S2413)
987
988 BRACK3=2.*T243*PD143/QMP2+3.*T123
989 % +T213*(2.*PD413/QMP2+1.)+T341*(PD241/QMP3+1.)
990 % +T342*(PD142/QMP3+3.)
991 % -3./4.*(3.*S2314+S1324+S1423+3.*S2413)
992
993 BRACK4A=2.*T143*(PD243/QQ*(QP1/QMP1+1.)+PD143/QQ)
994 % +2.*T243*(PD143/QQ*(QP2/QMP2+1.)+PD243/QQ)
995 % +T123+T213
996 % +2.*T123*(PD423/QQ*(QP1/QMP1+1.)+PD123/QQ)
997 % +2.*T213*(PD413/QQ*(QP2/QMP2+1.)+PD213/QQ)
998 % +T341*(PD241/QMP3+1.-2.*PD241/QQ*(QP3/QMP3+1.)
999 % -2.*PD341/QQ)
1000 % +T342*(PD142/QMP3+1.-2.*PD142/QQ*(QP3/QMP3+1.)
1001 % -2.*PD342/QQ)
1002
1003 BRACK4B=-3./4.*(S2314*(2.*(QP2-QP3)/QQ+1.)
1004 % +S1324*(2.*(QP1-QP3)/QQ+1.)
1005 % +S1423*(2.*(QP1-QP4)/QQ+1.)
1006 % +S2413*(2.*(QP2-QP4)/QQ+1.)
1007 % +4.*S34*(QP4-QP3)/QQ)
1008
1009 BRACK4=BRACK4A+BRACK4B
1010
1011 DO 208 K=1,4
1012
1013 HCOMP1(K)=(PIM3(K)-PIM4(K))*BRACK1
1014 HCOMP2(K)=PIM1(K)*BRACK2
1015 HCOMP3(K)=PIM2(K)*BRACK3
1016 HCOMP4(K)=(PIM1(K)+PIM2(K)+PIM3(K)+PIM4(K))*BRACK4
1017
1018 208 CONTINUE
1019
1020 DO 209 I=1,4
1021
1022 HADCUR(I)=HCOMP1(I)-HCOMP2(I)-HCOMP3(I)+HCOMP4(I)
1023 HADCUR(I)=-COEF1*FRHO4(SQRT(QQ),AMPI,AMPI)*HADCUR(I)
1024
1025 209 CONTINUE
1026
1027
1028C --- END OF THE NON OMEGA CURRENT (3 POSSIBILITIES)
1029 201 CONTINUE
1030C
1031C
1032C --- THERE ARE TWO POSSIBILITIES FOR OMEGA CURRENT
1033C --- PA PB ARE CORRESPONDING FIRST AND SECOND PI-S
1034 DO 301 KK=1,2
1035 DO 302 I=1,4
1036 PA(I)=PP(KK,I)
1037 PB(I)=PP(3-KK,I)
1038 302 CONTINUE
1039C --- LORENTZ INVARIANTS
1040 QQA=0.0
1041 SS23=0.0
1042 SS24=0.0
1043 SS34=0.0
1044 QP1P2=0.0
1045 QP1P3=0.0
1046 QP1P4=0.0
1047 P1P2 =0.0
1048 P1P3 =0.0
1049 P1P4 =0.0
1050 DO 303 K=1,4
1051 SIGN=-1.0
1052.EQ. IF (K4) SIGN= 1.0
1053 QQA=QQA+SIGN*(PAA(K)-PA(K))**2
1054 SS23=SS23+SIGN*(PB(K) +PIM3(K))**2
1055 SS24=SS24+SIGN*(PB(K) +PIM4(K))**2
1056 SS34=SS34+SIGN*(PIM3(K)+PIM4(K))**2
1057 QP1P2=QP1P2+SIGN*(PAA(K)-PA(K))*PB(K)
1058 QP1P3=QP1P3+SIGN*(PAA(K)-PA(K))*PIM3(K)
1059 QP1P4=QP1P4+SIGN*(PAA(K)-PA(K))*PIM4(K)
1060 P1P2=P1P2+SIGN*PA(K)*PB(K)
1061 P1P3=P1P3+SIGN*PA(K)*PIM3(K)
1062 P1P4=P1P4+SIGN*PA(K)*PIM4(K)
1063 303 CONTINUE
1064C
1065 FORM2=COEF2*(BWIGN(QQ,AMRO,GAMRO)+ELPHA*BWIGN(QQ,AMROP,GAMROP))
1066C FORM3=BWIGN(QQA,AMOM,GAMOM)*(BWIGN(SS23,AMRO,GAMRO)+
1067C $ BWIGN(SS24,AMRO,GAMRO)+BWIGN(SS34,AMRO,GAMRO))
1068 FORM3=BWIGN(QQA,AMOM,GAMOM)
1069C
1070 DO 304 K=1,4
1071 HADCUR(K)=HADCUR(K)+FORM2*FORM3*(
1072 $ PB (K)*(QP1P3*P1P4-QP1P4*P1P3)
1073 $ +PIM3(K)*(QP1P4*P1P2-QP1P2*P1P4)
1074 $ +PIM4(K)*(QP1P2*P1P3-QP1P3*P1P2) )
1075 304 CONTINUE
1076 301 CONTINUE
1077C
1078 ELSE
1079C ===================================================================
1080C PI0 PI0 P0 PI- CASE ====
1081C ===================================================================
1082 QQ=PAA(4)**2-PAA(3)**2-PAA(2)**2-PAA(1)**2
1083
1084
1085C FIRST DEFINITION OF SCALAR PRODUCTS OF MOMENTUM VECTORS
1086
1087C DEFINE (Q-PI)**2 AS QPI:
1088
1089 QMP1=(PIM2(4)+PIM3(4)+PIM4(4))**2-(PIM2(3)+PIM3(3)+PIM4(3))**2
1090 % -(PIM2(2)+PIM3(2)+PIM4(2))**2-(PIM2(1)+PIM3(1)+PIM4(1))**2
1091
1092 QMP2=(PIM1(4)+PIM3(4)+PIM4(4))**2-(PIM1(3)+PIM3(3)+PIM4(3))**2
1093 % -(PIM1(2)+PIM3(2)+PIM4(2))**2-(PIM1(1)+PIM3(1)+PIM4(1))**2
1094
1095 QMP3=(PIM1(4)+PIM2(4)+PIM4(4))**2-(PIM1(3)+PIM2(3)+PIM4(3))**2
1096 % -(PIM1(2)+PIM2(2)+PIM4(2))**2-(PIM1(1)+PIM2(1)+PIM4(1))**2
1097
1098 QMP4=(PIM1(4)+PIM2(4)+PIM3(4))**2-(PIM1(3)+PIM2(3)+PIM3(3))**2
1099 % -(PIM1(2)+PIM2(2)+PIM3(2))**2-(PIM1(1)+PIM2(1)+PIM3(1))**2
1100
1101
1102C DEFINE (PI+PK)**2 AS PSIK:
1103
1104 PS14=(PIM1(4)+PIM4(4))**2-(PIM1(3)+PIM4(3))**2
1105 % -(PIM1(2)+PIM4(2))**2-(PIM1(1)+PIM4(1))**2
1106
1107 PS21=(PIM2(4)+PIM1(4))**2-(PIM2(3)+PIM1(3))**2
1108 % -(PIM2(2)+PIM1(2))**2-(PIM2(1)+PIM1(1))**2
1109
1110 PS23=(PIM2(4)+PIM3(4))**2-(PIM2(3)+PIM3(3))**2
1111 % -(PIM2(2)+PIM3(2))**2-(PIM2(1)+PIM3(1))**2
1112
1113 PS24=(PIM2(4)+PIM4(4))**2-(PIM2(3)+PIM4(3))**2
1114 % -(PIM2(2)+PIM4(2))**2-(PIM2(1)+PIM4(1))**2
1115
1116 PS31=(PIM3(4)+PIM1(4))**2-(PIM3(3)+PIM1(3))**2
1117 % -(PIM3(2)+PIM1(2))**2-(PIM3(1)+PIM1(1))**2
1118
1119 PS34=(PIM3(4)+PIM4(4))**2-(PIM3(3)+PIM4(3))**2
1120 % -(PIM3(2)+PIM4(2))**2-(PIM3(1)+PIM4(1))**2
1121
1122
1123
1124 PD324=PIM3(4)*(PIM2(4)-PIM4(4))-PIM3(3)*(PIM2(3)-PIM4(3))
1125 % -PIM3(2)*(PIM2(2)-PIM4(2))-PIM3(1)*(PIM2(1)-PIM4(1))
1126
1127 PD314=PIM3(4)*(PIM1(4)-PIM4(4))-PIM3(3)*(PIM1(3)-PIM4(3))
1128 % -PIM3(2)*(PIM1(2)-PIM4(2))-PIM3(1)*(PIM1(1)-PIM4(1))
1129
1130 PD234=PIM2(4)*(PIM3(4)-PIM4(4))-PIM2(3)*(PIM3(3)-PIM4(3))
1131 % -PIM2(2)*(PIM3(2)-PIM4(2))-PIM2(1)*(PIM3(1)-PIM4(1))
1132
1133 PD214=PIM2(4)*(PIM1(4)-PIM4(4))-PIM2(3)*(PIM1(3)-PIM4(3))
1134 % -PIM2(2)*(PIM1(2)-PIM4(2))-PIM2(1)*(PIM1(1)-PIM4(1))
1135
1136 PD134=PIM1(4)*(PIM3(4)-PIM4(4))-PIM1(3)*(PIM3(3)-PIM4(3))
1137 % -PIM1(2)*(PIM3(2)-PIM4(2))-PIM1(1)*(PIM3(1)-PIM4(1))
1138
1139 PD124=PIM1(4)*(PIM2(4)-PIM4(4))-PIM1(3)*(PIM2(3)-PIM4(3))
1140 % -PIM1(2)*(PIM2(2)-PIM4(2))-PIM1(1)*(PIM2(1)-PIM4(1))
1141
1142C DEFINE Q*PI = QPI:
1143
1144 QP1=PIM1(4)*(PIM1(4)+PIM2(4)+PIM3(4)+PIM4(4))
1145 % -PIM1(3)*(PIM1(3)+PIM2(3)+PIM3(3)+PIM4(3))
1146 % -PIM1(2)*(PIM1(2)+PIM2(2)+PIM3(2)+PIM4(2))
1147 % -PIM1(1)*(PIM1(1)+PIM2(1)+PIM3(1)+PIM4(1))
1148
1149 QP2=PIM2(4)*(PIM1(4)+PIM2(4)+PIM3(4)+PIM4(4))
1150 % -PIM2(3)*(PIM1(3)+PIM2(3)+PIM3(3)+PIM4(3))
1151 % -PIM2(2)*(PIM1(2)+PIM2(2)+PIM3(2)+PIM4(2))
1152 % -PIM2(1)*(PIM1(1)+PIM2(1)+PIM3(1)+PIM4(1))
1153
1154 QP3=PIM3(4)*(PIM1(4)+PIM2(4)+PIM3(4)+PIM4(4))
1155 % -PIM3(3)*(PIM1(3)+PIM2(3)+PIM3(3)+PIM4(3))
1156 % -PIM3(2)*(PIM1(2)+PIM2(2)+PIM3(2)+PIM4(2))
1157 % -PIM3(1)*(PIM1(1)+PIM2(1)+PIM3(1)+PIM4(1))
1158
1159 QP4=PIM4(4)*(PIM1(4)+PIM2(4)+PIM3(4)+PIM4(4))
1160 % -PIM4(3)*(PIM1(3)+PIM2(3)+PIM3(3)+PIM4(3))
1161 % -PIM4(2)*(PIM1(2)+PIM2(2)+PIM3(2)+PIM4(2))
1162 % -PIM4(1)*(PIM1(1)+PIM2(1)+PIM3(1)+PIM4(1))
1163
1164
1165C DEFINE T(PI;PJ,PK)= TIJK:
1166
1167 T324=BWIGA1(QMP3)*FPIKM(SQRT(PS24),AMPI,AMPI)*GAMMA1
1168 T314=BWIGA1(QMP3)*FPIKM(SQRT(PS14),AMPI,AMPI)*GAMMA1
1169 T234=BWIGA1(QMP2)*FPIKM(SQRT(PS34),AMPI,AMPI)*GAMMA1
1170 T214=BWIGA1(QMP2)*FPIKM(SQRT(PS14),AMPI,AMPI)*GAMMA1
1171 T134=BWIGA1(QMP1)*FPIKM(SQRT(PS34),AMPI,AMPI)*GAMMA1
1172 T124=BWIGA1(QMP1)*FPIKM(SQRT(PS24),AMPI,AMPI)*GAMMA1
1173
1174C DEFINE S(I,J;K,L)= SIJKL:
1175
1176 S1423=FRHO4(SQRT(PS14),AMPI,AMPI)*BWIGEPS(PS23)*GAMMA2
1177 S2431=FRHO4(SQRT(PS24),AMPI,AMPI)*BWIGEPS(PS31)*GAMMA2
1178 S3421=FRHO4(SQRT(PS34),AMPI,AMPI)*BWIGEPS(PS21)*GAMMA2
1179
1180
1181C DEFINITION OF AMPLITUDE, FIRST THE [] BRACKETS:
1182
1183 BRACK1=T234+T324+2.*T314+T134+2.*T214+T124
1184 % +T134*PD234/QMP1+T124*PD324/QMP1
1185 % -3./2.*(S3421+S2431+2.*S1423)
1186
1187
1188 BRACK2=T234*(1.+2.*PD134/QMP2)+3.*T324+3.*T124
1189 % +T134*(1.-PD234/QMP1)+2.*T214*PD314/QMP2
1190 % -T124*PD324/QMP1
1191 % -3./2.*(S3421+3.*S2431)
1192
1193 BRACK3=T324*(1.+2.*PD124/QMP3)+3.*T234+3.*T134
1194 % +T124*(1.-PD324/QMP1)+2.*T314*PD214/QMP3
1195 % -T134*PD234/QMP1
1196 % -3./2.*(3.*S3421+S2431)
1197
1198 BRACK4A=2.*T234*(1./2.+PD134/QQ*(QP2/QMP2+1.)+PD234/QQ)
1199 % +2.*T324*(1./2.+PD124/QQ*(QP3/QMP3+1.)+PD324/QQ)
1200 % +2.*T134*(1./2.+PD234/QQ*(QP1/QMP1+1.)
1201 % -1./2.*PD234/QMP1+PD134/QQ)
1202 % +2.*T124*(1./2.+PD324/QQ*(QP1/QMP1+1.)
1203 % -1./2.*PD324/QMP1+PD124/QQ)
1204 % +2.*T214*(PD314/QQ*(QP2/QMP2+1.)+PD214/QQ)
1205 % +2.*T314*(PD214/QQ*(QP3/QMP3+1.)+PD314/QQ)
1206
1207 BRACK4B=-3./2.*(S3421*(2.*(QP3-QP4)/QQ+1.)
1208 % +S2431*(2.*(QP2-QP4)/QQ+1.)
1209 % +S1423*2.*(QP1-QP4)/QQ)
1210
1211
1212 BRACK4=BRACK4A+BRACK4B
1213
1214 DO 308 K=1,4
1215
1216 HCOMP1(K)=(PIM1(K)-PIM4(K))*BRACK1
1217 HCOMP2(K)=PIM2(K)*BRACK2
1218 HCOMP3(K)=PIM3(K)*BRACK3
1219 HCOMP4(K)=(PIM1(K)+PIM2(K)+PIM3(K)+PIM4(K))*BRACK4
1220
1221 308 CONTINUE
1222
1223 DO 309 I=1,4
1224
1225 HADCUR(I)=HCOMP1(I)+HCOMP2(I)+HCOMP3(I)-HCOMP4(I)
1226 HADCUR(I)=COEF1*FRHO4(SQRT(QQ),AMPI,AMPI)*HADCUR(I)
1227
1228 309 CONTINUE
1229
1230 101 CONTINUE
1231 ENDIF
1232C M. Finkemeier et al. END
1233 END
STL class.