CDECK ID>, HWDCLE. *CMZ :- -28/01/92 12.34.44 by Mike Seymour *-- Author : Luca Stanco C----------------------------------------------------------------------- SUBROUTINE HWDCLE(IHEP) C----------------------------------------------------------------------- C INTERFACE TO QQ-CLEO MONTE CARLO (LS 11/12/91) C----------------------------------------------------------------------- INCLUDE 'HERWIG61.INC' INTEGER IHEP,IIHEP,NHEPHF,QQLMAT LOGICAL QQLERR CHARACTER*8 NAME EXTERNAL QQLMAT C---QQ-CLEO COMMON'S C*** MCPARS.INC INTEGER MCTRK, NTRKS, MCVRTX, NVTXS, MCHANS, MCDTRS, MPOLQQ INTEGER MCNUM, MCSTBL, MCSTAB, MCTLQQ, MDECQQ INTEGER MHLPRB, MHLLST, MHLANG, MCPLST, MFDECA PARAMETER (MCTRK = 512) PARAMETER (NTRKS = MCTRK) PARAMETER (MCVRTX = 256) PARAMETER (NVTXS = MCVRTX) PARAMETER (MCHANS = 4000) PARAMETER (MCDTRS = 8000) PARAMETER (MPOLQQ = 300) PARAMETER (MCNUM = 500) PARAMETER (MCSTBL = 40) PARAMETER (MCSTAB = 512) PARAMETER (MCTLQQ = 100) PARAMETER (MDECQQ = 300) PARAMETER (MHLPRB = 500) PARAMETER (MHLLST = 1000) PARAMETER (MHLANG = 500) PARAMETER (MCPLST = 200) PARAMETER (MFDECA = 5) C*** MCPROP.INC REAL AMASS, CHARGE, CTAU, SPIN, RWIDTH, RMASMN, RMASMX REAL RMIXPP, RCPMIX INTEGER NPMNQQ, NPMXQQ, IDMC, INVMC, LPARTY, CPARTY INTEGER IMIXPP, ICPMIX COMMON/MCMAS1/ * NPMNQQ, NPMXQQ, * AMASS(-20:MCNUM), CHARGE(-20:MCNUM), CTAU(-20:MCNUM), * IDMC(-20:MCNUM), SPIN(-20:MCNUM), * RWIDTH(-20:MCNUM), RMASMN(-20:MCNUM), RMASMX(-20:MCNUM), * LPARTY(-20:MCNUM), CPARTY(-20:MCNUM), * IMIXPP(-20:MCNUM), RMIXPP(-20:MCNUM), * ICPMIX(-20:MCNUM), RCPMIX(-20:MCNUM), * INVMC(0:MCSTBL) C INTEGER NPOLQQ, IPOLQQ COMMON/MCPOL1/ * NPOLQQ, IPOLQQ(5,MPOLQQ) C CHARACTER QNAME*10, PNAME*10 COMMON/MCNAMS/ * QNAME(37), PNAME(-20:MCNUM) C C*** MCCOMS.INC INTEGER NCTLQQ, NDECQQ, IVRSQQ, IORGQQ, IRS1QQ INTEGER IEVTQQ, IRUNQQ, IBMRAD INTEGER NTRKMC, QQNTRK, NSTBMC, NSTBQQ, NCHGMC, NCHGQQ INTEGER IRANQQ, IRANMC, IRANCC, IRS2QQ INTEGER IPFTQQ, IPCDQQ, IPRNTV, ITYPEV, IDECSV, IDAUTV INTEGER ISTBMC, NDAUTV INTEGER IVPROD, IVDECA REAL BFLDQQ REAL ENERQQ, BEAMQQ, BMPSQQ, BMNGQQ, EWIDQQ, BWPSQQ, BWNGQQ REAL BPOSQQ, BSIZQQ REAL ECM, P4CMQQ, P4PHQQ, ENERNW, BEAMNW, BEAMP, BEAMN REAL PSAV, P4QQ, HELCQQ CHARACTER DATEQQ*20, TIMEQQ*20, FOUTQQ*80, FCTLQQ*80, FDECQQ*80 CHARACTER FGEOQQ*80 CHARACTER CCTLQQ*80, CDECQQ*80 C COMMON/MCCM1A/ * NCTLQQ, NDECQQ, IVRSQQ, IORGQQ, IRS1QQ(3), BFLDQQ, * ENERQQ, BEAMQQ, BMPSQQ, BMNGQQ, EWIDQQ, BWPSQQ, BWNGQQ, * BPOSQQ(3), BSIZQQ(3), * IEVTQQ, IRUNQQ, * IBMRAD, ECM, P4CMQQ(4), P4PHQQ(4), * ENERNW, BEAMNW, BEAMP, BEAMN, * NTRKMC, QQNTRK, NSTBMC, NSTBQQ, NCHGMC, NCHGQQ, * IRANQQ(2), IRANMC(2), IRANCC(2), IRS2QQ(5), * IPFTQQ(MCTRK), IPCDQQ(MCTRK), IPRNTV(MCTRK), ITYPEV(MCTRK,2), * IDECSV(MCTRK), IDAUTV(MCTRK), ISTBMC(MCTRK), NDAUTV(MCTRK), * IVPROD(MCTRK), IVDECA(MCTRK), * PSAV(MCTRK,4), HELCQQ(MCTRK), P4QQ(4,MCTRK) C COMMON/MCCM1B/ * DATEQQ, TIMEQQ, FOUTQQ, FCTLQQ, FDECQQ, FGEOQQ, * CCTLQQ(MCTLQQ), CDECQQ(MDECQQ) INTEGER IDSTBL COMMON/MCCM1C/ * IDSTBL(MCSTAB) C INTEGER IFINAL(MCTRK), IFINSV(MCSTAB), NFINAL EQUIVALENCE (IFINAL,ISTBMC), (IFINSV,IDSTBL), (NFINAL,NSTBMC) C INTEGER NVRTX, ITRKIN, NTRKOU, ITRKOU, IVKODE REAL XVTX, TVTX, RVTX COMMON/MCCM2/ * NVRTX, XVTX(MCVRTX,3), TVTX(MCVRTX), RVTX(MCVRTX), * ITRKIN(MCVRTX), NTRKOU(MCVRTX), ITRKOU(MCVRTX), * IVKODE(MCVRTX) C*** MCGEN.INC INTEGER QQIST,QQIFR,QQN,QQK,QQMESO,QQNC,QQKC,QQLASTN REAL QQPUD,QQPS1,QQSIGM,QQMAS,QQPAR,QQCMIX,QQCND,QQBSPI,QQBSYM,QQP REAL QQPC,QQCZF C COMMON/DATA1/QQIST,QQIFR,QQPUD,QQPS1,QQSIGM,QQMAS(15),QQPAR(25) COMMON/DATA2/QQCZF(15),QQMESO(36),QQCMIX(6,2) COMMON/DATA3/QQCND(3) COMMON/DATA5/QQBSPI(5),QQBSYM(3) COMMON/JET/QQN,QQK(250,2),QQP(250,5),QQNC,QQKC(10),QQPC(10,4), * QQLASTN C--- IF(FSTEVT) THEN C---INITIALIZE QQ-CLEO CALL QQINIT(QQLERR) IF(QQLERR) CALL HWWARN('HWDEUR',500,*999) ENDIF C---CONSTRUCT THE HADRON FOR QQ-CLEO C NOTE: THE IDPDG CODE IS PROVIDED THROUGH THE QQLMAT ROUTINE C FROM THE CLEO PACKAGE (QQ-CLEO <--> IDPDG CODE TRANSFORMATION) QQN=1 IDHEP(IHEP)=IDPDG(IDHW(IHEP)) QQK(1,1)=0 QQK(1,2)=QQLMAT(IDHEP(IHEP),1) QQP(1,1)=PHEP(1,IHEP) QQP(1,2)=PHEP(2,IHEP) QQP(1,3)=PHEP(3,IHEP) QQP(1,5)=AMASS(QQK(1,2)) QQP(1,4)=SQRT(QQP(1,5)**2+QQP(1,1)**2+QQP(1,2)**2+QQP(1,3)**2) C---LET QQ-CLEO DO THE JOB QQNTRK=0 NVRTX=0 CALL DECADD(.FALSE.) C---UPDATE THE HERWIG TABLE : LOOP OVER QQN-CLEO FINAL PARTICLES DO 40 IIHEP=1,QQN NHEP=NHEP+1 ISTHEP(NHEP)=198 IF(ITYPEV(IIHEP,2).GE.0) ISTHEP(NHEP)=1 IDHEP(NHEP)=QQLMAT(ITYPEV(IIHEP,1),2) CALL HWUIDT(1,IDHEP(NHEP),IDHW(NHEP),NAME) IF(IIHEP.EQ.1) THEN ISTHEP(IHEP)=199 JDAHEP(1,IHEP)=NHEP JDAHEP(2,IHEP)=NHEP ISTHEP(NHEP)=199 NHEPHF=NHEP JMOHEP(1,NHEP)=IHEP JMOHEP(2,NHEP)=IHEP ELSE JMOHEP(1,NHEP)=IPRNTV(IIHEP)+NHEPHF-1 JMOHEP(2,NHEP)=NHEPHF ENDIF JDAHEP(1,NHEP)=0 JDAHEP(2,NHEP)=0 IF(NDAUTV(IIHEP).GT.0) THEN JDAHEP(1,NHEP)=IDAUTV(IIHEP)+NHEPHF-1 JDAHEP(2,NHEP)=JDAHEP(1,NHEP)+NDAUTV(IIHEP)-1 ENDIF PHEP(1,NHEP)=QQP(IIHEP,1) PHEP(2,NHEP)=QQP(IIHEP,2) PHEP(3,NHEP)=QQP(IIHEP,3) PHEP(4,NHEP)=QQP(IIHEP,4) PHEP(5,NHEP)=QQP(IIHEP,5) VHEP(1,NHEP)=XVTX(IVPROD(IIHEP),1) VHEP(2,NHEP)=XVTX(IVPROD(IIHEP),2) VHEP(3,NHEP)=XVTX(IVPROD(IIHEP),3) VHEP(4,NHEP)=0. 40 CONTINUE 999 END CDECK ID>, HWDEUR. *CMZ :- -28/01/92 12.34.44 by Mike Seymour *-- Author : Luca Stanco C----------------------------------------------------------------------- SUBROUTINE HWDEUR(IHEP) C----------------------------------------------------------------------- C INTERFACE TO EURODEC PACKAGE (LS 10/29/91) C----------------------------------------------------------------------- INCLUDE 'HERWIG61.INC' INTEGER IHEP,IIHEP,NHEPHF,IEUPDG,IPDGEU CHARACTER*8 NAME C---EURODEC COMMON'S : INITIAL INPUT INTEGER EULUN0,EULUN1,EULUN2,EURUN,EUEVNT CHARACTER*4 EUDATD,EUTIT REAL AMINIE(12),EUWEI COMMON/INPOUT/EULUN0,EULUN1,EULUN2 COMMON/FILNAM/EUDATD,EUTIT COMMON/HVYINI/AMINIE COMMON/RUNINF/EURUN,EUEVNT,EUWEI C---EURODEC WORKING COMMON'S INTEGER NPMAX,NTMAX PARAMETER (NPMAX=18,NTMAX=2000) INTEGER EUNP,EUIP(NPMAX),EUPHEL(NPMAX),EUTEIL,EUINDX(NTMAX), & EUORIG(NTMAX),EUDCAY(NTMAX),EUTHEL(NTMAX) REAL EUAPM(NPMAX),EUPCM(5,NPMAX),EUPVTX(3,NPMAX),EUPTEI(5,NTMAX), & EUSECV(3,NTMAX) COMMON/MOMGEN/EUNP,EUIP,EUAPM,EUPCM,EUPHEL,EUPVTX COMMON/RESULT/EUTEIL,EUPTEI,EUINDX,EUORIG,EUDCAY,EUTHEL,EUSECV C---EURODEC COMMON'S FOR DECAY PROPERTIES INTEGER NGMAX,NCMAX PARAMETER (NGMAX=400,NCMAX=9000) INTEGER EUNPA,EUIPC(NGMAX),EUIPDG(NGMAX),EUIDP(NGMAX), & EUCONV(NCMAX) REAL EUPM(NGMAX),EUPLT(NGMAX) COMMON/PCTABL/EUNPA,EUIPC,EUIPDG,EUPM,EUPLT,EUIDP COMMON/CONVRT/EUCONV C--- IF(FSTEVT) THEN C---CHANGE HERE THE DEFAULT VALUES OF EURODEC COMMON'S C C---INITIALIZE EURODEC COMMON'S CC CALL EUDCIN C---INITIALIZE EURODEC CALL EUDINI ENDIF C---CONSTRUCT THE HADRON FOR EURODEC FROM ID1,ID2 EUNP=1 IDHEP(IHEP)=IDPDG(IDHW(IHEP)) EUIP(1)=IPDGEU(IDHEP(IHEP)) EUAPM(1)=EUPM(EUCONV(IABS(EUIP(1)))) EUPCM(1,1)=PHEP(1,IHEP) EUPCM(2,1)=PHEP(2,IHEP) EUPCM(3,1)=PHEP(3,IHEP) EUPCM(5,1)=SQRT(PHEP(1,IHEP)**2+PHEP(2,IHEP)**2+PHEP(3,IHEP)**2) EUPCM(4,1)=SQRT(EUPCM(5,1)**2+EUAPM(1)**2) C NOT POLARIZED HADRONS EUPHEL(1)=0 C HADRONS START FROM PRIMARY VERTEX EUPVTX(1,1)=0. EUPVTX(2,1)=0. EUPVTX(3,1)=0. C---LET EURODEC DO THE JOB EUTEIL=0 CALL FRAGMT(1,1,0) C---UPDATE THE HERWIG TABLE : LOOP OVER N-EURODEC FINAL PARTICLES DO 40 IIHEP=1,EUTEIL NHEP=NHEP+1 ISTHEP(NHEP)=198 IF(EUDCAY(IIHEP).EQ.0) ISTHEP(NHEP)=1 IDHEP(NHEP)=IEUPDG(EUINDX(IIHEP)) CALL HWUIDT(1,IDHEP(NHEP),IDHW(NHEP),NAME) IF(IIHEP.EQ.1) THEN ISTHEP(IHEP)=199 JDAHEP(1,IHEP)=NHEP JDAHEP(2,IHEP)=NHEP ISTHEP(NHEP)=199 NHEPHF=NHEP JMOHEP(1,NHEP)=IHEP JMOHEP(2,NHEP)=IHEP JDAHEP(1,NHEP)=EUDCAY(IIHEP)/10000+NHEPHF-1 JDAHEP(2,NHEP)=MOD(EUDCAY(IIHEP),10000)+NHEPHF-1 ELSE JMOHEP(1,NHEP)=MOD(EUORIG(IIHEP),10000)+NHEPHF-1 JMOHEP(2,NHEP)=NHEPHF JDAHEP(1,NHEP)=EUDCAY(IIHEP)/10000+NHEPHF-1 JDAHEP(2,NHEP)=MOD(EUDCAY(IIHEP),10000)+NHEPHF-1 ENDIF PHEP(1,NHEP)=EUPTEI(1,IIHEP) PHEP(2,NHEP)=EUPTEI(2,IIHEP) PHEP(3,NHEP)=EUPTEI(3,IIHEP) PHEP(4,NHEP)=EUPTEI(4,IIHEP) PHEP(5,NHEP)=EUPTEI(5,IIHEP) VHEP(1,NHEP)=EUSECV(1,IIHEP) VHEP(2,NHEP)=EUSECV(2,IIHEP) VHEP(3,NHEP)=EUSECV(3,IIHEP) VHEP(4,NHEP)=0. IF (IIHEP.GT.NTMAX) CALL HWWARN('HWDEUR',99,*999) 40 CONTINUE 999 END CDECK ID>, HWDFOR. *CMZ :- -01/04/99 19.52.44 by Mike Seymour *-- Author : Ian Knowles C----------------------------------------------------------------------- SUBROUTINE HWDFOR(P0,P1,P2,P3,P4) C----------------------------------------------------------------------- C Generates 4-body decay 0->1+2+3+4 using pure phase space C----------------------------------------------------------------------- IMPLICIT NONE DOUBLE PRECISION HWR,P0(5),P1(5),P2(5),P3(5),P4(5),B,C,AA,BB, & CC,DD,EE,TT,S1,RS1,FF,S2,PP,QQ,RR,P1CM,P234(5),P2CM,P34(5),P3CM DOUBLE PRECISION TWO PARAMETER (TWO=2.D0) EXTERNAL HWR B=P0(5)-P1(5) C=P2(5)+P3(5)+P4(5) IF (B.LT.C) CALL HWWARN('HWDFOR',100,*999) AA=(P0(5)+P1(5))**2 BB=B**2 CC=C**2 DD=(P3(5)+P4(5))**2 EE=(P3(5)-P4(5))**2 TT=(B-C)*P0(5)**7/16 C Select squared masses S1 and S2 of 234 and 34 subsystems 10 S1=BB+HWR()*(CC-BB) RS1=SQRT(S1) FF=(RS1-P2(5))**2 S2=DD+HWR()*(FF-DD) PP=(AA-S1)*(BB-S1) QQ=((RS1+P2(5))**2-S2)*(FF-S2)/S1 RR=(S2-DD)*(S2-EE)/S2 IF (PP*QQ*RR*(FF-DD)**2.LT.TT*S1*S2*HWR()**2) GOTO 10 C Do two body decays: 0-->1+234, 234-->2+34 and 34-->3+4 P1CM=SQRT(PP/4)/P0(5) P234(5)=RS1 P2CM=SQRT(QQ/4) P34(5)=SQRT(S2) P3CM=SQRT(RR/4) CALL HWDTWO(P0 ,P1,P234,P1CM,TWO,.TRUE.) CALL HWDTWO(P234,P2,P34 ,P2CM,TWO,.TRUE.) CALL HWDTWO(P34 ,P3,P4 ,P3CM,TWO,.TRUE.) 999 END CDECK ID>, HWDFIV. *CMZ :- -01/04/99 19.52.44 by Mike Seymour *-- Author : Ian Knowles C----------------------------------------------------------------------- SUBROUTINE HWDFIV(P0,P1,P2,P3,P4,P5) C----------------------------------------------------------------------- C Generates 5-body decay 0->1+2+3+4+5 using pure phase space C----------------------------------------------------------------------- IMPLICIT NONE DOUBLE PRECISION HWR,P0(5),P1(5),P2(5),P3(5),P4(5),P5(5),B,C, & AA,BB,CC,DD,EE,FF,TT,S1,RS1,GG,S2,RS2,HH,S3,PP,QQ,RR,SS,P1CM, & P2345(5),P2CM,P345(5),P3CM,P45(5),P4CM DOUBLE PRECISION TWO PARAMETER (TWO=2.D0) EXTERNAL HWR B=P0(5)-P1(5) C=P2(5)+P3(5)+P4(5)+P5(5) IF (B.LT.C) CALL HWWARN('HWDFIV',100,*999) AA=(P0(5)+P1(5))**2 BB=B**2 CC=C**2 DD=(P3(5)+P4(5)+P5(5))**2 EE=(P4(5)+P5(5))**2 FF=(P4(5)-P5(5))**2 TT=(B-C)*P0(5)**11/729 C Select squared masses S1, S2 and S3 of 2345, 345 and 45 subsystems 10 S1=BB+HWR()*(CC-BB) RS1=SQRT(S1) GG=(RS1-P2(5))**2 S2=DD+HWR()*(GG-DD) RS2=SQRT(S2) HH=(RS2-P3(5))**2 S3=EE+HWR()*(HH-EE) PP=(AA-S1)*(BB-S1) QQ=((RS1+P2(5))**2-S2)*(GG-S2)/S1 RR=((RS2+P3(5))**2-S3)*(HH-S3)/S2 SS=(S3-EE)*(S3-FF)/S3 IF (PP*QQ*RR*SS*((GG-DD)*(HH-EE))**2.LT.TT*S1*S2*S3*HWR()**2) & GOTO 10 C Do two body decays: 0-->1+2345, 2345-->2+345, 345-->3+45 and 45-->4+5 P1CM=SQRT(PP/4)/P0(5) P2345(5)=RS1 P2CM=SQRT(QQ/4) P345(5)=RS2 P3CM=SQRT(RR/4) P45(5)=SQRT(S3) P4CM=SQRT(SS/4) CALL HWDTWO(P0 ,P1,P2345,P1CM,TWO,.TRUE.) CALL HWDTWO(P2345,P2,P345 ,P2CM,TWO,.TRUE.) CALL HWDTWO(P345 ,P3,P45 ,P3CM,TWO,.TRUE.) CALL HWDTWO(P45 ,P4,P5 ,P4CM,TWO,.TRUE.) 999 END CDECK ID>, HWDHAD. *CMZ :- -26/04/91 14.01.26 by Federico Carminati *-- Author : Ian Knowles, Bryan Webber & Mike Seymour C----------------------------------------------------------------------- SUBROUTINE HWDHAD C----------------------------------------------------------------------- C GENERATES DECAYS OF UNSTABLE HADRONS AND LEPTONS C----------------------------------------------------------------------- INCLUDE 'HERWIG61.INC' DOUBLE PRECISION HWR,HWULDO,RN,BF,COSANG,RSUM,DIST(4),VERTX(4), & PMIX,WTMX,WTMX2,XS,DOT1,DOT2,HWDPWT,HWDWWT,XXX,YYY INTEGER IHEP,ID,MHEP,IDM,I,IDS,IM,MO,IPDG LOGICAL STABLE EXTERNAL HWR,HWDPWT,HWDWWT,HWULDO IF (IERROR.NE.0) RETURN DO 100 IHEP=1,NMXHEP IF (IHEP.GT.NHEP) THEN ISTAT=90 RETURN ELSEIF (ISTHEP(IHEP).EQ.120 .AND. & JDAHEP(1,IHEP).EQ.IHEP.AND.JDAHEP(2,IHEP).EQ.IHEP) THEN C---COPY COLOUR SINGLET CMF NHEP=NHEP+1 IF (NHEP.GT.NMXHEP) CALL HWWARN('HWDHAD',100,*999) CALL HWVEQU(5,PHEP(1,IHEP),PHEP(1,NHEP)) CALL HWVEQU(4,VHEP(1,IHEP),VHEP(1,NHEP)) IDHW(NHEP)=IDHW(IHEP) IDHEP(NHEP)=IDHEP(IHEP) ISTHEP(NHEP)=190 JMOHEP(1,NHEP)=IHEP JMOHEP(2,NHEP)=NHEP JDAHEP(2,NHEP)=NHEP JDAHEP(1,IHEP)=NHEP JDAHEP(2,IHEP)=NHEP ELSEIF (ISTHEP(IHEP).GE.190.AND.ISTHEP(IHEP).LE.193) THEN C---FIRST CHECK FOR STABILITY ID=IDHW(IHEP) IF (RSTAB(ID)) THEN ISTHEP(IHEP)=1 JDAHEP(1,IHEP)=0 JDAHEP(2,IHEP)=0 C---SPECIAL FOR GAUGE BOSON DECAY IF (ID.GE.198.AND.ID.LE.200) CALL HWDBOS(IHEP) C---SPECIAL FOR HIGGS BOSON DECAY IF (ID.EQ.201) CALL HWDHIG(ZERO) ELSE C---UNSTABLE. C Calculate position of decay vertex IF (DKLTM(ID).EQ.ZERO) THEN CALL HWVEQU(4,VHEP(1,IHEP),VERTX) MHEP=IHEP IDM=ID ELSE CALL HWUDKL(ID,PHEP(1,IHEP),DIST) CALL HWVSUM(4,VHEP(1,IHEP),DIST,VERTX) IF (MAXDKL) THEN CALL HWDXLM(VERTX,STABLE) IF (STABLE) THEN ISTHEP(IHEP)=1 JDAHEP(1,IHEP)=0 JDAHEP(2,IHEP)=0 GOTO 100 ENDIF ENDIF IF (MIXING.AND.(ID.EQ.221.OR.ID.EQ.223.OR. & ID.EQ.245.OR.ID.EQ.247)) THEN C Select flavour of decaying b-meson allowing for flavour oscillation IDS=MOD(ID,3) XXX=XMRCT(IDS)*DIST(4)/PHEP(4,IHEP) YYY=YMRCT(IDS)*DIST(4)/PHEP(4,IHEP) IF (ABS(YYY).LT.10) THEN PMIX=HALF*(ONE-COS(XXX)/COSH(YYY)) ELSE PMIX=HALF ENDIF IF (HWR().LE.PMIX) THEN IF (ID.LE.223) THEN IDM=ID+24 ELSE IDM=ID-24 ENDIF ELSE IDM=ID ENDIF C Introduce a decaying neutral b-meson IF (NHEP+1.GT.NMXHEP) CALL HWWARN('HWDHAD',101,*999) MHEP=NHEP+1 ISTHEP(MHEP)=ISTHEP(IHEP) ISTHEP(IHEP)=200 JDAHEP(1,IHEP)=MHEP JDAHEP(2,IHEP)=MHEP IDHW(MHEP)=IDM IDHEP(MHEP)=IDPDG(IDM) JMOHEP(1,MHEP)=IHEP JMOHEP(2,MHEP)=JMOHEP(2,IHEP) CALL HWVEQU(5,PHEP(1,IHEP),PHEP(1,MHEP)) CALL HWVEQU(4,VERTX,VHEP(1,MHEP)) NHEP=NHEP+1 ELSE MHEP=IHEP IDM=ID ENDIF ENDIF C Use CLEO/EURODEC packages for b-hadrons if requested IF ((IDM.GE.221.AND.IDM.LE.231).OR. & (IDM.GE.245.AND.IDM.LE.254)) THEN IF (BDECAY.EQ.'CLEO') THEN CALL HWDCLE(MHEP) GOTO 100 ELSEIF (BDECAY.EQ.'EURO') THEN CALL HWDEUR(MHEP) GOTO 100 ENDIF ENDIF C Choose decay mode ISTHEP(MHEP)=ISTHEP(MHEP)+5 RN=HWR() BF=0. IM=LSTRT(IDM) DO 10 I=1,NMODES(IDM) BF=BF+BRFRAC(IM) IF (BF.GE.RN) GOTO 20 10 IM=LNEXT(IM) CALL HWWARN('HWDHAD',50,*20) 20 IF ((IDKPRD(1,IM).GE.1.AND.IDKPRD(1,IM).LE.13).OR. & (IDKPRD(3,IM).GE.1.AND.IDKPRD(3,IM).LE.13)) THEN C Partonic decay of a heavy-(b,c)-hadron, store details NQDK=NQDK+1 IF (NQDK.GT.NMXQDK) CALL HWWARN('HWDHAD',102,*999) LOCQ(NQDK)=MHEP IMQDK(NQDK)=IM CALL HWVEQU(4,VERTX,VTXQDK(1,NQDK)) GOTO 100 ELSE C Exclusive decay, add decay products to event record IF (NHEP+NPRODS(IM).GT.NMXHEP) & CALL HWWARN('HWDHAD',103,*999) JDAHEP(1,MHEP)=NHEP+1 DO 30 I=1,NPRODS(IM) NHEP=NHEP+1 IDHW(NHEP)=IDKPRD(I,IM) IDHEP(NHEP)=IDPDG(IDKPRD(I,IM)) ISTHEP(NHEP)=193 JMOHEP(1,NHEP)=MHEP JMOHEP(2,NHEP)=JMOHEP(2,MHEP) PHEP(5,NHEP)=RMASS(IDKPRD(I,IM)) 30 CALL HWVEQU(4,VERTX,VHEP(1,NHEP)) JDAHEP(2,MHEP)=NHEP ENDIF C Next choose momenta: IF (NPRODS(IM).EQ.1) THEN C 1-body decay: K0(BR) --> K0S,K0L CALL HWVEQU(4,PHEP(1,MHEP),PHEP(1,NHEP)) ELSEIF (NPRODS(IM).EQ.2) THEN C 2-body decay C---SPECIAL TREATMENT OF POLARIZED MESONS COSANG=TWO IF (ID.EQ.IDHW(JMOHEP(1,MHEP))) THEN MO=JMOHEP(1,MHEP) RSUM=0 DO 40 I=1,3 40 RSUM=RSUM+RHOHEP(I,MO) IF (RSUM.GT.ZERO) THEN RSUM=RSUM*HWR() IF (RSUM.LT.RHOHEP(1,MO)) THEN C---(1+COSANG)**2 COSANG=MAX(HWR(),HWR(),HWR())*TWO-ONE ELSEIF (RSUM.LT.RHOHEP(1,MO)+RHOHEP(2,MO)) THEN C---1-COSANG**2 COSANG=2*COS((ACOS(HWR()*TWO-ONE)+PIFAC)/THREE) ELSE C---(1-COSANG)**2 COSANG=MIN(HWR(),HWR(),HWR())*TWO-ONE ENDIF ENDIF ENDIF CALL HWDTWO(PHEP(1,MHEP),PHEP(1,NHEP-1), & PHEP(1,NHEP),CMMOM(IM),COSANG,.FALSE.) ELSEIF (NPRODS(IM).EQ.3) THEN C 3-body decay IF (NME(IM).EQ.100) THEN C Use free massless (V-A)*(V-A) Matrix Element CALL HWDTHR(PHEP(1,MHEP),PHEP(1,NHEP-1),PHEP(1,NHEP-2), & PHEP(1,NHEP),HWDWWT) ELSEIF (NME(IM).EQ.101) THEN C Use bound massless (V-A)*(V-A) Matrix Element WTMX=((PHEP(5,MHEP)-PHEP(5,NHEP)) & *(PHEP(5,MHEP)+PHEP(5,NHEP)) & +(PHEP(5,NHEP-1)-PHEP(5,NHEP-2)) & *(PHEP(5,NHEP-1)+PHEP(5,NHEP-2)))/TWO WTMX2=WTMX**2 IPDG=ABS(IDHEP(MHEP)) XS=ONE-MAX(RMASS(MOD(IPDG/1000,10)), & RMASS(MOD(IPDG/100,10)),RMASS(MOD(IPDG/10,10))) & /(RMASS(MOD(IPDG/1000,10))+RMASS(MOD(IPDG/100,10)) & +RMASS(MOD(IPDG/10,10))) 50 CALL HWDTHR(PHEP(1,MHEP),PHEP(1,NHEP-1),PHEP(1,NHEP-2), & PHEP(1,NHEP),HWDWWT) DOT1=HWULDO(PHEP(1,MHEP),PHEP(1,NHEP-1)) DOT2=HWULDO(PHEP(1,MHEP),PHEP(1,NHEP-2)) IF (DOT1*(WTMX-DOT1-XS*DOT2).LT.HWR()*WTMX2) GOTO 50 ELSE CALL HWDTHR(PHEP(1,MHEP),PHEP(1,NHEP-2),PHEP(1,NHEP-1), & PHEP(1,NHEP),HWDPWT) ENDIF ELSEIF (NPRODS(IM).EQ.4) THEN C 4-body decay CALL HWDFOR(PHEP(1,MHEP ),PHEP(1,NHEP-3),PHEP(1,NHEP-2), & PHEP(1,NHEP-1),PHEP(1,NHEP)) ELSEIF (NPRODS(IM).EQ.5) THEN C 5-body decay CALL HWDFIV(PHEP(1,MHEP ),PHEP(1,NHEP-4),PHEP(1,NHEP-3), & PHEP(1,NHEP-2),PHEP(1,NHEP-1),PHEP(1,NHEP)) ELSE CALL HWWARN('HWDHAD',104,*999) ENDIF ENDIF ENDIF 100 CONTINUE C---MAY HAVE OVERFLOWED /HEPEVT/ CALL HWWARN('HWDHAD',105,*999) 999 END CDECK ID>, HWDHGC. *CMZ :- -26/04/91 11.11.55 by Bryan Webber *-- Author : Mike Seymour C----------------------------------------------------------------------- SUBROUTINE HWDHGC(TAU,FNREAL,FNIMAG) C----------------------------------------------------------------------- C CALCULATE THE COMPLEX FUNCTION F OF HHG eq 2.18 C FOR USE IN H-->GAMMGAMM DECAYS C----------------------------------------------------------------------- INCLUDE 'HERWIG61.INC' DOUBLE PRECISION TAU,FNREAL,FNIMAG,FNLOG,FNSQR IF (TAU.GT.ONE) THEN FNREAL=(ASIN(1/SQRT(TAU)))**2 FNIMAG=0 ELSEIF (TAU.LT.ONE) THEN FNSQR=SQRT(1-TAU) FNLOG=LOG((1+FNSQR)/(1-FNSQR)) FNREAL=-0.25 * (FNLOG**2 - PIFAC**2) FNIMAG= 0.5 * PIFAC*FNLOG ELSE FNREAL=0.25*PIFAC**2 FNIMAG=0 ENDIF END