lepto.F

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

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

      SUBROUTINE lepto

      IMPLICIT NONE

C...Administer the generation of an event.
C...Note: if error flag LST(21) is non-zero, no proper event generated.

      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/

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

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


* if ARIADNE is used to simulate the parton shower evolution, the
* ARDAT1 Common block is neccessary for a proper interface.
* if ARIADNE is used to simulate the parton shower evolution, the
* ARDAT1 Common block is neccessary for a proper interface.
      COMMON /ardat1/ para(40),msta(40)
      REAL para
      INTEGER msta
      SAVE /ardat1/ 


      INTEGER nummis,nwarn,i,j,ns,l17
      REAL ulmass,ulangl,rlu,ulalps
      REAL qg,qqb,spq,srlu,plu,pt,phi,energy,p2,atan
      DOUBLE PRECISION detot,dari29,dari30
      dimension spq(17)
      DATA nummis,nwarn/0,10/,dari29,dari30/2*0.d0/

      l17=0
    1 lst(21)=0
      DO 10 i=1,10
      DO 10 j=1,5
      k(i,j)=0
   10 v(i,j)=0.
      DO 15 i=1,4
      k(i,1)=21
   15 k(i,2)=ksave(i)
      k(4,1)=1
      n=2

      IF(lst(17).NE.0.AND.lst(2).GT.0) THEN
C...Lepton and/or nucleon energy may vary from event to event,
        IF(l17.EQ.0) THEN
C...Momentum vectors from P(i,j) i=1,2 j=1,2,3 on entry in LEPTO
          DO 20 i=1,2
          p(i,5)=ulmass(k(i,2))
          p(i,4)=sqrt(p(i,1)**2+p(i,2)**2+p(i,3)**2+p(i,5)**2)
          DO 20 j=1,5
   20     psave(3,i,j)=p(i,j)
        ELSE
C...Momentum vectors from PSAVE if new try, i.e. jump back to 1
          DO 25 i=1,2
          DO 25 j=1,5
   25     p(i,j)=psave(3,i,j)
        ENDIF
        l17=1
C...Transform to cms of incoming particles, lepton along +z axis.
        DO 30 j=1,3
   30   dbeta(1,j)=(dble(p(1,j))+dble(p(2,j)))/
     &             (dble(p(1,4))+dble(p(2,4)))
        CALL ludbrb(0,0,0.,0.,-dbeta(1,1),-dbeta(1,2),-dbeta(1,3))
        sphi(1)=ulangl(p(1,1),p(1,2))
        CALL ludbrb(0,0,0.,-sphi(1),0.d0,0.d0,0.d0)
        stheta(1)=ulangl(p(1,3),p(1,1))
        CALL ludbrb(0,0,-stheta(1),0.,0.d0,0.d0,0.d0)
        lst(28)=2
        parl(21)=2.*(p(1,4)*p(2,4)-p(1,3)*p(2,3))
      ELSE
C...Initial state momenta fixed from LINIT call.
        DO 42 i=1,2
        DO 40 j=1,5
   40   p(i,j)=psave(3,i,j)
   42   IF(psave(3,1,3).LT.0.) p(i,3)=-psave(3,i,3)
        lst(28)=3
      ENDIF
      CALL leptox
C...Return if error or if no event to be generated.
      IF(lst(21).NE.0.OR.lst(2).LE.0.OR.lst(7).EQ.-1) RETURN

      IF(pari(29).LT.0.5) THEN
C...For first call, reset double precision counters.
        dari29=0.d0
        dari30=0.d0
      ENDIF
      dari29=dari29+1.d0
      pari(29)=dari29
C     CALL GULIST(-3,2)
C...Scattered lepton and exchanged boson added to event record in LKINEM
C...Transform to lepton-nucleon cms if not made earlier
      IF(lst(17).EQ.0) THEN
        DO 46 i=3,4
        DO 45 j=1,5
   45   psave(3,i,j)=p(i,j)
   46   IF(psave(3,1,3).LT.0.) psave(3,i,3)=-p(i,3)
        CALL ludbrb(0,0,0.,0.,0.d0,0.d0,-dbeta(1,3))
        lst(28)=2
      ENDIF
      DO 50 i=1,4
      DO 50 j=1,5
   50 psave(2,i,j)=p(i,j)
C     CALL GULIST(-2,2)

C...Prepare for parton cascade.
      IF(lst(8).GE.2.AND.mod(lst(8),10).NE.9) CALL lshowr(0)

C...Transform to hadronic cms, boost parameters in double precision.
      detot=dble(p(1,4))-dble(p(4,4))+dble(p(2,4))
      dbeta(2,1)=-dble(p(4,1))/detot
      dbeta(2,2)=-dble(p(4,2))/detot
      dbeta(2,3)=(dble(p(1,3))-dble(p(4,3))+dble(p(2,3)))/detot
      CALL ludbrb(0,0,0.,0.,-dbeta(2,1),-dbeta(2,2),-dbeta(2,3))
      sphi(2)=0.
      stheta(2)=ulangl(p(3,3),p(3,1))
      CALL ludbrb(0,0,-stheta(2),0.,0.d0,0.d0,0.d0)
      lst(28)=1
      DO 60 i=1,4
      DO 60 j=1,5
   60 psave(1,i,j)=p(i,j)
C...Save momentum of exchanged boson (used in subroutine LFRAME).
      DO 70 j=1,5
   70 pb(j)=p(3,j)
C     CALL GULIST(-1,2)

   90 n=4
      mstu(1)=n+1
      lst(26)=n+1
      lst(27)=0
      parl(25)=ulalps(q2)
      IF(lst(8).EQ.1.OR.lst(8)/10.EQ.1.OR.mod(lst(8),10).EQ.9) THEN
C...Probabilities for hard, first order QCD events.
CAE...Corrected what to do when LQGEV or LQQBEV fail. Now make LQEV.
         CALL lqcdpr(qg,qqb)

        DO 100 i=1,17
  100   spq(i)=pq(i)
  200   srlu=rlu(0)

        IF(srlu.GT.qqb+qg) THEN
          CALL lqev
        ELSEIF(srlu.GT.qqb) THEN
          IF(lst(8).EQ.9) THEN
            DO 211 i=1,17
  211         pq(i)=spq(i)
            CALL lqev
          ELSE
            CALL lqgev
            DO 212 i=1,17
              pq(i)=spq(i)
  212       CONTINUE
          ENDIF
        ELSE
          CALL lqqbev
          DO 213 i=1,17
            pq(i)=spq(i)
  213     CONTINUE
          IF(lst(8).EQ.9.AND.lst(21).EQ.0) THEN
            IF(plu(5,11).LT.q2*para(20)) THEN
              DO 220 i=1,17
  220         pq(i)=spq(i)
              CALL lqevar(k(5,2),k(7,2))
            ENDIF
          ENDIF
        ENDIF
        IF(lst(21).NE.0) THEN
  230     CALL lqev
          IF(lst(21).NE.0) goto 230
        ENDIF   
      ELSE
C...QPM model without QCD corrections (cascade applied later).
  300   CALL lqev
        IF(lst(21).NE.0) goto 300
      ENDIF

      ns=mstu(1)
      mstu(1)=0
C     CALL GULIST(-3,2)
C     WRITE(6,*) ' LST(24)=',LST(24)
CJR--       no preclustering of small systems
          mstj(14)=-1
CJR--            
      IF(lst(8).LE.1.OR.mod(lst(8),10).EQ.9) THEN
C...No parton cascade, introduce primordial kt.
        IF(parl(3).GT.1.e-03) THEN
          CALL lprikt(parl(3),pt,phi)
          CALL ludbrb(ns,n,0.,-phi,0.d0,0.d0,0.d0)
          CALL ludbrb(ns,n,atan(2.*pt/sqrt(w2)),phi,0.d0,0.d0,0.d0)
        ENDIF
        IF(mod(lst(8),10).NE.9) THEN
C...Check system against fragmentation cuts.
          mstu(24)=0
          CALL luprep(0)
          IF(mstu(24).NE.0) THEN
            IF(lst(3).GE.1) WRITE(6,*)'LUPREP error MSTU(24)=',mstu(24),
     &                                ', New event generated'
            lst(21)=1
            goto 1
          ENDIF
        ENDIF
      ELSEIF(lst(24).EQ.1) THEN
C...Include parton cascades (+ remnant & kt) on q-event
        CALL lshowr(1)
      ELSE
C...Include parton cascades (+ remnant & kt) on qg- or qqbar-event
        CALL lmeps
      ENDIF
      IF(lst(21).NE.0) THEN
C        IF(LST(3).GE.1)
C     &     WRITE(6,*)'Cascade error LST(21)= ',LST(21),
C     &               ', New event generated'
        goto 1
      ENDIF

CJR--       Soft colour interactions
      IF(lst(34).EQ.1) CALL lsci(parl(7))
      IF(lst(21).NE.0) goto 1
CJR--       take care of small systems
      CALL lsmall
      IF(lst(21).NE.0) THEN
        IF(lst(3).GE.1) WRITE(6,*)' LSMALL error LST(21)= ',lst(21),
     &                            ', New event generated'
        goto 1
      ENDIF
      mstj(14)=1
      CALL luprep(0)
      IF(mstu(24).NE.0) THEN
         IF(lst(3).GE.1) WRITE(6,*)' LUPREP error MSTU(24)= ',mstu(24),
     &                             ', New event generated'
         lst(21)=1
      ENDIF
CJR--            
      IF(lst(21).NE.0) goto 1

      DO 400 i=1,n
C...Correct energy-momentum-mass mismatch for real particle
      IF(p(i,5).LT.0.) goto 400
      energy=sqrt(dble(p(i,5))**2+dble(p(i,1))**2+dble(p(i,2))**2+
     &dble(p(i,3))**2)
      p2=dble(p(i,4))**2-dble(p(i,1))**2-dble(p(i,2))**2-dble(p(i,3))**2
      IF(abs(energy-p(i,4))/(psave(3,1,4)+psave(3,2,4)).GT.paru(11))THEN
        nummis=nummis+1
C...For testing purposes
C       IF(LST(3).GE.1.AND.NUMMIS.LE.NWARN) THEN
C         WRITE(6,1000) I,(K(I,J),J=1,2),(P(I,J),J=1,5),
C    &    SIGN(SQRT(ABS(P2)),P2),ENERGY,INT(DARI29),NWARN
C         IF(ABS(P2-P(I,5)**2).GT.400.) CALL LULIST(2)
C       ENDIF
CAE        WRITE(6,*) 'Energy mismatch',LST(24),PARL(28),PARL(29),NUMMIS
        goto 90
      ENDIF
      p(i,4)=energy
  400 CONTINUE

      dari30=dari30+1.d0
      pari(30)=dari30
Ctest IF(LST(23).EQ.2) PARL(24)=PARL(24)*DARI30/DARI29

      DO 500 i=1,n
      DO 500 j=1,5
  500 v(i,j)=0.
      IF(lst(7).EQ.1) THEN
        CALL luexec
        IF(mstu(24).NE.0) THEN
          WRITE(6,*) ' Error from JETSET, new event made'
          goto 90
        ENDIF
      ENDIF

C     CALL GULIST(-1,2)
C...Transform to desired frame
C     LST(28)=1
      lst(29)=0
      phir=6.2832*rlu(0)
      IF(lst(17).EQ.0) THEN
        IF(lst(5).GE.2) CALL lframe(lst(5),0)
C...Restore momenta (e,p,boson,L) due to numerical errors from boosts
C...But only in frames 1-3
        IF(lst(28).GT.0 .AND. lst(28).LT.4) THEN
           DO 600 i=1,4
           DO 600 j=1,5
  600      p(i,j)=psave(lst(28),i,j)
        ENDIF 
        IF(lst(6).EQ.1.AND.lst(28).GE.2) THEN
C...Random rotation in azimuthal angle
          CALL ludbrb(0,0,0.,phir,0.d0,0.d0,0.d0)
          lst(29)=1
        ENDIF
      ELSE
        IF(lst(5).GE.2) CALL lframe(lst(5),lst(6))
      ENDIF
C...Deactivate scattered lepton
      IF(mod(lst(4),10).EQ.0) k(4,1)=21
C     CALL GULIST(0,2)

      RETURN
 1000 FORMAT(' Warning: too large numerical mismatch in ',
     &'particle energy-momentum-mass',
     &/,3x,'I K(I,1) ..2)  P(I,1)  P(I,2)  P(I,3)',
     &'  P(I,4)  P(I,5)    mass  energy',/,i4,2i6,7f8.3,/,
     &' Event no.',i8,' regenerated. Only first',i5,' warnings printed')
      END