program GENERATE_DECAYS

c  Glen Cowan
c  12 December, 1999

c  Generates momentum vectors for 3 particles.
c  Uses the subroutine RAMBO by R. Kleiss, W.J. Stirling
c  and S.D. Ellis, Comp. Phys. Commun. 40 (1986) 359

       implicit        NONE

c  Constants

       integer           num_particles
       parameter         (num_particles = 3)

       real              m_pi
       parameter         (m_pi = 0.139)      ! GeV

       real              m2_ave_0_minus
       parameter         (m2_ave_0_minus = 3.6380252E-06)
       real              m2_ave_1_minus
       parameter         (m2_ave_1_minus = 2.7610781E-02)

       real              m2_max_0_minus
       parameter         (m2_max_0_minus = 2.9060666E-05)
       real              m2_max_1_minus
       parameter         (m2_max_1_minus = 5.6871582E-02)

       double precision  E_cm
       parameter         (E_cm = 1.0)        ! GeV
       double precision  mass (num_particles)  / num_particles*m_pi /

c  Needed for HBOOK routines

       integer           hsize
       parameter         (hsize = 1000000)
       integer           hmemor (hsize)
       common  /pawc/    hmemor

c  Common block for ntuple

       real              p1, p2, p3
       common  /evtblk/  p1(3), p2(3), p3(3)

c  Local variables

       character*2       jp                / ' ' /
       character*80      outfile

       integer           i
       integer           icycle
       integer           id
       integer           j
       integer           k
       integer           istat
       integer           lw
       integer           num_0_minus
       integer           num_events
       integer           option

       logical           accept_event

       real              r
       real              rvec(1)
       real              theta_0            !  BR for 0- decay
       real              m2_ave
       real              m2_max
       real              M2_0_MINUS
       real              M2_1_MINUS
       real              weight

       double precision  p (4, num_particles)
       double precision  wt

c  Initialize HBOOK

       call HLIMIT (hsize)
       write (*, *) 'output file name'
       read (*, fmt='(a80)') outfile
       call HROPEN (20, 'ntuple', outfile, 'N', 1024, istat)

c  Set up decay properties

       do i = 1, num_particles
         mass(i) = m_pi
       end do
       lw = 1             ! requests unweighted events

       write (*, *) 'Integrate phase space, find m2_max       1'
       write (*, *) 'Generate events of specific Jp           2'
       write (*, *) 'Generate events with mixture of Jp       3'
       write (*, *) 'Enter option:'
       read (*, *) option

       if ( option .eq. 1 ) then

         write (*, *) 'enter spin, parity of particle, e.g. 0-, 1-, 1+'
         read (*, fmt='(a2)') jp
         call HBOOK1 (100, 'w', 100, 0., 0.1, 0.)
         call HIDOPT (100, 'STAT')
         call HBOOK1 (101, '-log10 w', 100, 0., 50., 0.)
         call HIDOPT (101, 'STAT')

       elseif ( option .eq. 2 ) then

         write (*, *) 'enter spin, parity of particle, e.g. 0-, 1-, 1+'
         read (*, fmt='(a2)') jp
         id = 1
         call HBNT (id, 'p1(3), p2(3), p3(3)', ' ')
         call HBNAME (id, 'evtblk', p1, 'p1(3), p2(3), p3(3)')

         if ( jp .eq. '0-' ) then
           m2_max = m2_max_0_minus
         elseif ( jp .eq. '1-' ) then
           m2_max = m2_max_1_minus
         else
           write (*, *) 'spin-parity combination not recognized'
         endif

       elseif ( option .eq. 3 ) then

         write (*, *) 'Enter BR for 0- decay'
         read (*, *) theta_0
         id = 1
         call HBNT (id, 'p1(3), p2(3), p3(3)', ' ')
         call HBNAME (id, 'evtblk', p1, 'p1(3), p2(3), p3(3)')

       endif

c  Either integrate phase space or generate events

       if ( option .eq. 1 ) then           ! integrate phase space

         write (*, *) 'enter number of events'
         read (*, *) num_events

         m2_ave = 0.
         do i = 1, num_events

           call RAMBO (num_particles, E_cm, mass, p, wt, lw )
           do j = 1, 3
             p1(j) = p(j,1)
             p2(j) = p(j,2)
             p3(j) = p(j,3)
           end do

           if ( jp .eq. '0-' ) then
             weight = M2_0_MINUS (p1, p2, p3)
           elseif ( jp .eq. '1-' ) then
             weight = M2_1_MINUS (p1, p2, p3)
           endif
           m2_ave = m2_ave + weight
           call HF1 (100, weight, 1.)
           call HF1 (101, -LOG10 (weight), 1.)
           m2_max = MAX (weight, m2_max)

         end do
         m2_ave = m2_ave / FLOAT (num_events)

         write (*, *) 'm2_max = ', m2_max
         write (*, *) 'm2_ave = ', m2_ave

       elseif ( option .eq. 2 .or. option .eq. 3 ) then 

         num_0_minus = 0
         write (*, *) 'enter number of events to generate'
         read (*, *) num_events
         do i = 1, num_events

c  if option 3, flip coin to decide whether decay is 0- or 1-

           if ( option .eq. 3 ) then
             call RANMAR (rvec, 1)
             if ( rvec(1) .le. theta_0 ) then
               jp = '0-'
               m2_max = m2_max_0_minus
               num_0_minus = num_0_minus + 1
             else
               jp = '1-'
               m2_max = m2_max_1_minus
             endif
           endif

c  generate events according to phase space, accept with probability
c  proportional to the weight.

           accept_event = .false.
           do while ( .not. accept_event )

             call RAMBO (num_particles, E_cm, mass, p, wt, lw )
             do j = 1, 3
               p1(j) = p(j,1)
               p2(j) = p(j,2)
               p3(j) = p(j,3)
             end do

             if ( jp .eq. '0-' ) then
               weight = M2_0_MINUS (p1, p2, p3)
             elseif ( jp .eq. '1-' ) then
               weight = M2_1_MINUS (p1, p2, p3)
             endif

             call RANMAR (rvec, 1)
             r = rvec(1) * m2_max
             accept_event = r .lt. weight

           end do                         ! .not. accept_event
           call HFNT (id)

         end do                           ! i = 1, num_events

       endif                              ! option

       if ( option .eq. 3 ) then
         write (*, *) 'number of 0- events generated = ', num_0_minus
       endif

c  create subroutine needed later to read in ntuple (only needed once)

cc       open (30, file = 'read_cwn_setup.f', status = 'unknown')
cc       call HUWFUN (30, id, 'READ_CWN_SETUP', 0, 'B')
cc       close (30)

c  Store histogram and close.
     
       call HROUT (0, icycle, ' ')
       call HREND ('ntuple')

       stop
       END