*------------------------------------------------------------------------------
*
* Filename : TEST_HYDJET.F
*
*==============================================================================
*
* Description : Example program to simulate hadron spectra in AA collisions
* with HYDJET generator. It should be compiled with hydjet1_9.f,
* pyquen1_5.f and latest pythia (e.g. pythia-6.4.24.f),
* and also jetset_73.f with the extended array size of 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).
*
*==============================================================================
IMPLICIT DOUBLE PRECISION(A-H, O-Z)
double precision nbcol,npart
real p,v,plu
integer pycomp
external ludata,pydata
common /lujets/ n,k(150000,5),p(150000,5),v(150000,5)
common /hyjets/ nhj,nhp,khj(150000,5),phj(150000,5),vhj(150000,5)
common /hyfpar/ bgen,nbcol,npart,npyt,nhyd
common /hyflow/ ytfl,ylfl,Tf,fpart
common /hyjpar/ ptmin,sigin,sigjet,nhsel,ishad,njet
common /pydat1/ mstu(200),paru(200),mstj(200),parj(200)
common /pydat3/ mdcy(500,3),mdme(8000,2),brat(8000),kfdp(8000,5)
common /pysubs/ msel,mselpd,msub(500),kfin(2,-40:40),ckin(200)
common /pypars/ mstp(200),parp(200),msti(200),pari(200)
common /pyqpar/ T0,tau0,nf,ienglu,ianglu
save /lujets/,/hyjets/,/hyflow/,/hyjpar/,/hyfpar/,/pyqpar/,
> /pysubs/,/pypars/,/pydat1/
* set beam parameters:
energy=2760. ! c.m.s energy per nucleon pair in GeV
A=207. ! atomic weigth
nh=18500 ! mean soft hadron multiplicity for central Pb+Pb
ifb=0 ! fixed impact parameter
bfix=0. ! in nucleus radius units
c ifb=1 ! distribution over impact parameter
c bmin=0. ! from 'bmin'
c bmax=3. ! to 'bmax'
* set of input HYDJET parameters:
* nhsel=0 - hydro (no jets), nhsel=1 - hydro + pythia jets, nhsel=2 - hydro +
* pyquen jets, nhsel=3 - pythia jets (no hydro), nhsel=4 - pyquen jets (no hydro)
nhsel=2 ! flag to include hard scatterings
* ishad=0 - no shadowing, ishad=1 - shadowing is included
ishad=1 ! flag to include "nuclear shadowing"
ylfl=4.5 ! maximum longitudinal flow rapidity
ytfl=1.15 ! maximum transverse flow rapidity
Tf=0.125 ! freeze-out temperature in GeV
fpart=1. ! fraction of soft multiplicity proportional
! # of nucleons-participants
sigin=64. ! inelastic NN cross-section in mb
* set of input PYQUEN parameters:
* ienglu=0 - radiative and collisional loss, ienglu=1 - only radiative loss,
* ienglu=2 - only collisional loss;
* ianglu=0 - small-angular radiation, ianglu=1 - wide angular radiation,
* inanglu=2 - collinear radiation
ienglu=0 ! set type of partonic energy loss
ianglu=1 ! set angular spectrum of gluon radiation
T0=1. ! initial QGP temperature
tau0=0.1 ! proper time of QGP formation
nf=0 ! number of active quark flavours in QGP
* set input PYTHIA parameters:
msel=1 ! QCD-dijet production
mstj(22)=2 ! decay those unstable particles
parj(71)=10. ! for which ctau < 10 mm
ckin(3)=7.5 ! minimum pt in initial hard scattering, GeV
mstu(21)=1 ! avoid stopping run
paru(14)=1. ! tolerance parameter to adjust fragmentation
* Pro-Q2O tune
mstp(51)=7 ! PDF set
mstp(3)=2 ! QCD switch for choice of LambdaQCD
parp(62)=2.9 ! ISR IR cutoff
parp(64)=0.14 ! ISR renormalization scale prefactor
parp(67)=2.65 ! ISR Q2max factor
mstp(68)=3 ! ISR phase space choice & ME corrections
parp(71)=4. ! FSR Q2max factor for non-s-channel procs
parj(81)=0.29 ! FSR Lambda_QCD scale
parj(82)=1.65 ! FSR IR cutoff
mstp(33)=0 ! "K" switch for K-factor on/off & type
mstp(81)=1 ! UE model
parp(82)=1.9 ! UE IR cutoff at reference ecm
parp(89)=1800. ! UE IR cutoff reference ecm
parp(90)=0.22 ! UE IR cutoff ecm scaling power
mstp(82)=4 ! UE hadron transverse mass distribution
parp(83)=0.83 ! UE mass distribution parameter
parp(84)=0.6 ! UE mass distribution parameter
parp(85)=0.86 ! UE gg colour correlated fraction
parp(86)=0.93 ! UE total gg fraction
mstp(91)=1 ! BR primordial kT distribution
parp(91)=2.1 ! BR primordial kT width <|kT|>
parp(93)=5. ! BR primordial kT UV cutoff
mstj(11)=5 ! HAD choice of fragmentation function(s)
parj(1)=0.073 ! HAD diquark suppression
parj(2)=0.2 ! HAD strangeness suppression
parj(3)=0.94 ! HAD strange diquark suppression
parj(4)=0.032 ! HAD vector diquark suppression
parj(11)=0.31 ! HAD P(vector meson), u and d only
parj(12)=0.4 ! HAD P(vector meson), contains s
parj(13)=0.54 ! HAD P(vector meson), heavy quarks
parj(21)=0.325 ! HAD fragmentation pT
parj(25)=0.63 ! HAD eta0 suppression
parj(26)=0.12 ! HAD eta0' suppression
parj(41)=0.5 ! HAD string parameter a
parj(42)=0.6 ! HAD string parameter b
parj(46)=1. ! HAD Lund(=0)-Bowler(=1) rQ (rc)
parj(47)=0.67 ! HAD Lund(=0)-Bowler(=1) rb
* set original test values of mean hadron transverse momentum and event multiplicity (rounded)
pta0=0.5
dna0=45640
* set initial test values and their rms
ptam=0.d0
ptrms=0.d0
dnam=0.d0
dnrms=0.d0
* initialize HYDJET at given input parameters
call hyinit(energy,A,ifb,bmin,bmax,bfix,nh)
* set number of generated events
ntot=100
do ne=1,ntot ! cycle on events
call hyevnt(bfix) ! single event generation
call luedit(2) ! remove unstable particles and partons
* reset current test value of pt
ptamc=0.d0
if(n.ge.1) then
do ip=1,n ! cycle on particles
pt=plu(ip,10) ! transverse momentum...
* add current test value of pt and its rms
ptamc=ptamc+pt
ptrms=ptrms+(pt-pta0)**2
end do
* add current test values of mean pt and event multiplicity and their rms
ptam=ptam+ptamc/n
end if
dnam=dnam+n
dnrms=dnrms+(n-dna0)**2
write(6,*) 'Event #',ne
write(6,*) 'Impact parameter',bgen,'*RA',' Total multiplicity',n
write(6,*) 'Pt hard min',ptmin,' GeV',' Ndijets',njet
write(6,*) '***************************************************'
end do
* test calculation and printout of original "true" numbers
* and generated one's (with statistical errors)
ptam=ptam/ntot
ptrms=dsqrt(ptrms)/dnam
dnam=dnam/ntot
dnrms=dsqrt(dnrms)/ntot
write(6,1) dna0
1 format(2x,'True (rounded) mean multiplicity =',d12.5)
write(6,2) dnam, dnrms
2 format(2x,'Generated mean multiplicity =',d12.5,3x,
> '+- ',d11.4)
write(6,5) pta0
5 format(2x,'True mean transverse momentum =',d9.2)
write(6,6) ptam, ptrms
6 format(2x,'Generated mean transverse momentum =',d9.2,3x,
> '+- ',d9.2)
end
*******************************************************************************