      PROGRAM SUSANA
C ======================================================================
C Main program for standalone event generation.              JMB
C
C Error codes -> ierr < 0 - skip an event
C                ierr > 0 - End run (fatal error)
C                ierr = 0 - all ok, carry on. 
C ======================================================================
C
      IMPLICIT NONE

      CHARACTER*16 OP_NAME
      INTEGER I,LENOCC

C --  Set default parameters (access block data).
      EXTERNAL SUSDAT
      EXTERNAL LUDATA
      EXTERNAL PYDATA
C
      INCLUDE 'SUSANA.INC'
      INCLUDE 'PYTHIA57.INC'

      INTEGER NEV_NT
      REAL
     &     DISTANCE_NT, PPHINO_NT, SM_NT, TM_NT, UM_NT, XJB_NT,YJB_NT,
     &     XJ_NT, YJ_NT, XT_NT, YT_NT, EP_NT, RM_NT, ET_NT, PH_NT
      COMMON/NTUP /
     &     NEV_NT, DISTANCE_NT, PPHINO_NT(4), SM_NT, TM_NT, UM_NT, 
     &     XJB_NT, YJB_NT, XJ_NT, YJ_NT, XT_NT, YT_NT, EP_NT, RM_NT, 
     &     ET_NT, PH_NT


      INTEGER part, npos

      REAL
     &     pld1, pld2, pld3, pld4, plm, pldm,
     &     costh, costhe, cossqby2, sinsqby2,
     &     phipld, polpld, etapld, rmin, phip, etap, polp,
     &     dphi, deta, rad, ecut,sume, sumpz, sumpy,sumpx, 
     &     xjb, yjb, qsqjb

      PARAMETER (ecut = 1.0)

C ----------------------> End of common.

      INTEGER IErr
      INTEGER Evt

      INTEGER Nwds_HBOOK
      PARAMETER (Nwds_HBOOK=10000000)
      REAL          HMEM
      COMMON /PAWC/ HMEM(NWDS_HBOOK)

C     Get Event no
      WRITE( isustr(1),*) 'SUSANA: Number of events... ?'
      READ( isustr(2),*) NEvent

C --  Initialize generator.
      CALL SUSINT( ierr )
      IF (ierr.NE.0) THEN
C       Any error is fatal
        GOTO 7900
      ENDIF

C     Book histograms.
      CALL HLIMIT (NWDS_HBOOK)

C     Create output file.
      OP_NAME='test.his'
      WRITE(*,*) 'HROPEN OPENS ',OP_NAME(1:LENOCC(OP_NAME))
      CALL HROPEN(69,'OUTPUT',OP_NAME(1:LENOCC(OP_NAME)),'NQ',1024,IERR)
      
C     Go to the output directory.
      CALL HCDIR('//PAWC',' ')
      CALL HCDIR('//OUTPUT',' ')

C --  Book column-wise ntuple
      CALL HBNT(9999,'SUSY ntuple','D')
      
      CALL HBNAME(9999,'Susana',NEV_NT,
     & 'NEVT,DISTANCE,PPHINO(4),SM,TM,UM,XJB,YJB,XJ,YJ,XT,'
     & //'YT,EP,RM,ET,PH')

C     Event loop.
      DO 1 evt = 1, nevent
C     
         WRITE(isustr(1),*) 'Start generating event ',event

C     Generate an event.
         CALL SUSGEN( ierr )
         IF (ierr.GT.0) THEN
C     Fatal error. Terminate run.
            GOTO 7900
         ENDIF

C     Branch here for final states containing neutrinos

         IF (isustr(6).ne.1.and.isustr(6). ne.10) go to 30

C      Go Through LUJET array to find positron from photino decay
 
         pld1=1.
         pld2=1.
         pld3=1
         pld4=1.
         npos = 0

         Do 2 part=1,n
         IF(k(part,1).ne.1) go to 2
         If(k(part,2).eq.-11) then
                  npos = part
                  pld1=p(part,1)
                  pld2=p(part,2)
                  pld3=p(part,3)
                  pld4=p(part,4)
              endif
 2       continue
         If(npos.eq.0) Write(*,*) 'Positron not found'
         if(npos.eq.0) go to 1

C        Determine phi and eta of positron

         phipld = atan2(pld2,pld1)
         if(pld2.lt.0.0) phipld = phipld + 2.0*pi
         if(pld3.eq.0.0) then
              etapld = 0.0
         else
              polpld = atan(sqrt(pld1**2 +pld2**2)/pld3)
              if(polpld.lt.0.0) polpld = polpld + pi 
              etapld = -log(tan(polpld/2.0))
         endif

c        Now determine the minimum etaphi radius called rad

         rmin = 1.0E9
         Do 3 part =1,n
          IF(k(part,1).ne.1)go to 3
             if(part.eq.npos) go to 3
C        Determine phi and eta of particle
             if(p(part,4).lt.ecut) go to 3
             if(p(part,1).eq.0.0.and.p(part,2).eq.0.0) go to 3
             phip = atan2(p(part,2),p(part,1))
             if(p(part,2).lt.0.0) phip = phip + 2.0*pi
             if(p(part,3).eq.0.0) then
               etap = 0.0
             else
               polp = atan(sqrt(p(part,1)**2 +p(part,2)**2)/p(part,3))
               if(polp.lt.0.0) polp = polp + pi
               etap = -log(tan(polp/2.0))
             endif
             dphi = phip - phipld
             if(dphi.gt.pi) dphi = 2.0*pi - dphi
             deta = etap - etapld
             rad  = sqrt(dphi**2 + deta**2)
             if(rad.lt.rmin)  rmin = rad
 3       continue
         NEV_NT=EVT
         DISTANCE_NT=DISTANCE 
         DO I=1,4
            PPHINO_NT(I)=PPHINO(I)
         ENDDO
         SM_NT=SM
         TM_NT=TM 
         UM_NT=UM
         XJ_NT=XJ
         YJ_NT=YJ
         EP_NT=pld4
         RM_NT=rmin
         ET_NT=etapld
         PH_NT=phipld

C     Calculate x and y from the incoming electron and the positron
C     from the photino.

C     Calculate the cosine of the apparent scattering angle.

         plm = sqrt(pl(1)**2 +pl(2)**2 + pl(3)**2)
         pldm = sqrt(pld1**2 +pld2**2 +pld3**2)                  
         costh=(pl(1)*pld1 + pl(2)*pld2 +pl(3)*pld3)/(plm*pldm)
         costhe = -costh
         cossqby2 = 0.5*(1.0 +costhe)
         sinsqby2 = 0.5*(1.0 -costhe)
         xt_nt= (pld4/pp(4))*cossqby2 / (1.0 - (pld4/pl(4))*sinsqby2)
         yt_nt= 1.0 - (pld4/pl(4))*sinsqby2

C     Photino modes with neutrinos com here for Jacquet-Blondel Kinematics.

 30      Continue

         sume = 0.0
         sumpz = 0.0
         sumpy = 0.0
         sumpx = 0.0

         Do 20 part = 1,n
         If(k(part,1).ne.1) go to 20  
         If(k(part,2).eq.12) go to 20
         sume = sume + p(part,4)
         sumpz = sumpz + p(part,3)
         sumpy = sumpy + p(part,2)
         sumpx = sumpx + p(part,1)
 20      continue

C     Now calculate via the Jacquet-Blondel method.

      yjb = (sume - sumpz)/(2.0*pl(4))
      if(yjb. eq. 1.0) go to 1
      qsqjb = (sumpx**2 +sumpy**2)/(1.0 - yjb)
      xjb = qsqjb/(2.0*pp(4)*(sume - sumpz))
      XJB_NT = xjb
      YJB_NT = yjb
 
C     Fill ntuple
         CALL HFNT(9999)
         
C     Print out first event.
         IF (Evt.LT.10) CALL LULIST (1)
C     
 1    CONTINUE
C     
 7900 CONTINUE
      WRITE( isustr(1),1002) Evt 

      CALL HCDIR('//OUTPUT',' ')
      CALL HISTDO
      CALL HROUT( 0,I, ' ' )
      CALL HREND( 'OUTPUT' )
      CLOSE (69)
C     
      STOP
C     
 1002 FORMAT(/,18X,'Susana run generated ', I6, ' events', /)
      END
