C++ Interface to Tauola
tauola/demo-standalone/taumain.f
1/* copyright(c) 1991-2023 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 PROGRAM TAUDEM
44C **************
45C NOTE THAT THE ROUTINES ARE NOT LIKE IN CPC DECK THIS IS HISTORICAL !!
46C=======================================================================
47C====================== DECTES : TEST OF TAU DECAY LIBRARY===========
48C====================== KTORY = 1 : INTERFACE OF KORAL-Z TYPE ==========
49C====================== KTORY = 2 : INTERFACE OF KORAL-B TYPE =========
50C=======================================================================
51C 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)
68C ************************
69 REAL POL(4)
70 DOUBLE PRECISION HH(4)
71C SWITCHES FOR TAUOLA;
72 COMMON / JAKI / JAK1,JAK2,JAKP,JAKM,KTOM
73 COMMON / IDFC / IDFF
74C I/O UNITS NUMBERS
75 COMMON / INOUT / INUT,IOUT
76C LUND TYPE IDENTIFIER FOR A1
77 COMMON / IDPART / IA1
78C /PTAU/ IS USED IN ROUTINE TRALO4
79 COMMON /PTAU/ PTAU
80 COMMON / TAURAD / XK0DEC,ITDKRC
81 REAL*8 XK0DEC
82 COMMON /TESTA1/ KEYA1
83C special switch for tests of dGamma/dQ**2 in a1 decay
84C KEYA1=1 constant width of a1 and rho
85C KEYA1=2 free choice of rho propagator (defined in function FPIK)
86C and free choice of a1 mass and width. function g(Q**2)
87C (see formula 3.48 in Comp. Phys. Comm. 64 (1991) 275)
88C hard coded both in Monte Carlo and in testing distribution.
89C KEYA1=3 function g(Q**2) hardcoded in the Monte Carlo
90C (it is timy to calculate!), but appropriately adjusted in
91C testing distribution.
92C-----------------------------------------------------------------------
93C INITIALIZATION
94C-----------------------------------------------------------------------
95C======================================
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
109C======================================
110C 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
120C======================================
121 JAK=0
122C JAK1=5
123C JAK2=5
124C LUND IDENTIFIER (FOR TAU+) -15
125.EQ. IF (KTORY1) THEN
126 IDFF=-15
127 ELSE
128 IDFF= 15
129 ENDIF
130C KTO=1 DENOTES TAU DEFINED BY IDFF (I.E. TAU+)
131C 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
138C TAU POLARIZATION IN ITS RESTFRAME;
139 POL(1)=0.
140 POL(2)=0.
141 POL(3)=.9
142C TAU MOMENTUM IN GEV;
143C PTAU=CMSENE/2.D0
144C NUMBER OF EVENTS TO BE GENERATED;
145 NEVTES=10
146 NEVTES=NEVT
147 PRINT *, 'nevtes= ',NEVTES
148 WRITE(IOUT,7011) KEYA1
149C
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
155C INITIALISATION OF TAU DECAY PACKAGE TAUOLA
156C ******************************************
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
167C-----------------------------------------------------------------------
168C GENERATION
169C-----------------------------------------------------------------------
170 NEV=0
171 DO 300 IEV=1,NEVTES
172 NEV=NEV+1
173C RESLU INITIALISE THE LUND RECORD
174 CALL TAUFIL
175C 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
188C 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
195C-----------------------------------------------------------------------
196C POSTGENERATION
197C-----------------------------------------------------------------------
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
235C
236 REAL*4 AMTAU,AMNUTA,AMEL,AMNUE,AMMU,AMNUMU
237 * ,AMPIZ,AMPI,AMRO,GAMRO,AMA1,GAMA1
238 * ,AMK,AMKZ,AMKST,GAMKST
239C
240 AMROP=1.1
241 GAMROP=0.36
242 AMOM=.782
243 GAMOM=0.0084
244C 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
354C
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
371C
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
406CAM normalised to e nu nutau channel
407CAM enu munu pinu rhonu A1nu Knu K*nu pi
408CAM 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*
421C
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,
439C 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)
624C THIS ROUTINE CAN BE CALLED BEFORE ANY TAU+ OR TAU- EVENT IS GENERATED
625C IT CAN BE USED TO GENERATE TAU+ AND TAU- SAMPLES OF DIFFERENT
626C 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
631C ==================
632C 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
638C ====
639C 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
645C =====
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
664C ----------------------------------------------------------------------
665C INITIALISATION OF MASSES
666C
667C called by : KORALZ
668C ----------------------------------------------------------------------
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
676C
677C IN-COMING / OUT-GOING FERMION MASSES
678 AMTAU = 1.7842
679C --- 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
699C
700C
701C IN-COMING / OUT-GOING FERMION MASSES
702!! AMNUTA = PKORB(1,2)
703!! AMNUE = PKORB(1,4)
704!! AMNUMU = PKORB(1,6)
705C
706C 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)
717C
718
719 RETURN
720 END
721 SUBROUTINE TAUFIL
722C *****************
723C SUBSITUTE OF tau PRODUCTION GENERATOR
724C
725 COMMON / PARMAS / AMTAU,AMNUTA,AMEL,AMNUE,AMMU,AMNUMU
726 * ,AMPIZ,AMPI,AMRO,GAMRO,AMA1,GAMA1
727 * ,AMK,AMKZ,AMKST,GAMKST
728C
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
733C positions of taus in the LUND common block
734C it will be used by TAUOLA output routines.
735 COMMON /TAUPOS / NPA,NPB
736 DIMENSION XPB1(4),XPB2(4),AQF1(4),AQF2(4)
737C
738C --- 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
747C --- TAU MOMENTA
748 CALL TRALO4(1,AQF1,AQF1,AM)
749 CALL TRALO4(2,AQF2,AQF2,AM)
750C --- 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))
761C --- Position of first and second tau in LUND common
762 NPA=3
763 NPB=4
764C --- 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)
771C **************************
772C SUBSITUTE OF TRALO4
773 REAL P(4),Q(4)
774C
775 COMMON / PARMAS / AMTAU,AMNUTA,AMEL,AMNUE,AMMU,AMNUMU
776 * ,AMPIZ,AMPI,AMRO,GAMRO,AMA1,GAMA1
777 * ,AMK,AMKZ,AMKST,GAMKST
778C
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)
788C ======================================================================
789C END OF THE TEST JOB
790C ======================================================================
791 END
792 SUBROUTINE FILHEP(N,IST,ID,JMO1,JMO2,JDA1,JDA2,P4,PINV,PHFLAG)
793C ----------------------------------------------------------------------
794C this subroutine fills one entry into the HEPEVT common
795C and updates the information for affected mother entries
796C
797C written by Martin W. Gruenewald (91/01/28)
798C
799C called by : ZTOHEP,BTOHEP,DWLUxy
800C ----------------------------------------------------------------------
801C
802C 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
823C PARAMETER (NMXHEP=2000)
824C COMMON/HEPEVT/NEVHEP,NHEP,ISTHEP(NMXHEP),IDHEP(NMXHEP),
825C &JMOHEP(2,NMXHEP),JDAHEP(2,NMXHEP),PHEP(5,NMXHEP),VHEP(4,NMXHEP)
826C SAVE /HEPEVT/
827C COMMON/PHOQED/QEDRAD(NMXHEP)
828C LOGICAL QEDRAD
829C SAVE /PHOQED/
830 LOGICAL PHFLAG
831C
832 REAL*4 P4(4)
833C
834C check address mode
835.EQ. IF (N0) THEN
836C
837C append mode
838 IHEP=NHEP+1
839.GT. ELSE IF (N0) THEN
840C
841C absolute position
842 IHEP=N
843 ELSE
844C
845C relative position
846 IHEP=NHEP+N
847 END IF
848C
849C check on IHEP
850.LE..OR..GT. IF ((IHEP0)(IHEPNMXHEP)) RETURN
851C
852C 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
862C
863 DO I=1,4
864 PHEP(I,IHEP)=P4(I)
865C
866C 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
870C FLAG FOR PHOTOS...
871 QEDRAD(IHEP)=PHFLAG
872C
873C update process:
874 DO IP=JMOHEP(1,IHEP),JMOHEP(2,IHEP)
875.GT. IF(IP0)THEN
876C
877C if there is a daughter at IHEP, mother entry at IP has decayed
878.EQ. IF(ISTHEP(IP)1)ISTHEP(IP)=2
879C
880C 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
889C
890 RETURN
891 END