Merging the VirtualMC branch to the main development branch (HEAD)
[u/mrichter/AliRoot.git] / HIJING / hipyset1_35 / lutabu_hijing.F
CommitLineData
e74335a4 1* $Id$
2
3C*********************************************************************
4
5 SUBROUTINE LUTABU_HIJING(MTABU)
6
7C...Purpose: to evaluate various properties of an event, with
8C...statistics accumulated during the course of the run and
9C...printed at the end.
10#include "lujets_hijing.inc"
11#include "ludat1_hijing.inc"
12#include "ludat2_hijing.inc"
13#include "ludat3_hijing.inc"
14 DIMENSION KFIS(100,2),NPIS(100,0:10),KFFS(400),NPFS(400,4),
15 &FEVFM(10,4),FM1FM(3,10,4),FM2FM(3,10,4),FMOMA(4),FMOMS(4),
16 &FEVEE(50),FE1EC(50),FE2EC(50),FE1EA(25),FE2EA(25),
17 &KFDM(8),KFDC(200,0:8),NPDC(200)
18 SAVE NEVIS,NKFIS,KFIS,NPIS,NEVFS,NPRFS,NFIFS,NCHFS,NKFFS,
19 &KFFS,NPFS,NEVFM,NMUFM,FM1FM,FM2FM,NEVEE,FE1EC,FE2EC,FE1EA,
20 &FE2EA,NEVDC,NKFDC,NREDC,KFDC,NPDC
21 CHARACTER CHAU*16,CHIS(2)*12,CHDC(8)*12
22 DATA NEVIS/0/,NKFIS/0/,NEVFS/0/,NPRFS/0/,NFIFS/0/,NCHFS/0/,
23 &NKFFS/0/,NEVFM/0/,NMUFM/0/,FM1FM/120*0./,FM2FM/120*0./,
24 &NEVEE/0/,FE1EC/50*0./,FE2EC/50*0./,FE1EA/25*0./,FE2EA/25*0./,
25 &NEVDC/0/,NKFDC/0/,NREDC/0/
26
27C...Reset statistics on initial parton state.
28 IF(MTABU.EQ.10) THEN
29 NEVIS=0
30 NKFIS=0
31
32C...Identify and order flavour content of initial state.
33 ELSEIF(MTABU.EQ.11) THEN
34 NEVIS=NEVIS+1
35 KFM1=2*IABS(MSTU(161))
36 IF(MSTU(161).GT.0) KFM1=KFM1-1
37 KFM2=2*IABS(MSTU(162))
38 IF(MSTU(162).GT.0) KFM2=KFM2-1
39 KFMN=MIN(KFM1,KFM2)
40 KFMX=MAX(KFM1,KFM2)
41 DO 100 I=1,NKFIS
42 IF(KFMN.EQ.KFIS(I,1).AND.KFMX.EQ.KFIS(I,2)) THEN
43 IKFIS=-I
44 GOTO 110
45 ELSEIF(KFMN.LT.KFIS(I,1).OR.(KFMN.EQ.KFIS(I,1).AND.
46 & KFMX.LT.KFIS(I,2))) THEN
47 IKFIS=I
48 GOTO 110
49 ENDIF
50 100 CONTINUE
51 IKFIS=NKFIS+1
52 110 IF(IKFIS.LT.0) THEN
53 IKFIS=-IKFIS
54 ELSE
55 IF(NKFIS.GE.100) RETURN
56 DO 120 I=NKFIS,IKFIS,-1
57 KFIS(I+1,1)=KFIS(I,1)
58 KFIS(I+1,2)=KFIS(I,2)
59 DO 120 J=0,10
60 120 NPIS(I+1,J)=NPIS(I,J)
61 NKFIS=NKFIS+1
62 KFIS(IKFIS,1)=KFMN
63 KFIS(IKFIS,2)=KFMX
64 DO 130 J=0,10
65 130 NPIS(IKFIS,J)=0
66 ENDIF
67 NPIS(IKFIS,0)=NPIS(IKFIS,0)+1
68
69C...Count number of partons in initial state.
70 NP=0
71 DO 150 I=1,N
72 IF(K(I,1).LE.0.OR.K(I,1).GT.12) THEN
73 ELSEIF(IABS(K(I,2)).GT.80.AND.IABS(K(I,2)).LE.100) THEN
74 ELSEIF(IABS(K(I,2)).GT.100.AND.MOD(IABS(K(I,2))/10,10).NE.0)
75 & THEN
76 ELSE
77 IM=I
78 140 IM=K(IM,3)
79 IF(IM.LE.0.OR.IM.GT.N) THEN
80 NP=NP+1
81 ELSEIF(K(IM,1).LE.0.OR.K(IM,1).GT.20) THEN
82 NP=NP+1
83 ELSEIF(IABS(K(IM,2)).GT.80.AND.IABS(K(IM,2)).LE.100) THEN
84 ELSEIF(IABS(K(IM,2)).GT.100.AND.MOD(IABS(K(IM,2))/10,10).NE.0)
85 & THEN
86 ELSE
87 GOTO 140
88 ENDIF
89 ENDIF
90 150 CONTINUE
91 NPCO=MAX(NP,1)
92 IF(NP.GE.6) NPCO=6
93 IF(NP.GE.8) NPCO=7
94 IF(NP.GE.11) NPCO=8
95 IF(NP.GE.16) NPCO=9
96 IF(NP.GE.26) NPCO=10
97 NPIS(IKFIS,NPCO)=NPIS(IKFIS,NPCO)+1
98 MSTU(62)=NP
99
100C...Write statistics on initial parton state.
101 ELSEIF(MTABU.EQ.12) THEN
102 FAC=1./MAX(1,NEVIS)
103 WRITE(MSTU(11),1000) NEVIS
104 DO 160 I=1,NKFIS
105 KFMN=KFIS(I,1)
106 IF(KFMN.EQ.0) KFMN=KFIS(I,2)
107 KFM1=(KFMN+1)/2
108 IF(2*KFM1.EQ.KFMN) KFM1=-KFM1
109 CALL LUNAME_HIJING(KFM1,CHAU)
110 CHIS(1)=CHAU(1:12)
111 IF(CHAU(13:13).NE.' ') CHIS(1)(12:12)='?'
112 KFMX=KFIS(I,2)
113 IF(KFIS(I,1).EQ.0) KFMX=0
114 KFM2=(KFMX+1)/2
115 IF(2*KFM2.EQ.KFMX) KFM2=-KFM2
116 CALL LUNAME_HIJING(KFM2,CHAU)
117 CHIS(2)=CHAU(1:12)
118 IF(CHAU(13:13).NE.' ') CHIS(2)(12:12)='?'
119 160 WRITE(MSTU(11),1100) CHIS(1),CHIS(2),FAC*NPIS(I,0),
120 & (NPIS(I,J)/FLOAT(NPIS(I,0)),J=1,10)
121
122C...Copy statistics on initial parton state into /LUJETS_HIJING/.
123 ELSEIF(MTABU.EQ.13) THEN
124 FAC=1./MAX(1,NEVIS)
125 DO 170 I=1,NKFIS
126 KFMN=KFIS(I,1)
127 IF(KFMN.EQ.0) KFMN=KFIS(I,2)
128 KFM1=(KFMN+1)/2
129 IF(2*KFM1.EQ.KFMN) KFM1=-KFM1
130 KFMX=KFIS(I,2)
131 IF(KFIS(I,1).EQ.0) KFMX=0
132 KFM2=(KFMX+1)/2
133 IF(2*KFM2.EQ.KFMX) KFM2=-KFM2
134 K(I,1)=32
135 K(I,2)=99
136 K(I,3)=KFM1
137 K(I,4)=KFM2
138 K(I,5)=NPIS(I,0)
139 DO 170 J=1,5
140 P(I,J)=FAC*NPIS(I,J)
141 170 V(I,J)=FAC*NPIS(I,J+5)
142 N=NKFIS
143 DO 180 J=1,5
144 K(N+1,J)=0
145 P(N+1,J)=0.
146 180 V(N+1,J)=0.
147 K(N+1,1)=32
148 K(N+1,2)=99
149 K(N+1,5)=NEVIS
150 MSTU(3)=1
151
152C...Reset statistics on number of particles/partons.
153 ELSEIF(MTABU.EQ.20) THEN
154 NEVFS=0
155 NPRFS=0
156 NFIFS=0
157 NCHFS=0
158 NKFFS=0
159
160C...Identify whether particle/parton is primary or not.
161 ELSEIF(MTABU.EQ.21) THEN
162 NEVFS=NEVFS+1
163 MSTU(62)=0
164 DO 230 I=1,N
165 IF(K(I,1).LE.0.OR.K(I,1).GT.20.OR.K(I,1).EQ.13) GOTO 230
166 MSTU(62)=MSTU(62)+1
167 KC=LUCOMP_HIJING(K(I,2))
168 MPRI=0
169 IF(K(I,3).LE.0.OR.K(I,3).GT.N) THEN
170 MPRI=1
171 ELSEIF(K(K(I,3),1).LE.0.OR.K(K(I,3),1).GT.20) THEN
172 MPRI=1
173 ELSEIF(KC.EQ.0) THEN
174 ELSEIF(K(K(I,3),1).EQ.13) THEN
175 IM=K(K(I,3),3)
176 IF(IM.LE.0.OR.IM.GT.N) THEN
177 MPRI=1
178 ELSEIF(K(IM,1).LE.0.OR.K(IM,1).GT.20) THEN
179 MPRI=1
180 ENDIF
181 ELSEIF(KCHG(KC,2).EQ.0) THEN
182 KCM=LUCOMP_HIJING(K(K(I,3),2))
183 IF(KCM.NE.0) THEN
184 IF(KCHG(KCM,2).NE.0) MPRI=1
185 ENDIF
186 ENDIF
187 IF(KC.NE.0.AND.MPRI.EQ.1) THEN
188 IF(KCHG(KC,2).EQ.0) NPRFS=NPRFS+1
189 ENDIF
190 IF(K(I,1).LE.10) THEN
191 NFIFS=NFIFS+1
192 IF(LUCHGE_HIJING(K(I,2)).NE.0) NCHFS=NCHFS+1
193 ENDIF
194
195C...Fill statistics on number of particles/partons in event.
196 KFA=IABS(K(I,2))
197 KFS=3-ISIGN(1,K(I,2))-MPRI
198 DO 190 IP=1,NKFFS
199 IF(KFA.EQ.KFFS(IP)) THEN
200 IKFFS=-IP
201 GOTO 200
202 ELSEIF(KFA.LT.KFFS(IP)) THEN
203 IKFFS=IP
204 GOTO 200
205 ENDIF
206 190 CONTINUE
207 IKFFS=NKFFS+1
208 200 IF(IKFFS.LT.0) THEN
209 IKFFS=-IKFFS
210 ELSE
211 IF(NKFFS.GE.400) RETURN
212 DO 210 IP=NKFFS,IKFFS,-1
213 KFFS(IP+1)=KFFS(IP)
214 DO 210 J=1,4
215 210 NPFS(IP+1,J)=NPFS(IP,J)
216 NKFFS=NKFFS+1
217 KFFS(IKFFS)=KFA
218 DO 220 J=1,4
219 220 NPFS(IKFFS,J)=0
220 ENDIF
221 NPFS(IKFFS,KFS)=NPFS(IKFFS,KFS)+1
222 230 CONTINUE
223
224C...Write statistics on particle/parton composition of events.
225 ELSEIF(MTABU.EQ.22) THEN
226 FAC=1./MAX(1,NEVFS)
227 WRITE(MSTU(11),1200) NEVFS,FAC*NPRFS,FAC*NFIFS,FAC*NCHFS
228 DO 240 I=1,NKFFS
229 CALL LUNAME_HIJING(KFFS(I),CHAU)
230 KC=LUCOMP_HIJING(KFFS(I))
231 MDCYF=0
232 IF(KC.NE.0) MDCYF=MDCY(KC,1)
233 240 WRITE(MSTU(11),1300) KFFS(I),CHAU,MDCYF,(FAC*NPFS(I,J),J=1,4),
234 & FAC*(NPFS(I,1)+NPFS(I,2)+NPFS(I,3)+NPFS(I,4))
235
236C...Copy particle/parton composition information into /LUJETS_HIJING/.
237 ELSEIF(MTABU.EQ.23) THEN
238 FAC=1./MAX(1,NEVFS)
239 DO 260 I=1,NKFFS
240 K(I,1)=32
241 K(I,2)=99
242 K(I,3)=KFFS(I)
243 K(I,4)=0
244 K(I,5)=NPFS(I,1)+NPFS(I,2)+NPFS(I,3)+NPFS(I,4)
245 DO 250 J=1,4
246 P(I,J)=FAC*NPFS(I,J)
247 250 V(I,J)=0.
248 P(I,5)=FAC*K(I,5)
249 260 V(I,5)=0.
250 N=NKFFS
251 DO 270 J=1,5
252 K(N+1,J)=0
253 P(N+1,J)=0.
254 270 V(N+1,J)=0.
255 K(N+1,1)=32
256 K(N+1,2)=99
257 K(N+1,5)=NEVFS
258 P(N+1,1)=FAC*NPRFS
259 P(N+1,2)=FAC*NFIFS
260 P(N+1,3)=FAC*NCHFS
261 MSTU(3)=1
262
263C...Reset factorial moments statistics.
264 ELSEIF(MTABU.EQ.30) THEN
265 NEVFM=0
266 NMUFM=0
267 DO 280 IM=1,3
268 DO 280 IB=1,10
269 DO 280 IP=1,4
270 FM1FM(IM,IB,IP)=0.
271 280 FM2FM(IM,IB,IP)=0.
272
273C...Find particles to include, with (pion,pseudo)rapidity and azimuth.
274 ELSEIF(MTABU.EQ.31) THEN
275 NEVFM=NEVFM+1
276 NLOW=N+MSTU(3)
277 NUPP=NLOW
278 DO 360 I=1,N
279 IF(K(I,1).LE.0.OR.K(I,1).GT.10) GOTO 360
280 IF(MSTU(41).GE.2) THEN
281 KC=LUCOMP_HIJING(K(I,2))
282 IF(KC.EQ.0.OR.KC.EQ.12.OR.KC.EQ.14.OR.KC.EQ.16.OR.
283 & KC.EQ.18) GOTO 360
284 IF(MSTU(41).GE.3.AND.KCHG(KC,2).EQ.0.AND.LUCHGE_HIJING(K(I,2))
285 $ .EQ.0)GOTO 360
286 ENDIF
287 PMR=0.
288 IF(MSTU(42).EQ.1.AND.K(I,2).NE.22) PMR=ULMASS_HIJING(211)
289 IF(MSTU(42).GE.2) PMR=P(I,5)
290 PR=MAX(1E-20,PMR**2+P(I,1)**2+P(I,2)**2)
291 YETA=SIGN(LOG(MIN((SQRT(PR+P(I,3)**2)+ABS(P(I,3)))/SQRT(PR),
292 & 1E20)),P(I,3))
293 IF(ABS(YETA).GT.PARU(57)) GOTO 360
294 PHI=ULANGL_HIJING(P(I,1),P(I,2))
295 IYETA=512.*(YETA+PARU(57))/(2.*PARU(57))
296 IYETA=MAX(0,MIN(511,IYETA))
297 IPHI=512.*(PHI+PARU(1))/PARU(2)
298 IPHI=MAX(0,MIN(511,IPHI))
299 IYEP=0
300 DO 290 IB=0,9
301 290 IYEP=IYEP+4**IB*(2*MOD(IYETA/2**IB,2)+MOD(IPHI/2**IB,2))
302
303C...Order particles in (pseudo)rapidity and/or azimuth.
304 IF(NUPP.GT.MSTU(4)-5-MSTU(32)) THEN
305 CALL LUERRM_HIJING(11
306 $ ,'(LUTABU_HIJING:) no more memory left in LUJETS_HIJING'
307 $ )
308 RETURN
309 ENDIF
310 NUPP=NUPP+1
311 IF(NUPP.EQ.NLOW+1) THEN
312 K(NUPP,1)=IYETA
313 K(NUPP,2)=IPHI
314 K(NUPP,3)=IYEP
315 ELSE
316 DO 300 I1=NUPP-1,NLOW+1,-1
317 IF(IYETA.GE.K(I1,1)) GOTO 310
318 300 K(I1+1,1)=K(I1,1)
319 310 K(I1+1,1)=IYETA
320 DO 320 I1=NUPP-1,NLOW+1,-1
321 IF(IPHI.GE.K(I1,2)) GOTO 330
322 320 K(I1+1,2)=K(I1,2)
323 330 K(I1+1,2)=IPHI
324 DO 340 I1=NUPP-1,NLOW+1,-1
325 IF(IYEP.GE.K(I1,3)) GOTO 350
326 340 K(I1+1,3)=K(I1,3)
327 350 K(I1+1,3)=IYEP
328 ENDIF
329 360 CONTINUE
330 K(NUPP+1,1)=2**10
331 K(NUPP+1,2)=2**10
332 K(NUPP+1,3)=4**10
333
334C...Calculate sum of factorial moments in event.
335 DO 400 IM=1,3
336 DO 370 IB=1,10
337 DO 370 IP=1,4
338 370 FEVFM(IB,IP)=0.
339 DO 380 IB=1,10
340 IF(IM.LE.2) IBIN=2**(10-IB)
341 IF(IM.EQ.3) IBIN=4**(10-IB)
342 IAGR=K(NLOW+1,IM)/IBIN
343 NAGR=1
344 DO 380 I=NLOW+2,NUPP+1
345 ICUT=K(I,IM)/IBIN
346 IF(ICUT.EQ.IAGR) THEN
347 NAGR=NAGR+1
348 ELSE
349 IF(NAGR.EQ.1) THEN
350 ELSEIF(NAGR.EQ.2) THEN
351 FEVFM(IB,1)=FEVFM(IB,1)+2.
352 ELSEIF(NAGR.EQ.3) THEN
353 FEVFM(IB,1)=FEVFM(IB,1)+6.
354 FEVFM(IB,2)=FEVFM(IB,2)+6.
355 ELSEIF(NAGR.EQ.4) THEN
356 FEVFM(IB,1)=FEVFM(IB,1)+12.
357 FEVFM(IB,2)=FEVFM(IB,2)+24.
358 FEVFM(IB,3)=FEVFM(IB,3)+24.
359 ELSE
360 FEVFM(IB,1)=FEVFM(IB,1)+NAGR*(NAGR-1.)
361 FEVFM(IB,2)=FEVFM(IB,2)+NAGR*(NAGR-1.)*(NAGR-2.)
362 FEVFM(IB,3)=FEVFM(IB,3)+NAGR*(NAGR-1.)*(NAGR-2.)*(NAGR-3.)
363 FEVFM(IB,4)=FEVFM(IB,4)+NAGR*(NAGR-1.)*(NAGR-2.)*(NAGR-3.)*
364 & (NAGR-4.)
365 ENDIF
366 IAGR=ICUT
367 NAGR=1
368 ENDIF
369 380 CONTINUE
370
371C...Add results to total statistics.
372 DO 390 IB=10,1,-1
373 DO 390 IP=1,4
374 IF(FEVFM(1,IP).LT.0.5) THEN
375 FEVFM(IB,IP)=0.
376 ELSEIF(IM.LE.2) THEN
377 FEVFM(IB,IP)=2**((IB-1)*IP)*FEVFM(IB,IP)/FEVFM(1,IP)
378 ELSE
379 FEVFM(IB,IP)=4**((IB-1)*IP)*FEVFM(IB,IP)/FEVFM(1,IP)
380 ENDIF
381 FM1FM(IM,IB,IP)=FM1FM(IM,IB,IP)+FEVFM(IB,IP)
382 390 FM2FM(IM,IB,IP)=FM2FM(IM,IB,IP)+FEVFM(IB,IP)**2
383 400 CONTINUE
384 NMUFM=NMUFM+(NUPP-NLOW)
385 MSTU(62)=NUPP-NLOW
386
387C...Write accumulated statistics on factorial moments.
388 ELSEIF(MTABU.EQ.32) THEN
389 FAC=1./MAX(1,NEVFM)
390 IF(MSTU(42).LE.0) WRITE(MSTU(11),1400) NEVFM,'eta'
391 IF(MSTU(42).EQ.1) WRITE(MSTU(11),1400) NEVFM,'ypi'
392 IF(MSTU(42).GE.2) WRITE(MSTU(11),1400) NEVFM,'y '
393 DO 420 IM=1,3
394 WRITE(MSTU(11),1500)
395 DO 420 IB=1,10
396 BYETA=2.*PARU(57)
397 IF(IM.NE.2) BYETA=BYETA/2**(IB-1)
398 BPHI=PARU(2)
399 IF(IM.NE.1) BPHI=BPHI/2**(IB-1)
400 IF(IM.LE.2) BNAVE=FAC*NMUFM/FLOAT(2**(IB-1))
401 IF(IM.EQ.3) BNAVE=FAC*NMUFM/FLOAT(4**(IB-1))
402 DO 410 IP=1,4
403 FMOMA(IP)=FAC*FM1FM(IM,IB,IP)
404 410 FMOMS(IP)=SQRT(MAX(0.,FAC*(FAC*FM2FM(IM,IB,IP)-FMOMA(IP)**2)))
405 420 WRITE(MSTU(11),1600) BYETA,BPHI,BNAVE,(FMOMA(IP),FMOMS(IP),
406 & IP=1,4)
407
408C...Copy statistics on factorial moments into /LUJETS_HIJING/.
409 ELSEIF(MTABU.EQ.33) THEN
410 FAC=1./MAX(1,NEVFM)
411 DO 430 IM=1,3
412 DO 430 IB=1,10
413 I=10*(IM-1)+IB
414 K(I,1)=32
415 K(I,2)=99
416 K(I,3)=1
417 IF(IM.NE.2) K(I,3)=2**(IB-1)
418 K(I,4)=1
419 IF(IM.NE.1) K(I,4)=2**(IB-1)
420 K(I,5)=0
421 P(I,1)=2.*PARU(57)/K(I,3)
422 V(I,1)=PARU(2)/K(I,4)
423 DO 430 IP=1,4
424 P(I,IP+1)=FAC*FM1FM(IM,IB,IP)
425 430 V(I,IP+1)=SQRT(MAX(0.,FAC*(FAC*FM2FM(IM,IB,IP)-P(I,IP+1)**2)))
426 N=30
427 DO 440 J=1,5
428 K(N+1,J)=0
429 P(N+1,J)=0.
430 440 V(N+1,J)=0.
431 K(N+1,1)=32
432 K(N+1,2)=99
433 K(N+1,5)=NEVFM
434 MSTU(3)=1
435
436C...Reset statistics on Energy-Energy Correlation.
437 ELSEIF(MTABU.EQ.40) THEN
438 NEVEE=0
439 DO 450 J=1,25
440 FE1EC(J)=0.
441 FE2EC(J)=0.
442 FE1EC(51-J)=0.
443 FE2EC(51-J)=0.
444 FE1EA(J)=0.
445 450 FE2EA(J)=0.
446
447C...Find particles to include, with proper assumed mass.
448 ELSEIF(MTABU.EQ.41) THEN
449 NEVEE=NEVEE+1
450 NLOW=N+MSTU(3)
451 NUPP=NLOW
452 ECM=0.
453 DO 460 I=1,N
454 IF(K(I,1).LE.0.OR.K(I,1).GT.10) GOTO 460
455 IF(MSTU(41).GE.2) THEN
456 KC=LUCOMP_HIJING(K(I,2))
457 IF(KC.EQ.0.OR.KC.EQ.12.OR.KC.EQ.14.OR.KC.EQ.16.OR.
458 & KC.EQ.18) GOTO 460
459 IF(MSTU(41).GE.3.AND.KCHG(KC,2).EQ.0.AND.LUCHGE_HIJING(K(I,2))
460 $ .EQ.0)GOTO 460
461 ENDIF
462 PMR=0.
463 IF(MSTU(42).EQ.1.AND.K(I,2).NE.22) PMR=ULMASS_HIJING(211)
464 IF(MSTU(42).GE.2) PMR=P(I,5)
465 IF(NUPP.GT.MSTU(4)-5-MSTU(32)) THEN
466 CALL LUERRM_HIJING(11
467 $ ,'(LUTABU_HIJING:) no more memory left in LUJETS_HIJING'
468 $ )
469 RETURN
470 ENDIF
471 NUPP=NUPP+1
472 P(NUPP,1)=P(I,1)
473 P(NUPP,2)=P(I,2)
474 P(NUPP,3)=P(I,3)
475 P(NUPP,4)=SQRT(PMR**2+P(I,1)**2+P(I,2)**2+P(I,3)**2)
476 P(NUPP,5)=MAX(1E-10,SQRT(P(I,1)**2+P(I,2)**2+P(I,3)**2))
477 ECM=ECM+P(NUPP,4)
478 460 CONTINUE
479 IF(NUPP.EQ.NLOW) RETURN
480
481C...Analyze Energy-Energy Correlation in event.
482 FAC=(2./ECM**2)*50./PARU(1)
483 DO 470 J=1,50
484 470 FEVEE(J)=0.
485 DO 480 I1=NLOW+2,NUPP
486 DO 480 I2=NLOW+1,I1-1
487 CTHE=(P(I1,1)*P(I2,1)+P(I1,2)*P(I2,2)+P(I1,3)*P(I2,3))/
488 & (P(I1,5)*P(I2,5))
489 THE=ACOS(MAX(-1.,MIN(1.,CTHE)))
490 ITHE=MAX(1,MIN(50,1+INT(50.*THE/PARU(1))))
491 480 FEVEE(ITHE)=FEVEE(ITHE)+FAC*P(I1,4)*P(I2,4)
492 DO 490 J=1,25
493 FE1EC(J)=FE1EC(J)+FEVEE(J)
494 FE2EC(J)=FE2EC(J)+FEVEE(J)**2
495 FE1EC(51-J)=FE1EC(51-J)+FEVEE(51-J)
496 FE2EC(51-J)=FE2EC(51-J)+FEVEE(51-J)**2
497 FE1EA(J)=FE1EA(J)+(FEVEE(51-J)-FEVEE(J))
498 490 FE2EA(J)=FE2EA(J)+(FEVEE(51-J)-FEVEE(J))**2
499 MSTU(62)=NUPP-NLOW
500
501C...Write statistics on Energy-Energy Correlation.
502 ELSEIF(MTABU.EQ.42) THEN
503 FAC=1./MAX(1,NEVEE)
504 WRITE(MSTU(11),1700) NEVEE
505 DO 500 J=1,25
506 FEEC1=FAC*FE1EC(J)
507 FEES1=SQRT(MAX(0.,FAC*(FAC*FE2EC(J)-FEEC1**2)))
508 FEEC2=FAC*FE1EC(51-J)
509 FEES2=SQRT(MAX(0.,FAC*(FAC*FE2EC(51-J)-FEEC2**2)))
510 FEECA=FAC*FE1EA(J)
511 FEESA=SQRT(MAX(0.,FAC*(FAC*FE2EA(J)-FEECA**2)))
512 500 WRITE(MSTU(11),1800) 3.6*(J-1),3.6*J,FEEC1,FEES1,FEEC2,FEES2,
513 & FEECA,FEESA
514
515C...Copy statistics on Energy-Energy Correlation into /LUJETS_HIJING/.
516 ELSEIF(MTABU.EQ.43) THEN
517 FAC=1./MAX(1,NEVEE)
518 DO 510 I=1,25
519 K(I,1)=32
520 K(I,2)=99
521 K(I,3)=0
522 K(I,4)=0
523 K(I,5)=0
524 P(I,1)=FAC*FE1EC(I)
525 V(I,1)=SQRT(MAX(0.,FAC*(FAC*FE2EC(I)-P(I,1)**2)))
526 P(I,2)=FAC*FE1EC(51-I)
527 V(I,2)=SQRT(MAX(0.,FAC*(FAC*FE2EC(51-I)-P(I,2)**2)))
528 P(I,3)=FAC*FE1EA(I)
529 V(I,3)=SQRT(MAX(0.,FAC*(FAC*FE2EA(I)-P(I,3)**2)))
530 P(I,4)=PARU(1)*(I-1)/50.
531 P(I,5)=PARU(1)*I/50.
532 V(I,4)=3.6*(I-1)
533 510 V(I,5)=3.6*I
534 N=25
535 DO 520 J=1,5
536 K(N+1,J)=0
537 P(N+1,J)=0.
538 520 V(N+1,J)=0.
539 K(N+1,1)=32
540 K(N+1,2)=99
541 K(N+1,5)=NEVEE
542 MSTU(3)=1
543
544C...Reset statistics on decay channels.
545 ELSEIF(MTABU.EQ.50) THEN
546 NEVDC=0
547 NKFDC=0
548 NREDC=0
549
550C...Identify and order flavour content of final state.
551 ELSEIF(MTABU.EQ.51) THEN
552 NEVDC=NEVDC+1
553 NDS=0
554 DO 550 I=1,N
555 IF(K(I,1).LE.0.OR.K(I,1).GE.6) GOTO 550
556 NDS=NDS+1
557 IF(NDS.GT.8) THEN
558 NREDC=NREDC+1
559 RETURN
560 ENDIF
561 KFM=2*IABS(K(I,2))
562 IF(K(I,2).LT.0) KFM=KFM-1
563 DO 530 IDS=NDS-1,1,-1
564 IIN=IDS+1
565 IF(KFM.LT.KFDM(IDS)) GOTO 540
566 530 KFDM(IDS+1)=KFDM(IDS)
567 IIN=1
568 540 KFDM(IIN)=KFM
569 550 CONTINUE
570
571C...Find whether old or new final state.
572 DO 570 IDC=1,NKFDC
573 IF(NDS.LT.KFDC(IDC,0)) THEN
574 IKFDC=IDC
575 GOTO 580
576 ELSEIF(NDS.EQ.KFDC(IDC,0)) THEN
577 DO 560 I=1,NDS
578 IF(KFDM(I).LT.KFDC(IDC,I)) THEN
579 IKFDC=IDC
580 GOTO 580
581 ELSEIF(KFDM(I).GT.KFDC(IDC,I)) THEN
582 GOTO 570
583 ENDIF
584 560 CONTINUE
585 IKFDC=-IDC
586 GOTO 580
587 ENDIF
588 570 CONTINUE
589 IKFDC=NKFDC+1
590 580 IF(IKFDC.LT.0) THEN
591 IKFDC=-IKFDC
592 ELSEIF(NKFDC.GE.200) THEN
593 NREDC=NREDC+1
594 RETURN
595 ELSE
596 DO 590 IDC=NKFDC,IKFDC,-1
597 NPDC(IDC+1)=NPDC(IDC)
598 DO 590 I=0,8
599 590 KFDC(IDC+1,I)=KFDC(IDC,I)
600 NKFDC=NKFDC+1
601 KFDC(IKFDC,0)=NDS
602 DO 600 I=1,NDS
603 600 KFDC(IKFDC,I)=KFDM(I)
604 NPDC(IKFDC)=0
605 ENDIF
606 NPDC(IKFDC)=NPDC(IKFDC)+1
607
608C...Write statistics on decay channels.
609 ELSEIF(MTABU.EQ.52) THEN
610 FAC=1./MAX(1,NEVDC)
611 WRITE(MSTU(11),1900) NEVDC
612 DO 620 IDC=1,NKFDC
613 DO 610 I=1,KFDC(IDC,0)
614 KFM=KFDC(IDC,I)
615 KF=(KFM+1)/2
616 IF(2*KF.NE.KFM) KF=-KF
617 CALL LUNAME_HIJING(KF,CHAU)
618 CHDC(I)=CHAU(1:12)
619 610 IF(CHAU(13:13).NE.' ') CHDC(I)(12:12)='?'
620 620 WRITE(MSTU(11),2000) FAC*NPDC(IDC),(CHDC(I),I=1,KFDC(IDC,0))
621 IF(NREDC.NE.0) WRITE(MSTU(11),2100) FAC*NREDC
622
623C...Copy statistics on decay channels into /LUJETS_HIJING/.
624 ELSEIF(MTABU.EQ.53) THEN
625 FAC=1./MAX(1,NEVDC)
626 DO 650 IDC=1,NKFDC
627 K(IDC,1)=32
628 K(IDC,2)=99
629 K(IDC,3)=0
630 K(IDC,4)=0
631 K(IDC,5)=KFDC(IDC,0)
632 DO 630 J=1,5
633 P(IDC,J)=0.
634 630 V(IDC,J)=0.
635 DO 640 I=1,KFDC(IDC,0)
636 KFM=KFDC(IDC,I)
637 KF=(KFM+1)/2
638 IF(2*KF.NE.KFM) KF=-KF
639 IF(I.LE.5) P(IDC,I)=KF
640 640 IF(I.GE.6) V(IDC,I-5)=KF
641 650 V(IDC,5)=FAC*NPDC(IDC)
642 N=NKFDC
643 DO 660 J=1,5
644 K(N+1,J)=0
645 P(N+1,J)=0.
646 660 V(N+1,J)=0.
647 K(N+1,1)=32
648 K(N+1,2)=99
649 K(N+1,5)=NEVDC
650 V(N+1,5)=FAC*NREDC
651 MSTU(3)=1
652 ENDIF
653
654C...Format statements for output on unit MSTU(11) (default 6).
655 1000 FORMAT(///20X,'Event statistics - initial state'/
656 &20X,'based on an analysis of ',I6,' events'//
657 &3X,'Main flavours after',8X,'Fraction',4X,'Subfractions ',
658 &'according to fragmenting system multiplicity'/
659 &4X,'hard interaction',24X,'1',7X,'2',7X,'3',7X,'4',7X,'5',
660 &6X,'6-7',5X,'8-10',3X,'11-15',3X,'16-25',4X,'>25'/)
661 1100 FORMAT(3X,A12,1X,A12,F10.5,1X,10F8.4)
662 1200 FORMAT(///20X,'Event statistics - final state'/
663 &20X,'based on an analysis of ',I6,' events'//
664 &5X,'Mean primary multiplicity =',F8.3/
665 &5X,'Mean final multiplicity =',F8.3/
666 &5X,'Mean charged multiplicity =',F8.3//
667 &5X,'Number of particles produced per event (directly and via ',
668 &'decays/branchings)'/
669 &5X,'KF Particle/jet MDCY',8X,'Particles',9X,'Antiparticles',
670 &5X,'Total'/34X,'prim seco prim seco'/)
671 1300 FORMAT(1X,I6,4X,A16,I2,5(1X,F9.4))
672 1400 FORMAT(///20X,'Factorial moments analysis of multiplicity'/
673 &20X,'based on an analysis of ',I6,' events'//
674 &3X,'delta-',A3,' delta-phi <n>/bin',10X,'<F2>',18X,'<F3>',
675 &18X,'<F4>',18X,'<F5>'/35X,4(' value error '))
676 1500 FORMAT(10X)
677 1600 FORMAT(2X,2F10.4,F12.4,4(F12.4,F10.4))
678 1700 FORMAT(///20X,'Energy-Energy Correlation and Asymmetry'/
679 &20X,'based on an analysis of ',I6,' events'//
680 &2X,'theta range',8X,'EEC(theta)',8X,'EEC(180-theta)',7X,
681 &'EECA(theta)'/2X,'in degrees ',3(' value error')/)
682 1800 FORMAT(2X,F4.1,' - ',F4.1,3(F11.4,F9.4))
683 1900 FORMAT(///20X,'Decay channel analysis - final state'/
684 &20X,'based on an analysis of ',I6,' events'//
685 &2X,'Probability',10X,'Complete final state'/)
686 2000 FORMAT(2X,F9.5,5X,8(A12,1X))
687 2100 FORMAT(2X,F9.5,5X,'into other channels (more than 8 particles ',
688 &'or table overflow)')
689
690 RETURN
691 END