SUBROUTINE SUSBCH( when, iway, ierr )
C ----------------------------------------------------------------------------
C  Purpose: Calculate photino decay branching fractions and decide which way 
C           the photino has decayed for a given event.
C
C    Notes: If when = 1, the partial widths and branching ratios are 
C           reintegrated and iway new sets of energies are generated 
C           for each of the included decay paths.
C
C           If when = 2, a decision is made as to which of the decays actually
C           happened in this event.
C
C   Author: JMB 
C   Last change 18/2/92
C ----------------------------------------------------------------------------
C
      IMPLICIT NONE
C
      INCLUDE 'SUSANA.INC'
C
      INTEGER when, iway, ierr, i, j, temp, randum
      INTEGER maxnum
      REAL rno, prob, e12(2), ulmass, gnerr, geerr
      REAL ma, mb, mc, ecms, highe(2), lowe(2), rn32
C
      IF (when.EQ.1) THEN
C
        regen = .TRUE.
C --    Calculate included photino branching fractions.
C
        temp = isustr(6)
        DO 1, i = 1, 4
          dk(i) = MOD( temp, 10)
          IF ((dk(i).NE.1).AND.(dk(i).NE.0)) THEN
            WRITE(isustr(1),*)
     &     'SUSINT:Illegal photino branching selection!',isustr(6)
            ierr = 1
            GOTO 8900
          ENDIF
          temp = temp/10
 1      CONTINUE
C
        IF (isustr(3).GT.0) THEN
          WRITE(isustr(1),9001) iway, event
        ENDIF
C
C --    Integrate gamma_nu   --  --  --  --  
C
        IF ( msus(8)**2 .LT.(mdd2+msin2 + 2.0*sqrt(mdd2*msin2)) ) THEN
          gamma_nu = 0.0
        ELSE
C
C -       Set the kinematic bounds
          ma = 0
          mb = SQRT(msin2)
          mc = mdd2
          ecms = msus(8)
          lowe(1)  = ma
          lowe(2)  = mb
C
          highe(1) = max( ((ecms-mb)**2-mc**2+ma**2)/(2*(ecms-mb)),
     &                  ((ecms-mc)**2-mb**2+ma**2)/(2*(ecms-mc)))
          highe(2) = max( ((ecms-ma)**2-mc**2+mb**2)/(2*(ecms-ma)),
     &                  ((ecms-mc)**2-ma**2+mb**2)/(2*(ecms-mc)))
C
          function_type = 101
          maxnum = 100000
          CALL PARTN( 2, lowe, highe, 1.0, maxnum)
          CALL INTGRL( 2, 0, 200, gamma_nu, gnerr )
          IF (dk(3)+dk(4).GT.0) THEN
C           Generate a new set of neutrino decay energies
            DO 2, i = 1, iway
              CALL RANGEN( 2, E12 )
              e_nu(2*i-1) = E12(1)
              e_nu(2*i)   = E12(2)
 2          CONTINUE       
          ENDIF
          function_type = 0
        ENDIF
C
C --    Integrate gamma_e    --  --  --  
C
        IF ( msus(8)**2 .LT.(msin2+mdu2 + 2.0*sqrt(msin2*mdu2)) ) THEN
          gamma_e = 0.0
        ELSE
C -       Set the kinematic bounds
          ma = ulmass( 11 )
          mb = SQRT(msin2)
          mc = SQRT(mdu2)
          ecms = msus(8)
          lowe(1)  = ma
          lowe(2)  = mb
C
          highe(1) = max( ((ecms-mb)**2-mb**2+ma**2)/(2*(ecms-mb)),
     &                  ((ecms-mc)**2-mc**2+ma**2)/(2*(ecms-mc)))
          highe(2) = max( ((ecms-ma)**2-ma**2+mb**2)/(2*(ecms-ma)),
     &                  ((ecms-mc)**2-mc**2+mb**2)/(2*(ecms-mc)))
C
          function_type = 102
          maxnum = 100000
          CALL PARTN( 2, lowe, highe, 1.0, maxnum)
          CALL INTGRL( 2, 0, 200, gamma_e, geerr )
          IF (dk(1)+dk(2).GT.0) THEN
C           Generate a new set of electron decay energies
            DO 3, i = 1, iway
              CALL RANGEN( 2, E12 )
              e_e(2*i-1) = E12(1)
              e_e(2*i)   = E12(2)
 3          CONTINUE       
          ENDIF
          function_type = 0
C
        ENDIF
C -     Calculate total width (remember charge conjugates)
        gamma_tot = 2.*gamma_nu + 2.*gamma_e
C
C -     Calculate branching fractions
        br_nu = gamma_nu/gamma_tot
        br_e  = gamma_e/gamma_tot
C
C -     Calculate total included branching fraction for this run.
C       ( This multiplies the final cross section )
        branch = 0.0
        DO 4 i = 1, 2
          branch = branch + dk(i)*br_e
 4      CONTINUE
        DO 5 i = 3, 4
          branch = branch + dk(i)*br_nu
 5      CONTINUE
C
        IF (branch. GT. 1) THEN
          ierr = 1
          WRITE (isustr(1),*) 'SUSBCH: Branching fraction > 1!'
          GOTO 8900
        ENDIF
C
C -     Calculate photino lifetime (seconds).
        lphino = 6.582E-25/(gamma_tot*branch)
C
C       Zero counters
        DO 6 i = 1, 2
          e_used(i) = 0
 6      CONTINUE
C
      ELSEIF (when.EQ.2) THEN
C
C --    Decide which decay the  photino made.
C        
        rno = RN32(randum)
        prob = 0.0
        DO 7 i = 1, 4
C
          IF (i.LE.2) THEN
            prob = prob + dk(i)*br_e/branch
          ELSE
            prob = prob + dk(i)*br_nu/branch
          ENDIF
C
          IF (rno.LT.prob) THEN
            iway = i
            GOTO 7900
          ENDIF
C
 7      CONTINUE   
C
 7900   CONTINUE
C
      ENDIF
C
 8900 CONTINUE
      RETURN
C
 9001 FORMAT (4X,'SUSBCH:Generating ',I3,' more energies at event'
     &,           I6)
C
      END