HIJING/hipyset1_35/pyremn_hijing.F

   1 * $Id$
   2
   3 C*********************************************************************
   4
   5       SUBROUTINE PYREMN_HIJING(IPU1,IPU2)
   6
   7 C...Adds on target remnants (one or two from each side) and
   8 C...includes primordial kT.
   9 #include "hiparnt.inc"
  10 #include "histrng.inc"
  11 C...COMMON BLOCK FROM HIJING
  12 #include "lujets_hijing.inc"
  13 #include "ludat1_hijing.inc"
  14 #include "ludat2_hijing.inc"
  15 #include "pypars_hijing.inc"
  16 #include "pyint1_hijing.inc"
  17       DIMENSION KFLCH(2),KFLSP(2),CHI(2),PMS(6),IS(2),ROBO(5)
  18
  19 C...Special case for lepton-lepton interaction.
  20       IF(MINT(43).EQ.1) THEN
  21         DO 100 JT=1,2
  22         I=MINT(83)+JT+2
  23         K(I,1)=21
  24         K(I,2)=K(I-2,2)
  25         K(I,3)=I-2
  26         DO 100 J=1,5
  27   100   P(I,J)=P(I-2,J)
  28       ENDIF
  29
  30 C...Find event type, set pointers.
  31       IF(IPU1.EQ.0.AND.IPU2.EQ.0) RETURN
  32       ISUB=MINT(1)
  33       ILEP=0
  34       IF(IPU1.EQ.0) ILEP=1
  35       IF(IPU2.EQ.0) ILEP=2
  36       IF(ISUB.EQ.95) ILEP=-1
  37       IF(ILEP.EQ.1) IQ=MINT(84)+1
  38       IF(ILEP.EQ.2) IQ=MINT(84)+2
  39       IP=MAX(IPU1,IPU2)
  40       ILEPR=MINT(83)+5-ILEP
  41       NS=N
  42
  43 C...Define initial partons, including primordial kT.
  44   110 DO 130 JT=1,2
  45       I=MINT(83)+JT+2
  46       IF(JT.EQ.1) IPU=IPU1
  47       IF(JT.EQ.2) IPU=IPU2
  48       K(I,1)=21
  49       K(I,3)=I-2
  50       IF(ISUB.EQ.95) THEN
  51         K(I,2)=21
  52         SHS=0.
  53       ELSEIF(MINT(40+JT).EQ.1.AND.IPU.NE.0) THEN
  54         K(I,2)=K(IPU,2)
  55         P(I,5)=P(IPU,5)
  56         P(I,1)=0.
  57         P(I,2)=0.
  58         PMS(JT)=P(I,5)**2
  59       ELSEIF(IPU.NE.0) THEN
  60         K(I,2)=K(IPU,2)
  61         P(I,5)=P(IPU,5)
  62 C...No primordial kT or chosen according to truncated Gaussian or
  63 C...exponential.
  64 C
  65 c     X.N. Wang (7.22.97)
  66 c
  67         RPT1=0.0
  68         RPT2=0.0
  69         SS_W2=(PP(IHNT2(11),4)+PT(IHNT2(12),4))**2
  70      &       -(PP(IHNT2(11),1)+PT(IHNT2(12),1))**2
  71      &       -(PP(IHNT2(11),2)+PT(IHNT2(12),2))**2
  72      &       -(PP(IHNT2(11),3)+PT(IHNT2(12),3))**2
  73 C
  74 C********this is s of the current NN collision
  75         IF(SS_W2.LE.4.0*PARP(93)**2) GOTO 1211
  76 c
  77         IF(IHPR2(5).LE.0) THEN
  78 120          IF(MSTP(91).LE.0) THEN
  79                PTL=0.
  80              ELSEIF(MSTP(91).EQ.1) THEN
  81                PTL=PARP(91)*SQRT(-LOG(RLU_HIJING(0)))
  82              ELSE
  83                RPT1=RLU_HIJING(0)
  84                RPT2=RLU_HIJING(0)
  85                PTL=-PARP(92)*LOG(RPT1*RPT2)
  86              ENDIF
  87              IF(PTL.GT.PARP(93)) GOTO 120
  88              PHI=PARU(2)*RLU_HIJING(0)
  89              RPT1=PTL*COS(PHI)
  90              RPT2=PTL*SIN(PHI)
  91         ELSE IF(IHPR2(5).EQ.1) THEN
  92              IF(JT.EQ.1) JPT=NFP(IHNT2(11),11)
  93              IF(JT.EQ.2) JPT=NFT(IHNT2(12),11)
  94 1205         PTGS=PARP(91)*SQRT(-LOG(RLU_HIJING(0)))
  95              IF(PTGS.GT.PARP(93)) GO TO 1205
  96              PHI=2.0*HIPR1(40)*RLU_HIJING(0)
  97              RPT1=PTGS*COS(PHI)
  98              RPT2=PTGS*SIN(PHI)
  99              DO 1210 I_INT=1,JPT-1
 100                 PKCSQ=PARP(91)*SQRT(-LOG(RLU_HIJING(0)))
 101                 PHI=2.0*HIPR1(40)*RLU_HIJING(0)
 102                 RPT1=RPT1+PKCSQ*COS(PHI)
 103                 RPT2=RPT2+PKCSQ*SIN(PHI)
 104 1210         CONTINUE
 105              IF(RPT1**2+RPT2**2.GE.SS_W2/4.0) GO TO 1205
 106         ENDIF
 107 C     X.N. Wang
 108 C                       ********When initial interaction among soft partons is
 109 C                               assumed the primordial pt comes from the sum of
 110 C                               pt of JPT-1 number of initial interaction, JPT
 111 C                               is the number of interaction including present
 112 C                               one that nucleon hassuffered
 113 1211    P(I,1)=RPT1
 114         P(I,2)=RPT2
 115         PMS(JT)=P(I,5)**2+P(I,1)**2+P(I,2)**2
 116       ELSE
 117         K(I,2)=K(IQ,2)
 118         Q2=VINT(52)
 119         P(I,5)=-SQRT(Q2)
 120         PMS(JT)=-Q2
 121         SHS=(1.-VINT(43-JT))*Q2/VINT(43-JT)+VINT(5-JT)**2
 122       ENDIF
 123   130 CONTINUE
 124
 125 C...Kinematics construction for initial partons.
 126       I1=MINT(83)+3
 127       I2=MINT(83)+4
 128       IF(ILEP.EQ.0) SHS=VINT(141)*VINT(142)*VINT(2)+
 129      &(P(I1,1)+P(I2,1))**2+(P(I1,2)+P(I2,2))**2
 130       SHR=SQRT(MAX(0.,SHS))
 131       IF(ILEP.EQ.0) THEN
 132         IF((SHS-PMS(1)-PMS(2))**2-4.*PMS(1)*PMS(2).LE.0.) GOTO 110
 133         P(I1,4)=0.5*(SHR+(PMS(1)-PMS(2))/SHR)
 134         P(I1,3)=SQRT(MAX(0.,P(I1,4)**2-PMS(1)))
 135         P(I2,4)=SHR-P(I1,4)
 136         P(I2,3)=-P(I1,3)
 137       ELSEIF(ILEP.EQ.1) THEN
 138         P(I1,4)=P(IQ,4)
 139         P(I1,3)=P(IQ,3)
 140         P(I2,4)=P(IP,4)
 141         P(I2,3)=P(IP,3)
 142       ELSEIF(ILEP.EQ.2) THEN
 143         P(I1,4)=P(IP,4)
 144         P(I1,3)=P(IP,3)
 145         P(I2,4)=P(IQ,4)
 146         P(I2,3)=P(IQ,3)
 147       ENDIF
 148       IF(MINT(43).EQ.1) RETURN
 149
 150 C...Transform partons to overall CM-frame (not for leptoproduction).
 151       IF(ILEP.EQ.0) THEN
 152         ROBO(3)=(P(I1,1)+P(I2,1))/SHR
 153         ROBO(4)=(P(I1,2)+P(I2,2))/SHR
 154         CALL LUDBRB_HIJING(I1,I2,0.,0.,-DBLE(ROBO(3)),-DBLE(ROBO(4)),0D0
 155      $       )
 156         ROBO(2)=ULANGL_HIJING(P(I1,1),P(I1,2))
 157         CALL LUDBRB_HIJING(I1,I2,0.,-ROBO(2),0D0,0D0,0D0)
 158         ROBO(1)=ULANGL_HIJING(P(I1,3),P(I1,1))
 159         CALL LUDBRB_HIJING(I1,I2,-ROBO(1),0.,0D0,0D0,0D0)
 160         NMAX=MAX(MINT(52),IPU1,IPU2)
 161         CALL LUDBRB_HIJING(I1,NMAX,ROBO(1),ROBO(2),DBLE(ROBO(3))
 162      $       ,DBLE(ROBO(4)),0D0)
 163         ROBO(5)=MAX(-0.999999,MIN(0.999999,(VINT(141)-VINT(142))/
 164      &  (VINT(141)+VINT(142))))
 165         CALL LUDBRB_HIJING(I1,NMAX,0.,0.,0D0,0D0,DBLE(ROBO(5)))
 166       ENDIF
 167
 168 C...Check invariant mass of remnant system:
 169 C...hadronic events or leptoproduction.
 170       IF(ILEP.LE.0) THEN
 171         IF(MSTP(81).LE.0.OR.MSTP(82).LE.0.OR.ISUB.EQ.95) THEN
 172           VINT(151)=0.
 173           VINT(152)=0.
 174         ENDIF
 175         PEH=P(I1,4)+P(I2,4)+0.5*VINT(1)*(VINT(151)+VINT(152))
 176         PZH=P(I1,3)+P(I2,3)+0.5*VINT(1)*(VINT(151)-VINT(152))
 177         SHH=(VINT(1)-PEH)**2-(P(I1,1)+P(I2,1))**2-(P(I1,2)+P(I2,2))**2-
 178      &  PZH**2
 179         PMMIN=P(MINT(83)+1,5)+P(MINT(83)+2,5)+ULMASS_HIJING(K(I1,2))+
 180      &  ULMASS_HIJING(K(I2,2))
 181         IF(SHR.GE.VINT(1).OR.SHH.LE.(PMMIN+PARP(111))**2) THEN
 182           MINT(51)=1
 183           RETURN
 184         ENDIF
 185         SHR=SQRT(SHH+(P(I1,1)+P(I2,1))**2+(P(I1,2)+P(I2,2))**2)
 186       ELSE
 187         PEI=P(IQ,4)+P(IP,4)
 188         PZI=P(IQ,3)+P(IP,3)
 189         PMS(ILEP)=MAX(0.,PEI**2-PZI**2)
 190         PMMIN=P(ILEPR-2,5)+ULMASS_HIJING(K(ILEPR,2))+SQRT(PMS(ILEP))
 191         IF(SHR.LE.PMMIN+PARP(111)) THEN
 192           MINT(51)=1
 193           RETURN
 194         ENDIF
 195       ENDIF
 196
 197 C...Subdivide remnant if necessary, store first parton.
 198   140 I=NS
 199       DO 190 JT=1,2
 200       IF(JT.EQ.ILEP) GOTO 190
 201       IF(JT.EQ.1) IPU=IPU1
 202       IF(JT.EQ.2) IPU=IPU2
 203       CALL PYSPLI_HIJING(MINT(10+JT),MINT(12+JT),KFLCH(JT),KFLSP(JT))
 204       I=I+1
 205       IS(JT)=I
 206       DO 150 J=1,5
 207       K(I,J)=0
 208       P(I,J)=0.
 209   150 V(I,J)=0.
 210       K(I,1)=3
 211       K(I,2)=KFLSP(JT)
 212       K(I,3)=MINT(83)+JT
 213       P(I,5)=ULMASS_HIJING(K(I,2))
 214
 215 C...First parton colour connections and transverse mass.
 216       KFLS=(3-KCHG(LUCOMP_HIJING(KFLSP(JT)),2)*ISIGN(1,KFLSP(JT)))/2
 217       K(I,KFLS+3)=IPU
 218       K(IPU,6-KFLS)=MOD(K(IPU,6-KFLS),MSTU(5))+MSTU(5)*I
 219       IF(KFLCH(JT).EQ.0) THEN
 220         P(I,1)=-P(MINT(83)+JT+2,1)
 221         P(I,2)=-P(MINT(83)+JT+2,2)
 222         PMS(JT)=P(I,5)**2+P(I,1)**2+P(I,2)**2
 223
 224 C...When extra remnant parton or hadron: find relative pT, store.
 225       ELSE
 226         CALL LUPTDI_HIJING(1,P(I,1),P(I,2))
 227         PMS(JT+2)=P(I,5)**2+P(I,1)**2+P(I,2)**2
 228         I=I+1
 229         DO 160 J=1,5
 230         K(I,J)=0
 231         P(I,J)=0.
 232   160   V(I,J)=0.
 233         K(I,1)=1
 234         K(I,2)=KFLCH(JT)
 235         K(I,3)=MINT(83)+JT
 236         P(I,5)=ULMASS_HIJING(K(I,2))
 237         P(I,1)=-P(MINT(83)+JT+2,1)-P(I-1,1)
 238         P(I,2)=-P(MINT(83)+JT+2,2)-P(I-1,2)
 239         PMS(JT+4)=P(I,5)**2+P(I,1)**2+P(I,2)**2
 240 C...Relative distribution of energy for particle into two jets.
 241         IMB=1
 242         IF(MOD(MINT(10+JT)/1000,10).NE.0) IMB=2
 243         IF(IABS(KFLCH(JT)).LE.10.OR.KFLCH(JT).EQ.21) THEN
 244           CHIK=PARP(92+2*IMB)
 245           IF(MSTP(92).LE.1) THEN
 246             IF(IMB.EQ.1) CHI(JT)=RLU_HIJING(0)
 247             IF(IMB.EQ.2) CHI(JT)=1.-SQRT(RLU_HIJING(0))
 248           ELSEIF(MSTP(92).EQ.2) THEN
 249             CHI(JT)=1.-RLU_HIJING(0)**(1./(1.+CHIK))
 250           ELSEIF(MSTP(92).EQ.3) THEN
 251             CUT=2.*0.3/VINT(1)
 252   170       CHI(JT)=RLU_HIJING(0)**2
 253             IF((CHI(JT)**2/(CHI(JT)**2+CUT**2))**0.25*(1.-CHI(JT))**CHIK
 254      &      .LT.RLU_HIJING(0)) GOTO 170
 255           ELSE
 256             CUT=2.*0.3/VINT(1)
 257             CUTR=(1.+SQRT(1.+CUT**2))/CUT
 258   180       CHIR=CUT*CUTR**RLU_HIJING(0)
 259             CHI(JT)=(CHIR**2-CUT**2)/(2.*CHIR)
 260             IF((1.-CHI(JT))**CHIK.LT.RLU_HIJING(0)) GOTO 180
 261           ENDIF
 262 C...Relative distribution of energy for particle into jet plus particle.
 263         ELSE
 264   111     CONTINUE
 265           IF(MSTP(92).LE.1) THEN
 266             IF(IMB.EQ.1) CHI(JT)=RLU_HIJING(0)
 267             IF(IMB.EQ.2) CHI(JT)=1.-SQRT(RLU_HIJING(0))
 268           ELSE
 269             CHI(JT)=1.-RLU_HIJING(0)**(1./(1.+PARP(93+2*IMB)))
 270           ENDIF
 271           IF(MOD(KFLCH(JT)/1000,10).NE.0) CHI(JT)=1.-CHI(JT)
 272           IF (CHI(JT) .EQ. 1. .OR. CHI(JT) .EQ. 0.) GOTO 111
 273         ENDIF
 274         PMS(JT)=PMS(JT+4)/CHI(JT)+PMS(JT+2)/(1.-CHI(JT))
 275         KFLS=KCHG(LUCOMP_HIJING(KFLCH(JT)),2)*ISIGN(1,KFLCH(JT))
 276         IF(KFLS.NE.0) THEN
 277           K(I,1)=3
 278           KFLS=(3-KFLS)/2
 279           K(I,KFLS+3)=IPU
 280           K(IPU,6-KFLS)=MOD(K(IPU,6-KFLS),MSTU(5))+MSTU(5)*I
 281         ENDIF
 282       ENDIF
 283   190 CONTINUE
 284       IF(SHR.LE.SQRT(PMS(1))+SQRT(PMS(2))) GOTO 140
 285       N=I
 286
 287 C...Reconstruct kinematics of remnants.
 288       DO 200 JT=1,2
 289       IF(JT.EQ.ILEP) GOTO 200
 290       PE=0.5*(SHR+(PMS(JT)-PMS(3-JT))/SHR)
 291       PZ=SQRT(PE**2-PMS(JT))
 292       IF(KFLCH(JT).EQ.0) THEN
 293         P(IS(JT),4)=PE
 294         P(IS(JT),3)=PZ*(-1)**(JT-1)
 295       ELSE
 296         PW1=CHI(JT)*(PE+PZ)
 297         P(IS(JT)+1,4)=0.5*(PW1+PMS(JT+4)/PW1)
 298         P(IS(JT)+1,3)=0.5*(PW1-PMS(JT+4)/PW1)*(-1)**(JT-1)
 299         P(IS(JT),4)=PE-P(IS(JT)+1,4)
 300         P(IS(JT),3)=PZ*(-1)**(JT-1)-P(IS(JT)+1,3)
 301       ENDIF
 302   200 CONTINUE
 303
 304 C...Hadronic events: boost remnants to correct longitudinal frame.
 305       IF(ILEP.LE.0) THEN
 306          CALL LUDBRB_HIJING(NS+1,N,0.,0.,0D0,0D0,-DBLE(PZH/(VINT(1)-PEH)
 307      $        ))
 308 C...Leptoproduction events: boost colliding subsystem.
 309       ELSE
 310         NMAX=MAX(IP,MINT(52))
 311         PEF=SHR-PE
 312         PZF=PZ*(-1)**(ILEP-1)
 313         PT2=P(ILEPR,1)**2+P(ILEPR,2)**2
 314         PHIPT=ULANGL_HIJING(P(ILEPR,1),P(ILEPR,2))
 315         CALL LUDBRB_HIJING(MINT(84)+1,NMAX,0.,-PHIPT,0D0,0D0,0D0)
 316         RQP=P(IQ,3)*(PT2+PEI**2)-P(IQ,4)*PEI*PZI
 317         SINTH=P(IQ,4)*SQRT(PT2*(PT2+PEI**2)/(RQP**2+PT2*
 318      &  P(IQ,4)**2*PZI**2))*SIGN(1.,-RQP)
 319         CALL LUDBRB_HIJING(MINT(84)+1,NMAX,ASIN(SINTH),0.,0D0,0D0,0D0)
 320         BETAX=(-PEI*PZI*SINTH+SQRT(PT2*(PT2+PEI**2-(PZI*SINTH)**2)))/
 321      &  (PT2+PEI**2)
 322         CALL LUDBRB_HIJING(MINT(84)+1,NMAX,0.,0.,DBLE(BETAX),0D0,0D0)
 323         CALL LUDBRB_HIJING(MINT(84)+1,NMAX,0.,PHIPT,0D0,0D0,0D0)
 324         PEM=P(IQ,4)+P(IP,4)
 325         PZM=P(IQ,3)+P(IP,3)
 326         BETAZ=(-PEM*PZM+PZF*SQRT(PZF**2+PEM**2-PZM**2))/(PZF**2+PEM**2)
 327         CALL LUDBRB_HIJING(MINT(84)+1,NMAX,0.,0.,0D0,0D0,DBLE(BETAZ))
 328         CALL LUDBRB_HIJING(I1,I2,ASIN(SINTH),0.,DBLE(BETAX),0D0,0D0)
 329         CALL LUDBRB_HIJING(I1,I2,0.,PHIPT,0D0,0D0,DBLE(BETAZ))
 330       ENDIF
 331
 332       RETURN
 333       END