C...All real arithmetic in double precision. IMPLICIT DOUBLE PRECISION(A-H, O-Z) C...Three Pythia functions return integers, so need declaring. INTEGER PYK,PYCHGE,PYCOMP C...EXTERNAL statement links PYDATA on most machines. EXTERNAL PYDATA C...Commonblocks. C...The event record. COMMON/PYJETS/N,NPAD,K(4000,5),P(4000,5),V(4000,5) C...Parameters. COMMON/PYDAT1/MSTU(200),PARU(200),MSTJ(200),PARJ(200) C...Particle properties + some flavour parameters. COMMON/PYDAT2/KCHG(500,4),PMAS(500,4),PARF(2000),VCKM(4,4) C...Decay information. COMMON/PYDAT3/MDCY(500,3),MDME(4000,2),BRAT(4000),KFDP(4000,5) C...Selection of hard scattering subprocesses. COMMON/PYSUBS/MSEL,MSELPD,MSUB(500),KFIN(2,-40:40),CKIN(200) C...Parameters. COMMON/PYPARS/MSTP(200),PARP(200),MSTI(200),PARI(200) common/pyint1/mint(400),vint(400) C...Counter for number of generated events of each type. dimension ncount(2) data ncount/2*0/ C...Proton Variable dimension xpp(-25:25),xppbar(-25:25),term(20) C Compile C g77 -o test test.o /home/cdfsoft/products/pdf/v8_04/Linux+2/lib/libpdf.a /home/tsuno/lib/libpythia6203.a -L/usr/products/cern/2000/lib/ -lmathlib -lpacklib -lkernlib -L/usr/lib -lnsl -lm C...Number of events. nev = 100000 ! Number of events C...Switch off unnecessary aspects: initial and final state C...showers, multiple interactions, hadronization. mstp(61) = 0 ! Initial state radiation OFF mstp(71) = 0 ! Final state radiation OFF mstp(81) = 0 ! Multiple interaction OFF mstp(111) = 0 ! Hadronization OFF C...Set PDF and Coupling. mstp(52) = 1 ! internal/external : 1/2 swich mstp(51) = 7 ! CTEQ5L mstp(58) = 4 ! Nr. of flavor in PDF mstu(101) = 1 ! running alpha_em mstu(111) = 1 ! first-order running alpha_s mstu(112) = 5 ! Nr. of flavors assumed in alpha_s paru(112) = 0.146d0 ! Lambda used in alpha_s C...Set the weak mixing angle C paru(102) = 0.22395 C...Switch on process. CCCC msel = 19 CCC msel = 0 msub(102) = 1 pmas(25,1) = 120.0d0 mstp(49) = 0 PARU(102) = 0.232 do i = 210,288 mdme(i,1) = 0 enddo mdme(214,1) = 1 CCCCCC pmas(24,2) = 5.0d0 CCC CCC CCCC msub(102) = 1 CCCC msub(152) = 1 CCC msub(157) = 1 CCC CCC PARP(48) = 10000.0 CCC mstp(49) = 0 CCC CCC mstp(4) = 2 CCC pmas(25,1) = 90.9d0 CCC CCC ckin(1) = 60 CCC CCCC paru(141) = 45.0d0 CCCC paru(141) = 15.0d0 CCC paru(141) = 60.0d0 CCCC paru(141) = 35.0d0 CCCC paru(141) = 45.0d0 CCC CCC pmas(5,1) = 4.8d0 CCC CCCCCC do i = 210,288 CCCCCC mdme(i,1) = 0 CCCCCC enddo CCCCCC mdme(214,1) = 1 C...Initialize. CCC call pyinit('CMS','p','pbar',2000.0d0) call pyinit('USER',' ',' ',0,0d0) write(*,*) 'mass',pmas(25,1),pmas(36,1),pmas(35,1) write(*,*) 'width',pmas(25,2),pmas(36,2),pmas(35,2) C...Event loop. do iev = 1,nev call pyevnt CCC CCCC write(*,*) pmas(25,1),pmas(25,2) CCC C as = pyalps(120.0d0**2) C write(*,*) as isub = msti(1) icase = 1 if (isub.eq.200) icase = 2 ncount(icase) = ncount(icase) + 1 if (ncount(icase).le.3) then write(6,*) ' Following event is subprocess',isub call pylist(1) endif C...Data Store if (mod(iev,2000).eq.0) then write(*,*) iev endif C...End of loop over events. enddo C----------------------------------------------------------------- C...Third section: produce output and end. C...Cross section table. CALL PYSTAT(1) write(*,*) 'mass',pmas(25,1),pmas(36,1),pmas(35,1) write(*,*) 'width',pmas(25,2),pmas(36,2),pmas(35,2) END C########################################################### SUBROUTINE UPINIT C...All real arithmetic in double precision. IMPLICIT NONE C...Include HEP common parameter. 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) SAVE /HEPRUP/ C...PDF common parameter. CHARACTER*20 PARM(20) DOUBLE PRECISION VALUE(20) C...Set incoming beams: Tevatron Run II. IDBMUP(1) = 2212 ! proton IDBMUP(2) = -2212 ! anti-proton EBMUP(1) = 980D0 ! Beam energy EBMUP(2) = 980D0 C...Initialization of PDFLIB. PARM(1) = 'Init0' VALUE(1) = 0.0d0 CALL PDFSET(PARM, VALUE) C...Set PDF's of incoming beams: CTEQ 4L. PARM(1) = 'Nptype' VALUE(1) = 1 PARM(2) = 'Ngroup' VALUE(2) = 4 PARM(3) = 'Nset' VALUE(3) = 46 CALL PDFSET(PARM, VALUE) C...Set PDF's of incoming beams: CTEQ 5L. C...Note that Pythia will not look at PDFGUP and PDFSUP. PDFGUP(1) = VALUE(2) ! PDF group PDFSUP(1) = VALUE(3) ! PDF set PDFGUP(2) = PDFGUP(1) PDFSUP(2) = PDFSUP(1) C...Decide on weighting strategy: unweighted on input. IDWTUP = 1 C...Number of external processes. NPRUP = 1 ! Nr. of process LPRUP(1) = 200 ! proc. ID C XMAXUP(1) = 9.0D-1 ! maximum weight XMAXUP(1) = 6.0D-1 ! maximum weight C...Done. RETURN END SUBROUTINE UPEVNT C...All real arithmetic in double precision. C IMPLICIT NONE IMPLICIT DOUBLE PRECISION(A-H, O-Z) C...Include HEP common parameter. 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) SAVE /HEPRUP/ C...Include HEP common parameter. 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) SAVE /HEPEUP/ REAL*8 PI,GF,WCM,BETAF,Z1,Z2,PYR,ALPHAS2,PYMRUN REAL*8 SHAT,WHAT,AJACOB,X1,X2,BETA,GAMM,FACT,SUM REAL*8 AMH,AGH,MB,MT,MBRUN REAL*8 TAUMIN,TAUMAX,TAURE,GAMRE,THEMIN,THEMAX,THE,TAU REAL*8 YSTMIN,YSTMAX,YST REAL*8 COS3MIN,COS3MAX,COS3,PHI3 REAL*8 UPV1,DNV1,USEA1,DSEA1,STR1,CHM1,BOT1,TOP1,GL1 REAL*8 UPV2,DNV2,USEA2,DSEA2,STR2,CHM2,BOT2,TOP2,GL2 REAL*8 TQ,RQ,LQ,PQ1,PQ2,AQ1,AQ2,AP,AQ REAL*8 SIGIN,BWIG,RDCOR,SIGOUT REAL*8 E3,P3,SIN3 COMMON/PYDAT1/MSTU(200),PARU(200),MSTJ(200),PARJ(200) REAL*8 PYALPS,PYALEM,AMW REAL*8 XPP(-25:25),XPPBAR(-25:25) C...beta function. BETAF(Z1,Z2)=SQRT(1.0D0-2.0D0*(Z1+Z2)+(Z1-Z2)**2) C...Constant. PI = ACOS(-1.0D0) GF = 1.16639D-5 C...Set Mass and width AMH = 120.0D0 ! Higgs mass AGH = 0.00405 ! Width of Higgs MB = 4.8D0 ! b-quark mass MT = 175.0D0 ! top mass WCM = EBMUP(1) + EBMUP(2) ! tot. Beam energy C...Phase space Jacobian. AJACOB = 1.0D0 C...Set kinematic variables. TAUMIN = 4.0D0*MB**2/WCM**2 TAUMAX = 1.0D0 TAURE = AMH**2/WCM**2 GAMRE = AMH*AGH/WCM**2 THEMIN = ATAN((TAUMIN - TAURE)/GAMRE) THEMAX = ATAN((TAUMAX - TAURE)/GAMRE) THE = THEMIN + (THEMAX - THEMIN)*PYR(0) TAU = GAMRE*TAN(THE) + TAURE AJACOB = AJACOB*(THEMAX - THEMIN) $ *((TAU - TAURE)**2 + GAMRE**2)/(GAMRE*TAU) YSTMIN = LOG(SQRT(TAU)) YSTMAX = -YSTMIN YST = YSTMIN + (YSTMAX - YSTMIN)*PYR(0) AJACOB = AJACOB*(YSTMAX - YSTMIN) COS3MIN = -1.0D0 COS3MAX = 1.0D0 COS3 = COS3MIN + (COS3MAX - COS3MIN)*PYR(0) AJACOB = AJACOB*(COS3MAX - COS3MIN) PHI3 = PYR(0)*2.0D0*PI AJACOB = AJACOB*2.0D0*PI C...Set x1,x2,s^,w^ X1 = SQRT(TAU)*EXP(YST) X2 = SQRT(TAU)*EXP(-YST) SHAT = TAU*WCM**2 WHAT = SQRT(TAU)*WCM C...Set beta,gamma BETA = TANH(YST) GAMM = COSH(YST) FACT = 0.38937966D9/2.0D0/SHAT CCCC...External PDF CCC CALL STRUCTM(X1,WHAT, CCC $ UPV1,DNV1,USEA1,DSEA1,STR1,CHM1,BOT1,TOP1,GL1) CCC CALL STRUCTM(X2,WHAT, CCC $ UPV2,DNV2,USEA2,DSEA2,STR2,CHM2,BOT2,TOP2,GL2) CCCC...gluon from p, gluon from pbar CCC SUM = GL1*GL2 CALL PYPDFU( 2212,X1,SHAT,XPP) CALL PYPDFU(-2212,X2,SHAT,XPPBAR) SUM = XPP(0)*XPPBAR(0) AMW = 80.45D0 GF = 4.0D0*PI*PYALEM(SHAT)*SQRT(2.0D0)/8.0D0/AMW**2/PARU(102) AJACOB = AJACOB*SUM C...Phase space factor AJACOB = AJACOB*FACT/2.D0/PI**2*BETAF(MB**2/SHAT,MB**2/SHAT) C...Calculate Gamma(gg -> H). TQ = 4.0D0*MT**2/SHAT IF (TQ.LT.1) THEN RQ = SQRT(1.0D0-TQ) LQ = LOG((1.0D0+RQ)/(1.0D0-RQ)) PQ1 = -(LQ**2-PI**2)/4.0 PQ2 = PI/2.0*LQ AQ1 = -TQ/2.0D0*(1.0D0+(1.0D0-TQ)*PQ1) AQ2 = -TQ/2.0D0*(1.0D0-TQ)*PQ2 AQ = AQ1**2+AQ2**2 ELSE AP = (ASIN(1.0D0/SQRT(TQ)))**2 AQ = 9.0D0*(0.5D0*TQ*(1.0D0+(1.0D0-TQ)*AP))**2 ENDIF CCC SIGIN = GF/288.0D0/SQRT(2.0D0)/PI*ALPHAS2(WHAT)**2*SHAT**2*AQ SIGIN = GF/288.0D0/SQRT(2.0D0)/PI*PYALPS(SHAT)**2*SHAT**2*AQ C...Boson Propagator C BWIG = WHAT/((SHAT-AMH**2)**2+AGH**2*AMH**2) C...For s-dependence width. BWIG = WHAT/((SHAT-AMH**2)**2+AGH**2*SHAT) C...Calculate Gamma(H -> bb). MBRUN = PYMRUN(5,SHAT) ! MSbar b-mass. CCC RDCOR = 1.0D0 + ALPHAS2(WHAT)/PI ! Radiative correction RDCOR = 1.0D0 + PYALPS(SHAT)/PI SIGOUT = 3.0D0*GF/4.0D0/SQRT(2.0D0)/PI*MBRUN**2*WHAT*RDCOR C...Differencial cross section XWGTUP = AJACOB*SIGIN*BWIG*SIGOUT C...Flavour codes for entries. NUP = 4 ! Nr. of particles IDUP(1) = 21 ! gluon IDUP(2) = 21 ! gluon IDUP(3) = 5 ! b IDUP(4) = -5 ! b C...Status codes for entries. ISTUP(1) = -1 ISTUP(2) = -1 ISTUP(3) = 1 ISTUP(4) = 1 C...Position of the first and last mother particles for entries. MOTHUP(1,3) = 1 MOTHUP(2,3) = 2 MOTHUP(1,4) = 1 MOTHUP(2,4) = 2 C...Possible color strings for entries. ICOLUP(1,1) = 501 ICOLUP(2,2) = 501 ICOLUP(1,2) = 502 ICOLUP(2,1) = 502 ICOLUP(1,3) = 503 ICOLUP(2,4) = 503 C...Particle 1 PUP(1,1) = 0.0D0 PUP(2,1) = 0.0D0 PUP(3,1) = 0.5D0*X1*WCM PUP(4,1) = 0.5D0*X1*WCM PUP(5,1) = 0.0D0 C...Particle 2 PUP(1,2) = 0.0D0 PUP(2,2) = 0.0D0 PUP(3,2) = -0.5D0*X2*WCM PUP(4,2) = 0.5D0*X2*WCM PUP(5,2) = 0.0D0 C...Particle 3 E3 = WHAT/2.D0 P3 = SQRT((E3-MB)*(E3+MB)) SIN3 = SQRT((1.0D0 - COS3)*(1.0D0 + COS3)) PUP(1,3) = P3*SIN3*COS(PHI3) PUP(2,3) = P3*SIN3*SIN(PHI3) PUP(3,3) = GAMM*(P3*COS3 + BETA*E3) PUP(4,3) = GAMM*(E3 + BETA*P3*COS3) PUP(5,3) = MB C...Particle 4 PUP(1,4) = -PUP(1,3) PUP(2,4) = -PUP(2,3) PUP(3,4) = GAMM*(-P3*COS3 + BETA*E3) PUP(4,4) = GAMM*(E3 - BETA*P3*COS3) PUP(5,4) = MB C...Set Q^2 and couplings. (whatever...) SCALUP = WHAT AQEDUP = 1.0/127.0 AQCDUP = ALPHAS2(WHAT) C...Ensure proper normalization of weights for MODE=+-2 IDPRUP = 200 C...Done. RETURN END