TEvtGen/Tauola/tauola_extras.f

   1
   2       SUBROUTINE CHOICE(MNUM,RR,ICHAN,PROB1,PROB2,PROB3,
   3      $            AMRX,GAMRX,AMRA,GAMRA,AMRB,GAMRB)
   4       COMMON / PARMAS / AMTAU,AMNUTA,AMEL,AMNUE,AMMU,AMNUMU
   5      *                 ,AMPIZ,AMPI,AMRO,GAMRO,AMA1,GAMA1
   6      *                 ,AMK,AMKZ,AMKST,GAMKST
   7 C
   8       REAL*4            AMTAU,AMNUTA,AMEL,AMNUE,AMMU,AMNUMU
   9      *                 ,AMPIZ,AMPI,AMRO,GAMRO,AMA1,GAMA1
  10      *                 ,AMK,AMKZ,AMKST,GAMKST
  11 C
  12       AMROP=1.1
  13       GAMROP=0.36
  14       AMOM=.782
  15       GAMOM=0.0084
  16 C     XXXXA CORRESPOND TO S2 CHANNEL !
  17       IF(MNUM.EQ.0) THEN
  18        PROB1=0.5
  19        PROB2=0.5
  20        AMRX =AMA1
  21        GAMRX=GAMA1
  22        AMRA =AMRO
  23        GAMRA=GAMRO
  24        AMRB =AMRO
  25        GAMRB=GAMRO
  26       ELSEIF(MNUM.EQ.1) THEN
  27        PROB1=0.5
  28        PROB2=0.5
  29        AMRX =1.57
  30        GAMRX=0.9
  31        AMRB =AMKST
  32        GAMRB=GAMKST
  33        AMRA =AMRO
  34        GAMRA=GAMRO
  35       ELSEIF(MNUM.EQ.2) THEN
  36        PROB1=0.5
  37        PROB2=0.5
  38        AMRX =1.57
  39        GAMRX=0.9
  40        AMRB =AMKST
  41        GAMRB=GAMKST
  42        AMRA =AMRO
  43        GAMRA=GAMRO
  44       ELSEIF(MNUM.EQ.3) THEN
  45        PROB1=0.5
  46        PROB2=0.5
  47        AMRX =1.27
  48        GAMRX=0.3
  49        AMRA =AMKST
  50        GAMRA=GAMKST
  51        AMRB =AMKST
  52        GAMRB=GAMKST
  53       ELSEIF(MNUM.EQ.4) THEN
  54        PROB1=0.5
  55        PROB2=0.5
  56        AMRX =1.27
  57        GAMRX=0.3
  58        AMRA =AMKST
  59        GAMRA=GAMKST
  60        AMRB =AMKST
  61        GAMRB=GAMKST
  62       ELSEIF(MNUM.EQ.5) THEN
  63        PROB1=0.5
  64        PROB2=0.5
  65        AMRX =1.27
  66        GAMRX=0.3
  67        AMRA =AMKST
  68        GAMRA=GAMKST
  69        AMRB =AMRO
  70        GAMRB=GAMRO
  71       ELSEIF(MNUM.EQ.6) THEN
  72        PROB1=0.4
  73        PROB2=0.4
  74        AMRX =1.27
  75        GAMRX=0.3
  76        AMRA =AMRO
  77        GAMRA=GAMRO
  78        AMRB =AMKST
  79        GAMRB=GAMKST
  80       ELSEIF(MNUM.EQ.7) THEN
  81        PROB1=0.0
  82        PROB2=1.0
  83        AMRX =1.27
  84        GAMRX=0.9
  85        AMRA =AMRO
  86        GAMRA=GAMRO
  87        AMRB =AMRO
  88        GAMRB=GAMRO
  89       ELSEIF(MNUM.EQ.8) THEN
  90        PROB1=0.0
  91        PROB2=1.0
  92        AMRX =AMROP
  93        GAMRX=GAMROP
  94        AMRB =AMOM
  95        GAMRB=GAMOM
  96        AMRA =AMRO
  97        GAMRA=GAMRO
  98       ELSEIF(MNUM.EQ.101) THEN
  99        PROB1=.35
 100        PROB2=.35
 101        AMRX =1.2
 102        GAMRX=.46
 103        AMRB =AMOM
 104        GAMRB=GAMOM
 105        AMRA =AMOM
 106        GAMRA=GAMOM
 107       ELSEIF(MNUM.EQ.102) THEN
 108        PROB1=0.0
 109        PROB2=0.0
 110        AMRX =1.4
 111        GAMRX=.6
 112        AMRB =AMOM
 113        GAMRB=GAMOM
 114        AMRA =AMOM
 115        GAMRA=GAMOM
 116       ELSE
 117        PROB1=0.0
 118        PROB2=0.0
 119        AMRX =AMA1
 120        GAMRX=GAMA1
 121        AMRA =AMRO
 122        GAMRA=GAMRO
 123        AMRB =AMRO
 124        GAMRB=GAMRO
 125       ENDIF
 126 C
 127       IF    (RR.LE.PROB1) THEN
 128        ICHAN=1
 129       ELSEIF(RR.LE.(PROB1+PROB2)) THEN
 130        ICHAN=2
 131         AX   =AMRA
 132         GX   =GAMRA
 133         AMRA =AMRB
 134         GAMRA=GAMRB
 135         AMRB =AX
 136         GAMRB=GX
 137         PX   =PROB1
 138         PROB1=PROB2
 139         PROB2=PX
 140       ELSE
 141        ICHAN=3
 142       ENDIF
 143 C
 144       PROB3=1.0-PROB1-PROB2
 145       END
 146       SUBROUTINE INITDK
 147 * ----------------------------------------------------------------------
 148 *     INITIALISATION OF TAU DECAY PARAMETERS  and routines
 149 *
 150 *     called by : KORALZ
 151 * ----------------------------------------------------------------------
 152
 153       COMMON / DECPAR / GFERMI,GV,GA,CCABIB,SCABIB,GAMEL
 154       REAL*4            GFERMI,GV,GA,CCABIB,SCABIB,GAMEL
 155       COMMON / PARMAS / AMTAU,AMNUTA,AMEL,AMNUE,AMMU,AMNUMU
 156      *                 ,AMPIZ,AMPI,AMRO,GAMRO,AMA1,GAMA1
 157      *                 ,AMK,AMKZ,AMKST,GAMKST
 158 *
 159       REAL*4            AMTAU,AMNUTA,AMEL,AMNUE,AMMU,AMNUMU
 160      *                 ,AMPIZ,AMPI,AMRO,GAMRO,AMA1,GAMA1
 161      *                 ,AMK,AMKZ,AMKST,GAMKST
 162       COMMON / TAUBRA / GAMPRT(30),JLIST(30),NCHAN
 163       COMMON / TAUKLE / BRA1,BRK0,BRK0B,BRKS
 164       REAL*4            BRA1,BRK0,BRK0B,BRKS
 165
 166
 167
 168
 169
 170
 171       PARAMETER (NMODE=15,NM1=0,NM2=1,NM3=8,NM4=2,NM5=1,NM6=3)
 172       COMMON / TAUDCD /IDFFIN(9,NMODE),MULPIK(NMODE)
 173      &                ,NAMES
 174       CHARACTER NAMES(NMODE)*31
 175
 176       CHARACTER OLDNAMES(7)*31
 177       CHARACTER*80 bxINIT
 178       PARAMETER (
 179      $  bxINIT ='(1x,1h*,g17.8,            16x, a31,a4,a4, 1x,1h*)'
 180      $ )
 181       REAL*4 PI
 182 *
 183 *
 184 * LIST OF BRANCHING RATIOS
 185 CAM normalised to e nu nutau channel
 186 CAM                  enu   munu   pinu  rhonu   A1nu   Knu    K*nu   pi
 187 CAM   DATA JLIST  /    1,     2,     3,     4,     5,     6,     7,
 188 *AM   DATA GAMPRT /1.000,0.9730,0.6054,1.2432,0.8432,0.0432,O.O811,0.616
 189 *AM
 190 *AM  multipion decays
 191 *
 192 *    conventions of particles names
 193 *                 K-,P-,K+,  K0,P-,KB,  K-,P0,K0
 194 *                  3, 1,-3  , 4, 1,-4  , 3, 2, 4  ,
 195 *                 P0,P0,K-,  K-,P-,P+,  P-,KB,P0
 196 *                  2, 2, 3  , 3, 1,-1  , 1,-4, 2  ,
 197 *                 ET,P-,P0   P-,P0,GM
 198 *                  9, 1, 2  , 1, 2, 8
 199 *
 200 C
 201       DIMENSION NOPIK(6,NMODE),NPIK(NMODE)
 202 *AM   outgoing multiplicity and flavors of multi-pion /multi-K modes
 203       DATA   NPIK  /                4,                    4,
 204      1                              5,                    5,
 205      2                              6,                    6,
 206      3                              3,                    3,
 207      4                              3,                    3,
 208      5                              3,                    3,
 209      6                              3,                    3,
 210      7                              2                         /
 211       DATA  NOPIK / -1,-1, 1, 2, 0, 0,     2, 2, 2,-1, 0, 0,
 212      1              -1,-1, 1, 2, 2, 0,    -1,-1,-1, 1, 1, 0,
 213      2              -1,-1,-1, 1, 1, 2,    -1,-1, 1, 2, 2, 2,
 214      3              -3,-1, 3, 0, 0, 0,    -4,-1, 4, 0, 0, 0,
 215      4              -3, 2,-4, 0, 0, 0,     2, 2,-3, 0, 0, 0,
 216      5              -3,-1, 1, 0, 0, 0,    -1, 4, 2, 0, 0, 0,
 217      6               9,-1, 2, 0, 0, 0,    -1, 2, 8, 0, 0, 0,
 218 C AJWMOD fix sign bug, 2/22/99
 219      7              -3,-4, 0, 0, 0, 0                         /
 220 * LIST OF BRANCHING RATIOS
 221       NCHAN = NMODE + 7
 222       DO 1 I = 1,30
 223       IF (I.LE.NCHAN) THEN
 224         JLIST(I) = I
 225         IF(I.EQ. 1) GAMPRT(I) =0.1800
 226         IF(I.EQ. 2) GAMPRT(I) =0.1751
 227         IF(I.EQ. 3) GAMPRT(I) =0.1110
 228         IF(I.EQ. 4) GAMPRT(I) =0.2515
 229         IF(I.EQ. 5) GAMPRT(I) =0.1790
 230         IF(I.EQ. 6) GAMPRT(I) =0.0071
 231         IF(I.EQ. 7) GAMPRT(I) =0.0134
 232         IF(I.EQ. 8) GAMPRT(I) =0.0450
 233         IF(I.EQ. 9) GAMPRT(I) =0.0100
 234         IF(I.EQ.10) GAMPRT(I) =0.0009
 235         IF(I.EQ.11) GAMPRT(I) =0.0004
 236         IF(I.EQ.12) GAMPRT(I) =0.0003
 237         IF(I.EQ.13) GAMPRT(I) =0.0005
 238         IF(I.EQ.14) GAMPRT(I) =0.0015
 239         IF(I.EQ.15) GAMPRT(I) =0.0015
 240         IF(I.EQ.16) GAMPRT(I) =0.0015
 241         IF(I.EQ.17) GAMPRT(I) =0.0005
 242         IF(I.EQ.18) GAMPRT(I) =0.0050
 243         IF(I.EQ.19) GAMPRT(I) =0.0055
 244         IF(I.EQ.20) GAMPRT(I) =0.0017
 245         IF(I.EQ.21) GAMPRT(I) =0.0013
 246         IF(I.EQ.22) GAMPRT(I) =0.0010
 247         IF(I.EQ. 1) OLDNAMES(I)='  TAU-  -->   E-               '
 248         IF(I.EQ. 2) OLDNAMES(I)='  TAU-  -->  MU-               '
 249         IF(I.EQ. 3) OLDNAMES(I)='  TAU-  -->  PI-               '
 250         IF(I.EQ. 4) OLDNAMES(I)='  TAU-  -->  PI-, PI0          '
 251         IF(I.EQ. 5) OLDNAMES(I)='  TAU-  -->  A1- (two subch)   '
 252         IF(I.EQ. 6) OLDNAMES(I)='  TAU-  -->   K-               '
 253         IF(I.EQ. 7) OLDNAMES(I)='  TAU-  -->  K*- (two subch)   '
 254         IF(I.EQ. 8) NAMES(I-7)='  TAU-  --> 2PI-,  PI0,  PI+   '
 255         IF(I.EQ. 9) NAMES(I-7)='  TAU-  --> 3PI0,        PI-   '
 256         IF(I.EQ.10) NAMES(I-7)='  TAU-  --> 2PI-,  PI+, 2PI0   '
 257         IF(I.EQ.11) NAMES(I-7)='  TAU-  --> 3PI-, 2PI+,        '
 258         IF(I.EQ.12) NAMES(I-7)='  TAU-  --> 3PI-, 2PI+,  PI0   '
 259         IF(I.EQ.13) NAMES(I-7)='  TAU-  --> 2PI-,  PI+, 3PI0   '
 260         IF(I.EQ.14) NAMES(I-7)='  TAU-  -->  K-, PI-,  K+      '
 261         IF(I.EQ.15) NAMES(I-7)='  TAU-  -->  K0, PI-, K0B      '
 262         IF(I.EQ.16) NAMES(I-7)='  TAU-  -->  K-,  K0, PI0      '
 263         IF(I.EQ.17) NAMES(I-7)='  TAU-  --> PI0  PI0   K-      '
 264         IF(I.EQ.18) NAMES(I-7)='  TAU-  -->  K-  PI-  PI+      '
 265         IF(I.EQ.19) NAMES(I-7)='  TAU-  --> PI-  K0B  PI0      '
 266         IF(I.EQ.20) NAMES(I-7)='  TAU-  --> ETA  PI-  PI0      '
 267         IF(I.EQ.21) NAMES(I-7)='  TAU-  --> PI-  PI0  GAM      '
 268         IF(I.EQ.22) NAMES(I-7)='  TAU-  -->  K-  K0            '
 269       ELSE
 270         JLIST(I) = 0
 271         GAMPRT(I) = 0.
 272       ENDIF
 273    1  CONTINUE
 274       DO I=1,NMODE
 275         MULPIK(I)=NPIK(I)
 276         DO J=1,MULPIK(I)
 277          IDFFIN(J,I)=NOPIK(J,I)
 278         ENDDO
 279       ENDDO
 280 *
 281 *
 282 * --- COEFFICIENTS TO FIX RATIO OF:
 283 * --- A1 3CHARGED/ A1 1CHARGED 2 NEUTRALS MATRIX ELEMENTS (MASLESS LIM.)
 284 * --- PROBABILITY OF K0 TO BE KS
 285 * --- PROBABILITY OF K0B TO BE KS
 286 * --- RATIO OF COEFFICIENTS FOR K*--> K0 PI-
 287 * --- ALL COEFFICENTS SHOULD BE IN THE RANGE (0.0,1.0)
 288 * --- THEY MEANING IS PROBABILITY OF THE FIRST CHOICE ONLY IF ONE
 289 * --- NEGLECTS MASS-PHASE SPACE EFFECTS
 290       BRA1=0.5
 291       BRK0=0.5
 292       BRK0B=0.5
 293       BRKS=0.6667
 294 *
 295
 296       GFERMI = 1.16637E-5
 297       CCABIB = 0.975
 298       GV     = 1.0
 299       GA     =-1.0
 300
 301
 302
 303 * ZW 13.04.89 HERE WAS AN ERROR
 304       SCABIB = SQRT(1.-CCABIB**2)
 305       PI =4.*ATAN(1.)
 306       GAMEL  = GFERMI**2*AMTAU**5/(192*PI**3)
 307 *
 308 *      CALL DEXAY(-1,pol1)
 309 *
 310       RETURN
 311       END
 312       FUNCTION DCDMAS(IDENT)
 313       COMMON / PARMAS / AMTAU,AMNUTA,AMEL,AMNUE,AMMU,AMNUMU
 314      *                 ,AMPIZ,AMPI,AMRO,GAMRO,AMA1,GAMA1
 315      *                 ,AMK,AMKZ,AMKST,GAMKST
 316 *
 317       REAL*4            AMTAU,AMNUTA,AMEL,AMNUE,AMMU,AMNUMU
 318      *                 ,AMPIZ,AMPI,AMRO,GAMRO,AMA1,GAMA1
 319      *                 ,AMK,AMKZ,AMKST,GAMKST
 320       IF      (IDENT.EQ. 1) THEN
 321         APKMAS=AMPI
 322       ELSEIF  (IDENT.EQ.-1) THEN
 323         APKMAS=AMPI
 324       ELSEIF  (IDENT.EQ. 2) THEN
 325         APKMAS=AMPIZ
 326       ELSEIF  (IDENT.EQ.-2) THEN
 327         APKMAS=AMPIZ
 328       ELSEIF  (IDENT.EQ. 3) THEN
 329         APKMAS=AMK
 330       ELSEIF  (IDENT.EQ.-3) THEN
 331         APKMAS=AMK
 332       ELSEIF  (IDENT.EQ. 4) THEN
 333         APKMAS=AMKZ
 334       ELSEIF  (IDENT.EQ.-4) THEN
 335         APKMAS=AMKZ
 336       ELSEIF  (IDENT.EQ. 8) THEN
 337         APKMAS=0.0001
 338       ELSEIF  (IDENT.EQ.-8) THEN
 339         APKMAS=0.0001
 340       ELSEIF  (IDENT.EQ. 9) THEN
 341         APKMAS=0.5488
 342       ELSEIF  (IDENT.EQ.-9) THEN
 343         APKMAS=0.5488
 344       ELSE
 345         PRINT *, 'STOP IN APKMAS, WRONG IDENT=',IDENT
 346         STOP
 347       ENDIF
 348       DCDMAS=APKMAS
 349       END
 350       FUNCTION LUNPIK(ID,ISGN)
 351       COMMON / TAUKLE / BRA1,BRK0,BRK0B,BRKS
 352       REAL*4            BRA1,BRK0,BRK0B,BRKS
 353       REAL*4 XIO(1)
 354       IDENT=ID*ISGN
 355       IF      (IDENT.EQ. 1) THEN
 356         IPKDEF=-211
 357       ELSEIF  (IDENT.EQ.-1) THEN
 358         IPKDEF= 211
 359       ELSEIF  (IDENT.EQ. 2) THEN
 360         IPKDEF=111
 361       ELSEIF  (IDENT.EQ.-2) THEN
 362         IPKDEF=111
 363       ELSEIF  (IDENT.EQ. 3) THEN
 364         IPKDEF=-321
 365       ELSEIF  (IDENT.EQ.-3) THEN
 366         IPKDEF= 321
 367       ELSEIF  (IDENT.EQ. 4) THEN
 368 *
 369 * K0 --> K0_LONG (IS 130) / K0_SHORT (IS 310) = 1/1
 370         CALL RANMAR(XIO,1)
 371         IF (XIO(1).GT.BRK0) THEN
 372           IPKDEF= 130
 373         ELSE
 374           IPKDEF= 310
 375         ENDIF
 376       ELSEIF  (IDENT.EQ.-4) THEN
 377 *
 378 * K0B--> K0_LONG (IS 130) / K0_SHORT (IS 310) = 1/1
 379         CALL RANMAR(XIO,1)
 380         IF (XIO(1).GT.BRK0B) THEN
 381           IPKDEF= 130
 382         ELSE
 383           IPKDEF= 310
 384         ENDIF
 385       ELSEIF  (IDENT.EQ. 8) THEN
 386         IPKDEF= 22
 387       ELSEIF  (IDENT.EQ.-8) THEN
 388         IPKDEF= 22
 389       ELSEIF  (IDENT.EQ. 9) THEN
 390         IPKDEF= 221
 391       ELSEIF  (IDENT.EQ.-9) THEN
 392         IPKDEF= 221
 393       ELSE
 394         PRINT *, 'STOP IN IPKDEF, WRONG IDENT=',IDENT
 395         STOP
 396       ENDIF
 397       LUNPIK=IPKDEF
 398       END
 399
 400
 401
 402       SUBROUTINE TAURDF(KTO)
 403 C THIS ROUTINE CAN BE CALLED BEFORE ANY TAU+ OR TAU- EVENT IS GENERATED
 404 C IT CAN BE USED TO GENERATE TAU+ AND TAU- SAMPLES OF DIFFERENT
 405 C CONTENTS
 406       COMMON / TAUKLE / BRA1,BRK0,BRK0B,BRKS
 407       REAL*4            BRA1,BRK0,BRK0B,BRKS
 408       COMMON / TAUBRA / GAMPRT(30),JLIST(30),NCHAN
 409
 410 C Subroutine TAURDF is disabled
 411       RETURN
 412
 413
 414       IF (KTO.EQ.1) THEN
 415 C     ==================
 416 C AJWMOD: Set the BRs for (A1+ -> rho+ pi0) and (K*+ -> K0 pi+)
 417       BRA1 = PKORB(4,1)
 418       BRKS = PKORB(4,3)
 419       BRK0  = PKORB(4,5)
 420       BRK0B  = PKORB(4,6)
 421       ELSE
 422 C     ====
 423 C AJWMOD: Set the BRs for (A1+ -> rho+ pi0) and (K*+ -> K0 pi+)
 424       BRA1 = PKORB(4,2)
 425       BRKS = PKORB(4,4)
 426       BRK0  = PKORB(4,5)
 427       BRK0B  = PKORB(4,6)
 428       ENDIF
 429 C     =====
 430       END
 431
 432       SUBROUTINE INIPHY(XK00)
 433 * ----------------------------------------------------------------------
 434 *     INITIALISATION OF PARAMETERS
 435 *     USED IN QED and/or GSW ROUTINES
 436 * ----------------------------------------------------------------------
 437       COMMON / QEDPRM /ALFINV,ALFPI,XK0
 438       REAL*8           ALFINV,ALFPI,XK0
 439       REAL*8 PI8,XK00
 440 *
 441       PI8    = 4.D0*DATAN(1.D0)
 442       ALFINV = 137.03604D0
 443       ALFPI  = 1D0/(ALFINV*PI8)
 444       XK0=XK00
 445       END
 446
 447       SUBROUTINE INIMAS
 448 C ----------------------------------------------------------------------
 449 C     INITIALISATION OF MASSES
 450 C
 451 C     called by : KORALZ
 452 C ----------------------------------------------------------------------
 453       COMMON / PARMAS / AMTAU,AMNUTA,AMEL,AMNUE,AMMU,AMNUMU
 454      *                 ,AMPIZ,AMPI,AMRO,GAMRO,AMA1,GAMA1
 455      *                 ,AMK,AMKZ,AMKST,GAMKST
 456 *
 457       REAL*4            AMTAU,AMNUTA,AMEL,AMNUE,AMMU,AMNUMU
 458      *                 ,AMPIZ,AMPI,AMRO,GAMRO,AMA1,GAMA1
 459      *                 ,AMK,AMKZ,AMKST,GAMKST
 460 C
 461 C IN-COMING / OUT-GOING  FERMION MASSES
 462       AMTAU  = 1.7842
 463 C --- let us update tau mass ...
 464       AMTAU  = 1.777
 465       AMNUTA = 0.010
 466       AMEL   = 0.0005111
 467       AMNUE  = 0.0
 468       AMMU   = 0.105659
 469       AMNUMU = 0.0
 470 *
 471 * MASSES USED IN TAU DECAYS
 472       AMPIZ  = 0.134964
 473       AMPI   = 0.139568
 474       AMRO   = 0.773
 475       GAMRO  = 0.145
 476 *C    GAMRO  = 0.666
 477       AMA1   = 1.251
 478       GAMA1  = 0.599
 479       AMK    = 0.493667
 480       AMKZ   = 0.49772
 481       AMKST  = 0.8921
 482       GAMKST = 0.0513
 483 C
 484 C
 485 C IN-COMING / OUT-GOING  FERMION MASSES
 486 !!      AMNUTA = PKORB(1,2)
 487 !!      AMNUE  = PKORB(1,4)
 488 !!      AMNUMU = PKORB(1,6)
 489 C
 490 C MASSES USED IN TAU DECAYS  Cleo settings
 491 !!      AMPIZ  = PKORB(1,7)
 492 !!      AMPI   = PKORB(1,8)
 493 !!      AMRO   = PKORB(1,9)
 494 !!      GAMRO  = PKORB(2,9)
 495       AMA1   = 1.275   !! PKORB(1,10)
 496       GAMA1  = 0.615   !! PKORB(2,10)
 497 !!      AMK    = PKORB(1,11)
 498 !!      AMKZ   = PKORB(1,12)
 499 !!      AMKST  = PKORB(1,13)
 500 !!      GAMKST = PKORB(2,13)
 501 C
 502
 503       RETURN
 504       END
 505       SUBROUTINE ANGULU(PD1,PD2,Q1,Q2,COSTHE)
 506       REAL*8 PD1(4),PD2(4),Q1(4),Q2(4),COSTHE,P(4),QQ(4),QT(4)
 507 C take effective beam which is less massive, it should be irrelevant
 508 C but in case HEPEVT is particulary dirty may help.
 509 C this routine calculate reduced system transver and cosine of scattering
 510 C angle.
 511
 512       XM1=ABS(PD1(4)**2-PD1(3)**2-PD1(2)**2-PD1(1)**2)
 513       XM2=ABS(PD2(4)**2-PD2(3)**2-PD2(2)**2-PD2(1)**2)
 514       IF (XM1.LT.XM2) THEN
 515         SIGN=1D0
 516         DO K=1,4
 517           P(K)=PD1(K)
 518         ENDDO
 519       ELSE
 520         SIGN=-1D0
 521         DO K=1,4
 522           P(K)=PD2(K)
 523         ENDDO
 524       ENDIF
 525 C calculate space like part of P (in Z restframe)
 526       DO K=1,4
 527        QQ(K)=Q1(k)+Q2(K)
 528        QT(K)=Q1(K)-Q2(K)
 529       ENDDO
 530
 531        XMQQ=SQRT(QQ(4)**2-QQ(3)**2-QQ(2)**2-QQ(1)**2)
 532
 533        QTXQQ=QT(4)*QQ(4)-QT(3)*QQ(3)-QT(2)*QQ(2)-QT(1)*QQ(1)
 534       DO K=1,4
 535        QT(K)=QT(K)-QQ(K)*QTXQQ/XMQQ**2
 536       ENDDO
 537
 538        PXQQ=P(4)*QQ(4)-P(3)*QQ(3)-P(2)*QQ(2)-P(1)*QQ(1)
 539       DO K=1,4
 540        P(K)=P(K)-QQ(K)*PXQQ/XMQQ**2
 541       ENDDO
 542 C calculate costhe
 543        PXP  =SQRT(p(1)**2+p(2)**2+p(3)**2-p(4)**2)
 544        QTXQT=SQRT(QT(3)**2+QT(2)**2+QT(1)**2-QT(4)**2)
 545        PXQT =P(3)*QT(3)+P(2)*QT(2)+P(1)*QT(1)-P(4)*QT(4)
 546        COSTHE=PXQT/PXP/QTXQT
 547        COSTHE=COSTHE*SIGN
 548       END
 549
 550       FUNCTION PLZAP0(IDE,IDF,SVAR,COSTH0)
 551 C this function calculates probability for the helicity +1 +1 configuration
 552 C of taus for given Z/gamma transfer and COSTH0 cosine of scattering angle
 553       REAL*8 PLZAP0,SVAR,COSTHE,COSTH0,T_BORN
 554
 555       COSTHE=COSTH0
 556 C >>>>>      IF (IDE*IDF.LT.0) COSTHE=-COSTH0 ! this is probably not needed ID
 557 C >>>>>      of first beam is used by T_GIVIZ0 including sign
 558
 559       IF (IDF.GT.0) THEN
 560         CALL INITWK(IDE,IDF,SVAR)
 561       ELSE
 562         CALL INITWK(-IDE,-IDF,SVAR)
 563       ENDIF
 564       PLZAP0=T_BORN(0,SVAR,COSTHE,1D0,1D0)
 565      $  /(T_BORN(0,SVAR,COSTHE,1D0,1D0)+T_BORN(0,SVAR,COSTHE,-1D0,-1D0))
 566
 567
 568 C      write(*,*) 'svar=  ',  svar
 569 C      write(*,*) 'COSTHE=',  COSTHE
 570 C      write(*,*) ide,'  ',idf
 571 C      write(*,*) 'PLZAP0=', PLZAP0
 572 C      COSTHE=0.999
 573 C      write(*,*) 'TBORN+=', T_BORN(0,SVAR,COSTHE,1D0,1D0)
 574
 575 C      COSTHE=-0.999
 576 C      write(*,*) 'TBORN-=', T_BORN(0,SVAR,COSTHE,-1D0,-1D0)
 577 C 100  format (A,E8.4)
 578
 579 !      PLZAP0=0.5
 580       END
 581       FUNCTION T_BORN(MODE,SVAR,COSTHE,TA,TB)
 582 C ----------------------------------------------------------------------
 583 C THIS ROUTINE PROVIDES BORN CROSS SECTION. IT HAS THE SAME
 584 C STRUCTURE AS FUNTIS AND FUNTIH, THUS CAN BE USED AS SIMPLER
 585 C EXAMPLE OF THE METHOD APPLIED THERE
 586 C INPUT PARAMETERS ARE: SVAR    -- transfer
 587 C                       COSTHE  -- cosine of angle between tau+ and 1st beam
 588 C                       TA,TB   -- helicity states of tau+ tau-
 589 C
 590 C     called by : BORNY, BORAS, BORNV, WAGA, WEIGHT
 591 C ----------------------------------------------------------------------
 592       IMPLICIT REAL*8(A-H,O-Z)
 593       COMMON / T_BEAMPM / ENE ,AMIN,AMFIN,IDE,IDF
 594       REAL*8              ENE ,AMIN,AMFIN
 595       COMMON / T_GAUSPM /SS,POLN,T3E,QE,T3F,QF
 596      &                  ,XUPGI   ,XUPZI   ,XUPGF   ,XUPZF
 597      &                  ,NDIAG0,NDIAGA,KEYA,KEYZ
 598      &                  ,ITCE,JTCE,ITCF,JTCF,KOLOR
 599       REAL*8             SS,POLN,T3E,QE,T3F,QF
 600      &                  ,XUPGI(2),XUPZI(2),XUPGF(2),XUPZF(2)
 601       REAL*8            SEPS1,SEPS2
 602 C=====================================================================
 603       COMMON / T_GSWPRM /SWSQ,AMW,AMZ,AMH,AMTOP,GAMMZ
 604       REAL*8             SWSQ,AMW,AMZ,AMH,AMTOP,GAMMZ
 605 C     SWSQ        = sin2 (theta Weinberg)
 606 C     AMW,AMZ     = W & Z boson masses respectively
 607 C     AMH         = the Higgs mass
 608 C     AMTOP       = the top mass
 609 C     GAMMZ       = Z0 width
 610       COMPLEX*16 ABORN(2,2),APHOT(2,2),AZETT(2,2)
 611       COMPLEX*16 XUPZFP(2),XUPZIP(2)
 612       COMPLEX*16 ABORNM(2,2),APHOTM(2,2),AZETTM(2,2)
 613       COMPLEX*16 PROPA,PROPZ
 614       COMPLEX*16 XR,XI
 615       COMPLEX*16 XUPF,XUPI
 616       COMPLEX*16 XTHING
 617       DATA XI/(0.D0,1.D0)/,XR/(1.D0,0.D0)/
 618       DATA MODE0 /-5/
 619       DATA IDE0 /-55/
 620       DATA SVAR0,COST0 /-5.D0,-6.D0/
 621       DATA PI /3.141592653589793238462643D0/
 622       DATA SEPS1,SEPS2 /0D0,0D0/
 623 C
 624 C MEMORIZATION =========================================================
 625       IF ( MODE.NE.MODE0.OR.SVAR.NE.SVAR0.OR.COSTHE.NE.COST0
 626      $    .OR.IDE0.NE.IDE)THEN
 627 C
 628         KEYGSW=1
 629 C ** PROPAGATORS
 630         IDE0=IDE
 631         MODE0=MODE
 632         SVAR0=SVAR
 633         COST0=COSTHE
 634         SINTHE=SQRT(1.D0-COSTHE**2)
 635         BETA=SQRT(MAX(0D0,1D0-4D0*AMFIN**2/SVAR))
 636 C I MULTIPLY AXIAL COUPLING BY BETA FACTOR.
 637         XUPZFP(1)=0.5D0*(XUPZF(1)+XUPZF(2))+0.5*BETA*(XUPZF(1)-XUPZF(2))
 638         XUPZFP(2)=0.5D0*(XUPZF(1)+XUPZF(2))-0.5*BETA*(XUPZF(1)-XUPZF(2))
 639         XUPZIP(1)=0.5D0*(XUPZI(1)+XUPZI(2))+0.5*(XUPZI(1)-XUPZI(2))
 640         XUPZIP(2)=0.5D0*(XUPZI(1)+XUPZI(2))-0.5*(XUPZI(1)-XUPZI(2))
 641 C FINAL STATE VECTOR COUPLING
 642         XUPF     =0.5D0*(XUPZF(1)+XUPZF(2))
 643         XUPI     =0.5D0*(XUPZI(1)+XUPZI(2))
 644         XTHING   =0D0
 645
 646         PROPA =1D0/SVAR
 647         PROPZ =1D0/DCMPLX(SVAR-AMZ**2,SVAR/AMZ*GAMMZ)
 648         IF (KEYGSW.EQ.0) PROPZ=0.D0
 649         DO 50 I=1,2
 650          DO 50 J=1,2
 651           REGULA= (3-2*I)*(3-2*J) + COSTHE
 652           REGULM=-(3-2*I)*(3-2*J) * SINTHE *2.D0*AMFIN/SQRT(SVAR)
 653           APHOT(I,J)=PROPA*(XUPGI(I)*XUPGF(J)*REGULA)
 654           AZETT(I,J)=PROPZ*(XUPZIP(I)*XUPZFP(J)+XTHING)*REGULA
 655           ABORN(I,J)=APHOT(I,J)+AZETT(I,J)
 656           APHOTM(I,J)=PROPA*DCMPLX(0D0,1D0)*XUPGI(I)*XUPGF(J)*REGULM
 657           AZETTM(I,J)=PROPZ*DCMPLX(0D0,1D0)*(XUPZIP(I)*XUPF+XTHING)*REGULM
 658           ABORNM(I,J)=APHOTM(I,J)+AZETTM(I,J)
 659    50   CONTINUE
 660       ENDIF
 661 C
 662 C******************
 663 C* IN CALCULATING CROSS SECTION ONLY DIAGONAL ELEMENTS
 664 C* OF THE SPIN DENSITY MATRICES ENTER (LONGITUD. POL. ONLY.)
 665 C* HELICITY CONSERVATION EXPLICITLY OBEYED
 666       POLAR1=  (SEPS1)
 667       POLAR2= (-SEPS2)
 668       BORN=0D0
 669       DO 150 I=1,2
 670        HELIC= 3-2*I
 671        DO 150 J=1,2
 672         HELIT=3-2*J
 673         FACTOR=KOLOR*(1D0+HELIC*POLAR1)*(1D0-HELIC*POLAR2)/4D0
 674         FACTOM=FACTOR*(1+HELIT*TA)*(1-HELIT*TB)
 675         FACTOR=FACTOR*(1+HELIT*TA)*(1+HELIT*TB)
 676
 677         BORN=BORN+CDABS(ABORN(I,J))**2*FACTOR
 678 C      MASS TERM IN BORN
 679         IF (MODE.GE.1) THEN
 680          BORN=BORN+CDABS(ABORNM(I,J))**2*FACTOM
 681         ENDIF
 682
 683   150 CONTINUE
 684 C************
 685       FUNT=BORN
 686       IF(FUNT.LT.0.D0)  FUNT=BORN
 687
 688 C
 689       IF (SVAR.GT.4D0*AMFIN**2) THEN
 690 C PHASE SPACE THRESHOLD FACTOR
 691         THRESH=SQRT(1-4D0*AMFIN**2/SVAR)
 692         T_BORN= FUNT*SVAR**2*THRESH
 693       ELSE
 694         THRESH=0.D0
 695         T_BORN=0.D0
 696       ENDIF
 697 C ZW HERE WAS AN ERROR 19. 05. 1989
 698 !      write(*,*) 'KKKK ',PROPA,PROPZ,XUPGI,XUPGF,XUPZI,XUPZF
 699 !      write(*,*) 'KKKK X',svar,costhe,TA,TB,T_BORN
 700       END
 701
 702       SUBROUTINE INITWK(IDEX,IDFX,SVAR)
 703 ! initialization routine coupling masses etc.
 704       IMPLICIT REAL*8 (A-H,O-Z)
 705       COMMON / T_BEAMPM / ENE ,AMIN,AMFIN,IDE,IDF
 706       REAL*8              ENE ,AMIN,AMFIN
 707       COMMON / T_GAUSPM /SS,POLN,T3E,QE,T3F,QF
 708      &                  ,XUPGI   ,XUPZI   ,XUPGF   ,XUPZF
 709      &                  ,NDIAG0,NDIAGA,KEYA,KEYZ
 710      &                  ,ITCE,JTCE,ITCF,JTCF,KOLOR
 711       REAL*8             SS,POLN,T3E,QE,T3F,QF
 712      &                  ,XUPGI(2),XUPZI(2),XUPGF(2),XUPZF(2)
 713       COMMON / T_GSWPRM /SWSQ,AMW,AMZ,AMH,AMTOP,GAMMZ
 714       REAL*8             SWSQ,AMW,AMZ,AMH,AMTOP,GAMMZ
 715 C     SWSQ        = sin2 (theta Weinberg)
 716 C     AMW,AMZ     = W & Z boson masses respectively
 717 C     AMH         = the Higgs mass
 718 C     AMTOP       = the top mass
 719 C     GAMMZ       = Z0 width
 720 C
 721       ENE=SQRT(SVAR)/2
 722       AMIN=0.511D-3
 723       SWSQ=0.23147
 724       AMZ=91.1882
 725       GAMMZ=2.4952
 726       IF     (IDFX.EQ. 15) then
 727         IDF=2  ! denotes tau +2 tau-
 728         AMFIN=1.77703 !this mass is irrelevant if small, used in ME only
 729       ELSEIF (IDFX.EQ.-15) then
 730         IDF=-2  ! denotes tau -2 tau-
 731         AMFIN=1.77703 !this mass is irrelevant if small, used in ME only
 732       ELSE
 733         WRITE(*,*) 'INITWK: WRONG IDFX'
 734         STOP
 735       ENDIF
 736
 737       IF     (IDEX.EQ. 11) then      !electron
 738         IDE= 2
 739         AMIN=0.511D-3
 740       ELSEIF (IDEX.EQ.-11) then      !positron
 741         IDE=-2
 742         AMIN=0.511D-3
 743       ELSEIF (IDEX.EQ. 13) then      !mu+
 744         IDE= 2
 745         AMIN=0.105659
 746       ELSEIF (IDEX.EQ.-13) then      !mu-
 747         IDE=-2
 748         AMIN=0.105659
 749       ELSEIF (IDEX.EQ.  1) then      !d
 750         IDE= 4
 751         AMIN=0.05
 752       ELSEIF (IDEX.EQ.- 1) then      !d~
 753         IDE=-4
 754         AMIN=0.05
 755       ELSEIF (IDEX.EQ.  2) then      !u
 756         IDE= 3
 757         AMIN=0.02
 758       ELSEIF (IDEX.EQ.- 2) then      !u~
 759         IDE=-3
 760         AMIN=0.02
 761       ELSEIF (IDEX.EQ.  3) then      !s
 762         IDE= 4
 763         AMIN=0.3
 764       ELSEIF (IDEX.EQ.- 3) then      !s~
 765         IDE=-4
 766         AMIN=0.3
 767       ELSEIF (IDEX.EQ.  4) then      !c
 768         IDE= 3
 769         AMIN=1.3
 770       ELSEIF (IDEX.EQ.- 4) then      !c~
 771         IDE=-3
 772         AMIN=1.3
 773       ELSEIF (IDEX.EQ.  5) then      !b
 774         IDE= 4
 775         AMIN=4.5
 776       ELSEIF (IDEX.EQ.- 5) then      !b~
 777         IDE=-4
 778         AMIN=4.5
 779       ELSEIF (IDEX.EQ.  12) then     !nu_e
 780         IDE= 1
 781         AMIN=0.1D-3
 782       ELSEIF (IDEX.EQ.- 12) then     !nu_e~
 783         IDE=-1
 784         AMIN=0.1D-3
 785       ELSEIF (IDEX.EQ.  14) then     !nu_mu
 786         IDE= 1
 787         AMIN=0.1D-3
 788       ELSEIF (IDEX.EQ.- 14) then     !nu_mu~
 789         IDE=-1
 790         AMIN=0.1D-3
 791       ELSEIF (IDEX.EQ.  16) then     !nu_tau
 792         IDE= 1
 793         AMIN=0.1D-3
 794       ELSEIF (IDEX.EQ.- 16) then     !nu_tau~
 795         IDE=-1
 796         AMIN=0.1D-3
 797
 798       ELSE
 799         WRITE(*,*) 'INITWK: WRONG IDEX'
 800         STOP
 801       ENDIF
 802
 803 C ----------------------------------------------------------------------
 804 C
 805 C     INITIALISATION OF COUPLING CONSTANTS AND FERMION-GAMMA / Z0 VERTEX
 806 C
 807 C     called by : KORALZ
 808 C ----------------------------------------------------------------------
 809       ITCE=IDE/IABS(IDE)
 810       JTCE=(1-ITCE)/2
 811       ITCF=IDF/IABS(IDF)
 812       JTCF=(1-ITCF)/2
 813       CALL T_GIVIZO( IDE, 1,AIZOR,QE,KDUMM)
 814       CALL T_GIVIZO( IDE,-1,AIZOL,QE,KDUMM)
 815       XUPGI(1)=QE
 816       XUPGI(2)=QE
 817       T3E    = AIZOL+AIZOR
 818       XUPZI(1)=(AIZOR-QE*SWSQ)/SQRT(SWSQ*(1-SWSQ))
 819       XUPZI(2)=(AIZOL-QE*SWSQ)/SQRT(SWSQ*(1-SWSQ))
 820       CALL T_GIVIZO( IDF, 1,AIZOR,QF,KOLOR)
 821       CALL T_GIVIZO( IDF,-1,AIZOL,QF,KOLOR)
 822       XUPGF(1)=QF
 823       XUPGF(2)=QF
 824       T3F    =  AIZOL+AIZOR
 825       XUPZF(1)=(AIZOR-QF*SWSQ)/SQRT(SWSQ*(1-SWSQ))
 826       XUPZF(2)=(AIZOL-QF*SWSQ)/SQRT(SWSQ*(1-SWSQ))
 827 C
 828       NDIAG0=2
 829       NDIAGA=11
 830       KEYA  = 1
 831       KEYZ  = 1
 832 C
 833 C
 834       RETURN
 835       END
 836
 837       SUBROUTINE T_GIVIZO(IDFERM,IHELIC,SIZO3,CHARGE,KOLOR)
 838 C ----------------------------------------------------------------------
 839 C PROVIDES ELECTRIC CHARGE AND WEAK IZOSPIN OF A FAMILY FERMION
 840 C IDFERM=1,2,3,4 DENOTES NEUTRINO, LEPTON, UP AND DOWN QUARK
 841 C NEGATIVE IDFERM=-1,-2,-3,-4, DENOTES ANTIPARTICLE
 842 C IHELIC=+1,-1 DENOTES RIGHT AND LEFT HANDEDNES ( CHIRALITY)
 843 C SIZO3 IS THIRD PROJECTION OF WEAK IZOSPIN (PLUS MINUS HALF)
 844 C AND CHARGE IS ELECTRIC CHARGE IN UNITS OF ELECTRON CHARGE
 845 C KOLOR IS A QCD COLOUR, 1 FOR LEPTON, 3 FOR QUARKS
 846 C
 847 C     called by : EVENTE, EVENTM, FUNTIH, .....
 848 C ----------------------------------------------------------------------
 849       IMPLICIT REAL*8(A-H,O-Z)
 850 C
 851       IF(IDFERM.EQ.0.OR.IABS(IDFERM).GT.4) GOTO 901
 852       IF(IABS(IHELIC).NE.1)                GOTO 901
 853       IH  =IHELIC
 854       IDTYPE =IABS(IDFERM)
 855       IC  =IDFERM/IDTYPE
 856       LEPQUA=INT(IDTYPE*0.4999999D0)
 857       IUPDOW=IDTYPE-2*LEPQUA-1
 858       CHARGE  =(-IUPDOW+2D0/3D0*LEPQUA)*IC
 859       SIZO3   =0.25D0*(IC-IH)*(1-2*IUPDOW)
 860       KOLOR=1+2*LEPQUA
 861 C** NOTE THAT CONVENTIONALY Z0 COUPLING IS
 862 C** XOUPZ=(SIZO3-CHARGE*SWSQ)/SQRT(SWSQ*(1-SWSQ))
 863       RETURN
 864  901  PRINT *,' STOP IN GIVIZO: WRONG PARAMS.'
 865       STOP
 866       END
 867       SUBROUTINE PHYFIX(NSTOP,NSTART)
 868       COMMON/LUJETS/N,K(4000,5),P(4000,5),V(4000,5)
 869       SAVE /LUJETS/
 870 C NSTOP NSTART : when PHYTIA history ends and event starts.
 871       NSTOP=0
 872       NSTART=1
 873       DO I=1, N
 874        IF(K(I,1).NE.21) THEN
 875            NSTOP = I-1
 876            NSTART= I
 877            GOTO 500
 878        ENDIF
 879       ENDDO
 880  500  CONTINUE
 881       END
 882
 883
 884       SUBROUTINE TAUPI0(PI0,K)
 885 C no initialization required. Must be called once after every:
 886 C   1)    CALL DEKAY(1+10,...)
 887 C   2)    CALL DEKAY(2+10,...)
 888 C   3)    CALL DEXAY(1,...)
 889 C   4)    CALL DEXAY(2,...)
 890 C subroutine to decay originating from TAUOLA's taus:
 891 C 1) etas (with CALL TAUETA(JAK))
 892 C 2) later pi0's from taus.
 893 C 3) extensions to other applications possible.
 894 C this routine belongs to >tauola universal interface<, but uses
 895 C routines from >tauola< utilities as well.  25.08.2005
 896 C this is the hepevt class in old style. No d_h_ class pre-name
 897
 898 C position of taus, must be defined by host program:
 899       COMMON /TAUPOS/ NP1,NP2
 900 c
 901       REAL  PHOT1(4),PHOT2(4)
 902       REAL*8  R,X(4),Y(4),PI0(4)
 903       INTEGER JEZELI(3),K
 904       DATA JEZELI /0,0,0/
 905       SAVE JEZELI
 906
 907 ! random 3 vector on the sphere, masless
 908         R=SQRT(PI0(4)**2-PI0(3)**2-PI0(2)**2-PI0(1)**2)/2D0
 909         CALL SPHERD(R,X)
 910         X(4)=R
 911         Y(4)=R
 912
 913         Y(1)=-X(1)
 914         Y(2)=-X(2)
 915         Y(3)=-X(3)
 916 ! boost to lab and to real*4
 917         CALL bostdq(-1,PI0,X,X)
 918         CALL bostdq(-1,PI0,Y,Y)
 919         DO L=1,4
 920          PHOT1(L)=X(L)
 921          PHOT2(L)=Y(L)
 922         ENDDO
 923 C to hepevt
 924         CALL FILHEP(0,1,22,K,K,0,0,PHOT1,0.0,.TRUE.)
 925         CALL FILHEP(0,1,22,K,K,0,0,PHOT2,0.0,.TRUE.)
 926
 927 C
 928       END
 929       SUBROUTINE TAUETA(PETA,K)
 930 C subroutine to decay etas's from taus.
 931 C this routine belongs to tauola universal interface, but uses
 932 C routines from tauola utilities. Just flat phase space, but 4 channels.
 933 C it is called at the beginning of SUBR. TAUPI0(JAK)
 934 C and as far as hepevt search it is basically the same as TAUPI0.  25.08.2005
 935 C this is the hepevt class in old style. No d_h_ class pre-name
 936
 937       COMMON / PARMAS / AMTAU,AMNUTA,AMEL,AMNUE,AMMU,AMNUMU
 938      *                 ,AMPIZ,AMPI,AMRO,GAMRO,AMA1,GAMA1
 939      *                 ,AMK,AMKZ,AMKST,GAMKST
 940 *
 941       REAL*4            AMTAU,AMNUTA,AMEL,AMNUE,AMMU,AMNUMU
 942      *                 ,AMPIZ,AMPI,AMRO,GAMRO,AMA1,GAMA1
 943      *                 ,AMK,AMKZ,AMKST,GAMKST
 944
 945 C position of taus, must be defined by host program:
 946 c
 947       REAL  RRR(1),BRSUM(3), RR(2)
 948       REAL  PHOT1(4),PHOT2(4),PHOT3(4)
 949       REAL*8    X(4),    Y(4),    Z(4)
 950       REAL                                YM1,YM2,YM3
 951       REAL*8  R,RU,PETA(4),XM1,XM2,XM3,XLAM
 952       REAL*8 a,b,c
 953       INTEGER K
 954       XLAM(a,b,c)=SQRT(ABS((a-b-c)**2-4.0*b*c))
 955 C position of decaying particle:
 956 C        DO L=1,4
 957 C          PETA(L)= phep(L,K)  ! eta 4 momentum
 958 C        ENDDO
 959 C       eta cumulated branching ratios:
 960         BRSUM(1)=0.389  ! gamma gamma
 961         BRSUM(2)=BRSUM(1)+0.319  ! 3 pi0
 962         BRSUM(3)=BRSUM(2)+0.237  ! pi+ pi- pi0 rest is thus pi+pi-gamma
 963         CALL RANMAR(RRR,1)
 964
 965         IF (RRR(1).LT.BRSUM(1)) THEN ! gamma gamma channel exactly like pi0
 966 ! random 3 vector on the sphere, masless
 967          R=SQRT(PETA(4)**2-PETA(3)**2-PETA(2)**2-PETA(1)**2)/2D0
 968          CALL SPHERD(R,X)
 969          X(4)=R
 970          Y(4)=R
 971
 972          Y(1)=-X(1)
 973          Y(2)=-X(2)
 974          Y(3)=-X(3)
 975 ! boost to lab and to real*4
 976          CALL bostdq(-1,PETA,X,X)
 977          CALL bostdq(-1,PETA,Y,Y)
 978          DO L=1,4
 979           PHOT1(L)=X(L)
 980           PHOT2(L)=Y(L)
 981          ENDDO
 982 C to hepevt
 983          CALL FILHEP(0,1,22,K,K,0,0,PHOT1,0.0,.TRUE.)
 984          CALL FILHEP(0,1,22,K,K,0,0,PHOT2,0.0,.TRUE.)
 985         ELSE ! 3 body channels
 986          IF(RRR(1).LT.BRSUM(2)) THEN  ! 3 pi0
 987           ID1= 111
 988           ID2= 111
 989           ID3= 111
 990           XM1=AMPIZ ! masses
 991           XM2=AMPIZ
 992           XM3=AMPIZ
 993          ELSEIF(RRR(1).LT.BRSUM(3)) THEN ! pi+ pi- pi0
 994           ID1= 211
 995           ID2=-211
 996           ID3= 111
 997           XM1=AMPI ! masses
 998           XM2=AMPI
 999           XM3=AMPIZ
1000          ELSE                            ! pi+ pi- gamma
1001           ID1= 211
1002           ID2=-211
1003           ID3=  22
1004           XM1=AMPI ! masses
1005           XM2=AMPI
1006           XM3=0.0
1007          ENDIF
1008  7       CONTINUE  ! we generate mass of the first pair:
1009           CALL RANMAR(RR,2)
1010           R=SQRT(PETA(4)**2-PETA(3)**2-PETA(2)**2-PETA(1)**2)
1011           AMIN=XM1+XM2
1012           AMAX=R-XM3
1013           AM2=SQRT(AMIN**2+RR(1)*(AMAX**2-AMIN**2))
1014 C         weight for flat phase space
1015           WT=XLAM(1D0*R**2,1D0*AM2**2,1D0*XM3**2)
1016      &      *XLAM(1D0*AM2**2,1D0*XM1**2,1D0*XM2**2)
1017      &           /R**2                    /AM2**2
1018          IF (RR(2).GT.WT) GOTO 7
1019
1020          RU=XLAM(1D0*AM2**2,1D0*XM1**2,1D0*XM2**2)/AM2/2  ! momenta of the
1021                                               ! first two products
1022                                               ! in the rest frame of that pair
1023          CALL SPHERD(RU,X)
1024          X(4)=SQRT(RU**2+XM1**2)
1025          Y(4)=SQRT(RU**2+XM2**2)
1026
1027          Y(1)=-X(1)
1028          Y(2)=-X(2)
1029          Y(3)=-X(3)
1030 C generate momentum of that pair in rest frame of eta:
1031          RU=XLAM(1D0*R**2,1D0*AM2**2,1D0*XM3**2)/R/2
1032          CALL SPHERD(RU,Z)
1033          Z(4)=SQRT(RU**2+AM2**2)
1034 C and boost first two decay products to rest frame of eta.
1035          CALL bostdq(-1,Z,X,X)
1036          CALL bostdq(-1,Z,Y,Y)
1037 C redefine Z(4) to 4-momentum of the last decay product:
1038          Z(1)=-Z(1)
1039          Z(2)=-Z(2)
1040          Z(3)=-Z(3)
1041          Z(4)=SQRT(RU**2+XM3**2)
1042 C boost all to lab and move to real*4; also masses
1043          CALL bostdq(-1,PETA,X,X)
1044          CALL bostdq(-1,PETA,Y,Y)
1045          CALL bostdq(-1,PETA,Z,Z)
1046          DO L=1,4
1047           PHOT1(L)=X(L)
1048           PHOT2(L)=Y(L)
1049           PHOT3(L)=Z(L)
1050          ENDDO
1051          YM1=XM1
1052          YM2=XM2
1053          YM3=XM3
1054 C to hepevt
1055          CALL FILHEP(0,1,ID1,K,K,0,0,PHOT1,YM1,.TRUE.)
1056          CALL FILHEP(0,1,ID2,K,K,0,0,PHOT2,YM2,.TRUE.)
1057          CALL FILHEP(0,1,ID3,K,K,0,0,PHOT3,YM3,.TRUE.)
1058         ENDIF
1059
1060
1061 C
1062       END
1063       SUBROUTINE TAUK0S(PETA,K)
1064 C subroutine to decay K0S's from taus.
1065 C this routine belongs to tauola universal interface, but uses
1066 C routines from tauola utilities. Just flat phase space, but 4 channels.
1067 C it is called at the beginning of SUBR. TAUPI0(JAK)
1068 C and as far as hepevt search it is basically the same as TAUPI0.  25.08.2005
1069
1070       COMMON / PARMAS / AMTAU,AMNUTA,AMEL,AMNUE,AMMU,AMNUMU
1071      *                 ,AMPIZ,AMPI,AMRO,GAMRO,AMA1,GAMA1
1072      *                 ,AMK,AMKZ,AMKST,GAMKST
1073 *
1074       REAL*4            AMTAU,AMNUTA,AMEL,AMNUE,AMMU,AMNUMU
1075      *                 ,AMPIZ,AMPI,AMRO,GAMRO,AMA1,GAMA1
1076      *                 ,AMK,AMKZ,AMKST,GAMKST
1077
1078 C position of taus, must be defined by host program:
1079       COMMON /TAUPOS/ NP1,NP2
1080 c
1081       REAL  RRR(1),BRSUM(3)
1082       REAL  PHOT1(4),PHOT2(4)
1083       REAL*8    X(4),    Y(4)
1084       REAL                                YM1,YM2
1085       REAL*8  R,PETA(4),XM1,XM2,XLAM
1086       REAL*8 a,b,c
1087       INTEGER K
1088       XLAM(a,b,c)=SQRT(ABS((a-b-c)**2-4.0*b*c))
1089 C position of decaying particle:
1090
1091
1092 !        DO L=1,4
1093 !          PETA(L)= phep(L,K)  ! K0S 4 momentum  (this is cloned from eta decay)
1094 !        ENDDO
1095 C       K0S cumulated branching ratios:
1096         BRSUM(1)=0.313  ! 2 PI0
1097         BRSUM(2)=1.0 ! BRSUM(1)+0.319  ! Pi+ PI-
1098         BRSUM(3)=BRSUM(2)+0.237  ! pi+ pi- pi0 rest is thus pi+pi-gamma
1099         CALL RANMAR(RRR,1)
1100
1101          IF(RRR(1).LT.BRSUM(1)) THEN  ! 2 pi0
1102           ID1= 111
1103           ID2= 111
1104           XM1=AMPIZ ! masses
1105           XM2=AMPIZ
1106          ELSEIF(RRR(1).LT.BRSUM(2)) THEN ! pi+ pi-
1107           ID1= 211
1108           ID2=-211
1109           XM1=AMPI ! masses
1110           XM2=AMPI
1111          ELSE                            ! gamma gamma unused !!!
1112           ID1= 22
1113           ID2= 22
1114           XM1= 0.0 ! masses
1115           XM2= 0.0
1116          ENDIF
1117
1118 ! random 3 vector on the sphere, of equal mass !!
1119          R=SQRT(PETA(4)**2-PETA(3)**2-PETA(2)**2-PETA(1)**2)/2D0
1120          R4=R
1121          R=SQRT(ABS(R**2-XM1**2))
1122          CALL SPHERD(R,X)
1123          X(4)=R4
1124          Y(4)=R4
1125
1126          Y(1)=-X(1)
1127          Y(2)=-X(2)
1128          Y(3)=-X(3)
1129 ! boost to lab and to real*4
1130          CALL bostdq(-1,PETA,X,X)
1131          CALL bostdq(-1,PETA,Y,Y)
1132          DO L=1,4
1133           PHOT1(L)=X(L)
1134           PHOT2(L)=Y(L)
1135          ENDDO
1136
1137          YM1=XM1
1138          YM2=XM2
1139 C to hepevt
1140          CALL FILHEP(0,1,ID1,K,K,0,0,PHOT1,YM1,.TRUE.)
1141          CALL FILHEP(0,1,ID2,K,K,0,0,PHOT2,YM2,.TRUE.)
1142
1143 C
1144       END
1145
1146       subroutine bostdq(idir,vv,pp,q)
1147 *     *******************************
1148 c Boost along arbitrary vector v (see eg. J.D. Jacson, Classical
1149 c Electrodynamics).
1150 c Four-vector pp is boosted from an actual frame to the rest frame
1151 c of the four-vector v (for idir=1) or back (for idir=-1).
1152 c q is a resulting four-vector.
1153 c Note: v must be time-like, pp may be arbitrary.
1154 c
1155 c Written by: Wieslaw Placzek            date: 22.07.1994
1156 c Last update: 3/29/95                     by: M.S.
1157 c
1158       implicit DOUBLE PRECISION (a-h,o-z)
1159       parameter (nout=6)
1160       DOUBLE PRECISION v(4),p(4),q(4),pp(4),vv(4)
1161       save
1162 !
1163       do 1 i=1,4
1164       v(i)=vv(i)
1165  1    p(i)=pp(i)
1166       amv=(v(4)**2-v(1)**2-v(2)**2-v(3)**2)
1167       if (amv.le.0d0) then
1168         write(6,*) 'bosstv: warning amv**2=',amv
1169       endif
1170       amv=sqrt(abs(amv))
1171       if (idir.eq.-1) then
1172         q(4)=( p(1)*v(1)+p(2)*v(2)+p(3)*v(3)+p(4)*v(4))/amv
1173         wsp =(q(4)+p(4))/(v(4)+amv)
1174       elseif (idir.eq.1) then
1175         q(4)=(-p(1)*v(1)-p(2)*v(2)-p(3)*v(3)+p(4)*v(4))/amv
1176         wsp =-(q(4)+p(4))/(v(4)+amv)
1177       else
1178         write(nout,*)' >>> boostv: wrong value of idir = ',idir
1179       endif
1180       q(1)=p(1)+wsp*v(1)
1181       q(2)=p(2)+wsp*v(2)
1182       q(3)=p(3)+wsp*v(3)
1183       end
1184
1185
1186
1187