pyofsh.f

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

 
C***********************************************************************
 
C...PYOFSH
C...Calculates partial width and differential cross-section maxima
C...of channels/processes not allowed on mass-shell, and selects
C...masses in such channels/processes.
 
      SUBROUTINE pyofsh(MOFSH,KFMO,KFD1,KFD2,PMMO,RET1,RET2)
 
C...Double precision and integer declarations.
      IMPLICIT DOUBLE PRECISION(a-h, o-z)
      IMPLICIT INTEGER(i-n)
      INTEGER pyk,pychge,pycomp
C...Commonblocks.
      common/pydat1/mstu(200),paru(200),mstj(200),parj(200)
      common/pydat2/kchg(500,4),pmas(500,4),parf(2000),vckm(4,4)
      common/pydat3/mdcy(500,3),mdme(8000,2),brat(8000),kfdp(8000,5)
      common/pysubs/msel,mselpd,msub(500),kfin(2,-40:40),ckin(200)
      common/pypars/mstp(200),parp(200),msti(200),pari(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 /pydat1/,/pydat2/,/pydat3/,/pysubs/,/pypars/,/pyint1/,
     &/pyint2/,/pyint5/
C...Local arrays.
      dimension kfd(2),mbw(2),pmd(2),pgd(2),pmg(2),pml(2),pmu(2),
     &pmh(2),atl(2),atu(2),ath(2),rmg(2),inx1(100),xpt1(100),
     &fpt1(100),inx2(100),xpt2(100),fpt2(100),wdtp(0:400),
     &wdte(0:400,0:5)
 
C...Find if particles equal, maximum mass, matrix elements, etc.
      mint(51)=0
      isub=mint(1)
      kfd(1)=iabs(kfd1)
      kfd(2)=iabs(kfd2)
      meql=0
      IF(kfd(1).EQ.kfd(2)) meql=1
      mlm=0
      IF(mofsh.GE.2.AND.meql.EQ.1) mlm=int(1.5d0+pyr(0))
      IF(mofsh.LE.2.OR.mofsh.EQ.5) THEN
        noff=44
        pmmx=pmmo
      ELSE
        noff=40
        pmmx=vint(1)
        IF(ckin(2).GT.ckin(1)) pmmx=min(ckin(2),vint(1))
      ENDIF
      mmed=0
      IF((kfmo.EQ.25.OR.kfmo.EQ.35.OR.kfmo.EQ.36).AND.meql.EQ.1.AND.
     &(kfd(1).EQ.23.OR.kfd(1).EQ.24)) mmed=1
      IF((kfmo.EQ.32.OR.iabs(kfmo).EQ.34).AND.(kfd(1).EQ.23.OR.
     &kfd(1).EQ.24).AND.(kfd(2).EQ.23.OR.kfd(2).EQ.24)) mmed=2
      IF((kfmo.EQ.32.OR.iabs(kfmo).EQ.34).AND.(kfd(2).EQ.25.OR.
     &kfd(2).EQ.35.OR.kfd(2).EQ.36)) mmed=3
      loop=1
 
C...Find where Breit-Wigners are required, else select discrete masses.
  100 DO 110 i=1,2
        kfca=pycomp(kfd(i))
        IF(kfca.GT.0) THEN
          pmd(i)=pmas(kfca,1)
          pgd(i)=pmas(kfca,2)
        ELSE
          pmd(i)=0d0
          pgd(i)=0d0
        ENDIF
        IF(mstp(42).LE.0.OR.pgd(i).LT.parp(41)) THEN
          mbw(i)=0
          pmg(i)=pmd(i)
          rmg(i)=(pmg(i)/pmmx)**2
        ELSE
          mbw(i)=1
        ENDIF
  110 CONTINUE
 
C...Find allowed mass range and Breit-Wigner parameters.
      DO 120 i=1,2
        IF(mofsh.EQ.1.AND.loop.EQ.1.AND.mbw(i).EQ.1) THEN
          pml(i)=parp(42)
          pmu(i)=pmmx-parp(42)
          IF(mbw(3-i).EQ.0) pmu(i)=min(pmu(i),pmmx-pmd(3-i))
          IF(pmu(i).LT.pml(i)+parj(64)) mbw(i)=-1
        ELSEIF(mbw(i).EQ.1.AND.mofsh.NE.5) THEN
          ilm=i
          IF(mlm.EQ.2) ilm=3-i
          pml(i)=max(ckin(noff+2*ilm-1),parp(42))
          IF(mbw(3-i).EQ.0) THEN
            pmu(i)=pmmx-pmd(3-i)
          ELSE
            pmu(i)=pmmx-max(ckin(noff+5-2*ilm),parp(42))
          ENDIF
          IF(ckin(noff+2*ilm).GT.ckin(noff+2*ilm-1)) pmu(i)=
     &    min(pmu(i),ckin(noff+2*ilm))
          IF(i.EQ.mlm) pmu(i)=min(pmu(i),0.5d0*pmmx)
          IF(meql.EQ.0) pmh(i)=min(pmu(i),0.5d0*pmmx)
          IF(pmu(i).LT.pml(i)+parj(64)) mbw(i)=-1
          IF(mbw(i).EQ.1) THEN
            atl(i)=atan((pml(i)**2-pmd(i)**2)/(pmd(i)*pgd(i)))
            atu(i)=atan((pmu(i)**2-pmd(i)**2)/(pmd(i)*pgd(i)))
            IF(meql.EQ.0) ath(i)=atan((pmh(i)**2-pmd(i)**2)/(pmd(i)*
     &      pgd(i)))
          ENDIF
        ELSEIF(mbw(i).EQ.1.AND.mofsh.EQ.5) THEN
          ilm=i
          IF(mlm.EQ.2) ilm=3-i
          pml(i)=max(ckin(48+i),parp(42))
          pmu(i)=pmmx-max(ckin(51-i),parp(42))
          IF(mbw(3-i).EQ.0) pmu(i)=min(pmu(i),pmmx-pmd(3-i))
          IF(i.EQ.mlm) pmu(i)=min(pmu(i),0.5d0*pmmx)
          IF(meql.EQ.0) pmh(i)=min(pmu(i),0.5d0*pmmx)
          IF(pmu(i).LT.pml(i)+parj(64)) mbw(i)=-1
          IF(mbw(i).EQ.1) THEN
            atl(i)=atan((pml(i)**2-pmd(i)**2)/(pmd(i)*pgd(i)))
            atu(i)=atan((pmu(i)**2-pmd(i)**2)/(pmd(i)*pgd(i)))
            IF(meql.EQ.0) ath(i)=atan((pmh(i)**2-pmd(i)**2)/(pmd(i)*
     &      pgd(i)))
          ENDIF
        ENDIF
  120 CONTINUE
      IF(mbw(1).LT.0.OR.mbw(2).LT.0.OR.(mbw(1).EQ.0.AND.mbw(2).EQ.0))
     &THEN
        CALL pyerrm(3,'(PYOFSH:) no allowed decay product masses')
        mint(51)=1
        RETURN
      ENDIF
 
C...Calculation of partial width of resonance.
      IF(mofsh.EQ.1) THEN
 
C..If only one integration, pick that to be the inner.
        IF(mbw(1).EQ.0) THEN
          pm2=pmd(1)
          pmd(1)=pmd(2)
          pgd(1)=pgd(2)
          pml(1)=pml(2)
          pmu(1)=pmu(2)
        ELSEIF(mbw(2).EQ.0) THEN
          pm2=pmd(2)
        ENDIF
 
C...Start outer loop of integration.
        IF(mbw(1).EQ.1.AND.mbw(2).EQ.1) THEN
          atl2=atan((pml(2)**2-pmd(2)**2)/(pmd(2)*pgd(2)))
          atu2=atan((pmu(2)**2-pmd(2)**2)/(pmd(2)*pgd(2)))
          npt2=1
          xpt2(1)=1d0
          inx2(1)=0
          fmax2=0d0
        ENDIF
  130   IF(mbw(1).EQ.1.AND.mbw(2).EQ.1) THEN
          pm2s=pmd(2)**2+pmd(2)*pgd(2)*tan(atl2+xpt2(npt2)*(atu2-atl2))
          pm2=min(pmu(2),max(pml(2),sqrt(max(0d0,pm2s))))
        ENDIF
        rm2=(pm2/pmmx)**2
 
C...Start inner loop of integration.
        pml1=pml(1)
        pmu1=min(pmu(1),pmmx-pm2)
        IF(meql.EQ.1) pmu1=min(pmu1,pm2)
        atl1=atan((pml1**2-pmd(1)**2)/(pmd(1)*pgd(1)))
        atu1=atan((pmu1**2-pmd(1)**2)/(pmd(1)*pgd(1)))
        IF(pml1+parj(64).GE.pmu1.OR.atl1+1d-7.GE.atu1) THEN
          func2=0d0
          goto 180
        ENDIF
        npt1=1
        xpt1(1)=1d0
        inx1(1)=0
        fmax1=0d0
  140   pm1s=pmd(1)**2+pmd(1)*pgd(1)*tan(atl1+xpt1(npt1)*(atu1-atl1))
        pm1=min(pmu1,max(pml1,sqrt(max(0d0,pm1s))))
        rm1=(pm1/pmmx)**2
 
C...Evaluate function value - inner loop.
        func1=sqrt(max(0d0,(1d0-rm1-rm2)**2-4d0*rm1*rm2))
        IF(mmed.EQ.1) func1=func1*((1d0-rm1-rm2)**2+8d0*rm1*rm2)
        IF(mmed.EQ.2) func1=func1**3*(1d0+10d0*rm1+10d0*rm2+rm1**2+
     &  rm2**2+10d0*rm1*rm2)
        IF(func1.GT.fmax1) fmax1=func1
        fpt1(npt1)=func1
 
C...Go to next position in inner loop.
        IF(npt1.EQ.1) THEN
          npt1=npt1+1
          xpt1(npt1)=0d0
          inx1(npt1)=1
          goto 140
        ELSEIF(npt1.LE.8) THEN
          npt1=npt1+1
          IF(npt1.LE.4.OR.npt1.EQ.6) ish1=1
          ish1=ish1+1
          xpt1(npt1)=0.5d0*(xpt1(ish1)+xpt1(inx1(ish1)))
          inx1(npt1)=inx1(ish1)
          inx1(ish1)=npt1
          goto 140
        ELSEIF(npt1.LT.100) THEN
          isn1=ish1
  150     ish1=ish1+1
          IF(ish1.GT.npt1) ish1=2
          IF(ish1.EQ.isn1) goto 160
          dfpt1=abs(fpt1(ish1)-fpt1(inx1(ish1)))
          IF(dfpt1.LT.parp(43)*fmax1) goto 150
          npt1=npt1+1
          xpt1(npt1)=0.5d0*(xpt1(ish1)+xpt1(inx1(ish1)))
          inx1(npt1)=inx1(ish1)
          inx1(ish1)=npt1
          goto 140
        ENDIF
 
C...Calculate integral over inner loop.
  160   fsum1=0d0
        DO 170 ipt1=2,npt1
          fsum1=fsum1+0.5d0*(fpt1(ipt1)+fpt1(inx1(ipt1)))*
     &    (xpt1(inx1(ipt1))-xpt1(ipt1))
  170   CONTINUE
        func2=fsum1*(atu1-atl1)/paru(1)
  180   IF(mbw(1).EQ.1.AND.mbw(2).EQ.1) THEN
          IF(func2.GT.fmax2) fmax2=func2
          fpt2(npt2)=func2
 
C...Go to next position in outer loop.
          IF(npt2.EQ.1) THEN
            npt2=npt2+1
            xpt2(npt2)=0d0
            inx2(npt2)=1
            goto 130
          ELSEIF(npt2.LE.8) THEN
            npt2=npt2+1
            IF(npt2.LE.4.OR.npt2.EQ.6) ish2=1
            ish2=ish2+1
            xpt2(npt2)=0.5d0*(xpt2(ish2)+xpt2(inx2(ish2)))
            inx2(npt2)=inx2(ish2)
            inx2(ish2)=npt2
            goto 130
          ELSEIF(npt2.LT.100) THEN
            isn2=ish2
  190       ish2=ish2+1
            IF(ish2.GT.npt2) ish2=2
            IF(ish2.EQ.isn2) goto 200
            dfpt2=abs(fpt2(ish2)-fpt2(inx2(ish2)))
            IF(dfpt2.LT.parp(43)*fmax2) goto 190
            npt2=npt2+1
            xpt2(npt2)=0.5d0*(xpt2(ish2)+xpt2(inx2(ish2)))
            inx2(npt2)=inx2(ish2)
            inx2(ish2)=npt2
            goto 130
          ENDIF
 
C...Calculate integral over outer loop.
  200     fsum2=0d0
          DO 210 ipt2=2,npt2
            fsum2=fsum2+0.5d0*(fpt2(ipt2)+fpt2(inx2(ipt2)))*
     &      (xpt2(inx2(ipt2))-xpt2(ipt2))
  210     CONTINUE
          fsum2=fsum2*(atu2-atl2)/paru(1)
          IF(meql.EQ.1) fsum2=2d0*fsum2
        ELSE
          fsum2=func2
        ENDIF
 
C...Save result; second integration for user-selected mass range.
        IF(loop.EQ.1) widw=fsum2
        wid2=fsum2
        IF(loop.EQ.1.AND.(ckin(46).GE.ckin(45).OR.ckin(48).GE.ckin(47)
     &  .OR.max(ckin(45),ckin(47)).GE.1.01d0*parp(42))) THEN
          loop=2
          goto 100
        ENDIF
        ret1=widw
        ret2=wid2/widw
 
C...Select two decay product masses of a resonance.
      ELSEIF(mofsh.EQ.2.OR.mofsh.EQ.5) THEN
  220   DO 230 i=1,2
          IF(mbw(i).EQ.0) goto 230
          pmbw=pmd(i)**2+pmd(i)*pgd(i)*tan(atl(i)+pyr(0)*
     &    (atu(i)-atl(i)))
          pmg(i)=min(pmu(i),max(pml(i),sqrt(max(0d0,pmbw))))
          rmg(i)=(pmg(i)/pmmx)**2
  230   CONTINUE
        IF((meql.EQ.1.AND.pmg(max(1,mlm)).GT.pmg(min(2,3-mlm))).OR.
     &  pmg(1)+pmg(2)+parj(64).GT.pmmx) goto 220
 
C...Weight with matrix element (if none known, use beta factor).
        flam=sqrt(max(0d0,(1d0-rmg(1)-rmg(2))**2-4d0*rmg(1)*rmg(2)))
        IF(mmed.EQ.1) THEN
          wtbe=flam*((1d0-rmg(1)-rmg(2))**2+8d0*rmg(1)*rmg(2))
        ELSEIF(mmed.EQ.2) THEN
          wtbe=flam**3*(1d0+10d0*rmg(1)+10d0*rmg(2)+rmg(1)**2+
     &    rmg(2)**2+10d0*rmg(1)*rmg(2))
        ELSEIF(mmed.EQ.3) THEN
          wtbe=flam*(rmg(1)+flam**2/12d0)
        ELSE
          wtbe=flam
        ENDIF
        IF(wtbe.LT.pyr(0)) goto 220
        ret1=pmg(1)
        ret2=pmg(2)
 
C...Find suitable set of masses for initialization of 2 -> 2 processes.
      ELSEIF(mofsh.EQ.3) THEN
        IF(mbw(1).NE.0.AND.mbw(2).EQ.0) THEN
          pmg(1)=min(pmd(1),0.5d0*(pml(1)+pmu(1)))
          pmg(2)=pmd(2)
        ELSEIF(mbw(2).NE.0.AND.mbw(1).EQ.0) THEN
          pmg(1)=pmd(1)
          pmg(2)=min(pmd(2),0.5d0*(pml(2)+pmu(2)))
        ELSE
          idiv=-1
  240     idiv=idiv+1
          pmg(1)=min(pmd(1),0.1d0*(idiv*pml(1)+(10-idiv)*pmu(1)))
          pmg(2)=min(pmd(2),0.1d0*(idiv*pml(2)+(10-idiv)*pmu(2)))
          IF(idiv.LE.9.AND.pmg(1)+pmg(2).GT.0.9d0*pmmx) goto 240
        ENDIF
        ret1=pmg(1)
        ret2=pmg(2)
 
C...Evaluate importance of excluded tails of Breit-Wigners.
        IF(meql.EQ.0.AND.mbw(1).EQ.1.AND.mbw(2).EQ.1.AND.pmd(1)+pmd(2)
     &  .GT.pmmx.AND.pmh(1).GT.pml(1).AND.pmh(2).GT.pml(2)) meql=2
        IF(meql.LE.1) THEN
          vint(80)=1d0
          DO 250 i=1,2
            IF(mbw(i).NE.0) vint(80)=vint(80)*1.25d0*(atu(i)-atl(i))/
     &      paru(1)
  250     CONTINUE
        ELSE
          vint(80)=(1.25d0/paru(1))**2*max((atu(1)-atl(1))*
     &    (ath(2)-atl(2)),(ath(1)-atl(1))*(atu(2)-atl(2)))
        ENDIF
        IF((isub.EQ.15.OR.isub.EQ.19.OR.isub.EQ.30.OR.isub.EQ.35).AND.
     &  mstp(43).NE.2) vint(80)=2d0*vint(80)
        IF(isub.EQ.22.AND.mstp(43).NE.2) vint(80)=4d0*vint(80)
        IF(meql.GE.1) vint(80)=2d0*vint(80)
 
C...Pick one particle to be the lighter (if improves efficiency).
      ELSEIF(mofsh.EQ.4) THEN
        IF(meql.EQ.0.AND.mbw(1).EQ.1.AND.mbw(2).EQ.1.AND.pmd(1)+pmd(2)
     &  .GT.pmmx.AND.pmh(1).GT.pml(1).AND.pmh(2).GT.pml(2)) meql=2
  260   IF(meql.EQ.2) mlm=int(1.5d0+pyr(0))
 
C...Select two masses according to Breit-Wigner + flat in s + 1/s.
        DO 270 i=1,2
          IF(mbw(i).EQ.0) goto 270
          pmv=pmu(i)
          IF(meql.EQ.2.AND.i.EQ.mlm) pmv=pmh(i)
          atv=atu(i)
          IF(meql.EQ.2.AND.i.EQ.mlm) atv=ath(i)
          rbr=pyr(0)
          IF((isub.EQ.15.OR.isub.EQ.19.OR.isub.EQ.22.OR.isub.EQ.30.OR.
     &    isub.EQ.35).AND.mstp(43).NE.2) rbr=2d0*rbr
          IF(rbr.LT.0.8d0) THEN
            pmsr=pmd(i)**2+pmd(i)*pgd(i)*tan(atl(i)+pyr(0)*(atv-atl(i)))
            pmg(i)=min(pmv,max(pml(i),sqrt(max(0d0,pmsr))))
          ELSEIF(rbr.LT.0.9d0) THEN
            pmg(i)=sqrt(max(0d0,pml(i)**2+pyr(0)*(pmv**2-pml(i)**2)))
          ELSEIF(rbr.LT.1.5d0) THEN
            pmg(i)=pml(i)*(pmv/pml(i))**pyr(0)
          ELSE
            pmg(i)=sqrt(max(0d0,pml(i)**2*pmv**2/(pml(i)**2+pyr(0)*
     &      (pmv**2-pml(i)**2))))
          ENDIF
  270   CONTINUE
        IF((meql.GE.1.AND.pmg(max(1,mlm)).GT.pmg(min(2,3-mlm))).OR.
     &  pmg(1)+pmg(2)+parj(64).GT.pmmx) THEN
          IF(mint(48).EQ.1.AND.mstp(171).EQ.0) THEN
            ngen(0,1)=ngen(0,1)+1
            ngen(mint(1),1)=ngen(mint(1),1)+1
            goto 260
          ELSE
            mint(51)=1
            RETURN
          ENDIF
        ENDIF
        ret1=pmg(1)
        ret2=pmg(2)
 
C...Give weight for selected mass distribution.
        vint(80)=1d0
        DO 280 i=1,2
          IF(mbw(i).EQ.0) goto 280
          pmv=pmu(i)
          IF(meql.EQ.2.AND.i.EQ.mlm) pmv=pmh(i)
          atv=atu(i)
          IF(meql.EQ.2.AND.i.EQ.mlm) atv=ath(i)
          f0=pmd(i)*pgd(i)/((pmg(i)**2-pmd(i)**2)**2+
     &    (pmd(i)*pgd(i))**2)/paru(1)
          f1=1d0
          f2=1d0/pmg(i)**2
          f3=1d0/pmg(i)**4
          fi0=(atv-atl(i))/paru(1)
          fi1=pmv**2-pml(i)**2
          fi2=2d0*log(pmv/pml(i))
          fi3=1d0/pml(i)**2-1d0/pmv**2
          IF((isub.EQ.15.OR.isub.EQ.19.OR.isub.EQ.22.OR.isub.EQ.30.OR.
     &    isub.EQ.35).AND.mstp(43).NE.2) THEN
            vint(80)=vint(80)*20d0/(8d0+(fi0/f0)*(f1/fi1+6d0*f2/fi2+
     &      5d0*f3/fi3))
          ELSE
            vint(80)=vint(80)*10d0/(8d0+(fi0/f0)*(f1/fi1+f2/fi2))
          ENDIF
          vint(80)=vint(80)*fi0
  280   CONTINUE
        IF(meql.GE.1) vint(80)=2d0*vint(80)
      ENDIF
 
      RETURN
      END