Coding violation fixes
[u/mrichter/AliRoot.git] / THydjet / hydjet1_1 / test_hydjet.f
CommitLineData
cb220f83 1*------------------------------------------------------------------------------
2*
3* Filename : TEST_HYDJET.F
4*
5*==============================================================================
6*
7* Description : Example program to simulate hadron spectra in AA collisions
8* at LHC with HYDJET-code (should be compiled with object files
9* obtained with hydjet1_1.f, pyquen1_1.f, pythia6401.f (or later
10* pythia versions) and jetset_73.f with extended array size of
11* common block LUJETS)
12*
13*==============================================================================
14
15 double precision ckin,parp,pari
16 real A,bmin,bmax,bfix
17 external ludata,pydata
18 common /lujets/ n,k(150000,5),p(150000,5),v(150000,5)
19 common /hyfpar/ bgen,nbcol,npart,npyt,nhyd
20 common /hyflow/ ytfl,ylfl,fpart
21 common /hyjpar/ nhsel,ptmin,njet
22 common /pydat1/ mstu(200),paru(200),mstj(200),parj(200)
23 common /pysubs/ msel,mselpd,msub(500),kfin(2,-40:40),ckin(200)
24 common /pypars/ mstp(200),parp(200),msti(200),pari(200)
25c common /pawc/ hmemor(20000)
26 save /lujets/,/hyflow/,/hyjpar/,/hyfpar/
27 save /pysubs/,/pypars/,/pydat1/
28
29* prepare hbook memory
30c call hlimit(20000)
31* open hrout file to write histograms
32c call HROPEN(1,'HISTO','hydjet.hrout','N',1024,ISTAT)
33* prepare hbook histograms
34c call hbook1(1,'dN/dy $',100,-10.,10.,0.) ! rapidity
35c call hbook1(2,'dN/deta $',100,-10.,10.,0.) ! pseudorapidity
36c call hbook1(3,'dN/dpt $',100,0.,10.,0.) ! transverse momentum
37c call hbook1(4,'dN/dphi $',100,-3.15,3.15,0.) ! azimuthal angle
38c call hbarx(0)
39
40* set initial beam parameters
41 A=207. ! atomic weigth
42 nh=20000 ! mean soft multiplicity for central Pb-Pb
43 ifb=0 ! fixed impact parameter
44 bfix=0. ! in nucleus radius units
45c ifb=1 ! distribution over impact parameter
46c bmin=0. ! from 'bmin'
47c bmax=1. ! to 'bmax'
48
49* set hydro parameters
50* nhsel=0 - hydro (no jets), nhsel=1 - hydro + pythia jets, nhsel=2 - hydro +
51* pyquen jets, nhsel=3 - pythia jets (no hydro), nhsel=4 - pyquen jets (no hydro)
52 nhsel=2 ! flag to include hard scatterings
53 ylfl=5. ! maximum longitudinal flow rapidity
54 ytfl=1. ! maximum transverse flow rapidity
55 fpart=1. ! fraction of soft multiplicity proportional
56 ! # of nucleons-participants
57* set input PYTHIA parameters
58 msel=1 ! QCD-dijet production
59 ptmin=10.
60 ckin(3)=dble(ptmin) ! minimum pt in initial hard sub-process
61 mstp(51)=7 ! CTEQ5M pdf
62 mstp(81)=0 ! pp multiple scattering off
63 mstu(21)=1 ! avoid stopping run
64 paru(14)=1.d0 ! tolerance parameter to adjust fragmentation
65
66* set original (rounded) test values and its rms for current model parameters
67 pta0=0.56
68 eta0=0.
69 dna0=29850.
70* set initial test values and its rms
71 ptam=0.
72 ptrms=0.
73 etam=0.
74 etrms=0.
75 dnam=0.
76 dnrms=0.
77
78* initialize PYTHIA for hard parton-parton scattering
79 if(nhsel.ne.0) call pyinit('CMS','p','p',5500.D0)
80
81* set number of generated events
82 ntot=10
83
84 do ne=1,ntot ! cycle on events
85 call hydro(A,ifb,bmin,bmax,bfix,nh)! single event generation
86
87 call luedit(2) ! remove unstable particles and partons
88
89 do ip=1,n ! cycle on particles
90 pt=plu(ip,10) ! transverse momentum...
91 ycm=plu(ip,17) ! rapidity...
92 eta=plu(ip,19) ! pseudorapidity...
93 phi=plu(ip,15) ! azimuthal angle...
94 charge=plu(ip,6) ! electric charge...
95
96* add current test values of eta, pt and its rms
97 etam=etam+eta
98 etrms=etrms+(eta-eta0)**2
99 ptam=ptam+pt
100 ptrms=ptrms+(pt-pta0)**2
101
102* fill histograms for charged particles
103c if(abs(charge).gt.0.) then
104c call hfill(1,ycm,0.,1.) ! rapidity
105c call hfill(2,eta,0.,1.) ! pseudorapidity
106c call hfill(3,pt,0.,1.) ! transverse momentum
107c call hfill(4,phi,0.,1.) ! azimuthal angle
108c end if
109 end do
110 write(6,*) 'Event #',ne
111 write(6,*) 'Impact parameter',bgen,'*RA',' Total multiplicity',n
112 write(6,*) 'Pt hard min',ptmin,' GeV',' Ndijets',njet
113 write(6,*) '***************************************************'
114
115* add current test value of event multiplicity and its rms
116 dnam=dnam+n
117 dnrms=dnrms+(n-dna0)**2
118 end do
119
120* test calculating and printing of original "true" (rounded) numbers
121* and generated one's (with statistical errors)
122 etam=etam/dnam
123 etrms=sqrt(etrms)/dnam
124 ptam=ptam/dnam
125 ptrms=sqrt(ptrms)/dnam
126 dnam=dnam/float(ntot)
127 dnrms=sqrt(dnrms)/float(ntot)
128 write(6,1) dna0
1291 format(2x,'True (rounded) mean multiplicity =',f7.0)
130 write(6,2) dnam, dnrms
1312 format(2x,'Generated mean multiplicity =',f7.0,3x,
132 > '+- ',f6.0)
133 write(6,3) eta0
1343 format(2x,'True (rounded) mean pseudorapidity =',f5.2)
135 write(6,4) etam, etrms
1364 format(2x,'Generated mean pseudorapidity =',f5.2,3x,
137 > '+- ',f5.2)
138 write(6,5) pta0
1395 format(2x,'True (rounded) mean transverse momentum =',f5.2)
140 write(6,6) ptam, ptrms
1416 format(2x,'Generated mean transverse momentum =',f5.2,3x,
142 > '+- ',f5.2)
143
144* finish histograms writing procedure
145c call hidopt(0,'ERRO')
146c call histdo
147c CALL HROUT(0,ICYCLE,' ')
148c CALL HREND('HISTO')
149
150 end
151*******************************************************************************