hijing/src/pyhiremn.f

0001
0002 C*********************************************************************
0003
0004       SUBROUTINE PYHIREMN(IPU1,IPU2)
0005
0006 C...Adds on target remnants (one or two from each side) and
0007 C...includes primordial kT.
0008       COMMON/HIPARNT/HIPR1(100),IHPR2(50),HINT1(100),IHNT2(50)
0009       SAVE  /HIPARNT/
0010       COMMON/HISTRNG/NFP(300,15),PPHI(300,15),NFT(300,15),PTHI(300,15)
0011       SAVE  /HISTRNG/
0012 C...COMMON BLOCK FROM HIJING
0013       COMMON/LUJETS/N,K(9000,5),P(9000,5),V(9000,5)
0014       SAVE /LUJETS/
0015       COMMON/LUDAT1/MSTU(200),PARU(200),MSTJ(200),PARJ(200)
0016       SAVE /LUDAT1/
0017       COMMON/LUDAT2/KCHG(500,3),PMAS(500,4),PARF(2000),VCKM(4,4)
0018       SAVE /LUDAT2/
0019       COMMON/PYHIPARS/MSTP(200),PARP(200),MSTI(200),PARI(200)
0020       SAVE /PYHIPARS/
0021       COMMON/PYHIINT1/MINT(400),VINT(400)
0022       SAVE /PYHIINT1/
0023       DIMENSION KFLCH(2),KFLSP(2),CHI(2),PMS(6),IS(2),ROBO(5)
0024
0025 C...Special case for lepton-lepton interaction.
0026       IF(MINT(43).EQ.1) THEN
0027         DO 100 JT=1,2
0028         I=MINT(83)+JT+2
0029         K(I,1)=21
0030         K(I,2)=K(I-2,2)
0031         K(I,3)=I-2
0032         DO 100 J=1,5
0033   100   P(I,J)=P(I-2,J)
0034       ENDIF
0035
0036 C...Find event type, set pointers.
0037       IF(IPU1.EQ.0.AND.IPU2.EQ.0) RETURN
0038       ISUB=MINT(1)
0039       ILEP=0
0040       IF(IPU1.EQ.0) ILEP=1
0041       IF(IPU2.EQ.0) ILEP=2
0042       IF(ISUB.EQ.95) ILEP=-1
0043       IF(ILEP.EQ.1) IQ=MINT(84)+1
0044       IF(ILEP.EQ.2) IQ=MINT(84)+2
0045       IP=MAX(IPU1,IPU2)
0046       ILEPR=MINT(83)+5-ILEP
0047       NS=N
0048
0049 C...Define initial partons, including primordial kT.
0050   110 DO 130 JT=1,2
0051       I=MINT(83)+JT+2
0052       IF(JT.EQ.1) IPU=IPU1
0053       IF(JT.EQ.2) IPU=IPU2
0054       K(I,1)=21
0055       K(I,3)=I-2
0056       IF(ISUB.EQ.95) THEN
0057         K(I,2)=21
0058         SHS=0.
0059       ELSEIF(MINT(40+JT).EQ.1.AND.IPU.NE.0) THEN
0060         K(I,2)=K(IPU,2)
0061         P(I,5)=P(IPU,5)
0062         P(I,1)=0.
0063         P(I,2)=0.
0064         PMS(JT)=P(I,5)**2
0065       ELSEIF(IPU.NE.0) THEN
0066         K(I,2)=K(IPU,2)
0067         P(I,5)=P(IPU,5)
0068 C...No primordial kT or chosen according to truncated Gaussian or
0069 C...exponential.
0070 C
0071 c     X.N. Wang (7.22.97)
0072 c
0073         RPT1=0.0
0074         RPT2=0.0
0075         SS_W2=(PPHI(IHNT2(11),4)+PTHI(IHNT2(12),4))**2
0076      &       -(PPHI(IHNT2(11),1)+PTHI(IHNT2(12),1))**2
0077      &       -(PPHI(IHNT2(11),2)+PTHI(IHNT2(12),2))**2
0078      &       -(PPHI(IHNT2(11),3)+PTHI(IHNT2(12),3))**2
0079 C
0080 C********this is s of the current NN collision
0081         IF(SS_W2.LE.4.0*PARP(93)**2) GOTO 1211
0082 c
0083         IF(IHPR2(5).LE.0) THEN
0084 120          IF(MSTP(91).LE.0) THEN
0085                PT=0.
0086              ELSEIF(MSTP(91).EQ.1) THEN
0087                PT=PARP(91)*SQRT(-LOG(RLU(0)))
0088              ELSE
0089                RPT1=RLU(0)
0090                RPT2=RLU(0)
0091                PT=-PARP(92)*LOG(RPT1*RPT2)
0092              ENDIF
0093              IF(PT.GT.PARP(93)) GOTO 120
0094              PHI=PARU(2)*RLU(0)
0095              RPT1=PT*COS(PHI)
0096              RPT2=PT*SIN(PHI)
0097         ELSE IF(IHPR2(5).EQ.1) THEN
0098              IF(JT.EQ.1) JPT=NFP(IHNT2(11),11)
0099              IF(JT.EQ.2) JPT=NFT(IHNT2(12),11)
0100 1205         PTGS=PARP(91)*SQRT(-LOG(RLU(0)))
0101              IF(PTGS.GT.PARP(93)) GO TO 1205
0102              PHI=2.0*HIPR1(40)*RLU(0)
0103              RPT1=PTGS*COS(PHI)
0104              RPT2=PTGS*SIN(PHI)
0105              DO 1210 I_INT=1,JPT-1
0106                 PKCSQ=PARP(91)*SQRT(-LOG(RLU(0)))
0107                 PHI=2.0*HIPR1(40)*RLU(0)
0108                 RPT1=RPT1+PKCSQ*COS(PHI)
0109                 RPT2=RPT2+PKCSQ*SIN(PHI)
0110 1210         CONTINUE
0111              IF(RPT1**2+RPT2**2.GE.SS_W2/4.0) GO TO 1205
0112         ENDIF
0113 C     X.N. Wang
0114 C                       ********When initial interaction among soft partons is
0115 C                               assumed the primordial pt comes from the sum of
0116 C                               pt of JPT-1 number of initial interaction, JPT
0117 C                               is the number of interaction including present
0118 C                               one that nucleon hassuffered
0119 1211    P(I,1)=RPT1
0120         P(I,2)=RPT2
0121         PMS(JT)=P(I,5)**2+P(I,1)**2+P(I,2)**2
0122       ELSE
0123         K(I,2)=K(IQ,2)
0124         Q2=VINT(52)
0125         P(I,5)=-SQRT(Q2)
0126         PMS(JT)=-Q2
0127         SHS=(1.-VINT(43-JT))*Q2/VINT(43-JT)+VINT(5-JT)**2
0128       ENDIF
0129   130 CONTINUE
0130
0131 C...Kinematics construction for initial partons.
0132       I1=MINT(83)+3
0133       I2=MINT(83)+4
0134       IF(ILEP.EQ.0) SHS=VINT(141)*VINT(142)*VINT(2)+
0135      &(P(I1,1)+P(I2,1))**2+(P(I1,2)+P(I2,2))**2
0136       SHR=SQRT(MAX(0.,SHS))
0137       IF(ILEP.EQ.0) THEN
0138         IF((SHS-PMS(1)-PMS(2))**2-4.*PMS(1)*PMS(2).LE.0.) GOTO 110
0139         P(I1,4)=0.5*(SHR+(PMS(1)-PMS(2))/SHR)
0140         P(I1,3)=SQRT(MAX(0.,P(I1,4)**2-PMS(1)))
0141         P(I2,4)=SHR-P(I1,4)
0142         P(I2,3)=-P(I1,3)
0143       ELSEIF(ILEP.EQ.1) THEN
0144         P(I1,4)=P(IQ,4)
0145         P(I1,3)=P(IQ,3)
0146         P(I2,4)=P(IP,4)
0147         P(I2,3)=P(IP,3)
0148       ELSEIF(ILEP.EQ.2) THEN
0149         P(I1,4)=P(IP,4)
0150         P(I1,3)=P(IP,3)
0151         P(I2,4)=P(IQ,4)
0152         P(I2,3)=P(IQ,3)
0153       ENDIF
0154       IF(MINT(43).EQ.1) RETURN
0155
0156 C...Transform partons to overall CM-frame (not for leptoproduction).
0157       IF(ILEP.EQ.0) THEN
0158         ROBO(3)=(P(I1,1)+P(I2,1))/SHR
0159         ROBO(4)=(P(I1,2)+P(I2,2))/SHR
0160         CALL LUDBRB(I1,I2,0.,0.,-DBLE(ROBO(3)),-DBLE(ROBO(4)),0D0)
0161         ROBO(2)=ULANGL(P(I1,1),P(I1,2))
0162         CALL LUDBRB(I1,I2,0.,-ROBO(2),0D0,0D0,0D0)
0163         ROBO(1)=ULANGL(P(I1,3),P(I1,1))
0164         CALL LUDBRB(I1,I2,-ROBO(1),0.,0D0,0D0,0D0)
0165         NMAX=MAX(MINT(52),IPU1,IPU2)
0166         CALL LUDBRB(I1,NMAX,ROBO(1),ROBO(2),DBLE(ROBO(3)),DBLE(ROBO(4)),
0167      &  0D0)
0168         ROBO(5)=MAX(-0.999999,MIN(0.999999,(VINT(141)-VINT(142))/
0169      &  (VINT(141)+VINT(142))))
0170         CALL LUDBRB(I1,NMAX,0.,0.,0D0,0D0,DBLE(ROBO(5)))
0171       ENDIF
0172
0173 C...Check invariant mass of remnant system:
0174 C...hadronic events or leptoproduction.
0175       IF(ILEP.LE.0) THEN
0176         IF(MSTP(81).LE.0.OR.MSTP(82).LE.0.OR.ISUB.EQ.95) THEN
0177           VINT(151)=0.
0178           VINT(152)=0.
0179         ENDIF
0180         PEH=P(I1,4)+P(I2,4)+0.5*VINT(1)*(VINT(151)+VINT(152))
0181         PZH=P(I1,3)+P(I2,3)+0.5*VINT(1)*(VINT(151)-VINT(152))
0182         SHH=(VINT(1)-PEH)**2-(P(I1,1)+P(I2,1))**2-(P(I1,2)+P(I2,2))**2-
0183      &  PZH**2
0184         PMMIN=P(MINT(83)+1,5)+P(MINT(83)+2,5)+ULMASS(K(I1,2))+
0185      &  ULMASS(K(I2,2))
0186         IF(SHR.GE.VINT(1).OR.SHH.LE.(PMMIN+PARP(111))**2) THEN
0187           MINT(51)=1
0188           RETURN
0189         ENDIF
0190         SHR=SQRT(SHH+(P(I1,1)+P(I2,1))**2+(P(I1,2)+P(I2,2))**2)
0191       ELSE
0192         PEI=P(IQ,4)+P(IP,4)
0193         PZI=P(IQ,3)+P(IP,3)
0194         PMS(ILEP)=MAX(0.,PEI**2-PZI**2)
0195         PMMIN=P(ILEPR-2,5)+ULMASS(K(ILEPR,2))+SQRT(PMS(ILEP))
0196         IF(SHR.LE.PMMIN+PARP(111)) THEN
0197           MINT(51)=1
0198           RETURN
0199         ENDIF
0200       ENDIF
0201
0202 C...Subdivide remnant if necessary, store first parton.
0203   140 I=NS
0204       DO 190 JT=1,2
0205       IF(JT.EQ.ILEP) GOTO 190
0206       IF(JT.EQ.1) IPU=IPU1
0207       IF(JT.EQ.2) IPU=IPU2
0208       CALL PYHISPLI(MINT(10+JT),MINT(12+JT),KFLCH(JT),KFLSP(JT))
0209       I=I+1
0210       IS(JT)=I
0211       DO 150 J=1,5
0212       K(I,J)=0
0213       P(I,J)=0.
0214   150 V(I,J)=0.
0215       K(I,1)=3
0216       K(I,2)=KFLSP(JT)
0217       K(I,3)=MINT(83)+JT
0218       P(I,5)=ULMASS(K(I,2))
0219
0220 C...First parton colour connections and transverse mass.
0221       KFLS=(3-KCHG(LUCOMP(KFLSP(JT)),2)*ISIGN(1,KFLSP(JT)))/2
0222       K(I,KFLS+3)=IPU
0223       K(IPU,6-KFLS)=MOD(K(IPU,6-KFLS),MSTU(5))+MSTU(5)*I
0224       IF(KFLCH(JT).EQ.0) THEN
0225         P(I,1)=-P(MINT(83)+JT+2,1)
0226         P(I,2)=-P(MINT(83)+JT+2,2)
0227         PMS(JT)=P(I,5)**2+P(I,1)**2+P(I,2)**2
0228
0229 C...When extra remnant parton or hadron: find relative pT, store.
0230       ELSE
0231         CALL LUPTDI(1,P(I,1),P(I,2))
0232         PMS(JT+2)=P(I,5)**2+P(I,1)**2+P(I,2)**2
0233         I=I+1
0234         DO 160 J=1,5
0235         K(I,J)=0
0236         P(I,J)=0.
0237   160   V(I,J)=0.
0238         K(I,1)=1
0239         K(I,2)=KFLCH(JT)
0240         K(I,3)=MINT(83)+JT
0241         P(I,5)=ULMASS(K(I,2))
0242         P(I,1)=-P(MINT(83)+JT+2,1)-P(I-1,1)
0243         P(I,2)=-P(MINT(83)+JT+2,2)-P(I-1,2)
0244         PMS(JT+4)=P(I,5)**2+P(I,1)**2+P(I,2)**2
0245 C...Relative distribution of energy for particle into two jets.
0246         IMB=1
0247         IF(MOD(MINT(10+JT)/1000,10).NE.0) IMB=2
0248         IF(IABS(KFLCH(JT)).LE.10.OR.KFLCH(JT).EQ.21) THEN
0249           CHIK=PARP(92+2*IMB)
0250           IF(MSTP(92).LE.1) THEN
0251             IF(IMB.EQ.1) CHI(JT)=RLU(0)
0252             IF(IMB.EQ.2) CHI(JT)=1.-SQRT(RLU(0))
0253           ELSEIF(MSTP(92).EQ.2) THEN
0254             CHI(JT)=1.-RLU(0)**(1./(1.+CHIK))
0255           ELSEIF(MSTP(92).EQ.3) THEN
0256             CUT=2.*0.3/VINT(1)
0257   170       CHI(JT)=RLU(0)**2
0258             IF((CHI(JT)**2/(CHI(JT)**2+CUT**2))**0.25*(1.-CHI(JT))**CHIK
0259      &      .LT.RLU(0)) GOTO 170
0260           ELSE
0261             CUT=2.*0.3/VINT(1)
0262             CUTR=(1.+SQRT(1.+CUT**2))/CUT
0263   180       CHIR=CUT*CUTR**RLU(0)
0264             CHI(JT)=(CHIR**2-CUT**2)/(2.*CHIR)
0265             IF((1.-CHI(JT))**CHIK.LT.RLU(0)) GOTO 180
0266           ENDIF
0267 C...Relative distribution of energy for particle into jet plus particle.
0268         ELSE
0269           IF(MSTP(92).LE.1) THEN
0270             IF(IMB.EQ.1) CHI(JT)=RLU(0)
0271             IF(IMB.EQ.2) CHI(JT)=1.-SQRT(RLU(0))
0272           ELSE
0273             CHI(JT)=1.-RLU(0)**(1./(1.+PARP(93+2*IMB)))
0274           ENDIF
0275           IF(MOD(KFLCH(JT)/1000,10).NE.0) CHI(JT)=1.-CHI(JT)
0276         ENDIF
0277         PMS(JT)=PMS(JT+4)/CHI(JT)+PMS(JT+2)/(1.-CHI(JT))
0278         KFLS=KCHG(LUCOMP(KFLCH(JT)),2)*ISIGN(1,KFLCH(JT))
0279         IF(KFLS.NE.0) THEN
0280           K(I,1)=3
0281           KFLS=(3-KFLS)/2
0282           K(I,KFLS+3)=IPU
0283           K(IPU,6-KFLS)=MOD(K(IPU,6-KFLS),MSTU(5))+MSTU(5)*I
0284         ENDIF
0285       ENDIF
0286   190 CONTINUE
0287       IF(SHR.LE.SQRT(PMS(1))+SQRT(PMS(2))) GOTO 140
0288       N=I
0289
0290 C...Reconstruct kinematics of remnants.
0291       DO 200 JT=1,2
0292       IF(JT.EQ.ILEP) GOTO 200
0293       PE=0.5*(SHR+(PMS(JT)-PMS(3-JT))/SHR)
0294       PZ=SQRT(PE**2-PMS(JT))
0295       IF(KFLCH(JT).EQ.0) THEN
0296         P(IS(JT),4)=PE
0297         P(IS(JT),3)=PZ*(-1)**(JT-1)
0298       ELSE
0299         PW1=CHI(JT)*(PE+PZ)
0300         P(IS(JT)+1,4)=0.5*(PW1+PMS(JT+4)/PW1)
0301         P(IS(JT)+1,3)=0.5*(PW1-PMS(JT+4)/PW1)*(-1)**(JT-1)
0302         P(IS(JT),4)=PE-P(IS(JT)+1,4)
0303         P(IS(JT),3)=PZ*(-1)**(JT-1)-P(IS(JT)+1,3)
0304       ENDIF
0305   200 CONTINUE
0306
0307 C...Hadronic events: boost remnants to correct longitudinal frame.
0308       IF(ILEP.LE.0) THEN
0309         CALL LUDBRB(NS+1,N,0.,0.,0D0,0D0,-DBLE(PZH/(VINT(1)-PEH)))
0310 C...Leptoproduction events: boost colliding subsystem.
0311       ELSE
0312         NMAX=MAX(IP,MINT(52))
0313         PEF=SHR-PE
0314         PZF=PZ*(-1)**(ILEP-1)
0315         PT2=P(ILEPR,1)**2+P(ILEPR,2)**2
0316         PHIPT=ULANGL(P(ILEPR,1),P(ILEPR,2))
0317         CALL LUDBRB(MINT(84)+1,NMAX,0.,-PHIPT,0D0,0D0,0D0)
0318         RQP=P(IQ,3)*(PT2+PEI**2)-P(IQ,4)*PEI*PZI
0319         SINTH=P(IQ,4)*SQRT(PT2*(PT2+PEI**2)/(RQP**2+PT2*
0320      &  P(IQ,4)**2*PZI**2))*SIGN(1.,-RQP)
0321         CALL LUDBRB(MINT(84)+1,NMAX,ASIN(SINTH),0.,0D0,0D0,0D0)
0322         BETAX=(-PEI*PZI*SINTH+SQRT(PT2*(PT2+PEI**2-(PZI*SINTH)**2)))/
0323      &  (PT2+PEI**2)
0324         CALL LUDBRB(MINT(84)+1,NMAX,0.,0.,DBLE(BETAX),0D0,0D0)
0325         CALL LUDBRB(MINT(84)+1,NMAX,0.,PHIPT,0D0,0D0,0D0)
0326         PEM=P(IQ,4)+P(IP,4)
0327         PZM=P(IQ,3)+P(IP,3)
0328         BETAZ=(-PEM*PZM+PZF*SQRT(PZF**2+PEM**2-PZM**2))/(PZF**2+PEM**2)
0329         CALL LUDBRB(MINT(84)+1,NMAX,0.,0.,0D0,0D0,DBLE(BETAZ))
0330         CALL LUDBRB(I1,I2,ASIN(SINTH),0.,DBLE(BETAX),0D0,0D0)
0331         CALL LUDBRB(I1,I2,0.,PHIPT,0D0,0D0,DBLE(BETAZ))
0332       ENDIF
0333
0334       RETURN
0335       END