pypile.f

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

 
C*********************************************************************
 
C...PYPILE
C...Initializes multiplicity distribution and selects mutliplicity
C...of pileup events, i.e. several events occuring at the same
C...beam crossing.
 
      SUBROUTINE pypile(MPILE)
 
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/pypars/mstp(200),parp(200),msti(200),pari(200)
      common/pyint1/mint(400),vint(400)
      common/pyint7/sigt(0:6,0:6,0:5)
      SAVE /pydat1/,/pypars/,/pyint1/,/pyint7/
C...Local arrays and saved variables.
      dimension wti(0:200)
      SAVE imin,imax,wti,wts
 
C...Sum of allowed cross-sections for pileup events.
      IF(mpile.EQ.1) THEN
        vint(131)=sigt(0,0,5)
        IF(mstp(132).GE.2) vint(131)=vint(131)+sigt(0,0,4)
        IF(mstp(132).GE.3) vint(131)=vint(131)+sigt(0,0,2)+sigt(0,0,3)
        IF(mstp(132).GE.4) vint(131)=vint(131)+sigt(0,0,1)
        IF(mstp(133).LE.0) RETURN
 
C...Initialize multiplicity distribution at maximum.
        xnave=vint(131)*parp(131)
        IF(xnave.GT.120d0) WRITE(mstu(11),5000) xnave
        inave=max(1,min(200,nint(xnave)))
        wti(inave)=1d0
        wts=wti(inave)
        wtn=wti(inave)*inave
 
C...Find shape of multiplicity distribution below maximum.
        imin=inave
        DO 100 i=inave-1,1,-1
          IF(mstp(133).EQ.1) wti(i)=wti(i+1)*(i+1)/xnave
          IF(mstp(133).GE.2) wti(i)=wti(i+1)*i/xnave
          IF(wti(i).LT.1d-6) goto 110
          wts=wts+wti(i)
          wtn=wtn+wti(i)*i
          imin=i
  100   CONTINUE
 
C...Find shape of multiplicity distribution above maximum.
  110   imax=inave
        DO 120 i=inave+1,200
          IF(mstp(133).EQ.1) wti(i)=wti(i-1)*xnave/i
          IF(mstp(133).GE.2) wti(i)=wti(i-1)*xnave/(i-1)
          IF(wti(i).LT.1d-6) goto 130
          wts=wts+wti(i)
          wtn=wtn+wti(i)*i
          imax=i
  120   CONTINUE
  130   vint(132)=xnave
        vint(133)=wtn/wts
        IF(mstp(133).EQ.1.AND.imin.EQ.1) vint(134)=
     &  wts/(wts+wti(1)/xnave)
        IF(mstp(133).EQ.1.AND.imin.GT.1) vint(134)=1d0
        IF(mstp(133).GE.2) vint(134)=xnave
 
C...Pick multiplicity of pileup events.
      ELSE
        IF(mstp(133).LE.0) THEN
          mint(81)=max(1,mstp(134))
        ELSE
          wtr=wts*pyr(0)
          DO 140 i=imin,imax
            mint(81)=i
            wtr=wtr-wti(i)
            IF(wtr.LE.0d0) goto 150
  140     CONTINUE
  150     CONTINUE
        ENDIF
      ENDIF
 
C...Format statement for error message.
 5000 FORMAT(1x,'Warning: requested average number of events per bunch',
     &'crossing too large, ',1p,d12.4)
 
      RETURN
      END