dpmjetint.f

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

*
*===program crint======================================================*
*
C      OPTIONS/ EXTEND_SOURCE
C      SUBROUTINE CRINT
      SUBROUTINE dt_produceevent(ENERGY_SL, NPARTICLES)


      IMPLICIT DOUBLE PRECISION (a-h,o-z)
      REAL energy_sl
      INTEGER init
      REAL ne,etest,prob,slump
      SAVE

* Call the init sub routine in the first event
      DATA init /0/

      parameter(nmxhkk=200000)

      COMMON /dtiont/ linp,lout,ldat

      COMMON /dtevt1/ nhkk,nevhkk,isthkk(nmxhkk),idhkk(nmxhkk),
     &                jmohkk(2,nmxhkk),jdahkk(2,nmxhkk),
     &                phkk(5,nmxhkk),vhkk(4,nmxhkk),whkk(4,nmxhkk)

*     event flag
      COMMON /dtevno/ nevent, icasca

      IF(init.EQ.0) THEN
         OPEN (unit = 50, file = "my.input")    
   linp = 50
         CALL dt_dtuini(nevts,epn,npmass,npchar,ntmass,ntchar,idp,iemu)
*        Init called, make sure it's not called again
         init = 1
      ENDIF
*-----------------------------------------------------------------------
*     generation of one event
      nevent = 1
      kkmat = -1

*   If an energy-range has been defined with the ENERGY input-card the
*   laboratory energy ELAB can be set to any value within that range,..
C        ELAB = DT_RNDM(EPN)*(EPN-0.5D7)+0.5D7

*   ..otherwise it has to coincide with EPN.
C        ELAB = EPN

      elab = energy_sl

*   sampling of one event

*     TEST

      CALL dt_kkinc(npmass,npchar,ntmass,ntchar,idp,elab,kkmat,irej)

      IF (irej.NE.0) RETURN

c     Return the number of particles produced
      
c     Fill the particle info 
      CALL dt_getparticles(nparticles)

      END


      SUBROUTINE dt_getparticles(NPARTICLES)

      IMPLICIT DOUBLE PRECISION (a-h,o-z)
      INTEGER pid,qch,q_sum,ntpc,nfinal,naccept,ipart,res
      DOUBLE PRECISION yrap,pt,mass,mt,etot
      DOUBLE PRECISION pt_cut_tpc
      parameter(pt_cut_tpc=0.050)

      SAVE
*
* COMMON /DTEVT1/ :
*                   NHKK         number of entries in common block
*                   NEVHKK       number of the event
*                   ISTHKK(i)    status code for entry i
*                   IDHKK(i)     identifier for the entry
*                                (for particles: identifier according
*                                 to the PDG numbering scheme)
*                   JMOHKK(1,i)  pointer to the entry of the first mother
*                                of entry i
*                   JMOHKK(2,i)  pointer to the entry of the second mother
*                                of entry i
*                   JDAHKK(1,i)  pointer to the entry of the first daughter
*                                of entry i
*                   JDAHKK(2,i)  pointer to the entry of the second daughter
*                                of entry i
*                   PHKK(1..3,i) 3-momentum
*                   PHKK(4,i)    energy
*                   PHKK(5,i)    mass
*
* event history

      parameter(nmxhkk=200000)

      COMMON /dtevt1/ nhkk,nevhkk,isthkk(nmxhkk),idhkk(nmxhkk),
     &                jmohkk(2,nmxhkk),jdahkk(2,nmxhkk),
     &                phkk(5,nmxhkk),vhkk(4,nmxhkk),whkk(4,nmxhkk)

* extended event history
      COMMON /dtevt2/ idres(nmxhkk),idxres(nmxhkk),nobam(nmxhkk),
     &                idbam(nmxhkk),idch(nmxhkk),npoint(10),
     &                ihist(2,nmxhkk)

      DOUBLE PRECISION slpx, slpy, slpz, sle, slm
      INTEGER slpid, slcharge
      COMMON /dpmjetparticle/ slpx(nmxhkk), slpy(nmxhkk), slpz(nmxhkk),
     &       sle(nmxhkk), slm(nmxhkk), slpid(nmxhkk), slcharge(nmxhkk)


C     >> Set Counter to Zero

      nfinal=0
      
      DO 42 i=1, nhkk
c      I = IPART

CC       >> Remove all non-final-state particles
        IF(.not.(isthkk(i).eq.1.or.isthkk(i).eq.-1.or.
     $isthkk(i).eq.1001)) goto 42

C >> Find Particle Charge, qch
        IF((abs(isthkk(i)).eq.1).and.(idhkk(i).ne.80000))THEN
C         >> final state ptcles except nuclei

          qch=ipho_chr3(idhkk(i),1)/3
        ELSEIF(idhkk(i).eq.80000)THEN
C         >> final state nuclei
          qch=idxres(i)
        ELSE
C         >> not a final state particle, qch not interesting
          qch=-999
        ENDIF

  nfinal = nfinal + 1
  slpx(nfinal) = phkk(1,i)
        slpy(nfinal) = phkk(2,i)
        slpz(nfinal) = phkk(3,i)
        sle(nfinal) = phkk(4,i)
        slm(nfinal) = phkk(5,i)
        slpid(nfinal) = idhkk(i)
        slcharge(nfinal) = qch

 42     CONTINUE
        nparticles = nfinal
  
      END

      SUBROUTINE dt_usrhis(MODE)
c Dummy to make the linker happy
      END