C...Program to compare event properties of total cross section
C...(including elastic and diffractive, unless changed) of
C...gamma-gamma, gamma-p and p-p at some desired energy.
C...Requires a second separate run to produce HPLOT output.

      COMMON/LUJETS/N,K(4000,5),P(4000,5),V(4000,5)
      COMMON/LUDAT1/MSTU(200),PARU(200),MSTJ(200),PARJ(200)
      COMMON/LUDAT2/KCHG(500,3),PMAS(500,4),PARF(2000),VCKM(4,4)
      COMMON/LUDAT3/MDCY(500,3),MDME(2000,2),BRAT(2000),KFDP(2000,5)
      COMMON/PYSUBS/MSEL,MSUB(200),KFIN(2,-40:40),CKIN(200)
      COMMON/PYPARS/MSTP(200),PARP(200),MSTI(200),PARI(200)
      COMMON/PYINT1/MINT(400),VINT(400)
      DIMENSION NTR(10)
      DATA PI/3.1415926/, NTR/10*0/
 
C...Number of events, CM energy, all processes, photon treatment.
      NEV=50000
      ECM=20.
      MSEL=2
      MSTP(14)=10

C***Special run without QCD radiation or multiple interactions.
C      MSTP(61)=0
C      MSTP(71)=0
C      MSTP(81)=0

C...Set grid for LUCELL. PARU(53) = min jet energy. 
      MSTU(51)=100
      MSTU(52)=64
      PARU(51)=5.
      PARU(52)=1.
      PARU(53)=5.
      PARU(54)=1.      
 
C...Book histograms.
      DO 100 ICLA=1,10
      CALL GBOOK1(ICLA,'dN(jet)/d(pT) hardest',40,0.,20.)
      CALL GBOOK1(ICLA+10,'sum ET',50,0.,100.)
      CALL GBOOK1(ICLA+20,'dET/dy',100,-10.,10.)
      CALL GBOOK1(ICLA+30,'N(charged)',50,0.,100.)
      CALL GBOOK1(ICLA+40,'dN(charged)/dy',100,-10.,10.)
      CALL GBOOK1(ICLA+50,'dN(charged)/d(pT)',100,0.,5.)
      CALL GBOOK1(ICLA+60,'maximum rapidity gap',80,0.,20.)
      CALL GBOOK1(ICLA+70,'ET spectrum of jets',50,0.,25.)
      CALL GBOOK1(ICLA+80,'eta spectrum of jets',80,-8.,8.)
      CALL GBOOK1(ICLA+90,'pT beam remnant',80,0.,4.)
  100 CONTINUE

C...Loop over gamma-gamma, gamma-p and p-p.
      DO 250 IBM=1,3
      ISM=6+IBM
      
C...Initialize.
      MSTP(122)=0
      IF(IBM.EQ.1) THEN
        CALL PYINIT('CMS','GAMMA','GAMMA',ECM)
      ELSEIF(IBM.EQ.2) THEN
        CALL PYINIT('CMS','GAMMA','P',ECM)
      ELSEIF(IBM.EQ.3) THEN
C...Here we have to set pTmin as in the other cases.
        PARP(81)=1.30+0.15*LOG(ECM/200.)/LOG(900./200.)
        CKIN(1)=0. 
        CKIN(3)=0.    
        CALL PYINIT('CMS','P','P',ECM)
      ENDIF
 
C...Generate events; look at first few.
      DO 200 IEV=1,NEV
      CALL PYEVNT
c      IF(IEV.LE.2) CALL LULIST(1)
      ITY=MINT(122)
      IF(IBM.GE.2) ITY=10
      NTR(ISM)=NTR(ISM)+1
      NTR(ITY)=NTR(ITY)+1

C...Hard process pT spectrum.
      PTHARD=PARI(17)
      CALL GFILL1(ISM,PTHARD,1.)
      CALL GFILL1(ITY,PTHARD,1.)

C...pT of photon beam remnant.
      IF(MINT(11).EQ.22.AND.MINT(107).EQ.0) THEN
        PTREM=PARI(17)
      ELSEIF(MINT(107).EQ.2.OR.MINT(11).NE.22) THEN
        PTREM=VINT(157)
        IF(MSTI(1).EQ.95) THEN
          DO 110 I=1,N
  110     IF(K(I,3).EQ.1) PTREM=SQRT(P(I,1)**2+P(I,2)**2)
        ENDIF
      ELSE
        PTREM=SQRT(MAX(0.,VINT(283)))
      ENDIF
      CALL GFILL1(ISM+90,PTREM,1.)
      CALL GFILL1(ITY+90,PTREM,1.)
      IF(MINT(12).EQ.22.AND.MINT(108).EQ.0) THEN
        PTREM=PARI(17)
      ELSEIF(MINT(108).EQ.2.OR.MINT(12).NE.22) THEN
        PTREM=VINT(158)
        IF(MSTI(1).EQ.95) THEN
          DO 115 I=1,N
  115     IF(K(I,3).EQ.2) PTREM=SQRT(P(I,1)**2+P(I,2)**2)
        ENDIF
      ELSE
        PTREM=SQRT(MAX(0.,VINT(284)))
      ENDIF
      CALL GFILL1(ISM+90,PTREM,1.)
      CALL GFILL1(ITY+90,PTREM,1.)

C...ET flow as function of rapidity and total ET.
      CALL LUEDIT(2)
      ETSUM=0.
      DO 120 I=1,N
      Y=PLU(I,17)
      ET=P(I,4)*SQRT((P(I,1)**2+P(I,2)**2)/
     &(P(I,1)**2+P(I,2)**2+P(I,3)**2))
      CALL GFILL1(ISM+20,Y,ET)
      CALL GFILL1(ITY+20,Y,ET)
      ETSUM=ETSUM+ET
  120 CONTINUE
      CALL GFILL1(ISM+10,ETSUM,1.)
      CALL GFILL1(ITY+10,ETSUM,1.)

C...Jet properties.   
      CALL LUCELL(NJET)
      DO 140 IJET=1,NJET
      CALL GFILL1(ISM+70,P(N+IJET,5),1.)
      CALL GFILL1(ITY+70,P(N+IJET,5),1.)
      CALL GFILL1(ISM+80,P(N+IJET,1),1.)
      CALL GFILL1(ITY+80,P(N+IJET,1),1.)
  140 CONTINUE

C...Charged multiplicity.
      CALL LUEDIT(3)
      XNCH=N
      CALL GFILL1(ISM+30,XNCH,1.)
      CALL GFILL1(ITY+30,XNCH,1.)

C...Charged rapidity and pT spectra.
      DO 150 I=1,N
      Y=PLU(I,17)
      CALL GFILL1(ISM+40,Y,1.)
      CALL GFILL1(ITY+40,Y,1.)
      PT=SQRT(P(I,1)**2+P(I,2)**2)
      CALL GFILL1(ISM+50,PT,1.)
      CALL GFILL1(ITY+50,PT,1.)
  150 CONTINUE

C...Find size of maximum gap.
      CALL MAXGAP(YMXGAP)
      CALL GFILL1(ISM+60,YMXGAP,1.)
      CALL GFILL1(ITY+60,YMXGAP,1.)

C...End of event loop.
  200 CONTINUE
      CALL PYSTAT(1)
  250 CONTINUE

C...Normalize and print histograms.
      DO 300 ICLA=1,10
      FREV=1./MAX(1,NTR(ICLA))
      CALL GSCALE(ICLA,2.*FREV)
      CALL GSCALE(ICLA+10,2.*FREV)
      CALL GSCALE(ICLA+20,5.*FREV)
      CALL GSCALE(ICLA+30,FREV)
      CALL GSCALE(ICLA+40,5.*FREV)
      CALL GSCALE(ICLA+50,20.*FREV)
      CALL GSCALE(ICLA+60,4.*FREV)
      CALL GSCALE(ICLA+70,2.*FREV)
      CALL GSCALE(ICLA+80,5.*FREV)
      CALL GSCALE(ICLA+90,20.*FREV)
  300 CONTINUE
      CALL GDUMP('gout.33') 
c      CALL GCLEAR
      
      END

C*******************************************************************

      SUBROUTINE MAXGAP(YMX)
      COMMON/LUJETS/N,K(4000,5),P(4000,5),V(4000,5)
      DIMENSION YY(100)

      NY=MIN(100,N)
      DO 110 I=1,NY
      Y=PLU(I,17)
      II=I
  100 II=II-1
      IF(II.EQ.0) THEN
      ELSEIF(Y.GT.YY(II)) THEN
        YY(II+1)=YY(II)
        GOTO 100
      ENDIF
      YY(II+1)=Y
  110 CONTINUE

      YMX=0.
      DO 120 I=1,NY-1
      DY=YY(I)-YY(I+1)
      IF(DY.GT.YMX) YMX=DY
  120 CONTINUE

      RETURN
      END