jewel-2.2.0.f

Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file jewel-2.2.0.f

C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C++ Copyright (C) 2017 Korinna C. Zapp [Korinna.Zapp@cern.ch]       ++
C++                                                                 ++
C++ This file is part of JEWEL 2.2.0                                ++
C++                                                                 ++
C++ The JEWEL homepage is jewel.hepforge.org                        ++
C++                                                                 ++
C++ The medium model was partly implemented by Jochen Klein.        ++
C++ Raghav Kunnawalkam Elayavalli helped with the implementation    ++
C++ of the V+jet processes.                                         ++
C++                                                                 ++
C++ Please follow the MCnet GUIDELINES and cite Eur.Phys.J. C74     ++
C++ (2014) no.2, 2762 [arXiv:1311.0048] for the code and            ++
C++ JHEP 1303 (2013) 080 [arXiv:1212.1599] and                      ++
C++ optionally EPJC 60 (2009) 617 [arXiv:0804.3568] for the         ++
C++ physics. The reference for V+jet processes is EPJC 76 (2016)    ++
C++ no.12 695 [arXiv:1608.03099] and for recoil effects it is       ++
C++ arXiv:1707.01539.
C++                                                                 ++
C++ JEWEL relies heavily on PYTHIA 6 for the event generation. The  ++
C++ modified version of PYTHIA 6.4.25 that is distributed with      ++
C++ JEWEL is, however, not an official PYTHIA release and must not  ++
C++ be used for anything else. Please refer to results as           ++
C++ "JEWEL+PYTHIA".                                                 ++
C++                                                                 ++
C++ JEWEL also uses code provided by S. Zhang and J. M. Jing        ++
C++ (Computation of Special Functions, John Wiley & Sons, New York, ++
C++ 1996 and http://jin.ece.illinois.edu) for computing the         ++
C++ exponential integral Ei(x).                                     ++
C++                                                                 ++
C++                                                                 ++
C++ JEWEL  is free software; you can redistribute it and/or         ++
C++ modify it under the terms of the GNU General Public License     ++
C++ as published by the Free Software Foundation; either version 2  ++
C++ of the License, or (at your option) any later version.          ++
C++                                                                 ++
C++ JEWEL is distributed in the hope that it will be useful,        ++
C++ but WITHOUT ANY WARRANTY; without even the implied warranty of  ++
C++ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the    ++
C++ GNU General Public License for more details.                    ++
C++                                                                 ++
C++ You should have received a copy of the GNU General Public       ++  
C++ License along with this program; if not, write to the Free      ++
C++ Software Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, ++
C++ MA 02110-1301 USA                                               ++
C++                                                                 ++
C++ Linking JEWEL statically or dynamically with other modules is   ++
C++ making a combined work based on JEWEL. Thus, the terms and      ++
C++ conditions of the GNU General Public License cover the whole    ++
C++ combination.                                                    ++
C++                                                                 ++
C++ In addition, as a special exception, I give you permission to   ++
C++ combine JEWEL with the code for the computation of special      ++
C++ functions provided by S. Zhang and J. M. Jing. You may copy and ++
C++ distribute such a system following the terms of the GNU GPL for ++
C++ JEWEL and the licenses of the other code concerned, provided    ++
C++ that you include the source code of that other code when and as ++
C++ the GNU GPL requires distribution of source code.               ++
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

      PROGRAM jewel
  IMPLICIT NONE
C--Common block of Pythia
      common/pyjets/n,npad,k(23000,5),p(23000,5),v(23000,5)
  INTEGER n,npad,k
  DOUBLE PRECISION p,v
      common/pydat1/mstu(200),paru(200),mstj(200),parj(200)
  INTEGER mstu,mstj
  DOUBLE PRECISION paru,parj
      common/pydat3/mdcy(500,3),mdme(8000,2),brat(8000),kfdp(8000,5)
  INTEGER mdcy,mdme,kfdp
  DOUBLE PRECISION brat
      common/pysubs/msel,mselpd,msub(500),kfin(2,-40:40),ckin(200)
  INTEGER msel,mselpd,msub,kfin
  DOUBLE PRECISION ckin 
      common/pypars/mstp(200),parp(200),msti(200),pari(200)
  INTEGER mstp,msti
  DOUBLE PRECISION parp,pari
      common/pydatr/mrpy(6),rrpy(100)
  INTEGER mrpy
  DOUBLE PRECISION rrpy
C--identifier of file for hepmc output and logfile
  common/hepmcid/hpmcfid,logfid
  integer hpmcfid,logfid
C--use nuclear pdf?      
      common/npdf/mass,nset,eps09,initstr
      INTEGER nset
      DOUBLE PRECISION mass
      LOGICAL eps09
      CHARACTER*10 initstr
C--number of protons
  common/np/nproton
  integer nproton
C--organisation of event record
  common/evrecord/nsim,npart,offset,hadrotype,sqrts,collider,hadro,
     &shorthepmc,channel,isochannel
  integer nsim,npart,offset,hadrotype
  double precision sqrts
  character*4 collider,channel
  character*2 isochannel
  logical hadro,shorthepmc
C--discard event flag
  common/disc/ndisc,nstrange,ngood,errcount,wdisc,discard
  LOGICAL discard
  INTEGER ndisc,nstrange,ngood,errcount
  double precision wdisc
C--event weight
  common/weight/evweight,sumofweights
  double precision evweight,sumofweights
C--number of scattering events
  common/check/nscat,nscateff,nsplit
  DOUBLE PRECISION nscat,nscateff,nsplit
C--number of extrapolations in tables
  common/extrapolations/ntotspliti,noverspliti,ntotpdf,noverpdf,
     &ntotxsec,noverxsec,ntotsuda,noversuda
  integer ntotspliti,noverspliti,ntotpdf,noverpdf,
     &ntotxsec,noverxsec,ntotsuda,noversuda
C--local variables
  integer j,i,kk,poissonian
      integer nsimpp,nsimpn,nsimnp,nsimnn,nsimsum,nsimchn
  double precision sumofweightstot,wdisctot,scalefac
  double precision gettemp,r,tau
  character*2 b1,b2

  call init()

  sumofweightstot=0.d0
      wdisctot=0.d0

C--e+ + e- event generation
  if (collider.eq.'EEJJ') then
    b1 = 'e+'
    b2 = 'e-'
    write(logfid,*)
    write(logfid,*)
     &'####################################################'
    write(logfid,*)
    write(logfid,*)'generating ',nsim,' events in ',b1,' + ',b2,
     &' channel'
    write(logfid,*)
    write(logfid,*)
     &'####################################################'
    write(logfid,*)
    sumofweights=0.d0
        wdisc=0.d0
    call initpythia(b1,b2)
      write(logfid,*)
C--e+ + e- event loop
    DO 100 j=1,nsim
      call genevent(j,b1,b2)
 100    CONTINUE
    sumofweightstot = sumofweightstot+sumofweights
    wdisctot = wdisctot + wdisc
    write(logfid,*)
    write(logfid,*)'cross section in e+ + e- channel:',pari(1),'mb'
    write(logfid,*)'sum of event weights in e+ + e- channel:',
     &  sumofweights-wdisc
    write(logfid,*)

  else
C--hadronic event generation
    if (isochannel.eq.'PP') then
      nsimpp = nsim
      nsimpn = 0
      nsimnp = 0
      nsimnn = 0
    elseif (isochannel.eq.'PN') then
      nsimpp = 0
      nsimpn = nsim
      nsimnp = 0
      nsimnn = 0
    elseif (isochannel.eq.'NP') then
      nsimpp = 0
      nsimpn = 0
      nsimnp = nsim
      nsimnn = 0
    elseif (isochannel.eq.'NN') then
      nsimpp = 0
      nsimpn = 0
      nsimnp = 0
      nsimnn = nsim
    else
      nsimpp = poissonian(nsim*nproton**2/mass**2)
      nsimpn = poissonian(nsim*nproton*(mass-nproton*1.d0)/mass**2)
      nsimnp = poissonian(nsim*nproton*(mass-nproton*1.d0)/mass**2)
      nsimnn = poissonian(nsim*(mass-nproton*1.d0)**2/mass**2)
      nsimsum = nsimpp + nsimpn + nsimnp + nsimnn
      scalefac = nsim*1.d0/(nsimsum*1.d0)
      nsimpp = int(nsimpp*scalefac)
      nsimpn = int(nsimpn*scalefac)
      nsimnp = int(nsimnp*scalefac)
      nsimnn = int(nsimnn*scalefac)
      nsimsum = nsimpp + nsimpn + nsimnp + nsimnn
    endif
C--loop over channels
    do 101 kk=1,4
      if (kk.eq.1) then
        b1 = 'p+'
        b2 = 'p+'
        nsimchn = nsimpp
      elseif (kk.eq.2) then
        b1 = 'p+'
        b2 = 'n0'
        nsimchn = nsimpn
      elseif (kk.eq.3) then
        b1 = 'n0'
        b2 = 'p+'
        nsimchn = nsimnp
      else
        b1 = 'n0'
        b2 = 'n0'
        nsimchn = nsimnn
      endif
      write(logfid,*)
      write(logfid,*)
     &'####################################################'
      write(logfid,*)
      write(logfid,*)'generating ',nsimchn,' events in ',
     &b1,' + ',b2,' channel'
      write(logfid,*)
      write(logfid,*)
     &'####################################################'
      write(logfid,*)
      sumofweights=0.d0
          wdisc=0.d0
      call initpythia(b1,b2)
      write(logfid,*)
C--event loop
      DO 102 j=1,nsimchn
        call genevent(j,b1,b2)
 102      CONTINUE
      sumofweightstot = sumofweightstot+sumofweights
      wdisctot = wdisctot + wdisc
      write(logfid,*)
      write(logfid,*)'cross section in ',b1,' + ',b2,' channel:',
     &  pari(1),'mb'
      write(logfid,*)'sum of event weights in ',b1,' + ',b2,
     &  ' channel:',sumofweights-wdisc
      write(logfid,*)
 101    continue
  endif
 
C--finish
  WRITE(hpmcfid,'(A)')'HepMC::IO_GenEvent-END_EVENT_LISTING'
  WRITE(hpmcfid,*)
  CLOSE(hpmcfid,status='keep')

  write(logfid,*)
  write(logfid,*)'mean number of scatterings:',
     &      nscat/(sumofweightstot-wdisctot)
  write(logfid,*)'mean number of effective scatterings:',
     &      nscateff/(sumofweightstot-wdisctot)
  write(logfid,*)'mean number of splittings:',
     &      nsplit/(sumofweightstot-wdisctot)
  write(logfid,*)
  write(logfid,*)'number of extrapolations in splitting integral: ',
     &  noverspliti,' (',(noverspliti*1.d0)/(ntotspliti*1.d0),'%)'
  write(logfid,*)
     &  'number of extrapolations in splitting partonic PDFs: ',
     &  noverpdf,' (',(noverpdf*1.d0)/(ntotpdf*1.d0),'%)'
  write(logfid,*)
     &  'number of extrapolations in splitting cross sections: ',
     &  noverxsec,' (',(noverxsec*1.d0)/(ntotxsec*1.d0),'%)'
  write(logfid,*)
     &  'number of extrapolations in Sudakov form factor: ',
     &  noversuda,' (',(noversuda*1.d0)/(ntotsuda*1.d0),'%)'
  write(logfid,*)
  write(logfid,*)'number of good events: ',ngood
  write(logfid,*)'total number of discarded events: ',ndisc
  write(logfid,*)'number of events for which conversion '//
     &'to hepmc failed: ',nstrange
  call printtime

  close(logfid,status='keep')

  END



***********************************************************************
***********************************************************************
***   END OF MAIN PROGRAM - NOW COME THE SUBROUTINES   ****************
***********************************************************************
***********************************************************************


***********************************************************************
***   subroutine init
***********************************************************************
  subroutine init()
  implicit none
  INTEGER pycomp
  INTEGER nmxhep
C--Common block of Pythia
      common/pyjets/n,npad,k(23000,5),p(23000,5),v(23000,5)
  INTEGER n,npad,k
  DOUBLE PRECISION p,v
      common/pydat1/mstu(200),paru(200),mstj(200),parj(200)
  INTEGER mstu,mstj
  DOUBLE PRECISION paru,parj
      common/pydat3/mdcy(500,3),mdme(8000,2),brat(8000),kfdp(8000,5)
  INTEGER mdcy,mdme,kfdp
  DOUBLE PRECISION brat
      common/pysubs/msel,mselpd,msub(500),kfin(2,-40:40),ckin(200)
  INTEGER msel,mselpd,msub,kfin
  DOUBLE PRECISION ckin 
      common/pypars/mstp(200),parp(200),msti(200),pari(200)
  INTEGER mstp,msti
  DOUBLE PRECISION parp,pari
      common/pydatr/mrpy(6),rrpy(100)
  INTEGER mrpy
  DOUBLE PRECISION rrpy
C--use nuclear pdf?      
      common/npdf/mass,nset,eps09,initstr
      INTEGER nset
      DOUBLE PRECISION mass
      LOGICAL eps09
      CHARACTER*10 initstr
C--pdfset
  common/pdf/pdfset
  integer pdfset
C--number of protons
  common/np/nproton
  integer nproton
C--Parameter common block
  common/param/q0,lps,lqcd,ltime,scalefacm,angord,scatrecoil,
     &allhad,compress,nf
      INTEGER nf
  DOUBLE PRECISION q0,lqcd,ltime,lps,scalefacm
      LOGICAL angord,scatrecoil,allhad,compress
C--splitting integral
      common/splitint/splitiggv(1000,1000),splitiqqv(1000,1000),
     &splitiqgv(1000,1000),qval(1000),zmval(1000),qmax,zmmin,npoint
      INTEGER npoint
      DOUBLE PRECISION splitiggv,splitiqqv,splitiqgv,
     &qval,zmval,qmax,zmmin
C--pdf common block
  common/pdfs/qinqx(2,1000),ginqx(2,1000),qingx(2,1000),
     &gingx(2,1000)
  DOUBLE PRECISION qinqx,ginqx,qingx,gingx
C--cross secttion common block
  common/xsecs/intq1(1001,101),intq2(1001,101),
     &intg1(1001,101),intg2(1001,101)
  DOUBLE PRECISION intq1,intq2,intg1,intg2
C--Sudakov common block
  common/insuda/sudaqq(1000,2),sudaqg(1000,2),sudagg(1000,2)
     &,sudagc(1000,2)
  DOUBLE PRECISION sudaqq,sudaqg,sudagg,sudagc
C--exponential integral for negative arguments
      common/expint/eix(3,1000),valmax,nval
      INTEGER nval
      DOUBLE PRECISION eix,valmax
C--discard event flag
  common/disc/ndisc,nstrange,ngood,errcount,wdisc,discard
  LOGICAL discard
  INTEGER ndisc,nstrange,ngood,errcount
  double precision wdisc
C--factor in front of formation times
  common/ftimefac/ftfac
  DOUBLE PRECISION ftfac
C--factor in front of alphas argument
  common/alphasfac/ptfac
  DOUBLE PRECISION ptfac
C--number of scattering events
  common/check/nscat,nscateff,nsplit
  DOUBLE PRECISION nscat,nscateff,nsplit
C--number of extrapolations in tables
  common/extrapolations/ntotspliti,noverspliti,ntotpdf,noverpdf,
     &ntotxsec,noverxsec,ntotsuda,noversuda
  integer ntotspliti,noverspliti,ntotpdf,noverpdf,
     &ntotxsec,noverxsec,ntotsuda,noversuda
C--event weight
  common/weight/evweight,sumofweights
  double precision evweight,sumofweights
C--event weight exponent
  common/wexpo/weightex
  DOUBLE PRECISION weightex
C--identifier of file for hepmc output and logfile
  common/hepmcid/hpmcfid,logfid
  integer hpmcfid,logfid
C--max rapidity
  common/rapmax/etamax
  double precision etamax
C--memory for error message from getdeltat
  common/errline/errl
  integer errl
C--organisation of event record
  common/evrecord/nsim,npart,offset,hadrotype,sqrts,collider,hadro,
     &shorthepmc,channel,isochannel
  integer nsim,npart,offset,hadrotype
  double precision sqrts
  character*4 collider,channel
  character*2 isochannel
  logical hadro,shorthepmc
C--extra storage for scattering centres before interactions
      common/storescatcen/nscatcen,maxnscatcen,scatflav(10000),
     &scatcen(10000,5),writescatcen,writedummies
  integer nscatcen,maxnscatcen,scatflav
  double precision scatcen
  logical writescatcen,writedummies
C--Pythia parameters
  common/pythiaparams/ptmin,ptmax,weighted
  double precision ptmin,ptmax
  LOGICAL weighted

C--Variables local to this program
  INTEGER njob,ios,pos,i,j,jj,intmass
  DOUBLE PRECISION getltimemax,eovest,r,pyr
  character firstchar
  CHARACTER*2 snset
      CHARACTER*80 pdffile,xsecfile,filemed,filesplit,buffer,
     &label,value
      CHARACTER*100 hepmcfile,logfile,filename2
  CHARACTER(LEN=100) filename
  LOGICAL pdfexist,splitiexist,xsecexist

  data maxnscatcen/10000/

      hpmcfid = 4
  logfid = 3

C--default settings
  nsim = 10
  njob = 0
  logfile = 'out.log'
  hepmcfile = 'out.hepmc'
  filesplit = 'splitint.dat'
  pdffile = 'pdfs.dat'
  xsecfile = 'xsecs.dat'
  filemed = 'medium-params.dat'
  nf = 3
  lqcd = 0.4
  q0 = 1.5
  ptmin = 5.
  ptmax = 350.
  etamax = 3.1
  collider = 'PPJJ'
  isochannel = 'XX'
  channel = 'MUON'
  sqrts = 2760
  pdfset = 10042
  nset = 1
  mass = 208.
      nproton = 82
  weighted = .true.
  weightex = 5.
  angord = .true.
  allhad = .false.
  hadro = .true.
  hadrotype = 0
  shorthepmc = .true.
  compress = .true.
  writescatcen = .false.
  writedummies = .false.
  
  lps = lqcd
  scatrecoil = .false.
  if (.not.hadro) shorthepmc = .true.

  scalefacm=1.
  ptfac=1.
  ftfac=1.d0

  if (iargc().eq.0) then
    write(*,*)'No parameter file given, '// 
     &'will run with default settings.'
  else
    call getarg(1,filename)
    write(*,*)'Reading parameters from ',filename
    open(unit=1,file=filename,status='old',err=110)
    do 120 i=1,1000
          read(1, '(A)', iostat=ios) buffer
      if(ios.ne.0) goto 130
      firstchar = buffer(1:1)
      if (firstchar.eq.'#') goto 120
          pos=scan(buffer,' ')
          label=buffer(1:pos)
          value=buffer(pos+1:)
          if(label.eq."NEVENT")then
            read(value,*,iostat=ios) nsim
          elseif(label.eq."NJOB")then
            read(value,*,iostat=ios) njob
          elseif(label.eq."LOGFILE")then
            read(value,'(a)',iostat=ios) logfile
          elseif(label.eq."HEPMCFILE")then
            read(value,'(a)',iostat=ios) hepmcfile
          elseif(label.eq."SPLITINTFILE")then
            read(value,'(a)',iostat=ios) filesplit
          elseif(label.eq."PDFFILE")then
            read(value,'(a)',iostat=ios) pdffile
          elseif(label.eq."XSECFILE")then
            read(value,'(a)',iostat=ios) xsecfile
          elseif(label.eq."MEDIUMPARAMS")then
            read(value,'(a)',iostat=ios) filemed
          elseif(label.eq."NF")then
            read(value,*,iostat=ios) nf
          elseif(label.eq."LAMBDAQCD")then
            read(value,*,iostat=ios) lqcd
          elseif(label.eq."Q0")then
            read(value,*,iostat=ios) q0
          elseif(label.eq."PTMIN")then
            read(value,*,iostat=ios) ptmin
          elseif(label.eq."PTMAX")then
            read(value,*,iostat=ios) ptmax
          elseif(label.eq."ETAMAX")then
            read(value,*,iostat=ios) etamax
          elseif(label.eq."PROCESS")then
            read(value,*,iostat=ios) collider
          elseif(label.eq."ISOCHANNEL")then
            read(value,*,iostat=ios) isochannel
      elseif(label.eq."CHANNEL")then
      read(value,*,iostat=ios) channel
          elseif(label.eq."SQRTS")then
            read(value,*,iostat=ios) sqrts
          elseif(label.eq."PDFSET")then
            read(value,*,iostat=ios) pdfset
          elseif(label.eq."NSET")then
            read(value,*,iostat=ios) nset
          elseif(label.eq."MASS")then
            read(value,*,iostat=ios) mass
          elseif(label.eq."NPROTON")then
            read(value,*,iostat=ios) nproton
          elseif(label.eq."WEIGHTED")then
            read(value,*,iostat=ios) weighted
          elseif(label.eq."WEXPO")then
            read(value,*,iostat=ios) weightex
          elseif(label.eq."ANGORD")then
            read(value,*,iostat=ios) angord
          elseif(label.eq."KEEPRECOILS")then
            read(value,*,iostat=ios) allhad
          elseif(label.eq."HADRO")then
            read(value,*,iostat=ios) hadro
          elseif(label.eq."HADROTYPE")then
            read(value,*,iostat=ios) hadrotype
          elseif(label.eq."SHORTHEPMC")then
            read(value,*,iostat=ios) shorthepmc
          elseif(label.eq."COMPRESS")then
            read(value,*,iostat=ios) compress
          elseif(label.eq."WRITESCATCEN")then
            read(value,*,iostat=ios) writescatcen
          elseif(label.eq."WRITEDUMMIES")then
            read(value,*,iostat=ios) writedummies
      else
        write(*,*)'unknown label ',label
      endif
 120    continue


 110    write(*,*)
     &    'Unable to open parameter file, will exit the run.'
    call exit(1)

 130    close(1,status='keep')
    write(*,*)'...done'
  endif

  if (ptmin.lt.3.d0) ptmin = 3.d0
  if (.not.writescatcen) writedummies = .false.

  OPEN(unit=logfid,file=logfile,status='unknown')
  mstu(11)=logfid

  call printtime
  call printlogo(logfid)


  write(logfid,*)
  write(logfid,*)'parameters of the run:'
  write(logfid,*)'NEVENT       = ',nsim
  write(logfid,*)'NJOB         = ',njob
  write(logfid,*)'LOGFILE      = ',logfile
  write(logfid,*)'HEPMCFILE    = ',hepmcfile
  write(logfid,*)'SPLITINTFILE = ',filesplit
  write(logfid,*)'PDFFILE      = ',pdffile
  write(logfid,*)'XSECFILE     = ',xsecfile
  write(logfid,*)'MEDIUMPARAMS = ',filemed
  write(logfid,*)'NF           = ',nf
  write(logfid,*)'LAMBDAQCD    = ',lqcd
  write(logfid,*)'Q0           = ',q0
  write(logfid,*)'PTMIN        = ',ptmin
  write(logfid,*)'PTMAX        = ',ptmax
  write(logfid,*)'ETAMAX       = ',etamax
  write(logfid,*)'PROCESS      = ',collider
  write(logfid,*)'ISOCHANNEL   = ',isochannel
  write(logfid,*)'CHANNEL      = ',channel
  write(logfid,*)'SQRTS        = ',sqrts
  write(logfid,*)'PDFSET       = ',pdfset
  write(logfid,*)'NSET         = ',nset
  write(logfid,*)'MASS         = ',mass
  write(logfid,*)'NPROTON      = ',nproton
  write(logfid,*)'WEIGHTED     = ',weighted
  write(logfid,*)'WEXPO        = ',weightex
  write(logfid,*)'ANGORD       = ',angord
  write(logfid,*)'KEEPRECOILS  = ',allhad
  write(logfid,*)'HADRO        = ',hadro
  write(logfid,*)'HADROTYPE    = ',hadrotype
  write(logfid,*)'SHORTHEPMC   = ',shorthepmc
  write(logfid,*)'COMPRESS     = ',compress
  write(logfid,*)'WRITESCATCEN = ',writescatcen
  write(logfid,*)'WRITEDUMMIES = ',writedummies
  write(logfid,*)
  call flush(logfid)

  if ((collider.ne.'PPJJ').and.(collider.ne.'EEJJ')
     &  .and.(collider.ne.'PPYJ').and.(collider.ne.'PPYQ')
     &  .and.(collider.ne.'PPYG')
     &  .and.(collider.ne.'PPZJ').and.(collider.ne.'PPZQ')
     &  .and.(collider.ne.'PPZG').and.(collider.ne.'PPWJ')
     &  .and.(collider.ne.'PPWQ').and.(collider.ne.'PPWG')
     &      .and.(collider.ne.'PPDY')) then
    write(logfid,*)'Fatal error: colliding system unknown, '//
     &  'will exit now'
    call exit(1)
  endif

C--initialize medium
  intmass = int(mass)
      CALL medinit(filemed,logfid,etamax,intmass)
      CALL mednextevt

  OPEN(unit=hpmcfid,file=hepmcfile,status='unknown')
  WRITE(hpmcfid,*)
  WRITE(hpmcfid,'(A)')'HepMC::Version 2.06.05'
  WRITE(hpmcfid,'(A)')'HepMC::IO_GenEvent-START_EVENT_LISTING'

  npart=2
  
  if(ptmax.gt.0.)then
    eovest=min(1.5*(ptmax+50.)*cosh(etamax),sqrts/2.)
  else
    eovest=sqrts/2.
  endif

  
  CALL eixint
  CALL insudaint(eovest)

  write(logfid,*)
   INQUIRE(file=filesplit,exist=splitiexist)
   IF(splitiexist)THEN
    write(logfid,*)'read splitting integrals from ',filesplit
    OPEN(unit=10,file=filesplit,status='old')
    READ(10,*)qmax,zmmin,npoint
    DO 893 i=1,npoint+1
     READ(10,*) qval(i),zmval(i)
 893    CONTINUE   
    DO 891 i=1,npoint+1
     DO 892 j=1,npoint+1
      READ(10,*)splitiggv(i,j),splitiqqv(i,j),splitiqgv(i,j)
 892     CONTINUE
 891    CONTINUE
    CLOSE(10,status='keep')
   ELSE
    write(logfid,*)'have to integrate splitting functions, '// 
     &'this may take some time'
    CALL splitfncint(eovest)
    INQUIRE(file=filesplit,exist=splitiexist)
    IF(.NOT.splitiexist)THEN
     write(logfid,*)'write splitting integrals to ',filesplit
     OPEN(unit=10,file=filesplit,status='new')
     WRITE(10,*)qmax,zmmin,npoint
     DO 896 i=1,npoint+1
      WRITE(10,*) qval(i),zmval(i)
 896     CONTINUE  
     DO 897 i=1,npoint+1
      DO 898 j=1,npoint+1
       WRITE(10,*)splitiggv(i,j),splitiqqv(i,j),splitiqgv(i,j)
 898      CONTINUE
 897     CONTINUE
     CLOSE(10,status='keep')
    ENDIF 
   ENDIF
  write(logfid,*)

  INQUIRE(file=pdffile,exist=pdfexist)
  IF(pdfexist)THEN
  write(logfid,*)'read pdfs from ',pdffile
   OPEN(unit=10,file=pdffile,status='old')
   DO 872 i=1,2
    DO 873 j=1,1000
     READ(10,*)qinqx(i,j),ginqx(i,j),qingx(i,j),gingx(i,j)
 873    CONTINUE
 872   CONTINUE
   CLOSE(10,status='keep')
  ELSE
   write(logfid,*)'have to integrate pdfs, this may take some time'
   CALL pdfint(eovest)
   INQUIRE(file=pdffile,exist=pdfexist)
   IF(.NOT.pdfexist)THEN
    write(logfid,*)'write pdfs to ',pdffile
    OPEN(unit=10,file=pdffile,status='new')
    DO 876 i=1,2
     DO 877 j=1,1000
      WRITE(10,*)qinqx(i,j),ginqx(i,j),qingx(i,j),gingx(i,j)
 877     CONTINUE
 876    CONTINUE
    CLOSE(10,status='keep')
   ENDIF
  ENDIF 
  write(logfid,*)

  INQUIRE(file=xsecfile,exist=xsecexist)
  IF(xsecexist)THEN
  write(logfid,*)'read cross sections from ',xsecfile
   OPEN(unit=10,file=xsecfile,status='old')
    DO 881 j=1,1001
         DO 885 jj=1,101
     READ(10,*)intq1(j,jj),intq2(j,jj),
     &intg1(j,jj),intg2(j,jj)
 885     CONTINUE
 881    CONTINUE
   CLOSE(10,status='keep')
  ELSE
   write(logfid,*)'have to integrate cross sections, '//
     &'this may take some time'
   CALL xsecint(eovest)
   INQUIRE(file=xsecfile,exist=xsecexist)
   IF(.NOT.xsecexist)THEN
    write(logfid,*)'write cross sections to ',xsecfile
    OPEN(unit=10,file=xsecfile,status='new')
     DO 883 j=1,1001
          DO 884 jj=1,101
      WRITE(10,*)intq1(j,jj),intq2(j,jj),
     &intg1(j,jj),intg2(j,jj)
 884      CONTINUE
 883     CONTINUE
    CLOSE(10,status='keep')
   ENDIF 
  ENDIF
  write(logfid,*)
  CALL flush(3)



C--initialise random number generator status
      IF(njob.GT.0)THEN
       mrpy(1)=njob*1000
       mrpy(2)=0
      ENDIF

C--Call PYR once for initialization
  r=pyr(0)

  ndisc=0
      ngood=0
      nstrange=0
      
  errcount=0
  errl = 0

  nscat=0.d0
  nscateff=0.d0
  nsplit=0.d0

  ntotspliti=0
  noverspliti=0
  ntotpdf=0
  noverpdf=0
  ntotxsec=0
  noverxsec=0
  ntotsuda=0
  noversuda=0

  IF(nset.EQ.0)THEN
   eps09=.false.
  ELSE
   eps09=.true.
   IF(nset.LT.10)THEN
    WRITE(snset,'(i1)') nset
   ELSE
    WRITE(snset,'(i2)') nset
   ENDIF
    initstr='EPS09LO,'//snset
  ENDIF 

  end



***********************************************************************
***   subroutine initpythia
***********************************************************************
  subroutine initpythia(beam1,beam2)
  implicit none
  INTEGER pycomp
  INTEGER nmxhep
C--Common block of Pythia
      common/pyjets/n,npad,k(23000,5),p(23000,5),v(23000,5)
  INTEGER n,npad,k
  DOUBLE PRECISION p,v
      common/pydat1/mstu(200),paru(200),mstj(200),parj(200)
  INTEGER mstu,mstj
  DOUBLE PRECISION paru,parj
      common/pydat3/mdcy(500,3),mdme(8000,2),brat(8000),kfdp(8000,5)
  INTEGER mdcy,mdme,kfdp
  DOUBLE PRECISION brat
      common/pysubs/msel,mselpd,msub(500),kfin(2,-40:40),ckin(200)
  INTEGER msel,mselpd,msub,kfin
  DOUBLE PRECISION ckin 
      common/pypars/mstp(200),parp(200),msti(200),pari(200)
  INTEGER mstp,msti
  DOUBLE PRECISION parp,pari
      common/pydatr/mrpy(6),rrpy(100)
  INTEGER mrpy
  DOUBLE PRECISION rrpy
C--use nuclear pdf?      
      common/npdf/mass,nset,eps09,initstr
      INTEGER nset
      DOUBLE PRECISION mass
      LOGICAL eps09
      CHARACTER*10 initstr
C--pdfset
  common/pdf/pdfset
  integer pdfset
C--Parameter common block
  common/param/q0,lps,lqcd,ltime,scalefacm,angord,scatrecoil,
     &allhad,compress,nf
      INTEGER nf
  DOUBLE PRECISION q0,lqcd,ltime,lps,scalefacm
      LOGICAL angord,scatrecoil,allhad,compress
C--discard event flag
  common/disc/ndisc,nstrange,ngood,errcount,wdisc,discard
  LOGICAL discard
  INTEGER ndisc,nstrange,ngood,errcount
  double precision wdisc
C--event weight
  common/weight/evweight,sumofweights
  double precision evweight,sumofweights
C--event weight exponent
  common/wexpo/weightex
  DOUBLE PRECISION weightex
C--memory for error message from getdeltat
  common/errline/errl
  integer errl
C--organisation of event record
  common/evrecord/nsim,npart,offset,hadrotype,sqrts,collider,hadro,
     &shorthepmc,channel,isochannel
  integer nsim,npart,offset,hadrotype
  double precision sqrts
  character*4 collider,channel
  character*2 isochannel
  logical hadro,shorthepmc
C--Pythia parameters
  common/pythiaparams/ptmin,ptmax,weighted
  double precision ptmin,ptmax
  LOGICAL weighted

C--Variables local to this program
  character*2 beam1,beam2


C--initialise PYTHIA
C--no multiple interactions
   mstp(81) = 0
C--initial state radiation
   mstp(61)=1
C--switch off final state radiation
   mstp(71)=0
C--No hadronisation (yet)
       mstp(111)=0
C--parameter affecting treatment of string corners
       paru(14)=1.
C--Min shat in simulation
       ckin(1)=2.      
C--pT-cut
       ckin(3)=ptmin
       ckin(4)=ptmax
C--use LHAPDF
   mstp(52)=2
C--choose pdf: CTEQ6ll (LO fit/LO alphas) - 10042
C          MSTW2008 (LO central) - 21000
   mstp(51)=pdfset
   IF(collider.EQ.'PPYQ')THEN
    msel=0
    msub(29)=1
   ELSEIF(collider.EQ.'PPYG')THEN
    msel=0
    msub(14)=1
    msub(115)=1
   ELSEIF(collider.EQ.'PPYJ')THEN
    msel=0
    msub(14)=1
    msub(29)=1
    msub(115)=1
   ELSEIF((collider.EQ.'PPZJ').or.(collider.EQ.'PPZQ')
     &  .or.(collider.EQ.'PPZG')
     &      .or.(collider.eq.'PPDY'))THEN
    msel=0
    IF((collider.EQ.'PPZJ').or.(collider.EQ.'PPZQ')) msub(30)=1
    IF((collider.EQ.'PPZJ').or.(collider.EQ.'PPZG')) msub(15)=1
    IF(collider.EQ.'PPDY') msub(1)=1
    mdme(174,1)=0          !Z decay into d dbar', 
    mdme(175,1)=0          !Z decay into u ubar', 
    mdme(176,1)=0          !Z decay into s sbar', 
    mdme(177,1)=0          !Z decay into c cbar', 
    mdme(178,1)=0          !Z decay into b bbar', 
    mdme(179,1)=0          !Z decay into t tbar', 
    mdme(182,1)=0          !Z decay into e- e+', 
    mdme(183,1)=0          !Z decay into nu_e nu_ebar', 
    mdme(184,1)=0          !Z decay into mu- mu+', 
    mdme(185,1)=0          !Z decay into nu_mu nu_mubar', 
    mdme(186,1)=0          !Z decay into tau- tau+', 
    mdme(187,1)=0          !Z decay into nu_tau nu_taubar',
    if (channel.EQ.'ELEC')THEN
      mdme(182,1)=1
    ELSEIF(channel.EQ.'MUON')THEN
      mdme(184,1)=1
    ENDIF
   ELSEIF((collider.EQ.'PPWJ').or.(collider.EQ.'PPWQ')
     &  .or.(collider.EQ.'PPWG'))THEN
    msel=0
    IF((collider.EQ.'PPWJ').or.(collider.EQ.'PPWQ')) msub(31)=1
    IF((collider.EQ.'PPWJ').or.(collider.EQ.'PPWG')) msub(16)=1
    mdme(190,1)=0          ! W+ decay into dbar u,
    mdme(191,1)=0          ! W+ decay into dbar c,
    mdme(192,1)=0          ! W+ decay into dbar t,
    mdme(194,1)=0          ! W+ decay into sbar u,
    mdme(195,1)=0          ! W+ decay into sbar c,
    mdme(196,1)=0          ! W+ decay into sbar t,
    mdme(198,1)=0          ! W+ decay into bbar u,
    mdme(199,1)=0          ! W+ decay into bbar c,
    mdme(200,1)=0          ! W+ decay into bbar t,
    mdme(202,1)=0          ! W+ decay into b'bar u,
    mdme(203,1)=0          ! W+ decay into b'bar c,
    mdme(204,1)=0          ! W+ decay into b'bar t,
    mdme(206,1)=0          ! W+ decay into e+ nu_e,
    mdme(207,1)=0          ! W+ decay into mu+ nu_mu,
    mdme(208,1)=0          ! W+ decay into tau+ nu_tau,
    mdme(209,1)=0      ! W+ decay into tau'+ nu'_tau,
    if (channel.EQ.'ELEC')THEN
     mdme(206,1)=1
    ELSEIF(channel.EQ.'MUON')THEN
     mdme(207,1)=1
    ENDIF
   ELSE
C--All QCD processes are active
        msel=1
   ENDIF
!  MSEL=0
!  MSUB(11)=1
!  MSUB(12)=1
!  MSUB(53)=1
!  MSUB(13)=1
!  MSUB(68)=1
!  MSUB(28)=1

C--weighted events
       IF(weighted) mstp(142)=1

C--number of errors to be printed
   mstu(22)=max(10,int(5.*nsim/100.))

C--number of lines in event record
  mstu(4)=23000
  mstu(5)=23000

C--switch off pi0 decay
      mdcy(pycomp(111),1)=0
C--initialisation call
   IF(collider.EQ.'EEJJ')THEN
    offset=9
    CALL pyinit('CMS',beam1,beam2,sqrts)
   ELSEIF((collider.EQ.'PPJJ').OR.(collider.EQ.'PPYJ').OR.
     &    (collider.EQ.'PPYG').OR.(collider.EQ.'PPYQ'))THEN
    offset=8
    CALL pyinit('CMS',beam1,beam2,sqrts)
   ELSEIF((collider.EQ.'PPWJ').OR.(collider.EQ.'PPZJ').or.
     &  (collider.EQ.'PPWQ').OR.(collider.EQ.'PPZQ').or.
     &  (collider.EQ.'PPWG').OR.(collider.EQ.'PPZG'))THEN
    offset=10
    CALL pyinit('CMS',beam1,beam2,sqrts)
   elseif (collider.eq.'PPDY') then
    CALL pyinit('CMS',beam1,beam2,sqrts)
   ENDIF

  end



***********************************************************************
***   subroutine genevent
***********************************************************************
  subroutine genevent(j,b1,b2)
  implicit none
C--identifier of file for hepmc output and logfile
  common/hepmcid/hpmcfid,logfid
  integer hpmcfid,logfid
  INTEGER pycomp
  INTEGER nmxhep
C--Common block of Pythia
      common/pyjets/n,npad,k(23000,5),p(23000,5),v(23000,5)
  INTEGER n,npad,k
  DOUBLE PRECISION p,v
      common/pydat1/mstu(200),paru(200),mstj(200),parj(200)
  INTEGER mstu,mstj
  DOUBLE PRECISION paru,parj
      common/pydat3/mdcy(500,3),mdme(8000,2),brat(8000),kfdp(8000,5)
  INTEGER mdcy,mdme,kfdp
  DOUBLE PRECISION brat
      common/pysubs/msel,mselpd,msub(500),kfin(2,-40:40),ckin(200)
  INTEGER msel,mselpd,msub,kfin
  DOUBLE PRECISION ckin 
      common/pypars/mstp(200),parp(200),msti(200),pari(200)
  INTEGER mstp,msti
  DOUBLE PRECISION parp,pari
      common/pydatr/mrpy(6),rrpy(100)
  INTEGER mrpy
  DOUBLE PRECISION rrpy
C--Parameter common block
  common/param/q0,lps,lqcd,ltime,scalefacm,angord,scatrecoil,
     &allhad,compress,nf
      INTEGER nf
  DOUBLE PRECISION q0,lqcd,ltime,lps,scalefacm
      LOGICAL angord,scatrecoil,allhad,compress
C--discard event flag
  common/disc/ndisc,nstrange,ngood,errcount,wdisc,discard
  LOGICAL discard
  INTEGER ndisc,nstrange,ngood,errcount
  double precision wdisc
C--variables for angular ordering
      common/angor/za(23000),zd(23000),thetaa(23000),qqbard(23000)
  DOUBLE PRECISION za,zd,thetaa
      LOGICAL qqbard
C--factor in front of formation times
  common/ftimefac/ftfac
  DOUBLE PRECISION ftfac
C--time common block
      common/time/mv(23000,5)
      DOUBLE PRECISION mv
C--colour index common block
  common/colour/trip(23000),anti(23000),colmax
  INTEGER trip,anti,colmax
C--number of scattering events
  common/check/nscat,nscateff,nsplit
  DOUBLE PRECISION nscat,nscateff,nsplit
C--event weight
  common/weight/evweight,sumofweights
  double precision evweight,sumofweights
C--event weight exponent
  common/wexpo/weightex
  DOUBLE PRECISION weightex
C--max rapidity
  common/rapmax/etamax
  double precision etamax
C--production point
  common/jetpoint/x0,y0
  double precision x0,y0
C--organisation of event record
  common/evrecord/nsim,npart,offset,hadrotype,sqrts,collider,hadro,
     &shorthepmc,channel,isochannel
  integer nsim,npart,offset,hadrotype
  double precision sqrts
  character*4 collider,channel
  character*2 isochannel
  logical hadro,shorthepmc
C--extra storage for scattering centres before interactions
      common/storescatcen/nscatcen,maxnscatcen,scatflav(10000),
     &scatcen(10000,5),writescatcen,writedummies
  integer nscatcen,maxnscatcen,scatflav
  double precision scatcen
  logical writescatcen,writedummies

C--Variables local to this program
  INTEGER nold,pid,ipart,lme1,lme2,j,i,lme1orig,lme2orig,llep1,
     &llep2,lv
  DOUBLE PRECISION pyr,eni,qmax1,r,getmass,pyp,q1,q2,p21,p22,etot,
     &qmax2,pold,en1,en2,beta(3),enew1,enew2,emax,lambda,x1,x2,x3,
     &meweight,psweight,weight,eps1,eps2,theta1,theta2,z1,z2,
     &getltimemax,pi,m1,m2
  character*2 b1,b2
  CHARACTER*2 type1,type2
  LOGICAL firsttrip,which1,which2,isdiquark
  DATA pi/3.141592653589793d0/

   n=0
   colmax=600
   discard=.false.
       DO 91 i=1,23000
        mv(i,1)=0.d0
        mv(i,2)=0.d0
        mv(i,3)=0.d0
        mv(i,4)=0.d0
        mv(i,5)=0.d0
 91    CONTINUE
   nscatcen = 0

       CALL mednextevt

C--initialisation with matrix element  
C--production vertex
        CALL pickvtx(x0,y0)
        ltime=getltimemax()

 99   CALL pyevnt
        npart=n-offset
        evweight=pari(10)
    sumofweights=sumofweights+evweight
    IF((collider.EQ.'EEJJ').AND.(abs(k(8,2)).GT.6))THEN
     wdisc=wdisc+evweight
     ndisc=ndisc+1
     goto 102
    ELSE
     ngood=ngood+1
    ENDIF 

C--DY: don't have to do anything
    if (collider.eq.'PPDY') then
      CALL pyexec
      call converttohepmc(hpmcfid,ngood,pid,b1,b2)
      goto 102
    endif


C--   prepare event record
    if((collider.EQ.'PPZJ').OR.(collider.EQ.'PPZQ').or.
     &  (collider.EQ.'PPZG').or.(collider.EQ.'PPWJ').or.
     &  (collider.EQ.'PPWQ').or.(collider.EQ.'PPWG'))THEN 
             lme1orig=7
             lme2orig=8
         if(abs(k(7,2)).gt.21) then
           lv=7
     else
           lv=8
         endif
          ELSE
             lme1orig=offset-1
             lme2orig=offset
          ENDIF
        DO 180 ipart=offset+1, offset+npart
C--find decay leptons in V+jet events
    if((collider.EQ.'PPZJ').OR.(collider.EQ.'PPZQ').or.
     &  (collider.EQ.'PPZG').or.(collider.EQ.'PPWJ').or.
     &  (collider.EQ.'PPWQ').or.(collider.EQ.'PPWG'))THEN 
       if(k(ipart,3).eq.offset-1) llep1=ipart
       if(k(ipart,3).eq.offset) llep2=ipart
     endif
         IF(k(ipart,3).EQ.(lme1orig))THEN
          lme1=ipart
      IF(k(ipart,2).EQ.21)THEN
       type1='GC'
      ELSE
       type1='QQ'
      ENDIF
         ELSEIF(k(ipart,3).EQ.lme2orig)THEN
          lme2=ipart        
      IF(k(ipart,2).EQ.21)THEN
       type2='GC'
      ELSE
       type2='QQ'
      ENDIF
     ELSE
      trip(ipart)=0
      anti(ipart)=0
      zd(ipart)=0.d0
      thetaa(ipart)=0.d0
     ENDIF 
C--assign colour indices
         IF(k(ipart,1).EQ.2)THEN
      IF(k(ipart-1,1).EQ.2)THEN
C--in middle of colour singlet
       IF(firsttrip)THEN
        trip(ipart)=colmax+1
        anti(ipart)=trip(ipart-1)
       ELSE
        trip(ipart)=anti(ipart-1)
        anti(ipart)=colmax+1
       ENDIF
       colmax=colmax+1
      ELSE
C--beginning of colour singlet
       IF(((abs(k(ipart,2)).LT.10).AND.(k(ipart,2).GT.0))
     &      .OR.(isdiquark(k(ipart,2)).AND.(k(ipart,2).LT.0)))THEN
        trip(ipart)=colmax+1
        anti(ipart)=0
        firsttrip=.true.
       ELSE
        trip(ipart)=0
        anti(ipart)=colmax+1
        firsttrip=.false.
       ENDIF
       colmax=colmax+1
      ENDIF
     ENDIF 
         IF(k(ipart,1).EQ.1)THEN
C--end of colour singlet
      IF(firsttrip)THEN
       trip(ipart)=0
       anti(ipart)=trip(ipart-1)
      ELSE
       trip(ipart)=anti(ipart-1)
       anti(ipart)=0
      ENDIF
     ENDIF
 180    CONTINUE
    if (k(lme1,1).lt.11) k(lme1,1)=1
    if (k(lme2,1).lt.11) k(lme2,1)=1
    pid=k(lme1,2)
    eni=max(p(lme1,4),p(lme2,4))
    DO 183 ipart=offset+1, offset+npart
     IF((ipart.NE.lme1).AND.(ipart.NE.lme2).AND.(k(ipart,1).LT.11))
     &     k(ipart,1)=4
     if (k(ipart,2).eq.22) k(ipart,1)=4
 183    CONTINUE    

C--find virtualities and adapt four-vectors
    if((collider.EQ.'PPZJ').OR.(collider.EQ.'PPZQ').or.
     &  (collider.EQ.'PPZG').or.(collider.EQ.'PPWJ').or.
     &  (collider.EQ.'PPWQ').or.(collider.EQ.'PPWG'))THEN 
      if (abs(k(lme1,2)).gt.21) then
           qmax1=0.d0
           qmax2=sqrt(pari(18)+p(lme1,5)**2)
      else
           qmax1=sqrt(pari(18)+p(lme2,5)**2)
           qmax2=0.d0
      endif
           emax=p(lme1,4)+p(lme2,4)
           theta1=-1.d0
           theta2=-1.d0
        ELSEIF(collider.EQ.'PPJJ'.OR.collider.EQ.'PPYJ'
     &          .OR.collider.EQ.'PPYQ'.OR.collider.EQ.'PPYG')THEN
       if (k(lme1,1).eq.4) then
         qmax1 = 0.d0
       else
             qmax1=pari(17)
       endif
       if (k(lme2,1).eq.4) then
         qmax2 = 0.d0
       else
             qmax2=pari(17)
       endif
!        QMAX1=PYP(LME1,10)*exp(0.3*abs(pyp(lme1,17)-pyp(lme2,17))/2.)/2.
!        QMAX2=PYP(LME2,10)*exp(0.3*abs(pyp(lme1,17)-pyp(lme2,17))/2.)/2.
         emax=p(lme1,4)+p(lme2,4)
         theta1=-1.d0
         theta2=-1.d0
        ENDIF 
        en1=p(lme1,4)
        en2=p(lme2,4)
        beta(1)=(p(lme1,1)+p(lme2,1))/(p(lme1,4)+p(lme2,4))
        beta(2)=(p(lme1,2)+p(lme2,2))/(p(lme1,4)+p(lme2,4))
        beta(3)=(p(lme1,3)+p(lme2,3))/(p(lme1,4)+p(lme2,4))
        CALL pyrobo(lme1,lme1,0d0,0d0,-beta(1),-beta(2),-beta(3))
        CALL pyrobo(lme2,lme2,0d0,0d0,-beta(1),-beta(2),-beta(3))
    etot=p(lme1,4)+p(lme2,4)
    IF(collider.EQ.'EEJJ')THEN
         qmax1=etot
         qmax2=etot
     emax=p(lme1,4)+p(lme2,4)
     theta1=-1.d0
     theta2=-1.d0
        ENDIF
C--   find virtuality
        q1=getmass(0.d0,qmax1,theta1,emax,type1,emax,.false.,
     &       z1,which1)
        q2=getmass(0.d0,qmax2,theta2,emax,type2,emax,.false.,
     &       z2,which2)
 182    if (abs(k(lme1,2)).gt.21) then
      m1=p(lme1,5)
    else
      m1=q1
    endif
    if (abs(k(lme2,2)).gt.21) then
      m2=p(lme2,5)
    else
      m2=q2
    endif
        enew1=etot/2.d0 + (m1**2-m2**2)/(2.*etot)
        enew2=etot/2.d0 - (m1**2-m2**2)/(2.*etot)
    p21 = (etot/2.d0 + (m1**2-m2**2)/(2.*etot))**2 - m1**2
    p22 = (etot/2.d0 - (m1**2-m2**2)/(2.*etot))**2 - m2**2
    weight=1.d0
    IF((pyr(0).GT.weight).OR.(p21.LT.0.d0).OR.(p22.LT.0.d0)
     &  .OR.(enew1.LT.0.d0).OR.(enew2.LT.0.d0)
     &  )THEN
     IF(q1.GT.q2)THEN
          q1=getmass(0.d0,q1,theta1,emax,type1,emax,.false.,
     &  z1,which1)
     ELSE
          q2=getmass(0.d0,q2,theta2,emax,type2,emax,.false.,
     &  z2,which2)
     ENDIF
     goto 182
    ENDIF
        pold=pyp(lme1,8)
    p(lme1,1)=p(lme1,1)*sqrt(p21)/pold
    p(lme1,2)=p(lme1,2)*sqrt(p21)/pold
    p(lme1,3)=p(lme1,3)*sqrt(p21)/pold
    p(lme1,4)=enew1
    p(lme1,5)=m1
        pold=pyp(lme2,8)
    p(lme2,1)=p(lme2,1)*sqrt(p22)/pold
    p(lme2,2)=p(lme2,2)*sqrt(p22)/pold
    p(lme2,3)=p(lme2,3)*sqrt(p22)/pold
    p(lme2,4)=enew2
    p(lme2,5)=m2
        CALL pyrobo(lme1,lme1,0d0,0d0,beta(1),beta(2),beta(3))
        CALL pyrobo(lme2,lme2,0d0,0d0,beta(1),beta(2),beta(3))
C--correct for overestimated energy
    IF(q1.GT.0.d0)THEN
     eps1=0.5-0.5*sqrt(1.-q0**2/q1**2)
     &     *sqrt(1.-q1**2/p(lme1,4)**2)
     IF((z1.LT.eps1).OR.(z1.GT.(1.-eps1)))THEN
          q1=getmass(0.d0,q1,theta1,emax,type1,emax,.false.,
     &  z1,which1)
          CALL pyrobo(lme1,lme1,0d0,0d0,-beta(1),-beta(2),-beta(3))
          CALL pyrobo(lme2,lme2,0d0,0d0,-beta(1),-beta(2),-beta(3))
      goto 182
     ENDIF
    ENDIF 
    IF(q2.GT.0.d0)THEN
     eps2=0.5-0.5*sqrt(1.-q0**2/q2**2)
     &     *sqrt(1.-q2**2/p(lme2,4)**2)
         IF((z2.LT.eps2).OR.(z2.GT.(1.-eps2)))THEN
          q2=getmass(0.d0,q2,theta2,emax,type2,emax,.false.,
     &  z2,which2)
          CALL pyrobo(lme1,lme1,0d0,0d0,-beta(1),-beta(2),-beta(3))
          CALL pyrobo(lme2,lme2,0d0,0d0,-beta(1),-beta(2),-beta(3))
      goto 182
         ENDIF
        ENDIF
        
C--correct to ME for first parton
    IF(collider.EQ.'EEJJ')THEN
         beta(1)=(p(lme1,1)+p(lme2,1))/(p(lme1,4)+p(lme2,4))
         beta(2)=(p(lme1,2)+p(lme2,2))/(p(lme1,4)+p(lme2,4))
         beta(3)=(p(lme1,3)+p(lme2,3))/(p(lme1,4)+p(lme2,4))
         CALL pyrobo(lme1,lme1,0d0,0d0,-beta(1),-beta(2),-beta(3))
         CALL pyrobo(lme2,lme2,0d0,0d0,-beta(1),-beta(2),-beta(3))
         IF(q1.GT.0.d0)THEN
C--generate z value      
      x1=z1*(etot**2+q1**2)/etot**2
      x2=(etot**2-q1**2)/etot**2
      x3=(1.-z1)*(etot**2+q1**2)/etot**2
      psweight=(1.-x1)*(1.+(x1/(2.-x2))**2)/x3
     &  + (1.-x2)*(1.+(x2/(2.-x1))**2)/x3 
      meweight=x1**2+x2**2
      weight=meweight/psweight
      IF(pyr(0).GT.weight)THEN
 184       q1=getmass(0.d0,q1,theta1,emax,type1,emax,.false.,
     &  z1,which1)
      ENDIF
     ENDIF 
C--correct to ME for second parton
     IF(q2.GT.0.d0)THEN
C--generate z value      
      x1=(etot**2-q2**2)/etot**2
      x2=z2*(etot**2+q2**2)/etot**2
      x3=(1.-z2)*(etot**2+q2**2)/etot**2
      psweight=(1.-x1)*(1.+(x1/(2.-x2))**2)/x3
     &  + (1.-x2)*(1.+(x2/(2.-x1))**2)/x3 
      meweight=x1**2+x2**2
      weight=meweight/psweight
      IF(pyr(0).GT.weight)THEN
 185       q2=getmass(0.d0,q2,theta2,emax,type2,emax,.false.,
     &  z2,which2)
      ENDIF
     ENDIF
 186     enew1=etot/2.d0 + (q1**2-q2**2)/(2.*etot)
         enew2=etot/2.d0 - (q1**2-q2**2)/(2.*etot)
     p21 = (etot/2.d0 + (q1**2-q2**2)/(2.*etot))**2 - q1**2
     p22 = (etot/2.d0 - (q1**2-q2**2)/(2.*etot))**2 - q2**2
         pold=pyp(lme1,8)
     p(lme1,1)=p(lme1,1)*sqrt(p21)/pold
     p(lme1,2)=p(lme1,2)*sqrt(p21)/pold
     p(lme1,3)=p(lme1,3)*sqrt(p21)/pold
     p(lme1,4)=enew1
     p(lme1,5)=q1
         pold=pyp(lme2,8)
     p(lme2,1)=p(lme2,1)*sqrt(p22)/pold
     p(lme2,2)=p(lme2,2)*sqrt(p22)/pold
     p(lme2,3)=p(lme2,3)*sqrt(p22)/pold
     p(lme2,4)=enew2
     p(lme2,5)=q2
         CALL pyrobo(lme1,lme1,0d0,0d0,beta(1),beta(2),beta(3))
         CALL pyrobo(lme2,lme2,0d0,0d0,beta(1),beta(2),beta(3))
C--correct for overestimated energy
     IF(q1.GT.0.d0)THEN
     eps1=0.5-0.5*sqrt(1.-q0**2/q1**2)
     &     *sqrt(1.-q1**2/p(lme1,4)**2)
      IF((z1.LT.eps1).OR.(z1.GT.(1.-eps1)))THEN
           q1=getmass(0.d0,q1,theta1,emax,type1,emax,.false.,
     &  z1,which1)
           CALL pyrobo(lme1,lme1,0d0,0d0,-beta(1),-beta(2),-beta(3))
           CALL pyrobo(lme2,lme2,0d0,0d0,-beta(1),-beta(2),-beta(3))
       goto 186
      ENDIF
     ENDIF 
     IF(q2.GT.0.d0)THEN
     eps2=0.5-0.5*sqrt(1.-q0**2/q2**2)
     &     *sqrt(1.-q2**2/p(lme2,4)**2)
          IF((z2.LT.eps2).OR.(z2.GT.(1.-eps2)))THEN
           q2=getmass(0.d0,q2,theta2,emax,type2,emax,.false.,
     &  z2,which2)
           CALL pyrobo(lme1,lme1,0d0,0d0,-beta(1),-beta(2),-beta(3))
           CALL pyrobo(lme2,lme2,0d0,0d0,-beta(1),-beta(2),-beta(3))
       goto 186
          ENDIF
         ENDIF 
    ENDIF

C--transfer recoil to decay leptons in V+jet
    if((collider.EQ.'PPZJ').OR.(collider.EQ.'PPZQ').or.
     &  (collider.EQ.'PPZG').or.(collider.EQ.'PPWJ').or.
     &  (collider.EQ.'PPWQ').or.(collider.EQ.'PPWG'))THEN 
      beta(1)=p(lv,1)/p(lv,4)
      beta(2)=p(lv,2)/p(lv,4)
      beta(3)=p(lv,3)/p(lv,4)
          CALL pyrobo(llep1,llep1,0d0,0d0,-beta(1),-beta(2),-beta(3))
          CALL pyrobo(llep2,llep2,0d0,0d0,-beta(1),-beta(2),-beta(3))
      if (abs(k(lme1,2)).gt.21) then
        beta(1)=p(lme1,1)/p(lme1,4)
        beta(2)=p(lme1,2)/p(lme1,4)
        beta(3)=p(lme1,3)/p(lme1,4)
      else
        beta(1)=p(lme2,1)/p(lme2,4)
        beta(2)=p(lme2,2)/p(lme2,4)
        beta(3)=p(lme2,3)/p(lme2,4)
      endif
          CALL pyrobo(llep1,llep1,0d0,0d0,beta(1),beta(2),beta(3))
          CALL pyrobo(llep2,llep2,0d0,0d0,beta(1),beta(2),beta(3))
    endif

  
        za(lme1)=1.d0
        za(lme2)=1.d0
    thetaa(lme1)=p(lme1,5)/(sqrt(z1*(1.-z1))*p(lme1,4))
    thetaa(lme2)=p(lme2,5)/(sqrt(z2*(1.-z2))*p(lme2,4))
    zd(lme1)=z1
    zd(lme2)=z2
    qqbard(lme1)=which1
    qqbard(lme2)=which2

        mv(lme1,1)=x0
        mv(lme1,2)=y0
        mv(lme1,3)=0.d0
        mv(lme1,4)=0.d0
        IF(p(lme1,5).GT.0.d0)THEN
         lambda=1.d0/(ftfac*p(lme1,4)*0.2/q1**2)
          mv(lme1,5)=-log(1.d0-pyr(0))/lambda
        ELSE
         mv(lme1,5)=ltime
        ENDIF
         
        mv(lme2,1)=x0
        mv(lme2,2)=y0
        mv(lme2,3)=0.d0
        mv(lme2,4)=0.d0
        IF(p(lme2,5).GT.0.d0)THEN
         lambda=1.d0/(ftfac*p(lme2,4)*0.2/q2**2)
          mv(lme2,5)=-log(1.d0-pyr(0))/lambda
        ELSE
         mv(lme2,5)=ltime
        ENDIF

C--develop parton shower
   CALL makecascade
   IF(discard) THEN
    ngood=ngood-1
    wdisc=wdisc+evweight
    ndisc=ndisc+1
        write(logfid,*)'discard event',j
    goto 102
   ENDIF

       IF(.NOT.allhad)THEN
        DO 86 i=1,n
         IF(k(i,1).EQ.3) k(i,1)=22
 86     CONTINUE
       ENDIF
       IF(hadro)THEN
        CALL makestrings(hadrotype)
    IF(discard) THEN
         write(logfid,*)'discard event',j
     wdisc=wdisc+evweight
     ndisc=ndisc+1
     ngood=ngood-1
     goto 102
    ENDIF
        CALL pyexec
    IF(mstu(30).NE.errcount)THEN
         write(logfid,*)'PYTHIA discards event',j,
     &  '  (error number',mstu(30),')'
     errcount=mstu(30)
     wdisc=wdisc+evweight
     ndisc=ndisc+1
     ngood=ngood-1
     goto 102
    ENDIF
       ENDIF

   IF(mstu(30).NE.errcount)THEN
    errcount=mstu(30)
   ELSE 
    CALL converttohepmc(hpmcfid,ngood,pid,b1,b2)
   ENDIF

C--write message to log-file
 102  IF(nsim.GT.100)THEN
       IF(mod(j,nsim/100).EQ.0)THEN
    write(logfid,*) 'done with event number ',j
   ENDIF
  else
    write(logfid,*) 'done with event number ',j
      ENDIF
  call flush(logfid)
  end



***********************************************************************
***   subroutine makestrings
***********************************************************************
  SUBROUTINE makestrings(WHICH)
  IMPLICIT NONE
C--identifier of file for hepmc output and logfile
  common/hepmcid/hpmcfid,logfid
  integer hpmcfid,logfid
  INTEGER which
  IF(which.EQ.0)THEN
   CALL makestrings_vac
  ELSEIF(which.EQ.1)THEN
   CALL makestrings_minl
  ELSE
  WRITE(logfid,*)'error: unknown hadronisation type in MAKESTRINGS'
  ENDIF
  END


***********************************************************************
***   subroutine makestrings_vac
***********************************************************************
      SUBROUTINE makestrings_vac
      IMPLICIT NONE
C--identifier of file for hepmc output and logfile
  common/hepmcid/hpmcfid,logfid
  integer hpmcfid,logfid
C--Common block of Pythia
      common/pyjets/n,npad,k(23000,5),p(23000,5),v(23000,5)
  INTEGER n,npad,k
  DOUBLE PRECISION p,v
C--Parameter common block
  common/param/q0,lps,lqcd,ltime,scalefacm,angord,scatrecoil,
     &allhad,compress,nf
      INTEGER nf
  DOUBLE PRECISION q0,lqcd,ltime,lps,scalefacm
      LOGICAL angord,scatrecoil,allhad,compress
C--colour index common block
  common/colour/trip(23000),anti(23000),colmax
  INTEGER trip,anti,colmax
C--discard event flag
  common/disc/ndisc,nstrange,ngood,errcount,wdisc,discard
  LOGICAL discard
  INTEGER ndisc,nstrange,ngood,errcount
  double precision wdisc
C--local variables
      INTEGER nold,i,j,lquark,lmatch,lloose,nold1
      DOUBLE PRECISION eaddend,pyr,dir
      LOGICAL isdiquark,compressevent,roomleft
      DATA eaddend/10.d0/
  
  i = 0
  if (compress) roomleft = compressevent(i)
      nold1=n
C--remove all active lines that are leptons, gammas, hadrons etc.
  DO 52 i=1,nold1
   IF((k(i,1).EQ.4).AND.(trip(i).EQ.0).AND.(anti(i).EQ.0))THEN
C--copy line to end of event record
        n=n+1
        IF(n.GT.22990) THEN
         write(logfid,*)'event too long for event record'
         discard=.true.
         RETURN
        ENDIF
        k(n,1)=11
        k(n,2)=k(i,2)
        k(n,3)=i
        k(n,4)=0
        k(n,5)=0
        p(n,1)=p(i,1)
        p(n,2)=p(i,2)
        p(n,3)=p(i,3)
        p(n,4)=p(i,4)
        p(n,5)=p(i,5)
        k(i,1)=17
        k(i,4)=n
        k(i,5)=n
    trip(n)=trip(i)
    anti(n)=anti(i)
   ENDIF
 52 CONTINUE
      nold=n
C--first do strings with existing (anti)triplets
C--find string end (=quark or antiquark)
 43   lquark=0
      DO 40 i=1,nold
       IF((k(i,1).EQ.11).OR.(k(i,1).EQ.12).OR.(k(i,1).EQ.13)
     &            .OR.(k(i,1).EQ.14)) k(i,1)=17
       IF(((k(i,1).EQ.1).OR.(k(i,1).EQ.3).OR.(k(i,1).EQ.4).OR.
     &   (k(i,1).EQ.5)).AND.((k(i,2).LT.6).OR.isdiquark(k(i,2))))THEN
        lquark=i
    goto 41
       ENDIF
 40   CONTINUE
  goto 50
 41 CONTINUE
C--copy string end to end of event record
      n=n+1
      IF(n.GT.22990) THEN
       write(logfid,*)'event too long for event record'
       discard=.true.
       RETURN
      ENDIF
      k(n,1)=2
      k(n,2)=k(lquark,2)
      k(n,3)=lquark
      k(n,4)=0
      k(n,5)=0
      p(n,1)=p(lquark,1)
      p(n,2)=p(lquark,2)
      p(n,3)=p(lquark,3)
      p(n,4)=p(lquark,4)
      p(n,5)=p(lquark,5)
      k(lquark,1)=16
      k(lquark,4)=n
      k(lquark,5)=n
  trip(n)=trip(lquark)
  anti(n)=anti(lquark)
C--append matching colour partner
  lmatch=0
  DO 44 j=1,10000000
   DO 42 i=1,nold
    IF(((k(i,1).EQ.1).OR.(k(i,1).EQ.3).OR.(k(i,1).EQ.4)
     &            .OR.(k(i,1).EQ.5))
     &      .AND.(((trip(i).EQ.anti(n)).AND.(trip(i).NE.0))
     &    .OR.((anti(i).EQ.trip(n)).AND.(anti(i).NE.0))))THEN
         n=n+1
         IF(n.GT.22990) THEN
          write(logfid,*)'event too long for event record'
          discard=.true.
          RETURN
         ENDIF
         k(n,2)=k(i,2)
         k(n,3)=i
         k(n,4)=0
         k(n,5)=0
         p(n,1)=p(i,1)
         p(n,2)=p(i,2)
         p(n,3)=p(i,3)
         p(n,4)=p(i,4)
         p(n,5)=p(i,5)
     trip(n)=trip(i)
     anti(n)=anti(i)
         k(i,1)=16
         k(i,4)=n
         k(i,5)=n
         IF(k(i,2).EQ.21)THEN
          k(n,1)=2
          goto 44
         ELSE
          k(n,1)=1
          goto 43
         ENDIF
    ENDIF
 42  CONTINUE
C--no matching colour partner found
   write(logfid,*)'Error in MAKESTRINGS_VAC: failed to reconstruct '//
     &'colour singlet system, will discard event'
   discard = .true.
   return
 44 CONTINUE
C--now take care of purely gluonic remainder system
C-----------------------------------------
C--find gluon where anti-triplet is not matched
 50   lloose=0
      DO 45 i=1,nold
       IF(((k(i,1).EQ.1).OR.(k(i,1).EQ.3).OR.(k(i,1).EQ.4)
     &          .OR.(k(i,1).EQ.5)))THEN
    DO 46 j=1,nold
     IF(((k(i,1).EQ.1).OR.(k(i,1).EQ.3).OR.(k(i,1).EQ.4)
     &          .OR.(k(i,1).EQ.5)))THEN
      IF(anti(i).EQ.trip(j)) goto 45
     ENDIF
 46   CONTINUE
        lloose=i
    goto 47
       ENDIF
 45   CONTINUE
  goto 51
 47 CONTINUE
C--generate artificial triplet end
   write(logfid,*)'Error in MAKESTRINGS_VAC: failed to reconstruct '//
     &'colour singlet system, will discard event'
   discard = .true.
   return
C--copy loose gluon to end of event record
      n=n+1
      IF(n.GT.22990) THEN
       write(logfid,*)'event too long for event record'
       discard=.true.
       RETURN
      ENDIF
      k(n,1)=2
      k(n,2)=k(lloose,2)
      k(n,3)=lloose
      k(n,4)=0
      k(n,5)=0
      p(n,1)=p(lloose,1)
      p(n,2)=p(lloose,2)
      p(n,3)=p(lloose,3)
      p(n,4)=p(lloose,4)
      p(n,5)=p(lloose,5)
      k(lloose,1)=16
      k(lloose,4)=n
      k(lloose,5)=n
  trip(n)=trip(lloose)
  anti(n)=anti(lloose)
C--append matching colour partner
  lmatch=0
  DO 48 j=1,10000000
   DO 49 i=1,nold
    IF(((k(i,1).EQ.1).OR.(k(i,1).EQ.3).OR.(k(i,1).EQ.4)
     &        .OR.(k(i,1).EQ.5))
     &    .AND.(anti(i).EQ.trip(n)))THEN
         n=n+1
         IF(n.GT.22990) THEN
          write(logfid,*)'event too long for event record'
          discard=.true.
          RETURN
         ENDIF
         k(n,2)=k(i,2)
         k(n,3)=i
         k(n,4)=0
         k(n,5)=0
         p(n,1)=p(i,1)
         p(n,2)=p(i,2)
         p(n,3)=p(i,3)
         p(n,4)=p(i,4)
         p(n,5)=p(i,5)
     trip(n)=trip(i)
     anti(n)=anti(i)
         k(i,1)=16
         k(i,4)=n
         k(i,5)=n
         k(n,1)=2
         goto 48
    ENDIF
 49  CONTINUE
C--no matching colour partner found, add artificial end point
   write(logfid,*)'Error in MAKESTRINGS_VAC: failed to reconstruct '//
     &'colour singlet system, will discard event'
   discard = .true.
   return
 48 CONTINUE
 51 CONTINUE
  CALL cleanup(nold1)
  END


***********************************************************************
***   subroutine makestrings_minl
***********************************************************************
      SUBROUTINE makestrings_minl
      IMPLICIT NONE
C--Common block of Pythia
      common/pyjets/n,npad,k(23000,5),p(23000,5),v(23000,5)
  INTEGER n,npad,k
  DOUBLE PRECISION p,v
C--Parameter common block
  common/param/q0,lps,lqcd,ltime,scalefacm,angord,scatrecoil,
     &allhad,compress,nf
      INTEGER nf
  DOUBLE PRECISION q0,lqcd,ltime,lps,scalefacm
      LOGICAL angord,scatrecoil,allhad,compress
C--colour index common block
  common/colour/trip(23000),anti(23000),colmax
  INTEGER trip,anti,colmax
C--local variables
      INTEGER nold,i,j,lmax,lmin,lend,nold1
      DOUBLE PRECISION emax,minv,mmin,z,generatez,mcut,eaddend,pyr,dir,
     &pyp
      DATA mcut/1.d8/
      DATA eaddend/10.d0/
C--identifier of file for hepmc output and logfile
  common/hepmcid/hpmcfid,logfid
  integer hpmcfid,logfid
C--discard event flag
  common/disc/ndisc,nstrange,ngood,errcount,wdisc,discard
  LOGICAL discard
  INTEGER ndisc,nstrange,ngood,errcount
  double precision wdisc
  logical compressevent,roomleft

   i = 0
   if (compress) roomleft = compressevent(i)
      nold1=n
C--remove all active lines that are leptons, gammas, hadrons etc.
  DO 52 i=1,nold1
   IF((k(i,1).EQ.4).AND.(trip(i).EQ.0).AND.(anti(i).EQ.0))THEN
C--copy line to end of event record
        n=n+1
        IF(n.GT.22990) THEN
         write(logfid,*)'event too long for event record'
         discard=.true.
         RETURN
        ENDIF
        k(n,1)=11
        k(n,2)=k(i,2)
        k(n,3)=i
        k(n,4)=0
        k(n,5)=0
        p(n,1)=p(i,1)
        p(n,2)=p(i,2)
        p(n,3)=p(i,3)
        p(n,4)=p(i,4)
        p(n,5)=p(i,5)
        k(i,1)=17
        k(i,4)=n
        k(i,5)=n
    trip(n)=trip(i)
    anti(n)=anti(i)
   ENDIF
 52 CONTINUE
       nold=n
C--find most energetic unfragmented parton in event
 43    emax=0
       lmax=0
       DO 40 i=1,nold
        IF((k(i,1).EQ.11).OR.(k(i,1).EQ.12).OR.(k(i,1).EQ.13)
     &            .OR.(k(i,1).EQ.14)) k(i,1)=17
        if (abs(pyp(i,17)).gt.4.d0) k(i,1)=17
        IF(((k(i,1).EQ.1).OR.(k(i,1).EQ.3).OR.(k(i,1).EQ.4)
     &  .OR.(k(i,1).EQ.5)).AND.(p(i,4).GT.emax))THEN
         emax=p(i,4)
         lmax=i
        ENDIF
 40    CONTINUE
C--if there is non, we are done
       IF(lmax.EQ.0) goto 50
C--check if highest energy parton is (anti)quark or gluon
       IF(k(lmax,2).EQ.21)THEN
C--split gluon in qqbar pair and store one temporarily in line 1
C--make new line in event record for string end
        n=n+2
        IF(n.GT.22990) THEN
         write(logfid,*)'event too long for event record'
         discard=.true.
         RETURN
        ENDIF
    IF((n-2).GT.nold)THEN
         DO 47 j=nold,n-3
          k(n+nold-j,1)=k(n+nold-j-2,1)
          k(n+nold-j,2)=k(n+nold-j-2,2)
          IF(k(n+nold-j-2,3).GT.nold) THEN
           k(n+nold-j,3)=k(n+nold-j-2,3)+2
          ELSE
           k(n+nold-j,3)=k(n+nold-j-2,3)
          ENDIF
          k(n+nold-j,4)=0
          k(n+nold-j,5)=0
          p(n+nold-j,1)=p(n+nold-j-2,1)
          p(n+nold-j,2)=p(n+nold-j-2,2)
          p(n+nold-j,3)=p(n+nold-j-2,3)
          p(n+nold-j,4)=p(n+nold-j-2,4)
          p(n+nold-j,5)=p(n+nold-j-2,5)
          k(k(n+nold-j-2,3),4)=k(k(n+nold-j-2,3),4)+2
          k(k(n+nold-j-2,3),5)=k(k(n+nold-j-2,3),5)+2
 47      CONTINUE
    ENDIF
        nold=nold+2
        k(lmax,1)=18
        z=generatez(0.d0,0.d0,1.d-3,'QG')
        IF(z.GT.0.5)THEN
         k(nold-1,2)=1
         k(nold,2)=-1
        ELSE
         z=1.-z
         k(nold-1,2)=-1
         k(nold,2)=1
        ENDIF
        k(nold-1,1)=1
        k(nold-1,3)=lmax
        k(nold-1,4)=0
        k(nold-1,5)=0
        p(nold-1,1)=(1.-z)*p(lmax,1)
        p(nold-1,2)=(1.-z)*p(lmax,2)
        p(nold-1,3)=(1.-z)*p(lmax,3)
        p(nold-1,4)=(1.-z)*p(lmax,4)
        p(nold-1,5)=p(lmax,5)
        k(nold,1)=1
        k(nold,3)=lmax
        k(nold,4)=0
        k(nold,5)=0
        p(nold,1)=z*p(lmax,1)
        p(nold,2)=z*p(lmax,2)
        p(nold,3)=z*p(lmax,3)
        p(nold,4)=z*p(lmax,4)
        p(nold,5)=p(lmax,5)
        k(lmax,1)=18
        k(lmax,4)=nold-1
        k(lmax,5)=nold
        lmax=nold
       ENDIF
       n=n+1
       IF(n.GT.22990) THEN
        write(logfid,*)'event too long for event record'
        discard=.true.
        RETURN
       ENDIF
       k(n,1)=2
       k(n,2)=k(lmax,2)
       k(n,3)=lmax
       k(n,4)=0
       k(n,5)=0
       p(n,1)=p(lmax,1)
       p(n,2)=p(lmax,2)
       p(n,3)=p(lmax,3)
       p(n,4)=p(lmax,4)
       p(n,5)=p(lmax,5)
       k(lmax,1)=16
       k(lmax,4)=n
       k(lmax,5)=n
       lend=lmax
C--find closest partner
 42    mmin=1.d10
       lmin=0
       DO 41 i=1,nold
        IF(((k(i,1).EQ.1).OR.(k(i,1).EQ.3).OR.(k(i,1)
     &      .EQ.4).OR.(k(i,1).EQ.5))
     &      .AND.((k(i,2).EQ.21).OR.((k(i,2)*k(lend,2).LT.0.d0).AND.
     &    (k(i,3).NE.k(lend,3))))
     &      .AND.(p(i,1)*p(lend,1).GT.0.d0))THEN
         minv=p(i,4)*p(lmax,4)-p(i,1)*p(lmax,1)-p(i,2)*p(lmax,2)
     &            -p(i,3)*p(lmax,3)
         IF((minv.LT.mmin).AND.(minv.GT.0.d0).AND.(minv.LT.mcut))THEN
          mmin=minv
          lmin=i
         ENDIF
        ENDIF
 41    CONTINUE
C--if no closest partner can be found, generate artificial end point for string
       IF(lmin.EQ.0)THEN
        n=n+1
        IF(n.GT.22990) THEN
         write(logfid,*)'event too long for event record'
         discard=.true.
         RETURN
        ENDIF
        k(n,1)=1
        k(n,2)=-k(lend,2)
        k(n,3)=0
        k(n,4)=0
        k(n,5)=0
        p(n,1)=0.d0
        p(n,2)=0.d0
        IF(pyr(0).LT.0.5)THEN
         dir=1.d0
        ELSE
         dir=-1.d0
        ENDIF
        p(n,3)=dir*eaddend
        p(n,4)=eaddend
        p(n,5)=0.d0
        goto 43
       ELSE
C--else build closest partner in string
        n=n+1
        IF(n.GT.22990) THEN
         write(logfid,*)'event too long for event record'
         discard=.true.
         RETURN
        ENDIF
        k(n,2)=k(lmin,2)
        k(n,3)=lmin
        k(n,4)=0
        k(n,5)=0
        p(n,1)=p(lmin,1)
        p(n,2)=p(lmin,2)
        p(n,3)=p(lmin,3)
        p(n,4)=p(lmin,4)
        p(n,5)=p(lmin,5)
        k(lmin,1)=16
        k(lmin,4)=n
        k(lmin,5)=n
        IF(k(lmin,2).EQ.21)THEN
         k(n,1)=2
         lmax=lmin
         goto 42
        ELSE
         k(n,1)=1
         goto 43
        ENDIF
       ENDIF
 50    CONTINUE
       CALL cleanup(nold)
      END


***********************************************************************
***   subroutine cleanup
***********************************************************************
  SUBROUTINE cleanup(NFIRST)
  IMPLICIT NONE
C--Common block of Pythia
      common/pyjets/n,npad,k(23000,5),p(23000,5),v(23000,5)
  INTEGER n,npad,k
  DOUBLE PRECISION p,v
C--local variables
  INTEGER nfirst,nlast,i,j
  
  nlast=n
  DO 21 i=1,nlast-nfirst
   DO 22 j=1,5
    k(i,j)=k(nfirst+i,j)
    p(i,j)=p(nfirst+i,j)
    v(i,j)=v(nfirst+i,j)
 22  CONTINUE
   k(i,3)=0  
 21 CONTINUE
      n=nlast-nfirst
  END


***********************************************************************
***   subroutine makecascade
***********************************************************************
  SUBROUTINE makecascade
      IMPLICIT NONE
C--Common block of Pythia
      common/pyjets/n,npad,k(23000,5),p(23000,5),v(23000,5)
  INTEGER n,npad,k
  DOUBLE PRECISION p,v
C--time common block
      common/time/mv(23000,5)
      DOUBLE PRECISION mv
C--Parameter common block
  common/param/q0,lps,lqcd,ltime,scalefacm,angord,scatrecoil,
     &allhad,compress,nf
      INTEGER nf
  DOUBLE PRECISION q0,lqcd,ltime,lps,scalefacm
      LOGICAL angord,scatrecoil,allhad,compress
C--discard event flag
  common/disc/ndisc,nstrange,ngood,errcount,wdisc,discard
  LOGICAL discard
  INTEGER ndisc,nstrange,ngood,errcount
  double precision wdisc

C--local variables
  INTEGER nold,i
  LOGICAL cont

 10 nold=n
  cont=.false.
  DO 11 i=2,nold
   if (i.gt.n) goto 10
C--check if parton may evolve, i.e. do splitting or scattering
   IF((k(i,1).EQ.1).OR.(k(i,1).EQ.2))THEN
    cont=.true.
    CALL makebranch(i)
    IF(discard) goto 12
   ENDIF
 11 CONTINUE
  IF(cont) goto 10
 12 END


***********************************************************************
***   subroutine makebranch
***********************************************************************
      SUBROUTINE makebranch(L)
      IMPLICIT NONE
C--Common block of Pythia
      common/pyjets/n,npad,k(23000,5),p(23000,5),v(23000,5)
  INTEGER n,npad,k
  DOUBLE PRECISION p,v
C--time common block
      common/time/mv(23000,5)
      DOUBLE PRECISION mv
C--Parameter common block
  common/param/q0,lps,lqcd,ltime,scalefacm,angord,scatrecoil,
     &allhad,compress,nf
      INTEGER nf
  DOUBLE PRECISION q0,lqcd,ltime,lps,scalefacm
      LOGICAL angord,scatrecoil,allhad,compress
C--discard event flag
  common/disc/ndisc,nstrange,ngood,errcount,wdisc,discard
  LOGICAL discard
  INTEGER ndisc,nstrange,ngood,errcount
  double precision wdisc
C--variables for angular ordering
      common/angor/za(23000),zd(23000),thetaa(23000),qqbard(23000)
  DOUBLE PRECISION za,zd,thetaa
      LOGICAL qqbard
C--number of scattering events
  common/check/nscat,nscateff,nsplit
  DOUBLE PRECISION nscat,nscateff,nsplit
C--variables for coherent scattering
  common/coherent/nstart,nend,allqs(10000,6),scatcentres(10000,10),
     &qsumvec(4),qsum2
  INTEGER nstart,nend
  DOUBLE PRECISION allqs,scatcentres,qsumvec,qsum2
C--event weight
  common/weight/evweight,sumofweights
  double precision evweight,sumofweights
C--identifier of file for hepmc output and logfile
  common/hepmcid/hpmcfid,logfid
  integer hpmcfid,logfid
C--extra storage for scattering centres before interactions
       common/storescatcen/nscatcen,maxnscatcen,scatflav(10000),
     & scatcen(10000,5),writescatcen,writedummies
   integer nscatcen,maxnscatcen,scatflav
   double precision scatcen
   logical writescatcen,writedummies
C--local variables
      INTEGER l,line,nold,typi,lineold,lkine,nendold,nscatcenold
      DOUBLE PRECISION theta,phi,pyp,formtime,starttime,tleft,
     &tsum,deltat,newmass,getmass,q,getms,zdec,x,dtcorr
  LOGICAL overq0,qqbardec
  CHARACTER typ
  LOGICAL radiation,retrysplit,medind,roomleft,compressevent

  line=l
  nstart=0
  nend=0
  starttime=mv(line,4)
  tsum=0.d0
  qsum2=0.d0
  qsumvec(1)=0.d0
  qsumvec(2)=0.d0
  qsumvec(3)=0.d0
  qsumvec(4)=0.d0
  retrysplit=.false.
      medind=.false.
  x=0.d0
  q=0.d0
  typi=0

      IF ((n.GT.20000).and.compress) roomleft = compressevent(line)

20  IF(discard) RETURN
  IF(((k(line,1).EQ.1).AND.(p(line,5).GT.0.d0))
     &  .OR.((k(line,1).EQ.2).AND.(zd(line).gt.0.d0)))THEN
       IF(medind)THEN
        formtime=starttime
       ELSE 
    formtime=min(mv(line,5),ltime)
   ENDIF
   radiation=.true.
  ELSE
   formtime=ltime
   radiation=.false.
  ENDIF
  tleft=formtime-starttime
      IF(k(line,2).EQ.21)THEN
       typ='G'
      ELSE
       typ='Q'
      ENDIF
      medind=.false.

      IF(tleft.LE.1.d-10)THEN
C--no scattering
   IF(radiation)THEN
C--if there is radiation associated with the parton then form it now
C--rotate such that momentum points in z-direction
        nold=n
        nscatcenold=nscatcen
        theta=pyp(line,13)
        phi=pyp(line,15)
        CALL pyrobo(line,line,0d0,-phi,0d0,0d0,0d0)
        CALL pyrobo(line,line,-theta,0d0,0d0,0d0,0d0)
        CALL makesplitting(line)
C--rotate back
        CALL pyrobo(line,line,theta,0d0,0d0,0d0,0d0)
        CALL pyrobo(line,line,0d0,phi,0d0,0d0,0d0)
        IF(discard) RETURN
        CALL pyrobo(n-1,n,theta,0d0,0d0,0d0,0d0)
        CALL pyrobo(n-1,n,0d0,phi,0d0,0d0,0d0)
C--set the production vertices: x_mother + (tprod - tprod_mother) * beta_mother
        mv(n-1,1)=mv(line,1)
     &  +(mv(n-1,4)-mv(line,4))*p(line,1)/max(pyp(line,8),p(line,4))
        mv(n-1,2)=mv(line,2)
     &  +(mv(n-1,4)-mv(line,4))*p(line,2)/max(pyp(line,8),p(line,4))
        mv(n-1,3)=mv(line,3)
     &  +(mv(n-1,4)-mv(line,4))*p(line,3)/max(pyp(line,8),p(line,4))
        mv(n,  1)=mv(line,1)
     &  +(mv(n,  4)-mv(line,4))*p(line,1)/max(pyp(line,8),p(line,4))
        mv(n,  2)=mv(line,2)
     &  +(mv(n,  4)-mv(line,4))*p(line,2)/max(pyp(line,8),p(line,4))
        mv(n,  3)=mv(line,3)
     &  +(mv(n,  4)-mv(line,4))*p(line,3)/max(pyp(line,8),p(line,4))

    line=n
    nstart=0
    nend=0
    starttime=mv(n,4)
    qsumvec(1)=0.d0
    qsumvec(2)=0.d0
    qsumvec(3)=0.d0
    qsumvec(4)=0.d0
    qsum2=0.d0
    tsum=0.d0
    goto 21
   ELSE
    nstart=0
    nend=0
    starttime=formtime
    qsumvec(1)=0.d0
    qsumvec(2)=0.d0
    qsumvec(3)=0.d0
    qsumvec(4)=0.d0
    qsum2=0.d0
    tsum=0.d0
    goto 21
   ENDIF
  ELSE
C--do scattering
C--find delta t for the scattering
   deltat=tleft
   overq0=.false.
   CALL doinstatescat(line,x,typi,q,starttime+tsum,deltat,
     &    overq0,.false.)
   tsum=tsum+deltat
   tleft=tleft-deltat
C--do initial state splitting if there is one
   nold=n
   lineold=line
   zdec=zd(line)
   qqbardec=qqbard(line)
        nscatcenold=nscatcen
 25  IF(x.LT.1.d0) THEN
    CALL makeinsplit(line,x,qsum2,q,typi,starttime+tsum,deltat)
        IF(discard) RETURN
    IF(x.LT.1.d0)THEN
     line=n
     lkine=n
     IF(k(line,2).EQ.21)THEN
      newmass=getmass(0.d0,scalefacm*sqrt(-qsum2),-1.d0,p(line,4),
     &      'GC',sqrt(-qsum2),.false.,zdec,qqbardec)
          IF(zdec.GT.0.d0)THEN
           thetaa(line)=newmass/(sqrt(zdec*(1.-zdec))*p(line,4))
          ELSE
           thetaa(line)=0.d0
          ENDIF 
      zd(line)=zdec
      qqbard(line)=qqbardec
     ELSE 
      newmass=getmass(0.d0,scalefacm*sqrt(-qsum2),-1.d0,p(line,4),
     &      'QQ',sqrt(-qsum2),.false.,zdec,qqbardec)
      IF(zdec.GT.0.d0)THEN
           thetaa(line)=newmass/(sqrt(zdec*(1.-zdec))*p(line,4))
          ELSE
           thetaa(line)=0.d0
          ENDIF 
      zd(line)=zdec
      qqbard(line)=qqbardec
     ENDIF
     zdec=zd(line)
     qqbardec=qqbard(line)
    ELSE
     lkine=line
     nend=nstart
     qsum2=allqs(nend,1)
     qsumvec(1)=allqs(nend,2)
     qsumvec(2)=allqs(nend,3)
     qsumvec(3)=allqs(nend,4)
     qsumvec(4)=allqs(nend,5)
     IF(-allqs(nend,1).GT.q0**2/scalefacm**2)THEN
      overq0=.true.
     ELSE
      overq0=.false.
     ENDIF
     tleft = starttime+tsum+tleft-allqs(1,6)
     tsum = allqs(1,6)-starttime
    ENDIF 
   ENDIF
   IF(x.EQ.1.d0)THEN
    newmass=0.d0
    IF(nend.GT.0)THEN
     CALL dofistatescat(line,starttime+tsum,tleft,deltat,
     &    newmass,overq0,zdec,qqbardec)
     IF(newmass.GT.(p(line,5)*(1.d0+1.d-6)))THEN
      medind=.true.
     ELSE
      medind=.false.
      zdec=zd(line)
      qqbardec=qqbard(line)
     ENDIF 
     tsum=tsum+deltat
     tleft=tleft-deltat
     lkine=line
    ENDIF
   ENDIF
C--do kinematics
   retrysplit=.false.
   IF(nend.GT.0) THEN
    nendold=nend
    CALL dokinematics(lkine,lineold,nstart,nend,newmass,retrysplit,
     &    starttime+tsum,x,zdec,qqbardec)
    IF(retrysplit) THEN
     tleft = starttime+tsum+tleft-allqs(1,6)
     tsum = allqs(1,6)-starttime
     if (x.lt.1.d0) then
       nend=nstart
       qsum2=allqs(nend,1)
       qsumvec(1)=allqs(nend,2)
       qsumvec(2)=allqs(nend,3)
       qsumvec(3)=allqs(nend,4)
       qsumvec(4)=allqs(nend,5)
       typi=k(l,2)
       IF(-allqs(nend,1).GT.q0**2/scalefacm**2)THEN
         overq0=.true.
       ELSE
         overq0=.false.
       ENDIF
       n=nold
       line=lineold
       x=1.d0
       k(line,1)=1
       nscatcen=nscatcenold
       nsplit=nsplit-evweight
       goto 25
     else
       line=n
       starttime=starttime+tsum
       tsum=0.d0
     endif
    ELSE
     line=n
     starttime=starttime+tsum
     tsum=0.d0
    ENDIF
   ELSE
    starttime=starttime+tsum
    tsum=0.d0
   ENDIF
   IF(p(line,5).GT.0.d0) radiation=.true.
  ENDIF

 21   IF(((k(line,1).EQ.1).AND.(p(line,5).GT.0.d0))
     &  .OR.((k(line,1).EQ.2).AND.(zd(line).gt.0.d0))
     &  .OR.(starttime.LT.ltime))THEN
   goto 20
  ENDIF
  IF((k(line,1).EQ.1).AND.(p(line,5).EQ.0.d0)) k(line,1)=4
  IF((k(line,1).EQ.2).AND.(zd(line).lt.0.d0)) k(line,1)=5
      END



***********************************************************************
***   subroutine makesplitting
***********************************************************************
  SUBROUTINE makesplitting(L)
  IMPLICIT NONE
C--identifier of file for hepmc output and logfile
  common/hepmcid/hpmcfid,logfid
  integer hpmcfid,logfid
C--Common block of Pythia
      common/pyjets/n,npad,k(23000,5),p(23000,5),v(23000,5)
  INTEGER n,npad,k
  DOUBLE PRECISION p,v
C--time common block
      common/time/mv(23000,5)
      DOUBLE PRECISION mv
C--factor in front of formation times
  common/ftimefac/ftfac
  DOUBLE PRECISION ftfac
C--colour index common block
  common/colour/trip(23000),anti(23000),colmax
  INTEGER trip,anti,colmax
C--Parameter common block
  common/param/q0,lps,lqcd,ltime,scalefacm,angord,scatrecoil,
     &allhad,compress,nf
      INTEGER nf
  DOUBLE PRECISION q0,lqcd,ltime,lps,scalefacm
      LOGICAL angord,scatrecoil,allhad,compress
C--discard event flag
  common/disc/ndisc,nstrange,ngood,errcount,wdisc,discard
  LOGICAL discard
  INTEGER ndisc,nstrange,ngood,errcount
  double precision wdisc
C--variables for angular ordering
      common/angor/za(23000),zd(23000),thetaa(23000),qqbard(23000)
  DOUBLE PRECISION za,zd,thetaa
      LOGICAL qqbard
C--number of scattering events
  common/check/nscat,nscateff,nsplit
  DOUBLE PRECISION nscat,nscateff,nsplit
C--event weight
  common/weight/evweight,sumofweights
  double precision evweight,sumofweights

C--local variables
  INTEGER l,dir
  DOUBLE PRECISION phiq,pyr,pi,generatez,bmax1,cmax1,pts,mb,mc,
     &getmass,pz,eps,qh,z,r,lambda,weight,zdecb,zdecc,xdec(3),theta,
     &gettemp
      LOGICAL quark,qqbar,qqbardecb,qqbardecc
  integer bin
  DATA pi/3.141592653589793d0/

      IF((n+2).GT.22990) THEN
       write(logfid,*)'event too long for event record'
       discard=.true.
       RETURN
      ENDIF

      xdec(1)=mv(l,1)+(mv(l,5)-mv(l,4))*p(l,1)/p(l,4)
      xdec(2)=mv(l,2)+(mv(l,5)-mv(l,4))*p(l,2)/p(l,4)
      xdec(3)=mv(l,3)+(mv(l,5)-mv(l,4))*p(l,3)/p(l,4)
  IF(gettemp(xdec(1),xdec(2),xdec(3),mv(l,5)).GT.0.d0)THEN
   theta=-1.d0
  ELSE
   theta=thetaa(l)
  ENDIF 

C--on-shell partons cannot split
  IF((p(l,5).EQ.0d0).OR.(k(l,1).EQ.11).OR.(k(l,1).EQ.12)
     &  .OR.(k(l,1).EQ.13).OR.(k(l,1).EQ.14).OR.(k(l,1).EQ.3)
     &  .or.(zd(l).lt.0.d0)) goto 31
C--quark or gluon?
  IF(k(l,2).EQ.21)THEN
   quark=.false.
  ELSE
   quark=.true.
   qqbar=.false.
  ENDIF
C--if gluon decide on kind of splitting
  qqbar=qqbard(l)
C--if g->gg splitting decide on colour order
  IF(quark.OR.qqbar)THEN
   dir=0
  ELSE
   IF(pyr(0).LT.0.5)THEN
    dir=1
   ELSE
    dir=-1
   ENDIF
  ENDIF
  z=zd(l)
  IF(z.EQ.0.d0)THEN
   write(logfid,*)'makesplitting: z=0',l
   goto 36
  ENDIF  
  goto 35
C--generate z value
 36 IF(angord.AND.(za(l).NE.1.d0))THEN
C--additional z constraint due to angular ordering
   qh=4.*p(l,5)**2*(1.-za(l))/(za(l)*p(k(l,3),5)**2)
   IF(qh.GT.1)THEN
    write(logfid,*)l,': reject event: angular ordering
     &      conflict in medium'
    CALL pylist(3)
    discard=.true.
    goto 31
   ENDIF
   eps=0.5-0.5*sqrt(1.-qh)
  ELSE
   eps=0d0
  ENDIF
  IF(quark)THEN
   z=generatez(p(l,5)**2,p(l,4),eps,'QQ')
  ELSE
   IF(qqbar)THEN
    z=generatez(p(l,5)**2,p(l,4),eps,'QG')
   ELSE
    z=generatez(p(l,5)**2,p(l,4),eps,'GG')
   ENDIF
  ENDIF
 35 CONTINUE
C--maximum virtualities for daughters
  bmax1=min(p(l,5),z*p(l,4))
      cmax1=min(p(l,5),(1.-z)*p(l,4))
C--generate mass of quark or gluon (particle b) from Sudakov FF
 30 IF(quark.OR.qqbar)THEN
   mb=getmass(0.d0,bmax1,theta,z*p(l,4),'QQ',
     &      bmax1,.false.,zdecb,qqbardecb)
  ELSE
   mb=getmass(0.d0,bmax1,theta,z*p(l,4),'GC',
     &      bmax1,.false.,zdecb,qqbardecb)
  ENDIF
C--generate mass gluon (particle c) from Sudakov FF
  IF(quark.OR.(.NOT.qqbar))THEN
       mc=getmass(0.d0,cmax1,theta,(1.-z)*p(l,4),'GC',
     &  cmax1,.false.,zdecc,qqbardecc)
  ELSE
       mc=getmass(0.d0,cmax1,theta,(1.-z)*p(l,4),'QQ',
     &  cmax1,.false.,zdecc,qqbardecc)
  ENDIF
C--quark (parton b) momentum
 182  pz=(2.*z*p(l,4)**2-p(l,5)**2-mb**2+mc**2)/(2.*p(l,3))
  pts=z**2*(p(l,4)**2)-pz**2-mb**2
C--if kinematics doesn't work out, generate new virtualities
C     for daughters
C--massive phase space weight 
      IF((mb.EQ.0.d0).AND.(mc.EQ.0.d0).AND.(pts.LT.0.d0)) goto 36
  weight=1.d0
  IF((pyr(0).GT.weight).OR.(pts.LT.0.d0)
     &  .OR.((mb+mc).GT.p(l,5)))THEN
   IF(mb.GT.mc)THEN
    IF(quark.OR.qqbar)THEN
     mb=getmass(0.d0,mb,theta,z*p(l,4),'QQ',
     &      bmax1,.false.,zdecb,qqbardecb)
    ELSE
     mb=getmass(0.d0,mb,theta,z*p(l,4),'GC',
     &      bmax1,.false.,zdecb,qqbardecb)
    ENDIF
   ELSE
    IF(quark.OR.(.NOT.qqbar))THEN
         mc=getmass(0.d0,mc,theta,(1.-z)*p(l,4),'GC',
     &  cmax1,.false.,zdecc,qqbardecc)
    ELSE
         mc=getmass(0.d0,mc,theta,(1.-z)*p(l,4),'QQ',
     &  cmax1,.false.,zdecc,qqbardecc)
    ENDIF
   ENDIF
   goto 182
  ENDIF
  n=n+2
C--take care of first daughter (radiated gluon or antiquark)
  k(n-1,1)=k(l,1)
  IF(qqbar)THEN
   k(n-1,2)=-1
   trip(n-1)=0
   anti(n-1)=anti(l)
  ELSE
   k(n-1,2)=21
   IF((k(l,2).GT.0).AND.(dir.GE.0))THEN
    trip(n-1)=trip(l)
    anti(n-1)=colmax+1
   ELSE
    trip(n-1)=colmax+1
    anti(n-1)=anti(l)
   ENDIF
   colmax=colmax+1
  ENDIF
  k(n-1,3)=l
  k(n-1,4)=0
  k(n-1,5)=0
  p(n-1,4)=(1-z)*p(l,4)
  p(n-1,5)=mc
  za(n-1)=1.-z
  IF(zdecc.GT.0.d0)THEN
   thetaa(n-1)=p(n-1,5)/(sqrt(zdecc*(1.-zdecc))*p(n-1,4))
  ELSE
   thetaa(n-1)=0.d0
  ENDIF 
  zd(n-1)=zdecc
  qqbard(n-1)=qqbardecc
C--take care of second daughter (final quark or gluon or quark from 
C  gluon splitting)
  k(n,1)=k(l,1)
  IF(quark)THEN
   k(n,2)=k(l,2)
   IF(k(n,2).GT.0)THEN
    trip(n)=anti(n-1)
    anti(n)=0
   ELSE
    trip(n)=0
    anti(n)=trip(n-1)
   ENDIF
  ELSEIF(qqbar)THEN
   k(n,2)=1
   trip(n)=trip(l)
   anti(n)=0
  ELSE
   k(n,2)=21
   IF(dir.EQ.1)THEN
    trip(n)=anti(n-1)
    anti(n)=anti(l)
   ELSE
    trip(n)=trip(l)
    anti(n)=trip(n-1)
   ENDIF
  ENDIF
  k(n,3)=l
  k(n,4)=0
  k(n,5)=0
  p(n,3)=pz
  p(n,4)=z*p(l,4)
  p(n,5)=mb
  za(n)=z
  IF(zdecb.GT.0.d0)THEN
   thetaa(n)=p(n,5)/(sqrt(zdecb*(1.-zdecb))*p(n,4))
  ELSE 
   thetaa(n)=0.d0
  ENDIF 
  zd(n)=zdecb
  qqbard(n)=qqbardecb
C--azimuthal angle
  phiq=2*pi*pyr(0)
  p(n,1)=sqrt(pts)*cos(phiq)
  p(n,2)=sqrt(pts)*sin(phiq)
C--gluon momentum
  p(n-1,1)=p(l,1)-p(n,1)
  p(n-1,2)=p(l,2)-p(n,2)
  p(n-1,3)=p(l,3)-p(n,3)
      mv(n-1,4)=mv(l,5)
      IF(p(n-1,5).GT.0.d0)THEN
       lambda=1.d0/(ftfac*p(n-1,4)*0.2/p(n-1,5)**2)
   mv(n-1,5)=mv(l,5)-log(1.d0-pyr(0))/lambda
      ELSE
      mv(n-1,5)=0.d0
      ENDIF
      mv(n,4)=mv(l,5)
      IF(p(n,5).GT.0.d0)THEN
       lambda=1.d0/(ftfac*p(n,4)*0.2/p(n,5)**2)
   mv(n,5)=mv(l,5)-log(1.d0-pyr(0))/lambda
      ELSE
       mv(n,5)=0.d0
      ENDIF
C--take care of initial quark (or gluon)
      IF(k(l,1).EQ.2)THEN
       k(l,1)=13
      ELSE
   k(l,1)=11
      ENDIF
  k(l,4)=n-1
  k(l,5)=n
  nsplit=nsplit+evweight
 31 CONTINUE
  END


***********************************************************************
***   subroutine makeinsplit
***********************************************************************
  SUBROUTINE makeinsplit(L,X,TSUM,VIRT,TYPI,TIME,TAURAD)
  IMPLICIT NONE
C--identifier of file for hepmc output and logfile
  common/hepmcid/hpmcfid,logfid
  integer hpmcfid,logfid
C--Common block of Pythia
      common/pyjets/n,npad,k(23000,5),p(23000,5),v(23000,5)
  INTEGER n,npad,k
  DOUBLE PRECISION p,v
C--time common block
      common/time/mv(23000,5)
      DOUBLE PRECISION mv
C--factor in front of formation times
  common/ftimefac/ftfac
  DOUBLE PRECISION ftfac
C--colour index common block
  common/colour/trip(23000),anti(23000),colmax
  INTEGER trip,anti,colmax
C--variables for angular ordering
      common/angor/za(23000),zd(23000),thetaa(23000),qqbard(23000)
  DOUBLE PRECISION za,zd,thetaa
      LOGICAL qqbard
C--discard event flag
  common/disc/ndisc,nstrange,ngood,errcount,wdisc,discard
  LOGICAL discard
  INTEGER ndisc,nstrange,ngood,errcount
  double precision wdisc
C--Parameter common block
  common/param/q0,lps,lqcd,ltime,scalefacm,angord,scatrecoil,
     &allhad,compress,nf
      INTEGER nf
  DOUBLE PRECISION q0,lqcd,ltime,lps,scalefacm
      LOGICAL angord,scatrecoil,allhad,compress
C--number of scattering events
  common/check/nscat,nscateff,nsplit
  DOUBLE PRECISION nscat,nscateff,nsplit
C--event weight
  common/weight/evweight,sumofweights
  double precision evweight,sumofweights

C--local variables
  INTEGER l,typi,nold,dir
  DOUBLE PRECISION x,virt,mb2,mc2,getmass,pz,kt2,theta,phi,pi,
     &phiq,pyp,pyr,r,time,tsum,taurad,lambda,zdec
      LOGICAL qqbardec
  CHARACTER*2 typ2,typc
  integer bin
  DATA pi/3.141592653589793d0/

      IF((n+2).GT.22990) THEN
       write(logfid,*)'event too long for event record'
       discard=.true.
       RETURN
      ENDIF

  IF(k(l,2).EQ.21)THEN
   IF(typi.EQ.21)THEN
    typ2='GG'
    typc='GC'
   ELSE
    typ2='QG'
    typc='QQ'
   ENDIF
  ELSE
   IF(typi.EQ.21)THEN
    typ2='GQ'
    typc='QQ'
   ELSE
    typ2='QQ'
    typc='GC'
   ENDIF
  ENDIF

C--if g->gg decide on colour configuration
  IF(typ2.EQ.'GG')THEN
   IF(pyr(0).LT.0.5)THEN
    dir=1
   ELSE
    dir=-1
   ENDIF
  ELSE
   dir=0
  ENDIF

  mb2=virt**2
  mb2=p(l,5)**2-mb2
  mc2=getmass(0.d0,scalefacm*sqrt(-tsum),-1.d0,
     &  (1.-x)*p(l,4),typc,(1.-x)*p(l,4),
     &      .false.,zdec,qqbardec)**2

C--rotate such that momentum points in z-direction
      nold=n
      theta=pyp(l,13)
      phi=pyp(l,15)
      CALL pyrobo(l,l,0d0,-phi,0d0,0d0,0d0)
      CALL pyrobo(l,l,-theta,0d0,0d0,0d0,0d0)
  pz=(2*x*p(l,4)**2-p(l,5)**2-mb2+mc2)/(2*p(l,3))
  kt2=x**2*(p(l,4)**2)-pz**2-mb2
  IF(kt2.LT.0.d0)THEN
   mc2=0.d0
   pz=(2*x*p(l,4)**2-p(l,5)**2-mb2+mc2)/(2*p(l,3))
   kt2=x**2*(p(l,4)**2)-pz**2-mb2
   IF(kt2.LT.0.d0)THEN
        CALL pyrobo(l,l,theta,0d0,0d0,0d0,0d0)
        CALL pyrobo(l,l,0d0,phi,0d0,0d0,0d0)
        x=1.d0
    RETURN
   ENDIF
  ENDIF 
  n=n+2
C--take care of first daughter (radiated gluon or antiquark)
  k(n-1,1)=k(l,1)
  IF(typ2.EQ.'QG')THEN
   k(n-1,2)=-typi
   IF(k(n-1,2).GT.0)THEN
    trip(n-1)=trip(l)
    anti(n-1)=0
   ELSE
    trip(n-1)=0
    anti(n-1)=anti(l)
   ENDIF
  ELSEIF(typ2.EQ.'GQ')THEN
   k(n-1,2)=k(l,2)
       IF(k(n-1,2).GT.0)THEN
    trip(n-1)=colmax+1
    anti(n-1)=0
   ELSE
    trip(n-1)=0
    anti(n-1)=colmax+1
   ENDIF
   colmax=colmax+1
  ELSE
   k(n-1,2)=21
   IF((k(l,2).GT.0).AND.(dir.GE.0))THEN
    trip(n-1)=trip(l)
    anti(n-1)=colmax+1
   ELSE
    trip(n-1)=colmax+1
    anti(n-1)=anti(l)
   ENDIF
   colmax=colmax+1
  ENDIF
  k(n-1,3)=l
  k(n-1,4)=0
  k(n-1,5)=0
  p(n-1,4)=(1.-x)*p(l,4)
  p(n-1,5)=sqrt(mc2)
C--take care of second daughter (final quark or gluon or quark from 
C  gluon splitting)
  k(n,1)=k(l,1)
  IF(typ2.EQ.'QG')THEN
   k(n,2)=typi
   IF(k(n,2).GT.0)THEN
    trip(n)=trip(l)
    anti(n)=0
   ELSE
    trip(n)=0
    anti(n)=anti(l)
   ENDIF
  ELSEIF(typi.NE.21)THEN
   k(n,2)=k(l,2)
       IF(k(n,2).GT.0)THEN
    trip(n)=anti(n-1)
    anti(n)=0
   ELSE
    trip(n)=0
    anti(n)=trip(n-1)
   ENDIF
  ELSE
   k(n,2)=21
   IF(k(n-1,2).EQ.21)THEN
    IF(dir.EQ.1)THEN
     trip(n)=anti(n-1)
     anti(n)=anti(l)
    ELSE
     trip(n)=trip(l)
     anti(n)=trip(n-1)
    ENDIF
   ELSEIF(k(n-1,2).GT.0)THEN
    trip(n)=trip(l)
    anti(n)=trip(n-1)
   ELSE
    trip(n)=anti(n-1)
    anti(n)=anti(l)
   ENDIF
  ENDIF
  k(n,3)=l
  k(n,4)=0
  k(n,5)=0
  p(n,3)=pz
  p(n,4)=x*p(l,4)
  IF(mb2.LT.0.d0)THEN
   p(n,5)=-sqrt(-mb2)
  ELSE
   p(n,5)=sqrt(mb2)
  ENDIF
C--azimuthal angle
  phiq=2*pi*pyr(0)
  p(n,1)=sqrt(kt2)*cos(phiq)
  p(n,2)=sqrt(kt2)*sin(phiq)
C--gluon momentum
  p(n-1,1)=p(l,1)-p(n,1)
  p(n-1,2)=p(l,2)-p(n,2)
  p(n-1,3)=p(l,3)-p(n,3)
  mv(l,5)=time-taurad
      mv(n-1,4)=mv(l,5)
      IF(p(n-1,5).GT.0.d0)THEN
       lambda=1.d0/(ftfac*p(n-1,4)*0.2/p(n-1,5)**2)
   mv(n-1,5)=mv(l,5)-log(1.d0-pyr(0))/lambda
      ELSE
       mv(n-1,5)=0.d0
      ENDIF
      mv(n,4)=mv(l,5)
      IF(p(n,5).GT.0.d0)THEN
   mv(n,5)=time
      ELSE
       mv(n,5)=0.d0
      ENDIF
  za(n-1)=1.d0
      thetaa(n-1)=-1.d0
  zd(n-1)=zdec
  qqbard(n-1)=qqbardec
  za(n)=1.d0
  thetaa(n)=-1.d0
  zd(n)=0.d0
  qqbard(n)=.false.
C--take care of initial quark (or gluon)
      IF(k(l,1).EQ.2)THEN
       k(l,1)=13
      ELSE
   k(l,1)=11
      ENDIF
  k(l,4)=n-1
  k(l,5)=n
  nsplit=nsplit+evweight
      CALL pyrobo(l,l,theta,0d0,0d0,0d0,0d0)
      CALL pyrobo(n-1,n,theta,0d0,0d0,0d0,0d0)
      CALL pyrobo(l,l,0d0,phi,0d0,0d0,0d0)
      CALL pyrobo(n-1,n,0d0,phi,0d0,0d0,0d0)

C--set the production vertices: x_mother + (tprod - tprod_mother) * beta_mother
      mv(n-1,1)=mv(l,1)+(mv(n-1,4)-mv(l,4))*p(l,1)/max(pyp(l,8),p(l,4))
      mv(n-1,2)=mv(l,2)+(mv(n-1,4)-mv(l,4))*p(l,2)/max(pyp(l,8),p(l,4))
      mv(n-1,3)=mv(l,3)+(mv(n-1,4)-mv(l,4))*p(l,3)/max(pyp(l,8),p(l,4))
      mv(n,  1)=mv(l,1)+(mv(n,  4)-mv(l,4))*p(l,1)/max(pyp(l,8),p(l,4))
      mv(n,  2)=mv(l,2)+(mv(n,  4)-mv(l,4))*p(l,2)/max(pyp(l,8),p(l,4))
      mv(n,  3)=mv(l,3)+(mv(n,  4)-mv(l,4))*p(l,3)/max(pyp(l,8),p(l,4))

  END


***********************************************************************
***   subroutine doinstatescat
***********************************************************************
  SUBROUTINE doinstatescat(L,X,TYPI,Q,TSTART,DELTAT,OVERQ0,
     &        retrysplit)
  IMPLICIT NONE
C--Common block of Pythia
      common/pyjets/n,npad,k(23000,5),p(23000,5),v(23000,5)
  INTEGER n,npad,k
  DOUBLE PRECISION p,v
C--time common block
      common/time/mv(23000,5)
      DOUBLE PRECISION mv
C--factor in front of formation times
  common/ftimefac/ftfac
  DOUBLE PRECISION ftfac
C--Parameter common block
  common/param/q0,lps,lqcd,ltime,scalefacm,angord,scatrecoil,
     &allhad,compress,nf
      INTEGER nf
  DOUBLE PRECISION q0,lqcd,ltime,lps,scalefacm
      LOGICAL angord,scatrecoil,allhad,compress
C--discard event flag
  common/disc/ndisc,nstrange,ngood,errcount,wdisc,discard
  LOGICAL discard
  INTEGER ndisc,nstrange,ngood,errcount
  double precision wdisc
C--variables for coherent scattering
  common/coherent/nstart,nend,allqs(10000,6),scatcentres(10000,10),
     &qsumvec(4),qsum2
  INTEGER nstart,nend
  DOUBLE PRECISION allqs,scatcentres,qsumvec,qsum2
C--identifier of file for hepmc output and logfile
  common/hepmcid/hpmcfid,logfid
  integer hpmcfid,logfid
C--local variables
  INTEGER l,typi,counter,countmax,count2
  DOUBLE PRECISION x,deltat,deltal,pyr,r,pnorad,getpnorad1,getnoscat,
     &weight,low,fmax,getpdf,sigmatot,getsscat,pfchange,pi,tnow,tleft,
     &xmax,pqq,pqg,pgq,pgg,alphas,tstart,tsum,q,qold,q2old,getnewmass,
     &generatez,tmax,tmaxnew,dt,xsc,ysc,zsc,tsc,ms1,md1,getms,getmd,
     &gettemp,getneff,lambda,rtau,phi,tauest,qsumvecold(4),zdum,weight,
     &pyp
  LOGICAL fchange,norad,overq0,noscat,getdeltat,retrysplit,
     &qqbardum  
  CHARACTER typ
  CHARACTER*2 typ2
  DATA pi/3.141592653589793d0/
  DATA countmax/10000/

  counter=0
  
      xsc=mv(l,1)+(tstart-mv(l,4))*p(l,1)/p(l,4)
      ysc=mv(l,2)+(tstart-mv(l,4))*p(l,2)/p(l,4)
      zsc=mv(l,3)+(tstart-mv(l,4))*p(l,3)/p(l,4)
      tsc=tstart
      md1=getmd(xsc,ysc,zsc,tsc)
      ms1=getms(xsc,ysc,zsc,tsc)

      IF(md1.LE.1.d-4.OR.ms1.LE.1.d-4)THEN
       write(logfid,*)'problem!',gettemp(xsc,ysc,zsc,tsc),
     &getneff(xsc,ysc,zsc,tsc)
      ENDIF

C--check for scattering
      noscat=.NOT.getdeltat(l,tstart,deltat,dt)
  IF(noscat.AND.(.NOT.retrysplit)) goto 116

C--decide whether there will be radiation
  pnorad=getpnorad1(l,xsc,ysc,zsc,tsc)
  IF((pyr(0).LT.pnorad).OR.(p(l,4).LT.1.001*q0))THEN
   norad=.true.
  ELSE
   norad=.false.
  ENDIF

C--decide whether q or g is to be scattered
      IF(k(l,2).EQ.21)THEN
       typ='G'
       typ2='GC'
   sigmatot=getsscat(p(l,4),p(l,1),p(l,2),p(l,3),p(l,5),
     &  q0,'G','C',xsc,ysc,zsc,tsc,0)
   IF((sigmatot.EQ.0.d0).OR.(pnorad.EQ.1.d0))THEN
    pfchange=0.d0
   ELSE
    pfchange=getsscat(p(l,4),p(l,1),p(l,2),p(l,3),p(l,5),
     &  q0,'G','Q',xsc,ysc,zsc,tsc,0)
     &  /sigmatot
   ENDIF
   sigmatot=getsscat(p(l,4),p(l,1),p(l,2),p(l,3),p(l,5),
     &  0.d0,'G','C',xsc,ysc,zsc,tsc,0)
      ELSE
       typ='Q'
       typ2='QQ'
   sigmatot=getsscat(p(l,4),p(l,1),p(l,2),p(l,3),p(l,5),
     &  q0,'Q','C',xsc,ysc,zsc,tsc,0)
   IF((sigmatot.EQ.0.d0).OR.(pnorad.EQ.1.d0))THEN
    pfchange=0.d0
   ELSE
    pfchange=getsscat(p(l,4),p(l,1),p(l,2),p(l,3),p(l,5),
     &  q0,'Q','G',xsc,ysc,zsc,tsc,0)
     &  /sigmatot
   ENDIF
   sigmatot=getsscat(p(l,4),p(l,1),p(l,2),p(l,3),p(l,5),
     &  0.d0,'Q','C',xsc,ysc,zsc,tsc,0)
      ENDIF
  IF((pfchange.LT.-1.d-4).OR.(pfchange.GT.1.d0+1.d-4)) THEN
      write(logfid,*)'error: flavour change probability=',
     &  pfchange,'for ',typ
  ENDIF
  IF(pyr(0).LT.pfchange)THEN
   fchange=.true.
  ELSE
   fchange=.false.
  ENDIF
      IF (norad) fchange=.false.
C--set TYPI
  IF(typ.EQ.'G')THEN
   IF(fchange)THEN
    typi=int(sign(2.d0,pyr(0)-0.5))
   ELSE
    typi=k(l,2)
   ENDIF
  ELSE
   IF(fchange)THEN
    typi=21
   ELSE
    typi=k(l,2)
   ENDIF
  ENDIF
  low=q0**2/scalefacm**2
  tmax=4.*(p(l,4)**2-p(l,5)**2)
  xmax=1.-q0**2/(scalefacm**2*4.*tmax)

  IF(sigmatot.EQ.0.d0) goto 116

  rtau=pyr(0)

C--generate a trial emission
C--pick a x value from splitting function
 112  counter=counter+1
  IF(typ.EQ.'G')THEN
   IF(fchange)THEN
    x=generatez(0.d0,0.d0,1.-xmax,'QG')
   ELSE
    x=generatez(0.d0,0.d0,1.-xmax,'GG')
   ENDIF
  ELSE
   IF(fchange)THEN
    x=1.-generatez(0.d0,0.d0,1.-xmax,'QQ')
   ELSE
    x=generatez(0.d0,0.d0,1.-xmax,'QQ')
   ENDIF
  ENDIF
      IF(norad) x=1.d0
C--initialisation
      tmaxnew=(x*p(l,4))**2
  phi=0.d0
  tleft=deltat
  tnow=tstart
  qsumvec(1)=0.d0
  qsumvec(2)=0.d0
  qsumvec(3)=0.d0
  qsumvec(4)=0.d0
  qsum2=-1.d-10
  overq0=.false.
  q=p(l,5)
  qold=p(l,5)
      tauest=deltat
C--generate first momentum transfer
  deltal=dt
  nstart=1
  nend=1
  tnow=tnow+deltal
  tsum=deltal
  tleft=tleft-deltal
  allqs(nend,6)=tnow
  q2old=qsum2
C--get new momentum transfer
  count2=0
 118  CALL getqvec(l,nend,tnow-mv(l,4),x)
  IF(-qsum2.GT.p(l,4)**2)THEN
   qsumvec(1)=0.d0
   qsumvec(2)=0.d0
   qsumvec(3)=0.d0
   qsumvec(4)=0.d0
   qsum2=q2old
   IF(count2.LT.100)THEN
    count2=count2+1
    goto 118
   ELSE
    allqs(nend,1)=0.d0
    allqs(nend,2)=0.d0
    allqs(nend,3)=0.d0
    allqs(nend,4)=0.d0
    allqs(nend,5)=0.d0
   ENDIF
  ENDIF
C--update OVERQ0
  IF(-allqs(nend,1).GT.low) overq0=.true.
C--get new virtuality
   IF(overq0.AND.(.NOT.norad))THEN
    q=getnewmass(l,scalefacm**2*qsum2,scalefacm**2*q2old,0.d0,
     &    .true.,x,zdum,qqbardum)
   ELSE
    q=0.d0
   ENDIF

C--estimate formation time
 111  IF((q.EQ.0.d0).OR.(q.EQ.p(l,5)))THEN
   tauest=deltat
  ELSE
   tauest=ftfac*(1.-phi)*0.2*x*p(l,4)/q**2
  ENDIF
  lambda=1.d0/tauest
  tauest=-log(1.d0-rtau)/lambda

C--find number, position and momentum transfers of further scatterings
  noscat=.NOT.getdeltat(l,tnow,min(tleft,tauest),deltal)
  IF((.NOT.noscat).AND.(.NOT.norad))THEN
C--add a momentum transfer
   nend=nend+1
   IF(nend.GE.100)THEN
    nend=nend-1
    goto 114
   ENDIF
   tnow=tnow+deltal
   tsum=tsum+deltal
   tleft=tleft-deltal
C--update phase
   IF((q.NE.0.d0).AND.(q.NE.p(l,5)))THEN
    phi=phi+5.*deltal*q**2/(1.*x*p(l,4))
   ENDIF
C--get new momentum transfer
   allqs(nend,6)=tnow
   q2old=qsum2
   qsumvecold(1)=qsumvec(1)
   qsumvecold(2)=qsumvec(2)
   qsumvecold(3)=qsumvec(3)
   qsumvecold(4)=qsumvec(4)
   count2=0
 119   CALL getqvec(l,nend,tnow-mv(l,4),x)
   IF(-qsum2.GT.p(l,4)**2)THEN
    qsumvec(1)=qsumvecold(1)
    qsumvec(2)=qsumvecold(2)
    qsumvec(3)=qsumvecold(3)
    qsumvec(4)=qsumvecold(4)
    qsum2=q2old
    IF(count2.LT.100)THEN
     count2=count2+1
     goto 119
    ELSE
     allqs(nend,1)=0.d0
     allqs(nend,2)=0.d0
     allqs(nend,3)=0.d0
     allqs(nend,4)=0.d0
     allqs(nend,5)=0.d0
    ENDIF
   ENDIF
C--update OVERQ0
   IF((-qsum2.GT.low)
     &  .OR.(-allqs(nend,1).GT.low)) overq0=.true.
C--get new virtuality
   qold=q
   IF(overq0.AND.(.NOT.norad))THEN
    q=getnewmass(l,scalefacm**2*qsum2,scalefacm**2*q2old,0.d0,
     &    .true.,x,zdum,qqbardum)
   ELSE
    q=0.d0
   ENDIF
   goto 111
  ENDIF

C--do reweighting
 114  tmaxnew=x**2*p(l,4)**2
  IF(norad)THEN
   weight=1.d0
   q=0.d0
   x=1.d0
  ELSEIF((-qsum2.LT.low).OR.(q.EQ.0.d0))THEN
   weight=0.d0
  ELSEIF(-qsum2.GT.p(l,4)**2)THEN
   weight=0.d0
  ELSE   
   IF(typ.EQ.'G')THEN
    fmax=2.*log(-scalefacm**2*qsum2/q0**2)
     &    *alphas(q0**2/4.,lps)/(2.*pi)
    IF(qsum2.EQ.0.d0)THEN
     weight=0.d0
     norad=.true.
    ELSE
     IF(fchange)THEN
      weight=2.*getpdf(x,scalefacm*sqrt(-qsum2),'QG')/(pqg(x)*fmax)
      IF((weight.GT.1.d0+1.d-4).OR.(weight.LT.-1.d-4))THEN
        write(logfid,*)'x,sqrt(qsum^2),getpdf,fmax:',x,
     &  sqrt(-qsum2),getpdf(x,scalefacm*sqrt(-qsum2),'QG'),'qg',
     &  fmax
          ENDIF
     ELSE
      weight=getpdf(x,scalefacm*sqrt(-qsum2),'GG')/(pgg(x)*fmax)
      IF((weight.GT.1.d0+1.d-4).OR.(weight.LT.-1.d-4))THEN
        write(logfid,*)'x,sqrt(qsum^2),getpdf,fmax:',x,
     &  sqrt(-qsum2),getpdf(x,scalefacm*sqrt(-qsum2),'GG'),'gg',
     &  fmax
          ENDIF
     ENDIF
    ENDIF
   ELSE
    fmax=log(-scalefacm**2*qsum2/q0**2)
     &    *alphas(q0**2/4.,lps)/(2.*pi)
    IF(qsum2.EQ.0.d0)THEN
     weight=0.d0
     norad=.true.
    ELSE
     IF(fchange)THEN
      weight=getpdf(x,scalefacm*sqrt(-qsum2),'GQ')/(pgq(x)*fmax)
      IF((weight.GT.1.d0+1.d-4).OR.(weight.LT.-1.d-4))THEN
       write(logfid,*)'x,sqrt(qsum^2),getpdf:,fmax',x,
     &  sqrt(-qsum2),getpdf(x,scalefacm*sqrt(-qsum2),'GQ'),'gq',
     &  fmax
          ENDIF
     ELSE
      weight=getpdf(x,scalefacm*sqrt(-qsum2),'QQ')/(pqq(x)*fmax)
      IF((weight.GT.1.d0+1.d-4).OR.(weight.LT.-1.d-4))THEN
       write(logfid,*)'x,sqrt(qsum^2),getpdf,fmax:',x,
     &  sqrt(-qsum2),getpdf(x,scalefacm*sqrt(-qsum2),'QQ'),'qq',
     &  fmax
          ENDIF
     ENDIF
    ENDIF
   ENDIF
  ENDIF
  IF((weight.GT.1.d0+1.d-4).OR.(weight.LT.-1.d-4))
     &  write(logfid,*)'error: weight=',weight
 115  IF(pyr(0).GT.weight)THEN
   IF(counter.LT.countmax)THEN
    goto 112
   ELSE
    q=0.d0
    x=1.d0
    nend=nstart
    qsum2=allqs(nend,1)
    qsumvec(1)=allqs(nend,2)
    qsumvec(2)=allqs(nend,3)
    qsumvec(3)=allqs(nend,4)
    qsumvec(4)=allqs(nend,5)
    typi=k(l,2)
    IF(-allqs(nend,1).GT.low)THEN
     overq0=.true.
    ELSE
     overq0=.false.
    ENDIF
        deltat=allqs(nend,6)-tstart
    tnow=allqs(1,6)
    RETURN
   ENDIF
  ENDIF
C--found meaningful configuration, now do final checks
C--check if phase is unity and weight with 1/Nscat
      IF(((tleft.LT.tauest).OR.(pyr(0).GT.1.d0/(nend*1.d0)))
     &      .AND.(.NOT.norad))THEN
   q=0.d0
   x=1.d0
   nend=nstart
   qsum2=allqs(nend,1)
   qsumvec(1)=allqs(nend,2)
   qsumvec(2)=allqs(nend,3)
   qsumvec(3)=allqs(nend,4)
   qsumvec(4)=allqs(nend,5)
   typi=k(l,2)
   IF(-allqs(nend,1).GT.low)THEN
    overq0=.true.
   ELSE
    overq0=.false.
   ENDIF
       deltat=allqs(nend,6)-tstart
   tnow=allqs(1,6)
  ELSE
       IF(.NOT.norad)THEN
    tleft=tleft-tauest
    tnow=tnow+tauest
    tsum=tsum+tauest
   ENDIF
       deltat=tsum
  ENDIF
  RETURN
C--exit in case of failure
 116  q=0.d0
  x=1.d0
  nstart=0
  nend=0
  qsumvec(1)=0.d0
  qsumvec(2)=0.d0
  qsumvec(3)=0.d0
  qsumvec(4)=0.d0
  qsum2=0.d0
  overq0=.false.
  typi=k(l,2)
  RETURN
  END


***********************************************************************
***   subroutine dofistatescat
***********************************************************************
  SUBROUTINE dofistatescat(L,TNOW,DTLEFT,DELTAT,NEWMASS,
     &    overq0,z,qqbar)
  IMPLICIT NONE
C--identifier of file for hepmc output and logfile
  common/hepmcid/hpmcfid,logfid
  integer hpmcfid,logfid
C--Common block of Pythia
      common/pyjets/n,npad,k(23000,5),p(23000,5),v(23000,5)
  INTEGER n,npad,k
  DOUBLE PRECISION p,v
C--time common block
      common/time/mv(23000,5)
      DOUBLE PRECISION mv
C--factor in front of formation times
  common/ftimefac/ftfac
  DOUBLE PRECISION ftfac
C--Parameter common block
  common/param/q0,lps,lqcd,ltime,scalefacm,angord,scatrecoil,
     &allhad,compress,nf
      INTEGER nf
  DOUBLE PRECISION q0,lqcd,ltime,lps,scalefacm
      LOGICAL angord,scatrecoil,allhad,compress
C--discard event flag
  common/disc/ndisc,nstrange,ngood,errcount,wdisc,discard
  LOGICAL discard
  INTEGER ndisc,nstrange,ngood,errcount
  double precision wdisc
C--variables for coherent scattering
  common/coherent/nstart,nend,allqs(10000,6),scatcentres(10000,10),
     &qsumvec(4),qsum2
  INTEGER nstart,nend
  DOUBLE PRECISION allqs,scatcentres,qsumvec,qsum2
C--local variables
  INTEGER l,counter,countmax,count2
  DOUBLE PRECISION tnow,deltat,newmass,tleft,deltal,q2old,
     &getnewmass,pyr,tsum,qsumvecold(4),rtau,lambda,dtleft,phi,
     &tauest,low,z,pyp
  LOGICAL overq0,noscat,getdeltat,qqbar
  CHARACTER typ
  DATA countmax/100/
  deltal=0.d0

  IF(-qsum2.GT.p(l,4)**2)
     & write(logfid,*) 'DOFISTATESCAT has a problem:',-qsum2,p(l,4)**2

      IF(k(l,2).EQ.21)THEN
       typ='G'
  ELSE
   typ='Q'
  ENDIF
  low=q0**2/scalefacm**2

  tsum=0.d0
  phi=0.d0
  deltat=0.d0

C--check for radiation with first (given) momentum transfer
  q2old=0.d0
  IF(overq0.OR.(-qsum2.GT.low))THEN
   newmass=getnewmass(l,scalefacm**2*qsum2,scalefacm**2*q2old,
     &  newmass,.false.,1.d0,z,qqbar)
   overq0=.true.
  ELSE
   newmass=p(l,5)
  ENDIF

  rtau=pyr(0)

  tleft=dtleft
 222  IF((newmass.EQ.0.d0).OR.(newmass.EQ.p(l,5)))THEN
   tauest=tleft
  ELSE
   tauest=ftfac*(1.-phi)*0.2*p(l,4)/newmass**2
  ENDIF
  lambda=1.d0/tauest
  tauest=-log(1.d0-rtau)/lambda
      noscat=.NOT.getdeltat(l,tnow+tsum,min(tauest,tleft),deltal)
  IF(.NOT.noscat)THEN
C--do scattering
   nend=nend+1
   IF(nend.gt.countmax)THEN
    nend=nend-1
    goto 218
   ENDIF
   IF(nstart.EQ.0) nstart=1
   tsum=tsum+deltal
   tleft=tleft-deltal
   IF((newmass.NE.0.d0).AND.(newmass.NE.p(l,5)))THEN
    phi=phi+5.*deltal*newmass**2/(1.*p(l,4))
   ENDIF
   allqs(nend,6)=tnow+tsum
   qsumvecold(1)=qsumvec(1)
   qsumvecold(2)=qsumvec(2)
   qsumvecold(3)=qsumvec(3)
   qsumvecold(4)=qsumvec(4)
   q2old=qsum2
C--get new momentum transfer
   count2=0
 219   CALL getqvec(l,nend,tnow+tsum-mv(l,4),1.d0)
   IF(-qsum2.GT.p(l,4)**2)THEN
    qsumvec(1)=qsumvecold(1)
    qsumvec(2)=qsumvecold(2)
    qsumvec(3)=qsumvecold(3)
    qsumvec(4)=qsumvecold(4)
    qsum2=q2old
    IF(count2.LT.100)THEN
     count2=count2+1
     goto 219
    ELSE
     allqs(nend,1)=0.d0
     allqs(nend,2)=0.d0
     allqs(nend,3)=0.d0
     allqs(nend,4)=0.d0
     allqs(nend,5)=0.d0
    ENDIF
   ENDIF
C--figure out new virtuality
   IF(overq0.OR.(-qsum2.GT.low))THEN
    newmass=getnewmass(l,scalefacm**2*qsum2,scalefacm**2*q2old,
     &    newmass,.false.,1.d0,z,qqbar)
    overq0=.true.
   ENDIF
   goto 222
  ENDIF
C--no more scattering
 218  if ((newmass**2.gt.low).and.(newmass.ne.p(l,5))) then
    if ((tleft.LT.tauest).OR.(pyr(0).GT.1.d0/(nend*1.d0))) then
      if (nend.eq.countmax) then
        deltat=tsum
      else if (tleft.LT.tauest) then
        deltat=tsum+tleft
      else
        deltat=tsum+tauest
      endif
      newmass=p(l,5)
    ELSE
      deltat=tsum+tauest
    ENDIF
  else  
    deltat=0.d0
    nstart=1
    nend=1
    qsum2=allqs(nend,1)
    qsumvec(1)=allqs(nend,2)
    qsumvec(2)=allqs(nend,3)
    qsumvec(3)=allqs(nend,4)
    qsumvec(4)=allqs(nend,5)
    IF(-allqs(nend,1).GT.low)THEN
      overq0=.true.
    ELSE
      overq0=.false.
    ENDIF
    newmass=p(l,5)
  endif
  return
  END


***********************************************************************
***   function getnewmass
***********************************************************************
  DOUBLE PRECISION FUNCTION getnewmass(L,Q2,QOLD2,MASS,IN,X,
     &  zdec,qqbardec)
  IMPLICIT NONE
C--Common block of Pythia
      common/pyjets/n,npad,k(23000,5),p(23000,5),v(23000,5)
  INTEGER n,npad,k
  DOUBLE PRECISION p,v
C--time common block
      common/time/mv(23000,5)
      DOUBLE PRECISION mv
C--variables for angular ordering
      common/angor/za(23000),zd(23000),thetaa(23000),qqbard(23000)
  DOUBLE PRECISION za,zd,thetaa
      LOGICAL qqbard
C--Parameter common block
  common/param/q0,lps,lqcd,ltime,scalefacm,angord,scatrecoil,
     &allhad,compress,nf
      INTEGER nf
  DOUBLE PRECISION q0,lqcd,ltime,lps,scalefacm
      LOGICAL angord,scatrecoil,allhad,compress
C--local variables
  INTEGER l
  DOUBLE PRECISION q2,qold2,r,pyr,pnosplit1,pnosplit2,z,qa,
     &getsudakov,getmass,pkeep,x,mass,zdec,qtmp,zold
  LOGICAL in,qqbardec,qqbarold
  CHARACTER*2 typ 

  IF(x*p(l,4).LT.q0)THEN
   getnewmass=0.d0
   zdec=0.d0
   qqbardec=.false.
   RETURN
  ENDIF
  IF (-q2.LT.q0**2)THEN
   getnewmass=0.d0
   RETURN
  ENDIF
      IF(k(l,2).EQ.21)THEN
       typ='GC'
      ELSE
       typ='QQ'
      ENDIF
  IF(sqrt(-qold2).LE.q0)THEN
   IF(in)THEN
    getnewmass=getmass(0.d0,sqrt(-q2),-1.d0,
     &  x*p(l,4),typ,x*p(l,4),in,zdec,qqbardec)
   ELSE
    getnewmass=getmass(0.d0,sqrt(-q2),-1.d0,p(l,4),typ,
     &  sqrt(-q2),in,zdec,qqbardec)
   ENDIF
   getnewmass=min(getnewmass,x*p(l,4))
   RETURN
  ENDIF
  z=1.d0
  qa=1.d0 
  IF(max(p(l,5),mass).GT.0.d0)THEN
     IF(-q2.GT.-qold2)THEN
        zold=zdec
        qqbarold=qqbardec
        qtmp=getmass(0.d0,sqrt(-q2),-1.d0,x*p(l,4),typ,
     &    sqrt(-q2),in,zdec,qqbardec)
        IF(qtmp.LT.sqrt(-qold2))THEN
          getnewmass=mass
          zdec=zold
              qqbardec=qqbarold
        ELSE
           getnewmass=qtmp
        ENDIF
     ELSE
       pnosplit1=getsudakov(sqrt(-qold2),qa,q0,z,x*p(l,4),
     &      typ,mv(l,4),in)
       pnosplit2=getsudakov(sqrt(-q2),qa,q0,z,x*p(l,4),
     &      typ,mv(l,4),in)
       pkeep=(1.-pnosplit2)/(1.-pnosplit1)
       IF(pyr(0).LT.pkeep)THEN
         IF(p(l,5).LT.sqrt(-q2))THEN
       getnewmass=mass
     ELSE
 55      getnewmass=getmass(q0,sqrt(-q2),-1.d0,x*p(l,4),typ,
     &    sqrt(-q2),in,zdec,qqbardec)
       IF((getnewmass.EQ.0.d0).AND.(x*p(l,4).GT.q0)) goto 55
     ENDIF
       ELSE
         getnewmass=0.d0
         zdec=0.d0
         qqbardec=.false.
       ENDIF
     ENDIF
   ELSE
     IF(-q2.GT.-qold2)THEN
       getnewmass=getmass(0.d0,sqrt(-q2),-1.d0,
     &        x*p(l,4),typ,x*p(l,4),in,zdec,qqbardec)
           if(getnewmass.lt.sqrt(-qold2))then
         getnewmass=0.d0
         zdec=0.d0
         qqbardec=.false.
           endif
     ELSE
       getnewmass=0.d0
       zdec=0.d0
       qqbardec=.false.
     ENDIF
   ENDIF
   getnewmass=min(getnewmass,x*p(l,4))
  END 


***********************************************************************
***   function getpnorad1
***********************************************************************
  DOUBLE PRECISION FUNCTION getpnorad1(LINE,x,y,z,t)
  IMPLICIT NONE
C--identifier of file for hepmc output and logfile
  common/hepmcid/hpmcfid,logfid
  integer hpmcfid,logfid
C--Common block of Pythia
      common/pyjets/n,npad,k(23000,5),p(23000,5),v(23000,5)
  INTEGER n,npad,k
  DOUBLE PRECISION p,v
C--Parameter common block
  common/param/q0,lps,lqcd,ltime,scalefacm,angord,scatrecoil,
     &allhad,compress,nf
      INTEGER nf
  DOUBLE PRECISION q0,lqcd,ltime,lps,scalefacm
      LOGICAL angord,scatrecoil,allhad,compress
C--local variables
  INTEGER line
  DOUBLE PRECISION up,low,ccol,sigmatot,getsscat,getxsecint,
     &scatprimfunc,ms1,md1,shat,pcms2,avmom(5),x,y,z,t,getmd
  
  md1 = getmd(x,y,z,t)
  call avscatcen(x,y,z,t,
     &avmom(1),avmom(2),avmom(3),avmom(4),avmom(5))
  ms1 = avmom(5)
  shat = avmom(5)**2 + p(line,5)**2 + 2.*(avmom(4)*p(line,4)
     &       -avmom(1)*p(line,1)-avmom(2)*p(line,2)-avmom(3)*p(line,3))
  pcms2 = (shat+p(line,5)**2-ms1**2)**2/(4.*shat)-p(line,5)**2
  up = 4.*pcms2
   low=q0**2/scalefacm**2
   IF((up.LE.low).OR.(p(line,4).LT.q0/scalefacm))THEN
    getpnorad1=1.d0
    RETURN
   ENDIF
   IF(k(line,2).EQ.21)THEN
    ccol=3./2.
C--probability for no initial state radiation
    sigmatot=getsscat(p(line,4),p(line,1),p(line,2),p(line,3),
     &    p(line,5),0.d0,'G','C',x,y,z,t,0)
    IF(sigmatot.EQ.0.d0)THEN
     getpnorad1=-1.d0
     RETURN
    ENDIF
     getpnorad1=(ccol*(scatprimfunc(low,md1)-
     &scatprimfunc(0.d0,md1))
     &    + getxsecint(up,md1,'GB'))/sigmatot
   ELSE
    ccol=2./3.
C--probability for no initial state radiation
    sigmatot=getsscat(p(line,4),p(line,1),p(line,2),p(line,3),
     &    p(line,5),0.d0,'Q','C',x,y,z,t,0)
    IF(sigmatot.EQ.0.d0)THEN
     getpnorad1=1.d0
     RETURN
    ENDIF
     getpnorad1=(ccol*(scatprimfunc(low,md1)-
     &scatprimfunc(0.d0,md1))
     &    + getxsecint(up,md1,'QB'))/sigmatot
   ENDIF
  IF((getpnorad1.LT.-1.d-4).OR.(getpnorad1.GT.1.d0+1.d-4))THEN
       write(logfid,*)'error: P_norad=',getpnorad1,
     &  p(line,4),p(line,5),low,up,k(line,2),md1
  ENDIF
  END


***********************************************************************
***   subroutine getqvec
***********************************************************************
  SUBROUTINE getqvec(L,J,DT,X)
  IMPLICIT NONE
C--identifier of file for hepmc output and logfile
  common/hepmcid/hpmcfid,logfid
  integer hpmcfid,logfid
C--Common block of Pythia
      common/pyjets/n,npad,k(23000,5),p(23000,5),v(23000,5)
  INTEGER n,npad,k
  DOUBLE PRECISION p,v
C--time common block
      common/time/mv(23000,5)
      DOUBLE PRECISION mv
C--variables for coherent scattering
  common/coherent/nstart,nend,allqs(10000,6),scatcentres(10000,10),
     &qsumvec(4),qsum2
  INTEGER nstart,nend
  DOUBLE PRECISION allqs,scatcentres,qsumvec,qsum2
C--discard event flag
  common/disc/ndisc,nstrange,ngood,errcount,wdisc,discard
  LOGICAL discard
  INTEGER ndisc,nstrange,ngood,errcount
  double precision wdisc
C--Parameter common block
  common/param/q0,lps,lqcd,ltime,scalefacm,angord,scatrecoil,
     &allhad,compress,nf
      INTEGER nf
  DOUBLE PRECISION q0,lqcd,ltime,lps,scalefacm
      LOGICAL angord,scatrecoil,allhad,compress
C--local variables
  INTEGER l,j,counter,countmax,count2,i
      DOUBLE PRECISION xsc,ysc,zsc,tsc,getmd,gettemp,dt,x,pyr,newmom(4),
     &t,pt,maxt,phi2,beta(3),phi,theta,gett,pyp,pi,pt2,getms,
     &savemom(5),theta2,mb2,pz,kt2,phiq,maxt2,xi,md,shat,pcms2,
     &avmom(5)
  CHARACTER typs
  DATA pi/3.141592653589793d0/
  DATA countmax/1000/

      IF (j.GT.10000)THEN
       discard = .true.
   return
      ENDIF

  counter=0
  count2=0

      xsc=mv(l,1)+dt*p(l,1)/p(l,4)
      ysc=mv(l,2)+dt*p(l,2)/p(l,4)
      zsc=mv(l,3)+dt*p(l,3)/p(l,4)
      tsc=mv(l,4)+dt
  md = getmd(xsc,ysc,zsc,tsc)

  call avscatcen(xsc,ysc,zsc,tsc,
     &avmom(1),avmom(2),avmom(3),avmom(4),avmom(5))

  do 210 i=1,5
    savemom(i) = p(l,i)
 210  continue

  xi = sqrt(max(x**2*p(l,4)**2,p(l,5)**2) - p(l,5)**2)/pyp(l,8)
  p(l,1) = xi*p(l,1)
  p(l,2) = xi*p(l,2)
  p(l,3) = xi*p(l,3)
  p(l,4) = max(x*p(l,4),p(l,5))


 444  CALL getscatterer(xsc,ysc,zsc,tsc,
     &k(1,2),p(1,1),p(1,2),p(1,3),p(1,4),p(1,5))
      mv(1,1)=xsc
      mv(1,2)=ysc
      mv(1,3)=zsc
      mv(1,4)=tsc
      typs='Q'
      IF(k(1,2).EQ.21)typs='G'

  shat = avmom(5)**2 + savemom(5)**2 + 2.*(avmom(4)*savemom(4)
     &    -avmom(1)*savemom(1)-avmom(2)*savemom(2)-avmom(3)*savemom(3))
  pcms2 = (shat+savemom(5)**2-avmom(5)**2)**2/(4.*shat)
     &  -savemom(5)**2
  maxt = 4.*pcms2

      k(1,1)=13
  scatcentres(j,1)=k(1,2)
  scatcentres(j,2)=p(1,1)
  scatcentres(j,3)=p(1,2)
  scatcentres(j,4)=p(1,3)
  scatcentres(j,5)=p(1,4)
  scatcentres(j,6)=p(1,5)
  scatcentres(j,7)=mv(1,1)
  scatcentres(j,8)=mv(1,2)
  scatcentres(j,9)=mv(1,3)
  scatcentres(j,10)=mv(1,4)
C--transform to scattering centre's rest frame and rotate such that parton momentum is in z-direction
      beta(1)=p(1,1)/p(1,4)
      beta(2)=p(1,2)/p(1,4)
      beta(3)=p(1,3)/p(1,4)
      CALL pyrobo(l,l,0d0,0d0,-beta(1),-beta(2),-beta(3))
      CALL pyrobo(1,1,0d0,0d0,-beta(1),-beta(2),-beta(3))
      theta=pyp(l,13)
      phi=pyp(l,15)
      CALL pyrobo(l,l,0d0,-phi,0d0,0d0,0d0)
      CALL pyrobo(1,1,0d0,-phi,0d0,0d0,0d0)
      CALL pyrobo(l,l,-theta,0d0,0d0,0d0,0d0)
      CALL pyrobo(1,1,-theta,0d0,0d0,0d0,0d0)
C--pick a t from differential scattering cross section
 204  t=-gett(0.d0,maxt,md)
 202  newmom(4)=p(l,4)+t/(2.*p(1,5))
  newmom(3)=(t-2.*p(l,5)**2+2.*p(l,4)*newmom(4))/(2.*p(l,3))
  pt2=newmom(4)**2-newmom(3)**2-p(l,5)**2
  IF(dabs(pt2).LT.1.d-10) pt2=0.d0  
  IF(t.EQ.0.d0) pt2=0.d0
  IF(pt2.LT.0.d0)THEN
   t=0.d0
   goto 202
  ENDIF
  pt=sqrt(pt2)
      phi2=pyr(0)*2*pi
  newmom(1)=pt*cos(phi2)
  newmom(2)=pt*sin(phi2)
  p(1,1)=newmom(1)-p(l,1)
  p(1,2)=newmom(2)-p(l,2)
  p(1,3)=newmom(3)-p(l,3)
  p(1,4)=newmom(4)-p(l,4)
  p(1,5)=0.d0
C--transformation to lab
      CALL pyrobo(l,l,theta,0d0,0d0,0d0,0d0)
      CALL pyrobo(1,1,theta,0d0,0d0,0d0,0d0)
      CALL pyrobo(l,l,0d0,phi,0d0,0d0,0d0)
      CALL pyrobo(1,1,0d0,phi,0d0,0d0,0d0)
      CALL pyrobo(l,l,0d0,0d0,beta(1),beta(2),beta(3))
      CALL pyrobo(1,1,0d0,0d0,beta(1),beta(2),beta(3))
  allqs(j,1)=t
  allqs(j,2)=p(1,1)
  allqs(j,3)=p(1,2)
  allqs(j,4)=p(1,3)
  allqs(j,5)=p(1,4)
  qsumvec(1)=qsumvec(1)+allqs(nend,2)
  qsumvec(2)=qsumvec(2)+allqs(nend,3)
  qsumvec(3)=qsumvec(3)+allqs(nend,4)
  qsumvec(4)=qsumvec(4)+allqs(nend,5)
  qsum2=qsumvec(4)**2-qsumvec(1)**2-qsumvec(2)**2-qsumvec(3)**2
  IF(qsum2.GT.0.d0)THEN
   qsumvec(1)=qsumvec(1)-allqs(nend,2)
   qsumvec(2)=qsumvec(2)-allqs(nend,3)
   qsumvec(3)=qsumvec(3)-allqs(nend,4)
   qsumvec(4)=qsumvec(4)-allqs(nend,5)
   qsum2=qsumvec(4)**2-qsumvec(1)**2-qsumvec(2)**2-qsumvec(3)**2
   IF(counter.GT.countmax)THEN
    write(logfid,*)'GETQVEC unable to find q vector'
    allqs(j,1)=0.d0
    allqs(j,2)=0.d0
    allqs(j,3)=0.d0
    allqs(j,4)=0.d0
    allqs(j,5)=0.d0
   ELSE
    counter=counter+1
    goto 444
   ENDIF
  ENDIF
  do 211 i=1,5
    p(l,i) = savemom(i)
 211  continue
  END

***********************************************************************
***   subroutine dokinematics
***********************************************************************
      SUBROUTINE dokinematics(L,lold,N1,N2,NEWM,RETRYSPLIT,
     &  time,x,z,qqbar)
      IMPLICIT NONE
C--identifier of file for hepmc output and logfile
  common/hepmcid/hpmcfid,logfid
  integer hpmcfid,logfid
C--Common block of Pythia
      common/pyjets/n,npad,k(23000,5),p(23000,5),v(23000,5)
  INTEGER n,npad,k
  DOUBLE PRECISION p,v
C--time common block
      common/time/mv(23000,5)
      DOUBLE PRECISION mv
C--factor in front of formation times
  common/ftimefac/ftfac
  DOUBLE PRECISION ftfac
C--colour index common block
  common/colour/trip(23000),anti(23000),colmax
  INTEGER trip,anti,colmax
C--Parameter common block
  common/param/q0,lps,lqcd,ltime,scalefacm,angord,scatrecoil,
     &allhad,compress,nf
      INTEGER nf
  DOUBLE PRECISION q0,lqcd,ltime,lps,scalefacm
      LOGICAL angord,scatrecoil,allhad,compress
C--discard event flag
  common/disc/ndisc,nstrange,ngood,errcount,wdisc,discard
  LOGICAL discard
  INTEGER ndisc,nstrange,ngood,errcount
  double precision wdisc
C--variables for angular ordering
      common/angor/za(23000),zd(23000),thetaa(23000),qqbard(23000)
  DOUBLE PRECISION za,zd,thetaa
      LOGICAL qqbard
C--variables for coherent scattering
  common/coherent/nstart,nend,allqs(10000,6),scatcentres(10000,10),
     &qsumvec(4),qsum2
  INTEGER nstart,nend
  DOUBLE PRECISION allqs,scatcentres,qsumvec,qsum2
C--number of scattering events
  common/check/nscat,nscateff,nsplit
  DOUBLE PRECISION nscat,nscateff,nsplit
C--event weight
  common/weight/evweight,sumofweights
  double precision evweight,sumofweights
C--extra storage for scattering centres before interactions
      common/storescatcen/nscatcen,maxnscatcen,scatflav(10000),
     &scatcen(10000,5),writescatcen,writedummies
  integer nscatcen,maxnscatcen,scatflav
  double precision scatcen
  logical writescatcen,writedummies
C--local variables
      INTEGER l,line,n1,n2,j,dir,lold,nold,colmaxold,statold,nscatcenold
      DOUBLE PRECISION pyr,pi,beta(3),theta,phi,pyp,phi2,maxt,t,
     &newmass,deltam,dm,ttot,dmleft,lambda,time,endtime,x,tmp,
     &m32,newm2,shat,theta2,z,gettemp,e3new,e4new,p32,p42,p3old,
     &newm,mass2,enew,pt2,pt,pl,m12,firsttime,pcms2
      CHARACTER*2 typ
  LOGICAL retrysplit,qqbar,qqbardec,rejectt,redokin,reshuffle
  DATA pi/3.141592653589793d0/

      IF((n+2*(n2-n1+1)).GT.22990)THEN
        write(logfid,*)'event too long for event record'
        discard=.true.
        RETURN
      ENDIF

  firsttime = mv(l,5)

  redokin = .false.

  newm2=newm
  nold=n
  colmaxold=colmax
  statold=k(l,1)
 204  deltam=newm2-p(l,5)
  dmleft=deltam

  ttot=0.d0
  DO 220 j=n1,n2
   ttot=ttot+allqs(j,1)
 220  CONTINUE

  line=l

  DO 222 j=n1,n2
  
C--projectile type
   IF(k(line,2).EQ.21)THEN
    typ='GC'
    IF(pyr(0).LT.0.5)THEN
     dir=1
    ELSE
     dir=-1
    ENDIF
   ELSE
    typ='QQ'
    dir=0
   ENDIF
       k(1,1)=6
   k(1,2)=scatcentres(j,1)
   p(1,1)=scatcentres(j,2)
   p(1,2)=scatcentres(j,3)
   p(1,3)=scatcentres(j,4)
   p(1,4)=scatcentres(j,5)
   p(1,5)=scatcentres(j,6)
       mv(1,1)=scatcentres(j,7)
       mv(1,2)=scatcentres(j,8)
       mv(1,3)=scatcentres(j,9)
       mv(1,4)=scatcentres(j,10)
   t=allqs(j,1)
   if (t.eq.0.d0) then
     rejectt = .true.
   else 
     rejectt = .false.
   endif

C--transform to c.m.s. and rotate such that parton momentum is in z-direction
       beta(1)=(p(1,1)+p(line,1))/(p(1,4)+p(line,4))
       beta(2)=(p(1,2)+p(line,2))/(p(1,4)+p(line,4))
       beta(3)=(p(1,3)+p(line,3))/(p(1,4)+p(line,4))
       IF ((beta(1).GT.1.d0).OR.(beta(2).GT.1.d0).OR.(beta(3).GT.1.d0)
     &  .or.(sqrt(beta(1)**2+beta(2)**2+beta(3)**2).gt.1.d0))THEN
     reshuffle = .false.
   else 
     reshuffle = .true.
   endif
 205   if (.not.reshuffle) then
         beta(1)=p(1,1)/p(1,4)
         beta(2)=p(1,2)/p(1,4)
         beta(3)=p(1,3)/p(1,4)
         CALL pyrobo(line,line,0d0,0d0,-beta(1),-beta(2),-beta(3))
         CALL pyrobo(1,1,0d0,0d0,-beta(1),-beta(2),-beta(3))
         theta=pyp(line,13)
         phi=pyp(line,15)
         CALL pyrobo(line,line,0d0,-phi,0d0,0d0,0d0)
         CALL pyrobo(1,1,0d0,-phi,0d0,0d0,0d0)
         CALL pyrobo(line,line,-theta,0d0,0d0,0d0,0d0)
         CALL pyrobo(1,1,-theta,0d0,0d0,0d0,0d0)

     maxt = -2.*p(1,5)*p(line,4)
     if (t.lt.maxt) then
       t=0.d0
       rejectt = .true.
     endif
     m12 = -p(line,5)**2
 203     enew = p(line,4)+t/(2.*p(1,5))
     pl = (t+2.*p(line,4)*enew-2.*m12)/(2.*p(line,3))
     pt2 = enew**2-pl**2-m12
     if (t.eq.0.d0) pt2 = 0.d0
     if (dabs(pt2).lt.1.d-8) pt2 = 0.d0
     if (pt2.lt.0.d0) then
       write(logfid,*)' This should not have happened: pt^2<0!'
       write(logfid,*)t,enew,pl,pt2
       t = 0.d0
       rejectt = .true.
       goto 203
     endif
     pt = sqrt(pt2)
     phi2 = pyr(0)*2.*pi
     n=n+2
     p(n,1)=pt*cos(phi2)
     p(n,2)=pt*sin(phi2)
     p(n,3)=pl
     p(n,4)=enew
     p(n,5)=p(line,5)
!---------------------------------       
         p(n-1,1)=p(1,1)+p(line,1)-p(n,1)
         p(n-1,2)=p(1,2)+p(line,2)-p(n,2)
         p(n-1,3)=p(1,3)+p(line,3)-p(n,3)
         p(n-1,4)=p(1,4)+p(line,4)-p(n,4)
     mass2 = p(n-1,4)**2-p(n-1,1)**2-p(n-1,2)**2-p(n-1,3)**2
     if ((mass2.lt.0.d0).and.(mass2.gt.-1.-6))  mass2=0.d0
         if (mass2.lt.0.d0)  
     &  write(logfid,*)'messed up scattering centres mass^2: ',
     &  mass2,p(1,5)**2
         p(n-1,5)=sqrt(mass2)
     if (abs(p(n-1,5)-p(1,5)).gt.1.d-6)
     &  write(logfid,*)'messed up scattering centres mass: ',
     &  p(n-1,5),p(1,5),p(l,5)
     call flush(logfid)
!---------------------------------       
!        P(N-1,1)=P(1,1)
!        P(N-1,2)=P(1,2)
!        P(N-1,3)=P(1,3)
!        P(N-1,4)=P(1,4)
!        P(N-1,5)=P(1,5)
!---------------------------------       
   else 
         CALL pyrobo(line,line,0d0,0d0,-beta(1),-beta(2),-beta(3))
         CALL pyrobo(1,1,0d0,0d0,-beta(1),-beta(2),-beta(3))
     if ((p(1,4).lt.0.d0).or.(p(line,4).lt.0.d0)) then
           CALL pyrobo(1,1,0d0,0d0,beta(1),beta(2),beta(3))
           CALL pyrobo(line,line,0d0,0d0,beta(1),beta(2),beta(3))
       reshuffle = .false.
       goto 205
     endif
         theta=pyp(line,13)
         phi=pyp(line,15)
         CALL pyrobo(line,line,0d0,-phi,0d0,0d0,0d0)
         CALL pyrobo(1,1,0d0,-phi,0d0,0d0,0d0)
         CALL pyrobo(line,line,-theta,0d0,0d0,0d0,0d0)
         CALL pyrobo(1,1,-theta,0d0,0d0,0d0,0d0)
     shat = (p(1,4)+p(line,4))**2
     p3old = p(line,3)

     maxt = -4.*p(line,3)**2
     if (t.lt.maxt) then
       t=0.d0
       rejectt = .true.
     endif
     theta2 = acos(1.d0+t/(2.*p(line,3)**2))
     phi2 = pyr(0)*2.*pi
     n=n+2
     p(n,1)=p(line,3)*sin(theta2)*cos(phi2)
     p(n,2)=p(line,3)*sin(theta2)*sin(phi2)
     p(n,3)=p(line,3)*cos(theta2)
     p(n,4)=p(line,4)
     p(n,5)=p(line,5)
!---------------------------------       
         p(n-1,1)=p(1,1)+p(line,1)-p(n,1)
         p(n-1,2)=p(1,2)+p(line,2)-p(n,2)
         p(n-1,3)=p(1,3)+p(line,3)-p(n,3)
         p(n-1,4)=p(1,4)+p(line,4)-p(n,4)
     mass2 = p(n-1,4)**2-p(n-1,1)**2-p(n-1,2)**2-p(n-1,3)**2
     if ((mass2.lt.0.d0).and.(mass2.gt.-1.-6))  mass2=0.d0
         if (mass2.lt.0.d0)  
     &  write(logfid,*)'messed up scattering centres mass^2: ',
     &  mass2,p(1,5)**2
         p(n-1,5)=sqrt(mass2)
     if (abs(p(n-1,5)-p(1,5)).gt.1.d-6)
     &  write(logfid,*)'messed up scattering centres mass: ',
     &  p(n-1,5),p(1,5),p(l,5)
     call flush(logfid)
!---------------------------------       
!        P(N-1,1)=P(1,1)
!        P(N-1,2)=P(1,2)
!        P(N-1,3)=P(1,3)
!        P(N-1,4)=P(1,4)
!        P(N-1,5)=P(1,5)
!---------------------------------       
   endif
C--outgoing projectile
       za(n)=1.d0
   thetaa(n)=-1.d0
       zd(n)=z
       qqbard(n)=qqbar
       k(n,1)=k(line,1)
       k(n,2)=k(line,2)
   k(n,3)=l
   k(n,4)=0
   k(n,5)=0
   IF(allhad.and.(.not.rejectt))THEN
    IF(k(n,2).EQ.21)THEN
     IF(dir.EQ.1)THEN
      trip(n)=colmax+1
      anti(n)=anti(line)
     ELSE
      trip(n)=trip(line)
      anti(n)=colmax+1
     ENDIF
    ELSEIF(k(n,2).GT.0)THEN
     trip(n)=colmax+1 
     anti(n)=0
    ELSE
     trip(n)=0
     anti(n)=colmax+1
    ENDIF
    colmax=colmax+1
   ELSE
    trip(n)=trip(line)
    anti(n)=anti(line)
   ENDIF
C--take care of incoming projectile
       IF(k(line,1).EQ.1)THEN
    k(line,1)=12
       ELSE
        k(line,1)=14
       ENDIF
   k(line,4)=n-1
   k(line,5)=n
C--outgoing scattering centre
       za(n-1)=1.d0
   thetaa(n-1)=-1.d0
       zd(n-1)=-1.d0
       qqbard(n-1)=.false.
C--temporary status code, will be overwritten later
       k(n-1,1)=3
   k(n-1,2)=21
   k(n-1,3)=0
   k(n-1,4)=0
   k(n-1,5)=0
   IF(allhad.and.(.not.rejectt))THEN
    IF((k(n,2).GT.0).AND.(dir.GE.0))THEN
     trip(n-1)=trip(line)
     anti(n-1)=trip(n)
    ELSE
     trip(n-1)=anti(n)
     anti(n-1)=anti(line)
    ENDIF
   ELSE
    trip(n-1)=0
    anti(n-1)=0
   ENDIF

   if (reshuffle.and.(dm.gt.0.d0)) then
C--adjust mass and re-shuffle momenta

     IF(ttot.EQ.0.d0)THEN
      dm=0.d0
     ELSE
      if (dmleft.lt.0.d0) then
        dm=max(dmleft*t/ttot*1.5d0,dmleft)
      else
        dm=min(dmleft*t/ttot*1.5d0,dmleft)
      endif
     ENDIF
     ttot=ttot-allqs(j,1)

     newmass = p(n,5)+dm
     if (newmass.lt.0.d0) then
       m32 = -newmass**2
     else
       m32 = newmass**2
     endif
     e3new = (shat + m32 - p(1,5)**2)/(2.d0*sqrt(shat))
     e4new = (shat - m32 + p(1,5)**2)/(2.d0*sqrt(shat))
     p32 = e3new**2 - m32
     p42 = e4new**2 - p(1,5)**2
     if ((p32.lt.0.d0).or.(p42.lt.0.d0).or.
     &       (e3new.lt.0.d0).or.(e4new.lt.0.d0)) then
       p32 = 0.d0
       p42 = 0.d0
       e4new = p(n-1,5)
       e3new = sqrt(shat) - e4new
       m32 = e3new**2
       if ((e3new.lt.0.d0).or.(e4new.lt.0.d0)) then
         e3new = p(n,4)
         e4new = p(n-1,4)
         p32 = p3old**2
         p42 = p3old**2
       if (p(n,5).lt.0.d0) then
           m32 = -p(n,5)**2
       else
           m32 = p(n,5)**2
       endif 
       endif
     endif
     p(n,1) = sqrt(p32)*p(n,1)/p3old
     p(n,2) = sqrt(p32)*p(n,2)/p3old
     p(n,3) = sqrt(p32)*p(n,3)/p3old
     p(n,4) = e3new
     p(n,5) = sign(sqrt(abs(m32)),newmass)
     tmp = p(n,4)**2-p(n,1)**2-p(n,2)**2-p(n,3)**2
     if (abs(tmp-m32).gt.1.d-6) 
     &  write(logfid,*) 'Oups, messed up projectiles mass:',
     &  tmp,m32,p(n,5)
!---------------------------------       
     p(n-1,1) = sqrt(p42)*p(n-1,1)/p3old
     p(n-1,2) = sqrt(p42)*p(n-1,2)/p3old
     p(n-1,3) = sqrt(p42)*p(n-1,3)/p3old
     p(n-1,4) = e4new
     tmp = p(n-1,4)**2-p(n-1,1)**2-p(n-1,2)**2-p(n-1,3)**2
     &  -p(n-1,5)**2
     if (abs(tmp).gt.1.d-6) 
     &  write(logfid,*) 'Oups, messed up scattering centres mass:',
     &  tmp,p3old,p(n-1,1),p(n-1,2),p(n-1,3),p(n-1,4),p(n-1,5)
     if ((abs(p(n,1)+p(n-1,1)).gt.1.d-6).or.
     &     (abs(p(n,2)+p(n-1,2)).gt.1.d-6).or.
     &     (abs(p(n,3)+p(n-1,3)).gt.1.d-6)) 
     &  write(logfid,*) 'Oups, momentum not conserved', 
     &  p(n,1)+p(n-1,1),p(n,2)+p(n-1,2),p(n,3)+p(n-1,3)
!---------------------------------       
!        P(N-1,1)=P(1,1)
!        P(N-1,2)=P(1,2)
!        P(N-1,3)=P(1,3)
!        P(N-1,4)=P(1,4)
!        P(N-1,5)=P(1,5)
!---------------------------------       
   endif

C--transformation to lab
       CALL pyrobo(n-1,n,theta,0d0,0d0,0d0,0d0)
       CALL pyrobo(line,line,theta,0d0,0d0,0d0,0d0)
       CALL pyrobo(n-1,n,0d0,phi,0d0,0d0,0d0)
       CALL pyrobo(line,line,0d0,phi,0d0,0d0,0d0)
       CALL pyrobo(n-1,n,0d0,0d0,beta(1),beta(2),beta(3))
       CALL pyrobo(line,line,0d0,0d0,beta(1),beta(2),beta(3))
       CALL pyrobo(1,1,theta,0d0,0d0,0d0,0d0)
       CALL pyrobo(1,1,0d0,phi,0d0,0d0,0d0)
       CALL pyrobo(1,1,0d0,0d0,beta(1),beta(2),beta(3))
      if (.not.allhad) then
    k(n-1,1)=13
   else
        IF(scatrecoil.AND.(p(n-1,4).GT.(10.*3.*
     &gettemp(mv(1,1),mv(1,2),mv(1,3),mv(1,4)))))THEN
         k(n-1,1)=2
        ELSE
         k(n-1,1)=3
        ENDIF
   endif
   if (rejectt) k(n-1,1)=11
       mv(n,4)=mv(1,4)
       mv(n-1,4)=mv(1,4)
C--set the production vertices: x_mother + (tprod - tprod_mother) * beta_mother
       mv(n-1,1)=mv(line,1)
     &  +(mv(n-1,4)-mv(line,4))*p(line,1)/max(pyp(line,8),p(line,4))
       mv(n-1,2)=mv(line,2)
     &  +(mv(n-1,4)-mv(line,4))*p(line,2)/max(pyp(line,8),p(line,4))
       mv(n-1,3)=mv(line,3)
     &  +(mv(n-1,4)-mv(line,4))*p(line,3)/max(pyp(line,8),p(line,4))
       mv(n,  1)=mv(line,1)
     &  +(mv(n,  4)-mv(line,4))*p(line,1)/max(pyp(line,8),p(line,4))
       mv(n,  2)=mv(line,2)
     &  +(mv(n,  4)-mv(line,4))*p(line,2)/max(pyp(line,8),p(line,4))
       mv(n,  3)=mv(line,3)
     &  +(mv(n,  4)-mv(line,4))*p(line,3)/max(pyp(line,8),p(line,4))
   IF(p(n-1,5).GT.p(1,5))THEN
     lambda=1.d0/(ftfac*0.2*p(n-1,4)/p(n-1,5)**2)
     mv(n-1,5)=mv(n-1,4)-log(1.d0-pyr(0))/lambda
   ELSE
        mv(n-1,5)=0.d0
   ENDIF
   IF(j.LT.n2)THEN
        mv(n,5)=scatcentres(j+1,10)
   ELSE
    IF(p(n,5).GT.0.d0)THEN
     IF(deltam.EQ.0.d0)THEN
      endtime=firsttime
     ELSE
      IF(x.LT.1.d0)THEN
           lambda=1.d0/(ftfac*p(n,4)*0.2/p(n,5)**2)
       endtime=scatcentres(j,10)-log(1.d0-pyr(0))/lambda
      ELSE
       endtime=time
      ENDIF
     ENDIF
     mv(n,5)=endtime
    ELSE
         mv(n,5)=0.d0
    ENDIF
   ENDIF
   mv(line,5)=allqs(j,6)


C--store scattering centre before interaction in separate common block
   if (writescatcen.and.(.not.rejectt).and.
     &    (nscatcen.lt.maxnscatcen)) then
    nscatcen = nscatcen+1
    if (nscatcen.le.maxnscatcen) then
     scatflav(nscatcen) = k(1,2)
     scatcen(nscatcen,1) = p(1,1)
     scatcen(nscatcen,2) = p(1,2)
     scatcen(nscatcen,3) = p(1,3)
     scatcen(nscatcen,4) = p(1,4)
     scatcen(nscatcen,5) = p(1,5)
    else
     write(logfid,*) 
     &'WARNING: no room left to store further scattering centres'
    endif
   endif

! if ((p(line,4).gt.100.d0).and.(p(n,4)-p(line,4).gt.1.d0)) then
!   write(*,*)p(line,1),p(line,2),p(line,3),p(line,4),p(line,5)
!   write(*,*)p(n,1),p(n,2),p(n,3),p(n,4),p(n,5)
!   write(*,*)p(1,1),p(1,2),p(1,3),p(1,4),p(1,5)
!   write(*,*)p(n-1,1),p(n-1,2),p(n-1,3),p(n-1,4),p(n-1,5)
!   write(*,*)t
!   write(*,*)GETTEMP(MV(1,1),MV(1,2),MV(1,3),MV(1,4))
!   write(*,*)
! endif

   dmleft=dmleft-(p(n,5)-p(line,5))
   line=n
   tmp = abs(p(n,4)**2-p(n,1)**2-p(n,2)**2-p(n,3)**2)-p(n,5)**2
   if (abs(tmp).ge.1.d-6) 
     &  write(logfid,*)tmp,j,p(l,5),p(line,5),p(n,5)
 222  CONTINUE
  if (p(n,5).lt.0.d0) then
    retrysplit=.true.
    return
  endif
  if (p(n,5).ne.newm2) then
    retrysplit=.true.
    redokin = .true.
    n=nold
    colmax=colmaxold
    k(l,1)=statold
    if (p(l,5).le.0.d0) then
      newm2 = 0.d0
    else
          if (p(l,5).lt.q0) then
            if ((newm2.eq.newm).and.(newm.ne.q0+1.d-6)) then
              newm2=q0+1.d-6
            else    
              retrysplit=.true.
              return
            endif
          else
            newm2=p(l,5)
          endif
          n2=n1
        endif
    goto 204
  endif
  if ((k(n,1).eq.1).and.
     &  ((p(n,5).lt.0.d0).or.((p(n,5).gt.0.d0).and.(p(n,5).lt.q0))))
     &write(logfid,*)'dokinematics did not reach sensible mass: ',
     &p(n,5),newm,p(l,5),newm2
  nscateff=nscateff+evweight
      END



***********************************************************************
***   function getproba
***********************************************************************
  DOUBLE PRECISION FUNCTION getproba(QI,QF,QAA,ZAA,EBB,TYPE,
     &  t1,ins2)
  IMPLICIT NONE
C--variables for Sudakov integration
  common/sudaint/qa,za2,eb,t,instate,typ
  DOUBLE PRECISION qa,za2,eb,t
  CHARACTER*2 typ
  LOGICAL instate
C--local variables
  DOUBLE PRECISION qi,qf,qaa,zaa,ebb,getsudakov,deriv,t1
  CHARACTER*2 type
  LOGICAL ins2

  qa=qaa
  za2=zaa
  eb=ebb
  typ=TYPE
  t=t1
  instate=ins2
  getproba=getsudakov(qi,qaa,qf,zaa,ebb,type,t1,ins2)
     &      *deriv(qf,1)
  END


***********************************************************************
***   function getsudakov
***********************************************************************
  DOUBLE PRECISION FUNCTION getsudakov(QMAX1,QA1,QB1,ZA1,EB1,
     &                                                type3,t2,ins)
  IMPLICIT NONE
C--identifier of file for hepmc output and logfile
  common/hepmcid/hpmcfid,logfid
  integer hpmcfid,logfid
C--Parameter common block
  common/param/q0,lps,lqcd,ltime,scalefacm,angord,scatrecoil,
     &allhad,compress,nf
      INTEGER nf
  DOUBLE PRECISION q0,lqcd,ltime,lps,scalefacm
      LOGICAL angord,scatrecoil,allhad,compress
C--variables for Sudakov integration
  common/sudaint/qa,za2,eb,t,instate,typ
  DOUBLE PRECISION qa,za2,eb,t
  CHARACTER*2 typ
  LOGICAL instate
C--local variables
  DOUBLE PRECISION qmax1,qa1,qb1,za1,eb1,tmax,tb,ystart,epsi,
     &hfirst,t2,getinsudafast,qb2
  CHARACTER*2 type3
  LOGICAL ins
      DATA epsi/1.d-4/

  qb2=qb1
  IF(ins)THEN
       IF(qb2.LT.q0) write(logfid,*) 'error: Q < Q0',qb2,qmax1
       IF(qb2.LT.(q0+1.d-10)) qb2=qb2+1.d-10
      ELSE 
       IF(qb2.LT.q0) write(logfid,*) 'error: Q < min',qb2,qmax1
       IF(qb2.LT.(q0+1.d-10)) qb2=qb2+1.d-10
      ENDIF 
      IF(qb2.GE.(qmax1-1.d-10)) THEN
       getsudakov=1.d0
      ELSE
   IF(ins)THEN
    getsudakov=getinsudafast(qb1,qmax1,type3)
   ELSE
    qa=qa1
    za2=za1
    eb=eb1
    typ=type3
    t=t2
    instate=.false.
        hfirst=0.01*(qmax1-qb1)
        ystart=0.d0
        CALL odeint(ystart,qb2,qmax1,epsi,hfirst,0.d0,1)
        getsudakov=exp(-ystart)
   ENDIF
      ENDIF
  END


***********************************************************************
***   function getinsudakov
***********************************************************************
  DOUBLE PRECISION FUNCTION getinsudakov(QB,QMAX1,TYPE3)
  IMPLICIT NONE
C--identifier of file for hepmc output and logfile
  common/hepmcid/hpmcfid,logfid
  integer hpmcfid,logfid
C--Parameter common block
  common/param/q0,lps,lqcd,ltime,scalefacm,angord,scatrecoil,
     &allhad,compress,nf
      INTEGER nf
  DOUBLE PRECISION q0,lqcd,ltime,lps,scalefacm
      LOGICAL angord,scatrecoil,allhad,compress
C--variables for Sudakov integration
  common/sudaint/qa,za2,eb,t,instate,typ
  DOUBLE PRECISION qa,za2,eb,t
  CHARACTER*2 typ
  LOGICAL instate
C--local variables
  DOUBLE PRECISION qmax1,qb,qb1,za1,ea1,ystart,epsi,
     &hfirst
  CHARACTER*2 type3
      DATA epsi/1.d-4/

      qb1=qb
      IF(qb1.LT.q0) write(logfid,*) 'error: Q < Q0',qb1,qmax1
      IF(qb1.LT.(q0+1.d-12)) qb1=qb1+1.d-12
      IF(qb1.GE.(qmax1-1.d-12)) THEN
       getinsudakov=1.d0
      ELSE
   typ=type3
       hfirst=0.01*(qmax1-qb1)
       ystart=0.d0
       CALL odeint(ystart,qb1,qmax1,epsi,hfirst,0.d0,6)
       getinsudakov=exp(-ystart)
      ENDIF
  END


***********************************************************************
***   function deriv
***********************************************************************
      DOUBLE PRECISION FUNCTION deriv(XVAL,W4)
      IMPLICIT NONE
C--Parameter common block
  common/param/q0,lps,lqcd,ltime,scalefacm,angord,scatrecoil,
     &allhad,compress,nf
      INTEGER nf
  DOUBLE PRECISION q0,lqcd,ltime,lps,scalefacm
      LOGICAL angord,scatrecoil,allhad,compress
C--variables for splitting function integration
  common/intsplitf/qquad,fm
  DOUBLE PRECISION qquad,fm
C--variables for Sudakov integration
  common/sudaint/qa,za2,eb,t,instate,typ
  DOUBLE PRECISION qa,za2,eb,t
  CHARACTER*2 typ
  LOGICAL instate
C--variables for pdf integration
  common/pdfintv/xmax,z
  DOUBLE PRECISION xmax,z
C--variables for cross section integration 
  common/xsecv/qlow,mdx
  DOUBLE PRECISION qlow,mdx
C--local variables
  INTEGER w4
      DOUBLE PRECISION xval,getspliti,pi,alphas,getinspliti,
     &getinsudafast,scatprimfunc,pqq,pqg,pgg,pgq,
     &medderiv
  DATA pi/3.141592653589793d0/

  IF(w4.EQ.1)THEN
C--Sudakov integration
   IF(instate)THEN
        deriv=2.*getinspliti(xval,typ)/xval
   ELSE
        deriv=2.*getspliti(qa,xval,za2,eb,typ)/xval
   ENDIF
  ELSEIF(w4.EQ.2)THEN
C--P(q->qg) integration
   deriv=(1.+fm)*alphas(xval*(1.-xval)*qquad/1.,lps)*
     &    pqq(xval)/(2.*pi)
  ELSEIF(w4.EQ.3)THEN
C--P(g->gg) integration
       deriv=(1.+fm)*alphas(xval*(1.-xval)*qquad/1.,lps)
     &           *pgg(xval)/(2.*pi)
  ELSEIF(w4.EQ.4)THEN
C--P(g->qq) integration
   deriv=(1.+fm)*alphas(xval*(1-xval)*qquad/1.,lps)*
     &  pqg(xval)/(2.*pi) 
  ELSEIF(w4.EQ.5)THEN
   deriv=exp(-xval)/xval
  ELSEIF(w4.EQ.6)THEN
       deriv=2.*getinspliti(xval,typ)/xval
  ELSEIF(w4.EQ.7)THEN
   deriv=2.*getinsudafast(xval,xmax,'QQ')
     &  *alphas((1.-z)*xval**2/1.,lps)
     &  *pqq(z)/(2.*pi*xval)
  ELSEIF(w4.EQ.8)THEN
   deriv=2.*getinsudafast(xval,xmax,'GC')
     &  *alphas((1.-z)*xval**2/1.,lps)
     &  *pgq(z)/(2.*pi*xval)
  ELSEIF(w4.EQ.9)THEN
   deriv=2.*getinsudafast(xval,xmax,'QQ')
     &  *alphas((1.-z)*xval**2/1.,lps)
     &  *pqg(z)/(2.*pi*xval)  
  ELSEIF(w4.EQ.10)THEN
   deriv=2.*getinsudafast(xval,xmax,'GC')
     &  *alphas((1.-z)*xval**2/1.,lps)*
     &      *2.*pgg(z)/(2.*pi*xval)
  ELSEIF(w4.EQ.11)THEN
   deriv=3.*getinspliti(scalefacm*sqrt(xval),'GQ')
     &  *scatprimfunc(xval,mdx)/(2.*xval)
  ELSEIF(w4.EQ.12)THEN
   deriv=2.*getinspliti(scalefacm*sqrt(xval),'QG')
     &  *scatprimfunc(xval,mdx)/(3.*xval)
  ELSEIF(w4.EQ.13)THEN
   deriv=getinsudafast(qlow,scalefacm*sqrt(xval),'GC')
     &  *3.*2.*pi*alphas(xval+mdx**2,lqcd)**2/(2.*(xval+mdx**2)**2)
  ELSEIF(w4.EQ.14)THEN
   deriv=getinsudafast(qlow,scalefacm*sqrt(xval),'QQ')
     &  *2.*2.*pi*alphas(xval+mdx**2,lqcd)**2/(3.*(xval+mdx**2)**2)
  ELSEIF(w4.EQ.21)THEN
   deriv=2.*getinsudafast(xval,xmax,'QQ')*getinspliti(xval,'QQ')
     &  /xval
  ELSEIF(w4.EQ.22)THEN
   deriv=2.*getinsudafast(xval,xmax,'GC')*getinspliti(xval,'GQ')
     &  /xval
  ELSEIF(w4.EQ.23)THEN
   deriv=2.*getinsudafast(xval,xmax,'QQ')*getinspliti(xval,'QG')
     &  /xval
  ELSEIF(w4.EQ.24)THEN
   deriv=2.*getinsudafast(xval,xmax,'GC')*2.
     &  *getinspliti(xval,'GG')/xval
      ELSE
       deriv=medderiv(xval,w4-100)
      ENDIF
      END


***********************************************************************
***   function getspliti
***********************************************************************
  DOUBLE PRECISION FUNCTION getspliti(QA,QB,ZETA,EB,TYPE1)
  IMPLICIT NONE
C--identifier of file for hepmc output and logfile
  common/hepmcid/hpmcfid,logfid
  integer hpmcfid,logfid
C--Parameter common block
  common/param/q0,lps,lqcd,ltime,scalefacm,angord,scatrecoil,
     &allhad,compress,nf
      INTEGER nf
  DOUBLE PRECISION q0,lqcd,ltime,lps,scalefacm
      LOGICAL angord,scatrecoil,allhad,compress
C--splitting integral
      common/splitint/splitiggv(1000,1000),splitiqqv(1000,1000),
     &splitiqgv(1000,1000),qval(1000),zmval(1000),qmax,zmmin,npoint
      INTEGER npoint
      DOUBLE PRECISION splitiggv,splitiqqv,splitiqgv,
     &qval,zmval,qmax,zmmin
C--variables for splitting function integration
  common/intsplitf/qquad,fm
  DOUBLE PRECISION qquad,fm
C--number of extrapolations in tables
  common/extrapolations/ntotspliti,noverspliti,ntotpdf,noverpdf,
     &ntotxsec,noverxsec,ntotsuda,noversuda
  integer ntotspliti,noverspliti,ntotpdf,noverpdf,
     &ntotxsec,noverxsec,ntotsuda,noversuda
C--local variables
  INTEGER i,j,lt,qlmax,zlmax,qline,zline
  DOUBLE PRECISION qa,qb,zeta,eb,low,x1a(2),x2a(2),ya(2,2),y,
     &splitintgg,splitintqg,a,b,yb(2)
  CHARACTER*2 type1 

  ntotspliti=ntotspliti+1
  if (qb.gt.qmax) then
    noverspliti=noverspliti+1
    if (noverspliti.le.25) 
     &  write(logfid,*)'WARNING in getspliti: need to extrapolate: ',
     &  qb,qmax
  endif

C--find boundaries for z integration
      IF(angord.AND.(zeta.NE.1.d0))THEN
       low=max(0.5-0.5*sqrt(1.-q0**2/qb**2)
     &  *sqrt(1.-qb**2/eb**2),
     &     0.5-0.5*sqrt(1.-4.*qb**2*(1.-zeta)/(zeta*qa**2)))
      ELSE
       low=0.5-0.5*sqrt(1.-q0**2/qb**2)
     &  *sqrt(1.-qb**2/eb**2)
      ENDIF
C--find values in array
        qlmax=int((qb-qval(1))*npoint/(qval(1000)-qval(1))+1)
        qline=max(qlmax,1)
        qline=min(qline,npoint)
        zlmax=int((log(low)-log(zmval(1)))*npoint/
     &        (log(zmval(1000))-log(zmval(1)))+1)
        zline=max(zlmax,1)
        zline=min(zline,npoint)
    IF((qline.GT.999).OR.(zline.GT.999).OR.
     &  (qline.LT.1).OR.(zline.LT.1))THEN 
         write(logfid,*)'ERROR in GETSPLITI: line number out of bound',
     &  qline,zline
    ENDIF
        IF((type1.EQ.'GG').OR.(type1.EQ.'GC'))THEN
         DO 17 i=1,2
          x1a(i)=qval(qline-1+i)
          x2a(i)=zmval(zline-1+i)
          DO 16 j=1,2
           ya(i,j)=splitiggv(qline-1+i,zline-1+j)
 16       CONTINUE
 17      CONTINUE
     DO 30 i=1,2
      a=(ya(i,2)-ya(i,1))/(x2a(2)-x2a(1))
      b=ya(i,1)-a*x2a(1)
      yb(i)=a*low+b
 30    CONTINUE
     IF(x1a(1).EQ.x1a(2))THEN
      y=(yb(1)+yb(2))/2.
     ELSE
      a=(yb(2)-yb(1))/(x1a(2)-x1a(1))
      b=yb(1)-a*x1a(1)
      y=a*qb+b
     ENDIF
         IF(type1.EQ.'GG')THEN
          getspliti=min(y,10.d0)
         ELSE
          splitintgg=min(y,10.d0)
         ENDIF
        ENDIF
        IF((type1.EQ.'QG').OR.(type1.EQ.'GC'))THEN
         DO 19 i=1,2
          x1a(i)=qval(qline-1+i)
          x2a(i)=zmval(zline-1+i)
          DO 18 j=1,2
           ya(i,j)=splitiqgv(qline-1+i,zline-1+j)
 18       CONTINUE
 19      CONTINUE
     DO 31 i=1,2
      a=(ya(i,2)-ya(i,1))/(x2a(2)-x2a(1))
      b=ya(i,1)-a*x2a(1)
      yb(i)=a*low+b
 31    CONTINUE
     IF(x1a(1).EQ.x1a(2))THEN
      y=(yb(1)+yb(2))/2.
     ELSE
      a=(yb(2)-yb(1))/(x1a(2)-x1a(1))
      b=yb(1)-a*x1a(1)
      y=a*qb+b
     ENDIF
         IF(type1.EQ.'QG')THEN
          getspliti=nf*min(y,10.d0)
         ELSE
          splitintqg=nf*min(y,10.d0)
         ENDIF
        ENDIF
        IF(type1.EQ.'QQ')THEN
         DO 21 i=1,2
          x1a(i)=qval(qline-1+i)
          x2a(i)=zmval(zline-1+i)
          DO 20 j=1,2
           ya(i,j)=splitiqqv(qline-1+i,zline-1+j)
 20       CONTINUE
 21      CONTINUE
     DO 32 i=1,2
      a=(ya(i,2)-ya(i,1))/(x2a(2)-x2a(1))
      b=ya(i,1)-a*x2a(1)
      yb(i)=a*low+b
 32    CONTINUE
     IF(x1a(1).EQ.x1a(2))THEN
      y=(yb(1)+yb(2))/2.
     ELSE
      a=(yb(2)-yb(1))/(x1a(2)-x1a(1))
      b=yb(1)-a*x1a(1)
      y=a*qb+b
     ENDIF
         getspliti=min(y,10.d0)
        ENDIF
        IF(type1.EQ.'GC') getspliti=splitintgg+splitintqg
      END


***********************************************************************
***   function getinspliti
***********************************************************************
  DOUBLE PRECISION FUNCTION getinspliti(QB,TYPE1)
  IMPLICIT NONE
C--Parameter common block
  common/param/q0,lps,lqcd,ltime,scalefacm,angord,scatrecoil,
     &allhad,compress,nf
      INTEGER nf
  DOUBLE PRECISION q0,lqcd,ltime,lps,scalefacm
      LOGICAL angord,scatrecoil,allhad,compress
C--local variables
  DOUBLE PRECISION qb,low,pi,y,splitintgg,splitintqg,up,ei
  CHARACTER*2 type1 
  DATA pi/3.141592653589793d0/

C--find boundaries for z integration
   up = 1. - q0**2/(4.*qb**2)
       IF((type1.EQ.'GG').OR.(type1.EQ.'GC'))THEN
    low=1.d0-up
    IF (up.LE.low) THEN
     getinspliti=0.d0
     RETURN
    ENDIF
    y = 2.* ( log(log((1.-low)*qb**2/lps**2))
     &  - lps**2*ei(log((1.-low)*qb**2/lps**2))/qb**2
     &  + lps**4*ei(2.*log((1.-low)*qb**2/lps**2))/qb**4
     &  - lps**6*ei(3.*log((1.-low)*qb**2/lps**2))/qb**6
     &      - log(log((1.-up)*qb**2/lps**2))
     &  + lps**2*ei(log((1.-up)*qb**2/lps**2))/qb**2
     &  - lps**4*ei(2.*log((1.-up)*qb**2/lps**2))/qb**4
     &  + lps**6*ei(3.*log((1.-up)*qb**2/lps**2))/qb**6
     &  + low - log(low) - up + log(up) )
     &  *3.*12.*pi/(2.*pi*(33.-2.*nf))
        IF(type1.EQ.'GG')THEN
         getinspliti=y
        ELSE
         splitintgg=y
        ENDIF
       ENDIF
       IF((type1.EQ.'QG').OR.(type1.EQ.'GC'))THEN
    low=0.d0
    IF (up.LE.low) THEN
     getinspliti=0.d0
     RETURN
    ENDIF
    y = ( 2.*lps**6*ei(3.*log((1.-low)*qb**2/lps**2))/qb**6
     &  - 2.*lps**4*ei(2.*log((1.-low)*qb**2/lps**2))/qb**4
     &  + 2.*lps**2*ei(log((1.-low)*qb**2/lps**2))/qb**2
     &  - 2.*lps**6*ei(3.*log((1.-up)*qb**2/lps**2))/qb**6
     &  + 2.*lps**4*ei(2.*log((1.-up)*qb**2/lps**2))/qb**4
     &  - 2.*lps**2*ei(log((1.-up)*qb**2/lps**2))/qb**2 )
     &  *12.*pi/(2.*2.*pi*(33.-2.*nf))
        IF(type1.EQ.'QG')THEN
         getinspliti=nf*y
        ELSE
         splitintqg=nf*y
        ENDIF
       ENDIF
       IF(type1.EQ.'QQ')THEN
    low=0.d0
    IF (up.LE.low) THEN
     getinspliti=0.d0
     RETURN
    ENDIF
    y = ( 2.*log(log((1.-low)*qb**2/lps**2))
     &  - 2.*lps**2*ei(log((1.-low)*qb**2/lps**2))/qb**2
     &  + lps**4*ei(2.*log((1.-low)*qb**2/lps**2))/qb**4
     &  - 2.*log(log((1.-up)*qb**2/lps**2))
     &  + 2.*lps**2*ei(log((1.-up)*qb**2/lps**2))/qb**2
     &  - lps**4*ei(2.*log((1.-up)*qb**2/lps**2))/qb**4 ) 
     &  *4.*12.*pi/(3.*2.*pi*(33.-2.*nf))
        getinspliti=y
       ENDIF
       IF(type1.EQ.'GQ')THEN
    low=1.d0-up
    IF (up.LE.low) THEN
     getinspliti=0.d0
     RETURN
    ENDIF
    y = (up**2/2.-2.*up+2.*log(up)-low**2/2.+2.*low- 2.*log(low)) 
     &  *4.*12.*pi/(3.*2.*pi*(33.-2.*nf)*log(qb**2/lps**2))
        getinspliti=y
       ENDIF
       IF(type1.EQ.'GC') getinspliti=splitintgg+splitintqg
      END


***********************************************************************
***   function getpdf
***********************************************************************
  DOUBLE PRECISION FUNCTION getpdf(X,Q,TYP)
  IMPLICIT NONE
C--identifier of file for hepmc output and logfile
  common/hepmcid/hpmcfid,logfid
  integer hpmcfid,logfid
C--pdf common block
  common/pdfs/qinqx(2,1000),ginqx(2,1000),qingx(2,1000),
     &gingx(2,1000)
  DOUBLE PRECISION qinqx,ginqx,qingx,gingx
C--Parameter common block
  common/param/q0,lps,lqcd,ltime,scalefacm,angord,scatrecoil,
     &allhad,compress,nf
      INTEGER nf
  DOUBLE PRECISION q0,lqcd,ltime,lps,scalefacm
      LOGICAL angord,scatrecoil,allhad,compress
C--variables for pdf integration
  common/pdfintv/xmax,z
  DOUBLE PRECISION xmax,z
C--local variables
  DOUBLE PRECISION x,q,qlow,qhigh,ystart,epsi,hfirst
  CHARACTER*2 typ
  DATA epsi/1.d-4/  

  IF((x.LT.0.d0).OR.(x.GT.1.d0).OR.(q.LT.q0))THEN
   write(logfid,*)'error in GETPDF: parameter out of bound',x,q
   getpdf=0.d0
   RETURN
  ENDIF

  IF(typ.EQ.'QQ')THEN
    z=x
    xmax=q
C--f_q^q
    qlow=max(q0,q0/(2.*sqrt(1.-x)))
    qhigh=q
    IF((qlow.GE.qhigh*(1.d0-1.d-10)).OR.(x.GT.1.d0-1.d-10))THEN
     ystart=0.d0
    ELSE
         hfirst=0.01*(qhigh-qlow)
         ystart=0.d0
         CALL odeint(ystart,qlow,qhigh,epsi,hfirst,0.d0,7)
    ENDIF
    getpdf=ystart
  ELSEIF(typ.EQ.'GQ')THEN
    z=x
    xmax=q
C--f_q^g
    qlow=max(q0,max(q0/(2.*sqrt(x)),q0/(2.*sqrt(1.-x))))
    qhigh=q
    IF((qlow.GE.qhigh*(1.d0-1.d-10)).OR.(x.LT.0.d0+1.d-10)
     &  .OR.(x.GT.1.d0-1.d-10))THEN
     ystart=0.d0
    ELSE
         hfirst=0.01*(qhigh-qlow)
         ystart=0.d0
         CALL odeint(ystart,qlow,qhigh,epsi,hfirst,0.d0,8)
    ENDIF
    getpdf=ystart
  ELSEIF(typ.EQ.'QG')THEN
    z=x
    xmax=q
C--f_q^g
    qlow=max(q0,q0/(2.*sqrt(1.-x)))
    qhigh=q
    IF((qlow.GE.qhigh*(1.d0-1.d-10)).OR.(x.GT.1.d0-1.d-10))THEN
     ystart=0.d0
    ELSE
         hfirst=0.01*(qhigh-qlow)
         ystart=0.d0
         CALL odeint(ystart,qlow,qhigh,epsi,hfirst,0.d0,9)
    ENDIF
    getpdf=ystart
  ELSEIF(typ.EQ.'GG')THEN
    z=x
    xmax=q
C--f_q^q
  qlow=max(q0,max(q0/(2.*sqrt(x)),q0/(2.*sqrt(1.-x))))
    qhigh=q
    IF((qlow.GE.qhigh*(1.d0-1.d-10)).OR.(x.LT.0.d0+1.d-10)
     &  .OR.(x.GT.1.d0-1d-10))THEN
     ystart=0.d0
    ELSE
         hfirst=0.01*(qhigh-qlow)
         ystart=0.d0
         CALL odeint(ystart,qlow,qhigh,epsi,hfirst,0.d0,10)
    ENDIF
    getpdf=ystart
  ELSE
   write(logfid,*)'error: pdf-type ',typ,' does not exist'
   getpdf=0.d0
  ENDIF
  END

***********************************************************************
***   function getpdfxint
***********************************************************************
  DOUBLE PRECISION FUNCTION getpdfxint(Q,TYP)
  IMPLICIT NONE
C--identifier of file for hepmc output and logfile
  common/hepmcid/hpmcfid,logfid
  integer hpmcfid,logfid
C--pdf common block
  common/pdfs/qinqx(2,1000),ginqx(2,1000),qingx(2,1000),
     &gingx(2,1000)
  DOUBLE PRECISION qinqx,ginqx,qingx,gingx
C--number of extrapolations in tables
  common/extrapolations/ntotspliti,noverspliti,ntotpdf,noverpdf,
     &ntotxsec,noverxsec,ntotsuda,noversuda
  integer ntotspliti,noverspliti,ntotpdf,noverpdf,
     &ntotxsec,noverxsec,ntotsuda,noversuda
C--local variables
  INTEGER j,q2close,q2line
  DOUBLE PRECISION q,xa(2),ya(2),y,a,b
  CHARACTER*2 typ

  ntotpdf=ntotpdf+1
  if (q**2.gt.qinqx(1,1000)) then
    noverpdf=noverpdf+1
    if (noverpdf.le.25) 
     &  write(logfid,*)'WARNING in getpdfxint: need to extrapolate: ',
     &  q**2,qinqx(1,1000)
  endif

      q2close=int((log(q**2)-log(qinqx(1,1)))*999.d0/
     &  (log(qinqx(1,1000))-log(qinqx(1,1)))+1)
      q2line=max(q2close,1)
      q2line=min(q2line,999)
  IF((q2line.GT.999).OR.(q2line.LT.1))THEN
       write(logfid,*)'ERROR in GETPDFXINT: line number out of bound',
     &  q2line
  ENDIF

      IF(typ.EQ.'QQ')THEN
       DO 11 j=1,2
        xa(j)=qinqx(1,q2line-1+j)
        ya(j)=qinqx(2,q2line-1+j)
 11    CONTINUE
      ELSEIF(typ.EQ.'GQ')THEN
       DO 13 j=1,2
        xa(j)=ginqx(1,q2line-1+j)
        ya(j)=ginqx(2,q2line-1+j)
 13    CONTINUE
      ELSEIF(typ.EQ.'QG')THEN
       DO 15 j=1,2
        xa(j)=qingx(1,q2line-1+j)
        ya(j)=qingx(2,q2line-1+j)
 15    CONTINUE
      ELSEIF(typ.EQ.'GG')THEN
       DO 17 j=1,2
        xa(j)=gingx(1,q2line-1+j)
        ya(j)=gingx(2,q2line-1+j)
 17    CONTINUE
  ELSE
   write(logfid,*)'error in GETPDFXINT: unknown integral type ',typ
  ENDIF
  a=(ya(2)-ya(1))/(xa(2)-xa(1))
  b=ya(1)-a*xa(1)
  y=a*q**2+b
  getpdfxint=y
  END


***********************************************************************
***   subroutine getpdfxintexact
***********************************************************************
  DOUBLE PRECISION FUNCTION getpdfxintexact(Q,TYP)
  IMPLICIT NONE
C--Parameter common block
  common/param/q0,lps,lqcd,ltime,scalefacm,angord,scatrecoil,
     &allhad,compress,nf
      INTEGER nf
  DOUBLE PRECISION q0,lqcd,ltime,lps,scalefacm
      LOGICAL angord,scatrecoil,allhad,compress
C--variables for pdf integration
  common/pdfintv/xmax,z
  DOUBLE PRECISION xmax,z
C--local variables
  DOUBLE PRECISION q,epsi,ystart,hfirst
  CHARACTER*2 typ
  DATA epsi/1.d-4/
  
      hfirst=0.01d0
      ystart=0.d0
  xmax=q
  z=0.d0
  IF(typ.EQ.'QQ')THEN
       CALL odeint(ystart,q0,q,epsi,hfirst,0.d0,21)
  ELSEIF(typ.EQ.'QG')THEN
       CALL odeint(ystart,q0,q,epsi,hfirst,0.d0,23)
  ELSEIF(typ.EQ.'GQ')THEN
       CALL odeint(ystart,q0,q,epsi,hfirst,0.d0,22)
  ELSEIF(typ.EQ.'GG')THEN
       CALL odeint(ystart,q0,q,epsi,hfirst,0.d0,24)
  ENDIF
  getpdfxintexact=ystart 
  END


***********************************************************************
***   function getxsecint
***********************************************************************
  DOUBLE PRECISION FUNCTION getxsecint(TM,MD,TYP2)
  IMPLICIT NONE
C--identifier of file for hepmc output and logfile
  common/hepmcid/hpmcfid,logfid
  integer hpmcfid,logfid
C--Parameter common block
  common/param/q0,lps,lqcd,ltime,scalefacm,angord,scatrecoil,
     &allhad,compress,nf
      INTEGER nf
  DOUBLE PRECISION q0,lqcd,ltime,lps,scalefacm
      LOGICAL angord,scatrecoil,allhad,compress
C--cross secttion common block
  common/xsecs/intq1(1001,101),intq2(1001,101),
     &intg1(1001,101),intg2(1001,101)
  DOUBLE PRECISION intq1,intq2,intg1,intg2
C--variables for cross section integration 
  common/xsecv/qlow,mdx
  DOUBLE PRECISION qlow,mdx
C--number of extrapolations in tables
  common/extrapolations/ntotspliti,noverspliti,ntotpdf,noverpdf,
     &ntotxsec,noverxsec,ntotsuda,noversuda
  integer ntotspliti,noverspliti,ntotpdf,noverpdf,
     &ntotxsec,noverxsec,ntotsuda,noversuda
C--local variables
  INTEGER tline,tclose,mdclose,mdline,i,j
  DOUBLE PRECISION tm,x1a(2),x2a(2),ya(2,2),y,md,yb(2),a,b
  CHARACTER*2 typ2

  ntotxsec=ntotxsec+1
  if (tm.gt.intq1(1000,101)) then
    noverxsec=noverxsec+1
    if (noverpdf.le.25) 
     &  write(logfid,*)'WARNING in getxsecint: need to extrapolate: ',
     &  tm,intq1(1000,101)
  endif

       tclose=int((log(tm)-log(intq1(1,101)))*999.d0/
     &  (log(intq1(1000,101))-log(intq1(1,101)))+1)
       tline=max(tclose,1)
       tline=min(tline,999)
       mdclose=int((md-intq1(1001,1))*99.d0/
     &(intq1(1001,100)-intq1(1001,1))+1)
       mdline=max(mdclose,1)
       mdline=min(mdline,99)
   IF((tline.GT.999).OR.(mdline.GT.99)
     &  .OR.(tline.LT.1).OR.(mdline.LT.1)) THEN
      write(logfid,*)'ERROR in GETXSECINT: line number out of bound',
     &  tline,mdline
   ENDIF

       IF(typ2.EQ.'QA')THEN
C--first quark integral
        DO 12 i=1,2
         x1a(i)=intq1(1001,mdline-1+i)
         x2a(i)=intq1(tline-1+i,101)
         DO 11 j=1,2
          ya(i,j)=intq1(tline-1+j,mdline-1+i)
 11      CONTINUE
 12     CONTINUE
   ELSEIF(typ2.EQ.'QB')THEN
C--second quark integral
        DO 18 i=1,2
         x1a(i)=intq2(1001,mdline-1+i)
         x2a(i)=intq2(tline-1+i,101)
         DO 17 j=1,2
          ya(i,j)=intq2(tline-1+j,mdline-1+i)
 17      CONTINUE
 18     CONTINUE
   ELSEIF(typ2.EQ.'GA')THEN
C--first gluon integral
        DO 14 i=1,2
         x1a(i)=intg1(1001,mdline-1+i)
         x2a(i)=intg1(tline-1+i,101)
         DO 13 j=1,2
          ya(i,j)=intg1(tline-1+j,mdline-1+i)
 13      CONTINUE
 14     CONTINUE
   ELSEIF(typ2.EQ.'GB')THEN
C--second gluon integral
        DO 16 i=1,2
         x1a(i)=intg2(1001,mdline-1+i)
         x2a(i)=intg2(tline-1+i,101)
         DO 15 j=1,2
          ya(i,j)=intg2(tline-1+j,mdline-1+i)
 15      CONTINUE
 16     CONTINUE
   ELSE
    write(logfid,*)'error in GETXSECINT: unknown integral type ',
     &                    typ2
   ENDIF
   DO 19 i=1,2
    a=(ya(i,2)-ya(i,1))/(x2a(2)-x2a(1))
    b=ya(i,1)-a*x2a(1)
    yb(i)=a*tm+b
 19  CONTINUE
   IF(x1a(1).EQ.x1a(2))THEN
    y=yb(1)
   ELSE
    a=(yb(2)-yb(1))/(x1a(2)-x1a(1))
    b=yb(1)-a*x1a(1)
    y=a*md+b
   ENDIF
   getxsecint=y
  END


***********************************************************************
***   function getinsudafast
***********************************************************************
  DOUBLE PRECISION FUNCTION getinsudafast(Q1,Q2,TYP)
  IMPLICIT NONE
C--identifier of file for hepmc output and logfile
  common/hepmcid/hpmcfid,logfid
  integer hpmcfid,logfid
C--Parameter common block
  common/param/q0,lps,lqcd,ltime,scalefacm,angord,scatrecoil,
     &allhad,compress,nf
      INTEGER nf
  DOUBLE PRECISION q0,lqcd,ltime,lps,scalefacm
      LOGICAL angord,scatrecoil,allhad,compress
C--local variables
  DOUBLE PRECISION q1,q2,getinsudared
  CHARACTER*2 typ
  
  IF(q2.LE.q1)THEN
   getinsudafast=1.d0
  ELSEIF(q1.LE.q0)THEN
   getinsudafast=getinsudared(q2,typ)
  ELSE
   getinsudafast=getinsudared(q2,typ)/getinsudared(q1,typ)
  ENDIF
      IF(getinsudafast.GT.1.d0) getinsudafast=1.d0
  IF(getinsudafast.LT.(-1.d-10))THEN
   write(logfid,*)'ERROR: GETINSUDAFAST < 0:',
     &  getinsudafast,' for',q1,' ',q2,' ',typ
  ENDIF
  if (getinsudafast.lt.0.d0) getinsudafast = 0.d0
  END


***********************************************************************
***   function getinsudared
***********************************************************************
  DOUBLE PRECISION FUNCTION getinsudared(Q,TYP2)
  IMPLICIT NONE
C--identifier of file for hepmc output and logfile
  common/hepmcid/hpmcfid,logfid
  integer hpmcfid,logfid
C--Parameter common block
  common/param/q0,lps,lqcd,ltime,scalefacm,angord,scatrecoil,
     &allhad,compress,nf
      INTEGER nf
  DOUBLE PRECISION q0,lqcd,ltime,lps,scalefacm
      LOGICAL angord,scatrecoil,allhad,compress
C--Sudakov common block
  common/insuda/sudaqq(1000,2),sudaqg(1000,2),sudagg(1000,2),
     &sudagc(1000,2)
  DOUBLE PRECISION sudaqq,sudaqg,sudagg,sudagc
C--number of extrapolations in tables
  common/extrapolations/ntotspliti,noverspliti,ntotpdf,noverpdf,
     &ntotxsec,noverxsec,ntotsuda,noversuda
  integer ntotspliti,noverspliti,ntotpdf,noverpdf,
     &ntotxsec,noverxsec,ntotsuda,noversuda
C--local variables
  INTEGER qclose,qbin,i
  DOUBLE PRECISION q,xa(2),ya(2),y,a,b
  CHARACTER*2 typ2

  ntotsuda=ntotsuda+1
  if (q.gt.sudaqq(1000,1)) then
    noversuda=noversuda+1
    if (noversuda.le.25) 
     &  write(logfid,*)'WARNING in getinsudared: need to extrapolate: ',
     &  q,sudaqq(1000,1)
  endif

      qclose=int((log(q)-log(sudaqq(1,1)))*999.d0
     &  /(log(sudaqq(1000,1))-log(sudaqq(1,1)))+1)
      qbin=max(qclose,1)
      qbin=min(qbin,999)
  IF((qbin.GT.999).OR.(qbin.LT.1)) THEN
       write(logfid,*)
     &  'ERROR in GETINSUDARED: line number out of bound',qbin
  ENDIF
  IF(typ2.EQ.'QQ')THEN
       DO 16 i=1,2
        xa(i)=sudaqq(qbin-1+i,1)
        ya(i)=sudaqq(qbin-1+i,2)
 16    CONTINUE
  ELSEIF(typ2.EQ.'QG')THEN
       DO 17 i=1,2
        xa(i)=sudaqg(qbin-1+i,1)
        ya(i)=sudaqg(qbin-1+i,2)
 17    CONTINUE
  ELSEIF(typ2.EQ.'GG')THEN
       DO 18 i=1,2
        xa(i)=sudagg(qbin-1+i,1)
        ya(i)=sudagg(qbin-1+i,2)
 18    CONTINUE
  ELSEIF(typ2.EQ.'GC')THEN
       DO 19 i=1,2
        xa(i)=sudagc(qbin-1+i,1)
        ya(i)=sudagc(qbin-1+i,2)
 19    CONTINUE
  ELSE
   write(logfid,*)'error in GETINSUDARED: unknown type ',typ2
  ENDIF
  a=(ya(2)-ya(1))/(xa(2)-xa(1))
  b=ya(1)-a*xa(1)
  y=a*q+b
  getinsudared=y
  IF(getinsudared.LT.(-1.d-10))THEN
   write(logfid,*) 'ERROR: GETINSUDARED < 0:',getinsudared,q,typ2
  ENDIF
  if (getinsudared.lt.0.d0) getinsudared = 0.d0
  END


***********************************************************************
***   function getsscat
***********************************************************************
      DOUBLE PRECISION FUNCTION getsscat(EN,px,py,PZ,MP,LW,TYPE1,TYPE2,
     &  x,y,z,t,mode)
      IMPLICIT NONE
C--identifier of file for hepmc output and logfile
  common/hepmcid/hpmcfid,logfid
  integer hpmcfid,logfid
C--Parameter common block
  common/param/q0,lps,lqcd,ltime,scalefacm,angord,scatrecoil,
     &allhad,compress,nf
      INTEGER nf
  DOUBLE PRECISION q0,lqcd,ltime,lps,scalefacm
      LOGICAL angord,scatrecoil,allhad,compress
C--variables for cross section integration 
  common/xsecv/qlow,mdx
  DOUBLE PRECISION qlow,mdx
C--local variables
  integer mode
      DOUBLE PRECISION up,en,lw,scatprimfunc,ccol,mp,
     &low,getpdfxint,getxsecint,mdeb,pz,pcms2,shat,gettemp,
     &x,y,z,t,getmd,avmom(5),px,py,getmdmin,getmdmax,pproj,psct
      CHARACTER type1,type2

       IF(type1.EQ.'Q')THEN
        ccol=2./3.
       ELSE
        ccol=3./2.
       ENDIF 
   if (mode.eq.0) then
     mdeb = getmd(x,y,z,t)
     call avscatcen(x,y,z,t,
     &  avmom(1),avmom(2),avmom(3),avmom(4),avmom(5))
     shat = avmom(5)**2 + mp**2 + 
     &  2.*(avmom(4)*en - avmom(1)*px - avmom(2)*py - avmom(3)*pz)
     pcms2 = (shat+mp**2-avmom(5)**2)**2/(4.*shat)-mp**2
     up = 4.*pcms2
   else
     if (mode.eq.1) then
       mdeb = getmdmin()
     else 
       mdeb = getmdmax()
     endif 
     call maxscatcen(avmom(1),avmom(2),avmom(3),avmom(4),avmom(5))
     psct = sqrt(avmom(1)**2+avmom(2)**2+avmom(3)**2)
     pproj = sqrt(px**2+py**2+pz**2)
     shat = avmom(5)**2 + mp**2 + 2.*(en*avmom(4) + pproj*psct)
     pcms2 = (shat+mp**2-avmom(5)**2)**2/(4.*shat)-mp**2
     up = 4.*pcms2
   endif
   low=lw**2
   IF(low.GT.up)THEN
    getsscat=0.d0
    RETURN
   ENDIF
   IF((type2.EQ.'C').OR.
     &  ((type1.EQ.'Q').AND.(type2.EQ.'Q')).OR.
     &    ((type1.EQ.'G').AND.(type2.EQ.'G')))THEN
        getsscat=ccol*(scatprimfunc(up,mdeb)-scatprimfunc(low,mdeb))
   ELSE
    getsscat=0.d0
   ENDIF
   low=q0**2/scalefacm**2
   IF(up.GT.low)THEN
        IF(type1.EQ.'Q')THEN
     IF((type2.EQ.'C').OR.(type2.EQ.'G'))THEN
      getsscat=getsscat+getpdfxint(scalefacm*sqrt(up),'GQ')
     &  *3.*scatprimfunc(up,mdeb)/2.
      getsscat=getsscat-getxsecint(up,mdeb,'QA')
     ENDIF
    ELSE
     IF((type2.EQ.'C').OR.(type2.EQ.'G'))THEN
      getsscat=getsscat+ccol*(scatprimfunc(up,mdeb)-
     &      scatprimfunc(low,mdeb))
     &    - getxsecint(up,mdeb,'GB')
     ENDIF
     IF((type2.EQ.'C').OR.(type2.EQ.'Q'))THEN
      getsscat=getsscat+2.*getpdfxint(scalefacm*sqrt(up),'QG')
     &  *2.*scatprimfunc(up,mdeb)/3.
      getsscat=getsscat-2.*getxsecint(up,mdeb,'GA')
     ENDIF
    ENDIF
   ENDIF
  IF(getsscat.LT.-1.d-4)
     &    write(logfid,*) 'error: cross section < 0',getsscat,'for',
     &  en,mp,lw,type1,type2,lw**2,up
  getsscat=max(getsscat,0.d0)
      END



***********************************************************************
***   function getmass
***********************************************************************
  DOUBLE PRECISION FUNCTION getmass(QBMIN,QBMAX,THETA,EP,TYPE,
     &                                   max2,ins,zdec,qqbardec)
  IMPLICIT NONE
C--identifier of file for hepmc output and logfile
  common/hepmcid/hpmcfid,logfid
  integer hpmcfid,logfid
C--Common block of Pythia
      common/pyjets/n,npad,k(23000,5),p(23000,5),v(23000,5)
  INTEGER n,npad,k
  DOUBLE PRECISION p,v
      common/pydat1/mstu(200),paru(200),mstj(200),parj(200)
  INTEGER mstu,mstj
  DOUBLE PRECISION paru,parj
      common/pydat3/mdcy(500,3),mdme(8000,2),brat(8000),kfdp(8000,5)
  INTEGER mdcy,mdme,kfdp
  DOUBLE PRECISION brat
C--Parameter common block
  common/param/q0,lps,lqcd,ltime,scalefacm,angord,scatrecoil,
     &allhad,compress,nf
      INTEGER nf
  DOUBLE PRECISION q0,lqcd,ltime,lps,scalefacm
      LOGICAL angord,scatrecoil,allhad,compress
C--time common block
      common/time/mv(23000,5)
      DOUBLE PRECISION mv
C--factor in front of alphas argument
  common/alphasfac/ptfac
  DOUBLE PRECISION ptfac
C--local variables
  DOUBLE PRECISION qbmin,qbmax,theta,ep,max2,zdec,
     &q2min,alphmax,alphas,log14,pref,q2max,sudaover,gmin,
     &gmax,arg,cand,eps,trueeps,trueval,oest,weight,getinspliti,
     &r,pyr,z,rz,thetanew,r2,pi,pqq,pgg,pqg,rmin
      CHARACTER*2 type
  LOGICAL ins,qqbardec
      DATA pi/3.141592653589793d0/
  
  q2min = q0**2

  alphmax = alphas(3.*ptfac*q2min/16.,lps)
  log14 = log(0.25)

      IF(type.EQ.'QQ')THEN
   pref=4.*alphmax/(3.*2.*pi)
      ELSE
   pref=29.*alphmax/(8.*2.*pi)
      ENDIF

C--check if phase space available, return 0.d0 otherwise
  IF((qbmax.LE.qbmin).OR.(ep.LT.qbmin)) THEN
   getmass=0.d0
   zdec=0.d0
   qqbardec=.false.
   RETURN
  ENDIF

      q2max = qbmax**2
! 21  sudaover = exp(-pref*(log(q2min/(4.*q2max))**2 - log14**2))
! IF(pyr(0).LE.sudaover)THEN
 21   if (q2max-qbmin**2.lt.1e-4)then
      getmass=qbmin
      zdec=0.5
      IF(type.EQ.'QQ')THEN
        qqbardec=.false.
      ELSE
        IF(pyr(0).LT.pqg(0.5d0)/(pqg(0.5d0)+pgg(0.5d0)))THEN
          qqbardec=.true.
        ELSE 
          qqbardec=.false.
        ENDIF
      endif
      return
        endif
        gmax = pref*log(q2min/(4.*q2max))**2
        if (qbmin.gt.0.d0) then
          rmin = exp(pref*log(q2min/(4.*qbmin**2))**2-gmax)
        else
      rmin = 0.d0
    endif  
    
       r=pyr(0)*(1.d0-rmin)+rmin
       arg=gmax+log(r)
       if(arg.lt.0.d0)then
   getmass=0.d0
   zdec=0.d0
   qqbardec=.false.
   RETURN
  endif
! r=pyr(0)
! gmin = pref*log14**2
! gmax = pref*log(q2min/(4.*q2max))**2
! arg = log(r*exp(gmax)+(1.-r)*exp(gmin))
  cand = q2min*exp(sqrt(arg/pref))/4.
  eps = q2min/(4.*cand)

  if ((cand.lt.q2min).or.(cand.lt.qbmin**2)) then
   getmass=0.d0
   zdec=0.d0
   qqbardec=.false.
   RETURN
  endif

  IF((cand.GT.max2**2).OR.(cand.GT.ep**2))THEN
   q2max=cand
   goto 21
  ENDIF

  if (ins) then
    trueval=getinspliti(sqrt(cand),type)
    oest = -2.*pref*log(eps)
        weight = trueval/oest
  else
C--find true z interval
        trueeps=0.5-0.5*sqrt(1.-q2min/cand)
     &  *sqrt(1.-cand/ep**2)
        IF(trueeps.LT.eps)
     &  WRITE(logfid,*)'error in getmass: true eps < eps',trueeps,eps
    rz=pyr(0)
    z = 1.-eps**rz
    if ((z.lt.trueeps).or.(z.gt.(1.-trueeps))) then
      weight = 0.
    else
      if (type.eq.'QQ')then
!       if (ins) then
!                trueval = alphas(ptfac*(1.-z)*cand,lps)*pqq(z)/(2.*pi)
!              else
          trueval = alphas(ptfac*z*(1.-z)*cand,lps)*pqq(z)/(2.*pi)
!              endif
        oest = 2.*pref/(1.-z)
        weight = trueval/oest
      else
        if (pyr(0).lt.(17./29.)) z = 1.-z
!       if (ins)then
!         trueval = alphas(ptfac*(1.-z)*cand,lps)
!     &     *(pgg(z)+pqg(z))/(2.*pi)
!              else
          trueval = alphas(ptfac*z*(1.-z)*cand,lps)
     &      *(pgg(z)+pqg(z))/(2.*pi)
!              endif
        oest = alphmax*(17./(4.*z)+3./(1.-z))/(2.*pi)
        weight = trueval/oest
      endif
      thetanew = sqrt(cand/(z*(1.-z)))/ep
      if (angord.and.(theta.gt.0.).and.(thetanew.gt.theta)) 
     &                weight = 0.d0
    endif
  endif
  IF (weight.GT.1.d0) WRITE(logfid,*) 
     &  'problem in getmass: weight> 1',
     &    weight,type,eps,trueeps,z,cand
  r2=pyr(0)
  IF(r2.GT.weight)THEN
   q2max=cand
   goto 21
  ELSE
   getmass=sqrt(cand)
   if (.not.ins) then
     zdec=z
     IF(type.EQ.'QQ')THEN
       qqbardec=.false.
     ELSE
       IF(pyr(0).LT.pqg(z)/(pqg(z)+pgg(z)))THEN
         qqbardec=.true.
       ELSE 
         qqbardec=.false.
       ENDIF
     ENDIF
    endif
  ENDIF
  END



***********************************************************************
***   function generatez
***********************************************************************
  DOUBLE PRECISION FUNCTION generatez(TI,EA,EPSI,TYPE)
      IMPLICIT NONE
C--Parameter common block
  common/param/q0,lps,lqcd,ltime,scalefacm,angord,scatrecoil,
     &allhad,compress,nf
      INTEGER nf
  DOUBLE PRECISION q0,lqcd,ltime,lps,scalefacm
      LOGICAL angord,scatrecoil,allhad,compress
C--local variables
      DOUBLE PRECISION ti,ea,eps,pyr,x,r,help,r1,epsi
  CHARACTER*2 type

      IF(ti.EQ.0.d0)THEN
       eps=epsi
      ELSE
       eps=max(0.5-0.5*sqrt(1.-q0**2/ti)
     &      *sqrt(1.-ti/ea**2),epsi)
      ENDIF
      IF(eps.GT.0.5)THEN
       generatez=0.5
       goto 61
      ENDIF
 60   r=pyr(0)
  IF(type.EQ.'QQ')THEN
       x=1.-(1.-eps)*(eps/(1.-eps))**r
       r=pyr(0)
       IF(r.LT.((1.+x**2)/2.))THEN
        generatez=x
       ELSE
        goto 60
       ENDIF
  ELSEIF(type.EQ.'GG')THEN
       x=1./(1.+((1.-eps)/eps)**(1.-2.*r))
       r=pyr(0)
   help=((1.-x)/x+x/(1.-x)+x*(1.-x))/(1./(1.-x)+1./x)
       IF(r.LT.help)THEN
        generatez=x
       ELSE
        goto 60
       ENDIF
  ELSE
   r=pyr(0)*(1.-2.*eps)+eps
   r1=pyr(0)/2.
   help=0.5*(r**2+(1.-r)**2)
   IF(r1.LT.help)THEN
    generatez=r
   ELSE
    goto 60
   ENDIF
  ENDIF
 61 END



***********************************************************************
***   function scatprimfunc
***********************************************************************
      DOUBLE PRECISION FUNCTION scatprimfunc(T,MDEB)
      IMPLICIT NONE
C--Parameter common block
  common/param/q0,lps,lqcd,ltime,scalefacm,angord,scatrecoil,
     &allhad,compress,nf
      INTEGER nf
  DOUBLE PRECISION q0,lqcd,ltime,lps,scalefacm
      LOGICAL angord,scatrecoil,allhad,compress
C--local variables
      DOUBLE PRECISION t,pi,s,ei,alphas,t1,mdeb
      DATA pi/3.141592653589793d0/

   scatprimfunc = 2.*pi*(12.*pi)**2*(
     &  - ei(-log((t+mdeb**2)/lqcd**2))/lqcd**2
     &  - 1./((t+mdeb**2)*log((t+mdeb**2)/lqcd**2)))/(33.-2.*nf)**2
      END



***********************************************************************
***   function intpqq
***********************************************************************
  DOUBLE PRECISION FUNCTION intpqq(Z,Q)
  IMPLICIT NONE
C--Parameter common block
  common/param/q0,lps,lqcd,ltime,scalefacm,angord,scatrecoil,
     &allhad,compress,nf
      INTEGER nf
  DOUBLE PRECISION q0,lqcd,ltime,lps,scalefacm
      LOGICAL angord,scatrecoil,allhad,compress
C--local variables
  DOUBLE PRECISION z,q

  intpqq=6.*4.*(-2.*log(log(q**2/lps**2)
     &  +log(1.-z)))/((33.-2.*nf)*3.)
  END



***********************************************************************
***   function intpgglow
***********************************************************************
  DOUBLE PRECISION FUNCTION intpgglow(Z,Q)
  IMPLICIT NONE
C--Parameter common block
  common/param/q0,lps,lqcd,ltime,scalefacm,angord,scatrecoil,
     &allhad,compress,nf
      INTEGER nf
  DOUBLE PRECISION q0,lqcd,ltime,lps,scalefacm
      LOGICAL angord,scatrecoil,allhad,compress
C--local variables
  DOUBLE PRECISION z,q

  intpgglow=6.*3.*(log(log(q**2/lps**2)+log(z)))/(33.-2.*nf)
  END
  


***********************************************************************
***   function intpgghigh
***********************************************************************
  DOUBLE PRECISION FUNCTION intpgghigh(Z,Q)
  IMPLICIT NONE
C--Parameter common block
  common/param/q0,lps,lqcd,ltime,scalefacm,angord,scatrecoil,
     &allhad,compress,nf
      INTEGER nf
  DOUBLE PRECISION q0,lqcd,ltime,lps,scalefacm
      LOGICAL angord,scatrecoil,allhad,compress
C--local variables
  DOUBLE PRECISION z,q

  intpgghigh=-6.*3.*(log(log(q**2/lps**2)+log(1.-z)))/(33.-2.*nf)
  END
  


***********************************************************************
***   function intpqglow
***********************************************************************
  DOUBLE PRECISION FUNCTION intpqglow(Z,Q)
  IMPLICIT NONE
C--Parameter common block
  common/param/q0,lps,lqcd,ltime,scalefacm,angord,scatrecoil,
     &allhad,compress,nf
      INTEGER nf
  DOUBLE PRECISION q0,lqcd,ltime,lps,scalefacm
      LOGICAL angord,scatrecoil,allhad,compress
C--local variables
  DOUBLE PRECISION z,q,ei

  intpqglow=6.*(lps**2*ei(log(q**2/lps**2)+log(z))/q**2 
     & - 2.*lps**4*ei(2.*(log(q**2/lps**2)+log(z)))/q**4
     & + 2.*lps**6*ei(3.*(log(q**2/lps**2)+log(z)))/q**6)/
     &((33.-2.*nf)*2.)
  END
  


***********************************************************************
***   function intpqghigh
***********************************************************************
  DOUBLE PRECISION FUNCTION intpqghigh(Z,Q)
  IMPLICIT NONE
C--Parameter common block
  common/param/q0,lps,lqcd,ltime,scalefacm,angord,scatrecoil,
     &allhad,compress,nf
      INTEGER nf
  DOUBLE PRECISION q0,lqcd,ltime,lps,scalefacm
      LOGICAL angord,scatrecoil,allhad,compress
C--local variables
  DOUBLE PRECISION z,q,ei

  intpqghigh=-6.*(lps**2*ei(log(q**2/lps**2)+log(1.-z))/q**2 
     & - 2.*lps**4*ei(2.*(log(q**2/lps**2)+log(1.-z)))/q**4
     & + 2.*lps**6*ei(3.*(log(q**2/lps**2)+log(1.-z)))/q**6)/
     &((33.-2.*nf)*2.)
  END



***********************************************************************
***   function gett
***********************************************************************
  DOUBLE PRECISION FUNCTION gett(MINT,MAXT,MDEB)
  IMPLICIT NONE
C--Parameter common block
  common/param/q0,lps,lqcd,ltime,scalefacm,angord,scatrecoil,
     &allhad,compress,nf
      INTEGER nf
  DOUBLE PRECISION q0,lqcd,ltime,lps,scalefacm
      LOGICAL angord,scatrecoil,allhad,compress
C--local variables
  DOUBLE PRECISION tmin,tmax,maxi,pyr,r1,r2,alphas,pi,y,maxt,
     &mdeb,mint,t
  DATA pi/3.141592653589793d0/

  tmax=maxt+mdeb**2
  tmin=mint+mdeb**2
  IF(tmin.GT.tmax) THEN
   gett=0.d0
   RETURN
  ENDIF
 20 r1=pyr(0)
  t=tmax*tmin/(tmax+r1*(tmin-tmax))
  r2=pyr(0)
  IF(r2.LT.alphas(t,lqcd)**2/alphas(tmin,lqcd)**2)THEN
   gett=t-mdeb**2
  ELSE
   goto 20
  ENDIF

  END



***********************************************************************
***   function ei
***********************************************************************
      DOUBLE PRECISION FUNCTION ei(X)
      IMPLICIT NONE
C--identifier of file for hepmc output and logfile
  common/hepmcid/hpmcfid,logfid
  integer hpmcfid,logfid
C--exponential integral for negative arguments
      common/expint/eixs(3,1000),valmax,nval
      INTEGER nval
      DOUBLE PRECISION eixs,valmax
C--local variables
      INTEGER k,line,lmax
      DOUBLE PRECISION x,r,ga,xa(2),ya(2),y,dy,a,b
  DOUBLE PRECISION ystart,epsi,hfirst
  DATA epsi/1.e-5/
  
  IF(dabs(x).GT.valmax)
     &  write(logfid,*)'warning: value out of array in Ei(x)',x,valmax

      IF(x.GE.0.d0)THEN
       lmax=int(x*nval/valmax)
       line=max(lmax,1)
       line=min(line,999)
   IF((line.GT.999).OR.(line.LT.1)) THEN
        write(logfid,*)'ERROR in EI: line number out of bound',line
   ENDIF
       DO 26 k=1,2
        xa(k)=eixs(1,line-1+k)
        ya(k)=eixs(3,line-1+k)
 26    CONTINUE
   a=(ya(2)-ya(1))/(xa(2)-xa(1))
   b=ya(1)-a*xa(1)
   y=a*x+b
      ELSE
       lmax=int(-x*nval/valmax)
       line=max(lmax,1)
       line=min(line,999)
   IF((line.GT.999).OR.(line.LT.1)) THEN
        write(logfid,*)'ERROR in EI: line number out of bound',line
   ENDIF
       DO 27 k=1,2
        xa(k)=eixs(1,line-1+k)
        ya(k)=eixs(2,line-1+k)
 27    CONTINUE
   a=(ya(2)-ya(1))/(xa(2)-xa(1))
   b=ya(1)-a*xa(1)
   y=-a*x+b
      ENDIF
      ei=y
      END



***********************************************************************
***   function pqq
***********************************************************************
  DOUBLE PRECISION FUNCTION pqq(Z)
  IMPLICIT NONE
  DOUBLE PRECISION z
  pqq=4.*(1.+z**2)/(3.*(1.-z))
  END



***********************************************************************
***   function pgq
***********************************************************************
  DOUBLE PRECISION FUNCTION pgq(Z)
  IMPLICIT NONE
  DOUBLE PRECISION z
  pgq=4.*(1.+(1.-z)**2)/(3.*z)
  END



***********************************************************************
***   function pgg
***********************************************************************
  DOUBLE PRECISION FUNCTION pgg(Z)
  IMPLICIT NONE
  DOUBLE PRECISION z
  pgg=3.*((1.-z)/z + z/(1.-z) + z*(1.-z))
  END



***********************************************************************
***   function pqg
***********************************************************************
  DOUBLE PRECISION FUNCTION pqg(Z)
  IMPLICIT NONE
  DOUBLE PRECISION z
  pqg=0.5*(z**2 + (1.-z)**2)
  END



***********************************************************************
***   function alphas
***********************************************************************
  DOUBLE PRECISION FUNCTION alphas(T,LAMBDA)
  IMPLICIT NONE
C--Parameter common block
  common/param/q0,lps,lqcd,ltime,scalefacm,angord,scatrecoil,
     &allhad,compress,nf
      INTEGER nf
  DOUBLE PRECISION q0,lqcd,ltime,lps,scalefacm
      LOGICAL angord,scatrecoil,allhad,compress
C--local variables
  DOUBLE PRECISION t,l0,pi,lambda
  DATA pi/3.141592653589793d0/

   alphas=4.*pi/((11.-2.*nf/3.)*log(t/lambda**2))
  END



***********************************************************************
***   subroutine splitfncint
***********************************************************************
  SUBROUTINE splitfncint(EMAX)
  IMPLICIT NONE
C--Parameter common block
  common/param/q0,lps,lqcd,ltime,scalefacm,angord,scatrecoil,
     &allhad,compress,nf
      INTEGER nf
  DOUBLE PRECISION q0,lqcd,ltime,lps,scalefacm
      LOGICAL angord,scatrecoil,allhad,compress
C--splitting integral
      common/splitint/splitiggv(1000,1000),splitiqqv(1000,1000),
     &splitiqgv(1000,1000),qval(1000),zmval(1000),qmax,zmmin,npoint
      INTEGER npoint
      DOUBLE PRECISION splitiggv,splitiqqv,splitiqgv,
     &qval,zmval,qmax,zmmin
C--variables for splitting function integration
  common/intsplitf/qquad,fm
  DOUBLE PRECISION qquad,fm
C--max rapidity
  common/rapmax/etamax
  double precision etamax
C--local variables
  INTEGER nstep,i,j
  DOUBLE PRECISION emax,zmmax,epsi,hfirst,ystart,lnzmmin,
     &lnzmmax,zm,zm2,q,getmsmax,avmom(5),shat,pcms2
      DATA zmmax/0.5/
      DATA nstep/999/
  DATA epsi/1.d-5/

  call maxscatcen(avmom(1),avmom(2),avmom(3),avmom(4),avmom(5))
  shat = avmom(5)**2 +
     &    2.*emax*(avmom(4)+sqrt(avmom(1)**2+avmom(2)**2+avmom(3)**2))
  pcms2 = (shat-avmom(5)**2)**2/(4.*shat)
  qmax = sqrt(scalefacm*4.*pcms2)

  zmmin=q0/emax

      lnzmmin=log(zmmin)
      lnzmmax=log(zmmax)

  npoint=nstep

  DO 100 i=1,nstep+1
   q=(i-1)*(qmax-q0)/nstep+q0
       qval(i)=q
   qquad=q**2
       DO 110 j=1,nstep+1
        zm=exp((j-1)*(lnzmmax-lnzmmin)/nstep+lnzmmin)
        zmval(j)=zm
    IF(q**2.LT.q0**2)THEN
     zm2=0.5
    ELSE 
     zm2=0.5-0.5*sqrt(1.-q0**2/q**2)
    ENDIF 
    zm=max(zm,zm2)
    IF(zm.EQ.0.5)THEN 
     splitiqqv(i,j)=0.d0
     splitiggv(i,j)=0.d0
     splitiqgv(i,j)=0.d0
    ELSE
     ystart=0d0
     hfirst=0.01
     fm=0.d0
     CALL odeint(ystart,zm,1.-zm,epsi,hfirst,0d0,2)
     splitiqqv(i,j)=ystart
     ystart=0d0
     hfirst=0.01
     fm=0.d0
     CALL odeint(ystart,zm,1.-zm,epsi,hfirst,0d0,3)
     splitiggv(i,j)=ystart
     ystart=0d0
     hfirst=0.01
     fm=0.d0
     CALL odeint(ystart,zm,1.-zm,epsi,hfirst,0d0,4)
     splitiqgv(i,j)=ystart
    ENDIF
 110   CONTINUE
 100  CONTINUE

  END



***********************************************************************
***   subroutine pdfint
***********************************************************************
  SUBROUTINE pdfint(EMAX)
  IMPLICIT NONE
C--Parameter common block
  common/param/q0,lps,lqcd,ltime,scalefacm,angord,scatrecoil,
     &allhad,compress,nf
      INTEGER nf
  DOUBLE PRECISION q0,lqcd,ltime,lps,scalefacm
      LOGICAL angord,scatrecoil,allhad,compress
C--pdf common block
  common/pdfs/qinqx(2,1000),ginqx(2,1000),qingx(2,1000),
     &gingx(2,1000)
  DOUBLE PRECISION qinqx,ginqx,qingx,gingx
C--variables for pdf integration
  common/pdfintv/xmax,z
  DOUBLE PRECISION xmax,z
C--max rapidity
  common/rapmax/etamax
  double precision etamax
C--local variables
  INTEGER i,j
  DOUBLE PRECISION emax,q2,getpdfxintexact,ystart,hfirst,epsi,
     &q2max,deltaq2,avmom(5),shat,pcms2
  DATA epsi/1.d-4/

  call maxscatcen(avmom(1),avmom(2),avmom(3),avmom(4),avmom(5))
  shat = avmom(5)**2 +
     &    2.*emax*(avmom(4)+sqrt(avmom(1)**2+avmom(2)**2+avmom(3)**2))
  pcms2 = (shat-avmom(5)**2)**2/(4.*shat)
  q2max = scalefacm*4.*pcms2

  deltaq2=log(q2max)-log(q0**2)
  qinqx(1,1)=q0**2
  ginqx(1,1)=q0**2
  qingx(1,1)=q0**2
  gingx(1,1)=q0**2
  qinqx(2,1)=0.d0
  ginqx(2,1)=0.d0
  qingx(2,1)=0.d0
  gingx(2,1)=0.d0
   DO 12 j=2,1000
    q2 = exp((j-1)*deltaq2/999.d0 + log(q0**2))
    qinqx(1,j)=q2
    ginqx(1,j)=q2
    qingx(1,j)=q2
    gingx(1,j)=q2
    qinqx(2,j)=getpdfxintexact(sqrt(q2),'QQ')
    ginqx(2,j)=getpdfxintexact(sqrt(q2),'GQ')
    qingx(2,j)=getpdfxintexact(sqrt(q2),'QG')
    gingx(2,j)=getpdfxintexact(sqrt(q2),'GG')
 12  CONTINUE
  END



***********************************************************************
***   subroutine xsecint
***********************************************************************
  SUBROUTINE xsecint(EMAX)
  IMPLICIT NONE
C--Parameter common block
  common/param/q0,lps,lqcd,ltime,scalefacm,angord,scatrecoil,
     &allhad,compress,nf
      INTEGER nf
  DOUBLE PRECISION q0,lqcd,ltime,lps,scalefacm
      LOGICAL angord,scatrecoil,allhad,compress
C--cross secttion common block
  common/xsecs/intq1(1001,101),intq2(1001,101),
     &intg1(1001,101),intg2(1001,101)
  DOUBLE PRECISION intq1,intq2,intg1,intg2
C--variables for cross section integration 
  common/xsecv/qlow,mdx
  DOUBLE PRECISION qlow,mdx
C--max rapidity
  common/rapmax/etamax
  double precision etamax
C--local variables
  INTEGER j,k
  DOUBLE PRECISION emax,tmax,tmaxmax,deltatmax,ystart,hfirst,epsi,
     &getmsmax,getmdmax,mdmin,mdmax,deltamd,getmdmin,avmom(5),shat,pcms2
  DATA epsi/1.d-4/

  call maxscatcen(avmom(1),avmom(2),avmom(3),avmom(4),avmom(5))
  shat = avmom(5)**2 +
     &    2.*emax*(avmom(4)+sqrt(avmom(1)**2+avmom(2)**2+avmom(3)**2))
  pcms2 = (shat-avmom(5)**2)**2/(4.*shat)
  tmaxmax = scalefacm*4.*pcms2
  deltatmax=(log(tmaxmax)-
     &  log(q0**2*(1.d0+1.d-6)/scalefacm**2))/999.d0
      mdmin=getmdmin()
      mdmax=max(mdmin,getmdmax())
      deltamd=(mdmax-mdmin)/99.d0

   DO 12 j=1,1000
    tmax = exp((j-1)*deltatmax
     &    + log(q0**2*(1.d0+1.d-6)/scalefacm**2))
    intq1(j,101)=tmax
    intq2(j,101)=tmax
    intg1(j,101)=tmax
    intg2(j,101)=tmax
        DO 13 k=1,100
         mdx=mdmin+(k-1)*deltamd
         intq1(1001,k)=mdx
         intq2(1001,k)=mdx
         intg1(1001,k)=mdx
         intg2(1001,k)=mdx
    IF(tmax.LT.q0**2/scalefacm**2)THEN
     intq1(j,k)=0.d0
     intq2(j,k)=0.d0
     intg1(j,k)=0.d0
     intg2(j,k)=0.d0
    ELSE
C--first quark integral
     qlow=q0
       hfirst=0.01*(tmax-q0**2/scalefacm**2)
         ystart=0.d0
        CALL odeint(ystart,q0**2/scalefacm**2,tmax,epsi,hfirst
     &        ,0.d0,11)
     intq1(j,k)=ystart
C--second quark integral
     qlow=q0
       hfirst=0.01*(tmax-q0**2/scalefacm**2)
         ystart=0.d0
        CALL odeint(ystart,q0**2/scalefacm**2,tmax,epsi,hfirst
     &        ,0.d0,14)
     intq2(j,k)=ystart
C--first gluon integral
     qlow=q0
         ystart=0.d0
        CALL odeint(ystart,q0**2/scalefacm**2,tmax,epsi,hfirst
     &        ,0.d0,12)
     intg1(j,k)=ystart
C--second gluon integral
     qlow=q0
         ystart=0.d0
        CALL odeint(ystart,q0**2/scalefacm**2,tmax,epsi,hfirst
     &        ,0.d0,13)
     intg2(j,k)=ystart
    ENDIF
 13     CONTINUE
 12  CONTINUE
  END



***********************************************************************
***   function insudaint
***********************************************************************
  SUBROUTINE insudaint(EMAX)
  IMPLICIT NONE
C--Parameter common block
  common/param/q0,lps,lqcd,ltime,scalefacm,angord,scatrecoil,
     &allhad,compress,nf
      INTEGER nf
  DOUBLE PRECISION q0,lqcd,ltime,lps,scalefacm
      LOGICAL angord,scatrecoil,allhad,compress
C--Sudakov common block
  common/insuda/sudaqq(1000,2),sudaqg(1000,2),sudagg(1000,2),
     &sudagc(1000,2)
  DOUBLE PRECISION sudaqq,sudaqg,sudagg,sudagc
C--max rapidity
  common/rapmax/etamax
  double precision etamax
C--local variables
  INTEGER i
  DOUBLE PRECISION qmax,q,getinsudakov,delta,emax,avmom(5),
     &shat,pcms2
  
  call maxscatcen(avmom(1),avmom(2),avmom(3),avmom(4),avmom(5))
  shat = avmom(5)**2 +
     &    2.*emax*(avmom(4)+sqrt(avmom(1)**2+avmom(2)**2+avmom(3)**2))
  pcms2 = (shat-avmom(5)**2)**2/(4.*shat)
  qmax = sqrt(scalefacm*4.*pcms2)
  delta=(log(3.*qmax)-log(q0**2*(1.d0+1.d-6)))/999.d0
  DO 22 i=1,1000
   q = exp((i-1)*delta + log(q0**2*(1.d0+1.d-6)))
   sudaqq(i,1)=q
   sudaqg(i,1)=q
   sudagg(i,1)=q
   sudagc(i,1)=q
   sudaqq(i,2)=getinsudakov(q0,q,'QQ')
   sudaqg(i,2)=getinsudakov(q0,q,'QG')
   sudagg(i,2)=getinsudakov(q0,q,'GG')
   sudagc(i,2)=getinsudakov(q0,q,'GC')
 22 CONTINUE
  END



***********************************************************************
***   function eixint
***********************************************************************
  SUBROUTINE eixint
  IMPLICIT NONE
C--exponential integral for negative arguments
      common/expint/eixs(3,1000),valmax,nval
      INTEGER nval
      DOUBLE PRECISION eixs,valmax
C-local variables
  INTEGER i,k
  DOUBLE PRECISION x,epsi,hfirst,ystart,ei,ga,r 
  DATA  epsi/1.d-5/

  nval=1000
  valmax=55.

      DO 10 i=1,nval
       x=i*valmax/(nval*1.d0)
       eixs(1,i)=x
C--do negative arguments first
   ystart=0d0
   hfirst=0.01
   CALL odeint(ystart,x,1000.d0,epsi,hfirst,0.d0,5)
       eixs(2,i)=-ystart
C--now do the positive arguments
       IF (x.EQ.0.0) THEN
        ei=-1.0d+300
       ELSE IF (x.LE.40.0) THEN
        ei=1.0d0
        r=1.0d0
        DO 15 k=1,100
         r=r*k*x/(k+1.0d0)**2
         ei=ei+r
         IF (dabs(r/ei).LE.1.0d-15) go to 20
15      CONTINUE
20      ga=0.5772156649015328d0
        ei=ga+dlog(x)+x*ei
       ELSE
        ei=1.0d0
        r=1.0d0
        DO 25 k=1,20
         r=r*k/x
25       ei=ei+r
         ei=dexp(x)/x*ei
       ENDIF
   eixs(3,i)=ei
 10   CONTINUE
  END



***********************************************************************
***   function odeint
***********************************************************************
  subroutine odeint(ystart,a,b,eps,h1,hmin,w1)
  implicit none
C--identifier of file for hepmc output and logfile
  common/hepmcid/hpmcfid,logfid
  integer hpmcfid,logfid
C--local variables
  integer nmax,nstep,w1
  double precision ystart,a,b,eps,h1,hmin,x,h,y,dydx,
     &deriv,yscale,hdid,hnew
  data nmax/100000/

  x = a
  y = ystart
  h = sign(h1,b-a)
  do 20 nstep=1,nmax
    dydx = deriv(x,w1)
    yscale = abs(y) + abs(h*dydx) + 1.e-25
    if (((x + h - b)*h).gt.0.) h = b-x
    call rkstepper(x,y,dydx,h,hdid,hnew,yscale,eps,w1)
    if ((x - b)*h.ge.0) then
      ystart = y
      return
    endif
    h = hnew
    if (abs(h).lt.abs(hmin)) then
      write(logfid,*)'Error in odeint: stepsize too small',w1
     &  ,ystart,a,b,h1
      return
    endif   
 20 continue
  write(logfid,*)'Error in odeint: too many steps',w1
     &  ,ystart,a,b,h1
  end



***********************************************************************
***   function rkstepper
***********************************************************************
  subroutine rkstepper(x,y,dydx,htest,hdid,hnew,yscale,eps,w1)
  implicit none
C--identifier of file for hepmc output and logfile
  common/hepmcid/hpmcfid,logfid
  integer hpmcfid,logfid
C--local variables
  integer w1
  double precision x,y,dydx,htest,hdid,hnew,yscale,eps,
     &yhalf,y1,y2,rk4step,dydxhalf,xnew,delta,err,h,safety, powerdown,
     &powerup,maxup,maxdown,deriv,fac
  logical reject
  data powerdown/0.25/
  data powerup/0.2/
  data safety/0.9/
  data maxdown/10./
  data maxup/5./

  reject = .false.
  h = htest
 10 xnew = x + h
  if (x.eq.xnew) then
    write(logfid,*)'Error in rkstepper: step size not significant'
    return
  endif
  yhalf = rk4step(x,y,dydx,h/2.,w1)
  dydxhalf = deriv(x+h/2.,w1)
  y2 = rk4step(x+h/2.,yhalf,dydxhalf,h/2.,w1)
  y1 = rk4step(x,y,dydx,h,w1)
  delta = y2-y1
  err = abs(delta)/(yscale*eps)
  if (err.gt.1.) then
    reject = .true.
    fac = max(1./maxdown,safety/err**powerdown)
    h = h*fac
    goto 10 
  else
    if (reject) then
      hnew = h
    else
      fac = min(maxup,safety/err**powerup)
      hnew = fac*h
    endif
    x = xnew
    y = y2 + delta/15.
    hdid = h
  endif
  end



***********************************************************************
***   function rk4step
***********************************************************************
  double precision function rk4step(x,y,dydx,h,w1)
  implicit none
  integer w1
  double precision x,y,dydx,h,k1,k2,k4,yout,deriv
  k1 = h*dydx
  k2 = h*deriv(x+h/2.,w1)
  k4 = h*deriv(x+h,w1)
  yout = y+k1/6.+2.*k2/3.+k4/6.
  rk4step = yout
  end



***********************************************************************
***   function getdeltat
***********************************************************************
      LOGICAL FUNCTION getdeltat(LINE,TSTART,DTMAX1,DELTAT)
      IMPLICIT NONE
C--identifier of file for hepmc output and logfile
  common/hepmcid/hpmcfid,logfid
  integer hpmcfid,logfid
C--pythia common block
      common/pyjets/n,npad,k(23000,5),p(23000,5),v(23000,5)
  INTEGER n,npad,k
  DOUBLE PRECISION p,v
C--Parameter common block
  common/param/q0,lps,lqcd,ltime,scalefacm,angord,scatrecoil,
     &allhad,compress,nf
      INTEGER nf
  DOUBLE PRECISION q0,lqcd,ltime,lps,scalefacm
      LOGICAL angord,scatrecoil,allhad,compress
C--time common block
      common/time/mv(23000,5)
      DOUBLE PRECISION mv
C--max rapidity
  common/rapmax/etamax
  double precision etamax
C--memory for error message from getdeltat
  common/errline/errl
  integer errl
C--local variables
      INTEGER line,i,nnull
      DOUBLE PRECISION dtmax,sigmamax,neffmax,linvmax,pyr,
     &r,toff,xs,ys,zs,ts,getsscat,getmsmax,getmdmin,msmax,mdmin,
     &xstart,ystart,zstart,weight,ms,md,neff,sigma,getneff,
     &getneffmax,getms,getmd,tau,mdmax,getmdmax,getnatmdmin,
     &sigmamin,neffmin,tstart,dtmax1,deltat
  CHARACTER ptype
  LOGICAL stopnow

C--initialization
  getdeltat=.false.
      deltat=0.d0
  dtmax=dtmax1
  IF(k(line,2).EQ.21)THEN
   ptype='G'
  ELSE
   ptype='Q'
  ENDIF

  nnull=0
  stopnow=.false.

C--check for upper bound from plasma lifetime
      IF((tstart+dtmax).GE.ltime)dtmax=ltime-tstart
      IF(dtmax.LT.0.d0) RETURN
  
C--calculate time relative to production of the considered parton
      toff=tstart-mv(line,4)
  xstart=mv(line,1)+toff*p(line,1)/p(line,4)
  ystart=mv(line,2)+toff*p(line,2)/p(line,4)
  zstart=mv(line,3)+toff*p(line,3)/p(line,4)

C--calculate upper limit for density*cross section
  sigmamax=getsscat(p(line,4),p(line,1),p(line,2),p(line,3),
!     & xstart,ystart,-sign(abs(zstart),p(line,3)),zstart+1.d-6)
     &  p(line,5),0.d0,ptype,'C',xstart,ystart,zstart,tstart,1)
  sigmamin=getsscat(p(line,4),p(line,1),p(line,2),p(line,3),
!     & xstart,ystart,-sign(abs(zstart),p(line,3)),zstart+1.d-6)
     &  p(line,5),0.d0,ptype,'C',xstart,ystart,zstart,tstart,2)
  neffmax=getneffmax()
  neffmin=getnatmdmin()
  linvmax=5.d0*max(neffmin*sigmamax,neffmax*sigmamin)
  if(linvmax.eq.0.d0) return

  DO 333 i=1,1000000
   deltat=deltat-log(pyr(0))/linvmax
   xs=xstart+deltat*p(line,1)/p(line,4)
   ys=ystart+deltat*p(line,2)/p(line,4)
   zs=zstart+deltat*p(line,3)/p(line,4)
   ts=tstart+deltat
   IF(ts.LT.zs)THEN
    tau=-1.d0
   ELSE
    tau=sqrt(ts**2-zs**2)
   ENDIF
   neff=getneff(xs,ys,zs,ts)
   IF((tau.GT.1.d0).AND.(neff.EQ.0.d0))THEN
    IF(nnull.GT.4)THEN
     stopnow=.true.
    ELSE 
     nnull=nnull+1
    ENDIF
   ELSE
    nnull=0
   ENDIF
   IF((deltat.GT.dtmax).OR.stopnow) THEN
    deltat=dtmax
    RETURN
   ENDIF
   IF(neff.GT.0.d0)THEN
    sigma=getsscat(p(line,4),p(line,1),p(line,2),p(line,3),
     &  p(line,5),0.d0,ptype,'C',xs,ys,zs,ts,0)
   ELSE
    sigma=0.d0
   ENDIF
   weight=5.d0*neff*sigma/linvmax
   IF(weight.GT.1.d0+1d-6) then
     if (line.ne.errl) then
           write(logfid,*)'error in GETDELTAT: weight > 1',weight,
     &   neff*sigma/(neffmax*sigmamin),neff*sigma/(neffmin*sigmamax),
     &       p(line,4)
       errl=line
     endif
   endif
       r=pyr(0)
   IF(r.LT.weight)THEN
    getdeltat=.true.
    RETURN
   ENDIF
 333  CONTINUE
  END


  integer function poissonian(lambda)
  implicit none
  integer n
  double precision lambda,disc,p,pyr,u,v,pi
  data pi/3.141592653589793d0/
  
  if (lambda.gt.745.d0) then
    u = pyr(0);
    v = pyr(0);
    poissonian = 
     &  int(sqrt(lambda)*sqrt(-2.*log(u))*cos(2.*pi*v)+lambda)
  else
   disc=exp(-lambda)
   p=1.d0
   n=0  
 800   p = p*pyr(0)
   if (p.gt.disc) then
     n = n+1
     goto 800
   endif
   poissonian=n
  endif
  end


***********************************************************************
***   function ishadron
***********************************************************************
  LOGICAL FUNCTION ishadron(ID)
  IMPLICIT NONE
C--local variables
  INTEGER id  
  IF(abs(id).LT.100) THEN
   ishadron=.false.
  ELSE
   IF(mod(int(abs(id)/10.),10).EQ.0) THEN
    ishadron = .false.
   ELSE
    ishadron = .true.
       ENDIF
      ENDIF
      END



***********************************************************************
***   function isdiquark
***********************************************************************
  LOGICAL FUNCTION isdiquark(ID)
  IMPLICIT NONE
C--local variables
  INTEGER id  
  IF(abs(id).LT.1000) THEN
   isdiquark=.false.
  ELSE 
   IF(mod(int(id/10),10).EQ.0) THEN
    isdiquark = .true.
   ELSE
    isdiquark = .false.
       ENDIF
      ENDIF 
      END

***********************************************************************
***   function islepton
***********************************************************************
      LOGICAL FUNCTION islepton(ID)
      IMPLICIT NONE
C--   local variables
      INTEGER id
      IF((abs(id).EQ.11).OR.(abs(id).EQ.13).OR.(abs(id).EQ.15)) THEN
         islepton=.true.
      ELSE
         islepton=.false.
      ENDIF
      END
      
***********************************************************************
***   function isparton
***********************************************************************
  LOGICAL FUNCTION isparton(ID)
  IMPLICIT NONE
C--local variables
  INTEGER id  
  LOGICAL isdiquark
  IF((abs(id).LT.6).OR.(id.EQ.21).OR.isdiquark(id)) THEN
   isparton=.true.
  ELSE 
   isparton=.false.
      ENDIF 
      END      



***********************************************************************
***   function isprimstring
***********************************************************************
      logical function isprimstring(l)
      implicit none
      common/pyjets/n,npad,k(23000,5),p(23000,5),v(23000,5)
  INTEGER n,npad,k
  DOUBLE PRECISION p,v
C--local variables
  integer l
  logical isparton
  if ((k(l,2).ne.91).and.(k(l,2).ne.92)) then
    isprimstring=.false.
    return
  endif
  if ((k(k(l,3),3).eq.0).or.(isparton(k(k(k(l,3),3),2)))) then
        isprimstring=.true.
  else 
        isprimstring=.false.
  endif
  end



***********************************************************************
***   function issecstring
***********************************************************************
      logical function issecstring(l)
      implicit none
      common/pyjets/n,npad,k(23000,5),p(23000,5),v(23000,5)
  INTEGER n,npad,k
  DOUBLE PRECISION p,v
C--local variables
  integer l
  logical isparton,isprimstring
  if ((k(l,2).ne.91).and.(k(l,2).ne.92)) then
    issecstring = .false.
    return
  endif
  if (isprimstring(l)) then
    issecstring = .false.
    return
  endif
  if (isparton(k(k(k(l,3),3),2))) then 
    issecstring = .false.
  else
    issecstring = .true.
  endif
  end



***********************************************************************
***   function isprimhadron
***********************************************************************
      logical function isprimhadron(l)
      implicit none
      common/pyjets/n,npad,k(23000,5),p(23000,5),v(23000,5)
  INTEGER n,npad,k
  DOUBLE PRECISION p,v
C--local variables
  integer l
  logical isprimstring,isparton
  if (((k(k(l,3),2).EQ.91).OR.(k(k(l,3),2).EQ.92))
     &  .and.isprimstring(k(l,3))
     &  .and.(.not.isparton(k(l,2)))) then
    isprimhadron=.true.
  else 
        isprimhadron=.false.
  endif
  if (k(l,1).eq.17) isprimhadron=.true.
  end



***********************************************************************
***   function compressevent
***********************************************************************
  logical function compressevent(l1)
  implicit none
C--identifier of file for hepmc output and logfile
  common/hepmcid/hpmcfid,logfid
  integer hpmcfid,logfid
      common/pyjets/n,npad,k(23000,5),p(23000,5),v(23000,5)
  INTEGER n,npad,k
  DOUBLE PRECISION p,v
C--variables for angular ordering
      common/angor/za(23000),zd(23000),thetaa(23000),qqbard(23000)
  DOUBLE PRECISION za,zd,thetaa
      LOGICAL qqbard
C--time common block
      common/time/mv(23000,5)
      DOUBLE PRECISION mv
C--colour index common block
  common/colour/trip(23000),anti(23000),colmax
  INTEGER trip,anti,colmax
C--local variables
  integer l1,i,j,nold,nnew,nstart
  
  nold = n

  do 777 i=2,nold
    if (((k(i,1).eq.11).or.(k(i,1).eq.12).or.(k(i,1).eq.13)).and.
     &  (i.ne.l1)) then
      nnew = i
      goto 778
    endif
 777  continue
  compressevent = .false.
  return
 778  continue
  nstart = nnew
  do 779 i=nstart,nold
    if (((k(i,1).ne.11).and.(k(i,1).ne.12).and.(k(i,1).ne.13)).or.
     &  (i.eq.l1)) then
      do 780 j=1,5
        p(nnew,j)=p(i,j)
        v(nnew,j)=v(i,j)
        mv(nnew,j)=mv(i,j)
 780      continue
      trip(nnew)=trip(i)
      anti(nnew)=anti(i)
      za(nnew)=za(i)
      zd(nnew)=zd(i)
      thetaa(nnew)=thetaa(i)
      qqbard(nnew)=qqbard(i)
      k(nnew,1)=k(i,1)
      k(nnew,2)=k(i,2)
      k(nnew,3)=0
      k(nnew,4)=0
      k(nnew,5)=0
      if (l1.eq.i) l1=nnew
      nnew=nnew+1
    endif
 779  continue
  n=nnew-1
  if ((nold-n).le.10) then
    compressevent = .false.
  else
    compressevent = .true.
  endif
  do 781 i=nnew,nold
    do 782 j=1,5
      k(i,j)=0
      p(i,j)=0.d0
      v(i,j)=0.d0
      mv(i,j)=0.d0
 782    continue
    trip(i)=0
    anti(i)=0
    za(i)=0.d0
    zd(i)=0.d0
    thetaa(i)=0.d0
    qqbard(i)=.false.
 781  continue
  if (n.gt.23000) write(logfid,*)'Error in compressevent: n = ',n 
  if (l1.gt.n) write(logfid,*)'Error in compressevent: l1 = ',l1  
  call flush(logfid)
  return
  end



***********************************************************************
***   subroutine pevrec
***********************************************************************
      SUBROUTINE pevrec(NUM,COL)
C--identifier of file for hepmc output and logfile
  implicit none
  common/hepmcid/hpmcfid,logfid
  integer hpmcfid,logfid
      common/pyjets/n,npad,k(23000,5),p(23000,5),v(23000,5)
  INTEGER n,npad,k
  DOUBLE PRECISION p,v
C--variables for angular ordering
      common/angor/za(23000),zd(23000),thetaa(23000),qqbard(23000)
  DOUBLE PRECISION za,zd,thetaa
      LOGICAL qqbard
C--time common block
      common/time/mv(23000,5)
      DOUBLE PRECISION mv
C--colour index common block
  common/colour/trip(23000),anti(23000),colmax
  INTEGER trip,anti,colmax
  INTEGER num,i
  LOGICAL col

      DO 202 i=1,n
       v(i,1)=mv(i,1)
       v(i,2)=mv(i,2)
       v(i,3)=mv(i,3)
       v(i,4)=mv(i,4)
       v(i,5)=mv(i,5)
   IF(col) write(logfid,*)i,' (',trip(i),',',anti(i),')    [',
     &k(i,3),k(i,4),k(i,5),' ]  {',k(i,2),k(i,1),' } ',  
     &zd(i),thetaa(i)
 202  CONTINUE
      CALL pylist(num)

      END



***********************************************************************
***   subroutine converttohepmc
***********************************************************************
  SUBROUTINE converttohepmc(J,EVNUM,PID,beam1,beam2)
  IMPLICIT NONE
      common/pyjets/n,npad,k(23000,5),p(23000,5),v(23000,5)
  INTEGER n,npad,k
  DOUBLE PRECISION p,v
      common/pypars/mstp(200),parp(200),msti(200),pari(200)
  INTEGER mstp,msti
  DOUBLE PRECISION parp,pari
C--Parameter common block
  common/param/q0,lps,lqcd,ltime,scalefacm,angord,scatrecoil,
     &allhad,compress,nf
      INTEGER nf
  DOUBLE PRECISION q0,lqcd,ltime,lps,scalefacm
      LOGICAL angord,scatrecoil,allhad,compress
C--organisation of event record
  common/evrecord/nsim,npart,offset,hadrotype,sqrts,collider,hadro,
     &shorthepmc,channel,isochannel
  integer nsim,npart,offset,hadrotype
  double precision sqrts
  character*4 collider,channel
  character*2 isochannel
  logical hadro,shorthepmc
C--extra storage for scattering centres before interactions
      common/storescatcen/nscatcen,maxnscatcen,scatflav(10000),
     &scatcen(10000,5),writescatcen,writedummies
  integer nscatcen,maxnscatcen,scatflav
  double precision scatcen
  logical writescatcen,writedummies
C--local variables
  INTEGER evnum,pbarcode,vbarcode,codelist(25000),i,pid,nstart,
     &nfirst,nvertex,ntot,j,codefirst
      DOUBLE PRECISION mproton,mneutron,pdummy,pscatcen
      LOGICAL ishadron,isdiquark,isparton,isprimhadron,isprimstring,
     &issecstring
  character*2 beam1,beam2
  data mproton/0.9383/
  data mneutron/0.9396/
  data pdummy/1.d-6/  
  
 5000 FORMAT(a2,i10,i3,3e14.6,2i2,i6,4i2,e14.6)
 5100 FORMAT(a2,2e14.6)
 5200 FORMAT(a2,6i7,2i2,1i7,4e14.6)
 5300 FORMAT(a2,2i2,5e14.6,2i2)
 5400 FORMAT(a2,i6,6i2,i6,i2)
 5500 FORMAT(a2,i6,i6,5e14.6,3i2,i6,i2)

  pbarcode=0
  vbarcode=0

  if (shorthepmc) then
C--short output
        IF(collider.EQ.'EEJJ')THEN
          nvertex=3
      pbarcode=5
        ELSE
          nvertex=1
      pbarcode=2
        ENDIF
    nfirst = 0
    do 131 i=1,n
      if (((k(i,1).lt.6).or.(k(i,1).eq.17)))
     &  nfirst = nfirst+1
 131    continue
    if(writescatcen) nfirst=nfirst+nscatcen
    if(writedummies) nfirst=nfirst+nscatcen

    WRITE(j,5000)'E ',evnum,-1,0.d0,0.d0,0.d0,0,0,nvertex,1,2,0,1,
     &pari(10)
    WRITE(j,'(A2,I2,A5)')'N ',1,'"0"' 
    WRITE(j,'(A)')'U GEV MM'
    WRITE(j,5100)'C ',pari(1)*1.d9,0.d0
      WRITE(j,5200)'H ',0,0,0,0,0,0,0,0,0,0.d0,0.d0,0.d0,0.d0
    WRITE(j,5300)'F ',0,0,-1.d0,-1.d0,-1.d0,-1.d0,-1.d0,0,0
C--write out vertex line    
    IF(collider.EQ.'EEJJ')THEN
      WRITE(j,5400)'V ',-1,0,0,0,0,0,2,1,0
      WRITE(j,5500)'P ',1,-11,0.d0,0.d0,sqrts/2.,sqrts/2.,
     &  0.00051,2,0,0,-1,0
      WRITE(j,5500)'P ',2,11,0.d0,0.d0,-sqrts/2.,sqrts/2.,
     &  0.00051,2,0,0,-1,0
      WRITE(j,5500)'P ',3,23,0.d0,0.d0,0.d0,sqrts,
     &  91.2,2,0,0,-2,0
      WRITE(j,5400)'V ',-2,0,0,0,0,0,0,2,0
      WRITE(j,5500)'P ',4,pid,sqrts/2.,0.d0,0.d0,sqrts/2.,
     &  0.000,2,0,0,-3,0
      WRITE(j,5500)'P ',5,-pid,-sqrts/2.,0.d0,0.d0,sqrts/2.,
     &  0.000,2,0,0,-3,0
      WRITE(j,5400)'V ',-3,0,0,0,0,0,0,nfirst,0
        ELSE
      WRITE(j,5400)'V ',-1,0,0,0,0,0,2,nfirst,0
      if (beam1.eq.'p+') then
      WRITE(j,5500)'P ',1,2212,0.d0,0.d0,
     &  sqrt(sqrts**2/4.-mproton**2),sqrts/2.,mproton,2,0,0,-1,0
      else
      WRITE(j,5500)'P ',1,2112,0.d0,0.d0,
     &  sqrt(sqrts**2/4.-mneutron**2),sqrts/2.,mneutron,2,0,0,-1,0
      endif
      if (beam2.eq.'p+') then
        WRITE(j,5500)'P ',2,2212,0.d0,0.d0,
     &  -sqrt(sqrts**2/4.-mproton**2),sqrts/2.,mproton,2,0,0,-1,0
      else
        WRITE(j,5500)'P ',2,2112,0.d0,0.d0,
     &  -sqrt(sqrts**2/4.-mneutron**2),sqrts/2.,mneutron,2,0,0,-1,0
      endif
    ENDIF
C--write out scattering centres
  if(writescatcen) then
      do 133 i=1,nscatcen
        pbarcode=pbarcode+1
        WRITE(j,5500)'P ',pbarcode,scatflav(i),scatcen(i,1),
     &    scatcen(i,2),scatcen(i,3),scatcen(i,4),scatcen(i,5),
     &    3,0,0,0,0
 133      continue
    endif   
C--write out dummy particles
    if(writedummies) then
      do 135 i=1,nscatcen
        pbarcode=pbarcode+1
        pscatcen=sqrt(scatcen(i,1)**2+scatcen(i,2)**2+
     &    scatcen(i,3)**2)
        WRITE(j,5500)'P ',pbarcode,111,pdummy*scatcen(i,1)/pscatcen,
     &    pdummy*scatcen(i,2)/pscatcen,pdummy*scatcen(i,3)/pscatcen,
     &    pdummy,0.d0,1,0,0,0,0
 135      continue
    endif   
C--write out particle lines
    do 132 i=1,n
      if(((k(i,1).lt.6).or.(k(i,1).eq.17))) then
        pbarcode=pbarcode+1
    if((k(i,1).eq.3).or.(k(i,1).eq.5)) then
          WRITE(j,5500)'P ',pbarcode,k(i,2),p(i,1),p(i,2),p(i,3),
     &    p(i,4),p(i,5),4,0,0,0,0
        else
          WRITE(j,5500)'P ',pbarcode,k(i,2),p(i,1),p(i,2),p(i,3),
     &    p(i,4),p(i,5),1,0,0,0,0
    endif
      endif
 132    continue

  else
C--long output
    if (hadro) then
C--hadronised events
      nfirst=0
          IF(collider.EQ.'EEJJ')THEN
            nvertex=3
          ELSE
            nvertex=1
          ENDIF
      DO 123 i=1,n
        IF(k(i,3).ne.0)THEN
          nstart=i
          goto 124
        ENDIF
 123      CONTINUE   
 124      CONTINUE   
      nstart=0

          DO 126 i=nstart+1,n
        IF(isprimhadron(i)) nfirst=nfirst+1
        IF((ishadron(k(i,2)).OR.(abs(k(i,2)).EQ.15))
     &    .AND.(k(i,4).NE.0)) nvertex=nvertex+1
 126      CONTINUE   
 127      CONTINUE   

      if(writescatcen) nfirst=nfirst+nscatcen
      if(writedummies) nfirst=nfirst+nscatcen

      WRITE(j,5000)'E ',evnum,-1,0.d0,0.d0,0.d0,0,0,nvertex,
     &1,2,0,1,pari(10)
      WRITE(j,'(A2,I2,A5)')'N ',1,'"0"' 
      WRITE(j,'(A)')'U GEV MM'
      WRITE(j,5100)'C ',pari(1)*1.d9,0.d0
      WRITE(j,5200)'H ',0,0,0,0,0,0,0,0,0,0.d0,0.d0,0.d0,0.d0
      WRITE(j,5300)'F ',0,0,-1.d0,-1.d0,-1.d0,-1.d0,-1.d0,0,0

C--write out vertex line    
          IF(collider.EQ.'EEJJ')THEN
        vbarcode=-3
        pbarcode=5
      ELSE
        vbarcode=-1
        pbarcode=2
      ENDIF
      IF(collider.EQ.'EEJJ')THEN
        WRITE(j,5400)'V ',-1,0,0,0,0,0,2,1,0
        WRITE(j,5500)'P ',1,-11,0.d0,0.d0,sqrts/2.,sqrts/2.,
     &  0.00051,2,0,0,-1,0
        WRITE(j,5500)'P ',2,11,0.d0,0.d0,-sqrts/2.,sqrts/2.,
     &  0.00051,2,0,0,-1,0
        WRITE(j,5500)'P ',3,23,0.d0,0.d0,0.d0,sqrts,
     &  91.2,2,0,0,-2,0
        WRITE(j,5400)'V ',-2,0,0,0,0,0,0,2,0
        WRITE(j,5500)'P ',4,pid,sqrts/2.,0.d0,0.d0,sqrts/2.,
     &  0.000,2,0,0,-3,0
    WRITE(j,5500)'P ',5,-pid,-sqrts/2.,0.d0,0.d0,sqrts/2.,
     &  0.000,2,0,0,-3,0
    WRITE(j,5400)'V ',vbarcode,0,0,0,0,0,0,nfirst,0
          ELSE
        WRITE(j,5400)'V ',-1,0,0,0,0,0,2,nfirst,0
      if (beam1.eq.'p+') then
      WRITE(j,5500)'P ',1,2212,0.d0,0.d0,
     &  sqrt(sqrts**2/4.-mproton**2),sqrts/2.,mproton,2,0,0,-1,0
      else
      WRITE(j,5500)'P ',1,2112,0.d0,0.d0,
     &  sqrt(sqrts**2/4.-mneutron**2),sqrts/2.,mneutron,2,0,0,-1,0
      endif
      if (beam2.eq.'p+') then
        WRITE(j,5500)'P ',2,2212,0.d0,0.d0,
     &  -sqrt(sqrts**2/4.-mproton**2),sqrts/2.,mproton,2,0,0,-1,0
      else
        WRITE(j,5500)'P ',2,2112,0.d0,0.d0,
     &  -sqrt(sqrts**2/4.-mneutron**2),sqrts/2.,mneutron,2,0,0,-1,0
      endif
      ENDIF
       
      codefirst=nfirst+pbarcode

C--write out scattering centres
    if(writescatcen) then
      do 134 i=1,nscatcen
        pbarcode=pbarcode+1
        WRITE(j,5500)'P ',pbarcode,scatflav(i),scatcen(i,1),
     &    scatcen(i,2),scatcen(i,3),scatcen(i,4),scatcen(i,5),
     &    3,0,0,0,0
 134      continue
    endif   
C--write out dummy particles
    if(writedummies) then
      do 136 i=1,nscatcen
        pbarcode=pbarcode+1
        pscatcen=sqrt(scatcen(i,1)**2+scatcen(i,2)**2+
     &    scatcen(i,3)**2)
        WRITE(j,5500)'P ',pbarcode,111,pdummy*scatcen(i,1)/pscatcen,
     &    pdummy*scatcen(i,2)/pscatcen,pdummy*scatcen(i,3)/pscatcen,
     &    pdummy,0.d0,1,0,0,0,0
 136      continue
    endif   

C--first write out all particles coming directly from string or cluster decays
       DO 125 i=nstart+1,n
         IF(.not.isprimhadron(i))THEN
           goto 125
         ELSE
           IF (pbarcode.EQ.codefirst) goto 130
           pbarcode=pbarcode+1
C--write out particle line    
           IF(k(i,4).GT.0)THEN
             vbarcode=vbarcode-1
             codelist(i)=vbarcode
            WRITE(j,5500)'P ',pbarcode,k(i,2),p(i,1),p(i,2),p(i,3),
     &       p(i,4),p(i,5),2,0,0,vbarcode,0
           ELSE 
            WRITE(j,5500)'P ',pbarcode,k(i,2),p(i,1),p(i,2),p(i,3),
     &       p(i,4),p(i,5),1,0,0,0,0
           ENDIF      
         ENDIF   
 125       CONTINUE    
 130       CONTINUE 
C--now write out all other particles and vertices 
       DO 129 i=nstart+1,n
         if (isprimhadron(i).or.isprimstring(i)) goto 129
         if (isparton(k(i,2))) then
           if (ishadron(k(k(i,3),2))) codelist(i)=codelist(k(i,3))
           goto 129
         endif
         if (issecstring(i)) then
           codelist(i)=codelist(k(i,3))
           goto 129
         endif
         pbarcode=pbarcode+1
         IF((k(i,3).NE.k(i-1,3)))THEN
C--write out vertex line    
           WRITE(j,5400)'V ',codelist(k(i,3)),0,0,0,0,0,0,
     &        k(k(i,3),5)-k(k(i,3),4)+1,0
         ENDIF 
C--write out particle line    
         IF(k(i,4).GT.0)THEN
           vbarcode=vbarcode-1
           codelist(i)=vbarcode
           WRITE(j,5500)'P ',pbarcode,k(i,2),p(i,1),p(i,2),p(i,3),
     &    p(i,4),p(i,5),2,0,0,vbarcode,0
         ELSE 
           WRITE(j,5500)'P ',pbarcode,k(i,2),p(i,1),p(i,2),p(i,3),
     &    p(i,4),p(i,5),1,0,0,0,0
         ENDIF      
 129       CONTINUE

    else
C--partonic events
    endif
  endif
  call flush(j)
  END
  


***********************************************************************
***   subroutine printlogo
***********************************************************************
  subroutine printlogo(fid)
  implicit none
  integer fid

  write(fid,*)
  write(fid,*)'                   _______________'//
     &'__________________________                  '
  write(fid,*)'                  |               '//
     &'                          |                 '
  write(fid,*)'                  |  JJJJJ  EEEEE '//
     &' W       W  EEEEE  L      |                  '
  write(fid,*)'                  |      J  E     '//
     &' W       W  E      L      |                  '
  write(fid,*)' _________________|      J  EEE   '//
     &'  W  W  W   EEE    L      |_________________ '
  write(fid,*)'|                 |  J   J  E     '//
     &'  W W W W   E      L      |                 |'
  write(fid,*)'|                 |   JJJ   EEEEE '//
     &'   W   W    EEEEE  LLLLL  |                 |'
  write(fid,*)'|                 |_______________'//
     &'__________________________|                 |'
  write(fid,*)'|                                 '//
     &'                                            |'
  write(fid,*)'|                            '//
     &'this is JEWEL 2.1.0                              |'
  write(fid,*)'|                                 '//
     &'                                            |'
  write(fid,*)'| Copyright Korinna C. Zapp (2016)'//
     &'  [Korinna.Zapp@cern.ch]                    |'
  write(fid,*)'|                                 '//
     &'                                            |'
  write(fid,*)'| The JEWEL homepage is jewel.hepforge.org '//
     &'                                   |'
  write(fid,*)'|                                 '//
     &'                                            |'
  write(fid,*)'| The medium model was partly '//
     &'implemented by Jochen Klein                     |'
  write(fid,*)'| [Jochen.Klein@cern.ch]. Raghav '//
     &'Kunnawalkam Elayavalli helped with the       |'
  write(fid,*)'| implementation of the V+jet processes '//
     &'[raghav.k.e@cern.ch].                 |'
  write(fid,*)'|                                 '//
     &'                                            |'
  write(fid,*)'| Please cite JHEP 1303 (2013) '//
     &'080 [arXiv:1212.1599] and optionally           |'
  write(fid,*)'| EPJC C60 (2009) 617 [arXiv:0804.3568] '//
     &'for the physics and arXiv:1311.0048   |'
  write(fid,*)'| for the code. The reference for '//
     &'V+jet processes is EPJC 76 (2016) no.12 695 |'
       write(fid,*)'| [arXiv:1608.03099] and for recoil effects'//
     &' it is arXiv:1707.01539.           |'
  write(fid,*)'|                                 '//
     &'                                            |'
  write(fid,*)'| JEWEL contains code provided by '//
     &'S. Zhang and J. M. Jing                     |'
  write(fid,*)'| (Computation of Special Functions, '//
     &'John Wiley & Sons, New York, 1996 and    |'
  write(fid,*)'| http://jin.ece.illinois.edu) for '//
     &'computing the exponential integral Ei(x).  |'
  write(fid,*)'|                                 '//
     &'                                            |'
  write(fid,*)'| JEWEL relies heavily on PYTHIA 6'//
     &' for the event generation. The modified     |'
  write(fid,*)'| version of PYTHIA 6.4.25 that is'//
     &' shipped with JEWEL is, however, not an     |'
  write(fid,*)'| official PYTHIA release and must'//
     &' not be used for anything else. Please      |'
  write(fid,*)'| refer to results as "JEWEL+PYTHIA".'//
     &'                                         |'
  write(fid,*)'|                                 '//
     &'                                            |'
  write(fid,*)'|_________________________________'//
     &'____________________________________________|'
  write(fid,*)
  write(fid,*)
  end


***********************************************************************
***   subroutine printtime
***********************************************************************
  subroutine printtime
  implicit none
C--identifier of file for hepmc output and logfile
  common/hepmcid/hpmcfid,logfid
  integer hpmcfid,logfid
C--local variables
  integer*4 date(3),time(3)

 1000 format (i2.2, '.', i2.2, '.', i4.4, ', ',
     &         i2.2, ':', i2.2, ':', i2.2 )
  call idate(date)
  call itime(time)
  write(logfid,1000)date,time
  end