c =======================================================================
c     A sample program that uses ojf_015 (optimal jet finder ver. 015)
c     ( works also with the obsolete version ojf_014 )
c     by D.Yu.Grigoriev, E.Jankowski and F.V.Tkachov.
c
c     For more comments see ww160.f and ww160a.f
c
c     This is for playing with Q_minimize.
c     An event taken from a file 
c     (slightly different format of event file than ww160.in: 
c     here number of particles is not specified).
c
c     A random initial jet configuration with 3 jets is specified,
c     the minimization is performed, and the results printed out.
c
        PROGRAM  demo

        INCLUDE 'ojf_kin.fh'

        CHARACTER*2000  name
        CHARACTER*1     event_type

        DOUBLE PRECISION  Radius, e, theta, phi
        DOUBLE PRECISION  o_fin, y_fin, e_fin
        INTEGER  a, seed, nparts, njets, kinematics
        LOGICAL  success

c  defaults for this demo:
        name  = 'demo.in'
        OPEN(6,  FILE = 'demo.out')
        njets = 3
        seed  = 13
        Radius = 1.0

c begin dialog to update defaults:

100     WRITE(0,*) 'default filename = ', name(1:index(name,' ')),
     &             '   Any updates?'
        READ(5,*, end=101, err=100) name
101     CONTINUE
        WRITE(0,*) 'filename : ', name(1:index(name,' '))

200     WRITE(0,'(1x,a,i2,a)') 'default no. of jets = ', njets, 
     &                      '   Any updates?'
        READ(5,*, end=201, err=200) njets
201     CONTINUE
        WRITE(0,'(1x,a,i2)') 'number of jets =', njets

299     FORMAT(1x,a,1p,g10.2,a)
300     WRITE(0,299) 'default Radius = ', Radius, '   Any updates?'
        READ(5,*, end=301, err=300) Radius
301     CONTINUE
        WRITE(0,299) 'Radius =', Radius

400     WRITE(0,'(1x,a,i7,a)') 'default seed = ', seed,
     &                         '   Any updates?'
        READ(5,*, end=401, err=400) seed
401     CONTINUE
        WRITE(0,'(1x,a,i7)') 'seed = ', seed
c
c end dialog

c Set up event from the specified file:
c
        OPEN(10, FILE = name, form = 'formatted', status = 'old')
c
c Event setup must know the kinematics 
c (sphere [c.m.s. collisions] or cylinder[hadron collisions]).
c It is assumed in this demo 
c that the file starts with a single word on the first line.
c Its first character determines kinematics -- sphere or cylinder.
c
        READ(10,*) event_type

        IF ( (event_type .EQ. 's') .OR. (event_type .EQ. 'S') )
     &  THEN
            kinematics = sphere

        ELSE IF ( (event_type .EQ. 'c') .OR. (event_type .EQ. 'C') ) 
     &  THEN
            kinematics = cylinder
        ELSE
            STOP 'Wrong event type in input file'
        END IF
c
c Start event initialization: set kinematics; reset internal stuff.
c
        CALL event_setup_begin ( kinematics )
c
c read particles
c
        nparts = 0
        DO a = 1, 1999
            READ(10,*, end=1000, err=1000) e, theta, phi
            CALL add_particle ( e, theta, phi )
            nparts = nparts + 1
        END DO
1000    CLOSE(10)
c       nparts is now the total number of particles in the event.
c
        CALL event_setup_end

        WRITE(0,'(1x,i5,a)') nparts, ' particles read'
        WRITE(6,'(1x,i5,a)') nparts, ' particles read'

        IF (nparts .LE. 0) STOP 'file without particles'

c
c set up jet configuration
c
        CALL jets_setup_begin ( njets, Radius )
c
c Set the seed of random generator (not mandatory, but recommended 
c -- different seeds mean different order of found jets):
c
        CALL set_seed ( seed )
c
c One must initialize the recombination matrix. 
c The simplest way is a complete randomization:
c
        CALL init_z_random_all
c
c Finish set up of initial jet configuration:
c
        CALL jets_setup_end

c
c And now perform minimization:
c
        CALL Q_minimize ( success )

        IF (.NOT. success) THEN
            WRITE(6,*) 'W A R N I N G :   minimum NOT FOUND.'
c           this is very unlikely unless ojf_015 was changed.
        ELSE
            WRITE(6,*) 'minimum found with'

            CALL get_criterion ( o_fin, y_fin, e_fin )
c              see comments in ww160a.f about these parameters
c
            WRITE(0,'(1x,3(a,1p,g13.5))')
     &       ' Omega =', o_fin, '  Y =', y_fin, '  Esoft =', e_fin
            WRITE(6,*)
            WRITE(6,'(1x,3(a,1p,g13.5))')
     &       ' Omega =', o_fin, '  Y =', y_fin, '  Esoft =', e_fin
        END IF
c
c Optional: sort jets in decreasing energy
c (normal/transverse energy for spherical/cylindrical kinematics):
c
        CALL sort_jets
c
c Different ways to print results. 
c Code in the procedures demonstrates different ways to get info
c about the resulting configuration of jets:
c
        WRITE(6,*)
        WRITE(6,*) '==============================================='
        WRITE(6,*) 'A primitive output of recombination matrix (z):'

        CALL print_z_raw
c
        WRITE(6,*)
        WRITE(6,*) '==============================================='
        WRITE(6,*) 'A more readable output of z:'

        CALL print_z_nice
c
        WRITE(6,*)
        WRITE(6,*) '==============================================='
        WRITE(6,*) 'Output by jet:'
        WRITE(6,*)

        CALL print_jets
c
        WRITE(6,*)
        WRITE(6,*) '==============================================='
        WRITE(6,*) 'Output by particle:'
        WRITE(6,*)

        CALL print_particles
c
        END
c
c END OF DEMO
c ==================================================================