lprwts.F

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

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

      SUBROUTINE lprwts(NSTEP)

      IMPLICIT NONE

C...Prints probabilities for q-, qg- and qqbar-events using the present
C...QCD weights stored in common block LGRID.


      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 /lgrid/ nxx,nww,xx(31),ww(21),pqg(31,21,3),pqqb(31,21,2), 
     &qgmax(31,21,3),qqbmax(31,21,2),ycut(31,21),xtot(31,21),np
      REAL    xx,ww,pqg,pqqb,qgmax,qqbmax,ycut,xtot
      INTEGER nxx,nww,np
      SAVE /lgrid/


      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 nstep,lw,iw,lx,ix,ip
      REAL wmax,w,rqg,rqqb,qtot,rq,ulalps

      wmax=sqrt(parl(21)+psave(3,1,5)**2+psave(3,2,5)**2)
      WRITE(6,1000) parl(11),lst(13),mstu(112),paru(112),
     &parl(8),parl(9),parl(12),parl(13)
      IF(np.EQ.1) THEN
        WRITE(6,1010)
      ELSE
        IF(lst(19).LT.10) WRITE(6,1020)
        IF(lst(19).GE.10) WRITE(6,2020)
      ENDIF
      IF(lst(19).LT.10) THEN
        WRITE(6,1030) lst(19),nxx,nww,xx,ww
        IF(wmax.GT.ww(nww)) WRITE(6,1040) wmax,ww(nww)
        WRITE(6,1100)
      ELSEIF(lst(19).GE.10) THEN
        WRITE(6,2030) lst(19),nxx,nww,xx,ww
        WRITE(6,2100)
      ENDIF

      lw=0
      DO 500 iw=1,nww,max(1,nstep)
      w=ww(iw)
      y=ww(iw)
      IF(lw.GT.0) goto 600
      IF(lst(19).LT.10.AND.w.GT.wmax) lw=lw+1
      w2=w**2
      lx=0
      DO 400 ix=1,nxx,max(1,nstep)
      x=xx(ix)
      IF(lx.GT.0) goto 500
      IF(lst(19).LT.10) THEN
C...x,W2 given.
        u=(w2-psave(3,2,5)**2)/(2.*psave(3,2,5)*(1.-x))
        q2=2.*psave(3,2,5)*u*x
        y=q2/(parl(21)*x)
      ELSEIF(lst(19).GE.10) THEN
C...x,y given.
        parl(22)=y*parl(21)
        q2=x*parl(22)
        u=parl(22)/(2.*psave(3,2,5))
        w2=parl(22)*(1.-x)+psave(3,2,5)**2
      ENDIF
      pari(24)=(1.+(1.-y)**2)/2.
      pari(25)=1.-y
      pari(26)=(1.-(1.-y)**2)/2.
      parl(25)=ulalps(q2)
      IF(y.GT.1.) lx=lx+1
      rqg=0.
      rqqb=0.
      DO 200 ip=1,np
      IF(np.EQ.1) THEN
        rqg=pqg(ix,iw,ip)
        rqqb=pqqb(ix,iw,ip)
      ELSE
        rqg=rqg+pqg(ix,iw,ip)*pari(23+ip)/xtot(ix,iw)
        IF(ip.LT.3) rqqb=rqqb+pqqb(ix,iw,ip)*pari(23+ip)/xtot(ix,iw)
      ENDIF
  200 CONTINUE
C...Include alpha-strong in weight.
      rqg=rqg*parl(25)
      rqqb=rqqb*parl(25)
      IF(lst(33).EQ.-91) THEN
C...Include 3-jet cross section in denominator
        qtot=1.+rqg+rqqb
        rqg =rqg/qtot
        rqqb=rqqb/qtot
      ENDIF
      rq=1.-rqg-rqqb
      IF(lst(19).LT.10) THEN
        WRITE(6,1200) w,x,y,q2,parl(25),ycut(ix,iw),rq,rqg,rqqb
      ELSEIF(lst(19).GE.10) THEN
        WRITE(6,2200) x,y,q2,w,parl(25),ycut(ix,iw),rq,rqg,rqqb
      ENDIF
  400 CONTINUE
  500 CONTINUE
  600 CONTINUE
      RETURN

 1000 FORMAT('1',/,5x,'SUMMARY OF QCD MATRIX ELEMENT INTEGRATION',
     &           /,5x,'-----------------------------------------',//,
     &/,' for gluon radiation (qg-event) and boson-gluon fusion ',
     &'(qq-event) probability.',
     &//,' Required precision in integration, PARL(11) =',f8.4,
     &//,' Heaviest flavour produced in boson-gluon fusion, LST(13) =',
     &i5,//,' Alpha-strong parameters: # flavours, MSTU(112) =',i3,
     &'  QCD lambda, PARU(112) =',f6.3,' GeV',
     &//,' Cuts on matrix elements:',
     &/,' PARL(8), PARL(9), PARL(12), PARL(13) =',4f8.4,/)
 1010 FORMAT(' Lepton energy not allowed to vary in simulation.',/)
 1020 FORMAT(' Lepton energy allowed to vary in simulation, ',/,
     &' y in table below calculated assuming max energy.',/)
 1030 FORMAT(' Grid choice, LST(19) =',i3,5x,'# grid points in x, W =',
     &2i5,/,' x-values in array XX:',/,10f8.5,/,10f8.5,/,11f8.5,
     &    /,' W-values in array WW:',/,10f7.1,/,11f7.1,/)
 1040 FORMAT(' Max W outside grid, execution stopped ! Wmax, grid-max ='
     &,2f12.1)
 1100 FORMAT( 
     &5x,'cut',2x,'q-event',1x,'qg-event',1x,'qq-event',
     &/,1x,77(1h-),/)
 1200 FORMAT(f7.1,2f8.4,1pg12.3,0pf8.2,f8.4,3f9.4)
 2020 FORMAT(' Lepton energy allowed to vary in simulation, ',/,
     &' W in table below calculated assuming max energy.',/)
 2030 FORMAT(' Grid choice, LST(19) =',i3,5x,'# grid points in x, y =',
     &2i5,/,' x-values in array XX:',/,10f8.5,/,10f8.5,/,11f8.5,
     &    /,' y-values in array WW:',/,10f7.4,/,11f7.4,/)
 2100 FORMAT( 
     &5x,'cut',2x,'q-event',1x,'qg-event',1x,'qq-event',
     &/,1x,77(1h-),/)
 2200 FORMAT(2f8.5,1pg12.3,0pf7.1,f8.2,f8.4,3f9.4)
      END