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 jets are searched from random initial jet configurations. c c See print_* routines in ojf_015.f for how jets' parameters are accessed. c PROGRAM ww160 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_cut, & o_fin, y_fin, e_fin, & e_scale INTEGER kinematics, event, a, nparts, njets, ntries, event4j INTEGER seed c defaults for this demo: name = 'ww160.in' OPEN(6, FILE='ww160.out') c c ww160.in contains 100 events generated at 160 GeV. c c The two key control parameters of the new version of the algorithm: Radius = 1 c Radius is very similar to the R_cone c of cone algorithms. c Think about it as the maximal distance of infinitesimally soft particle from c any jet's axis above which the particle is left out of jet formation c (in the terminology of the new version: "is relegated to soft energy"). c c Quantitatively, I'd roughly relate R_cone ~ Radius/SQRT(2) c (see hep-ph/9901444, "second edition"), so the above value c roughly corresponds to the popular choice R_cone~0.7. o_cut = 0.1 c o_cut (o from omega) is related to both y_cut and the 'f'-cut. c Think about it as an upper bound on the [transverse] energy fraction c allowed to be left out of jet formation -- but this is only c for the purposes of numerical estimates c because its second (or may be first) role is as a cumulative version c of the usual y_cut. ntries = 10 c The larger ntries, the lower the probability that the algorithm c misses the global minimum. c c Changing ntries from 10 to 5 causes ~1% of events in these samples c to change the number of found jets. 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 The jet resolution parameter: c 400 WRITE(0,299) 'default o_cut = ', o_cut, ' Any updates?' READ(5,*, end=401, err=400) o_cut 401 CONTINUE WRITE(0,299) 'o_cut =', o_cut c No. of tries to repeat search c from different initial jet configurations: c 500 WRITE(0,*) 'default ntries = ', ntries, ' Any updates?' READ(5,*, end=501, err=500) ntries 501 CONTINUE WRITE(0,299) 'ntries =', ntries seed = 13 c A seed for random number generator. c Change seed arbitrarily to verify stability of results against c changes of random initial conditions. 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') c c The algorithm must know at event setup about the kinematics c -- sphere [c.m.s. collisions] or cylinder [hadron collisions]. c It is assumed in this demo that the first line of the file c contains a single word "sphere" or "cylinder" (case-insensitive). c READ(10,*) event_type c 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 For tests it is convenient to have the event file c contain info about format of events c -- energy, theta, phi or px, py, pz. c This is specified by a single word on the second line in the file c -- "angles" or "vectors" (case-insensitive). c 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 c c Once again: larger ntries --> higher probability that the minimum found c is a global minimum: CALL set_ntries ( ntries ) c c different seeds generate different random initial configurations: CALL set_seed ( seed ) c c set event counter: event = 0 c set the counter of events with 4 jets found: event4j = 0 c loop over events: DO WHILE(.TRUE.) event = event + 1 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. Read it: c READ(10,*,END=1000,err=1000) nparts 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 the same units as used for input) c algorithm deals internally with energy fractions c and knows nothing about physical units like GeV. 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 c Remember the seed used for this event -- just in case c something interesting happens so you could reproduce the results: CALL get_seed ( seed ) WRITE(6,*) 'random number generator seed for this event is:', & seed c c The jet search proper: CALL Q_search ( o_cut, Radius, njets ) c njets is the found no. of jets. CALL get_criterion ( o_fin, y_fin, e_fin ) c These are final values of o_fin as well as its components c y_fin = the value of Y, c e_fin = the value of E_soft. c This is useful for producing the Y-E_soft distribution, c which is one of the interesting new things in OJF c -- for a meaningful plot one should minimize with a fixed no. of jets. c See ww160a.f. WRITE(6,*) WRITE(6,*) 'final values: ' WRITE(6,'(1x,3(a,1p,g11.3))') ' Omega =', o_fin, & ' Y =', y_fin, ' E_soft =', e_fin IF (njets .EQ. 0) THEN WRITE(0,*) 'W A R N I N G : jets not found' WRITE(6,*) 'W A R N I N G : jets not found' c this means algorithm could not achieve Omega < o_cut c using the default maximal no. of jets (20) built into algorithm c (the parameter njets_max in params.fh) c ELSE WRITE(0,'(50x,i3,a)') njets, ' jets found' WRITE(6,'(1x,i3,a)') njets, ' jets found' c different ways to output results -- see code in the corresponding subroutines c (in ojf_015.f) for how data are accessed. c c Optional: sort jets in decreasing energy c (normal/transverse energy for spherical/cylindrical kinematics): c CALL sort_jets c c A readable output of z:' CALL print_z_nice c c Output by jet: CALL print_jets c c Output by particle: CALL print_particles IF (njets .EQ. 4) THEN event4j = event4j + 1 END IF END IF END DO 1000 CONTINUE WRITE(0,*) event4j, ' events with 4 jets found' WRITE(6,*) event4j, ' events with 4 jets found' STOP END c