HIJING/hipyset1_35/lutabu_hijing.F

   1 * $Id$
   2
   3 C*********************************************************************
   4
   5       SUBROUTINE LUTABU_HIJING(MTABU)
   6
   7 C...Purpose: to evaluate various properties of an event, with
   8 C...statistics accumulated during the course of the run and
   9 C...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
  27 C...Reset statistics on initial parton state.
  28       IF(MTABU.EQ.10) THEN
  29         NEVIS=0
  30         NKFIS=0
  31
  32 C...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
  69 C...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
 100 C...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
 122 C...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
 152 C...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
 160 C...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
 195 C...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
 224 C...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
 236 C...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
 263 C...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
 273 C...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
 303 C...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
 334 C...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
 371 C...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
 387 C...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
 408 C...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
 436 C...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
 447 C...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
 481 C...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
 501 C...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
 515 C...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
 544 C...Reset statistics on decay channels.
 545       ELSEIF(MTABU.EQ.50) THEN
 546         NEVDC=0
 547         NKFDC=0
 548         NREDC=0
 549
 550 C...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
 571 C...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
 608 C...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
 623 C...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
 654 C...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