      SUBROUTINE GUKINE
      implicit none
c
c     Authors: S.R. Tonse and C.F. Maguire (FOPI relic)
c
C     Sets up subevent filling of primary particles into GEANT common blocks
C
c
c     Revision History:
c
C         Modification: 8/4/93 SRTonse added option KINE=14 to call a
C                       MU_SPGEN particle generator used by muon people.
C         JPSullivan Oct 4, 1993 --- redirected output to LOUT for non-error
C                       messages and to "*" for error messages.
C
C         C. F. Maguire: November 22, 1994, changes for GEN_EVT generic events
C
c         H.van Hecke Jan 96: allow 3d spread of vertex
c
c         C. F. Maguire: August 27, 1997, change FORMAT 800 I4 ---> I7
c
c         C. F. Maguire: May 30, 1998, add primary particle ZEBRA bank;
c                                      new call to mzwipe (was in gudigi)
c
c
c     Global variables
c
      include 'gcflag.inc'
      include 'gckine.inc'
      include 'guphnx.inc'
      include 'guevgen.inc'
      include 'sublink.inc'
      include 'fpxlink.inc'
      include 'fstore.inc'
      include 'subevt.inc'
      include 'udst.inc'
      include 'gcunit.inc'
      real ekin_init_jps
      common / jpstesting / ekin_init_jps
 
c
c
c     Local variables
c     ===============
 
      INTEGER   IDENT, IER, I
      Integer   EVT_STATUS, IDXPRT
      Integer   NTTRCK, NREMAIN, NSTORE, NVTX, NUBUF
 
c
c     KINE and VERT user space
c         local variables moved into GUKINE
c     ========================
      Integer NUKINE, NUVERT
      Parameter (NUKINE = 0)
      Parameter (NUVERT = 0)
      Real UVERT (10), UKINE (10)

c
c     FPRI ZEBRA bank additions
c
      character*20 chform            ! ZEBRA formatting
      integer ifirst /0/             ! initialization key
      integer ioixin /0/             ! MZFORM return value
      integer nmbook                 ! number of words to book
      integer jpoint                 ! first position
      integer ipoint                 ! counts number of records
      integer itrack_sub             ! track number in subevent
c
c     jps tests:
c
C
C     external functions:
C
      REAL      F_DNDE           ! function which returns shape
                                 ! of beta decay spectrum
      EXTERNAL  F_DNDE
C
      REAL      XRAND(2)         ! for random numbers
      INTEGER   J                ! loop index
      REAL      E                ! electron total energy
      REAL      DNDEMAX          ! max value of dN/dE
      SAVE      DNDEMAX
      REAL      PROB_E           ! probability for given E
      REAL      XYZ_VERT_JPS(3)  ! for vertex position
      REAL      PLAB_JPS(3)      ! for lab momentum into gskine
      REAL      PINIT_JPS        ! total momentum (MeV/c) initial
      REAL      COSTH_JPS        ! random costh
      REAL      SINTH_JPS        ! random sinth
      REAL      PHI_JPS          ! random phi
c
C     various parameters
c
      REAL      ENDPT            ! endpoint kinetic energy (MeV)
      REAL      E0               ! endpoint total energy (MeV)
      REAL      MASS             ! electron mass (MeV)
      REAL      EPS0             ! E0/MASS
      INTEGER   NHST_FUN         ! hist # to store beta dist. shape
      INTEGER   NTRY             ! # random entries in spectrum
      INTEGER   NTEST            ! max number of tests for a point
      PARAMETER (MASS=0.511)
      PARAMETER (ENDPT=3.54)
      PARAMETER (E0=ENDPT+MASS)
      PARAMETER (EPS0= E0/MASS)
      PARAMETER (NHST_FUN=7998)
      PARAMETER (NTRY=10000)
      PARAMETER (NTEST=1000)
      REAL      TWOPI            ! 2*pi
      PARAMETER (TWOPI=6.283185)



c
c     Executable code
c     ===============
 
      ikine = 99        ! my own code
c
      if(ifirst.eq.0)then
         ifirst = 1
c
c        FPRI data bank as 5 words per entry
c
         chform = '1I / 2I 3F'
         call mzform('FPRI', chform, ioixin)
c
c        jps tests
c
         CALL HBOOK1 ( 7999,'dN/dE chosen',100,1.0,EPS0,0. )
         CALL HBOOK1 ( 7997,'Ekin hit scint',100,0.,5.,0. )
         CALL HBOOK1 ( 7996,'Ekin_init hit scint',100,0.,5.,0. )
         CALL HBOOK1 ( 7995,'Ekin_init all',100,0.,5.,0. )

c        initialize beta distribution:
c
         CALL BETA_INI ( NHST_FUN, EPS0, E0, MASS, DNDEMAX )
C
         WRITE ( 6,* ) ' DNDEMAX=',DNDEMAX
c
      endif  ! initialization check

      if (ikine.eq.0) then          ! read data from file
   10 call gget (lun_dinp,cin_keys,nin_keys,ident,ier)
      if (ident.eq.0) then
         write(LOUT,*)'  <I> GUKINE: Read INIT structure'
         goto 10
         end if ! check on ident = 0
         if(ier.ne.0) then
            WRITE(*,*)' <E> GUKINE: Error code IER=',IER
            WRITE(*,*)' <E> GUKINE: terminate input at event #: ',
     +                IEVENT
            ieotri = 1  ! terminate this event
            ieorun = 1  ! terminate this run
         endif ! check on ier <> 0
      Else if (ikine .eq. 5) then
         Call LW_EVENT
C  option KINE 14 is an event generator similar to LW_LEAK, in that
C  each particle entered via GSKINE has a distinct vertex (origin). Code
C  is from the muon folks (S.Saini) who run this so as not to have to swim thru
C  magnet iron all the time. SRTonse 4th Aug 93
      Else if (ikine .eq. 14) then
         Call MU_SPGEN
      ELSE IF (IKINE .EQ. 99 ) THEN
c
c       jps test
c
         call mzwipe(ixdiv_fe)          ! clear event division
 101     CONTINUE
         DO J=1,NTEST
           CALL RANMAR(XRAND,2)
           E      = MASS + XRAND(1)*(E0-MASS)
           PROB_E = F_DNDE(EPS0,E)/DNDEMAX
           IF ( XRAND(2).LE.PROB_E ) GO TO 100
         END DO
 100     CONTINUE
         EKIN_INIT_JPS = E - MASS
         CALL HFILL(7999,E/MASS,0.,1.)
         CALL HFILL(7995,EKIN_INIT_JPS,0.,1.)
         IF ( EKIN_INIT_JPS.LE.0.5 ) GO TO 101  ! below threshold
         IPART = 2
         XYZ_VERT_JPS(1) = 0.
         XYZ_VERT_JPS(2) = 0.
         XYZ_VERT_JPS(3) = 0.
         CALL GSVERT ( XYZ_VERT_JPS, 0, 0, UVERT, 0, NVTX )
         CALL RANMAR(XRAND,2)
         PINIT_JPS   = SQRT(E**2 - MASS**2)/1000.  ! MeV to GeV
         CALL RANMAR(XRAND,2)
         COSTH_JPS   = -1. + 2.0*XRAND(1)
         COSTH_JPS   = MIN ( COSTH_JPS, 1.0 )
         COSTH_JPS   = MAX ( COSTH_JPS,-1.0 )
         SINTH_JPS   = SQRT(1.0 - COSTH_JPS**2)
         PHI_JPS     = TWOPI*XRAND(2)
         PLAB_JPS(1) = PINIT_JPS*SINTH_JPS*COS(PHI_JPS)
         PLAB_JPS(2) = PINIT_JPS*SINTH_JPS*SIN(PHI_JPS)
         PLAB_JPS(3) = PINIT_JPS*COSTH_JPS
c        PLAB_JPS(1) = PINIT_JPS
c        PLAB_JPS(2) = 0.
c        PLAB_JPS(3) = 0.
         CALL GSKINE ( PLAB_JPS, IPART, NVTX, UKINE,0,NTTRCK) 
      else
c
c       Handle Pisa standard event input. GUEVGEN reads complete events.
c         Breaking of events into sub-events is handled uniformly here for
c         all generators. GUEVGEN also produces vertex positions for the
c         complete event which must be used to create the vertex for
c         each sub-event.
c       ==================================================================
 
         If (END_EVTFLG) Then
c
c         If last event complete, call event input routine to fetch
c           complete event from user-chosen source of events.
c         =========================================================
            CALL GUEVGEN (EVT_STATUS)
c
c      EVT_STATUS has to be set in by the user event routines
c
            If (EVT_STATUS .NE. -1) Then
               If (EVT_STATUS .EQ. 2) Then
                  Write (*, 900)
                  IEOTRI = 1 ! terminate this event
               Else
                  Write (*, 901) EVT_STATUS
                  IEOTRI = 1 ! terminate this event
                  IEORUN = 1 ! terminate this run
               End if
               Go to 9999
            End if
c
c     Check that the predicted number of subevents will be less than the
c     limit of MAXSUB which can be handled by PISORP
c
            IF(IPOPSUB.GT.0)THEN
               IF((1 + MXTOT/IPOPSUB) .GT. MAXSUB)THEN
                  Write(*,902)(1 + MXTOT/IPOPSUB), MAXSUB
                  Write(*,903)IPOPSUB, 1 + MXTOT/MAXSUB
                  IEOTRI = 1 ! terminate this event
                  Go to 9999
               ENDIF  ! Safety check on predicted number of subevents
           ENDIF ! check on IPOPSUB > 0
 
c
c         Start new event, FZOUT will write "start of new evt" marker
c         ============================================================
            END_EVTFLG = .FALSE.
            BUDST_NEW_EVT = .TRUE.
            IEVPRT = 1
            NSUB_EVT = 0
            NTRU_EVT = NUMEVT
            CALL EVENT_HIST   ! Event Histogramming Array (J. Sullivan request)
         End if
 
c
c       Determine how many particles for this sub-event
c       ===============================================
         NSUB_EVT = NSUB_EVT + 1
         NREMAIN = MXTOT - IEVPRT + 1
         If (NREMAIN .LE. IPOPSUB) Then
            NSTORE = NREMAIN
            END_EVTFLG = .TRUE.
         Else
            NSTORE = IPOPSUB
         End if
c
c     Book the FPRI bank
c
         call mzwipe(ixdiv_fe)          ! clear event division
         nmbook = (nstore + 1)*mfx_kine
         lfx_kine(1) = 0                    !safety set
         call MZBOOK(ixdiv_fe,
     $               LFX_KINE(1),           !return bank address
     $               LFX_KINE(1),           !supporting link
     $               1,                     !JBIAS=1 ---> top level bank
     $               'FPRI',                !bank name
     $               0,                     !# of links
     $               0,                     !# of down links
     $               NMBOOK,                !# of data words
     $               IOIXIN,                !I/O characteristics
     $               -1)                    !do not clear memory contents
         IPOINT = LFX_KINE(1) + 1           ! counts the number of records
         JPOINT = IPOINT                    ! first position
         iqf(jpoint) = nstore

         if(nsub_evt.eq.1 .and. (kvtx_par(1,2).ne.0.0. or.
     &      kvtx_par(2,2).ne.0.0 .or. kvtx_par(3,2).ne.0.0) ) then
            write(LOUT,800)ntru_evt, (kvtx_par(i,2),i=1,3), xyz(3)
 800        format(/,' For event # ',i7,'  using X,Y,Z-RMS = ',3f6.2,
     +               '  obtain Z0 = ',f6.2,/)
         endif ! check on Z0 printout

         if(numevt.lt.3.or.(50*(numevt/50).eq.0))then
            write(LOUT, '(/,a,a, 4(/,1x,(a,i6)) )' ) '  GUKINE <I>:  '
     +      //'Event name   ',chevt_name, '  Total particles ',mxtot,
     +      '  Subevent number ',nsub_evt, '  Number stored this '
     +      //'subevent ',nstore, '  Number remaining ', nremain -
     +      nstore
         endif ! check on subevent printout
 
c
c     Randomize in the XY direction if requested
c
         IF(BUDST_NEW_EVT.AND.IXY_RNDM.GT.0)THEN
            CALL XY_RNDM
            IF(IXY_CHANGE.EQ.0)THEN
               IXY_CHANGE = 1
               WRITE(*, '(/,A,I2)' ) '   XY Randomization with IXYRND ='
     +         //' ', IXY_RNDM
            ENDIF ! check on notification
         ENDIF
 
c
c       Load VERT and KINE for next sub-event
c       =====================================
         IF(MULTI_XYZ.NE.1)THEN
            Call GSVERT (XYZ, 0, 0, UVERT, NUVERT, NVTX)
         ENDIF  ! GET A SINGLE VERTEX
 
c
c       Store tracks for this sub-event in KINE.
c       ========================================
         itrack_sub = 0
         Do IDXPRT = IEVPRT, IEVPRT + NSTORE - 1
            IF(MULTI_XYZ.EQ.1)THEN
               CALL GSVERT(XYZMV(1,IDXPRT), 0, 0, UVERT, NUVERT,
     +                     NVTX)
            ENDIF  ! CHECK ON DIFFERENT VERTEX FOR EACH PARTICLE
            IF(ID_PARENT(IDXPRT).GT.0)THEN
               NUBUF = 3     ! Will indicate an orphan track
               UKINE(1) = 0  ! Event generator does not furnish a parent track
               UKINE(2) = ID_PARENT(IDXPRT)  ! Parent ID from event generator
            ELSE
               NUBUF = NUKINE
            ENDIF  ! CHECK ON SENSIBLE PARENT ID NUMBER
            Call GSKINE (PPTOT (2, IDXPRT), IDTOT (IDXPRT), NVTX,
     +      UKINE, NUBUF, NTTRCK)
            iqf(ipoint + ofsx_tr) = idxprt
            iqf(ipoint + ofsx_id) = idtot(idxprt)
             qf(ipoint + ofsx_px) = pptot(2,idxprt)
             qf(ipoint + ofsx_py) = pptot(3,idxprt)
             qf(ipoint + ofsx_pz) = pptot(4,idxprt)
            itrack_sub = itrack_sub + 1
            iqf(ipoint + ofsx_trsub) = itrack_sub
            ipoint = ipoint + mfx_kine
         End do
 
         IEVPRT = IEVPRT + NSTORE
         MAXTRK = NTTRCK ! Last track stored
      End if
 
c
c     If user requests it, print out KINE bank
c     ========================================
      IF (IDEBUG .NE. 0) Then
         CALL GPKINE(0)
      End if
 
 9999 Continue
 
  900 Format (1H , 'GUKINE - Error ', I3, ' in GUEVGEN, aborting event')
  901 Format (1H , 'GUKINE - Error ', I3, ' in GUEVGEN, aborting run')
  902 Format (1H , 'GUKINE - Error:  Predicted number of subevents = ',
     1        i3, '  exceeds current limit in PISORP = ',i3)
  903 Format (1H , 'GUKINE <E>: You need to increase the ',
     1        'number of tracks per subevent from ', i4,
     2        ' to ', i4)
      Return
      End
      REAL FUNCTION F_DNDE ( EPS0, E )
c
c     Subroutine to calculate dN/dE (shape of beta decay
c     spectrum for electrons or positions). Most of the
c     details of the assumptions about the decay system
c     are contained in parameters inside this routine.
c     The gamow approximation is used for the coulomb
c     correction.
c     J.P.Sullivan 5-Jan-1999
c
      IMPLICIT NONE
C
C     input:
      REAL      EPS0             ! E0/MASS
      REAL      E                ! electron total energy (MeV)
C
C     the shape of the spectrum is controlled by these parameters:
C
      REAL      CHARGE           ! beta charge (-1 or +1)
      REAL      Z                ! charge of residual nucleus
      PARAMETER (CHARGE=-1.)
      PARAMETER (Z     =46.)

c     misc constants:
C
      REAL      MASS             ! electron mass (MeV)
      REAL      TWOPI            ! 2*pi
      REAL      ALPHA            ! fine structure constant
      PARAMETER (MASS  = 0.511  )
      PARAMETER (TWOPI=6.283185)
      PARAMETER (ALPHA = 1./137.036)
C
c     kinematic variables of the electron:
c
      REAL      BETA             ! electron velocity
      REAL      EPS              ! E/MASS
      REAL      P                ! electron momentum
c
c     variables related to the coulomb correction:
c
      REAL      ETA              ! Z*ALPHA/BETA
      REAL      ETA2PI           ! 2*PI*ETA
      REAL      GAMOW            ! coulomb correction (point approx)
C
Ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
C
c     write ( 6,* ) ' e=',e
      EPS     = E/MASS
      P       = SQRT(E**2 - MASS**2)
      BETA    = P/E
      ETA     = Z*ALPHA/BETA
      ETA2PI  = -CHARGE*TWOPI*ETA
      IF ( Z.NE.0 ) THEN
        GAMOW   = ETA2PI/(1.0-EXP(-ETA2PI))
      ELSE
        GAMOW   = 1.0
      END IF
      F_DNDE = ((EPS0-EPS)**2)*EPS*SQRT(EPS**2-1.0)*GAMOW

      RETURN
      END
      SUBROUTINE BETA_INI ( NHST, EPS0, E0, MASS, DNDEMAX )
c
c     Initialization routine -- calculate histogram with shape
c     if beta spectrum (in hist NHST) and find the max value
c     in the histogram.
c     J.P.Sullivan 5-Jan-1999
c
c
      IMPLICIT NONE
C
C     Input:
      INTEGER   NHST             ! histogram to fill
      REAL      EPS0             ! E0/MASS
      REAL      E0               ! endpoint total energy (MeV)
      REAL      MASS             ! electron mass (MeV)
C     Output:
      REAL      DNDEMAX          ! max value of dN/dE
C
C     external references:
C
      REAL      F_DNDE
      EXTERNAL  F_DNDE
C
C     local variables:
C
      INTEGER   I
      REAL      E                ! electron total energy
      REAL      STEP             ! binsize in dN/de histogram
      REAL      OFFSET           ! offset for calculating bins
C
C
      INTEGER   NSTEP            ! # bins in histogram of dN/de
                                 ! also sets granularity is search
                                 ! for max of dN/de function
      PARAMETER (NSTEP=1000)     ! bins in histogram
      REAL      DNDE(NSTEP)      ! dN/dE histogram value
      REAL      EPS(NSTEP)       ! eps value in bin
C
c     book the histogram which will be filled
c
      CALL HBOOK1 ( NHST,'F dN/dE',NSTEP,1.0,EPS0,0. )
c
      DNDEMAX = -1.
c
      STEP   = (E0-MASS)/FLOAT(NSTEP)
      OFFSET = -STEP/2. + MASS
c
      DO I=1,NSTEP
c
        E       = FLOAT(I)*STEP + OFFSET
        EPS(I)  = E/MASS
        DNDE(I) = F_DNDE ( EPS0, E )
        DNDEMAX = MAX(DNDE(I),DNDEMAX)
c
c       WRITE ( 6,* ) ' EPS(I)=',EPS(I),' DNDE(I)=',DNDE(I),
c    1                ' DNDEMAX=',DNDEMAX
      END DO
c
c     normalize the spectrum to 1 and put it in a histogram
c
      DO I=1,NSTEP
c
        DNDE(I) = DNDE(I)/DNDEMAX
        CALL HFILL(NHST,EPS(I),0.,DNDE(I))
c
      END DO

      RETURN
      END
