C++InterfacetoTauola
F/standalone-F/taumain.f
1 /* copyright(c) 1991-2018 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 /* 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  PROGRAM TAUDEM
44 C **************
45 C NOTE THAT THE ROUTINES ARE NOT LIKE IN CPC DECK THIS IS HISTORICAL !!
46 C=======================================================================
47 C====================== DECTES : TEST OF TAU DECAY LIBRARY===========
48 C====================== KTORY = 1 : INTERFACE OF KORAL-Z TYPE ==========
49 C====================== KTORY = 2 : INTERFACE OF KORAL-B TYPE =========
50 C=======================================================================
51 C COMMON /PAWC/ BLAN(10000)
52  COMMON / / BLAN(10000)
53  CHARACTER*7 DNAME
54  COMMON / INOUT / INUT,IOUT
55  DNAME='kkpi'
56 ! CALL GLIMIT(20000)
57 ! CALL GOUTPU(16)
58  INUT=5
59  IOUT=6
60  OPEN(IOUT,FILE="./tauola.output")
61  OPEN(INUT,FILE="./dane.dat")
62  KTORY=1
63  CALL DECTES(KTORY)
64  KTORY=2
65  CALL DECTES(KTORY)
66  END
67  SUBROUTINE DECTES(KTORY)
68 C ************************
69  REAL POL(4)
70  DOUBLE PRECISION HH(4)
71 C SWITCHES FOR TAUOLA;
72  COMMON / JAKI / JAK1,JAK2,JAKP,JAKM,KTOM
73  COMMON / IDFC / IDFF
74 C I/O UNITS NUMBERS
75  COMMON / INOUT / INUT,IOUT
76 C LUND TYPE IDENTIFIER FOR A1
77  COMMON / IDPART / IA1
78 C /PTAU/ IS USED IN ROUTINE TRALO4
79  COMMON /PTAU/ PTAU
80  COMMON / TAURAD / XK0DEC,ITDKRC
81  REAL*8 XK0DEC
82  COMMON /TESTA1/ KEYA1
83 C special switch for tests of dGamma/dQ**2 in a1 decay
84 C KEYA1=1 constant width of a1 and rho
85 C KEYA1=2 free choice of rho propagator (defined in function FPIK)
86 C and free choice of a1 mass and width. function g(Q**2)
87 C (see formula 3.48 in Comp. Phys. Comm. 64 (1991) 275)
88 C hard coded both in Monte Carlo and in testing distribution.
89 C KEYA1=3 function g(Q**2) hardcoded in the Monte Carlo
90 C (it is timy to calculate!), but appropriately adjusted in
91 C testing distribution.
92 C-----------------------------------------------------------------------
93 C INITIALIZATION
94 C-----------------------------------------------------------------------
95 C======================================
96  NINP=INUT
97  NOUT=IOUT
98  3000 FORMAT(A80)
99  3001 FORMAT(8I2)
100  3002 FORMAT(I10)
101  3003 FORMAT(F10.0)
102 .EQ. IF (KTORY1) THEN
103  READ( NINP,3000) TESTIT
104  WRITE(NOUT,3000) TESTIT
105  READ( NINP,3001) KAT1,KAT2,KAT3,KAT4,KAT5,KAT6
106  READ( NINP,3002) NEVT,JAK1,JAK2,ITDKRC
107  READ( NINP,3003) PTAU,XK0DEC
108  ENDIF
109 C======================================
110 C control output
111  WRITE(NOUT,'(6a6/6i6)')
112  $ 'kat1','kat2','kat3','kat4','kat5','kat6',
113  $ KAT1 , KAT2 , KAT3 , KAT4 , KAT5 , KAT6
114  WRITE(NOUT,'(4a12/4i12)')
115  $ 'nevt','jak1','jak2','itdkrc',
116  $ NEVT, JAK1 , JAK2 , ITDKRC
117  WRITE(NOUT,'(2a12/2f12.6)')
118  $ 'ptau','xk0dec',
119  $ PTAU , XK0DEC
120 C======================================
121  JAK=0
122 C JAK1=5
123 C JAK2=5
124 C LUND IDENTIFIER (FOR TAU+) -15
125 .EQ. IF (KTORY1) THEN
126  IDFF=-15
127  ELSE
128  IDFF= 15
129  ENDIF
130 C KTO=1 DENOTES TAU DEFINED BY IDFF (I.E. TAU+)
131 C KTO=2 DENOTES THE OPPOSITE (I.E. TAU-)
132  KTO=2
133 .NE. IF (KTO2) THEN
134  PRINT *, 'for the sake of these tests kto has to be 2'
135  PRINT *, 'to change tau- to tau+ change idff from -15 to 15'
136  STOP
137  ENDIF
138 C TAU POLARIZATION IN ITS RESTFRAME;
139  POL(1)=0.
140  POL(2)=0.
141  POL(3)=.9
142 C TAU MOMENTUM IN GEV;
143 C PTAU=CMSENE/2.D0
144 C NUMBER OF EVENTS TO BE GENERATED;
145  NEVTES=10
146  NEVTES=NEVT
147  PRINT *, 'nevtes= ',NEVTES
148  WRITE(IOUT,7011) KEYA1
149 C
150 .EQ. IF (KTORY1) THEN
151  WRITE(IOUT,7001) JAK,IDFF,POL(3),PTAU
152  ELSE
153  WRITE(IOUT,7004) JAK,IDFF,POL(3),PTAU
154  ENDIF
155 C INITIALISATION OF TAU DECAY PACKAGE TAUOLA
156 C ******************************************
157  CALL INIMAS
158  CALL INITDK
159 
160 
161  CALL INIPHY(0.1D0)
162 .EQ. IF (KTORY1) THEN
163  CALL DEXAY(-1,POL)
164  ELSE
165  CALL DEKAY(-1,HH)
166  ENDIF
167 C-----------------------------------------------------------------------
168 C GENERATION
169 C-----------------------------------------------------------------------
170  NEV=0
171  DO 300 IEV=1,NEVTES
172  NEV=NEV+1
173 C RESLU INITIALISE THE LUND RECORD
174  CALL TAUFIL
175 C DECAY....
176 .EQ. IF (KTORY1) THEN
177  CALL DEXAY(KTO,POL)
178  ELSE
179  CALL DEKAY(KTO,HH)
180  CALL DEKAY(KTO+10,HH)
181  ENDIF
182  CALL LUHEPC(2)
183 .LE. IF(IEV44) THEN
184  WRITE(IOUT,7002) IEV
185 .NE. IF (KTORY1) THEN
186  WRITE(IOUT,7003) HH
187  ENDIF
188 C CALL LULIST(11)
189  CALL LULIST(2)
190  ENDIF
191  IPRI=MOD(NEV,1000)
192 .EQ. IF(IPRI1) PRINT *, ' event no: ',NEV,' nevtes: ',NEVTES
193  300 CONTINUE
194  301 CONTINUE
195 C-----------------------------------------------------------------------
196 C POSTGENERATION
197 C-----------------------------------------------------------------------
198 .EQ. IF (KTORY1) THEN
199  CALL DEXAY(100,POL)
200  ELSE
201  CALL DEKAY(100,HH)
202  ENDIF
203  RETURN
204  7001 FORMAT(//4(/1X,15(5H=====))
205  $ /,' ', 19X,' test of rad. corr in electron decay ',9X,1H ,
206  $ /,' ', 19X,' tests of tau decay routines ',9X,1H ,
207  $ /,' ', 19X,' INTERFACE of the koral-z TYPE ',9X,1H ,
208  $ 2(/,1X,15(5H=====)),
209  $ /,5X ,'jak =',I7 ,' key defining decay TYPE ',9X,1H ,
210  $ /,5X ,'idff =',I7 ,' lund identifier for first tau ',9X,1H ,
211  $ /,5X ,'pol(3)=',F7.2,' third component of tau polariz. ',9X,1H ,
212  $ /,5X ,'ptau =',F7.2,' third component of tau mom. gev ',9X,1H ,
213  $ 2(/,1X,15(5H=====))/)
214  7002 FORMAT(///1X, '===== event no.',I4,1X,5H=====)
215  7003 FORMAT(5X,'polarimetric vector: ',
216  $ 7X,'hh(1)',7X,'hh(2)',7X,'hh(3)',7X,'hh(4)',
217  $ /, 5X,' ', 4(1X,F11.8) )
218  7004 FORMAT(//4(/1X,15(5H=====))
219  $ /,' ', 19X,' test of rad. corr in electron decay ',9X,1H ,
220  $ /,' ', 19X,' tests of tau decay routines ',9X,1H ,
221  $ /,' ', 19X,' INTERFACE of the koral-b TYPE ',9X,1H ,
222  $ 2(/,1X,15(5H=====)),
223  $ /,5X ,'jak =',I7 ,' key defining decay TYPE ',9X,1H ,
224  $ /,5X ,'idff =',I7 ,' lund identifier for first tau ',9X,1H ,
225  $ /,5X ,'pol(3)=',F7.2,' third component of tau polariz. ',9X,1H ,
226  $ /,5X ,'ptau =',F7.2,' third component of tau mom. gev ',9X,1H ,
227  $ 2(/,1X,15(5H=====))/)
228  7011 FORMAT(///1X, '===== TYPE of current',I4,1X,5H=====)
229  END
230  SUBROUTINE CHOICE(MNUM,RR,ICHAN,PROB1,PROB2,PROB3,
231  $ AMRX,GAMRX,AMRA,GAMRA,AMRB,GAMRB)
232  COMMON / PARMAS / AMTAU,AMNUTA,AMEL,AMNUE,AMMU,AMNUMU
233  * ,AMPIZ,AMPI,AMRO,GAMRO,AMA1,GAMA1
234  * ,AMK,AMKZ,AMKST,GAMKST
235 C
236  REAL*4 AMTAU,AMNUTA,AMEL,AMNUE,AMMU,AMNUMU
237  * ,AMPIZ,AMPI,AMRO,GAMRO,AMA1,GAMA1
238  * ,AMK,AMKZ,AMKST,GAMKST
239 C
240  AMROP=1.1
241  GAMROP=0.36
242  AMOM=.782
243  GAMOM=0.0084
244 C XXXXA CORRESPOND TO S2 CHANNEL !
245 .EQ. IF(MNUM0) THEN
246  PROB1=0.5
247  PROB2=0.5
248  AMRX =AMA1
249  GAMRX=GAMA1
250  AMRA =AMRO
251  GAMRA=GAMRO
252  AMRB =AMRO
253  GAMRB=GAMRO
254 .EQ. ELSEIF(MNUM1) 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 .EQ. ELSEIF(MNUM2) THEN
264  PROB1=0.5
265  PROB2=0.5
266  AMRX =1.57
267  GAMRX=0.9
268  AMRB =AMKST
269  GAMRB=GAMKST
270  AMRA =AMRO
271  GAMRA=GAMRO
272 .EQ. ELSEIF(MNUM3) 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 .EQ. ELSEIF(MNUM4) 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 =AMKST
289  GAMRB=GAMKST
290 .EQ. ELSEIF(MNUM5) THEN
291  PROB1=0.5
292  PROB2=0.5
293  AMRX =1.27
294  GAMRX=0.3
295  AMRA =AMKST
296  GAMRA=GAMKST
297  AMRB =AMRO
298  GAMRB=GAMRO
299 .EQ. ELSEIF(MNUM6) THEN
300  PROB1=0.4
301  PROB2=0.4
302  AMRX =1.27
303  GAMRX=0.3
304  AMRA =AMRO
305  GAMRA=GAMRO
306  AMRB =AMKST
307  GAMRB=GAMKST
308 .EQ. ELSEIF(MNUM7) THEN
309  PROB1=0.0
310  PROB2=1.0
311  AMRX =1.27
312  GAMRX=0.9
313  AMRA =AMRO
314  GAMRA=GAMRO
315  AMRB =AMRO
316  GAMRB=GAMRO
317 .EQ. ELSEIF(MNUM8) THEN
318  PROB1=0.0
319  PROB2=1.0
320  AMRX =AMROP
321  GAMRX=GAMROP
322  AMRB =AMOM
323  GAMRB=GAMOM
324  AMRA =AMRO
325  GAMRA=GAMRO
326 .EQ. ELSEIF(MNUM101) THEN
327  PROB1=.35
328  PROB2=.35
329  AMRX =1.2
330  GAMRX=.46
331  AMRB =AMOM
332  GAMRB=GAMOM
333  AMRA =AMOM
334  GAMRA=GAMOM
335 .EQ. ELSEIF(MNUM102) THEN
336  PROB1=0.0
337  PROB2=0.0
338  AMRX =1.4
339  GAMRX=.6
340  AMRB =AMOM
341  GAMRB=GAMOM
342  AMRA =AMOM
343  GAMRA=GAMOM
344  ELSE
345  PROB1=0.0
346  PROB2=0.0
347  AMRX =AMA1
348  GAMRX=GAMA1
349  AMRA =AMRO
350  GAMRA=GAMRO
351  AMRB =AMRO
352  GAMRB=GAMRO
353  ENDIF
354 C
355 .LE. IF (RRPROB1) THEN
356  ICHAN=1
357 .LE. ELSEIF(RR(PROB1+PROB2)) THEN
358  ICHAN=2
359  AX =AMRA
360  GX =GAMRA
361  AMRA =AMRB
362  GAMRA=GAMRB
363  AMRB =AX
364  GAMRB=GX
365  PX =PROB1
366  PROB1=PROB2
367  PROB2=PX
368  ELSE
369  ICHAN=3
370  ENDIF
371 C
372  PROB3=1.0-PROB1-PROB2
373  END
374  SUBROUTINE INITDK
375 * ----------------------------------------------------------------------
376 * INITIALISATION OF TAU DECAY PARAMETERS and routines
377 *
378 * called by : KORALZ
379 * ----------------------------------------------------------------------
380 
381  COMMON / DECPAR / GFERMI,GV,GA,CCABIB,SCABIB,GAMEL
382  REAL*4 GFERMI,GV,GA,CCABIB,SCABIB,GAMEL
383  COMMON / PARMAS / AMTAU,AMNUTA,AMEL,AMNUE,AMMU,AMNUMU
384  * ,AMPIZ,AMPI,AMRO,GAMRO,AMA1,GAMA1
385  * ,AMK,AMKZ,AMKST,GAMKST
386 *
387  REAL*4 AMTAU,AMNUTA,AMEL,AMNUE,AMMU,AMNUMU
388  * ,AMPIZ,AMPI,AMRO,GAMRO,AMA1,GAMA1
389  * ,AMK,AMKZ,AMKST,GAMKST
390  COMMON / TAUBRA / GAMPRT(30),JLIST(30),NCHAN
391  COMMON / TAUKLE / BRA1,BRK0,BRK0B,BRKS
392  REAL*4 BRA1,BRK0,BRK0B,BRKS
393  PARAMETER (NMODE=15,NM1=0,NM2=1,NM3=8,NM4=2,NM5=1,NM6=3)
394  COMMON / TAUDCD /IDFFIN(9,NMODE),MULPIK(NMODE)
395  & ,NAMES
396  CHARACTER NAMES(NMODE)*31
397  CHARACTER OLDNAMES(7)*31
398  CHARACTER*80 bxINIT
399  PARAMETER (
400  $ bxINIT ='(1x,1h*,g17.8, 16x, a31,a4,a4, 1x,1h*)'
401  $ )
402  REAL*4 PI,POL1(4)
403 *
404 *
405 * LIST OF BRANCHING RATIOS
406 CAM normalised to e nu nutau channel
407 CAM enu munu pinu rhonu A1nu Knu K*nu pi
408 CAM DATA JLIST / 1, 2, 3, 4, 5, 6, 7,
409 *AM DATA GAMPRT /1.000,0.9730,0.6054,1.2432,0.8432,0.0432,O.O811,0.616
410 *AM
411 *AM multipion decays
412 *
413 * conventions of particles names
414 * K-,P-,K+, K0,P-,KB, K-,P0,K0
415 * 3, 1,-3 , 4, 1,-4 , 3, 2, 4 ,
416 * P0,P0,K-, K-,P-,P+, P-,KB,P0
417 * 2, 2, 3 , 3, 1,-1 , 1,-4, 2 ,
418 * ET,P-,P0 P-,P0,GM
419 * 9, 1, 2 , 1, 2, 8
420 *
421 C
422  DIMENSION NOPIK(6,NMODE),NPIK(NMODE)
423 *AM outgoing multiplicity and flavors of multi-pion /multi-K modes
424  DATA NPIK / 4, 4,
425  1 5, 5,
426  2 6, 6,
427  3 3, 3,
428  4 3, 3,
429  5 3, 3,
430  6 3, 3,
431  7 2 /
432  DATA NOPIK / -1,-1, 1, 2, 0, 0, 2, 2, 2,-1, 0, 0,
433  1 -1,-1, 1, 2, 2, 0, -1,-1,-1, 1, 1, 0,
434  2 -1,-1,-1, 1, 1, 2, -1,-1, 1, 2, 2, 2,
435  3 -3,-1, 3, 0, 0, 0, -4,-1, 4, 0, 0, 0,
436  4 -3, 2,-4, 0, 0, 0, 2, 2,-3, 0, 0, 0,
437  5 -3,-1, 1, 0, 0, 0, -1, 4, 2, 0, 0, 0,
438  6 9,-1, 2, 0, 0, 0, -1, 2, 8, 0, 0, 0,
439 C AJWMOD fix sign bug, 2/22/99
440  7 -3,-4, 0, 0, 0, 0 /
441 * LIST OF BRANCHING RATIOS
442  NCHAN = NMODE + 7
443  DO 1 I = 1,30
444 .LE. IF (INCHAN) THEN
445  JLIST(I) = I
446 .EQ. IF(I 1) GAMPRT(I) =0.1800
447 .EQ. IF(I 2) GAMPRT(I) =0.1751
448 .EQ. IF(I 3) GAMPRT(I) =0.1110
449 .EQ. IF(I 4) GAMPRT(I) =0.2515
450 .EQ. IF(I 5) GAMPRT(I) =0.1790
451 .EQ. IF(I 6) GAMPRT(I) =0.0071
452 .EQ. IF(I 7) GAMPRT(I) =0.0134
453 .EQ. IF(I 8) GAMPRT(I) =0.0450
454 .EQ. IF(I 9) GAMPRT(I) =0.0100
455 .EQ. IF(I10) GAMPRT(I) =0.0009
456 .EQ. IF(I11) GAMPRT(I) =0.0004
457 .EQ. IF(I12) GAMPRT(I) =0.0003
458 .EQ. IF(I13) GAMPRT(I) =0.0005
459 .EQ. IF(I14) GAMPRT(I) =0.0015
460 .EQ. IF(I15) GAMPRT(I) =0.0015
461 .EQ. IF(I16) GAMPRT(I) =0.0015
462 .EQ. IF(I17) GAMPRT(I) =0.0005
463 .EQ. IF(I18) GAMPRT(I) =0.0050
464 .EQ. IF(I19) GAMPRT(I) =0.0055
465 .EQ. IF(I20) GAMPRT(I) =0.0017
466 .EQ. IF(I21) GAMPRT(I) =0.0013
467 .EQ. IF(I22) GAMPRT(I) =0.0010
468 .EQ. IF(I 1) OLDNAMES(I)=' tau- --> e- '
469 .EQ. IF(I 2) OLDNAMES(I)=' tau- --> mu- '
470 .EQ. IF(I 3) OLDNAMES(I)=' tau- --> pi- '
471 .EQ. IF(I 4) OLDNAMES(I)=' tau- --> pi-, pi0 '
472 .EQ. IF(I 5) OLDNAMES(I)=' tau- --> a1- (two subch) '
473 .EQ. IF(I 6) OLDNAMES(I)=' tau- --> k- '
474 .EQ. IF(I 7) OLDNAMES(I)=' tau- --> k*- (two subch) '
475 .EQ. IF(I 8) NAMES(I-7)=' tau- --> 2pi-, pi0, pi+ '
476 .EQ. IF(I 9) NAMES(I-7)=' tau- --> 3pi0, pi- '
477 .EQ. IF(I10) NAMES(I-7)=' tau- --> 2pi-, pi+, 2pi0 '
478 .EQ. IF(I11) NAMES(I-7)=' tau- --> 3pi-, 2pi+, '
479 .EQ. IF(I12) NAMES(I-7)=' tau- --> 3pi-, 2pi+, pi0 '
480 .EQ. IF(I13) NAMES(I-7)=' tau- --> 2pi-, pi+, 3pi0 '
481 .EQ. IF(I14) NAMES(I-7)=' tau- --> k-, pi-, k+ '
482 .EQ. IF(I15) NAMES(I-7)=' tau- --> k0, pi-, k0b '
483 .EQ. IF(I16) NAMES(I-7)=' tau- --> k-, k0, pi0 '
484 .EQ. IF(I17) NAMES(I-7)=' tau- --> pi0 pi0 k- '
485 .EQ. IF(I18) NAMES(I-7)=' tau- --> k- pi- pi+ '
486 .EQ. IF(I19) NAMES(I-7)=' tau- --> pi- k0b pi0 '
487 .EQ. IF(I20) NAMES(I-7)=' tau- --> eta pi- pi0 '
488 .EQ. IF(I21) NAMES(I-7)=' tau- --> pi- pi0 gam '
489 .EQ. IF(I22) NAMES(I-7)=' tau- --> k- k0 '
490  ELSE
491  JLIST(I) = 0
492  GAMPRT(I) = 0.
493  ENDIF
494  1 CONTINUE
495  DO I=1,NMODE
496  MULPIK(I)=NPIK(I)
497  DO J=1,MULPIK(I)
498  IDFFIN(J,I)=NOPIK(J,I)
499  ENDDO
500  ENDDO
501 *
502 *
503 * --- COEFFICIENTS TO FIX RATIO OF:
504 * --- A1 3CHARGED/ A1 1CHARGED 2 NEUTRALS MATRIX ELEMENTS (MASLESS LIM.)
505 * --- PROBABILITY OF K0 TO BE KS
506 * --- PROBABILITY OF K0B TO BE KS
507 * --- RATIO OF COEFFICIENTS FOR K*--> K0 PI-
508 * --- ALL COEFFICENTS SHOULD BE IN THE RANGE (0.0,1.0)
509 * --- THEY MEANING IS PROBABILITY OF THE FIRST CHOICE ONLY IF ONE
510 * --- NEGLECTS MASS-PHASE SPACE EFFECTS
511  BRA1=0.5
512  BRK0=0.5
513  BRK0B=0.5
514  BRKS=0.6667
515 *
516 
517  GFERMI = 1.16637E-5
518  CCABIB = 0.975
519  GV = 1.0
520  GA =-1.0
521 
522 
523 
524 * ZW 13.04.89 HERE WAS AN ERROR
525  SCABIB = SQRT(1.-CCABIB**2)
526  PI =4.*ATAN(1.)
527  GAMEL = GFERMI**2*AMTAU**5/(192*PI**3)
528 *
529 * CALL DEXAY(-1,pol1)
530 *
531  RETURN
532  END
533  FUNCTION DCDMAS(IDENT)
534  COMMON / PARMAS / AMTAU,AMNUTA,AMEL,AMNUE,AMMU,AMNUMU
535  * ,AMPIZ,AMPI,AMRO,GAMRO,AMA1,GAMA1
536  * ,AMK,AMKZ,AMKST,GAMKST
537 *
538  REAL*4 AMTAU,AMNUTA,AMEL,AMNUE,AMMU,AMNUMU
539  * ,AMPIZ,AMPI,AMRO,GAMRO,AMA1,GAMA1
540  * ,AMK,AMKZ,AMKST,GAMKST
541 .EQ. IF (IDENT 1) THEN
542  APKMAS=AMPI
543 .EQ. ELSEIF (IDENT-1) THEN
544  APKMAS=AMPI
545 .EQ. ELSEIF (IDENT 2) THEN
546  APKMAS=AMPIZ
547 .EQ. ELSEIF (IDENT-2) THEN
548  APKMAS=AMPIZ
549 .EQ. ELSEIF (IDENT 3) THEN
550  APKMAS=AMK
551 .EQ. ELSEIF (IDENT-3) THEN
552  APKMAS=AMK
553 .EQ. ELSEIF (IDENT 4) THEN
554  APKMAS=AMKZ
555 .EQ. ELSEIF (IDENT-4) THEN
556  APKMAS=AMKZ
557 .EQ. ELSEIF (IDENT 8) THEN
558  APKMAS=0.0001
559 .EQ. ELSEIF (IDENT-8) THEN
560  APKMAS=0.0001
561 .EQ. ELSEIF (IDENT 9) THEN
562  APKMAS=0.5488
563 .EQ. ELSEIF (IDENT-9) THEN
564  APKMAS=0.5488
565  ELSE
566  PRINT *, 'stop in apkmas, wrong ident=',IDENT
567  STOP
568  ENDIF
569  DCDMAS=APKMAS
570  END
571  FUNCTION LUNPIK(ID,ISGN)
572  COMMON / TAUKLE / BRA1,BRK0,BRK0B,BRKS
573  REAL*4 BRA1,BRK0,BRK0B,BRKS
574  REAL*4 XIO(1)
575  IDENT=ID*ISGN
576 .EQ. IF (IDENT 1) THEN
577  IPKDEF=-211
578 .EQ. ELSEIF (IDENT-1) THEN
579  IPKDEF= 211
580 .EQ. ELSEIF (IDENT 2) THEN
581  IPKDEF=111
582 .EQ. ELSEIF (IDENT-2) THEN
583  IPKDEF=111
584 .EQ. ELSEIF (IDENT 3) THEN
585  IPKDEF=-321
586 .EQ. ELSEIF (IDENT-3) THEN
587  IPKDEF= 321
588 .EQ. ELSEIF (IDENT 4) THEN
589 *
590 * K0 --> K0_LONG (IS 130) / K0_SHORT (IS 310) = 1/1
591  CALL RANMAR(XIO,1)
592 .GT. IF (XIO(1)BRK0) THEN
593  IPKDEF= 130
594  ELSE
595  IPKDEF= 310
596  ENDIF
597 .EQ. ELSEIF (IDENT-4) THEN
598 *
599 * K0B--> K0_LONG (IS 130) / K0_SHORT (IS 310) = 1/1
600  CALL RANMAR(XIO,1)
601 .GT. IF (XIO(1)BRK0B) THEN
602  IPKDEF= 130
603  ELSE
604  IPKDEF= 310
605  ENDIF
606 .EQ. ELSEIF (IDENT 8) THEN
607  IPKDEF= 22
608 .EQ. ELSEIF (IDENT-8) THEN
609  IPKDEF= 22
610 .EQ. ELSEIF (IDENT 9) THEN
611  IPKDEF= 221
612 .EQ. ELSEIF (IDENT-9) THEN
613  IPKDEF= 221
614  ELSE
615  PRINT *, 'stop in ipkdef, wrong ident=',IDENT
616  STOP
617  ENDIF
618  LUNPIK=IPKDEF
619  END
620 
621 
622 
623  SUBROUTINE TAURDF(KTO)
624 C THIS ROUTINE CAN BE CALLED BEFORE ANY TAU+ OR TAU- EVENT IS GENERATED
625 C IT CAN BE USED TO GENERATE TAU+ AND TAU- SAMPLES OF DIFFERENT
626 C CONTENTS
627  COMMON / TAUKLE / BRA1,BRK0,BRK0B,BRKS
628  REAL*4 BRA1,BRK0,BRK0B,BRKS
629  COMMON / TAUBRA / GAMPRT(30),JLIST(30),NCHAN
630 .EQ. IF (KTO1) THEN
631 C ==================
632 C AJWMOD: Set the BRs for (A1+ -> rho+ pi0) and (K*+ -> K0 pi+)
633  BRA1 = PKORB(4,1)
634  BRKS = PKORB(4,3)
635  BRK0 = PKORB(4,5)
636  BRK0B = PKORB(4,6)
637  ELSE
638 C ====
639 C AJWMOD: Set the BRs for (A1+ -> rho+ pi0) and (K*+ -> K0 pi+)
640  BRA1 = PKORB(4,2)
641  BRKS = PKORB(4,4)
642  BRK0 = PKORB(4,5)
643  BRK0B = PKORB(4,6)
644  ENDIF
645 C =====
646  END
647 
648  SUBROUTINE INIPHY(XK00)
649 * ----------------------------------------------------------------------
650 * INITIALISATION OF PARAMETERS
651 * USED IN QED and/or GSW ROUTINES
652 * ----------------------------------------------------------------------
653  COMMON / QEDPRM /ALFINV,ALFPI,XK0
654  REAL*8 ALFINV,ALFPI,XK0
655  REAL*8 PI8,XK00
656 *
657  PI8 = 4.D0*DATAN(1.D0)
658  ALFINV = 137.03604D0
659  ALFPI = 1D0/(ALFINV*PI8)
660  XK0=XK00
661  END
662 
663  SUBROUTINE INIMAS
664 C ----------------------------------------------------------------------
665 C INITIALISATION OF MASSES
666 C
667 C called by : KORALZ
668 C ----------------------------------------------------------------------
669  COMMON / PARMAS / AMTAU,AMNUTA,AMEL,AMNUE,AMMU,AMNUMU
670  * ,AMPIZ,AMPI,AMRO,GAMRO,AMA1,GAMA1
671  * ,AMK,AMKZ,AMKST,GAMKST
672 *
673  REAL*4 AMTAU,AMNUTA,AMEL,AMNUE,AMMU,AMNUMU
674  * ,AMPIZ,AMPI,AMRO,GAMRO,AMA1,GAMA1
675  * ,AMK,AMKZ,AMKST,GAMKST
676 C
677 C IN-COMING / OUT-GOING FERMION MASSES
678  AMTAU = 1.7842
679 C --- let us update tau mass ...
680  AMTAU = 1.777
681  AMNUTA = 0.010
682  AMEL = 0.0005111
683  AMNUE = 0.0
684  AMMU = 0.105659
685  AMNUMU = 0.0
686 *
687 * MASSES USED IN TAU DECAYS
688  AMPIZ = 0.134964
689  AMPI = 0.139568
690  AMRO = 0.773
691  GAMRO = 0.145
692 *C GAMRO = 0.666
693  AMA1 = 1.251
694  GAMA1 = 0.599
695  AMK = 0.493667
696  AMKZ = 0.49772
697  AMKST = 0.8921
698  GAMKST = 0.0513
699 C
700 C
701 C IN-COMING / OUT-GOING FERMION MASSES
702 !! AMNUTA = PKORB(1,2)
703 !! AMNUE = PKORB(1,4)
704 !! AMNUMU = PKORB(1,6)
705 C
706 C MASSES USED IN TAU DECAYS Cleo settings
707 !! AMPIZ = PKORB(1,7)
708 !! AMPI = PKORB(1,8)
709 !! AMRO = PKORB(1,9)
710 !! GAMRO = PKORB(2,9)
711  AMA1 = 1.275 !! PKORB(1,10)
712  GAMA1 = 0.615 !! PKORB(2,10)
713 !! AMK = PKORB(1,11)
714 !! AMKZ = PKORB(1,12)
715 !! AMKST = PKORB(1,13)
716 !! GAMKST = PKORB(2,13)
717 C
718 
719  RETURN
720  END
721  SUBROUTINE TAUFIL
722 C *****************
723 C SUBSITUTE OF tau PRODUCTION GENERATOR
724 C
725  COMMON / PARMAS / AMTAU,AMNUTA,AMEL,AMNUE,AMMU,AMNUMU
726  * ,AMPIZ,AMPI,AMRO,GAMRO,AMA1,GAMA1
727  * ,AMK,AMKZ,AMKST,GAMKST
728 C
729  REAL*4 AMTAU,AMNUTA,AMEL,AMNUE,AMMU,AMNUMU
730  * ,AMPIZ,AMPI,AMRO,GAMRO,AMA1,GAMA1
731  * ,AMK,AMKZ,AMKST,GAMKST
732  COMMON / IDFC / IDFF
733 C positions of taus in the LUND common block
734 C it will be used by TAUOLA output routines.
735  COMMON /TAUPOS / NPA,NPB
736  DIMENSION XPB1(4),XPB2(4),AQF1(4),AQF2(4)
737 C
738 C --- DEFINING DUMMY EVENTS MOMENTA
739  DO 4 K=1,3
740  XPB1(K)=0.0
741  XPB2(K)=0.0
742  AQF1(K)=0.0
743  AQF2(K)=0.0
744  4 CONTINUE
745  AQF1(4)=AMTAU
746  AQF2(4)=AMTAU
747 C --- TAU MOMENTA
748  CALL TRALO4(1,AQF1,AQF1,AM)
749  CALL TRALO4(2,AQF2,AQF2,AM)
750 C --- BEAMS MOMENTA AND IDENTIFIERS
751  KFB1= 11*IDFF/IABS(IDFF)
752  KFB2=-11*IDFF/IABS(IDFF)
753  XPB1(4)= AQF1(4)
754  XPB1(3)= AQF1(4)
755 .NE. IF(AQF1(3)0.0)
756  $ XPB1(3)= AQF1(4)*AQF1(3)/ABS(AQF1(3))
757  XPB2(4)= AQF2(4)
758  XPB2(3)=-AQF2(4)
759 .NE. IF(AQF2(3)0.0)
760  $ XPB2(3)= AQF2(4)*AQF2(3)/ABS(AQF2(3))
761 C --- Position of first and second tau in LUND common
762  NPA=3
763  NPB=4
764 C --- FILL TO LUND COMMON
765  CALL FILHEP( 1,3, KFB1,0,0,0,0,XPB1, AMEL,.TRUE.)
766  CALL FILHEP( 2,3, KFB2,0,0,0,0,XPB2, AMEL,.TRUE.)
767  CALL FILHEP(NPA,1, IDFF,1,2,0,0,AQF1,AMTAU,.TRUE.)
768  CALL FILHEP(NPB,1,-IDFF,1,2,0,0,AQF2,AMTAU,.TRUE.)
769  END
770  SUBROUTINE TRALO4(KTO,P,Q,AM)
771 C **************************
772 C SUBSITUTE OF TRALO4
773  REAL P(4),Q(4)
774 C
775  COMMON / PARMAS / AMTAU,AMNUTA,AMEL,AMNUE,AMMU,AMNUMU
776  * ,AMPIZ,AMPI,AMRO,GAMRO,AMA1,GAMA1
777  * ,AMK,AMKZ,AMKST,GAMKST
778 C
779  REAL*4 AMTAU,AMNUTA,AMEL,AMNUE,AMMU,AMNUMU
780  * ,AMPIZ,AMPI,AMRO,GAMRO,AMA1,GAMA1
781  * ,AMK,AMKZ,AMKST,GAMKST
782  COMMON /PTAU/ PTAU
783  AM=AMAS4(P)
784  ETAU=SQRT(PTAU**2+AMTAU**2)
785  EXE=(ETAU+PTAU)/AMTAU
786 .EQ. IF(KTO2) EXE=(ETAU-PTAU)/AMTAU
787  CALL BOSTR3(EXE,P,Q)
788 C ======================================================================
789 C END OF THE TEST JOB
790 C ======================================================================
791  END
792  SUBROUTINE FILHEP(N,IST,ID,JMO1,JMO2,JDA1,JDA2,P4,PINV,PHFLAG)
793 C ----------------------------------------------------------------------
794 C this subroutine fills one entry into the HEPEVT common
795 C and updates the information for affected mother entries
796 C
797 C written by Martin W. Gruenewald (91/01/28)
798 C
799 C called by : ZTOHEP,BTOHEP,DWLUxy
800 C ----------------------------------------------------------------------
801 C
802 C this is the hepevt class in old style. No d_h_ class pre-name
803  INTEGER NMXHEP
804  PARAMETER (NMXHEP=10000)
805  REAL*8 phep, vhep ! to be real*4/ *8 depending on host
806  INTEGER nevhep,nhep,isthep,idhep,jmohep,
807  $ jdahep
808  COMMON /hepevt/
809  $ nevhep, ! serial number
810  $ nhep, ! number of particles
811  $ isthep(nmxhep), ! status code
812  $ idhep(nmxhep), ! particle ident KF
813  $ jmohep(2,nmxhep), ! parent particles
814  $ jdahep(2,nmxhep), ! childreen particles
815  $ phep(5,nmxhep), ! four-momentum, mass [GeV]
816  $ vhep(4,nmxhep) ! vertex [mm]
817 * ----------------------------------------------------------------------
818  LOGICAL qedrad
819  COMMON /phoqed/
820  $ qedrad(nmxhep) ! Photos flag
821 * ----------------------------------------------------------------------
822  SAVE hepevt,phoqed
823 C PARAMETER (NMXHEP=2000)
824 C COMMON/HEPEVT/NEVHEP,NHEP,ISTHEP(NMXHEP),IDHEP(NMXHEP),
825 C &JMOHEP(2,NMXHEP),JDAHEP(2,NMXHEP),PHEP(5,NMXHEP),VHEP(4,NMXHEP)
826 C SAVE /HEPEVT/
827 C COMMON/PHOQED/QEDRAD(NMXHEP)
828 C LOGICAL QEDRAD
829 C SAVE /PHOQED/
830  LOGICAL PHFLAG
831 C
832  REAL*4 P4(4)
833 C
834 C check address mode
835 .EQ. IF (N0) THEN
836 C
837 C append mode
838  IHEP=NHEP+1
839 .GT. ELSE IF (N0) THEN
840 C
841 C absolute position
842  IHEP=N
843  ELSE
844 C
845 C relative position
846  IHEP=NHEP+N
847  END IF
848 C
849 C check on IHEP
850 .LE..OR..GT. IF ((IHEP0)(IHEPNMXHEP)) RETURN
851 C
852 C add entry
853  NHEP=IHEP
854  ISTHEP(IHEP)=IST
855  IDHEP(IHEP)=ID
856  JMOHEP(1,IHEP)=JMO1
857 .LT. IF(JMO10)JMOHEP(1,IHEP)=JMOHEP(1,IHEP)+IHEP
858  JMOHEP(2,IHEP)=JMO2
859 .LT. IF(JMO20)JMOHEP(2,IHEP)=JMOHEP(2,IHEP)+IHEP
860  JDAHEP(1,IHEP)=JDA1
861  JDAHEP(2,IHEP)=JDA2
862 C
863  DO I=1,4
864  PHEP(I,IHEP)=P4(I)
865 C
866 C KORAL-B and KORAL-Z do not provide vertex and/or lifetime informations
867  VHEP(I,IHEP)=0.0
868  END DO
869  PHEP(5,IHEP)=PINV
870 C FLAG FOR PHOTOS...
871  QEDRAD(IHEP)=PHFLAG
872 C
873 C update process:
874  DO IP=JMOHEP(1,IHEP),JMOHEP(2,IHEP)
875 .GT. IF(IP0)THEN
876 C
877 C if there is a daughter at IHEP, mother entry at IP has decayed
878 .EQ. IF(ISTHEP(IP)1)ISTHEP(IP)=2
879 C
880 C and daughter pointers of mother entry must be updated
881 .EQ. IF(JDAHEP(1,IP)0)THEN
882  JDAHEP(1,IP)=IHEP
883  JDAHEP(2,IP)=IHEP
884  ELSE
885  JDAHEP(2,IP)=MAX(IHEP,JDAHEP(2,IP))
886  END IF
887  END IF
888  END DO
889 C
890  RETURN
891  END