EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
luboei.f
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file luboei.f
1 
2 C*********************************************************************
3 
4  SUBROUTINE luboei(NSAV)
5 
6 C...Purpose: to modify event so as to approximately take into account
7 C...Bose-Einstein effects according to a simple phenomenological
8 C...parametrization.
9  IMPLICIT DOUBLE PRECISION(d)
10  common/lujets/n,k(9000,5),p(9000,5),v(9000,5)
11  SAVE /lujets/
12  common/ludat1/mstu(200),paru(200),mstj(200),parj(200)
13  SAVE /ludat1/
14  dimension dps(4),kfbe(9),nbe(0:9),bei(100)
15  DATA kfbe/211,-211,111,321,-321,130,310,221,331/
16 
17 C...Boost event to overall CM frame. Calculate CM energy.
18  IF((mstj(51).NE.1.AND.mstj(51).NE.2).OR.n-nsav.LE.1) RETURN
19  DO 100 j=1,4
20  100 dps(j)=0.
21  DO 120 i=1,n
22  IF(k(i,1).LE.0.OR.k(i,1).GT.10) goto 120
23  DO 110 j=1,4
24  110 dps(j)=dps(j)+p(i,j)
25  120 CONTINUE
26  CALL ludbrb(0,0,0.,0.,-dps(1)/dps(4),-dps(2)/dps(4),
27  &-dps(3)/dps(4))
28  pecm=0.
29  DO 130 i=1,n
30  130 IF(k(i,1).GE.1.AND.k(i,1).LE.10) pecm=pecm+p(i,4)
31 
32 C...Reserve copy of particles by species at end of record.
33  nbe(0)=n+mstu(3)
34  DO 160 ibe=1,min(9,mstj(51))
35  nbe(ibe)=nbe(ibe-1)
36  DO 150 i=nsav+1,n
37  IF(k(i,2).NE.kfbe(ibe)) goto 150
38  IF(k(i,1).LE.0.OR.k(i,1).GT.10) goto 150
39  IF(nbe(ibe).GE.mstu(4)-mstu(32)-5) THEN
40  CALL luerrm(11,'(LUBOEI:) no more memory left in LUJETS')
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(211)
59  IF(ibe.EQ.4) pmhq=2.*ulmass(321)
60  IF(ibe.EQ.8) pmhq=2.*ulmass(221)
61  IF(ibe.EQ.9) pmhq=2.*ulmass(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(0,0,0.,0.,dps(1)/dps(4),dps(2)/dps(4),dps(3)/dps(4))
140 
141  RETURN
142  END