pyhirand.f

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

    
C*********************************************************************  
    
      SUBROUTINE pyhirand 
    
C...Generates quantities characterizing the high-pT scattering at the   
C...parton level according to the matrix elements. Chooses incoming,    
C...reacting partons, their momentum fractions and one of the possible  
C...subprocesses.   
      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/pyhiint4/widp(21:40,0:40),wide(21:40,0:40),wids(21:40,3) 
      SAVE /pyhiint4/ 
      common/pyhiint5/ngen(0:200,3),xsec(0:200,3) 
      SAVE /pyhiint5/ 
    
C...Initial values, specifically for (first) semihard interaction.  
      mint(17)=0    
      mint(18)=0    
      vint(143)=1.  
      vint(144)=1.  
      IF(msub(95).EQ.1.OR.mint(82).GE.2) CALL pyhimult(2) 
      isub=0    
  100 mint(51)=0    
    
C...Choice of process type - first event of overlay.    
      IF(mint(82).EQ.1.AND.(isub.LE.90.OR.isub.GT.96)) THEN 
        rsub=xsec(0,1)*rlu(0)   
        DO 110 i=1,200  
        IF(msub(i).NE.1) goto 110   
        isub=i  
        rsub=rsub-xsec(i,1) 
        IF(rsub.LE.0.) goto 120 
  110   CONTINUE    
  120   IF(isub.EQ.95) isub=96  
    
C...Choice of inclusive process type - overlayed events.    
      ELSEIF(mint(82).GE.2.AND.isub.EQ.0) THEN  
        rsub=vint(131)*rlu(0)   
        isub=96 
        IF(rsub.GT.vint(106)) isub=93   
        IF(rsub.GT.vint(106)+vint(104)) isub=92 
        IF(rsub.GT.vint(106)+vint(104)+vint(103)) isub=91   
      ENDIF 
      IF(mint(82).EQ.1) ngen(0,1)=ngen(0,1)+1   
      IF(mint(82).EQ.1) ngen(isub,1)=ngen(isub,1)+1 
      mint(1)=isub  
    
C...Find resonances (explicit or implicit in cross-section).    
      mint(72)=0    
      kfr1=0    
      IF(iset(isub).EQ.1.OR.iset(isub).EQ.3) THEN   
        kfr1=kfpr(isub,1)   
      ELSEIF(isub.GE.71.AND.isub.LE.77) THEN    
        kfr1=25 
      ENDIF 
      IF(kfr1.NE.0) THEN    
        taur1=pmas(kfr1,1)**2/vint(2)   
        gamr1=pmas(kfr1,1)*pmas(kfr1,2)/vint(2) 
        mint(72)=1  
        mint(73)=kfr1   
        vint(73)=taur1  
        vint(74)=gamr1  
      ENDIF 
      IF(isub.EQ.141) THEN  
        kfr2=23 
        taur2=pmas(kfr2,1)**2/vint(2)   
        gamr2=pmas(kfr2,1)*pmas(kfr2,2)/vint(2) 
        mint(72)=2  
        mint(74)=kfr2   
        vint(75)=taur2  
        vint(76)=gamr2  
      ENDIF 
    
C...Find product masses and minimum pT of process,  
C...optionally with broadening according to a truncated Breit-Wigner.   
      vint(63)=0.   
      vint(64)=0.   
      mint(71)=0    
      vint(71)=ckin(3)  
      IF(mint(82).GE.2) vint(71)=0. 
      IF(iset(isub).EQ.2.OR.iset(isub).EQ.4) THEN   
        DO 130 i=1,2    
        IF(kfpr(isub,i).EQ.0) THEN  
        ELSEIF(mstp(42).LE.0) THEN  
          vint(62+i)=pmas(kfpr(isub,i),1)**2    
        ELSE    
          vint(62+i)=ulmass(kfpr(isub,i))**2    
        ENDIF   
  130   CONTINUE    
        IF(min(vint(63),vint(64)).LT.ckin(6)**2) mint(71)=1 
        IF(mint(71).EQ.1) vint(71)=max(ckin(3),ckin(5)) 
      ENDIF 
    
      IF(iset(isub).EQ.0) THEN  
C...Double or single diffractive, or elastic scattering:    
C...choose m^2 according to 1/m^2 (diffractive), constant (elastic) 
        is=int(1.5+rlu(0))  
        vint(63)=vint(3)**2 
        vint(64)=vint(4)**2 
        IF(isub.EQ.92.OR.isub.EQ.93) vint(62+is)=parp(111)**2   
        IF(isub.EQ.93) vint(65-is)=parp(111)**2 
        sh=vint(2)  
        sqm1=vint(3)**2 
        sqm2=vint(4)**2 
        sqm3=vint(63)   
        sqm4=vint(64)   
        sqla12=(sh-sqm1-sqm2)**2-4.*sqm1*sqm2   
        sqla34=(sh-sqm3-sqm4)**2-4.*sqm3*sqm4   
        thter1=sqm1+sqm2+sqm3+sqm4-(sqm1-sqm2)*(sqm3-sqm4)/sh-sh    
        thter2=sqrt(max(0.,sqla12))*sqrt(max(0.,sqla34))/sh 
        thl=0.5*(thter1-thter2) 
        thu=0.5*(thter1+thter2) 
        thm=min(max(thl,parp(101)),thu) 
        jtmax=0 
        IF(isub.EQ.92.OR.isub.EQ.93) jtmax=isub-91  
        DO 140 jt=1,jtmax   
        mint(13+3*jt-is*(2*jt-3))=1 
        sqmmin=vint(59+3*jt-is*(2*jt-3))    
        sqmi=vint(8-3*jt+is*(2*jt-3))**2    
        sqmj=vint(3*jt-1-is*(2*jt-3))**2    
        sqmf=vint(68-3*jt+is*(2*jt-3))  
        squa=0.5*sh/sqmi*((1.+(sqmi-sqmj)/sh)*thm+sqmi-sqmf-    
     &  sqmj**2/sh+(sqmi+sqmj)*sqmf/sh+(sqmi-sqmj)**2/sh**2*sqmf)   
        quar=sh/sqmi*(thm*(thm+sh-sqmi-sqmj-sqmf*(1.-(sqmi-sqmj)/sh))+  
     &  sqmi*sqmj-sqmj*sqmf*(1.+(sqmi-sqmj-sqmf)/sh))   
        sqmmax=squa+sqrt(max(0.,squa**2-quar))  
        IF(abs(quar/squa**2).LT.1.e-06) sqmmax=0.5*quar/squa    
        sqmmax=min(sqmmax,(vint(1)-sqrt(sqmf))**2)  
        vint(59+3*jt-is*(2*jt-3))=sqmmin*(sqmmax/sqmmin)**rlu(0)    
  140   CONTINUE    
C...Choose t-hat according to exp(B*t-hat+C*t-hat^2).   
        sqm3=vint(63)   
        sqm4=vint(64)   
        sqla34=(sh-sqm3-sqm4)**2-4.*sqm3*sqm4   
        thter1=sqm1+sqm2+sqm3+sqm4-(sqm1-sqm2)*(sqm3-sqm4)/sh-sh    
        thter2=sqrt(max(0.,sqla12))*sqrt(max(0.,sqla34))/sh 
        thl=0.5*(thter1-thter2) 
        thu=0.5*(thter1+thter2) 
        b=vint(121) 
        c=vint(122) 
        IF(isub.EQ.92.OR.isub.EQ.93) THEN   
          b=0.5*b   
          c=0.5*c   
        ENDIF   
        thm=min(max(thl,parp(101)),thu) 
        expth=0.    
        tharg=b*(thm-thu)   
        IF(tharg.GT.-20.) expth=exp(tharg)  
  150   th=thu+log(expth+(1.-expth)*rlu(0))/b   
        th=max(thm,min(thu,th)) 
        ratlog=min((b+c*(th+thm))*(th-thm),(b+c*(th+thu))*(th-thu)) 
        IF(ratlog.LT.log(rlu(0))) goto 150  
        vint(21)=1. 
        vint(22)=0. 
        vint(23)=min(1.,max(-1.,(2.*th-thter1)/thter2)) 
    
C...Note: in the following, by In is meant the integral over the    
C...quantity multiplying coefficient cn.    
C...Choose tau according to h1(tau)/tau, where  
C...h1(tau) = c0 + I0/I1*c1*1/tau + I0/I2*c2*1/(tau+tau_R) +    
C...I0/I3*c3*tau/((s*tau-m^2)^2+(m*Gamma)^2) +  
C...I0/I4*c4*1/(tau+tau_R') +   
C...I0/I5*c5*tau/((s*tau-m'^2)^2+(m'*Gamma')^2), and    
C...c0 + c1 + c2 + c3 + c4 + c5 = 1 
      ELSEIF(iset(isub).GE.1.AND.iset(isub).LE.4) THEN  
        CALL pyhiklim(1)  
        IF(mint(51).NE.0) goto 100  
        rtau=rlu(0) 
        mtau=1  
        IF(rtau.GT.coef(isub,1)) mtau=2 
        IF(rtau.GT.coef(isub,1)+coef(isub,2)) mtau=3    
        IF(rtau.GT.coef(isub,1)+coef(isub,2)+coef(isub,3)) mtau=4   
        IF(rtau.GT.coef(isub,1)+coef(isub,2)+coef(isub,3)+coef(isub,4)) 
     &  mtau=5  
        IF(rtau.GT.coef(isub,1)+coef(isub,2)+coef(isub,3)+coef(isub,4)+ 
     &  coef(isub,5)) mtau=6    
        CALL pyhikmap(1,mtau,rlu(0))  
    
C...2 -> 3, 4 processes:    
C...Choose tau' according to h4(tau,tau')/tau', where   
C...h4(tau,tau') = c0 + I0/I1*c1*(1 - tau/tau')^3/tau', and 
C...c0 + c1 = 1.    
        IF(iset(isub).EQ.3.OR.iset(isub).EQ.4) THEN 
          CALL pyhiklim(4)    
          IF(mint(51).NE.0) goto 100    
          rtaup=rlu(0)  
          mtaup=1   
          IF(rtaup.GT.coef(isub,15)) mtaup=2    
          CALL pyhikmap(4,mtaup,rlu(0))   
        ENDIF   
    
C...Choose y* according to h2(y*), where    
C...h2(y*) = I0/I1*c1*(y*-y*min) + I0/I2*c2*(y*max-y*) +    
C...I0/I3*c3*1/cosh(y*), I0 = y*max-y*min, and c1 + c2 + c3 = 1.    
        CALL pyhiklim(2)  
        IF(mint(51).NE.0) goto 100  
        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))  
    
C...2 -> 2 processes:   
C...Choose cos(theta-hat) (cth) according to h3(cth), where 
C...h3(cth) = c0 + I0/I1*c1*1/(A - cth) + I0/I2*c2*1/(A + cth) +    
C...I0/I3*c3*1/(A - cth)^2 + I0/I4*c4*1/(A + cth)^2,    
C...A = 1 + 2*(m3*m4/sh)^2 (= 1 for massless products), 
C...and c0 + c1 + c2 + c3 + c4 = 1. 
        CALL pyhiklim(3)  
        IF(mint(51).NE.0) goto 100  
        IF(iset(isub).EQ.2.OR.iset(isub).EQ.4) THEN 
          rcth=rlu(0)   
          mcth=1    
          IF(rcth.GT.coef(isub,10)) mcth=2  
          IF(rcth.GT.coef(isub,10)+coef(isub,11)) mcth=3    
          IF(rcth.GT.coef(isub,10)+coef(isub,11)+coef(isub,12)) mcth=4  
          IF(rcth.GT.coef(isub,10)+coef(isub,11)+coef(isub,12)+ 
     &    coef(isub,13)) mcth=5 
          CALL pyhikmap(3,mcth,rlu(0))    
        ENDIF   
    
C...Low-pT or multiple interactions (first semihard interaction).   
      ELSEIF(iset(isub).EQ.5) THEN  
        CALL pyhimult(3)  
        isub=mint(1)    
      ENDIF 
    
C...Choose azimuthal angle. 
      vint(24)=paru(2)*rlu(0)   
    
C...Check against user cuts on kinematics at parton level.  
      mint(51)=0    
      IF(isub.LE.90.OR.isub.GT.100) CALL pyhiklim(0)  
      IF(mint(51).NE.0) goto 100    
      IF(mint(82).EQ.1.AND.mstp(141).GE.1) THEN 
        mcut=0  
        IF(msub(91)+msub(92)+msub(93)+msub(94)+msub(95).EQ.0)   
     &  CALL pyhikcut(mcut)   
        IF(mcut.NE.0) goto 100  
      ENDIF 
    
C...Calculate differential cross-section for different subprocesses.    
      CALL pyhisigh(nchn,sigs)    
    
C...Calculations for Monte Carlo estimate of all cross-sections.    
      IF(mint(82).EQ.1.AND.isub.LE.90.OR.isub.GE.96) THEN   
        xsec(isub,2)=xsec(isub,2)+sigs  
      ELSEIF(mint(82).EQ.1) THEN    
        xsec(isub,2)=xsec(isub,2)+xsec(isub,1)  
      ENDIF 
    
C...Multiple interactions: store results of cross-section calculation.  
      IF(mint(43).EQ.4.AND.mstp(82).GE.3) THEN  
        vint(153)=sigs  
        CALL pyhimult(4)  
      ENDIF 
    
C...Weighting using estimate of maximum of differential cross-section.  
      viol=sigs/xsec(isub,1)    
      IF(viol.LT.rlu(0)) goto 100   
    
C...Check for possible violation of estimated maximum of differential   
C...cross-section used in weighting.    
      IF(mstp(123).LE.0) THEN   
        IF(viol.GT.1.) THEN 
          WRITE(mstu(11),1000) viol,ngen(0,3)+1 
          WRITE(mstu(11),1100) isub,vint(21),vint(22),vint(23),vint(26) 
          stop  
        ENDIF   
      ELSEIF(mstp(123).EQ.1) THEN   
        IF(viol.GT.vint(108)) THEN  
          vint(108)=viol    
C          IF(VIOL.GT.1.) THEN   
C            WRITE(MSTU(11),1200) VIOL,NGEN(0,3)+1   
C            WRITE(MSTU(11),1100) ISUB,VINT(21),VINT(22),VINT(23),   
C     &      VINT(26)    
C          ENDIF 
        ENDIF   
      ELSEIF(viol.GT.vint(108)) THEN    
        vint(108)=viol  
        IF(viol.GT.1.) THEN 
          xdif=xsec(isub,1)*(viol-1.)   
          xsec(isub,1)=xsec(isub,1)+xdif    
          IF(msub(isub).EQ.1.AND.(isub.LE.90.OR.isub.GT.96))    
     &    xsec(0,1)=xsec(0,1)+xdif  
C          WRITE(MSTU(11),1200) VIOL,NGEN(0,3)+1 
C          WRITE(MSTU(11),1100) ISUB,VINT(21),VINT(22),VINT(23),VINT(26) 
C          IF(ISUB.LE.9) THEN    
C            WRITE(MSTU(11),1300) ISUB,XSEC(ISUB,1)  
C          ELSEIF(ISUB.LE.99) THEN   
C            WRITE(MSTU(11),1400) ISUB,XSEC(ISUB,1)  
C          ELSE  
C            WRITE(MSTU(11),1500) ISUB,XSEC(ISUB,1)  
C          ENDIF 
          vint(108)=1.  
        ENDIF   
      ENDIF 
    
C...Multiple interactions: choose impact parameter. 
      vint(148)=1.  
      IF(mint(43).EQ.4.AND.(isub.LE.90.OR.isub.GE.96).AND.mstp(82).GE.3)    
     &THEN  
        CALL pyhimult(5)  
        IF(vint(150).LT.rlu(0)) goto 100    
      ENDIF 
      IF(mint(82).EQ.1.AND.msub(95).EQ.1) THEN  
        IF(isub.LE.90.OR.isub.GE.95) ngen(95,1)=ngen(95,1)+1    
        IF(isub.LE.90.OR.isub.GE.96) ngen(96,2)=ngen(96,2)+1    
      ENDIF 
      IF(isub.LE.90.OR.isub.GE.96) mint(31)=mint(31)+1  
    
C...Choose flavour of reacting partons (and subprocess).    
      rsigs=sigs*rlu(0) 
      qt2=vint(48)  
      rqqbar=parp(87)*(1.-(qt2/(qt2+(parp(88)*parp(82))**2))**2)    
      IF(isub.NE.95.AND.(isub.NE.96.OR.mstp(82).LE.1.OR.    
     &rlu(0).GT.rqqbar)) THEN   
        DO 190 ichn=1,nchn  
        kfl1=isig(ichn,1)   
        kfl2=isig(ichn,2)   
        mint(2)=isig(ichn,3)    
        rsigs=rsigs-sigh(ichn)  
        IF(rsigs.LE.0.) goto 210    
  190   CONTINUE    
    
C...Multiple interactions: choose qqbar preferentially at small pT. 
      ELSEIF(isub.EQ.96) THEN   
        CALL pyhispli(mint(11),21,kfl1,kfldum)    
        CALL pyhispli(mint(12),21,kfl2,kfldum)    
        mint(1)=11  
        mint(2)=1   
        IF(kfl1.EQ.kfl2.AND.rlu(0).LT.0.5) mint(2)=2    
    
C...Low-pT: choose string drawing configuration.    
      ELSE  
        kfl1=21 
        kfl2=21 
        rsigs=6.*rlu(0) 
        mint(2)=1   
        IF(rsigs.GT.1.) mint(2)=2   
        IF(rsigs.GT.2.) mint(2)=3   
      ENDIF 
    
C...Reassign QCD process. Partons before initial state radiation.   
  210 IF(mint(2).GT.10) THEN    
        mint(1)=mint(2)/10  
        mint(2)=mod(mint(2),10) 
      ENDIF 
      mint(15)=kfl1 
      mint(16)=kfl2 
      mint(13)=mint(15) 
      mint(14)=mint(16) 
      vint(141)=vint(41)    
      vint(142)=vint(42)    
    
C...Format statements for differential cross-section maximum violations.    
 1000 FORMAT(1x,'Error: maximum violated by',1p,e11.3,1x,   
     &'in event',1x,i7,'.'/1x,'Execution stopped!') 
 1100 FORMAT(1x,'ISUB = ',i3,'; Point of violation:'/1x,'tau =',1p, 
     &e11.3,', y* =',e11.3,', cthe = ',0p,f11.7,', tau'' =',1p,e11.3)   
 1200 FORMAT(1x,'Warning: maximum violated by',1p,e11.3,1x, 
     &'in event',1x,i7) 
 1300 FORMAT(1x,'XSEC(',i1,',1) increased to',1p,e11.3) 
 1400 FORMAT(1x,'XSEC(',i2,',1) increased to',1p,e11.3) 
 1500 FORMAT(1x,'XSEC(',i3,',1) increased to',1p,e11.3) 
    
      RETURN    
      END