C***********************************************************************
C
C This is a set of routines to perform a standard QCD analysis of 
C event shapes and inclusive charged particle distributions
C with KINGAL (HVFL05, based on JETSET7.4) (used for full LEP1 MC)
C Histogram-based, not Ntuple-based
C
C Author : Guenther Dissertori / 23-MAY-01
C
C Version: 1.0 First setup
C***********************************************************************

C-----------------------------------------------------------------------
C-----------------------------------------------------------------------


      SUBROUTINE USER_INIT
C-----------------------------------------------------------------------
C user initialization routine
C =================================
C
C Author: Guenther Dissertori  
C-----------------------------------------------------------------------

      Implicit NONE

C -- LUND COMMON BLOCKS :   K has an extra P here
C -- JETSET part
      Integer    N,KP(4000,5),MSTU(200),MDCY(500,3),MDME(2000,2)
      Integer    KFDP(2000,5),MRLU(6),KCHG(500,3),MSTJ(200)
      Real*4     P(4000,5),V(4000,5),PARU(200),PARJ(200)
      Real*4     PMAS(500,4),PARF(2000),VCKM(4,4)
      Real*4     BRAT(2000),CHAF(500),RRLU(100)

      Common /LUJETS/N,KP,P,V
      Common /LUDAT1/MSTU,PARU,MSTJ,PARJ


      Integer   i,j,LU
C -- output unit
      Parameter (LU=6)

C -- event counter
      Integer   Nhad
      Common /Norm/ Nhad

C -- the histo definitions
      Include 'histonaming.h'

C -- BOS stuff
      INTEGER LMHLEN,LMHROW,LMHCOL,LBCS,LCHAR,IW
      PARAMETER (LMHLEN=2, LMHROW=2, LMHCOL=1)
      PARAMETER (LBCS=100000,LCHAR=4)
      COMMON/BCS/ IW(LBCS)
      REAL RW(LBCS)
      EQUIVALENCE (RW(1),IW(1))
      INTEGER  NLINK,JGLUN,JDYM
      External NLINK

C -- book keeping
      Integer    NEV,NEVNT,LEVEL
      Real*4     PI,RADEG,ECM
      Common /MYBOOKING/ PI,RADEG,NEV,NEVNT,ECM,LEVEL

      SAVE
C-----------------------------------------------------------------------


C -- specific initalization
C -- number of events completed at hadron level
      Nhad = 0
      NEV  = 0

C -- to convert radians to degrees :
C -- ALPHA knows about PI and RADEG
      PI = ATAN2 (0., -1.) 
      RADEG = 180.0/PI

C -- HBOOK-output to unit 65
      call HOUTPU(65)
      call HERMES(65)
      
C -- get some settings from cards
      JGLUN = NLINK('GLUN',0)
      JDYM  = NLINK('GDYM',0)
      IF     (JGLUN.NE.0) THEN
        ECM = RW(JGLUN+2)
      ELSEIF (JDYM.NE.0) THEN
        ECM = RW(JDYM+2)*2.
      ELSE
        write(6,*)' ******> USER_INIT : NO ECM setting !! STOP'
        STOP
      ENDIF

      LEVEL = MSTU(200)

C -- reset particle statistics
      CALL LULIST (13)                                                  !
      CALL LUTABU (10)                                                  !
      CALL LUTABU (20)

      WRITE (LU,*) ' '
      WRITE (LU,*) '          -------------    BEGIN   --------------'  !
      WRITE (LU,*) '          ------------ QCD ANALYSIS -------------'
      WRITE (LU,*) '          ------------ with  HVFL05 -------------'
      WRITE (LU,*) ' '
      WRITE (LU,*) ' '


C -- book histos
C -- NOTE : BINNING MIGHT BE OPTIMIZED FOR SPECIFIC NEEDS

C -- event shapes
      call HBOOK1(idth,'Thrust       ',50,0.5,1.0,0.)
      call HBOOK1(idta,'Thrust Major ',50,0.0,1.0,0.)
      call HBOOK1(idtm,'Thrust Minor ',50,0.0,0.5,0.)
      call HBOOK1(idmh,'Mh2/s        ',50,0.0,0.5,0.) 
      call HBOOK1(idcp,'Cparameter   ',50,0.0,1.0,0.)
      call HBOOK1(idbt,'Btot         ',50,0.0,0.4,0.) 
      call HBOOK1(idbw,'Bwide        ',50,0.0,0.4,0.) 
      call HBOOK1(idly,'-lny3        ',50,0.0,10.,0.) 

C -- D parameter for different topologies
      call HBOOK1(iddi,'Dpar, incl   ',50,0.,1.,0.)   
      call HBOOK1(idd1,'Dpar, y=0.1  ',50,0.,1.,0.)   
      call HBOOK1(idd2,'Dpar, y=0.05 ',50,0.,1.,0.)   
      call HBOOK1(idd3,'Dpar, y=0.025',50,0.,1.,0.)   

C -- moments
      call HBOOK1(idmo,'Moments of Thrust',8,0.5,8.5,0.)

C -- three-jet rate for three ycut values
      call HBOOK1(idjr,'3 jets',3,0.5,3.5,0.)

C -- angular correlations in 4-jet events
      call HBOOK1(idbz,'|cos(BZ)| ',20,0.,1.,0.)   
      call HBOOK1(idnr,'|cos(NR)| ',20,0.,1.,0.)   
      call HBOOK1(idks,'cos(KSW)  ',20,-1.,1.,0.)   
      call HBOOK1(id34,'cos(a34)  ',20,-1.,1.,0.)   

C -- inclusive charged particle distributions
      call HBOOK1(idxp,'xp   ',50,0.,1.,0.)  ! = pch/pbeam
      call HBOOK1(idxi,'Xsi  ',50,0.,6.,0.)  ! = -ln(xp)
      call HBOOK1(idra,'y    ',50,0.,6.,0.)  ! rapidity w.r.t. thrust axis 
      call HBOOK1(idpi,'PTin ',50,0.,15.,0.) ! momentum in event plane 
      call HBOOK1(idpo,'PTout',50,0.,5.,0.)  ! momentum out of event plane

C -- angle of thrust axis
      call HBOOK1(idt1,'|cos(theta_thrust)| bef',100,0.,1.,0.) ! before selection
      call HBOOK1(idt2,'|cos(theta_thrust)| aft',100,0.,1.,0.) ! after selection
      call HBOOK1(idt3,'|cos(theta_thrust)| chg',100,0.,1.,0.) ! after selection,
C                                                                chg tracks only

C -- some info and couting
      call HBOOK1(idc1,'Iret',21,-20.5,0.5,0.)
      call HBOOK1(idc2,'counting',6,0.5,6.5,0.)

C -- some monitoring
      call HBOOK1(idet,'Etot sel.obj',150,0.0,150.0,0.)
      call HBOOK1(idnt,'Numb sel obj',100,0.0,100.0,0.)
      call HBOOK1(idfl,'flavour',4,-0.5,3.5,0.)

      CALL HMINIM (0,0.0)
      CALL HIDOPT (0,'STAT')
      CALL HIDOPT (0,'BLAC')
      CALL HIDOPT (0,'NPHI')
      CALL HSQUEZ ('YES')

      RETURN
      END


C-----------------------------------------------------------------------
C-----------------------------------------------------------------------


      SUBROUTINE USER_EVENT
C-----------------------------------------------------------------------
C Event Analysis Routine
C ============================
C
C Author: Guenther Dissertori 
C-----------------------------------------------------------------------

      Implicit NONE

C -- LUND COMMON BLOCKS :   K has an extra P here
C -- JETSET part
      Integer    N,KP(4000,5),MSTU(200),MDCY(500,3),MDME(2000,2)
      Integer    KFDP(2000,5),MRLU(6),KCHG(500,3),MSTJ(200)
      Real*4     P(4000,5),V(4000,5),PARU(200),PARJ(200)
      Real*4     PMAS(500,4),PARF(2000),VCKM(4,4)
      Real*4     BRAT(2000),CHAF(500),RRLU(100)

      Common /LUJETS/N,KP,P,V
      Common /LUDAT1/MSTU,PARU,MSTJ,PARJ


C -- Internal Workspace
      Real*4     Hmemo(200000)
      Common /PAWC/Hmemo

C -- book keeping
      Integer    NEV,NEVNT,LEVEL
      Real*4     PI,RADEG,ECM
      Common /MYBOOKING/ PI,RADEG,NEV,NEVNT,ECM,LEVEL

      Integer   NefoMax
      Parameter (NefoMax = 300)
      Real      Pefo(4,NefoMax)
      Real      etot
      Integer   i,itkl,ntrk,Nstring,flavor,mylevel,IretGlobal,iret

C -- event counting
      Integer   Nhad
      Common /Norm/ Nhad

C -- the histo definitions
      Include 'histonaming.h'

C-----------------------------------------------------------------------

C -- initialization
      NEV = NEV + 1
      IretGlobal = 0

C -- reset internal counters
      CALL LUTABU (11)
      CALL LUTABU (21)

C -- do counting of all events which are looked at
      call HFILL(idc2,1.,0.,1.)

C -- find flavor (do be done before any PYEDIT!)
C this works for jetset ps and herwig 
      flavor = 0
      do i=1,N
        if (KP(KP(i,3),2).eq.23) then ! origin is Z0
          if (abs(KP(i,2)).le.3) flavor = 1 
          if (abs(KP(i,2)).eq.4) flavor = 2 
          if (abs(KP(i,2)).eq.5) flavor = 3 
        endif
      enddo
C -- check for flavor
      if (flavor.eq.0) then
        call ErrMess('USER_EVENT, not able to determine flavor',
     >                float(flavor))
        call LULIST(1)
        IretGlobal = -1
        call HF1(idc1,float(IretGlobal),1.)        
        RETURN                  ! ----------->
      endif

C -- dump some event
      if (NEV.le.3) then
        write(MSTU(11),*)' '
        write(MSTU(11),*)'=====> NEV = ',NEV
        call LULIST(1)
      endif

C======================================================
C -- NOTE : by setting the variable "mylevel" you can choose
C -- which level should be analyzed : 
C    1=parton, 2=hadron direct from string, 3=hadron,final
C======================================================
      mylevel = LEVEL

C -- select stable final state particles
      if (mylevel.eq.3) CALL LUEDIT(1)

C -- for remembering
      if (NEV.le.2) then
        write(MSTU(11),*)' '
        write(MSTU(11),*)'=====> Level = ',mylevel
        call LULIST(1)
      endif

C -- strange event, should never happen
      if (N.le.1)  then
        write(MSTU(11),*)' '
        write(MSTU(11),*)'!!! =====> NEV = ',NEV
        write(MSTU(11),*)'!!! =====> N = 1 !!! '
        call LULIST(1)
        IretGlobal = -2
        call HF1(idc1,float(IretGlobal),1.)        
        RETURN                  ! ----------->
      endif

C -- protection against very large events
      if (N.gt.NefoMax)  then
        write(MSTU(11),*)' '
        write(MSTU(11),*)'!!! =====> NEV = ',NEV
        write(MSTU(11),*)'!!! =====> N > NefoMax !!! '
        call LULIST(1)
        IretGlobal = -3
        call HF1(idc1,float(IretGlobal),1.)        
        RETURN                  ! ----------->
      endif

C -- CoLU momenta to internal array
      itkl = 0
      etot = 0.

C -- search string
      Nstring=0
      if (mylevel.lt.3) then
        do i=1,N
          if (KP(i,2).eq.92.or.KP(i,2).eq.91) then
            Nstring=i
            goto 100
          endif
        enddo
      endif
 100  continue

      do i=1,N
C -- find partonlevel
        if (mylevel.eq.1) then
C -- check for quarks and gluons at the end of the shower
          if ((((KP(i,1).eq.11).or.(KP(i,1).eq.12))
     >        .and.((abs(KP(i,2)).le.5).or.(abs(KP(i,2)).eq.21))) .or.
C -- check for ISR and/or FSR photons
     >        ((KP(i,1).eq.1).and.(KP(i,2).eq.22).and.(i.lt.Nstring)))
     >    then 
            itkl = itkl+1
            Pefo(1,itkl) = P(i,1)
            Pefo(2,itkl) = P(i,2)
            Pefo(3,itkl) = P(i,3)
            Pefo(4,itkl) = P(i,4)
            etot         = etot + P(i,4)
            if (NEV.lt.5) then
              write(MSTU(11),*)'got a quark or gluon ... : ',
     >             i,KP(i,1),KP(i,2),etot
            endif
          endif

C -- else : you want the hadron level, search for direct hadrons from string 
        elseif (mylevel.eq.2) then
C -- check for particles for which origin is string or cluster
          if (((KP(i,3).lt.Nstring).and.(i.gt.Nstring).and.
     >         (KP(i,2).ne.91).and.(KP(i,2).ne.92).and.
     >        ((KP(KP(i,3),1).eq.11).or.(KP(KP(i,3),1).eq.12))) .or. 
C -- check for ISR and/or FSR photons
     >        ((KP(i,1).eq.1).and.(KP(i,2).eq.22).and.(i.lt.Nstring))) 
     >    then
            itkl = itkl+1
            Pefo(1,itkl) = P(i,1)
            Pefo(2,itkl) = P(i,2)
            Pefo(3,itkl) = P(i,3)
            Pefo(4,itkl) = P(i,4)
            etot         = etot + P(i,4)
          endif


C -- else : you want the hadron level, after all decays
        else
          itkl = itkl+1
          Pefo(1,itkl) = P(i,1)
          Pefo(2,itkl) = P(i,2)
          Pefo(3,itkl) = P(i,3)
          Pefo(4,itkl) = P(i,4)
          etot         = etot + P(i,4)
        endif
      enddo


      if (mylevel.eq.2.and.etot.lt.91.1) then
        write(MSTU(11),*)' ETOT problem ',itkl,etot
        call LULIST(1)
      endif

      ntrk = itkl

C -- some bookkeeping
      call HF1(idet,etot,1.)
      call HF1(idnt,float(ntrk),1.)
      call HF1(idfl,float(flavor),1.)

C -- count the events now
      Nhad=Nhad+1
      call HFILL(idc2,2.,0.,1.)

C========================================================================
C -- calculations of event shapes  + jet rates
C========================================================================
      call Calc_Evshapes(ntrk,Pefo,iret)
      
C -- bookkeeping of events for which all calculations ok
      if (iret.eq.0) call HFILL(idc2,5.,0.,1.)

C========================================================================
C -- calculations of inclusive charged particle distributions
C -- makes sense only for LEVEL=3
C========================================================================
      iret = 1
      if (mylevel.eq.3) then
        call Calc_Inclusives(iret)
      endif
C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
C -- NOTE: every analysis after this routine will have 
C          charged particles only!!
C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

C -- bookkeeping of events for which all incl. calculations ok
      if (iret.eq.0) call HFILL(idc2,6.,0.,1.)

C -- fill status code
      IretGlobal = 0
      call HFILL(idc1,float(IretGlobal),0.,1.)

C -- that's it....


 9999 RETURN 
      END


C----------------------------------------------------------------------
C----------------------------------------------------------------------


      SUBROUTINE USER_TERM
C-----------------------------------------------------------------------
C user termination routine
C =================================
C
C Author: Guenther Dissertori 
C-----------------------------------------------------------------------

      Integer      i,j

      Real         Ngood,N3j,N4j,N3e,N4e
      Real         sum1,sum2,mean,erro

      Integer   Nhad
      Common /Norm/ Nhad

      Real         Hi,Hie,Hstati,Hij
      Logical      Hexist
      External     Hi,Hie,Hexist,Hstati,Hij

C -- histo definition
      Include 'histonaming.h'

C-----------------------------------------------------------------------

      write(6,*)
      write(6,*)'*****************************************************'
      write(6,*)' ====================> final calculations ....'
      write(6,*)

C -- do statistics
      call LUTABU(22)

      if (Nhad.eq.0) then
        write(6,*)'No hadronic event found !!!!!'
        RETURN
      endif

C -- number of accepted hadronic events
      write(6,*)
      write(6,'(1X,A,I8)')'Nhad used in normalization = ',Nhad 
      write(6,*)


      write(6,'(2X,A,F10.0,A,F10.0,A)')'Out of ',HI(idc2,1),' events ',
     >      HI(idc2,2),' passed the basic criteria.'
      write(6,*)

      write(6,'(2X,A,F10.0,A,F10.0,A)')'Out of ',HI(idc2,2),' events ',
     >      HI(idc2,3),' had good event shape calculations.'
      write(6,*)

      write(6,'(2X,A,F10.0,A,F10.0,A)')'Out of ',HI(idc2,2),' events ',
     >      HI(idc2,4),' had good 4-jet calculations.'
      write(6,*)

      write(6,'(2X,A,F10.0,A,F10.0,A)')'Out of ',HI(idc2,2),' events ',
     >      HI(idc2,5),' had good event shape + jet calculations.'
      write(6,*)

      write(6,'(2X,A,F10.0,A,F10.0,A)')'Out of ',HI(idc2,2),' events ',
     >      HI(idc2,6),' had good inclusive chg tracks calculations.'
      write(6,*)
      
C -- compute the 3-jet rates
      Ngood = HI(idc2,2)
      if (Ngood.gt.0) then 
        do i=1,3
C -- get the numbers
          N3j = HI(idjr,i) 
          N3j = N3j/Ngood
          N3e = sqrt(N3j*(1-N3j)/Ngood)

          write(6,990)' 3-jet rate for ycut ',Ycut3j(i),' = ',N3j,N3e
        enddo
      endif


C -- compute the 4-jet rate
      Ngood = HI(idc2,2)
      if (Ngood.gt.0) then 
C -- get the numbers
        N4j = HI(idc2,4) 
        N4j = N4j/Ngood
        N4e = sqrt(N4j*(1-N4j)/Ngood)
        write(6,*)
        write(6,990)' 4-jet rate for ycut ',Ycut4j,' = ',N4j,N4e
      endif


C -- compute the first moment of thrust
      Ngood = HI(idc2,3)
      if (Ngood.gt.1) then 
C -- get the numbers
        sum1 = HI(idmo,1)
        sum2 = HI(idmo,5)
        mean = sum1/Ngood
        erro = Ngood/(Ngood-1.) * (sum2/Ngood - mean**2)
        erro = sqrt(max(1.e-7,erro/Ngood))
        write(6,*)
        write(6,991)' < 1 - Thrust > = ',mean,erro
      endif

      write(6,*)
      write(6,*)' ====================> end of final calculations ....'
      write(6,*)'*****************************************************'
      write(6,*)

 990  Format(2X,A,F6.3,A,F7.5," +- ",F7.5)
 991  Format(2X,A,F7.5," +- ",F7.5)

      RETURN
      END
      

C-----------------------------------------------------------------------
C-----------------------------------------------------------------------


      SUBROUTINE ErrMess(Message,Var)
C-----------------------------------------------------------------------
C Routine which prints out a message and correlated variable
C ==========================================================
C
C Author: Guenther Dissertori  / 5.1.95
C-----------------------------------------------------------------------


      Real*4         Var
      Character*(*) Message

      Integer    NEV,NEVNT,LEVEL
      Real*4     PI,RADEG,ECM
      Common /MYBOOKING/ PI,RADEG,NEV,NEVNT,ECM,LEVEL

      Integer   ErrOut
      Parameter (ErrOut = 33)

C-----------------------------------------------------------------------

      write(ErrOut,*)' '
      write(ErrOut,'(X,A,I8)')'NEV : ',NEV
      write(ErrOut,'(X,A,F10.2)')Message(1:),Var
      write(ErrOut,*)' '

      RETURN
      END

C-----------------------------------------------------------------------
C-----------------------------------------------------------------------


      SUBROUTINE Calc_Evshapes(Ntrk,Pefo,iret)
C-----------------------------------------------------------------------
C ALPHA user analysis of event shapes and jet rates
C =================================================
C
C Input : Integer  Ntrk         : Number of objects (particles)
C         Real     Pefo(4,Ntrk) : Array containing 4-momenta
C Output: Integer  iret         : 0 = OK
C
C Author: Guenther Dissertori  / 28.5.2001 
C-----------------------------------------------------------------------

      Implicit NONE

C -- LUND COMMON BLOCKS :   K has an extra P here
C -- JETSET part
      Integer    N,KP(4000,5),MSTU(200),MDCY(500,3),MDME(2000,2)
      Integer    KFDP(2000,5),MRLU(6),KCHG(500,3),MSTJ(200)
      Real*4     P(4000,5),V(4000,5),PARU(200),PARJ(200)
      Real*4     PMAS(500,4),PARF(2000),VCKM(4,4)
      Real*4     BRAT(2000),CHAF(500),RRLU(100)

      Common /LUJETS/N,KP,P,V
      Common /LUDAT1/MSTU,PARU,MSTJ,PARJ


      Integer   Ntrk,iret
      Real      Pefo(4,*)

      Logical   Miss
      Integer   i,j,IretGlobal,mylevel
      Integer   Scheme,Error

      Integer   NefoMax
      Parameter(NefoMax=300)
      Real      Pjet(4,NefoMax)
      Real      EvS(20),ThD(4,3),Yn(NefoMax),xth
      Real      Abz,Anr,Aks,A34
      
      Integer   ncut,icut,Njet,iddp
      Parameter(ncut=3)      
      Logical  Is3jet(ncut)
      Real     yc

C -- booking
      Integer    NEV,NEVNT,LEVEL
      Real*4     PI,RADEG,ECM
      Common /MYBOOKING/ PI,RADEG,NEV,NEVNT,ECM,LEVEL

C -- counting
      Integer   Nhad
      Common /Norm/ Nhad

C -- histo definitions
      Include 'histonaming.h'

C-----------------------------------------------------------------------

C -- initialization
      iret = 1
      mylevel = LEVEL

C-- do jet clustering
      Scheme = 1         ! Durham E-scheme
      do icut=1,ncut
        Is3jet(icut) = .False.
        yc = Ycut3j(icut)
        call DurhamJn(Pefo,ntrk,Scheme,yc,Njet,Pjet,Error)
        if (Error.ne.0) then
          call ErrMess('Calc_Evshapes, Error in clustering=',
     >         float(Error))
          IretGlobal = -4
          call HFILL(idc1,float(IretGlobal),0.,1.)
          RETURN  
        else
          if (Njet.eq.3) then
            Is3jet(icut) = .True.
            call HFILL(idjr,float(icut),0.,1.)
          endif
        endif
      enddo


C -- do not consider missing momentum
      Miss = .False.

C -- Calculate Eventshapes for all particles
      call EvShapes(Pefo,ntrk,Miss,EvS,ThD,Error)

      if (Error.ne.0) then
        call ErrMess('Calc_EvShapes: error occured,Error=',
     >                float(Error))
        IretGlobal = -5
        call HFILL(idc1,float(IretGlobal),0.,1.)
        RETURN   ! ------>
      endif

C -- store the angle of the axis     
      call HFILL(idt2,abs(ThD(3,1)),0.,1.)

C -- for -ln(y3) we use this routine
C -- Note : if Level=1 (parton level), it can happen 
C -- that there are only two particles ==> Y3 = 0
      if (ntrk.lt.3) then
        Yn(1) = 1.E-8
      else
        Scheme = 1              ! Durham E-scheme
        call DurhamYn(Pefo,ntrk,Scheme,Yn,Error)
        if (Error.ne.0) then
          write(MSTU(11),*)
          write(MSTU(11),*)' ===> YN problem: ',mylevel,ntrk
          call LULIST(1)
          call ErrMess('Calc_Evshapes, Error in Yn clustering=',
     >         float(Error))
          IretGlobal = -6
          call HFILL(idc1,float(IretGlobal),0.,1.)
          RETURN                ! ----->
        endif
        if (Yn(1).le.0.) then
          call ErrMess('CalcEvshapes, strange Y3 value=',Yn(1))
          IretGlobal = -7
          call HFILL(idc1,float(IretGlobal),0.,1.)
          RETURN                ! ----->
        endif
      endif

C -- bookkeeping of events for which eventshape calc. is ok
      call HFILL(idc2,3.,0.,1.)

C -- fill the histos 
C -- see routine EvShapes for definition of the array EvS
      call HFILL(idth,EvS(1),0.,1.)
      call HFILL(idta,EvS(2),0.,1.)
      call HFILL(idtm,EvS(3),0.,1.)
      call HFILL(idmh,EvS(6),0.,1.)
      call HFILL(idcp,EvS(5),0.,1.)
      call HFILL(idbt,EvS(9),0.,1.)
      call HFILL(idbw,EvS(10),0.,1.)
      call HFILL(idly,-log(Yn(1)),0.,1.)

C -- fill the Dparameter values
      call HFILL(iddi,EvS(11),0.,1.)
      do icut=1,ncut
        if     (icut.eq.1) then
          iddp = idd1
        elseif (icut.eq.2) then
          iddp = idd2
        elseif (icut.eq.3) then
          iddp = idd3
        endif
        if (Is3jet(icut)) then 
          call HFILL(iddp,EvS(11),0.,1.)
        endif
      enddo

C -- fill the moments of 1-Thrust
      xth = 1-EvS(1)
      do i=1,4
        call HFILL(idmo,float(i),0.,xth**i)
        call HFILL(idmo,float(i+4),0.,xth**(2*i))
      enddo

C -- angular correlations in 4-jet events
C -- check first if it is a 4-jet event
      Scheme = 1         ! Durham E-scheme
      yc = Ycut4j
      call DurhamJn(Pefo,ntrk,Scheme,yc,Njet,Pjet,Error)
      if (Error.ne.0) then
        call ErrMess('Calc_Evshapes, Error in 4-jet clustering=',
     >       float(Error))
        IretGlobal = -8
        call HFILL(idc1,float(IretGlobal),0.,1.)
        RETURN   ! ----->
      else
C -- if it is a 4-jet event, compute the angles and fill them
        if (Njet.eq.4) then
          call Fourjets(Pjet,Abz,Aks,Anr,A34,Error)
          if (Error.ne.0) then
            call ErrMess('Calc_Evshapes, Error in 4-jet angles=',
     >       float(Error))
            IretGlobal = -9
            call HFILL(idc1,float(IretGlobal),0.,1.)
            RETURN              ! ----->
          else
            call HFILL(idbz,abs(cos(Abz)),0.,1.)
            call HFILL(idnr,abs(cos(Anr)),0.,1.)
            call HFILL(idks,    cos(Aks) ,0.,1.)
            call HFILL(id34,    cos(A34) ,0.,1.)

C -- bookkeeping of events for which 4-jet angular corr. are OK
            call HFILL(idc2,4.,0.,1.)

          endif
        endif
      endif

C -- OK, everthing done
      iret = 0

 9999 RETURN 
      END

C-----------------------------------------------------------------------
C-----------------------------------------------------------------------


      SUBROUTINE Calc_Inclusives(iret)
C-----------------------------------------------------------------------
C ALPHA user analysis of inclusive charged particle distributions
C ===============================================================
C
C Output: Integer  iret         : 0 = OK
C
C Author: Guenther Dissertori  / 28.5.2001 
C-----------------------------------------------------------------------

      Implicit NONE

C -- LUND COMMON BLOCKS :   K has an extra P here
C -- JETSET part
      Integer    N,KP(4000,5),MSTU(200),MDCY(500,3),MDME(2000,2)
      Integer    KFDP(2000,5),MRLU(6),KCHG(500,3),MSTJ(200)
      Real*4     P(4000,5),V(4000,5),PARU(200),PARJ(200)
      Real*4     PMAS(500,4),PARF(2000),VCKM(4,4)
      Real*4     BRAT(2000),CHAF(500),RRLU(100)

      Common /LUJETS/N,KP,P,V
      Common /LUDAT1/MSTU,PARU,MSTJ,PARJ

      Integer   iret

      Real      eigval,THRUST,Saxis(3,3),Sphe,Apla
      Real      xp,xsi,ptth,psq,yra,pti,pto,xup,xdo
      Integer   i,j,IretGlobal,LUCHGE,itkl
      Logical   Miss
      External  LUCHGE

      Integer   NefoMax
      Parameter(NefoMax=300)
      Real      Pefo(4,NefoMax),ThD(4,3)

C -- booking
      Integer    NEV,NEVNT,LEVEL
      Real*4     PI,RADEG,ECM
      Common /MYBOOKING/ PI,RADEG,NEV,NEVNT,ECM,LEVEL

C -- counting
      Integer   Nhad
      Common /Norm/ Nhad

C -- histo definitions
      Include 'histonaming.h'

C-----------------------------------------------------------------------

C -- initialization
      iret = 1

C -- now look at charged only 
      do i = 1,N      
        if (LUCHGE(kp(i,2)).eq.0) then
          kp(i,1) = 0
        end if
      enddo

C -- select again stable final state particles
C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
C -- NOTE: every analysis after this routine will have 
C          charged particles only!!
C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
      CALL LUEDIT(1)

      IF (NEV.LE.3) CALL LULIST(1)

C -- reject pathological case
      if (N.le.1)  then
        call ErrMess('Calc_Inclusives, Only 1 chg particle! Nev =',
     >       float(NEV))
        IretGlobal = -10
        call HFILL(idc1,float(IretGlobal),0.,1.)
        RETURN   ! ----->
      endif
      
C -- fill them into internal array
      itkl = 0
      do i=1,N
        itkl = itkl+1
        Pefo(1,itkl) = P(i,1)
        Pefo(2,itkl) = P(i,2)
        Pefo(3,itkl) = P(i,3)
        Pefo(4,itkl) = P(i,4)        

C -- protection against abnormally large events
        if (itkl.eq.NefoMax) then
          call ErrMess('Calc_Inclusives, limit NefoMax reached!',
     >         float(itkl))
          IretGlobal = -11
          call HFILL(idc1,float(IretGlobal),0.,1.)
          RETURN ! ----->  
        endif 
      enddo

C -- compute event plane axes (from Sphericity tensor) and Thrust, 
C -- with charged tracks only !!
      MSTU(41) = 1
      call LUSPHE(Sphe,Apla)

      if (Sphe.eq.-1.) then
        call ErrMess('Calc_Inclusives, Sph = -1, Nev =',
     >       float(NEV))
        IretGlobal = -12
        call HFILL(idc1,float(IretGlobal),0.,1.)
        RETURN   ! ----->
      endif

C -- get the sphericity axis
      do i=1,3
        do j=1,3
          Saxis(i,j) = P(N+i,j) 
        enddo
      enddo

C -- compute thrust
      Miss = .False.
      call ThrustDG(Pefo,itkl,ThD,Miss)
      THRUST = ThD(4,1)

      if (THRUST.eq.-1.or.THRUST.eq.-2) then
        call ErrMess('Calc_Inclusives, problem with Thrust, Nev =',
     >       float(NEV))
        IretGlobal = -13
        call HFILL(idc1,float(IretGlobal),0.,1.)
        RETURN  ! ----->
      endif

C -- store the angle of the axis with charged tracks only      
      call HFILL(idt3,abs(ThD(3,1)),0.,1.)


C -- now the loop over charged eflow objects
      do i=1,N

        psq = sqrt(P(i,1)**2+P(i,2)**2+P(i,3)**2)
        xp    = 2*psq/Ecm
        xsi   = -log(xp)
            
C -- momentum component parallel to thrust
        ptth  = P(i,1)*ThD(1,1)+P(i,2)*ThD(2,1)+P(i,3)*ThD(3,1)
        ptth  = abs(ptth)
        
C -- rapidity w.r.t. thrust axis
        xup   = P(i,4) + ptth
        xdo   = P(i,4) - ptth
        if (xdo.eq.0..or.xup.eq.0.) then
          yra = -1.
        else
          yra   = 0.5*log(xup/xdo)
        endif

C -- momenta w.r.t. event plane
        pti = P(i,1)*Saxis(2,1)+P(i,2)*Saxis(2,2)+P(i,3)*Saxis(2,3)
        pto = P(i,1)*Saxis(3,1)+P(i,2)*Saxis(3,2)+P(i,3)*Saxis(3,3)
        pto = abs(pto)
        pti = abs(pti)

C -- fill the histos
        call HFILL(idxp,xp ,0.,1.)
        call HFILL(idxi,xsi,0.,1.)
        call HFILL(idra,yra,0.,1.)
        call HFILL(idpi,pti,0.,1.)
        call HFILL(idpo,pto,0.,1.)
      enddo


C -- OK, everthing done
      iret = 0

 9999 RETURN 
      END
      

C-----------------------------------------------------------------------
C-----------------------------------------------------------------------


C include additional code : the utilities package
C =======================
      Include 'qcd_utilities.F'

C-----------------------------------------------------------------------


