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 processing events generated by pythia for ee-->WW-->4 jets c (thanks to A.Skasyrskaya) c c for each event we here compute "fuzziness with 4 axes" c (see hep-ph/9901444, "second edition") c for the purpose of plotting the Y-Esoft distribution c which is one of the interesing by-products of the theory behind ojf. c c See also the accompanying sample program ww160.f. c PROGRAM ww160a c two parameters to set kinematics -- sphere and cylinder: c INCLUDE 'ojf_kin.fh' CHARACTER*2000 name CHARACTER*1 event_type, event_format DOUBLE PRECISION px, py, pz, & Radius, & o_fin, y_fin, e_fin, & e_scale, omega, y, esoft INTEGER ntries, seed, kinematics, event, a, nparts, njets INTEGER succs LOGICAL success name = 'ww160.in' OPEN(6, FILE='ww160a.out') njets = 4 c in this example we fix the number of jets and minimize c the criterion and look at the values of Y and E_soft. Radius = 1 c omega_cut is not set because we simply find the minimum. seed = 13 ntries = 10 c c begin dialog to update defaults: c c Input file name: c 100 WRITE(0,'(1x,3a)') 'default filename = ', & name(1:index(name,' ')), ' Any updates ?' READ(5,*, end=101, err=100) name 101 CONTINUE WRITE(0,*) 'filename : ', name(1:index(name,' ')) c The R parameter: c 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 c==c No. of tries to repeat search c==c from different initial jet configurations: c==c c==500 WRITE(0,*) 'default ntries = ', ntries, ' Any updates?' c== READ(5,*, end=501, err=500) ntries c==501 CONTINUE c== WRITE(0,299) 'ntries =', ntries c c Default seed (useful to vary initial configurations): c 600 WRITE(0,*) 'default seed = ', seed, ' Any updates?' READ(5,*, end=601, err=600) seed 601 CONTINUE WRITE(0,299) 'seed =', seed c c end dialog c Set up event from the specified file: c OPEN(10, file=name, form='formatted', status='old') 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 READ(10,*) event_format IF ((event_format .EQ. 'a') .OR. (event_format .EQ. 'A')) THEN event_format = 'a' c angles (energy, theta, phi) ELSE IF ((event_format .EQ. 'v') .OR. (event_format .EQ. 'V')) & THEN event_format = 'v' c vectors (px, py, pz) ELSE STOP 'Wrong event format in input file' END IF CALL set_ntries ( ntries ) CALL set_seed ( seed ) c c loop over all events in the file c event = 0 DO WHILE(.TRUE.) c c Start event initialization: set kinematics; reset internal stuff c CALL event_setup_begin ( kinematics ) c c each event in this demo file starts with the number of particles c on a separate line (the way pythia produced them). Read it: c READ(10,*,END=1000,err=1000) nparts event = event + 1 WRITE(0,*) 'event', event,': nparts =', nparts c c read particles c IF (nparts .LE. 0) STOP 'empty event' DO a = 1, nparts READ (10,*,END=1000,err=1000) px, py, pz IF (event_format .EQ. 'a') THEN CALL add_particle ( px, py, pz ) ELSE IF (event_format .EQ. 'v') THEN CALL add_particle_raw ( px, py, pz ) ELSE STOP 'invalid event format specified' END IF END DO c nparts is the label of this particle c CALL event_setup_end c "closing bracket" for event setup. CALL get_nparts ( nparts, e_scale ) c this returns: c nparts = no. of particles in the event c e_scale = total energy of the event c (in absolute units) c algorithm deals internally with energy fractions. WRITE(6,'(1x,79(''=''))') WRITE(6,*) WRITE(6,'(1x,(a,i5,a,i3,a))') 'Event', event,': ', & nparts, ' particles read' IF ( kinematics .EQ. sphere ) THEN WRITE(6,'(2x,a,1p,g13.5)') 'total energy =', e_scale ELSE IF ( kinematics .EQ. cylinder ) THEN WRITE(6,'(2x,a,1p,g13.5)') & 'total transverse energy =', e_scale ELSE STOP 'invalid kinematics' END IF cccccccccccccccccc c A loop is needed because Q_minimize does not set c initial configuration to start minimum search from. c So we need to randomly change it a few times. c o_fin = 10000. succs = 0 c counter for protection against no-local-minima-found condition c DO iter=1, ntries c CALL jets_setup_begin( njets, Radius ) c two parameters have to be specified. c c remember the seed used for this minimization attempt -- just in case c something interesting happens so we could repeat it: CALL get_seed ( seed ) c completely random initial configuration to start minimization from CALL init_z_random_all c it can also be set selectively; see comments in ojf_015.f CALL jets_setup_end c CALL Q_minimize ( success ) IF (success) THEN succs = succs + 1 CALL get_criterion ( omega, y, esoft ) c The found values of Omega (omega) and its components: c y = the value of Y, c esoft = the value of E_soft. IF (omega < o_fin) THEN o_fin = omega y_fin = y e_fin = esoft END IF ELSE c the algorithm failed to descend into a local minimum c from the specified starting point. extremely unlikely END IF END DO IF (succs .GT. 0) THEN WRITE(6,'(1x,3(a,1p,g11.3))') & ' Omega =', o_fin, ' Y =', y_fin, ' E_soft =', e_fin WRITE(0,'(1x,3(a,1p,g11.3))') & ' Omega =', o_fin, ' Y =', y_fin, ' E_soft =', e_fin ELSE c this branch is extremely unlikely WRITE(0,*) 'W A R N I N G : OJF failed to find any minima' WRITE(0,*) & 'Please send the event in text form to the authors!' WRITE(6,*) 'W A R N I N G : OJF failed to find any minima' WRITE(6,*) & 'Please send the event in text form to the authors!' END IF cccccccccccccccccc END DO 1000 CONTINUE c end the loop over events. STOP END c