C++InterfacetoTauola
demo-standalone/taumain.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  PROGRAM taudem
35 c **************
36 c note that the routines are not like in cpc deck this is historical !!
37 c=======================================================================
38 c====================== dectes : test of tau decay library===========
39 c====================== ktory = 1 : INTERFACE of koral-z TYPE ==========
40 c====================== ktory = 2 : INTERFACE of koral-b TYPE =========
41 c=======================================================================
42 c COMMON /pawc/ blan(10000)
43  COMMON / / blan(10000)
44  CHARACTER*7 dname
45  COMMON / inout / inut,iout
46  dname='KKPI'
47 ! CALL GLIMIT(20000)
48 ! CALL GOUTPU(16)
49  inut=5
50  iout=6
51  OPEN(iout,file="./tauola.output")
52  OPEN(inut,file="./dane.dat")
53  ktory=1
54  CALL dectes(ktory)
55  ktory=2
56  CALL dectes(ktory)
57  END
58  SUBROUTINE dectes(KTORY)
59 c ************************
60  REAL pol(4)
61  DOUBLE PRECISION hh(4)
62 c switches for tauola;
63  COMMON / jaki / jak1,jak2,jakp,jakm,ktom
64  COMMON / idfc / idff
65 c i/o units numbers
66  COMMON / inout / inut,iout
67 c lund TYPE identifier for a1
68  COMMON / idpart / ia1
69 c /ptau/ is used in routine tralo4
70  COMMON /ptau/ ptau
71  COMMON / taurad / xk0dec,itdkrc
72  REAL*8 xk0dec
73  COMMON /testa1/ keya1
74 c special switch for tests of dgamma/dq**2 in a1 decay
75 c keya1=1 constant width of a1 and rho
76 c keya1=2 free choice of rho propagator(defined in function fpik)
77 c and free choice of a1 mass and width. function g(Q**2)
78 c(see formula 3.48 in comp. phys. comm. 64 (1991) 275)
79 c hard coded both in monte carlo and in testing distribution.
80 c keya1=3 function g(Q**2) hardcoded in the Monte Carlo
81 c(it is timy to calculate!), but appropriately adjusted in
82 c testing distribution.
83 c-----------------------------------------------------------------------
84 c initialization
85 c-----------------------------------------------------------------------
86 c======================================
87  ninp=inut
88  nout=iout
89  3000 FORMAT(a80)
90  3001 FORMAT(8i2)
91  3002 FORMAT(i10)
92  3003 FORMAT(f10.0)
93  IF (ktory.EQ.1) THEN
94  READ( ninp,3000) testit
95  WRITE(nout,3000) testit
96  READ( ninp,3001) kat1,kat2,kat3,kat4,kat5,kat6
97  READ( ninp,3002) nevt,jak1,jak2,itdkrc
98  READ( ninp,3003) ptau,xk0dec
99  ENDIF
100 c======================================
101 c control output
102  WRITE(nout,'(6A6/6I6)')
103  $ 'KAT1','KAT2','KAT3','KAT4','KAT5','KAT6',
104  $ kat1 , kat2 , kat3 , kat4 , kat5 , kat6
105  WRITE(nout,'(4A12/4I12)')
106  $ 'NEVT','JAK1','JAK2','ITDKRC',
107  $ nevt, jak1 , jak2 , itdkrc
108  WRITE(nout,'(2A12/2F12.6)')
109  $ 'PTAU','XK0DEC',
110  $ ptau , xk0dec
111 c======================================
112  jak=0
113 c jak1=5
114 c jak2=5
115 c lund identifier(for tau+) -15
116  IF (ktory.EQ.1) THEN
117  idff=-15
118  ELSE
119  idff= 15
120  ENDIF
121 c kto=1 denotes tau defined by idff(i.e. tau+)
122 c kto=2 denotes the opposite(i.e. tau-)
123  kto=2
124  IF (kto.NE.2) THEN
125  print *, 'for the sake of these tests KTO has to be 2'
126  print *, 'to change tau- to tau+ change IDFF from -15 to 15'
127  stop
128  ENDIF
129 c tau polarization in its restframe;
130  pol(1)=0.
131  pol(2)=0.
132  pol(3)=.9
133 c tau momentum in gev;
134 c ptau=cmsene/2.d0
135 c number of events to be generated;
136  nevtes=10
137  nevtes=nevt
138  print *, 'NEVTES= ',nevtes
139  WRITE(iout,7011) keya1
140 c
141  IF (ktory.EQ.1) THEN
142  WRITE(iout,7001) jak,idff,pol(3),ptau
143  ELSE
144  WRITE(iout,7004) jak,idff,pol(3),ptau
145  ENDIF
146 c initialisation of tau decay package tauola
147 c ******************************************
148  CALL inimas
149  CALL initdk
150 
151 
152  CALL iniphy(0.1d0)
153  IF (ktory.EQ.1) THEN
154  CALL dexay(-1,pol)
155  ELSE
156  CALL dekay(-1,hh)
157  ENDIF
158 c-----------------------------------------------------------------------
159 c generation
160 c-----------------------------------------------------------------------
161  nev=0
162  DO 300 iev=1,nevtes
163  nev=nev+1
164 c reslu initialise the lund record
165  CALL taufil
166 c decay....
167  IF (ktory.EQ.1) THEN
168  CALL dexay(kto,pol)
169  ELSE
170  CALL dekay(kto,hh)
171  CALL dekay(kto+10,hh)
172  ENDIF
173  CALL luhepc(2)
174  IF(iev.LE.44) THEN
175  WRITE(iout,7002) iev
176  IF (ktory.NE.1) THEN
177  WRITE(iout,7003) hh
178  ENDIF
179 c CALL lulist(11)
180  CALL lulist(2)
181  ENDIF
182  ipri=mod(nev,1000)
183  IF(ipri.EQ.1) print *, ' event no: ',nev,' NEVTES: ',nevtes
184  300 CONTINUE
185  301 CONTINUE
186 c-----------------------------------------------------------------------
187 c postgeneration
188 c-----------------------------------------------------------------------
189  IF (ktory.EQ.1) THEN
190  CALL dexay(100,pol)
191  ELSE
192  CALL dekay(100,hh)
193  ENDIF
194  RETURN
195  7001 FORMAT(//4(/1x,15(5h=====))
196  $ /,' ', 19x,' TEST OF RAD. CORR IN ELECTRON DECAY ',9x,1h ,
197  $ /,' ', 19x,' TESTS OF TAU DECAY ROUTINES ',9x,1h ,
198  $ /,' ', 19x,' INTERFACE OF THE KORAL-Z TYPE ',9x,1h ,
199  $ 2(/,1x,15(5h=====)),
200  $ /,5x ,'JAK =',i7 ,' KEY DEFINING DECAY TYPE ',9x,1h ,
201  $ /,5x ,'IDFF =',i7 ,' LUND IDENTIFIER FOR FIRST TAU ',9x,1h ,
202  $ /,5x ,'POL(3)=',f7.2,' THIRD COMPONENT OF TAU POLARIZ. ',9x,1h ,
203  $ /,5x ,'PTAU =',f7.2,' THIRD COMPONENT OF TAU MOM. GEV ',9x,1h ,
204  $ 2(/,1x,15(5h=====))/)
205  7002 FORMAT(///1x, '===== EVENT NO.',i4,1x,5h=====)
206  7003 FORMAT(5x,'POLARIMETRIC VECTOR: ',
207  $ 7x,'HH(1)',7x,'HH(2)',7x,'HH(3)',7x,'HH(4)',
208  $ /, 5x,' ', 4(1x,f11.8) )
209  7004 FORMAT(//4(/1x,15(5h=====))
210  $ /,' ', 19x,' TEST OF RAD. CORR IN ELECTRON DECAY ',9x,1h ,
211  $ /,' ', 19x,' TESTS OF TAU DECAY ROUTINES ',9x,1h ,
212  $ /,' ', 19x,' INTERFACE OF THE KORAL-B TYPE ',9x,1h ,
213  $ 2(/,1x,15(5h=====)),
214  $ /,5x ,'JAK =',i7 ,' KEY DEFINING DECAY TYPE ',9x,1h ,
215  $ /,5x ,'IDFF =',i7 ,' LUND IDENTIFIER FOR FIRST TAU ',9x,1h ,
216  $ /,5x ,'POL(3)=',f7.2,' THIRD COMPONENT OF TAU POLARIZ. ',9x,1h ,
217  $ /,5x ,'PTAU =',f7.2,' THIRD COMPONENT OF TAU MOM. GEV ',9x,1h ,
218  $ 2(/,1x,15(5h=====))/)
219  7011 FORMAT(///1x, '===== TYPE OF CURRENT',i4,1x,5h=====)
220  END
221  SUBROUTINE choice(MNUM,RR,ICHAN,PROB1,PROB2,PROB3,
222  $ amrx,gamrx,amra,gamra,amrb,gamrb)
223  COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
224  * ,ampiz,ampi,amro,gamro,ama1,gama1
225  * ,amk,amkz,amkst,gamkst
226 c
227  REAL*4 amtau,amnuta,amel,amnue,ammu,amnumu
228  * ,ampiz,ampi,amro,gamro,ama1,gama1
229  * ,amk,amkz,amkst,gamkst
230 c
231  amrop=1.1
232  gamrop=0.36
233  amom=.782
234  gamom=0.0084
235 c xxxxa correspond to s2 channel !
236  IF(mnum.EQ.0) THEN
237  prob1=0.5
238  prob2=0.5
239  amrx =ama1
240  gamrx=gama1
241  amra =amro
242  gamra=gamro
243  amrb =amro
244  gamrb=gamro
245  ELSEIF(mnum.EQ.1) THEN
246  prob1=0.5
247  prob2=0.5
248  amrx =1.57
249  gamrx=0.9
250  amrb =amkst
251  gamrb=gamkst
252  amra =amro
253  gamra=gamro
254  ELSEIF(mnum.EQ.2) THEN
255  prob1=0.5
256  prob2=0.5
257  amrx =1.57
258  gamrx=0.9
259  amrb =amkst
260  gamrb=gamkst
261  amra =amro
262  gamra=gamro
263  ELSEIF(mnum.EQ.3) THEN
264  prob1=0.5
265  prob2=0.5
266  amrx =1.27
267  gamrx=0.3
268  amra =amkst
269  gamra=gamkst
270  amrb =amkst
271  gamrb=gamkst
272  ELSEIF(mnum.EQ.4) THEN
273  prob1=0.5
274  prob2=0.5
275  amrx =1.27
276  gamrx=0.3
277  amra =amkst
278  gamra=gamkst
279  amrb =amkst
280  gamrb=gamkst
281  ELSEIF(mnum.EQ.5) THEN
282  prob1=0.5
283  prob2=0.5
284  amrx =1.27
285  gamrx=0.3
286  amra =amkst
287  gamra=gamkst
288  amrb =amro
289  gamrb=gamro
290  ELSEIF(mnum.EQ.6) THEN
291  prob1=0.4
292  prob2=0.4
293  amrx =1.27
294  gamrx=0.3
295  amra =amro
296  gamra=gamro
297  amrb =amkst
298  gamrb=gamkst
299  ELSEIF(mnum.EQ.7) THEN
300  prob1=0.0
301  prob2=1.0
302  amrx =1.27
303  gamrx=0.9
304  amra =amro
305  gamra=gamro
306  amrb =amro
307  gamrb=gamro
308  ELSEIF(mnum.EQ.8) THEN
309  prob1=0.0
310  prob2=1.0
311  amrx =amrop
312  gamrx=gamrop
313  amrb =amom
314  gamrb=gamom
315  amra =amro
316  gamra=gamro
317  ELSEIF(mnum.EQ.101) THEN
318  prob1=.35
319  prob2=.35
320  amrx =1.2
321  gamrx=.46
322  amrb =amom
323  gamrb=gamom
324  amra =amom
325  gamra=gamom
326  ELSEIF(mnum.EQ.102) THEN
327  prob1=0.0
328  prob2=0.0
329  amrx =1.4
330  gamrx=.6
331  amrb =amom
332  gamrb=gamom
333  amra =amom
334  gamra=gamom
335  ELSE
336  prob1=0.0
337  prob2=0.0
338  amrx =ama1
339  gamrx=gama1
340  amra =amro
341  gamra=gamro
342  amrb =amro
343  gamrb=gamro
344  ENDIF
345 c
346  IF (rr.LE.prob1) THEN
347  ichan=1
348  ELSEIF(rr.LE.(prob1+prob2)) THEN
349  ichan=2
350  ax =amra
351  gx =gamra
352  amra =amrb
353  gamra=gamrb
354  amrb =ax
355  gamrb=gx
356  px =prob1
357  prob1=prob2
358  prob2=px
359  ELSE
360  ichan=3
361  ENDIF
362 c
363  prob3=1.0-prob1-prob2
364  END
365  SUBROUTINE initdk
366 * ----------------------------------------------------------------------
367 * initialisation of tau decay parameters and routines
368 *
369 * called by : koralz
370 * ----------------------------------------------------------------------
371 
372  COMMON / decpar / gfermi,gv,ga,ccabib,scabib,gamel
373  REAL*4 gfermi,gv,ga,ccabib,scabib,gamel
374  COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
375  * ,ampiz,ampi,amro,gamro,ama1,gama1
376  * ,amk,amkz,amkst,gamkst
377 *
378  REAL*4 amtau,amnuta,amel,amnue,ammu,amnumu
379  * ,ampiz,ampi,amro,gamro,ama1,gama1
380  * ,amk,amkz,amkst,gamkst
381  COMMON / taubra / gamprt(30),jlist(30),nchan
382  COMMON / taukle / bra1,brk0,brk0b,brks
383  REAL*4 bra1,brk0,brk0b,brks
384  parameter(nmode=15,nm1=0,nm2=1,nm3=8,nm4=2,nm5=1,nm6=3)
385  COMMON / taudcd /idffin(9,nmode),mulpik(nmode)
386  & ,names
387  CHARACTER names(nmode)*31
388  CHARACTER oldnames(7)*31
389  CHARACTER*80 bxinit
390  parameter(
391  $ bxinit ='(1x,1h*,g17.8, 16x, a31,a4,a4, 1x,1h*)'
392  $ )
393  REAL*4 pi,pol1(4)
394 *
395 *
396 * list of branching ratios
397 cam normalised to e nu nutau channel
398 cam enu munu pinu rhonu a1nu knu k*nu pi
399 cam DATA jlist / 1, 2, 3, 4, 5, 6, 7,
400 *am DATA gamprt /1.000,0.9730,0.6054,1.2432,0.8432,0.0432,o.o811,0.616
401 *am
402 *am multipion decays
403 *
404 * conventions of particles names
405 * k-,p-,k+, k0,p-,kb, k-,p0,k0
406 * 3, 1,-3 , 4, 1,-4 , 3, 2, 4 ,
407 * p0,p0,k-, k-,p-,p+, p-,kb,p0
408 * 2, 2, 3 , 3, 1,-1 , 1,-4, 2 ,
409 * et,p-,p0 p-,p0,gm
410 * 9, 1, 2 , 1, 2, 8
411 *
412 c
413  dimension nopik(6,nmode),npik(nmode)
414 *am outgoing multiplicity and flavors of multi-pion /multi-k modes
415  DATA npik / 4, 4,
416  1 5, 5,
417  2 6, 6,
418  3 3, 3,
419  4 3, 3,
420  5 3, 3,
421  6 3, 3,
422  7 2 /
423  DATA nopik / -1,-1, 1, 2, 0, 0, 2, 2, 2,-1, 0, 0,
424  1 -1,-1, 1, 2, 2, 0, -1,-1,-1, 1, 1, 0,
425  2 -1,-1,-1, 1, 1, 2, -1,-1, 1, 2, 2, 2,
426  3 -3,-1, 3, 0, 0, 0, -4,-1, 4, 0, 0, 0,
427  4 -3, 2,-4, 0, 0, 0, 2, 2,-3, 0, 0, 0,
428  5 -3,-1, 1, 0, 0, 0, -1, 4, 2, 0, 0, 0,
429  6 9,-1, 2, 0, 0, 0, -1, 2, 8, 0, 0, 0,
430 c ajwmod fix sign bug, 2/22/99
431  7 -3,-4, 0, 0, 0, 0 /
432 * list of branching ratios
433  nchan = nmode + 7
434  DO 1 i = 1,30
435  IF (i.LE.nchan) THEN
436  jlist(i) = i
437  IF(i.EQ. 1) gamprt(i) =0.1800
438  IF(i.EQ. 2) gamprt(i) =0.1751
439  IF(i.EQ. 3) gamprt(i) =0.1110
440  IF(i.EQ. 4) gamprt(i) =0.2515
441  IF(i.EQ. 5) gamprt(i) =0.1790
442  IF(i.EQ. 6) gamprt(i) =0.0071
443  IF(i.EQ. 7) gamprt(i) =0.0134
444  IF(i.EQ. 8) gamprt(i) =0.0450
445  IF(i.EQ. 9) gamprt(i) =0.0100
446  IF(i.EQ.10) gamprt(i) =0.0009
447  IF(i.EQ.11) gamprt(i) =0.0004
448  IF(i.EQ.12) gamprt(i) =0.0003
449  IF(i.EQ.13) gamprt(i) =0.0005
450  IF(i.EQ.14) gamprt(i) =0.0015
451  IF(i.EQ.15) gamprt(i) =0.0015
452  IF(i.EQ.16) gamprt(i) =0.0015
453  IF(i.EQ.17) gamprt(i) =0.0005
454  IF(i.EQ.18) gamprt(i) =0.0050
455  IF(i.EQ.19) gamprt(i) =0.0055
456  IF(i.EQ.20) gamprt(i) =0.0017
457  IF(i.EQ.21) gamprt(i) =0.0013
458  IF(i.EQ.22) gamprt(i) =0.0010
459  IF(i.EQ. 1) oldnames(i)=' TAU- --> E- '
460  IF(i.EQ. 2) oldnames(i)=' TAU- --> MU- '
461  IF(i.EQ. 3) oldnames(i)=' TAU- --> PI- '
462  IF(i.EQ. 4) oldnames(i)=' TAU- --> PI-, PI0 '
463  IF(i.EQ. 5) oldnames(i)=' TAU- --> A1- (two subch) '
464  IF(i.EQ. 6) oldnames(i)=' TAU- --> K- '
465  IF(i.EQ. 7) oldnames(i)=' TAU- --> K*- (two subch) '
466  IF(i.EQ. 8) names(i-7)=' TAU- --> 2PI-, PI0, PI+ '
467  IF(i.EQ. 9) names(i-7)=' TAU- --> 3PI0, PI- '
468  IF(i.EQ.10) names(i-7)=' TAU- --> 2PI-, PI+, 2PI0 '
469  IF(i.EQ.11) names(i-7)=' TAU- --> 3PI-, 2PI+, '
470  IF(i.EQ.12) names(i-7)=' TAU- --> 3PI-, 2PI+, PI0 '
471  IF(i.EQ.13) names(i-7)=' TAU- --> 2PI-, PI+, 3PI0 '
472  IF(i.EQ.14) names(i-7)=' TAU- --> K-, PI-, K+ '
473  IF(i.EQ.15) names(i-7)=' TAU- --> K0, PI-, K0B '
474  IF(i.EQ.16) names(i-7)=' TAU- --> K-, K0, PI0 '
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 (ident.EQ. 1) THEN
568  ipkdef=-211
569  ELSEIF (ident.EQ.-1) THEN
570  ipkdef= 211
571  ELSEIF (ident.EQ. 2) THEN
572  ipkdef=111
573  ELSEIF (ident.EQ.-2) THEN
574  ipkdef=111
575  ELSEIF (ident.EQ. 3) THEN
576  ipkdef=-321
577  ELSEIF (ident.EQ.-3) THEN
578  ipkdef= 321
579  ELSEIF (ident.EQ. 4) THEN
580 *
581 * k0 --> k0_long(is 130) / k0_short(is 310) = 1/1
582  CALL ranmar(xio,1)
583  IF (xio(1).GT.brk0) THEN
584  ipkdef= 130
585  ELSE
586  ipkdef= 310
587  ENDIF
588  ELSEIF (ident.EQ.-4) THEN
589 *
590 * k0b--> k0_long(is 130) / k0_short(is 310) = 1/1
591  CALL ranmar(xio,1)
592  IF (xio(1).GT.brk0b) THEN
593  ipkdef= 130
594  ELSE
595  ipkdef= 310
596  ENDIF
597  ELSEIF (ident.EQ. 8) THEN
598  ipkdef= 22
599  ELSEIF (ident.EQ.-8) THEN
600  ipkdef= 22
601  ELSEIF (ident.EQ. 9) THEN
602  ipkdef= 221
603  ELSEIF (ident.EQ.-9) THEN
604  ipkdef= 221
605  ELSE
606  print *, 'STOP IN IPKDEF, WRONG IDENT=',ident
607  stop
608  ENDIF
609  lunpik=ipkdef
610  END
611 
612 
613 
614  SUBROUTINE taurdf(KTO)
615 c this routine can be called before any tau+ or tau- event is generated
616 c it can be used to generate tau+ and tau- samples of different
617 c contents
618  COMMON / taukle / bra1,brk0,brk0b,brks
619  REAL*4 bra1,brk0,brk0b,brks
620  COMMON / taubra / gamprt(30),jlist(30),nchan
621  IF (kto.EQ.1) THEN
622 c ==================
623 c ajwmod: set the brs for(a1+ -> rho+ pi0) and(k*+ -> k0 pi+)
624  bra1 = pkorb(4,1)
625  brks = pkorb(4,3)
626  brk0 = pkorb(4,5)
627  brk0b = pkorb(4,6)
628  ELSE
629 c ====
630 c ajwmod: set the brs for(a1+ -> rho+ pi0) and(k*+ -> k0 pi+)
631  bra1 = pkorb(4,2)
632  brks = pkorb(4,4)
633  brk0 = pkorb(4,5)
634  brk0b = pkorb(4,6)
635  ENDIF
636 c =====
637  END
638 
639  SUBROUTINE iniphy(XK00)
640 * ----------------------------------------------------------------------
641 * initialisation of parameters
642 * used in qed and/or gsw routines
643 * ----------------------------------------------------------------------
644  COMMON / qedprm /alfinv,alfpi,xk0
645  REAL*8 alfinv,alfpi,xk0
646  REAL*8 pi8,xk00
647 *
648  pi8 = 4.d0*datan(1.d0)
649  alfinv = 137.03604d0
650  alfpi = 1d0/(alfinv*pi8)
651  xk0=xk00
652  END
653 
654  SUBROUTINE inimas
655 c ----------------------------------------------------------------------
656 c initialisation of masses
657 c
658 c called by : koralz
659 c ----------------------------------------------------------------------
660  COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
661  * ,ampiz,ampi,amro,gamro,ama1,gama1
662  * ,amk,amkz,amkst,gamkst
663 *
664  REAL*4 amtau,amnuta,amel,amnue,ammu,amnumu
665  * ,ampiz,ampi,amro,gamro,ama1,gama1
666  * ,amk,amkz,amkst,gamkst
667 c
668 c in-coming / out-going fermion masses
669  amtau = 1.7842
670 c --- let us update tau mass ...
671  amtau = 1.777
672  amnuta = 0.010
673  amel = 0.0005111
674  amnue = 0.0
675  ammu = 0.105659
676  amnumu = 0.0
677 *
678 * masses used in tau decays
679  ampiz = 0.134964
680  ampi = 0.139568
681  amro = 0.773
682  gamro = 0.145
683 *c gamro = 0.666
684  ama1 = 1.251
685  gama1 = 0.599
686  amk = 0.493667
687  amkz = 0.49772
688  amkst = 0.8921
689  gamkst = 0.0513
690 c
691 c
692 c in-coming / out-going fermion masses
693 !! AMNUTA = PKORB(1,2)
694 !! AMNUE = PKORB(1,4)
695 !! AMNUMU = PKORB(1,6)
696 c
697 c masses used in tau decays cleo settings
698 !! AMPIZ = PKORB(1,7)
699 !! AMPI = PKORB(1,8)
700 !! AMRO = PKORB(1,9)
701 !! GAMRO = PKORB(2,9)
702  ama1 = 1.275 !! PKORB(1,10)
703  gama1 = 0.615 !! PKORB(2,10)
704 !! AMK = PKORB(1,11)
705 !! AMKZ = PKORB(1,12)
706 !! AMKST = PKORB(1,13)
707 !! GAMKST = PKORB(2,13)
708 c
709 
710  RETURN
711  END
712  SUBROUTINE taufil
713 c *****************
714 c subsitute of tau production generator
715 c
716  COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
717  * ,ampiz,ampi,amro,gamro,ama1,gama1
718  * ,amk,amkz,amkst,gamkst
719 c
720  REAL*4 amtau,amnuta,amel,amnue,ammu,amnumu
721  * ,ampiz,ampi,amro,gamro,ama1,gama1
722  * ,amk,amkz,amkst,gamkst
723  COMMON / idfc / idff
724 c positions of taus in the lund common block
725 c it will be used by tauola output routines.
726  COMMON /taupos / npa,npb
727  dimension xpb1(4),xpb2(4),aqf1(4),aqf2(4)
728 c
729 c --- defining dummy events momenta
730  DO 4 k=1,3
731  xpb1(k)=0.0
732  xpb2(k)=0.0
733  aqf1(k)=0.0
734  aqf2(k)=0.0
735  4 CONTINUE
736  aqf1(4)=amtau
737  aqf2(4)=amtau
738 c --- tau momenta
739  CALL tralo4(1,aqf1,aqf1,am)
740  CALL tralo4(2,aqf2,aqf2,am)
741 c --- beams momenta and identifiers
742  kfb1= 11*idff/iabs(idff)
743  kfb2=-11*idff/iabs(idff)
744  xpb1(4)= aqf1(4)
745  xpb1(3)= aqf1(4)
746  IF(aqf1(3).NE.0.0)
747  $ xpb1(3)= aqf1(4)*aqf1(3)/abs(aqf1(3))
748  xpb2(4)= aqf2(4)
749  xpb2(3)=-aqf2(4)
750  IF(aqf2(3).NE.0.0)
751  $ xpb2(3)= aqf2(4)*aqf2(3)/abs(aqf2(3))
752 c --- position of first and second tau in lund common
753  npa=3
754  npb=4
755 c --- fill to lund COMMON
756  CALL filhep( 1,3, kfb1,0,0,0,0,xpb1, amel,.true.)
757  CALL filhep( 2,3, kfb2,0,0,0,0,xpb2, amel,.true.)
758  CALL filhep(npa,1, idff,1,2,0,0,aqf1,amtau,.true.)
759  CALL filhep(npb,1,-idff,1,2,0,0,aqf2,amtau,.true.)
760  END
761  SUBROUTINE tralo4(KTO,P,Q,AM)
762 c **************************
763 c subsitute of tralo4
764  REAL p(4),q(4)
765 c
766  COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
767  * ,ampiz,ampi,amro,gamro,ama1,gama1
768  * ,amk,amkz,amkst,gamkst
769 c
770  REAL*4 amtau,amnuta,amel,amnue,ammu,amnumu
771  * ,ampiz,ampi,amro,gamro,ama1,gama1
772  * ,amk,amkz,amkst,gamkst
773  COMMON /ptau/ ptau
774  am=amas4(p)
775  etau=sqrt(ptau**2+amtau**2)
776  exe=(etau+ptau)/amtau
777  IF(kto.EQ.2) exe=(etau-ptau)/amtau
778  CALL bostr3(exe,p,q)
779 c ======================================================================
780 c END of the test job
781 c ======================================================================
782  END
783  SUBROUTINE filhep(N,IST,ID,JMO1,JMO2,JDA1,JDA2,P4,PINV,PHFLAG)
784 c ----------------------------------------------------------------------
785 c this subroutine fills one entry into the hepevt common
786 c and updates the information for affected mother entries
787 c
788 c written by martin w. gruenewald(91/01/28)
789 c
790 c called by : ztohep,btohep,dwluxy
791 c ----------------------------------------------------------------------
792 c
793 c this is the hepevt class in old style. no d_h_ class pre-name
794  INTEGER nmxhep
795  parameter(nmxhep=10000)
796  REAL*8 phep, vhep ! to be real*4/ *8 depending on host
797  INTEGER nevhep,nhep,isthep,idhep,jmohep,
798  $ jdahep
799  COMMON /hepevt/
800  $ nevhep, ! serial number
801  $ nhep, ! number of particles
802  $ isthep(nmxhep), ! status code
803  $ idhep(nmxhep), ! particle ident KF
804  $ jmohep(2,nmxhep), ! parent particles
805  $ jdahep(2,nmxhep), ! childreen particles
806  $ phep(5,nmxhep), ! four-momentum, mass [GeV]
807  $ vhep(4,nmxhep) ! vertex [mm]
808 * ----------------------------------------------------------------------
809  LOGICAL qedrad
810  COMMON /phoqed/
811  $ qedrad(nmxhep) ! Photos flag
812 * ----------------------------------------------------------------------
813  SAVE hepevt,phoqed
814 c parameter(nmxhep=2000)
815 c common/hepevt/nevhep,nhep,isthep(nmxhep),idhep(nmxhep),
816 c &jmohep(2,nmxhep),jdahep(2,nmxhep),phep(5,nmxhep),vhep(4,nmxhep)
817 c SAVE /hepevt/
818 c common/phoqed/qedrad(nmxhep)
819 c LOGICAL qedrad
820 c SAVE /phoqed/
821  LOGICAL phflag
822 c
823  REAL*4 p4(4)
824 c
825 c check address mode
826  IF (n.EQ.0) THEN
827 c
828 c append mode
829  ihep=nhep+1
830  ELSE IF (n.GT.0) THEN
831 c
832 c absolute position
833  ihep=n
834  ELSE
835 c
836 c relative position
837  ihep=nhep+n
838  END IF
839 c
840 c check on ihep
841  IF ((ihep.LE.0).OR.(ihep.GT.nmxhep)) RETURN
842 c
843 c add entry
844  nhep=ihep
845  isthep(ihep)=ist
846  idhep(ihep)=id
847  jmohep(1,ihep)=jmo1
848  IF(jmo1.LT.0)jmohep(1,ihep)=jmohep(1,ihep)+ihep
849  jmohep(2,ihep)=jmo2
850  IF(jmo2.LT.0)jmohep(2,ihep)=jmohep(2,ihep)+ihep
851  jdahep(1,ihep)=jda1
852  jdahep(2,ihep)=jda2
853 c
854  DO i=1,4
855  phep(i,ihep)=p4(i)
856 c
857 c koral-b and koral-z do not provide vertex and/or lifetime informations
858  vhep(i,ihep)=0.0
859  END DO
860  phep(5,ihep)=pinv
861 c flag for photos...
862  qedrad(ihep)=phflag
863 c
864 c update process:
865  DO ip=jmohep(1,ihep),jmohep(2,ihep)
866  IF(ip.GT.0)THEN
867 c
868 c if there is a daughter at ihep, mother entry at ip has decayed
869  IF(isthep(ip).EQ.1)isthep(ip)=2
870 c
871 c and daughter pointers of mother entry must be updated
872  IF(jdahep(1,ip).EQ.0)THEN
873  jdahep(1,ip)=ihep
874  jdahep(2,ip)=ihep
875  ELSE
876  jdahep(2,ip)=max(ihep,jdahep(2,ip))
877  END IF
878  END IF
879  END DO
880 c
881  RETURN
882  END