pythiaeRHIC/pythia/pyveto.f

0001
0002 C*********************************************************************
0003
0004 C...PYVETO
0005 C...Interface to UPVETO, which allows user to veto event generation
0006 C...on the parton level, after parton showers but before multiple
0007 C...interactions, beam remnants and hadronization is added.
0008
0009       SUBROUTINE PYVETO(IVETO)
0010
0011 C...All real arithmetic in double precision.
0012       IMPLICIT DOUBLE PRECISION(A-H, O-Z)
0013 C...Three Pythia functions return integers, so need declaring.
0014       INTEGER PYK,PYCHGE,PYCOMP
0015
0016 C...PYTHIA commonblocks.
0017       COMMON/PYJETS/N,NPAD,K(4000,5),P(4000,5),V(4000,5)
0018       COMMON/PYPARS/MSTP(200),PARP(200),MSTI(200),PARI(200)
0019       COMMON/PYDAT1/MSTU(200),PARU(200),MSTJ(200),PARJ(200)
0020       COMMON/PYINT1/MINT(400),VINT(400)
0021       SAVE /PYJETS/,/PYPARS/,/PYINT1/
0022 C...HEPEVT commonblock.
0023       PARAMETER (NMXHEP=4000)
0024       COMMON/HEPEVT/NEVHEP,NHEP,ISTHEP(NMXHEP),IDHEP(NMXHEP),
0025      &JMOHEP(2,NMXHEP),JDAHEP(2,NMXHEP),PHEP(5,NMXHEP),VHEP(4,NMXHEP)
0026       DOUBLE PRECISION PHEP,VHEP
0027       SAVE /HEPEVT/
0028 C...Local array.
0029       DIMENSION IRESO(100)
0030
0031 C...Define longitudinal boost from initiator rest frame to cm frame.
0032       GAMMA=0.5D0*(VINT(141)+VINT(142))/SQRT(VINT(141)*VINT(142))
0033       GABEZ=0.5D0*(VINT(141)-VINT(142))/SQRT(VINT(141)*VINT(142))
0034
0035 C... Reset counters.
0036       NEVHEP=0
0037       NHEP=0
0038       NRESO=0
0039
0040 C...Oth pass: identify beam and incoming partons
0041       DO 140 I=MINT(83)+1,MINT(83)+6
0042         ISTORE=0
0043         IF(K(I,2).EQ.94) THEN
0044
0045         ELSE
0046           ISTORE=1
0047           NHEP=NHEP+1
0048           II=NHEP
0049           NRESO=NRESO+1
0050           IRESO(NRESO)=I
0051           IMOTH=K(I,3)
0052         ENDIF
0053         IF(ISTORE.EQ.1) THEN
0054 C...Copy parton info, boosting momenta along z axis to cm frame.
0055           ISTHEP(II)=2
0056           IDHEP(II)=K(I,2)
0057           PHEP(1,II)=P(I,1)
0058           PHEP(2,II)=P(I,2)
0059           PHEP(3,II)=GAMMA*P(I,3)+GABEZ*P(I,4)
0060           PHEP(4,II)=GAMMA*P(I,4)+GABEZ*P(I,3)
0061           PHEP(5,II)=P(I,5)
0062 C...Store one mother. Rest of history and vertex info zeroed.
0063           JMOHEP(1,II)=IMOTH
0064           JMOHEP(2,II)=0
0065           JDAHEP(1,II)=0
0066           JDAHEP(2,II)=0
0067           VHEP(1,II)=0D0
0068           VHEP(2,II)=0D0
0069           VHEP(3,II)=0D0
0070           VHEP(4,II)=0D0
0071         ENDIF
0072  140  CONTINUE
0073
0074 C...First pass: identify final locations of resonances
0075 C...and of their daughters before showering.
0076       DO 150 I=MINT(84)+3,N
0077         ISTORE=0
0078         IMOTH=0
0079
0080 C...Skip shower CM frame documentation lines.
0081         IF(K(I,2).EQ.94) THEN
0082
0083 C...  Store a new intermediate product, when mother in documentation.
0084         ELSEIF(MSTP(128).EQ.0.AND.K(I,3).GT.MINT(83)+6.AND.
0085      &  K(I,3).LE.MINT(84)) THEN
0086           ISTORE=1
0087           NHEP=NHEP+1
0088           II=NHEP
0089           NRESO=NRESO+1
0090           IRESO(NRESO)=I
0091           IMOTH=K(K(I,3),3)
0092
0093 C...  Store a new intermediate product, when mother in main section.
0094         ELSEIF(MSTP(128).EQ.1.AND.K(I-MINT(84)+MINT(83)+4,1).EQ.21.AND.
0095      &  K(I-MINT(84)+MINT(83)+4,2).EQ.K(I,2)) THEN
0096           ISTORE=1
0097           NHEP=NHEP+1
0098           II=NHEP
0099           NRESO=NRESO+1
0100           IRESO(NRESO)=I
0101           IMOTH=MAX(0,K(I-MINT(84)+MINT(83)+4,3))
0102         ENDIF
0103
0104         IF(ISTORE.EQ.1) THEN
0105 C...Copy parton info, boosting momenta along z axis to cm frame.
0106           ISTHEP(II)=2
0107           IDHEP(II)=K(I,2)
0108           PHEP(1,II)=P(I,1)
0109           PHEP(2,II)=P(I,2)
0110           PHEP(3,II)=GAMMA*P(I,3)+GABEZ*P(I,4)
0111           PHEP(4,II)=GAMMA*P(I,4)+GABEZ*P(I,3)
0112           PHEP(5,II)=P(I,5)
0113 C...Store one mother. Rest of history and vertex info zeroed.
0114           JMOHEP(1,II)=IMOTH
0115           JMOHEP(2,II)=0
0116           JDAHEP(1,II)=I
0117           JDAHEP(2,II)=0
0118           VHEP(1,II)=0D0
0119           VHEP(2,II)=0D0
0120           VHEP(3,II)=0D0
0121           VHEP(4,II)=0D0
0122         ENDIF
0123  150  CONTINUE
0124
0125 C...Second pass: identify current set of "final" partons.
0126       DO 200 I=MINT(84)+3,N
0127         ISTORE=0
0128         IMOTH=0
0129
0130 C...Store a final parton.
0131         IF(K(I,1).GE.1.AND.K(I,1).LE.10) THEN
0132           ISTORE=1
0133           NHEP=NHEP+1
0134           II=NHEP
0135 C..Trace it back through shower, to check if from documented particle.
0136           IHIST=I
0137           ISAVE=IHIST
0138   160     CONTINUE
0139           IF(IHIST.GT.MINT(84)) THEN
0140             IF(K(IHIST,2).EQ.94) IHIST=K(IHIST,3)+(ISAVE-1-IHIST)
0141             DO 170 IRI=1,NRESO
0142               IF(IHIST.EQ.IRESO(IRI)) IMOTH=IRI
0143   170       CONTINUE
0144             ISAVE=IHIST
0145             IHIST=K(IHIST,3)
0146             IF(IMOTH.EQ.0) GOTO 160
0147           ELSEIF(IHIST.LE.4) THEN
0148             IF(IHIST.EQ.1.OR.IHIST.EQ.2) THEN
0149               ISTORE=0
0150               NHEP=NHEP-1
0151             ELSE
0152               IMOTH=IHIST
0153             ENDIF
0154           ENDIF
0155         ENDIF
0156
0157         IF(ISTORE.EQ.1) THEN
0158 C...Copy parton info, boosting momenta along z axis to cm frame.
0159           ISTHEP(II)=1
0160           IDHEP(II)=K(I,2)
0161           PHEP(1,II)=P(I,1)
0162           PHEP(2,II)=P(I,2)
0163           PHEP(3,II)=GAMMA*P(I,3)+GABEZ*P(I,4)
0164           PHEP(4,II)=GAMMA*P(I,4)+GABEZ*P(I,3)
0165           PHEP(5,II)=P(I,5)
0166 C...Store one mother. Rest of history and vertex info zeroed.
0167           JMOHEP(1,II)=IMOTH
0168           JMOHEP(2,II)=0
0169           JDAHEP(1,II)=0
0170           JDAHEP(2,II)=0
0171           VHEP(1,II)=0D0
0172           VHEP(2,II)=0D0
0173           VHEP(3,II)=0D0
0174           VHEP(4,II)=0D0
0175         ENDIF
0176   200 CONTINUE
0177
0178 C...Call user-written routine to decide whether to keep events.
0179       CALL UPVETO(IVETO)
0180
0181       RETURN
0182       END