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 C...Local arrays. REAL*8 XLUM character*72 FNAME character*5 cgive character*30 cgive0 integer istr integer lok,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. OPEN(55,FILE='pyxtra.in',status='old') NEV=100 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) C...Initialize with external process. CALL PYINIT('USER',' ',' ',0D0) NEV = INT (XLUM * XSECUP(1)) IF(NEV.LE.0) NEV=1 write(*,*) ' requesting ',XLUM,' pbi of data ' write(*,*) ' cross section is ',xsecup(1),' pb ' write(*,*) ' events to be generated = ',nev C.....Opening stdhep file for writing call stdxwinit(FNAME,'StdHep/Pythia example', 1 nev,istr,lok) if(lok.ne.0) write(lnhout,*) 1 ' Problem opening file ' C Write Stdhep begin-run record call stdxwrt(100,istr,lok) if(lok.ne.0) write(lnhout,*) 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 CALL LUNHEP(1) IF(MSTI(51).EQ.1) GOTO 140 IF(IEV.LT.50) THEN IF(MOD(IEV-1,10).EQ.0) THEN PRINT*,' IEV = ',IEV CALL PYLIST(7) CALL PYLIST(2) ENDIF ENDIF C.......Write one event call stdxwrt(1,istr,lok) 130 CONTINUE C...Statistics and histograms. 140 CALL PYSTAT(1) C Fill Stdhep common block 1 with run information call stdflpyxsec(nev) C Write end-of-run record call stdxwrt(200,istr,lok) if(lok.ne.0) write(lnhout,*) ' Problem writing end run record' c...close event file 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*80 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,WTGRT INTEGER I,KNTEV,J,NPARTS,IMAX,NGRT C...Extra commonblock to transfer run info. INTEGER MODE,NLIM REAL*8 SCALEF,WTMAX COMMON/PRIV/MODE,NLIM,SCALEF,WTMAX SAVE/PRIV/ 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 WTGRT=0D0 NGRT=0D0 IMAX=0 DO I=1,10000000 READ(77,*,END=777) NPARTS,KNTEV,WTEVNT IF(WTEVNT.GT.WTMAX) WTMAX=WTEVNT XSTOT=XSTOT+WTEVNT DO J=1,5 READ(77,*) CHAR_READ ENDDO DO J=1,NPARTS READ(77,*) CHAR_READ ENDDO IMAX=I ENDDO 777 CONTINUE REWIND(77) IF(SCALEF.LT.1D0.AND.SCALEF.GT.0D0) THEN WTMAX=WTMAX*SCALEF DO I=1,IMAX READ(77,*,END=778) NPARTS,KNTEV,WTEVNT IF(WTEVNT.GT.WTMAX) THEN NGRT=NGRT+1 WTGRT=WTGRT+WTEVNT ENDIF DO J=1,5 READ(77,*) CHAR_READ ENDDO DO J=1,NPARTS READ(77,*) CHAR_READ ENDDO ENDDO 778 CONTINUE REWIND(77) WRITE(*,*) ' For unweighting, scaled maximum by ',SCALEF WRITE(*,*) ' Number of Events with wt > wtmax ',NGRT, $ ' out of ',IMAX WRITE(*,*) ' Total Cross Section with wt > wtmax ',WTGRT WRITE(*,*) ' Relative to total ',WTGRT/XSTOT ENDIF 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,ID,IDUMB,IC,J C...If PYTHIA is supposed to select event type, do not modify this choice. IF(IABS(MODE).LE.2) THEN IDPRUP=661 C...Else free hands to mix; here evenly. ELSE ENDIF 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 READ(77,*,END=888) NUP,ID,XWGTUP C.....Account for scale factor IF(XWGTUP.GT.WTMAX) XWGTUP=WTMAX READ(77,*,END=888) (IDUP(I),I=1,NUP) ISTUP(1)=-1 ISTUP(2)=-1 DO I=3,NUP ISTUP(I)=1 ENDDO READ(77,*,END=888) (MOTHUP(1,I),I=1,NUP) READ(77,*,END=888) (MOTHUP(2,I),I=1,NUP) READ(77,*,END=888) (ICOLUP(1,I),I=1,NUP) READ(77,*,END=888) (ICOLUP(2,I),I=1,NUP) 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...Some other compulsory quantities. SCALUP=-1D0 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)) C DO IJ=1,NUP IDC1=ICOLUP(1,IJ) IDC2=ICOLUP(2,IJ) IF(IDC1.EQ.0) IDC1=-1 IF(IDC2.EQ.0) IDC2=-2 DO IC=IJ+1,NUP IC1=ICOLUP(1,IC) IC2=ICOLUP(2,IC) IF(ISTUP(IC)*ISTUP(IJ).EQ.1) THEN IF(IDC1.EQ.IC2.OR.IDC2.EQ.IC1) THEN QTMP=SMDOT5(PUP(1,IJ),PUP(1,IC)) QMIN=MIN(QMIN,QTMP) c print*,ij,ic,qmin,idc1,ic2,idc2,ic1 ENDIF ELSEIF(ISTUP(IC)*ISTUP(IJ).EQ.-1) THEN IF(IDC1.EQ.IC1.OR.IDC2.EQ.IC2) THEN QTMP=SMDOT5(PUP(1,IJ),PUP(1,IC)) QMIN=MIN(QMIN,QTMP) c print*,ij,ic,qmin,idc1,ic1,idc2,ic2 ENDIF ENDIF ENDDO ENDDO SCALUP=QMIN c print*,' scalup = ',scalup RETURN 888 CONTINUE CC REWIND(77) CC GOTO 99 NUP=0 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(2D0*ABS(TEMP)) RETURN END