lqmcut.F

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

C **********************************************************************

      FUNCTION lqmcut(XP,ZP,AM1,AM2,AM3)

C...Apply cuts, if necessary, on the event configuration
C...obtained from QCD matrix elements.

*
* to avoid variable conflictions, a second keep element is necessary
* with the same common block name (see LPTOU2)
*
      COMMON /leptou/ cut(14),lst(40),parl(30),
     &                x,y,w2,q2,u
      REAL cut,parl,x,y,w2,q2,u
      INTEGER lst
      SAVE /leptou/

      COMMON /linter/ pari(50),ewqc(2,2,8),qc(8),zl(2,4),zq(2,8),pq(17)
      REAL  pari,ewqc,qc,zl,zq,pq
      SAVE /linter/


      INTEGER lqmcut
      REAL xp,zp,am1,am2,am3,s12,s23,s13
      REAL w,x1,x2,x3,pa1,pa2,pa3
      DATA s12,s23,s13/3*0./

      IF(lst(24).EQ.2) THEN
        s12=q2*(1.-xp)/xp
        s23=q2*(xp-x)*(1.-zp)/x/xp+am2**2+am3**2
        s13=q2*(xp-x)*zp/x/xp+am1**2+am3**2
      ELSEIF(lst(24).EQ.3) THEN
        s13=q2*(1.-xp)/xp
        s23=q2*(xp-x)*(1.-zp)/x/xp+am2**2+am3**2
        s12=q2*(xp-x)*zp/x/xp+am1**2+am2**2
        IF(s13.LT.(am1+am3)**2) goto 900
      ENDIF

      w=sqrt(w2)
      x1=1.-(s23-am1**2)/w2
      x3=1.-(s12-am3**2)/w2
      x2=2.-x1-x3
      pari(21)=x1
      pari(22)=x2
      pari(23)=x3
                
      IF(x1.GT.1..OR.x2.GT.1..OR.x3.GT.1.) goto 900
      IF(x1*w/2..LT.am1.OR.x2*w/2..LT.am2.OR.x3*w/2..LT.am3) goto 900
      pa1=sqrt((0.5*x1*w)**2-am1**2)
      pa2=sqrt((0.5*x2*w)**2-am2**2)
      pa3=sqrt((0.5*x3*w)**2-am3**2)
      IF(abs((pa3**2-pa1**2-pa2**2)/(2.*pa1*pa2)).GE.1.) goto 900
      IF(abs((pa2**2-pa1**2-pa3**2)/(2.*pa1*pa3)).GE.1.) goto 900
      lqmcut=0
      RETURN

  900 lqmcut=1
      RETURN
      END