EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
lsmall.F
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file lsmall.F
1 C ********************************************************************
2 
3  SUBROUTINE lsmall
4 C-- --C
5 C-- Created: 951031 --C
6 C-- Last update: 980722 --C
7 C-- Purpose: Take care of small mass systems with one diquark, --C
8 C-- charmonium or bottonium --C
9 
10  IMPLICIT NONE
11 
12 C-- arg for spec. version
13  INTEGER ilon
14  REAL weight
15 
16 C-- global variables
17  INTEGER nlupdm,nplbuf
18  parameter(nlupdm=4000,nplbuf=5)
19  common/lujets/n,k(nlupdm,5),p(nlupdm,nplbuf),v(nlupdm,5)
20  INTEGER n,k
21  REAL p,v
22  SAVE /lujets/
23 
24 *
25 * to avoid variable conflictions, a second keep element is necessary
26 * with the same common block name (see LPTOU2)
27 *
28  COMMON /leptou/ cut(14),lst(40),parl(30),
29  & x,y,w2,q2,u
30  REAL cut,parl,x,y,w2,q2,u
31  INTEGER lst
32  SAVE /leptou/
33 
34  common/ludat1/mstu(200),paru(200),mstj(200),parj(200)
35  INTEGER mstu,mstj
36  REAL paru,parj
37  SAVE /ludat1/
38 
39  common/ludat2/kchg(500,3),pmas(500,4),parf(2000),vckm(4,4)
40  INTEGER kchg
41  REAL pmas,parf,vckm
42  SAVE /ludat2/
43 
44 
45 C-- functions
46  REAL rlu,ulmass
47  INTEGER lucomp
48 C-- local variables
49  INTEGER i,j,ifirst,i1,i2,iaq,iq,ntry1,ntry2,kfb,kfh1,kfh2
50  INTEGER kiaq,kiq,kdummy,nparton
51  REAL*8 energy,maxenergy,invmass,threshold
52  REAL*8 psum(5)
53  REAL*8 tot2,m1,m2,rotarg,pabs,costhe,ptemp,phi,pi
54  REAL*8 pcps,pc2,pn2,ps2,a,b,c,eps2,eps1
55  REAL*8 w441,w443,w445,w10441,w20443,w30443,tmp
56  REAL*8 w551,w553,w555,w10551,w20553,w30553
57  LOGICAL first,charmonium,bottonium,diquark
58  DATA pi/3.14159265359d0/
59 
60 C-- do not use spec. version
61  ilon=0
62 
63 C-- find light singlet system and hadronize them
64  first=.true.
65  nparton=n
66  charmonium=.false.
67  bottonium=.false.
68  diquark=.false.
69  energy=2.d0*parj(32)
70  DO 20 i=1,nparton
71  IF (k(i,1).EQ.2) THEN
72  IF (first) THEN
73  psum(1)=p(i,1)
74  psum(2)=p(i,2)
75  psum(3)=p(i,3)
76  psum(4)=p(i,4)
77  mstj(93)=1
78  psum(5)=ulmass(k(i,2))
79  first=.false.
80  ifirst=i
81  ELSE
82  psum(1)=psum(1)+p(i,1)
83  psum(2)=psum(2)+p(i,2)
84  psum(3)=psum(3)+p(i,3)
85  psum(4)=psum(4)+p(i,4)
86  mstj(93)=1
87  psum(5)=psum(5)+ulmass(k(i,2))
88  ENDIF
89  ELSEIF (k(i,1).EQ.1 .AND. .NOT. first) THEN
90  psum(1)=psum(1)+p(i,1)
91  psum(2)=psum(2)+p(i,2)
92  psum(3)=psum(3)+p(i,3)
93  psum(4)=psum(4)+p(i,4)
94  mstj(93)=1
95  psum(5)=psum(5)+ulmass(k(i,2))
96  invmass=psum(4)**2-psum(1)**2-psum(2)**2-psum(3)**2
97  energy=sqrt(max(0.d0,invmass))-psum(5)
98  first=.true.
99 C-- only light systems with two c- or b-quarks or one diquark are of interest
100  IF (energy.LT.dble(parj(32)) .AND.
101  & ((abs(k(ifirst,2)).EQ.4 .AND. abs(k(i,2)).EQ.4
102  & .OR. abs(k(ifirst,2)).EQ.5 .AND. abs(k(i,2)).EQ.5) .OR.
103  / (lucomp(k(ifirst,2)).EQ.90 .OR. lucomp(k(i,2)).EQ.90))) THEN
104 
105  i1=ifirst
106  i2=i
107 C-- add system as cluster
108  n=n+1
109  DO 30 j=1,4
110  p(n,j)=psum(j)
111  30 CONTINUE
112  p(n,5)=sqrt(max(0.d0,invmass))
113  k(n,1)=11
114  k(n,2)=91
115  k(n,3)=i1
116  k(n,4)=n+1
117  k(n,5)=n+2
118  IF(abs(k(i1,2)).EQ.4) THEN
119  charmonium=.true.
120  threshold=2.*ulmass(421)
121  ELSEIF(abs(k(i1,2)).EQ.5) THEN
122  bottonium=.true.
123  threshold=2.*ulmass(521)
124  ELSEIF(lucomp(k(i1,2)).EQ.90 .OR.
125  / lucomp(k(i2,2)).EQ.90) THEN
126  diquark=.true.
127  threshold=ulmass(2212)+ulmass(111)
128  ENDIF
129 C-- inactivate old system
130  DO 40 j=i1,i2
131  k(j,1)=k(j,1)+10
132  k(j,4)=n
133  40 CONTINUE
134 
135 C-- try to make two particles
136  ntry1=0
137  ntry2=0
138  50 CONTINUE
139 C-- take quark or diquark end first
140  IF (k(i1,2).EQ.4 .OR. k(i1,2).EQ.5 .OR.
141  / lucomp(k(i1,2)).EQ.90) THEN
142  iq=i1
143  iaq=i2
144  ELSE
145  iq=i2
146  iaq=i1
147  ENDIF
148  kiq=k(iq,2)
149  kiaq=k(iaq,2)
150  kdummy=0
151  kfb=0
152  kfh1=0
153  kfh2=0
154  CALL lukfdi(kiq,kdummy,kfb,kfh1)
155  kdummy=0
156  CALL lukfdi(kiaq,-kfb,kdummy,kfh2)
157  IF (kfh1.EQ.0 .OR. kfh2.EQ.0) THEN
158  ntry1=ntry1+1
159  IF (ntry1.GE.100) THEN
160  lst(21)=200
161  RETURN
162  ENDIF
163  goto 50
164  ENDIF
165 C-- consistency checks
166  IF (lucomp(kfh1).EQ.0 .OR. lucomp(kfh2).EQ.0) THEN
167  lst(21)=201
168  RETURN
169  ENDIF
170  IF (kchg(lucomp(kfh1),2)*isign(1,kfh1)+
171  + kchg(lucomp(kfh2),2)*isign(1,kfh2) .NE.0) THEN
172  lst(21)=202
173  RETURN
174  ENDIF
175  p(n+1,5)=ulmass(kfh1)
176  p(n+2,5)=ulmass(kfh2)
177  IF (p(n,5).LE.p(n+1,5)+p(n+2,5)+parj(64) .AND.
178  & p(n,5).GE.threshold+parj(64) .AND.
179  & ntry2.LE.100) THEN
180  ntry2=ntry2+1
181  goto 50
182  ENDIF
183  IF (p(n,5).GE.p(n+1,5)+p(n+2,5)+parj(64)) THEN
184 C-- make two particles
185 C-- isotropic decay in cms (dcostheta*dphi)
186  tot2=invmass
187  m1=dble(p(n+1,5))
188  m2=dble(p(n+2,5))
189  rotarg=(tot2-m1**2-m2**2)**2-4.d0*m1**2*m2**2
190  IF (rotarg.LT.0.d0) THEN
191  lst(21)=203
192  RETURN
193  ENDIF
194  pabs=0.5d0*sqrt(rotarg/tot2)
195  costhe=-1.d0+2.d0*rlu(0)
196  ptemp=pabs*sqrt(1.d0-costhe**2)
197  phi=2.d0*pi*rlu(0)
198  p(n+1,4)=sqrt(pabs**2+m1**2)
199  p(n+1,3)=pabs*costhe
200  p(n+1,2)=ptemp*cos(phi)
201  p(n+1,1)=ptemp*sin(phi)
202  p(n+2,4)=sqrt(pabs**2+m2**2)
203  p(n+2,3)=-p(n+1,3)
204  p(n+2,2)=-p(n+1,2)
205  p(n+2,1)=-p(n+1,1)
206 C-- K-vector
207  k(n+1,1)=1
208  k(n+1,2)=kfh1
209  k(n+1,3)=n
210  k(n+1,4)=0
211  k(n+1,5)=0
212  k(n+2,1)=1
213  k(n+2,2)=kfh2
214  k(n+2,3)=n
215  k(n+2,4)=0
216  k(n+2,5)=0
217 C-- boost to lab
218  mstu(33)=1
219  CALL ludbrb(n+1,n+2,0.,0.,psum(1)/psum(4),
220  & psum(2)/psum(4),psum(3)/psum(4))
221 C-- V-vector
222  DO 60 j=1,5
223  v(n,j)=v(iq,j)
224  v(n+1,j)=v(iq,j)
225  v(n+2,j)=v(iaq,j)
226  60 CONTINUE
227  v(n,5)=0.
228  v(n+1,5)=0.
229  v(n+2,5)=0.
230  n=n+2
231  ELSE
232 
233 
234 C-- make one particle instead
235  IF(charmonium) THEN
236 
237  w441=1.d0
238  w443=3.d0
239  w445=5.d0
240  w10441=1.d0
241  w20443=3.d0
242 C-- extra suppression of psiprime from radial exitation
243  w30443=3.d0/2.d0
244 
245  tmp=rlu(0)*(w441+w443+w445+w10441+w20443+w30443)
246  IF (tmp .LT. w441) THEN
247  kfh1=441
248  ELSEIF (tmp .LT. w441+w443) THEN
249  kfh1=443
250  ELSEIF (tmp .LT. w441+w443+w445) THEN
251  kfh1=445
252  ELSEIF (tmp .LT. w441+w443+w445+w10441) THEN
253  kfh1=10441
254  ELSEIF (tmp .LT. w441+w443+w445+w10441+w20443) THEN
255  kfh1=20443
256  ELSE
257  kfh1=30443
258  ENDIF
259 
260 C-- code for spec. version
261  IF (ilon.EQ.1) THEN
262  kfh1=441
263  weight=w441/(w441+w443+w445+w10441+w20443+w30443)
264  ENDIF
265  IF (ilon.EQ.2) THEN
266  kfh1=443
267  weight=w443/(w441+w443+w445+w10441+w20443+w30443)
268  ENDIF
269  IF (ilon.EQ.3) THEN
270  kfh1=445
271  weight=w445/(w441+w443+w445+w10441+w20443+w30443)
272  ENDIF
273  IF (ilon.EQ.4) THEN
274  kfh1=10441
275  weight=w10441/(w441+w443+w445+w10441+w20443+w30443)
276  ENDIF
277  IF (ilon.EQ.5) THEN
278  kfh1=20443
279  weight=w20443/(w441+w443+w445+w10441+w20443+w30443)
280  ENDIF
281  IF (ilon.EQ.6) THEN
282  kfh1=30443
283  weight=w30443/(w441+w443+w445+w10441+w20443+w30443)
284  ENDIF
285 
286  ELSEIF(bottonium) THEN
287 
288  w551=1.d0
289  w553=3.d0
290  w555=5.d0
291  w10551=1.d0
292  w20553=3.d0
293  w30553=3.d0/2.d0
294 
295  tmp=rlu(0)*(w551+w553+w555+w10551+w20553+w30553)
296  IF (tmp .LT. w551) THEN
297  kfh1=551
298  ELSEIF (tmp .LT. w551+w553) THEN
299  kfh1=553
300  ELSEIF (tmp .LT. w551+w553+w555) THEN
301  kfh1=555
302  ELSEIF (tmp .LT. w551+w553+w555+w10551) THEN
303  kfh1=10551
304  ELSEIF (tmp .LT. w551+w553+w555+w10551+w20553) THEN
305  kfh1=20553
306  ELSE
307  kfh1=30553
308  ENDIF
309 
310 C-- code for spec. version
311  IF (ilon.EQ.1) THEN
312  kfh1=551
313  weight=w551/(w551+w553+w555+w10551+w20553+w30553)
314  ENDIF
315  IF (ilon.EQ.2) THEN
316  kfh1=553
317  weight=w553/(w551+w553+w555+w10551+w20553+w30553)
318  ENDIF
319  IF (ilon.EQ.3) THEN
320  kfh1=555
321  weight=w555/(w551+w553+w555+w10551+w20553+w30553)
322  ENDIF
323  IF (ilon.EQ.4) THEN
324  kfh1=10551
325  weight=w10551/(w551+w553+w555+w10551+w20553+w30553)
326  ENDIF
327  IF (ilon.EQ.5) THEN
328  kfh1=20553
329  weight=w20553/(w551+w553+w555+w10551+w20553+w30553)
330  ENDIF
331  IF (ilon.EQ.6) THEN
332  kfh1=30553
333  weight=w30553/(w551+w553+w555+w10551+w20553+w30553)
334  ENDIF
335 
336  ELSEIF(diquark) THEN
337 70 CONTINUE
338  kiq=k(iq,2)
339  kiaq=k(iaq,2)
340  kdummy=0
341  kfh1=0
342  CALL lukfdi(kiq,kiaq,kdummy,kfh1)
343  IF (kfh1.EQ.0) goto 70
344 C-- isospin conservation
345  IF (kfh1.EQ.2214) kfh1=2212
346  IF (kfh1.EQ.2114) kfh1=2112
347  IF (kfh1.EQ.-2214) kfh1=-2212
348  IF (kfh1.EQ.-2114) kfh1=-2112
349 
350  ENDIF
351 
352  k(n,5)=n+1
353  k(n+1,1)=1
354  k(n+1,2)=kfh1
355  k(n+1,3)=n
356  k(n+1,4)=0
357  k(n+1,5)=0
358  p(n+1,5)=ulmass(kfh1)
359 C-- find particle to shuffle energy & momentum to and from
360  maxenergy=0.d0
361  i1=0
362  DO 80 j=1,n-1
363  IF (0.LT.k(j,1) .AND. k(j,1).LT.10 .AND.
364  & lucomp(k(j,2)).NE.0 .AND.
365  & (k(j,2).EQ.21 .OR. abs(k(j,2)).LT.10 .OR.
366  & abs(k(j,2)).GT.100) ) THEN
367  invmass=(p(n,4)+p(j,4))**2-(p(n,1)+p(j,1))**2-
368  - (p(n,2)+p(j,2))**2-(p(n,3)+p(j,3))**2
369  energy=sqrt(max(0.d0,invmass))-p(n,5)-p(j,5)
370  IF (energy.GT.maxenergy ) THEN
371  i1=j
372  maxenergy=energy
373  ENDIF
374  ENDIF
375  80 CONTINUE
376 C-- shuffle energy & momentum
377  IF (i1.NE.0) THEN
378  pcps=dble(p(n,4))*dble(p(i1,4))
379  - -dble(p(n,1))*dble(p(i1,1))
380  - -dble(p(n,2))*dble(p(i1,2))
381  - -dble(p(n,3))*dble(p(i1,3))
382  pc2=dble(p(n,5))**2
383  pn2=dble(p(n+1,5))**2
384  ps2=dble(p(i1,5))**2
385  a=pc2+ps2+2.d0*pcps
386  b=pc2+pn2+2.d0*pcps
387  c=(pn2-pc2)*(4.d0*pcps*(pcps+pc2)-pc2*(pn2-pc2))/
388  / 4.d0/(pcps**2-pc2*ps2)
389  IF (b**2-4.d0*c*a.LT.0.d0) THEN
390  lst(21)=204
391  RETURN
392  ENDIF
393  eps2=(-b+sqrt(max(0.d0,b**2-4.d0*c*a)))/2.d0/a
394  eps1=(pn2-pc2+2.d0*eps2*(ps2+pcps))/2.d0/(pc2+pcps)
395  DO 90 j=1,4
396  p(n+1,j)=(1.+eps1)*p(n,j)-eps2*p(i1,j)
397  p(i1,j)=(1.+eps2)*p(i1,j)-eps1*p(n,j)
398  v(n,j)=v(i1,j)
399  v(n+1,j)=v(i1,j)
400  90 CONTINUE
401  v(n,5)=0.
402  v(n+1,5)=0.
403  n=n+1
404  ELSE
405  lst(21)=205
406  RETURN
407  ENDIF
408  ENDIF
409  charmonium=.false.
410  bottonium=.false.
411  diquark=.false.
412  ENDIF
413  ENDIF
414  20 CONTINUE
415 
416  RETURN
417 
418  END