pykcut.f

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

 
C*********************************************************************
 
C...PYKCUT
C...Dummy routine, which the user can replace in order to make cuts on
C...the kinematics on the parton level before the matrix elements are
C...evaluated and the event is generated. The cross-section estimates
C...will automatically take these cuts into account, so the given
C...values are for the allowed phase space region only. MCUT=0 means
C...that the event has passed the cuts, MCUT=1 that it has failed.
 
      SUBROUTINE pykcut(MCUT)
 
C...Double precision and integer declarations.
      IMPLICIT DOUBLE PRECISION(a-h, o-z)
      IMPLICIT INTEGER(i-n)
      INTEGER pyk,pychge,pycomp
C...Commonblocks.
      common/pydat1/mstu(200),paru(200),mstj(200),parj(200)
      common/pyint1/mint(400),vint(400)
      common/pyint2/iset(500),kfpr(500,2),coef(500,20),icol(40,4,2)
      SAVE /pydat1/,/pyint1/,/pyint2/
 
C...Set default value (accepting event) for MCUT.
      mcut=0
 
C...Read out subprocess number.
      isub=mint(1)
      istsb=iset(isub)
 
C...Read out tau, y*, cos(theta), tau' (where defined, else =0).
      tau=vint(21)
      yst=vint(22)
      cth=0d0
      IF(istsb.EQ.2.OR.istsb.EQ.4) cth=vint(23)
      taup=0d0
      IF(istsb.GE.3.AND.istsb.LE.5) taup=vint(26)
 
C...Calculate x_1, x_2, x_F.
      IF(istsb.LE.2.OR.istsb.GE.5) THEN
        x1=sqrt(tau)*exp(yst)
        x2=sqrt(tau)*exp(-yst)
      ELSE
        x1=sqrt(taup)*exp(yst)
        x2=sqrt(taup)*exp(-yst)
      ENDIF
      xf=x1-x2
 
C...Calculate shat, that, uhat, p_T^2.
      shat=tau*vint(2)
      sqm3=vint(63)
      sqm4=vint(64)
      rm3=sqm3/shat
      rm4=sqm4/shat
      be34=sqrt(max(0d0,(1d0-rm3-rm4)**2-4d0*rm3*rm4))
      rpts=4d0*vint(71)**2/shat
      be34l=sqrt(max(0d0,(1d0-rm3-rm4)**2-4d0*rm3*rm4-rpts))
      rm34=2d0*rm3*rm4
      rsqm=1d0+rm34
      rthm=(4d0*rm3*rm4+rpts)/(1d0-rm3-rm4+be34l)
      that=-0.5d0*shat*max(rthm,1d0-rm3-rm4-be34*cth)
      uhat=-0.5d0*shat*max(rthm,1d0-rm3-rm4+be34*cth)
      pt2=max(vint(71)**2,0.25d0*shat*be34**2*(1d0-cth**2))
 
C...Decisions by user to be put here.
 
C...Stop program if this routine is ever called.
C...You should not copy these lines to your own routine.
      WRITE(mstu(11),5000)
      CALL pystop(6)
 
C...Format for error printout.
 5000 FORMAT(1x,'Error: you did not link your PYKCUT routine ',
     &'correctly.'/1x,'Dummy routine in PYTHIA file called instead.'/
     &1x,'Execution stopped!')
 
      RETURN
      END