hijing/examples/example_7.f

0001
0002 C****************************************************************************
0003 C     The program for calculation of the general properties of
0004 C    neutral particles in the interactions by 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       COMMON/HIMAIN2/KATT(130000,4),PATT(130000,4)
0012
0013 C         ********information of produced particles
0014
0015       COMMON/HIPARNT/HIPR1(100),IHPR2(50),HINT1(100),IHNT2(50)
0016       SAVE  /HIPARNT/
0017 C
0018 C         ********Event Options and parameters
0019
0020       CHARACTER *1 KEY
0021       CHARACTER *12 FNAME
0022
0023       PARAMETER (NWPAWC=1000000)
0024       COMMON/PAWC/HMEMORY(NWPAWC)
0025
0026       COMMON/RANSEED/NSEED
0027       SAVE  /RANSEED/
0028       NSEED=0
0029
0030       write(6,*)'====================================================='
0031       write(6,*)'The program for calculation of general properties of '
0032       write(6,*)'   neutral particles (Pi_0, Lambda0, K0_S, Phi       '
0033       write(6,*)'             in the HIJING model.                    '
0034       write(6,*)'====================================================='
0035       write(6,*)
0036       write(6,*)'        You can work in Lab. or CM systems           '
0037       write(6,*)' Would you like to use CM system? Y - Yes, N - No?   '
0038
0039       read(5,1001) KEY
0040 1001  format(A1)
0041
0042       if(KEY.eq.'Y'.or.KEY.eq.'y') then
0043         write(6,*)'CM system is used --------------------------------'
0044         FRAME='CMS'
0045       else
0046         write(6,*)'LAB system is used -------------------------------'
0047         FRAME='LAB'
0048       endif
0049
0050       write(6,*)'Enter the corresponding energy per NN collisions (GeV)'
0051
0052       read(5,*) EFRM
0053
0054       write(6,*)
0055       write(6,*)'Enter a type of the projectile particle'
0056       write(6,*)
0057       write(6,*)' P proton,            PBAR anti-proton,'
0058       write(6,*)' N neutron,           NBAR anti-neutron,'
0059       write(6,*)' PI+ - positive pion, PI- negative pion,'
0060       write(6,*)' K+ positive kaon,    K- negative kaon'
0061       write(6,*)' A - nucleus --------------------------'
0062
0063       read(5,1002) PROJ
0064 1002  format(A8)
0065
0066       if(PROJ.ne.'A') then
0067         IAP=1
0068         if(PROJ.eq.'P'   ) IZP= 1
0069         if(PROJ.eq.'PBAR') IZP=-1
0070         if(PROJ.eq.'N'   ) IZP= 0
0071         if(PROJ.eq.'NBAR') IZP= 0
0072         if(PROJ.eq.'PI+' ) IZP= 1
0073         if(PROJ.eq.'PI-' ) IZP=-1
0074         if(PROJ.eq.'K+'  ) IZP= 1
0075         if(PROJ.eq.'K-'  ) IZP=-1
0076       else
0077         write(6,*)
0078         write(6,*)'Enter mass number and charge of the proj. nucleus'
0079         read(5,*)  IAP, IZP
0080       endif
0081
0082       write(6,*)
0083       write(6,*)'Enter a type of the target particle (same notations)'
0084       read(5,1002) TARG
0085
0086       if(TARG.ne.'A') then
0087         IAT=1
0088         if(TARG.eq.'P'   ) IZT= 1
0089         if(TARG.eq.'PBAR') IZT=-1
0090         if(TARG.eq.'N'   ) IZT= 0
0091         if(TARG.eq.'NBAR') IZT= 0
0092         if(TARG.eq.'PI+' ) IZT= 1
0093         if(TARG.eq.'PI-' ) IZT=-1
0094         if(TARG.eq.'K+'  ) IZT= 1
0095         if(TARG.eq.'K-'  ) IZT=-1
0096       else
0097         write(6,*)
0098         write(6,*)'Enter mass number and charge of the target nucleus'
0099         read(5,*)  IAT, IZT
0100       endif
0101
0102       write(6,*)
0103       write(6,*)'Enter number of events'
0104       read(5,*) N_events
0105
0106       write(6,*)'Enter FILENAME for HBOOK output'
0107       read(5,1003) FNAME
0108 1003  format(A12)
0109
0110 c------------------------------------------------------------------
0111       CALL HIJSET(EFRM,FRAME,PROJ,TARG,IAP,IZP,IAT,IZT)
0112 C                       ********Initialize HIJING
0113
0114       WRITE(6,*)' Simulation of interactions with'
0115       WRITE(6,*)
0116       WRITE(6,*)' Proj = ',PROJ,' and  Targ = ',TARG
0117       WRITE(6,*)' IAP  =',IAP  ,'            IAT  =',IAT
0118       WRITE(6,*)' IZP  =',IZP  ,'            IZT  =',IZT
0119       WRITE(6,*)
0120       WRITE(6,*)' Reference frame -   ',FRAME
0121       WRITE(6,*)' ENERGY            ',EFRM,' GeV'
0122       WRITE(6,*)' Number of generated events -',N_events
0123       WRITE(6,*)
0124
0125       BMIN=0.
0126       BMAX=HIPR1(34)+HIPR1(35)
0127
0128       WRITE(6,*)
0129
0130 c------------------------------------------------------------------
0131       CALL HLIMIT(NWPAWC)
0132 c------------------------------------------------------------------
0133 c------------------------ Pi0 mesons ------------------------------
0134 c------------------------------------------------------------------
0135
0136       Nu_max=MIN(IAP*IAT+1,2000)
0137
0138       if(PROJ.eq.'A'.or.TARG.eq.'A') then
0139         CH_max=5*Nu_max/2.
0140       else
0141         CH_max=99.5
0142       endif
0143
0144       CALL HBOOK1(10,' Pi0 meson multiplicity distribution',
0145      ,100,-0.5,CH_max,0.)
0146
0147       if(FRAME.eq.'CMS') then
0148         Eta_l=-15.
0149         Eta_h=+15.
0150       else
0151         Eta_l=-2.0
0152         Eta_h=Alog(2.*EFRM)+2.
0153       endif
0154       NbinEta=(Eta_h-Eta_l)/0.2
0155
0156       CALL HBOOK1(20,' Pi0 meson pseudo-rapidity distribution',
0157      ,NbinEta,Eta_l,Eta_h,0.)
0158
0159       if(FRAME.eq.'CMS') then
0160         Y_l=-Alog(EFRM)-2.
0161         Y_h=+Alog(EFRM)+2.
0162       else
0163         Y_l=-2.0
0164         Y_h=Alog(2.*EFRM)+2.
0165       endif
0166       NbinY=(Y_h-Y_l)/0.2
0167
0168       CALL HBOOK1(30,' Pi0 meson rapidity distribution',
0169      ,NbinY,Y_l,Y_h,0.)
0170
0171       CALL HBOOK1(40,' Pi0 meson Pt distribution',
0172      ,100,0.,10.,0.)
0173
0174       if(FRAME.eq.'CMS') then
0175         E_l= 0.
0176         E_h=EFRM/2.
0177       else
0178         E_l= 0.
0179         E_h=EFRM
0180       endif
0181       NbinE=(E_h-E_l)/0.5
0182
0183       CALL HBOOK1(50,' Pi0 meson energy distribution',
0184      ,NbinE,E_l,E_h,0.)
0185
0186       CALL HBOOK1(60,' Pi0 meson Cos(Theta) distribution',
0187      ,40,-1.,1.,0.)
0188
0189       CALL HBOOK1(70,' Pi0 meson Phi distribution',
0190      ,180,0.,180.,0.)
0191
0192       CALL HBOOK2(80,' Phi - Eta correlation of Pi0 meson',
0193      ,NbinEta,Eta_l,Eta_h,180,-180.,180.,0.)
0194
0195 C----------------------------------------------------------------
0196 c--------------------------- Lambda0 ---------------------------
0197 c----------------------------------------------------------------
0198       CALL HBOOK1(110,' Lambda0 multiplicity distribution',
0199      ,100,-0.5,float(IAP+IAT),0.)
0200
0201       CALL HBOOK1(120,' Lambda0 pseudo-rapidity distribution',
0202      ,NbinEta,Eta_l,Eta_h,0.)
0203
0204       CALL HBOOK1(130,' Lambda0 rapidity distribution',
0205      ,NbinY,Y_l,Y_h,0.)
0206
0207       CALL HBOOK1(140,' Lambda0 Pt distribution',
0208      ,100,0.,10.,0.)
0209
0210       CALL HBOOK1(150,' Lambda0 energy distribution',
0211      ,NbinE,E_l,E_h,0.)
0212
0213       CALL HBOOK1(160,' Lambda0 Cos(Theta) distribution',
0214      ,40,-1.,1.,0.)
0215
0216       CALL HBOOK1(170,' Lambda0 Phi distribution',
0217      ,180,0.,180.,0.)
0218
0219       CALL HBOOK2(180,' Phi - Eta correlation of Lambda0',
0220      ,NbinEta,Eta_l,Eta_h,180,-180.,180.,0.)
0221
0222 c----------------------------------------------------------------
0223 c------------------------------- K0_S ---------------------------
0224 c----------------------------------------------------------------
0225       CALL HBOOK1(210,' K0S multiplicity distribution',
0226      ,100,-0.5,CH_max,0.)
0227
0228       CALL HBOOK1(220,' K0S pseudo-rapidity distribution',
0229      ,NbinEta,Eta_l,Eta_h,0.)
0230
0231       CALL HBOOK1(230,' K0S rapidity distribution',
0232      ,NbinY,Y_l,Y_h,0.)
0233
0234       CALL HBOOK1(240,' K0S Pt distribution',
0235      ,100,0.,10.,0.)
0236
0237       CALL HBOOK1(250,' K0S energy distribution',
0238      ,NbinE,E_l,E_h,0.)
0239
0240       CALL HBOOK1(260,' K0S Cos(Theta) distribution',
0241      ,40,-1.,1.,0.)
0242
0243       CALL HBOOK1(270,' K0S Phi distribution',
0244      ,180,0.,180.,0.)
0245
0246       CALL HBOOK2(280,' Phi - Eta correlation of K0S',
0247      ,NbinEta,Eta_l,Eta_h,180,-180.,180.,0.)
0248
0249 c----------------------------------------------------------------
0250 c------------------------- Phi meson  ---------------------------
0251 c----------------------------------------------------------------
0252       CALL HBOOK1(310,' Phi meson multiplicity distribution',
0253      ,100,-0.5,CH_max/3.,0.)
0254
0255       CALL HBOOK1(320,' Phi meson pseudo-rapidity distribution',
0256      ,NbinEta,Eta_l,Eta_h,0.)
0257
0258       CALL HBOOK1(330,' Phi meson rapidity distribution',
0259      ,NbinY,Y_l,Y_h,0.)
0260
0261       CALL HBOOK1(340,' Phi meson Pt distribution',
0262      ,100,0.,10.,0.)
0263
0264       CALL HBOOK1(350,' Phi meson energy distribution',
0265      ,NbinE,E_l,E_h,0.)
0266
0267       CALL HBOOK1(360,' Phi meson Cos(Theta) distribution',
0268      ,40,-1.,1.,0.)
0269
0270       CALL HBOOK1(370,'  Phi meson Phi distribution',
0271      ,180,0.,180.,0.)
0272
0273       CALL HBOOK2(380,' Phi - Eta correlation of Phi meson',
0274      ,NbinEta,Eta_l,Eta_h,180,-180.,180.,0.)
0275
0276 c----------------------------------------------------------------
0277
0278       IHPR2(12)=1                   ! To suppress the particle decays
0279       CALL LUGIVE('MDCY(C333,1)=0') ! To suppress Phi decay
0280 c     CALL LUGIVE('MDCY(C1114,1)=0;MDCY(C-1114,1)=0') ! To suppress Delta- Decay
0281
0282 c----------------------------------------------------------------
0283       Pi=4.*ATAN(1.)                   ! 3.14159
0284
0285       DO 2000 I_event=1,N_events
0286
0287         WRITE(6,*)' Event # ',I_event,' ------------------'
0288 C
0289         CALL HIJING(FRAME,BMIN,BMAX)
0290 C
0291         N_pi0=0
0292         N_lam=0
0293         N_k0s=0
0294         N_phi=0
0295
0296 c        write(6,*)' NATT --------- ',NATT
0297         DO 3000 I=1,NATT
0298
0299 c          write(6,*)i,KATT(I,1)
0300
0301           if(KATT(I,2).eq. 0) go to 3000 ! reject non-interacting projectile
0302           if(KATT(I,2).eq. 1) go to 3000 ! reject elastic scattering
0303           if(KATT(I,2).eq.10) go to 3000 ! reject non-interacting target
0304           if(KATT(I,2).eq.11) go to 3000 ! reject elastic scattering
0305
0306           ID=KATT(i,1)
0307
0308           ICH=LUCHGE(KATT(I,1))/3
0309
0310           Amass=ULMASS(KATT(I,1))
0311
0312           E=PATT(i,4)
0313           Pz=PATT(i,3)
0314           Pt=sqrt(PATT(i,1)**2+PATT(i,2)**2)
0315           P=sqrt(Pz**2+Pt**2)
0316
0317           if(P-Pz.ge.1.0e-4) then
0318             Eta= 0.5*Alog((P+Pz)/(P-Pz))
0319           elseif(P+Pz.ge.1.0e-4) then
0320             Eta=-0.5*Alog((P-Pz)/(P+pz))
0321           else
0322             Eta=15.
0323           endif
0324
0325           if(E-Pz.ge.1.0e-4) then
0326             Y= 0.5*Alog((E+Pz)/(E-Pz))
0327           elseif(E+Pz.ge.1.0e-4) then
0328             Y=-0.5*Alog((E-Pz)/(E+pz))
0329           else
0330             Y=15.
0331           endif
0332
0333
0334           if(P.ge.1.0e-4) then
0335             Cos_Theta=Pz/P
0336           else
0337             Cos_Theta=1.
0338           endif
0339
0340           Phi=ULANGL(PATT(i,1),PATT(i,2))*180./Pi
0341 c==========================================================
0342
0343           if(ID.eq.111) then             ! Pi_0
0344             N_pi0=N_pi0+1
0345
0346             Call HF1(20,Eta,1.)
0347             Call HF1(30,  Y,1.)
0348             Call HF1(40, Pt,1.)
0349             Call HF1(50,  E,1.)
0350             Call HF1(60,Cos_Theta,1.)
0351             Call HF1(70,Phi,1.)
0352             Call HFILL(80,Eta,Phi,1.)
0353             Call Hf1(90,Float(Id),1.)
0354           endif
0355
0356           if(ID.eq.3122) then         ! Lambda0
0357             N_lam=N_lam+1
0358
0359             Call HF1(120,Eta,1.)
0360             Call HF1(130,  Y,1.)
0361             Call HF1(140, Pt,1.)
0362             Call HF1(150,  E,1.)
0363             Call HF1(160,Cos_Theta,1.)
0364             Call HF1(170,Phi,1.)
0365             Call HFILL(180,Eta,Phi,1.)
0366           endif
0367
0368           if(ID.eq.310) then          ! K0_S
0369             N_k0s=N_k0s+1
0370
0371             Call HF1(220,Eta,1.)
0372             Call HF1(230,  Y,1.)
0373             Call HF1(240, Pt,1.)
0374             Call HF1(250,  E,1.)
0375             Call HF1(260,Cos_Theta,1.)
0376             Call HF1(270,Phi,1.)
0377             Call HFILL(280,Eta,Phi,1.)
0378           endif
0379
0380           if(Id.eq.333) then          ! Phi
0381             N_phi=N_phi+1
0382
0383             Call HF1(320,Eta,1.)
0384             Call HF1(330,  Y,1.)
0385             Call HF1(340, Pt,1.)
0386             Call HF1(350,  E,1.)
0387             Call HF1(360,Cos_Theta,1.)
0388             Call HF1(370,Phi,1.)
0389             Call HFILL(380,Eta,Phi,1.)
0390           endif
0391
0392 3000    CONTINUE
0393
0394 c       PAUSE
0395         CALL HF1( 10,float(N_pi0),1.)
0396         CALL HF1(110,float(N_lam),1.)
0397         CALL HF1(210,float(N_k0s),1.)
0398         CALL HF1(310,float(N_phi),1.)
0399
0400 2000  CONTINUE
0401
0402 c================ Normalization ===========================
0403
0404       C1=1./float(N_events)                        ! Multiplicity distr.
0405       C2=0.
0406
0407       Call HOPERA( 10,'+', 10, 10,C1,C2)
0408       Call HOPERA(110,'+',110,110,C1,C2)
0409       Call HOPERA(210,'+',210,210,C1,C2)
0410       Call HOPERA(310,'+',310,310,C1,C2)
0411
0412       C1=1./float(N_events)/((Eta_h-Eta_l)/NbinEta)! Eta distr.
0413       C2=0.
0414
0415       Call HOPERA( 20,'+', 20, 20,C1,C2)
0416       Call HOPERA(120,'+',120,120,C1,C2)
0417       Call HOPERA(220,'+',220,220,C1,C2)
0418       Call HOPERA(320,'+',320,320,C1,C2)
0419
0420       C1=1./float(N_events)/((Y_h-Y_l)/NbinY)      ! Y distr.
0421       C2=0.
0422
0423       Call HOPERA( 30,'+', 30, 30,C1,C2)
0424       Call HOPERA(130,'+',130,130,C1,C2)
0425       Call HOPERA(230,'+',230,230,C1,C2)
0426       Call HOPERA(330,'+',330,330,C1,C2)
0427
0428       C1=1./float(N_events)/0.1                    ! Pt distr.
0429       C2=0.
0430
0431       Call HOPERA( 40,'+', 40, 40,C1,C2)
0432       Call HOPERA(140,'+',140,140,C1,C2)
0433       Call HOPERA(240,'+',240,240,C1,C2)
0434       Call HOPERA(340,'+',340,340,C1,C2)
0435
0436       C1=1./float(N_events)/((E_h-E_l)/NbinE)      ! E distr.
0437       C2=0.
0438
0439       Call HOPERA( 50,'+', 50, 50,C1,C2)
0440       Call HOPERA(150,'+',150,150,C1,C2)
0441       Call HOPERA(250,'+',250,250,C1,C2)
0442       Call HOPERA(350,'+',350,350,C1,C2)
0443
0444       C1=1./float(N_events)/0.05                   ! Cos Theta distr.
0445       C2=0.
0446
0447       Call HOPERA( 60,'+', 60, 60,C1,C2)
0448       Call HOPERA(160,'+',160,160,C1,C2)
0449       Call HOPERA(260,'+',260,260,C1,C2)
0450       Call HOPERA(360,'+',360,360,C1,C2)
0451
0452       C1=1./float(N_events)                        ! Phi distr.
0453       C2=0.
0454
0455       Call HOPERA( 70,'+', 70, 70,C1,C2)
0456       Call HOPERA(170,'+',170,170,C1,C2)
0457       Call HOPERA(270,'+',270,270,C1,C2)
0458       Call HOPERA(370,'+',370,370,C1,C2)
0459
0460 c================ Writing results =========================
0461
0462       CALL HRPUT(0,FNAME,'N')
0463       END
0464
0465       FUNCTION RAN(NSEED)
0466       RAN=RLU(NSEED)
0467       RETURN
0468       END