pykmap.f

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

 
C*********************************************************************
 
C...PYKMAP
C...Maps a uniform distribution into a distribution of a kinematical
C...variable according to one of the possibilities allowed. It is
C...assumed that kinematical limits have been set by a PYKLIM call.
 
      SUBROUTINE pykmap(IVAR,MVAR,VVAR)
 
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/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)
      SAVE /pydat1/,/pydat2/,/pysubs/,/pypars/,/pyint1/,/pyint2/
 
C...Convert VVAR to tau variable.
      isub=mint(1)
      istsb=iset(isub)
      IF(ivar.EQ.1) THEN
        taumin=vint(11)
        taumax=vint(31)
        IF(mvar.EQ.3.OR.mvar.EQ.4) THEN
          taure=vint(73)
          gamre=vint(74)
        ELSEIF(mvar.EQ.5.OR.mvar.EQ.6) THEN
          taure=vint(75)
          gamre=vint(76)
        ENDIF
        IF(mint(47).EQ.1.AND.(istsb.EQ.1.OR.istsb.EQ.2)) THEN
          tau=1d0
        ELSEIF(mvar.EQ.1) THEN
          tau=taumin*(taumax/taumin)**vvar
        ELSEIF(mvar.EQ.2) THEN
          tau=taumax*taumin/(taumin+(taumax-taumin)*vvar)
        ELSEIF(mvar.EQ.3.OR.mvar.EQ.5) THEN
          ratgen=(taure+taumax)/(taure+taumin)*taumin/taumax
          tau=taure*taumin/((taure+taumin)*ratgen**vvar-taumin)
        ELSEIF(mvar.EQ.4.OR.mvar.EQ.6) THEN
          aupp=atan((taumax-taure)/gamre)
          alow=atan((taumin-taure)/gamre)
          tau=taure+gamre*tan(alow+(aupp-alow)*vvar)
        ELSEIF(mint(47).EQ.5) THEN
          aupp=log(max(2d-10,1d0-taumax))
          alow=log(max(2d-10,1d0-taumin))
          tau=1d0-exp(aupp+vvar*(alow-aupp))
        ELSE
          aupp=log(max(1d-10,1d0-taumax))
          alow=log(max(1d-10,1d0-taumin))
          tau=1d0-exp(aupp+vvar*(alow-aupp))
        ENDIF
        vint(21)=min(taumax,max(taumin,tau))
 
C...Convert VVAR to y* variable.
      ELSEIF(ivar.EQ.2) THEN
        ystmin=vint(12)
        ystmax=vint(32)
        taue=vint(21)
        IF(istsb.GE.3.AND.istsb.LE.5) taue=vint(26)
        IF(mint(47).EQ.1) THEN
          yst=0d0
        ELSEIF(mint(47).EQ.2.OR.mint(47).EQ.6) THEN
          yst=-0.5d0*log(taue)
        ELSEIF(mint(47).EQ.3.OR.mint(47).EQ.7) THEN
          yst=0.5d0*log(taue)
        ELSEIF(mvar.EQ.1) THEN
          yst=ystmin+(ystmax-ystmin)*sqrt(vvar)
        ELSEIF(mvar.EQ.2) THEN
          yst=ystmax-(ystmax-ystmin)*sqrt(1d0-vvar)
        ELSEIF(mvar.EQ.3) THEN
          aupp=atan(exp(ystmax))
          alow=atan(exp(ystmin))
          yst=log(tan(alow+(aupp-alow)*vvar))
        ELSEIF(mvar.EQ.4) THEN
          yst0=-0.5d0*log(taue)
          aupp=log(max(1d-10,exp(yst0-ystmin)-1d0))
          alow=log(max(1d-10,exp(yst0-ystmax)-1d0))
          yst=yst0-log(1d0+exp(alow+vvar*(aupp-alow)))
        ELSE
          yst0=-0.5d0*log(taue)
          aupp=log(max(1d-10,exp(yst0+ystmin)-1d0))
          alow=log(max(1d-10,exp(yst0+ystmax)-1d0))
          yst=log(1d0+exp(aupp+vvar*(alow-aupp)))-yst0
        ENDIF
        vint(22)=min(ystmax,max(ystmin,yst))
 
C...Convert VVAR to cos(theta-hat) variable.
      ELSEIF(ivar.EQ.3) THEN
        rm34=max(1d-20,2d0*vint(63)*vint(64)/(vint(21)*vint(2))**2)
        rsqm=1d0+rm34
        IF(2d0*vint(71)**2/(vint(21)*vint(2)).LT.0.0001d0)
     &  rm34=max(rm34,2d0*vint(71)**2/(vint(21)*vint(2)))
        ctnmin=vint(13)
        ctnmax=vint(33)
        ctpmin=vint(14)
        ctpmax=vint(34)
        IF(mvar.EQ.1) THEN
          aneg=ctnmax-ctnmin
          apos=ctpmax-ctpmin
          IF(aneg.GT.0d0.AND.vvar*(aneg+apos).LE.aneg) THEN
            vctn=vvar*(aneg+apos)/aneg
            cth=ctnmin+(ctnmax-ctnmin)*vctn
          ELSE
            vctp=(vvar*(aneg+apos)-aneg)/apos
            cth=ctpmin+(ctpmax-ctpmin)*vctp
          ENDIF
        ELSEIF(mvar.EQ.2) THEN
          rmnmin=max(rm34,rsqm-ctnmin)
          rmnmax=max(rm34,rsqm-ctnmax)
          rmpmin=max(rm34,rsqm-ctpmin)
          rmpmax=max(rm34,rsqm-ctpmax)
          aneg=log(rmnmin/rmnmax)
          apos=log(rmpmin/rmpmax)
          IF(aneg.GT.0d0.AND.vvar*(aneg+apos).LE.aneg) THEN
            vctn=vvar*(aneg+apos)/aneg
            cth=rsqm-rmnmin*(rmnmax/rmnmin)**vctn
          ELSE
            vctp=(vvar*(aneg+apos)-aneg)/apos
            cth=rsqm-rmpmin*(rmpmax/rmpmin)**vctp
          ENDIF
        ELSEIF(mvar.EQ.3) THEN
          rmnmin=max(rm34,rsqm+ctnmin)
          rmnmax=max(rm34,rsqm+ctnmax)
          rmpmin=max(rm34,rsqm+ctpmin)
          rmpmax=max(rm34,rsqm+ctpmax)
          aneg=log(rmnmax/rmnmin)
          apos=log(rmpmax/rmpmin)
          IF(aneg.GT.0d0.AND.vvar*(aneg+apos).LE.aneg) THEN
            vctn=vvar*(aneg+apos)/aneg
            cth=rmnmin*(rmnmax/rmnmin)**vctn-rsqm
          ELSE
            vctp=(vvar*(aneg+apos)-aneg)/apos
            cth=rmpmin*(rmpmax/rmpmin)**vctp-rsqm
          ENDIF
        ELSEIF(mvar.EQ.4) THEN
          rmnmin=max(rm34,rsqm-ctnmin)
          rmnmax=max(rm34,rsqm-ctnmax)
          rmpmin=max(rm34,rsqm-ctpmin)
          rmpmax=max(rm34,rsqm-ctpmax)
          aneg=1d0/rmnmax-1d0/rmnmin
          apos=1d0/rmpmax-1d0/rmpmin
          IF(aneg.GT.0d0.AND.vvar*(aneg+apos).LE.aneg) THEN
            vctn=vvar*(aneg+apos)/aneg
            cth=rsqm-1d0/(1d0/rmnmin+aneg*vctn)
          ELSE
            vctp=(vvar*(aneg+apos)-aneg)/apos
            cth=rsqm-1d0/(1d0/rmpmin+apos*vctp)
          ENDIF
        ELSEIF(mvar.EQ.5) THEN
          rmnmin=max(rm34,rsqm+ctnmin)
          rmnmax=max(rm34,rsqm+ctnmax)
          rmpmin=max(rm34,rsqm+ctpmin)
          rmpmax=max(rm34,rsqm+ctpmax)
          aneg=1d0/rmnmin-1d0/rmnmax
          apos=1d0/rmpmin-1d0/rmpmax
          IF(aneg.GT.0d0.AND.vvar*(aneg+apos).LE.aneg) THEN
            vctn=vvar*(aneg+apos)/aneg
            cth=1d0/(1d0/rmnmin-aneg*vctn)-rsqm
          ELSE
            vctp=(vvar*(aneg+apos)-aneg)/apos
            cth=1d0/(1d0/rmpmin-apos*vctp)-rsqm
          ENDIF
        ENDIF
        IF(cth.LT.0d0) cth=min(ctnmax,max(ctnmin,cth))
        IF(cth.GT.0d0) cth=min(ctpmax,max(ctpmin,cth))
        vint(23)=cth
 
C...Convert VVAR to tau' variable.
      ELSEIF(ivar.EQ.4) THEN
        tau=vint(21)
        taupmn=vint(16)
        taupmx=vint(36)
        IF(mint(47).EQ.1) THEN
          taup=1d0
        ELSEIF(mvar.EQ.1) THEN
          taup=taupmn*(taupmx/taupmn)**vvar
        ELSEIF(mvar.EQ.2) THEN
          aupp=(1d0-tau/taupmx)**4
          alow=(1d0-tau/taupmn)**4
          taup=tau/max(1d-10,1d0-(alow+(aupp-alow)*vvar)**0.25d0)
        ELSEIF(mint(47).EQ.5) THEN
          aupp=log(max(2d-10,1d0-taupmx))
          alow=log(max(2d-10,1d0-taupmn))
          taup=1d0-exp(aupp+vvar*(alow-aupp))
        ELSE
          aupp=log(max(1d-10,1d0-taupmx))
          alow=log(max(1d-10,1d0-taupmn))
          taup=1d0-exp(aupp+vvar*(alow-aupp))
        ENDIF
        vint(26)=min(taupmx,max(taupmn,taup))
 
C...Selection of extra variables needed in 2 -> 3 process:
C...pT1, pT2, phi1, phi2, y3 for three outgoing particles.
C...Since no options are available, the functions of PYKLIM
C...and PYKMAP are joint for these choices.
      ELSEIF(ivar.EQ.5) THEN
 
C...Read out total energy and particle masses.
        mint(51)=0
        mptpk=1
        IF(isub.EQ.123.OR.isub.EQ.124.OR.isub.EQ.173.OR.isub.EQ.174
     &  .OR.isub.EQ.178.OR.isub.EQ.179.OR.isub.EQ.351.OR.isub.EQ.352)
     &  mptpk=2
        shp=vint(26)*vint(2)
        shpr=sqrt(shp)
        pm1=vint(201)
        pm2=vint(206)
        pm3=sqrt(vint(21))*vint(1)
        IF(pm1+pm2+pm3.GT.0.9999d0*shpr) THEN
          mint(51)=1
          RETURN
        ENDIF
        pmrs1=vint(204)**2
        pmrs2=vint(209)**2
 
C...Specify coefficients of pT choice; upper and lower limits.
        IF(mptpk.EQ.1) THEN
          hwt1=0.4d0
          hwt2=0.4d0
        ELSE
          hwt1=0.05d0
          hwt2=0.05d0
        ENDIF
        hwt3=1d0-hwt1-hwt2
        ptsmx1=((shp-pm1**2-(pm2+pm3)**2)**2-(2d0*pm1*(pm2+pm3))**2)/
     &  (4d0*shp)
        IF(ckin(52).GT.0d0) ptsmx1=min(ptsmx1,ckin(52)**2)
        ptsmn1=ckin(51)**2
        ptsmx2=((shp-pm2**2-(pm1+pm3)**2)**2-(2d0*pm2*(pm1+pm3))**2)/
     &  (4d0*shp)
        IF(ckin(54).GT.0d0) ptsmx2=min(ptsmx2,ckin(54)**2)
        ptsmn2=ckin(53)**2
 
C...Select transverse momenta according to
C...dp_T^2 * (a + b/(M^2 + p_T^2) + c/(M^2 + p_T^2)^2).
        hmx=pmrs1+ptsmx1
        hmn=pmrs1+ptsmn1
        IF(hmx.LT.1.0001d0*hmn) THEN
          mint(51)=1
          RETURN
        ENDIF
        hde=ptsmx1-ptsmn1
        rpt=pyr(0)
        IF(rpt.LT.hwt1) THEN
          pts1=ptsmn1+pyr(0)*hde
        ELSEIF(rpt.LT.hwt1+hwt2) THEN
          pts1=max(ptsmn1,hmn*(hmx/hmn)**pyr(0)-pmrs1)
        ELSE
          pts1=max(ptsmn1,hmn*hmx/(hmn+pyr(0)*hde)-pmrs1)
        ENDIF
        wtpts1=hde/(hwt1+hwt2*hde/(log(hmx/hmn)*(pmrs1+pts1))+
     &  hwt3*hmn*hmx/(pmrs1+pts1)**2)
        hmx=pmrs2+ptsmx2
        hmn=pmrs2+ptsmn2
        IF(hmx.LT.1.0001d0*hmn) THEN
          mint(51)=1
          RETURN
        ENDIF
        hde=ptsmx2-ptsmn2
        rpt=pyr(0)
        IF(rpt.LT.hwt1) THEN
          pts2=ptsmn2+pyr(0)*hde
        ELSEIF(rpt.LT.hwt1+hwt2) THEN
          pts2=max(ptsmn2,hmn*(hmx/hmn)**pyr(0)-pmrs2)
        ELSE
          pts2=max(ptsmn2,hmn*hmx/(hmn+pyr(0)*hde)-pmrs2)
        ENDIF
        wtpts2=hde/(hwt1+hwt2*hde/(log(hmx/hmn)*(pmrs2+pts2))+
     &  hwt3*hmn*hmx/(pmrs2+pts2)**2)
 
C...Select azimuthal angles and check pT choice.
        phi1=paru(2)*pyr(0)
        phi2=paru(2)*pyr(0)
        phir=phi2-phi1
        pts3=max(0d0,pts1+pts2+2d0*sqrt(pts1*pts2)*cos(phir))
        IF(pts3.LT.ckin(55)**2.OR.(ckin(56).GT.0d0.AND.pts3.GT.
     &  ckin(56)**2)) THEN
          mint(51)=1
          RETURN
        ENDIF
 
C...Calculate transverse masses and check phase space not closed.
        pms1=pm1**2+pts1
        pms2=pm2**2+pts2
        pms3=pm3**2+pts3
        pmt1=sqrt(pms1)
        pmt2=sqrt(pms2)
        pmt3=sqrt(pms3)
        pm12=(pmt1+pmt2)**2
        IF(pmt1+pmt2+pmt3.GT.0.9999d0*shpr) THEN
          mint(51)=1
          RETURN
        ENDIF
 
C...Select rapidity for particle 3 and check phase space not closed.
        y3max=log((shp+pms3-pm12+sqrt(max(0d0,(shp-pms3-pm12)**2-
     &  4d0*pms3*pm12)))/(2d0*shpr*pmt3))
        IF(y3max.LT.1d-6) THEN
          mint(51)=1
          RETURN
        ENDIF
        y3=(2d0*pyr(0)-1d0)*0.999999d0*y3max
        pz3=pmt3*sinh(y3)
        pe3=pmt3*cosh(y3)
 
C...Find momentum transfers in two mirror solutions (in 1-2 frame).
        pz12=-pz3
        pe12=shpr-pe3
        pms12=pe12**2-pz12**2
        sql12=sqrt(max(0d0,(pms12-pms1-pms2)**2-4d0*pms1*pms2))
        IF(sql12.LT.1d-6*shp) THEN
          mint(51)=1
          RETURN
        ENDIF
        pmm1=pms12+pms1-pms2
        pmm2=pms12+pms2-pms1
        tfac=-shpr/(2d0*pms12)
        t1p=tfac*(pe12-pz12)*(pmm1-sql12)
        t1n=tfac*(pe12-pz12)*(pmm1+sql12)
        t2p=tfac*(pe12+pz12)*(pmm2-sql12)
        t2n=tfac*(pe12+pz12)*(pmm2+sql12)
 
C...Construct relative mirror weights and make choice.
        IF(mptpk.EQ.1.OR.isub.EQ.351.OR.isub.EQ.352) THEN
          wtpu=1d0
          wtnu=1d0
        ELSE
          wtpu=1d0/((t1p-pmrs1)*(t2p-pmrs2))**2
          wtnu=1d0/((t1n-pmrs1)*(t2n-pmrs2))**2
        ENDIF
        wtp=wtpu/(wtpu+wtnu)
        wtn=wtnu/(wtpu+wtnu)
        eps=1d0
        IF(wtn.GT.pyr(0)) eps=-1d0
 
C...Store result of variable choice and associated weights.
        vint(202)=pts1
        vint(207)=pts2
        vint(203)=phi1
        vint(208)=phi2
        vint(205)=wtpts1
        vint(210)=wtpts2
        vint(211)=y3
        vint(212)=y3max
        vint(213)=eps
        IF(eps.GT.0d0) THEN
          vint(214)=1d0/wtp
          vint(215)=t1p
          vint(216)=t2p
        ELSE
          vint(214)=1d0/wtn
          vint(215)=t1n
          vint(216)=t2n
        ENDIF
        vint(217)=-0.5d0*tfac*(pe12-pz12)*(pmm2+eps*sql12)
        vint(218)=-0.5d0*tfac*(pe12+pz12)*(pmm1+eps*sql12)
        vint(219)=0.5d0*(pms12-pts3)
        vint(220)=sql12
      ENDIF
 
      RETURN
      END