Merging the VirtualMC branch to the main development branch (HEAD)
[u/mrichter/AliRoot.git] / 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           IF(MSTP(92).LE.1) THEN    
265             IF(IMB.EQ.1) CHI(JT)=RLU_HIJING(0) 
266             IF(IMB.EQ.2) CHI(JT)=1.-SQRT(RLU_HIJING(0))    
267           ELSE  
268             CHI(JT)=1.-RLU_HIJING(0)**(1./(1.+PARP(93+2*IMB))) 
269           ENDIF 
270           IF(MOD(KFLCH(JT)/1000,10).NE.0) CHI(JT)=1.-CHI(JT)    
271         ENDIF   
272         PMS(JT)=PMS(JT+4)/CHI(JT)+PMS(JT+2)/(1.-CHI(JT))    
273         KFLS=KCHG(LUCOMP_HIJING(KFLCH(JT)),2)*ISIGN(1,KFLCH(JT))   
274         IF(KFLS.NE.0) THEN  
275           K(I,1)=3  
276           KFLS=(3-KFLS)/2   
277           K(I,KFLS+3)=IPU   
278           K(IPU,6-KFLS)=MOD(K(IPU,6-KFLS),MSTU(5))+MSTU(5)*I    
279         ENDIF   
280       ENDIF 
281   190 CONTINUE  
282       IF(SHR.LE.SQRT(PMS(1))+SQRT(PMS(2))) GOTO 140 
283       N=I   
284     
285 C...Reconstruct kinematics of remnants. 
286       DO 200 JT=1,2 
287       IF(JT.EQ.ILEP) GOTO 200   
288       PE=0.5*(SHR+(PMS(JT)-PMS(3-JT))/SHR)  
289       PZ=SQRT(PE**2-PMS(JT))    
290       IF(KFLCH(JT).EQ.0) THEN   
291         P(IS(JT),4)=PE  
292         P(IS(JT),3)=PZ*(-1)**(JT-1) 
293       ELSE  
294         PW1=CHI(JT)*(PE+PZ) 
295         P(IS(JT)+1,4)=0.5*(PW1+PMS(JT+4)/PW1)   
296         P(IS(JT)+1,3)=0.5*(PW1-PMS(JT+4)/PW1)*(-1)**(JT-1)  
297         P(IS(JT),4)=PE-P(IS(JT)+1,4)    
298         P(IS(JT),3)=PZ*(-1)**(JT-1)-P(IS(JT)+1,3)   
299       ENDIF 
300   200 CONTINUE  
301     
302 C...Hadronic events: boost remnants to correct longitudinal frame.  
303       IF(ILEP.LE.0) THEN    
304          CALL LUDBRB_HIJING(NS+1,N,0.,0.,0D0,0D0,-DBLE(PZH/(VINT(1)-PEH)
305      $        ))  
306 C...Leptoproduction events: boost colliding subsystem.  
307       ELSE  
308         NMAX=MAX(IP,MINT(52))   
309         PEF=SHR-PE  
310         PZF=PZ*(-1)**(ILEP-1)   
311         PT2=P(ILEPR,1)**2+P(ILEPR,2)**2 
312         PHIPT=ULANGL_HIJING(P(ILEPR,1),P(ILEPR,2)) 
313         CALL LUDBRB_HIJING(MINT(84)+1,NMAX,0.,-PHIPT,0D0,0D0,0D0)  
314         RQP=P(IQ,3)*(PT2+PEI**2)-P(IQ,4)*PEI*PZI    
315         SINTH=P(IQ,4)*SQRT(PT2*(PT2+PEI**2)/(RQP**2+PT2*    
316      &  P(IQ,4)**2*PZI**2))*SIGN(1.,-RQP)   
317         CALL LUDBRB_HIJING(MINT(84)+1,NMAX,ASIN(SINTH),0.,0D0,0D0,0D0) 
318         BETAX=(-PEI*PZI*SINTH+SQRT(PT2*(PT2+PEI**2-(PZI*SINTH)**2)))/   
319      &  (PT2+PEI**2)    
320         CALL LUDBRB_HIJING(MINT(84)+1,NMAX,0.,0.,DBLE(BETAX),0D0,0D0)  
321         CALL LUDBRB_HIJING(MINT(84)+1,NMAX,0.,PHIPT,0D0,0D0,0D0)   
322         PEM=P(IQ,4)+P(IP,4) 
323         PZM=P(IQ,3)+P(IP,3) 
324         BETAZ=(-PEM*PZM+PZF*SQRT(PZF**2+PEM**2-PZM**2))/(PZF**2+PEM**2) 
325         CALL LUDBRB_HIJING(MINT(84)+1,NMAX,0.,0.,0D0,0D0,DBLE(BETAZ))  
326         CALL LUDBRB_HIJING(I1,I2,ASIN(SINTH),0.,DBLE(BETAX),0D0,0D0)   
327         CALL LUDBRB_HIJING(I1,I2,0.,PHIPT,0D0,0D0,DBLE(BETAZ)) 
328       ENDIF 
329     
330       RETURN    
331       END