lkinem.F

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

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

      FUNCTION lkinem(L)
      
      IMPLICIT NONE
      
C...Calculate kinematical variables and reject (optionally) if outside
C...required limits.

      COMMON /lintrl/ psave(3,4,5),ksave(4),xmin,xmax,ymin,ymax,
     &q2min,q2max,w2min,w2max,ilep,inu,ig,iz
      REAL psave,xmin,xmax,ymin,ymax,q2min,q2max,w2min,w2max
      INTEGER ksave,ilep,inu,ig,iz
      SAVE /lintrl/

      INTEGER  nlupdm,nplbuf
      parameter(nlupdm=4000,nplbuf=5)
      common/lujets/n,k(nlupdm,5),p(nlupdm,nplbuf),v(nlupdm,5)
      INTEGER  n,k
      REAL  p,v
      SAVE /lujets/

      common/ludat1/mstu(200),paru(200),mstj(200),parj(200)
      INTEGER  mstu,mstj
      REAL  paru,parj
      SAVE /ludat1/

*
* 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/

      COMMON /lboost/ dbeta(2,3),stheta(2),sphi(2),pb(5),phir
      DOUBLE PRECISION dbeta
      REAL stheta,sphi,pb,phir
      SAVE /lboost/


      INTEGER lkinem,l,j
      REAL thetal,plu
      REAL ulmass
      DOUBLE PRECISION de,dpz,dpt

      lkinem=1
      IF(l.EQ.-3) THEN
C...x,W known from LWEITS, no cuts applied.
        u=(w2-p(2,5)**2)/(2.*p(2,5)*(1.-x))
        q2=2.*p(2,5)*u*x
        y=q2/(parl(21)*x)
        goto 200
      ENDIF
C...x,y given.
      parl(22)=y*parl(21)
      q2=x*parl(22)
      u=parl(22)/(2.*p(2,5))
      w2=parl(22)*(1.-x)+p(2,5)**2
      p(4,5)=ulmass(k(4,2))
      IF(p(4,5)/sqrt(parl(21)).LT.0.001) THEN
C...Simpler formulae for effectively massless scattered lepton.
        de=dble(p(1,4))*(1.-dble(y))+dble(x)*dble(y)*dble(abs(p(2,3)))
        dpz=de-dble(x)*dble(y)*(dble(p(2,4))+dble(abs(p(2,3))))
      ELSE
C...Formulae for massive scattered lepton.
        de=dble(p(1,4))+(dble(abs(p(2,3)))*(dble(q2)+dble(p(1,5))**2+
     &  dble(p(4,5))**2)/(2.d0*dble(p(1,4)))-dble(parl(22))/2.d0)/
     &  (dble(p(2,4))+dble(abs(p(2,3))))
        dpz=dble(p(1,4))-(dble(p(2,4))*(dble(q2)+dble(p(1,5))**2+
     &  dble(p(4,5))**2)/(2.d0*dble(p(1,4)))+dble(parl(22))/2.d0)/
     &  (dble(p(2,4))+dble(abs(p(2,3))))
      ENDIF
      dpt=de**2-dpz**2-dble(p(4,5))**2
      IF(dpt.LT.0.d0) RETURN
      dpt=sqrt(dpt)
      p(4,1)=dpt
      p(4,2)=0.
      p(4,3)=dpz
      p(4,4)=de
      p(3,1)=-dpt
      p(3,2)=0.
      p(3,3)=dble(p(1,3))-dpz
      p(3,4)=dble(p(1,4))-de
      p(3,5)=-sqrt(q2)
      k(3,3)=1
      k(4,3)=1
      n=4
      IF(l.EQ.3) goto 200

      IF(x.LT.xmin.OR.x.GT.xmax) RETURN
      IF(y.LT.ymin.OR.y.GT.ymax) RETURN
      IF(q2.LT.q2min.OR.q2.GT.q2max) RETURN
      IF(w2.LT.w2min.OR.w2.GT.w2max) RETURN
C-check: CUT(9),CUT(10) --> UMIN,UMAX needs change in /LINTRL/ --> next update 
      IF(u.LT.cut(9).OR.u.GT.cut(10)) RETURN
      IF(lst(17).EQ.0) THEN
        IF(p(4,4).LT.cut(11).OR.p(4,4).GT.cut(12))  RETURN
        thetal=plu(4,13)
C       THETAL=ACOS((P(1,1)*P(4,1)+P(1,2)*P(4,2)+P(1,3)*P(4,3))
C    &  /SQRT(P(1,1)**2+P(1,2)**2+P(1,3)**2)/
C    &  SQRT(P(4,1)**2+P(4,2)**2+P(4,3)**2))
      ELSE
C...No cuts on energy, angle for initialisation of varying energy mode
        IF(lst(32).NE.0) goto 200
C...Transform scattered lepton back to lab system to make cut
C...in energy and angle (defined as space angle to incoming lepton).
        DO 110 j=1,5
        k(6,j)=k(4,j)
  110   p(6,j)=p(4,j)
        CALL ludbrb(6,6,stheta(1),sphi(1),0.d0,0.d0,0.d0)
        CALL ludbrb(6,6,0.,0.,dbeta(1,1),dbeta(1,2),dbeta(1,3))
        IF(p(6,4).LT.cut(11).OR.p(6,4).GT.cut(12))  RETURN
        thetal=acos((psave(3,1,1)*p(6,1)+psave(3,1,2)*p(6,2)+
     &  psave(3,1,3)*p(6,3))
     &  /sqrt(psave(3,1,1)**2+psave(3,1,2)**2+psave(3,1,3)**2)/
     &  sqrt(p(6,1)**2+p(6,2)**2+p(6,3)**2))
      ENDIF
      IF(thetal.LT.cut(13).OR.thetal.GT.cut(14))  RETURN
  200 lkinem=0
      RETURN
      END