      SUBROUTINE SUSINT( ierr )
C -------------------------------------------------------------------------
C  Purpose: Initialise common block for R-parity-violating SUSY
C           cross-section calculation.
C
C   Errors: ierr = 0        : All OK
C           ierr =  1 -> 10 : Error in SUSINT (fatal)
C           ierr = 11 -> 100: Error in SUSBCH (fatal)
C           ierr      >  100: Error in SUSARR (fatal)
C
C   Author: JMB 
C   Last change 18/2/92
C -------------------------------------------------------------------------
C
      IMPLICIT NONE
C
      INCLUDE 'SUSANA.INC'
      INCLUDE 'PYTHIA57.INC'
C
      EXTERNAL SUSDAT
      EXTERNAL LUDATA
      EXTERNAL PYDATA
C
      REAL time, ulmass, term, eq1, eq2
      INTEGER ierr, i, iway
      LOGICAL ok
      CHARACTER*20 pdfpar( 20 )
      DOUBLE PRECISION pdfval(20)
C
C -   Start timer.
      CALL TIMED (time)
C
C -   Some zeroing.
      Event = 0
      ierr = 0
      total_xsec = 0.0
      error_xsec = 0.0
      br_e  = 0.0
      br_nu = 0.0
C
C -   Print banner.
      WRITE( isustr(1), 1001 )
C
C -   Initialize random number generator.
      IF (suseed.NE.0) THEN
        CALL RN32IN(suseed)
      ENDIF
C
      pmas(6,1) = mtop
      mp2 = msus(8)**2
C
C --  Store beam momenta and CMS energy squared (GeV, GeV^2).
C     (massless proton, electron assumed)
      pl(4) = SQRT( pl(1)*pl(1) + pl(2)*pl(2) + pl(3)*pl(3) )
      pp(4) = SQRT( pp(1)*pp(1) + pp(2)*pp(2) + pp(3)*pp(3) )
      epcms = 2.0*(pl(4)*pp(4) - 
     &                pl(1)*pp(1) - pl(2)*pp(2) - pl(3)*pp(3))
C
C -   Print beam conditions.
C
C     If Electron is incident.
C
      If(isustr(8).eq.0) then
      WRITE( isustr(1),1002) pp(1), pp(2), pp(3), pl(1), pl(2), pl(3), 
     &      SQRT( Epcms)
      End If
C
C    If Positron is incident.
C
      If(isustr(8).eq.1) then
      WRITE( isustr(1),1022) pp(1), pp(2), pp(3), pl(1), pl(2), pl(3), 
     &      SQRT( Epcms)
      End If

C     
C -   Print kinematic limits.
      WRITE( isustr(1),1004) low(1), high(1), -high(2)/(epcms*high(1))
     &,     (epcms*high(1)-msus(8))/(epcms*high(1)),-high(2),-low(2)
C
C   If Electron incident.
C
      If(ISUSTR(8).eq.0) then
      eq1 = ed
      eq2 = eu
C
C --Supersymmetric  Singlet part - Specify outgoing quark in ss decay
C
      IF ( MOD(ISUSTR(5),3).EQ.1 )  THEN
        msin2 = ulmass(1)**2
        ssq = 1
      ELSE IF ( MOD(isustr(5),3).EQ.2 ) THEN
        msin2 = ulmass(3)**2
        ssq = 3
      ELSE IF ( MOD(ISUSTR(5),3).EQ.0 ) THEN
        msin2 = ulmass(5)**2
        ssq = 5
      ENDIF
C
C --  Doublet part.
      IF ( ISUSTR(5).LE.3 )  THEN
        mdd2 = ulmass(1)**2
        mdu2 = ulmass(2)**2
        dsq1 = 1
        dsq2 = 2
      ELSE IF ( isustr(5).LE.6 ) THEN
        mdd2 = ulmass(3)**2
        mdu2 = ulmass(4)**2
        dsq1 = 3
        dsq2 = 4
      ELSE IF ( isustr(5).LE.9 ) THEN
        mdd2 = ulmass(5)**2
        mdu2 = mtop**2
        dsq1 = 5
        dsq2 = 6
      ENDIF
      ENDIF
C
C   If Positron incident.
C
      If(ISUSTR(8).eq.1) then
      eq1 = eu
      eq2 = ed
C
C --Supersymmetric  Singlet part - Specify outgoing quark in ss decay
C
      IF ( MOD(ISUSTR(5),3).EQ.1 )  THEN
        msin2 = ulmass(2)**2
        ssq = 2
      ELSE IF ( MOD(isustr(5),3).EQ.2 ) THEN
        msin2 = ulmass(4)**2
        ssq = 4
      ELSE IF ( MOD(ISUSTR(5),3).EQ.0 ) THEN
        msin2 = ulmass(6)**2
        ssq = 6
      ENDIF
C
C --  Doublet part.
      IF ( ISUSTR(5).LE.3 )  THEN
        mdd2 = ulmass(2)**2
        mdu2 = ulmass(1)**2
        dsq1 = 2
        dsq2 = 1
      ELSE IF ( isustr(5).LE.6 ) THEN
        mdd2 = ulmass(4)**2
        mdu2 = ulmass(3)**2
        dsq1 = 4
        dsq2 = 3
      ELSE IF ( isustr(5).LE.9 ) THEN
        mdd2 = mtop**2
        mdu2 = ulmass(6)**2
        dsq1 = 6
        dsq2 = 5
      ENDIF
      ENDIF
C
      CALL SUSCHK(ok)
      IF (.NOT.ok) THEN
        ierr = 1
        GOTO 8900
      ENDIF
C
C -   Calculate d-like squark (electron) or u-like (positron) width (GeV).
      IF (msus(ssq)**2.LE.(mp2+msin2 + 2.0*sqrt(mp2*msin2))) THEN
        term = 0.0
      ELSE
        term = e2*eq1**2*msus(ssq)*(1.-(mp2+msin2)/msus(ssq)**2) 
     &     * SQRT( (1.-(mp2+msin2)/msus(ssq)**2) 
     &     *       (1.-(mp2-msin2)/msus(ssq)**2) )/(8.*pi)
      ENDIF
      wssq = term 
     &     + lambda**2 * msus(ssq) / (16.*pi) 
     &     *( SQRT(1.-mdu2/msus(ssq)**2) + SQRT(1.-mdd2/msus(ssq)**2)) 
C
C -   Calculate u-like squark (electron) or d-like (positron) width (GeV).
      IF (msus(dsq2)**2.LE.(mp2+mdu2 + 2.0*sqrt(mp2*mdu2))) THEN
        term = 0.0
      ELSE
        term = e2*eq2**2*msus(dsq2)*(1.-(mp2+mdu2)/msus(dsq2)**2) 
     &     * SQRT( (1.-(mp2+mdu2)/msus(dsq2)**2) 
     &     *       (1.-(mp2-mdu2)/msus(dsq2)**2) )/(8.*pi)
      ENDIF
      wdsq2  = term 
     &     + lambda**2 * msus(dsq2) / (16.*pi) 
     &     *SQRT(1.-msin2/msus(dsq2)**2) 
C
C -   Calculate inlcuded photino branching fractions and store energies.
      IF (nevent.GT.100) THEN
        iway = 100
      ELSE
        iway = nevent
      ENDIF
      CALL SUSBCH( 1, iway, ierr )
      IF (ierr.NE.0) THEN
        ierr = ierr + 10
        GOTO 8900
      ENDIF
C
C --  Check kinematics and print out initial calculations.
      WRITE (isustr(1),1111) (msus(i), i=1, 8), lphino, mtop, lambda
     &, isustr(6), branch
C
C   Do Electron Write Outs Here
C
      IF (     ISUSTR(8).EQ. 0) THEN
C
      IF ( MOD(ISUSTR(5),3).EQ.1 )  THEN
C
        IF (ISUSTR(5).EQ.1) THEN
          WRITE(ISUSTR(1),1009)
          WRITE(isustr(1), 1005) wssq 
          WRITE(isustr(1), 1025) wdsq2
        ELSE IF (ISUSTR(5).EQ.4) THEN
          WRITE(ISUSTR(1),1012)
          WRITE(isustr(1), 1005) wssq 
          WRITE(isustr(1), 1026) wdsq2
        ELSE
          WRITE(ISUSTR(1),1015)
          WRITE(isustr(1), 1005) wssq 
          WRITE(isustr(1), 1027) wdsq2
        ENDIF
C
      ELSE IF ( MOD(isustr(5),3).EQ.2 ) THEN
        WRITE(isustr(1), 1006) wssq, wdsq2
        IF (ISUSTR(5).EQ.2) THEN
          WRITE(ISUSTR(1),1010)
          WRITE(isustr(1), 1006) wssq 
          WRITE(isustr(1), 1025) wdsq2
        ELSE IF (ISUSTR(5).EQ.5) THEN
          WRITE(ISUSTR(1),1013)
          WRITE(isustr(1), 1006) wssq 
          WRITE(isustr(1), 1026) wdsq2
        ELSE
          WRITE(isustr(1),1016)
          WRITE(isustr(1), 1006) wssq 
          WRITE(isustr(1), 1027) wdsq2
        ENDIF
C
      ELSE IF ( MOD(ISUSTR(5),3).EQ.0 ) THEN
        WRITE(isustr(1),1007) wssq, wdsq2
        IF (ISUSTR(5).EQ.3) THEN
          WRITE(ISUSTR(1),1011)
          WRITE(isustr(1), 1007) wssq 
          WRITE(isustr(1), 1025) wdsq2
        ELSE IF (ISUSTR(5).EQ.6) THEN
          WRITE(ISUSTR(1),1014)
          WRITE(isustr(1), 1007) wssq 
          WRITE(isustr(1), 1026) wdsq2
        ELSE
          WRITE(ISUSTR(1),1017)
          WRITE(isustr(1), 1007) wssq 
          WRITE(isustr(1), 1027) wdsq2
        ENDIF
      ENDIF
      ENDIF
C
C   Do Positron Write Outs Here
C
      IF (     ISUSTR(8).EQ. 1) THEN
C
      IF ( MOD(ISUSTR(5),3).EQ.1 )  THEN
C
        IF (ISUSTR(5).EQ.1) THEN
          WRITE(ISUSTR(1),1209)
          WRITE(isustr(1), 1205) wssq 
          WRITE(isustr(1), 1225) wdsq2
        ELSE IF (ISUSTR(5).EQ.4) THEN
          WRITE(ISUSTR(1),1212)
          WRITE(isustr(1), 1205) wssq 
          WRITE(isustr(1), 1226) wdsq2
        ELSE
          WRITE(ISUSTR(1),1215)
          WRITE(isustr(1), 1205) wssq 
          WRITE(isustr(1), 1227) wdsq2
        ENDIF
C
      ELSE IF ( MOD(isustr(5),3).EQ.2 ) THEN
        WRITE(isustr(1), 1206) wssq, wdsq2
        IF (ISUSTR(5).EQ.2) THEN
          WRITE(ISUSTR(1),1210)
          WRITE(isustr(1), 1206) wssq 
          WRITE(isustr(1), 1225) wdsq2
        ELSE IF (ISUSTR(5).EQ.5) THEN
          WRITE(ISUSTR(1),1213)
          WRITE(isustr(1), 1206) wssq 
          WRITE(isustr(1), 1226) wdsq2
        ELSE
          WRITE(isustr(1),1216)
          WRITE(isustr(1), 1206) wssq 
          WRITE(isustr(1), 1227) wdsq2
        ENDIF
C
      ELSE IF ( MOD(ISUSTR(5),3).EQ.0 ) THEN
        WRITE(isustr(1),1007) wssq, wdsq2
        IF (ISUSTR(5).EQ.3) THEN
          WRITE(ISUSTR(1),1211)
          WRITE(isustr(1), 1207) wssq 
          WRITE(isustr(1), 1225) wdsq2
        ELSE IF (ISUSTR(5).EQ.6) THEN
          WRITE(ISUSTR(1),1214)
          WRITE(isustr(1), 1207) wssq 
          WRITE(isustr(1), 1226) wdsq2
        ELSE
          WRITE(ISUSTR(1),1217)
          WRITE(isustr(1), 1207) wssq 
          WRITE(isustr(1), 1227) wdsq2
        ENDIF
      ENDIF
      ENDIF
C
C --  Set up chosen Structure function.
      pdfpar(1) = 'NGROUP'
      pdfval(1) = 4
      pdfpar(2) = 'NSET'
      pdfval(2) = 30
      pdfpar(3) = 'NPTYPE'
      pdfval(3) = 1
      pdfpar(4) = 'TMAS'
      pdfval(4) = mtop
      pdfpar(5) = 'NFL'
      pdfval(5) = 6
      CALL PDFSET(pdfpar,pdfval)
C
C --  Calculate total cross section and add beam particles to listing.
      CALL SUSARR( ierr )
      IF (ierr.NE.0) THEN
        ierr = ierr + 100
        GOTO 8900
      ENDIF
C
C     Stop timer.
      CALL TIMED (time)
      WRITE( isustr(1), 1008 ) time
      WRITE(isustr(1),1000) NEvent
C
 8900 CONTINUE
      RETURN
C
 1000 FORMAT(21X,'Susana run generating ', I6, ' events', /)
 1001 FORMAT(///
     & 16X, 'o----------------------------------------------o', /
     &,16X, '| R-Violating SUSY Monte Carlo Event Generator |', /
     &,16X, '|                                              |', /
     &,16X, '|       SUSANA version 1.4, 23-March-1992      |', /
     &,16X, '|                                              |', /
     &,16X, '|        Jon Butterworth & Herbi Dreiner       |', /
     &,16X, '|        Oxford HEP        Oxford Theory       |', /
     &,16X, 'o----------------------------------------------o')
 1002 FORMAT(
     & 20X,'Beam momenta in GeV:', /
     &,20X,'Proton     ', 2F6.2, F8.2, /
     &,20X,'Electron   ', 2F6.2, F8.2, /
     &,20X,'CMS energy ',  F7.2,  /)
 1022 FORMAT(
     & 20X,'Beam momenta in GeV:', /
     &,20X,'Proton     ', 2F6.2, F8.2, /
     &,20X,'Positron   ', 2F6.2, F8.2, /
     &,20X,'CMS energy ',  F7.2,  /)
C
 1004 FORMAT(
     & 20X,'Kinematic limits:', /
     &,20X,'Minimum value of x       = ', F7.4, /
     &,20X,'Maximum value of x       = ', F7.4, /
     &,20X,'Minimum value of y       = ', F7.4, /
     &,20X,'Maximum value of y       = ', F7.4, /
     &,20X,'Minimum value of Q2      = ', F7.0, /
     &,20X,'Maximum value of Q2      = ', F7.0, /)
 1111 FORMAT(
     & 18X,'SUSANA will run with the following parameters:',/
     &,20X,'Down squark mass         =',f9.4,' GeV',/
     &,20X,'Up-squark mass           =',f9.4,' GeV',/
     &,20X,'Strange-squark mass      =',f9.4,' GeV',/
     &,20X,'Charm-squark mass        =',f9.4,' GeV',/
     &,20X,'Bottom-squark mass       =',f9.4,' GeV',/
     &,20X,'Top-squark mass          =',f9.4,' GeV',/
     &,20X,'Selectron mass           =',f9.4,' GeV',/
     &,20X,'Photino mass             =',f9.4,' GeV',/
     &,20X,'Photino lifetime         =',E9.4,' s',/
     &,20X,'Top quark mass           =',f9.4,' GeV',/
     &,20X,'Coupling strength        =',f9.4,/
     &,20X,'Photino decays allowed   =',I4,/
     &,20X,'(Branching fraction      =',f9.4,')'/)
C
C    Formats for Electron
C 
 1005 FORMAT(
     & 20X,'Down squark width        =',f9.4,' GeV')
 1025 FORMAT(
     & 20X,'Up squark width          =',f9.4,' GeV')
 1006 FORMAT(
     & 20X,'Strange squark width     =',f9.4,' GeV')
 1026 FORMAT(
     & 20X,'Charm squark width       =',f9.4,' GeV')
 1007 FORMAT(
     & 20X,'Bottom squark width      =',f9.4,' GeV')
 1027 FORMAT(
     & 20X,'Top squark width         =',f9.4,' GeV')
 1008 FORMAT(20X,'Time initializing generator = ', F6.2, ' s', /)
 1009 FORMAT(20X,'None zero lambda = 111. e- + u -> d squark'/)
 1010 FORMAT(20X,'None zero lambda = 112. e- + u -> s squark'/)
 1011 FORMAT(20X,'None zero lambda = 113. e- + u -> b squark'/)
 1012 FORMAT(20X,'None zero lambda = 121. e- + c -> d squark'/)
 1013 FORMAT(20X,'None zero lambda = 122. e- + c -> s squark'/)
 1014 FORMAT(20X,'None zero lambda = 123. e- + c -> b squark'/)
 1015 FORMAT(20X,'None zero lambda = 131. e- + t -> d squark'/)
 1016 FORMAT(20X,'None zero lambda = 132. e- + t -> s squark'/)
 1017 FORMAT(20X,'None zero lambda = 133. e- + t -> b squark'/)
C
C     Formats for Positron
C
 1205 FORMAT(
     & 20X,'Up squark width        =',f9.4,' GeV')
 1225 FORMAT(
     & 20X,'Down squark width          =',f9.4,' GeV')
 1206 FORMAT(
     & 20X,'Charm  squark width     =',f9.4,' GeV')
 1226 FORMAT(
     & 20X,'Strange squark width       =',f9.4,' GeV')
 1207 FORMAT(
     & 20X,'Top squark width      =',f9.4,' GeV')
 1227 FORMAT(
     & 20X,'Bottom squark width         =',f9.4,' GeV')
 1208 FORMAT(20X,'Time initializing generator = ', F6.2, ' s', /)
 1209 FORMAT(20X,'None zero lambda = 111. e+ + d -> u squark'/)
 1210 FORMAT(20X,'None zero lambda = 112. e+ + d -> c squark'/)
 1211 FORMAT(20X,'None zero lambda = 113. e+ + d -> t squark'/)
 1212 FORMAT(20X,'None zero lambda = 121. e+ + s -> u squark'/)
 1213 FORMAT(20X,'None zero lambda = 122. e+ + s -> c squark'/)
 1214 FORMAT(20X,'None zero lambda = 123. e+ + s -> t squark'/)
 1215 FORMAT(20X,'None zero lambda = 131. e+ + b -> u squark'/)
 1216 FORMAT(20X,'None zero lambda = 132. e+ + b -> c squark'/)
 1217 FORMAT(20X,'None zero lambda = 133. e+ + b -> t squark'/)
      END
