hijing/examples/example_6.f

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