pyhithia.f

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

    
C*********************************************************************  
    
      SUBROUTINE pyhithia 
    
C...Administers the generation of a high-pt event via calls to a number 
C...of subroutines; also computes cross-sections.   
      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/pyhiint5/ngen(0:200,3),xsec(0:200,3) 
      SAVE /pyhiint5/ 
    
C...Loop over desired number of overlayed events (normally 1).  
      mint(7)=0 
      mint(8)=0 
      novl=1    
      IF(mstp(131).NE.0) CALL pyhiovly(2) 
      IF(mstp(131).NE.0) novl=mint(81)  
      mint(83)=0    
      mint(84)=mstp(126)    
      mstu(70)=0    
      DO 190 iovl=1,novl    
      IF(mint(84)+100.GE.mstu(4)) THEN  
        CALL luerrm(11, 
     &  '(PYHITHIA:) no more space in LUJETS for overlayed events')   
        IF(mstu(21).GE.1) goto 200  
      ENDIF 
      mint(82)=iovl 
    
C...Generate variables of hard scattering.  
  100 CONTINUE  
      IF(iovl.EQ.1) ngen(0,2)=ngen(0,2)+1   
      mint(31)=0    
      mint(51)=0    
      CALL pyhirand   
      isub=mint(1)  
      IF(iovl.EQ.1) THEN    
        ngen(isub,2)=ngen(isub,2)+1 
    
C...Store information on hard interaction.  
        DO 110 j=1,200  
        msti(j)=0   
  110   pari(j)=0.  
        msti(1)=mint(1) 
        msti(2)=mint(2) 
        msti(11)=mint(11)   
        msti(12)=mint(12)   
        msti(15)=mint(15)   
        msti(16)=mint(16)   
        msti(17)=mint(17)   
        msti(18)=mint(18)   
        pari(11)=vint(1)    
        pari(12)=vint(2)    
        IF(isub.NE.95) THEN 
          DO 120 j=13,22    
  120     pari(j)=vint(30+j)    
          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(41)=vint(23) 
        ENDIF   
      ENDIF 
    
      IF(mstp(111).EQ.-1) goto 160  
      IF(isub.LE.90.OR.isub.GE.95) THEN 
C...Hard scattering (including low-pT): 
C...reconstruct kinematics and colour flow of hard scattering.  
        CALL pyhiscat 
        IF(mint(51).EQ.1) goto 100  
    
C...Showering of initial state partons (optional).  
        ipu1=mint(84)+1 
        ipu2=mint(84)+2 
        IF(mstp(61).GE.1.AND.mint(43).NE.1.AND.isub.NE.95)  
     &  CALL pyhisspa(ipu1,ipu2)  
        nsav1=n 
    
C...Multiple interactions.  
        IF(mstp(81).GE.1.AND.mint(43).EQ.4.AND.isub.NE.95)  
     &  CALL pyhimult(6)  
        mint(1)=isub    
        nsav2=n 
    
C...Hadron remnants and primordial kT.  
        CALL pyhiremn(ipu1,ipu2)  
        IF(mint(51).EQ.1) goto 100  
        nsav3=n 
    
C...Showering of final state partons (optional).    
        ipu3=mint(84)+3 
        ipu4=mint(84)+4 
        IF(mstp(71).GE.1.AND.isub.NE.95.AND.k(ipu3,1).GT.0.AND. 
     &  k(ipu3,1).LE.10.AND.k(ipu4,1).GT.0.AND.k(ipu4,1).LE.10) THEN    
          qmax=sqrt(parp(71)*vint(52))  
          IF(isub.EQ.5) qmax=sqrt(pmas(23,1)**2)    
          IF(isub.EQ.8) qmax=sqrt(pmas(24,1)**2)    
          CALL lushow(ipu3,ipu4,qmax)   
        ENDIF   
    
C...Sum up transverse and longitudinal momenta. 
        IF(iovl.EQ.1) THEN  
          pari(65)=2.*pari(17)  
          DO 130 i=mstp(126)+1,n    
          IF(k(i,1).LE.0.OR.k(i,1).GT.10) goto 130  
          pt=sqrt(p(i,1)**2+p(i,2)**2)  
          pari(69)=pari(69)+pt  
          IF(i.LE.nsav1.OR.i.GT.nsav3) pari(66)=pari(66)+pt 
          IF(i.GT.nsav1.AND.i.LE.nsav2) pari(68)=pari(68)+pt    
  130     CONTINUE  
          pari(67)=pari(68) 
          pari(71)=vint(151)    
          pari(72)=vint(152)    
          pari(73)=vint(151)    
          pari(74)=vint(152)    
        ENDIF   
    
C...Decay of final state resonances.    
        IF(mstp(41).GE.1.AND.isub.NE.95) CALL pyhiresd    
    
      ELSE  
C...Diffractive and elastic scattering. 
        CALL pyhidiff 
        IF(iovl.EQ.1) THEN  
          pari(65)=2.*pari(17)  
          pari(66)=pari(65) 
          pari(69)=pari(65) 
        ENDIF   
      ENDIF 
    
C...Recalculate energies from momenta and masses (if desired).  
      IF(mstp(113).GE.1) THEN   
        DO 140 i=mint(83)+1,n   
  140   IF(k(i,1).GT.0.AND.k(i,1).LE.10) p(i,4)=sqrt(p(i,1)**2+ 
     &  p(i,2)**2+p(i,3)**2+p(i,5)**2)  
      ENDIF 
    
C...Rearrange partons along strings, check invariant mass cuts. 
      mstu(28)=0    
      CALL luprep(mint(84)+1)   
      IF(mstp(112).EQ.1.AND.mstu(28).EQ.3) goto 100 
      IF(mstp(125).EQ.0.OR.mstp(125).EQ.1) THEN 
        DO 150 i=mint(84)+1,n   
        IF(k(i,2).NE.94) goto 150   
        k(i+1,3)=mod(k(i+1,4)/mstu(5),mstu(5))  
        k(i+2,3)=mod(k(i+2,4)/mstu(5),mstu(5))  
  150   CONTINUE    
        CALL luedit(12) 
        CALL luedit(14) 
        IF(mstp(125).EQ.0) CALL luedit(15)  
        IF(mstp(125).EQ.0) mint(4)=0    
      ENDIF 
    
C...Introduce separators between sections in LULIST event listing.  
      IF(iovl.EQ.1.AND.mstp(125).LE.0) THEN 
        mstu(70)=1  
        mstu(71)=n  
      ELSEIF(iovl.EQ.1) THEN    
        mstu(70)=3  
        mstu(71)=2  
        mstu(72)=mint(4)    
        mstu(73)=n  
      ENDIF 
    
C...Perform hadronization (if desired). 
      IF(mstp(111).GE.1) CALL luexec    
      IF(mstp(125).EQ.0.OR.mstp(125).EQ.1) CALL luedit(14)  
    
C...Calculate Monte Carlo estimates of cross-sections.  
  160 IF(iovl.EQ.1) THEN    
        IF(mstp(111).NE.-1) ngen(isub,3)=ngen(isub,3)+1 
        ngen(0,3)=ngen(0,3)+1   
        xsec(0,3)=0.    
        DO 170 i=1,200  
        IF(i.EQ.96) THEN    
          xsec(i,3)=0.  
        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(1.,float(ngen(96,1))*  
     &    float(ngen(96,2)))    
        ELSEIF(ngen(i,1).EQ.0) THEN 
          xsec(i,3)=0.  
        ELSEIF(ngen(i,2).EQ.0) THEN 
          xsec(i,3)=xsec(i,2)*ngen(0,3)/(float(ngen(i,1))*  
     &    float(ngen(0,2))) 
        ELSE    
          xsec(i,3)=xsec(i,2)*ngen(i,3)/(float(ngen(i,1))*  
     &    float(ngen(i,2))) 
        ENDIF   
  170   xsec(0,3)=xsec(0,3)+xsec(i,3)   
        IF(msub(95).EQ.1) THEN  
          ngens=ngen(91,3)+ngen(92,3)+ngen(93,3)+ngen(94,3)+ngen(95,3)  
          xsecs=xsec(91,3)+xsec(92,3)+xsec(93,3)+xsec(94,3)+xsec(95,3)  
          xmaxs=xsec(95,1)  
          IF(msub(91).EQ.1) xmaxs=xmaxs+xsec(91,1)  
          IF(msub(92).EQ.1) xmaxs=xmaxs+xsec(92,1)  
          IF(msub(93).EQ.1) xmaxs=xmaxs+xsec(93,1)  
          IF(msub(94).EQ.1) xmaxs=xmaxs+xsec(94,1)  
          fac=1.    
          IF(ngens.LT.ngen(0,3)) fac=(xmaxs-xsecs)/(xsec(0,3)-xsecs)    
          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(0,3)=xsec(91,3)+xsec(92,3)+xsec(93,3)+xsec(94,3)+    
     &    xsec(95,1)    
        ENDIF   
    
C...Store final information.    
        mint(5)=mint(5)+1   
        msti(3)=mint(3) 
        msti(4)=mint(4) 
        msti(5)=mint(5) 
        msti(6)=mint(6) 
        msti(7)=mint(7) 
        msti(8)=mint(8) 
        msti(13)=mint(13)   
        msti(14)=mint(14)   
        msti(21)=mint(21)   
        msti(22)=mint(22)   
        msti(23)=mint(23)   
        msti(24)=mint(24)   
        msti(25)=mint(25)   
        msti(26)=mint(26)   
        msti(31)=mint(31)   
        pari(1)=xsec(0,3)   
        pari(2)=xsec(0,3)/mint(5)   
        pari(31)=vint(141)  
        pari(32)=vint(142)  
        IF(isub.NE.95.AND.mint(7)*mint(8).NE.0) THEN    
          pari(42)=2.*vint(47)/vint(1)  
          DO 180 is=7,8 
          pari(36+is)=p(mint(is),3)/vint(1) 
          pari(38+is)=p(mint(is),4)/vint(1) 
          i=mint(is)    
          pr=max(1e-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),1e20)),p(i,3))   
          pr=max(1e-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),1e20)),p(i,3))   
          pari(44+is)=p(i,3)/sqrt(p(i,1)**2+p(i,2)**2+p(i,3)**2)    
          pari(46+is)=ulangl(p(i,3),sqrt(p(i,1)**2+p(i,2)**2))  
          pari(48+is)=ulangl(p(i,1),p(i,2)) 
  180     CONTINUE  
        ENDIF   
        pari(61)=vint(148)  
        IF(iset(isub).EQ.1.OR.iset(isub).EQ.3) THEN 
          mstu(161)=mint(21)    
          mstu(162)=0   
        ELSE    
          mstu(161)=mint(21)    
          mstu(162)=mint(22)    
        ENDIF   
      ENDIF 
    
C...Prepare to go to next overlayed event.  
      msti(41)=iovl 
      IF(iovl.GE.2.AND.iovl.LE.10) msti(40+iovl)=isub   
      IF(mstu(70).LT.10) THEN   
        mstu(70)=mstu(70)+1 
        mstu(70+mstu(70))=n 
      ENDIF 
      mint(83)=n    
      mint(84)=n+mstp(126)  
  190 CONTINUE  
    
C...Information on overlayed events.    
      IF(mstp(131).EQ.1.AND.mstp(133).GE.1) THEN    
        pari(91)=vint(132)  
        pari(92)=vint(133)  
        pari(93)=vint(134)  
        IF(mstp(133).EQ.2) pari(93)=pari(93)*xsec(0,3)/vint(131)    
      ENDIF 
    
C...Transform to the desired coordinate frame.  
  200 CALL pyhifram(mstp(124))    
    
      RETURN    
      END