hijing/examples/example_5.f

0001
0002 C****************************************************************************
0003 C       The program for checking conservation laws fulfillment in
0004 C                         the HIJING model.
0005 C                 Written by V.Uzhinsky, CERN, Oct. 2003
0006 C***************************************************************************
0007
0008       CHARACTER FRAME*8,PROJ*8,TARG*8
0009
0010       COMMON/HIMAIN1/ NATT,EATT,JATT,NT,NP,N0,N01,N10,N11
0011       SAVE  /HIMAIN1/
0012
0013       COMMON/HIMAIN2/KATT(130000,4),PATT(130000,4)
0014       SAVE  /HIMAIN2/
0015
0016 C         ********information of produced particles
0017 C
0018       CHARACTER *1 KEY
0019       CHARACTER *12 FNAME
0020
0021       PARAMETER (NWPAWC=50000)
0022       COMMON/PAWC/HMEMORY(NWPAWC)
0023
0024       COMMON/RANSEED/NSEED
0025       SAVE  /RANSEED/
0026       NSEED=0
0027
0028       write(6,*)'====================================================='
0029       write(6,*)'The program for checking conservation laws fulfillment'
0030       write(6,*)'              in the HIJING model                    '
0031       write(6,*)'  Only hadronic (hh) interactions are considered     '
0032       write(6,*)'====================================================='
0033       write(6,*)'        You can work in Lab. or CM systems           '
0034       write(6,*)' Would you like to use CM system? Y - Yes, N - No?   '
0035
0036       read(5,1001) KEY
0037 1001  format(A1)
0038
0039       if(KEY.eq.'Y'.or.KEY.eq.'y') then
0040         write(6,*)'CM system is used --------------------------------'
0041         FRAME='CMS'
0042       else
0043         write(6,*)'LAB system is used -------------------------------'
0044         FRAME='LAB'
0045       endif
0046
0047       write(6,*)'Enter the corresponding energy (GeV)'
0048
0049       read(5,*) EFRM
0050
0051       write(6,*)
0052       write(6,*)'Enter a type of the projectile particle'
0053       write(6,*)
0054       write(6,*)' P proton,            PBAR anti-proton,'
0055       write(6,*)' N neutron,           NBAR anti-neutron,'
0056       write(6,*)' PI+ - positive pion, PI- negative pion,'
0057       write(6,*)' K+ positive kaon,    K- negative kaon'
0058
0059       read(5,1002) PROJ
0060 1002  format(A8)
0061
0062       if(PROJ.ne.'A') then
0063         IAP=1
0064         if(PROJ.eq.'P'   ) IZP= 1
0065         if(PROJ.eq.'PBAR') IZP=-1
0066         if(PROJ.eq.'N'   ) IZP= 0
0067         if(PROJ.eq.'NBAR') IZP= 0
0068         if(PROJ.eq.'PI+' ) IZP= 1
0069         if(PROJ.eq.'PI-' ) IZP=-1
0070         if(PROJ.eq.'K+'  ) IZP= 1
0071         if(PROJ.eq.'K-'  ) IZP=-1
0072       else
0073         write(6,*)
0074         write(6,*)'Enter mass number and charge of the proj. nucleus'
0075         read(5,*)  IAP, IZP
0076       endif
0077
0078       write(6,*)
0079       write(6,*)'Enter a type of the target particle (same notations)'
0080       read(5,1002) TARG
0081
0082       if(TARG.ne.'A') then
0083         IAT=1
0084         if(TARG.eq.'P'   ) IZT= 1
0085         if(TARG.eq.'PBAR') IZT=-1
0086         if(TARG.eq.'N'   ) IZT= 0
0087         if(TARG.eq.'NBAR') IZT= 0
0088         if(TARG.eq.'PI+' ) IZT= 1
0089         if(TARG.eq.'PI-' ) IZT=-1
0090         if(TARG.eq.'K+'  ) IZT= 1
0091         if(TARG.eq.'K-'  ) IZT=-1
0092       else
0093         write(6,*)
0094         write(6,*)'Enter mass number and charge of the target nucleus'
0095         read(5,*)  IAT, IZT
0096       endif
0097
0098       write(6,*)
0099       write(6,*)'Enter number of events'
0100       read(5,*) N_events
0101
0102       write(6,*)'Enter FILENAME for HBOOK output'
0103       read(5,1003) FNAME
0104 1003  format(A12)
0105
0106 c------------------------------------------------------------------
0107       CALL HLIMIT(NWPAWC)
0108
0109       EFRMl=EFRM*(9./10.)
0110       EFRMh=EFRM*(11./10.)
0111       CALL HBOOK1(10,' Energy conservation',100,EFRMl,EFRMh,0.)
0112
0113       Pxl=-5.
0114       Pxh= 5.
0115       CALL HBOOK1(20,' Px momentum conservation',100,Pxl,Pxh,0.)
0116
0117       Pyl=-5.
0118       Pyh= 5.
0119       CALL HBOOK1(30,' Py momentum conservation',100,Pyl,Pyh,0.)
0120
0121       if(FRAME.eq.'CMS') then
0122         Pzl=-5.
0123         Pzh= 5.
0124       else
0125         if(PROJ.eq.'P'.or.PROJ.eq.'PBAR'.or.
0126      &     PROJ.eq.'N'.or.PROJ.eq.'NBAR'.or.PROJ.eq.'A') then
0127          Plab=sqrt(EFRM**2-0.88)
0128         endif
0129
0130         if(PROJ.eq.'PI+'.or.PROJ.eq.'PI-') then
0131          Plab=sqrt(EFRM**2-0.0196)
0132         endif
0133
0134         if(PROJ.eq.'K+'.or.PROJ.eq.'K-') then
0135          Plab=sqrt(EFRM-0.25)
0136         endif
0137
0138         Pzl=Plab*(3./4.)
0139         Pzh=Plab*(5./4.)
0140       endif
0141
0142       CALL HBOOK1(40,' Pz momentum conservation',100,Pzl,Pzh,0.)
0143
0144       CALL HBOOK1(50,' Charge conservation'       ,100,-5.,5.,0.)
0145
0146       CALL HBOOK1(60,' Baryon number conservation',100,-5.,5.,0.)
0147
0148       CALL HBOOK1(70,' Lepton number conservation',100,-5.,5.,0.)
0149
0150       CALL HBOOK1(80,' Azimuthal isotropy'      ,100,-3.2,3.2,0.)
0151 C----------------------------------------------------------------
0152
0153       CALL HIJSET(EFRM,FRAME,PROJ,TARG,IAP,IZP,IAT,IZT)
0154 C                       ********Initialize HIJING
0155
0156       WRITE(6,*)' Simulation of interactions with'
0157       WRITE(6,*)
0158       WRITE(6,*)' Proj = ',PROJ,' and Targ = ',TARG
0159       WRITE(6,*)' IAP  =',IAP  ,'            IAT  =',IAT
0160       WRITE(6,*)' IZP  =',IZP  ,'            IZT  =',IZT
0161       WRITE(6,*)
0162       WRITE(6,*)' Reference frame -   ',FRAME
0163       WRITE(6,*)' ENERGY            ',EFRM,' GeV'
0164       WRITE(6,*)' Number of generated events -',N_events
0165       WRITE(6,*)
0166
0167
0168
0169       BMIN=0.0
0170       BMAX=0.0
0171
0172       DO 2000 I_event=1,N_events
0173 C
0174         WRITE(6,*)' Event # ',I_event,' ------------------'
0175 C
0176         CALL HIJING(FRAME,BMIN,BMAX)
0177 C
0178         Sum_E =0.
0179         Sum_Px=0.
0180         Sum_Py=0.
0181         Sum_Pz=0.
0182         Sum_Q =0.
0183         Sum_B =0.
0184         Sum_L =0.
0185
0186         DO 3000 I=1,NATT
0187
0188           if(KATT(I,2).eq. 0) go to 3000 ! reject non-interacting projectile
0189           if(KATT(I,2).eq. 1) go to 3000 ! reject elastic scattering
0190           if(KATT(I,2).eq.10) go to 3000 ! reject non-interacting target
0191           if(KATT(I,2).eq.11) go to 3000 ! reject elastic scattering
0192
0193           ICH=LUCHGE(KATT(I,1))/3
0194
0195           Amass=ULMASS(KATT(I,1))
0196
0197           Sum_E =Sum_E + PATT(i,4)
0198           Sum_Px=Sum_Px+ PATT(i,1)
0199           Sum_Py=Sum_Py+ PATT(i,2)
0200           Sum_Pz=Sum_Pz+ PATT(i,3)
0201           Sum_Q =Sum_Q + ICH
0202
0203           if(IABS(KATT(i,1)).gt.2000) then
0204             Sum_B=Sum_B+IABS(KATT(i,1))/KATT(i,1)
0205           endif
0206
0207           if((IABS(KATT(i,1)).gt.10.).and.
0208      &       (IABS(KATT(i,1)).lt.20.))           then
0209             Sum_L=Sum_L+IABS(KATT(i,1))/KATT(i,1)
0210           endif
0211
0212           Phi=ULANGL(PATT(i,1),PATT(i,2))
0213           CALL HF1(80,Phi,1.)
0214
0215 3000    CONTINUE
0216
0217         CALL HF1(10,Sum_E ,1.)
0218         CALL HF1(20,Sum_Px,1.)
0219         CALL HF1(30,Sum_Py,1.)
0220         CALL HF1(40,Sum_Pz,1.)
0221         CALL HF1(50,Sum_Q ,1.)
0222         CALL HF1(60,Sum_B ,1.)
0223         CALL HF1(70,Sum_L ,1.)
0224 2000  CONTINUE
0225
0226       CALL HRPUT(0,FNAME,'N')
0227       END
0228
0229       FUNCTION RAN(NSEED)
0230       RAN=RLU(NSEED)
0231       RETURN
0232       END