C++InterfacetoTauola
demo-factory/prod/tauola_photos_ini.F
1 C this file is created by hand from taumain.F
2 C actions: Remove routines: TAUDEM DECTES TAUFIL FILHEP
3 C add: INIETC will not necesarily work fine ...
4 C replace TRALO4
5 C rename INIPHY to INIPHX
6 
7  SUBROUTINE inietc(jakk1,jakk2,itd,ifpho)
8  COMMON / idfc / idff
9  COMMON / taurad / xk0dec,itdkrc
10  DOUBLE PRECISION XK0DEC
11  COMMON / jaki / jak1,jak2,jakp,jakm,ktom
12  COMMON /phoact/ ifphot
13  SAVE
14 C KTO=1 will denote tau+, thus :: IDFF=-15
15  idff=-15
16 C XK0 for tau decays.
17  xk0dec=0.01
18 C radiative correction switch in tau --> e (mu) decays !
19  itdkrc=itd
20 C switches of tau+ tau- decay modes !!
21  jak1=jakk1
22  jak2=jakk2
23 C photos activation switch
24  ifphot=ifpho
25  end
26 
27  SUBROUTINE tralo4(KTOS,PHOI,PHOF,AM)
28 !! Corrected 11.10.96 (ZW) tralor for KORALW.
29 !! better treatment is to cascade from tau rest-frame through W
30 !! restframe down to LAB.
31  COMMON / momdec / q1,q2,p1,p2,p3,p4
32  COMMON /tralid/ idtra
33  double precision Q1(4),Q2(4),P1(4),P2(4),P3(4),P4(4),P1QQ(4),P2QQ(4)
34  double precision PIN(4),POUT(4),PBST(4),PBS1(4),QQ(4),PI
35  double precision THET,PHI,EXE
36  REAL*4 PHOI(4),PHOF(4)
37  SAVE
38  DATA pi /3.141592653589793238462643d0/
39  am=sqrt(abs
40  $ (phoi(4)**2-phoi(3)**2-phoi(2)**2-phoi(1)**2))
41  idtra=ktos
42  DO k=1,4
43  pin(k)=phoi(k)
44  phof(k)=phoi(k)
45  ENDDO
46 ! write(*,*) idtra
47  IF (idtra.EQ.1) THEN
48  DO k=1,4
49  pbst(k)=p1(k)
50  qq(k)=q1(k)
51  ENDDO
52  ELSEIF(idtra.EQ.2) THEN
53  DO k=1,4
54  pbst(k)=p2(k)
55  qq(k)=q1(k)
56  ENDDO
57  ELSEIF(idtra.EQ.3) THEN
58  DO k=1,4
59  pbst(k)=p3(k)
60  qq(k)=q2(k)
61  ENDDO
62  ELSE
63  DO k=1,4
64  pbst(k)=p4(k)
65  qq(k)=q2(k)
66  ENDDO
67  ENDIF
68 
69 
70 
71  CALL bostdq(1,qq,pbst,pbst)
72  CALL bostdq(1,qq,p1,p1qq)
73  CALL bostdq(1,qq,p2,p2qq)
74  pbs1(4)=pbst(4)
75  pbs1(3)=sqrt(pbst(3)**2+pbst(2)**2+pbst(1)**2)
76  pbs1(2)=0d0
77  pbs1(1)=0d0
78  exe=(pbs1(4)+pbs1(3))/dsqrt(pbs1(4)**2-pbs1(3)**2)
79 C for KTOS=1 boost is antiparallel to 4-momentum of P2.
80 C restframes of tau+ tau- and 'first' frame of 'higgs' are all connected
81 C by boosts along z axis
82  IF(ktos.EQ.1) exe=(pbs1(4)-pbs1(3))/dsqrt(pbs1(4)**2-pbs1(3)**2)
83  CALL bostd3(exe,pin,pout)
84 
85 C once in Z/gamma/Higgs rest frame we control further kinematics by P2QQ for KTOS=1,2
86  thet=acos(p2qq(3)/sqrt(p2qq(3)**2+p2qq(2)**2+p2qq(1)**2))
87  phi=0d0
88  phi=acos(p2qq(1)/sqrt(p2qq(2)**2+p2qq(1)**2))
89  IF(p2qq(2).LT.0d0) phi=2*pi-phi
90 
91  CALL rotpox(thet,phi,pout)
92  CALL bostdq(-1,qq,pout,pout)
93  DO k=1,4
94  phof(k)=pout(k)
95  ENDDO
96  END
97 
98  SUBROUTINE choice(MNUM,RR,ICHAN,PROB1,PROB2,PROB3,
99  $ AMRX,GAMRX,AMRA,GAMRA,AMRB,GAMRB)
100  COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
101  * ,ampiz,ampi,amro,gamro,ama1,gama1
102  * ,amk,amkz,amkst,gamkst
103 C
104  REAL*4 AMTAU,AMNUTA,AMEL,AMNUE,AMMU,AMNUMU
105  * ,ampiz,ampi,amro,gamro,ama1,gama1
106  * ,amk,amkz,amkst,gamkst
107 C
108  amrop=1.1
109  gamrop=0.36
110  amom=.782
111  gamom=0.0084
112 C XXXXA CORRESPOND TO S2 CHANNEL !
113  IF(mnum.EQ.0) THEN
114  prob1=0.5
115  prob2=0.5
116  amrx =ama1
117  gamrx=gama1
118  amra =amro
119  gamra=gamro
120  amrb =amro
121  gamrb=gamro
122  ELSEIF(mnum.EQ.1) THEN
123  prob1=0.5
124  prob2=0.5
125  amrx =1.57
126  gamrx=0.9
127  amrb =amkst
128  gamrb=gamkst
129  amra =amro
130  gamra=gamro
131  ELSEIF(mnum.EQ.2) THEN
132  prob1=0.5
133  prob2=0.5
134  amrx =1.57
135  gamrx=0.9
136  amrb =amkst
137  gamrb=gamkst
138  amra =amro
139  gamra=gamro
140  ELSEIF(mnum.EQ.3) THEN
141  prob1=0.5
142  prob2=0.5
143  amrx =1.27
144  gamrx=0.3
145  amra =amkst
146  gamra=gamkst
147  amrb =amkst
148  gamrb=gamkst
149  ELSEIF(mnum.EQ.4) THEN
150  prob1=0.5
151  prob2=0.5
152  amrx =1.27
153  gamrx=0.3
154  amra =amkst
155  gamra=gamkst
156  amrb =amkst
157  gamrb=gamkst
158  ELSEIF(mnum.EQ.5) THEN
159  prob1=0.5
160  prob2=0.5
161  amrx =1.27
162  gamrx=0.3
163  amra =amkst
164  gamra=gamkst
165  amrb =amro
166  gamrb=gamro
167  ELSEIF(mnum.EQ.6) THEN
168  prob1=0.4
169  prob2=0.4
170  amrx =1.27
171  gamrx=0.3
172  amra =amro
173  gamra=gamro
174  amrb =amkst
175  gamrb=gamkst
176  ELSEIF(mnum.EQ.7) THEN
177  prob1=0.0
178  prob2=1.0
179  amrx =1.27
180  gamrx=0.9
181  amra =amro
182  gamra=gamro
183  amrb =amro
184  gamrb=gamro
185  ELSEIF(mnum.EQ.8) THEN
186  prob1=0.0
187  prob2=1.0
188  amrx =amrop
189  gamrx=gamrop
190  amrb =amom
191  gamrb=gamom
192  amra =amro
193  gamra=gamro
194  ELSEIF(mnum.EQ.101) THEN
195  prob1=.35
196  prob2=.35
197  amrx =1.2
198  gamrx=.46
199  amrb =amom
200  gamrb=gamom
201  amra =amom
202  gamra=gamom
203  ELSEIF(mnum.EQ.102) THEN
204  prob1=0.0
205  prob2=0.0
206  amrx =1.4
207  gamrx=.6
208  amrb =amom
209  gamrb=gamom
210  amra =amom
211  gamra=gamom
212  ELSE
213  prob1=0.0
214  prob2=0.0
215  amrx =ama1
216  gamrx=gama1
217  amra =amro
218  gamra=gamro
219  amrb =amro
220  gamrb=gamro
221  ENDIF
222 C
223  IF (rr.LE.prob1) THEN
224  ichan=1
225  ELSEIF(rr.LE.(prob1+prob2)) THEN
226  ichan=2
227  ax =amra
228  gx =gamra
229  amra =amrb
230  gamra=gamrb
231  amrb =ax
232  gamrb=gx
233  px =prob1
234  prob1=prob2
235  prob2=px
236  ELSE
237  ichan=3
238  ENDIF
239 C
240  prob3=1.0-prob1-prob2
241  END
242  SUBROUTINE initdk
243 * ----------------------------------------------------------------------
244 * INITIALISATION OF TAU DECAY PARAMETERS and routines
245 *
246 * called by : KORALZ
247 * ----------------------------------------------------------------------
248 
249  COMMON / decpar / gfermi,gv,ga,ccabib,scabib,gamel
250  REAL*4 GFERMI,GV,GA,CCABIB,SCABIB,GAMEL
251  COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
252  * ,ampiz,ampi,amro,gamro,ama1,gama1
253  * ,amk,amkz,amkst,gamkst
254 *
255  REAL*4 AMTAU,AMNUTA,AMEL,AMNUE,AMMU,AMNUMU
256  * ,ampiz,ampi,amro,gamro,ama1,gama1
257  * ,amk,amkz,amkst,gamkst
258  COMMON / taubra / gamprt(30),jlist(30),nchan
259  COMMON / taukle / bra1,brk0,brk0b,brks
260  REAL*4 BRA1,BRK0,BRK0B,BRKS
261 #if defined (ALEPH)
262  parameter(nmode=15,nm1=0,nm2=1,nm3=8,nm4=2,nm5=1,nm6=3)
263  COMMON / taudcd /idffin(9,nmode),mulpik(nmode)
264  & ,names
265  CHARACTER NAMES(NMODE)*31
266 #else
267  parameter(nmode=15,nm1=0,nm2=1,nm3=8,nm4=2,nm5=1,nm6=3)
268  COMMON / taudcd /idffin(9,nmode),mulpik(nmode)
269  & ,names
270  CHARACTER NAMES(NMODE)*31
271 #endif
272  CHARACTER OLDNAMES(7)*31
273  CHARACTER*80 bxINIT
274  parameter(
275  $ bxinit ='(1x,1h*,g17.8, 16x, a31,a4,a4, 1x,1h*)'
276  $ )
277  REAL*4 PI,POL1(4)
278 *
279 *
280 * LIST OF BRANCHING RATIOS
281 CAM normalised to e nu nutau channel
282 CAM enu munu pinu rhonu A1nu Knu K*nu pi
283 CAM DATA JLIST / 1, 2, 3, 4, 5, 6, 7,
284 #if defined (ALEPH)
285 CAM /0.1779,0.1731,0.1106,0.2530,0.1811,0.0072,0.0139
286 CAM DATA GAMPRT / 1.000,0.9732,0.6217,1.4221,1.0180,0.0405,0.0781
287 CAM DATA GAMPRT /1.000,0.9676,0.6154,1.3503,1.0225,0.0368,O.O758
288 CAM
289 C
290 C conventions of particles names
291 c
292 cam mode (JAK) 8 9
293 CAM channel pi- pi- pi0 pi+ 3pi0 pi-
294 cam particle code -1,-1, 2, 1, 0, 0, 2, 2, 2,-1, 0, 0,
295 CAM BR relative to electron .2414, .0601,
296 c
297 * 10 11
298 * 1 3pi+- 2pi0 5pi+-
299 * 1 -1,-1, 1, 2, 2, 0, -1,-1,-1, 1, 1, 0,
300 * 1 .0281, .0045,
301 
302 * 12 13
303 * 2 5pi+- pi0 3pi+- 3pi0
304 * 2 -1,-1,-1, 1, 1, 2, -1,-1, 1, 2, 2, 2,
305 * 2 .0010, .0062,
306 
307 * 14 15
308 * 3 K- pi- K+ K0 pi- KB
309 * 3 -3,-1, 3, 0, 0, 0, 4,-1,-4, 0, 0, 0,
310 * 3 .0096, .0169,
311 
312 * 16 17
313 * 4 K- pi0 K0 2pi0 K-
314 * 4 -3, 2, 4, 0, 0, 0, 2, 2,-3, 0, 0, 0,
315 * 4 .0056, .0045,
316 
317 * 18 19
318 * 5 K- pi- pi+ pi- KB pi0
319 * 5 -3,-1, 1, 0, 0, 0, -1,-4, 2, 0, 0, 0,
320 * 5 .0219, .0180,
321 
322 * 20 21
323 * 6 eta pi- pi0 pi- pi0 gamma
324 * 6 9,-1, 2, 0, 0, 0, -1, 2, 8, 0, 0, 0,
325 * 6 .0096, .0088,
326 
327 * 22 /
328 * 7 K- K0 /
329 * 7 -3, 4 /
330 * 7 .0146 /
331 #else
332 *AM DATA GAMPRT /1.000,0.9730,0.6054,1.2432,0.8432,0.0432,O.O811,0.616
333 *AM
334 *AM multipion decays
335 *
336 * conventions of particles names
337 * K-,P-,K+, K0,P-,KB, K-,P0,K0
338 * 3, 1,-3 , 4, 1,-4 , 3, 2, 4 ,
339 * P0,P0,K-, K-,P-,P+, P-,KB,P0
340 * 2, 2, 3 , 3, 1,-1 , 1,-4, 2 ,
341 * ET,P-,P0 P-,P0,GM
342 * 9, 1, 2 , 1, 2, 8
343 *
344 #endif
345 C
346  dimension nopik(6,nmode),npik(nmode)
347 *AM outgoing multiplicity and flavors of multi-pion /multi-K modes
348  DATA npik / 4, 4,
349  1 5, 5,
350  2 6, 6,
351  3 3, 3,
352  4 3, 3,
353  5 3, 3,
354  6 3, 3,
355  7 2 /
356 #if defined (ALEPH)
357  DATA nopik / -1,-1, 2, 1, 0, 0, 2, 2, 2,-1, 0, 0,
358  1 -1,-1, 1, 2, 2, 0, -1,-1,-1, 1, 1, 0,
359  2 -1,-1,-1, 1, 1, 2, -1,-1, 1, 2, 2, 2,
360  3 -3,-1, 3, 0, 0, 0, 4,-1,-4, 0, 0, 0,
361  4 -3, 2, 4, 0, 0, 0, 2, 2,-3, 0, 0, 0,
362  5 -3,-1, 1, 0, 0, 0, -1,-4, 2, 0, 0, 0,
363  6 9,-1, 2, 0, 0, 0, -1, 2, 8, 0, 0, 0,
364 #else
365  DATA nopik / -1,-1, 1, 2, 0, 0, 2, 2, 2,-1, 0, 0,
366  1 -1,-1, 1, 2, 2, 0, -1,-1,-1, 1, 1, 0,
367  2 -1,-1,-1, 1, 1, 2, -1,-1, 1, 2, 2, 2,
368  3 -3,-1, 3, 0, 0, 0, -4,-1, 4, 0, 0, 0,
369  4 -3, 2,-4, 0, 0, 0, 2, 2,-3, 0, 0, 0,
370  5 -3,-1, 1, 0, 0, 0, -1, 4, 2, 0, 0, 0,
371  6 9,-1, 2, 0, 0, 0, -1, 2, 8, 0, 0, 0,
372 #endif
373 #if defined (CLEO)
374 C AJWMOD fix sign bug, 2/22/99
375  7 -3,-4, 0, 0, 0, 0 /
376 #else
377  7 -3, 4, 0, 0, 0, 0 /
378 #endif
379 * LIST OF BRANCHING RATIOS
380  nchan = nmode + 7
381  DO 1 i = 1,30
382  IF (i.LE.nchan) THEN
383  jlist(i) = i
384 #if defined (CePeCe)
385  IF(i.EQ. 1) gamprt(i) = 1.0000
386  IF(i.EQ. 2) gamprt(i) = 1.0000
387  IF(i.EQ. 3) gamprt(i) = 1.0000
388  IF(i.EQ. 4) gamprt(i) = 1.0000
389  IF(i.EQ. 5) gamprt(i) = 1.0000
390  IF(i.EQ. 6) gamprt(i) = 1.0000
391  IF(i.EQ. 7) gamprt(i) = 1.0000
392  IF(i.EQ. 8) gamprt(i) = 1.0000
393  IF(i.EQ. 9) gamprt(i) = 1.0000
394  IF(i.EQ.10) gamprt(i) = 1.0000
395  IF(i.EQ.11) gamprt(i) = 1.0000
396  IF(i.EQ.12) gamprt(i) = 1.0000
397  IF(i.EQ.13) gamprt(i) = 1.0000
398  IF(i.EQ.14) gamprt(i) = 1.0000
399  IF(i.EQ.15) gamprt(i) = 1.0000
400  IF(i.EQ.16) gamprt(i) = 1.0000
401  IF(i.EQ.17) gamprt(i) = 1.0000
402  IF(i.EQ.18) gamprt(i) = 1.0000
403  IF(i.EQ.19) gamprt(i) = 1.0000
404  IF(i.EQ.20) gamprt(i) = 1.0000
405  IF(i.EQ.21) gamprt(i) = 1.0000
406  IF(i.EQ.22) gamprt(i) = 1.0000
407 #elif defined (CLEO)
408  IF(i.EQ. 1) gamprt(i) =0.1800
409  IF(i.EQ. 2) gamprt(i) =0.1751
410  IF(i.EQ. 3) gamprt(i) =0.1110
411  IF(i.EQ. 4) gamprt(i) =0.2515
412  IF(i.EQ. 5) gamprt(i) =0.1790
413  IF(i.EQ. 6) gamprt(i) =0.0071
414  IF(i.EQ. 7) gamprt(i) =0.0134
415  IF(i.EQ. 8) gamprt(i) =0.0450
416  IF(i.EQ. 9) gamprt(i) =0.0100
417  IF(i.EQ.10) gamprt(i) =0.0009
418  IF(i.EQ.11) gamprt(i) =0.0004
419  IF(i.EQ.12) gamprt(i) =0.0003
420  IF(i.EQ.13) gamprt(i) =0.0005
421  IF(i.EQ.14) gamprt(i) =0.0015
422  IF(i.EQ.15) gamprt(i) =0.0015
423  IF(i.EQ.16) gamprt(i) =0.0015
424  IF(i.EQ.17) gamprt(i) =0.0005
425  IF(i.EQ.18) gamprt(i) =0.0050
426  IF(i.EQ.19) gamprt(i) =0.0055
427  IF(i.EQ.20) gamprt(i) =0.0017
428  IF(i.EQ.21) gamprt(i) =0.0013
429  IF(i.EQ.22) gamprt(i) =0.0010
430 #elif defined (ALEPH)
431  IF(i.EQ. 1) gamprt(i) = 1.0000
432  IF(i.EQ. 2) gamprt(i) = .9732
433  IF(i.EQ. 3) gamprt(i) = .6217
434  IF(i.EQ. 4) gamprt(i) = 1.4221
435  IF(i.EQ. 5) gamprt(i) = 1.0180
436  IF(i.EQ. 6) gamprt(i) = .0405
437  IF(i.EQ. 7) gamprt(i) = .0781
438  IF(i.EQ. 8) gamprt(i) = .2414
439  IF(i.EQ. 9) gamprt(i) = .0601
440  IF(i.EQ.10) gamprt(i) = .0281
441  IF(i.EQ.11) gamprt(i) = .0045
442  IF(i.EQ.12) gamprt(i) = .0010
443  IF(i.EQ.13) gamprt(i) = .0062
444  IF(i.EQ.14) gamprt(i) = .0096
445  IF(i.EQ.15) gamprt(i) = .0169
446  IF(i.EQ.16) gamprt(i) = .0056
447  IF(i.EQ.17) gamprt(i) = .0045
448  IF(i.EQ.18) gamprt(i) = .0219
449  IF(i.EQ.19) gamprt(i) = .0180
450  IF(i.EQ.20) gamprt(i) = .0096
451  IF(i.EQ.21) gamprt(i) = .0088
452  IF(i.EQ.22) gamprt(i) = .0146
453 #else
454 #endif
455  IF(i.EQ. 1) oldnames(i)=' TAU- --> E- '
456  IF(i.EQ. 2) oldnames(i)=' TAU- --> MU- '
457  IF(i.EQ. 3) oldnames(i)=' TAU- --> PI- '
458  IF(i.EQ. 4) oldnames(i)=' TAU- --> PI-, PI0 '
459  IF(i.EQ. 5) oldnames(i)=' TAU- --> A1- (two subch) '
460  IF(i.EQ. 6) oldnames(i)=' TAU- --> K- '
461  IF(i.EQ. 7) oldnames(i)=' TAU- --> K*- (two subch) '
462  IF(i.EQ. 8) names(i-7)=' TAU- --> 2PI-, PI0, PI+ '
463  IF(i.EQ. 9) names(i-7)=' TAU- --> 3PI0, PI- '
464  IF(i.EQ.10) names(i-7)=' TAU- --> 2PI-, PI+, 2PI0 '
465  IF(i.EQ.11) names(i-7)=' TAU- --> 3PI-, 2PI+, '
466  IF(i.EQ.12) names(i-7)=' TAU- --> 3PI-, 2PI+, PI0 '
467  IF(i.EQ.13) names(i-7)=' TAU- --> 2PI-, PI+, 3PI0 '
468  IF(i.EQ.14) names(i-7)=' TAU- --> K-, PI-, K+ '
469  IF(i.EQ.15) names(i-7)=' TAU- --> K0, PI-, K0B '
470 #if defined (ALEPH)
471  IF(i.EQ.16) names(i-7)=' TAU- --> K- PI0 K0 '
472 #else
473  IF(i.EQ.16) names(i-7)=' TAU- --> K-, K0, PI0 '
474 #endif
475  IF(i.EQ.17) names(i-7)=' TAU- --> PI0 PI0 K- '
476  IF(i.EQ.18) names(i-7)=' TAU- --> K- PI- PI+ '
477  IF(i.EQ.19) names(i-7)=' TAU- --> PI- K0B PI0 '
478  IF(i.EQ.20) names(i-7)=' TAU- --> ETA PI- PI0 '
479  IF(i.EQ.21) names(i-7)=' TAU- --> PI- PI0 GAM '
480  IF(i.EQ.22) names(i-7)=' TAU- --> K- K0 '
481  ELSE
482  jlist(i) = 0
483  gamprt(i) = 0.
484  ENDIF
485  1 CONTINUE
486  DO i=1,nmode
487  mulpik(i)=npik(i)
488  DO j=1,mulpik(i)
489  idffin(j,i)=nopik(j,i)
490  ENDDO
491  ENDDO
492 *
493 *
494 * --- COEFFICIENTS TO FIX RATIO OF:
495 * --- A1 3CHARGED/ A1 1CHARGED 2 NEUTRALS MATRIX ELEMENTS (MASLESS LIM.)
496 * --- PROBABILITY OF K0 TO BE KS
497 * --- PROBABILITY OF K0B TO BE KS
498 * --- RATIO OF COEFFICIENTS FOR K*--> K0 PI-
499 * --- ALL COEFFICENTS SHOULD BE IN THE RANGE (0.0,1.0)
500 * --- THEY MEANING IS PROBABILITY OF THE FIRST CHOICE ONLY IF ONE
501 * --- NEGLECTS MASS-PHASE SPACE EFFECTS
502  bra1=0.5
503  brk0=0.5
504  brk0b=0.5
505  brks=0.6667
506 *
507 
508  gfermi = 1.16637e-5
509  ccabib = 0.975
510  gv = 1.0
511  ga =-1.0
512 
513 
514 
515 * ZW 13.04.89 HERE WAS AN ERROR
516  scabib = sqrt(1.-ccabib**2)
517  pi =4.*atan(1.)
518  gamel = gfermi**2*amtau**5/(192*pi**3)
519 *
520  CALL dexay(-1,pol1)
521 *
522  RETURN
523  END
524  FUNCTION dcdmas(IDENT)
525  COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
526  * ,ampiz,ampi,amro,gamro,ama1,gama1
527  * ,amk,amkz,amkst,gamkst
528 *
529  REAL*4 AMTAU,AMNUTA,AMEL,AMNUE,AMMU,AMNUMU
530  * ,ampiz,ampi,amro,gamro,ama1,gama1
531  * ,amk,amkz,amkst,gamkst
532  IF (ident.EQ. 1) THEN
533  apkmas=ampi
534  ELSEIF (ident.EQ.-1) THEN
535  apkmas=ampi
536  ELSEIF (ident.EQ. 2) THEN
537  apkmas=ampiz
538  ELSEIF (ident.EQ.-2) THEN
539  apkmas=ampiz
540  ELSEIF (ident.EQ. 3) THEN
541  apkmas=amk
542  ELSEIF (ident.EQ.-3) THEN
543  apkmas=amk
544  ELSEIF (ident.EQ. 4) THEN
545  apkmas=amkz
546  ELSEIF (ident.EQ.-4) THEN
547  apkmas=amkz
548  ELSEIF (ident.EQ. 8) THEN
549  apkmas=0.0001
550  ELSEIF (ident.EQ.-8) THEN
551  apkmas=0.0001
552  ELSEIF (ident.EQ. 9) THEN
553  apkmas=0.5488
554  ELSEIF (ident.EQ.-9) THEN
555  apkmas=0.5488
556  ELSE
557  print *, 'STOP IN APKMAS, WRONG IDENT=',ident
558  stop
559  ENDIF
560  dcdmas=apkmas
561  END
562  FUNCTION lunpik(ID,ISGN)
563  COMMON / taukle / bra1,brk0,brk0b,brks
564  REAL*4 BRA1,BRK0,BRK0B,BRKS
565  REAL*4 XIO(1)
566  ident=id*isgn
567 #if defined (ALEPH)
568  IF (ident.EQ. 1) THEN
569  ipkdef= 211
570  ELSEIF (ident.EQ.-1) THEN
571  ipkdef=-211
572  ELSEIF (ident.EQ. 2) THEN
573  ipkdef= 111
574  ELSEIF (ident.EQ.-2) THEN
575  ipkdef= 111
576  ELSEIF (ident.EQ. 3) THEN
577  ipkdef= 321
578  ELSEIF (ident.EQ.-3) THEN
579  ipkdef=-321
580 #else
581  IF (ident.EQ. 1) THEN
582  ipkdef=-211
583  ELSEIF (ident.EQ.-1) THEN
584  ipkdef= 211
585  ELSEIF (ident.EQ. 2) THEN
586  ipkdef=111
587  ELSEIF (ident.EQ.-2) THEN
588  ipkdef=111
589  ELSEIF (ident.EQ. 3) THEN
590  ipkdef=-321
591  ELSEIF (ident.EQ.-3) THEN
592  ipkdef= 321
593 #endif
594  ELSEIF (ident.EQ. 4) THEN
595 *
596 * K0 --> K0_LONG (IS 130) / K0_SHORT (IS 310) = 1/1
597  CALL ranmar(xio,1)
598  IF (xio(1).GT.brk0) THEN
599  ipkdef= 130
600  ELSE
601  ipkdef= 310
602  ENDIF
603  ELSEIF (ident.EQ.-4) THEN
604 *
605 * K0B--> K0_LONG (IS 130) / K0_SHORT (IS 310) = 1/1
606  CALL ranmar(xio,1)
607  IF (xio(1).GT.brk0b) THEN
608  ipkdef= 130
609  ELSE
610  ipkdef= 310
611  ENDIF
612  ELSEIF (ident.EQ. 8) THEN
613  ipkdef= 22
614  ELSEIF (ident.EQ.-8) THEN
615  ipkdef= 22
616  ELSEIF (ident.EQ. 9) THEN
617  ipkdef= 221
618  ELSEIF (ident.EQ.-9) THEN
619  ipkdef= 221
620  ELSE
621  print *, 'STOP IN IPKDEF, WRONG IDENT=',ident
622  stop
623  ENDIF
624  lunpik=ipkdef
625  END
626 
627 
628 #if defined (CLEO)
629 
630  SUBROUTINE taurdf(KTO)
631 C THIS ROUTINE CAN BE CALLED BEFORE ANY TAU+ OR TAU- EVENT IS GENERATED
632 C IT CAN BE USED TO GENERATE TAU+ AND TAU- SAMPLES OF DIFFERENT
633 C CONTENTS
634  COMMON / taukle / bra1,brk0,brk0b,brks
635  REAL*4 BRA1,BRK0,BRK0B,BRKS
636  COMMON / taubra / gamprt(30),jlist(30),nchan
637  IF (kto.EQ.1) THEN
638 C ==================
639 C AJWMOD: Set the BRs for (A1+ -> rho+ pi0) and (K*+ -> K0 pi+)
640  bra1 = pkorb(4,1)
641  brks = pkorb(4,3)
642  brk0 = pkorb(4,5)
643  brk0b = pkorb(4,6)
644  ELSE
645 C ====
646 C AJWMOD: Set the BRs for (A1+ -> rho+ pi0) and (K*+ -> K0 pi+)
647  bra1 = pkorb(4,2)
648  brks = pkorb(4,4)
649  brk0 = pkorb(4,5)
650  brk0b = pkorb(4,6)
651  ENDIF
652 C =====
653  END
654 #else
655 
656  SUBROUTINE taurdf(KTO)
657 * THIS ROUTINE CAN BE CALLED BEFORE ANY TAU+ OR TAU- EVENT IS GENERATED
658 * IT CAN BE USED TO GENERATE TAU+ AND TAU- SAMPLES OF DIFFERENT
659 * CONTENTS
660  COMMON / taukle / bra1,brk0,brk0b,brks
661  REAL*4 BRA1,BRK0,BRK0B,BRKS
662  COMMON / taubra / gamprt(30),jlist(30),nchan
663  IF (kto.EQ.1) THEN
664 * ==================
665 * LIST OF BRANCHING RATIOS
666  nchan = 19
667  DO 1 i = 1,30
668  IF (i.LE.nchan) THEN
669  jlist(i) = i
670  IF(i.EQ. 1) gamprt(i) = .0000
671  IF(i.EQ. 2) gamprt(i) = .0000
672  IF(i.EQ. 3) gamprt(i) = .0000
673  IF(i.EQ. 4) gamprt(i) = .0000
674  IF(i.EQ. 5) gamprt(i) = .0000
675  IF(i.EQ. 6) gamprt(i) = .0000
676  IF(i.EQ. 7) gamprt(i) = .0000
677  IF(i.EQ. 8) gamprt(i) = 1.0000
678  IF(i.EQ. 9) gamprt(i) = 1.0000
679  IF(i.EQ.10) gamprt(i) = 1.0000
680  IF(i.EQ.11) gamprt(i) = 1.0000
681  IF(i.EQ.12) gamprt(i) = 1.0000
682  IF(i.EQ.13) gamprt(i) = 1.0000
683  IF(i.EQ.14) gamprt(i) = 1.0000
684  IF(i.EQ.15) gamprt(i) = 1.0000
685  IF(i.EQ.16) gamprt(i) = 1.0000
686  IF(i.EQ.17) gamprt(i) = 1.0000
687  IF(i.EQ.18) gamprt(i) = 1.0000
688  IF(i.EQ.19) gamprt(i) = 1.0000
689  ELSE
690  jlist(i) = 0
691  gamprt(i) = 0.
692  ENDIF
693  1 CONTINUE
694 * --- COEFFICIENTS TO FIX RATIO OF:
695 * --- A1 3CHARGED/ A1 1CHARGED 2 NEUTRALS MATRIX ELEMENTS (MASLESS LIM.)
696 * --- PROBABILITY OF K0 TO BE KS
697 * --- PROBABILITY OF K0B TO BE KS
698 * --- RATIO OF COEFFICIENTS FOR K*--> K0 PI-
699 * --- ALL COEFFICENTS SHOULD BE IN THE RANGE (0.0,1.0)
700 * --- THEY MEANING IS PROBABILITY OF THE FIRST CHOICE ONLY IF ONE
701 * --- NEGLECTS MASS-PHASE SPACE EFFECTS
702  bra1=0.5
703  brk0=0.5
704  brk0b=0.5
705  brks=0.6667
706  ELSE
707 * ====
708 * LIST OF BRANCHING RATIOS
709  nchan = 19
710  DO 2 i = 1,30
711  IF (i.LE.nchan) THEN
712  jlist(i) = i
713  IF(i.EQ. 1) gamprt(i) = .0000
714  IF(i.EQ. 2) gamprt(i) = .0000
715  IF(i.EQ. 3) gamprt(i) = .0000
716  IF(i.EQ. 4) gamprt(i) = .0000
717  IF(i.EQ. 5) gamprt(i) = .0000
718  IF(i.EQ. 6) gamprt(i) = .0000
719  IF(i.EQ. 7) gamprt(i) = .0000
720  IF(i.EQ. 8) gamprt(i) = 1.0000
721  IF(i.EQ. 9) gamprt(i) = 1.0000
722  IF(i.EQ.10) gamprt(i) = 1.0000
723  IF(i.EQ.11) gamprt(i) = 1.0000
724  IF(i.EQ.12) gamprt(i) = 1.0000
725  IF(i.EQ.13) gamprt(i) = 1.0000
726  IF(i.EQ.14) gamprt(i) = 1.0000
727  IF(i.EQ.15) gamprt(i) = 1.0000
728  IF(i.EQ.16) gamprt(i) = 1.0000
729  IF(i.EQ.17) gamprt(i) = 1.0000
730  IF(i.EQ.18) gamprt(i) = 1.0000
731  IF(i.EQ.19) gamprt(i) = 1.0000
732  ELSE
733  jlist(i) = 0
734  gamprt(i) = 0.
735  ENDIF
736  2 CONTINUE
737 * --- COEFFICIENTS TO FIX RATIO OF:
738 * --- A1 3CHARGED/ A1 1CHARGED 2 NEUTRALS MATRIX ELEMENTS (MASLESS LIM.)
739 * --- PROBABILITY OF K0 TO BE KS
740 * --- PROBABILITY OF K0B TO BE KS
741 * --- RATIO OF COEFFICIENTS FOR K*--> K0 PI-
742 * --- ALL COEFFICENTS SHOULD BE IN THE RANGE (0.0,1.0)
743 * --- THEY MEANING IS PROBABILITY OF THE FIRST CHOICE ONLY IF ONE
744 * --- NEGLECTS MASS-PHASE SPACE EFFECTS
745  bra1=0.5
746  brk0=0.5
747  brk0b=0.5
748  brks=0.6667
749  ENDIF
750 * =====
751  END
752 #endif
753 
754  SUBROUTINE iniphx(XK00)
755 * ----------------------------------------------------------------------
756 * INITIALISATION OF PARAMETERS
757 * USED IN QED and/or GSW ROUTINES
758 * ----------------------------------------------------------------------
759  COMMON / qedprm /alfinv,alfpi,xk0
760  REAL*8 ALFINV,ALFPI,XK0
761  REAL*8 PI8,XK00
762 *
763  pi8 = 4.d0*datan(1.d0)
764  alfinv = 137.03604d0
765  alfpi = 1d0/(alfinv*pi8)
766  xk0=xk00
767  END
768 
769  SUBROUTINE inimas
770 C ----------------------------------------------------------------------
771 C INITIALISATION OF MASSES
772 C
773 C called by : KORALZ
774 C ----------------------------------------------------------------------
775  COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
776  * ,ampiz,ampi,amro,gamro,ama1,gama1
777  * ,amk,amkz,amkst,gamkst
778 *
779  REAL*4 AMTAU,AMNUTA,AMEL,AMNUE,AMMU,AMNUMU
780  * ,ampiz,ampi,amro,gamro,ama1,gama1
781  * ,amk,amkz,amkst,gamkst
782 C
783 C IN-COMING / OUT-GOING FERMION MASSES
784  amtau = 1.7842
785 C --- tau mass must be the same as in the host program, what-so-ever
786  amtau = 1.777
787  amnuta = 0.010
788  amel = 0.0005111
789  amnue = 0.0
790  ammu = 0.105659
791  amnumu = 0.0
792 *
793 * MASSES USED IN TAU DECAYS
794 #if defined (CePeCe)
795  ampiz = 0.134964
796  ampi = 0.139568
797  amro = 0.773
798  gamro = 0.145
799 *C GAMRO = 0.666
800  ama1 = 1.251
801  gama1 = 0.599
802  amk = 0.493667
803  amkz = 0.49772
804  amkst = 0.8921
805  gamkst = 0.0513
806 #elif defined (CLEO)
807  ampiz = 0.134964
808  ampi = 0.139568
809  amro = 0.773
810  gamro = 0.145
811 *C GAMRO = 0.666
812  ama1 = 1.251
813  gama1 = 0.599
814  amk = 0.493667
815  amkz = 0.49772
816  amkst = 0.8921
817  gamkst = 0.0513
818 C
819 C
820 C IN-COMING / OUT-GOING FERMION MASSES
821 !! AMNUTA = PKORB(1,2)
822 !! AMNUE = PKORB(1,4)
823 !! AMNUMU = PKORB(1,6)
824 C
825 C MASSES USED IN TAU DECAYS Cleo settings
826 !! AMPIZ = PKORB(1,7)
827 !! AMPI = PKORB(1,8)
828 !! AMRO = PKORB(1,9)
829 !! GAMRO = PKORB(2,9)
830  ama1 = 1.275 !! PKORB(1,10)
831  gama1 = 0.615 !! PKORB(2,10)
832 !! AMK = PKORB(1,11)
833 !! AMKZ = PKORB(1,12)
834 !! AMKST = PKORB(1,13)
835 !! GAMKST = PKORB(2,13)
836 C
837 #elif defined (ALEPH)
838  ampiz = 0.134964
839  ampi = 0.139568
840  amro = 0.7714
841  gamro = 0.153
842 cam AMRO = 0.773
843 cam GAMRO = 0.145
844  ama1 = 1.251! PMAS(LUCOMP(ia1),1) ! AMA1 = 1.251
845  gama1 = 0.599! PMAS(LUCOMP(ia1),2) ! GAMA1 = 0.599
846  print *,'INIMAS a1 mass= ',ama1,gama1
847  amk = 0.493667
848  amkz = 0.49772
849  amkst = 0.8921
850  gamkst = 0.0513
851 #else
852 #endif
853 
854  RETURN
855  END
856  subroutine bostdq(idir,vv,pp,q)
857 * *******************************
858 c Boost along arbitrary vector v (see eg. J.D. Jacson, Classical
859 c Electrodynamics).
860 c Four-vector pp is boosted from an actual frame to the rest frame
861 c of the four-vector v (for idir=1) or back (for idir=-1).
862 c q is a resulting four-vector.
863 c Note: v must be time-like, pp may be arbitrary.
864 c
865 c Written by: Wieslaw Placzek date: 22.07.1994
866 c Last update: 3/29/95 by: M.S.
867 c
868  implicit DOUBLE PRECISION (a-h,o-z)
869  parameter(nout=6)
870  DOUBLE PRECISION v(4),p(4),q(4),pp(4),vv(4)
871  save
872 !
873  do 1 i=1,4
874  v(i)=vv(i)
875  1 p(i)=pp(i)
876  amv=(v(4)**2-v(1)**2-v(2)**2-v(3)**2)
877  if (amv.le.0d0) then
878  write(6,*) 'bosstv: warning amv**2=',amv
879  endif
880  amv=sqrt(abs(amv))
881  if (idir.eq.-1) then
882  q(4)=( p(1)*v(1)+p(2)*v(2)+p(3)*v(3)+p(4)*v(4))/amv
883  wsp =(q(4)+p(4))/(v(4)+amv)
884  elseif (idir.eq.1) then
885  q(4)=(-p(1)*v(1)-p(2)*v(2)-p(3)*v(3)+p(4)*v(4))/amv
886  wsp =-(q(4)+p(4))/(v(4)+amv)
887  else
888  write(nout,*)' >>> boostv: wrong value of idir = ',idir
889  endif
890  q(1)=p(1)+wsp*v(1)
891  q(2)=p(2)+wsp*v(2)
892  q(3)=p(3)+wsp*v(3)
893  end
894 
895 
896 
897 
898 
899 
900 
901 
902 
903 #if defined (ALEPH)
904  FUNCTION dilogy(X)
905 C *****************
906  IMPLICIT REAL*8(a-h,o-z)
907 CERN C304 VERSION 29/07/71 DILOG 59 C
908  z=-1.64493406684822
909  IF(x .LT.-1.0) GO TO 1
910  IF(x .LE. 0.5) GO TO 2
911  IF(x .EQ. 1.0) GO TO 3
912  IF(x .LE. 2.0) GO TO 4
913  z=3.2898681336964
914  1 t=1.0/x
915  s=-0.5
916  z=z-0.5* log(abs(x))**2
917  GO TO 5
918  2 t=x
919  s=0.5
920  z=0.
921  GO TO 5
922  3 dilogy=1.64493406684822
923  RETURN
924  4 t=1.0-x
925  s=-0.5
926  z=1.64493406684822 - log(x)* log(abs(t))
927  5 y=2.66666666666666 *t+0.66666666666666
928  b= 0.00000 00000 00001
929  a=y*b +0.00000 00000 00004
930  b=y*a-b+0.00000 00000 00011
931  a=y*b-a+0.00000 00000 00037
932  b=y*a-b+0.00000 00000 00121
933  a=y*b-a+0.00000 00000 00398
934  b=y*a-b+0.00000 00000 01312
935  a=y*b-a+0.00000 00000 04342
936  b=y*a-b+0.00000 00000 14437
937  a=y*b-a+0.00000 00000 48274
938  b=y*a-b+0.00000 00001 62421
939  a=y*b-a+0.00000 00005 50291
940  b=y*a-b+0.00000 00018 79117
941  a=y*b-a+0.00000 00064 74338
942  b=y*a-b+0.00000 00225 36705
943  a=y*b-a+0.00000 00793 87055
944  b=y*a-b+0.00000 02835 75385
945  a=y*b-a+0.00000 10299 04264
946  b=y*a-b+0.00000 38163 29463
947  a=y*b-a+0.00001 44963 00557
948  b=y*a-b+0.00005 68178 22718
949  a=y*b-a+0.00023 20021 96094
950  b=y*a-b+0.00100 16274 96164
951  a=y*b-a+0.00468 63619 59447
952  b=y*a-b+0.02487 93229 24228
953  a=y*b-a+0.16607 30329 27855
954  a=y*a-b+1.93506 43008 6996
955  dilogy=s*t*(a-b)+z
956  RETURN
957 C=======================================================================
958 C===================END OF CPC PART ====================================
959 C=======================================================================
960  END
961 #endif