tauola_extras.f

00001 
00002       SUBROUTINE CHOICE(MNUM,RR,ICHAN,PROB1,PROB2,PROB3,
00003      $            AMRX,GAMRX,AMRA,GAMRA,AMRB,GAMRB)
00004       COMMON / PARMAS / AMTAU,AMNUTA,AMEL,AMNUE,AMMU,AMNUMU
00005      *                 ,AMPIZ,AMPI,AMRO,GAMRO,AMA1,GAMA1
00006      *                 ,AMK,AMKZ,AMKST,GAMKST
00007 C
00008       REAL*4            AMTAU,AMNUTA,AMEL,AMNUE,AMMU,AMNUMU
00009      *                 ,AMPIZ,AMPI,AMRO,GAMRO,AMA1,GAMA1
00010      *                 ,AMK,AMKZ,AMKST,GAMKST
00011 C
00012       AMROP=1.1
00013       GAMROP=0.36
00014       AMOM=.782
00015       GAMOM=0.0084
00016 C     XXXXA CORRESPOND TO S2 CHANNEL !
00017       IF(MNUM.EQ.0) THEN
00018        PROB1=0.5
00019        PROB2=0.5
00020        AMRX =AMA1
00021        GAMRX=GAMA1
00022        AMRA =AMRO
00023        GAMRA=GAMRO
00024        AMRB =AMRO
00025        GAMRB=GAMRO
00026       ELSEIF(MNUM.EQ.1) THEN
00027        PROB1=0.5
00028        PROB2=0.5
00029        AMRX =1.57
00030        GAMRX=0.9
00031        AMRB =AMKST
00032        GAMRB=GAMKST
00033        AMRA =AMRO
00034        GAMRA=GAMRO
00035       ELSEIF(MNUM.EQ.2) THEN
00036        PROB1=0.5
00037        PROB2=0.5
00038        AMRX =1.57
00039        GAMRX=0.9
00040        AMRB =AMKST
00041        GAMRB=GAMKST
00042        AMRA =AMRO
00043        GAMRA=GAMRO
00044       ELSEIF(MNUM.EQ.3) THEN
00045        PROB1=0.5
00046        PROB2=0.5
00047        AMRX =1.27
00048        GAMRX=0.3
00049        AMRA =AMKST
00050        GAMRA=GAMKST
00051        AMRB =AMKST
00052        GAMRB=GAMKST
00053       ELSEIF(MNUM.EQ.4) THEN
00054        PROB1=0.5
00055        PROB2=0.5
00056        AMRX =1.27
00057        GAMRX=0.3
00058        AMRA =AMKST
00059        GAMRA=GAMKST
00060        AMRB =AMKST
00061        GAMRB=GAMKST
00062       ELSEIF(MNUM.EQ.5) THEN
00063        PROB1=0.5
00064        PROB2=0.5
00065        AMRX =1.27
00066        GAMRX=0.3
00067        AMRA =AMKST
00068        GAMRA=GAMKST
00069        AMRB =AMRO
00070        GAMRB=GAMRO
00071       ELSEIF(MNUM.EQ.6) THEN
00072        PROB1=0.4
00073        PROB2=0.4
00074        AMRX =1.27
00075        GAMRX=0.3
00076        AMRA =AMRO
00077        GAMRA=GAMRO
00078        AMRB =AMKST
00079        GAMRB=GAMKST
00080       ELSEIF(MNUM.EQ.7) THEN
00081        PROB1=0.0
00082        PROB2=1.0
00083        AMRX =1.27
00084        GAMRX=0.9
00085        AMRA =AMRO
00086        GAMRA=GAMRO
00087        AMRB =AMRO
00088        GAMRB=GAMRO
00089       ELSEIF(MNUM.EQ.8) THEN
00090        PROB1=0.0
00091        PROB2=1.0
00092        AMRX =AMROP
00093        GAMRX=GAMROP
00094        AMRB =AMOM
00095        GAMRB=GAMOM
00096        AMRA =AMRO
00097        GAMRA=GAMRO
00098       ELSEIF(MNUM.EQ.101) THEN
00099        PROB1=.35
00100        PROB2=.35
00101        AMRX =1.2
00102        GAMRX=.46
00103        AMRB =AMOM
00104        GAMRB=GAMOM
00105        AMRA =AMOM
00106        GAMRA=GAMOM
00107       ELSEIF(MNUM.EQ.102) THEN
00108        PROB1=0.0
00109        PROB2=0.0
00110        AMRX =1.4
00111        GAMRX=.6
00112        AMRB =AMOM
00113        GAMRB=GAMOM
00114        AMRA =AMOM
00115        GAMRA=GAMOM
00116       ELSE
00117        PROB1=0.0
00118        PROB2=0.0
00119        AMRX =AMA1
00120        GAMRX=GAMA1
00121        AMRA =AMRO
00122        GAMRA=GAMRO
00123        AMRB =AMRO
00124        GAMRB=GAMRO
00125       ENDIF
00126 C
00127       IF    (RR.LE.PROB1) THEN
00128        ICHAN=1
00129       ELSEIF(RR.LE.(PROB1+PROB2)) THEN
00130        ICHAN=2
00131         AX   =AMRA
00132         GX   =GAMRA
00133         AMRA =AMRB
00134         GAMRA=GAMRB
00135         AMRB =AX
00136         GAMRB=GX
00137         PX   =PROB1
00138         PROB1=PROB2
00139         PROB2=PX
00140       ELSE
00141        ICHAN=3
00142       ENDIF
00143 C
00144       PROB3=1.0-PROB1-PROB2
00145       END
00146       SUBROUTINE INITDK
00147 * ----------------------------------------------------------------------
00148 *     INITIALISATION OF TAU DECAY PARAMETERS  and routines
00149 *
00150 *     called by : KORALZ
00151 * ----------------------------------------------------------------------
00152 
00153       COMMON / DECPAR / GFERMI,GV,GA,CCABIB,SCABIB,GAMEL
00154       REAL*4            GFERMI,GV,GA,CCABIB,SCABIB,GAMEL
00155       COMMON / PARMAS / AMTAU,AMNUTA,AMEL,AMNUE,AMMU,AMNUMU
00156      *                 ,AMPIZ,AMPI,AMRO,GAMRO,AMA1,GAMA1
00157      *                 ,AMK,AMKZ,AMKST,GAMKST
00158 *
00159       REAL*4            AMTAU,AMNUTA,AMEL,AMNUE,AMMU,AMNUMU
00160      *                 ,AMPIZ,AMPI,AMRO,GAMRO,AMA1,GAMA1
00161      *                 ,AMK,AMKZ,AMKST,GAMKST
00162       COMMON / TAUBRA / GAMPRT(30),JLIST(30),NCHAN
00163       COMMON / TAUKLE / BRA1,BRK0,BRK0B,BRKS
00164       REAL*4            BRA1,BRK0,BRK0B,BRKS
00165 
00166 
00167 
00168 
00169 
00170 
00171       PARAMETER (NMODE=15,NM1=0,NM2=1,NM3=8,NM4=2,NM5=1,NM6=3)
00172       COMMON / TAUDCD /IDFFIN(9,NMODE),MULPIK(NMODE)
00173      &                ,NAMES
00174       CHARACTER NAMES(NMODE)*31
00175 
00176       CHARACTER OLDNAMES(7)*31
00177       CHARACTER*80 bxINIT
00178       PARAMETER (
00179      $  bxINIT ='(1x,1h*,g17.8,            16x, a31,a4,a4, 1x,1h*)'
00180      $ )
00181       REAL*4 PI
00182 *
00183 *
00184 * LIST OF BRANCHING RATIOS
00185 CAM normalised to e nu nutau channel
00186 CAM                  enu   munu   pinu  rhonu   A1nu   Knu    K*nu   pi
00187 CAM   DATA JLIST  /    1,     2,     3,     4,     5,     6,     7,
00188 *AM   DATA GAMPRT /1.000,0.9730,0.6054,1.2432,0.8432,0.0432,O.O811,0.616
00189 *AM
00190 *AM  multipion decays
00191 *
00192 *    conventions of particles names
00193 *                 K-,P-,K+,  K0,P-,KB,  K-,P0,K0
00194 *                  3, 1,-3  , 4, 1,-4  , 3, 2, 4  ,
00195 *                 P0,P0,K-,  K-,P-,P+,  P-,KB,P0
00196 *                  2, 2, 3  , 3, 1,-1  , 1,-4, 2  ,
00197 *                 ET,P-,P0   P-,P0,GM
00198 *                  9, 1, 2  , 1, 2, 8
00199 *
00200 C
00201       DIMENSION NOPIK(6,NMODE),NPIK(NMODE)
00202 *AM   outgoing multiplicity and flavors of multi-pion /multi-K modes    
00203       DATA   NPIK  /                4,                    4,  
00204      1                              5,                    5,
00205      2                              6,                    6,
00206      3                              3,                    3,            
00207      4                              3,                    3,            
00208      5                              3,                    3,            
00209      6                              3,                    3,  
00210      7                              2                         /         
00211       DATA  NOPIK / -1,-1, 1, 2, 0, 0,     2, 2, 2,-1, 0, 0,  
00212      1              -1,-1, 1, 2, 2, 0,    -1,-1,-1, 1, 1, 0,  
00213      2              -1,-1,-1, 1, 1, 2,    -1,-1, 1, 2, 2, 2, 
00214      3              -3,-1, 3, 0, 0, 0,    -4,-1, 4, 0, 0, 0,  
00215      4              -3, 2,-4, 0, 0, 0,     2, 2,-3, 0, 0, 0,  
00216      5              -3,-1, 1, 0, 0, 0,    -1, 4, 2, 0, 0, 0,  
00217      6               9,-1, 2, 0, 0, 0,    -1, 2, 8, 0, 0, 0,
00218 C AJWMOD fix sign bug, 2/22/99
00219      7              -3,-4, 0, 0, 0, 0                         /
00220 * LIST OF BRANCHING RATIOS
00221       NCHAN = NMODE + 7
00222       DO 1 I = 1,30
00223       IF (I.LE.NCHAN) THEN
00224         JLIST(I) = I
00225         IF(I.EQ. 1) GAMPRT(I) =0.1800 
00226         IF(I.EQ. 2) GAMPRT(I) =0.1751 
00227         IF(I.EQ. 3) GAMPRT(I) =0.1110 
00228         IF(I.EQ. 4) GAMPRT(I) =0.2515 
00229         IF(I.EQ. 5) GAMPRT(I) =0.1790 
00230         IF(I.EQ. 6) GAMPRT(I) =0.0071 
00231         IF(I.EQ. 7) GAMPRT(I) =0.0134
00232         IF(I.EQ. 8) GAMPRT(I) =0.0450
00233         IF(I.EQ. 9) GAMPRT(I) =0.0100
00234         IF(I.EQ.10) GAMPRT(I) =0.0009
00235         IF(I.EQ.11) GAMPRT(I) =0.0004 
00236         IF(I.EQ.12) GAMPRT(I) =0.0003 
00237         IF(I.EQ.13) GAMPRT(I) =0.0005 
00238         IF(I.EQ.14) GAMPRT(I) =0.0015 
00239         IF(I.EQ.15) GAMPRT(I) =0.0015 
00240         IF(I.EQ.16) GAMPRT(I) =0.0015 
00241         IF(I.EQ.17) GAMPRT(I) =0.0005
00242         IF(I.EQ.18) GAMPRT(I) =0.0050
00243         IF(I.EQ.19) GAMPRT(I) =0.0055
00244         IF(I.EQ.20) GAMPRT(I) =0.0017 
00245         IF(I.EQ.21) GAMPRT(I) =0.0013 
00246         IF(I.EQ.22) GAMPRT(I) =0.0010 
00247         IF(I.EQ. 1) OLDNAMES(I)='  TAU-  -->   E-               '
00248         IF(I.EQ. 2) OLDNAMES(I)='  TAU-  -->  MU-               '
00249         IF(I.EQ. 3) OLDNAMES(I)='  TAU-  -->  PI-               '
00250         IF(I.EQ. 4) OLDNAMES(I)='  TAU-  -->  PI-, PI0          '
00251         IF(I.EQ. 5) OLDNAMES(I)='  TAU-  -->  A1- (two subch)   '
00252         IF(I.EQ. 6) OLDNAMES(I)='  TAU-  -->   K-               '
00253         IF(I.EQ. 7) OLDNAMES(I)='  TAU-  -->  K*- (two subch)   '
00254         IF(I.EQ. 8) NAMES(I-7)='  TAU-  --> 2PI-,  PI0,  PI+   '
00255         IF(I.EQ. 9) NAMES(I-7)='  TAU-  --> 3PI0,        PI-   '
00256         IF(I.EQ.10) NAMES(I-7)='  TAU-  --> 2PI-,  PI+, 2PI0   '
00257         IF(I.EQ.11) NAMES(I-7)='  TAU-  --> 3PI-, 2PI+,        '
00258         IF(I.EQ.12) NAMES(I-7)='  TAU-  --> 3PI-, 2PI+,  PI0   '
00259         IF(I.EQ.13) NAMES(I-7)='  TAU-  --> 2PI-,  PI+, 3PI0   '
00260         IF(I.EQ.14) NAMES(I-7)='  TAU-  -->  K-, PI-,  K+      '
00261         IF(I.EQ.15) NAMES(I-7)='  TAU-  -->  K0, PI-, K0B      '
00262         IF(I.EQ.16) NAMES(I-7)='  TAU-  -->  K-,  K0, PI0      '
00263         IF(I.EQ.17) NAMES(I-7)='  TAU-  --> PI0  PI0   K-      '
00264         IF(I.EQ.18) NAMES(I-7)='  TAU-  -->  K-  PI-  PI+      '
00265         IF(I.EQ.19) NAMES(I-7)='  TAU-  --> PI-  K0B  PI0      '
00266         IF(I.EQ.20) NAMES(I-7)='  TAU-  --> ETA  PI-  PI0      '
00267         IF(I.EQ.21) NAMES(I-7)='  TAU-  --> PI-  PI0  GAM      '
00268         IF(I.EQ.22) NAMES(I-7)='  TAU-  -->  K-  K0            '
00269       ELSE
00270         JLIST(I) = 0
00271         GAMPRT(I) = 0.
00272       ENDIF
00273    1  CONTINUE
00274       DO I=1,NMODE
00275         MULPIK(I)=NPIK(I)
00276         DO J=1,MULPIK(I)
00277          IDFFIN(J,I)=NOPIK(J,I)
00278         ENDDO
00279       ENDDO
00280 *
00281 *
00282 * --- COEFFICIENTS TO FIX RATIO OF:
00283 * --- A1 3CHARGED/ A1 1CHARGED 2 NEUTRALS MATRIX ELEMENTS (MASLESS LIM.)
00284 * --- PROBABILITY OF K0 TO BE KS
00285 * --- PROBABILITY OF K0B TO BE KS
00286 * --- RATIO OF COEFFICIENTS FOR K*--> K0 PI-
00287 * --- ALL COEFFICENTS SHOULD BE IN THE RANGE (0.0,1.0)
00288 * --- THEY MEANING IS PROBABILITY OF THE FIRST CHOICE ONLY IF ONE
00289 * --- NEGLECTS MASS-PHASE SPACE EFFECTS
00290       BRA1=0.5
00291       BRK0=0.5
00292       BRK0B=0.5
00293       BRKS=0.6667
00294 *
00295 
00296       GFERMI = 1.16637E-5
00297       CCABIB = 0.975
00298       GV     = 1.0
00299       GA     =-1.0
00300 
00301 
00302 
00303 * ZW 13.04.89 HERE WAS AN ERROR
00304       SCABIB = SQRT(1.-CCABIB**2)
00305       PI =4.*ATAN(1.)
00306       GAMEL  = GFERMI**2*AMTAU**5/(192*PI**3)
00307 *
00308 *      CALL DEXAY(-1,pol1)
00309 *
00310       RETURN
00311       END
00312       FUNCTION DCDMAS(IDENT)
00313       COMMON / PARMAS / AMTAU,AMNUTA,AMEL,AMNUE,AMMU,AMNUMU
00314      *                 ,AMPIZ,AMPI,AMRO,GAMRO,AMA1,GAMA1
00315      *                 ,AMK,AMKZ,AMKST,GAMKST
00316 *
00317       REAL*4            AMTAU,AMNUTA,AMEL,AMNUE,AMMU,AMNUMU
00318      *                 ,AMPIZ,AMPI,AMRO,GAMRO,AMA1,GAMA1
00319      *                 ,AMK,AMKZ,AMKST,GAMKST
00320       IF      (IDENT.EQ. 1) THEN
00321         APKMAS=AMPI
00322       ELSEIF  (IDENT.EQ.-1) THEN
00323         APKMAS=AMPI
00324       ELSEIF  (IDENT.EQ. 2) THEN
00325         APKMAS=AMPIZ
00326       ELSEIF  (IDENT.EQ.-2) THEN
00327         APKMAS=AMPIZ
00328       ELSEIF  (IDENT.EQ. 3) THEN
00329         APKMAS=AMK
00330       ELSEIF  (IDENT.EQ.-3) THEN
00331         APKMAS=AMK
00332       ELSEIF  (IDENT.EQ. 4) THEN
00333         APKMAS=AMKZ
00334       ELSEIF  (IDENT.EQ.-4) THEN
00335         APKMAS=AMKZ
00336       ELSEIF  (IDENT.EQ. 8) THEN
00337         APKMAS=0.0001
00338       ELSEIF  (IDENT.EQ.-8) THEN
00339         APKMAS=0.0001
00340       ELSEIF  (IDENT.EQ. 9) THEN
00341         APKMAS=0.5488
00342       ELSEIF  (IDENT.EQ.-9) THEN
00343         APKMAS=0.5488
00344       ELSE
00345         PRINT *, 'STOP IN APKMAS, WRONG IDENT=',IDENT
00346         STOP
00347       ENDIF
00348       DCDMAS=APKMAS
00349       END
00350       FUNCTION LUNPIK(ID,ISGN)
00351       COMMON / TAUKLE / BRA1,BRK0,BRK0B,BRKS
00352       REAL*4            BRA1,BRK0,BRK0B,BRKS
00353       REAL*4 XIO(1)
00354       IDENT=ID*ISGN
00355       IF      (IDENT.EQ. 1) THEN
00356         IPKDEF=-211
00357       ELSEIF  (IDENT.EQ.-1) THEN
00358         IPKDEF= 211
00359       ELSEIF  (IDENT.EQ. 2) THEN
00360         IPKDEF=111
00361       ELSEIF  (IDENT.EQ.-2) THEN
00362         IPKDEF=111
00363       ELSEIF  (IDENT.EQ. 3) THEN
00364         IPKDEF=-321
00365       ELSEIF  (IDENT.EQ.-3) THEN
00366         IPKDEF= 321
00367       ELSEIF  (IDENT.EQ. 4) THEN
00368 *
00369 * K0 --> K0_LONG (IS 130) / K0_SHORT (IS 310) = 1/1
00370         CALL RANMAR(XIO,1)
00371         IF (XIO(1).GT.BRK0) THEN
00372           IPKDEF= 130
00373         ELSE
00374           IPKDEF= 310
00375         ENDIF
00376       ELSEIF  (IDENT.EQ.-4) THEN
00377 *
00378 * K0B--> K0_LONG (IS 130) / K0_SHORT (IS 310) = 1/1
00379         CALL RANMAR(XIO,1)
00380         IF (XIO(1).GT.BRK0B) THEN
00381           IPKDEF= 130
00382         ELSE
00383           IPKDEF= 310
00384         ENDIF
00385       ELSEIF  (IDENT.EQ. 8) THEN
00386         IPKDEF= 22
00387       ELSEIF  (IDENT.EQ.-8) THEN
00388         IPKDEF= 22
00389       ELSEIF  (IDENT.EQ. 9) THEN
00390         IPKDEF= 221
00391       ELSEIF  (IDENT.EQ.-9) THEN
00392         IPKDEF= 221
00393       ELSE
00394         PRINT *, 'STOP IN IPKDEF, WRONG IDENT=',IDENT
00395         STOP
00396       ENDIF
00397       LUNPIK=IPKDEF
00398       END
00399 
00400 
00401 
00402       SUBROUTINE TAURDF(KTO)
00403 C THIS ROUTINE CAN BE CALLED BEFORE ANY TAU+ OR TAU- EVENT IS GENERATED
00404 C IT CAN BE USED TO GENERATE TAU+ AND TAU- SAMPLES OF DIFFERENT
00405 C CONTENTS
00406       COMMON / TAUKLE / BRA1,BRK0,BRK0B,BRKS
00407       REAL*4            BRA1,BRK0,BRK0B,BRKS
00408       COMMON / TAUBRA / GAMPRT(30),JLIST(30),NCHAN
00409 
00410 C Subroutine TAURDF is disabled
00411       RETURN
00412 
00413 
00414       IF (KTO.EQ.1) THEN
00415 C     ==================
00416 C AJWMOD: Set the BRs for (A1+ -> rho+ pi0) and (K*+ -> K0 pi+)
00417       BRA1 = PKORB(4,1)
00418       BRKS = PKORB(4,3)
00419       BRK0  = PKORB(4,5)
00420       BRK0B  = PKORB(4,6)
00421       ELSE
00422 C     ====
00423 C AJWMOD: Set the BRs for (A1+ -> rho+ pi0) and (K*+ -> K0 pi+)
00424       BRA1 = PKORB(4,2)
00425       BRKS = PKORB(4,4)
00426       BRK0  = PKORB(4,5)
00427       BRK0B  = PKORB(4,6)
00428       ENDIF
00429 C     =====
00430       END
00431 
00432       SUBROUTINE INIPHY(XK00)
00433 * ----------------------------------------------------------------------
00434 *     INITIALISATION OF PARAMETERS
00435 *     USED IN QED and/or GSW ROUTINES
00436 * ----------------------------------------------------------------------
00437       COMMON / QEDPRM /ALFINV,ALFPI,XK0
00438       REAL*8           ALFINV,ALFPI,XK0
00439       REAL*8 PI8,XK00
00440 *
00441       PI8    = 4.D0*DATAN(1.D0)
00442       ALFINV = 137.03604D0
00443       ALFPI  = 1D0/(ALFINV*PI8)
00444       XK0=XK00
00445       END
00446 
00447       SUBROUTINE INIMAS
00448 C ----------------------------------------------------------------------
00449 C     INITIALISATION OF MASSES
00450 C
00451 C     called by : KORALZ
00452 C ----------------------------------------------------------------------
00453       COMMON / PARMAS / AMTAU,AMNUTA,AMEL,AMNUE,AMMU,AMNUMU
00454      *                 ,AMPIZ,AMPI,AMRO,GAMRO,AMA1,GAMA1
00455      *                 ,AMK,AMKZ,AMKST,GAMKST
00456 *
00457       REAL*4            AMTAU,AMNUTA,AMEL,AMNUE,AMMU,AMNUMU
00458      *                 ,AMPIZ,AMPI,AMRO,GAMRO,AMA1,GAMA1
00459      *                 ,AMK,AMKZ,AMKST,GAMKST
00460 C
00461 C IN-COMING / OUT-GOING  FERMION MASSES
00462       AMTAU  = 1.7842
00463 C --- let us update tau mass ...
00464       AMTAU  = 1.777
00465       AMNUTA = 0.010
00466       AMEL   = 0.0005111
00467       AMNUE  = 0.0
00468       AMMU   = 0.105659 
00469       AMNUMU = 0.0
00470 *
00471 * MASSES USED IN TAU DECAYS
00472       AMPIZ  = 0.134964
00473       AMPI   = 0.139568
00474       AMRO   = 0.773
00475       GAMRO  = 0.145
00476 *C    GAMRO  = 0.666
00477       AMA1   = 1.251
00478       GAMA1  = 0.599
00479       AMK    = 0.493667
00480       AMKZ   = 0.49772
00481       AMKST  = 0.8921
00482       GAMKST = 0.0513
00483 C
00484 C
00485 C IN-COMING / OUT-GOING  FERMION MASSES
00486 !!      AMNUTA = PKORB(1,2)
00487 !!      AMNUE  = PKORB(1,4)
00488 !!      AMNUMU = PKORB(1,6)
00489 C
00490 C MASSES USED IN TAU DECAYS  Cleo settings
00491 !!      AMPIZ  = PKORB(1,7)
00492 !!      AMPI   = PKORB(1,8)
00493 !!      AMRO   = PKORB(1,9)
00494 !!      GAMRO  = PKORB(2,9)
00495       AMA1   = 1.275   !! PKORB(1,10)
00496       GAMA1  = 0.615   !! PKORB(2,10)
00497 !!      AMK    = PKORB(1,11)
00498 !!      AMKZ   = PKORB(1,12)
00499 !!      AMKST  = PKORB(1,13)
00500 !!      GAMKST = PKORB(2,13)
00501 C
00502 
00503       RETURN
00504       END
00505       SUBROUTINE ANGULU(PD1,PD2,Q1,Q2,COSTHE)
00506       REAL*8 PD1(4),PD2(4),Q1(4),Q2(4),COSTHE,P(4),QQ(4),QT(4)
00507 C take effective beam which is less massive, it should be irrelevant
00508 C but in case HEPEVT is particulary dirty may help.
00509 C this routine calculate reduced system transver and cosine of scattering 
00510 C angle.
00511 
00512       XM1=ABS(PD1(4)**2-PD1(3)**2-PD1(2)**2-PD1(1)**2)
00513       XM2=ABS(PD2(4)**2-PD2(3)**2-PD2(2)**2-PD2(1)**2)
00514       IF (XM1.LT.XM2) THEN
00515         SIGN=1D0
00516         DO K=1,4
00517           P(K)=PD1(K)
00518         ENDDO
00519       ELSE
00520         SIGN=-1D0
00521         DO K=1,4
00522           P(K)=PD2(K)
00523         ENDDO
00524       ENDIF
00525 C calculate space like part of P (in Z restframe)
00526       DO K=1,4
00527        QQ(K)=Q1(k)+Q2(K)
00528        QT(K)=Q1(K)-Q2(K)
00529       ENDDO
00530 
00531        XMQQ=SQRT(QQ(4)**2-QQ(3)**2-QQ(2)**2-QQ(1)**2)
00532 
00533        QTXQQ=QT(4)*QQ(4)-QT(3)*QQ(3)-QT(2)*QQ(2)-QT(1)*QQ(1)
00534       DO K=1,4
00535        QT(K)=QT(K)-QQ(K)*QTXQQ/XMQQ**2
00536       ENDDO
00537 
00538        PXQQ=P(4)*QQ(4)-P(3)*QQ(3)-P(2)*QQ(2)-P(1)*QQ(1)
00539       DO K=1,4
00540        P(K)=P(K)-QQ(K)*PXQQ/XMQQ**2
00541       ENDDO
00542 C calculate costhe
00543        PXP  =SQRT(p(1)**2+p(2)**2+p(3)**2-p(4)**2)
00544        QTXQT=SQRT(QT(3)**2+QT(2)**2+QT(1)**2-QT(4)**2)
00545        PXQT =P(3)*QT(3)+P(2)*QT(2)+P(1)*QT(1)-P(4)*QT(4)
00546        COSTHE=PXQT/PXP/QTXQT
00547        COSTHE=COSTHE*SIGN
00548       END
00549 
00550       FUNCTION PLZAP0(IDE,IDF,SVAR,COSTH0)
00551 C this function calculates probability for the helicity +1 +1 configuration
00552 C of taus for given Z/gamma transfer and COSTH0 cosine of scattering angle
00553       REAL*8 PLZAP0,SVAR,COSTHE,COSTH0,T_BORN
00554 
00555       COSTHE=COSTH0
00556 C >>>>>      IF (IDE*IDF.LT.0) COSTHE=-COSTH0 ! this is probably not needed ID
00557 C >>>>>      of first beam is used by T_GIVIZ0 including sign
00558 
00559       IF (IDF.GT.0) THEN
00560         CALL INITWK(IDE,IDF,SVAR)
00561       ELSE
00562         CALL INITWK(-IDE,-IDF,SVAR)
00563       ENDIF
00564       PLZAP0=T_BORN(0,SVAR,COSTHE,1D0,1D0)
00565      $  /(T_BORN(0,SVAR,COSTHE,1D0,1D0)+T_BORN(0,SVAR,COSTHE,-1D0,-1D0))
00566 
00567 
00568 C      write(*,*) 'svar=  ',  svar
00569 C      write(*,*) 'COSTHE=',  COSTHE
00570 C      write(*,*) ide,'  ',idf
00571 C      write(*,*) 'PLZAP0=', PLZAP0
00572 C      COSTHE=0.999
00573 C      write(*,*) 'TBORN+=', T_BORN(0,SVAR,COSTHE,1D0,1D0)
00574 
00575 C      COSTHE=-0.999
00576 C      write(*,*) 'TBORN-=', T_BORN(0,SVAR,COSTHE,-1D0,-1D0)
00577 C 100  format (A,E8.4)
00578 
00579 !      PLZAP0=0.5
00580       END
00581       FUNCTION T_BORN(MODE,SVAR,COSTHE,TA,TB)
00582 C ----------------------------------------------------------------------
00583 C THIS ROUTINE PROVIDES BORN CROSS SECTION. IT HAS THE SAME         
00584 C STRUCTURE AS FUNTIS AND FUNTIH, THUS CAN BE USED AS SIMPLER       
00585 C EXAMPLE OF THE METHOD APPLIED THERE                               
00586 C INPUT PARAMETERS ARE: SVAR    -- transfer
00587 C                       COSTHE  -- cosine of angle between tau+ and 1st beam
00588 C                       TA,TB   -- helicity states of tau+ tau-
00589 C
00590 C     called by : BORNY, BORAS, BORNV, WAGA, WEIGHT
00591 C ----------------------------------------------------------------------
00592       IMPLICIT REAL*8(A-H,O-Z)
00593       COMMON / T_BEAMPM / ENE ,AMIN,AMFIN,IDE,IDF
00594       REAL*8              ENE ,AMIN,AMFIN
00595       COMMON / T_GAUSPM /SS,POLN,T3E,QE,T3F,QF
00596      &                  ,XUPGI   ,XUPZI   ,XUPGF   ,XUPZF
00597      &                  ,NDIAG0,NDIAGA,KEYA,KEYZ
00598      &                  ,ITCE,JTCE,ITCF,JTCF,KOLOR
00599       REAL*8             SS,POLN,T3E,QE,T3F,QF
00600      &                  ,XUPGI(2),XUPZI(2),XUPGF(2),XUPZF(2)
00601       REAL*8            SEPS1,SEPS2
00602 C=====================================================================
00603       COMMON / T_GSWPRM /SWSQ,AMW,AMZ,AMH,AMTOP,GAMMZ
00604       REAL*8             SWSQ,AMW,AMZ,AMH,AMTOP,GAMMZ
00605 C     SWSQ        = sin2 (theta Weinberg)
00606 C     AMW,AMZ     = W & Z boson masses respectively
00607 C     AMH         = the Higgs mass
00608 C     AMTOP       = the top mass
00609 C     GAMMZ       = Z0 width
00610       COMPLEX*16 ABORN(2,2),APHOT(2,2),AZETT(2,2)
00611       COMPLEX*16 XUPZFP(2),XUPZIP(2)
00612       COMPLEX*16 ABORNM(2,2),APHOTM(2,2),AZETTM(2,2)
00613       COMPLEX*16 PROPA,PROPZ
00614       COMPLEX*16 XR,XI
00615       COMPLEX*16 XUPF,XUPI
00616       COMPLEX*16 XTHING
00617       DATA XI/(0.D0,1.D0)/,XR/(1.D0,0.D0)/
00618       DATA MODE0 /-5/
00619       DATA IDE0 /-55/
00620       DATA SVAR0,COST0 /-5.D0,-6.D0/
00621       DATA PI /3.141592653589793238462643D0/
00622       DATA SEPS1,SEPS2 /0D0,0D0/
00623 C
00624 C MEMORIZATION =========================================================
00625       IF ( MODE.NE.MODE0.OR.SVAR.NE.SVAR0.OR.COSTHE.NE.COST0
00626      $    .OR.IDE0.NE.IDE)THEN
00627 C
00628         KEYGSW=1
00629 C ** PROPAGATORS
00630         IDE0=IDE
00631         MODE0=MODE
00632         SVAR0=SVAR
00633         COST0=COSTHE
00634         SINTHE=SQRT(1.D0-COSTHE**2)
00635         BETA=SQRT(MAX(0D0,1D0-4D0*AMFIN**2/SVAR))
00636 C I MULTIPLY AXIAL COUPLING BY BETA FACTOR.
00637         XUPZFP(1)=0.5D0*(XUPZF(1)+XUPZF(2))+0.5*BETA*(XUPZF(1)-XUPZF(2))
00638         XUPZFP(2)=0.5D0*(XUPZF(1)+XUPZF(2))-0.5*BETA*(XUPZF(1)-XUPZF(2))
00639         XUPZIP(1)=0.5D0*(XUPZI(1)+XUPZI(2))+0.5*(XUPZI(1)-XUPZI(2))
00640         XUPZIP(2)=0.5D0*(XUPZI(1)+XUPZI(2))-0.5*(XUPZI(1)-XUPZI(2))
00641 C FINAL STATE VECTOR COUPLING
00642         XUPF     =0.5D0*(XUPZF(1)+XUPZF(2))
00643         XUPI     =0.5D0*(XUPZI(1)+XUPZI(2))
00644         XTHING   =0D0
00645 
00646         PROPA =1D0/SVAR
00647         PROPZ =1D0/DCMPLX(SVAR-AMZ**2,SVAR/AMZ*GAMMZ)
00648         IF (KEYGSW.EQ.0) PROPZ=0.D0
00649         DO 50 I=1,2
00650          DO 50 J=1,2
00651           REGULA= (3-2*I)*(3-2*J) + COSTHE
00652           REGULM=-(3-2*I)*(3-2*J) * SINTHE *2.D0*AMFIN/SQRT(SVAR)
00653           APHOT(I,J)=PROPA*(XUPGI(I)*XUPGF(J)*REGULA)
00654           AZETT(I,J)=PROPZ*(XUPZIP(I)*XUPZFP(J)+XTHING)*REGULA
00655           ABORN(I,J)=APHOT(I,J)+AZETT(I,J)
00656           APHOTM(I,J)=PROPA*DCMPLX(0D0,1D0)*XUPGI(I)*XUPGF(J)*REGULM
00657           AZETTM(I,J)=PROPZ*DCMPLX(0D0,1D0)*(XUPZIP(I)*XUPF+XTHING)*REGULM
00658           ABORNM(I,J)=APHOTM(I,J)+AZETTM(I,J)
00659    50   CONTINUE
00660       ENDIF
00661 C
00662 C******************
00663 C* IN CALCULATING CROSS SECTION ONLY DIAGONAL ELEMENTS
00664 C* OF THE SPIN DENSITY MATRICES ENTER (LONGITUD. POL. ONLY.)
00665 C* HELICITY CONSERVATION EXPLICITLY OBEYED
00666       POLAR1=  (SEPS1)
00667       POLAR2= (-SEPS2)
00668       BORN=0D0
00669       DO 150 I=1,2
00670        HELIC= 3-2*I
00671        DO 150 J=1,2
00672         HELIT=3-2*J
00673         FACTOR=KOLOR*(1D0+HELIC*POLAR1)*(1D0-HELIC*POLAR2)/4D0
00674         FACTOM=FACTOR*(1+HELIT*TA)*(1-HELIT*TB)
00675         FACTOR=FACTOR*(1+HELIT*TA)*(1+HELIT*TB)
00676 
00677         BORN=BORN+CDABS(ABORN(I,J))**2*FACTOR
00678 C      MASS TERM IN BORN
00679         IF (MODE.GE.1) THEN
00680          BORN=BORN+CDABS(ABORNM(I,J))**2*FACTOM
00681         ENDIF
00682 
00683   150 CONTINUE
00684 C************
00685       FUNT=BORN
00686       IF(FUNT.LT.0.D0)  FUNT=BORN
00687 
00688 C
00689       IF (SVAR.GT.4D0*AMFIN**2) THEN
00690 C PHASE SPACE THRESHOLD FACTOR
00691         THRESH=SQRT(1-4D0*AMFIN**2/SVAR)
00692         T_BORN= FUNT*SVAR**2*THRESH
00693       ELSE
00694         THRESH=0.D0
00695         T_BORN=0.D0
00696       ENDIF
00697 C ZW HERE WAS AN ERROR 19. 05. 1989
00698 !      write(*,*) 'KKKK ',PROPA,PROPZ,XUPGI,XUPGF,XUPZI,XUPZF
00699 !      write(*,*) 'KKKK X',svar,costhe,TA,TB,T_BORN
00700       END
00701 
00702       SUBROUTINE INITWK(IDEX,IDFX,SVAR)
00703 ! initialization routine coupling masses etc.
00704       IMPLICIT REAL*8 (A-H,O-Z)
00705       COMMON / T_BEAMPM / ENE ,AMIN,AMFIN,IDE,IDF
00706       REAL*8              ENE ,AMIN,AMFIN
00707       COMMON / T_GAUSPM /SS,POLN,T3E,QE,T3F,QF
00708      &                  ,XUPGI   ,XUPZI   ,XUPGF   ,XUPZF
00709      &                  ,NDIAG0,NDIAGA,KEYA,KEYZ
00710      &                  ,ITCE,JTCE,ITCF,JTCF,KOLOR
00711       REAL*8             SS,POLN,T3E,QE,T3F,QF
00712      &                  ,XUPGI(2),XUPZI(2),XUPGF(2),XUPZF(2)
00713       COMMON / T_GSWPRM /SWSQ,AMW,AMZ,AMH,AMTOP,GAMMZ
00714       REAL*8             SWSQ,AMW,AMZ,AMH,AMTOP,GAMMZ
00715 C     SWSQ        = sin2 (theta Weinberg)
00716 C     AMW,AMZ     = W & Z boson masses respectively
00717 C     AMH         = the Higgs mass
00718 C     AMTOP       = the top mass
00719 C     GAMMZ       = Z0 width
00720 C
00721       ENE=SQRT(SVAR)/2
00722       AMIN=0.511D-3
00723       SWSQ=0.23147
00724       AMZ=91.1882
00725       GAMMZ=2.4952
00726       IF     (IDFX.EQ. 15) then       
00727         IDF=2  ! denotes tau +2 tau-
00728         AMFIN=1.77703 !this mass is irrelevant if small, used in ME only
00729       ELSEIF (IDFX.EQ.-15) then
00730         IDF=-2  ! denotes tau -2 tau-
00731         AMFIN=1.77703 !this mass is irrelevant if small, used in ME only
00732       ELSE
00733         WRITE(*,*) 'INITWK: WRONG IDFX'
00734         STOP
00735       ENDIF
00736 
00737       IF     (IDEX.EQ. 11) then      !electron
00738         IDE= 2
00739         AMIN=0.511D-3
00740       ELSEIF (IDEX.EQ.-11) then      !positron
00741         IDE=-2
00742         AMIN=0.511D-3
00743       ELSEIF (IDEX.EQ. 13) then      !mu+
00744         IDE= 2
00745         AMIN=0.105659
00746       ELSEIF (IDEX.EQ.-13) then      !mu-
00747         IDE=-2
00748         AMIN=0.105659
00749       ELSEIF (IDEX.EQ.  1) then      !d
00750         IDE= 4
00751         AMIN=0.05
00752       ELSEIF (IDEX.EQ.- 1) then      !d~
00753         IDE=-4
00754         AMIN=0.05
00755       ELSEIF (IDEX.EQ.  2) then      !u
00756         IDE= 3
00757         AMIN=0.02
00758       ELSEIF (IDEX.EQ.- 2) then      !u~
00759         IDE=-3
00760         AMIN=0.02
00761       ELSEIF (IDEX.EQ.  3) then      !s
00762         IDE= 4
00763         AMIN=0.3
00764       ELSEIF (IDEX.EQ.- 3) then      !s~
00765         IDE=-4
00766         AMIN=0.3
00767       ELSEIF (IDEX.EQ.  4) then      !c
00768         IDE= 3
00769         AMIN=1.3
00770       ELSEIF (IDEX.EQ.- 4) then      !c~
00771         IDE=-3
00772         AMIN=1.3
00773       ELSEIF (IDEX.EQ.  5) then      !b
00774         IDE= 4
00775         AMIN=4.5
00776       ELSEIF (IDEX.EQ.- 5) then      !b~
00777         IDE=-4
00778         AMIN=4.5
00779       ELSEIF (IDEX.EQ.  12) then     !nu_e
00780         IDE= 1
00781         AMIN=0.1D-3
00782       ELSEIF (IDEX.EQ.- 12) then     !nu_e~
00783         IDE=-1
00784         AMIN=0.1D-3
00785       ELSEIF (IDEX.EQ.  14) then     !nu_mu
00786         IDE= 1
00787         AMIN=0.1D-3
00788       ELSEIF (IDEX.EQ.- 14) then     !nu_mu~
00789         IDE=-1
00790         AMIN=0.1D-3
00791       ELSEIF (IDEX.EQ.  16) then     !nu_tau
00792         IDE= 1
00793         AMIN=0.1D-3
00794       ELSEIF (IDEX.EQ.- 16) then     !nu_tau~
00795         IDE=-1
00796         AMIN=0.1D-3
00797 
00798       ELSE
00799         WRITE(*,*) 'INITWK: WRONG IDEX'
00800         STOP
00801       ENDIF
00802 
00803 C ----------------------------------------------------------------------
00804 C
00805 C     INITIALISATION OF COUPLING CONSTANTS AND FERMION-GAMMA / Z0 VERTEX
00806 C
00807 C     called by : KORALZ
00808 C ----------------------------------------------------------------------
00809       ITCE=IDE/IABS(IDE)
00810       JTCE=(1-ITCE)/2
00811       ITCF=IDF/IABS(IDF)
00812       JTCF=(1-ITCF)/2
00813       CALL T_GIVIZO( IDE, 1,AIZOR,QE,KDUMM)
00814       CALL T_GIVIZO( IDE,-1,AIZOL,QE,KDUMM)
00815       XUPGI(1)=QE
00816       XUPGI(2)=QE
00817       T3E    = AIZOL+AIZOR
00818       XUPZI(1)=(AIZOR-QE*SWSQ)/SQRT(SWSQ*(1-SWSQ))
00819       XUPZI(2)=(AIZOL-QE*SWSQ)/SQRT(SWSQ*(1-SWSQ))
00820       CALL T_GIVIZO( IDF, 1,AIZOR,QF,KOLOR)
00821       CALL T_GIVIZO( IDF,-1,AIZOL,QF,KOLOR)
00822       XUPGF(1)=QF
00823       XUPGF(2)=QF
00824       T3F    =  AIZOL+AIZOR
00825       XUPZF(1)=(AIZOR-QF*SWSQ)/SQRT(SWSQ*(1-SWSQ))
00826       XUPZF(2)=(AIZOL-QF*SWSQ)/SQRT(SWSQ*(1-SWSQ))
00827 C
00828       NDIAG0=2
00829       NDIAGA=11
00830       KEYA  = 1
00831       KEYZ  = 1
00832 C
00833 C
00834       RETURN
00835       END
00836 
00837       SUBROUTINE T_GIVIZO(IDFERM,IHELIC,SIZO3,CHARGE,KOLOR)
00838 C ----------------------------------------------------------------------
00839 C PROVIDES ELECTRIC CHARGE AND WEAK IZOSPIN OF A FAMILY FERMION
00840 C IDFERM=1,2,3,4 DENOTES NEUTRINO, LEPTON, UP AND DOWN QUARK
00841 C NEGATIVE IDFERM=-1,-2,-3,-4, DENOTES ANTIPARTICLE
00842 C IHELIC=+1,-1 DENOTES RIGHT AND LEFT HANDEDNES ( CHIRALITY)
00843 C SIZO3 IS THIRD PROJECTION OF WEAK IZOSPIN (PLUS MINUS HALF)
00844 C AND CHARGE IS ELECTRIC CHARGE IN UNITS OF ELECTRON CHARGE
00845 C KOLOR IS A QCD COLOUR, 1 FOR LEPTON, 3 FOR QUARKS
00846 C
00847 C     called by : EVENTE, EVENTM, FUNTIH, .....
00848 C ----------------------------------------------------------------------
00849       IMPLICIT REAL*8(A-H,O-Z)
00850 C
00851       IF(IDFERM.EQ.0.OR.IABS(IDFERM).GT.4) GOTO 901
00852       IF(IABS(IHELIC).NE.1)                GOTO 901
00853       IH  =IHELIC
00854       IDTYPE =IABS(IDFERM)
00855       IC  =IDFERM/IDTYPE
00856       LEPQUA=INT(IDTYPE*0.4999999D0)
00857       IUPDOW=IDTYPE-2*LEPQUA-1
00858       CHARGE  =(-IUPDOW+2D0/3D0*LEPQUA)*IC
00859       SIZO3   =0.25D0*(IC-IH)*(1-2*IUPDOW)
00860       KOLOR=1+2*LEPQUA
00861 C** NOTE THAT CONVENTIONALY Z0 COUPLING IS
00862 C** XOUPZ=(SIZO3-CHARGE*SWSQ)/SQRT(SWSQ*(1-SWSQ))
00863       RETURN
00864  901  PRINT *,' STOP IN GIVIZO: WRONG PARAMS.'
00865       STOP
00866       END
00867       SUBROUTINE PHYFIX(NSTOP,NSTART)
00868       COMMON/LUJETS/N,K(4000,5),P(4000,5),V(4000,5) 
00869       SAVE /LUJETS/ 
00870 C NSTOP NSTART : when PHYTIA history ends and event starts.
00871       NSTOP=0
00872       NSTART=1
00873       DO I=1, N
00874        IF(K(I,1).NE.21) THEN
00875            NSTOP = I-1
00876            NSTART= I
00877            GOTO 500
00878        ENDIF
00879       ENDDO
00880  500  CONTINUE
00881       END
00882  
00883 
00884       SUBROUTINE TAUPI0(PI0,K)
00885 C no initialization required. Must be called once after every:
00886 C   1)    CALL DEKAY(1+10,...)
00887 C   2)    CALL DEKAY(2+10,...)
00888 C   3)    CALL DEXAY(1,...)
00889 C   4)    CALL DEXAY(2,...)
00890 C subroutine to decay originating from TAUOLA's taus: 
00891 C 1) etas (with CALL TAUETA(JAK))
00892 C 2) later pi0's from taus.
00893 C 3) extensions to other applications possible. 
00894 C this routine belongs to >tauola universal interface<, but uses 
00895 C routines from >tauola< utilities as well.  25.08.2005      
00896 C this is the hepevt class in old style. No d_h_ class pre-name
00897  
00898 C position of taus, must be defined by host program:
00899       COMMON /TAUPOS/ NP1,NP2
00900 c
00901       REAL  PHOT1(4),PHOT2(4)
00902       REAL*8  R,X(4),Y(4),PI0(4)
00903       INTEGER JEZELI(3),K
00904       DATA JEZELI /0,0,0/
00905       SAVE JEZELI
00906 
00907 ! random 3 vector on the sphere, masless
00908         R=SQRT(PI0(4)**2-PI0(3)**2-PI0(2)**2-PI0(1)**2)/2D0
00909         CALL SPHERD(R,X)
00910         X(4)=R
00911         Y(4)=R
00912         
00913         Y(1)=-X(1)
00914         Y(2)=-X(2)
00915         Y(3)=-X(3)
00916 ! boost to lab and to real*4
00917         CALL bostdq(-1,PI0,X,X)
00918         CALL bostdq(-1,PI0,Y,Y)
00919         DO L=1,4
00920          PHOT1(L)=X(L)
00921          PHOT2(L)=Y(L)
00922         ENDDO
00923 C to hepevt
00924         CALL FILHEP(0,1,22,K,K,0,0,PHOT1,0.0,.TRUE.)
00925         CALL FILHEP(0,1,22,K,K,0,0,PHOT2,0.0,.TRUE.)
00926 
00927 C
00928       END
00929       SUBROUTINE TAUETA(PETA,K)
00930 C subroutine to decay etas's from taus. 
00931 C this routine belongs to tauola universal interface, but uses 
00932 C routines from tauola utilities. Just flat phase space, but 4 channels.
00933 C it is called at the beginning of SUBR. TAUPI0(JAK)
00934 C and as far as hepevt search it is basically the same as TAUPI0.  25.08.2005    
00935 C this is the hepevt class in old style. No d_h_ class pre-name
00936  
00937       COMMON / PARMAS / AMTAU,AMNUTA,AMEL,AMNUE,AMMU,AMNUMU
00938      *                 ,AMPIZ,AMPI,AMRO,GAMRO,AMA1,GAMA1
00939      *                 ,AMK,AMKZ,AMKST,GAMKST
00940 *
00941       REAL*4            AMTAU,AMNUTA,AMEL,AMNUE,AMMU,AMNUMU
00942      *                 ,AMPIZ,AMPI,AMRO,GAMRO,AMA1,GAMA1
00943      *                 ,AMK,AMKZ,AMKST,GAMKST
00944 
00945 C position of taus, must be defined by host program:
00946 c
00947       REAL  RRR(1),BRSUM(3), RR(2)
00948       REAL  PHOT1(4),PHOT2(4),PHOT3(4)
00949       REAL*8    X(4),    Y(4),    Z(4)
00950       REAL                                YM1,YM2,YM3
00951       REAL*8  R,RU,PETA(4),XM1,XM2,XM3,XLAM
00952       REAL*8 a,b,c
00953       INTEGER K
00954       XLAM(a,b,c)=SQRT(ABS((a-b-c)**2-4.0*b*c))
00955 C position of decaying particle:
00956 C        DO L=1,4
00957 C          PETA(L)= phep(L,K)  ! eta 4 momentum
00958 C        ENDDO
00959 C       eta cumulated branching ratios:
00960         BRSUM(1)=0.389  ! gamma gamma
00961         BRSUM(2)=BRSUM(1)+0.319  ! 3 pi0
00962         BRSUM(3)=BRSUM(2)+0.237  ! pi+ pi- pi0 rest is thus pi+pi-gamma
00963         CALL RANMAR(RRR,1) 
00964         
00965         IF (RRR(1).LT.BRSUM(1)) THEN ! gamma gamma channel exactly like pi0
00966 ! random 3 vector on the sphere, masless   
00967          R=SQRT(PETA(4)**2-PETA(3)**2-PETA(2)**2-PETA(1)**2)/2D0
00968          CALL SPHERD(R,X) 
00969          X(4)=R
00970          Y(4)=R
00971         
00972          Y(1)=-X(1)
00973          Y(2)=-X(2)
00974          Y(3)=-X(3)
00975 ! boost to lab and to real*4
00976          CALL bostdq(-1,PETA,X,X)  
00977          CALL bostdq(-1,PETA,Y,Y)
00978          DO L=1,4
00979           PHOT1(L)=X(L)
00980           PHOT2(L)=Y(L)
00981          ENDDO
00982 C to hepevt
00983          CALL FILHEP(0,1,22,K,K,0,0,PHOT1,0.0,.TRUE.)
00984          CALL FILHEP(0,1,22,K,K,0,0,PHOT2,0.0,.TRUE.)
00985         ELSE ! 3 body channels
00986          IF(RRR(1).LT.BRSUM(2)) THEN  ! 3 pi0
00987           ID1= 111
00988           ID2= 111
00989           ID3= 111
00990           XM1=AMPIZ ! masses
00991           XM2=AMPIZ
00992           XM3=AMPIZ
00993          ELSEIF(RRR(1).LT.BRSUM(3)) THEN ! pi+ pi- pi0
00994           ID1= 211
00995           ID2=-211
00996           ID3= 111
00997           XM1=AMPI ! masses
00998           XM2=AMPI
00999           XM3=AMPIZ
01000          ELSE                            ! pi+ pi- gamma 
01001           ID1= 211
01002           ID2=-211
01003           ID3=  22
01004           XM1=AMPI ! masses
01005           XM2=AMPI
01006           XM3=0.0
01007          ENDIF
01008  7       CONTINUE  ! we generate mass of the first pair:
01009           CALL RANMAR(RR,2)
01010           R=SQRT(PETA(4)**2-PETA(3)**2-PETA(2)**2-PETA(1)**2)
01011           AMIN=XM1+XM2
01012           AMAX=R-XM3
01013           AM2=SQRT(AMIN**2+RR(1)*(AMAX**2-AMIN**2))
01014 C         weight for flat phase space
01015           WT=XLAM(1D0*R**2,1D0*AM2**2,1D0*XM3**2)
01016      &      *XLAM(1D0*AM2**2,1D0*XM1**2,1D0*XM2**2)
01017      &           /R**2                    /AM2**2
01018          IF (RR(2).GT.WT) GOTO 7
01019 
01020          RU=XLAM(1D0*AM2**2,1D0*XM1**2,1D0*XM2**2)/AM2/2  ! momenta of the
01021                                               ! first two products
01022                                               ! in the rest frame of that pair
01023          CALL SPHERD(RU,X)
01024          X(4)=SQRT(RU**2+XM1**2)
01025          Y(4)=SQRT(RU**2+XM2**2)
01026         
01027          Y(1)=-X(1)
01028          Y(2)=-X(2)
01029          Y(3)=-X(3)
01030 C generate momentum of that pair in rest frame of eta:
01031          RU=XLAM(1D0*R**2,1D0*AM2**2,1D0*XM3**2)/R/2
01032          CALL SPHERD(RU,Z)
01033          Z(4)=SQRT(RU**2+AM2**2)
01034 C and boost first two decay products to rest frame of eta.
01035          CALL bostdq(-1,Z,X,X)
01036          CALL bostdq(-1,Z,Y,Y)
01037 C redefine Z(4) to 4-momentum of the last decay product: 
01038          Z(1)=-Z(1)
01039          Z(2)=-Z(2)
01040          Z(3)=-Z(3)
01041          Z(4)=SQRT(RU**2+XM3**2)
01042 C boost all to lab and move to real*4; also masses
01043          CALL bostdq(-1,PETA,X,X)
01044          CALL bostdq(-1,PETA,Y,Y)
01045          CALL bostdq(-1,PETA,Z,Z)
01046          DO L=1,4
01047           PHOT1(L)=X(L)
01048           PHOT2(L)=Y(L)
01049           PHOT3(L)=Z(L)
01050          ENDDO
01051          YM1=XM1
01052          YM2=XM2
01053          YM3=XM3
01054 C to hepevt
01055          CALL FILHEP(0,1,ID1,K,K,0,0,PHOT1,YM1,.TRUE.)
01056          CALL FILHEP(0,1,ID2,K,K,0,0,PHOT2,YM2,.TRUE.)
01057          CALL FILHEP(0,1,ID3,K,K,0,0,PHOT3,YM3,.TRUE.)
01058         ENDIF
01059 
01060 
01061 C
01062       END
01063       SUBROUTINE TAUK0S(PETA,K)
01064 C subroutine to decay K0S's from taus. 
01065 C this routine belongs to tauola universal interface, but uses 
01066 C routines from tauola utilities. Just flat phase space, but 4 channels.
01067 C it is called at the beginning of SUBR. TAUPI0(JAK)
01068 C and as far as hepevt search it is basically the same as TAUPI0.  25.08.2005   
01069 
01070       COMMON / PARMAS / AMTAU,AMNUTA,AMEL,AMNUE,AMMU,AMNUMU
01071      *                 ,AMPIZ,AMPI,AMRO,GAMRO,AMA1,GAMA1
01072      *                 ,AMK,AMKZ,AMKST,GAMKST
01073 *
01074       REAL*4            AMTAU,AMNUTA,AMEL,AMNUE,AMMU,AMNUMU 
01075      *                 ,AMPIZ,AMPI,AMRO,GAMRO,AMA1,GAMA1
01076      *                 ,AMK,AMKZ,AMKST,GAMKST
01077 
01078 C position of taus, must be defined by host program:
01079       COMMON /TAUPOS/ NP1,NP2
01080 c
01081       REAL  RRR(1),BRSUM(3)
01082       REAL  PHOT1(4),PHOT2(4)
01083       REAL*8    X(4),    Y(4)
01084       REAL                                YM1,YM2
01085       REAL*8  R,PETA(4),XM1,XM2,XLAM
01086       REAL*8 a,b,c
01087       INTEGER K
01088       XLAM(a,b,c)=SQRT(ABS((a-b-c)**2-4.0*b*c))
01089 C position of decaying particle:
01090 
01091       
01092 !        DO L=1,4
01093 !          PETA(L)= phep(L,K)  ! K0S 4 momentum  (this is cloned from eta decay)
01094 !        ENDDO
01095 C       K0S cumulated branching ratios:
01096         BRSUM(1)=0.313  ! 2 PI0
01097         BRSUM(2)=1.0 ! BRSUM(1)+0.319  ! Pi+ PI-
01098         BRSUM(3)=BRSUM(2)+0.237  ! pi+ pi- pi0 rest is thus pi+pi-gamma
01099         CALL RANMAR(RRR,1) 
01100 
01101          IF(RRR(1).LT.BRSUM(1)) THEN  ! 2 pi0
01102           ID1= 111
01103           ID2= 111
01104           XM1=AMPIZ ! masses
01105           XM2=AMPIZ
01106          ELSEIF(RRR(1).LT.BRSUM(2)) THEN ! pi+ pi- 
01107           ID1= 211
01108           ID2=-211
01109           XM1=AMPI ! masses
01110           XM2=AMPI
01111          ELSE                            ! gamma gamma unused !!!
01112           ID1= 22
01113           ID2= 22
01114           XM1= 0.0 ! masses
01115           XM2= 0.0
01116          ENDIF
01117         
01118 ! random 3 vector on the sphere, of equal mass !!  
01119          R=SQRT(PETA(4)**2-PETA(3)**2-PETA(2)**2-PETA(1)**2)/2D0
01120          R4=R
01121          R=SQRT(ABS(R**2-XM1**2))
01122          CALL SPHERD(R,X) 
01123          X(4)=R4
01124          Y(4)=R4
01125         
01126          Y(1)=-X(1)
01127          Y(2)=-X(2)
01128          Y(3)=-X(3)
01129 ! boost to lab and to real*4
01130          CALL bostdq(-1,PETA,X,X)  
01131          CALL bostdq(-1,PETA,Y,Y)
01132          DO L=1,4
01133           PHOT1(L)=X(L)
01134           PHOT2(L)=Y(L)
01135          ENDDO
01136 
01137          YM1=XM1
01138          YM2=XM2
01139 C to hepevt
01140          CALL FILHEP(0,1,ID1,K,K,0,0,PHOT1,YM1,.TRUE.)
01141          CALL FILHEP(0,1,ID2,K,K,0,0,PHOT2,YM2,.TRUE.)
01142 
01143 C
01144       END
01145 
01146       subroutine bostdq(idir,vv,pp,q)
01147 *     *******************************
01148 c Boost along arbitrary vector v (see eg. J.D. Jacson, Classical 
01149 c Electrodynamics).
01150 c Four-vector pp is boosted from an actual frame to the rest frame 
01151 c of the four-vector v (for idir=1) or back (for idir=-1). 
01152 c q is a resulting four-vector.
01153 c Note: v must be time-like, pp may be arbitrary.
01154 c
01155 c Written by: Wieslaw Placzek            date: 22.07.1994
01156 c Last update: 3/29/95                     by: M.S.
01157 c 
01158       implicit DOUBLE PRECISION (a-h,o-z)
01159       parameter (nout=6)
01160       DOUBLE PRECISION v(4),p(4),q(4),pp(4),vv(4)  
01161       save
01162 !
01163       do 1 i=1,4
01164       v(i)=vv(i)
01165  1    p(i)=pp(i)
01166       amv=(v(4)**2-v(1)**2-v(2)**2-v(3)**2)
01167       if (amv.le.0d0) then
01168         write(6,*) 'bosstv: warning amv**2=',amv
01169       endif
01170       amv=sqrt(abs(amv))
01171       if (idir.eq.-1) then
01172         q(4)=( p(1)*v(1)+p(2)*v(2)+p(3)*v(3)+p(4)*v(4))/amv
01173         wsp =(q(4)+p(4))/(v(4)+amv)
01174       elseif (idir.eq.1) then
01175         q(4)=(-p(1)*v(1)-p(2)*v(2)-p(3)*v(3)+p(4)*v(4))/amv
01176         wsp =-(q(4)+p(4))/(v(4)+amv)
01177       else
01178         write(nout,*)' >>> boostv: wrong value of idir = ',idir
01179       endif
01180       q(1)=p(1)+wsp*v(1)
01181       q(2)=p(2)+wsp*v(2)
01182       q(3)=p(3)+wsp*v(3)
01183       end
01184         
01185 
01186 
01187
Generated on Sun Oct 20 20:24:11 2013 for C++InterfacetoTauola by  doxygen 1.6.3