pyhiinit.f

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

      SUBROUTINE pyhiinit(FRAME,BEAM,TARGET,WIN)  
    
C...Initializes the generation procedure; finds maxima of the   
C...differential cross-sections to be used for weighting.   
      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/ludat3/mdcy(500,3),mdme(2000,2),brat(2000),kfdp(2000,5)    
      SAVE /ludat3/ 
      common/ludat4/chaf(500)   
      CHARACTER chaf*8  
      SAVE /ludat4/ 
      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/ 
      CHARACTER*(*) frame,beam,TARGET   
      CHARACTER chfram*8,chbeam*8,chtarg*8,chmo(12)*3,chlh(2)*6 
      DATA chmo/'Jan','Feb','Mar','Apr','May','Jun','Jul','Aug','Sep',  
     &'Oct','Nov','Dec'/, chlh/'lepton','hadron'/   
    
C...Write headers.  
C      IF(MSTP(122).GE.1) WRITE(MSTU(11),1000) MSTP(181),MSTP(182),  
C     &MSTP(185),CHMO(MSTP(184)),MSTP(183)   
      CALL lulist(0)
C      IF(MSTP(122).GE.1) WRITE(MSTU(11),1100)  
    
C...Identify beam and target particles and initialize kinematics.   
      chfram=frame//' ' 
      chbeam=beam//' '  
      chtarg=TARGET//' '    
      CALL pyhiinki(chfram,chbeam,chtarg,win) 
    
C...Select partonic subprocesses to be included in the simulation.  
      IF(msel.NE.0) THEN    
        DO 100 i=1,200  
  100   msub(i)=0   
      ENDIF 
      IF(mint(43).EQ.1.AND.(msel.EQ.1.OR.msel.EQ.2)) THEN   
C...Lepton+lepton -> gamma/Z0 or W. 
        IF(mint(11)+mint(12).EQ.0) msub(1)=1    
        IF(mint(11)+mint(12).NE.0) msub(2)=1    
      ELSEIF(msel.EQ.1) THEN    
C...High-pT QCD processes:  
        msub(11)=1  
        msub(12)=1  
        msub(13)=1  
        msub(28)=1  
        msub(53)=1  
        msub(68)=1  
        IF(mstp(82).LE.1.AND.ckin(3).LT.parp(81)) msub(95)=1    
        IF(mstp(82).GE.2.AND.ckin(3).LT.parp(82)) msub(95)=1    
      ELSEIF(msel.EQ.2) THEN    
C...All QCD processes:  
        msub(11)=1  
        msub(12)=1  
        msub(13)=1  
        msub(28)=1  
        msub(53)=1  
        msub(68)=1  
        msub(91)=1  
        msub(92)=1  
        msub(93)=1  
        msub(95)=1  
      ELSEIF(msel.GE.4.AND.msel.LE.8) THEN  
C...Heavy quark production. 
        msub(81)=1  
        msub(82)=1  
        DO 110 j=1,min(8,mdcy(21,3))    
  110   mdme(mdcy(21,2)+j-1,1)=0    
        mdme(mdcy(21,2)+msel-1,1)=1 
      ELSEIF(msel.EQ.10) THEN   
C...Prompt photon production:   
        msub(14)=1  
        msub(18)=1  
        msub(29)=1  
      ELSEIF(msel.EQ.11) THEN   
C...Z0/gamma* production:   
        msub(1)=1   
      ELSEIF(msel.EQ.12) THEN   
C...W+/- production:    
        msub(2)=1   
      ELSEIF(msel.EQ.13) THEN   
C...Z0 + jet:   
        msub(15)=1  
        msub(30)=1  
      ELSEIF(msel.EQ.14) THEN   
C...W+/- + jet: 
        msub(16)=1  
        msub(31)=1  
      ELSEIF(msel.EQ.15) THEN   
C...Z0 & W+/- pair production:  
        msub(19)=1  
        msub(20)=1  
        msub(22)=1  
        msub(23)=1  
        msub(25)=1  
      ELSEIF(msel.EQ.16) THEN   
C...H0 production:  
        msub(3)=1   
        msub(5)=1   
        msub(8)=1   
        msub(102)=1 
      ELSEIF(msel.EQ.17) THEN   
C...H0 & Z0 or W+/- pair production:    
        msub(24)=1  
        msub(26)=1  
      ELSEIF(msel.EQ.21) THEN   
C...Z'0 production: 
        msub(141)=1 
      ELSEIF(msel.EQ.22) THEN   
C...H+/- production:    
        msub(142)=1 
      ELSEIF(msel.EQ.23) THEN   
C...R production:   
        msub(143)=1 
      ENDIF 
    
C...Count number of subprocesses on.    
      mint(44)=0    
      DO 120 isub=1,200 
      IF(mint(43).LT.4.AND.isub.GE.91.AND.isub.LE.96.AND.   
     &msub(isub).EQ.1) THEN 
        WRITE(mstu(11),1200) isub,chlh(mint(41)),chlh(mint(42)) 
        stop    
      ELSEIF(msub(isub).EQ.1.AND.iset(isub).EQ.-1) THEN 
        WRITE(mstu(11),1300) isub   
        stop    
      ELSEIF(msub(isub).EQ.1.AND.iset(isub).LE.-2) THEN 
        WRITE(mstu(11),1400) isub   
        stop    
      ELSEIF(msub(isub).EQ.1) THEN  
        mint(44)=mint(44)+1 
      ENDIF 
  120 CONTINUE  
      IF(mint(44).EQ.0) THEN    
        WRITE(mstu(11),1500)    
        stop    
      ENDIF 
      mint(45)=mint(44)-msub(91)-msub(92)-msub(93)-msub(94) 
    
C...Maximum 4 generations; set maximum number of allowed flavours.  
      mstp(1)=min(4,mstp(1))    
      mstu(114)=min(mstu(114),2*mstp(1))    
      mstp(54)=min(mstp(54),2*mstp(1))  
    
C...Sum up Cabibbo-Kobayashi-Maskawa factors for each quark/lepton. 
      DO 140 i=-20,20   
      vint(180+i)=0.    
      ia=iabs(i)    
      IF(ia.GE.1.AND.ia.LE.2*mstp(1)) THEN  
        DO 130 j=1,mstp(1)  
        ib=2*j-1+mod(ia,2)  
        ipm=(5-isign(1,i))/2    
        idc=j+mdcy(ia,2)+2  
  130   IF(mdme(idc,1).EQ.1.OR.mdme(idc,1).EQ.ipm) vint(180+i)= 
     &  vint(180+i)+vckm((ia+1)/2,(ib+1)/2) 
      ELSEIF(ia.GE.11.AND.ia.LE.10+2*mstp(1)) THEN  
        vint(180+i)=1.  
      ENDIF 
  140 CONTINUE  
    
C...Choose Lambda value to use in alpha-strong. 
      mstu(111)=mstp(2) 
      IF(mstp(3).GE.1) THEN 
        alam=parp(1)    
        IF(mstp(51).EQ.1) alam=0.2  
        IF(mstp(51).EQ.2) alam=0.29 
        IF(mstp(51).EQ.3) alam=0.2  
        IF(mstp(51).EQ.4) alam=0.4  
        IF(mstp(51).EQ.11) alam=0.16    
        IF(mstp(51).EQ.12) alam=0.26    
        IF(mstp(51).EQ.13) alam=0.36    
        parp(1)=alam    
        parp(61)=alam   
        paru(112)=alam  
        parj(81)=alam   
      ENDIF 
    
C...Initialize widths and partial widths for resonances.    
      CALL pyhiinre   
    
C...Reset variables for cross-section calculation.  
      DO 150 i=0,200    
      DO 150 j=1,3  
      ngen(i,j)=0   
  150 xsec(i,j)=0.  
      vint(108)=0.  
    
C...Find parametrized total cross-sections. 
      IF(mint(43).EQ.4) CALL pyhixtot 
    
C...Maxima of differential cross-sections.  
      IF(mstp(121).LE.0) CALL pyhimaxi    
    
C...Initialize possibility of overlayed events. 
      IF(mstp(131).NE.0) CALL pyhiovly(1) 
    
C...Initialize multiple interactions with variable impact parameter.    
      IF(mint(43).EQ.4.AND.(mint(45).NE.0.OR.mstp(131).NE.0).AND.   
     &mstp(82).GE.2) CALL pyhimult(1) 
C      IF(MSTP(122).GE.1) WRITE(MSTU(11),1600)  
    
C...Formats for initialization information. 
 1000 FORMAT(///20x,'The Lund Monte Carlo - PYHITHIA version ',i1,
     &'.',i1/ 
     &20x,'**  Last date of change:  ',i2,1x,a3,1x,i4,'  **'/)  
 1100 FORMAT('1',18('*'),1x,'PYHIINIT: initialization of PYHITHIA ',    
     &'(hijing pythia) routines',1x,17('*'))    
 1200 FORMAT(1x,'Error: process number ',i3,' not meaningful for ',a6,  
     &'-',a6,' interactions.'/1x,'Execution stopped!')  
 1300 FORMAT(1x,'Error: requested subprocess',i4,' not implemented.'/   
     &1x,'Execution stopped!')  
 1400 FORMAT(1x,'Error: requested subprocess',i4,' not existing.'/  
     &1x,'Execution stopped!')  
 1500 FORMAT(1x,'Error: no subprocess switched on.'/    
     &1x,'Execution stopped.')  
 1600 FORMAT(/1x,22('*'),1x,'PYHIINIT: initialization completed',1x,  
     &22('*'))  
    
      RETURN    
      END