hijing/examples/example_8.f

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