pyveto.f

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

 
C*********************************************************************
 
C...PYVETO
C...Interface to UPVETO, which allows user to veto event generation
C...on the parton level, after parton showers but before multiple
C...interactions, beam remnants and hadronization is added.
 
      SUBROUTINE pyveto(IVETO)
 
C...All real arithmetic in double precision.
      IMPLICIT DOUBLE PRECISION(a-h, o-z)
C...Three Pythia functions return integers, so need declaring.
      INTEGER pyk,pychge,pycomp
 
C...PYTHIA commonblocks.
      common/pyjets/n,npad,k(4000,5),p(4000,5),v(4000,5)
      common/pypars/mstp(200),parp(200),msti(200),pari(200)
      common/pydat1/mstu(200),paru(200),mstj(200),parj(200)
      common/pyint1/mint(400),vint(400)
      SAVE /pyjets/,/pypars/,/pyint1/
C...HEPEVT commonblock.
      parameter(nmxhep=4000)
      common/hepevt/nevhep,nhep,isthep(nmxhep),idhep(nmxhep),
     &jmohep(2,nmxhep),jdahep(2,nmxhep),phep(5,nmxhep),vhep(4,nmxhep)
      DOUBLE PRECISION phep,vhep
      SAVE /hepevt/
C...Local array.
      dimension ireso(100)
 
C...Define longitudinal boost from initiator rest frame to cm frame.
      gamma=0.5d0*(vint(141)+vint(142))/sqrt(vint(141)*vint(142))
      gabez=0.5d0*(vint(141)-vint(142))/sqrt(vint(141)*vint(142))
 
C... Reset counters.
      nevhep=0
      nhep=0
      nreso=0
      
C...Oth pass: identify beam and incoming partons
      DO 140 i=mint(83)+1,mint(83)+6
        istore=0
        IF(k(i,2).EQ.94) THEN

        ELSE
          istore=1
          nhep=nhep+1
          ii=nhep
          nreso=nreso+1
          ireso(nreso)=i
          imoth=k(i,3)
        ENDIF
        IF(istore.EQ.1) THEN
C...Copy parton info, boosting momenta along z axis to cm frame.
          isthep(ii)=2
          idhep(ii)=k(i,2)
          phep(1,ii)=p(i,1)
          phep(2,ii)=p(i,2)
          phep(3,ii)=gamma*p(i,3)+gabez*p(i,4)
          phep(4,ii)=gamma*p(i,4)+gabez*p(i,3)
          phep(5,ii)=p(i,5)
C...Store one mother. Rest of history and vertex info zeroed.
          jmohep(1,ii)=imoth
          jmohep(2,ii)=0
          jdahep(1,ii)=0
          jdahep(2,ii)=0
          vhep(1,ii)=0d0
          vhep(2,ii)=0d0
          vhep(3,ii)=0d0
          vhep(4,ii)=0d0
        ENDIF
 140  CONTINUE

C...First pass: identify final locations of resonances
C...and of their daughters before showering.
      DO 150 i=mint(84)+3,n
        istore=0
        imoth=0
 
C...Skip shower CM frame documentation lines.
        IF(k(i,2).EQ.94) THEN
 
C...  Store a new intermediate product, when mother in documentation.
        ELSEIF(mstp(128).EQ.0.AND.k(i,3).GT.mint(83)+6.AND.
     &  k(i,3).LE.mint(84)) THEN
          istore=1
          nhep=nhep+1
          ii=nhep
          nreso=nreso+1
          ireso(nreso)=i
          imoth=k(k(i,3),3)
 
C...  Store a new intermediate product, when mother in main section.
        ELSEIF(mstp(128).EQ.1.AND.k(i-mint(84)+mint(83)+4,1).EQ.21.AND.
     &  k(i-mint(84)+mint(83)+4,2).EQ.k(i,2)) THEN
          istore=1
          nhep=nhep+1
          ii=nhep
          nreso=nreso+1
          ireso(nreso)=i
          imoth=max(0,k(i-mint(84)+mint(83)+4,3))
        ENDIF
  
        IF(istore.EQ.1) THEN
C...Copy parton info, boosting momenta along z axis to cm frame.
          isthep(ii)=2
          idhep(ii)=k(i,2)
          phep(1,ii)=p(i,1)
          phep(2,ii)=p(i,2)
          phep(3,ii)=gamma*p(i,3)+gabez*p(i,4)
          phep(4,ii)=gamma*p(i,4)+gabez*p(i,3)
          phep(5,ii)=p(i,5)
C...Store one mother. Rest of history and vertex info zeroed.
          jmohep(1,ii)=imoth
          jmohep(2,ii)=0
          jdahep(1,ii)=i
          jdahep(2,ii)=0
          vhep(1,ii)=0d0
          vhep(2,ii)=0d0
          vhep(3,ii)=0d0
          vhep(4,ii)=0d0
        ENDIF
 150  CONTINUE

C...Second pass: identify current set of "final" partons.
      DO 200 i=mint(84)+3,n
        istore=0
        imoth=0
 
C...Store a final parton.
        IF(k(i,1).GE.1.AND.k(i,1).LE.10) THEN
          istore=1
          nhep=nhep+1
          ii=nhep
C..Trace it back through shower, to check if from documented particle.
          ihist=i
          isave=ihist
  160     CONTINUE
          IF(ihist.GT.mint(84)) THEN
            IF(k(ihist,2).EQ.94) ihist=k(ihist,3)+(isave-1-ihist)
            DO 170 iri=1,nreso
              IF(ihist.EQ.ireso(iri)) imoth=iri
  170       CONTINUE
            isave=ihist
            ihist=k(ihist,3)
            IF(imoth.EQ.0) goto 160
          ELSEIF(ihist.LE.4) THEN
            IF(ihist.EQ.1.OR.ihist.EQ.2) THEN
              istore=0
              nhep=nhep-1
            ELSE
              imoth=ihist
            ENDIF
          ENDIF
        ENDIF
 
        IF(istore.EQ.1) THEN
C...Copy parton info, boosting momenta along z axis to cm frame.
          isthep(ii)=1
          idhep(ii)=k(i,2)
          phep(1,ii)=p(i,1)
          phep(2,ii)=p(i,2)
          phep(3,ii)=gamma*p(i,3)+gabez*p(i,4)
          phep(4,ii)=gamma*p(i,4)+gabez*p(i,3)
          phep(5,ii)=p(i,5)
C...Store one mother. Rest of history and vertex info zeroed.
          jmohep(1,ii)=imoth
          jmohep(2,ii)=0
          jdahep(1,ii)=0
          jdahep(2,ii)=0
          vhep(1,ii)=0d0
          vhep(2,ii)=0d0
          vhep(3,ii)=0d0
          vhep(4,ii)=0d0
        ENDIF
  200 CONTINUE

C...Call user-written routine to decide whether to keep events.
      CALL upveto(iveto)
 
      RETURN
      END