c Write the showered events in stdhep format. Requires the stdhep (and
c mcfio) libraries. The event file will be identical to the Les Houches
c parton-level input event file with ".hep" appended.
C----------------------------------------------------------------------
      SUBROUTINE RCLOS()
C     DUMMY IF HBOOK IS USED
C----------------------------------------------------------------------
      END

C----------------------------------------------------------------------
      SUBROUTINE HWABEG
C     USER'S ROUTINE FOR INITIALIZATION
C----------------------------------------------------------------------
      implicit none
      integer lok,niter,maxevv
      character*80 output_file
      integer istrstd
      common/canalysis/istrstd
      CHARACTER*50 QQIN         ! LH event file
      COMMON/VVJIN/QQIN
      COMMON/CMAXEVV/MAXEVV
c For the stdhep event file, use the input file with ".hep" appended
      output_file=QQIN
      write(output_file(index(QQIN,' '):),'(a)') '.hep'

      NITER = MAXEVV !number of events
C...Open STDHEP output file
      call stdxwinit(output_file,'HERWIG file',NITER,istrstd,lok)
      if (lok.ne.0) then
         print *,'Failed to open output file ',output_file,
     &        '. Quitting!'
         stop
      endif
C   Store Herwig run info in STDHEP common block
      call stdflhwxsec
C   Write STDHEP begin run record
      call stdxwrt(100,istrstd,lok)
      if (lok.ne.0) then
         print '('' Write failure on'', '' begin run record: '',i5)',lok
         stop
      endif
      PRINT *,'stdhep initialized'
 999  END


C----------------------------------------------------------------------
      SUBROUTINE HWAEND
C     USER'S ROUTINE FOR TERMINAL CALCULATIONS, HISTOGRAM OUTPUT, ETC
C----------------------------------------------------------------------
      implicit none
      integer lok
      integer istrstd
      common/canalysis/istrstd
C...Finalize STDHEP file
      call stdflhwxsec
      call stdxwrt(200,istrstd,lok) ! write STDHEP end run record
      if (lok.ne.0) then
         print '('' Write failure on'', '' end run record: '',i5)',lok
      endif
      call stdxend(istrstd)
 999  END

C----------------------------------------------------------------------
      SUBROUTINE HWANAL
C     USER'S ROUTINE TO ANALYSE DATA FROM EVENT
C----------------------------------------------------------------------
      INCLUDE 'HERWIG65.INC'
      integer lok
      integer istrstd
      common/canalysis/istrstd
      INTEGER NMAX
      PARAMETER (NMAX=2000)
      INTEGER i,ist,id,IHEP,NS,n_match,IHARD,NHARD,JMO1,JMO2
      DOUBLE PRECISION PS(4,NMAX),matching_scale
      LOGICAL max_multiplicity,matched
      integer nevhep_save,icount
      data icount /0/
c matching_scale: this is the \mu_Q of the paper (arXiv:1209.6215)
      parameter (matching_scale=30d0)
c n_match identified the i-parton sample of the paper: that's the number
c of partons AT THE BORN LEVEL in a given multiplicity
      parameter (n_match=0)
c set max_multiplicity=.true. if the current sample is that of largest
c multiplicity, and equal to .false. otherwise
      parameter (max_multiplicity=.false.)
c
      IF (IERROR.NE.0) RETURN
      NS=0
      IHARD=-1
      NHARD=0
      DO IHEP=1,min(NHEP,10)
        IF(ISTHEP(IHEP).EQ.120)THEN
          IHARD=IHEP
          NHARD=NHARD+1
        ENDIF
      ENDDO
      IF( (IHARD.EQ.-1.OR.NHARD.NE.1).AND.
     #    IERROR.EQ.0) THEN
         CALL HWUEPR
         WRITE(*,*)IHARD,NHARD
         CALL HWWARN('HWANAL',503)
      ENDIF
      DO 100 IHEP=1,NHEP
         IST=ISTHEP(IHEP)  ! status code of the current particle
         ID =IDHEP(IHEP)   ! PDG-code of the current particle
c Keep QCD partons from primary shower. This already excludes radiation
c from particles that are decayed in Herwig, but not decays that have
c been included via MadSpin.
         IF( IST.EQ.2 .AND. ( (ABS(ID).GE.1.AND.ABS(ID).LE.5) .OR.
     $        ID.EQ.21 ) ) THEN
            JMO1=JMOHEP(1,JMOHEP(1,IHEP))
            JMO2=JMOHEP(2,JMOHEP(1,IHEP))
            IF(ISTHEP(JMO1).GE.121.AND.ISTHEP(JMO1).LE.124.AND.
     $           JMOHEP(1,JMO1).EQ.IHARD.AND.JMO2.EQ.IHARD)THEN
               NS=NS+1
               IF (NS.GT.NMAX) STOP 'Too many particles [merging]!'
               DO I=1,4
                  PS(I,NS)=PHEP(I,IHEP)
               ENDDO
            ENDIF
         ENDIF
  100 CONTINUE
C MLM-type rejection
      call MLM_matching(ps,ns,matching_scale,n_match
     #     ,max_multiplicity, matched)
c If jets after and before shower are not matched, the event has to be
c rejected. Simply exit this subroutine so that the StdHEP file will not
c be filled.
      if (.not.matched) return

c Keep a counter for the number of events passing the MLM-rejection
c criteria.  Update the 'nevlep' accordingly such that the events in the
c StdHEP file will have the correct ID.
      icount=icount+1
      nevhep_save=nevhep
      nevhep=icount
C...Convert herwig event to stdhep format
      call hwghep(1)
C...Write event to STDHEP output file
      call stdxwrt(4,istrstd,lok)
      if (lok.ne.0) then
         print '('' Failure writing event '',i5)',lok
      endif
      nevhep=nevhep_save
 999  END



      subroutine MLM_matching(p_shower,n_shower,matching_scale,n_match
     &     ,max_multiplicity, matched)
      implicit none
C Hard event information:
      INTEGER MAXNUP
      PARAMETER (MAXNUP=500)
      INTEGER NUP,IDPRUP,IDUP,ISTUP,MOTHUP,ICOLUP
      DOUBLE PRECISION XWGTUP,SCALUP,AQEDUP,AQCDUP,PUP,VTIMUP,SPINUP
      COMMON/HEPEUP/NUP,IDPRUP,XWGTUP,SCALUP,AQEDUP,AQCDUP,
     &              IDUP(MAXNUP),ISTUP(MAXNUP),MOTHUP(2,MAXNUP),
     &              ICOLUP(2,MAXNUP),PUP(5,MAXNUP),VTIMUP(MAXNUP),
     &              SPINUP(MAXNUP)
      integer n_hard,i,j
      double precision p_hard(4,MAXNUP)
      INTEGER njet_hard,jet_hard(MAXNUP),n_match,n_shower,njet_shower
     &     ,jet_shower(n_shower)
      DOUBLE PRECISION palg,rfj,sycut,pjet_hard(4,MAXNUP),matching_scale
     &     ,p_shower(4,n_shower),pjet_shower(4,n_shower),getdrv
      logical matched,max_multiplicity,jet_shower_matched
      double precision SHSCALE
      common/cSHSCALE/SHSCALE
      matched=.true.
      palg=1d0
      sycut=matching_scale
      rfj=1.0d0
      call fastjetppgenkt(p_shower,n_shower,rfj,sycut,palg,pjet_shower
     &     ,njet_shower,jet_shower)
c First check that the number of jets after showering is okay if the
c number of jets required in the matching
      If (max_multiplicity) then
         if (njet_shower.lt.n_match) then
            matched=.false.
            return
         endif
      else
         if (njet_shower.ne.n_match) then
            matched=.false.
            return
         endif
      endif
      n_hard=0
      do i=1,nup
         if (ISTUP(i).eq.1.and.(abs(IDUP(i)).lt.6.or.IDUP(i).eq.21))
     &        then
c Final state QCD parton            
            n_hard=n_hard+1
            do j=1,4
               p_hard(j,n_hard)=PUP(j,i)
            enddo
         endif
      enddo
      sycut=0d0
      call fastjetppgenkt(p_hard,n_hard,rfj,sycut,palg,pjet_hard
     &     ,njet_hard,jet_hard)
c Check that the hardest jet after showering 'match' with the hardest
c jets before showering
      do i=1,n_match
         jet_shower_matched=.false.
         do j=1,njet_hard
            if (getdrv(pjet_shower(1,i),pjet_hard(1,j)).lt.1.5d0*rfj)
     &           then
               jet_shower_matched=.true.
               exit
            endif
         enddo
         if (.not.jet_shower_matched) then
            matched=.false.
            return
            exit
         endif
      enddo
c
      return
      end


      function getpseudorap(en,ptx,pty,pl)
      implicit none
      real*8 getpseudorap,en,ptx,pty,pl,tiny,pt,eta,th
      parameter (tiny=1.d-5)
c
      pt=sqrt(ptx**2+pty**2)
      if(pt.lt.tiny.and.abs(pl).lt.tiny)then
        eta=sign(1.d0,pl)*1.d8
      else
        th=atan2(pt,pl)
        eta=-log(tan(th/2.d0))
      endif
      getpseudorap=eta
      return
      end


      function getdelphi(ptx1,pty1,ptx2,pty2)
      implicit none
      real*8 getdelphi,ptx1,pty1,ptx2,pty2,tiny,pt1,pt2,tmp
      parameter (tiny=1.d-5)
c
      pt1=sqrt(ptx1**2+pty1**2)
      pt2=sqrt(ptx2**2+pty2**2)
      if(pt1.ne.0.d0.and.pt2.ne.0.d0)then
        tmp=ptx1*ptx2+pty1*pty2
        tmp=tmp/(pt1*pt2)
        if(abs(tmp).gt.1.d0+tiny)then
          write(*,*)'Cosine larger than 1'
          stop
        elseif(abs(tmp).ge.1.d0)then
          tmp=sign(1.d0,tmp)
        endif
        tmp=acos(tmp)
      else
        tmp=1.d8
      endif
      getdelphi=tmp
      return
      end


      function getdr(en1,ptx1,pty1,pl1,en2,ptx2,pty2,pl2)
      implicit none
      real*8 getdr,en1,ptx1,pty1,pl1,en2,ptx2,pty2,pl2,deta,dphi,
     # getpseudorap,getdelphi
c
      deta=getpseudorap(en1,ptx1,pty1,pl1)-
     #     getpseudorap(en2,ptx2,pty2,pl2)
      dphi=getdelphi(ptx1,pty1,ptx2,pty2)
      getdr=sqrt(dphi**2+deta**2)
      return
      end


      function getdrv(p1,p2)
      implicit none
      real*8 getdrv,p1(4),p2(4)
      real*8 getdr
c
      getdrv=getdr(p1(4),p1(1),p1(2),p1(3),
     #             p2(4),p2(1),p2(2),p2(3))
      return
      end


