CDECK ID>, HWHV1J. *CMZ :- -18/05/99 14.37.45 by Mike Seymour *-- Author : Mike Seymour C----------------------------------------------------------------------- SUBROUTINE HWHV1J C----------------------------------------------------------------------- C V + 1 JET PRODUCTION, WHERE V=W (IHPRO.LT.5) OR Z (IHPRO.GE.5). C USES CROSS-SECTIONS OF EHLQ FOR ANNIHILATION AND COMPTON SCATTERING C IHPRO=0 FOR BOTH, 1 FOR ANNIHILATION, AND 2 FOR COMPTON. C----------------------------------------------------------------------- INCLUDE 'HERWIG61.INC' DOUBLE PRECISION HWR,HWRUNI,DISFAC(2,12,2),EMV2,DISMAX,S,T,U, & SHAT,THAT,UHAT,Z,HWUALF,PT,EMT,GFACTR,SIGANN,SIGCOM(2),CSFAC,ET, & EJ,YMIN,YMAX,VYMIN,VYMAX,EMAX,CV,CA,BR,EMV,GAMV INTEGER HWRINT,IDINIT(2,12,2),ICOFLO(4,2),I,J,K,L,M,ID1,ID2, $ IDV,IDI,IDM EXTERNAL HWRINT SAVE DISFAC,SHAT,THAT,EMV,EMV2,IDV,IDI C---IDINIT HOLDS THE INITIAL STATES FOR ANNIHILATION PROCESSES DATA IDINIT/1,8,2,7,3,10,4,9,5,12,6,11,1,10,2,9,3,8,4,7,5,12,6,11, $ 1,7,2,8,3,9,4,10,5,11,6,12,1,7,2,8,3,9,4,10,5,11,6,12/ C---ICOFLO HOLDS THE COLOR FLOW FOR EACH PROCESS C---DISFAC HOLDS THE DISTRIBUTION FUNCTION*CROSS-SECTION FOR EACH C POSSIBLE SUB-PROCESS. C INDEX1=INITIAL STATE PERMUTATION (1=AS IDINIT/QG;2=OPPOSITE/GQ), C 2=QUARK (FOR ANNIHILATION, >6 IMPLIES CABIBBO ROTATED PAIR), C 3=PROCESS (1=ANNIHILATION, 2=COMPTON) DATA ICOFLO,DISFAC/2,4,3,1,4,1,3,2,48*0./ IF (GENEV) THEN DISMAX=0 DO 110 I=1,2 DO 110 J=1,12 DO 110 K=1,2 110 DISMAX=MAX(DISFAC(K,J,I),DISMAX) 120 I=HWRINT(1,2) J=HWRINT(1,12) K=HWRINT(1,2) IF (HWR()*DISMAX.GT.DISFAC(K,J,I)) GOTO 120 IF (I.EQ.1) THEN C---ANNIHILATION IDN(1)=IDINIT(K,J,IDI) IDN(2)=IDINIT(3-K,J,IDI) IDN(4)=13 ELSE C---COMPTON SCATTERING IDN(1)=J IDN(2)=13 IF (IDV.EQ.200) THEN IDN(4)=J ELSE IF (J.EQ.5.OR.J.EQ.6.OR.J.GE.11.OR.HWR().GT.SCABI) THEN C---CHANGE QUARKS (1->2,2->1,3->4,4->3,...) IDN(4)=4*INT((J-1)/2)-J+3 ELSE C---CHANGE AND CABIBBO ROTATE QUARKS (1->4,2->3,3->2,...) IDN(4)=12*INT((J-1)/6)-J+5 ENDIF ENDIF IF ((SQRT(EMV2)+RMASS(IDN(4)))**2.GT.SHAT) GOTO 120 IF (K.EQ.2) THEN C---SWAP INITIAL STATES IDN(3)=IDN(1) IDN(1)=IDN(2) IDN(2)=IDN(3) ENDIF ENDIF IF (IDV.EQ.200) THEN IDN(3)=200 ELSE C---W+ OR W-? USE CHARGE CONSERVATION TO WORK OUT IDN(3)=NINT(198.5-.1667*FLOAT(ICHRG(IDN(1))+ICHRG(IDN(2)))) ENDIF M=K IF (I.EQ.2.AND.J.LE.6) M=3-K DO 130 L=1,4 130 ICO(L)=ICOFLO(L,M) IDCMF=15 COSTH=(SHAT+2*THAT-EMV2)/(SHAT-EMV2) C---TRICK HWETWO INTO USING THE OFF-SHELL V MASS RMASS(IDN(3))=SQRT(EMV2) CALL HWETWO RMASS(IDN(3))=EMV RHOHEP(1,NHEP-1)=0.5 RHOHEP(2,NHEP-1)=0.0 RHOHEP(3,NHEP-1)=0.5 ELSE EVWGT=0. IHPRO=MOD(IPROC,100)/10 IF (IHPRO.LT.5) THEN IDV=198 IDI=1 IDM=10 GAMV=GAMW ELSE IDV=200 IDI=2 IDM=6 GAMV=GAMZ IHPRO=IHPRO-5 ENDIF EMV=RMASS(IDV) EMV2=EMV*(EMV+GAMV*TAN(HWRUNI(0,-ONE-HALF,ONE+HALF))) IF (EMV2.LE.ZERO) RETURN CALL HWRPOW(ET,EJ) PT=0.5*ET EMT=SQRT(PT**2+EMV2) EMAX=0.5*(PHEP(5,3)+EMV2/PHEP(5,3)) IF (EMAX.LE.EMT) RETURN VYMAX=0.5*LOG((EMAX+SQRT(EMAX**2-EMT**2)) & /(EMAX-SQRT(EMAX**2-EMT**2))) VYMIN=-VYMAX IF (VYMAX.LE.VYMIN) RETURN Z=EXP(HWRUNI(0,VYMIN,VYMAX)) S= PHEP(5,3)**2 T=-PHEP(5,3)*EMT/Z+EMV2 U=-PHEP(5,3)*EMT*Z+EMV2 XXMIN=-U/(S+T-EMV2) IF (XXMIN.LT.ZERO.OR.XXMIN.GT.ONE) RETURN YMIN=MAX(LOG((XXMIN*PHEP(5,3)-EMT*Z)/PT),YJMIN) YMAX=MIN(LOG((PHEP(5,3)-EMT*Z)/PT),YJMAX) IF (YMAX.LE.YMIN) RETURN XX(1)=(Z*EMT+EXP(HWRUNI(2,YMIN,YMAX))*PT)/PHEP(5,3) IF (XX(1).LE.ZERO.OR.XX(1).GT.ONE) RETURN THAT =XX(1)*T+(1.-XX(1))*EMV2 XX(2)=-THAT / (XX(1)*S+U-EMV2) IF (XX(2).LT.ZERO.OR.XX(2).GT.ONE) RETURN UHAT =XX(2)*U+(1.-XX(2))*EMV2 SHAT =XX(1)*XX(2)*S EMSCA=EMT CALL HWSGEN(.FALSE.) GFACTR=GEV2NB*2.*PIFAC*ALPHEM*HWUALF(1,EMSCA)/(9.*SWEIN) SIGANN=GFACTR*((THAT-EMV2)**2+(UHAT-EMV2)**2) & /(SHAT**2*THAT*UHAT) SIGCOM(2)=.375*GFACTR*(SHAT**2+UHAT**2+2*EMV2*THAT) & /(-UHAT*SHAT**3) SIGCOM(1)=.375*GFACTR*(SHAT**2+THAT**2+2*EMV2*UHAT) & /(-THAT*SHAT**3) C---IF USER SPECIFIED A SUB-PROCESS, ZERO THE OTHER IF (IHPRO.EQ.1) THEN SIGCOM(1)=0. SIGCOM(2)=0. ENDIF IF (IHPRO.EQ.2) SIGANN=0. DO 210 I=1,IDM IF (IDV.EQ.200) THEN J=I IF(I.GT.6) J=I-6 DISFAC(1,I,1)=4*SWEIN*(VFCH(J,1)**2+AFCH(J,1)**2) ELSE IF (I.LE.4) THEN DISFAC(1,I,1)=1-SCABI ELSEIF (I.GE.7) THEN DISFAC(1,I,1)=SCABI ELSE DISFAC(1,I,1)=1. ENDIF ENDIF DISFAC(2,I,1)=DISFAC(1,I,1) * & SIGANN*DISF(IDINIT(1,I,IDI),2)*DISF(IDINIT(2,I,IDI),1) DISFAC(1,I,1)=DISFAC(1,I,1) * & SIGANN*DISF(IDINIT(1,I,IDI),1)*DISF(IDINIT(2,I,IDI),2) 210 CONTINUE DO 211 I=IDM+1,12 DISFAC(1,I,1)=0 DISFAC(2,I,1)=0 211 CONTINUE DO 220 I=1,12 IF (IDV.EQ.200) THEN J=I IF(I.GT.6) J=I-6 DISFAC(1,I,2)=4*SWEIN*(VFCH(J,1)**2+AFCH(J,1)**2) ELSE DISFAC(1,I,2)=1. ENDIF DISFAC(2,I,2)=DISFAC(1,I,2)*SIGCOM(2)*DISF(I,2)*DISF(13,1) DISFAC(1,I,2)=DISFAC(1,I,2)*SIGCOM(1)*DISF(I,1)*DISF(13,2) 220 CONTINUE DO 230 I=1,2 DO 230 J=1,12 DO 230 K=1,2 230 EVWGT=EVWGT+DISFAC(K,J,I) CSFAC=PT*EJ*(YMAX-YMIN)*(VYMAX-VYMIN)*THREE/PIFAC C---INCLUDE BRANCHING RATIO OF V CALL HWDBOZ(IDV,ID1,ID2,CV,CA,BR,0) EVWGT=EVWGT*CSFAC*BR ENDIF 999 END CDECK ID>, HWHWEX. *CMZ :- -26/04/91 14.55.45 by Federico Carminati *-- Author : Mike Seymour C----------------------------------------------------------------------- SUBROUTINE HWHWEX C----------------------------------------------------------------------- C TOP QUARK PRODUCTION VIA W EXCHANGE: MEAN EVWGT=TOP PROD C-S IN NB C C-S IS SUM OF: C UbarBbar, DBbar, DbarB, UB, CbarBbar, SBbar, SbarB, AND CB C UNLESS USER SPECIFIES OTHERWISE BY MOD(IPROC,100)=1-8 RESPECTIVELY C---DSDCOS HOLDS THE CROSS-SECTIONS FOR THE PROCESSES LISTED ABOVE C (1-8) ARE WITH B FROM BEAM 1, (9-16) ARE WITH B FROM BEAM 2. C----------------------------------------------------------------------- INCLUDE 'HERWIG61.INC' DOUBLE PRECISION HWR,HWRUNI,DSDCOS(16),EMT2,EMT,EMW2,EMW, & CMFMIN,TAUMIN,TAUMLN,S,T,U,ROOTS,DSMAX INTEGER HWRINT,IDHWEX(2,16),I EXTERNAL HWR,HWRUNI,HWRINT SAVE DSDCOS,DSMAX EQUIVALENCE (EMW,RMASS(198)),(EMT,RMASS(6)) C---IDHWEX HOLDS THE IDs OF THE INCOMING PARTICLES FOR EACH SUB-PROCESS DATA IDHWEX/11,8,11,1,5,7,5,2,11,10,11,3,5,9,5,4, & 8,11,1,11,7,5,2,5,10,11,3,11,9,5,4,5/ EMT2=EMT**2 EMW2=EMW**2 IF (GENEV) THEN 300 IHPRO=HWRINT(1,16) IF (HWR().GT.DSDCOS(IHPRO)/DSMAX) GOTO 300 DO 10 I=1,2 IDN(I)=IDHWEX(I,IHPRO) IF (IDN(I).EQ.5 .OR. IDN(I).EQ.11) THEN C---CHANGE B QUARK INTO T QUARK IDN(I+2)=IDN(I)+1 ELSEIF (HWR().GT.SCABI) THEN C---CHANGE QUARKS (1->2,2->1,3->4,4->3,7->8,8->7,...) IDN(I+2)=4*INT((IDN(I)-1)/2)-IDN(I)+3 ELSE C---CHANGE AND CABIBBO ROTATE QUARKS (1->4,2->3,3->2,4->1,7->10,...) IDN(I+2)=12*INT((IDN(I)-1)/6)-IDN(I)+5 ENDIF ICO(I)=I+2 ICO(I+2)=I 10 CONTINUE IDCMF=15 CALL HWETWO ELSE EVWGT=0. CMFMIN=EMT TAUMIN=(CMFMIN/PHEP(5,3))**2 TAUMLN=LOG(TAUMIN) ROOTS=PHEP(5,3)*SQRT(EXP(HWRUNI(0,ZERO,TAUMLN))) XXMIN=(ROOTS/PHEP(5,3))**2 XLMIN=LOG(XXMIN) COSTH=HWRUNI(0,-ONE, ONE) S=ROOTS**2 T=-0.5*S*(1-COSTH) U=-0.5*S*(1+COSTH) EMSCA=SQRT(2*S*T*U/(S*S+T*T+U*U)) DSDCOS(1)=GEV2NB*PIFAC*.125*(ALPHEM/SWEIN)**2 & *(S-EMT2)**2 / S / (EMW2 + 0.5*(S-EMT2)*(1-COSTH))**2 DSDCOS(2)=DSDCOS(1) / 4 & * (1 + EMT2/S + 2*COSTH + (1-EMT2/S)*COSTH**2) DSDCOS(3)=DSDCOS(2) DSDCOS(4)=DSDCOS(1) C---IF USER SPECIFIED SUB-PROCESS THEN ZERO ALL THE OTHERS IHPRO=MOD(IPROC,100) IF (IHPRO.GT.8) THEN CALL HWWARN('HWHWEX',1,*999) IHPRO=0 ENDIF DO 100 I=1,8 IF (I.LE.4) DSDCOS(I+4)=DSDCOS(I) IF (IHPRO.NE.0 .AND. IHPRO.NE.I) DSDCOS(I)=0 DSDCOS(I+8)=DSDCOS(I) 100 CONTINUE CALL HWSGEN(.TRUE.) DSMAX=0 DO 200 I=1,16 DSDCOS(I)=DSDCOS(I)*DISF(IDHWEX(1,I),1)*DISF(IDHWEX(2,I),2) EVWGT=EVWGT + 2*TAUMLN*XLMIN*DSDCOS(I) IF (DSDCOS(I).GT.DSMAX) DSMAX=DSDCOS(I) 200 CONTINUE ENDIF 999 END CDECK ID>, HWHWPR. *CMZ :- -18/05/99 14.22.13 by Mike Seymour *-- Author : Bryan Webber C----------------------------------------------------------------------- SUBROUTINE HWHWPR C----------------------------------------------------------------------- C W+/- PRODUCTION AND DECAY VIA DRELL-YAN PROCESS C MEAN EVWGT IS SIG(W+/-)*(BRANCHING FRACTION) IN NB C----------------------------------------------------------------------- INCLUDE 'HERWIG61.INC' DOUBLE PRECISION HWR,HWRUNI,HWUPCM,PRAN,PROB,COEF,CSFAC,EMW, & FTQK,PTOP,ETOP,EBOT,PMAX,FHAD,FTOT,BRAF,FLEP,TMIN,HWUAEM INTEGER HWRINT,ICH,IC,IL,ID,IDEC,JDEC,IWP(2,16) LOGICAL HWRLOG EXTERNAL HWR,HWRUNI,HWUPCM,HWRINT,HWRLOG SAVE CSFAC,IDEC,FLEP,FTQK,ETOP,PTOP,EBOT,PMAX,PROB DATA IWP/2,7,1,8,7,2,8,1,4,9,3,10,9,4,10,3, & 2,9,3,8,9,2,8,3,4,7,1,10,7,4,10,1/ IF (GENEV) THEN C---GENERATE EVENT (X'S AND STRUCTURE FUNCTIONS ALREADY FOUND) PRAN=PROB*HWR() C---LOOP OVER PARTON FLAVOURS PROB=0. COEF=1.-SCABI DO 10 IC=1,16 IF (IC.EQ.9) COEF=SCABI PROB=PROB+DISF(IWP(1,IC),1)*DISF(IWP(2,IC),2)*COEF IF (PROB.GE.PRAN) GOTO 20 10 CONTINUE C---STORE INCOMING PARTONS 20 IDN(1)=IWP(1,IC) IDN(2)=IWP(2,IC) ICO(1)=2 ICO(2)=1 C---ICH=1/2 FOR W+/- ICH=2-MOD(IC,2) IF ((IDEC.GT.49.AND.IDEC.LT.54).OR. & (IDEC.EQ.99.AND.HWRLOG(FLEP))) THEN C---LEPTONIC DECAY IL=IDEC-50 IF (IL.EQ.0.OR.IL.GT.3) IL=HWRINT(1,3) IDN(3)=2*IL+121-ICH IDN(4)=2*IL+124+ICH C---W DECAY ANGLE (1+COSTH)**2 COSTH=2.*HWR()**0.3333-1. ELSEIF (IDEC.EQ.5.OR.IDEC.EQ.6.OR. & ((IDEC.EQ.0.OR.IDEC.EQ.99).AND.HWRLOG(FTQK))) THEN C---W -> TOP + BOTTOM DECAY IDN(3)=7-ICH IDN(4)=10+ICH 21 COSTH=HWRUNI(1,-ONE, ONE) IF ((ETOP+(PTOP*COSTH))*(EBOT+(PTOP*COSTH)).LT. & PMAX*HWR()) GOTO 21 ELSE C---OTHER HADRONIC DECAY 25 PROB=0. PRAN=2.*HWR() COEF=1.-SCABI DO 30 ID=ICH,16,4 IF (ID.GT.8) COEF=SCABI PROB=PROB+COEF IF (PROB.GE.PRAN) THEN IDN(3)=IWP(1,ID) IDN(4)=IWP(2,ID) GOTO 40 ENDIF 30 CONTINUE 40 CONTINUE IF (IDEC.GT.0.AND.IDEC.LT.5) THEN JDEC=IDEC+6 IF (IDN(3).NE.IDEC.AND.IDN(4).NE.IDEC & .AND.IDN(3).NE.JDEC.AND.IDN(4).NE.JDEC) GOTO 25 ENDIF COSTH=2.*HWR()**0.3333-1. ENDIF IDCMF=197+ICH IF (IDN(1).GT.6) COSTH=-COSTH ICO(3)=4 ICO(4)=3 CALL HWETWO ELSE IDEC=MOD(IPROC,100) IF (IDEC.EQ.5.OR.IDEC.EQ.6) THEN TMIN=ATAN((RMASS(6)**2-RMASS(199)**2)/(GAMW*RMASS(199))) ELSE TMIN=-ATAN(RMASS(199)/GAMW) ENDIF EVWGT=0. EMW=GAMW*TAN(HWRUNI(0,TMIN,PIFAC/2.))+RMASS(199) IF (EMW.LE.ZERO) RETURN EMW=SQRT(EMW*RMASS(199)) IF (EMW.LE.QSPAC.OR.EMW.GE.PHEP(5,3)) RETURN EMSCA=EMW IF (EMLST.NE.EMW) THEN EMLST=EMW XXMIN=(EMW/PHEP(5,3))**2 XLMIN=LOG(XXMIN) CSFAC=-GEV2NB*PIFAC**2*HWUAEM(EMSCA**2) & /(3.*SWEIN*RMASS(199)**2)*XLMIN C---COMPUTE TOP AND LEPTONIC FRACTIONS FTQK=0. IF (NFLAV.GT.5) THEN PTOP=HWUPCM(EMW,RMASS(5),RMASS(6)) IF (PTOP.GT.ZERO) THEN ETOP=SQRT(PTOP**2+RMASS(6)**2) EBOT=EMW-ETOP FTQK=2.*PTOP*(3.*ETOP*EBOT+PTOP**2)/EMW**3 PMAX=(ETOP+PTOP)*(EBOT+PTOP) ENDIF ENDIF FHAD=FTQK+2. FTOT=FTQK+3. C---MULTIPLY WEIGHT BY BRANCHING FRACTION IF (IDEC.EQ.0) THEN BRAF=FHAD ELSEIF (IDEC.LT.5.OR.IDEC.EQ.50) THEN BRAF=1. ELSEIF (IDEC.LT.7) THEN BRAF=FTQK ELSEIF (IDEC.EQ.99) THEN BRAF=FTOT ELSE BRAF=1/THREE ENDIF CSFAC=CSFAC*BRAF/FTOT*(0.5-TMIN/PIFAC) FTQK=FTQK/FHAD FLEP=1./FTOT ENDIF CALL HWSGEN(.TRUE.) C---LOOP OVER PARTON FLAVOURS PROB=0. COEF=1.-SCABI DO 100 IC=1,16 IF (IC.EQ.9) COEF=SCABI PROB=PROB+DISF(IWP(1,IC),1)*DISF(IWP(2,IC),2)*COEF 100 CONTINUE EVWGT=PROB*CSFAC ENDIF 999 END CDECK ID>, HWIGIN. *CMZ :- -01/04/99 19.44.55 by Mike Seymour *-- Author : Bryan Webber C---------------------------------------------------------------------- SUBROUTINE HWIGIN C----------------------------------------------------------------------- C SETS INPUT PARAMETERS C---------------------------------------------------------------------- INCLUDE 'HERWIG61.INC' DOUBLE PRECISION FAC,ANGLE INTEGER I,J,N,L CHARACTER*28 TITLE DATA TITLE/'HERWIG 6.100 December 1999'/ WRITE (6,10) TITLE 10 FORMAT(//10X,A28//, & 10X,'Please reference: G. Marchesini, B.R. Webber,',/, & 10X,'G.Abbiendi, I.G.Knowles, M.H.Seymour & L.Stanco',/, & 10X,'Computer Physics Communications 67 (1992) 465') C---PRINT OPTIONS: C IPRINT=0 NO PRINTOUT C 1 PRINT SELECTED INPUT PARAMETERS C 2 1 + TABLE OF PARTICLE CODES AND PROPERTIES C 3 2 + TABLES OF SUDAKOV FORM FACTORS IPRINT=1 C Format for track numbers in event listing C PRNDEC=.TRUE. use decimal C .FALSE. use hexadecimal PRNDEC=(NMXHEP.LE.9999) C Number of significant figures to print out in event listing C NPRFMT (< 2) compact 80 character stout and A4-long tex output, C (= 2) 2 decimal places in stout, (> 2) - 5 decimal places in stout NPRFMT=1 C Print out vertex information PRVTX=.TRUE. C Print out particle properties/event record to stout, tex or web PRNDEF=.TRUE. PRNTEX=.FALSE. PRNWEB=.FALSE. C---MAX NO OF EVENTS TO PRINT MAXPR=0 C---UNIT FOR READING SUDAKOV FORM FACTORS (IF ZERO THEN COMPUTE THEM) LRSUD=0 C---UNIT FOR WRITING SUDAKOV FORM FACTORS (IF ZERO THEN NOT WRITTEN) LWSUD=77 C---UNIT FOR WRITING EVENT DATA IN HWANAL (IF ZERO THEN NOT WRITTEN) LWEVT=0 C---SEEDS FOR RANDOM NUMBER GENERATOR (CALLED HWR) NRN(1)= 17673 NRN(2)= 63565 C---AZIMUTHAL CORRELATIONS? C THESE INCLUDE SOFT GLUON (INSIDE CONE) AZSOFT=.TRUE. C AND NEAREST-NEIGHBOUR SPIN CORRELATIONS AZSPIN=.TRUE. C---MATRIX-ELEMENT MATCHING FOR E+E-, DIS, DRELL-YAN AND TOP DECAY C---HARD EMISSION HARDME=.TRUE. C---SOFT EMISSION SOFTME=.TRUE. C---GLUON ENERGY CUT FOR TOP DECAY CASE GCUTME=2 C Electromagnetic fine structure constant: Thomson limit ALPHEM=.0072993 C---QCD LAMBDA: CORRESPONDS TO 5-FLAVOUR LAMBDA-MS-BAR AT LARGE X ONLY QCDLAM=0.18 C---NUMBER OF COLOURS NCOLO=3 C---NUMBER OF FLAVOURS NFLAV=6 C---QUARK, GLUON AND PHOTON VIRTUAL MASS CUTOFFS IN C PARTON SHOWER (ADDED TO MASSES GIVEN BELOW) VQCUT=0.48 VGCUT=0.10 VPCUT=0.40 ALPFAC=1 C---D,U,S,C,B,T QUARK AND GLUON MASSES (IN THAT ORDER) RMASS(1)=0.32 RMASS(2)=0.32 RMASS(3)=0.5 RMASS(4)=1.55 RMASS(5)=4.95 RMASS(6)=175. RMASS(13)=0.75 C---W+/- AND Z0 MASSES RMASS(198)=80.4 RMASS(199)=80.4 RMASS(200)=91.2 C---HIGGS BOSON MASS RMASS(201)=150. C---WIDTHS OF W, Z, HIGGS GAMW=2.0 GAMZ=2.5 GAMH=0.02 C Include additional neutral, massive vector boson (Z') ZPRIME=.FALSE. C Z' mass and width RMASS(202)=500. GAMZP=5. C Lepton (EPOLN) and anti-lepton (PPOLN) beam polarisations used in: C e+e- --> ffbar/qqbar g; and l/lbar N DIS. C Cpts. 1,2 Transverse polarisation; cpt. 3 longitudinal polarisation. C Note require POLN(1)**2+POLN(2)**2+POLN(3)**2 < 1. DO 20 I=1,3 EPOLN(I)=0. 20 PPOLN(I)=0. C----------------------------------------------------------------------- C Specify couplings of weak vector bosons to fermions: C C electric current: QFCH(I)*e*G_mu (electric charge, e>0) C weak neutral current: [VFCH(I,J).1+AFCH(I,J).G_5]*e*G_mu C weak charged current: SQRT(VCKM(K,L)/2.)*g*(1+G_5)*G_mu C C I= 1- 6: d,u,s,c,b,t (quarks) C =11-16: e,nu_e,mu,nu_mu,tau,nu_tau (leptons) (`I=IDHW-110') C J=1 for minimal SM: C =2 for Z' couplings (ZPRIME=.TRUE.) C K=1,2,3 for u,c,t; L=1,2,3 for d,s,b C----------------------------------------------------------------------- C Minimal standard model neutral vector boson couplings C VFCH(I,1)=(T3/2-Q*S^2_W)/(C_W*S_W); AFCH(I,1)=T3/(2*C_W*S_W) C sin**2 Weinberg angle (PDG '94) SWEIN=.2319 FAC=1./SQRT(SWEIN*(1.-SWEIN)) DO 30 I=1,3 C Down-type quarks J=2*I-1 QFCH(J)=-1./3. VFCH(J,1)=(-0.25+SWEIN/3.)*FAC AFCH(J,1)= -0.25*FAC C Up-type quarks J=2*I QFCH(J)=+2./3. VFCH(J,1)=(+0.25-2.*SWEIN/3.)*FAC AFCH(J,1)= +0.25*FAC C Charged leptons J=2*I+9 QFCH(J)=-1. VFCH(J,1)=(-0.25+SWEIN)*FAC AFCH(J,1)= -0.25*FAC C Neutrinos J=2*I+10 QFCH(J)=0. VFCH(J,1)=+0.25*FAC AFCH(J,1)=+0.25*FAC 30 CONTINUE C Additional Z' couplings (To be set by the user) IF (.NOT.ZPRIME) THEN DO 40 I=1,6 AFCH(I,2)=0. AFCH(10+I,2)=0. VFCH(I,2)=0. VFCH(10+I,2)=0. 40 CONTINUE ENDIF C Cabibbo-Kobayashi-Maskawa matrix elements squared (PDG '92): C sin**2 of Cabibbo angle SCABI=.0488 C u ---> d,s,b VCKM(1,1)=1.-SCABI VCKM(1,2)=SCABI VCKM(1,3)=0.0 C c ---> d,s,b VCKM(2,1)=SCABI VCKM(2,2)=1.-SCABI-.002 VCKM(2,3)=0.002 C t ---> d,b,s VCKM(3,1)=0.0 VCKM(3,2)=0.002 VCKM(3,3)=0.998 C---GAUGE BOSON DECAYS DO 50 I=1,12 BRHIG(I)=1.D0/12 ENHANC(I)=1.D0 50 IF (I.LE.MODMAX) MODBOS(I)=0 C C THE iTH GAUGE BOSON DECAY PER EVENT IS CONTROLLED BY MODBOS AS FOLLOWS C MODBOS(i) W DECAY Z DECAY C 0 all all C 1 qqbar qqbar C 2 enu e+e- C 3 munu mu+mu- C 4 taunu tau+tau- C 5 enu & munu ee & mumu C 6 all nunu C 7 all bbbar C >7 all all C BOSON PAIRS (eg FROM HIGGS DECAY) ARE CHOSEN FROM MODBOS(i),MODBOS(i+1 C C---CONTROL OF LARGE EMH BEHAVIOUR (SEE HWHIGM FOR DETAILS) IOPHIG=3 GAMMAX=10. C Specicify approximation used in HWHIGA IAPHIG=1 C---MASSES OF HYPOTHETICAL NEW QUARKS GO C INTO 209-214 (ANTIQUARKS IN 215-220) C ID = 209,210 ARE B',T' WITH DECAYS T'->B'->C C 211,212 ARE B',T' WITH DECAYS T'->B'->T C 215-218 ARE THEIR ANTIQUARKS RMASS(209)=200. RMASS(215)=200. C---MAXIMUM CLUSTER MASS PARAMETERS C N.B. LIMIT FOR Q1-Q2BAR CLUSTER MASS C IS (CLMAX**CLPOW + (QM1+QM2)**CLPOW)**(1/CLPOW) CLMAX=3.35 CLPOW=2.0 C For PSPLT(I), CLDIR(I) & CLSMR(I): I=1 light u,d,s,c cluster C =2 heavy b cluster C---MASS SPECTRUM OF PRODUCTS IN CLUSTER C SPLITTING ABOVE CLMAX - FLAT IN M**PSPLT(*) PSPLT(1)=1.0 PSPLT(2)=PSPLT(1) C---KINEMATIC TREATMENT OF CLUSTER DECAY C 0=ISOTROPIC, 1=REMEMBER DIRECTION OF PERTURBATIVELY PRODUCED QUARKS CLDIR(1)=1 CLDIR(2)=CLDIR(1) C IF CLDIR(*)=1, DO GAUSSIAN SMEARING OF DIRECTION: C ACTUALLY EXPONENTIAL IN 1-COS(THETA) WITH MEAN CLSMR(*) CLSMR(1)=0.0 CLSMR(2)=CLSMR(1) C---OPTION FOR TREATMENT OF REMNANT CLUSTERS: C 0=BOTH CHILDREN ARE SOFT, (EQUIVALENT TO PREVIOUS VERSIONS) C 1=REMNANT CHILD IS SOFT, BUT PERTURBATIVE CHILD IS NORMAL IOPREM=1 C---TREATMENT OF LOWER LIMIT FOR SPACELIKE EVOLUTION C 0=EVOLUTION STOPS AT QSPAC, BUT STRUCT FUNS CAN GET CALLED AT C SMALLER SCALES IN FORCED EMISSION (EQUIVALENT TO V5.7 AND EARLIER) C 1=EVOLUTION STOPS AT QSPAC, STRUCTURE FUNCTIONS FREEZE AT QSPAC C 2=EVOLUTION CONTINUES TO INFRARED CUT, BUT S.F.S FREEZE AT QSPAC ISPAC=0 C---LOWER LIMIT FOR SPACELIKE EVOLUTION QSPAC=2.5 C---SWITCH OFF SPACE-LIKE SHOWERS NOSPAC=.FALSE. C---INTRINSIC PT OF SPACELIKE PARTONS (RMS) PTRMS=0.0 C---MASS PARAMETER IN REMNANT FRAGMENTATION BTCLM=1.0 C---STRUCTURE FUNCTION SET: C SET MODPDF(I)=MODE AND AUTPDF='AUTHOR GROUP' TO USE CERN LIBRARY C PDFLIB PACKAGE FOR STRUCTURE FUNCTIONS IN BEAM I MODPDF(1)=-1 MODPDF(2)=-1 AUTPDF(1)='MRS' AUTPDF(2)='MRS' C OR SET MODPDF(I)=-1 TO USE BUILT-IN STRUCTURE FUNCTION SET: C 1,2 FOR DUKE+OWENS SETS 1,2 (SOFT/HARD GLUE) C 3,4 FOR EICHTEN+AL SETS 1,2 (NUCLEONS ONLY) C 5 FOR OWENS SET 1.1 (SOFT GLUE ONLY) NSTRU=5 C PARAMETER FOR B CLUSTER DECAY TO 1 HADRON. IF MCL IS CLUSTER MASS C AND MTH IS THRESHOLD FOR 2-HADRON DECAY, THEN PROBABILITY IS C 1 IF MCL(1+B1LIM)*MTH, WITH LINEAR INTERPOLATION, B1LIM=0.0 C---B DECAY PACKAGE ('HERW'=>HERWIG, 'EURO'=>EURODEC, 'CLEO'=>CLEO) BDECAY='HERW' C---HARD SUBPROCESS SCALE TO BE USED IN BOSON-GLUON FUSION C IF (BGSHAT) THEN SCALE=SHAT C ELSE SCALE=2.*SHAT*THAT*UHAT/(SHAT**2+THAT**2+UHAT**2) BGSHAT=.FALSE. C---RECONSTRUCT DIS EVENTS IN BREIT FRAME BREIT=.TRUE. C---TREAT ALL EVENTS IN THEIR CMF (ELSE USE LAB FRAME) USECMF=.TRUE. C---PROBABILITY OF UNDERLYING SOFT EVENT: PRSOF=1. C---SOFT UNDERLYING OR MIN BIAS EVENT PARAMETERS C DEFAULT VALUES ARE FROM UA5 COLLAB, NPB291(1987)445 C NCH_PPBAR(SQRT(S)) = PMBN1*S**PMBN2+PMBN3 PMBN1= 9.11 PMBN2= 0.115 PMBN3=-9.50 C 1/K (IN NEG BINOMIAL) = PMBK1*LN(S)+PMBK2 PMBK1= 0.029 PMBK2=-0.104 C SOFT CLUSTER MASS SPECTRUM (M-M1-M2-PMBM1)*EXP(-PMBM2*M) PMBM1= 0.4 PMBM2= 2.0 C SOFT CLUSTER PT SPECTRUM PT*EXP(-B*SQRT(PT**2+M**2)) C B=PMBP1 FOR D,U, PMBP2 FOR S,C, PMBP3 FOR DIQUARKS PMBP1= 5.2 PMBP2= 3.0 PMBP3= 5.2 C---MULTIPLICITY ENHANCEMENT FOR UNDERLYING SOFT EVENT: C NCH = NCH_PPBAR(ENSOF*SQRT(S)) ENSOF=1. C PARAMETERS FOR MUELLER TANG FORMUA: IPROC=2400 C---THE VALUE TO USE FOR FIXED ALPHA_S IN DENOMINATOR ASFIXD=0.25 C---OMEGA0=12*LOG(2)*ALPHA_S/PI, BUT NOT NECESSARILY THE SAME ALPHA_S OMEGA0=0.3 C---MIN AND MAX JET RAPIDITIES IN QCD 2->2, C HEAVY FLAVOUR, SUSY AND DIRECT PHOTON PROCESSES YJMAX=8. YJMIN=-YJMAX C---MIN AND MAX PARTON TRANSVERSE MOMENTUM C IN ELEMENTARY 2 -> 2 SUBPROCESSES PTMIN=1D1 PTMAX=1D8 C---UPPER LIMIT ON HARD PROCESS SCALE QLIM=1D8 C---MAX PARTON THRUST IN 2->3 HARD PROCESSES THMAX=0.9 C Set parameters for 2->4 hard process C Choose inter-jet metric (else JADE) and minimum y-cut DURHAM=.TRUE. Y4JT=0.01 C---TREATMENT OF COLOUR INTERFERENCE IN E+E- -> 4 JETS: C qqbar-gg case: C IOP4JT(1)=0 neglect, =1 extreme 2341; =2 extreme 3421 C qqbar-qqbar (identical quark flavour) case: C IOP4JT(2)=0 neglect, =1 extreme 4123; =2 extreme 2143 IOP4JT(1)=0 IOP4JT(2)=0 C---MIN AND MAX DILEPTON INVARIANT MASS IN DRELL-YAN PROCESS EMMIN=0D0 EMMAX=1D8 C---MIN AND MAX ABS(Q**2) IN DEEP INELASTIC LEPTON SCATTERING Q2MIN=0D0 Q2MAX=1D10 C---MIN AND MAX ABS(Q**2) IN WEISZACKER-WILLIAMS APPROXIMATION Q2WWMN=0. Q2WWMX=4. C---MIN AND MAX ENERGY FRACTION IN WEISZACKER-WILLIAMS APPROXIMATION YWWMIN=0. YWWMAX=1. C---MINIMUM HADRONIC MASS FOR PHOTON-INDUCED PROCESSES (INCLUDING DIS) WHMIN=0. C---IF PHOMAS IS NON-ZERO, PARTON DISTRIBUTION FUNCTIONS FOR OFF-SHELL C PHOTONS IS DAMPED, WITH MASS PARAMETER = PHOMAS PHOMAS=0. C---MIN AND MAX FLAVOURS GENERATED BY IPROC=9100,9110,9130 IFLMIN=1 IFLMAX=5 C---MAX Z IN J/PSI PHOTO- AND ELECTRO- PRODUCTION ZJMAX=0.9 C---MIN AND MAX BJORKEN-Y YBMIN=0. YBMAX=1. C---MAX COS(THETA) FOR W'S IN E+E- -> W+W- CTMAX=0.9999 C Minimum virtuality^2 of partons to use in calculating distances VMIN2=0.1 C Exageration factor for lifetimes of weakly decaying heavy particles EXAG=1. C Include colour rearrangement in cluster formation CLRECO=.FALSE. C Probability for colour rearrangement to occur PRECO=1./9. C Minimum lifetime for particle to be considered stable PLTCUT=1.D-8 C Incude neutral B-meson mixing MIXING=.TRUE. C Set B_s and B_d mixing parameters: X=Delta m/Gamma XMIX(1)=10.0 XMIX(2)=0.70 C Y=Delta Gamma/2*Gamma YMIX(1)=0.2 YMIX(2)=0.0 C Include a cut on particle decay lengths MAXDKL=.FALSE. C Set option for decay length cut (see HWDXLM) IOPDKL=1 C Smear the primary interaction vertex: see HWRPIP for details PIPSMR=.TRUE. DO 60 I=0,NMXRES C Veto cluster decays into particle type I VTOCDK(I)=.FALSE. C Veto unstable particle decays into modes involving particle type I 60 VTORDK(I)=.FALSE. C Veto f_0(980) and a_0(980) production in cluster decays VTOCDK(290)=.TRUE. VTOCDK(291)=.TRUE. VTOCDK(292)=.TRUE. VTOCDK(293)=.TRUE. C---MINIMUM AND MAXIMUM S-HAT/S RANGE FOR PHOTON ISR TMNISR=1D-4 ZMXISR=1-1D-6 C---COLISR IS .TRUE. TO MAKE ISR PHOTONS COLLINEAR WITH BEAMS COLISR=.FALSE. C A Priori weights for mesons w.r.t. pionic n=1, 0-(+) states: C old VECWT=REPWT(0,1,0) & TENWT=REPWT(0,2,0) DO 70 N=0,4 DO 70 J=0,4 DO 70 L=0,3 70 REPWT(L,J,N)=1. C and singlet (Lambda-like) and decuplet barons SNGWT=1. DECWT=1. C---A PRIORI WEIGHTS FOR D,U,S,C,B,T QUARKS AND DIQUARKS (IN THAT ORDER) PWT(1)=1. PWT(2)=1. PWT(3)=1. PWT(4)=1. PWT(5)=1. PWT(6)=1. PWT(7)=1. C Octet-Singlet isoscalar mixing angles in degrees C (use ANGLE for ideal mixing, recommended for F0MIX & OMHMIX) ANGLE=ATAN(ONE/SQRT(TWO))*180./ACOS(-ONE) C eta - eta' ETAMIX=-23. C phi - omega PHIMIX=+36. C h_1(1380) - h_1(1170) H1MIX=ANGLE C MISSING - f_0(1370) F0MIX=ANGLE C f_1(1420) - f_1(1285) F1MIX=ANGLE C f'_2 - f_2 F2MIX=+26. C MISSING - omega(1600) OMHMIX=ANGLE C eta_2(1645) - eta_2(1870) ET2MIX=ANGLE C phi_3 - omega_3 PH3MIX=+28. C---PARAMETERS FOR NON-PERTURBATIVE SPLITTING OF GLUONS INTO C DIQUARK-ANTIDIQUARK PAIRS: C SCALE AT WHICH GLUONS CAN BE SPLIT INTO DIQUARKS C (0.0 FOR NO SPLITTING) QDIQK=0.0 C PROBABILITY (PER UNIT LOG SCALE) OF DIQUARK SPLITTING PDIQK=5.0 C---PARAMETERS FOR IMPORTANCE SAMPLING C ASSUME QCD 2->2 DSIG/DET FALLS LIKE ET**(-PTPOW) C WHERE ET=SQRT(MQ**2+PT**2) FOR HEAVY FLAVOURS PTPOW=4. C DEFAULT PTPOW=2 FOR SUSY PROCESSES IF (MOD(IPROC/100,100).EQ.30) PTPOW=2. C ASSUME DRELL-YAN DSIG/DEM FALLS LIKE EM**(-EMPOW) EMPOW=4. C ASSUME DEEP INELASTIC DSIG/DQ**2 FALLS LIKE (Q**2)**(-Q2POW) Q2POW=2.5 C---GENERATE UNWEIGHTED EVENTS (EVWGT=AVWGT)? NOWGT=.TRUE. C---DEFAULT MEAN EVENT WEIGHT AVWGT=1. C---ASSUMED MAXIMUM WEIGHT (ZERO TO RECOMPUTE) WGTMAX=0. C---MINIMUM ACCEPTABLE EVENT GENERATION EFFICIENCY EFFMIN=1E-3 C---MAX NO OF (CODE.GE.100) ERRORS MAXER=10 C---TIME (SEC) NEEDED TO TERMINATE GRACEFULLY TLOUT=5. C---CURRENT NO OF EVENTS NEVHEP=0 C---CURRENT NO OF ENTRIES IN /HEPEVT/ NHEP=0 C---ISTAT IS STATUS OF EVENT (I.E. STAGE IN PROCESSING) ISTAT=0 C---IERROR IS ERROR CODE IERROR=0 C---MORE TECHNICAL PARAMETERS - SHOULDN'T NEED ADJUSTMENT C---PI PIFAC=ACOS(-1.D0) C Speed of light (mm/s) CSPEED=2.99792D11 C Cross-section conversion factor (hbar.c/e)**2 GEV2NB=389380 C---NUMBER OF SHOTS FOR INITIAL MAX WEIGHT SEARCH IBSH=10000 C---RANDOM NO. SEEDS FOR INITIAL MAX WEIGHT SEARCH IBRN(1)=1246579 IBRN(2)=8447766 C---NUMBER OF ENTRIES IN LOOKUP TABLES OF SUDAKOV FORM FACTORS NQEV=1024 C---MAXIMUM BIN SIZE IN Z FOR SPACELIKE BRANCHING ZBINM=0.05 C---MAXIMUM NUMBER OF Z BINS FOR SPACELIKE BRANCHING NZBIN=100 C---MAXIMUM NUMBER OF BRANCH REJECTIONS (TO AVOID INFINITE LOOPS) NBTRY=200 C---MAXIMUM NUMBER OF TRIES TO GENERATE CLUSTER DECAY NCTRY=200 C---MAXIMUM NUMBER OF TRIES TO GENERATE MASS REQUESTED NETRY=200 C---MAXIMUM NUMBER OF TRIES TO GENERATE SOFT SUBPROCESS NSTRY=200 C---PRECISION FOR GAUSSIAN INTEGRATION ACCUR=1.D-6 C---ORDER OF INTERPOLATION IN SUDAKOV TABLES INTER=3 C---ORDER TO USE FOR ALPHAS IN SUDAKOV TABLES SUDORD=1 C--CONSERVATION OF RPARITY RPARTY = .TRUE. C--CHECK WHETHER SUSY DATA INPUTTED SUSYIN=.FALSE. 999 END CDECK ID>, HWIODK. *CMZ :- -27/07/99 13.33.03 by Mike Seymour *-- Author : Ian Knowles C----------------------------------------------------------------------- SUBROUTINE HWIODK(IUNIT,IOPT,IME) C----------------------------------------------------------------------- C If IUNIT > 0 writes out present HERWIG decay tables to unit IUNIT C < 0 reads in decay tables from unit IUNIT C The format used during the read/write is specified by IOPT C =1 PDG; =2 HERWIG numeric; =3 HERWIG character name. C When reading in if IME =1 matrix element codes >= 100 are accepted C 0 are set zero. C----------------------------------------------------------------------- INCLUDE 'HERWIG61.INC' INTEGER IUNIT,IOPT,IME,JUNIT,I,J,K,L,IDKY,ITMP(5),IDUM CHARACTER*8 CDK(NMXDKS),CDKPRD(5,NMXDKS),CDUM JUNIT=ABS(IUNIT) OPEN(UNIT=JUNIT,FORM='FORMATTED',STATUS='UNKNOWN') IF (IUNIT.GT.0) THEN C Write out the decay table WRITE(JUNIT,100) NDKYS IF (IOPT.EQ.1) THEN DO 20 I=1,NRES IF (NMODES(I).EQ.0) GOTO 20 K=LSTRT(I) DO 10 J=1,NMODES(I) WRITE(JUNIT,110) IDPDG(I),BRFRAC(K),NME(K), & (IDPDG(IDKPRD(L,K)),L=1,5) 10 K=LNEXT(K) 20 CONTINUE ELSEIF (IOPT.EQ.2) THEN DO 40 I=1,NRES IF (NMODES(I).EQ.0) GOTO 40 K=LSTRT(I) DO 30 J=1,NMODES(I) WRITE(JUNIT,120) I,BRFRAC(K),NME(K),(IDKPRD(L,K),L=1,5) 30 K=LNEXT(K) 40 CONTINUE ELSEIF (IOPT.EQ.3) THEN DO 60 I=1,NRES IF (NMODES(I).EQ.0) GOTO 60 K=LSTRT(I) DO 50 J=1,NMODES(I) WRITE(JUNIT,130) RNAME(I),BRFRAC(K),NME(K), & (RNAME(IDKPRD(L,K)),L=1,5) 50 K=LNEXT(K) 60 CONTINUE ENDIF ELSEIF (IUNIT.LT.0) THEN C Read in the decay table and convert to HERWIG numeric format READ(JUNIT,100) NDKYS IF (NDKYS.GT.NMXDKS) CALL HWWARN('HWIODK',100,*999) IF (IOPT.EQ.1) THEN DO 70 I=1,NDKYS READ(JUNIT,110) IDKY,BRFRAC(I),NME(I),ITMP IF (IME.EQ.0.AND.NME(I).GE.100) NME(I)=0 CALL HWUIDT(1,IDKY,IDK(I),CDUM) DO 70 J=1,5 70 CALL HWUIDT(1,ITMP(J),IDKPRD(J,I),CDUM) ELSEIF (IOPT.EQ.2) THEN DO 80 I=1,NDKYS READ(JUNIT,120) IDK(I),BRFRAC(I),NME(I),(IDKPRD(J,I),J=1,5) IF (IDK(I).LT.0.OR.IDK(I).GT.NRES) IDK(I)=20 80 IF (IME.EQ.0.AND.NME(I).GE.100) NME(I)=0 ELSEIF (IOPT.EQ.3) THEN DO 90 I=1,NDKYS READ(JUNIT,130) CDK(I),BRFRAC(I),NME(I),(CDKPRD(J,I),J=1,5) IF (IME.EQ.0.AND.NME(I).GE.100) NME(I)=0 CALL HWUIDT(3,IDUM,IDK(I),CDK(I)) DO 90 J=1,5 90 CALL HWUIDT(3,IDUM,IDKPRD(J,I),CDKPRD(J,I)) ELSE CALL HWWARN('HWIODK',101,*999) ENDIF ENDIF CLOSE(UNIT=JUNIT) 100 FORMAT(1X,I4) 110 FORMAT(1X,I7,1X,F7.5,1X,I3,5(1X,I7)) 120 FORMAT(1X,I3,1X,F7.5,6(1X,I3)) 130 FORMAT(1X,A8,1X,F7.5,1X,I3,5(1X,A8)) 999 RETURN END