pythiaeRHIC/pythia/pytabu.f

0001
0002 C*********************************************************************
0003
0004 C...PYTABU
0005 C...Evaluates various properties of an event, with statistics
0006 C...accumulated during the course of the run and
0007 C...printed at the end.
0008
0009       SUBROUTINE PYTABU(MTABU)
0010
0011 C...Double precision and integer declarations.
0012       IMPLICIT DOUBLE PRECISION(A-H, O-Z)
0013       IMPLICIT INTEGER(I-N)
0014       INTEGER PYK,PYCHGE,PYCOMP
0015 C...Parameter statement to help give large particle numbers.
0016       PARAMETER (KSUSY1=1000000,KSUSY2=2000000,KTECHN=3000000,
0017      &KEXCIT=4000000,KDIMEN=5000000)
0018 C...Commonblocks.
0019       COMMON/PYJETS/N,NPAD,K(4000,5),P(4000,5),V(4000,5)
0020       COMMON/PYDAT1/MSTU(200),PARU(200),MSTJ(200),PARJ(200)
0021       COMMON/PYDAT2/KCHG(500,4),PMAS(500,4),PARF(2000),VCKM(4,4)
0022       COMMON/PYDAT3/MDCY(500,3),MDME(8000,2),BRAT(8000),KFDP(8000,5)
0023       SAVE /PYJETS/,/PYDAT1/,/PYDAT2/,/PYDAT3/
0024 C...Local arrays, character variables, saved variables and data.
0025       DIMENSION KFIS(100,2),NPIS(100,0:10),KFFS(400),NPFS(400,4),
0026      &FEVFM(10,4),FM1FM(3,10,4),FM2FM(3,10,4),FMOMA(4),FMOMS(4),
0027      &FEVEE(50),FE1EC(50),FE2EC(50),FE1EA(25),FE2EA(25),
0028      &KFDM(8),KFDC(200,0:8),NPDC(200)
0029       SAVE NEVIS,NKFIS,KFIS,NPIS,NEVFS,NPRFS,NFIFS,NCHFS,NKFFS,
0030      &KFFS,NPFS,NEVFM,NMUFM,FM1FM,FM2FM,NEVEE,FE1EC,FE2EC,FE1EA,
0031      &FE2EA,NEVDC,NKFDC,NREDC,KFDC,NPDC
0032       CHARACTER CHAU*16,CHIS(2)*12,CHDC(8)*12
0033       DATA NEVIS/0/,NKFIS/0/,NEVFS/0/,NPRFS/0/,NFIFS/0/,NCHFS/0/,
0034      &NKFFS/0/,NEVFM/0/,NMUFM/0/,FM1FM/120*0D0/,FM2FM/120*0D0/,
0035      &NEVEE/0/,FE1EC/50*0D0/,FE2EC/50*0D0/,FE1EA/25*0D0/,FE2EA/25*0D0/,
0036      &NEVDC/0/,NKFDC/0/,NREDC/0/
0037
0038 C...Reset statistics on initial parton state.
0039       IF(MTABU.EQ.10) THEN
0040         NEVIS=0
0041         NKFIS=0
0042
0043 C...Identify and order flavour content of initial state.
0044       ELSEIF(MTABU.EQ.11) THEN
0045         NEVIS=NEVIS+1
0046         KFM1=2*IABS(MSTU(161))
0047         IF(MSTU(161).GT.0) KFM1=KFM1-1
0048         KFM2=2*IABS(MSTU(162))
0049         IF(MSTU(162).GT.0) KFM2=KFM2-1
0050         KFMN=MIN(KFM1,KFM2)
0051         KFMX=MAX(KFM1,KFM2)
0052         DO 100 I=1,NKFIS
0053           IF(KFMN.EQ.KFIS(I,1).AND.KFMX.EQ.KFIS(I,2)) THEN
0054             IKFIS=-I
0055             GOTO 110
0056           ELSEIF(KFMN.LT.KFIS(I,1).OR.(KFMN.EQ.KFIS(I,1).AND.
0057      &      KFMX.LT.KFIS(I,2))) THEN
0058             IKFIS=I
0059             GOTO 110
0060           ENDIF
0061   100   CONTINUE
0062         IKFIS=NKFIS+1
0063   110   IF(IKFIS.LT.0) THEN
0064           IKFIS=-IKFIS
0065         ELSE
0066           IF(NKFIS.GE.100) RETURN
0067           DO 130 I=NKFIS,IKFIS,-1
0068             KFIS(I+1,1)=KFIS(I,1)
0069             KFIS(I+1,2)=KFIS(I,2)
0070             DO 120 J=0,10
0071               NPIS(I+1,J)=NPIS(I,J)
0072   120       CONTINUE
0073   130     CONTINUE
0074           NKFIS=NKFIS+1
0075           KFIS(IKFIS,1)=KFMN
0076           KFIS(IKFIS,2)=KFMX
0077           DO 140 J=0,10
0078             NPIS(IKFIS,J)=0
0079   140     CONTINUE
0080         ENDIF
0081         NPIS(IKFIS,0)=NPIS(IKFIS,0)+1
0082
0083 C...Count number of partons in initial state.
0084         NP=0
0085         DO 160 I=1,N
0086           IF(K(I,1).LE.0.OR.K(I,1).GT.12) THEN
0087           ELSEIF(IABS(K(I,2)).GT.80.AND.IABS(K(I,2)).LE.100) THEN
0088           ELSEIF(IABS(K(I,2)).GT.100.AND.MOD(IABS(K(I,2))/10,10).NE.0)
0089      &      THEN
0090           ELSE
0091             IM=I
0092   150       IM=K(IM,3)
0093             IF(IM.LE.0.OR.IM.GT.N) THEN
0094               NP=NP+1
0095             ELSEIF(K(IM,1).LE.0.OR.K(IM,1).GT.20) THEN
0096               NP=NP+1
0097             ELSEIF(IABS(K(IM,2)).GT.80.AND.IABS(K(IM,2)).LE.100) THEN
0098             ELSEIF(IABS(K(IM,2)).GT.100.AND.MOD(IABS(K(IM,2))/10,10)
0099      &        .NE.0) THEN
0100             ELSE
0101               GOTO 150
0102             ENDIF
0103           ENDIF
0104   160   CONTINUE
0105         NPCO=MAX(NP,1)
0106         IF(NP.GE.6) NPCO=6
0107         IF(NP.GE.8) NPCO=7
0108         IF(NP.GE.11) NPCO=8
0109         IF(NP.GE.16) NPCO=9
0110         IF(NP.GE.26) NPCO=10
0111         NPIS(IKFIS,NPCO)=NPIS(IKFIS,NPCO)+1
0112         MSTU(62)=NP
0113
0114 C...Write statistics on initial parton state.
0115       ELSEIF(MTABU.EQ.12) THEN
0116         FAC=1D0/MAX(1,NEVIS)
0117         WRITE(MSTU(11),5000) NEVIS
0118         DO 170 I=1,NKFIS
0119           KFMN=KFIS(I,1)
0120           IF(KFMN.EQ.0) KFMN=KFIS(I,2)
0121           KFM1=(KFMN+1)/2
0122           IF(2*KFM1.EQ.KFMN) KFM1=-KFM1
0123           CALL PYNAME(KFM1,CHAU)
0124           CHIS(1)=CHAU(1:12)
0125           IF(CHAU(13:13).NE.' ') CHIS(1)(12:12)='?'
0126           KFMX=KFIS(I,2)
0127           IF(KFIS(I,1).EQ.0) KFMX=0
0128           KFM2=(KFMX+1)/2
0129           IF(2*KFM2.EQ.KFMX) KFM2=-KFM2
0130           CALL PYNAME(KFM2,CHAU)
0131           CHIS(2)=CHAU(1:12)
0132           IF(CHAU(13:13).NE.' ') CHIS(2)(12:12)='?'
0133           WRITE(MSTU(11),5100) CHIS(1),CHIS(2),FAC*NPIS(I,0),
0134      &    (NPIS(I,J)/DBLE(NPIS(I,0)),J=1,10)
0135   170   CONTINUE
0136
0137 C...Copy statistics on initial parton state into /PYJETS/.
0138       ELSEIF(MTABU.EQ.13) THEN
0139         FAC=1D0/MAX(1,NEVIS)
0140         DO 190 I=1,NKFIS
0141           KFMN=KFIS(I,1)
0142           IF(KFMN.EQ.0) KFMN=KFIS(I,2)
0143           KFM1=(KFMN+1)/2
0144           IF(2*KFM1.EQ.KFMN) KFM1=-KFM1
0145           KFMX=KFIS(I,2)
0146           IF(KFIS(I,1).EQ.0) KFMX=0
0147           KFM2=(KFMX+1)/2
0148           IF(2*KFM2.EQ.KFMX) KFM2=-KFM2
0149           K(I,1)=32
0150           K(I,2)=99
0151           K(I,3)=KFM1
0152           K(I,4)=KFM2
0153           K(I,5)=NPIS(I,0)
0154           DO 180 J=1,5
0155             P(I,J)=FAC*NPIS(I,J)
0156             V(I,J)=FAC*NPIS(I,J+5)
0157   180     CONTINUE
0158   190   CONTINUE
0159         N=NKFIS
0160         DO 200 J=1,5
0161           K(N+1,J)=0
0162           P(N+1,J)=0D0
0163           V(N+1,J)=0D0
0164   200   CONTINUE
0165         K(N+1,1)=32
0166         K(N+1,2)=99
0167         K(N+1,5)=NEVIS
0168         MSTU(3)=1
0169
0170 C...Reset statistics on number of particles/partons.
0171       ELSEIF(MTABU.EQ.20) THEN
0172         NEVFS=0
0173         NPRFS=0
0174         NFIFS=0
0175         NCHFS=0
0176         NKFFS=0
0177
0178 C...Identify whether particle/parton is primary or not.
0179       ELSEIF(MTABU.EQ.21) THEN
0180         NEVFS=NEVFS+1
0181         MSTU(62)=0
0182         DO 260 I=1,N
0183           IF(K(I,1).LE.0.OR.K(I,1).GT.20.OR.K(I,1).EQ.13) GOTO 260
0184           MSTU(62)=MSTU(62)+1
0185           KC=PYCOMP(K(I,2))
0186           MPRI=0
0187           IF(K(I,3).LE.0.OR.K(I,3).GT.N) THEN
0188             MPRI=1
0189           ELSEIF(K(K(I,3),1).LE.0.OR.K(K(I,3),1).GT.20) THEN
0190             MPRI=1
0191           ELSEIF(K(K(I,3),2).GE.91.AND.K(K(I,3),2).LE.93) THEN
0192             MPRI=1
0193           ELSEIF(KC.EQ.0) THEN
0194           ELSEIF(K(K(I,3),1).EQ.13) THEN
0195             IM=K(K(I,3),3)
0196             IF(IM.LE.0.OR.IM.GT.N) THEN
0197               MPRI=1
0198             ELSEIF(K(IM,1).LE.0.OR.K(IM,1).GT.20) THEN
0199               MPRI=1
0200             ENDIF
0201           ELSEIF(KCHG(KC,2).EQ.0) THEN
0202             KCM=PYCOMP(K(K(I,3),2))
0203             IF(KCM.NE.0) THEN
0204               IF(KCHG(KCM,2).NE.0) MPRI=1
0205             ENDIF
0206           ENDIF
0207           IF(KC.NE.0.AND.MPRI.EQ.1) THEN
0208             IF(KCHG(KC,2).EQ.0) NPRFS=NPRFS+1
0209           ENDIF
0210           IF(K(I,1).LE.10) THEN
0211             NFIFS=NFIFS+1
0212             IF(PYCHGE(K(I,2)).NE.0) NCHFS=NCHFS+1
0213           ENDIF
0214
0215 C...Fill statistics on number of particles/partons in event.
0216           KFA=IABS(K(I,2))
0217           KFS=3-ISIGN(1,K(I,2))-MPRI
0218           DO 210 IP=1,NKFFS
0219             IF(KFA.EQ.KFFS(IP)) THEN
0220               IKFFS=-IP
0221               GOTO 220
0222             ELSEIF(KFA.LT.KFFS(IP)) THEN
0223               IKFFS=IP
0224               GOTO 220
0225             ENDIF
0226   210     CONTINUE
0227           IKFFS=NKFFS+1
0228   220     IF(IKFFS.LT.0) THEN
0229             IKFFS=-IKFFS
0230           ELSE
0231             IF(NKFFS.GE.400) RETURN
0232             DO 240 IP=NKFFS,IKFFS,-1
0233               KFFS(IP+1)=KFFS(IP)
0234               DO 230 J=1,4
0235                 NPFS(IP+1,J)=NPFS(IP,J)
0236   230         CONTINUE
0237   240       CONTINUE
0238             NKFFS=NKFFS+1
0239             KFFS(IKFFS)=KFA
0240             DO 250 J=1,4
0241               NPFS(IKFFS,J)=0
0242   250       CONTINUE
0243           ENDIF
0244           NPFS(IKFFS,KFS)=NPFS(IKFFS,KFS)+1
0245   260   CONTINUE
0246
0247 C...Write statistics on particle/parton composition of events.
0248       ELSEIF(MTABU.EQ.22) THEN
0249         FAC=1D0/MAX(1,NEVFS)
0250         WRITE(MSTU(11),5200) NEVFS,FAC*NPRFS,FAC*NFIFS,FAC*NCHFS
0251         DO 270 I=1,NKFFS
0252           CALL PYNAME(KFFS(I),CHAU)
0253           KC=PYCOMP(KFFS(I))
0254           MDCYF=0
0255           IF(KC.NE.0) MDCYF=MDCY(KC,1)
0256           WRITE(MSTU(11),5300) KFFS(I),CHAU,MDCYF,(FAC*NPFS(I,J),J=1,4),
0257      &    FAC*(NPFS(I,1)+NPFS(I,2)+NPFS(I,3)+NPFS(I,4))
0258   270   CONTINUE
0259
0260 C...Copy particle/parton composition information into /PYJETS/.
0261       ELSEIF(MTABU.EQ.23) THEN
0262         FAC=1D0/MAX(1,NEVFS)
0263         DO 290 I=1,NKFFS
0264           K(I,1)=32
0265           K(I,2)=99
0266           K(I,3)=KFFS(I)
0267           K(I,4)=0
0268           K(I,5)=NPFS(I,1)+NPFS(I,2)+NPFS(I,3)+NPFS(I,4)
0269           DO 280 J=1,4
0270             P(I,J)=FAC*NPFS(I,J)
0271             V(I,J)=0D0
0272   280     CONTINUE
0273           P(I,5)=FAC*K(I,5)
0274           V(I,5)=0D0
0275   290   CONTINUE
0276         N=NKFFS
0277         DO 300 J=1,5
0278           K(N+1,J)=0
0279           P(N+1,J)=0D0
0280           V(N+1,J)=0D0
0281   300   CONTINUE
0282         K(N+1,1)=32
0283         K(N+1,2)=99
0284         K(N+1,5)=NEVFS
0285         P(N+1,1)=FAC*NPRFS
0286         P(N+1,2)=FAC*NFIFS
0287         P(N+1,3)=FAC*NCHFS
0288         MSTU(3)=1
0289
0290 C...Reset factorial moments statistics.
0291       ELSEIF(MTABU.EQ.30) THEN
0292         NEVFM=0
0293         NMUFM=0
0294         DO 330 IM=1,3
0295           DO 320 IB=1,10
0296             DO 310 IP=1,4
0297               FM1FM(IM,IB,IP)=0D0
0298               FM2FM(IM,IB,IP)=0D0
0299   310       CONTINUE
0300   320     CONTINUE
0301   330   CONTINUE
0302
0303 C...Find particles to include, with (pion,pseudo)rapidity and azimuth.
0304       ELSEIF(MTABU.EQ.31) THEN
0305         NEVFM=NEVFM+1
0306         NLOW=N+MSTU(3)
0307         NUPP=NLOW
0308         DO 410 I=1,N
0309           IF(K(I,1).LE.0.OR.K(I,1).GT.10) GOTO 410
0310           IF(MSTU(41).GE.2) THEN
0311             KC=PYCOMP(K(I,2))
0312             IF(KC.EQ.0.OR.KC.EQ.12.OR.KC.EQ.14.OR.KC.EQ.16.OR.
0313      &      KC.EQ.18.OR.K(I,2).EQ.KSUSY1+22.OR.K(I,2).EQ.39.OR.
0314      &      K(I,2).EQ.KSUSY1+39) GOTO 410
0315             IF(MSTU(41).GE.3.AND.KCHG(KC,2).EQ.0.AND.
0316      &      PYCHGE(K(I,2)).EQ.0) GOTO 410
0317           ENDIF
0318           PMR=0D0
0319           IF(MSTU(42).EQ.1.AND.K(I,2).NE.22) PMR=PYMASS(211)
0320           IF(MSTU(42).GE.2) PMR=P(I,5)
0321           PR=MAX(1D-20,PMR**2+P(I,1)**2+P(I,2)**2)
0322           YETA=SIGN(LOG(MIN((SQRT(PR+P(I,3)**2)+ABS(P(I,3)))/SQRT(PR),
0323      &    1D20)),P(I,3))
0324           IF(ABS(YETA).GT.PARU(57)) GOTO 410
0325           PHI=PYANGL(P(I,1),P(I,2))
0326           IYETA=512D0*(YETA+PARU(57))/(2D0*PARU(57))
0327           IYETA=MAX(0,MIN(511,IYETA))
0328           IPHI=512D0*(PHI+PARU(1))/PARU(2)
0329           IPHI=MAX(0,MIN(511,IPHI))
0330           IYEP=0
0331           DO 340 IB=0,9
0332             IYEP=IYEP+4**IB*(2*MOD(IYETA/2**IB,2)+MOD(IPHI/2**IB,2))
0333   340     CONTINUE
0334
0335 C...Order particles in (pseudo)rapidity and/or azimuth.
0336           IF(NUPP.GT.MSTU(4)-5-MSTU(32)) THEN
0337             CALL PYERRM(11,'(PYTABU:) no more memory left in PYJETS')
0338             RETURN
0339           ENDIF
0340           NUPP=NUPP+1
0341           IF(NUPP.EQ.NLOW+1) THEN
0342             K(NUPP,1)=IYETA
0343             K(NUPP,2)=IPHI
0344             K(NUPP,3)=IYEP
0345           ELSE
0346             DO 350 I1=NUPP-1,NLOW+1,-1
0347               IF(IYETA.GE.K(I1,1)) GOTO 360
0348               K(I1+1,1)=K(I1,1)
0349   350       CONTINUE
0350   360       K(I1+1,1)=IYETA
0351             DO 370 I1=NUPP-1,NLOW+1,-1
0352               IF(IPHI.GE.K(I1,2)) GOTO 380
0353               K(I1+1,2)=K(I1,2)
0354   370       CONTINUE
0355   380       K(I1+1,2)=IPHI
0356             DO 390 I1=NUPP-1,NLOW+1,-1
0357               IF(IYEP.GE.K(I1,3)) GOTO 400
0358               K(I1+1,3)=K(I1,3)
0359   390       CONTINUE
0360   400       K(I1+1,3)=IYEP
0361           ENDIF
0362   410   CONTINUE
0363         K(NUPP+1,1)=2**10
0364         K(NUPP+1,2)=2**10
0365         K(NUPP+1,3)=4**10
0366
0367 C...Calculate sum of factorial moments in event.
0368         DO 480 IM=1,3
0369           DO 430 IB=1,10
0370             DO 420 IP=1,4
0371               FEVFM(IB,IP)=0D0
0372   420       CONTINUE
0373   430     CONTINUE
0374           DO 450 IB=1,10
0375             IF(IM.LE.2) IBIN=2**(10-IB)
0376             IF(IM.EQ.3) IBIN=4**(10-IB)
0377             IAGR=K(NLOW+1,IM)/IBIN
0378             NAGR=1
0379             DO 440 I=NLOW+2,NUPP+1
0380               ICUT=K(I,IM)/IBIN
0381               IF(ICUT.EQ.IAGR) THEN
0382                 NAGR=NAGR+1
0383               ELSE
0384                 IF(NAGR.EQ.1) THEN
0385                 ELSEIF(NAGR.EQ.2) THEN
0386                   FEVFM(IB,1)=FEVFM(IB,1)+2D0
0387                 ELSEIF(NAGR.EQ.3) THEN
0388                   FEVFM(IB,1)=FEVFM(IB,1)+6D0
0389                   FEVFM(IB,2)=FEVFM(IB,2)+6D0
0390                 ELSEIF(NAGR.EQ.4) THEN
0391                   FEVFM(IB,1)=FEVFM(IB,1)+12D0
0392                   FEVFM(IB,2)=FEVFM(IB,2)+24D0
0393                   FEVFM(IB,3)=FEVFM(IB,3)+24D0
0394                 ELSE
0395                   FEVFM(IB,1)=FEVFM(IB,1)+NAGR*(NAGR-1D0)
0396                   FEVFM(IB,2)=FEVFM(IB,2)+NAGR*(NAGR-1D0)*(NAGR-2D0)
0397                   FEVFM(IB,3)=FEVFM(IB,3)+NAGR*(NAGR-1D0)*(NAGR-2D0)*
0398      &            (NAGR-3D0)
0399                   FEVFM(IB,4)=FEVFM(IB,4)+NAGR*(NAGR-1D0)*(NAGR-2D0)*
0400      &            (NAGR-3D0)*(NAGR-4D0)
0401                 ENDIF
0402                 IAGR=ICUT
0403                 NAGR=1
0404               ENDIF
0405   440       CONTINUE
0406   450     CONTINUE
0407
0408 C...Add results to total statistics.
0409           DO 470 IB=10,1,-1
0410             DO 460 IP=1,4
0411               IF(FEVFM(1,IP).LT.0.5D0) THEN
0412                 FEVFM(IB,IP)=0D0
0413               ELSEIF(IM.LE.2) THEN
0414                 FEVFM(IB,IP)=2D0**((IB-1)*IP)*FEVFM(IB,IP)/FEVFM(1,IP)
0415               ELSE
0416                 FEVFM(IB,IP)=4D0**((IB-1)*IP)*FEVFM(IB,IP)/FEVFM(1,IP)
0417               ENDIF
0418               FM1FM(IM,IB,IP)=FM1FM(IM,IB,IP)+FEVFM(IB,IP)
0419               FM2FM(IM,IB,IP)=FM2FM(IM,IB,IP)+FEVFM(IB,IP)**2
0420   460       CONTINUE
0421   470     CONTINUE
0422   480   CONTINUE
0423         NMUFM=NMUFM+(NUPP-NLOW)
0424         MSTU(62)=NUPP-NLOW
0425
0426 C...Write accumulated statistics on factorial moments.
0427       ELSEIF(MTABU.EQ.32) THEN
0428         FAC=1D0/MAX(1,NEVFM)
0429         IF(MSTU(42).LE.0) WRITE(MSTU(11),5400) NEVFM,'eta'
0430         IF(MSTU(42).EQ.1) WRITE(MSTU(11),5400) NEVFM,'ypi'
0431         IF(MSTU(42).GE.2) WRITE(MSTU(11),5400) NEVFM,'y  '
0432         DO 510 IM=1,3
0433           WRITE(MSTU(11),5500)
0434           DO 500 IB=1,10
0435             BYETA=2D0*PARU(57)
0436             IF(IM.NE.2) BYETA=BYETA/2**(IB-1)
0437             BPHI=PARU(2)
0438             IF(IM.NE.1) BPHI=BPHI/2**(IB-1)
0439             IF(IM.LE.2) BNAVE=FAC*NMUFM/DBLE(2**(IB-1))
0440             IF(IM.EQ.3) BNAVE=FAC*NMUFM/DBLE(4**(IB-1))
0441             DO 490 IP=1,4
0442               FMOMA(IP)=FAC*FM1FM(IM,IB,IP)
0443               FMOMS(IP)=SQRT(MAX(0D0,FAC*(FAC*FM2FM(IM,IB,IP)-
0444      &        FMOMA(IP)**2)))
0445   490       CONTINUE
0446             WRITE(MSTU(11),5600) BYETA,BPHI,BNAVE,(FMOMA(IP),FMOMS(IP),
0447      &      IP=1,4)
0448   500     CONTINUE
0449   510   CONTINUE
0450
0451 C...Copy statistics on factorial moments into /PYJETS/.
0452       ELSEIF(MTABU.EQ.33) THEN
0453         FAC=1D0/MAX(1,NEVFM)
0454         DO 540 IM=1,3
0455           DO 530 IB=1,10
0456             I=10*(IM-1)+IB
0457             K(I,1)=32
0458             K(I,2)=99
0459             K(I,3)=1
0460             IF(IM.NE.2) K(I,3)=2**(IB-1)
0461             K(I,4)=1
0462             IF(IM.NE.1) K(I,4)=2**(IB-1)
0463             K(I,5)=0
0464             P(I,1)=2D0*PARU(57)/K(I,3)
0465             V(I,1)=PARU(2)/K(I,4)
0466             DO 520 IP=1,4
0467               P(I,IP+1)=FAC*FM1FM(IM,IB,IP)
0468               V(I,IP+1)=SQRT(MAX(0D0,FAC*(FAC*FM2FM(IM,IB,IP)-
0469      &        P(I,IP+1)**2)))
0470   520       CONTINUE
0471   530     CONTINUE
0472   540   CONTINUE
0473         N=30
0474         DO 550 J=1,5
0475           K(N+1,J)=0
0476           P(N+1,J)=0D0
0477           V(N+1,J)=0D0
0478   550   CONTINUE
0479         K(N+1,1)=32
0480         K(N+1,2)=99
0481         K(N+1,5)=NEVFM
0482         MSTU(3)=1
0483
0484 C...Reset statistics on Energy-Energy Correlation.
0485       ELSEIF(MTABU.EQ.40) THEN
0486         NEVEE=0
0487         DO 560 J=1,25
0488           FE1EC(J)=0D0
0489           FE2EC(J)=0D0
0490           FE1EC(51-J)=0D0
0491           FE2EC(51-J)=0D0
0492           FE1EA(J)=0D0
0493           FE2EA(J)=0D0
0494   560   CONTINUE
0495
0496 C...Find particles to include, with proper assumed mass.
0497       ELSEIF(MTABU.EQ.41) THEN
0498         NEVEE=NEVEE+1
0499         NLOW=N+MSTU(3)
0500         NUPP=NLOW
0501         ECM=0D0
0502         DO 570 I=1,N
0503           IF(K(I,1).LE.0.OR.K(I,1).GT.10) GOTO 570
0504           IF(MSTU(41).GE.2) THEN
0505             KC=PYCOMP(K(I,2))
0506             IF(KC.EQ.0.OR.KC.EQ.12.OR.KC.EQ.14.OR.KC.EQ.16.OR.
0507      &      KC.EQ.18.OR.K(I,2).EQ.KSUSY1+22.OR.K(I,2).EQ.39.OR.
0508      &      K(I,2).EQ.KSUSY1+39) GOTO 570
0509             IF(MSTU(41).GE.3.AND.KCHG(KC,2).EQ.0.AND.
0510      &      PYCHGE(K(I,2)).EQ.0) GOTO 570
0511           ENDIF
0512           PMR=0D0
0513           IF(MSTU(42).EQ.1.AND.K(I,2).NE.22) PMR=PYMASS(211)
0514           IF(MSTU(42).GE.2) PMR=P(I,5)
0515           IF(NUPP.GT.MSTU(4)-5-MSTU(32)) THEN
0516             CALL PYERRM(11,'(PYTABU:) no more memory left in PYJETS')
0517             RETURN
0518           ENDIF
0519           NUPP=NUPP+1
0520           P(NUPP,1)=P(I,1)
0521           P(NUPP,2)=P(I,2)
0522           P(NUPP,3)=P(I,3)
0523           P(NUPP,4)=SQRT(PMR**2+P(I,1)**2+P(I,2)**2+P(I,3)**2)
0524           P(NUPP,5)=MAX(1D-10,SQRT(P(I,1)**2+P(I,2)**2+P(I,3)**2))
0525           ECM=ECM+P(NUPP,4)
0526   570   CONTINUE
0527         IF(NUPP.EQ.NLOW) RETURN
0528
0529 C...Analyze Energy-Energy Correlation in event.
0530         FAC=(2D0/ECM**2)*50D0/PARU(1)
0531         DO 580 J=1,50
0532           FEVEE(J)=0D0
0533   580   CONTINUE
0534         DO 600 I1=NLOW+2,NUPP
0535           DO 590 I2=NLOW+1,I1-1
0536             CTHE=(P(I1,1)*P(I2,1)+P(I1,2)*P(I2,2)+P(I1,3)*P(I2,3))/
0537      &      (P(I1,5)*P(I2,5))
0538             THE=ACOS(MAX(-1D0,MIN(1D0,CTHE)))
0539             ITHE=MAX(1,MIN(50,1+INT(50D0*THE/PARU(1))))
0540             FEVEE(ITHE)=FEVEE(ITHE)+FAC*P(I1,4)*P(I2,4)
0541   590     CONTINUE
0542   600   CONTINUE
0543         DO 610 J=1,25
0544           FE1EC(J)=FE1EC(J)+FEVEE(J)
0545           FE2EC(J)=FE2EC(J)+FEVEE(J)**2
0546           FE1EC(51-J)=FE1EC(51-J)+FEVEE(51-J)
0547           FE2EC(51-J)=FE2EC(51-J)+FEVEE(51-J)**2
0548           FE1EA(J)=FE1EA(J)+(FEVEE(51-J)-FEVEE(J))
0549           FE2EA(J)=FE2EA(J)+(FEVEE(51-J)-FEVEE(J))**2
0550   610   CONTINUE
0551         MSTU(62)=NUPP-NLOW
0552
0553 C...Write statistics on Energy-Energy Correlation.
0554       ELSEIF(MTABU.EQ.42) THEN
0555         FAC=1D0/MAX(1,NEVEE)
0556         WRITE(MSTU(11),5700) NEVEE
0557         DO 620 J=1,25
0558           FEEC1=FAC*FE1EC(J)
0559           FEES1=SQRT(MAX(0D0,FAC*(FAC*FE2EC(J)-FEEC1**2)))
0560           FEEC2=FAC*FE1EC(51-J)
0561           FEES2=SQRT(MAX(0D0,FAC*(FAC*FE2EC(51-J)-FEEC2**2)))
0562           FEECA=FAC*FE1EA(J)
0563           FEESA=SQRT(MAX(0D0,FAC*(FAC*FE2EA(J)-FEECA**2)))
0564           WRITE(MSTU(11),5800) 3.6D0*(J-1),3.6D0*J,FEEC1,FEES1,
0565      &    FEEC2,FEES2,FEECA,FEESA
0566   620   CONTINUE
0567
0568 C...Copy statistics on Energy-Energy Correlation into /PYJETS/.
0569       ELSEIF(MTABU.EQ.43) THEN
0570         FAC=1D0/MAX(1,NEVEE)
0571         DO 630 I=1,25
0572           K(I,1)=32
0573           K(I,2)=99
0574           K(I,3)=0
0575           K(I,4)=0
0576           K(I,5)=0
0577           P(I,1)=FAC*FE1EC(I)
0578           V(I,1)=SQRT(MAX(0D0,FAC*(FAC*FE2EC(I)-P(I,1)**2)))
0579           P(I,2)=FAC*FE1EC(51-I)
0580           V(I,2)=SQRT(MAX(0D0,FAC*(FAC*FE2EC(51-I)-P(I,2)**2)))
0581           P(I,3)=FAC*FE1EA(I)
0582           V(I,3)=SQRT(MAX(0D0,FAC*(FAC*FE2EA(I)-P(I,3)**2)))
0583           P(I,4)=PARU(1)*(I-1)/50D0
0584           P(I,5)=PARU(1)*I/50D0
0585           V(I,4)=3.6D0*(I-1)
0586           V(I,5)=3.6D0*I
0587   630   CONTINUE
0588         N=25
0589         DO 640 J=1,5
0590           K(N+1,J)=0
0591           P(N+1,J)=0D0
0592           V(N+1,J)=0D0
0593   640   CONTINUE
0594         K(N+1,1)=32
0595         K(N+1,2)=99
0596         K(N+1,5)=NEVEE
0597         MSTU(3)=1
0598
0599 C...Reset statistics on decay channels.
0600       ELSEIF(MTABU.EQ.50) THEN
0601         NEVDC=0
0602         NKFDC=0
0603         NREDC=0
0604
0605 C...Identify and order flavour content of final state.
0606       ELSEIF(MTABU.EQ.51) THEN
0607         NEVDC=NEVDC+1
0608         NDS=0
0609         DO 670 I=1,N
0610           IF(K(I,1).LE.0.OR.K(I,1).GE.6) GOTO 670
0611           NDS=NDS+1
0612           IF(NDS.GT.8) THEN
0613             NREDC=NREDC+1
0614             RETURN
0615           ENDIF
0616           KFM=2*IABS(K(I,2))
0617           IF(K(I,2).LT.0) KFM=KFM-1
0618           DO 650 IDS=NDS-1,1,-1
0619             IIN=IDS+1
0620             IF(KFM.LT.KFDM(IDS)) GOTO 660
0621             KFDM(IDS+1)=KFDM(IDS)
0622   650     CONTINUE
0623           IIN=1
0624   660     KFDM(IIN)=KFM
0625   670   CONTINUE
0626
0627 C...Find whether old or new final state.
0628         DO 690 IDC=1,NKFDC
0629           IF(NDS.LT.KFDC(IDC,0)) THEN
0630             IKFDC=IDC
0631             GOTO 700
0632           ELSEIF(NDS.EQ.KFDC(IDC,0)) THEN
0633             DO 680 I=1,NDS
0634               IF(KFDM(I).LT.KFDC(IDC,I)) THEN
0635                 IKFDC=IDC
0636                 GOTO 700
0637               ELSEIF(KFDM(I).GT.KFDC(IDC,I)) THEN
0638                 GOTO 690
0639               ENDIF
0640   680       CONTINUE
0641             IKFDC=-IDC
0642             GOTO 700
0643           ENDIF
0644   690   CONTINUE
0645         IKFDC=NKFDC+1
0646   700   IF(IKFDC.LT.0) THEN
0647           IKFDC=-IKFDC
0648         ELSEIF(NKFDC.GE.200) THEN
0649           NREDC=NREDC+1
0650           RETURN
0651         ELSE
0652           DO 720 IDC=NKFDC,IKFDC,-1
0653             NPDC(IDC+1)=NPDC(IDC)
0654             DO 710 I=0,8
0655               KFDC(IDC+1,I)=KFDC(IDC,I)
0656   710       CONTINUE
0657   720     CONTINUE
0658           NKFDC=NKFDC+1
0659           KFDC(IKFDC,0)=NDS
0660           DO 730 I=1,NDS
0661             KFDC(IKFDC,I)=KFDM(I)
0662   730     CONTINUE
0663           NPDC(IKFDC)=0
0664         ENDIF
0665         NPDC(IKFDC)=NPDC(IKFDC)+1
0666
0667 C...Write statistics on decay channels.
0668       ELSEIF(MTABU.EQ.52) THEN
0669         FAC=1D0/MAX(1,NEVDC)
0670         WRITE(MSTU(11),5900) NEVDC
0671         DO 750 IDC=1,NKFDC
0672           DO 740 I=1,KFDC(IDC,0)
0673             KFM=KFDC(IDC,I)
0674             KF=(KFM+1)/2
0675             IF(2*KF.NE.KFM) KF=-KF
0676             CALL PYNAME(KF,CHAU)
0677             CHDC(I)=CHAU(1:12)
0678             IF(CHAU(13:13).NE.' ') CHDC(I)(12:12)='?'
0679   740     CONTINUE
0680           WRITE(MSTU(11),6000) FAC*NPDC(IDC),(CHDC(I),I=1,KFDC(IDC,0))
0681   750   CONTINUE
0682         IF(NREDC.NE.0) WRITE(MSTU(11),6100) FAC*NREDC
0683
0684 C...Copy statistics on decay channels into /PYJETS/.
0685       ELSEIF(MTABU.EQ.53) THEN
0686         FAC=1D0/MAX(1,NEVDC)
0687         DO 780 IDC=1,NKFDC
0688           K(IDC,1)=32
0689           K(IDC,2)=99
0690           K(IDC,3)=0
0691           K(IDC,4)=0
0692           K(IDC,5)=KFDC(IDC,0)
0693           DO 760 J=1,5
0694             P(IDC,J)=0D0
0695             V(IDC,J)=0D0
0696   760     CONTINUE
0697           DO 770 I=1,KFDC(IDC,0)
0698             KFM=KFDC(IDC,I)
0699             KF=(KFM+1)/2
0700             IF(2*KF.NE.KFM) KF=-KF
0701             IF(I.LE.5) P(IDC,I)=KF
0702             IF(I.GE.6) V(IDC,I-5)=KF
0703   770     CONTINUE
0704           V(IDC,5)=FAC*NPDC(IDC)
0705   780   CONTINUE
0706         N=NKFDC
0707         DO 790 J=1,5
0708           K(N+1,J)=0
0709           P(N+1,J)=0D0
0710           V(N+1,J)=0D0
0711   790   CONTINUE
0712         K(N+1,1)=32
0713         K(N+1,2)=99
0714         K(N+1,5)=NEVDC
0715         V(N+1,5)=FAC*NREDC
0716         MSTU(3)=1
0717       ENDIF
0718
0719 C...Format statements for output on unit MSTU(11) (default 6).
0720  5000 FORMAT(///20X,'Event statistics - initial state'/
0721      &20X,'based on an analysis of ',I6,' events'//
0722      &3X,'Main flavours after',8X,'Fraction',4X,'Subfractions ',
0723      &'according to fragmenting system multiplicity'/
0724      &4X,'hard interaction',24X,'1',7X,'2',7X,'3',7X,'4',7X,'5',
0725      &6X,'6-7',5X,'8-10',3X,'11-15',3X,'16-25',4X,'>25'/)
0726  5100 FORMAT(3X,A12,1X,A12,F10.5,1X,10F8.4)
0727  5200 FORMAT(///20X,'Event statistics - final state'/
0728      &20X,'based on an analysis of ',I7,' events'//
0729      &5X,'Mean primary multiplicity =',F10.4/
0730      &5X,'Mean final   multiplicity =',F10.4/
0731      &5X,'Mean charged multiplicity =',F10.4//
0732      &5X,'Number of particles produced per event (directly and via ',
0733      &'decays/branchings)'/
0734      &8X,'KF    Particle/jet  MDCY',10X,'Particles',13X,'Antiparticles',
0735      &8X,'Total'/35X,'prim        seco        prim        seco'/)
0736  5300 FORMAT(1X,I9,4X,A16,I2,5(1X,F11.6))
0737  5400 FORMAT(///20X,'Factorial moments analysis of multiplicity'/
0738      &20X,'based on an analysis of ',I6,' events'//
0739      &3X,'delta-',A3,' delta-phi     <n>/bin',10X,'<F2>',18X,'<F3>',
0740      &18X,'<F4>',18X,'<F5>'/35X,4('     value     error  '))
0741  5500 FORMAT(10X)
0742  5600 FORMAT(2X,2F10.4,F12.4,4(F12.4,F10.4))
0743  5700 FORMAT(///20X,'Energy-Energy Correlation and Asymmetry'/
0744      &20X,'based on an analysis of ',I6,' events'//
0745      &2X,'theta range',8X,'EEC(theta)',8X,'EEC(180-theta)',7X,
0746      &'EECA(theta)'/2X,'in degrees ',3('      value    error')/)
0747  5800 FORMAT(2X,F4.1,' - ',F4.1,3(F11.4,F9.4))
0748  5900 FORMAT(///20X,'Decay channel analysis - final state'/
0749      &20X,'based on an analysis of ',I6,' events'//
0750      &2X,'Probability',10X,'Complete final state'/)
0751  6000 FORMAT(2X,F9.5,5X,8(A12,1X))
0752  6100 FORMAT(2X,F9.5,5X,'into other channels (more than 8 particles ',
0753      &'or table overflow)')
0754
0755       RETURN
0756       END