[u/mrichter/AliRoot.git] / THydjet / hydjet1_1 / test_hydjet.f

*------------------------------------------------------------------------------
*
* Filename             : TEST_HYDJET.F
*
*==============================================================================
*
* Description : Example program to simulate hadron spectra in AA collisions 
*               at LHC with HYDJET-code (should be compiled with object files
*               obtained with hydjet1_1.f, pyquen1_1.f, pythia6401.f (or later 
*               pythia versions) and jetset_73.f with extended array size of 
*               common block LUJETS) 
*
*==============================================================================

      double precision ckin,parp,pari 
      real A,bmin,bmax,bfix 
      external ludata,pydata 
      common /lujets/ n,k(150000,5),p(150000,5),v(150000,5)
      common /hyfpar/ bgen,nbcol,npart,npyt,nhyd
      common /hyflow/ ytfl,ylfl,fpart 
      common /hyjpar/ nhsel,ptmin,njet       
      common /pydat1/ mstu(200),paru(200),mstj(200),parj(200)
      common /pysubs/ msel,mselpd,msub(500),kfin(2,-40:40),ckin(200) 
      common /pypars/ mstp(200),parp(200),msti(200),pari(200)
c      common /pawc/ hmemor(20000)
      save /lujets/,/hyflow/,/hyjpar/,/hyfpar/
      save /pysubs/,/pypars/,/pydat1/ 

* prepare hbook memory 
c      call hlimit(20000)
* open hrout file to write histograms 
c      call HROPEN(1,'HISTO','hydjet.hrout','N',1024,ISTAT)
* prepare hbook histograms 
c      call hbook1(1,'dN/dy $',100,-10.,10.,0.)     ! rapidity 
c      call hbook1(2,'dN/deta $',100,-10.,10.,0.)   ! pseudorapidity 
c      call hbook1(3,'dN/dpt $',100,0.,10.,0.)      ! transverse momentum  
c      call hbook1(4,'dN/dphi $',100,-3.15,3.15,0.) ! azimuthal angle 
c      call hbarx(0)

* set initial beam parameters 
      A=207.                        ! atomic weigth        
      nh=20000                      ! mean soft 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=1.                      ! to 'bmax'  

* set hydro 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 
      ylfl=5.                        ! maximum longitudinal flow rapidity
      ytfl=1.                        ! maximum transverse flow rapidity
      fpart=1.                       ! fraction of soft multiplicity proportional 
                                      ! # of nucleons-participants
* set input PYTHIA parameters 
      msel=1                        ! QCD-dijet production 
      ptmin=10.
      ckin(3)=dble(ptmin)           ! minimum pt in initial hard sub-process 
      mstp(51)=7                    ! CTEQ5M pdf 
      mstp(81)=0                    ! pp multiple scattering off 
      mstu(21)=1                    ! avoid stopping run 
      paru(14)=1.d0                 ! tolerance parameter to adjust fragmentation 

* set original (rounded) test values and its rms for current model parameters 
      pta0=0.56
      eta0=0.         
      dna0=29850.   
* set initial test values and its rms 
      ptam=0. 
      ptrms=0. 
      etam=0.  
      etrms=0.        
      dnam=0.  
      dnrms=0.  
       
* initialize PYTHIA for hard parton-parton scattering 
      if(nhsel.ne.0) call pyinit('CMS','p','p',5500.D0)       
  
* set number of generated events 
      ntot=10

      do ne=1,ntot                        ! cycle on events 
       call hydro(A,ifb,bmin,bmax,bfix,nh)! single event generation

       call luedit(2)                     ! remove unstable particles and partons 
                      
       do ip=1,n                          ! cycle on particles        
        pt=plu(ip,10)                     ! transverse momentum... 
        ycm=plu(ip,17)                    ! rapidity...  
        eta=plu(ip,19)                    ! pseudorapidity...  
        phi=plu(ip,15)                    ! azimuthal angle...
        charge=plu(ip,6)                  ! electric charge...

* add current test values of eta, pt and its rms 
        etam=etam+eta 
	etrms=etrms+(eta-eta0)**2 
	ptam=ptam+pt  
	ptrms=ptrms+(pt-pta0)**2
           	 
* fill histograms for charged particles 
c        if(abs(charge).gt.0.) then 
c         call hfill(1,ycm,0.,1.)          ! rapidity    
c         call hfill(2,eta,0.,1.)          ! pseudorapidity 
c         call hfill(3,pt,0.,1.)           ! transverse momentum 
c         call hfill(4,phi,0.,1.)          ! azimuthal angle  
c        end if 
       end do 
       write(6,*) 'Event #',ne
       write(6,*) 'Impact parameter',bgen,'*RA',' Total multiplicity',n 
       write(6,*) 'Pt hard min',ptmin,' GeV','   Ndijets',njet 
       write(6,*) '***************************************************'

* add current test value of event multiplicity and its rms 
       dnam=dnam+n          
       dnrms=dnrms+(n-dna0)**2 
      end do 

* test calculating and printing of original "true" (rounded) numbers 
* and generated one's (with statistical errors) 
      etam=etam/dnam
      etrms=sqrt(etrms)/dnam
      ptam=ptam/dnam 
      ptrms=sqrt(ptrms)/dnam
      dnam=dnam/float(ntot)
      dnrms=sqrt(dnrms)/float(ntot) 
      write(6,1) dna0
1     format(2x,'True (rounded) mean multiplicity =',f7.0) 
      write(6,2) dnam, dnrms 
2     format(2x,'Generated mean multiplicity      =',f7.0,3x,
     > '+-  ',f6.0) 
      write(6,3) eta0
3     format(2x,'True (rounded) mean pseudorapidity =',f5.2) 
       write(6,4) etam, etrms 
4     format(2x,'Generated mean pseudorapidity      =',f5.2,3x,
     > '+-  ',f5.2) 
      write(6,5) pta0
5     format(2x,'True (rounded) mean transverse momentum =',f5.2)   
      write(6,6) ptam, ptrms 
6     format(2x,'Generated mean transverse momentum      =',f5.2,3x,
     > '+-  ',f5.2)    
 
* finish histograms writing procedure        
c      call hidopt(0,'ERRO')
c      call histdo
c      CALL HROUT(0,ICYCLE,' ')
c      CALL HREND('HISTO')
       
      end
*******************************************************************************
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	*******************************************************************************