pyhimult.f

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

    
C*********************************************************************  
    
      SUBROUTINE pyhimult(MMUL)   
    
C...Initializes treatment of multiple interactions, selects kinematics  
C...of hardest interaction if low-pT physics included in run, and   
C...generates all non-hardest interactions. 
      common/lujets/n,k(9000,5),p(9000,5),v(9000,5)
      SAVE /lujets/ 
      common/ludat1/mstu(200),paru(200),mstj(200),parj(200) 
      SAVE /ludat1/ 
      common/ludat2/kchg(500,3),pmas(500,4),parf(2000),vckm(4,4)    
      SAVE /ludat2/ 
      common/pyhisubs/msel,msub(200),kfin(2,-40:40),ckin(200) 
      SAVE /pyhisubs/ 
      common/pyhipars/mstp(200),parp(200),msti(200),pari(200) 
      SAVE /pyhipars/ 
      common/pyhiint1/mint(400),vint(400) 
      SAVE /pyhiint1/ 
      common/pyhiint2/iset(200),kfpr(200,2),coef(200,20),icol(40,4,2) 
      SAVE /pyhiint2/ 
      common/pyhiint3/xsfx(2,-40:40),isig(1000,3),sigh(1000)  
      SAVE /pyhiint3/ 
      common/pyhiint5/ngen(0:200,3),xsec(0:200,3) 
      SAVE /pyhiint5/ 
      dimension nmul(20),sigm(20),kstr(500,2)   
      SAVE xt2,xt2fac,xc2,xts,irbin,rbin,nmul,sigm  
    
C...Initialization of multiple interaction treatment.   
      IF(mmul.EQ.1) THEN    
        IF(mstp(122).GE.1) WRITE(mstu(11),1000) mstp(82)    
        isub=96 
        mint(1)=96  
        vint(63)=0. 
        vint(64)=0. 
        vint(143)=1.    
        vint(144)=1.    
    
C...Loop over phase space points: xT2 choice in 20 bins.    
  100   sigsum=0.   
        DO 120 ixt2=1,20    
        nmul(ixt2)=mstp(83) 
        sigm(ixt2)=0.   
        DO 110 itry=1,mstp(83)  
        rsca=0.05*((21-ixt2)-rlu(0))    
        xt2=vint(149)*(1.+vint(149))/(vint(149)+rsca)-vint(149) 
        xt2=max(0.01*vint(149),xt2) 
        vint(25)=xt2    
    
C...Choose tau and y*. Calculate cos(theta-hat).    
        IF(rlu(0).LE.coef(isub,1)) THEN 
          taup=(2.*(1.+sqrt(1.-xt2))/xt2-1.)**rlu(0)    
          tau=xt2*(1.+taup)**2/(4.*taup)    
        ELSE    
          tau=xt2*(1.+tan(rlu(0)*atan(sqrt(1./xt2-1.)))**2) 
        ENDIF   
        vint(21)=tau    
        CALL pyhiklim(2)  
        ryst=rlu(0) 
        myst=1  
        IF(ryst.GT.coef(isub,7)) myst=2 
        IF(ryst.GT.coef(isub,7)+coef(isub,8)) myst=3    
        CALL pyhikmap(2,myst,rlu(0))  
        vint(23)=sqrt(max(0.,1.-xt2/tau))*(-1)**int(1.5+rlu(0)) 
    
C...Calculate differential cross-section.   
        vint(71)=0.5*vint(1)*sqrt(xt2)  
        CALL pyhisigh(nchn,sigs)  
  110   sigm(ixt2)=sigm(ixt2)+sigs  
  120   sigsum=sigsum+sigm(ixt2)    
        sigsum=sigsum/(20.*mstp(83))    
    
C...Reject result if sigma(parton-parton) is smaller than hadronic one. 
        IF(sigsum.LT.1.1*vint(106)) THEN    
          IF(mstp(122).GE.1) WRITE(mstu(11),1100) parp(82),sigsum   
          parp(82)=0.9*parp(82) 
          vint(149)=4.*parp(82)**2/vint(2)  
          goto 100  
        ENDIF   
        IF(mstp(122).GE.1) WRITE(mstu(11),1200) parp(82), sigsum    
    
C...Start iteration to find k factor.   
        yke=sigsum/vint(106)    
        so=0.5  
        xi=0.   
        yi=0.   
        xk=0.5  
        iit=0   
  130   IF(iit.EQ.0) THEN   
          xk=2.*xk  
        ELSEIF(iit.EQ.1) THEN   
          xk=0.5*xk 
        ELSE    
          xk=xi+(yke-yi)*(xf-xi)/(yf-yi)    
        ENDIF   
    
C...Evaluate overlap integrals. 
        IF(mstp(82).EQ.2) THEN  
          sp=0.5*paru(1)*(1.-exp(-xk))  
          sop=sp/paru(1)    
        ELSE    
          IF(mstp(82).EQ.3) deltab=0.02 
          IF(mstp(82).EQ.4) deltab=min(0.01,0.05*parp(84))  
          sp=0. 
          sop=0.    
          b=-0.5*deltab 
  140     b=b+deltab    
          IF(mstp(82).EQ.3) THEN    
            ov=exp(-b**2)/paru(2)   
          ELSE  
            cq2=parp(84)**2 
            ov=((1.-parp(83))**2*exp(-min(100.,b**2))+2.*parp(83)*  
     &      (1.-parp(83))*2./(1.+cq2)*exp(-min(100.,b**2*2./(1.+cq2)))+ 
     &      parp(83)**2/cq2*exp(-min(100.,b**2/cq2)))/paru(2)   
          ENDIF 
          pacc=1.-exp(-min(100.,paru(1)*xk*ov)) 
          sp=sp+paru(2)*b*deltab*pacc   
          sop=sop+paru(2)*b*deltab*ov*pacc  
          IF(b.LT.1..OR.b*pacc.GT.1e-6) goto 140    
        ENDIF   
        yk=paru(1)*xk*so/sp 
    
C...Continue iteration until convergence.   
        IF(yk.LT.yke) THEN  
          xi=xk 
          yi=yk 
          IF(iit.EQ.1) iit=2    
        ELSE    
          xf=xk 
          yf=yk 
          IF(iit.EQ.0) iit=1    
        ENDIF   
        IF(abs(yk-yke).GE.1e-5*yke) goto 130    
    
C...Store some results for subsequent use.  
        vint(145)=sigsum    
        vint(146)=sop/so    
        vint(147)=sop/sp    
    
C...Initialize iteration in xT2 for hardest interaction.    
      ELSEIF(mmul.EQ.2) THEN    
        IF(mstp(82).LE.0) THEN  
        ELSEIF(mstp(82).EQ.1) THEN  
          xt2=1.    
          xt2fac=xsec(96,1)/vint(106)*vint(149)/(1.-vint(149))  
        ELSEIF(mstp(82).EQ.2) THEN  
          xt2=1.    
          xt2fac=vint(146)*xsec(96,1)/vint(106)*vint(149)*(1.+vint(149))    
        ELSE    
          xc2=4.*ckin(3)**2/vint(2) 
          IF(ckin(3).LE.ckin(5).OR.mint(82).GE.2) xc2=0.    
        ENDIF   
    
      ELSEIF(mmul.EQ.3) THEN    
C...Low-pT or multiple interactions (first semihard interaction):   
C...choose xT2 according to dpT2/pT2**2*exp(-(sigma above pT2)/norm)    
C...or (MSTP(82)>=2) dpT2/(pT2+pT0**2)**2*exp(-....).   
        isub=mint(1)    
        IF(mstp(82).LE.0) THEN  
          xt2=0.    
        ELSEIF(mstp(82).EQ.1) THEN  
          xt2=xt2fac*xt2/(xt2fac-xt2*log(rlu(0)))   
        ELSEIF(mstp(82).EQ.2) THEN  
          IF(xt2.LT.1..AND.exp(-xt2fac*xt2/(vint(149)*(xt2+ 
     &    vint(149)))).GT.rlu(0)) xt2=1.    
          IF(xt2.GE.1.) THEN    
            xt2=(1.+vint(149))*xt2fac/(xt2fac-(1.+vint(149))*log(1.-    
     &      rlu(0)*(1.-exp(-xt2fac/(vint(149)*(1.+vint(149)))))))-  
     &      vint(149)   
          ELSE  
            xt2=-xt2fac/log(exp(-xt2fac/(xt2+vint(149)))+rlu(0)*    
     &      (exp(-xt2fac/vint(149))-exp(-xt2fac/(xt2+vint(149)))))- 
     &      vint(149)   
          ENDIF 
          xt2=max(0.01*vint(149),xt2)   
        ELSE    
          xt2=(xc2+vint(149))*(1.+vint(149))/(1.+vint(149)- 
     &    rlu(0)*(1.-xc2))-vint(149)    
          xt2=max(0.01*vint(149),xt2)   
        ENDIF   
        vint(25)=xt2    
    
C...Low-pT: choose xT2, tau, y* and cos(theta-hat) fixed.   
        IF(mstp(82).LE.1.AND.xt2.LT.vint(149)) THEN 
          IF(mint(82).EQ.1) ngen(0,1)=ngen(0,1)-1   
          IF(mint(82).EQ.1) ngen(isub,1)=ngen(isub,1)-1 
          isub=95   
          mint(1)=isub  
          vint(21)=0.01*vint(149)   
          vint(22)=0.   
          vint(23)=0.   
          vint(25)=0.01*vint(149)   
    
        ELSE    
C...Multiple interactions (first semihard interaction). 
C...Choose tau and y*. Calculate cos(theta-hat).    
          IF(rlu(0).LE.coef(isub,1)) THEN   
            taup=(2.*(1.+sqrt(1.-xt2))/xt2-1.)**rlu(0)  
            tau=xt2*(1.+taup)**2/(4.*taup)  
          ELSE  
            tau=xt2*(1.+tan(rlu(0)*atan(sqrt(1./xt2-1.)))**2)   
          ENDIF 
          vint(21)=tau  
          CALL pyhiklim(2)    
          ryst=rlu(0)   
          myst=1    
          IF(ryst.GT.coef(isub,7)) myst=2   
          IF(ryst.GT.coef(isub,7)+coef(isub,8)) myst=3  
          CALL pyhikmap(2,myst,rlu(0))    
          vint(23)=sqrt(max(0.,1.-xt2/tau))*(-1)**int(1.5+rlu(0))   
        ENDIF   
        vint(71)=0.5*vint(1)*sqrt(vint(25)) 
    
C...Store results of cross-section calculation. 
      ELSEIF(mmul.EQ.4) THEN    
        isub=mint(1)    
        xts=vint(25)    
        IF(iset(isub).EQ.1) xts=vint(21)    
        IF(iset(isub).EQ.2) xts=(4.*vint(48)+2.*vint(63)+2.*vint(64))/  
     &  vint(2) 
        IF(iset(isub).EQ.3.OR.iset(isub).EQ.4) xts=vint(26) 
        rbin=max(0.000001,min(0.999999,xts*(1.+vint(149))/  
     &  (xts+vint(149))))   
        irbin=int(1.+20.*rbin)  
        IF(isub.EQ.96) nmul(irbin)=nmul(irbin)+1    
        IF(isub.EQ.96) sigm(irbin)=sigm(irbin)+vint(153)    
    
C...Choose impact parameter.    
      ELSEIF(mmul.EQ.5) THEN    
        IF(mstp(82).EQ.3) THEN  
          vint(148)=rlu(0)/(paru(2)*vint(147))  
        ELSE    
          rtype=rlu(0)  
          cq2=parp(84)**2   
          IF(rtype.LT.(1.-parp(83))**2) THEN    
            b2=-log(rlu(0)) 
          ELSEIF(rtype.LT.1.-parp(83)**2) THEN  
            b2=-0.5*(1.+cq2)*log(rlu(0))    
          ELSE  
            b2=-cq2*log(rlu(0)) 
          ENDIF 
          vint(148)=((1.-parp(83))**2*exp(-min(100.,b2))+2.*parp(83)*   
     &    (1.-parp(83))*2./(1.+cq2)*exp(-min(100.,b2*2./(1.+cq2)))+ 
     &    parp(83)**2/cq2*exp(-min(100.,b2/cq2)))/(paru(2)*vint(147))   
        ENDIF   
    
C...Multiple interactions (variable impact parameter) : reject with 
C...probability exp(-overlap*cross-section above pT/normalization). 
        rncor=(irbin-20.*rbin)*nmul(irbin)  
        sigcor=(irbin-20.*rbin)*sigm(irbin) 
        DO 150 ibin=irbin+1,20  
        rncor=rncor+nmul(ibin)  
  150   sigcor=sigcor+sigm(ibin)    
        sigabv=(sigcor/rncor)*vint(149)*(1.-xts)/(xts+vint(149))    
        vint(150)=exp(-min(100.,vint(146)*vint(148)*sigabv/vint(106)))  
    
C...Generate additional multiple semihard interactions. 
      ELSEIF(mmul.EQ.6) THEN    
    
C...Reconstruct strings in hard scattering. 
        isub=mint(1)    
        nmax=mint(84)+4 
        IF(iset(isub).EQ.1) nmax=mint(84)+2 
        nstr=0  
        DO 170 i=mint(84)+1,nmax    
        kcs=kchg(lucomp(k(i,2)),2)*isign(1,k(i,2))  
        IF(kcs.EQ.0) goto 170   
        DO 160 j=1,4    
        IF(kcs.EQ.1.AND.(j.EQ.2.OR.j.EQ.4)) goto 160    
        IF(kcs.EQ.-1.AND.(j.EQ.1.OR.j.EQ.3)) goto 160   
        IF(j.LE.2) THEN 
          ist=mod(k(i,j+3)/mstu(5),mstu(5)) 
        ELSE    
          ist=mod(k(i,j+1),mstu(5)) 
        ENDIF   
        IF(ist.LT.mint(84).OR.ist.GT.i) goto 160    
        IF(kchg(lucomp(k(ist,2)),2).EQ.0) goto 160  
        nstr=nstr+1 
        IF(j.EQ.1.OR.j.EQ.4) THEN   
          kstr(nstr,1)=i    
          kstr(nstr,2)=ist  
        ELSE    
          kstr(nstr,1)=ist  
          kstr(nstr,2)=i    
        ENDIF   
  160   CONTINUE    
  170   CONTINUE    
    
C...Set up starting values for iteration in xT2.    
        xt2=vint(25)    
        IF(iset(isub).EQ.1) xt2=vint(21)    
        IF(iset(isub).EQ.2) xt2=(4.*vint(48)+2.*vint(63)+2.*vint(64))/  
     &  vint(2) 
        IF(iset(isub).EQ.3.OR.iset(isub).EQ.4) xt2=vint(26) 
        isub=96 
        mint(1)=96  
        IF(mstp(82).LE.1) THEN  
          xt2fac=xsec(isub,1)*vint(149)/((1.-vint(149))*vint(106))  
        ELSE    
          xt2fac=vint(146)*vint(148)*xsec(isub,1)/vint(106)*    
     &    vint(149)*(1.+vint(149))  
        ENDIF   
        vint(63)=0. 
        vint(64)=0. 
        vint(151)=0.    
        vint(152)=0.    
        vint(143)=1.-vint(141)  
        vint(144)=1.-vint(142)  
    
C...Iterate downwards in xT2.   
  180   IF(mstp(82).LE.1) THEN  
          xt2=xt2fac*xt2/(xt2fac-xt2*log(rlu(0)))   
          IF(xt2.LT.vint(149)) goto 220 
        ELSE    
          IF(xt2.LE.0.01*vint(149)) goto 220    
          xt2=xt2fac*(xt2+vint(149))/(xt2fac-(xt2+vint(149))*   
     &    log(rlu(0)))-vint(149)    
          IF(xt2.LE.0.) goto 220    
          xt2=max(0.01*vint(149),xt2)   
        ENDIF   
        vint(25)=xt2    
    
C...Choose tau and y*. Calculate cos(theta-hat).    
        IF(rlu(0).LE.coef(isub,1)) THEN 
          taup=(2.*(1.+sqrt(1.-xt2))/xt2-1.)**rlu(0)    
          tau=xt2*(1.+taup)**2/(4.*taup)    
        ELSE    
          tau=xt2*(1.+tan(rlu(0)*atan(sqrt(1./xt2-1.)))**2) 
        ENDIF   
        vint(21)=tau    
        CALL pyhiklim(2)  
        ryst=rlu(0) 
        myst=1  
        IF(ryst.GT.coef(isub,7)) myst=2 
        IF(ryst.GT.coef(isub,7)+coef(isub,8)) myst=3    
        CALL pyhikmap(2,myst,rlu(0))  
        vint(23)=sqrt(max(0.,1.-xt2/tau))*(-1)**int(1.5+rlu(0)) 
    
C...Check that x not used up. Accept or reject kinematical variables.   
        x1m=sqrt(tau)*exp(vint(22)) 
        x2m=sqrt(tau)*exp(-vint(22))    
        IF(vint(143)-x1m.LT.0.01.OR.vint(144)-x2m.LT.0.01) goto 180 
        vint(71)=0.5*vint(1)*sqrt(xt2)  
        CALL pyhisigh(nchn,sigs)  
        IF(sigs.LT.xsec(isub,1)*rlu(0)) goto 180    
    
C...Reset K, P and V vectors. Select some variables.    
        DO 190 i=n+1,n+2    
        DO 190 j=1,5    
        k(i,j)=0    
        p(i,j)=0.   
  190   v(i,j)=0.   
        rflav=rlu(0)    
        pt=0.5*vint(1)*sqrt(xt2)    
        phi=paru(2)*rlu(0)  
        cth=vint(23)    
    
C...Add first parton to event record.   
        k(n+1,1)=3  
        k(n+1,2)=21 
        IF(rflav.GE.max(parp(85),parp(86))) k(n+1,2)=   
     &  1+int((2.+parj(2))*rlu(0))  
        p(n+1,1)=pt*cos(phi)    
        p(n+1,2)=pt*sin(phi)    
        p(n+1,3)=0.25*vint(1)*(vint(41)*(1.+cth)-vint(42)*(1.-cth)) 
        p(n+1,4)=0.25*vint(1)*(vint(41)*(1.+cth)+vint(42)*(1.-cth)) 
        p(n+1,5)=0. 
    
C...Add second parton to event record.  
        k(n+2,1)=3  
        k(n+2,2)=21 
        IF(k(n+1,2).NE.21) k(n+2,2)=-k(n+1,2)   
        p(n+2,1)=-p(n+1,1)  
        p(n+2,2)=-p(n+1,2)  
        p(n+2,3)=0.25*vint(1)*(vint(41)*(1.-cth)-vint(42)*(1.+cth)) 
        p(n+2,4)=0.25*vint(1)*(vint(41)*(1.-cth)+vint(42)*(1.+cth)) 
        p(n+2,5)=0. 
    
        IF(rflav.LT.parp(85).AND.nstr.GE.1) THEN    
C....Choose relevant string pieces to place gluons on.  
          DO 210 i=n+1,n+2  
          dmin=1e8  
          DO 200 istr=1,nstr    
          i1=kstr(istr,1)   
          i2=kstr(istr,2)   
          dist=(p(i,4)*p(i1,4)-p(i,1)*p(i1,1)-p(i,2)*p(i1,2)-   
     &    p(i,3)*p(i1,3))*(p(i,4)*p(i2,4)-p(i,1)*p(i2,1)-   
     &    p(i,2)*p(i2,2)-p(i,3)*p(i2,3))/max(1.,p(i1,4)*p(i2,4)-    
     &    p(i1,1)*p(i2,1)-p(i1,2)*p(i2,2)-p(i1,3)*p(i2,3))  
          IF(istr.EQ.1.OR.dist.LT.dmin) THEN    
            dmin=dist   
            ist1=i1 
            ist2=i2 
            istm=istr   
          ENDIF 
  200     CONTINUE  
    
C....Colour flow adjustments, new string pieces.    
          IF(k(ist1,4)/mstu(5).EQ.ist2) k(ist1,4)=mstu(5)*i+    
     &    mod(k(ist1,4),mstu(5))    
          IF(mod(k(ist1,5),mstu(5)).EQ.ist2) k(ist1,5)= 
     &    mstu(5)*(k(ist1,5)/mstu(5))+i 
          k(i,5)=mstu(5)*ist1   
          k(i,4)=mstu(5)*ist2   
          IF(k(ist2,5)/mstu(5).EQ.ist1) k(ist2,5)=mstu(5)*i+    
     &    mod(k(ist2,5),mstu(5))    
          IF(mod(k(ist2,4),mstu(5)).EQ.ist1) k(ist2,4)= 
     &    mstu(5)*(k(ist2,4)/mstu(5))+i 
          kstr(istm,2)=i    
          kstr(nstr+1,1)=i  
          kstr(nstr+1,2)=ist2   
  210     nstr=nstr+1   
    
C...String drawing and colour flow for gluon loop.  
        ELSEIF(k(n+1,2).EQ.21) THEN 
          k(n+1,4)=mstu(5)*(n+2)    
          k(n+1,5)=mstu(5)*(n+2)    
          k(n+2,4)=mstu(5)*(n+1)    
          k(n+2,5)=mstu(5)*(n+1)    
          kstr(nstr+1,1)=n+1    
          kstr(nstr+1,2)=n+2    
          kstr(nstr+2,1)=n+2    
          kstr(nstr+2,2)=n+1    
          nstr=nstr+2   
    
C...String drawing and colour flow for q-qbar pair. 
        ELSE    
          k(n+1,4)=mstu(5)*(n+2)    
          k(n+2,5)=mstu(5)*(n+1)    
          kstr(nstr+1,1)=n+1    
          kstr(nstr+1,2)=n+2    
          nstr=nstr+1   
        ENDIF   
    
C...Update remaining energy; iterate.   
        n=n+2   
        IF(n.GT.mstu(4)-mstu(32)-10) THEN   
          CALL luerrm(11,'(PYHIMULT:) no more memory left in LUJETS') 
          IF(mstu(21).GE.1) RETURN  
        ENDIF   
        mint(31)=mint(31)+1 
        vint(151)=vint(151)+vint(41)    
        vint(152)=vint(152)+vint(42)    
        vint(143)=vint(143)-vint(41)    
        vint(144)=vint(144)-vint(42)    
        IF(mint(31).LT.240) goto 180    
  220   CONTINUE    
      ENDIF 
    
C...Format statements for printout. 
 1000 FORMAT(/1x,'****** PYHIMULT: initialization of multiple inter', 
     &'actions for MSTP(82) =',i2,' ******')    
 1100 FORMAT(8x,'pT0 =',f5.2,' GeV gives sigma(parton-parton) =',1p,    
     &e9.2,' mb: rejected') 
 1200 FORMAT(8x,'pT0 =',f5.2,' GeV gives sigma(parton-parton) =',1p,    
     &e9.2,' mb: accepted') 
    
      RETURN    
      END