Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
lutabu.f
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file lutabu.f
1 
2 C*********************************************************************
3 
4  SUBROUTINE lutabu(MTABU)
5 
6 C...Purpose: to evaluate various properties of an event, with
7 C...statistics accumulated during the course of the run and
8 C...printed at the end.
9  common/lujets/n,k(9000,5),p(9000,5),v(9000,5)
10  SAVE /lujets/
11  common/ludat1/mstu(200),paru(200),mstj(200),parj(200)
12  SAVE /ludat1/
13  common/ludat2/kchg(500,3),pmas(500,4),parf(2000),vckm(4,4)
14  SAVE /ludat2/
15  common/ludat3/mdcy(500,3),mdme(2000,2),brat(2000),kfdp(2000,5)
16  SAVE /ludat3/
17  dimension kfis(100,2),npis(100,0:10),kffs(400),npfs(400,4),
18  &fevfm(10,4),fm1fm(3,10,4),fm2fm(3,10,4),fmoma(4),fmoms(4),
19  &fevee(50),fe1ec(50),fe2ec(50),fe1ea(25),fe2ea(25),
20  &kfdm(8),kfdc(200,0:8),npdc(200)
21  SAVE nevis,nkfis,kfis,npis,nevfs,nprfs,nfifs,nchfs,nkffs,
22  &kffs,npfs,nevfm,nmufm,fm1fm,fm2fm,nevee,fe1ec,fe2ec,fe1ea,
23  &fe2ea,nevdc,nkfdc,nredc,kfdc,npdc
24  CHARACTER chau*16,chis(2)*12,chdc(8)*12
25  DATA nevis/0/,nkfis/0/,nevfs/0/,nprfs/0/,nfifs/0/,nchfs/0/,
26  &nkffs/0/,nevfm/0/,nmufm/0/,fm1fm/120*0./,fm2fm/120*0./,
27  &nevee/0/,fe1ec/50*0./,fe2ec/50*0./,fe1ea/25*0./,fe2ea/25*0./,
28  &nevdc/0/,nkfdc/0/,nredc/0/
29 
30 C...Reset statistics on initial parton state.
31  IF(mtabu.EQ.10) THEN
32  nevis=0
33  nkfis=0
34 
35 C...Identify and order flavour content of initial state.
36  ELSEIF(mtabu.EQ.11) THEN
37  nevis=nevis+1
38  kfm1=2*iabs(mstu(161))
39  IF(mstu(161).GT.0) kfm1=kfm1-1
40  kfm2=2*iabs(mstu(162))
41  IF(mstu(162).GT.0) kfm2=kfm2-1
42  kfmn=min(kfm1,kfm2)
43  kfmx=max(kfm1,kfm2)
44  DO 100 i=1,nkfis
45  IF(kfmn.EQ.kfis(i,1).AND.kfmx.EQ.kfis(i,2)) THEN
46  ikfis=-i
47  goto 110
48  ELSEIF(kfmn.LT.kfis(i,1).OR.(kfmn.EQ.kfis(i,1).AND.
49  & kfmx.LT.kfis(i,2))) THEN
50  ikfis=i
51  goto 110
52  ENDIF
53  100 CONTINUE
54  ikfis=nkfis+1
55  110 IF(ikfis.LT.0) THEN
56  ikfis=-ikfis
57  ELSE
58  IF(nkfis.GE.100) RETURN
59  DO 120 i=nkfis,ikfis,-1
60  kfis(i+1,1)=kfis(i,1)
61  kfis(i+1,2)=kfis(i,2)
62  DO 120 j=0,10
63  120 npis(i+1,j)=npis(i,j)
64  nkfis=nkfis+1
65  kfis(ikfis,1)=kfmn
66  kfis(ikfis,2)=kfmx
67  DO 130 j=0,10
68  130 npis(ikfis,j)=0
69  ENDIF
70  npis(ikfis,0)=npis(ikfis,0)+1
71 
72 C...Count number of partons in initial state.
73  np=0
74  DO 150 i=1,n
75  IF(k(i,1).LE.0.OR.k(i,1).GT.12) THEN
76  ELSEIF(iabs(k(i,2)).GT.80.AND.iabs(k(i,2)).LE.100) THEN
77  ELSEIF(iabs(k(i,2)).GT.100.AND.mod(iabs(k(i,2))/10,10).NE.0)
78  & THEN
79  ELSE
80  im=i
81  140 im=k(im,3)
82  IF(im.LE.0.OR.im.GT.n) THEN
83  np=np+1
84  ELSEIF(k(im,1).LE.0.OR.k(im,1).GT.20) THEN
85  np=np+1
86  ELSEIF(iabs(k(im,2)).GT.80.AND.iabs(k(im,2)).LE.100) THEN
87  ELSEIF(iabs(k(im,2)).GT.100.AND.mod(iabs(k(im,2))/10,10).NE.0)
88  & THEN
89  ELSE
90  goto 140
91  ENDIF
92  ENDIF
93  150 CONTINUE
94  npco=max(np,1)
95  IF(np.GE.6) npco=6
96  IF(np.GE.8) npco=7
97  IF(np.GE.11) npco=8
98  IF(np.GE.16) npco=9
99  IF(np.GE.26) npco=10
100  npis(ikfis,npco)=npis(ikfis,npco)+1
101  mstu(62)=np
102 
103 C...Write statistics on initial parton state.
104  ELSEIF(mtabu.EQ.12) THEN
105  fac=1./max(1,nevis)
106  WRITE(mstu(11),1000) nevis
107  DO 160 i=1,nkfis
108  kfmn=kfis(i,1)
109  IF(kfmn.EQ.0) kfmn=kfis(i,2)
110  kfm1=(kfmn+1)/2
111  IF(2*kfm1.EQ.kfmn) kfm1=-kfm1
112  CALL luname(kfm1,chau)
113  chis(1)=chau(1:12)
114  IF(chau(13:13).NE.' ') chis(1)(12:12)='?'
115  kfmx=kfis(i,2)
116  IF(kfis(i,1).EQ.0) kfmx=0
117  kfm2=(kfmx+1)/2
118  IF(2*kfm2.EQ.kfmx) kfm2=-kfm2
119  CALL luname(kfm2,chau)
120  chis(2)=chau(1:12)
121  IF(chau(13:13).NE.' ') chis(2)(12:12)='?'
122  160 WRITE(mstu(11),1100) chis(1),chis(2),fac*npis(i,0),
123  & (npis(i,j)/float(npis(i,0)),j=1,10)
124 
125 C...Copy statistics on initial parton state into /LUJETS/.
126  ELSEIF(mtabu.EQ.13) THEN
127  fac=1./max(1,nevis)
128  DO 170 i=1,nkfis
129  kfmn=kfis(i,1)
130  IF(kfmn.EQ.0) kfmn=kfis(i,2)
131  kfm1=(kfmn+1)/2
132  IF(2*kfm1.EQ.kfmn) kfm1=-kfm1
133  kfmx=kfis(i,2)
134  IF(kfis(i,1).EQ.0) kfmx=0
135  kfm2=(kfmx+1)/2
136  IF(2*kfm2.EQ.kfmx) kfm2=-kfm2
137  k(i,1)=32
138  k(i,2)=99
139  k(i,3)=kfm1
140  k(i,4)=kfm2
141  k(i,5)=npis(i,0)
142  DO 170 j=1,5
143  p(i,j)=fac*npis(i,j)
144  170 v(i,j)=fac*npis(i,j+5)
145  n=nkfis
146  DO 180 j=1,5
147  k(n+1,j)=0
148  p(n+1,j)=0.
149  180 v(n+1,j)=0.
150  k(n+1,1)=32
151  k(n+1,2)=99
152  k(n+1,5)=nevis
153  mstu(3)=1
154 
155 C...Reset statistics on number of particles/partons.
156  ELSEIF(mtabu.EQ.20) THEN
157  nevfs=0
158  nprfs=0
159  nfifs=0
160  nchfs=0
161  nkffs=0
162 
163 C...Identify whether particle/parton is primary or not.
164  ELSEIF(mtabu.EQ.21) THEN
165  nevfs=nevfs+1
166  mstu(62)=0
167  DO 230 i=1,n
168  IF(k(i,1).LE.0.OR.k(i,1).GT.20.OR.k(i,1).EQ.13) goto 230
169  mstu(62)=mstu(62)+1
170  kc=lucomp(k(i,2))
171  mpri=0
172  IF(k(i,3).LE.0.OR.k(i,3).GT.n) THEN
173  mpri=1
174  ELSEIF(k(k(i,3),1).LE.0.OR.k(k(i,3),1).GT.20) THEN
175  mpri=1
176  ELSEIF(kc.EQ.0) THEN
177  ELSEIF(k(k(i,3),1).EQ.13) THEN
178  im=k(k(i,3),3)
179  IF(im.LE.0.OR.im.GT.n) THEN
180  mpri=1
181  ELSEIF(k(im,1).LE.0.OR.k(im,1).GT.20) THEN
182  mpri=1
183  ENDIF
184  ELSEIF(kchg(kc,2).EQ.0) THEN
185  kcm=lucomp(k(k(i,3),2))
186  IF(kcm.NE.0) THEN
187  IF(kchg(kcm,2).NE.0) mpri=1
188  ENDIF
189  ENDIF
190  IF(kc.NE.0.AND.mpri.EQ.1) THEN
191  IF(kchg(kc,2).EQ.0) nprfs=nprfs+1
192  ENDIF
193  IF(k(i,1).LE.10) THEN
194  nfifs=nfifs+1
195  IF(luchge(k(i,2)).NE.0) nchfs=nchfs+1
196  ENDIF
197 
198 C...Fill statistics on number of particles/partons in event.
199  kfa=iabs(k(i,2))
200  kfs=3-isign(1,k(i,2))-mpri
201  DO 190 ip=1,nkffs
202  IF(kfa.EQ.kffs(ip)) THEN
203  ikffs=-ip
204  goto 200
205  ELSEIF(kfa.LT.kffs(ip)) THEN
206  ikffs=ip
207  goto 200
208  ENDIF
209  190 CONTINUE
210  ikffs=nkffs+1
211  200 IF(ikffs.LT.0) THEN
212  ikffs=-ikffs
213  ELSE
214  IF(nkffs.GE.400) RETURN
215  DO 210 ip=nkffs,ikffs,-1
216  kffs(ip+1)=kffs(ip)
217  DO 210 j=1,4
218  210 npfs(ip+1,j)=npfs(ip,j)
219  nkffs=nkffs+1
220  kffs(ikffs)=kfa
221  DO 220 j=1,4
222  220 npfs(ikffs,j)=0
223  ENDIF
224  npfs(ikffs,kfs)=npfs(ikffs,kfs)+1
225  230 CONTINUE
226 
227 C...Write statistics on particle/parton composition of events.
228  ELSEIF(mtabu.EQ.22) THEN
229  fac=1./max(1,nevfs)
230  WRITE(mstu(11),1200) nevfs,fac*nprfs,fac*nfifs,fac*nchfs
231  DO 240 i=1,nkffs
232  CALL luname(kffs(i),chau)
233  kc=lucomp(kffs(i))
234  mdcyf=0
235  IF(kc.NE.0) mdcyf=mdcy(kc,1)
236  240 WRITE(mstu(11),1300) kffs(i),chau,mdcyf,(fac*npfs(i,j),j=1,4),
237  & fac*(npfs(i,1)+npfs(i,2)+npfs(i,3)+npfs(i,4))
238 
239 C...Copy particle/parton composition information into /LUJETS/.
240  ELSEIF(mtabu.EQ.23) THEN
241  fac=1./max(1,nevfs)
242  DO 260 i=1,nkffs
243  k(i,1)=32
244  k(i,2)=99
245  k(i,3)=kffs(i)
246  k(i,4)=0
247  k(i,5)=npfs(i,1)+npfs(i,2)+npfs(i,3)+npfs(i,4)
248  DO 250 j=1,4
249  p(i,j)=fac*npfs(i,j)
250  250 v(i,j)=0.
251  p(i,5)=fac*k(i,5)
252  260 v(i,5)=0.
253  n=nkffs
254  DO 270 j=1,5
255  k(n+1,j)=0
256  p(n+1,j)=0.
257  270 v(n+1,j)=0.
258  k(n+1,1)=32
259  k(n+1,2)=99
260  k(n+1,5)=nevfs
261  p(n+1,1)=fac*nprfs
262  p(n+1,2)=fac*nfifs
263  p(n+1,3)=fac*nchfs
264  mstu(3)=1
265 
266 C...Reset factorial moments statistics.
267  ELSEIF(mtabu.EQ.30) THEN
268  nevfm=0
269  nmufm=0
270  DO 280 im=1,3
271  DO 280 ib=1,10
272  DO 280 ip=1,4
273  fm1fm(im,ib,ip)=0.
274  280 fm2fm(im,ib,ip)=0.
275 
276 C...Find particles to include, with (pion,pseudo)rapidity and azimuth.
277  ELSEIF(mtabu.EQ.31) THEN
278  nevfm=nevfm+1
279  nlow=n+mstu(3)
280  nupp=nlow
281  DO 360 i=1,n
282  IF(k(i,1).LE.0.OR.k(i,1).GT.10) goto 360
283  IF(mstu(41).GE.2) THEN
284  kc=lucomp(k(i,2))
285  IF(kc.EQ.0.OR.kc.EQ.12.OR.kc.EQ.14.OR.kc.EQ.16.OR.
286  & kc.EQ.18) goto 360
287  IF(mstu(41).GE.3.AND.kchg(kc,2).EQ.0.AND.luchge(k(i,2)).EQ.0)
288  & goto 360
289  ENDIF
290  pmr=0.
291  IF(mstu(42).EQ.1.AND.k(i,2).NE.22) pmr=ulmass(211)
292  IF(mstu(42).GE.2) pmr=p(i,5)
293  pr=max(1e-20,pmr**2+p(i,1)**2+p(i,2)**2)
294  yeta=sign(log(min((sqrt(pr+p(i,3)**2)+abs(p(i,3)))/sqrt(pr),
295  & 1e20)),p(i,3))
296  IF(abs(yeta).GT.paru(57)) goto 360
297  phi=ulangl(p(i,1),p(i,2))
298  iyeta=512.*(yeta+paru(57))/(2.*paru(57))
299  iyeta=max(0,min(511,iyeta))
300  iphi=512.*(phi+paru(1))/paru(2)
301  iphi=max(0,min(511,iphi))
302  iyep=0
303  DO 290 ib=0,9
304  290 iyep=iyep+4**ib*(2*mod(iyeta/2**ib,2)+mod(iphi/2**ib,2))
305 
306 C...Order particles in (pseudo)rapidity and/or azimuth.
307  IF(nupp.GT.mstu(4)-5-mstu(32)) THEN
308  CALL luerrm(11,'(LUTABU:) no more memory left in LUJETS')
309  RETURN
310  ENDIF
311  nupp=nupp+1
312  IF(nupp.EQ.nlow+1) THEN
313  k(nupp,1)=iyeta
314  k(nupp,2)=iphi
315  k(nupp,3)=iyep
316  ELSE
317  DO 300 i1=nupp-1,nlow+1,-1
318  IF(iyeta.GE.k(i1,1)) goto 310
319  300 k(i1+1,1)=k(i1,1)
320  310 k(i1+1,1)=iyeta
321  DO 320 i1=nupp-1,nlow+1,-1
322  IF(iphi.GE.k(i1,2)) goto 330
323  320 k(i1+1,2)=k(i1,2)
324  330 k(i1+1,2)=iphi
325  DO 340 i1=nupp-1,nlow+1,-1
326  IF(iyep.GE.k(i1,3)) goto 350
327  340 k(i1+1,3)=k(i1,3)
328  350 k(i1+1,3)=iyep
329  ENDIF
330  360 CONTINUE
331  k(nupp+1,1)=2**10
332  k(nupp+1,2)=2**10
333  k(nupp+1,3)=4**10
334 
335 C...Calculate sum of factorial moments in event.
336  DO 400 im=1,3
337  DO 370 ib=1,10
338  DO 370 ip=1,4
339  370 fevfm(ib,ip)=0.
340  DO 380 ib=1,10
341  IF(im.LE.2) ibin=2**(10-ib)
342  IF(im.EQ.3) ibin=4**(10-ib)
343  iagr=k(nlow+1,im)/ibin
344  nagr=1
345  DO 380 i=nlow+2,nupp+1
346  icut=k(i,im)/ibin
347  IF(icut.EQ.iagr) THEN
348  nagr=nagr+1
349  ELSE
350  IF(nagr.EQ.1) THEN
351  ELSEIF(nagr.EQ.2) THEN
352  fevfm(ib,1)=fevfm(ib,1)+2.
353  ELSEIF(nagr.EQ.3) THEN
354  fevfm(ib,1)=fevfm(ib,1)+6.
355  fevfm(ib,2)=fevfm(ib,2)+6.
356  ELSEIF(nagr.EQ.4) THEN
357  fevfm(ib,1)=fevfm(ib,1)+12.
358  fevfm(ib,2)=fevfm(ib,2)+24.
359  fevfm(ib,3)=fevfm(ib,3)+24.
360  ELSE
361  fevfm(ib,1)=fevfm(ib,1)+nagr*(nagr-1.)
362  fevfm(ib,2)=fevfm(ib,2)+nagr*(nagr-1.)*(nagr-2.)
363  fevfm(ib,3)=fevfm(ib,3)+nagr*(nagr-1.)*(nagr-2.)*(nagr-3.)
364  fevfm(ib,4)=fevfm(ib,4)+nagr*(nagr-1.)*(nagr-2.)*(nagr-3.)*
365  & (nagr-4.)
366  ENDIF
367  iagr=icut
368  nagr=1
369  ENDIF
370  380 CONTINUE
371 
372 C...Add results to total statistics.
373  DO 390 ib=10,1,-1
374  DO 390 ip=1,4
375  IF(fevfm(1,ip).LT.0.5) THEN
376  fevfm(ib,ip)=0.
377  ELSEIF(im.LE.2) THEN
378  fevfm(ib,ip)=2**((ib-1)*ip)*fevfm(ib,ip)/fevfm(1,ip)
379  ELSE
380  fevfm(ib,ip)=4**((ib-1)*ip)*fevfm(ib,ip)/fevfm(1,ip)
381  ENDIF
382  fm1fm(im,ib,ip)=fm1fm(im,ib,ip)+fevfm(ib,ip)
383  390 fm2fm(im,ib,ip)=fm2fm(im,ib,ip)+fevfm(ib,ip)**2
384  400 CONTINUE
385  nmufm=nmufm+(nupp-nlow)
386  mstu(62)=nupp-nlow
387 
388 C...Write accumulated statistics on factorial moments.
389  ELSEIF(mtabu.EQ.32) THEN
390  fac=1./max(1,nevfm)
391  IF(mstu(42).LE.0) WRITE(mstu(11),1400) nevfm,'eta'
392  IF(mstu(42).EQ.1) WRITE(mstu(11),1400) nevfm,'ypi'
393  IF(mstu(42).GE.2) WRITE(mstu(11),1400) nevfm,'y '
394  DO 420 im=1,3
395  WRITE(mstu(11),1500)
396  DO 420 ib=1,10
397  byeta=2.*paru(57)
398  IF(im.NE.2) byeta=byeta/2**(ib-1)
399  bphi=paru(2)
400  IF(im.NE.1) bphi=bphi/2**(ib-1)
401  IF(im.LE.2) bnave=fac*nmufm/float(2**(ib-1))
402  IF(im.EQ.3) bnave=fac*nmufm/float(4**(ib-1))
403  DO 410 ip=1,4
404  fmoma(ip)=fac*fm1fm(im,ib,ip)
405  410 fmoms(ip)=sqrt(max(0.,fac*(fac*fm2fm(im,ib,ip)-fmoma(ip)**2)))
406  420 WRITE(mstu(11),1600) byeta,bphi,bnave,(fmoma(ip),fmoms(ip),
407  & ip=1,4)
408 
409 C...Copy statistics on factorial moments into /LUJETS/.
410  ELSEIF(mtabu.EQ.33) THEN
411  fac=1./max(1,nevfm)
412  DO 430 im=1,3
413  DO 430 ib=1,10
414  i=10*(im-1)+ib
415  k(i,1)=32
416  k(i,2)=99
417  k(i,3)=1
418  IF(im.NE.2) k(i,3)=2**(ib-1)
419  k(i,4)=1
420  IF(im.NE.1) k(i,4)=2**(ib-1)
421  k(i,5)=0
422  p(i,1)=2.*paru(57)/k(i,3)
423  v(i,1)=paru(2)/k(i,4)
424  DO 430 ip=1,4
425  p(i,ip+1)=fac*fm1fm(im,ib,ip)
426  430 v(i,ip+1)=sqrt(max(0.,fac*(fac*fm2fm(im,ib,ip)-p(i,ip+1)**2)))
427  n=30
428  DO 440 j=1,5
429  k(n+1,j)=0
430  p(n+1,j)=0.
431  440 v(n+1,j)=0.
432  k(n+1,1)=32
433  k(n+1,2)=99
434  k(n+1,5)=nevfm
435  mstu(3)=1
436 
437 C...Reset statistics on Energy-Energy Correlation.
438  ELSEIF(mtabu.EQ.40) THEN
439  nevee=0
440  DO 450 j=1,25
441  fe1ec(j)=0.
442  fe2ec(j)=0.
443  fe1ec(51-j)=0.
444  fe2ec(51-j)=0.
445  fe1ea(j)=0.
446  450 fe2ea(j)=0.
447 
448 C...Find particles to include, with proper assumed mass.
449  ELSEIF(mtabu.EQ.41) THEN
450  nevee=nevee+1
451  nlow=n+mstu(3)
452  nupp=nlow
453  ecm=0.
454  DO 460 i=1,n
455  IF(k(i,1).LE.0.OR.k(i,1).GT.10) goto 460
456  IF(mstu(41).GE.2) THEN
457  kc=lucomp(k(i,2))
458  IF(kc.EQ.0.OR.kc.EQ.12.OR.kc.EQ.14.OR.kc.EQ.16.OR.
459  & kc.EQ.18) goto 460
460  IF(mstu(41).GE.3.AND.kchg(kc,2).EQ.0.AND.luchge(k(i,2)).EQ.0)
461  & goto 460
462  ENDIF
463  pmr=0.
464  IF(mstu(42).EQ.1.AND.k(i,2).NE.22) pmr=ulmass(211)
465  IF(mstu(42).GE.2) pmr=p(i,5)
466  IF(nupp.GT.mstu(4)-5-mstu(32)) THEN
467  CALL luerrm(11,'(LUTABU:) no more memory left in LUJETS')
468  RETURN
469  ENDIF
470  nupp=nupp+1
471  p(nupp,1)=p(i,1)
472  p(nupp,2)=p(i,2)
473  p(nupp,3)=p(i,3)
474  p(nupp,4)=sqrt(pmr**2+p(i,1)**2+p(i,2)**2+p(i,3)**2)
475  p(nupp,5)=max(1e-10,sqrt(p(i,1)**2+p(i,2)**2+p(i,3)**2))
476  ecm=ecm+p(nupp,4)
477  460 CONTINUE
478  IF(nupp.EQ.nlow) RETURN
479 
480 C...Analyze Energy-Energy Correlation in event.
481  fac=(2./ecm**2)*50./paru(1)
482  DO 470 j=1,50
483  470 fevee(j)=0.
484  DO 480 i1=nlow+2,nupp
485  DO 480 i2=nlow+1,i1-1
486  cthe=(p(i1,1)*p(i2,1)+p(i1,2)*p(i2,2)+p(i1,3)*p(i2,3))/
487  & (p(i1,5)*p(i2,5))
488  the=acos(max(-1.,min(1.,cthe)))
489  ithe=max(1,min(50,1+int(50.*the/paru(1))))
490  480 fevee(ithe)=fevee(ithe)+fac*p(i1,4)*p(i2,4)
491  DO 490 j=1,25
492  fe1ec(j)=fe1ec(j)+fevee(j)
493  fe2ec(j)=fe2ec(j)+fevee(j)**2
494  fe1ec(51-j)=fe1ec(51-j)+fevee(51-j)
495  fe2ec(51-j)=fe2ec(51-j)+fevee(51-j)**2
496  fe1ea(j)=fe1ea(j)+(fevee(51-j)-fevee(j))
497  490 fe2ea(j)=fe2ea(j)+(fevee(51-j)-fevee(j))**2
498  mstu(62)=nupp-nlow
499 
500 C...Write statistics on Energy-Energy Correlation.
501  ELSEIF(mtabu.EQ.42) THEN
502  fac=1./max(1,nevee)
503  WRITE(mstu(11),1700) nevee
504  DO 500 j=1,25
505  feec1=fac*fe1ec(j)
506  fees1=sqrt(max(0.,fac*(fac*fe2ec(j)-feec1**2)))
507  feec2=fac*fe1ec(51-j)
508  fees2=sqrt(max(0.,fac*(fac*fe2ec(51-j)-feec2**2)))
509  feeca=fac*fe1ea(j)
510  feesa=sqrt(max(0.,fac*(fac*fe2ea(j)-feeca**2)))
511  500 WRITE(mstu(11),1800) 3.6*(j-1),3.6*j,feec1,fees1,feec2,fees2,
512  & feeca,feesa
513 
514 C...Copy statistics on Energy-Energy Correlation into /LUJETS/.
515  ELSEIF(mtabu.EQ.43) THEN
516  fac=1./max(1,nevee)
517  DO 510 i=1,25
518  k(i,1)=32
519  k(i,2)=99
520  k(i,3)=0
521  k(i,4)=0
522  k(i,5)=0
523  p(i,1)=fac*fe1ec(i)
524  v(i,1)=sqrt(max(0.,fac*(fac*fe2ec(i)-p(i,1)**2)))
525  p(i,2)=fac*fe1ec(51-i)
526  v(i,2)=sqrt(max(0.,fac*(fac*fe2ec(51-i)-p(i,2)**2)))
527  p(i,3)=fac*fe1ea(i)
528  v(i,3)=sqrt(max(0.,fac*(fac*fe2ea(i)-p(i,3)**2)))
529  p(i,4)=paru(1)*(i-1)/50.
530  p(i,5)=paru(1)*i/50.
531  v(i,4)=3.6*(i-1)
532  510 v(i,5)=3.6*i
533  n=25
534  DO 520 j=1,5
535  k(n+1,j)=0
536  p(n+1,j)=0.
537  520 v(n+1,j)=0.
538  k(n+1,1)=32
539  k(n+1,2)=99
540  k(n+1,5)=nevee
541  mstu(3)=1
542 
543 C...Reset statistics on decay channels.
544  ELSEIF(mtabu.EQ.50) THEN
545  nevdc=0
546  nkfdc=0
547  nredc=0
548 
549 C...Identify and order flavour content of final state.
550  ELSEIF(mtabu.EQ.51) THEN
551  nevdc=nevdc+1
552  nds=0
553  DO 550 i=1,n
554  IF(k(i,1).LE.0.OR.k(i,1).GE.6) goto 550
555  nds=nds+1
556  IF(nds.GT.8) THEN
557  nredc=nredc+1
558  RETURN
559  ENDIF
560  kfm=2*iabs(k(i,2))
561  IF(k(i,2).LT.0) kfm=kfm-1
562  DO 530 ids=nds-1,1,-1
563  iin=ids+1
564  IF(kfm.LT.kfdm(ids)) goto 540
565  530 kfdm(ids+1)=kfdm(ids)
566  iin=1
567  540 kfdm(iin)=kfm
568  550 CONTINUE
569 
570 C...Find whether old or new final state.
571  DO 570 idc=1,nkfdc
572  IF(nds.LT.kfdc(idc,0)) THEN
573  ikfdc=idc
574  goto 580
575  ELSEIF(nds.EQ.kfdc(idc,0)) THEN
576  DO 560 i=1,nds
577  IF(kfdm(i).LT.kfdc(idc,i)) THEN
578  ikfdc=idc
579  goto 580
580  ELSEIF(kfdm(i).GT.kfdc(idc,i)) THEN
581  goto 570
582  ENDIF
583  560 CONTINUE
584  ikfdc=-idc
585  goto 580
586  ENDIF
587  570 CONTINUE
588  ikfdc=nkfdc+1
589  580 IF(ikfdc.LT.0) THEN
590  ikfdc=-ikfdc
591  ELSEIF(nkfdc.GE.200) THEN
592  nredc=nredc+1
593  RETURN
594  ELSE
595  DO 590 idc=nkfdc,ikfdc,-1
596  npdc(idc+1)=npdc(idc)
597  DO 590 i=0,8
598  590 kfdc(idc+1,i)=kfdc(idc,i)
599  nkfdc=nkfdc+1
600  kfdc(ikfdc,0)=nds
601  DO 600 i=1,nds
602  600 kfdc(ikfdc,i)=kfdm(i)
603  npdc(ikfdc)=0
604  ENDIF
605  npdc(ikfdc)=npdc(ikfdc)+1
606 
607 C...Write statistics on decay channels.
608  ELSEIF(mtabu.EQ.52) THEN
609  fac=1./max(1,nevdc)
610  WRITE(mstu(11),1900) nevdc
611  DO 620 idc=1,nkfdc
612  DO 610 i=1,kfdc(idc,0)
613  kfm=kfdc(idc,i)
614  kf=(kfm+1)/2
615  IF(2*kf.NE.kfm) kf=-kf
616  CALL luname(kf,chau)
617  chdc(i)=chau(1:12)
618  610 IF(chau(13:13).NE.' ') chdc(i)(12:12)='?'
619  620 WRITE(mstu(11),2000) fac*npdc(idc),(chdc(i),i=1,kfdc(idc,0))
620  IF(nredc.NE.0) WRITE(mstu(11),2100) fac*nredc
621 
622 C...Copy statistics on decay channels into /LUJETS/.
623  ELSEIF(mtabu.EQ.53) THEN
624  fac=1./max(1,nevdc)
625  DO 650 idc=1,nkfdc
626  k(idc,1)=32
627  k(idc,2)=99
628  k(idc,3)=0
629  k(idc,4)=0
630  k(idc,5)=kfdc(idc,0)
631  DO 630 j=1,5
632  p(idc,j)=0.
633  630 v(idc,j)=0.
634  DO 640 i=1,kfdc(idc,0)
635  kfm=kfdc(idc,i)
636  kf=(kfm+1)/2
637  IF(2*kf.NE.kfm) kf=-kf
638  IF(i.LE.5) p(idc,i)=kf
639  640 IF(i.GE.6) v(idc,i-5)=kf
640  650 v(idc,5)=fac*npdc(idc)
641  n=nkfdc
642  DO 660 j=1,5
643  k(n+1,j)=0
644  p(n+1,j)=0.
645  660 v(n+1,j)=0.
646  k(n+1,1)=32
647  k(n+1,2)=99
648  k(n+1,5)=nevdc
649  v(n+1,5)=fac*nredc
650  mstu(3)=1
651  ENDIF
652 
653 C...Format statements for output on unit MSTU(11) (default 6).
654  1000 FORMAT(///20x,'Event statistics - initial state'/
655  &20x,'based on an analysis of ',i6,' events'//
656  &3x,'Main flavours after',8x,'Fraction',4x,'Subfractions ',
657  &'according to fragmenting system multiplicity'/
658  &4x,'hard interaction',24x,'1',7x,'2',7x,'3',7x,'4',7x,'5',
659  &6x,'6-7',5x,'8-10',3x,'11-15',3x,'16-25',4x,'>25'/)
660  1100 FORMAT(3x,a12,1x,a12,f10.5,1x,10f8.4)
661  1200 FORMAT(///20x,'Event statistics - final state'/
662  &20x,'based on an analysis of ',i6,' events'//
663  &5x,'Mean primary multiplicity =',f8.3/
664  &5x,'Mean final multiplicity =',f8.3/
665  &5x,'Mean charged multiplicity =',f8.3//
666  &5x,'Number of particles produced per event (directly and via ',
667  &'decays/branchings)'/
668  &5x,'KF Particle/jet MDCY',8x,'Particles',9x,'Antiparticles',
669  &5x,'Total'/34x,'prim seco prim seco'/)
670  1300 FORMAT(1x,i6,4x,a16,i2,5(1x,f9.4))
671  1400 FORMAT(///20x,'Factorial moments analysis of multiplicity'/
672  &20x,'based on an analysis of ',i6,' events'//
673  &3x,'delta-',a3,' delta-phi <n>/bin',10x,'<F2>',18x,'<F3>',
674  &18x,'<F4>',18x,'<F5>'/35x,4(' value error '))
675  1500 FORMAT(10x)
676  1600 FORMAT(2x,2f10.4,f12.4,4(f12.4,f10.4))
677  1700 FORMAT(///20x,'Energy-Energy Correlation and Asymmetry'/
678  &20x,'based on an analysis of ',i6,' events'//
679  &2x,'theta range',8x,'EEC(theta)',8x,'EEC(180-theta)',7x,
680  &'EECA(theta)'/2x,'in degrees ',3(' value error')/)
681  1800 FORMAT(2x,f4.1,' - ',f4.1,3(f11.4,f9.4))
682  1900 FORMAT(///20x,'Decay channel analysis - final state'/
683  &20x,'based on an analysis of ',i6,' events'//
684  &2x,'Probability',10x,'Complete final state'/)
685  2000 FORMAT(2x,f9.5,5x,8(a12,1x))
686  2100 FORMAT(2x,f9.5,5x,'into other channels (more than 8 particles ',
687  &'or table overflow)')
688 
689  RETURN
690  END