PROGRAM BETA_SPECT
c
c     Program to calculate the shape of a beta decay spectrum
c     and to select a set of initial energies from that 
c     distribution
c
      IMPLICIT NONE
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   I                ! loop index
      INTEGER   J                ! loop index
      REAL      E                ! electron total energy
      REAL      DNDEMAX          ! max value of dN/dE
      REAL      PROB_E           ! probability for given E
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=1)
      PARAMETER (NTRY=10000)
      PARAMETER (NTEST=1000)
C
C     Hbook stuff:
C
      REAL HMEMOR
      INTEGER NWPAW
      PARAMETER (NWPAW=100000)
      COMMON / PAWC / HMEMOR(NWPAW)
c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
      CALL HLIMIT ( NWPAW )                ! initialize hbook
      CALL HBOOK1 ( 2,'dN/dE chosen ',100,1.0,EPS0,0. )

c
c     initialize beta distribution:
c
      CALL BETA_INI ( NHST_FUN, EPS0, E0, MASS, DNDEMAX )
C
      WRITE ( 6,* ) ' DNDEMAX=',DNDEMAX
c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c     Now, choose a series of points from this distribution
c
      DO I=1,NTRY
        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
        CALL HFILL(2,E/MASS,0.,1.)
      END DO
c
c     save the histograms in a file:
c
      write ( 6,* ) ' save hists'
      CALL HRPUT ( 0,'beta.hst','N')
C
      STOP
      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
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
      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
      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