Back to home page

sPhenix code displayed by LXR
 
 

    

File indexing completed on 2025-08-05 08:21:08 

0001  
0002 C*********************************************************************
0003  
0004 C...PY4JET
0005 C...An interface from a four-parton generator to include
0006 C...parton showers and hadronization.
0007  
0008       SUBROUTINE PY4JET(PMAX,IRAD,ICOM)
0009  
0010 C...Double precision and integer declarations.
0011       IMPLICIT DOUBLE PRECISION(A-H, O-Z)
0012       IMPLICIT INTEGER(I-N)
0013       INTEGER PYK,PYCHGE,PYCOMP
0014 C...Commonblocks.
0015       COMMON/PYJETS/N,NPAD,K(4000,5),P(4000,5),V(4000,5)
0016       COMMON/PYDAT1/MSTU(200),PARU(200),MSTJ(200),PARJ(200)
0017       SAVE /PYJETS/,/PYDAT1/
0018 C...Local arrays.
0019       DIMENSION IJOIN(2),PTOT(4),BETA(3)
0020  
0021 C...Call PYHEPC to convert input from HEPEVT to PYJETS common.
0022       IF(ICOM.EQ.0) THEN
0023         MSTU(28)=0
0024         CALL PYHEPC(2)
0025       ENDIF
0026  
0027 C...Loop through entries and pick up all final partons.
0028       I1=0
0029       I2=0
0030       I3=0
0031       I4=0
0032       DO 100 I=1,N
0033       IF(K(I,1).LE.0.OR.K(I,1).GT.10) GOTO 100
0034       KFA=IABS(K(I,2))
0035       IF((KFA.GE.1.AND.KFA.LE.6).OR.KFA.EQ.21) THEN
0036         IF(K(I,2).GT.0.AND.K(I,2).LE.6) THEN
0037           IF(I1.EQ.0) THEN
0038             I1=I
0039           ELSEIF(I3.EQ.0) THEN
0040             I3=I
0041           ELSE
0042             CALL PYERRM(16,'(PY4JET:) more than two quarks')
0043           ENDIF
0044         ELSEIF(K(I,2).LT.0) THEN
0045           IF(I2.EQ.0) THEN
0046             I2=I
0047           ELSEIF(I4.EQ.0) THEN
0048             I4=I
0049           ELSE
0050             CALL PYERRM(16,'(PY4JET:) more than two antiquarks')
0051           ENDIF
0052         ELSE
0053           IF(I3.EQ.0) THEN
0054             I3=I
0055           ELSEIF(I4.EQ.0) THEN
0056             I4=I
0057           ELSE
0058             CALL PYERRM(16,'(PY4JET:) more than two gluons')
0059           ENDIF
0060         ENDIF
0061       ENDIF
0062   100 CONTINUE
0063  
0064 C...Check that event is arranged according to conventions.
0065       IF(I1.EQ.0.OR.I2.EQ.0.OR.I3.EQ.0.OR.I4.EQ.0) THEN
0066         CALL PYERRM(16,'(PY4JET:) event contains too few partons')
0067       ENDIF
0068       IF(I2.LT.I1.OR.I3.LT.I2.OR.I4.LT.I3) THEN
0069         CALL PYERRM(6,'(PY4JET:) partons arranged in wrong order')
0070       ENDIF
0071  
0072 C...Check whether second pair are quarks or gluons.
0073       IF(IABS(K(I3,2)).LT.10.AND.IABS(K(I4,2)).LT.10) THEN
0074         IQG34=1
0075       ELSEIF(K(I3,2).EQ.21.AND.K(I4,2).EQ.21) THEN
0076         IQG34=2
0077       ELSE
0078         CALL PYERRM(16,'(PY4JET:) second parton pair inconsistent')
0079       ENDIF
0080  
0081 C...Boost partons to their cm frame.
0082       DO 110 J=1,4
0083         PTOT(J)=P(I1,J)+P(I2,J)+P(I3,J)+P(I4,J)
0084   110 CONTINUE
0085       ECM=SQRT(MAX(0D0,PTOT(4)**2-PTOT(1)**2-PTOT(2)**2-PTOT(3)**2))
0086       DO 120 J=1,3
0087         BETA(J)=PTOT(J)/PTOT(4)
0088   120 CONTINUE
0089       CALL PYROBO(I1,I1,0D0,0D0,-BETA(1),-BETA(2),-BETA(3))
0090       CALL PYROBO(I2,I2,0D0,0D0,-BETA(1),-BETA(2),-BETA(3))
0091       CALL PYROBO(I3,I3,0D0,0D0,-BETA(1),-BETA(2),-BETA(3))
0092       CALL PYROBO(I4,I4,0D0,0D0,-BETA(1),-BETA(2),-BETA(3))
0093       NSAV=N
0094  
0095 C...Decide and set up shower history for q qbar q' qbar' events.
0096       IF(IQG34.EQ.1) THEN
0097         W1=PY4JTW(0,I1,I3,I4)
0098         W2=PY4JTW(0,I2,I3,I4)
0099         IF(W1.GT.PYR(0)*(W1+W2)) THEN
0100           CALL PY4JTS(0,I1,I3,I4,I2,QMAX)
0101         ELSE
0102           CALL PY4JTS(0,I2,I3,I4,I1,QMAX)
0103         ENDIF
0104  
0105 C...Decide and set up shower history for q qbar g g events.
0106       ELSE
0107         W1=PY4JTW(I1,I3,I2,I4)
0108         W2=PY4JTW(I1,I4,I2,I3)
0109         W3=PY4JTW(0,I3,I1,I4)
0110         W4=PY4JTW(0,I4,I1,I3)
0111         W5=PY4JTW(0,I3,I2,I4)
0112         W6=PY4JTW(0,I4,I2,I3)
0113         W7=PY4JTW(0,I1,I3,I4)
0114         W8=PY4JTW(0,I2,I3,I4)
0115         WR=(W1+W2+W3+W4+W5+W6+W7+W8)*PYR(0)
0116         IF(W1.GT.WR) THEN
0117           CALL PY4JTS(I1,I3,I2,I4,0,QMAX)
0118         ELSEIF(W1+W2.GT.WR) THEN
0119           CALL PY4JTS(I1,I4,I2,I3,0,QMAX)
0120         ELSEIF(W1+W2+W3.GT.WR) THEN
0121           CALL PY4JTS(0,I3,I1,I4,I2,QMAX)
0122         ELSEIF(W1+W2+W3+W4.GT.WR) THEN
0123           CALL PY4JTS(0,I4,I1,I3,I2,QMAX)
0124         ELSEIF(W1+W2+W3+W4+W5.GT.WR) THEN
0125           CALL PY4JTS(0,I3,I2,I4,I1,QMAX)
0126         ELSEIF(W1+W2+W3+W4+W5+W6.GT.WR) THEN
0127           CALL PY4JTS(0,I4,I2,I3,I1,QMAX)
0128         ELSEIF(W1+W2+W3+W4+W5+W6+W7.GT.WR) THEN
0129           CALL PY4JTS(0,I1,I3,I4,I2,QMAX)
0130         ELSE
0131           CALL PY4JTS(0,I2,I3,I4,I1,QMAX)
0132         ENDIF
0133       ENDIF
0134  
0135 C...Boost back original partons and mark them as deleted.
0136       CALL PYROBO(I1,I1,0D0,0D0,BETA(1),BETA(2),BETA(3))
0137       CALL PYROBO(I2,I2,0D0,0D0,BETA(1),BETA(2),BETA(3))
0138       CALL PYROBO(I3,I3,0D0,0D0,BETA(1),BETA(2),BETA(3))
0139       CALL PYROBO(I4,I4,0D0,0D0,BETA(1),BETA(2),BETA(3))
0140       K(I1,1)=K(I1,1)+10
0141       K(I2,1)=K(I2,1)+10
0142       K(I3,1)=K(I3,1)+10
0143       K(I4,1)=K(I4,1)+10
0144  
0145 C...Rotate shower initiating partons to be along z axis.
0146       PHI=PYANGL(P(NSAV+1,1),P(NSAV+1,2))
0147       CALL PYROBO(NSAV+1,NSAV+6,0D0,-PHI,0D0,0D0,0D0)
0148       THE=PYANGL(P(NSAV+1,3),P(NSAV+1,1))
0149       CALL PYROBO(NSAV+1,NSAV+6,-THE,0D0,0D0,0D0,0D0)
0150  
0151 C...Set up copy of shower initiating partons as on mass shell.
0152       DO 140 I=N+1,N+2
0153         DO 130 J=1,5
0154           K(I,J)=0
0155           P(I,J)=0D0
0156           V(I,J)=V(I1,J)
0157   130   CONTINUE
0158         K(I,1)=1
0159         K(I,2)=K(I-6,2)
0160   140 CONTINUE
0161       IF(K(NSAV+1,2).EQ.K(I1,2)) THEN
0162         K(N+1,3)=I1
0163         P(N+1,5)=P(I1,5)
0164         K(N+2,3)=I2
0165         P(N+2,5)=P(I2,5)
0166       ELSE
0167         K(N+1,3)=I2
0168         P(N+1,5)=P(I2,5)
0169         K(N+2,3)=I1
0170         P(N+2,5)=P(I1,5)
0171       ENDIF
0172       PABS=SQRT(MAX(0D0,(ECM**2-P(N+1,5)**2-P(N+2,5)**2)**2-
0173      &(2D0*P(N+1,5)*P(N+2,5))**2))/(2D0*ECM)
0174       P(N+1,3)=PABS
0175       P(N+1,4)=SQRT(PABS**2+P(N+1,5)**2)
0176       P(N+2,3)=-PABS
0177       P(N+2,4)=SQRT(PABS**2+P(N+2,5)**2)
0178       N=N+2
0179  
0180 C...Decide whether to allow or not photon radiation in showers.
0181 C...Connect up colours.
0182       MSTJ(41)=2
0183       IF(IRAD.EQ.0) MSTJ(41)=1
0184       IJOIN(1)=N-1
0185       IJOIN(2)=N
0186       CALL PYJOIN(2,IJOIN)
0187  
0188 C...Decide on maximum virtuality and do parton shower.
0189       IF(PMAX.LT.PARJ(82)) THEN
0190         PQMAX=QMAX
0191       ELSE
0192         PQMAX=PMAX
0193       ENDIF
0194       CALL PYSHOW(NSAV+1,-100,PQMAX)
0195  
0196 C...Rotate and boost back system.
0197       CALL PYROBO(NSAV+1,N,THE,PHI,BETA(1),BETA(2),BETA(3))
0198  
0199 C...Do fragmentation and decays.
0200       CALL PYEXEC
0201  
0202 C...Call PYHEPC to convert output from PYJETS to HEPEVT common.
0203       IF(ICOM.EQ.0) THEN
0204         MSTU(28)=0
0205         CALL PYHEPC(1)
0206       ENDIF
0207  
0208       RETURN
0209       END
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.3.7 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!