C-----------------------------------------------------------------------
      PROGRAM ZGEN
C-----------------------------------------------------------------------
C     Example main program
C-----------------------------------------------------------------------
      IMPLICIT NONE
C--Les Houches run common block
      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)
C--common block to pass the parameters
      DOUBLE PRECISION SMIN,SMAX
      INTEGER IPR
      COMMON /LPARAM/SMIN,SMAX,IPR
C--local variables
      INTEGER I,NPOINT
      DOUBLE PRECISION WGT,WGTSM,WGTSQ,MAXWGT
      PARAMETER(NPOINT=100000)
C--beam particles and energies (only proton=2212 and antiproton=-2212)
      IDBMUP(1) = 2212
      IDBMUP(2) = 2212
      EBMUP(1)  = 7000.0D0
      EBMUP(2)  = 7000.0D0
C--pdf's for the beams (use CTEQ5L)
C--MUST BE THE SAME FOR BOTH BEAMS
      PDFGUP(1) = 4
      PDFGUP(2) = 4
      PDFSUP(1) = 46
      PDFSUP(2) = 46
C--integration parameters
      SMIN = 10.0D0**2
      SMAX = (EBMUP(1)+EBMUP(2))**2
      print *,'limits on sqrt(sh)',sqrt(smin),sqrt(smax)
C--which process the following are supported
C      100 all fermion-antifermion production
C      110 all lepton production
C      120 all quark production
C      130+I production of the fermion i=(1,6 =e,ve,mu,vmu,tau,vtau)
C                                        (7,12-d,u,s,c,b,t)
      IPR = 131
C--set up the parameters of the Little Higgs Model
      CALL LITTLE
      WGTSM = 0.0D0
      WGTSQ = 0.0D0
      MAXWGT = 0.0D0
C--code to compute the average weight
      DO I=1,NPOINT
         CALL ZWGHT(WGT)
         WGTSM = WGTSM+WGT
         WGTSQ = WGTSQ+WGT**2
         MAXWGT=MAX(MAXWGT,WGT)
      ENDDO
      WGTSM = WGTSM/DBLE(NPOINT)
      WGTSQ = MAX(0.0D0,WGTSQ/DBLE(NPOINT)-WGTSM**2)
      WGTSQ = SQRT(WGTSQ/DBLE(NPOINT))
      print *,'Cross section = ',wgtsm,'+/-',wgtsq,' pb'
      END
C-----------------------------------------------------------------------
      SUBROUTINE ZWGHT(WGT)
C-----------------------------------------------------------------------
C     Subroutine to return the weight for a given phase space point
C-----------------------------------------------------------------------
      IMPLICIT NONE
C--Les Houches  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--Les Houches run common block
      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)
C--common block to pass the parameters
      DOUBLE PRECISION SMIN,SMAX
      INTEGER IPR
      COMMON /LPARAM/SMIN,SMAX,IPR
C--little Higgs common block
C--common block with Liggle Higgs parameters
      DOUBLE PRECISION MW,MZ,MF,ALPHA,SW,MT,FH,THETA,THETAP,G,GP,E,PI,
     &     CW,VEV,VEVP,ST,CT,STP,CTP,LAM,MZH,MWH,MAH,GL,GR,WIDTH
      COMMON /LHIGGS/MW,MZ,MF(12),ALPHA,SW,MT,FH,THETA,THETAP,G,GP,E,PI,
     &     CW,VEV,VEVP,ST,CT,STP,CTP,
     &     LAM(2),MZH,MWH,MAH,GL(14,6),GR(14,6),WIDTH(6)
C--local variables
      INTEGER NCH
      PARAMETER(NCH=4)
      INTEGER ITYPE(NCH),I,IN,IOUT,IFERM(3),J
      DOUBLE PRECISION MCH(NCH),MCH2(NCH),MGAM(NCH),MGAM2(NCH),ALP(NCH),
     &     WI(NCH),LIM(2,NCH),RANGEN
      DOUBLE PRECISION WGT,CTHETA,PHI,PIFAC,RANUNI,MEWGT,
     &     GEV2PB,TEMP,SH,XX(2),FJAC,DISF(13,2),ME(2,12,6),UH,TH,PCM,
     &     STHETA,CPHI,SPHI
      DOUBLE COMPLEX GLL,GRR,GLR,GRL,ZI,PROP
      PARAMETER(ZI=(0.0D0,1.0D0))
      EXTERNAL RANGEN,RANUNI
      LOGICAL FIRST
      DATA FIRST/.TRUE./
      SAVE FIRST,PIFAC,ITYPE,WI,ALP,MCH,MCH2,MGAM,MGAM2,LIM,
     &     GEV2PB,IFERM
C--useful parameters and initilisation
      IF(FIRST) THEN
C--constants
         PIFAC = ACOS(-1.0D0)
         FIRST = .FALSE.
C--conversion to picobarn from GeV
         GEV2PB=389379.D3
C--which processes to generate (IFERM controls the loop over final state fermions)
C--(e,ve,mu,vmu,tau,nvtau,d,u,s,c,b,t)
         IF(IPR.EQ.100) THEN
            IFERM(1) = 1
            IFERM(2) = 12
         ELSEIF(IPR.EQ.110) THEN
            IFERM(1) = 1
            IFERM(2) = 6
         ELSEIF(IPR.EQ.120) THEN
            IFERM(1) = 7
            IFERM(2) = 12
         ELSEIF(IPR.GE.131.AND.IPR.LE.142) THEN
            IFERM(1) = IPR-130
            IFERM(2) = IPR-130
         ELSE
            print *,'unknown process stopping'
            STOP
         ENDIF
         print *,'selected process code =',ipr
         print *,'limits on fermion loop',iferm
C--initialisation of the multichannel integration
C--first a channel to smooth the photon
         ITYPE(1) = 1
         WI   (1) = 0.4D0
         ALP  (1) = 2.0D0
         MCH  (1) = 0.0D0
         MCH2 (1) = 0.0D0
         MGAM (1) = 0.0D0
         MGAM2(1) = 0.0D0
C--second a channel to smooth the Z
         ITYPE(2) = 2
         WI   (2) = 0.4D0
         ALP  (2) = 0.0D0
         MCH  (2) = MZ
         MCH2 (2) = MZ**2
         MGAM (2) = MZ*WIDTH(2)
         MGAM2(2) = MGAM(2)**2
C--third a channel to smooth the heavy photon
         ITYPE(3) = 2
         WI   (3) = 0.1D0
         ALP  (3) = 0.0D0
         MCH  (3) = MAH
         MCH2 (3) = MAH**2
         MGAM (3) = MAH*WIDTH(3)
         MGAM2(3) = MGAM(3)**2
C--fourth a channel to smooth the heavy Z boson
         ITYPE(4) = 2
         WI   (4) = 0.1D0
         ALP  (4) = 0.0D0
         MCH  (4) = MZH
         MCH2 (4) = MZH**2
         MGAM (4) = MZH*WIDTH(3)
         MGAM2(4) = MGAM(4)**2
C--compute the limits on the variables
         CTHETA = 0.0D0
         DO I=1,NCH
            CTHETA = CTHETA+WI(I)
            IF(ITYPE(I).EQ.1) THEN
               LIM(1,I) = SMAX**(1.0D0-ALP(I))
               LIM(2,I) = SMIN**(1.0D0-ALP(I))
            ELSEIF(ITYPE(I).EQ.2) THEN
               LIM(1,I) = ATAN((SMIN-MCH2(I))/MGAM(I))
               LIM(2,I) = ATAN((SMAX-MCH2(I))/MGAM(I))
            ENDIF
         ENDDO
C--ensure weights add up to one
         CTHETA = 1.0D0/CTHETA
         DO I=1,NCH
            WI(I)=WI(I)*CTHETA
         ENDDO
      ENDIF
C--calculation of the phase space weight and momenta
C--first the values of theta and phi
      CTHETA = RANUNI(1,-1.0D0,1.0D0)
      PHI    = RANUNI(2,0.0D0,2.0D0*PIFAC)
C--weight for this piece of the process
      WGT = 4.0D0*PIFAC
C--calculation of sh
C--first pick the channel
      TEMP = RANGEN(3)
      DO I=1,NCH
         IF(WI(I).GT.TEMP) GOTO 10
         TEMP = TEMP-WI(I)
      ENDDO
C--generate the transformed variable
 10   TEMP = RANUNI(4,LIM(1,I),LIM(2,I))
C--calculate sh from transformed varaible
      IF(ITYPE(I).EQ.1) THEN
         SH = TEMP**(1.0D0/(1.0D0-ALP(I)))
      ELSE
         SH = MGAM(I)*TAN(TEMP)+MCH2(I)
      ENDIF
C--compute the Jacobian function
      FJAC = 0.0D0
      DO I=1,NCH
         TEMP = LIM(2,I)-LIM(1,I)
         IF(ITYPE(I).EQ.1) THEN
            TEMP = 1.0D0/SH**ALP(I)/TEMP
         ELSE
            TEMP = MGAM(I)/TEMP/((SH-MCH2(I))**2+MGAM2(I))
         ENDIF
         FJAC = FJAC+WI(I)*TEMP
      ENDDO
C--include this in the phase-space weight
      WGT = WGT/FJAC
C--generate the value of x_1 and x_2
      TEMP = LOG(SH/SMAX)
      XX(1) = EXP(RANUNI(5,LOG(SH/SMAX),0.0D0))
      XX(2) = SH/SMAX/XX(1)
C--include this in the phase-space weight
      WGT = -WGT*TEMP
C--final factors for the phase space piece of the weight
      WGT = WGT/32.0D0/PIFAC**2/SH**2.5D0
C--call the structure functions
      CALL PDF(SQRT(SH),XX,DISF)
C--calculation of the matrix element 
      MEWGT = 0.0D0
      DO 100 IN  =1,2
      DO 100 IOUT=IFERM(1),IFERM(2)
C--final state particles kinematically allowed
      IF(SH.GT.4.0D0*MF(IOUT)**2) THEN
C--compute the centre-of-mass momentum
         PCM = 0.5D0*SQRT(SH-4.0D0*MF(IOUT)**2)
C--compute the mandelstam variables
C     ( as final state masses equal this is that-mf^2 and
C      uhat-mf^2 as this is what we need in the matrix element)
         UH  =-SQRT(SH)*(SQRT(PCM**2+MF(IOUT)**2)-PCM*CTHETA)
         TH  =-SQRT(SH)*(SQRT(PCM**2+MF(IOUT)**2)+PCM*CTHETA)
C--compute the couplings
         GLL = 0.0D0
         GRL = 0.0D0
         GRR = 0.0D0
         GLR = 0.0D0
         DO I=1,4
            PROP = 1.0D0/(SH-MCH2(I)+ZI*MGAM(I))
            GLL = GLL+GL(IN+6,I)*GL(IOUT,I)*PROP
            GRR = GRR+GR(IN+6,I)*GR(IOUT,I)*PROP
            GLR = GLR+GL(IN+6,I)*GR(IOUT,I)*PROP
            GRL = GRL+GR(IN+6,I)*GL(IOUT,I)*PROP
         ENDDO
         TEMP = DBLE( TH**2*(GRL*DCONJG(GRL)+GLR*DCONJG(GLR))
     &               +UH**2*(GRR*DCONJG(GRR)+GLL*DCONJG(GLL))
     &          +2.0D0*MF(IOUT)**2*(GRL*DCONJG(GRR)+GLR*DCONJG(GLL)))
         TEMP = TEMP*PCM
         IF(IOUT.LE.6) TEMP=TEMP/3.0D0
C--compute the matrix elements
         DO I=1,3
            ME(1,IOUT,2*I+IN-2) = DISF(2*I+IN-2,1)*DISF(2*I+IN+4,2)*TEMP
            ME(2,IOUT,2*I+IN-2) = DISF(2*I+IN-2,2)*DISF(2*I+IN+4,1)*TEMP
            MEWGT = MEWGT+ME(1,IOUT,2*I+IN-2)+ME(2,IOUT,2*I+IN-2)
         ENDDO
C--set to zero as not allowed
      ELSE
         DO I=1,3
            ME(1,IOUT,2*I+IN-2) = 0.0D0
            ME(2,IOUT,2*I+IN-2) = 0.0D0
         ENDDO
      ENDIF
 100  CONTINUE
C--multiply the two bits of the weight together
      WGT = WGT*MEWGT*GEV2PB
C--first workout what the process is
      MEWGT = MEWGT*RANGEN(5)
      DO IN=1,6
        DO IOUT=IFERM(1),IFERM(2)
           DO I=1,2
              IF(ME(I,IOUT,IN).GT.MEWGT) THEN
                 GOTO 200
              ELSE
                 MEWGT=MEWGT-ME(I,IOUT,IN)
              ENDIF
           ENDDO
        ENDDO
      ENDDO
 200  CONTINUE
C--set up identities and colour flow of the incoming particles
      IF(I.EQ.1) THEN
         IDUP(1) = IN
         IDUP(2) =-IN
         ICOLUP(1,1) = 501 
         ICOLUP(2,1) = 0
         ICOLUP(1,2) = 0
         ICOLUP(2,2) = 501
      ELSE
         IDUP(1) =-IN
         IDUP(2) = IN
         ICOLUP(1,1) = 0 
         ICOLUP(2,1) = 501
         ICOLUP(1,2) = 501
         ICOLUP(2,2) = 0
      ENDIF
C--status of the incoming particles
      ISTUP(1) = -1
      ISTUP(2) = -1
C--set up the momenta of the incoming particles
      DO J=1,2
         PUP(1,J) = 0.0D0
         PUP(2,J) = 0.0D0
         PUP(3,J) = EBMUP(J)*XX(J)
         PUP(4,J) = EBMUP(J)*XX(J)
         PUP(5,J) = 0.0D0
      ENDDO
      PUP(3,2) = -PUP(3,2)
C--add the CMF as a resonance
C--momentum of the centre-of-mass system
      PUP(1,3) = 0.0D0
      PUP(2,3) = 0.0D0
      PUP(3,3) = EBMUP(1)*XX(1)-EBMUP(2)*XX(2)
      PUP(4,3) = EBMUP(1)*XX(1)+EBMUP(2)*XX(2)
      PUP(5,3) = SQRT(SH)
C--status
      ISTUP(3) = 2
C--mothers
      MOTHUP(1,3) = 1
      MOTHUP(2,3) = 2
C--colours
      ICOLUP(1,3) = 0
      ICOLUP(2,3) = 0
C--id, call it a Z
      IDUP(3)     = 23
C--momentum of the outgoing particles
      PCM = 0.5D0*SQRT(SH-4.0D0*MF(IOUT)**2)
C--change sign of costheta if antiquark first
      IF(I.EQ.2) CTHETA=-CTHETA
C--momentum of outgoing particles in CMF frame
      STHETA   = SQRT(1.0D0-CTHETA**2)
      CPHI     = COS(PHI)
      SPHI     = SIN(PHI)
      PUP(1,4) = PCM*STHETA*CPHI
      PUP(2,4) = PCM*STHETA*SPHI
      PUP(3,4) = PCM*CTHETA
      PUP(4,4) = SQRT(PCM**2+MF(IOUT)**2)
      PUP(5,4) = MF(IOUT)
      PUP(1,5) =-PCM*STHETA*CPHI
      PUP(2,5) =-PCM*STHETA*SPHI
      PUP(3,5) =-PCM*CTHETA
      PUP(4,5) = SQRT(PCM**2+MF(IOUT)**2)
      PUP(5,5) = MF(IOUT)
C--boost these to nthe lab frame
      CALL BSTLOB(PUP(1,3),PUP(1,4),PUP(1,4))
      CALL BSTLOB(PUP(1,3),PUP(1,5),PUP(1,5))
C--ids of outgoing particles
      IF(IOUT.LE.6) THEN
         IDUP(4) = IOUT+10
      ELSE
         IDUP(4) = IOUT-6
      ENDIF
      IDUP(5) =-IDUP(4)
C--status of outgoing particle
      ISTUP(4) = 1
      ISTUP(5) = 1
C--mothers
      MOTHUP(1,4) = 3
      MOTHUP(2,4) = 3
      MOTHUP(1,5) = 3
      MOTHUP(2,5) = 3
C--colour flows
C--zero if leptons
      IF(IOUT.LE.6) THEN
         ICOLUP(1,4) = 0
         ICOLUP(2,4) = 0
         ICOLUP(1,5) = 0
         ICOLUP(2,5) = 0
C--normal for quarks
      ELSE
         ICOLUP(1,4) = 502
         ICOLUP(2,4) = 0
         ICOLUP(1,5) = 0
         ICOLUP(2,5) = 502
      ENDIF
C--number of particles
      NUP = 5
C--scale
      SCALUP=SQRT(SH)
      END