PROGRAM PYXTRA C...Long example how to interface user-defined processes to PYTHIA C...based on the Les Houches commonblock agreement. Generates events C...of several different kinds, to test the new code under varied C...conditions, but kinematics selection and cross sections are C...completely unphysical. C...Double precision and integer declarations. IMPLICIT NONE INTEGER PYK,PYCHGE,PYCOMP C...User process 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...User process initialization commonblock. 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...Standard PYTHIA commonblocks. INTEGER N,NPAD,K REAL*8 P,V COMMON/PYJETS/N,NPAD,K(4000,5),P(4000,5),V(4000,5) INTEGER MSTU,MSTJ REAL*8 PARU,PARJ COMMON/PYDAT1/MSTU(200),PARU(200),MSTJ(200),PARJ(200) INTEGER MDCY,MDME,KFDP REAL*8 BRAT COMMON/PYDAT3/MDCY(500,3),MDME(8000,2),BRAT(8000),KFDP(8000,5) INTEGER MSEL,MSELPD,MSUB,KFIN REAL*8 CKIN COMMON/PYSUBS/MSEL,MSELPD,MSUB(500),KFIN(2,-40:40),CKIN(200) INTEGER MSTP,MSTI REAL*8 PARP,PARI COMMON/PYPARS/MSTP(200),PARP(200),MSTI(200),PARI(200) C...Extra commonblock to transfer run info. INTEGER MODE,NLIM REAL*8 SCALEF,WTMAX COMMON/PRIV/MODE,NLIM,SCALEF,WTMAX INTEGER IREAD COMMON/INFO/IREAD C...Local arrays. REAL*8 XLUM character*72 FNAME character*5 cgive character*30 cgive0 integer NEV,iev INTEGER LNHWRT,LNHRD,LNHDCY,LNHOUT EXTERNAL PYDATA C...Switch process mode; agrees with IDWTUP code (+-1,+-2,+-3,+-4). MODE=2 C.....2 means weighted in and unweighted out C......Events selected according to xsection CCC MODE=1 C.....1 means weighted in and unweighted out C......Events selected according to maximum weight C...Maximum number of events to generate. C/* OPEN(55,FILE='pyxtra.in',status='old') C NEV=100 C READ(55,*) MODE,XLUM,FNAME,SCALEF C initialize HEP logical units lnhwrt=23 lnhrd=0 lnhdcy=0 lnhout=22 C WRITE(CGIVE,'(I5)') lnhout CGIVE0='MSTU(11)='//CGIVE CALL PYGIVE(CGIVE0) CALL PYGIVE('MDCY(C111,1)=0') CALL PYGIVE('MDCY(C211,1)=0') C...Initialize with external process. CALL PYINIT('USER',' ',' ',0D0) NEV = IREAD IF(NEV.LE.0) NEV=1 XLUM=NEV/XSECUP(1) write(lnhout,*) ' requesting ',XLUM,' pbi of data ' write(lnhout,*) ' cross section is ',xsecup(1),' pb ' write(lnhout,*) ' events to be generated = ',nev C.....Opening stdhep file for writing cc call stdxwinit(FNAME,'StdHep/Pythia example', cc 1 nev,istr,lok) cc if(lok.ne.0) write(lnhout,*) cc 1 ' Problem opening file ' C Write Stdhep begin-run record cc call stdxwrt(100,istr,lok) cc if(lok.ne.0) write(lnhout,*) cc 1 ' Problem writing stdhep begin run record' C...Event loop; generate event; check it was obtained or quit. DO 130 IEV=1,NEV CALL PYEVNT C/* CALL LUNHEP(1) IF(MSTI(51).EQ.1) THEN WRITE(lnhout,*) 'MSTI(51) = ',MSTI(51) GOTO 140 ENDIF IF(IEV.LT.10) THEN CALL PYLIST(7) CALL PYLIST(2) ELSEIF(IEV.LT.50) THEN IF(MOD(IEV-1,10).EQ.0) THEN WRITE(lnhout,*) ' IEV = ',IEV CALL PYLIST(7) CALL PYLIST(2) ENDIF ENDIF C.......Write one event cc call stdxwrt(1,istr,lok) 130 CONTINUE C...Statistics and histograms. 140 CALL PYSTAT(1) C Fill Stdhep common block 1 with run information cc call stdflpyxsec(nev) C Write end-of-run record cc call stdxwrt(200,istr,lok) cc if(lok.ne.0) write(lnhout,*) ' Problem writing end run record' c...close event file cc call stdxend(istr) END C********************************************************************* C...UPINIT C...Routine to be called by user to set up user-defined processes. C...Code below only intended as example, without any claim of realism. C...Especially it shows what info needs to be put in HEPRUP. SUBROUTINE UPINIT C...Double precision and integer declarations. IMPLICIT NONE CHARACTER*132 CHAR_READ C...User process initialization commonblock. 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/ REAL*8 XSTOT,WTEVNT INTEGER I,KNTEV,J,NPARTS,IMAX,IBEAM,ICOM1,IPND1,LEN character*5 cgive character*10 cgive0 C...Extra commonblock to transfer run info. INTEGER MODE,NLIM REAL*8 SCALEF,WTMAX COMMON/PRIV/MODE,NLIM,SCALEF,WTMAX SAVE/PRIV/ INTEGER IREAD COMMON/INFO/IREAD C...Set incoming beams: Tevatron Run II. IDBMUP(1)=2212 IDBMUP(2)=-2212 EBMUP(1)=980D0 EBMUP(2)=980D0 C...Set PDF's of incoming beams: CTEQ 5L. C...Note that Pythia will not look at PDFGUP and PDFSUP. PDFGUP(1)=4 PDFSUP(1)=46 PDFGUP(2)=PDFGUP(1) PDFSUP(2)=PDFSUP(1) C...If you want Pythia to use PDFLIB, you have to set it by hand. C...(You also have to ensure that the dummy routines C...PDFSET, STRUCTM and STRUCTP in Pythia are not linked.) C MSTP(52)=2 C MSTP(51)=1000*PDFGUP(1)+PDFSUP(1) C...Decide on weighting strategy: unweighted on input. IDWTUP=MODE C...Number of external processes. NPRUP=1 C.....read through data file to find maximum weight....only first time WTMAX=0D0 XSTOT=0D0 IMAX=0 IREAD=0 IBEAM=0 DO I=1,10000000 C-- ignore comment lines READ(77,'(A132)',END=777,ERR=777) CHAR_READ DO WHILE(INDEX(CHAR_READ,"#") .NE. 0) READ(77,'(A132)',END=777,ERR=777) CHAR_READ IPND1=INDEX(CHAR_READ,"#") IF(INDEX(CHAR_READ,"= lpp(1)").NE.0) THEN ICOM1=INDEX(CHAR_READ,"= lpp(1)") LEN=ICOM1-IPND1-1 WRITE(CGIVE,'(I5)') LEN CGIVE0='(I'//CGIVE(1:lnblnk(CGIVE))//')' READ(CHAR_READ(IPND1+1:ICOM1-1),CGIVE0) IBEAM IF(IBEAM.EQ.1) THEN IDBMUP(1)=2212 ELSEIF(IBEAM.EQ.-1) THEN IDBMUP(1)=-2212 ELSE IDBMUP(1)=IBEAM ENDIF ELSEIF(INDEX(CHAR_READ,"= lpp(2)").NE.0) THEN ICOM1=INDEX(CHAR_READ,"= lpp(2)") LEN=ICOM1-IPND1-1 WRITE(CGIVE,'(I5)') LEN CGIVE0='(I'//CGIVE(1:lnblnk(CGIVE))//')' READ(CHAR_READ(IPND1+1:ICOM1-1),CGIVE0) IBEAM IF(IBEAM.EQ.1) THEN IDBMUP(2)=2212 ELSEIF(IBEAM.EQ.-1) THEN IDBMUP(2)=-2212 ELSE IDBMUP(2)=IBEAM ENDIF ELSEIF(INDEX(CHAR_READ,"= ebeam(1)").NE.0) THEN ICOM1=INDEX(CHAR_READ,"= ebeam(1)") LEN=ICOM1-IPND1-1 WRITE(CGIVE,'(I5)') LEN CGIVE0='(F'//CGIVE(1:lnblnk(CGIVE))//')' READ(CHAR_READ(IPND1+1:ICOM1-1),CGIVE0) EBMUP(1) ELSEIF(INDEX(CHAR_READ,"= ebeam(2)").NE.0) THEN ICOM1=INDEX(CHAR_READ,"= ebeam(2)") LEN=ICOM1-IPND1-1 WRITE(CGIVE,'(I5)') LEN CGIVE0='(F'//CGIVE(1:lnblnk(CGIVE))//')' READ(CHAR_READ(IPND1+1:ICOM1-1),CGIVE0) EBMUP(2) ELSEIF(INDEX(CHAR_READ,"= tmass").NE.0) THEN ICOM1=INDEX(CHAR_READ,"= tmass") call pygive("pmas(6,1)="//CHAR_READ(IPND1+1:ICOM1-1)) ELSEIF(INDEX(CHAR_READ,"= bmass").NE.0) THEN ICOM1=INDEX(CHAR_READ,"= bmass") call pygive("pmas(5,1)="//CHAR_READ(IPND1+1:ICOM1-1)) ELSEIF(INDEX(CHAR_READ,"= cmass").NE.0) THEN ICOM1=INDEX(CHAR_READ,"= cmass") call pygive("pmas(4,1)="//CHAR_READ(IPND1+1:ICOM1-1)) ELSEIF(INDEX(CHAR_READ,"= lmass").NE.0) THEN ICOM1=INDEX(CHAR_READ,"= lmass") call pygive("pmas(15,1)="//CHAR_READ(IPND1+1:ICOM1-1)) ELSEIF(INDEX(CHAR_READ,"= hmass").NE.0) THEN ICOM1=INDEX(CHAR_READ,"= hmass") call pygive("pmas(25,1)="//CHAR_READ(IPND1+1:ICOM1-1)) ENDIF ENDDO C-- read number of particles and weight. Ignore the rest. READ(CHAR_READ,*,END=777) NPARTS,KNTEV,WTEVNT IF(WTEVNT.GT.WTMAX) WTMAX=WTEVNT XSTOT=XSTOT+WTEVNT IREAD=IREAD+1 c-- skip the rest of the event record DO J=1,7 !idup,mothup(2),icolup(2),istup,spinup READ(77,*) CHAR_READ ENDDO DO J=1,NPARTS !momenta READ(77,*) CHAR_READ ENDDO IMAX=I ENDDO 777 CONTINUE REWIND(77) XSTOT=XSTOT*1D3 WTMAX=WTMAX*1D3 WRITE(*,*) ' Number of Events Read:: ',Iread WRITE(*,*) ' Total Cross Section :: ',XSTOT WRITE(*,*) ' Maximum Weight :: ',WTMAX XSECUP(1)=XSTOT XMAXUP(1)=WTMAX LPRUP(1)=661 RETURN END C********************************************************************* C...UPEVNT C...Sample routine to generate events of various kinds. C...Not intended to be realistic, but rather to show in closed C...and understandable form what such a routine might look like. C...Especially it shows what info needs to be put in HEPEUP. SUBROUTINE UPEVNT C...Double precision and integer declarations. IMPLICIT NONE INTEGER MSTP,MSTI REAL*8 PARP,PARI COMMON/PYPARS/MSTP(200),PARP(200),MSTI(200),PARI(200) C...Extra commonblock to transfer run info. INTEGER MODE,NLIM REAL*8 SCALEF,WTMAX COMMON/PRIV/MODE,NLIM,SCALEF,WTMAX SAVE/PRIV/ C...User process 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) SAVE /HEPEUP/ REAL*8 TEST,QMIN,SMDOT5,QTMP INTEGER I,IJ,IC1,IC2,IDC1,IDC2,IC,J,IMA,ICC1,ICC2 INTEGER IDH,IDL,ID,IDP,IDA,IDPA,IMO,IS,IDUMB CHARACTER*(132) BUFF C...If PYTHIA is supposed to select event type, do not modify this choice. IDPRUP=661 99 CONTINUE C...Zero some arrays in common blocks to simplify filling. NUP=12 DO 100 I=1,NUP MOTHUP(1,I)=0 MOTHUP(2,I)=0 ICOLUP(1,I)=0 ICOLUP(2,I)=0 SPINUP(I)=0D0 PUP(1,I)=0D0 PUP(2,I)=0D0 PUP(5,I)=0D0 VTIMUP(I)=0D0 100 CONTINUE C--ignore comment lines READ(77,'(A132)',END=888,ERR=888) BUFF DO WHILE(INDEX(BUFF,"#") .NE. 0) READ(77,'(A132)',END=888,ERR=888) BUFF ENDDO C--read number of particles,id#,weight,scale,alpha_em,alpha_s READ(BUFF,*,END=888) NUP,ID,XWGTUP,SCALUP,AQEDUP,AQCDUP C.....Account for scale factor XWGTUP=XWGTUP*1D3 IF(XWGTUP.GT.WTMAX) XWGTUP=WTMAX C--read the particle id's READ(77,*,END=888) (IDUP(I),I=1,NUP) C--read the mothers READ(77,*,END=888) (MOTHUP(1,I),I=1,NUP) READ(77,*,END=888) (MOTHUP(2,I),I=1,NUP) C--read the colours READ(77,*,END=888) (ICOLUP(1,I),I=1,NUP) READ(77,*,END=888) (ICOLUP(2,I),I=1,NUP) c--read the status code READ(77,*,END=888,ERR=888) (ISTUP (I),I=1,NUP) c--read the spin READ(77,*,END=888,ERR=888) (SPINUP(I),I=1,NUP) C--read in the momenta DO I=1,NUP READ(77,*,END=888) IDUMB,PUP(4,I),(PUP(J,I),J=1,3) TEST=PUP(1,I)**2+PUP(2,I)**2+PUP(3,I)**2-PUP(4,I)**2 IF(TEST.LE.0D0) THEN PUP(5,I)=DSQRT(-TEST) ELSEIF(DABS(TEST).LT.1D-4) THEN PUP(5,I)=0D0 ELSE PUP(5,I)=-1D0 PRINT*,' NEGATIVE MASS ' ENDIF ENDDO C...Guesses for the correct scale C Assumptions: C (1) if the initial state is a color singlet, then C use s-hat for the scale C C (2) if color flow to the final state, use the minimum C of the dot products of color connected pairs C (times two for consistency with above) QMIN=SMDOT5(PUP(1,1),PUP(1,2)) ICC1=1 ICC2=2 C/* C/* For now, there is no generic way to guarantee the "right" C/* scale choice. Here, we take the HERWIG pt. of view and C/* choose the dot product of the colored connected "primary" C/* pairs. C/* DO 101 IJ=1,NUP IF(MOTHUP(2,IJ).GT.2) GOTO 101 IDC1=ICOLUP(1,IJ) IDC2=ICOLUP(2,IJ) IF(IDC1.EQ.0) IDC1=-1 IF(IDC2.EQ.0) IDC2=-2 DO 201 IC=IJ+1,NUP IF(MOTHUP(2,IC).GT.2) GOTO 201 IC1=ICOLUP(1,IC) IC2=ICOLUP(2,IC) IF(ISTUP(IC)*ISTUP(IJ).GE.1) THEN IF(IDC1.EQ.IC2.OR.IDC2.EQ.IC1) THEN QTMP=SMDOT5(PUP(1,IJ),PUP(1,IC)) IF(QTMP.LT.QMIN) THEN QMIN=QTMP ICC1=IJ ICC2=IC ENDIF ENDIF ELSEIF(ISTUP(IC)*ISTUP(IJ).LE.-1) THEN IF(IDC1.EQ.IC1.OR.IDC2.EQ.IC2) THEN QTMP=SMDOT5(PUP(1,IJ),PUP(1,IC)) IF(QTMP.LT.QMIN) THEN QMIN=QTMP ICC1=IJ ICC2=IC ENDIF ENDIF ENDIF 201 CONTINUE 101 CONTINUE SCALUP=QMIN C/* Add mothers for l+ l- and l nu_l pairs C/* but only if they don't already exist IJ=3 555 CONTINUE IS=ISTUP(IJ) IF(IS.NE.1) GOTO 998 ID=IDUP(IJ) IDA=ABS(ID) IMA=MOTHUP(1,IJ) IF(IDA.GE.11.AND.IDA.LE.16.AND.IJ.LT.NUP) THEN IF(ICOLUP(1,IJ).EQ.0 .AND. ICOLUP(2,IJ).EQ.0) THEN IDP=IDUP(IJ+1) IDPA=ABS(IDP) IMO=0 IF(IDP.EQ.-ID) THEN C.......gamma*/Z* mother IMO=22 ELSE IDL=MIN(IDA,IDPA) IDH=MAX(IDA,IDPA) IF(MOD(IDH,2).EQ.0.AND.(IDH-IDL.EQ.1)) THEN IF(ABS(ID).EQ.IDL) THEN IMO=-SIGN(24,ID) ELSE IMO=-SIGN(24,IDP) ENDIF ENDIF ENDIF IF(IMO.EQ.0) GOTO 998 C/* If motherhood already treated correctly IF(IMA.NE.0.AND.IDUP(IMA).GT.21) GOTO 998 NUP=NUP+1 DO IC=NUP,IJ+1,-1 IDUP(IC)=IDUP(IC-1) ISTUP(IC)=ISTUP(IC-1) MOTHUP(1,IC)=MOTHUP(1,IC-1) MOTHUP(2,IC)=MOTHUP(2,IC-1) ICOLUP(1,IC)=ICOLUP(1,IC-1) ICOLUP(2,IC)=ICOLUP(2,IC-1) DO J=1,5 PUP(J,IC)=PUP(J,IC-1) ENDDO ENDDO IDUP(IJ)=IMO ISTUP(IJ)=2 DO J=1,4 PUP(J,IJ)=PUP(J,IJ+1)+PUP(J,IJ+2) ENDDO MOTHUP(1,IJ+1)=IJ MOTHUP(2,IJ+1)=IJ MOTHUP(1,IJ+2)=IJ MOTHUP(2,IJ+2)=IJ TEST=PUP(1,IJ)**2+PUP(2,IJ)**2+PUP(3,IJ)**2-PUP(4,IJ)**2 IF(TEST.LE.0D0) THEN PUP(5,IJ)=DSQRT(-TEST) ELSEIF(DABS(TEST).LT.1D-4) THEN PUP(5,IJ)=0D0 ELSE PUP(5,IJ)=-1D0 PRINT*,' NEGATIVE MASS ' ENDIF IJ=IJ+2 ENDIF ENDIF 998 CONTINUE IJ=IJ+1 IF(IJ.GT.NUP) GOTO 999 GOTO 555 999 CONTINUE RETURN 888 CONTINUE NUP=0 WRITE(*,*) 'last event !' RETURN END C FUNCTION SMDOT5(V1,V2) IMPLICIT NONE REAL*8 SMDOT5,TEMP REAL*8 V1(5),V2(5) INTEGER I SMDOT5=0D0 TEMP=V1(4)*V2(4) DO I=1,3 TEMP=TEMP-V1(I)*V2(I) ENDDO SMDOT5=SQRT(ABS(TEMP)) RETURN END