-------------------------------------------------------------
      HYDJET, event generator for simulation of jet quenching and   
      flow effects in ultrarelativistic heavy ion AA collisions
      -------------------------------------------------------------
      This code is merging HYDRO (flow effects), PYTHIA6.4xx (hard 
      jet production) and PYQUEN (jet quenching)
      --------------------------------------------------------------

      Igor Lokhtin, SINP MSU, Moscow, RU
      e-mail: Igor.Lokhtin@cern.ch 

      Main reference for HYDJET and PYQUEN: 
      I.P. Lokhtin, A.M. Snigirev, Eur. Phys. J. C 45 (2006) 211.

      Reference for PYTHIA6.4:  
      T.Sjostrand, S. Mrenna and P. Skands, JHEP05 (2006) 026.
      
      Reference for JETSET event format:  
      T.Sjostrand, Comput.Phys.Commun. 82 (1994) 74.  
      
      Reference for nuclear shadowing model:
      K. Tywoniuk, I.C. Arsene, L. Bravina, A. Kaidalov and 
      E. Zabrodin, Phys. Lett. B 657 (2007) 170.
      
      -------------------------------------------------------------- 
      Web-page: 
      http://cern.ch/lokhtin/hydro/hydjet.html
      --------------------------------------------------------------

      Description of routine

HYDJET initialization should be done at the beginning of main user's routine: 

CALL hyinit(energy,A,ifb,bmin,bmax,bfix,nh) 
       
input parameters to specify event configuration: 
(set in main user's routine before call HYINIT)

  energy - c.m.s energy per nucleon pair in GeV;
  A    - beam and target nucleus atomic weight; 
  ifb  - flag of type of centrality generation 
         =0 impact parameter is fixed (bfix)  
	 >0 or <0 impact parameter is generated with standard Glauber geometry 
	          between minimum (bmin) and maximum (bmax) values; 
  bmin - minimum impact parameter in units of nucleus radius RA
         (i.e. minimum value in [fm] will be bmin*RA), 
	  valid only if ifb not equal to zero
         (allowed range is 0 < bmin < bmax);   
  bmax - maximum impact parameter in units of nucleus radius RA
         (i.e. maximum value in [fm] will be bmax*RA),
	  valid only if ifb not equal to zero 
         (allowed range is bmin < bmax < 3); 
  bfix - fixed impact parameter in units of nucleus radius RA 
         (i.e. fixed value in [fm] will be bfix*RA), 
	 valid only if ifb=0
         (allowed range is 0 < bfix < 3);
  nh   - mean primary hadron multiplicity for soft component in central 
         PbPb collisions (multiplicity for other centralities and atomic 
         numbers is calculated automatically).  
	 

HYINIT includes:    
- PYTHIA initialization at given c.m.s. energy per nucleon pair;
- calculation of total inelastic NN cross section at given energy;
- calculation of hard scattering NN cross section at given ptmin and energy;
- tabulation of nuclear thickness function and nuclear overlap function; 
- calculation of number of participants & binary collisions at Pb+Pb (b=0); 
- test of input parameters. 

CALL hyevnt(bfix) - generates one event

  bfix - fixed impact parameter in units of nucleus radius RA,
         valid only if ifb=0;
	 the same as for hyinit. 
     
      --------------------------------------------------------------

Parameters in COMMON BLOCKS which can be varied by user: 

PYQUEN:

COMMON /pyqpar/ T0,tau0,nf,ienglu,ianglu
   T0 - maximum initial temperature of quark-gluon plasma for central Pb+Pb 
   collisions at mid-rapidity (initial temperature for other centralities and 
   atomic numbers will be calculated automatically);  
   (allowed range is 0.2 GeV < T0 < 2 GeV, default value is T0=1 GeV);
   tau0 - proper time of quark-gluon plasma formation 
   (allowed range is 0.01 < tau0 < 10 fm/c, default value is tau0=0.1 fm/c);
   nf - number of active quark flavours in quark-gluon plasma
   (nf=0, 1, 2 or 3, default value is nf=0); 
   ienglu - flag to fix type of medium-induced partonic energy loss
   (ienglu=0 - radiative and collisional loss, 
    ienglu=1 - radiative loss only, ienglu=2 - collisional loss only,
    default value is ienglu=0);
   ianglu - flag to fix type of angular distribution of emitted gluons  
   (ianglu=0 - small-angular, ianglu=1 - wide-angular, ianglu=2 - collinear, 
    default value is ianglu-0).  

HYDJET:

COMMON /hyflow/ ytfl,ylfl,Tf,fpart 
ytfl - maximum transverse collective rapidity, controls slope of low-pt spectra 
(allowed range is 0.01<ytfl<3.0, default value is ytfl=1.5);
ylfl - maximum longitudinal collective rapidity, controls width of eta-spectra
(allowed range is 0.01<ylfl<7.0, default value is ylfl=4.);
Tf - hadron freeze-out temperature 
(allowed range is 0.08 GeV < Tf < 0.2 GeV, default value is Tf=0.1 GeV); 
fpart - fraction of soft multiplicity proportional to the number of nucleon
participants; then (1.-fpart) will be the fraction of soft multiplicity
proportional to the number of nucleon-nucleon binary sub-collisions  
(allowed range is 0.01<fpart<1.0, default value is fpart=1.). 

COMMON /hyjpar/ ptmin,sigin,sigjet,nhsel,ishad,njet
nhsel - flag to include jet production in hydro-type event: 
nhsel=0 - jet production off (pure HYDRO event),
nhsel=1 - jet production on, jet quenching off (HYDRO+njet*PYTHIA events),
nhsel=2 - jet production & jet quenching on (HYDRO+njet*PYQUEN events),
nhsel=3 - jet production on, jet quenching off, HYDRO off (njet*PYTHIA events),
nhsel=4 - jet production & jet quenching on, HYDRO off (njet*PYQUEN events);
ishad - flag to include impact parameter dependent nuclear shadowing 
correction for gluons and light sea quarks (u,d,s) in hard part of event: 
ishad=0 - no shadowing, ishad=1 - shadowing is included;
ptmin- minimal pt of parton-parton scattering in PYTHIA event in GeV, is equal 
to parameter ckin(3) in PYTHIA common block /pysubs/;
sigin - total inelastic NN cross section at given energy in mb, is calculated 
by Pythia, but also can be changed by user in main routine before call HYINIT 
(allowed range is 10 mb < sigin < 200 mb, default value is given by PYTHIA). 

NOTE! If specified by user value of such parameter is out of allowed range, 
      the default value is used in HYDJET run. 

NOTE! Default parameters of quark-gluon plasma (T0, tau0, nf) in PYQUEN were 
      selected as an estimation for LHC heavy ion beam energies. 

NOTE! At this moment nuclear shadowing parameterization is implemented only 
      for Pb, Au, Pd or Ca ions. 

      -------------------------------------------------------------- 

Output event parameters in COMMON BLOCKS: 

COMMON /hyjpar/ ptmin,sigin,sigjet,nhsel,ishad,njet
njet - number of hard parton-parton scatterings with pt>ptmin in event;
sigin - total inelastic NN cross section at given c.m.s. energy (in mb);
sigjet - hard scattering NN cross section at given ptmin and energy (in mb).

common /hyfpar/ bgen,nbcol,npart,npyt,nhyd
bgen - generated value of impact parameter in units of nucleus radius RA 
(i.e the value in [fm] will be bgen*RA).
nbcol - mean number of nucleon-nucleon binary sub-collisions at given 'bgen'. 
npart - mean number of nucleon participants at given 'bgen'. 
npyt - multiplicity of hard PYTHIA/PYQUEN-induced particles in event 
       (including full parton and resonance decay story).  
nhyd - multiplicity of soft hydro-induced particles in event
       (including resonance decay story).  

NOTE! NBCOL and NPART are DOUBLE PRECISION variables, and so they should be 
      described as 'double precision' in main user's routine. 

      -------------------------------------------------------------- 

Output particle information is stored using standard JETSET/PYTHIA event 
format in common blocks HYJETS (in HYDJET) and LUJETS (in JETSET). 

COMMON /lujets/ n,k(150000,5),p(150000,5),v(150000,5)
n - total event multiplicity    
k(i,1-5) - particle codes       
p(i,1-5) - particle four-momentum and mass 
v(i,1-5) - particle vertex, production time and lifetime   

First 'npyt' lines in event list correspond to PYTHIA/PYQUEN-induced particles
(including partons, strings and unstable and final hadrons), 
last 'nhyd' lines -- soft hydro-induced hadrons.  

NOTE! Since jetset_73.f operates on REAL variables, arrays 'p' and 'v', as 
well as other REAL parameters and functions of JETSET to be used, should be 
described as 'real' in main user's routine. Copying the original DOUBLE 
PRECISION output particle information from HYDJET arrays (common block HYJETS,
see below) to JETSET arrays allows the user easily to perform some useful 
manipulations with the event record using standard JETSET subroutines and 
functions (to exclude unstable or undetectable particles, to list an event, to 
provide various event data, etc.) 

COMMON /hyjets/ nhj,nhp,khj(150000,5),phj(150000,5),vhj(150000,5)  
contains original event record in the same format as /lujets/  
(plus additional parameter 'nhp' - dummy to ensure as even number of integers):
n=nhj, k(i,j)=khj(i,j), p(i,j)=phj(i,j), v(i,j)=vhj, 
but in DOUBLE PRECISION for 'phj' and 'vhj' arrays.  

      --------------------------------------------------------------  
    
NOTE! Main user's routine should be compiled with object files obtained by:  
      - latest pythia (pythia6.401.f or later versions);
      - pyquen1_5.f;  
      - hydjet1_7.f;
      - jetset_73.f with extended array size in common block /lujets/ 
        (if using the standard JETSET subroutines and functions to manipulate 
	 with the event record and to provide various event data is needed).  
                 
      -------------------------------------------------------------------------    
        
      Physics validity of the model 
      
1. Hydro-type approximation for heavy ion collisions is expected to be 
applicable for central and semi-central collisions. The result obtained for 
very peripheral collisions (b~2*RA) can be not adequate. 
2. We do not expect correct event description for very forward rapidities 
(|y|>3), where other mechanisms of particle production among hydro-flow and
jets can be important. 
3. Physics model for description of medium-induced gluon radiation is based on 
Baier-Dokshitzer-Mueller-Schiff formalism and valid for relatively high
transverse momenta of jet partons (>> 1 GeV/c). Thus setting reasonably high  
values of minimum pt in initial hard parton sub-processes (parameter 'ptmin') 
is preferable.