]>
Commit | Line | Data |
---|---|---|
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) | |
25 | c common /pawc/ hmemor(20000) | |
26 | save /lujets/,/hyflow/,/hyjpar/,/hyfpar/ | |
27 | save /pysubs/,/pypars/,/pydat1/ | |
28 | ||
29 | * prepare hbook memory | |
30 | c call hlimit(20000) | |
31 | * open hrout file to write histograms | |
32 | c call HROPEN(1,'HISTO','hydjet.hrout','N',1024,ISTAT) | |
33 | * prepare hbook histograms | |
34 | c call hbook1(1,'dN/dy $',100,-10.,10.,0.) ! rapidity | |
35 | c call hbook1(2,'dN/deta $',100,-10.,10.,0.) ! pseudorapidity | |
36 | c call hbook1(3,'dN/dpt $',100,0.,10.,0.) ! transverse momentum | |
37 | c call hbook1(4,'dN/dphi $',100,-3.15,3.15,0.) ! azimuthal angle | |
38 | c 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 | |
45 | c ifb=1 ! distribution over impact parameter | |
46 | c bmin=0. ! from 'bmin' | |
47 | c 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 | |
103 | c if(abs(charge).gt.0.) then | |
104 | c call hfill(1,ycm,0.,1.) ! rapidity | |
105 | c call hfill(2,eta,0.,1.) ! pseudorapidity | |
106 | c call hfill(3,pt,0.,1.) ! transverse momentum | |
107 | c call hfill(4,phi,0.,1.) ! azimuthal angle | |
108 | c 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 | |
129 | 1 format(2x,'True (rounded) mean multiplicity =',f7.0) | |
130 | write(6,2) dnam, dnrms | |
131 | 2 format(2x,'Generated mean multiplicity =',f7.0,3x, | |
132 | > '+- ',f6.0) | |
133 | write(6,3) eta0 | |
134 | 3 format(2x,'True (rounded) mean pseudorapidity =',f5.2) | |
135 | write(6,4) etam, etrms | |
136 | 4 format(2x,'Generated mean pseudorapidity =',f5.2,3x, | |
137 | > '+- ',f5.2) | |
138 | write(6,5) pta0 | |
139 | 5 format(2x,'True (rounded) mean transverse momentum =',f5.2) | |
140 | write(6,6) ptam, ptrms | |
141 | 6 format(2x,'Generated mean transverse momentum =',f5.2,3x, | |
142 | > '+- ',f5.2) | |
143 | ||
144 | * finish histograms writing procedure | |
145 | c call hidopt(0,'ERRO') | |
146 | c call histdo | |
147 | c CALL HROUT(0,ICYCLE,' ') | |
148 | c CALL HREND('HISTO') | |
149 | ||
150 | end | |
151 | ******************************************************************************* |