C++InterfacetoTauola
F/prod/pkorb.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  REAL FUNCTION pkorb(IF1,IF2)
35 **********************************************************************
36 *
37 * this function returns a real value
38 * needed in the 1 version of koralb/tauola
39 * corresponding to a mass, width, mixing amplitude, or branching fraction
40 * depending on whether if1 = 1, 2, 3, 4 respectively.
41 * the idea is to make minimal mods to the 3-rd party koralb/tauola code,
42 * so this function supplies all the 1-specific parameters.
43 *
44 * alan weinstein, ajw, 11/97
45 **********************************************************************
46 
47 * arguments:
48  INTEGER if1 ! input, flag for type of data required
49  INTEGER if2 ! input, flag for type of data required
50 
51 * mc info
52 *#include "seq/clinc/qqpars.inc"
53 *#include "seq/clinc/qqprop.inc"
54 *#include "qqlib/seq/qqbrat.inc"
55 
56  INTEGER jak1,jak2,jakp,jakm,ktom
57  COMMON / jaki / jak1,jak2,jakp,jakm,ktom
58  REAL*4 rrr(1)
59  REAL parm(4,100)
60  integer imixpp(300)
61  INTEGER init,i,j
62  REAL c1270,c1402,a1270_kspi,a1270_krho,a1402_kspi,a1402_krho
63  REAL cg1,cg2,r,bra1,brks
64  SAVE init,parm
65  DATA init/0/
66 
67 **********************************************************************
68 * initialize return variable:
69  pkorb = 0.
70 
71 **********************************************************************
72 * initialize:
73  IF (init.EQ.0) THEN
74  init = 1
75 
76 c warning: isospin symmetry enforced, cleo or babar were not enforcing it.
77 c simplification to be used for precision tau decay simulations.
78  bra1=0.0
79  brks=0.0
80 
81 c CALL vzero(parm,400)
82  DO i=1,4
83  DO j=1,100
84  parm(i,j) = 0
85  END DO
86  END DO
87 
88 c youd better be using korb.dec, not decay.dec!!!!
89 c masses(needed in dist/inimas, formf/form*, etc)
90  parm(1, 1) = 1.777000 ! TAU
91  parm(1, 2) = 0. ! NUTA
92  parm(1, 3) = 0.000511 ! EL
93  parm(1, 4) = 0. ! NUEL
94  parm(1, 5) = 0.105658 ! MU
95  parm(1, 6) = 0. ! NUMU
96  parm(1, 7) = 0.134976 ! PIZ
97  parm(1, 8) = 0.139570 ! PI+
98  parm(1, 9) = 0.769900 ! RHO+
99  parm(1,10) = 1.275000 ! A1+
100  parm(1,11) = 0.493677 ! K+
101  parm(1,12) = 0.497670 ! KZ
102  parm(1,13) = 0.891590 ! K*+
103  parm(1,14) = 0.781940 ! OMEG
104  parm(1,15) = 1.370000 ! RHOP+
105  parm(1,16) = 1.700000 ! K*P+
106  parm(1,17) = 1.461000 ! A1P+
107  parm(1,18) = 1.300000 ! PIP+
108  parm(1,19) = 1.270000 ! K1A+
109  parm(1,20) = 1.402000 ! K1B+
110  parm(1,21) = 1.465000 ! RHOPP+
111  parm(1,22) = 1.700000 ! RHOPPP+
112 
113 c widths(needed in dist/inimas, formf/form*, etc)
114  parm(2, 1) = 0. ! TAU
115  parm(2, 2) = 0. ! NUTA
116  parm(2, 3) = 0. ! EL
117  parm(2, 4) = 0. ! NUEL
118  parm(2, 5) = 0. ! MU
119  parm(2, 6) = 0. ! NUMU
120  parm(2, 7) = 0. ! PIZ
121  parm(2, 8) = 0. ! PI+
122  parm(2, 9) = 0.1512 ! RHO+
123  parm(2,10) = 0.700 ! A1+
124  parm(2,11) = 0. ! K+
125  parm(2,12) = 0. ! KZ
126  parm(2,13) = 0.0498 ! K*+
127  parm(2,14) = 0.00843 ! OMEG
128  parm(2,15) = 0.510 ! RHOP+
129  parm(2,16) = 0.235 ! K*P+
130  parm(2,17) = 0.250 ! A1P+
131  parm(2,18) = 0.400 ! PIP+
132  parm(2,19) = 0.090 ! K1A+
133  parm(2,20) = 0.174 ! K1B+
134  parm(2,21) = 0.310 ! RHOPP+
135  parm(2,22) = 0.235 ! RHOPPP+
136 
137 c now store mixing parameters for 2pi and 4pi ffs
138 c needed in tauola/fpik, tauola/bwigs, formf/form* , formf/curr :
139 
140  parm(3,15) = -0.145
141 
142  imixpp(205)=1
143  imixpp(207)=1
144  imixpp(209)=1
145  imixpp(211)=1
146  imixpp(201)=1
147  imixpp(203)=1
148  imixpp(213)=1
149  imixpp(215)=1
150 
151 
152  IF (imixpp(205).NE.0) parm(3,15) = -0.110
153  IF (imixpp(207).NE.0) parm(3,16) = -0.038
154  IF (imixpp(209).NE.0) parm(3,17) = 0.00
155  IF (imixpp(211).NE.0) parm(3,18) = 0.00
156  IF (imixpp(201).NE.0) parm(3,19) = 1.0
157  IF (imixpp(203).NE.0) parm(3,20) = 0.8
158  IF (imixpp(213).NE.0) parm(3,21) = -0.110
159  IF (imixpp(215).NE.0) parm(3,22) = -0.110
160 
161  print *,' KORB: rho/rhop -> pi-pi0 mixing:'
162  print *,' KORB: rho =',parm(1,9) ,parm(2,9)
163  print *,' KORB: rhop =',parm(1,15),parm(2,15),parm(3,15)
164  print *,' KORB: K*/K*prime -> Kpi mixing:'
165  print *,' KORB: kstp =',parm(1,16),parm(2,16),parm(3,16)
166  print *,' KORB: a1/a1prime -> 3pi, KKpi mixing:'
167  print *,' KORB: a1 =',parm(1,10),parm(2,10)
168  print *,' KORB: a1prim=',parm(1,17),parm(2,17),parm(3,17)
169  print *,' KORB: K1A/K1B -> Kpipi mixing:'
170  print *,' KORB: K1A =',parm(1,19),parm(2,19),parm(3,19)
171  print *,' KORB: K1B =',parm(1,20),parm(2,20),parm(3,20)
172  print *,' KORB: rho/rhop/rhopp -> 4pi mixing:'
173  print *,' KORB: rho =',parm(1,9) ,parm(2,9)
174  print *,' KORB: rhopp =',parm(1,21),parm(2,21),parm(3,21)
175  print *,' KORB: rhoppp=',parm(1,22),parm(2,22),parm(3,22)
176 
177 c amplitudes for curr_cleo.f:
178 c for(3pi)-pi0: 4pi phase space; rho0pi-pi0; rho-pi+pi-; rho+pi-pi-; pi-omega
179  parm(3,31) = 0.
180  parm(3,32) = 0.1242
181  parm(3,33) = 0.1604
182  parm(3,34) = 0.2711
183  parm(3,35) = 0.4443
184 c for pi-3pi0: 4pi phase space; rho-pi0pi0
185  parm(3,36) = 0.
186  parm(3,37) = 1.0
187 
188 c modify amplitudes for 4pi form-factor in formf/curr, from korb.dec:
189 ccc IF (iplist(2,282).EQ.5) THEN
190  iplist=0
191  IF (iplist.EQ.5) THEN
192  parm(3,31) = 0.0000
193  parm(3,32) = 0.1242
194  parm(3,33) = 0.1604
195  parm(3,34) = 0.2711
196  parm(3,35) = 0.4443
197  parm(3,36) = 0.0000
198  parm(3,37) = 1.0000
199  END IF
200 
201  print *,' KORB: 3PI-PI0 PARAMS:',(parm(3,i),i=31,35)
202  print *,' KORB: PI-3PI0 PARAMS:',(parm(3,i),i=36,37)
203 
204 c the 4pi models are the most complicated in tauola.
205 c If the user has not modified any parameters of the 4pi model,
206 c we can use the wtmax determined with many trials.
207  IF (abs(parm(3,31)-0.0000).GT.0.0001 .OR.
208  1 abs(parm(3,32)-0.1242).GT.0.0001 .OR.
209  1 abs(parm(3,33)-0.1604).GT.0.0001 .OR.
210  1 abs(parm(3,34)-0.2711).GT.0.0001 .OR.
211  1 abs(parm(3,35)-0.4443).GT.0.0001 ) THEN
212  parm(3,38) = -1.
213  ELSE
214  parm(3,38) = 6.9673671e-14
215  END IF
216 
217  IF (abs(parm(3,36)-0.0000).GT.0.0001 .OR.
218  1 abs(parm(3,37)-1.0000).GT.0.0001 ) THEN
219  parm(3,39) = -1.
220  ELSE
221  parm(3,39) = 3.5374880e-13
222  END IF
223 
224 
225 c phases for curr_cleo.f:
226  parm(3,42) = -0.40
227  parm(3,43) = 0.00
228  parm(3,44) = -0.20+3.1416
229  parm(3,45) = -1.50
230 
231 c rho' contributions to rho' -> pi-omega:
232  parm(3,51) = -0.10
233  parm(3,52) = 1.00
234  parm(3,53) = -0.10
235  parm(3,54) = -0.04
236 
237 c rho' contribtions to rho' -> rhopipi:
238  parm(3,55) = 1.00
239  parm(3,56) = 0.14
240  parm(3,57) = -0.05
241  parm(3,58) = -0.05
242 
243 c rho contributions to rhopipi, rho -> 2pi:
244  parm(3,59) = 1.000
245  parm(3,60) = -0.145
246  parm(3,61) = 0.000
247 
248 c set the brs for(a1+ -> rho+ pi0) and(k*+ -> k0 pi+)
249 c needed in dist/taurdf:
250  parm(4,1) = 0.4920 ! BRA1+
251  parm(4,2) = 0.4920 ! BRA1-
252  parm(4,3) = 0.6660 ! BRKS+
253  parm(4,4) = 0.6660 ! BRKS-
254  parm(4,5) = 0.5 ! BRK0
255  parm(4,6) = 0.5 ! BRK0B
256 
257 c amplitude coefficients for tau -> k1(1270) / k1(1402)
258  c1270 = parm(3,19)
259  c1402 = parm(3,20)
260  IF (c1270.EQ.0.AND.c1402.EQ.0.) THEN
261  c1270 = 1.
262  c1402 = 0.6
263  END IF
264 c from pdg96, square roots of branching fractions:
265  a1270_kspi = sqrt(0.16)
266  a1270_krho = sqrt(0.42)
267  a1402_kspi = sqrt(0.94)
268  a1402_krho = sqrt(0.03)
269 c c-g coefficients for k1- -> cg1 * |k- pi0> + cg2 * |k0bar pi->
270  cg1 = -sqrt(2./3.)
271  cg2 = sqrt(1./3.)
272 c and the resulting amplitudes(times normalized ff):
273  parm(3,81) = c1270*a1270_kspi*cg1 ! K1270 -> K*0B pi-
274  parm(3,82) = c1402*a1402_kspi*cg1 ! K1402 -> K*0B pi-
275  parm(3,83) = c1270*a1270_krho*cg1 ! K1270 -> K0B rho-
276  parm(3,84) = c1402*a1402_krho*cg1 ! K1402 -> K0B rho-
277  parm(3,85) = c1270*a1270_kspi*cg2 ! K1270 -> K*- pi0
278  parm(3,86) = c1402*a1402_kspi*cg2 ! K1402 -> K*- pi0
279  parm(3,87) = c1270*a1270_krho*cg2 ! K1270 -> K- rho0
280  parm(3,88) = c1402*a1402_krho*cg2 ! K1402 -> K- rho0
281 
282  END IF
283 **********************************************************************
284 
285  r = 0.
286  IF (if1.GE.1 .AND. if1.LE.4 .AND. if2.GE.1 .AND. if2.LE.100) THEN
287  r = parm(if1,if2)
288 
289 cajw 4/4/94 better to decide on a1 br now, avoid dadmaa/dphsaa problem.
290  IF (if1.EQ.4.AND.jak1.EQ.5) THEN
291  IF (if2.EQ.11) THEN
292 c Return the br used in the last call:
293  r = bra1
294  ELSE IF (if2.EQ.1) THEN
295  bra1 = r
296  CALL ranmar(rrr,1)
297  IF (rrr(1).LT.bra1) THEN
298  r = 1. ! 3pi
299  ELSE
300  r = 0. ! pi-2pi0
301  END IF
302  bra1 = r
303  END IF
304  ELSEIF (if1.EQ.4.AND.jak1.EQ.7) THEN
305  IF (if2.EQ.13) THEN
306 c Return the br used in the last call:
307  r = brks
308  ELSE IF (if2.EQ.3) THEN
309  brks = r
310  CALL ranmar(rrr,1)
311  IF (rrr(1).LT.brks) THEN
312  r = 1. ! K0 pi-
313  ELSE
314  r = 0. ! K- pi0
315  END IF
316  brks = r
317  END IF
318  END IF
319 
320  END IF
321 
322  pkorb = r
323  RETURN
324  END