4 *CMZ :- -20/10/99 09:46:43 by Peter Richardson
6 *-- Author : Bryan Webber, modified by Kosuke Odagiri
8 C-----------------------------------------------------------------------
10 SUBROUTINE HWISSP(LRSUSY)
12 C-----------------------------------------------------------------------
14 C Reads in SUSY particle properties and decays,
16 C in format generated by ISAWIG
18 C-----------------------------------------------------------------------
20 INCLUDE 'HERWIG61.INC'
22 INTEGER I,J,K,IH,IHW,NSSP,NDEC,MDKYS,LRSUSY
24 DOUBLE PRECISION BETAH, WEINCOS,WEINSIN, MW,MZ, RMMAX
26 DOUBLE PRECISION FTM,FTMUU(4),FTMDD(4),FTMTT(4),FTMBB(4),FTMU,FTMD
28 DOUBLE PRECISION YTM,YTM1,DTERM(4), SQHF,SNBCSB,MZSW2
32 EQUIVALENCE (MW,RMASS(198)), (MZ,RMASS(200))
50 C--reset susy input flag
52 IF (LRSUSY.LT.0) CALL HWWARN('HWISSP',500,*999)
58 C Input SUSY particle + top quark table
64 10 FORMAT (10X,'Reading in SUSY data from unit',I3)
66 READ (LRSUSY,'(I4)') NSSP
70 RMMAX=SQRT(HALF*(EBEAM1*EBEAM2+PBEAM1*PBEAM2))
76 READ (LRSUSY,1) IHW,RMASS(IHW),RLTIM(IHW)
78 C Negative gaugino mass means physical field is gamma_5*psi
82 IF ((IHW.GE.450).AND.(IHW.LE.457)) THEN
88 ZSGNSS(J)=RMASS(IHW)/ABS(RMASS(IHW))
90 ELSEIF (IHW.LE.455) THEN
94 WSGNSS(J)=RMASS(IHW)/ABS(RMASS(IHW))
98 RMASS(IHW)=ABS(RMASS(IHW))
102 IF (ABS(IDPDG(IHW)).GT.1000000.AND.(RMASS(IHW).NE.ZERO))
104 & RMMNSS=MIN(RMMNSS,RMASS(IHW))
106 IF (IHW.GT.NRES) THEN
108 IF (IHW.GT.NMXRES) CALL HWWARN('HWISSP',501,*999)
116 XLMNSS=TWO*LOG(RMMNSS/RMMAX)
118 1 FORMAT(I5,F12.4,E15.5)
128 READ (LRSUSY,'(I4)') NDEC
136 IF (NDKYS.GT.NMXDKS) CALL HWWARN('HWISSP',100,*999)
138 READ (LRSUSY,11) IDK(NDKYS),BRFRAC(NDKYS),NME(NDKYS),
140 & (IDKPRD(K,NDKYS),K=1,5)
142 11 FORMAT(I6,F16.8,6I6)
152 C Mixings and other SUSY parameters
156 READ (LRSUSY,'(2F16.8)') TANB,ALPHAH
160 READ (LRSUSY,13) ZMXNSS(I,1),ZMXNSS(I,2),ZMXNSS(I,3),ZMXNSS(I,4)
164 WEINSIN = SQRT(SWEIN)
166 WEINCOS = SQRT(1.-SWEIN)
170 ZMIXSS(I,1) = WEINCOS*ZMXNSS(I,1)+WEINSIN*ZMXNSS(I,2)
172 ZMIXSS(I,2) = -WEINSIN*ZMXNSS(I,1)+WEINCOS*ZMXNSS(I,2)
174 ZMIXSS(I,3) = ZMXNSS(I,3)
176 ZMIXSS(I,4) = ZMXNSS(I,4)
182 IF ((J.LE.6).OR.(J.GE.11)) THEN
184 LFCH(J)=VFCH(J,1)+AFCH(J,1)
186 RFCH(J)=VFCH(J,1)-AFCH(J,1)
190 SLFCH(J,I)= ZMIXSS(I,1)*QFCH(J)+ZMIXSS(I,2)*LFCH(J)
192 SRFCH(J,I)=-ZMIXSS(I,1)*QFCH(J)-ZMIXSS(I,2)*RFCH(J)
200 READ (LRSUSY,13) WMXVSS(1,1),WMXVSS(1,2), WMXVSS(2,1),WMXVSS(2,2)
202 READ (LRSUSY,13) WMXUSS(1,1),WMXUSS(1,2), WMXUSS(2,1),WMXUSS(2,2)
204 READ (LRSUSY,'(3F16.8)') THETAT,THETAB,THETAL
206 READ (LRSUSY,'(3F16.8)') ATSS,ABSS,ALSS
208 READ (LRSUSY,'( F16.8)') MUSS
230 QMIXSS(6,1,1)= COS(THETAT)
232 QMIXSS(6,1,2)= SIN(THETAT)
234 QMIXSS(6,2,1)=-QMIXSS(6,1,2)
236 QMIXSS(6,2,2)= QMIXSS(6,1,1)
238 QMIXSS(5,1,1)= COS(THETAB)
240 QMIXSS(5,1,2)= SIN(THETAB)
242 QMIXSS(5,2,1)=-QMIXSS(5,1,2)
244 QMIXSS(5,2,2)= QMIXSS(5,1,1)
246 LMIXSS(5,1,1)= COS(THETAL)
248 LMIXSS(5,1,2)= SIN(THETAL)
250 LMIXSS(5,2,1)=-LMIXSS(5,1,2)
252 LMIXSS(5,2,2)= LMIXSS(5,1,1)
254 C--Evaluating Higgs parameters and couplings
260 COSBPA=COS(BETAH+ALPHAH)
262 SINBPA=SIN(BETAH+ALPHAH)
264 COSBMA=COS(BETAH-ALPHAH)
266 SINBMA=SIN(BETAH-ALPHAH)
306 MZSW2 = MZ**2 * SQRT(SWEIN*(ONE-SWEIN))
308 DTERM(1) =-SINBPA*MZSW2
310 DTERM(2) = COSBPA*MZSW2
314 FTMUU(1) =-MUSS*SINA/SINB
316 FTMUU(2) = MUSS*COSA/SINB
322 FTMTT(1) =-ATSS*COSA/SINB
324 FTMTT(2) =-ATSS*SINA/SINB
330 FTMDD(1) = MUSS*COSA/COSB
332 FTMDD(2) = MUSS*SINA/COSB
338 FTMBB(1) = ABSS*SINA/COSB
340 FTMBB(2) =-ABSS*COSA/COSB
354 IF (I.EQ.5) FTMU=FTMU+FTMTT(IH)
356 IF (I.EQ.5) FTMD=FTMD+FTMBB(IH)
358 IF (MOD(I,2).EQ.0) THEN
374 GHSQSS(IH,I,1,1) = ZERO
376 GHSQSS(IH,I,2,2) = ZERO
378 GHSQSS(IH,I,1,2) = FTM*HALF*RMASS(I)/MW
380 GHSQSS(IH,I,2,1) = - GHSQSS(IH,I,1,2)
384 ELSEIF (IH.EQ.4) THEN
394 IF (MOD(I,2).EQ.1) THEN
396 GHSQSS(IH,I,J,K)=SQHF*(
398 & RMASS(I )*FTMD*QMIXSS(I,2,J)*QMIXSS(I+1,1,K)
400 & +RMASS(I+1)*FTMU*QMIXSS(I,1,J)*QMIXSS(I+1,2,K)
402 & +( MW**2*TWO*SNBCSB-RMASS(I+1)**2*COTB
404 & -RMASS(I )**2*TANB )*QMIXSS(I,1,J)*QMIXSS(I+1,1,K)
406 & -RMASS(I)*RMASS(I+1)/SNBCSB
408 & *QMIXSS(I,2,J)*QMIXSS(I+1,2,K) ) / MW
412 GHSQSS(IH,I,J,K)=GHSQSS(IH,I-1,K,J)
428 IF (J.EQ.K) YTM1=YTM*RMASS(I)**2
430 GHSQSS(IH,I,J,K)=( YTM1
432 & +( LFCH(I)*QMIXSS(I,1,J)*QMIXSS(I,1,K)
434 & -RFCH(I)*QMIXSS(I,2,J)*QMIXSS(I,2,K) )*DTERM(IH)
436 & +FTM*HALF*RMASS(I)*(QMIXSS(I,1,J)*QMIXSS(I,2,K)
438 & +QMIXSS(I,2,J)*QMIXSS(I,1,K)) ) / MW
452 READ (LRSUSY,'(L5)') RPARTY
456 READ(LRSUSY,20) (((LAMDA1(I,J,K),K=1,3),J=1,3),I=1,3)
458 READ(LRSUSY,20) (((LAMDA2(I,J,K),K=1,3),J=1,3),I=1,3)
460 READ(LRSUSY,20) (((LAMDA3(I,J,K),K=1,3),J=1,3),I=1,3)
474 *CMZ :- -04/05/99 14.28.59 by Bryan Webber
476 *-- Author : Bryan Webber
478 C-----------------------------------------------------------------------
482 C-----------------------------------------------------------------------
484 C IPROC = 1000,... ADDS SOFT UNDERLYING EVENT
486 C = 8000: CREATES MINIMUM-BIAS EVENT
488 C SUPPRESSED BY ADDING 10000 TO IPROC
490 C-----------------------------------------------------------------------
492 INCLUDE 'HERWIG61.INC'
494 DOUBLE PRECISION HWREXP,ENFAC,TECM,SECM,SUMM,EMCL,BMP(5),BMR(3,3)
496 INTEGER HWRINT,NETC,IBT,IDBT,ID1,ID2,ID3,KHEP,LHEP,NTRY,ICMS,
498 & NPPBAR,MCHT,JCL,JD1,JD2,JD3,ICH,MODC,NCHT,INHEP(2),
502 EXTERNAL HWREXP,HWRINT
504 IF (IERROR.NE.0) RETURN
506 IF (.NOT.GENSOF) GOTO 990
508 IF (IPROC.EQ.8000) THEN
510 C---SET UP BEAM AND TARGET CLUSTERS
518 IF (JDAHEP(1,IBT).NE.0) JBT=JDAHEP(1,IBT)
522 IF (IDBT.EQ.73.OR.IDBT.EQ.75) THEN
524 INID(1,IBT)=HWRINT(1,2)
528 ELSEIF (IDBT.EQ.91.OR.IDBT.EQ.93) THEN
532 INID(2,IBT)=HWRINT(7,8)
534 ELSEIF (IDBT.EQ.30) THEN
536 INID(1,IBT)=HWRINT(1,2)
540 ELSEIF (IDBT.EQ.38) THEN
544 INID(2,IBT)=HWRINT(7,8)
546 ELSEIF (IDBT.EQ.34) THEN
550 INID(2,IBT)=HWRINT(7,8)
552 ELSEIF (IDBT.EQ.46) THEN
554 INID(1,IBT)=HWRINT(1,2)
558 ELSEIF (IDBT.EQ.59) THEN
560 INID(1,IBT)=HWRINT(1,2)
562 INID(2,IBT)=HWRINT(7,8)
566 CALL HWWARN('HWMEVT',100,*999)
570 NETC=NETC+ICHRG(IDBT)
572 & -(ICHRG(INID(1,IBT))+ICHRG(INID(2,IBT)))/3
580 ISTHEP(NHEP+IBT)=163+IBT
582 JMOHEP(1,NHEP+IBT)=JBT
590 ELSEIF (NETC.EQ.-1) THEN
594 ELSEIF (NETC.EQ.1) THEN
610 IF (JDAHEP(1,IBT).NE.0) JBT=JDAHEP(1,IBT)
612 CALL HWVEQU(5,PHEP(1,JBT),PHEP(1,NHEP))
618 C---FIND BEAM AND TARGET CLUSTERS
624 IF (ISTHEP(KHEP).EQ.163+IBT) THEN
628 INID(1,IBT)=IDHW(JMOHEP(1,KHEP))
630 INID(2,IBT)=IDHW(JMOHEP(2,KHEP))
638 C---COULDN'T FIND ONE
646 C---TEST FOR BOTH FOUND
648 IF (INHEP(1).EQ.0) JCL=INHEP(2)
650 IF (INHEP(2).EQ.0) JCL=INHEP(1)
652 IF (JCL.EQ.0) CALL HWWARN('HWMEVT',101,*999)
678 C---FIND SOFT CM MOMENTUM AND MULTIPLICITY
684 IF (NHEP.GT.NMXHEP) CALL HWWARN('HWMEVT',102,*999)
694 CALL HWVSUM(4,PHEP(1,INHEP(1)),PHEP(1,INHEP(2)),PHEP(1,NHEP))
696 CALL HWUMAS(PHEP(1,NHEP))
700 IF (IPROC/1000.EQ.9.OR.IPROC/1000.EQ.5) THEN
710 C---CHOOSE MULTIPLICITY
712 25 CALL HWMULT(SECM,NPPBAR)
732 IF (NHEP.GT.NMXHEP) CALL HWWARN('HWMEVT',103,*999)
752 IF (NCL.EQ.3) ID1=ID3
764 CALL HWVZRO(3,PHEP(1,JCL))
766 PHEP(4,JCL)=RMASS(ID1)+RMASS(ID2)+PMBM1+HWREXP(TWO/PMBM2)
768 PHEP(5,JCL)=PHEP(4,JCL)
770 C---HADRONIZE AND DECAY CLUSTERS
772 CALL HWCFLA(ID1,ID2,JD1,JD2)
774 CALL HWCHAD(JCL,JD1,JD2,JD3)
776 IF (IERROR.NE.0) RETURN
780 EMCL=RMASS(IDHW(NHEP))
782 IF (PHEP(4,JCL).NE.EMCL) THEN
810 IF (IERROR.NE.0) RETURN
812 C---CHECK CHARGED MULTIPLICITY
818 IF (ISTHEP(KHEP).EQ.1) THEN
820 ICH=ICHRG(IDHW(KHEP))
836 NCHT=NPPBAR+NETC+ABS(MODC)
840 ELSEIF (NCL.EQ.2) THEN
844 IF (NCHT.LT.0) NCHT=NCHT+2
848 IF (MCHT.LT.NCHT) THEN
852 ELSEIF (MCHT.GT.NCHT) THEN
854 IF (MOD(NTRY,50).EQ.0) GOTO 25
856 IF (NTRY.LT.NSTRY) GOTO 30
858 C---NO PHASE SPACE FOR SOFT EVENT
862 IF (IPROC.EQ.8000) THEN
864 C---MINIMUM BIAS: RELABEL BEAM AND TARGET CLUSTERS
874 IDHEP(KHEP)=IDHEP(LHEP)
876 IDHW(KHEP)=IDHW(LHEP)
882 C---UNDERLYING EVENT: DECAY THEM
902 C---GENERATE CLUSTER MOMENTA IN CLUSTER CM
904 C FRAME. N.B. SECOND CLUSTER IS TARGET
906 IF (SUMM.GT.TECM) GOTO 25
910 IF (NCL.EQ.0) GOTO 25
914 C---ROTATE & BOOST CLUSTERS & DECAY PRODUCTS
916 CALL HWULOF(PHEP(1,ICMS),PHEP(1,INHEP(1)),BMP)
918 CALL HWUROT(BMP, ONE,ZERO,BMR)
920 C---BMR PUTS BEAM ALONG Z AXIS (WE WANT INVERSE)
922 DO 70 KHEP=ICMS+1,NHEP
924 IF (ISTHEP(KHEP).GT.180.AND.ISTHEP(KHEP).LT.190
926 $ .AND.JMOHEP(1,KHEP).EQ.ICMS) THEN
928 ISTHEP(KHEP)=ISTHEP(KHEP)+3
934 CALL HWUROB(BMR,PPCL(1,JCL),PPCL(1,JCL))
936 CALL HWULOB(PHEP(1,ICMS),PPCL(1,JCL),PPCL(1,JCL))
938 C---NOW PPCL(*,JCL) IS LAB MOMENTUM OF JTH CLUSTER
942 CALL HWULOB(PPCL(1,JCL),PHEP(1,KHEP),PHEP(1,KHEP))
950 JMOHEP(1,ICMS)=INHEP(1)
952 JMOHEP(2,ICMS)=INHEP(2)
954 JDAHEP(1,INHEP(1))=ICMS
958 JDAHEP(1,INHEP(2))=ICMS
962 JDAHEP(1,ICMS)=ICMS+1
974 *CMZ :- -04/05/99 14.17.04 by Bryan Webber
976 *-- Author : David Ward, modified by Bryan Webber
978 C-----------------------------------------------------------------------
980 SUBROUTINE HWMLPS(TECM)
982 C-----------------------------------------------------------------------
984 C GENERATES CYLINDRICAL PHASE SPACE USING THE METHOD OF JADACH
986 C RETURNS WITH NCL=0 IF UNSUCCESSFUL
988 C-----------------------------------------------------------------------
990 INCLUDE 'HERWIG61.INC'
992 DOUBLE PRECISION HWREXT,HWRUNG,HWUSQR,TECM,ESS,ALOGS,EPS,SUMX,
994 & SUMY,PT,PX,PY,PT2,SUMPT2,SUMTM,XIMIN,XIMAX,YY,SUM1,SUM2,SUM3,
996 & SUM4,EX,FY,DD,DYY,ZZ,E1,TM,SLOP,XI(NMXCL)
998 INTEGER NTRY,I,NIT,IY(NMXCL),IDP
1000 EXTERNAL HWREXT,HWRUNG,HWUSQR
1002 IF (NCL.GT.NMXCL) THEN
1004 CALL HWWARN('HWMLPS',1,*999)
1020 IF (NTRY.GT.NSTRY) THEN
1034 C---Pt distribution of form exp(-b*Mt)
1036 C---Factors for pt slopes to fit data. IDCL contains the type of
1038 C q-qbar pair produced in this cluster (0 if 1-particle cluster).
1046 ELSEIF(IDP.EQ.3.OR.IDP.EQ.10) THEN
1050 ELSEIF(IDP.GT.3.AND.IDP.LE.9) THEN
1056 CALL HWWARN('HWMLPS',200,*999)
1060 PT=HWREXT(PPCL(5,I),SLOP)
1062 PT=HWUSQR(PT**2-PPCL(5,I)**2)
1064 CALL HWRAZM(PT,PX,PY)
1072 12 SUMY=SUMY+PPCL(2,I)
1084 PPCL(1,I)=PPCL(1,I)-SUMX
1086 PPCL(2,I)=PPCL(2,I)-SUMY
1088 PT2=PPCL(1,I)**2+PPCL(2,I)**2
1092 C---STORE TRANSVERSE MASS IN PPCL(3,I) TEMPORARILY
1094 PPCL(3,I)=SQRT(PT2+PPCL(5,I)**2)
1096 13 SUMTM=SUMTM+PPCL(3,I)
1098 IF (SUMTM.GT.TECM) GOTO 11
1102 C---Form of "reduced rapidity" distribution
1104 XI(I)=HWRUNG(0.6*ONE,ONE)
1108 CALL HWUSOR(XI,NCL,IY,1)
1114 C---N.B. TARGET CLUSTER IS SECOND
1120 XI(I+1)=(XI(I)-XIMIN)/XIMAX
1126 YY=LOG(ESS/(PPCL(3,1)*PPCL(3,2)))
1148 SUM3=SUM3+(TM*EX)*XI(I)
1150 19 SUM4=SUM4+(TM/EX)*XI(I)
1152 FY=ALOGS-LOG(SUM1*SUM2)
1154 DD=(SUM3*SUM2-SUM1*SUM4)/(SUM1*SUM2)
1158 IF(ABS(DYY/YY).LT.EPS) GOTO 20
1162 C---Y ITERATIONS EXCEEDED - TRY AGAIN
1164 IF (NTRY.LT.100) GOTO 11
1168 IF (EPS.GT.ONE) CALL HWWARN('HWMLPS',100,*999)
1170 CALL HWWARN('HWMLPS',50,*11)
1182 PPCL(3,I)=(0.5*TM)*((1./E1)-E1)
1184 PPCL(4,I)=(0.5*TM)*((1./E1)+E1)
1192 *CMZ :- -26/04/91 11.11.55 by Bryan Webber
1194 *-- Author : David Ward, modified by Bryan Webber
1196 C-----------------------------------------------------------------------
1198 FUNCTION HWMNBI(N,AVNCH,EK)
1200 C-----------------------------------------------------------------------
1202 C---Computes negative binomial probability
1204 C-----------------------------------------------------------------------
1206 DOUBLE PRECISION HWMNBI,AVNCH,EK,R
1218 HWMNBI=(1.+R)**(-EK)
1224 HWMNBI=HWMNBI*R*(EK+I-1)/I