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.....V3.1 ... intended to read CompHEP data 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=10 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.10) THEN CALL PYLIST(7) CALL PYLIST(2) ELSEIF(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,nproc 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,1000 READ(78,'(A80)') CHAR_READ WRITE(*,*) CHAR_READ(1:26) IF (CHAR_READ(1:26).EQ.'#Number_of_subprocesses = ') THEN READ(CHAR_READ,FMT='(26x,i4)') nproc print*,' nproc = ',nproc GOTO 1 ENDIF ENDDO 1 REWIND(78) READ(78,*) CHAR_READ DO I=1,NPROC DO J=1,11 READ(78,*) CHAR_READ ENDDO ENDDO DO J=1,4 READ(78,*) CHAR_READ ENDDO 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/ character*60 chain REAL*8 TEST,QMIN,SMDOT5,QTMP INTEGER I,IJ,IC1,IC2,IDC1,IDC2,IDUMB,IC,J INTEGER MUP,IDH,IDL,ID,IDP,IDA,IDPA,IMO,IS 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 NUP=4 cc 1001 FORMAT(i6,8(G18.10),G12.3,1x,A) c READ(78,'(A80)') chain c print*,' chain ',chain READ(78,*,END=888) ID,PUP(3,1),PUP(3,2), $((PUP(I,J),I=1,3),J=3,NUP),SCALUP,chain ISTUP(1)=-1 ISTUP(2)=-1 ISTUP(3)=1 ISTUP(4)=1 IDUP(3)=11 IDUP(4)=-11 IF(ID.EQ.1) THEN IDUP(1)=2 ELSEIF(ID.EQ.2) THEN IDUP(1)=-2 ELSEIF(ID.EQ.3) THEN IDUP(1)=1 ELSE IDUP(1)=-1 ENDIF IDUP(2)=-IDUP(1) PUP(4,1)=ABS(PUP(3,1)) PUP(4,2)=ABS(PUP(3,2)) if(idup(1).lt.0) then icolup(2,1)=501 icolup(1,2)=501 else icolup(1,1)=501 icolup(2,2)=501 endif DO IJ=3,NUP TEST=PUP(1,IJ)**2+PUP(2,IJ)**2+PUP(3,IJ)**2 PUP(4,IJ)=SQRT(TEST) 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)) 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) 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) ENDIF ENDIF ENDDO ENDDO SCALUP=QMIN c print*,' scalup = ',scalup C.....Add mothers for l+ l- and l nu_l pairs IJ=3 555 CONTINUE IS=ISTUP(IJ) IF(IS.NE.1) GOTO 998 ID=IDUP(IJ) IDA=ABS(ID) 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 MUP=NUP 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 MOTHUP(2,IJ)=2 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 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