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