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 *******************************************************************************