Merging the VirtualMC branch to the main development branch (HEAD)
[u/mrichter/AliRoot.git] / HIJING / hipyset1_35 / luboei_hijing.F
CommitLineData
e74335a4 1* $Id$
2
3C*********************************************************************
4
5 SUBROUTINE LUBOEI_HIJING(NSAV)
6
7C...Purpose: to modify event so as to approximately take into account
8C...Bose-Einstein effects according to a simple phenomenological
9C...parametrization.
10 IMPLICIT DOUBLE PRECISION(D)
11#include "lujets_hijing.inc"
12#include "ludat1_hijing.inc"
13 DIMENSION DPS(4),KFBE(9),NBE(0:9),BEI(100)
14 DATA KFBE/211,-211,111,321,-321,130,310,221,331/
15
16C...Boost event to overall CM frame. Calculate CM energy.
17 IF((MSTJ(51).NE.1.AND.MSTJ(51).NE.2).OR.N-NSAV.LE.1) RETURN
18 DO 100 J=1,4
19 100 DPS(J)=0.
20 DO 120 I=1,N
21 IF(K(I,1).LE.0.OR.K(I,1).GT.10) GOTO 120
22 DO 110 J=1,4
23 110 DPS(J)=DPS(J)+P(I,J)
24 120 CONTINUE
25 CALL LUDBRB_HIJING(0,0,0.,0.,-DPS(1)/DPS(4),-DPS(2)/DPS(4),
26 &-DPS(3)/DPS(4))
27 PECM=0.
28 DO 130 I=1,N
29 130 IF(K(I,1).GE.1.AND.K(I,1).LE.10) PECM=PECM+P(I,4)
30
31C...Reserve copy of particles by species at end of record.
32 NBE(0)=N+MSTU(3)
33 DO 160 IBE=1,MIN(9,MSTJ(51))
34 NBE(IBE)=NBE(IBE-1)
35 DO 150 I=NSAV+1,N
36 IF(K(I,2).NE.KFBE(IBE)) GOTO 150
37 IF(K(I,1).LE.0.OR.K(I,1).GT.10) GOTO 150
38 IF(NBE(IBE).GE.MSTU(4)-MSTU(32)-5) THEN
39 CALL LUERRM_HIJING(11
40 $ ,'(LUBOEI_HIJING:) no more memory left in LUJETS_HIJING')
41 RETURN
42 ENDIF
43 NBE(IBE)=NBE(IBE)+1
44 K(NBE(IBE),1)=I
45 DO 140 J=1,3
46 140 P(NBE(IBE),J)=0.
47 150 CONTINUE
48 160 CONTINUE
49
50C...Tabulate integral for subsequent momentum shift.
51 DO 210 IBE=1,MIN(9,MSTJ(51))
52 IF(IBE.NE.1.AND.IBE.NE.4.AND.IBE.LE.7) GOTO 180
53 IF(IBE.EQ.1.AND.MAX(NBE(1)-NBE(0),NBE(2)-NBE(1),NBE(3)-NBE(2)).
54 &LE.1) GOTO 180
55 IF(IBE.EQ.4.AND.MAX(NBE(4)-NBE(3),NBE(5)-NBE(4),NBE(6)-NBE(5),
56 &NBE(7)-NBE(6)).LE.1) GOTO 180
57 IF(IBE.GE.8.AND.NBE(IBE)-NBE(IBE-1).LE.1) GOTO 180
58 IF(IBE.EQ.1) PMHQ=2.*ULMASS_HIJING(211)
59 IF(IBE.EQ.4) PMHQ=2.*ULMASS_HIJING(321)
60 IF(IBE.EQ.8) PMHQ=2.*ULMASS_HIJING(221)
61 IF(IBE.EQ.9) PMHQ=2.*ULMASS_HIJING(331)
62 QDEL=0.1*MIN(PMHQ,PARJ(93))
63 IF(MSTJ(51).EQ.1) THEN
64 NBIN=MIN(100,NINT(9.*PARJ(93)/QDEL))
65 BEEX=EXP(0.5*QDEL/PARJ(93))
66 BERT=EXP(-QDEL/PARJ(93))
67 ELSE
68 NBIN=MIN(100,NINT(3.*PARJ(93)/QDEL))
69 ENDIF
70 DO 170 IBIN=1,NBIN
71 QBIN=QDEL*(IBIN-0.5)
72 BEI(IBIN)=QDEL*(QBIN**2+QDEL**2/12.)/SQRT(QBIN**2+PMHQ**2)
73 IF(MSTJ(51).EQ.1) THEN
74 BEEX=BEEX*BERT
75 BEI(IBIN)=BEI(IBIN)*BEEX
76 ELSE
77 BEI(IBIN)=BEI(IBIN)*EXP(-(QBIN/PARJ(93))**2)
78 ENDIF
79 170 IF(IBIN.GE.2) BEI(IBIN)=BEI(IBIN)+BEI(IBIN-1)
80
81C...Loop through particle pairs and find old relative momentum.
82 180 DO 200 I1M=NBE(IBE-1)+1,NBE(IBE)-1
83 I1=K(I1M,1)
84 DO 200 I2M=I1M+1,NBE(IBE)
85 I2=K(I2M,1)
86 Q2OLD=MAX(0.,(P(I1,4)+P(I2,4))**2-(P(I1,1)+P(I2,1))**2-(P(I1,2)+
87 &P(I2,2))**2-(P(I1,3)+P(I2,3))**2-(P(I1,5)+P(I2,5))**2)
88 QOLD=SQRT(Q2OLD)
89
90C...Calculate new relative momentum.
91 IF(QOLD.LT.0.5*QDEL) THEN
92 QMOV=QOLD/3.
93 ELSEIF(QOLD.LT.(NBIN-0.1)*QDEL) THEN
94 RBIN=QOLD/QDEL
95 IBIN=RBIN
96 RINP=(RBIN**3-IBIN**3)/(3*IBIN*(IBIN+1)+1)
97 QMOV=(BEI(IBIN)+RINP*(BEI(IBIN+1)-BEI(IBIN)))*
98 & SQRT(Q2OLD+PMHQ**2)/Q2OLD
99 ELSE
100 QMOV=BEI(NBIN)*SQRT(Q2OLD+PMHQ**2)/Q2OLD
101 ENDIF
102 Q2NEW=Q2OLD*(QOLD/(QOLD+3.*PARJ(92)*QMOV))**(2./3.)
103
104C...Calculate and save shift to be performed on three-momenta.
105 HC1=(P(I1,4)+P(I2,4))**2-(Q2OLD-Q2NEW)
106 HC2=(Q2OLD-Q2NEW)*(P(I1,4)-P(I2,4))**2
107 HA=0.5*(1.-SQRT(HC1*Q2NEW/(HC1*Q2OLD-HC2)))
108 DO 190 J=1,3
109 PD=HA*(P(I2,J)-P(I1,J))
110 P(I1M,J)=P(I1M,J)+PD
111 190 P(I2M,J)=P(I2M,J)-PD
112 200 CONTINUE
113 210 CONTINUE
114
115C...Shift momenta and recalculate energies.
116 DO 230 IM=NBE(0)+1,NBE(MIN(9,MSTJ(51)))
117 I=K(IM,1)
118 DO 220 J=1,3
119 220 P(I,J)=P(I,J)+P(IM,J)
120 230 P(I,4)=SQRT(P(I,5)**2+P(I,1)**2+P(I,2)**2+P(I,3)**2)
121
122C...Rescale all momenta for energy conservation.
123 PES=0.
124 PQS=0.
125 DO 240 I=1,N
126 IF(K(I,1).LE.0.OR.K(I,1).GT.10) GOTO 240
127 PES=PES+P(I,4)
128 PQS=PQS+P(I,5)**2/P(I,4)
129 240 CONTINUE
130 FAC=(PECM-PQS)/(PES-PQS)
131 DO 260 I=1,N
132 IF(K(I,1).LE.0.OR.K(I,1).GT.10) GOTO 260
133 DO 250 J=1,3
134 250 P(I,J)=FAC*P(I,J)
135 P(I,4)=SQRT(P(I,5)**2+P(I,1)**2+P(I,2)**2+P(I,3)**2)
136 260 CONTINUE
137
138C...Boost back to correct reference frame.
139 CALL LUDBRB_HIJING(0,0,0.,0.,DPS(1)/DPS(4),DPS(2)/DPS(4),DPS(3)
140 $ /DPS(4))
141
142 RETURN
143 END