anchor.diff:  *     generate 'njet' number of hard parton-parton scatterings with PYTHIA 
         IMPLICIT DOUBLE PRECISION(A-H, O-Z)
         double precision npar0,nbco0,npart,nbcol,npart0
  -      double precision gaus
  +      double precision qpgaus
         
         INTEGER PYK,PYCHGE,PYCOMP
         CHARACTER beam*2,targ*2
         external pydata 
  -      external pyp,pyr,pyk,pyquen,shad1,gaus
  +      external pyp,pyr,pyk,pyquen,shad1,qpgaus
         
         common /pyjets/ n,npad,k(4000,5),p(4000,5),v(4000,5)
         common /hyjets/ nhj,nhp,khj(150000,5),phj(150000,5),vhj(150000,5)
  
positive.code: 
  c       write(*,*)'in hyhard tau',tau
  
  	etaLj=gaus(rm,sig) !fm
  
  c       write(*,*)'in hyhard after gauss',tau
  	z=tau*sinh(etaLj)
  
negative.code:       v2=2.d0*u2-1.d0 
        s=v1**2+v2**2
        if(s.gt.1) go to 41
        qpgaus=v1*dsqrt(-2.d0*dlog(s)/s)*sig+x0
        return
        end    
  