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
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
CALL GBOOK1(1,'SHAT',100,10.0D0,200.0D0)
CALL GBOOK1(2,'SHAT',100,10.0D0,2000.0D0)
WGTSM = 0.0D0
WGTSQ = 0.0D0
C--code to compute the average weight
DO I=1,NPOINT
CALL ZWGHT(WGT)
WGTSM = WGTSM+WGT
WGTSQ = WGTSQ+WGT**2
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'
C--final manipulation of the histograms
CALL GSCALE(1,1.0D0/DBLE(NPOINT)/(200.0D0-10.0D0)*100.0D0)
CALL GSCALE(2,1.0D0/DBLE(NPOINT)/(2000.0D0-10.0D0)*100.0D0)
CALL GTOPDR(1,1,1,1)
CALL GTOPDR(2,1,1,1)
END
C-----------------------------------------------------------------------
SUBROUTINE ZWGHT(WGT)
C-----------------------------------------------------------------------
C Subroutine to return the weight for a given phase space point
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--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)
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
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--plot a couple of histograms
CALL GFILL1(1,SQRT(SH),WGT)
CALL GFILL1(2,SQRT(SH),WGT)
END