bugfix: boundery check for static hit array
[u/mrichter/AliRoot.git] / THydjet / hydjet1_1 / test_hydjet.f
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 *******************************************************************************