PROGRAM PYXTRA C...Long example how to interface user-defined processes to PYTHIA C...based on the Les Houches commonblock agreement. Generates events C...of several different kinds, to test the new code under varied C...conditions, but kinematics selection and cross sections are C...completely unphysical. C...Double precision and integer declarations. IMPLICIT DOUBLE PRECISION(A-H, O-Z) INTEGER PYK,PYCHGE,PYCOMP C...User process event common block. INTEGER MAXNUP PARAMETER (MAXNUP=500) INTEGER NUP,IDPRUP,IDUP,ISTUP,MOTHUP,ICOLUP DOUBLE PRECISION XWGTUP,SCALUP,AQEDUP,AQCDUP,PUP,VTIMUP,SPINUP COMMON/HEPEUP/NUP,IDPRUP,XWGTUP,SCALUP,AQEDUP,AQCDUP,IDUP(MAXNUP), &ISTUP(MAXNUP),MOTHUP(2,MAXNUP),ICOLUP(2,MAXNUP),PUP(5,MAXNUP), &VTIMUP(MAXNUP),SPINUP(MAXNUP) C...User process initialization commonblock. INTEGER MAXPUP PARAMETER (MAXPUP=100) INTEGER IDBMUP,PDFGUP,PDFSUP,IDWTUP,NPRUP,LPRUP DOUBLE PRECISION EBMUP,XSECUP,XERRUP,XMAXUP COMMON/HEPRUP/IDBMUP(2),EBMUP(2),PDFGUP(2),PDFSUP(2), &IDWTUP,NPRUP,XSECUP(MAXPUP),XERRUP(MAXPUP),XMAXUP(MAXPUP), &LPRUP(MAXPUP) SAVE /HEPRUP/ C...Standard PYTHIA commonblocks. INTEGER N,NPAD,K REAL*8 P,V COMMON/PYJETS/N,NPAD,K(4000,5),P(4000,5),V(4000,5) INTEGER MSTU,MSTJ REAL*8 PARU,PARJ COMMON/PYDAT1/MSTU(200),PARU(200),MSTJ(200),PARJ(200) INTEGER MDCY,MDME,KFDP REAL*8 BRAT COMMON/PYDAT3/MDCY(500,3),MDME(8000,2),BRAT(8000),KFDP(8000,5) INTEGER MSEL,MSELPD,MSUB,KFIN REAL*8 CKIN COMMON/PYSUBS/MSEL,MSELPD,MSUB(500),KFIN(2,-40:40),CKIN(200) INTEGER MSTP,MSTI REAL*8 PARP,PARI COMMON/PYPARS/MSTP(200),PARP(200),MSTI(200),PARI(200) C...Extra commonblock to transfer run info. INTEGER MODE,NLIM REAL*8 SCALEF,WTMAX COMMON/PRIV/MODE,NLIM,SCALEF,WTMAX C...Local arrays. DIMENSION IPRID(100) REAL*8 XLUM character*72 FNAME character*5 cgive character*30 cgive0 integer istr integer lok,NEV,iev INTEGER LNHWRT,LNHRD,LNHDCY,LNHOUT EXTERNAL PYDATA C...Switch process mode; agrees with IDWTUP code (+-1,+-2,+-3,+-4). MODE=2 C.....2 means weighted in and unweighted out C......Events selected according to xsection CCC MODE=1 C.....1 means weighted in and unweighted out C......Events selected according to maximum weight C...Maximum number of events to generate. OPEN(55,FILE='pyxtra.in',status='old') NEV=100 READ(55,*) MODE,XLUM,FNAME,SCALEF C initialize HEP logical units lnhwrt=23 lnhrd=0 lnhdcy=0 lnhout=22 C WRITE(CGIVE,'(I5)') lnhout CGIVE0='MSTU(11)='//CGIVE CALL PYGIVE(CGIVE0) C...Initialize with external process. CALL PYINIT('USER',' ',' ',0D0) NEV = INT (XLUM * XSECUP(1)) IF(NEV.LE.0) NEV=1 write(*,*) ' requesting ',XLUM,' pbi of data ' write(*,*) ' cross section is ',xsecup(1),' pb ' write(*,*) ' events to be generated = ',nev C.....Opening stdhep file for writing call stdxwinit(FNAME,'StdHep/Pythia example', 1 nev,istr,lok) if(lok.ne.0) write(lnhout,*) 1 ' Problem opening file ' C Write Stdhep begin-run record call stdxwrt(100,istr,lok) if(lok.ne.0) write(lnhout,*) 1 ' Problem writing stdhep begin run record' C...Event loop; generate event; check it was obtained or quit. DO 130 IEV=1,NEV CALL PYEVNT IF(MSTI(51).EQ.1) GOTO 140 CALL LUNHEP(1) IF(IEV.LT.10) THEN CALL PYLIST(7) CALL PYLIST(2) ELSEIF(IEV.LT.50) THEN IF(MOD(IEV-1,10).EQ.0) THEN PRINT*,' IEV = ',IEV CALL PYLIST(7) CALL PYLIST(2) ENDIF ENDIF C.......Write one event call stdxwrt(1,istr,lok) 130 CONTINUE C...Statistics and histograms. 140 CALL PYSTAT(1) C Fill Stdhep common block 1 with run information call stdflpyxsec(nev) C Write end-of-run record call stdxwrt(200,istr,lok) if(lok.ne.0) write(lnhout,*) ' Problem writing end run record' c...close event file call stdxend(istr) END C********************************************************************* C...UPINIT C...Routine to be called by user to set up user-defined processes. C...Code below only intended as example, without any claim of realism. C...Especially it shows what info needs to be put in HEPRUP. SUBROUTINE UPINIT C...Double precision and integer declarations. IMPLICIT DOUBLE PRECISION(A-H, O-Z) IMPLICIT INTEGER(I-N) include 'tauola.inc' CHARACTER*80 CHAR_READ C...User process initialization commonblock. INTEGER MAXPUP PARAMETER (MAXPUP=100) INTEGER IDBMUP,PDFGUP,PDFSUP,IDWTUP,NPRUP,LPRUP DOUBLE PRECISION EBMUP,XSECUP,XERRUP,XMAXUP COMMON/HEPRUP/IDBMUP(2),EBMUP(2),PDFGUP(2),PDFSUP(2), &IDWTUP,NPRUP,XSECUP(MAXPUP),XERRUP(MAXPUP),XMAXUP(MAXPUP), &LPRUP(MAXPUP) SAVE /HEPRUP/ REAL*8 XSTOT,WTEVNT,WTGRT INTEGER I,KNTEV,J,NPARTS,IMAX,NGRT C...Extra commonblock to transfer run info. INTEGER MODE,NLIM REAL*8 SCALEF,WTMAX COMMON/PRIV/MODE,NLIM,SCALEF,WTMAX SAVE/PRIV/ C....Pythia commonblock - needed for setting PDF's; see below. C COMMON/PYPARS/MSTP(200),PARP(200),MSTI(200),PARI(200) C SAVE /PYPARS/ C...Set incoming beams: Tevatron Run II. IDBMUP(1)=2212 IDBMUP(2)=-2212 EBMUP(1)=980D0 EBMUP(2)=980D0 C...Set PDF's of incoming beams: CTEQ 5L. C...Note that Pythia will not look at PDFGUP and PDFSUP. PDFGUP(1)=4 PDFSUP(1)=46 PDFGUP(2)=PDFGUP(1) PDFSUP(2)=PDFSUP(1) C...If you want Pythia to use PDFLIB, you have to set it by hand. C...(You also have to ensure that the dummy routines C...PDFSET, STRUCTM and STRUCTP in Pythia are not linked.) C MSTP(52)=2 C MSTP(51)=1000*PDFGUP(1)+PDFSUP(1) C...Decide on weighting strategy: unweighted on input. IDWTUP=MODE C...Number of external processes. NPRUP=1 open(unit=77,file='bk.77',status='old') C.....read through data file to find maximum weight....only first time WTMAX=0D0 XSTOT=0D0 DO I=1,10000000 READ(77,*,END=777) NPARTS,KNTEV,WTEVNT IF(WTEVNT.GT.WTMAX) WTMAX=WTEVNT XSTOT=XSTOT+WTEVNT DO J=1,6 READ(77,*) CHAR_READ ENDDO DO J=1,NPARTS READ(77,*) CHAR_READ ENDDO ENDDO 777 CONTINUE REWIND(77) XSECUP(1)=XSTOT XMAXUP(1)=WTMAX LPRUP(1)=661 cc print*,' calling dexay ' call tauola_init cc print*,' returned from dexay ' RETURN END C********************************************************************* C...UPEVNT C...Sample routine to generate events of various kinds. C...Not intended to be realistic, but rather to show in closed C...and understandable form what such a routine might look like. C...Especially it shows what info needs to be put in HEPEUP. SUBROUTINE UPEVNT C...Double precision and integer declarations. IMPLICIT DOUBLE PRECISION(A-H, O-Z) IMPLICIT INTEGER(I-N) COMMON/PYPARS/MSTP(200),PARP(200),MSTI(200),PARI(200) C...Extra commonblock to transfer run info. INTEGER MODE,NLIM REAL*8 SCALEF,WTMAX COMMON/PRIV/MODE,NLIM,SCALEF,WTMAX SAVE/PRIV/ C...User process event common block. INTEGER MAXNUP PARAMETER (MAXNUP=500) INTEGER NUP,IDPRUP,IDUP,ISTUP,MOTHUP,ICOLUP DOUBLE PRECISION XWGTUP,SCALUP,AQEDUP,AQCDUP,PUP,VTIMUP,SPINUP COMMON/HEPEUP/NUP,IDPRUP,XWGTUP,SCALUP,AQEDUP,AQCDUP,IDUP(MAXNUP), &ISTUP(MAXNUP),MOTHUP(2,MAXNUP),ICOLUP(2,MAXNUP),PUP(5,MAXNUP), &VTIMUP(MAXNUP),SPINUP(MAXNUP) SAVE /HEPEUP/ C.....For adding in tau polarization INTEGER NMXHEP PARAMETER (NMXHEP=2000) INTEGER NEVHEP,NHEP,ISTHEP,IDHEP,JMOHEP,JDAHEP REAL PHEP,VHEP COMMON/TDEVNT/NEVHEP,NHEP,ISTHEP(NMXHEP),IDHEP(NMXHEP), &JMOHEP(2,NMXHEP),JDAHEP(2,NMXHEP),PHEP(5,NMXHEP),VHEP(4,NMXHEP) SAVE /TDEVNT/ real*8 pin(5),pout(5),ppair(5),qpair(5),v1(5),v2(5) real*4 pout4(4),pmag,pol(4),amtau C.....For choosing the parton shower scale using ktclustering integer njmax parameter(njmax=20) real*8 ecut,y(njmax),pp(4,njmax) integer nn C...If PYTHIA is supposed to select event type, do not modify this choice. IF(IABS(MODE).LE.2) THEN IDPRUP=661 C...Else free hands to mix; here evenly. ELSE ENDIF 99 CONTINUE C...Zero some arrays in common blocks to simplify filling. NUP=12 DO 100 I=1,NUP MOTHUP(1,I)=0 MOTHUP(2,I)=0 ICOLUP(1,I)=0 ICOLUP(2,I)=0 SPINUP(I)=0D0 PUP(1,I)=0D0 PUP(2,I)=0D0 PUP(5,I)=0D0 VTIMUP(I)=0D0 100 CONTINUE READ(77,*,END=888) NUP,ID,XWGTUP READ(77,*,END=888) (IDUP(I),I=1,NUP) ISTUP(1)=-1 ISTUP(2)=-1 DO I=3,NUP ISTUP(I)=1 ENDDO READ(77,*,END=888) (MOTHUP(1,I),I=1,NUP) READ(77,*,END=888) (MOTHUP(2,I),I=1,NUP) READ(77,*,END=888) (ICOLUP(1,I),I=1,NUP) READ(77,*,END=888) (ICOLUP(2,I),I=1,NUP) C.....Read in the "new" helicity information READ(77,*,END=888) (SPINUP(I),I=1,NUP) DO I=1,NUP READ(77,*,END=888) IDUMB,PUP(4,I),(PUP(J,I),J=1,3) TEST=PUP(1,I)**2+PUP(2,I)**2+PUP(3,I)**2-PUP(4,I)**2 IF(TEST.LE.0D0) THEN PUP(5,I)=DSQRT(-TEST) ELSEIF(DABS(TEST).LT.1D-6) THEN PUP(5,I)=0D0 ELSE PUP(5,I)=-1D0 PRINT*,' NEGATIVE MASS ' ENDIF ENDDO C.....Add mothers for l+ l- and l nu_l pairs.... C.....Assumption that consecutive leptons make a pair. IJ=3 555 CONTINUE IS=ISTUP(IJ) IF(IS.NE.1) GOTO 998 ID=IDUP(IJ) IDA=ABS(ID) IF(IDA.GE.11.AND.IDA.LE.16.AND.IJ.LT.NUP) THEN IF(ICOLUP(1,IJ).EQ.0 .AND. ICOLUP(2,IJ).EQ.0) THEN IDP=IDUP(IJ+1) IDPA=ABS(IDP) IMO=0 IF(IDP.EQ.-ID) THEN C.......gamma*/Z* mother IMO=23 ELSE IDL=MIN(IDA,IDPA) IDH=MAX(IDA,IDPA) IF(MOD(IDH,2).EQ.0.AND.(IDH-IDL.EQ.1)) THEN IF(ABS(ID).EQ.IDL) THEN IMO=-SIGN(24,ID) ELSE IMO=-SIGN(24,IDP) ENDIF ENDIF ENDIF IF(IMO.EQ.0) GOTO 998 MUP=NUP NUP=NUP+1 DO IC=NUP,IJ+1,-1 IDUP(IC)=IDUP(IC-1) ISTUP(IC)=ISTUP(IC-1) MOTHUP(1,IC)=MOTHUP(1,IC-1) MOTHUP(2,IC)=MOTHUP(2,IC-1) ICOLUP(1,IC)=ICOLUP(1,IC-1) ICOLUP(2,IC)=ICOLUP(2,IC-1) DO J=1,5 PUP(J,IC)=PUP(J,IC-1) ENDDO ENDDO C....... IDUP(IJ)=IMO ISTUP(IJ)=2 DO J=1,4 PUP(J,IJ)=PUP(J,IJ+1)+PUP(J,IJ+2) ENDDO MOTHUP(1,IJ+1)=IJ MOTHUP(2,IJ+1)=IJ MOTHUP(1,IJ+2)=IJ MOTHUP(2,IJ+2)=IJ TEST=PUP(1,IJ)**2+PUP(2,IJ)**2+PUP(3,IJ)**2-PUP(4,IJ)**2 IF(TEST.LE.0D0) THEN PUP(5,IJ)=DSQRT(-TEST) ELSEIF(DABS(TEST).LT.1D-4) THEN PUP(5,IJ)=0D0 ELSE PUP(5,IJ)=-1D0 PRINT*,' NEGATIVE MASS ' ENDIF IJ=IJ+2 ENDIF ENDIF 998 CONTINUE IJ=IJ+1 IF(IJ.GT.NUP) GOTO 1000 GOTO 555 C...Some other compulsory quantities. 1000 SCALUP=-1D0 C... DO J=1,4 QPAIR(J)=PUP(J,1)+PUP(J,2) ENDDO IF(ABS(IDUP(1)).LE.6.and.ABS(IDUP(2)).LE.6) THEN IF(IDUP(1).GT.0.AND.IDUP(2).LT.0) THEN if(icolup(1,1).eq.icolup(2,2)) then scalup=sqrt(qpair(4)**2-qpair(3)**2-qpair(2)**2 $ -qpair(1)**2) endif elseif(idup(1).lt.0.and.idup(2).gt.0) then if(icolup(2,1).eq.icolup(1,2)) then scalup=sqrt(qpair(4)**2-qpair(3)**2-qpair(2)**2 $ -qpair(1)**2) endif endif endif if(scalup.lt.0D0) then CCCCCC....call to ktclustering subroutine NN=0 do 301 i=3,nup if(istup(i).ne.1) goto 301 if(abs(idup(i)).eq.11.or.abs(idup(i)).eq.12) goto 301 if(abs(idup(i)).eq.13.or.abs(idup(i)).eq.14) goto 301 if(abs(idup(i)).eq.15.or.abs(idup(i)).eq.16) goto 301 if(icolup(1,i).eq.0.and.icolup(2,i).eq.0) goto 301 NN=NN+1 DO J=1,4 PP(J,NN)=PUP(J,I) ENDDO 301 continue IF(NN.NE.0) THEN ECUT=1D0 CALL KTCLUS(4223,PP,NN,ECUT,Y,*999) scalup=sqrt(y(NN)) ENDIF 999 continue endif C... C.....decay taus with tauola DO 10 I=1,NUP id=idup(i) if(abs(id).ne.15) goto 10 istup(i)=2 C...Have TAUOLA decay the tau. DO J=1,4 QPAIR(J)=PUP(J,1)+PUP(J,2) PPAIR(J)=-QPAIR(J) PIN(J)=PUP(J,I) ENDDO PPAIR(4)=QPAIR(4) c print*,sqrt(ppair(4)**2-ppair(1)**2-ppair(2)**2-ppair(3)**2) c print*,' pin ',pin call boost_5(pin,ppair,pout) c print*,' pout ',pout do j=1,4 pout4(j)=sngl(pout(j)) enddo pmag=sqrt(pout4(1)**2+pout4(2)**2+pout4(3)**2+1D-10)*spinup(i) do j=1,3 pol(j)=pout4(j)/pmag enddo pol(4)=0D0 amtau=pup(5,i) if (id.eq.15) then ind=NUP call filhep(ind,1,15,0,0,0,0,pout4,amtau,.true.) call dexay(2,pol) idin=4 else if (id.eq.-15) then pol(1) = -pol(1) pol(2) = -pol(2) pol(3) = -pol(3) ind = NUP call filhep(ind,1,-15,0,0,0,0,pout4,amtau,.true.) call dexay(1,pol) idin=3 endif C........add in tau decay products do j=nup+1,nhep nup=nup+1 do k=1,4 v1(k)=phep(k,j) v2(k)=0D0 enddo call boost_5(v1,pout,v2) call boost_5(v2,qpair,v1) do k=1,4 pup(k,nup)=v1(k) enddo TEST=PUP(1,nup)**2+PUP(2,nup)**2+PUP(3,nup)**2 $ -PUP(4,nup)**2 IF(TEST.LE.0D0) THEN PUP(5,nup)=DSQRT(-TEST) ELSEIF(DABS(TEST).LT.1D-6) THEN PUP(5,nup)=0D0 ENDIF do k=1,2 mothup(k,nup)=jmohep(k,j) icolup(k,nup)=0 if(mothup(k,nup).eq.idin) then mothup(k,nup)=i endif enddo c vtimup(nup)=sqrt(vhep(4,k)**2-vhep(1,k)**2 $ -vhep(2,k)**2-vhep(3,k)**2) idup(nup)=idhep(j) istup(nup)=isthep(j) enddo 10 continue RETURN 888 CONTINUE REWIND(77) GOTO 99 NUP=0 RETURN END SUBROUTINE BOOST_5(P,R,Q) IMPLICIT DOUBLE PRECISION (A-H,J-Z) DIMENSION P(5),R(5),Q(5),BETA(3) X = 0D0 Y = 0D0 DO 10 I = 1,3 BETA(I) = R(I)/R(4) X = X + BETA(I)**2 10 Y = Y + BETA(I)*P(I) IF (X.LT.1D-16.OR.X.GE.(1D0-1D-12)) GOTO 30 GAMMA = 1D0/DSQRT(1D0-X) DO 20 I = 1,3 20 Q(I) = P(I)+BETA(I)*(Y*(GAMMA-1D0)/X + GAMMA*P(4)) Q(4) = GAMMA*(P(4) + Y) RETURN 30 CONTINUE DO 40 I = 1,4 40 Q(I) = P(I) IF(X.GE.(1D0-1D-12)) WRITE(*,1000) R RETURN 1000 FORMAT (' THE REFERENCE VECTOR ',4F10.3,' IS NOT TIMELIKE.') END C SUBROUTINE BOSTD3(EXE,PVEC,QVEC) IMPLICIT NONE DOUBLE PRECISION EXE,PVEC(4),QVEC(4) C ---------------------------------------------------------------------- C BOOST ALONG Z AXIS, EXE=EXP(ETA), ETA= HIPERBOLIC VELOCITY. C C USED BY : KORALZ RADKOR C ---------------------------------------------------------------------- INTEGER I DOUBLE PRECISION RVEC(4) DOUBLE PRECISION RPL,QPL,RMI,QMI C DO 10 I=1,4 10 RVEC(I)=PVEC(I) RPL=RVEC(4)+RVEC(3) RMI=RVEC(4)-RVEC(3) QPL=RPL*EXE QMI=RMI/EXE QVEC(1)=RVEC(1) QVEC(2)=RVEC(2) QVEC(3)=(QPL-QMI)/2 QVEC(4)=(QPL+QMI)/2 RETURN END SUBROUTINE BOSTR3(EXE,PVEC,QVEC) IMPLICIT NONE REAL EXE,PVEC(4),QVEC(4) C ---------------------------------------------------------------------- C BOOST ALONG Z AXIS, EXE=EXP(ETA), ETA= HIPERBOLIC VELOCITY. C C USED BY : TAUOLA KORALZ (?) C ---------------------------------------------------------------------- INTEGER I REAL RPL,QPL,RMI,QMI REAL RVEC(4) C DO 10 I=1,4 10 RVEC(I)=PVEC(I) RPL=RVEC(4)+RVEC(3) RMI=RVEC(4)-RVEC(3) QPL=RPL*EXE QMI=RMI/EXE QVEC(1)=RVEC(1) QVEC(2)=RVEC(2) QVEC(3)=(QPL-QMI)/2 QVEC(4)=(QPL+QMI)/2 END COMPLEX FUNCTION BWIG(S,M,G) IMPLICIT NONE C ********************************************************** C P-WAVE BREIT-WIGNER FOR RHO C ********************************************************** REAL S,M,G REAL QS,QM,W,GS REAL PI,PIM PARAMETER (PI=3.141592654) PARAMETER (PIM=0.139) C ------- BREIT-WIGNER ----------------------- IF (S.GT.4.*PIM**2) THEN QS=SQRT(ABS(ABS(S/4.-PIM**2)+(S/4.-PIM**2))/2.0) QM=SQRT(M**2/4.-PIM**2) W=SQRT(S) GS=G*(M/W)*(QS/QM)**3 ELSE GS=0.0 ENDIF BWIG=M**2/CMPLX(M**2-S,-M*GS) RETURN END COMPLEX FUNCTION BWIGS(S,M,G) IMPLICIT NONE C ********************************************************** C P-WAVE BREIT-WIGNER FOR K* C ********************************************************** REAL S,M,G REAL QS,QM,W,GS REAL PI,PIM,MK PARAMETER (PI=3.141592654) PARAMETER (PIM=0.139) PARAMETER (MK=0.493667) REAL P,A,B,C P(A,B,C)=SQRT(ABS(ABS(((A+B-C)**2-4.*A*B)/4./A) $ +(((A+B-C)**2-4.*A*B)/4./A))/2.0) C ------- BREIT-WIGNER ----------------------- QS=P(S,PIM**2,MK**2) QM=P(M**2,PIM**2,MK**2) W=SQRT(S) GS=G*(M/W)*(QS/QM)**3 BWIGS=M**2/CMPLX(M**2-S,-M*GS) RETURN END SUBROUTINE CLAXI(HJ,PN,PIA) IMPLICIT NONE INCLUDE 'tauola.inc' COMPLEX HJ(4) REAL PN(4),PIA(4) C ---------------------------------------------------------------------- * CALCULATES THE "AXIAL TYPE" PI-VECTOR PIA * NOTE THAT THE NEUTRINO MOM. PN IS ASSUMED TO BE ALONG Z-AXIS C SIGN is chosen +/- for decay of TAU +/- respectively C called by : DAMPAA, CLNUT C ---------------------------------------------------------------------- INTEGER I,J REAL SIGN COMPLEX HJC(4) C DET2(I,J)=AIMAG(HJ(I)*HJC(J)-HJ(J)*HJC(I)) C -- here was an error (ZW, 21.11.1991) REAL DET2 DET2(I,J)=AIMAG(HJC(I)*HJ(J)-HJC(J)*HJ(I)) C -- it was affecting sign of A_LR asymmetry in a1 decay. C -- note also collision of notation of gamma_va as defined in C -- TAUOLA paper and J.H. Kuhn and Santamaria Z. Phys C 48 (1990) 445 * ----------------------------------- IF (KTOM.EQ.1.OR.KTOM.EQ.-1) THEN SIGN= IDFF/ABS(IDFF) ELSEIF (KTOM.EQ.2) THEN SIGN=-IDFF/ABS(IDFF) ELSE PRINT *, 'STOP IN CLAXI: KTOM=',KTOM STOP ENDIF C DO 10 I=1,4 10 HJC(I)=CONJG(HJ(I)) PIA(1)= -2.*PN(3)*DET2(2,4)+2.*PN(4)*DET2(2,3) PIA(2)= -2.*PN(4)*DET2(1,3)+2.*PN(3)*DET2(1,4) PIA(3)= 2.*PN(4)*DET2(1,2) PIA(4)= 2.*PN(3)*DET2(1,2) C ALL FOUR INDICES ARE UP SO PIA(3) AND PIA(4) HAVE SAME SIGN DO 20 I=1,4 20 PIA(I)=PIA(I)*SIGN END SUBROUTINE CLNUT(HJ,B,HV) IMPLICIT NONE COMPLEX HJ(4) REAL B,HV(4) C ---------------------------------------------------------------------- * CALCULATES THE CONTRIBUTION BY NEUTRINO MASS * NOTE THE TAU IS ASSUMED TO BE AT REST C C called by : DAMPAA C ---------------------------------------------------------------------- REAL P(4) DATA P /3*0.,1.0/ C CALL CLAXI(HJ,P,HV) B=REAL( HJ(4)*AIMAG(HJ(4)) - HJ(3)*AIMAG(HJ(3)) & - HJ(2)*AIMAG(HJ(2)) - HJ(1)*AIMAG(HJ(1)) ) RETURN END SUBROUTINE CLVEC(HJ,PN,PIV) IMPLICIT NONE COMPLEX HJ(4) REAL PN(4),PIV(4) C ---------------------------------------------------------------------- * CALCULATES THE "VECTOR TYPE" PI-VECTOR PIV * NOTE THAT THE NEUTRINO MOM. PN IS ASSUMED TO BE ALONG Z-AXIS C C called by : DAMPAA C ---------------------------------------------------------------------- INTEGER I REAL HH COMPLEX HN C HN= HJ(4)*CMPLX(PN(4))-HJ(3)*CMPLX(PN(3)) HH= REAL(HJ(4)*CONJG(HJ(4))-HJ(3)*CONJG(HJ(3)) $ -HJ(2)*CONJG(HJ(2))-HJ(1)*CONJG(HJ(1))) DO 10 I=1,4 10 PIV(I)=4.*REAL(HN*CONJG(HJ(I)))-2.*HH*PN(I) RETURN END SUBROUTINE DADMAA(MODE,ISGN,HHV,PNU,PAA,PIM1,PIM2,PIPL,JAA) IMPLICIT NONE INCLUDE 'tauola.inc' INTEGER MODE,ISGN REAL HHV(4),PNU(4),PAA(4),PIM1(4),PIM2(4),PIPL(4) INTEGER JAA C ---------------------------------------------------------------------- * A1 DECAY UNWEIGHTED EVENTS C ---------------------------------------------------------------------- INTEGER I REAL HV(4) REAL PDUM1(4),PDUM2(4),PDUM3(4),PDUM4(4),PDUM5(4) REAL RRR(3) DOUBLE PRECISION SWT, SSWT REAL PI /3.141592653589793238462643/ REAL RN,COSTHE,THET,PHI,PARGAM,ERROR,RAT INTEGER IWARM/0/ INTEGER NEVRAW,NEVACC,NEVOVR REAL WT,WTMAX SAVE IWARM SAVE NEVRAW,NEVACC,NEVOVR SAVE SWT,SSWT,WT,WTMAX C IF(MODE.EQ.-1) THEN C =================== IWARM=1 NEVRAW=0 NEVACC=0 NEVOVR=0 SWT=0 SSWT=0 WTMAX=1E-20 DO 15 I=1,500 CALL DPHSAA(WT,HV,PDUM1,PDUM2,PDUM3,PDUM4,PDUM5,JAA) IF(WT.GT.WTMAX/1.2) WTMAX=WT*1.2 15 CONTINUE CC CALL HBOOK1(801,'WEIGHT DISTRIBUTION DADMAA $',100,0,2) C ELSEIF(MODE.EQ. 0) THEN C ======================= 300 CONTINUE IF(IWARM.EQ.0) GOTO 902 CALL DPHSAA(WT,HV,PNU,PAA,PIM1,PIM2,PIPL,JAA) CC CALL HFILL(801,WT/WTMAX) NEVRAW=NEVRAW+1 SWT=SWT+WT SSWT=SSWT+WT**2 CALL TDRAND(RRR,3) RN=RRR(1) IF(WT.GT.WTMAX) NEVOVR=NEVOVR+1 IF(RN*WTMAX.GT.WT) GOTO 300 C ROTATIONS TO BASIC TAU REST FRAME COSTHE=-1.+2.*RRR(2) THET=ACOS(COSTHE) PHI =2*PI*RRR(3) CALL ROTPOL(THET,PHI,PNU) CALL ROTPOL(THET,PHI,PAA) CALL ROTPOL(THET,PHI,PIM1) CALL ROTPOL(THET,PHI,PIM2) CALL ROTPOL(THET,PHI,PIPL) CALL ROTPOL(THET,PHI,HV) DO 44 I=1,3 44 HHV(I)=-ISGN*HV(I) NEVACC=NEVACC+1 C ELSEIF(MODE.EQ. 1) THEN C ======================= IF(NEVRAW.EQ.0) RETURN PARGAM=SWT/FLOAT(NEVRAW+1) ERROR=0 IF(NEVRAW.NE.0) ERROR=SQRT(SSWT/SWT**2-1./FLOAT(NEVRAW)) RAT=PARGAM/GAMEL WRITE(IOUT, 7010) NEVRAW,NEVACC,NEVOVR,PARGAM,RAT,ERROR CC CALL HPRINT(801) GAMPMC(5)=RAT GAMPER(5)=ERROR CAM NEVDEC(5)=NEVACC ENDIF C ===== RETURN 7003 FORMAT(///1X,15(5H*****) $ /,' *', 25X,'******** DADMAA INITIALISATION ********',9X,1H* $ /,' *',E20.5,5X,'WTMAX = MAXIMUM WEIGHT ',9X,1H* $ /,1X,15(5H*****)/) 7010 FORMAT(///1X,15(5H*****) $ /,' *', 25X,'******** DADMAA FINAL REPORT ******** ',9X,1H* $ /,' *',I20 ,5X,'NEVRAW = NO. OF A1 DECAYS TOTAL ',9X,1H* $ /,' *',I20 ,5X,'NEVACC = NO. OF A1 DECS. ACCEPTED ',9X,1H* $ /,' *',I20 ,5X,'NEVOVR = NO. OF OVERWEIGHTED EVENTS ',9X,1H* $ /,' *',E20.5,5X,'PARTIAL WTDTH (A1 DECAY) IN GEV UNITS ',9X,1H* $ /,' *',F20.9,5X,'IN UNITS GFERMI**2*MASS**5/192/PI**3 ',9X,1H* $ /,' *',F20.8,5X,'RELATIVE ERROR OF PARTIAL WIDTH ',9X,1H* $ /,1X,15(5H*****)/) 902 WRITE(IOUT, 9020) 9020 FORMAT(' ----- DADMAA: LACK OF INITIALISATION') STOP END SUBROUTINE DADMEL(MODE,ISGN,HHV,PNU,PWB,Q1,Q2,PHX) IMPLICIT NONE INCLUDE 'tauola.inc' INTEGER MODE,ISGN REAL HHV(4),PNU(4),PWB(4),Q1(4),Q2(4),PHX(4) C ---------------------------------------------------------------------- C C called by : DEXEL,(DEKAY,DEKAY1) C ---------------------------------------------------------------------- INTEGER I REAL RN,RR1,RR2,RR3,COSTHE,THET,PHI,PARGAM,ERROR,RAT REAL HV(4) REAL PDUM1(4),PDUM2(4),PDUM3(4),PDUM4(4),PDUM5(4) REAL RRR(3) REAL PI /3.141592653589793238462643/ INTEGER IWARM/0/ INTEGER NEVRAW,NEVACC,NEVOVR DOUBLE PRECISION SWT, SSWT REAL WT,WTMAX SAVE IWARM SAVE NEVRAW,NEVACC,NEVOVR,SWT,SSWT,WT,WTMAX C IF(MODE.EQ.-1) THEN C =================== IWARM=1 NEVRAW=0 NEVACC=0 NEVOVR=0 SWT=0 SSWT=0 WTMAX=1E-20 DO 15 I=1,500 CALL DPHSEL(WT,HV,PDUM1,PDUM2,PDUM3,PDUM4,PDUM5) IF(WT.GT.WTMAX/1.2) WTMAX=WT*1.2 15 CONTINUE CC CALL HBOOK1(803,'WEIGHT DISTRIBUTION DADMEL $',100,0,2) C ELSEIF(MODE.EQ. 0) THEN C ======================= 300 CONTINUE IF(IWARM.EQ.0) GOTO 902 NEVRAW=NEVRAW+1 CALL DPHSEL(WT,HV,PNU,PWB,Q1,Q2,PHX) CC CALL HFILL(803,WT/WTMAX) SWT=SWT+WT SSWT=SSWT+WT**2 CALL TDRAND(RRR,3) RN=RRR(1) IF(WT.GT.WTMAX) NEVOVR=NEVOVR+1 IF(RN*WTMAX.GT.WT) GOTO 300 C ROTATIONS TO BASIC TAU REST FRAME RR2=RRR(2) COSTHE=-1.+2.*RR2 THET=ACOS(COSTHE) RR3=RRR(3) PHI =2*PI*RR3 CALL ROTOR2(THET,PNU,PNU) CALL ROTOR3( PHI,PNU,PNU) CALL ROTOR2(THET,PWB,PWB) CALL ROTOR3( PHI,PWB,PWB) CALL ROTOR2(THET,Q1,Q1) CALL ROTOR3( PHI,Q1,Q1) CALL ROTOR2(THET,Q2,Q2) CALL ROTOR3( PHI,Q2,Q2) CALL ROTOR2(THET,HV,HV) CALL ROTOR3( PHI,HV,HV) CALL ROTOR2(THET,PHX,PHX) CALL ROTOR3( PHI,PHX,PHX) DO 44,I=1,3 44 HHV(I)=-ISGN*HV(I) NEVACC=NEVACC+1 C ELSEIF(MODE.EQ. 1) THEN C ======================= IF(NEVRAW.EQ.0) RETURN PARGAM=SWT/FLOAT(NEVRAW+1) ERROR=0 IF(NEVRAW.NE.0) ERROR=SQRT(SSWT/SWT**2-1./FLOAT(NEVRAW)) RAT=PARGAM/GAMEL WRITE(IOUT, 7010) NEVRAW,NEVACC,NEVOVR,PARGAM,RAT,ERROR CC CALL HPRINT(803) GAMPMC(1)=RAT GAMPER(1)=ERROR CAM NEVDEC(1)=NEVACC ENDIF C ===== RETURN 7010 FORMAT(///1X,15(5H*****) $ /,' *', 25X,'******** DADMEL FINAL REPORT ******** ',9X,1H* $ /,' *',I20 ,5X,'NEVRAW = NO. OF EL DECAYS TOTAL ',9X,1H* $ /,' *',I20 ,5X,'NEVACC = NO. OF EL DECS. ACCEPTED ',9X,1H* $ /,' *',I20 ,5X,'NEVOVR = NO. OF OVERWEIGHTED EVENTS ',9X,1H* $ /,' *',E20.5,5X,'PARTIAL WTDTH ( ELECTRON) IN GEV UNITS ',9X,1H* $ /,' *',F20.9,5X,'IN UNITS GFERMI**2*MASS**5/192/PI**3 ',9X,1H* $ /,' *',F20.9,5X,'RELATIVE ERROR OF PARTIAL WIDTH ',9X,1H* $ /,' *',25X, 'COMPLETE QED CORRECTIONS INCLUDED ',9X,1H* $ /,' *',25X, 'BUT ONLY V-A CUPLINGS ',9X,1H* $ /,1X,15(5H*****)/) 902 WRITE(IOUT, 9020) 9020 FORMAT(' ----- DADMEL: LACK OF INITIALISATION') STOP END SUBROUTINE DADMKK(MODE,ISGN,HV,PKK,PNU) IMPLICIT NONE INCLUDE 'tauola.inc' INTEGER MODE,ISGN REAL HV(4),PKK(4),PNU(4) C ---------------------------------------------------------------------- C FZ INTEGER I REAL EKK,ENU,XKK,PXQ,PXN,QXN,BRAK,FKK,GAMM,ERROR,RAT SAVE BRAK REAL PI /3.141592653589793238462643/ INTEGER NEVTOT SAVE NEVTOT C IF(MODE.EQ.-1) THEN C =================== NEVTOT=0 ELSEIF(MODE.EQ. 0) THEN C ======================= NEVTOT=NEVTOT+1 EKK= (AMTAU**2+AMK**2-AMNUTA**2)/(2*AMTAU) ENU= (AMTAU**2-AMK**2+AMNUTA**2)/(2*AMTAU) XKK= SQRT(EKK**2-AMK**2) C K MOMENTUM CALL SPHERA(XKK,PKK) PKK(4)=EKK C TAU-NEUTRINO MOMENTUM DO 30 I=1,3 30 PNU(I)=-PKK(I) PNU(4)=ENU PXQ=AMTAU*EKK PXN=AMTAU*ENU QXN=PKK(4)*PNU(4)-PKK(1)*PNU(1)-PKK(2)*PNU(2)-PKK(3)*PNU(3) BRAK=(GV**2+GA**2)*(2*PXQ*QXN-AMK**2*PXN) & +(GV**2-GA**2)*AMTAU*AMNUTA*AMK**2 DO 40 I=1,3 40 HV(I)=-ISGN*2*GA*GV*AMTAU*(2*PKK(I)*QXN-PNU(I)*AMK**2)/BRAK HV(4)=1 C ELSEIF(MODE.EQ. 1) THEN C ======================= IF(NEVTOT.EQ.0) RETURN FKK=0.0354 CFZ THERE WAS BRAK/AMTAU**4 BEFORE C GAMM=(GFERMI*FKK)**2/(16.*PI)*AMTAU**3* C * (BRAK/AMTAU**4)**2 CZW 7.02.93 here was an error affecting non standard model C configurations only GAMM=(GFERMI*FKK)**2/(16.*PI)*AMTAU**3* $ (BRAK/AMTAU**4)* $ SQRT((AMTAU**2-AMK**2-AMNUTA**2)**2 $ -4*AMK**2*AMNUTA**2 )/AMTAU**2 ERROR=0 ERROR=0 RAT=GAMM/GAMEL WRITE(IOUT, 7010) NEVTOT,GAMM,RAT,ERROR GAMPMC(6)=RAT GAMPER(6)=ERROR CAM NEVDEC(6)=NEVTOT ENDIF C ===== RETURN 7010 FORMAT(///1X,15(5H*****) $ /,' *', 25X,'******** DADMKK FINAL REPORT ********',9X,1H* $ /,' *',I20 ,5X,'NEVTOT = NO. OF K DECAYS TOTAL ',9X,1H*, $ /,' *',E20.5,5X,'PARTIAL WTDTH ( K DECAY) IN GEV UNITS ',9X,1H*, $ /,' *',F20.9,5X,'IN UNITS GFERMI**2*MASS**5/192/PI**3 ',9X,1H* $ /,' *',F20.8,5X,'RELATIVE ERROR OF PARTIAL WIDTH (STAT.)',9X,1H* $ /,1X,15(5H*****)/) END SUBROUTINE DADMKS(MODE,ISGN,HHV,PNU,PKS,PKK,PPI,JKST) IMPLICIT NONE INCLUDE 'tauola.inc' INTEGER MODE,ISGN REAL HHV(4),PNU(4),PKS(4),PKK(4),PPI(4) INTEGER JKST C ---------------------------------------------------------------------- INTEGER I REAL DEC1,RMOD,RN,COSTHE,THET,PHI,PARGAM,ERROR,RAT REAL HV(4) REAL PDUM1(4),PDUM2(4),PDUM3(4),PDUM4(4) REAL*4 RRR(3) REAL PI /3.141592653589793238462643/ INTEGER IWARM/0/ INTEGER NEVRAW,NEVACC,NEVOVR DOUBLE PRECISION SWT, SSWT REAL WT,WTMAX SAVE IWARM SAVE NEVRAW,NEVACC,NEVOVR,SWT,SSWT,WT,WTMAX C IF(MODE.EQ.-1) THEN C =================== IWARM=1 NEVRAW=0 NEVACC=0 NEVOVR=0 SWT=0 SSWT=0 WTMAX=1E-20 DO 15 I=1,500 C THE INITIALISATION IS DONE WITH THE 66.7% MODE JKST=10 CALL DPHSKS(WT,HV,PDUM1,PDUM2,PDUM3,PDUM4,JKST) IF(WT.GT.WTMAX/1.2) WTMAX=WT*1.2 15 CONTINUE CC CALL HBOOK1(801,'WEIGHT DISTRIBUTION DADMKS $',100,0,2) CC PRINT 7003,WTMAX CC CALL HBOOK1(112,'-------- K* MASS -------- $',100,0.,2.) ELSEIF(MODE.EQ. 0) THEN C ===================================== IF(IWARM.EQ.0) GOTO 902 C HERE WE CHOOSE RANDOMLY BETWEEN K0 PI+_ (66.7%) C AND K+_ PI0 (33.3%) DEC1=BRKS 400 CONTINUE CALL TDRAND(RMOD,1) IF(RMOD.LT.DEC1) THEN JKST=10 ELSE JKST=20 ENDIF CALL DPHSKS(WT,HV,PNU,PKS,PKK,PPI,JKST) CALL TDRAND(RRR,3) RN=RRR(1) IF(WT.GT.WTMAX) NEVOVR=NEVOVR+1 NEVRAW=NEVRAW+1 SWT=SWT+WT SSWT=SSWT+WT**2 IF(RN*WTMAX.GT.WT) GOTO 400 C ROTATIONS TO BASIC TAU REST FRAME COSTHE=-1.+2.*RRR(2) THET=ACOS(COSTHE) PHI =2*PI*RRR(3) CALL ROTOR2(THET,PNU,PNU) CALL ROTOR3( PHI,PNU,PNU) CALL ROTOR2(THET,PKS,PKS) CALL ROTOR3( PHI,PKS,PKS) CALL ROTOR2(THET,PKK,PKK) CALL ROTOR3(PHI,PKK,PKK) CALL ROTOR2(THET,PPI,PPI) CALL ROTOR3( PHI,PPI,PPI) CALL ROTOR2(THET,HV,HV) CALL ROTOR3( PHI,HV,HV) DO 44 I=1,3 44 HHV(I)=-ISGN*HV(I) NEVACC=NEVACC+1 C ELSEIF(MODE.EQ. 1) THEN C ======================= IF(NEVRAW.EQ.0) RETURN PARGAM=SWT/FLOAT(NEVRAW+1) ERROR=0 IF(NEVRAW.NE.0) ERROR=SQRT(SSWT/SWT**2-1./FLOAT(NEVRAW)) RAT=PARGAM/GAMEL WRITE(IOUT, 7010) NEVRAW,NEVACC,NEVOVR,PARGAM,RAT,ERROR CC CALL HPRINT(801) GAMPMC(7)=RAT GAMPER(7)=ERROR CAM NEVDEC(7)=NEVACC ENDIF C ===== RETURN 7003 FORMAT(///1X,15(5H*****) $ /,' *', 25X,'******** DADMKS INITIALISATION ********',9X,1H* $ /,' *',E20.5,5X,'WTMAX = MAXIMUM WEIGHT ',9X,1H* $ /,1X,15(5H*****)/) 7010 FORMAT(///1X,15(5H*****) $ /,' *', 25X,'******** DADMKS FINAL REPORT ********',9X,1H* $ /,' *',I20 ,5X,'NEVRAW = NO. OF K* DECAYS TOTAL ',9X,1H*, $ /,' *',I20 ,5X,'NEVACC = NO. OF K* DECS. ACCEPTED ',9X,1H*, $ /,' *',I20 ,5X,'NEVOVR = NO. OF OVERWEIGHTED EVENTS ',9X,1H* $ /,' *',E20.5,5X,'PARTIAL WTDTH (K* DECAY) IN GEV UNITS ',9X,1H*, $ /,' *',F20.9,5X,'IN UNITS GFERMI**2*MASS**5/192/PI**3 ',9X,1H* $ /,' *',F20.8,5X,'RELATIVE ERROR OF PARTIAL WIDTH ',9X,1H* $ /,1X,15(5H*****)/) 902 WRITE(IOUT, 9020) 9020 FORMAT(' ----- DADMKS: LACK OF INITIALISATION') STOP END SUBROUTINE DADMMU(MODE,ISGN,HHV,PNU,PWB,Q1,Q2,PHX) IMPLICIT NONE INCLUDE 'tauola.inc' INTEGER MODE,ISGN REAL HHV(4),PNU(4),PWB(4),Q1(4),Q2(4),PHX(4) C ---------------------------------------------------------------------- INTEGER I REAL RN,COSTHE,THET,PHI,PARGAM,ERROR,RAT REAL HV(4) REAL PDUM1(4),PDUM2(4),PDUM3(4),PDUM4(4),PDUM5(4) REAL RRR(3) REAL PI /3.141592653589793238462643/ INTEGER IWARM /0/ INTEGER NEVRAW,NEVACC,NEVOVR DOUBLE PRECISION SWT, SSWT REAL WT,WTMAX SAVE IWARM SAVE NEVRAW,NEVACC,NEVOVR,SWT,SSWT,WT,WTMAX C IF(MODE.EQ.-1) THEN C =================== IWARM=1 NEVRAW=0 NEVACC=0 NEVOVR=0 SWT=0 SSWT=0 WTMAX=1E-20 DO 15 I=1,500 CALL DPHSMU(WT,HV,PDUM1,PDUM2,PDUM3,PDUM4,PDUM5) IF(WT.GT.WTMAX/1.2) WTMAX=WT*1.2 15 CONTINUE CC CALL HBOOK1(802,'WEIGHT DISTRIBUTION DADMMU $',100,0,2) C ELSEIF(MODE.EQ. 0) THEN C ======================= 300 CONTINUE IF(IWARM.EQ.0) GOTO 902 NEVRAW=NEVRAW+1 CALL DPHSMU(WT,HV,PNU,PWB,Q1,Q2,PHX) CC CALL HFILL(802,WT/WTMAX) SWT=SWT+WT SSWT=SSWT+WT**2 CALL TDRAND(RRR,3) RN=RRR(1) IF(WT.GT.WTMAX) NEVOVR=NEVOVR+1 IF(RN*WTMAX.GT.WT) GOTO 300 C ROTATIONS TO BASIC TAU REST FRAME COSTHE=-1.+2.*RRR(2) THET=ACOS(COSTHE) PHI =2*PI*RRR(3) CALL ROTOR2(THET,PNU,PNU) CALL ROTOR3( PHI,PNU,PNU) CALL ROTOR2(THET,PWB,PWB) CALL ROTOR3( PHI,PWB,PWB) CALL ROTOR2(THET,Q1,Q1) CALL ROTOR3( PHI,Q1,Q1) CALL ROTOR2(THET,Q2,Q2) CALL ROTOR3( PHI,Q2,Q2) CALL ROTOR2(THET,HV,HV) CALL ROTOR3( PHI,HV,HV) CALL ROTOR2(THET,PHX,PHX) CALL ROTOR3( PHI,PHX,PHX) DO 44,I=1,3 44 HHV(I)=-ISGN*HV(I) NEVACC=NEVACC+1 C ELSEIF(MODE.EQ. 1) THEN C ======================= IF(NEVRAW.EQ.0) RETURN PARGAM=SWT/FLOAT(NEVRAW+1) ERROR=0 IF(NEVRAW.NE.0) ERROR=SQRT(SSWT/SWT**2-1./FLOAT(NEVRAW)) RAT=PARGAM/GAMEL WRITE(IOUT, 7010) NEVRAW,NEVACC,NEVOVR,PARGAM,RAT,ERROR CC CALL HPRINT(802) GAMPMC(2)=RAT GAMPER(2)=ERROR CAM NEVDEC(2)=NEVACC ENDIF C ===== RETURN 7010 FORMAT(///1X,15(5H*****) $ /,' *', 25X,'******** DADMMU FINAL REPORT ******** ',9X,1H* $ /,' *',I20 ,5X,'NEVRAW = NO. OF MU DECAYS TOTAL ',9X,1H* $ /,' *',I20 ,5X,'NEVACC = NO. OF MU DECS. ACCEPTED ',9X,1H* $ /,' *',I20 ,5X,'NEVOVR = NO. OF OVERWEIGHTED EVENTS ',9X,1H* $ /,' *',E20.5,5X,'PARTIAL WTDTH (MU DECAY) IN GEV UNITS ',9X,1H* $ /,' *',F20.9,5X,'IN UNITS GFERMI**2*MASS**5/192/PI**3 ',9X,1H* $ /,' *',F20.9,5X,'RELATIVE ERROR OF PARTIAL WIDTH ',9X,1H* $ /,' *',25X, 'COMPLETE QED CORRECTIONS INCLUDED ',9X,1H* $ /,' *',25X, 'BUT ONLY V-A CUPLINGS ',9X,1H* $ /,1X,15(5H*****)/) 902 WRITE(IOUT, 9020) 9020 FORMAT(' ----- DADMMU: LACK OF INITIALISATION') STOP END SUBROUTINE DADMPI(MODE,ISGN,HV,PPI,PNU) IMPLICIT NONE INCLUDE 'tauola.inc' INTEGER MODE,ISGN REAL HHV(4),PPI(4),PNU(4) C ---------------------------------------------------------------------- INTEGER I REAL EPI,ENU,XPI,PXQ,PXN,QXN,BRAK,FPI,GAMM,ERROR,RAT SAVE BRAK REAL HV(4) REAL PI /3.141592653589793238462643/ INTEGER NEVTOT SAVE NEVTOT C IF(MODE.EQ.-1) THEN C =================== NEVTOT=0 ELSEIF(MODE.EQ. 0) THEN C ======================= NEVTOT=NEVTOT+1 EPI= (AMTAU**2+AMPI**2-AMNUTA**2)/(2*AMTAU) ENU= (AMTAU**2-AMPI**2+AMNUTA**2)/(2*AMTAU) XPI= SQRT(EPI**2-AMPI**2) C PI MOMENTUM CALL SPHERA(XPI,PPI) PPI(4)=EPI C TAU-NEUTRINO MOMENTUM DO 30 I=1,3 30 PNU(I)=-PPI(I) PNU(4)=ENU PXQ=AMTAU*EPI PXN=AMTAU*ENU QXN=PPI(4)*PNU(4)-PPI(1)*PNU(1)-PPI(2)*PNU(2)-PPI(3)*PNU(3) BRAK=(GV**2+GA**2)*(2*PXQ*QXN-AMPI**2*PXN) & +(GV**2-GA**2)*AMTAU*AMNUTA*AMPI**2 DO 40 I=1,3 40 HV(I)=-ISGN*2*GA*GV*AMTAU*(2*PPI(I)*QXN-PNU(I)*AMPI**2)/BRAK HV(4)=1 C ELSEIF(MODE.EQ. 1) THEN C ======================= IF(NEVTOT.EQ.0) RETURN FPI=0.1284 C GAMM=(GFERMI*FPI)**2/(16.*PI)*AMTAU**3* C * (BRAK/AMTAU**4)**2 CZW 7.02.93 here was an error affecting non standard model C configurations only GAMM=(GFERMI*FPI)**2/(16.*PI)*AMTAU**3* $ (BRAK/AMTAU**4)* $ SQRT((AMTAU**2-AMPI**2-AMNUTA**2)**2 $ -4*AMPI**2*AMNUTA**2 )/AMTAU**2 ERROR=0 RAT=GAMM/GAMEL WRITE(IOUT, 7010) NEVTOT,GAMM,RAT,ERROR GAMPMC(3)=RAT GAMPER(3)=ERROR CAM NEVDEC(3)=NEVTOT ENDIF C ===== RETURN 7010 FORMAT(///1X,15(5H*****) $ /,' *', 25X,'******** DADMPI FINAL REPORT ******** ',9X,1H* $ /,' *',I20 ,5X,'NEVTOT = NO. OF PI DECAYS TOTAL ',9X,1H* $ /,' *',E20.5,5X,'PARTIAL WTDTH ( PI DECAY) IN GEV UNITS ',9X,1H* $ /,' *',F20.9,5X,'IN UNITS GFERMI**2*MASS**5/192/PI**3 ',9X,1H* $ /,' *',F20.8,5X,'RELATIVE ERROR OF PARTIAL WIDTH (STAT.)',9X,1H* $ /,1X,15(5H*****)/) END SUBROUTINE DADMRO(MODE,ISGN,HHV,PNU,PRO,PIC,PIZ) IMPLICIT NONE INCLUDE 'tauola.inc' INTEGER MODE,ISGN REAL HHV(4),PNU(4),PRO(4),PIC(4),PIZ(4) C ---------------------------------------------------------------------- INTEGER I REAL RN,COSTHE,THET,PHI,PARGAM,ERROR,RAT REAL HV(4) REAL PDUM1(4),PDUM2(4),PDUM3(4),PDUM4(4) REAL*4 RRR(3) REAL PI /3.141592653589793238462643/ INTEGER IWARM/0/ INTEGER NEVRAW,NEVACC,NEVOVR DOUBLE PRECISION SWT, SSWT REAL WT,WTMAX SAVE IWARM SAVE NEVRAW,NEVACC,NEVOVR,SWT,SSWT,WT,WTMAX C IF(MODE.EQ.-1) THEN C =================== IWARM=1 NEVRAW=0 NEVACC=0 NEVOVR=0 SWT=0 SSWT=0 WTMAX=1E-20 DO 15 I=1,500 CALL DPHSRO(WT,HV,PDUM1,PDUM2,PDUM3,PDUM4) IF(WT.GT.WTMAX/1.2) WTMAX=WT*1.2 15 CONTINUE CC CALL HBOOK1(801,'WEIGHT DISTRIBUTION DADMRO $',100,0,2) CC PRINT 7003,WTMAX C ELSEIF(MODE.EQ. 0) THEN C ======================= 300 CONTINUE IF(IWARM.EQ.0) GOTO 902 CALL DPHSRO(WT,HV,PNU,PRO,PIC,PIZ) CC CALL HFILL(801,WT/WTMAX) NEVRAW=NEVRAW+1 SWT=SWT+WT SSWT=SSWT+WT**2 CALL TDRAND(RRR,3) RN=RRR(1) IF(WT.GT.WTMAX) NEVOVR=NEVOVR+1 IF(RN*WTMAX.GT.WT) GOTO 300 C ROTATIONS TO BASIC TAU REST FRAME COSTHE=-1.+2.*RRR(2) THET=ACOS(COSTHE) PHI =2*PI*RRR(3) CALL ROTOR2(THET,PNU,PNU) CALL ROTOR3( PHI,PNU,PNU) CALL ROTOR2(THET,PRO,PRO) CALL ROTOR3( PHI,PRO,PRO) CALL ROTOR2(THET,PIC,PIC) CALL ROTOR3( PHI,PIC,PIC) CALL ROTOR2(THET,PIZ,PIZ) CALL ROTOR3( PHI,PIZ,PIZ) CALL ROTOR2(THET,HV,HV) CALL ROTOR3( PHI,HV,HV) DO 44 I=1,3 44 HHV(I)=-ISGN*HV(I) NEVACC=NEVACC+1 C ELSEIF(MODE.EQ. 1) THEN C ======================= IF(NEVRAW.EQ.0) RETURN PARGAM=SWT/FLOAT(NEVRAW+1) ERROR=0 IF(NEVRAW.NE.0) ERROR=SQRT(SSWT/SWT**2-1./FLOAT(NEVRAW)) RAT=PARGAM/GAMEL WRITE(IOUT, 7010) NEVRAW,NEVACC,NEVOVR,PARGAM,RAT,ERROR CC CALL HPRINT(801) GAMPMC(4)=RAT GAMPER(4)=ERROR CAM NEVDEC(4)=NEVACC ENDIF C ===== RETURN 7003 FORMAT(///1X,15(5H*****) $ /,' *', 25X,'******** DADMRO INITIALISATION ********',9X,1H* $ /,' *',E20.5,5X,'WTMAX = MAXIMUM WEIGHT ',9X,1H* $ /,1X,15(5H*****)/) 7010 FORMAT(///1X,15(5H*****) $ /,' *', 25X,'******** DADMRO FINAL REPORT ******** ',9X,1H* $ /,' *',I20 ,5X,'NEVRAW = NO. OF RHO DECAYS TOTAL ',9X,1H* $ /,' *',I20 ,5X,'NEVACC = NO. OF RHO DECS. ACCEPTED ',9X,1H* $ /,' *',I20 ,5X,'NEVOVR = NO. OF OVERWEIGHTED EVENTS ',9X,1H* $ /,' *',E20.5,5X,'PARTIAL WTDTH (RHO DECAY) IN GEV UNITS ',9X,1H* $ /,' *',F20.9,5X,'IN UNITS GFERMI**2*MASS**5/192/PI**3 ',9X,1H* $ /,' *',F20.8,5X,'RELATIVE ERROR OF PARTIAL WIDTH ',9X,1H* $ /,1X,15(5H*****)/) 902 WRITE(IOUT, 9020) 9020 FORMAT(' ----- DADMRO: LACK OF INITIALISATION') STOP END SUBROUTINE DADNEW(MODE,ISGN,HV,PNU,PWB,PNPI,JNPI) IMPLICIT NONE INCLUDE 'tauola.inc' INTEGER MODE,ISGN REAL HV(4),PNU(4),PWB(4),PNPI(4,9) INTEGER JNPI C ---------------------------------------------------------------------- INTEGER I,NMOD,INUM,ND SAVE NMOD REAL RN,COSTHE,THET,PHI,PARGAM,ERROR,RAT REAL HHV(4) REAL PDUM1(4),PDUM2(4),PDUMI(4,9) REAL RRR(3) C REAL PI /3.141592653589793238462643/ INTEGER IWARM/0/ INTEGER NEVRAW(NMODE),NEVOVR(NMODE),NEVACC(NMODE) DOUBLE PRECISION SWT(NMODE),SSWT(NMODE) REAL WT,WTMAX(NMODE) SAVE IWARM SAVE NEVRAW,NEVOVR,NEVACC,SWT,SSWT,WT,WTMAX C IF(MODE.EQ.-1) THEN C =================== C -- AT THE MOMENT ONLY TWO DECAY MODES OF MULTIPIONS HAVE M. ELEM NMOD=NMODE IWARM=1 C PRINT 7003 DO 1 JNPI=1,NMOD NEVRAW(JNPI)=0 NEVACC(JNPI)=0 NEVOVR(JNPI)=0 SWT(JNPI)=0 SSWT(JNPI)=0 WTMAX(JNPI)=-1. DO I=1,500 IF (JNPI.LE.0) THEN GOTO 903 ELSEIF(JNPI.LE.NM4) THEN CALL DPH4PI(WT,HV,PDUM1,PDUM2,PDUMI,JNPI) ELSEIF(JNPI.LE.NM4+NM5) THEN CALL DPH5PI(WT,HV,PDUM1,PDUM2,PDUMI,JNPI) ELSEIF(JNPI.LE.NM4+NM5+NM6) THEN CALL DPHNPI(WT,HV,PDUM1,PDUM2,PDUMI,JNPI) ELSEIF(JNPI.LE.NM4+NM5+NM6+NM3) THEN INUM=JNPI-NM4-NM5-NM6 CALL DPHSPK(WT,HV,PDUM1,PDUM2,PDUMI,INUM) ELSEIF(JNPI.LE.NM4+NM5+NM6+NM3+NM2) THEN INUM=JNPI-NM4-NM5-NM6-NM3 CALL DPHSRK(WT,HV,PDUM1,PDUM2,PDUMI,INUM) ELSE GOTO 903 ENDIF IF(WT.GT.WTMAX(JNPI)/1.2) WTMAX(JNPI)=WT*1.2 ENDDO C CALL HBOOK1(801,'WEIGHT DISTRIBUTION DADNPI $',100,0.,2.,.0) C PRINT 7004,WTMAX(JNPI) 1 CONTINUE WRITE(IOUT,7005) C ELSEIF(MODE.EQ. 0) THEN C ======================= IF(IWARM.EQ.0) GOTO 902 C 300 CONTINUE IF (JNPI.LE.0) THEN GOTO 903 ELSEIF(JNPI.LE.NM4) THEN CALL DPH4PI(WT,HHV,PNU,PWB,PNPI,JNPI) ELSEIF(JNPI.LE.NM4+NM5) THEN CALL DPH5PI(WT,HHV,PNU,PWB,PNPI,JNPI) ELSEIF(JNPI.LE.NM4+NM5+NM6) THEN CALL DPHNPI(WT,HHV,PNU,PWB,PNPI,JNPI) ELSEIF(JNPI.LE.NM4+NM5+NM6+NM3) THEN INUM=JNPI-NM4-NM5-NM6 CALL DPHSPK(WT,HHV,PNU,PWB,PNPI,INUM) ELSEIF(JNPI.LE.NM4+NM5+NM6+NM3+NM2) THEN INUM=JNPI-NM4-NM5-NM6-NM3 CALL DPHSRK(WT,HHV,PNU,PWB,PNPI,INUM) ELSE GOTO 903 ENDIF DO I=1,4 HV(I)=-ISGN*HHV(I) ENDDO C CALL HFILL(801,WT/WTMAX(JNPI)) NEVRAW(JNPI)=NEVRAW(JNPI)+1 SWT(JNPI)=SWT(JNPI)+WT SSWT(JNPI)=SSWT(JNPI)+WT**2 CALL TDRAND(RRR,3) RN=RRR(1) IF(WT.GT.WTMAX(JNPI)) NEVOVR(JNPI)=NEVOVR(JNPI)+1 IF(RN*WTMAX(JNPI).GT.WT) GOTO 300 C ROTATIONS TO BASIC TAU REST FRAME COSTHE=-1.+2.*RRR(2) THET=ACOS(COSTHE) PHI =2*PI*RRR(3) CALL ROTOR2(THET,PNU,PNU) CALL ROTOR3( PHI,PNU,PNU) CALL ROTOR2(THET,PWB,PWB) CALL ROTOR3( PHI,PWB,PWB) CALL ROTOR2(THET,HV,HV) CALL ROTOR3( PHI,HV,HV) ND=MULPIK(JNPI) DO 301 I=1,ND CALL ROTOR2(THET,PNPI(1,I),PNPI(1,I)) CALL ROTOR3( PHI,PNPI(1,I),PNPI(1,I)) 301 CONTINUE NEVACC(JNPI)=NEVACC(JNPI)+1 C ELSEIF(MODE.EQ. 1) THEN C ======================= DO 500 JNPI=1,NMOD IF(NEVRAW(JNPI).EQ.0) GOTO 500 PARGAM=SWT(JNPI)/FLOAT(NEVRAW(JNPI)+1) C...These lines were modified to avoid a square root of a negative number C...exception. This happens because of roundoff errors when there was C...only one event generated. Cal Loomis 7/14/94 c ERROR=0 c IF(NEVRAW(JNPI).NE.0) c & ERROR=SQRT(SSWT(JNPI)/SWT(JNPI)**2-1./FLOAT(NEVRAW(JNPI))) ERROR=0 IF(NEVRAW(JNPI).NE.0) & ERROR=SSWT(JNPI)/SWT(JNPI)**2-1./FLOAT(NEVRAW(JNPI)) error = sqrt(max(0.0,error)) C...End of modification. RAT=PARGAM/GAMEL WRITE(IOUT, 7010) NAMES(JNPI), & NEVRAW(JNPI),NEVACC(JNPI),NEVOVR(JNPI),PARGAM,RAT,ERROR CC CALL HPRINT(801) GAMPMC(8+JNPI-1)=RAT GAMPER(8+JNPI-1)=ERROR CAM NEVDEC(8+JNPI-1)=NEVACC(JNPI) 500 CONTINUE ENDIF C ===== RETURN 7003 FORMAT(///1X,15(5H*****) $ /,' *', 25X,'******** DADNEW INITIALISATION ********',9X,1H* $ ) 7004 FORMAT(' *',E20.5,5X,'WTMAX = MAXIMUM WEIGHT ',9X,1H*/) 7005 FORMAT( $ /,1X,15(5H*****)/) 7010 FORMAT(///1X,15(5H*****) $ /,' *', 25X,'******** DADNEW FINAL REPORT ******** ',9X,1H* $ /,' *', 25X,'CHANNEL:',A31 ,9X,1H* $ /,' *',I20 ,5X,'NEVRAW = NO. OF DECAYS TOTAL ',9X,1H* $ /,' *',I20 ,5X,'NEVACC = NO. OF DECAYS ACCEPTED ',9X,1H* $ /,' *',I20 ,5X,'NEVOVR = NO. OF OVERWEIGHTED EVENTS ',9X,1H* $ /,' *',E20.5,5X,'PARTIAL WTDTH IN GEV UNITS ',9X,1H* $ /,' *',F20.9,5X,'IN UNITS GFERMI**2*MASS**5/192/PI**3 ',9X,1H* $ /,' *',F20.8,5X,'RELATIVE ERROR OF PARTIAL WIDTH ',9X,1H* $ /,1X,15(5H*****)/) 902 WRITE(IOUT, 9020) 9020 FORMAT(' ----- DADNEW: LACK OF INITIALISATION') STOP 903 WRITE(IOUT, 9030) JNPI,MODE 9030 FORMAT(' ----- DADNEW: WRONG JNPI',2I5) STOP END SUBROUTINE DAM4PI(MNUM,PT,PN,PIM1,PIM2,PIM3,PIM4,AMPLIT,HV) IMPLICIT NONE INCLUDE 'tauola.inc' INTEGER MNUM REAL PT(4),PN(4),PIM1(4),PIM2(4),PIM3(4),PIM4(4),AMPLIT,HV(4) C ---------------------------------------------------------------------- * CALCULATES DIFFERENTIAL CROSS SECTION AND POLARIMETER VECTOR * FOR TAU DECAY INTO 4 PI MODES * ALL SPIN EFFECTS IN THE FULL DECAY CHAIN ARE TAKEN INTO ACCOUNT. * CALCULATIONS DONE IN TAU REST FRAME WITH Z-AXIS ALONG NEUTRINO MOMENT C MNUM DECAY MODE IDENTIFIER. C C called by : DPHSAA C ---------------------------------------------------------------------- INTEGER I REAL BRAK,BRAKM REAL PIVEC(4),PIAKS(4),HVM(4) COMPLEX HADCUR(4),FORM1,FORM2,FORM3,FORM4,FORM5 EXTERNAL FORM1,FORM2,FORM3,FORM4,FORM5 REAL PI /3.141592653589793238462643/ INTEGER ICONT /0/ C CALL CURR(MNUM,PIM1,PIM2,PIM3,PIM4,HADCUR) C * CALCULATE PI-VECTORS: VECTOR AND AXIAL CALL CLVEC(HADCUR,PN,PIVEC) CALL CLAXI(HADCUR,PN,PIAKS) CALL CLNUT(HADCUR,BRAKM,HVM) * SPIN INDEPENDENT PART OF DECAY DIFF-CROSS-SECT. IN TAU REST FRAME BRAK= (GV**2+GA**2)*PT(4)*PIVEC(4) +2.*GV*GA*PT(4)*PIAKS(4) & +2.*(GV**2-GA**2)*AMNUTA*AMTAU*BRAKM AMPLIT=(CCABIB*GFERMI)**2*BRAK/2. C POLARIMETER VECTOR IN TAU REST FRAME DO 90 I=1,3 HV(I)=-(AMTAU*((GV**2+GA**2)*PIAKS(I)+2.*GV*GA*PIVEC(I))) & +(GV**2-GA**2)*AMNUTA*AMTAU*HVM(I) C HV IS DEFINED FOR TAU- WITH GAMMA=B+HV*POL IF (BRAK.NE.0.0) &HV(I)=-HV(I)/BRAK 90 CONTINUE END SUBROUTINE DAMPAA(PT,PN,PIM1,PIM2,PIPL,AMPLIT,HV) IMPLICIT NONE INCLUDE 'tauola.inc' REAL PT(4),PN(4),PIM1(4),PIM2(4),PIPL(4),AMPLIT,HV(4) C ---------------------------------------------------------------------- * CALCULATES DIFFERENTIAL CROSS SECTION AND POLARIMETER VECTOR * FOR TAU DECAY INTO A1, A1 DECAYS NEXT INTO RHO+PI AND RHO INTO PI+PI. * ALL SPIN EFFECTS IN THE FULL DECAY CHAIN ARE TAKEN INTO ACCOUNT. * CALCULATIONS DONE IN TAU REST FRAME WITH Z-AXIS ALONG NEUTRINO MOMENT * THE ROUTINE IS WRITEN FOR ZERO NEUTRINO MASS. C C called by : DPHSAA C ---------------------------------------------------------------------- INTEGER I REAL XMAA,XMRO1,XMRO2,PROD1,PROD2,FA1 REAL FAROPI,FRO2PI,FNORM,FAMAX,BRAK,BRAKM REAL GAMAX REAL PAA(4),VEC1(4),VEC2(4) REAL PIVEC(4),PIAKS(4),HVM(4) COMPLEX BWIGN,HADCUR(4),FPIK INTEGER ICONT /1/ REAL GFUN EXTERNAL GFUN C * F CONSTANTS FOR A1, A1-RHO-PI, AND RHO-PI-PI * REAL FPI DATA FPI /93.3E-3/ * THIS INLINE FUNCT. CALCULATES THE SCALAR PART OF THE PROPAGATOR REAL XM,AM,GAMMA BWIGN(XM,AM,GAMMA)=1./CMPLX(XM**2-AM**2,GAMMA*AM) C * FOUR MOMENTUM OF A1 DO 10 I=1,4 10 PAA(I)=PIM1(I)+PIM2(I)+PIPL(I) * MASSES OF A1, AND OF TWO PI-PAIRS WHICH MAY FORM RHO XMAA =SQRT(ABS(PAA(4)**2-PAA(3)**2-PAA(2)**2-PAA(1)**2)) XMRO1 =SQRT(ABS((PIPL(4)+PIM1(4))**2-(PIPL(1)+PIM1(1))**2 $ -(PIPL(2)+PIM1(2))**2-(PIPL(3)+PIM1(3))**2)) XMRO2 =SQRT(ABS((PIPL(4)+PIM2(4))**2-(PIPL(1)+PIM2(1))**2 $ -(PIPL(2)+PIM2(2))**2-(PIPL(3)+PIM2(3))**2)) * ELEMENTS OF HADRON CURRENT PROD1 =PAA(4)*(PIM1(4)-PIPL(4))-PAA(1)*(PIM1(1)-PIPL(1)) $ -PAA(2)*(PIM1(2)-PIPL(2))-PAA(3)*(PIM1(3)-PIPL(3)) PROD2 =PAA(4)*(PIM2(4)-PIPL(4))-PAA(1)*(PIM2(1)-PIPL(1)) $ -PAA(2)*(PIM2(2)-PIPL(2))-PAA(3)*(PIM2(3)-PIPL(3)) DO 40 I=1,4 VEC1(I)= PIM1(I)-PIPL(I) -PAA(I)*PROD1/XMAA**2 40 VEC2(I)= PIM2(I)-PIPL(I) -PAA(I)*PROD2/XMAA**2 * HADRON CURRENT SATURATED WITH A1 AND RHO RESONANCES IF (KEYA1.EQ.1) THEN FA1=9.87 FAROPI=1.0 FRO2PI=1.0 FNORM=FA1/SQRT(2.)*FAROPI*FRO2PI DO 45 I=1,4 HADCUR(I)= CMPLX(FNORM) *AMA1**2*BWIGN(XMAA,AMA1,GAMA1) $ *(CMPLX(VEC1(I))*AMRO**2*BWIGN(XMRO1,AMRO,GAMRO) $ +CMPLX(VEC2(I))*AMRO**2*BWIGN(XMRO2,AMRO,GAMRO)) 45 CONTINUE ELSE FNORM=2.0*SQRT(2.)/3.0/FPI GAMAX=GAMA1*GFUN(XMAA**2)/GFUN(AMA1**2) DO 46 I=1,4 HADCUR(I)= CMPLX(FNORM) *AMA1**2*BWIGN(XMAA,AMA1,GAMAX) $ *(CMPLX(VEC1(I))*FPIK(XMRO1) $ +CMPLX(VEC2(I))*FPIK(XMRO2)) 46 CONTINUE ENDIF C * CALCULATE PI-VECTORS: VECTOR AND AXIAL CALL CLVEC(HADCUR,PN,PIVEC) CALL CLAXI(HADCUR,PN,PIAKS) CALL CLNUT(HADCUR,BRAKM,HVM) * SPIN INDEPENDENT PART OF DECAY DIFF-CROSS-SECT. IN TAU REST FRAME BRAK= (GV**2+GA**2)*PT(4)*PIVEC(4) +2.*GV*GA*PT(4)*PIAKS(4) & +2.*(GV**2-GA**2)*AMNUTA*AMTAU*BRAKM AMPLIT=(GFERMI*CCABIB)**2*BRAK/2. C THE STATISTICAL FACTOR FOR IDENTICAL PI'S WAS CANCELLED WITH C TWO, FOR TWO MODES OF A1 DECAY NAMELLY PI+PI-PI- AND PI-PI0PI0 C POLARIMETER VECTOR IN TAU REST FRAME DO 90 I=1,3 HV(I)=-(AMTAU*((GV**2+GA**2)*PIAKS(I)+2.*GV*GA*PIVEC(I))) & +(GV**2-GA**2)*AMNUTA*AMTAU*HVM(I) C HV IS DEFINED FOR TAU- WITH GAMMA=B+HV*POL HV(I)=-HV(I)/BRAK 90 CONTINUE END SUBROUTINE DAMPOG(PT,PN,PIM1,PIM2,PIPL,AMPLIT,HV) IMPLICIT NONE INCLUDE 'tauola.inc' REAL PT(4),PN(4),PIM1(4),PIM2(4),PIPL(4),AMPLIT,HV(4) C ---------------------------------------------------------------------- * CALCULATES DIFFERENTIAL CROSS SECTION AND POLARIMETER VECTOR * FOR TAU DECAY INTO A1, A1 DECAYS NEXT INTO RHO+PI AND RHO INTO PI+PI. * ALL SPIN EFFECTS IN THE FULL DECAY CHAIN ARE TAKEN INTO ACCOUNT. * CALCULATIONS DONE IN TAU REST FRAME WITH Z-AXIS ALONG NEUTRINO MOMENT * THE ROUTINE IS WRITEN FOR ZERO NEUTRINO MASS. C C called by : DPHSAA C ---------------------------------------------------------------------- INTEGER I,JJ REAL XMAA,XMOM,XMRO2,PROD1,PROD2,P12,P1PL,P2PL,GNORM REAL P1VEC1,P1VEC2,P2VEC1,P2VEC2,BRAK,BRAKM REAL PAA(4),VEC1(4),VEC2(4) REAL PIVEC(4),PIAKS(4),HVM(4) COMPLEX HADCUR(4),FNORM,FORMOM INTEGER ICONT /1/ C * FOUR MOMENTUM OF A1 DO 10 I=1,4 VEC1(I)=0.0 VEC2(I)=0.0 HV(I) =0.0 10 PAA(I)=PIM1(I)+PIM2(I)+PIPL(I) VEC1(1)=1.0 * MASSES OF A1, AND OF TWO PI-PAIRS WHICH MAY FORM RHO XMAA =SQRT(ABS(PAA(4)**2-PAA(3)**2-PAA(2)**2-PAA(1)**2)) XMOM =SQRT(ABS( (PIM2(4)+PIPL(4))**2-(PIM2(3)+PIPL(3))**2 $ -(PIM2(2)+PIPL(2))**2-(PIM2(1)+PIPL(1))**2 )) XMRO2 =(PIPL(1))**2 +(PIPL(2))**2 +(PIPL(3))**2 * ELEMENTS OF HADRON CURRENT PROD1 =VEC1(1)*PIPL(1) PROD2 =VEC2(2)*PIPL(2) P12 =PIM1(4)*PIM2(4)-PIM1(1)*PIM2(1) $ -PIM1(2)*PIM2(2)-PIM1(3)*PIM2(3) P1PL =PIM1(4)*PIPL(4)-PIM1(1)*PIPL(1) $ -PIM1(2)*PIPL(2)-PIM1(3)*PIPL(3) P2PL =PIPL(4)*PIM2(4)-PIPL(1)*PIM2(1) $ -PIPL(2)*PIM2(2)-PIPL(3)*PIM2(3) DO 40 I=1,3 VEC1(I)= (VEC1(I)-PROD1/XMRO2*PIPL(I)) 40 CONTINUE GNORM=SQRT(VEC1(1)**2+VEC1(2)**2+VEC1(3)**2) DO 41 I=1,3 VEC1(I)= VEC1(I)/GNORM 41 CONTINUE VEC2(1)=(VEC1(2)*PIPL(3)-VEC1(3)*PIPL(2))/SQRT(XMRO2) VEC2(2)=(VEC1(3)*PIPL(1)-VEC1(1)*PIPL(3))/SQRT(XMRO2) VEC2(3)=(VEC1(1)*PIPL(2)-VEC1(2)*PIPL(1))/SQRT(XMRO2) P1VEC1 =PIM1(4)*VEC1(4)-PIM1(1)*VEC1(1) $ -PIM1(2)*VEC1(2)-PIM1(3)*VEC1(3) P2VEC1 =VEC1(4)*PIM2(4)-VEC1(1)*PIM2(1) $ -VEC1(2)*PIM2(2)-VEC1(3)*PIM2(3) P1VEC2 =PIM1(4)*VEC2(4)-PIM1(1)*VEC2(1) $ -PIM1(2)*VEC2(2)-PIM1(3)*VEC2(3) P2VEC2 =VEC2(4)*PIM2(4)-VEC2(1)*PIM2(1) $ -VEC2(2)*PIM2(2)-VEC2(3)*PIM2(3) * HADRON CURRENT FNORM=FORMOM(XMAA,XMOM) BRAK=0.0 DO 120 JJ=1,2 DO 45 I=1,4 IF (JJ.EQ.1) THEN HADCUR(I) = FNORM *( $ VEC1(I)*(AMPI**2*P1PL-P2PL*(P12-P1PL)) $ -PIM2(I)*(P2VEC1*P1PL-P1VEC1*P2PL) $ +PIPL(I)*(P2VEC1*P12 -P1VEC1*(AMPI**2+P2PL)) ) ELSE HADCUR(I) = FNORM *( $ VEC2(I)*(AMPI**2*P1PL-P2PL*(P12-P1PL)) $ -PIM2(I)*(P2VEC2*P1PL-P1VEC2*P2PL) $ +PIPL(I)*(P2VEC2*P12 -P1VEC2*(AMPI**2+P2PL)) ) ENDIF 45 CONTINUE C * CALCULATE PI-VECTORS: VECTOR AND AXIAL CALL CLVEC(HADCUR,PN,PIVEC) CALL CLAXI(HADCUR,PN,PIAKS) CALL CLNUT(HADCUR,BRAKM,HVM) * SPIN INDEPENDENT PART OF DECAY DIFF-CROSS-SECT. IN TAU REST FRAME BRAK=BRAK+(GV**2+GA**2)*PT(4)*PIVEC(4) +2.*GV*GA*PT(4)*PIAKS(4) & +2.*(GV**2-GA**2)*AMNUTA*AMTAU*BRAKM DO 90 I=1,3 HV(I)=HV(I)-(AMTAU*((GV**2+GA**2)*PIAKS(I)+2.*GV*GA*PIVEC(I))) & +(GV**2-GA**2)*AMNUTA*AMTAU*HVM(I) 90 CONTINUE C HV IS DEFINED FOR TAU- WITH GAMMA=B+HV*POL 120 CONTINUE AMPLIT=(GFERMI*CCABIB)**2*BRAK/2. C THE STATISTICAL FACTOR FOR IDENTICAL PI'S WAS CANCELLED WITH C TWO, FOR TWO MODES OF A1 DECAY NAMELLY PI+PI-PI- AND PI-PI0PI0 C POLARIMETER VECTOR IN TAU REST FRAME DO 91 I=1,3 HV(I)=-HV(I)/BRAK 91 CONTINUE END SUBROUTINE DAMPPK(MNUM,PT,PN,PIM1,PIM2,PIM3,AMPLIT,HV) IMPLICIT NONE INCLUDE 'tauola.inc' INTEGER MNUM REAL PT(4),PN(4),PIM1(4),PIM2(4),PIM3(4),AMPLIT,HV(4) C ---------------------------------------------------------------------- * CALCULATES DIFFERENTIAL CROSS SECTION AND POLARIMETER VECTOR * FOR TAU DECAY INTO K K pi, K pi pi. * ALL SPIN EFFECTS IN THE FULL DECAY CHAIN ARE TAKEN INTO ACCOUNT. * CALCULATIONS DONE IN TAU REST FRAME WITH Z-AXIS ALONG NEUTRINO MOMENT C MNUM DECAY MODE IDENTIFIER. C C called by : DPHSAA C ---------------------------------------------------------------------- INTEGER I,K REAL DWAPI0,MXAA,MXRO1,XMAA,XMRO1,XMRO2,XMRO3 REAL PROD1,PROD2,PROD3,BRAK REAL BRAKM,ZNAK,XM1,XM2,XM3 REAL PAA(4),VEC1(4),VEC2(4),VEC3(4),VEC4(4),VEC5(4) REAL PIVEC(4),PIAKS(4),HVM(4) REAL FNORM(0:7),COEF(1:5,0:7) COMPLEX HADCUR(4),FORM1,FORM2,FORM3,FORM4,FORM5,UROJ EXTERNAL FORM1,FORM2,FORM3,FORM4,FORM5 REAL PI /3.141592653589793238462643/ INTEGER ICONT /0/ C REAL FPI DATA FPI /93.3E-3/ IF (ICONT.EQ.0) THEN ICONT=1 UROJ=CMPLX(0.0,1.0) DWAPI0=SQRT(2.0) FNORM(0)=CCABIB/FPI FNORM(1)=CCABIB/FPI FNORM(2)=CCABIB/FPI FNORM(3)=CCABIB/FPI FNORM(4)=SCABIB/FPI/DWAPI0 FNORM(5)=SCABIB/FPI FNORM(6)=SCABIB/FPI FNORM(7)=CCABIB/FPI C COEF(1,0)= 2.0*SQRT(2.)/3.0 COEF(2,0)=-2.0*SQRT(2.)/3.0 COEF(3,0)= 0.0 COEF(4,0)= FPI COEF(5,0)= 0.0 C COEF(1,1)=-SQRT(2.)/3.0 COEF(2,1)= SQRT(2.)/3.0 COEF(3,1)= 0.0 COEF(4,1)= FPI COEF(5,1)= SQRT(2.) C COEF(1,2)=-SQRT(2.)/3.0 COEF(2,2)= SQRT(2.)/3.0 COEF(3,2)= 0.0 COEF(4,2)= 0.0 COEF(5,2)=-SQRT(2.) C COEF(1,3)= 0.0 COEF(2,3)=-1.0 COEF(3,3)= 0.0 COEF(4,3)= 0.0 COEF(5,3)= 0.0 C COEF(1,4)= 1.0/SQRT(2.)/3.0 COEF(2,4)=-1.0/SQRT(2.)/3.0 COEF(3,4)= 0.0 COEF(4,4)= 0.0 COEF(5,4)= 0.0 C COEF(1,5)=-SQRT(2.)/3.0 COEF(2,5)= SQRT(2.)/3.0 COEF(3,5)= 0.0 COEF(4,5)= 0.0 COEF(5,5)=-SQRT(2.) C COEF(1,6)= 0.0 COEF(2,6)=-1.0 COEF(3,6)= 0.0 COEF(4,6)= 0.0 COEF(5,6)=-2.0 C COEF(1,7)= 0.0 COEF(2,7)= 0.0 COEF(3,7)= 0.0 COEF(4,7)= 0.0 COEF(5,7)=-SQRT(2.0/3.0) C ENDIF C DO 10 I=1,4 10 PAA(I)=PIM1(I)+PIM2(I)+PIM3(I) XMAA =SQRT(ABS(PAA(4)**2-PAA(3)**2-PAA(2)**2-PAA(1)**2)) XMRO1 =SQRT(ABS((PIM3(4)+PIM2(4))**2-(PIM3(1)+PIM2(1))**2 $ -(PIM3(2)+PIM2(2))**2-(PIM3(3)+PIM2(3))**2)) XMRO2 =SQRT(ABS((PIM3(4)+PIM1(4))**2-(PIM3(1)+PIM1(1))**2 $ -(PIM3(2)+PIM1(2))**2-(PIM3(3)+PIM1(3))**2)) XMRO3 =SQRT(ABS((PIM1(4)+PIM2(4))**2-(PIM1(1)+PIM2(1))**2 $ -(PIM1(2)+PIM2(2))**2-(PIM1(3)+PIM2(3))**2)) * ELEMENTS OF HADRON CURRENT PROD1 =PAA(4)*(PIM2(4)-PIM3(4))-PAA(1)*(PIM2(1)-PIM3(1)) $ -PAA(2)*(PIM2(2)-PIM3(2))-PAA(3)*(PIM2(3)-PIM3(3)) PROD2 =PAA(4)*(PIM3(4)-PIM1(4))-PAA(1)*(PIM3(1)-PIM1(1)) $ -PAA(2)*(PIM3(2)-PIM1(2))-PAA(3)*(PIM3(3)-PIM1(3)) PROD3 =PAA(4)*(PIM1(4)-PIM2(4))-PAA(1)*(PIM1(1)-PIM2(1)) $ -PAA(2)*(PIM1(2)-PIM2(2))-PAA(3)*(PIM1(3)-PIM2(3)) DO 40 I=1,4 VEC1(I)= PIM2(I)-PIM3(I) -PAA(I)*PROD1/XMAA**2 VEC2(I)= PIM3(I)-PIM1(I) -PAA(I)*PROD2/XMAA**2 VEC3(I)= PIM1(I)-PIM2(I) -PAA(I)*PROD3/XMAA**2 40 VEC4(I)= PIM1(I)+PIM2(I)+PIM3(I) CALL PROD5(PIM1,PIM2,PIM3,VEC5) * HADRON CURRENT C be aware that sign of vec2 is opposite to sign of vec1 in a1 case DO 45 I=1,4 HADCUR(I)= CMPLX(FNORM(MNUM)) * ( $CMPLX(VEC1(I)*COEF(1,MNUM))*FORM1(MNUM,XMAA**2,XMRO1**2,XMRO2**2)+ $CMPLX(VEC2(I)*COEF(2,MNUM))*FORM2(MNUM,XMAA**2,XMRO2**2,XMRO1**2)+ $CMPLX(VEC3(I)*COEF(3,MNUM))*FORM3(MNUM,XMAA**2,XMRO3**2,XMRO1**2)+ *(-1.0*UROJ)* $CMPLX(VEC4(I)*COEF(4,MNUM))*FORM4(MNUM,XMAA**2,XMRO1**2, $ XMRO2**2,XMRO3**2) + $(-1.0)*UROJ/4.0/PI**2/FPI**2* $CMPLX(VEC5(I)*COEF(5,MNUM))*FORM5(MNUM,XMAA**2,XMRO1**2,XMRO2**2)) 45 CONTINUE C * CALCULATE PI-VECTORS: VECTOR AND AXIAL CALL CLVEC(HADCUR,PN,PIVEC) CALL CLAXI(HADCUR,PN,PIAKS) CALL CLNUT(HADCUR,BRAKM,HVM) * SPIN INDEPENDENT PART OF DECAY DIFF-CROSS-SECT. IN TAU REST FRAME BRAK= (GV**2+GA**2)*PT(4)*PIVEC(4) +2.*GV*GA*PT(4)*PIAKS(4) & +2.*(GV**2-GA**2)*AMNUTA*AMTAU*BRAKM AMPLIT=(GFERMI)**2*BRAK/2. IF (MNUM.GE.9) THEN PRINT *, 'MNUM=',MNUM ZNAK=-1.0 XM1=0.0 XM2=0.0 XM3=0.0 DO 77 K=1,4 IF (K.EQ.4) ZNAK=1.0 XM1=ZNAK*PIM1(K)**2+XM1 XM2=ZNAK*PIM2(K)**2+XM2 XM3=ZNAK*PIM3(K)**2+XM3 77 PRINT *, 'PIM1=',PIM1(K),'PIM2=',PIM2(K),'PIM3=',PIM3(K) PRINT *, 'XM1=',SQRT(XM1),'XM2=',SQRT(XM2),'XM3=',SQRT(XM3) PRINT *, '************************************************' ENDIF C POLARIMETER VECTOR IN TAU REST FRAME DO 90 I=1,3 HV(I)=-(AMTAU*((GV**2+GA**2)*PIAKS(I)+2.*GV*GA*PIVEC(I))) & +(GV**2-GA**2)*AMNUTA*AMTAU*HVM(I) C HV IS DEFINED FOR TAU- WITH GAMMA=B+HV*POL HV(I)=-HV(I)/BRAK 90 CONTINUE END SUBROUTINE DAMPRY(ITDKRC_LOC,XK0DEC_LOC,XK,XA,QP,XN,AMPLIT,HV) IMPLICIT NONE INCLUDE 'tauola.inc' INTEGER ITDKRC_LOC DOUBLE PRECISION XK0DEC_LOC,XK(4),XA(4),QP(4),XN(4),AMPLIT,HV(4) C ---------------------------------------------------------------------- C IT CALCULATES MATRIX ELEMENT FOR THE C TAU --> MU(E) NU NUBAR DECAY MODE C INCLUDING COMPLETE ORDER ALPHA QED CORRECTIONS. C ---------------------------------------------------------------------- C DOUBLE PRECISION AK0 DOUBLE PRECISION THB,SQM2 EXTERNAL THB,SQM2 HV(4)=1.D0 AK0=XK0DEC_LOC*AMTAU IF(XK(4).LT.0.1D0*AK0) THEN AMPLIT=THB(ITDKRC_LOC,QP,XN,XA,AK0,HV) ELSE AMPLIT=SQM2(ITDKRC_LOC,QP,XN,XA,XK,AK0,HV) ENDIF RETURN END SUBROUTINE DEKAY(KTO,HX) IMPLICIT NONE INCLUDE 'tauola.inc' INTEGER KTO DOUBLE PRECISION HX(4) C *********************** C THIS DEKAY IS IN SPIRIT OF THE 'DECAY' WHICH C WAS INCLUDED IN KORAL-B PROGRAM, COMP. PHYS. COMMUN. C VOL. 36 (1985) 191, SEE COMMENTS ON GENERAL PHILOSOPHY THERE. C KTO=0 INITIALISATION (OBLIGATORY) C KTO=1,11 DENOTES TAU+ AND KTO=2,12 TAU- C DEKAY(1,H) AND DEKAY(2,H) IS CALLED INTERNALLY BY MC GENERATOR. C H DENOTES THE POLARIMETRIC VECTOR, USED BY THE HOST PROGRAM FOR C CALCULATION OF THE SPIN WEIGHT. C USER MAY OPTIONALLY CALL DEKAY(11,H) DEKAY(12,H) IN ORDER C TO TRANSFORM DECAY PRODUCTS TO CMS AND WRITE LUND RECORD IN /LUJETS/. C KTO=100, PRINT FINAL REPORT (OPTIONAL). C DECAY MODES: C JAK=1 ELECTRON DECAY C JAK=2 MU DECAY C JAK=3 PI DECAY C JAK=4 RHO DECAY C JAK=5 A1 DECAY C JAK=6 K DECAY C JAK=7 K* DECAY C JAK=8 NPI DECAY C JAK=0 INCLUSIVE: JAK=1,2,3,4,5,6,7,8 INTEGER I,K,ISGN,IDUM,JDUM,NEVTOT,NEV1,NEV2 SAVE NEVTOT,NEV1,NEV2 REAL X,PDUM REAL H(4) REAL PDUM1(4),PDUM2(4),PDUM3(4),PDUM4(4),PDUM5(4),HDUM(4) REAL PDUMX(4,9) INTEGER IWARM/0/ KTOM=KTO IF(KTO.EQ.-1) THEN C ================== C INITIALISATION OR REINITIALISATION KTOM=1 IF (IWARM.EQ.1) X=5/(IWARM-1) IWARM=1 WRITE(IOUT,7001) JAK1,JAK2 NEVTOT=0 NEV1=0 NEV2=0 IF(JAK1.NE.-1.OR.JAK2.NE.-1) THEN CALL DADMEL(-1,IDUM,HDUM,PDUM1,PDUM2,PDUM3,PDUM4,PDUM5) CALL DADMMU(-1,IDUM,HDUM,PDUM1,PDUM2,PDUM3,PDUM4,PDUM5) CALL DADMPI(-1,IDUM,PDUM,PDUM1,PDUM2) CALL DADMRO(-1,IDUM,HDUM,PDUM1,PDUM2,PDUM3,PDUM4) CALL DADMAA(-1,IDUM,HDUM,PDUM1,PDUM2,PDUM3,PDUM4,PDUM5,JDUM) CALL DADMKK(-1,IDUM,PDUM,PDUM1,PDUM2) CALL DADMKS(-1,IDUM,HDUM,PDUM1,PDUM2,PDUM3,PDUM4,JDUM) CALL DADNEW(-1,IDUM,HDUM,PDUM1,PDUM2,PDUMX,JDUM) ENDIF DO 21 I=1,30 NEVDEC(I)=0 GAMPMC(I)=0 21 GAMPER(I)=0 ELSEIF(KTO.EQ.1) THEN C ===================== C DECAY OF TAU+ IN THE TAU REST FRAME NEVTOT=NEVTOT+1 IF(IWARM.EQ.0) GOTO 902 ISGN= IDFF/IABS(IDFF) CALL DEKAY1(0,H,ISGN) ELSEIF(KTO.EQ.2) THEN C ================================= C DECAY OF TAU- IN THE TAU REST FRAME NEVTOT=NEVTOT+1 IF(IWARM.EQ.0) GOTO 902 ISGN=-IDFF/IABS(IDFF) CALL DEKAY2(0,H,ISGN) ELSEIF(KTO.EQ.11) THEN C ====================== C REST OF DECAY PROCEDURE FOR ACCEPTED TAU+ DECAY NEV1=NEV1+1 ISGN= IDFF/IABS(IDFF) CALL DEKAY1(1,H,ISGN) ELSEIF(KTO.EQ.12) THEN C ====================== C REST OF DECAY PROCEDURE FOR ACCEPTED TAU- DECAY NEV2=NEV2+1 ISGN=-IDFF/IABS(IDFF) CALL DEKAY2(1,H,ISGN) ELSEIF(KTO.EQ.100) THEN C ======================= IF(JAK1.NE.-1.OR.JAK2.NE.-1) THEN CALL DADMEL( 1,IDUM,HDUM,PDUM1,PDUM2,PDUM3,PDUM4,PDUM5) CALL DADMMU( 1,IDUM,HDUM,PDUM1,PDUM2,PDUM3,PDUM4,PDUM5) CALL DADMPI( 1,IDUM,PDUM,PDUM1,PDUM2) CALL DADMRO( 1,IDUM,HDUM,PDUM1,PDUM2,PDUM3,PDUM4) CALL DADMAA( 1,IDUM,HDUM,PDUM1,PDUM2,PDUM3,PDUM4,PDUM5,JDUM) CALL DADMKK( 1,IDUM,PDUM,PDUM1,PDUM2) CALL DADMKS( 1,IDUM,HDUM,PDUM1,PDUM2,PDUM3,PDUM4,JDUM) CALL DADNEW( 1,IDUM,HDUM,PDUM1,PDUM2,PDUMX,JDUM) WRITE(IOUT,7010) NEV1,NEV2,NEVTOT WRITE(IOUT,7011) (NEVDEC(I),GAMPMC(I),GAMPER(I),I= 1,7) WRITE(IOUT,7012) $ (NEVDEC(I),GAMPMC(I),GAMPER(I),NAMES(I-7),I=8,7+NMODE) WRITE(IOUT,7013) ENDIF ELSE C ==== GOTO 910 ENDIF C ===== DO 78 K=1,4 78 HX(K)=H(K) RETURN 7001 FORMAT(///1X,15(5H*****) $ /,' *', 25X,'*****TAUOLA LIBRARY: VERSION 2.5 ******',9X,1H*, $ /,' *', 25X,'***********JUNE 1994***************',9X,1H*, $ /,' *', 25X,'**AUTHORS: S.JADACH, Z.WAS*************',9X,1H*, $ /,' *', 25X,'**R. DECKER, M. JEZABEK, J.H.KUEHN*****',9X,1H*, $ /,' *', 25X,'**AVAILABLE FROM: WASM AT CERNVM ******',9X,1H*, $ /,' *', 25X,'***** PUBLISHED IN COMP. PHYS. COMM.***',9X,1H*, $ /,' *', 25X,'*******CERN-TH-5856 SEPTEMBER 1990*****',9X,1H*, $ /,' *', 25X,'*******CERN-TH-6195 SEPTEMBER 1991*****',9X,1H*, $ /,' *', 25X,'*******CERN TH-6793 NOVEMBER 1992*****',9X,1H*, $ /,' *', 25X,'**5 or more pi dec.: precision limited ',9X,1H*, $ /,' *', 25X,'****DEKAY ROUTINE: INITIALIZATION******',9X,1H*, $ /,' *',I20 ,5X,'JAK1 = DECAY MODE TAU+ ',9X,1H*, $ /,' *',I20 ,5X,'JAK2 = DECAY MODE TAU- ',9X,1H*, $ /,1X,15(5H*****)/) 7010 FORMAT(///1X,15(5H*****) $ /,' *', 25X,'*****TAUOLA LIBRARY: VERSION 2.5 ******',9X,1H*, $ /,' *', 25X,'***********JUNE 1994***************',9X,1H*, $ /,' *', 25X,'**AUTHORS: S.JADACH, Z.WAS*************',9X,1H*, $ /,' *', 25X,'**R. DECKER, M. JEZABEK, J.H.KUEHN*****',9X,1H*, $ /,' *', 25X,'**AVAILABLE FROM: WASM AT CERNVM ******',9X,1H*, $ /,' *', 25X,'***** PUBLISHED IN COMP. PHYS. COMM.***',9X,1H*, $ /,' *', 25X,'*******CERN-TH-5856 SEPTEMBER 1990*****',9X,1H*, $ /,' *', 25X,'*******CERN-TH-6195 SEPTEMBER 1991*****',9X,1H*, $ /,' *', 25X,'*******CERN TH-6793 NOVEMBER 1992*****',9X,1H*, $ /,' *', 25X,'*****DEKAY ROUTINE: FINAL REPORT*******',9X,1H*, $ /,' *',I20 ,5X,'NEV1 = NO. OF TAU+ DECS. ACCEPTED ',9X,1H*, $ /,' *',I20 ,5X,'NEV2 = NO. OF TAU- DECS. ACCEPTED ',9X,1H*, $ /,' *',I20 ,5X,'NEVTOT = SUM ',9X,1H*, $ /,' *',' NOEVTS ', $ ' PART.WIDTH ERROR ROUTINE DECAY MODE ',9X,1H*) 7011 FORMAT(1X,'*' $ ,I10,2F12.7 ,' DADMEL ELECTRON ',9X,1H* $ /,' *',I10,2F12.7 ,' DADMMU MUON ',9X,1H* $ /,' *',I10,2F12.7 ,' DADMPI PION ',9X,1H* $ /,' *',I10,2F12.7, ' DADMRO RHO (->2PI) ',9X,1H* $ /,' *',I10,2F12.7, ' DADMAA A1 (->3PI) ',9X,1H* $ /,' *',I10,2F12.7, ' DADMKK KAON ',9X,1H* $ /,' *',I10,2F12.7, ' DADMKS K* ',9X,1H*) 7012 FORMAT(1X,'*' $ ,I10,2F12.7,A31 ,8X,1H*) 7013 FORMAT(1X,'*' $ ,20X,'THE ERROR IS RELATIVE AND PART.WIDTH ',10X,1H* $ /,' *',20X,'IN UNITS GFERMI**2*MASS**5/192/PI**3 ',10X,1H* $ /,1X,15(5H*****)/) 902 PRINT 9020 9020 FORMAT(' ----- DEKAY: LACK OF INITIALISATION') STOP 910 PRINT 9100 9100 FORMAT(' ----- DEKAY: WRONG VALUE OF KTO ') STOP END SUBROUTINE DEKAY1(IMOD,HH,ISGN) IMPLICIT NONE INCLUDE 'tauola.inc' INTEGER IMOD REAL HH(4) INTEGER ISGN C ******************************* C THIS ROUTINE SIMULATES TAU+ DECAY INTEGER I,KTO,IMD,JAK,JAA,JKST SAVE JAK REAL HV(4),PNU(4),PPI(4) REAL PWB(4),PMU(4),PNM(4) REAL PRHO(4),PIC(4),PIZ(4) REAL PAA(4),PIM1(4),PIM2(4),PIPL(4) REAL PKK(4),PKS(4) REAL PNPI(4,9) REAL PHOT(4) REAL PDUM(4) INTEGER NEV/0/,NPRIN/10/ KTO=1 IF(JAK1.EQ.-1) RETURN IMD=IMOD IF(IMD.EQ.0) THEN C ================= JAK=JAK1 IF(JAK1.EQ.0) CALL JAKER(JAK) IF(JAK.EQ.1) THEN CALL DADMEL(0, ISGN,HV,PNU,PWB,PMU,PNM,PHOT) ELSEIF(JAK.EQ.2) THEN CALL DADMMU(0, ISGN,HV,PNU,PWB,PMU,PNM,PHOT) ELSEIF(JAK.EQ.3) THEN CALL DADMPI(0, ISGN,HV,PPI,PNU) ELSEIF(JAK.EQ.4) THEN CALL DADMRO(0, ISGN,HV,PNU,PRHO,PIC,PIZ) ELSEIF(JAK.EQ.5) THEN CALL DADMAA(0, ISGN,HV,PNU,PAA,PIM1,PIM2,PIPL,JAA) ELSEIF(JAK.EQ.6) THEN CALL DADMKK(0, ISGN,HV,PKK,PNU) ELSEIF(JAK.EQ.7) THEN CALL DADMKS(0, ISGN,HV,PNU,PKS ,PKK,PPI,JKST) ELSE CALL DADNEW(0, ISGN,HV,PNU,PWB,PNPI,JAK-7) ENDIF DO 33 I=1,3 33 HH(I)=HV(I) HH(4)=1.0 ELSEIF(IMD.EQ.1) THEN C ===================== NEV=NEV+1 IF (JAK.LT.31) THEN NEVDEC(JAK)=NEVDEC(JAK)+1 ENDIF DO 34 I=1,4 34 PDUM(I)=.0 IF(JAK.EQ.1) THEN CALL DWLUEL(1,ISGN,PNU,PWB,PMU,PNM) CALL DWRPH(KTOM,PHOT) DO 10 I=1,4 10 PP1(I)=PMU(I) ELSEIF(JAK.EQ.2) THEN CALL DWLUMU(1,ISGN,PNU,PWB,PMU,PNM) CALL DWRPH(KTOM,PHOT) DO 20 I=1,4 20 PP1(I)=PMU(I) ELSEIF(JAK.EQ.3) THEN CALL DWLUPI(1,ISGN,PPI,PNU) DO 30 I=1,4 30 PP1(I)=PPI(I) ELSEIF(JAK.EQ.4) THEN CALL DWLURO(1,ISGN,PNU,PRHO,PIC,PIZ) DO 40 I=1,4 40 PP1(I)=PRHO(I) ELSEIF(JAK.EQ.5) THEN CALL DWLUAA(1,ISGN,PNU,PAA,PIM1,PIM2,PIPL,JAA) DO 50 I=1,4 50 PP1(I)=PAA(I) ELSEIF(JAK.EQ.6) THEN CALL DWLUKK(1,ISGN,PKK,PNU) DO 60 I=1,4 60 PP1(I)=PKK(I) ELSEIF(JAK.EQ.7) THEN CALL DWLUKS(1,ISGN,PNU,PKS,PKK,PPI,JKST) DO 70 I=1,4 70 PP1(I)=PKS(I) ELSE CAM MULTIPION DECAY CALL DWLNEW(1,ISGN,PNU,PWB,PNPI,JAK) DO 80 I=1,4 80 PP1(I)=PWB(I) ENDIF ENDIF C ===== END SUBROUTINE DEKAY2(IMOD,HH,ISGN) IMPLICIT NONE INCLUDE 'tauola.inc' INTEGER IMOD REAL HH(4) INTEGER ISGN C ******************************* C THIS ROUTINE SIMULATES TAU- DECAY INTEGER I,KTO,IMD,JAK,JAA,JKST SAVE JAK REAL HV(4),PNU(4),PPI(4) REAL PWB(4),PMU(4),PNM(4) REAL PRHO(4),PIC(4),PIZ(4) REAL PAA(4),PIM1(4),PIM2(4),PIPL(4) REAL PKK(4),PKS(4) REAL PNPI(4,9) REAL PHOT(4) REAL PDUM(4) INTEGER NEV/0/,NPRIN/10/ KTO=2 IF(JAK2.EQ.-1) RETURN IMD=IMOD IF(IMD.EQ.0) THEN C ================= JAK=JAK2 IF(JAK2.EQ.0) CALL JAKER(JAK) IF(JAK.EQ.1) THEN CALL DADMEL(0, ISGN,HV,PNU,PWB,PMU,PNM,PHOT) ELSEIF(JAK.EQ.2) THEN CALL DADMMU(0, ISGN,HV,PNU,PWB,PMU,PNM,PHOT) ELSEIF(JAK.EQ.3) THEN CALL DADMPI(0, ISGN,HV,PPI,PNU) ELSEIF(JAK.EQ.4) THEN CALL DADMRO(0, ISGN,HV,PNU,PRHO,PIC,PIZ) ELSEIF(JAK.EQ.5) THEN CALL DADMAA(0, ISGN,HV,PNU,PAA,PIM1,PIM2,PIPL,JAA) ELSEIF(JAK.EQ.6) THEN CALL DADMKK(0, ISGN,HV,PKK,PNU) ELSEIF(JAK.EQ.7) THEN CALL DADMKS(0, ISGN,HV,PNU,PKS ,PKK,PPI,JKST) ELSE CALL DADNEW(0, ISGN,HV,PNU,PWB,PNPI,JAK-7) ENDIF DO 33 I=1,3 33 HH(I)=HV(I) HH(4)=1.0 ELSEIF(IMD.EQ.1) THEN C ===================== NEV=NEV+1 IF (JAK.LT.31) THEN NEVDEC(JAK)=NEVDEC(JAK)+1 ENDIF DO 34 I=1,4 34 PDUM(I)=.0 IF(JAK.EQ.1) THEN CALL DWLUEL(2,ISGN,PNU,PWB,PMU,PNM) CALL DWRPH(KTOM,PHOT) DO 10 I=1,4 10 PP2(I)=PMU(I) ELSEIF(JAK.EQ.2) THEN CALL DWLUMU(2,ISGN,PNU,PWB,PMU,PNM) CALL DWRPH(KTOM,PHOT) DO 20 I=1,4 20 PP2(I)=PMU(I) ELSEIF(JAK.EQ.3) THEN CALL DWLUPI(2,ISGN,PPI,PNU) DO 30 I=1,4 30 PP2(I)=PPI(I) ELSEIF(JAK.EQ.4) THEN CALL DWLURO(2,ISGN,PNU,PRHO,PIC,PIZ) DO 40 I=1,4 40 PP2(I)=PRHO(I) ELSEIF(JAK.EQ.5) THEN CALL DWLUAA(2,ISGN,PNU,PAA,PIM1,PIM2,PIPL,JAA) DO 50 I=1,4 50 PP2(I)=PAA(I) ELSEIF(JAK.EQ.6) THEN CALL DWLUKK(2,ISGN,PKK,PNU) DO 60 I=1,4 60 PP1(I)=PKK(I) ELSEIF(JAK.EQ.7) THEN CALL DWLUKS(2,ISGN,PNU,PKS,PKK,PPI,JKST) DO 70 I=1,4 70 PP1(I)=PKS(I) ELSE CAM MULTIPION DECAY CALL DWLNEW(2,ISGN,PNU,PWB,PNPI,JAK) DO 80 I=1,4 80 PP1(I)=PWB(I) ENDIF C ENDIF C ===== END SUBROUTINE DEXAA(MODE,ISGN,POL,PNU,PAA,PIM1,PIM2,PIPL,JAA) IMPLICIT NONE INCLUDE 'tauola.inc' INTEGER MODE,ISGN REAL POL(4),PNU(4),PAA(4),PIM1(4),PIM2(4),PIPL(4) INTEGER JAA C ---------------------------------------------------------------------- * THIS SIMULATES TAU DECAY IN TAU REST FRAME * INTO NU A1, NEXT A1 DECAYS INTO RHO PI AND FINALLY RHO INTO PI PI. * OUTPUT FOUR MOMENTA: PNU TAUNEUTRINO, * PAA A1 * PIM1 PION MINUS (OR PI0) 1 (FOR TAU MINUS) * PIM2 PION MINUS (OR PI0) 2 * PIPL PION PLUS (OR PI-) * (PIPL,PIM1) FORM A RHO C ---------------------------------------------------------------------- REAL HV(4) REAL WT,RN INTEGER IWARM/0/ C IF(MODE.EQ.-1) THEN C =================== IWARM=1 CALL DADMAA( -1,ISGN,HV,PNU,PAA,PIM1,PIM2,PIPL,JAA) CC CALL HBOOK1(816,'WEIGHT DISTRIBUTION DEXAA $',100,-2.,2.) C ELSEIF(MODE.EQ. 0) THEN * ======================= 300 CONTINUE IF(IWARM.EQ.0) GOTO 902 CALL DADMAA( 0,ISGN,HV,PNU,PAA,PIM1,PIM2,PIPL,JAA) WT=(1+POL(1)*HV(1)+POL(2)*HV(2)+POL(3)*HV(3))/2. CC CALL HFILL(816,WT) CALL TDRAND(RN,1) IF(RN.GT.WT) GOTO 300 C ELSEIF(MODE.EQ. 1) THEN * ======================= CALL DADMAA( 1,ISGN,HV,PNU,PAA,PIM1,PIM2,PIPL,JAA) CC CALL HPRINT(816) ENDIF C ===== RETURN 902 WRITE(IOUT, 9020) 9020 FORMAT(' ----- DEXAA: LACK OF INITIALISATION') STOP END SUBROUTINE DEXAY(KTO,POL) IMPLICIT NONE INCLUDE 'tauola.inc' INTEGER KTO REAL POL(4) C ---------------------------------------------------------------------- C THIS 'DEXAY' IS A ROUTINE WHICH GENERATES DECAY OF THE SINGLE C POLARIZED TAU, POL IS A POLARIZATION VECTOR (NOT A POLARIMETER C VECTOR AS IN DEKAY) OF THE TAU AND IT IS AN INPUT PARAMETER. C KTO=0 INITIALISATION (OBLIGATORY) C KTO=1 DENOTES TAU+ AND KTO=2 TAU- C DEXAY(1,POL) AND DEXAY(2,POL) ARE CALLED INTERNALLY BY MC GENERATOR. C DECAY PRODUCTS ARE TRANSFORMED READILY C TO CMS AND WRITEN IN THE LUND RECORD IN /LUJETS/ C KTO=100, PRINT FINAL REPORT (OPTIONAL). C C called by : KORALZ C ---------------------------------------------------------------------- INTEGER I,IDUM,ISGN,NEV1,NEV2,NEVTOT SAVE NEV1,NEV2,NEVTOT REAL PDUM1(4),PDUM2(4),PDUM3(4),PDUM4(4),PDUM5(4) REAL PDUM(4) REAL PDUMI(4,9) INTEGER IWARM/0/ KTOM=KTO C IF(KTO.EQ.-1) THEN C ================== C INITIALISATION OR REINITIALISATION IWARM=1 WRITE(IOUT, 7001) JAK1,JAK2 NEVTOT=0 NEV1=0 NEV2=0 IF(JAK1.NE.-1.OR.JAK2.NE.-1) THEN CALL DEXEL(-1,IDUM,PDUM,PDUM1,PDUM2,PDUM3,PDUM4,PDUM5) CALL DEXMU(-1,IDUM,PDUM,PDUM1,PDUM2,PDUM3,PDUM4,PDUM5) CALL DEXPI(-1,IDUM,PDUM,PDUM1,PDUM2) CALL DEXRO(-1,IDUM,PDUM,PDUM1,PDUM2,PDUM3,PDUM4) CALL DEXAA(-1,IDUM,PDUM,PDUM1,PDUM2,PDUM3,PDUM4,PDUM5,IDUM) CALL DEXKK(-1,IDUM,PDUM,PDUM1,PDUM2) CALL DEXKS(-1,IDUM,PDUM,PDUM1,PDUM2,PDUM3,PDUM4,IDUM) CALL DEXNEW(-1,IDUM,PDUM,PDUM1,PDUM2,PDUMI,IDUM) ENDIF DO 21 I=1,30 NEVDEC(I)=0 GAMPMC(I)=0 21 GAMPER(I)=0 ELSEIF(KTO.EQ.1) THEN C ===================== C DECAY OF TAU+ IN THE TAU REST FRAME NEVTOT=NEVTOT+1 NEV1=NEV1+1 IF(IWARM.EQ.0) GOTO 902 ISGN=IDFF/IABS(IDFF) CAM CALL DEXAY1(POL,ISGN) CALL DEXAY1(KTO,JAK1,JAKP,POL,ISGN) ELSEIF(KTO.EQ.2) THEN C ================================= C DECAY OF TAU- IN THE TAU REST FRAME NEVTOT=NEVTOT+1 NEV2=NEV2+1 IF(IWARM.EQ.0) GOTO 902 ISGN=-IDFF/IABS(IDFF) CAM CALL DEXAY2(POL,ISGN) CALL DEXAY1(KTO,JAK2,JAKM,POL,ISGN) ELSEIF(KTO.EQ.100) THEN C ======================= IF(JAK1.NE.-1.OR.JAK2.NE.-1) THEN CALL DEXEL( 1,IDUM,PDUM,PDUM1,PDUM2,PDUM3,PDUM4,PDUM5) CALL DEXMU( 1,IDUM,PDUM,PDUM1,PDUM2,PDUM3,PDUM4,PDUM5) CALL DEXPI( 1,IDUM,PDUM,PDUM1,PDUM2) CALL DEXRO( 1,IDUM,PDUM,PDUM1,PDUM2,PDUM3,PDUM4) CALL DEXAA( 1,IDUM,PDUM,PDUM1,PDUM2,PDUM3,PDUM4,PDUM5,IDUM) CALL DEXKK( 1,IDUM,PDUM,PDUM1,PDUM2) CALL DEXKS( 1,IDUM,PDUM,PDUM1,PDUM2,PDUM3,PDUM4,IDUM) CALL DEXNEW( 1,IDUM,PDUM,PDUM1,PDUM2,PDUMI,IDUM) WRITE(IOUT,7010) NEV1,NEV2,NEVTOT WRITE(IOUT,7011) (NEVDEC(I),GAMPMC(I),GAMPER(I),I= 1,7) WRITE(IOUT,7012) $ (NEVDEC(I),GAMPMC(I),GAMPER(I),NAMES(I-7),I=8,7+NMODE) WRITE(IOUT,7013) ENDIF ELSE GOTO 910 ENDIF RETURN 7001 FORMAT(///1X,15(5H*****) $ /,' *', 25X,'*****TAUOLA LIBRARY: VERSION 2.5 ******',9X,1H*, $ /,' *', 25X,'***********JUNE 1994***************',9X,1H*, $ /,' *', 25X,'**AUTHORS: S.JADACH, Z.WAS*************',9X,1H*, $ /,' *', 25X,'**R. DECKER, M. JEZABEK, J.H.KUEHN*****',9X,1H*, $ /,' *', 25X,'**AVAILABLE FROM: WASM AT CERNVM ******',9X,1H*, $ /,' *', 25X,'***** PUBLISHED IN COMP. PHYS. COMM.***',9X,1H*, $ /,' *', 25X,'*******CERN-TH-5856 SEPTEMBER 1990*****',9X,1H*, $ /,' *', 25X,'*******CERN-TH-6195 SEPTEMBER 1991*****',9X,1H*, $ /,' *', 25X,'*******CERN-TH-6793 NOVEMBER 1992*****',9X,1H*, $ /,' *', 25X,'**5 or more pi dec.: precision limited ',9X,1H*, $ /,' *', 25X,'******DEXAY ROUTINE: INITIALIZATION****',9X,1H* $ /,' *',I20 ,5X,'JAK1 = DECAY MODE FERMION1 (TAU+) ',9X,1H* $ /,' *',I20 ,5X,'JAK2 = DECAY MODE FERMION2 (TAU-) ',9X,1H* $ /,1X,15(5H*****)/) CHBU format 7010 had more than 19 continuation lines CHBU split into two 7010 FORMAT(///1X,15(5H*****) $ /,' *', 25X,'*****TAUOLA LIBRARY: VERSION 2.5 ******',9X,1H*, $ /,' *', 25X,'***********JUNE 1994***************',9X,1H*, $ /,' *', 25X,'**AUTHORS: S.JADACH, Z.WAS*************',9X,1H*, $ /,' *', 25X,'**R. DECKER, M. JEZABEK, J.H.KUEHN*****',9X,1H*, $ /,' *', 25X,'**AVAILABLE FROM: WASM AT CERNVM ******',9X,1H*, $ /,' *', 25X,'***** PUBLISHED IN COMP. PHYS. COMM.***',9X,1H*, $ /,' *', 25X,'*******CERN-TH-5856 SEPTEMBER 1990*****',9X,1H*, $ /,' *', 25X,'*******CERN-TH-6195 SEPTEMBER 1991*****',9X,1H*, $ /,' *', 25X,'*******CERN-TH-6793 NOVEMBER 1992*****',9X,1H*, $ /,' *', 25X,'******DEXAY ROUTINE: FINAL REPORT******',9X,1H* $ /,' *',I20 ,5X,'NEV1 = NO. OF TAU+ DECS. ACCEPTED ',9X,1H* $ /,' *',I20 ,5X,'NEV2 = NO. OF TAU- DECS. ACCEPTED ',9X,1H* $ /,' *',I20 ,5X,'NEVTOT = SUM ',9X,1H* $ /,' *',' NOEVTS ', $ ' PART.WIDTH ERROR ROUTINE DECAY MODE ',9X,1H*) 7011 FORMAT(1X,'*' $ ,I10,2F12.7 ,' DADMEL ELECTRON ',9X,1H* $ /,' *',I10,2F12.7 ,' DADMMU MUON ',9X,1H* $ /,' *',I10,2F12.7 ,' DADMPI PION ',9X,1H* $ /,' *',I10,2F12.7, ' DADMRO RHO (->2PI) ',9X,1H* $ /,' *',I10,2F12.7, ' DADMAA A1 (->3PI) ',9X,1H* $ /,' *',I10,2F12.7, ' DADMKK KAON ',9X,1H* $ /,' *',I10,2F12.7, ' DADMKS K* ',9X,1H*) 7012 FORMAT(1X,'*' $ ,I10,2F12.7,A31 ,8X,1H*) 7013 FORMAT(1X,'*' $ ,20X,'THE ERROR IS RELATIVE AND PART.WIDTH ',10X,1H* $ /,' *',20X,'IN UNITS GFERMI**2*MASS**5/192/PI**3 ',10X,1H* $ /,1X,15(5H*****)/) 902 WRITE(IOUT, 9020) 9020 FORMAT(' ----- DEXAY: LACK OF INITIALISATION') STOP 910 WRITE(IOUT, 9100) 9100 FORMAT(' ----- DEXAY: WRONG VALUE OF KTO ') STOP END SUBROUTINE DEXAY1(KTO,JAKIN,JAK,POL,ISGN) IMPLICIT NONE INCLUDE 'tauola.inc' INTEGER KTO,JAKIN,JAK REAL POL(4) INTEGER ISGN C --------------------------------------------------------------------- C THIS ROUTINE SIMULATES TAU+- DECAY C C called by : DEXAY C --------------------------------------------------------------------- INTEGER I,JAA,JKST,JNPI REAL POLAR(4) REAL PNU(4),PPI(4) REAL PRHO(4),PIC(4),PIZ(4) REAL PWB(4),PMU(4),PNM(4) REAL PAA(4),PIM1(4),PIM2(4),PIPL(4) REAL PKK(4),PKS(4) REAL PNPI(4,9) REAL PHOT(4) REAL PDUM(4) C IF(JAKIN.EQ.-1) RETURN DO 33 I=1,3 33 POLAR(I)=POL(I) POLAR(4)=0. DO 34 I=1,4 34 PDUM(I)=.0 JAK=JAKIN IF(JAK.EQ.0) CALL JAKER(JAK) CAM IF(JAK.EQ.1) THEN CALL DEXEL(0, ISGN,POLAR,PNU,PWB,PMU,PNM,PHOT) CALL DWLUEL(KTO,ISGN,PNU,PWB,PMU,PNM) CALL DWRPH(KTO,PHOT ) ELSEIF(JAK.EQ.2) THEN CALL DEXMU(0, ISGN,POLAR,PNU,PWB,PMU,PNM,PHOT) CALL DWLUMU(KTO,ISGN,PNU,PWB,PMU,PNM) CALL DWRPH(KTO,PHOT ) ELSEIF(JAK.EQ.3) THEN CALL DEXPI(0, ISGN,POLAR,PPI,PNU) CALL DWLUPI(KTO,ISGN,PPI,PNU) ELSEIF(JAK.EQ.4) THEN CALL DEXRO(0, ISGN,POLAR,PNU,PRHO,PIC,PIZ) CALL DWLURO(KTO,ISGN,PNU,PRHO,PIC,PIZ) ELSEIF(JAK.EQ.5) THEN CALL DEXAA(0, ISGN,POLAR,PNU,PAA,PIM1,PIM2,PIPL,JAA) CALL DWLUAA(KTO,ISGN,PNU,PAA,PIM1,PIM2,PIPL,JAA) ELSEIF(JAK.EQ.6) THEN CALL DEXKK(0, ISGN,POLAR,PKK,PNU) CALL DWLUKK(KTO,ISGN,PKK,PNU) ELSEIF(JAK.EQ.7) THEN CALL DEXKS(0, ISGN,POLAR,PNU,PKS,PKK,PPI,JKST) CALL DWLUKS(KTO,ISGN,PNU,PKS,PKK,PPI,JKST) ELSE JNPI=JAK-7 CALL DEXNEW(0, ISGN,POLAR,PNU,PWB,PNPI,JNPI) CALL DWLNEW(KTO,ISGN,PNU,PWB,PNPI,JAK) ENDIF NEVDEC(JAK)=NEVDEC(JAK)+1 END SUBROUTINE DEXEL(MODE,ISGN,POL,PNU,PWB,Q1,Q2,PH) IMPLICIT NONE INTEGER MODE,ISGN REAL POL(4),PNU(4),PWB(4),Q1(4),Q2(4),PH(4) C ---------------------------------------------------------------------- C THIS SIMULATES TAU DECAY IN TAU REST FRAME C INTO ELECTRON AND TWO NEUTRINOS C C called by : DEXAY,DEXAY1 C ---------------------------------------------------------------------- REAL WT,RN REAL HV(4) INTEGER IWARM/0/ C IF(MODE.EQ.-1) THEN C =================== IWARM=1 CALL DADMEL( -1,ISGN,HV,PNU,PWB,Q1,Q2,PH) CC CALL HBOOK1(813,'WEIGHT DISTRIBUTION DEXEL $',100,0,2) C ELSEIF(MODE.EQ. 0) THEN C ======================= 300 CONTINUE IF(IWARM.EQ.0) GOTO 902 CALL DADMEL( 0,ISGN,HV,PNU,PWB,Q1,Q2,PH) WT=(1+POL(1)*HV(1)+POL(2)*HV(2)+POL(3)*HV(3))/2. CC CALL HFILL(813,WT) CALL TDRAND(RN,1) IF(RN.GT.WT) GOTO 300 C ELSEIF(MODE.EQ. 1) THEN C ======================= CALL DADMEL( 1,ISGN,HV,PNU,PWB,Q1,Q2,PH) CC CALL HPRINT(813) ENDIF C ===== RETURN 902 PRINT 9020 9020 FORMAT(' ----- DEXEL: LACK OF INITIALISATION') STOP END SUBROUTINE DEXKK(MODE,ISGN,POL,PKK,PNU) IMPLICIT NONE INTEGER MODE,ISGN REAL POL(4),PKK(4),PNU(4) C ---------------------------------------------------------------------- C TAU DECAY INTO KAON AND TAU-NEUTRINO C IN TAU REST FRAME C OUTPUT FOUR MOMENTA: PNU TAUNEUTRINO, C PKK KAON CHARGED C ---------------------------------------------------------------------- REAL WT,RN REAL HV(4) C IF(MODE.EQ.-1) THEN C =================== CALL DADMKK(-1,ISGN,HV,PKK,PNU) CC CALL HBOOK1(815,'WEIGHT DISTRIBUTION DEXPI $',100,0,2) C ELSEIF(MODE.EQ. 0) THEN C ======================= 300 CONTINUE CALL DADMKK( 0,ISGN,HV,PKK,PNU) WT=(1+POL(1)*HV(1)+POL(2)*HV(2)+POL(3)*HV(3))/2. CC CALL HFILL(815,WT) CALL TDRAND(RN,1) IF(RN.GT.WT) GOTO 300 C ELSEIF(MODE.EQ. 1) THEN C ======================= CALL DADMKK( 1,ISGN,HV,PKK,PNU) CC CALL HPRINT(815) ENDIF C ===== RETURN END SUBROUTINE DEXKS(MODE,ISGN,POL,PNU,PKS,PKK,PPI,JKST) IMPLICIT NONE INCLUDE 'tauola.inc' INTEGER MODE,ISGN REAL POL(4),PNU(4),PKS(4),PKK(4),PPI(4) INTEGER JKST C ---------------------------------------------------------------------- C THIS SIMULATES TAU DECAY IN TAU REST FRAME C INTO NU K*, THEN K* DECAYS INTO PI0,K+-(JKST=20) C OR PI+-,K0(JKST=10). C OUTPUT FOUR MOMENTA: PNU TAUNEUTRINO, C PKS K* CHARGED C PK0 K ZERO C PKC K CHARGED C PIC PION CHARGED C PIZ PION ZERO C ---------------------------------------------------------------------- REAL WT,RN REAL HV(4) INTEGER IWARM/0/ C IF(MODE.EQ.-1) THEN C =================== IWARM=1 CFZ INITIALISATION DONE WITH THE GHARGED PION NEUTRAL KAON MODE(JKST=10 CALL DADMKS( -1,ISGN,HV,PNU,PKS,PKK,PPI,JKST) CC CALL HBOOK1(816,'WEIGHT DISTRIBUTION DEXKS $',100,0,2) CC CALL HBOOK1(916,'ABS2 OF HV IN ROUTINE DEXKS $',100,0,2) C ELSEIF(MODE.EQ. 0) THEN C ======================= 300 CONTINUE IF(IWARM.EQ.0) GOTO 902 CALL DADMKS( 0,ISGN,HV,PNU,PKS,PKK,PPI,JKST) WT=(1+POL(1)*HV(1)+POL(2)*HV(2)+POL(3)*HV(3))/2. CC CALL HFILL(816,WT) CC XHELP=HV(1)**2+HV(2)**2+HV(3)**2 CC CALL HFILL(916,XHELP) CALL TDRAND(RN,1) IF(RN.GT.WT) GOTO 300 C ELSEIF(MODE.EQ. 1) THEN C ====================================== CALL DADMKS( 1,ISGN,HV,PNU,PKS,PKK,PPI,JKST) CC CALL HPRINT(816) CC CALL HPRINT(916) ENDIF C ===== RETURN 902 WRITE(IOUT, 9020) 9020 FORMAT(' ----- DEXKS: LACK OF INITIALISATION') STOP END SUBROUTINE DEXMU(MODE,ISGN,POL,PNU,PWB,Q1,Q2,PH) IMPLICIT NONE INCLUDE 'tauola.inc' INTEGER MODE,ISGN REAL POL(4),PNU(4),PWB(4),Q1(4),Q2(4),PH(4) C ---------------------------------------------------------------------- C THIS SIMULATES TAU DECAY IN ITS REST FRAME C INTO MUON AND TWO NEUTRINOS C OUTPUT FOUR MOMENTA: PNU TAUNEUTRINO, C PWB W-BOSON C Q1 MUON C Q2 MUON-NEUTRINO C ---------------------------------------------------------------------- REAL WT,RN REAL HV(4) INTEGER IWARM/0/ C IF(MODE.EQ.-1) THEN C =================== IWARM=1 CALL DADMMU( -1,ISGN,HV,PNU,PWB,Q1,Q2,PH) CC CALL HBOOK1(814,'WEIGHT DISTRIBUTION DEXMU $',100,0,2) C ELSEIF(MODE.EQ. 0) THEN C ======================= 300 CONTINUE IF(IWARM.EQ.0) GOTO 902 CALL DADMMU( 0,ISGN,HV,PNU,PWB,Q1,Q2,PH) WT=(1+POL(1)*HV(1)+POL(2)*HV(2)+POL(3)*HV(3))/2. CC CALL HFILL(814,WT) CALL TDRAND(RN,1) IF(RN.GT.WT) GOTO 300 C ELSEIF(MODE.EQ. 1) THEN C ======================= CALL DADMMU( 1,ISGN,HV,PNU,PWB,Q1,Q2,PH) CC CALL HPRINT(814) ENDIF C ===== RETURN 902 WRITE(IOUT, 9020) 9020 FORMAT(' ----- DEXMU: LACK OF INITIALISATION') STOP END SUBROUTINE DEXNEW(MODE,ISGN,POL,PNU,PAA,PNPI,JNPI) IMPLICIT NONE INCLUDE 'tauola.inc' INTEGER MODE,ISGN REAL POL(4),PNU(4),PAA(4),PNPI(4,9) INTEGER JNPI C ---------------------------------------------------------------------- * THIS SIMULATES TAU DECAY IN TAU REST FRAME * INTO NU A1, NEXT A1 DECAYS INTO RHO PI AND FINALLY RHO INTO PI PI. * OUTPUT FOUR MOMENTA: PNU TAUNEUTRINO, * PAA A1 * PIM1 PION MINUS (OR PI0) 1 (FOR TAU MINUS) * PIM2 PION MINUS (OR PI0) 2 * PIPL PION PLUS (OR PI-) * (PIPL,PIM1) FORM A RHO C ---------------------------------------------------------------------- INTEGER JDUMM REAL WT,RN REAL HV(4) INTEGER IWARM/0/ C IF(MODE.EQ.-1) THEN C =================== IWARM=1 CALL DADNEW( -1,ISGN,HV,PNU,PAA,PNPI,JDUMM) CC CALL HBOOK1(816,'WEIGHT DISTRIBUTION DEXAA $',100,-2.,2.) C ELSEIF(MODE.EQ. 0) THEN * ======================= 300 CONTINUE IF(IWARM.EQ.0) GOTO 902 CALL DADNEW( 0,ISGN,HV,PNU,PAA,PNPI,JNPI) WT=(1+POL(1)*HV(1)+POL(2)*HV(2)+POL(3)*HV(3))/2. CC CALL HFILL(816,WT) CALL TDRAND(RN,1) IF(RN.GT.WT) GOTO 300 C ELSEIF(MODE.EQ. 1) THEN * ======================= CALL DADNEW( 1,ISGN,HV,PNU,PAA,PNPI,JDUMM) CC CALL HPRINT(816) ENDIF C ===== RETURN 902 WRITE(IOUT, 9020) 9020 FORMAT(' ----- DEXNEW: LACK OF INITIALISATION') STOP END SUBROUTINE DEXPI(MODE,ISGN,POL,PPI,PNU) IMPLICIT NONE INTEGER MODE,ISGN REAL POL(4),PPI(4),PNU(4) C ---------------------------------------------------------------------- C TAU DECAY INTO PION AND TAU-NEUTRINO C IN TAU REST FRAME C OUTPUT FOUR MOMENTA: PNU TAUNEUTRINO, C PPI PION CHARGED C ---------------------------------------------------------------------- REAL WT,RN REAL HV(4) CC IF(MODE.EQ.-1) THEN C =================== CALL DADMPI(-1,ISGN,HV,PPI,PNU) CC CALL HBOOK1(815,'WEIGHT DISTRIBUTION DEXPI $',100,0,2) ELSEIF(MODE.EQ. 0) THEN C ======================= 300 CONTINUE CALL DADMPI( 0,ISGN,HV,PPI,PNU) WT=(1+POL(1)*HV(1)+POL(2)*HV(2)+POL(3)*HV(3))/2. CC CALL HFILL(815,WT) CALL TDRAND(RN,1) IF(RN.GT.WT) GOTO 300 C ELSEIF(MODE.EQ. 1) THEN C ======================= CALL DADMPI( 1,ISGN,HV,PPI,PNU) CC CALL HPRINT(815) ENDIF C ===== RETURN END SUBROUTINE DEXRO(MODE,ISGN,POL,PNU,PRO,PIC,PIZ) IMPLICIT NONE INCLUDE 'tauola.inc' INTEGER MODE,ISGN REAL POL(4),PNU(4),PRO(4),PIC(4),PIZ(4) C ---------------------------------------------------------------------- C THIS SIMULATES TAU DECAY IN TAU REST FRAME C INTO NU RHO, NEXT RHO DECAYS INTO PION PAIR. C OUTPUT FOUR MOMENTA: PNU TAUNEUTRINO, C PRO RHO C PIC PION CHARGED C PIZ PION ZERO C ---------------------------------------------------------------------- REAL WT,RN REAL HV(4) INTEGER IWARM/0/ C IF(MODE.EQ.-1) THEN C =================== IWARM=1 CALL DADMRO( -1,ISGN,HV,PNU,PRO,PIC,PIZ) CC CALL HBOOK1(816,'WEIGHT DISTRIBUTION DEXRO $',100,0,2) CC CALL HBOOK1(916,'ABS2 OF HV IN ROUTINE DEXRO $',100,0,2) C ELSEIF(MODE.EQ. 0) THEN C ======================= 300 CONTINUE IF(IWARM.EQ.0) GOTO 902 CALL DADMRO( 0,ISGN,HV,PNU,PRO,PIC,PIZ) WT=(1+POL(1)*HV(1)+POL(2)*HV(2)+POL(3)*HV(3))/2. CC CALL HFILL(816,WT) CC XHELP=HV(1)**2+HV(2)**2+HV(3)**2 CC CALL HFILL(916,XHELP) CALL TDRAND(RN,1) IF(RN.GT.WT) GOTO 300 C ELSEIF(MODE.EQ. 1) THEN C ======================= CALL DADMRO( 1,ISGN,HV,PNU,PRO,PIC,PIZ) CC CALL HPRINT(816) CC CALL HPRINT(916) ENDIF C ===== RETURN 902 WRITE(IOUT, 9020) 9020 FORMAT(' ----- DEXRO: LACK OF INITIALISATION') STOP END SUBROUTINE DPH4PI(DGAMT,HV,PN,PAA,PMULT,JNPI) IMPLICIT NONE INCLUDE 'tauola.inc' REAL DGAMT,HV(4),PN(4),PAA(4),PMULT(4,9) INTEGER JNPI C ---------------------------------------------------------------------- * IT SIMULATES A1 DECAY IN TAU REST FRAME WITH * Z-AXIS ALONG A1 MOMENTUM C ---------------------------------------------------------------------- INTEGER I,K REAL PHSPAC,PHSP,PREZ,AMP1,AMP2,AMP3,AMP4,AMRX,GAMRX,AMROP REAL GAMROP,PROB1,PROB2,PROB3,AMRB,GAMRB REAL RR1,RR2,RR3,RR4,RR5,AMS1,AMS2 REAL ALP1,ALP2,ALP,AM4SQ,AM4,AM3SQ,AM3,AM2SQ,AM2 REAL THET,PHI,GG,AMPLIT,EXE,PPI,PPPI,ENQ1,ENQ2 REAL PT(4),PIM1(4),PIM2(4),PIPL(4) REAL PR(4),PIZ(4) REAL RRR(9) DOUBLE PRECISION UU,FF,FF1,FF2,FF3,FF4,GG1,GG2,GG3,GG4 REAL RR integer ichan REAL PI /3.141592653589793238462643/ INTEGER ICONT /0/ REAL XLAM,X,Y,Z XLAM(X,Y,Z)=SQRT(ABS((X-Y-Z)**2-4.0*Y*Z)) C AMRO, GAMRO IS ONLY A PARAMETER FOR GETING HIGHT EFFICIENCY C C THREE BODY PHASE SPACE NORMALISED AS IN BJORKEN-DRELL C D**3 P /2E/(2PI)**3 (2PI)**4 DELTA4(SUM P) PHSPAC=1./2**23/PI**11 PHSP=1./2**5/PI**2 IF (JNPI.EQ.1) THEN PREZ=0.7 AMP1=AMPI AMP2=AMPI AMP3=AMPI AMP4=AMPIZ AMRX=0.782 GAMRX=0.0084 AMROP =1.2 GAMROP=.46 ELSE PREZ=0.0 AMP1=AMPIZ AMP2=AMPIZ AMP3=AMPIZ AMP4=AMPI AMRX=1.4 GAMRX=.6 AMROP =AMRX GAMROP=GAMRX ENDIF RR=0.3 CALL CHOICE(100+JNPI,RR,ICHAN,PROB1,PROB2,PROB3, $ AMROP,GAMROP,AMRX,GAMRX,AMRB,GAMRB) PREZ=PROB1+PROB2 C TAU MOMENTUM PT(1)=0. PT(2)=0. PT(3)=0. PT(4)=AMTAU C CALL TDRAND(RRR,9) C * MASSES OF 4, 3 AND 2 PI SYSTEMS C 3 PI WITH SAMPLING FOR RESONANCE CAM RR1=RRR(6) AMS1=(AMP1+AMP2+AMP3+AMP4)**2 AMS2=(AMTAU-AMNUTA)**2 ALP1=ATAN((AMS1-AMROP**2)/AMROP/GAMROP) ALP2=ATAN((AMS2-AMROP**2)/AMROP/GAMROP) ALP=ALP1+RR1*(ALP2-ALP1) AM4SQ =AMROP**2+AMROP*GAMROP*TAN(ALP) AM4 =SQRT(AM4SQ) PHSPAC=PHSPAC* $ ((AM4SQ-AMROP**2)**2+(AMROP*GAMROP)**2)/(AMROP*GAMROP) PHSPAC=PHSPAC*(ALP2-ALP1) C RR1=RRR(1) AMS1=(AMP2+AMP3+AMP4)**2 AMS2=(AM4-AMP1)**2 IF (RRR(9).GT.PREZ) THEN AM3SQ=AMS1+ RR1*(AMS2-AMS1) AM3 =SQRT(AM3SQ) C --- this part of jacobian will be recovered later FF1=AMS2-AMS1 ELSE * PHASE SPACE WITH SAMPLING FOR OMEGA RESONANCE, ALP1=ATAN((AMS1-AMRX**2)/AMRX/GAMRX) ALP2=ATAN((AMS2-AMRX**2)/AMRX/GAMRX) ALP=ALP1+RR1*(ALP2-ALP1) AM3SQ =AMRX**2+AMRX*GAMRX*TAN(ALP) AM3 =SQRT(AM3SQ) C --- THIS PART OF THE JACOBIAN WILL BE RECOVERED LATER --------------- FF1=((AM3SQ-AMRX**2)**2+(AMRX*GAMRX)**2)/(AMRX*GAMRX) FF1=FF1*(ALP2-ALP1) ENDIF C MASS OF 2 RR2=RRR(2) AMS1=(AMP3+AMP4)**2 AMS2=(AM3-AMP2)**2 * FLAT PHASE SPACE; AM2SQ=AMS1+ RR2*(AMS2-AMS1) AM2 =SQRT(AM2SQ) C --- this part of jacobian will be recovered later FF2=(AMS2-AMS1) * 2 RESTFRAME, DEFINE PIZ AND PIPL ENQ1=(AM2SQ-AMP3**2+AMP4**2)/(2*AM2) ENQ2=(AM2SQ+AMP3**2-AMP4**2)/(2*AM2) PPI= ENQ1**2-AMP4**2 PPPI=SQRT(ABS(ENQ1**2-AMP4**2)) PHSPAC=PHSPAC*(4*PI)*(2*PPPI/AM2) * PIZ MOMENTUM IN 2 REST FRAME CALL SPHERA(PPPI,PIZ) PIZ(4)=ENQ1 * PIPL MOMENTUM IN 2 REST FRAME DO 30 I=1,3 30 PIPL(I)=-PIZ(I) PIPL(4)=ENQ2 * 3 REST FRAME, DEFINE PIM1 * PR MOMENTUM PR(1)=0 PR(2)=0 PR(4)=1./(2*AM3)*(AM3**2+AM2**2-AMP2**2) PR(3)= SQRT(ABS(PR(4)**2-AM2**2)) PPI = PR(4)**2-AM2**2 * PIM1 MOMENTUM PIM1(1)=0 PIM1(2)=0 PIM1(4)=1./(2*AM3)*(AM3**2-AM2**2+AMP2**2) PIM1(3)=-PR(3) C --- this part of jacobian will be recovered later FF3=(4*PI)*(2*PR(3)/AM3) * OLD PIONS BOOSTED FROM 2 REST FRAME TO 3 REST FRAME EXE=(PR(4)+PR(3))/AM2 CALL BOSTR3(EXE,PIZ,PIZ) CALL BOSTR3(EXE,PIPL,PIPL) RR3=RRR(3) RR4=RRR(4) THET =ACOS(-1.+2*RR3) PHI = 2*PI*RR4 CALL ROTPOL(THET,PHI,PIPL) CALL ROTPOL(THET,PHI,PIM1) CALL ROTPOL(THET,PHI,PIZ) CALL ROTPOL(THET,PHI,PR) * 4 REST FRAME, DEFINE PIM2 * PR MOMENTUM PR(1)=0 PR(2)=0 PR(4)=1./(2*AM4)*(AM4**2+AM3**2-AMP1**2) PR(3)= SQRT(ABS(PR(4)**2-AM3**2)) PPI = PR(4)**2-AM3**2 * PIM2 MOMENTUM PIM2(1)=0 PIM2(2)=0 PIM2(4)=1./(2*AM4)*(AM4**2-AM3**2+AMP1**2) PIM2(3)=-PR(3) C --- this part of jacobian will be recovered later FF4=(4*PI)*(2*PR(3)/AM4) * OLD PIONS BOOSTED FROM 3 REST FRAME TO 4 REST FRAME EXE=(PR(4)+PR(3))/AM3 CALL BOSTR3(EXE,PIZ,PIZ) CALL BOSTR3(EXE,PIPL,PIPL) CALL BOSTR3(EXE,PIM1,PIM1) RR3=RRR(7) RR4=RRR(8) THET =ACOS(-1.+2*RR3) PHI = 2*PI*RR4 CALL ROTPOL(THET,PHI,PIPL) CALL ROTPOL(THET,PHI,PIM1) CALL ROTPOL(THET,PHI,PIM2) CALL ROTPOL(THET,PHI,PIZ) CALL ROTPOL(THET,PHI,PR) C * NOW TO THE TAU REST FRAME, DEFINE PAA AND NEUTRINO MOMENTA * PAA MOMENTUM PAA(1)=0 PAA(2)=0 PAA(4)=1./(2*AMTAU)*(AMTAU**2-AMNUTA**2+AM4**2) PAA(3)= SQRT(ABS(PAA(4)**2-AM4**2)) PPI = PAA(4)**2-AM4**2 PHSPAC=PHSPAC*(4*PI)*(2*PAA(3)/AMTAU) PHSP=PHSP*(4*PI)*(2*PAA(3)/AMTAU) * TAU-NEUTRINO MOMENTUM PN(1)=0 PN(2)=0 PN(4)=1./(2*AMTAU)*(AMTAU**2+AMNUTA**2-AM4**2) PN(3)=-PAA(3) C WE INCLUDE REMAINING PART OF THE JACOBIAN C --- FLAT CHANNEL AM3SQ=(PIM1(4)+PIZ(4)+PIPL(4))**2-(PIM1(3)+PIZ(3)+PIPL(3))**2 $ -(PIM1(2)+PIZ(2)+PIPL(2))**2-(PIM1(1)+PIZ(1)+PIPL(1))**2 AMS2=(AM4-AMP2)**2 AMS1=(AMP1+AMP3+AMP4)**2 FF1=(AMS2-AMS1) AMS1=(AMP3+AMP4)**2 AMS2=(SQRT(AM3SQ)-AMP1)**2 FF2=AMS2-AMS1 FF3=(4*PI)*(XLAM(AM2**2,AMP1**2,AM3SQ)/AM3SQ) FF4=(4*PI)*(XLAM(AM3SQ,AMP2**2,AM4**2)/AM4**2) UU=FF1*FF2*FF3*FF4 C --- FIRST CHANNEL AM3SQ=(PIM1(4)+PIZ(4)+PIPL(4))**2-(PIM1(3)+PIZ(3)+PIPL(3))**2 $ -(PIM1(2)+PIZ(2)+PIPL(2))**2-(PIM1(1)+PIZ(1)+PIPL(1))**2 AMS2=(AM4-AMP2)**2 AMS1=(AMP1+AMP3+AMP4)**2 ALP1=ATAN((AMS1-AMRX**2)/AMRX/GAMRX) ALP2=ATAN((AMS2-AMRX**2)/AMRX/GAMRX) FF1=((AM3SQ-AMRX**2)**2+(AMRX*GAMRX)**2)/(AMRX*GAMRX) FF1=FF1*(ALP2-ALP1) AMS1=(AMP3+AMP4)**2 AMS2=(SQRT(AM3SQ)-AMP1)**2 FF2=AMS2-AMS1 FF3=(4*PI)*(XLAM(AM2**2,AMP1**2,AM3SQ)/AM3SQ) FF4=(4*PI)*(XLAM(AM3SQ,AMP2**2,AM4**2)/AM4**2) FF=FF1*FF2*FF3*FF4 C --- SECOND CHANNEL AM3SQ=(PIM2(4)+PIZ(4)+PIPL(4))**2-(PIM2(3)+PIZ(3)+PIPL(3))**2 $ -(PIM2(2)+PIZ(2)+PIPL(2))**2-(PIM2(1)+PIZ(1)+PIPL(1))**2 AMS2=(AM4-AMP1)**2 AMS1=(AMP2+AMP3+AMP4)**2 ALP1=ATAN((AMS1-AMRX**2)/AMRX/GAMRX) ALP2=ATAN((AMS2-AMRX**2)/AMRX/GAMRX) GG1=((AM3SQ-AMRX**2)**2+(AMRX*GAMRX)**2)/(AMRX*GAMRX) GG1=GG1*(ALP2-ALP1) AMS1=(AMP3+AMP4)**2 AMS2=(SQRT(AM3SQ)-AMP2)**2 GG2=AMS2-AMS1 GG3=(4*PI)*(XLAM(AM2**2,AMP2**2,AM3SQ)/AM3SQ) GG4=(4*PI)*(XLAM(AM3SQ,AMP1**2,AM4**2)/AM4**2) GG=GG1*GG2*GG3*GG4 C --- JACOBIAN AVERAGED OVER THE TWO IF ( ( (FF+GG)*UU+FF*GG ).GT.0.0D0) THEN RR=FF*GG*UU/(0.5*PREZ*(FF+GG)*UU+(1.0-PREZ)*FF*GG) PHSPAC=PHSPAC*RR ELSE PHSPAC=0.0 ENDIF * MOMENTA OF THE TWO PI-MINUS ARE RANDOMLY SYMMETRISED IF (JNPI.EQ.1) THEN RR5= RRR(5) IF(RR5.LE.0.5) THEN DO 70 I=1,4 X=PIM1(I) PIM1(I)=PIM2(I) 70 PIM2(I)=X ENDIF PHSPAC=PHSPAC/2. ELSE C MOMENTA OF PI0'S ARE GENERATED UNIFORMLY ONLY IF PREZ=0.0 RR5= RRR(5) IF(RR5.LE.0.5) THEN DO 71 I=1,4 X=PIM1(I) PIM1(I)=PIM2(I) 71 PIM2(I)=X ENDIF PHSPAC=PHSPAC/6. ENDIF * ALL PIONS BOOSTED FROM 4 REST FRAME TO TAU REST FRAME * Z-AXIS ANTIPARALLEL TO NEUTRINO MOMENTUM EXE=(PAA(4)+PAA(3))/AM4 CALL BOSTR3(EXE,PIZ,PIZ) CALL BOSTR3(EXE,PIPL,PIPL) CALL BOSTR3(EXE,PIM1,PIM1) CALL BOSTR3(EXE,PIM2,PIM2) CALL BOSTR3(EXE,PR,PR) C PARTIAL WIDTH CONSISTS OF PHASE SPACE AND AMPLITUDE C CHECK ON CONSISTENCY WITH DADNPI, THEN, CODE BREAKES UNIFORM PION C DISTRIBUTION IN HADRONIC SYSTEM CAM Assume neutrino mass=0. and sum over final polarisation C AMX2=AM4**2 C BRAK= 2*(AMTAU**2-AMX2) * (AMTAU**2+2.*AMX2) C AMPLIT=CCABIB**2*GFERMI**2/2. * BRAK * AMX2*SIGEE(AMX2,1) IF (JNPI.EQ.1) THEN CALL DAM4PI(JNPI,PT,PN,PIM1,PIM2,PIZ,PIPL,AMPLIT,HV) ELSEIF (JNPI.EQ.2) THEN CALL DAM4PI(JNPI,PT,PN,PIM1,PIM2,PIPL,PIZ,AMPLIT,HV) ENDIF DGAMT=1/(2.*AMTAU)*AMPLIT*PHSPAC C PHASE SPACE CHECK C DGAMT=PHSPAC DO 77 K=1,4 PMULT(K,1)=PIM1(K) PMULT(K,2)=PIM2(K) PMULT(K,3)=PIPL(K) PMULT(K,4)=PIZ (K) 77 CONTINUE END SUBROUTINE DPH5PI(DGAMT,HV,PN,PAA,PMULT,JNPI) IMPLICIT NONE INCLUDE 'tauola.inc' REAL DGAMT,HV(4),PN(4),PAA(4),PMULT(4,9) INTEGER JNPI C ---------------------------------------------------------------------- * IT SIMULATES 5pi DECAY IN TAU REST FRAME WITH * Z-AXIS ALONG 5pi MOMENTUM C ---------------------------------------------------------------------- INTEGER I,K REAL PHSPAC,RR1,RR2,RR3,RR4,AM2,ENQ1,ENQ2,PPI,PPPI,EXE,THET,PHI REAL PXQ,PXN,QXN,BRAK,FOMPP,FNORM,AMPLIT REAL PT(4) REAL PR(4),PI1(4),PI2(4),PI3(4),PI4(4),PI5(4) DOUBLE PRECISION AMP1,AMP2,AMP3,AMP4,AMP5,ams1,ams2,amom,gamom DOUBLE PRECISION AM5SQ,AM4SQ,AM3SQ,AM2SQ,AM5,AM4,AM3 REAL RRR(10) DOUBLE PRECISION gg1,gg2,gg3,ff1,ff2,ff3,ff4,alp,alp1,alp2 REAL PI /3.141592653589793238462643/ INTEGER ICONT /0/ REAL fpi /93.3e-3/ REAL DCDMAS EXTERNAL DCDMAS c COMPLEX BWIGN DOUBLE PRECISION XM,AM,GAMMA BWIGN(XM,AM,GAMMA)=XM**2/CMPLX(XM**2-AM**2,GAMMA*AM) C AMOM=.782 GAMOM=0.0085 c C 6 BODY PHASE SPACE NORMALISED AS IN BJORKEN-DRELL C D**3 P /2E/(2PI)**3 (2PI)**4 DELTA4(SUM P) PHSPAC=1./2**29/PI**14 c PHSPAC=1./2**5/PI**2 C init 5pi decay mode (JNPI) AMP1=DCDMAS(IDFFIN(1,JNPI)) AMP2=DCDMAS(IDFFIN(2,JNPI)) AMP3=DCDMAS(IDFFIN(3,JNPI)) AMP4=DCDMAS(IDFFIN(4,JNPI)) AMP5=DCDMAS(IDFFIN(5,JNPI)) c C TAU MOMENTUM PT(1)=0. PT(2)=0. PT(3)=0. PT(4)=AMTAU C CALL TDRAND(RRR,10) C c masses of 5, 4, 3 and 2 pi systems c 3 pi with sampling for omega resonance cam c mass of 5 (12345) rr1=rrr(10) ams1=(amp1+amp2+amp3+amp4+amp5)**2 ams2=(amtau-amnuta)**2 am5sq=ams1+ rr1*(ams2-ams1) am5 =sqrt(am5sq) phspac=phspac*(ams2-ams1) c c mass of 4 (2345) c flat phase space rr1=rrr(9) ams1=(amp2+amp3+amp4+amp5)**2 ams2=(am5-amp1)**2 am4sq=ams1+ rr1*(ams2-ams1) am4 =sqrt(am4sq) gg1=ams2-ams1 c c mass of 3 (234) C phase space with sampling for omega resonance rr1=rrr(1) ams1=(amp2+amp3+amp4)**2 ams2=(am4-amp5)**2 alp1=atan((ams1-amom**2)/amom/gamom) alp2=atan((ams2-amom**2)/amom/gamom) alp=alp1+rr1*(alp2-alp1) am3sq =amom**2+amom*gamom*tan(alp) am3 =sqrt(am3sq) c --- this part of the jacobian will be recovered later --------------- gg2=((am3sq-amom**2)**2+(amom*gamom)**2)/(amom*gamom) gg2=gg2*(alp2-alp1) c flat phase space; C am3sq=ams1+ rr1*(ams2-ams1) C am3 =sqrt(am3sq) c --- this part of jacobian will be recovered later C gg2=ams2-ams1 c C mass of 2 (34) rr2=rrr(2) ams1=(amp3+amp4)**2 ams2=(am3-amp2)**2 c flat phase space; am2sq=ams1+ rr2*(ams2-ams1) am2 =sqrt(am2sq) c --- this part of jacobian will be recovered later gg3=ams2-ams1 c c (34) restframe, define pi3 and pi4 enq1=(am2sq+amp3**2-amp4**2)/(2*am2) enq2=(am2sq-amp3**2+amp4**2)/(2*am2) ppi= enq1**2-amp3**2 pppi=sqrt(abs(enq1**2-amp3**2)) ff1=(4*pi)*(2*pppi/am2) c pi3 momentum in (34) rest frame call sphera(pppi,pi3) pi3(4)=enq1 c pi4 momentum in (34) rest frame do 30 i=1,3 30 pi4(i)=-pi3(i) pi4(4)=enq2 c c (234) rest frame, define pi2 c pr momentum pr(1)=0 pr(2)=0 pr(4)=1./(2*am3)*(am3**2+am2**2-amp2**2) pr(3)= sqrt(abs(pr(4)**2-am2**2)) ppi = pr(4)**2-am2**2 c pi2 momentum pi2(1)=0 pi2(2)=0 pi2(4)=1./(2*am3)*(am3**2-am2**2+amp2**2) pi2(3)=-pr(3) c --- this part of jacobian will be recovered later ff2=(4*pi)*(2*pr(3)/am3) c old pions boosted from 2 rest frame to 3 rest frame exe=(pr(4)+pr(3))/am2 call bostr3(exe,pi3,pi3) call bostr3(exe,pi4,pi4) rr3=rrr(3) rr4=rrr(4) thet =acos(-1.+2*rr3) phi = 2*pi*rr4 call rotpol(thet,phi,pi2) call rotpol(thet,phi,pi3) call rotpol(thet,phi,pi4) C C (2345) rest frame, define pi5 c pr momentum pr(1)=0 pr(2)=0 pr(4)=1./(2*am4)*(am4**2+am3**2-amp5**2) pr(3)= sqrt(abs(pr(4)**2-am3**2)) ppi = pr(4)**2-am3**2 c pi5 momentum pi5(1)=0 pi5(2)=0 pi5(4)=1./(2*am4)*(am4**2-am3**2+amp5**2) pi5(3)=-pr(3) c --- this part of jacobian will be recovered later ff3=(4*pi)*(2*pr(3)/am4) c old pions boosted from 3 rest frame to 4 rest frame exe=(pr(4)+pr(3))/am3 call bostr3(exe,pi2,pi2) call bostr3(exe,pi3,pi3) call bostr3(exe,pi4,pi4) rr3=rrr(5) rr4=rrr(6) thet =acos(-1.+2*rr3) phi = 2*pi*rr4 call rotpol(thet,phi,pi2) call rotpol(thet,phi,pi3) call rotpol(thet,phi,pi4) call rotpol(thet,phi,pi5) C C (12345) rest frame, define pi1 c pr momentum pr(1)=0 pr(2)=0 pr(4)=1./(2*am5)*(am5**2+am4**2-amp1**2) pr(3)= sqrt(abs(pr(4)**2-am4**2)) ppi = pr(4)**2-am4**2 c pi1 momentum pi1(1)=0 pi1(2)=0 pi1(4)=1./(2*am5)*(am5**2-am4**2+amp1**2) pi1(3)=-pr(3) c --- this part of jacobian will be recovered later ff4=(4*pi)*(2*pr(3)/am5) c old pions boosted from 4 rest frame to 5 rest frame exe=(pr(4)+pr(3))/am4 call bostr3(exe,pi2,pi2) call bostr3(exe,pi3,pi3) call bostr3(exe,pi4,pi4) call bostr3(exe,pi5,pi5) rr3=rrr(7) rr4=rrr(8) thet =acos(-1.+2*rr3) phi = 2*pi*rr4 call rotpol(thet,phi,pi1) call rotpol(thet,phi,pi2) call rotpol(thet,phi,pi3) call rotpol(thet,phi,pi4) call rotpol(thet,phi,pi5) c * now to the tau rest frame, define paa and neutrino momenta * paa momentum paa(1)=0 paa(2)=0 c paa(4)=1./(2*amtau)*(amtau**2-amnuta**2+am5**2) c paa(3)= sqrt(abs(paa(4)**2-am5**2)) c ppi = paa(4)**2-am5**2 paa(4)=1./(2*amtau)*(amtau**2-amnuta**2+am5sq) paa(3)= sqrt(abs(paa(4)**2-am5sq)) ppi = paa(4)**2-am5sq phspac=phspac*(4*pi)*(2*paa(3)/amtau) * tau-neutrino momentum pn(1)=0 pn(2)=0 pn(4)=1./(2*amtau)*(amtau**2+amnuta**2-am5**2) pn(3)=-paa(3) c phspac=phspac * gg1*gg2*gg3*ff1*ff2*ff3*ff4 c C all pions boosted from 5 rest frame to tau rest frame C z-axis antiparallel to neutrino momentum exe=(paa(4)+paa(3))/am5 call bostr3(exe,pi1,pi1) call bostr3(exe,pi2,pi2) call bostr3(exe,pi3,pi3) call bostr3(exe,pi4,pi4) call bostr3(exe,pi5,pi5) c C partial width consists of phase space and amplitude C AMPLITUDE (cf YS.Tsai Phys.Rev.D4,2821(1971) C or F.Gilman SH.Rhie Phys.Rev.D31,1066(1985) C PXQ=AMTAU*PAA(4) PXN=AMTAU*PN(4) QXN=PAA(4)*PN(4)-PAA(1)*PN(1)-PAA(2)*PN(2)-PAA(3)*PN(3) BRAK=2*(GV**2+GA**2)*(2*PXQ*QXN+AM5SQ*PXN) & -6*(GV**2-GA**2)*AMTAU*AMNUTA*AM5SQ fompp = cabs(bwign(am3,amom,gamom))**2 c normalisation factor (to some numerical undimensioned factor; c cf R.Fischer et al ZPhys C3, 313 (1980)) fnorm = 1/fpi**6 c AMPLIT=CCABIB**2*GFERMI**2/2. * BRAK * AM5SQ*SIGEE(AM5SQ,JNPI) AMPLIT=CCABIB**2*GFERMI**2/2. * BRAK amplit = amplit * fompp * fnorm c phase space test c amplit = amplit * fnorm DGAMT=1/(2.*AMTAU)*AMPLIT*PHSPAC c ignore spin terms DO 40 I=1,3 40 HV(I)=0. c do 77 k=1,4 pmult(k,1)=pi1(k) pmult(k,2)=pi2(k) pmult(k,3)=pi3(k) pmult(k,4)=pi4(k) pmult(k,5)=pi5(k) 77 continue return C missing: transposition of identical particles, startistical factors C for identical matrices, polarimetric vector. Matrix element rather naive. C flat phase space in pion system + with breit wigner for omega C anyway it is better than nothing, and code is improvable. end SUBROUTINE DPHNPI(DGAMT,HVX,PNX,PRX,PPIX,JNPI) IMPLICIT NONE INCLUDE 'tauola.inc' REAL DGAMT,HVX(4),PNX(4),PRX(4),PPIX(4,9) INTEGER JNPI C ---------------------------------------------------------------------- C IT SIMULATES MULTIPI DECAY IN TAU REST FRAME WITH C Z-AXIS OPPOSITE TO NEUTRINO MOMENTUM C ---------------------------------------------------------------------- INTEGER I,J,K,L,KK,IL,JL,NCONT,ND REAL PS,PHSPAC,RR1,AMX2,AMW,PXQ,PXN,QXN,BRAK,AMPLIT,RN,XNPI C DOUBLE PRECISION PN(4),PR(4),PPI(4,9),HV(4) DOUBLE PRECISION PV(5,9),PT(4),UE(3),BE(3) DOUBLE PRECISION PAWT,AMX,AMS1,AMS2,PA,PHS,PHSMAX,PMIN,PMAX DOUBLE PRECISION GAM,BEP,PHI,A,B,C DOUBLE PRECISION AMPIK REAL RRR(9),RRX(2) C REAL PI /3.141592653589793238462643/ DOUBLE PRECISION WETMAX(20) DATA WETMAX /20*1D-15/ REAL SIGEE,DCDMAS EXTERNAL SIGEE,DCDMAS C CC-- PAWT(A,B,C)=SQRT((A**2-(B+C)**2)*(A**2-(B-C)**2))/(2.*A) C PAWT(A,B,C)= $ SQRT(MAX(0.D0,(A**2-(B+C)**2)*(A**2-(B-C)**2)))/(2.D0*A) C AMPIK(I,J)=DCDMAS(IDFFIN(I,J)) C C IF ((JNPI.LE.0).OR.JNPI.GT.20) THEN WRITE(6,*) 'JNPI OUTSIDE RANGE DEFINED BY WETMAX; JNPI=',JNPI STOP ENDIF C TAU MOMENTUM PT(1)=0. PT(2)=0. PT(3)=0. PT(4)=AMTAU C 500 CONTINUE C MASS OF VIRTUAL W ND=MULPIK(JNPI) PS=0. PHSPAC = 1./2.**5 /PI**2 DO 4 I=1,ND 4 PS =PS+AMPIK(I,JNPI) CALL TDRAND(RR1,1) AMS1=PS**2 AMS2=(AMTAU-AMNUTA)**2 C C AMX2=AMS1+ RR1*(AMS2-AMS1) AMX =SQRT(AMX2) AMW =AMX PHSPAC=PHSPAC * (AMS2-AMS1) C C TAU-NEUTRINO MOMENTUM PN(1)=0 PN(2)=0 PN(4)=1./(2*AMTAU)*(AMTAU**2+AMNUTA**2-AMX2) PN(3)=-SQRT((PN(4)-AMNUTA)*(PN(4)+AMNUTA)) C W MOMENTUM PR(1)=0 PR(2)=0 PR(4)=1./(2*AMTAU)*(AMTAU**2-AMNUTA**2+AMX2) PR(3)=-PN(3) PHSPAC=PHSPAC * (4.*PI) * (2.*PR(3)/AMTAU) C C AMPLITUDE (cf YS.Tsai Phys.Rev.D4,2821(1971) C or F.Gilman SH.Rhie Phys.Rev.D31,1066(1985) C PXQ=AMTAU*PR(4) PXN=AMTAU*PN(4) QXN=PR(4)*PN(4)-PR(1)*PN(1)-PR(2)*PN(2)-PR(3)*PN(3) C HERE WAS AN ERROR. 20.10.91 (ZW) C BRAK=2*(GV**2+GA**2)*(2*PXQ*PXN+AMX2*QXN) BRAK=2*(GV**2+GA**2)*(2*PXQ*QXN+AMX2*PXN) & -6*(GV**2-GA**2)*AMTAU*AMNUTA*AMX2 CAM Assume neutrino mass=0. and sum over final polarisation C BRAK= 2*(AMTAU**2-AMX2) * (AMTAU**2+2.*AMX2) AMPLIT=CCABIB**2*GFERMI**2/2. * BRAK * AMX2*SIGEE(AMX2,JNPI) DGAMT=1./(2.*AMTAU)*AMPLIT*PHSPAC C C ISOTROPIC W DECAY IN W REST FRAME PHSMAX = 1. DO 200 I=1,4 200 PV(I,1)=PR(I) PV(5,1)=AMW PV(5,ND)=AMPIK(ND,JNPI) C COMPUTE MAX. PHASE SPACE FACTOR PMAX=AMW-PS+AMPIK(ND,JNPI) PMIN=.0 DO 220 IL=ND-1,1,-1 PMAX=PMAX+AMPIK(IL,JNPI) PMIN=PMIN+AMPIK(IL+1,JNPI) 220 PHSMAX=PHSMAX*PAWT(PMAX,PMIN,AMPIK(IL,JNPI))/PMAX C --- 2.02.94 ZW 9 lines AMX=AMW DO 222 IL=1,ND-2 AMS1=.0 DO 223 JL=IL+1,ND 223 AMS1=AMS1+AMPIK(JL,JNPI) AMS1=AMS1**2 AMX =(AMX-AMPIK(IL,JNPI)) AMS2=(AMX)**2 PHSMAX=PHSMAX * (AMS2-AMS1) 222 CONTINUE NCONT=0 100 CONTINUE NCONT=NCONT+1 CAM GENERATE ND-2 EFFECTIVE MASSES PHS=1.D0 PHSPAC = 1./2.**(6*ND-7) /PI**(3*ND-4) AMX=AMW CALL TDRAND(RRR,ND-2) DO 230 IL=1,ND-2 AMS1=.0D0 DO 231 JL=IL+1,ND 231 AMS1=AMS1+AMPIK(JL,JNPI) AMS1=AMS1**2 AMS2=(AMX-AMPIK(IL,JNPI))**2 RR1=RRR(IL) AMX2=AMS1+ RR1*(AMS2-AMS1) AMX=SQRT(AMX2) PV(5,IL+1)=AMX PHSPAC=PHSPAC * (AMS2-AMS1) C --- 2.02.94 ZW 1 line PHS=PHS* (AMS2-AMS1) PA=PAWT(PV(5,IL),PV(5,IL+1),AMPIK(IL,JNPI)) PHS =PHS *PA/PV(5,IL) 230 CONTINUE PA=PAWT(PV(5,ND-1),AMPIK(ND-1,JNPI),AMPIK(ND,JNPI)) PHS =PHS *PA/PV(5,ND-1) CALL TDRAND(RN,1) WETMAX(JNPI)=1.2D0*MAX(WETMAX(JNPI)/1.2D0,PHS/PHSMAX) IF (NCONT.EQ.500 000) THEN XNPI=0.0 DO KK=1,ND XNPI=XNPI+AMPIK(KK,JNPI) ENDDO WRITE(6,*) 'ROUNDING INSTABILITY IN DPHNPI ?' WRITE(6,*) 'AMW=',AMW,'XNPI=',XNPI WRITE(6,*) 'IF =AMW= IS NEARLY EQUAL =XNPI= THAT IS IT' WRITE(6,*) 'PHS=',PHS,'PHSMAX=',PHSMAX GOTO 500 ENDIF IF(RN*PHSMAX*WETMAX(JNPI).GT.PHS) GO TO 100 C...PERFORM SUCCESSIVE TWO-PARTICLE DECAYS IN RESPECTIVE CM FRAME 280 DO 300 IL=1,ND-1 PA=PAWT(PV(5,IL),PV(5,IL+1),AMPIK(IL,JNPI)) CALL TDRAND(RRX,2) UE(3)=2.*RRX(1)-1. PHI=2.*PI*RRX(2) UE(1)=SQRT(1.D0-UE(3)**2)*COS(PHI) UE(2)=SQRT(1.D0-UE(3)**2)*SIN(PHI) DO 290 J=1,3 PPI(J,IL)=PA*UE(J) 290 PV(J,IL+1)=-PA*UE(J) PPI(4,IL)=SQRT(PA**2+AMPIK(IL,JNPI)**2) PV(4,IL+1)=SQRT(PA**2+PV(5,IL+1)**2) PHSPAC=PHSPAC *(4.*PI)*(2.*PA/PV(5,IL)) 300 CONTINUE C...LORENTZ TRANSFORM DECAY PRODUCTS TO TAU FRAME DO 310 J=1,4 310 PPI(J,ND)=PV(J,ND) DO 340 IL=ND-1,1,-1 DO 320 J=1,3 320 BE(J)=PV(J,IL)/PV(4,IL) GAM=PV(4,IL)/PV(5,IL) DO 340 I=IL,ND BEP=BE(1)*PPI(1,I)+BE(2)*PPI(2,I)+BE(3)*PPI(3,I) DO 330 J=1,3 330 PPI(J,I)=PPI(J,I)+GAM*(GAM*BEP/(1.D0+GAM)+PPI(4,I))*BE(J) PPI(4,I)=GAM*(PPI(4,I)+BEP) 340 CONTINUE C HV(4)=1. HV(3)=0. HV(2)=0. HV(1)=0. DO K=1,4 PNX(K)=PN(K) PRX(K)=PR(K) HVX(K)=HV(K) DO L=1,ND PPIX(K,L)=PPI(K,L) ENDDO ENDDO RETURN END SUBROUTINE DPHSAA(DGAMT,HV,PN,PAA,PIM1,PIM2,PIPL,JAA) IMPLICIT NONE INCLUDE 'tauola.inc' REAL DGAMT,HV(4),PN(4),PAA(4),PIM1(4),PIM2(4),PIPL(4) INTEGER JAA C ---------------------------------------------------------------------- * IT SIMULATES A1 DECAY IN TAU REST FRAME WITH * Z-AXIS ALONG A1 MOMENTUM C ---------------------------------------------------------------------- INTEGER MNUM,KEYT REAL RMOD,AMP1,AMP2,AMP3 REAL*4 RRR(1) C MATRIX ELEMENT NUMBER: MNUM=0 C TYPE OF THE GENERATION: KEYT=1 CALL TDRAND(RRR,1) RMOD=RRR(1) IF (RMOD.LT.BRA1) THEN JAA=1 AMP1=AMPI AMP2=AMPI AMP3=AMPI ELSE JAA=2 AMP1=AMPIZ AMP2=AMPIZ AMP3=AMPI ENDIF CALL $ DPHTRE(DGAMT,HV,PN,PAA,PIM1,AMP1,PIM2,AMP2,PIPL,AMP3,KEYT,MNUM) END SUBROUTINE DPHSEL(DGAMX,HVX,XNX,PAAX,QPX,XAX,PHX) IMPLICIT NONE REAL DGAMX,HVX(4),XNX(4),PAAX(4),QPX(4),XAX(4),PHX(4) C XNX,XNA was flipped in parameters of dphsel and dphsmu C ********************************************************************* C * ELECTRON DECAY MODE * C ********************************************************************* INTEGER K,IELMU DOUBLE PRECISION HV(4),PH(4),PAA(4),XA(4),QP(4),XN(4) DOUBLE PRECISION DGAMT IELMU=1 CALL DRCMU(DGAMT,HV,PH,PAA,XA,QP,XN,IELMU) DO 7 K=1,4 HVX(K)=HV(K) PHX(K)=PH(K) PAAX(K)=PAA(K) XAX(K)=XA(K) QPX(K)=QP(K) XNX(K)=XN(K) 7 CONTINUE DGAMX=DGAMT END SUBROUTINE DPHSKS(DGAMT,HV,PN,PKS,PKK,PPI,JKST) IMPLICIT NONE INCLUDE 'tauola.inc' REAL DGAMT,HV(4),PN(4),PKS(4),PKK(4),PPI(4) INTEGER JKST C ---------------------------------------------------------------------- C IT SIMULATES KAON* DECAY IN TAU REST FRAME WITH C Z-AXIS ALONG KAON* MOMENTUM C JKST=10 FOR K* --->K0 + PI+- C JKST=20 FOR K* --->K+- + PI0 C ---------------------------------------------------------------------- INTEGER I REAL PHSPAC,RR1,AMS1,AMS2,ALP,ALP1,ALP2,AMX,AMX2,ENPI,PPPI REAL EXE,PKSD,QQPKS,PRODPQ,PRODNQ,PRODPN,QQ2,BRAK,FKS,AMPLIT REAL PT(4),QQ(4) COMPLEX BWIGS REAL PI /3.141592653589793238462643/ C INTEGER ICONT /0/ C THREE BODY PHASE SPACE NORMALISED AS IN BJORKEN-DRELL PHSPAC=1./2**11/PI**5 C TAU MOMENTUM PT(1)=0. PT(2)=0. PT(3)=0. PT(4)=AMTAU CALL TDRAND(RR1,1) C HERE BEGIN THE K0,PI+_ DECAY IF(JKST.EQ.10)THEN C ================== C MASS OF (REAL/VIRTUAL) K* AMS1=(AMPI+AMKZ)**2 AMS2=(AMTAU-AMNUTA)**2 C FLAT PHASE SPACE C AMX2=AMS1+ RR1*(AMS2-AMS1) C AMX=SQRT(AMX2) C PHSPAC=PHSPAC*(AMS2-AMS1) C PHASE SPACE WITH SAMPLING FOR K* RESONANCE ALP1=ATAN((AMS1-AMKST**2)/AMKST/GAMKST) ALP2=ATAN((AMS2-AMKST**2)/AMKST/GAMKST) ALP=ALP1+RR1*(ALP2-ALP1)