hijing/examples/example_3.f

0001 C****************************************************************************
0002 C          Program # 2 from Comp. Phys. Commun. 83 (1994) 307
0003 C                  by M. Gyulassy and X-.N. Wang
0004 C             Modified by V.Uzhinsky, CERN, Oct. 2003
0005 C***************************************************************************
0006
0007       CHARACTER FRAME*8,PROJ*8,TARG*8
0008
0009       COMMON/HIPARNT/HIPR1(100),IHPR2(50),HINT1(100),IHNT2(50)
0010       SAVE  /HIPARNT/
0011
0012       COMMON/HIMAIN1/NATT,EATT,JATT,NT,NP,N0,N01,N10,N11
0013       SAVE  /HIMAIN1/
0014
0015       COMMON/HIMAIN2/KATT(130000,4),PATT(130000,4)
0016       SAVE  /HIMAIN2/
0017 C
0018       DIMENSION GB(101), XB(101), DNDP(50)
0019
0020       COMMON/RANSEED/NSEED
0021       SAVE  /RANSEED/
0022
0023       NSEED=0
0024
0025       do i=1,50
0026          DNDP(i)=0.
0027       enddo
0028
0029       FRAME='CMS'
0030
0031       write(6,*)'===================================================='
0032       write(6,*)'  Calculation of transverse momentum distribution   '
0033       write(6,*)'        of charged pions in minimum bias            '
0034       write(6,*)'            hh-, hA- and AA-collisions              '
0035       write(6,*)'    Calculation will be performed in CM System      '
0036       write(6,*)'===================================================='
0037
0038       write(6,*)
0039       write(6,*)'Enter the energy per NN-collision (GeV)'
0040       read(5,*)  EFRM
0041
0042       write(6,*)
0043       write(6,*)'Enter a type of the "projectile" particle'
0044       write(6,*)
0045       write(6,*)' P proton,            PBAR anti-proton,'
0046       write(6,*)' N neutron,           NBAR anti-neutron,'
0047       write(6,*)' PI+ - positive pion, PI- negative pion,'
0048       write(6,*)' K+ positive kaon,    K- negative kaon'
0049       write(6,*)
0050       write(6,*)' A - nucleus --------------------------'
0051
0052       read(5,1) PROJ
0053 1     format(A8)
0054
0055       if(PROJ.ne.'A') then
0056         IAP=1
0057         if(PROJ.eq.'P'   ) IZP= 1
0058         if(PROJ.eq.'PBAR') IZP=-1
0059         if(PROJ.eq.'N'   ) IZP= 0
0060         if(PROJ.eq.'NBAR') IZP= 0
0061         if(PROJ.eq.'PI+' ) IZP= 1
0062         if(PROJ.eq.'PI-' ) IZP=-1
0063         if(PROJ.eq.'K+'  ) IZP= 1
0064         if(PROJ.eq.'K-'  ) IZP=-1
0065       else
0066         write(6,*)
0067         write(6,*)'Enter mass number and charge of the proj. nucleus'
0068         read(5,*)  IAP, IZP
0069       endif
0070
0071       write(6,*)
0072       write(6,*)'Enter a type of the "target" particle (same notations)'
0073       read(5,1) TARG
0074
0075       if(TARG.ne.'A') then
0076         IAT=1
0077         if(TARG.eq.'P'   ) IZT= 1
0078         if(TARG.eq.'PBAR') IZT=-1
0079         if(TARG.eq.'N'   ) IZT= 0
0080         if(TARG.eq.'NBAR') IZT= 0
0081         if(TARG.eq.'PI+' ) IZT= 1
0082         if(TARG.eq.'PI-' ) IZT=-1
0083         if(TARG.eq.'K+'  ) IZT= 1
0084         if(TARG.eq.'K-'  ) IZT=-1
0085       else
0086         write(6,*)
0087         write(6,*)'Enter mass number and charge of the target nucleus'
0088         read(5,*)  IAT, IZT
0089       endif
0090
0091       write(6,*)
0092       write(6,*)'Enter number of events per each value of impact',
0093      ,' parameter (e.a. 10)'
0094
0095       read(5,*)  N_EVENT
0096
0097 C....initialize HIJING for requested interactions
0098       CALL HIJSET(EFRM,FRAME,PROJ,TARG,IAP,IZP,IAT,IZT)
0099
0100       write(6,*)' Simulation of interactions with'
0101       write(6,*)
0102       write(6,*)' Proj = ',PROJ,' and  Targ = ',TARG
0103       write(6,*)' IAP  =',IAP  ,'            IAT  =',IAT
0104       write(6,*)' IZP  =',IZP  ,'            IZT  =',IZT
0105       write(6,*)
0106       write(6,*)' Reference frame -   ',FRAME
0107       write(6,*)' ENERGY            ',EFRM,' GeV'
0108       write(6,*)' Number of generated events per B interval -',N_event
0109       write(6,*)
0110
0111 C....set BMIN=0 and BMAX=R_A+R_B
0112
0113       BMIN=0.0
0114       BMAX=HIPR1(34)+HIPR1(35)
0115
0116 C....calculate the Glauber probability and its integrated value:
0117
0118       DIP=(BMAX-BMIN)/100.0
0119       GBTOT=0.0
0120       DO 100 I=1,101
0121          XB(I)=BMIN+(I-1)*DIP
0122          OV=PROFILE(XB(I))
0123          GB(I)=XB(I)*(1.0-exp(-HINT1(12)*OV))
0124          GBTOT=GBTOT+GB(I)
0125 100   CONTINUE
0126       write(6,*)'Inelastic X-section (mb) ',GBTOT*DIP*10.*6.28
0127
0128       PAUSE
0129 C....generating N_EVENT for each of 100 impact parameter intervals:
0130
0131       NONT=0
0132       GNORM=GBTOT
0133
0134       DO 300 IB=1,100
0135          B1=XB(iB)
0136          B2=XB(IB+1)
0137
0138 C....normalized Glauber probability:
0139
0140          W_GB=(GB(IB)+GB(IB+1))/2.0/GBTOT
0141
0142          DO 200 IE=1,N_EVENT
0143             CALL HIJING(FRAME,B1,B2)
0144
0145 C....count number of events without any interaction
0146 C....and renormalize the total Glauber prabability:
0147
0148             IF(NP+NT .EQ. 0) THEN
0149               NONT=NONT+1
0150               GNORM=GNORM-GB(IB)/FLOAT(N_EVENT)
0151               GO TO 200
0152             ENDIF
0153
0154 C....calculate pt distribution of charged pions:
0155
0156             DO 150 K=1,NATT
0157 C....select charged pions only:
0158                IF(ABS(KATT(K,1)) .NE. 211) GO TO 150
0159
0160 C....calculate pt:
0161                PTR=SQRT(PATT(K,1)**2+PATT(K,2)**2)
0162
0163 C....calculate pt distribution and weight with normalized
0164 C....Glauber probability to get minimum bias results:
0165
0166                IF(PTR .GE. 10.0) GO TO 150
0167                IPT=1+PTR/0.2
0168                DNDP(IPT)=DNDP(IPT)+W_GB/FLOAT(N_EVENT)/0.2
0169 150         CONTINUE
0170 200      CONTINUE
0171 300   CONTINUE
0172
0173 C....renormalize the distribution by the renormalized Glauber
0174 C....probability which excludes the events without any interaction:
0175
0176       IF(NONT.NE.0) THEN
0177         DO 400 I=1,50
0178            DNDP(I)=DNDP(I)*GBTOT/GNORM
0179            if(DNDP(I).ne.0.) write(6,350)0.2*(i-1),DNDP(I)
0180 350        format(1x,f5.1,2x,e11.4)
0181 400     CONTINUE
0182       ENDIF
0183
0184       STOP
0185       END
0186
0187       FUNCTION RAN(NSEED)
0188       RAN=RLU(NSEED)
0189       RETURN
0190       END