]>
Commit | Line | Data |
---|---|---|
e74335a4 | 1 | * $Id$ |
2 | ||
3 | C********************************************************************* | |
4 | ||
5 | SUBROUTINE LUBOEI_HIJING(NSAV) | |
6 | ||
7 | C...Purpose: to modify event so as to approximately take into account | |
8 | C...Bose-Einstein effects according to a simple phenomenological | |
9 | C...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 | ||
16 | C...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 | ||
31 | C...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 | ||
50 | C...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 | ||
81 | C...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 | ||
90 | C...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 | ||
104 | C...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 | ||
115 | C...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 | ||
122 | C...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 | ||
138 | C...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 |