pydocu.f

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

 
C*********************************************************************
 
C...PYDOCU
C...Handles the documentation of the process in MSTI and PARI,
C...and also computes cross-sections based on accumulated statistics.
 
      SUBROUTINE pydocu
 
C...Double precision and integer declarations.
      IMPLICIT DOUBLE PRECISION(a-h, o-z)
      IMPLICIT INTEGER(i-n)
      INTEGER pyk,pychge,pycomp
C...Commonblocks.
      common/pyjets/n,npad,k(4000,5),p(4000,5),v(4000,5)
      common/pydat1/mstu(200),paru(200),mstj(200),parj(200)
      common/pypars/mstp(200),parp(200),msti(200),pari(200)
      common/pysubs/msel,mselpd,msub(500),kfin(2,-40:40),ckin(200)
      common/pyint1/mint(400),vint(400)
      common/pyint2/iset(500),kfpr(500,2),coef(500,20),icol(40,4,2)
      common/pyint5/ngenpd,ngen(0:500,3),xsec(0:500,3)
      SAVE /pyjets/,/pydat1/,/pysubs/,/pypars/,/pyint1/,/pyint2/,
     &/pyint5/
 
C...Calculate Monte Carlo estimates of cross-sections.
      isub=mint(1)
      IF(mstp(111).NE.-1) ngen(isub,3)=ngen(isub,3)+1
      ngen(0,3)=ngen(0,3)+1
      xsec(0,3)=0d0
      DO 100 i=1,500
        IF(i.EQ.96.OR.i.EQ.97) THEN
          xsec(i,3)=0d0
        ELSEIF(msub(95).EQ.1.AND.(i.EQ.11.OR.i.EQ.12.OR.i.EQ.13.OR.
     &    i.EQ.28.OR.i.EQ.53.OR.i.EQ.68)) THEN
          xsec(i,3)=xsec(96,2)*ngen(i,3)/max(1d0,dble(ngen(96,1))*
     &    dble(ngen(96,2)))
        ELSEIF(msub(95).EQ.1.AND.i.GE.381.AND.i.LE.386) THEN
          xsec(i,3)=xsec(96,2)*ngen(i,3)/max(1d0,dble(ngen(96,1))*
     &    dble(ngen(96,2)))
        ELSEIF(msub(i).EQ.0.OR.ngen(i,1).EQ.0) THEN
          xsec(i,3)=0d0
        ELSEIF(ngen(i,2).EQ.0) THEN
          xsec(i,3)=xsec(i,2)*ngen(0,3)/(dble(ngen(i,1))*
     &    dble(ngen(0,2)))
        ELSE
          xsec(i,3)=xsec(i,2)*ngen(i,3)/(dble(ngen(i,1))*
     &    dble(ngen(i,2)))
        ENDIF
        xsec(0,3)=xsec(0,3)+xsec(i,3)
  100 CONTINUE
 
C...Rescale to known low-pT cross-section for standard QCD processes.
      IF(msub(95).EQ.1) THEN
        xsech=xsec(11,3)+xsec(12,3)+xsec(13,3)+xsec(28,3)+xsec(53,3)+
     &  xsec(68,3)+xsec(95,3)
        xsecw=xsec(97,2)/max(1d0,dble(ngen(97,1)))
        IF(xsech.GT.1d-20.AND.xsecw.GT.1d-20) THEN
          fac=xsecw/xsech
          xsec(11,3)=fac*xsec(11,3)
          xsec(12,3)=fac*xsec(12,3)
          xsec(13,3)=fac*xsec(13,3)
          xsec(28,3)=fac*xsec(28,3)
          xsec(53,3)=fac*xsec(53,3)
          xsec(68,3)=fac*xsec(68,3)
          xsec(95,3)=fac*xsec(95,3)
          xsec(0,3)=xsec(0,3)-xsech+xsecw
        ENDIF
      ENDIF
 
C...Save information for gamma-p and gamma-gamma.
      IF(mint(121).GT.1) THEN
        iga=mint(122)
        CALL pysave(2,iga)
        CALL pysave(5,0)
      ENDIF
 
C...Reset information on hard interaction.
      DO 110 j=1,200
        msti(j)=0
        pari(j)=0d0
  110 CONTINUE
 
C...Copy integer valued information from MINT into MSTI.
      DO 120 j=1,32
        msti(j)=mint(j)
  120 CONTINUE
      IF(mint(121).GT.1) msti(9)=mint(122)
 
C...Store cross-section variables in PARI.
      pari(1)=xsec(0,3)
      pari(2)=xsec(0,3)/mint(5)
      pari(7)=vint(97)
      pari(9)=vint(99)
      pari(10)=vint(100)
      vint(98)=vint(98)+vint(100)
      IF(mstp(142).EQ.1) pari(2)=xsec(0,3)/vint(98)
 
C...Store kinematics variables in PARI.
      pari(11)=vint(1)
      pari(12)=vint(2)
      IF(isub.NE.95) THEN
        DO 130 j=13,26
          pari(j)=vint(30+j)
  130   CONTINUE
        pari(29)=vint(39)
        pari(30)=vint(40)
        pari(31)=vint(141)
        pari(32)=vint(142)
        pari(33)=vint(41)
        pari(34)=vint(42)
        pari(35)=pari(33)-pari(34)
        pari(36)=vint(21)
        pari(37)=vint(22)
        pari(38)=vint(26)
        pari(39)=vint(157)
        pari(40)=vint(158)
        pari(41)=vint(23)
        pari(42)=2d0*vint(47)/vint(1)
      ENDIF
 
C...Store information on scattered partons in PARI.
      IF(isub.NE.95.AND.mint(7)*mint(8).NE.0) THEN
        DO 140 is=7,8
          i=mint(is)
          pari(36+is)=p(i,3)/vint(1)
          pari(38+is)=p(i,4)/vint(1)
          pr=max(1d-20,p(i,5)**2+p(i,1)**2+p(i,2)**2)
          pari(40+is)=sign(log(min((sqrt(pr+p(i,3)**2)+abs(p(i,3)))/
     &    sqrt(pr),1d20)),p(i,3))
          pr=max(1d-20,p(i,1)**2+p(i,2)**2)
          pari(42+is)=sign(log(min((sqrt(pr+p(i,3)**2)+abs(p(i,3)))/
     &    sqrt(pr),1d20)),p(i,3))
          pari(44+is)=p(i,3)/sqrt(1d-20+p(i,1)**2+p(i,2)**2+p(i,3)**2)
          pari(46+is)=pyangl(p(i,3),sqrt(p(i,1)**2+p(i,2)**2))
          pari(48+is)=pyangl(p(i,1),p(i,2))
  140   CONTINUE
      ENDIF
 
C...Store sum up transverse and longitudinal momenta.
      pari(65)=2d0*pari(17)
      IF(isub.LE.90.OR.isub.GE.95) THEN
        DO 150 i=mstp(126)+1,n
          IF(k(i,1).LE.0.OR.k(i,1).GT.10) goto 150
          pt=sqrt(p(i,1)**2+p(i,2)**2)
          pari(69)=pari(69)+pt
          IF(i.LE.mint(52)) pari(66)=pari(66)+pt
          IF(i.GT.mint(52).AND.i.LE.mint(53)) pari(68)=pari(68)+pt
  150   CONTINUE
        pari(67)=pari(68)
        pari(71)=vint(151)
        pari(72)=vint(152)
        pari(73)=vint(151)
        pari(74)=vint(152)
      ELSE
        pari(66)=pari(65)
        pari(69)=pari(65)
      ENDIF
 
C...Store various other pieces of information into PARI.
      pari(61)=vint(148)
      pari(75)=vint(155)
      pari(76)=vint(156)
      pari(77)=vint(159)
      pari(78)=vint(160)
      pari(81)=vint(138)
 
C...Store information on lepton -> lepton + gamma in PYGAGA.
      msti(71)=mint(141)
      msti(72)=mint(142)
      pari(101)=vint(301)
      pari(102)=vint(302)
      DO 160 i=103,114
        pari(i)=vint(i+202)
  160 CONTINUE
 
C...Set information for PYTABU.
      IF(iset(isub).EQ.1.OR.iset(isub).EQ.3) THEN
        mstu(161)=mint(21)
        mstu(162)=0
      ELSEIF(iset(isub).EQ.5) THEN
        mstu(161)=mint(23)
        mstu(162)=0
      ELSE
        mstu(161)=mint(21)
        mstu(162)=mint(22)
      ENDIF
 
      RETURN
      END