--------------------------------------- This is left just in case. Use file ojf_015 instead!!!! 2004-apr-13 --------------------------------------- --------------------------------------- --------------------------------------- C cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c c c OPTIMAL JET FINDER [ojf_014] c c c c (C) 1999-2001 D.Yu. Grigoriev and F.V.Tkachov, c c Institute for Nuclear Research of Russian Academy of c c Sciences, 60th October Ave., 7a, Moscow 117312, Russia. c c Permission is granted for anyone to use this software c c for any purpose on any computer. The authors are not c c responsible for the consequences of such use. c c This software is meant to be freely distributed. c cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c c The copyright formula is borrowed from G.P.Lepage's c VEGAS code. c c Physical background is described in: c F.V.Tkachov, "The definition of jets" hep-ph/9901444. c c NB The code is based on the revised version of the paper c (so-called "second edition", posted January, 2000) c which contains an important correction to the criterion c described in the Spring postings. c The correction is technically only a fine-tuning BUT c has important implications for both the definition c [direct connection to cone algorithms and thrust] c and for the resulting algorithm [which has become c much simpler, much more robust, and O(100) faster]. c c Implementation is described in: c D.Yu.Grigoriev and F.V.Tkachov, hep-ph/9912415. c See http://www.inr.ac.ru/~ftkachov/projects/jets/index.htm c cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c Authors thank Pablo Achard for beta-testing the code c c with real data and for useful suggestions not all of c c which we were able to implement. c cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c c ====================================================================== c USAGE c c See accompanying sample programs, sample subroutines below, c and comments after each subroutine heading. c c The key routine is Q_minimize. c None of the subroutines following it are user-callable. c c On the other hand, Q_minimize can be used in a variety of ways, c and Q_search implements only the simplest jet search. c c Regard Q_search as the simplest "driver" for Q_minimize. c Some special applications may require special such drivers. c c ====================================================================== c ACCESS TO DATA AND CONTROL PARAMETERS c c Normally, you will only INCLUDE ojf_kin.fh into your code. c c You should be able to do everything you may need c (short of modifications of the core minimization algorithm) c via special interface routines (set_... and get_...) c which ensure correctness of some important internal invariants. c c DO NOT write to the COMMON blocks directly. c c ====================================================================== c EXTERNAL OUTPUT c c All output is done via WRITE(6,... c c There are two kinds of such output: c sample printouts from subroutines print_..., c error messages from anywhere. c c If you need to redirect error messages, use text editor. c If you need a different (printed) output, write your own subroutines c using those provided here as examples. c c ====================================================================== c ERROR MESSAGES c c Message numbers 20..29 are due to user errors. c Numbers >= 30 are generated by failed internal checks, c so should you get such a message, please inform the authors; c please include the corresponding event in text form. c c ====================================================================== c DOUBLE PRECISION c c is used everywhere for floating point. c c We decided not to attempt implementing single precision in interfaces c because of lack of resources at our disposal. c c ====================================================================== c LIST OF USER-CALLABLE SUBROUTINES FOR OJF_014 c c Event setup: c event_setup_begin ( kinematics ) c add_particle ( energy, theta, phi ) c add_particle_raw ( px, py, pz ) c event_setup_end c c Setup of initial jet configuration: c jets_setup_begin ( njets, Radius ) c set_seed ( seed ) c reset_Radius ( Radius ) c init_z_random_all c assign_to_jet ( a, j ) c init_z_from ( a, z_in ) c init_z_random ( a ) c jets_setup_end c c Setting algorithm control parameters: c set_maxiter ( maxiter ) c set_njets_limits ( nstart, nstop ) c set_ntries ( n ) c set_trace_nmoved ( bool ) c c Access to parameters: c get_kinematics ( kinematics ) c get_nparts ( nparts, e_scale ) c get_particle ( a, e, xta, phi, p, ephys, pphys ) c get_njets ( njets ) c get_seed ( seed ) c get_Radius ( R ) c get_maxiter ( maxiter ) c get_njets_limits ( nstart, nstop ) c get_ntries ( n ) c c Access to resuls: c get_criterion ( omega, y, esoft ) c get_jet ( j, e, xta, phi, q, qtilde, ephys, qphys ) c get_z ( a, z_out ) c get_particle_split ( a, total_jets, jet, zj ) c get_jet_split ( j, nwhole, whole_a, nfract, fract_a, fract_z ) c c sort_jets c c Sample print routines, show how resulting jet configuration is accessed: c print_z_raw c print_z_nice c print_jets c print_particles c c Simplest example of straightforward jet search: c Q_search ( w_cut, Radius, njets ) c c Key minimization routine: c Q_minimize ( success ) c c Other routines are not user-callable. c c ====================================================================== c VERSION HISTORY c c Proof-of-concept algorithm in Component Pascal c by F.V.Tkachov, ftkachov@ms2.inr.ac.ru c c Fortran port, revisions 001-012 with algorithm fine-tunings c by D.Yu.Grigoriev, dgr@inr.ac.ru, Dmitry.Grigoriev@cern.ch c c Revision 013 by F.V.Tkachov. c c Beta tests of rev. 013 by Pablo.Achard@cern.ch c This public code does not contain any changes in the algorithm c or interfaces compared with the code used by P.Achard c except for: c -- a change in the value of an internal parameter c used for internal checks (this change does not affect results c in any way, only our expectations about typical size of c rounding errors in a sum which must be, ideally, zero, c and which therefore is used for an internal check); c -- replacement of the two corresponding STOP c statements by WRITE... 'W A R N I N G... c A few cosmetic changes have been done too. c c Revision 014 (official release) by F.V.Tkachov. c c ====================================================================== c version 0.14 28-oct-2001 --------------------------------------------- c c The third line of comments after Q_search now reads: c "It is likely..." instead of "It is not likely..." c c A correction at end of set_seed, forcing ojf_random c to be in the correct range. c c z_snap: check of input is moved out from the main loop. c c In order to fully agree with the main theory paper, c ojf_Radius2 is now Radius**2 c (instead of Radius**2/2 in version 013), c and eval_Y, grad_Y changed accordingly (now include factor 2). c These changes affect neither Omega nor the minima found c but only how Omega is decomposed into Y and Esoft. c The paper by Grigoried and Tkachov (hep-ph/9912415) c has been updated accordingly. c c eps_sum is now replaced by three parameters: c eps_sum, eps_sum0, eps_sum1 to better control rounding errors. c c tested: c MS PowerStation/Win98 c c version 0.13a 30-march-2000------------------------------------------- c c Microbugs in output routines removed (only manifest if Npart < Njets. c c version 0.13 10-sept-99 ---------------------------------------------- c c implemented the fine-tuned criterion as described in c the "second edition" of hep-ph/9901444 (with E_soft included additively); c the resulting multitudinous simplifications incorporated; c the corresponding changes of interface made; c c total zero-warning code clean-up [defies description]; c complete revision of every aspect of the code c -- numerics, interface routines, structure... [defies description]; c most notably, stricter and more differential control of rounding errors. c All within the fortran subset employed in previous versions by DG. c c tested: c c MS PowerStation/Win98 c large-scale real-life beta-tests by Pablo.Achard@cern.ch c c version 0.12 05-may-99 ----------------------------------------------- c c tested: c c g77/Linux on Intel, Alpha and Sparc platforms c Sun f77/Solaris (thwgs.cern.ch) c xlf/AIX (rsplus.cern.ch) c c z_assert is using eps_kin, not eps for accuracy checks; c major restructuring of interface; c interface modules now stop on all discrepancies instead of trying c to continue; c Q_minimize and Q_search no longer sort jets (explicit call c required [Dima, ty ne prav! no niqego, ja ispravil. -- FT]); c init_z_from corrected (now calls z_snap to set z_border array); c ojf_complete_search tries all jet numbers at least twice c c version 0.11 07-apr-99 ----------------------------------------------- c c eps returned to fixed value (10^(-3)); c deadlock processed by switching z_snap off; c z_assert split into full and light (for deadlock mode) versions; c consolidate_emiss introduced in deadlock mode to eliminate c artificially small coefficients z0; c abandoned attempts to predict step in minimize_wrt (appear to be c the reason for deadlocks) c c version 0.10 31-mar-99 ----------------------------------------------- c c more reorganizations; update_jets improved; sort_jets added; c g_eval and j_eval_nonlinear updated to correct processing of empty jets; c in cylindrical events E_miss now is a fraction of E_transverse, c not E_total; ojf_complete_search added; c value eps (z_snap, z_assert, d_assert) depends on prec_emiss c (removed in ver. 0.11) c c version 0.08 18-mar-99 ----------------------------------------------- c c modules add_raw_particle and p_set_raw added; c process of adding and normalizing particles reorganized; c code made reentrant (processing of several events is now possible) c c version 0.06 17-mar-99 ----------------------------------------------- c c tested: c c g77/Linux on Intel, Alpha and Sparc platforms c Sun f77/Solaris (thwgs.cern.ch) c xlf/AIX (rsplus.cern.ch) c MS Fortran PowerStation 4.0/Win98 c c interface subroutines added; c text rearranged; Pascal source removed; c correct algorithm used in set_z_random; c more optimization added in Q_minimize_wrt (effect ~ 20%, removed in ver.11) c optimizations in jets reevaluation c c version 0.01 09-mar-99 ----------------------------------------------- c c strictly follows the Component Pascal code by F.V.Tkachov c with one correction in CompleteSearch module (now in main c routine) -- Y_penalty should be reevaluated after increasing lambda c c ============================================================================= c c The code below consists of three parts (not counting the BLOCK DATA): c c Part 1. Sample subroutines. Simplest jet search, printing results c Part 2. Interface subroutines c Part 3. Core stuff c c ====================================================================== c ====================================================================== c c Part 1. SAMPLE SUBROUTINES. Simplest jet search, printing results c c ====================================================================== c The simplest jet search. c SUBROUTINE Q_search ( w_cut, Radius, njets ) c c This is only the simplest default. c It is likely that you may want to modify this c -- e.g. when trying to do something with local minima. c This subroutine uses only interface routines; c it does not access internal ojf data. c c It tries to find c the configuration of jets which minimizes Omega and ensures that c Omega < w_cut with a minimal njets starting from the number of jets c previously set via set_njets_start (usually the same for all events). c c For each njets, search is repeated ntries times, each time with a c different random initial z, and the configuration with the lowest c value of Omega is retained as a result. c c Failure of search is signaled by the condition njets = 0. c c Note that Q_search randomizes initial z, c so it is meaningless to use it if one wants to specify c the initial configuration for z. c c If you need to specify configuration of jets, c you probably only need Q_minimize. c c Other control options could be: c to continue attempts until a specified number of attempts fails to yield c a better configuration; c to stop search for new minimum if say three first random initial configurations c yielded the same configurations (the event has a single local minimum which is c automatically the global one; this is quite likely so if CPU time is an ussue..) c DOUBLE PRECISION w_cut, Radius INTEGER njets INCLUDE 'ojf_kin.fh' INCLUDE 'ojf_par.fh' INCLUDE 'ojf_com.fh' LOGICAL success INTEGER a, try, nstop, ntries DOUBLE PRECISION o_best, omega, dummy, & z_best( 0 : njets_max, 1 : nparts_max ) IF ( (w_cut .LE. 0) .OR. (Radius .LE. 0) ) THEN WRITE(6,*) 'Q_search: 21: invalid parameter values' WRITE(6,*) ' w_cut = ', w_cut, ', Radius = ', Radius STOP 'Q_search: 21' END IF o_best = inf CALL get_njets_limits ( njets, nstop ) CALL get_nparts ( nparts, dummy ) CALL get_ntries ( ntries ) DO WHILE ((njets .LE. nstop) .AND. (o_best .GE. w_cut)) DO try = 1, ntries CALL jets_setup_begin ( njets, Radius ) CALL init_z_random_all CALL jets_setup_end CALL Q_minimize ( success ) IF (success) THEN CALL get_criterion ( omega, dummy, dummy ) IF ( omega .LT. o_best ) THEN o_best = omega IF ( o_best .LT. w_cut ) THEN ccccccccc save the better configuration: DO a = 1, nparts CALL get_z ( a, z_best( 0, a ) ) END DO END IF END IF END IF END DO njets = njets + 1 END DO IF (njets .GT. nstop) THEN njets = 0 ELSE c c restore the best configuration c njets = njets - 1 CALL jets_setup_begin ( njets, Radius ) DO a = 1, nparts CALL init_z_from ( a, z_best( 0, a ) ) END DO CALL jets_setup_end END IF RETURN END c c ====================================================================== c Generic printing routines. All WRITE to unit=6 c SUBROUTINE print_z_raw INCLUDE 'ojf_par.fh' INTEGER a, j, nparts, njets DOUBLE PRECISION z( 0 : njets_max ), enorm CALL get_njets ( njets ) CALL get_nparts( nparts, enorm ) WRITE(6,*) WRITE(6,*) WRITE(6,*) 'Raw output of z by particle label (a):' WRITE(6,*) WRITE(6,'(26x,a)') 'j e t n u m b e r s' WRITE(6,'(1x,a,999(i7,5x))') ' a background',(j,j=1,njets) WRITE(6,'(1x,a,999(a))') '--- ',(' -----------',j=0,njets) DO a = 1, nparts c c get coefficients z(0:njets) for particle a; z(0) -- missing energy c CALL get_z( a, z ) c c ... and immediately print them ! c WRITE(6,'(1x,i3,2x,999(1p,g12.5))') a,(z(j),j=0,njets) END DO WRITE(6,*) RETURN END c c ---------------------------------------------------------------------- c SUBROUTINE print_z_nice INCLUDE 'ojf_par.fh' INTEGER a, j, nparts, njets DOUBLE PRECISION z( 0 : njets_max ), zj, enorm CHARACTER*12 place_holder( 0 : njets_max ) CALL get_njets ( njets ) WRITE(6,*) WRITE(6,*) WRITE(6,*) 'recombination matrix z by particle label a:' WRITE(6,*) WRITE(6,'(26x,a)') 'j e t n u m b e r s' WRITE(6,'(1x,a,999(i7,5x))') ' a background',(j,j=1,njets) WRITE(6,'(1x,a,999(a))') '--- ',(' -----------',j=0,njets) CALL get_nparts( nparts, enorm ) CALL get_njets( njets ) DO a = 1, nparts CALL get_z( a, z ) c fill out place_holder(..) for later print: DO j = 0, njets zj = z(j) IF (zj.EQ.0) THEN WRITE( place_holder(j), '(a)' ) ' - ' ELSE IF (zj.EQ.1) THEN WRITE( place_holder(j), '(a)' ) ' 1. ' ELSE WRITE( place_holder(j), '(1p,g12.5)' ) z(j) END IF END DO WRITE(6,'(1x,i3,2x,999(a))') a, (place_holder(j),j=0,njets) END DO WRITE(6,*) RETURN END c c ---------------------------------------------------------------------- c SUBROUTINE print_jets INCLUDE 'ojf_kin.fh' INCLUDE 'ojf_par.fh' INTEGER a, j, nparts, njets, nwhole, nfract, total_parts, kin DOUBLE PRECISION e, xta, phi, enorm, ephys DOUBLE PRECISION j_e_total DOUBLE PRECISION q(0:3), qtilde(0:3), qphys(0:3) DOUBLE PRECISION fract_z( 1 : nparts_max ) INTEGER fract_a( 1 : nparts_max ), whole_a( 1 : nparts_max ) CALL get_kinematics( kin ) CALL get_njets( njets ) CALL get_nparts( nparts, enorm ) WRITE(6,*) IF ( kin .EQ. sphere ) THEN WRITE(6,'(1x,a,i3,a)') 'SPHERE:', njets,' jets processed' ELSE IF ( kin .EQ. cylinder ) THEN WRITE(6,'(1x,a,i3,a)') 'CYLINDER:', njets,' jets processed' END IF WRITE(6,*) WRITE(6,*) 'Configuration by jet:' WRITE(6,*) IF ( kin .EQ. sphere ) THEN WRITE(6,'(1x,5a)') 'jet ',' E ', & ' E(%) ',' theta ', & ' phi ' ELSE IF ( kin .EQ. cylinder ) THEN WRITE(6,'(1x,5a)') 'jet ',' E_tr', & ' E_tr(%) ',' eta ', & ' phi ' END IF WRITE(6,'(1x,5a)') '--- ','-----------', & ' ------- -','------ ', & ' ------- ' j_e_total = 0 DO j = 1, njets c c get jet configuration (energy/tr_energy, theta/pseudorapidity c and phi; all angles measured in degrees) c CALL get_jet( j, e, xta, phi, q, qtilde, ephys, qphys ) j_e_total = j_e_total + e WRITE(6,'(1x,i3,3x,1p,g11.4,0p3f10.4)') & j, ephys, e * 100, xta, phi END DO WRITE(6,*)'---------------------------' WRITE(6,'(1x,a6,1p,g11.5,0pf10.4)')'TOTAL:', & j_e_total * enorm, j_e_total * 100 WRITE(6,*) WRITE(6,*) 'Particle content by jet:' DO j = 1, njets c c get jet kinematic parameters again c CALL get_jet( j, e, xta, phi, q, qtilde, ephys, qphys ) c c get particle distribution in jet (whole particles sorted c by label a, others -- by weight coefficient z) c CALL get_jet_split ( j, nwhole, whole_a, & nfract, fract_a, fract_z ) total_parts = nwhole + nfract c c now come to making printout c WRITE(6,*) WRITE(6,'(1x,2(a,i3),a)') 'jet label', j, & ' (',total_parts,' particle(s) ):' IF ( kin .EQ. sphere ) THEN WRITE(6,'(1x,a,1p,g12.4)') ' E(%) =', e * 100 WRITE(6,'(1x,a,1p,g12.4)') ' theta =', xta ELSE IF ( kin .EQ. cylinder ) THEN WRITE(6,'(1x,a,1p,g12.4)') ' E_tr(%) =', e * 100 WRITE(6,'(1x,a,1p,g12.4)') ' eta =', xta END IF WRITE(6,'(1x,a,1p,g12.4)') ' phi =', phi IF (nwhole .GT. 0) THEN WRITE(6,'(1x,i3,a,9999(i4))') nwhole, & ' particle(s) in jet as a whole:', & (whole_a(a), a=1, nwhole) END IF IF (nfract .GT. 0) THEN WRITE(6,'(1x,i3,a,9999(i4,a,1p,g11.4,a))') nfract, & ' particle(s) partially in jet:', (fract_a(a), ' [', & fract_z(a), '];', a=1, nfract) END IF END DO c c and the soft energy. It's just jet j=0 but get_jet_data( jet=0, ... ) won't work ! c CALL get_jet_split ( 0, nwhole, whole_a, & nfract, fract_a, fract_z ) total_parts = nwhole + nfract WRITE(6,*) WRITE(6,'(1x,a,i3,a)') & 'soft energy (', total_parts, ' particle(s) ):' IF (nwhole .GT. 0) THEN WRITE(6,'(1x,i3,a,9999(i4))') nwhole, & ' whole particle(s) in soft energy:', & (whole_a(a), a=1, nwhole) ELSE WRITE(6,*) 'no whole particles in soft energy' END IF IF (nfract .GT. 0) THEN WRITE(6,'(1x,i3,a,9999(i4,a,1p,g11.4,a))') nfract, & ' particle(s) partially in soft energy:', ( fract_a(a), & ' [', fract_z(a), '];', a = 1, nfract ) ELSE WRITE(6,*) 'no particles partially in soft energy' END IF RETURN END c c ---------------------------------------------------------------------- c SUBROUTINE print_particles INCLUDE 'ojf_kin.fh' INCLUDE 'ojf_par.fh' INTEGER a, j, nparts, kin, total_jets DOUBLE PRECISION e, xta, phi, p(0:3), p_e_total, enorm DOUBLE PRECISION ephys, pphys(0:3) INTEGER jet_labels ( njets_max + 1 ) DOUBLE PRECISION jet_weights( njets_max + 1 ) CALL get_kinematics( kin ) CALL get_nparts( nparts, enorm ) WRITE(6,*) WRITE(6,*) 'Configuration by particle:' WRITE(6,*) WRITE(6,'(40x,a)') '(soft energy is denoted as jet=0)' IF ( kin .EQ. sphere ) THEN WRITE(6,'(1x,6a)') ' a ',' E ', & ' E(%) ',' theta ', & ' phi ','jet [ fraction ]; ...' ELSE IF ( kin .EQ. cylinder ) THEN WRITE(6,'(1x,6a)') ' a ',' E_tr', & ' E_tr(%) ',' eta ', & ' phi ','jet [ fraction ]; ...' END IF WRITE(6,'(1x,6a)') '--- ','-----------', & ' ------- -','------ ', & ' ------- ','---' p_e_total = 0 DO a = 1, nparts c c getting particle parameters... c CALL get_particle ( a, e, xta, phi, p, ephys, pphys ) p_e_total = p_e_total + e c c and its split between jets (total_jets is total number of jets that c include this particle) c CALL get_particle_split ( a, total_jets, & jet_labels, jet_weights ) c c now a pretty otput c IF (total_jets .EQ. 1) THEN c c particle completely belongs to one of the jets (or soft energy) c WRITE(6,'(1x,i3,3x,1p,g11.4,0p,3f10.4,i6)') a, ephys, & e * 100, xta, phi, jet_labels(1) ELSE c c particle is split between several jets (and/or soft energy) c WRITE(6,'(1x,i3,3x,1p,g11.4,0p,3f10.4,3x,999(i3,a,1p,g11.4,a))') & a, ephys, e * 100, xta, phi,( jet_labels(j), ' [', & jet_weights(j), '];', j=1, total_jets ) END IF END DO WRITE(6,*)'---------------------------' WRITE(6,'(1x,a6,1p,g11.5,0p,f10.4)') 'TOTAL:', & p_e_total * enorm, p_e_total * 100 WRITE(6,*) RETURN END c c end of print routines c c ====================================================================== c c Part 2. INTERFACE MODULES c c ====================================================================== c Event setup c SUBROUTINE event_setup_begin ( kinematics ) c c Begins initialization of a new event c -- from the very beginning must know kinematics c (sphere [c.m.s. collisions] or cylinder [hadron collisions]). c c kinematics ought to be set once for all events in a job c but one must draw a line somewhere... c INTEGER kinematics INCLUDE 'ojf_kin.fh' INCLUDE 'ojf_par.fh' INCLUDE 'ojf_com.fh' ojf_pi180 = 4 * DATAN( one ) / 180 IF (ojf_event_begin) THEN WRITE(6,*) 'event_setup_begin: 20: wrong call sequence' STOP 'event_setup_begin: 20' END IF IF (.NOT.( ( kinematics .EQ. sphere ) .OR. & ( kinematics .EQ. cylinder ) )) THEN WRITE(6,*) 'event_setup_begin: 21: illegal kinematics' STOP 'event_setup_begin: 21' END IF ojf_kinematics = kinematics ojf_nparts = 0 ojf_e_scale = 0.0 ojf_event_begin = .TRUE. ojf_event_set = .FALSE. ojf_jets_begin = .FALSE. ojf_jets_set = .FALSE. RETURN END c c --------------------------------------------------------------------- c SUBROUTINE add_particle ( energy, theta, phi ) c c Adds a particle (=detector cell) to the event. c Energy in your favourite units; use GeV if in doubt. c theta and phi in degrees, c both completely standard: theta measured from beam axis, etc. c Must be called between event_setup_begin and event_setup_end. c DOUBLE PRECISION energy, theta, phi INCLUDE 'ojf_kin.fh' INCLUDE 'ojf_par.fh' INCLUDE 'ojf_com.fh' IF (.NOT. ojf_event_begin) THEN WRITE(6,*) 'add_particle: 20: wrong call sequence' WRITE(6,*) ' call event_setup_begin first' STOP 'add_particle: 20' END IF IF (energy .LT. 0) THEN WRITE(6,*) 'add_particle: 21: negative energy is illegal' STOP 'add_particle: 21' END IF ojf_nparts = ojf_nparts + 1 IF (ojf_nparts .GT. nparts_max) THEN WRITE(6,*) 'add_particle: 22: max no. of particles exceeded' STOP 'add_particle: 22' END IF CALL p_set( ojf_p( 0, ojf_nparts ), energy, theta, phi ) RETURN END c c --------------------------------------------------------------------- c SUBROUTINE add_particle_raw ( px, py, pz ) c c Adds a particle (=detector cell) to the event. c This is useful with output of event generators: c px, py, pz are 3-momentum components in the same units as energy in c add_particle. pz is the beam axis. c Must be called between event_setup_begin and event_setup_end. c c Can be freely mixed with add_particle. c DOUBLE PRECISION px, py, pz INCLUDE 'ojf_kin.fh' INCLUDE 'ojf_par.fh' INCLUDE 'ojf_com.fh' IF (.NOT. ojf_event_begin) THEN WRITE(6,*) 'add_particle_raw: 20: wrong call sequence' WRITE(6,*) ' call event_setup_begin first' STOP 'add_particle_raw: 20' END IF ojf_nparts = ojf_nparts + 1 IF (ojf_nparts .GT. nparts_max) THEN WRITE(6,*) & 'add_particle_raw: 22: max no. of particles exceeded' STOP 'add_particle_raw: 22' END IF CALL p_set_raw( ojf_p( 0, ojf_nparts ), px, py, pz ) RETURN END c c ---------------------------------------------------------------------- c SUBROUTINE event_setup_end c c Must be called before jet search can be undertaken. c No particles can be added to the event afterwards. c Needed for internal housekeeping, e.g. proper normalization of c particles' energies. c INCLUDE 'ojf_kin.fh' INCLUDE 'ojf_par.fh' INCLUDE 'ojf_com.fh' INTEGER a DOUBLE PRECISION e IF (.NOT. ojf_event_begin) THEN WRITE(6,*) 'event_setup_end: 20: wrong call sequence' WRITE(6,*) ' no event setup started' STOP 'event_setup_end: 20' END IF IF (ojf_event_set) THEN WRITE(6,*) 'event_setup_end: 21: wrong call sequence' WRITE(6,*) ' event alread set' STOP 'event_setup_end: 21' END IF IF (ojf_nparts.EQ.0) THEN WRITE(6,*) 'event_setup_end: 22: no particles entered' STOP 'event_setup_end: 22' END IF IF (ojf_nparts.LT.1 .OR. ojf_nparts.GT.nparts_max) THEN WRITE(6,*) 'event_setup_end: 100' STOP 'event_setup_end: 100' END IF c c normalize energy and momenta, in three steps: c c 1) fill in ojf_e(..): c DO a = 1, ojf_nparts IF (ojf_kinematics .EQ. sphere) THEN ojf_e( a ) = ojf_p( 0, a ) ELSE IF (ojf_kinematics .EQ. cylinder) THEN ojf_e( a ) = ojf_p( par_Et, a ) END IF END DO c c 2) evaluate ojf_e_scale: c e = 0 DO a = 1, ojf_nparts e = e + ojf_e( a ) END DO IF (e .LE. 1.0E-30) THEN WRITE(6,*) & 'event_setup_end: 29: event energy too small or negative' STOP 'event_setup_end: 29' END IF ojf_e_scale = e c c 3) normalization proper: c DO a = 1, ojf_nparts ojf_p( 0, a ) = ojf_p( 0, a ) / e ojf_p( 1, a ) = ojf_p( 1, a ) / e ojf_p( 2, a ) = ojf_p( 2, a ) / e ojf_p( 3, a ) = ojf_p( 3, a ) / e IF (ojf_kinematics .EQ. cylinder) THEN ojf_p( par_Et, a ) = ojf_p( par_Et, a ) / e ojf_p( par_Eteta, a ) = ojf_p( par_Eteta, a ) / e END IF ojf_e(a) = ojf_e(a) / e END DO ojf_event_begin = .FALSE. ojf_event_set = .TRUE. RETURN END c c ====================================================================== c Setup of initial jets configuration c SUBROUTINE jets_setup_begin ( njets, Radius ) c c Must be called to begin setup of the initial jet configuration for c iterative minimization of the criterion (Omega). c c Must be called explicitly if one intends to use Q_minimize. c c Is called automatically by Q_search. c c njets is the required number of jets. Must be > 0. c Radius is the parameter R of hep-ph/9901444. c Must be > 0 and not too close to zero. c c New configurations of jets can be set up any number of times for the c same event. c c The value of seed from which the random number generator will start c for this jet configuration is stored at this point. c (From this point till the first invocation of anything random, c seed can be reset by set_seed.) c c If you need only to change Radius and proceed with minimization c starting from the current configuration, use reset_Radius. c INTEGER njets DOUBLE PRECISION Radius INCLUDE 'ojf_kin.fh' INCLUDE 'ojf_par.fh' INCLUDE 'ojf_com.fh' INTEGER a IF (ojf_jets_begin) THEN WRITE(6,*) 'jets_setup_begin: 20: wrong call sequence' WRITE(6,*) ' cannot be called twice without jets_setup_end' STOP 'jets_setup_begin: 20' END IF IF (.NOT.( (njets .GE. 1) .OR. (njets .LE. njets_max) )) THEN WRITE(6,*) 'jets_setup_begin: 21: invalid number of jets' STOP 'jets_setup_begin: 21' END IF ojf_njets = njets IF ( Radius .LE. eps_radius ) THEN WRITE(6,*) 'jets_setup_begin: 22: invalid Radius too small' STOP 'jets_setup_begin: 22' END IF ojf_Radius = Radius c ojf_Radius2 = Radius**2 / 2 ojf_Radius2 = Radius**2 ojf_seed = ojf_random c -- default value: where previous search stopped. DO a = 1, nparts_max ojf_zb_set( a ) = .FALSE. END DO ojf_jets_begin = .TRUE. ojf_jets_set = .FALSE. ojf_ran_lock = .FALSE. RETURN END c c --------------------------------------------------------------------- c SUBROUTINE jets_setup_end c c Must be called prior to minimization. Does housekeeping like c initialization of the particles whose recombination matrix elements c have not been explicitly initialized by calls like init_z_random, etc. c INCLUDE 'ojf_kin.fh' INCLUDE 'ojf_par.fh' INCLUDE 'ojf_com.fh' INTEGER a IF (.NOT. ojf_jets_begin) THEN WRITE(6,*) 'jets_setup_end: 20: wrong call sequence' WRITE(6,*) ' must be preceded by jets_setup_begin' STOP 'jets_setup_end: 20' END IF DO a = 1, ojf_nparts IF (.NOT. ojf_zb_set(a)) THEN CALL z_force_to_simplex ( ojf_z( 0, a ), ojf_b( 0, a ) ) ojf_zb_set(a) = .TRUE. ELSE CALL z_assert ( ojf_z( 0, a ), ojf_b( 0, a ) ) END IF END DO ojf_jets_begin = .FALSE. ojf_jets_set = .TRUE. CALL reset_jets RETURN END c c --------------------------------------------------------------------- c SUBROUTINE set_seed ( seed ) c c This is to allow variation of random initial configuration for the search c in case there are several local minima. c May be called once for a whole sequence of events -- each event c starts with a seed set up by the internal random numbers generator. c c seed can be remembered (see get_seed) and used as a key to regenerate c the corresponding configuration of jets (i.e. local minimum; c so the local minimum is completely determined by its seed). c c Must be called after jets_setup_begin but cannot be called c after the first invocation of init_z_random or init_random_all or c jets_setup_end and until the next jets_setup_begin. c INTEGER seed INCLUDE 'ojf_kin.fh' INCLUDE 'ojf_par.fh' INCLUDE 'ojf_com.fh' IF (ojf_ran_lock) THEN WRITE(6,*) 'set_seed: 20: wrong call sequence' WRITE(6,*) ' move this call ahead of anything random' STOP 'set_seed: 20' END IF ojf_seed = seed ojf_random = MOD( seed, random_m ) RETURN END c c ---------------------------------------------------------------------- c SUBROUTINE reset_Radius ( Radius ) c c Radius cannot be too small (see above concerning Radius). c Can be called any time -- the current configuration of jets is not c affected (only Omega etc. is recalculated properly). c c This may be useful for setting up interesting variations of the c algorithm ("annealing") in which one starts from, c say, some small value of Radius and then changes it gradually, c fine-tuning the resulting jet configurations by calls to Q_minimize. c c With infinitesimal values of Radius, the global minimum c occurs on jet configurations with the most energetic particles c playing the role of jets so this can be used to obtain c the most energetic [narrow clusters of] particles. c DOUBLE PRECISION Radius INCLUDE 'ojf_kin.fh' INCLUDE 'ojf_par.fh' INCLUDE 'ojf_com.fh' IF ( Radius .LE. eps_radius ) THEN WRITE(6,*) 'reset_Radius: 21: argument value too small' STOP 'reset_Radius: 21' END IF ojf_Radius = Radius c ojf_Radius2 = Radius**2/2 ojf_Radius2 = Radius**2 IF (ojf_jets_set) THEN CALL reset_jets END IF RETURN END c c --------------------------------------------------------------------- c SUBROUTINE set_maxiter ( maxiter ) c c Usually can be ignored -- a pro forma protection against infinite c iteration loops. c Default value is worth ~1 sec CPU time on a modest workstation. c Can be called any time. c INTEGER maxiter INCLUDE 'ojf_kin.fh' INCLUDE 'ojf_par.fh' INCLUDE 'ojf_com.fh' IF ( maxiter .LT. 1 ) THEN WRITE(6,*) 'set_maxiter: 21: invalid argument ', maxiter STOP 'set_maxiter: 21' END IF ojf_max_iter = maxiter RETURN END c c ---------------------------------------------------------------------- c SUBROUTINE set_njets_limits ( nstart, nstop ) c c Only needed in conjunction with Q_search. c Set limiting values for no. of jets in the jet search. c Must be 1 <= nstart <= nstop <= njets_max (set in ojf_par.fh). c Can be called any time. c INTEGER nstart, nstop INCLUDE 'ojf_kin.fh' INCLUDE 'ojf_par.fh' INCLUDE 'ojf_com.fh' IF (.NOT.( (nstart .GE. 1) .AND. (nstart .LE. nstop) .AND. & (nstop .LE. njets_max) )) & THEN WRITE(6,*) 'set_njets_limits: 21: invalid arguments' STOP 'set_njets_limits: 21' END IF ojf_njets_start = nstart ojf_njets_stop = nstop RETURN END c c ---------------------------------------------------------------------- c SUBROUTINE set_ntries ( n ) c c Only needed in conjunction with Q_search. c Sets no. of attempts to find minimum with c different random initial configurations for each no. of jets. c Must be n > 0. c The larger n, the higher the probability that the found c configuration is the global minimum. c Note that no. of local minima correlates positively c with no. of hard partons. c Usually values ~10 should suffice. c Can be called any time. c INTEGER n INCLUDE 'ojf_kin.fh' INCLUDE 'ojf_par.fh' INCLUDE 'ojf_com.fh' IF (.NOT.(n .GE. 1)) THEN WRITE(6,*) 'set_ntries: 21: invalid argument' STOP 'set_ntries: 21' END IF ojf_ntries = n RETURN END c c ---------------------------------------------------------------------- c SUBROUTINE set_trace_nmoved ( bool ) c c If set to TRUE will print how many particles moved at each iteration c in Q_minimize. c Can be called any time. c LOGICAL bool INCLUDE 'ojf_kin.fh' INCLUDE 'ojf_par.fh' INCLUDE 'ojf_com.fh' ojf_trace_nmoved = bool RETURN END c c ---------------------------------------------------------------------- c SUBROUTINE init_z_random_all c c The simplest way to initialize recombination matrix: c completely and uniformly random in the direct product of all the c simplices corresponding to particles. c c If only specific particles need to be randomized, use init_z_random( a ). c If only specific particles need to be set non-randomly, c call init_z( a, z_in) or assign_to_jet( a, j ) for those particles; c then call this to randomize the remaining particles. c c If this is not called explicitly, the particles not explicitly initialized c are set to "neutral" positions (democratically shared between all jets c and the soft energy). c INCLUDE 'ojf_kin.fh' INCLUDE 'ojf_par.fh' INCLUDE 'ojf_com.fh' INTEGER a IF (.NOT. ojf_jets_begin) THEN WRITE(6,*) 'init_z_random_all: 20: wrong call sequence' WRITE(6,*) ' jets setup not begun' STOP 'init_z_random_all: 20' END IF DO a = 1, ojf_nparts IF (.NOT. ojf_zb_set( a )) THEN CALL init_z_random( a ) END IF END DO RETURN END c c --------------------------------------------------------------------- c SUBROUTINE assign_to_jet ( a, j ) c c Can only be called between jets_setup_begin and jets_setup_end c if one needs to explicitly set the initial configuration of jets c -- e.g. the output of another jet algorithm to be fine-tuned. c c Directly assigns the particle a to the jet j (0 <= j <= ojf_njets; c j = 0 corresponds to soft energy), i.e. sets z(a,j) = 1, and c z(a,j')=0 for all j'.NE.j. c c This sets only *initial* configuration. No element of recombination c matrix is protected from being changed by subsequent minimizations. c INTEGER a, j INCLUDE 'ojf_kin.fh' INCLUDE 'ojf_par.fh' INCLUDE 'ojf_com.fh' INTEGER nj IF (.NOT. ojf_jets_begin) THEN WRITE(6,*) 'assign_to_jet: 20: wrong call sequence' WRITE(6,*) ' jets setup not begun' STOP 'assign_to_jet: 20' END IF IF ( (a .LT. 1) .OR. (a .GT. ojf_nparts) ) THEN WRITE(6,*) & 'assign_to_jet: 21: illegal particle specified', a STOP 'assign_to_jet: 21' END IF IF ( j .LT. 0 .OR. j .GT. ojf_njets ) THEN WRITE(6,*) & 'assign_to_jet: 22: illegal jet specified', j STOP 'assign_to_jet: 22' END IF IF ( ojf_zb_set(a) ) THEN WRITE(6,*) & 'assign_to_jet: 23: z already set for this particle', a STOP 'assign_to_jet: 23' END IF DO nj = 0, ojf_njets ojf_z( nj, a ) = 0.0 END DO ojf_z( j, a ) = 1.0 c We shouldn't forget to set borders: CALL z_snap( ojf_z( 0, a ), ojf_b( 0, a ) ) ojf_zb_set(a) = .TRUE. RETURN END c c --------------------------------------------------------------------- c SUBROUTINE init_z_from ( a, z_in ) c c Can only be called between jets_setup_begin and jets_setup_end c if one needs to explicitly set the initial configuration of jets c -- e.g. to fine-tune the output of another jet algorithm. c c Performes explicit initialization of the a-th particle's c recomb. matrix elements z(a,j) by the array z_in. c z_in must be all non-negative but need not be normalized c correctly -- correct normalization will be imposed automatically. c z_in(0) is the particle's fraction relegated to soft energy. c c For instance, z_in(j) can be some measure of distance c between this particle and the jet. c c This sets only initial configuration. No element of recombination c matrix is protected from being changed by subsequent minimizations. c INCLUDE 'ojf_kin.fh' INCLUDE 'ojf_par.fh' INCLUDE 'ojf_com.fh' INTEGER a DOUBLE PRECISION z_in( 0 : ojf_njets ) INTEGER j, jz IF (.NOT.ojf_jets_begin) THEN WRITE(6,*) 'init_z_from: 20: wrong call sequence' WRITE(6,*) ' jets setup not begun' STOP 'init_z_from: 20' END IF IF ( (a .LT. 1) .OR. (a .GT. ojf_nparts) ) THEN WRITE(6,*) & 'init_z_from: 21: illegal particle label specified', a STOP 'init_z_from: 21' END IF IF ( ojf_zb_set(a) ) THEN WRITE(6,*) & 'init_z_from: 22: z already set for this particle', a STOP 'init_z_from: 22' END IF DO j = 0, ojf_njets IF ( z_in(j) .LT. 0 ) THEN WRITE(6,*) & 'init_z_from: 23: all z_in() must be non-negative' WRITE(6,*) ' z_in array is: ', (z_in(jz), jz=0, ojf_njets) STOP 'init_z_from: 23' END IF ojf_z( j, a ) = z_in(j) END DO CALL z_force_to_simplex( ojf_z( 0, a ), ojf_b( 0, a ) ) ojf_zb_set(a) = .TRUE. RETURN END c c --------------------------------------------------------------------- c SUBROUTINE init_z_random ( a ) c c Can only be called between jets_setup_begin and jets_setup_end c if one needs to explicitly set the initial configuration of jets c -- e.g. to fine-tune the output of another jet algorithm. c c Does random initialization of recomb. matrix elements for particle a c (uniform distribution in the simplex). c INTEGER a INCLUDE 'ojf_kin.fh' INCLUDE 'ojf_par.fh' INCLUDE 'ojf_com.fh' IF (.NOT. ojf_jets_begin) THEN WRITE(6,*) 'init_z_random: 20: wrong call sequence' WRITE(6,*) ' jets setup not begun' STOP 'init_z_random: 20' END IF IF ( (a .LT. 1) .OR. (a .GT. ojf_nparts) ) THEN WRITE(6,*) & 'init_z_random: 21: illegal particle label specified', a STOP 'init_z_random: 21' END IF IF ( ojf_zb_set( a ) ) THEN WRITE(6,*) & 'init_z_random: 22: z already set for this particle', a STOP 'init_z_random: 22' END IF IF ( ojf_e( a ) .GT. zero ) THEN CALL z_set_random( ojf_z( 0, a ), ojf_b( 0, a ) ) ELSE CALL assign_to_jet( a, 0 ) END IF ojf_zb_set(a) = .TRUE. RETURN END c c ===================================================================== c Routines for getting info: c SUBROUTINE get_kinematics ( kinematics ) c c kinematics becomes sphere or cylinder, the value set c for the current event. c Cannot be called prior to the very first event_setup_begin. c INTEGER kinematics INCLUDE 'ojf_kin.fh' INCLUDE 'ojf_par.fh' INCLUDE 'ojf_com.fh' IF (.NOT.( ojf_event_begin .OR. ojf_event_set )) THEN WRITE(6,*) 'get_kinematics: 20: wrong call sequence' WRITE(6,*) & ' cannot be called prior to 1st event_setup_begin' STOP 'get_kinematics: 20' END IF kinematics = ojf_kinematics RETURN END c c --------------------------------------------------------------------- c SUBROUTINE get_nparts ( nparts, e_scale ) c c Returns the number of particles (nparts) and the total energy c of the event -- therefore, can be called only after event_setup_end. c The total energy is the usual energy for spherical kinematics, c and transverse energy for cylindrical kinematics. c c This is necessary because all energy/momentum parameters returned c by other subroutines are normalized by this value. c c Must not be called between event_setup_begin and event_setup_end. c INTEGER nparts DOUBLE PRECISION e_scale INCLUDE 'ojf_kin.fh' INCLUDE 'ojf_par.fh' INCLUDE 'ojf_com.fh' IF (.NOT. ojf_event_set) THEN WRITE(6,*) 'get_nparts: 20: wrong call sequence' WRITE(6,*) ' finish event setup first (event_setup_end)' STOP 'get_nparts: 20' END IF nparts = ojf_nparts e_scale = ojf_e_scale RETURN END c c --------------------------------------------------------------------- c SUBROUTINE get_particle ( a, e, xta, phi, p, ephys, pphys ) c c Takes particle a (1 <= a <= ojf_nparts) and returns: c for spherical kinematics: c e = normalized energy, xta = theta, phi c for cylindrical kinematics: c e = normalized transverse energy, xta = rapidity, phi c c p becomes particle's normalized 4-momenum c [p(3) <-> pz <-> beam axis; metrics is Bjorken-Drell]. c c All angles are in degrees. c c Remember that energies, momenta, etc., are always normalized on event setup c so that total energy (transverse for cylindrical kinematics) c of the event is 1. c c ephys and pphys are the same as e and p but in physical units c (the same as used on input). c c Must not be called between event_setup_begin and event_setup_end. c INTEGER a DOUBLE PRECISION e, xta, phi, p(0:3), ephys, pphys(0:3) INCLUDE 'ojf_kin.fh' INCLUDE 'ojf_par.fh' INCLUDE 'ojf_com.fh' IF (.NOT. ojf_event_set) THEN WRITE(6,*) 'get_particle: 21: event not set' STOP 'get_particle: 21' END IF IF ( (a .LT. 1) .OR. (a .GT. ojf_nparts) ) THEN WRITE(6,*) & 'get_particle: 22: illegal particle label specified', a STOP 'get_particle: 22' END IF CALL p_get( e, xta, phi, p, ojf_p( 0, a ) ) ephys = e * ojf_e_scale pphys(0) = p(0) * ojf_e_scale pphys(1) = p(1) * ojf_e_scale pphys(2) = p(2) * ojf_e_scale pphys(3) = p(3) * ojf_e_scale RETURN END c c --------------------------------------------------------------------- c SUBROUTINE get_njets ( njets ) c c Returns the number of jets in the current configuration of jets. c Cannot be called before the first configuration of jets is setup. c INTEGER njets INCLUDE 'ojf_kin.fh' INCLUDE 'ojf_par.fh' INCLUDE 'ojf_com.fh' IF (.NOT.( ojf_jets_set .OR. ojf_jets_begin )) THEN WRITE(6,*) 'get_njets: 20: wrong call sequence' WRITE(6,*) ' first setup a jet configuration' STOP 'get_njets: 20' END IF njets = ojf_njets RETURN END c c --------------------------------------------------------------------- c SUBROUTINE get_seed ( seed ) c c Returns the seed for random generator used for setting up the current c jet configuration. c The value of seed is 'locked' (causing attempts to reset it to result c in STOP) by the first invocation of anything 'random' and retained c until 'unlocked' and reset by jets_setup_begin. c INTEGER seed INCLUDE 'ojf_kin.fh' INCLUDE 'ojf_par.fh' INCLUDE 'ojf_com.fh' seed = ojf_seed RETURN END c c --------------------------------------------------------------------- c SUBROUTINE get_Radius ( R ) DOUBLE PRECISION R INCLUDE 'ojf_kin.fh' INCLUDE 'ojf_par.fh' INCLUDE 'ojf_com.fh' R = ojf_Radius RETURN END c c --------------------------------------------------------------------- c SUBROUTINE get_maxiter ( maxiter ) INTEGER maxiter INCLUDE 'ojf_kin.fh' INCLUDE 'ojf_par.fh' INCLUDE 'ojf_com.fh' maxiter = ojf_max_iter RETURN END c c --------------------------------------------------------------------- c SUBROUTINE get_njets_limits ( nstart, nstop ) INTEGER nstart, nstop INCLUDE 'ojf_kin.fh' INCLUDE 'ojf_par.fh' INCLUDE 'ojf_com.fh' nstart = ojf_njets_start nstop = ojf_njets_stop RETURN END c c ---------------------------------------------------------------------- c SUBROUTINE get_ntries ( n ) INTEGER n INCLUDE 'ojf_kin.fh' INCLUDE 'ojf_par.fh' INCLUDE 'ojf_com.fh' n = ojf_ntries RETURN END c c ---------------------------------------------------------------------- c SUBROUTINE get_criterion ( omega, y, esoft ) c c Whenever a jet configuration is set up or modified, c the corresponding values of the criterion (Omega) and its c components Y and E_soft are recalculated and can be retrieved using c this. c INCLUDE 'ojf_kin.fh' INCLUDE 'ojf_par.fh' INCLUDE 'ojf_com.fh' DOUBLE PRECISION omega, y, esoft IF (.NOT. ojf_jets_set) THEN WRITE(6,*) 'get_criterion: 20: wrong call sequence' WRITE(6,*) ' setup a jet configuration first' STOP 'get_criterion: 20' END IF omega = ojf_Omega y = ojf_Y esoft = ojf_Esoft RETURN END c c --------------------------------------------------------------------- c SUBROUTINE get_jet ( j, e, xta, phi, q, qtilde, ephys, qphys ) c c Takes jet j (1 <= j <= ojf_njets) and returns: c for spherical kinematics: c e = normalized energy, xta = theta, phi c for cylindrical kinematics: c e = normalized transverse energy, xta = rapidity, phi c c q becomes jet's physical normalized 4-momenum (as defined in hep-ph/9901444). c qtilde becomes the light-like 4-vector denoted in the paper as $\tilde q_j$. c c All angles are in degrees. c c Remember that energies, momenta, etc., are always normalized on event setup c so that total energy (or total transverse energy in cylindrical kinematics) c of the event is 1, internally. c c ephys and qphys are the same as e and q but in physical units c (the same as used on input). c INTEGER j DOUBLE PRECISION e, xta, phi, q(0:3), qtilde(0:3) DOUBLE PRECISION ephys, qphys(0:3) INCLUDE 'ojf_kin.fh' INCLUDE 'ojf_par.fh' INCLUDE 'ojf_com.fh' INTEGER k IF (.NOT. ojf_jets_set) THEN WRITE(6,*) 'get_jet: 20: illegal call sequence' WRITE(6,*) ' jets not set' STOP 'get_jet: 20' END IF IF ( (j .LT. 1) .OR. (j .GT. ojf_njets) ) THEN WRITE(6,*) 'get_jet: 22: illegal jet tag specified', j STOP 'get_jet: 22' END IF CALL p_get( e, xta, phi, q, ojf_q( 0, j ) ) DO k = 0, 3 qtilde( k ) = ojf_q( k + tilde, j ) END DO ephys = e * ojf_e_scale qphys(0) = q(0) * ojf_e_scale qphys(1) = q(1) * ojf_e_scale qphys(2) = q(2) * ojf_e_scale qphys(3) = q(3) * ojf_e_scale RETURN END c c --------------------------------------------------------------------- c SUBROUTINE get_z ( a, z_out ) INTEGER a INCLUDE 'ojf_kin.fh' INCLUDE 'ojf_par.fh' INCLUDE 'ojf_com.fh' DOUBLE PRECISION z_out( 0 : ojf_njets ) c used to be ojf_nparts; corrected 2000-03-30 --FT c c Takes particle a (1 <= a <= ojf_nparts) and returns array of c coefficients z(0 : ojf_njets), c where z_out(j) is the particle's contribution to jet j. c j = 0 corresponds to soft energy. c Note: SUM( z_out ) = 1 c INTEGER j IF ( (a .LT. 1) .OR. (a .GT. ojf_nparts) ) THEN WRITE(6,*) 'get_z: 21: illegal particle label specified' STOP 'get_z: 21' END IF DO j = 0, ojf_njets z_out(j) = ojf_z( j, a ) END DO RETURN END c c --------------------------------------------------------------------- c SUBROUTINE get_particle_split ( a, total_jets, jet, zj ) c c Takes particle a (1 <= a <= ojf_nparts) and returns: c c total_jets -- no. of jets (including soft energy, j=0) which include c a non-zero fraction of the particle [i.e. have nonzero z(j,a)] c jet ( 0 : total_jets-1 ) -- their labels j; c zj ( 0 : total_jets-1 ) -- the corresponding z(a,j) c c The ordering of the two arrays is such that zj(j) >= zj(j+1). c INTEGER a, total_jets INCLUDE 'ojf_kin.fh' INCLUDE 'ojf_par.fh' INCLUDE 'ojf_com.fh' INTEGER jet ( 0 : ojf_njets ) DOUBLE PRECISION zj ( 0 : ojf_njets ) c used to be ojf_nparts; corrected 2000-03-30 --FT c INTEGER j DOUBLE PRECISION zz IF ((a .LT. 1) .OR. (a .GT. ojf_nparts)) THEN WRITE(6,*) & 'get_particle_split: 21: illegal particle label', a STOP 'get_particle_split: 21' END IF c c initialization: c total_jets = 0 DO j = 0, ojf_njets zz = ojf_z( j, a ) IF (zz .EQ. one) THEN IF (total_jets .NE. 0) THEN WRITE(6,*) & 'get_particle_split: 100: internal error' STOP 'get_particle_split: 100' END IF total_jets = 1 jet (0) = j zj (0) = 1 RETURN ELSE IF (zz .NE. zero) THEN jet ( total_jets ) = j zj ( total_jets ) = zz total_jets = total_jets + 1 END IF END DO CALL aux_sort( zj(0), jet(0), total_jets ) RETURN END c c --------------------------------------------------------------------- c SUBROUTINE get_jet_split ( j, nwhole, whole_a, & nfract, fract_a, fract_z ) INCLUDE 'ojf_kin.fh' INCLUDE 'ojf_par.fh' INCLUDE 'ojf_com.fh' INTEGER j, nwhole, nfract INTEGER whole_a( 1 : ojf_nparts ), fract_a( 1 : ojf_nparts ) DOUBLE PRECISION fract_z( 1 : ojf_nparts ) c c Takes jet j (0 <= j <= ojf_njets; j=0 is the soft energy) and returns: c c nwhole -- number of particles wholly in the jet [ z(j,a)=1 ]; c whole_a -- labels of such particles; c nfract -- number of particles partially in the jet [ 0 < z(j,a) < 1 ]; c fract_a -- tags of such particles; c fract_z -- the corresponding z(j,a) for such particles. c The latter two arrays are synchronously ordered so that fract_z(1) is largest. c INTEGER a DOUBLE PRECISION z IF ((j .LT. 0) .OR. (j .GT. ojf_njets)) THEN WRITE(6,*) & 'get_jet_split: 21: illegal jet label specified', j STOP 'get_jet_split: 21' END IF nwhole = 0 nfract = 0 DO a = 1, ojf_nparts z = ojf_z( j, a ) IF (z .EQ. zero) THEN c do nothing ELSE IF (z .EQ. one) THEN nwhole = nwhole + 1 whole_a( nwhole ) = a ELSE nfract = nfract + 1 fract_a( nfract ) = a fract_z( nfract ) = z END IF END DO CALL aux_sort( fract_z, fract_a, nfract ) RETURN END c c ------------------------------------------------------------------- c SUBROUTINE sort_jets c c Sort jets 1<=j<=ojf_njets so that first jet (j=1) had largest energy c [transverse energy for cylindrical geometry]. c All internal data are reordered correspondingly. c c The sort does not include j=0 (soft energy). c c This sort only affects jets' internal labels, not their physical content, c and one can proceed with things like changing Radius and applying c Q_minimize, etc. c INCLUDE 'ojf_kin.fh' INCLUDE 'ojf_par.fh' INCLUDE 'ojf_com.fh' DOUBLE PRECISION jet_e ( 1 : njets_max ) INTEGER order ( 1 : njets_max ) DOUBLE PRECISION buf_z ( 0 : njets_max ) LOGICAL buf_b ( 0 : njets_max ) INTEGER j, k, a IF (.NOT. ojf_jets_set) THEN WRITE(6,*) 'sort_jets: 20: wrong call sequence' STOP 'sort_jets: 20' END IF c c initialize jet_e and order: c DO j = 1, ojf_njets order(j) = j IF (ojf_kinematics .EQ. sphere) THEN jet_e(j) = ojf_q( 0, j ) ELSE IF (ojf_kinematics .EQ. cylinder) THEN jet_e(j) = ojf_q( par_Et, j ) END IF END DO c c sort jet_e and order in parallel: c CALL aux_sort( jet_e, order, ojf_njets ) c c reorder ojf_q, ojf_z, ojf_b according to order: c c (Dima, prosti. U menia nadpisj na monitore: "Keep It Simple, Stupid!" Poetomu.) c DO a = 1, ojf_nparts DO j = 1, ojf_njets buf_z( j ) = ojf_z( j, a ) END DO DO j = 1, ojf_njets ojf_z( j, a ) = buf_z( order( j ) ) END DO DO j = 1, ojf_njets buf_b( j ) = ojf_b( j, a ) END DO DO j = 1, ojf_njets ojf_b( j, a ) = buf_b( order( j ) ) END DO END DO DO k = 0, 12 DO j = 1, ojf_njets buf_z( j ) = ojf_q( k, j ) END DO DO j = 1, ojf_njets ojf_q( k, j ) = buf_z( order( j ) ) END DO END DO RETURN END c c end of interface modules --------------------------------------------- c c ====================================================================== c ====================================================================== c ====================================================================== c c Part 3. THE CORE STUFF c ************************************************************************ ************************************************************************ ************************************************************************ *** *** *** Like Shurik is fond of saying, THINK THRICE BEFORE... *** *** *** ************************************************************************ ************************************************************************ ************************************************************************ c BLOCK DATA ojf_lock INCLUDE 'ojf_kin.fh' INCLUDE 'ojf_par.fh' INCLUDE 'ojf_com.fh' DATA ojf_Radius/1.0/, ojf_ntries/10/ DATA ojf_max_iter/1000/, ojf_random/13/ DATA ojf_njets_start/1/, ojf_njets_stop/njets_max/ c DATA ojf_event_begin/.FALSE./, ojf_event_set/.FALSE./ DATA ojf_jets_begin/.FALSE./, ojf_jets_set/.FALSE./ DATA ojf_trace_nmoved/.FALSE./ DATA ojf_zb_set/nparts_max*.FALSE./, ojf_ran_lock/.FALSE./ END c c ---------------------------------------------------------------------- c SUBROUTINE Q_minimize ( success ) c c Minimizes Omega for a given no. of jets (ojf_njets) c and starting from the existing configuration of z. c c This is the only user callable subroutine here in Part 3 of the code. c LOGICAL success INCLUDE 'ojf_kin.fh' INCLUDE 'ojf_par.fh' INCLUDE 'ojf_com.fh' EXTERNAL ojf_lock INTEGER a, moved LOGICAL this_moved IF (.NOT. ojf_jets_set) THEN WRITE(6,*) 'Q_minimize: 20: wrong call sequence' WRITE(6,*) ' first set up a jet configuration' STOP 'Q_minimize: 20' END IF ojf_att_moves = 0 ojf_undone_moves = 0 ojf_a_old = -1 cccccccccccccccc IF (ojf_trace_nmoved) THEN WRITE(6,*) & 'tracing Q_minimize: no. of particles moved on each iteration' END IF cccccccccccccccc success = .FALSE. DO ojf_iter = 1, ojf_max_iter moved = 0 DO a = 1, ojf_nparts CALL Q_minimize_wrt( a, this_moved ) IF (this_moved) THEN moved = moved + 1 END IF END DO CALL reset_jets cccccccccccccccc IF (ojf_trace_nmoved) THEN WRITE(6,*) 'on iter =', ojf_iter, ' moved =', moved END IF cccccccccccccccc IF ( moved .EQ. 0 ) THEN success = .TRUE. GOTO 13 END IF END DO 13 CONTINUE RETURN END c c --------------------------------------------------------------------- c --------------------------------------------------------------------- c c NONE of the following routines are supposed to be called by the user. c c --------------------------------------------------------------------- c --------------------------------------------------------------------- c SUBROUTINE Q_minimize_wrt ( a, moved ) INTEGER a LOGICAL moved INCLUDE 'ojf_kin.fh' INCLUDE 'ojf_par.fh' INCLUDE 'ojf_com.fh' DOUBLE PRECISION d_norm EXTERNAL d_norm DOUBLE PRECISION dist, step, d( 0 : njets_max ), omega_old IF (.NOT. ojf_jets_set) THEN WRITE(6,*) 'Q_minimize_wrt: 30: wrong call sequence' STOP 'Q_minimize_wrt: 30' END IF CALL grad_Omega ( d, a ) CALL d_minus_snap( d, ojf_b( 0, a ) ) CALL d_eval_step ( d, ojf_z( 0, a ), ojf_b( 0, a ), step ) moved = .FALSE. IF (( step .GT. 0 ) .AND. (step .LT. inf_step )) THEN dist = step * d_norm( d ) DO WHILE ( dist .GT. eps_dist ) omega_old = ojf_Omega CALL move_particle( a, d, step ) moved = ( ojf_Omega .LT. omega_old ) IF (moved) GOTO 111 CALL undo_move_particle dist = dist * 0.3666667 step = step * 0.3666667 c Dima likes 0.4, Fyodor likes 1/3 c Dima + Fyodor = two scientific programmers END DO END IF 111 CONTINUE RETURN END c c --------------------------------------------------------------------- c SUBROUTINE move_particle ( a, d, step ) INCLUDE 'ojf_kin.fh' INCLUDE 'ojf_par.fh' INCLUDE 'ojf_com.fh' INTEGER a DOUBLE PRECISION d( 0 : ojf_njets ), step INTEGER j, k ojf_att_moves = ojf_att_moves + 1 c c save old configuration to make possible update_jets and undo: c DO j = 0, ojf_njets ojf_z_old(j) = ojf_z( j, a ) ojf_b_old(j) = ojf_b( j, a ) END DO DO j = 1, ojf_njets DO k = 0, 12 ojf_q_old( k, j ) = ojf_q( k, j ) END DO END DO ojf_a_old = a ojf_Esoft_old = ojf_Esoft ojf_Y_old = ojf_Y ojf_Omega_old = ojf_Omega CALL z_move_by( ojf_z( 0, a ), d, step, ojf_b( 0, a ) ) CALL update_jets( a ) RETURN END c c ---------------------------------------------------------------------- c SUBROUTINE update_jets ( a ) INTEGER a c c incremental reevaluation of jets' parameters ------------------------ c INCLUDE 'ojf_kin.fh' INCLUDE 'ojf_par.fh' INCLUDE 'ojf_com.fh' INTEGER j DOUBLE PRECISION z IF (.NOT. ojf_jets_set) THEN WRITE(6,*) 'update_jets: 30: wrong call sequence' STOP 'update_jets: 30' END IF c c update soft energy c z = ojf_z( 0, a ) - ojf_z_old( 0 ) IF (z.NE.0) THEN ojf_Esoft = ojf_Esoft + ojf_e( a ) * z END IF c c update Y c DO j = 1, ojf_njets z = ojf_z( j, a ) - ojf_z_old(j) IF (z.NE.0) THEN CALL j_add_mult( ojf_q( 0, j ), z, ojf_p( 0, a ) ) CALL j_eval_nonlinear( ojf_q( 0, j ) ) END IF END DO CALL eval_Y c c LASTLY, update Omega c CALL eval_Omega RETURN END c c ---------------------------------------------------------------------- c SUBROUTINE undo_move_particle INCLUDE 'ojf_kin.fh' INCLUDE 'ojf_par.fh' INCLUDE 'ojf_com.fh' INTEGER j, k ojf_undone_moves = ojf_undone_moves + 1 IF (.NOT.(ojf_a_old .GE. 0)) THEN WRITE(6,*) 'undo_move_particle: 30: wrong call sequence' STOP 'undo_move_particle: 30' END IF DO j = 0, ojf_njets ojf_z( j, ojf_a_old ) = ojf_z_old(j) ojf_b( j, ojf_a_old ) = ojf_b_old(j) END DO DO j = 1, ojf_njets DO k = 0, 12 ojf_q( k, j ) = ojf_q_old( k, j ) END DO END DO ojf_a_old = -1 ojf_Esoft = ojf_Esoft_old ojf_Y = ojf_Y_old ojf_Omega = ojf_Omega_old ojf_Omega_old = inf RETURN END c c --------------------------------------------------------------------- c SUBROUTINE reset_jets INCLUDE 'ojf_kin.fh' INCLUDE 'ojf_par.fh' INCLUDE 'ojf_com.fh' INTEGER a, j DOUBLE PRECISION z CALL eval_Esoft DO j = 1, ojf_njets CALL j_set_to_zero( ojf_q( 0, j ) ) DO a = 1, ojf_nparts z = ojf_z( j, a ) IF (z.NE.0) THEN CALL j_add_mult( ojf_q( 0, j ), z, ojf_p( 0, a ) ) END IF END DO CALL j_eval_nonlinear( ojf_q( 0, j ) ) END DO CALL eval_Y CALL eval_Omega RETURN END c c ====================================================================== c Criterion and its derivatives c SUBROUTINE eval_Omega INCLUDE 'ojf_kin.fh' INCLUDE 'ojf_par.fh' INCLUDE 'ojf_com.fh' IF (.NOT. ojf_jets_set) THEN WRITE(6,*) 'eval_Omega: 30: set up jets first' STOP 'eval_Omega: 30' END IF ojf_Omega = ojf_Y / ojf_Radius2 + ojf_Esoft RETURN END c c ---------------------------------------------------------------------- c SUBROUTINE eval_Y INCLUDE 'ojf_kin.fh' INCLUDE 'ojf_par.fh' INCLUDE 'ojf_com.fh' INTEGER j DOUBLE PRECISION s IF (.NOT. ojf_jets_set) THEN WRITE(6,*) 'eval_Y: 30: set up jets first' STOP 'eval_Y: 30' END IF s = 0 DO j = 1, ojf_njets s = s + ojf_q( par_Y, j ) END DO c ojf_Y = s ojf_Y = s * 2 RETURN END c c ---------------------------------------------------------------------- c SUBROUTINE eval_Esoft INCLUDE 'ojf_kin.fh' INCLUDE 'ojf_par.fh' INCLUDE 'ojf_com.fh' INTEGER a DOUBLE PRECISION s IF (.NOT. ojf_jets_set) THEN WRITE(6,*) 'eval_Esoft: 30: set up jets first' STOP 'eval_Esoft: 30' END IF s = 0 DO a = 1, ojf_nparts s = s + ojf_e( a ) * ojf_z( 0, a ) END DO IF ( .NOT.( s.GE.0 ) ) THEN WRITE(6,*) 'eval_Esoft: 60: invariant violated' STOP 'eval_Esoft: 60' END IF ojf_Esoft = s RETURN END c c ---------------------------------------------------------------------- c SUBROUTINE grad_Omega ( d, a ) INCLUDE 'ojf_kin.fh' INCLUDE 'ojf_par.fh' INCLUDE 'ojf_com.fh' DOUBLE PRECISION d( 0 : ojf_njets ), x INTEGER a, j IF (.NOT. ojf_jets_set) THEN WRITE(6,*) 'grad_Omega: 30: set up jets first' STOP 'grad_Omega: 30' END IF CALL grad_Y( d, a ) x = ojf_e( a ) DO j = 1, ojf_njets d(j) = d(j) / ojf_Radius2 - x END DO RETURN END c c ---------------------------------------------------------------------- c SUBROUTINE grad_Y ( d, a ) INCLUDE 'ojf_kin.fh' INCLUDE 'ojf_par.fh' INCLUDE 'ojf_com.fh' DOUBLE PRECISION d( 0 : ojf_njets ) INTEGER a, j EXTERNAL pos_prod DOUBLE PRECISION pos_prod DOUBLE PRECISION r024, res, x IF (.NOT. ojf_jets_set) THEN WRITE(6,*) 'grad_Y: 30: set up jets first' STOP 'grad_Y: 30' END IF d(0) = 0 DO j = 1, ojf_njets res = pos_prod( ojf_q( tilde, j ), ojf_p( 0, a ) ) IF (ojf_kinematics .EQ. cylinder) THEN x = MAX( ojf_q( par_Et, j ), eps_Et ) r024 = ojf_p( par_Et, a ) / x & * ( ojf_p( par_eta, a ) - ojf_q( par_eta, j ) ) res = res + r024 * ojf_q( p0shmpzch, j ) END IF c d(j) = res d(j) = res * 2 END DO RETURN END c c ====================================================================== c Doing things in a simplex c SUBROUTINE d_minus_snap ( d, b ) c c This is the most tricky subroutine in ojf. c d is supposed to be gradient [with d(0) undefined] c INCLUDE 'ojf_kin.fh' INCLUDE 'ojf_par.fh' INCLUDE 'ojf_com.fh' DOUBLE PRECISION d( 0 : ojf_njets ) LOGICAL b( 0 : ojf_njets ) EXTERNAL d_norm DOUBLE PRECISION d_norm INTEGER j, j_large LOGICAL b_large( 0 : njets_max ), bool DOUBLE PRECISION maxnorm, newnorm, s, & f( 0 : njets_max ), new( 0 : njets_max ) c c the minus part is done here: c f(0) = zero DO j = 1, ojf_njets f(j) = - d(j) END DO c c b_large models a set; initialize it: c DO j = 0, ojf_njets b_large(j) = b(j) END DO c c check: there are elements outside b_large c (indeed, a point cannot be on all boundaries of the simplex) c bool = .TRUE. DO j = 0, ojf_njets bool = bool .AND. b_large(j) END DO IF (bool) THEN WRITE(6,*) 'd_minus_snap: 100' STOP 'd_minus_snap: 100' END IF c c try all elements outside b_large and find which of them c yields what is presumably direction of fastest decrease c maxnorm = - 1 111 CONTINUE c find next j_large not belonging to b_large: j_large = - 1 222 CONTINUE j_large = j_large + 1 IF ( b_large( j_large ) ) GOTO 222 s = zero DO j = 0, ojf_njets new(j) = f(j) - f( j_large ) IF ( b(j) ) THEN new(j) = MAX( zero, new(j) ) END IF IF ( j .NE. j_large ) s = s + new(j) END DO new( j_large ) = - s c c check correctness of new: s = zero DO j = 0, ojf_njets s = s + new(j) END DO IF (.NOT.( DABS(s) .LT. eps_sum )) THEN WRITE(6,*) & 'W A R N I N G: d_minus_snap: 101: large rounding errors' WRITE(6,*) & 's must be zero but actually s = ', s END IF s = 0 DO j = 1, ojf_njets s = s + f(j) * new(j) END DO IF (.NOT.( s .GT. -eps_round )) THEN WRITE(6,*) & 'W A R N I N G: d_minus_snap: 102: large rounding errors' WRITE(6,*) & 's must be zero but actually s = ', s END IF c end of checks. c newnorm = d_norm( new ) IF ( newnorm .GT. maxnorm ) THEN maxnorm = newnorm DO j = 0, ojf_njets d(j) = new(j) END DO END IF b_large( j_large ) = .TRUE. bool = .TRUE. DO j = 0, ojf_njets bool = bool .AND. b_large(j) END DO IF (.NOT.bool) GOTO 111 RETURN END c c ---------------------------------------------------------------------- c SUBROUTINE d_eval_step ( d, z, b, step ) INCLUDE 'ojf_kin.fh' INCLUDE 'ojf_par.fh' INCLUDE 'ojf_com.fh' DOUBLE PRECISION d( 0 : ojf_njets ), z( 0 : ojf_njets ) LOGICAL b( 0 : ojf_njets ) DOUBLE PRECISION step DOUBLE PRECISION t INTEGER j step = inf_step DO j = 0, ojf_njets IF ( b(j) ) THEN IF (.NOT.( d(j) .GE. 0 )) THEN WRITE(6,*) 'd_eval_step: 31' STOP 'd_eval_step: 31' END IF t = inf_step ELSE IF ( d(j) .EQ. 0 ) THEN t = inf_step ELSE t = - z(j) / d(j) IF ( t .LT. 0 ) t = inf_step END IF END IF step = MIN( step, t ) END DO RETURN END c c ---------------------------------------------------------------------- c SUBROUTINE d_assert ( d, z, b ) INCLUDE 'ojf_kin.fh' INCLUDE 'ojf_par.fh' INCLUDE 'ojf_com.fh' DOUBLE PRECISION d( 0 : ojf_njets ), z( 0 : ojf_njets ) LOGICAL b( 0 : ojf_njets ) INTEGER j DOUBLE PRECISION s, dj CALL z_assert( z, b ) s = 0 DO j = 0, ojf_njets dj = d(j) IF ( b(j) .AND. (dj .LT. 0) ) THEN WRITE(6,*) 'd_assert: 31' STOP 'd_assert: 31' END IF s = s + dj END DO IF ( .NOT.(DABS(s) .LT. eps_sum) ) THEN WRITE(6,*) 'd_assert: 32' STOP 'd_assert: 32' END IF RETURN END c c ---------------------------------------------------------------------- c DOUBLE PRECISION FUNCTION d_norm ( d ) INCLUDE 'ojf_kin.fh' INCLUDE 'ojf_par.fh' INCLUDE 'ojf_com.fh' DOUBLE PRECISION d( 0 : ojf_njets ) INTEGER j DOUBLE PRECISION r r = 0 DO j = 0, ojf_njets r = MAX( r, DABS( d(j) ) ) END DO d_norm = r RETURN END c c ---------------------------------------------------------------------- c SUBROUTINE z_move_by ( z, d, step, b ) INCLUDE 'ojf_kin.fh' INCLUDE 'ojf_par.fh' INCLUDE 'ojf_com.fh' DOUBLE PRECISION z( 0 : ojf_njets ), d( 0 : ojf_njets ) LOGICAL b( 0 : ojf_njets ) DOUBLE PRECISION step INTEGER j CALL d_assert( d, z, b ) IF (.NOT.( ( step .GT. 0 ) . AND. ( step .LT. inf_step ) )) THEN WRITE(6,*) 'z_move_by: 31' STOP 'z_move_by: 31' END IF DO j = 0, ojf_njets z(j) = z(j) + step * d(j) END DO CALL z_snap( z, b ) RETURN END c c ---------------------------------------------------------------------- c SUBROUTINE z_set_random ( z, b ) INCLUDE 'ojf_kin.fh' INCLUDE 'ojf_par.fh' INCLUDE 'ojf_com.fh' DOUBLE PRECISION z( 0 : ojf_njets ) LOGICAL b( 0 : ojf_njets ) EXTERNAL random DOUBLE PRECISION random c c random point uniformly distributed in the simplex c INTEGER j, k DOUBLE PRECISION zj c c initially we need only ojf_njets random numbers c DO j = 0, ojf_njets - 1 z(j) = random() END DO z( ojf_njets ) = 1 c c sorting (ojf_njets is never large enough to justify fancy algorithms) c DO 10 j = ojf_njets - 2, 0, -1 zj = z(j) DO k = j + 1, ojf_njets - 1 IF ( z( k ) .LT. zj ) THEN z( k-1 ) = z( k ) ELSE z( k-1 ) = zj GOTO 10 END IF END DO z( ojf_njets-1 ) = zj 10 CONTINUE c c evaluating differences c DO j = ojf_njets, 1, -1 z(j) = z(j) - z( j-1 ) IF (.NOT.( z(j) .GE. 0 )) THEN WRITE(6,*) 'z_set_random: 60' STOP 'z_set_random: 60' END IF END DO CALL z_snap( z, b ) RETURN END c c ---------------------------------------------------------------------- c SUBROUTINE z_force_to_simplex ( z, b ) INCLUDE 'ojf_kin.fh' INCLUDE 'ojf_par.fh' INCLUDE 'ojf_com.fh' DOUBLE PRECISION z( 0 : ojf_njets ) LOGICAL b( 0 : ojf_njets ) INTEGER j DOUBLE PRECISION s s = 0 DO j = 0, ojf_njets z(j) = MAX( zero, z(j) ) s = s + z(j) END DO IF ( s .LE. eps_sum0 ) THEN s = one / ( ojf_njets + one ) DO j = 0, ojf_njets z(j) = s END DO ELSE DO j = 0, ojf_njets z(j) = z(j) / s END DO END IF CALL z_snap( z, b ) RETURN END c c ---------------------------------------------------------------------- c SUBROUTINE z_snap ( z, b ) INCLUDE 'ojf_kin.fh' INCLUDE 'ojf_par.fh' INCLUDE 'ojf_com.fh' DOUBLE PRECISION z( 0 : ojf_njets ) LOGICAL b( 0 : ojf_njets ) INTEGER j, nzeros, vertex DOUBLE PRECISION s, zj s = 0 DO j = 0, ojf_njets zj = z(j) IF (.NOT.( ( zj .GT. -eps_round ) .AND. & ( zj .LT. (1 + eps_round) ) )) THEN WRITE(6,*) 'z_snap: 31' STOP 'z_snap: 31' END IF s = s + zj END DO IF (.NOT.( DABS(s - 1) .LT. eps_sum1 )) & THEN WRITE(6,*) 'z_snap: 32' STOP 'z_snap: 32' END IF s = 0 nzeros = 0 DO j = 0, ojf_njets b(j) = .FALSE. END DO DO j = 0, ojf_njets zj = z(j) IF ( zj .LT. eps_snap ) THEN zj = 0 nzeros = nzeros + 1 b(j) = .TRUE. END IF z(j) = zj s = s + zj END DO IF (nzeros .EQ. ojf_njets) THEN j = 0 DO WHILE (b(j)) j = j + 1 END DO z(j) = 1 vertex = j ELSE vertex = - ( ojf_njets + 1 - nzeros ) s = 1 / s DO j = 0, ojf_njets z(j) = z(j) * s END DO END IF c vertex: c if >= 0 then vertex = jet to which the particle belongs wholly c else vertex = minus no. of jets which share the particle CALL z_assert( z, b ) RETURN END c c ---------------------------------------------------------------------- c SUBROUTINE z_assert ( z, b ) INCLUDE 'ojf_kin.fh' INCLUDE 'ojf_par.fh' INCLUDE 'ojf_com.fh' DOUBLE PRECISION z( 0 : ojf_njets ) LOGICAL b( 0 : ojf_njets ) INTEGER j DOUBLE PRECISION s, zj s = 0 DO j = 0, ojf_njets zj = z(j) IF (.NOT.( (zj .GE. eps_snap) .OR. & ( (zj .EQ. zero) .AND. b(j) ) )) THEN WRITE(6,*) 'z_assert: 31' STOP 'z_assert: 31' END IF s = s + zj END DO IF (.NOT.( DABS(s - 1) .LT. eps_sum )) THEN WRITE(6,*) 'z_assert: 32' STOP 'z_assert: 32' END IF RETURN END c c ====================================================================== c Operations on jets c SUBROUTINE j_eval_nonlinear ( qj ) DOUBLE PRECISION qj(0:12) INCLUDE 'ojf_kin.fh' INCLUDE 'ojf_par.fh' INCLUDE 'ojf_com.fh' EXTERNAL pos_prod, random DOUBLE PRECISION pos_prod, random DOUBLE PRECISION norm, x, y, z IF (ojf_kinematics .EQ. sphere) THEN qj( tilde ) = 1 norm = DSQRT( qj(1)**2 + qj(2)**2 + qj(3)**2 ) IF (norm .LT. eps_norm) THEN x = random() y = random() z = random() norm = DSQRT( x*x + y*y + z*z ) qj( tilde + 1 ) = x / norm qj( tilde + 2 ) = y / norm qj( tilde + 3 ) = z / norm c -- not quite uniform but will do. ELSE qj( tilde + 1 ) = qj( 1 ) / norm qj( tilde + 2 ) = qj( 2 ) / norm qj( tilde + 3 ) = qj( 3 ) / norm END IF ELSE IF (ojf_kinematics .EQ. cylinder) THEN IF ( ABS( qj( par_Et ) ) .LT. eps_Et ) THEN qj( par_eta ) = ( random() - 0.5 ) * 6 c -- random uniform between -3, +3; will do. ELSE qj( par_eta ) = qj( par_Eteta ) / qj( par_Et ) END IF qj( tilde + 0 ) = DCOSH( qj( par_eta ) ) qj( tilde + 3 ) = DSINH( qj( par_eta ) ) norm = DSQRT( qj(1)**2 + qj(2)**2 ) IF (norm .LT. eps_norm) THEN x = random() y = random() norm = DSQRT( x*x + y*y ) qj( tilde + 1 ) = x / norm qj( tilde + 2 ) = y / norm c -- not quite uniform but will do. ELSE qj( tilde + 1 ) = qj(1) / norm qj( tilde + 2 ) = qj(2) / norm END IF qj( p0shmpzch ) = qj(0) * qj( tilde + 3 ) & - qj(3) * qj( tilde + 0 ) END IF qj( par_Y ) = pos_prod( qj(0), qj( tilde ) ) RETURN END c c ---------------------------------------------------------------------- c SUBROUTINE j_add_mult ( qj, z, pa ) DOUBLE PRECISION qj(0:12), z, pa(0:6) c c accumulate contribution of particle pa to jet qj with weight z c INCLUDE 'ojf_kin.fh' INCLUDE 'ojf_par.fh' INCLUDE 'ojf_com.fh' CALL v_add_mult( qj, z, pa ) IF (ojf_kinematics .EQ. cylinder) THEN qj( par_Et ) = qj( par_Et ) + z * pa( par_Et ) qj( par_Eteta ) = qj( par_Eteta ) + z * pa( par_Eteta ) END IF RETURN END c c ---------------------------------------------------------------------- c SUBROUTINE j_set_to_zero ( qj ) DOUBLE PRECISION qj(0:12) c c Initialize to zero all components of jet qj c INTEGER k DO k = 0, 12 qj( k ) = 0 END DO RETURN END c c ====================================================================== c Operations on particles/jets c SUBROUTINE p_get ( e, xta, phi, p, p_from ) DOUBLE PRECISION e, xta, phi, p(0:3), p_from(0:6) INCLUDE 'ojf_kin.fh' INCLUDE 'ojf_par.fh' INCLUDE 'ojf_com.fh' INTEGER k DO k = 0, 3 p(k) = p_from(k) END DO IF (ojf_kinematics .EQ. sphere) THEN e = p_from(0) IF ( ( p(1)**2 + p(2)**2 + p(3)**2 ) .NE. 0 ) & THEN xta = DATAN2( DSQRT( p(1)**2 + p(2)**2 ), & p(3) ) & / ojf_pi180 ELSE xta = 0 END IF ELSE IF (ojf_kinematics .EQ. cylinder) THEN e = p_from( par_Et ) xta = p_from( par_eta ) END IF IF ( DSQRT( p(1)**2 + p(2)**2 ) .NE. 0 ) THEN phi = DATAN2( p(2), p(1) ) / ojf_pi180 ELSE phi = 0 END IF RETURN END c c ---------------------------------------------------------------------- c SUBROUTINE p_set ( p, e, theta, phi ) DOUBLE PRECISION p(0:6), e, theta, phi c c set particle c INCLUDE 'ojf_kin.fh' INCLUDE 'ojf_par.fh' INCLUDE 'ojf_com.fh' IF (.NOT.( e.GE.0 .AND. theta.GE.0 .AND. & theta.LE.180 .AND. DABS(phi).LE.360 )) THEN WRITE(6,*) 'p_set: 31' STOP 'p_set: 31' END IF theta = theta * ojf_pi180 phi = phi * ojf_pi180 p(0) = e p(1) = e * DSIN( theta ) * DCOS( phi ) p(2) = e * DSIN( theta ) * DSIN( phi ) p(3) = e * DCOS( theta ) CALL p_set_tr( p ) RETURN END c c ---------------------------------------------------------------------- c SUBROUTINE p_set_raw ( p, px, py, pz ) DOUBLE PRECISION p(0:6), px, py, pz c c set particle c INCLUDE 'ojf_kin.fh' INCLUDE 'ojf_par.fh' INCLUDE 'ojf_com.fh' p(0) = DSQRT( px**2 + py**2 + pz**2 ) p(1) = px p(2) = py p(3) = pz CALL p_set_tr( p ) RETURN END c c ---------------------------------------------------------------------- c SUBROUTINE p_set_tr ( p ) DOUBLE PRECISION p(0:6) c c set cylinder-specific components for particle pa c INCLUDE 'ojf_kin.fh' INCLUDE 'ojf_par.fh' INCLUDE 'ojf_com.fh' DOUBLE PRECISION e_tr DOUBLE PRECISION asinh, x asinh(x) = DSIGN( DLOG( DSQRT( x**2 + 1 ) + DABS( x ) ), x ) IF (ojf_kinematics .EQ. cylinder) THEN e_tr = DSQRT( p(1)**2 + p(2)**2 ) p( par_Et ) = e_tr p( par_eta ) = asinh( p(3) / e_tr ) p( par_Eteta ) = e_tr * p( par_eta ) END IF RETURN END c c ====================================================================== c Operations on Lorentz vectors c SUBROUTINE v_add_mult ( v, z, w ) DOUBLE PRECISION v(0:3), w(0:3), z c c summation of 4-vectors with scalar weights: v <-- v + z * w c INTEGER i DO i = 0, 3 v(i) = v(i) + z * w(i) END DO RETURN END c c ---------------------------------------------------------------------- c DOUBLE PRECISION FUNCTION pos_prod ( v, w ) DOUBLE PRECISION v(0:3), w(0:3) c c Lorentz scalar product of two vectors, always >=0 in our context c INCLUDE 'ojf_par.fh' DOUBLE PRECISION res INTEGER i res = v(0) * w(0) DO i = 1, 3 res = res - v(i) * w(i) END DO IF (.NOT.( res .GE. -eps_round )) THEN WRITE(6,*) 'pos_prod: 60: huge rounding errors' STOP 'pos_prod: 60' END IF pos_prod = MAX( zero, res ) RETURN END c c ====================================================================== c DOUBLE PRECISION FUNCTION random () c c standard uniform deviate: 0 < random < 1 c we only need it simple and fast c INCLUDE 'ojf_kin.fh' INCLUDE 'ojf_par.fh' INCLUDE 'ojf_com.fh' DOUBLE PRECISION res INTEGER m, b, c PARAMETER ( m = random_m, b = 7141, c = 54773 ) ojf_ran_lock = .TRUE. ojf_random = MOD( ojf_random * b + c, m ) res = ojf_random + 1 res = res/(m + 1) random = res RETURN END c c ====================================================================== c SUBROUTINE aux_sort ( w, a, l ) INTEGER l DOUBLE PRECISION w( 1:l ) INTEGER a( 1:l ) c c Sort w in decreasing order; a is reordered in parallel. c c This sort is supposed to be used only for short arrays. c INTEGER j, aj, k DOUBLE PRECISION wj DO 10 j = l - 1, 1, -1 wj = w(j) aj = a(j) DO k = j + 1, l IF (w( k ) .GT. wj) THEN w( k-1 ) = w( k ) a( k-1 ) = a( k ) ELSE w( k-1 ) = wj a( k-1 ) = aj GOTO 10 END IF END DO w( l ) = wj a( l ) = aj 10 CONTINUE RETURN END c c END OF OJF CODE c ======================================================================