example_2.f

Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file example_2.f

C****************************************************************************
C          Program # 1 from Comp. Phys. Commun. 83 (1994) 307
C            by M. Gyulassy and X-.N. Wang
C             Modified by V.Uzhinsky, CERN, Oct. 2003 
C***************************************************************************

      CHARACTER frame*8,proj*8,targ*8

      dimension dndpt(50), dndy(50)

      common/hiparnt/hipr1(100),ihpr2(50),hint1(100),ihnt2(50)
      SAVE  /hiparnt/

C....information of produced particles:

      common/himain1/natt,eatt,jatt,nt,np,n0,n01,n10,n11
      SAVE  /himain1/

      common/himain2/katt(130000,4),patt(130000,4)
      SAVE  /himain2/
C
C....information of produced partons:

      common/hijjet1/npj(300),kfpj(300,500),pjpx(300,500),pjpy(300,500)
     &         ,pjpz(300,500),pjpe(300,500),pjpm(300,500)
     &         ,ntj(300),kftj(300,500),pjtx(300,500),pjty(300,500)
     &         ,pjtz(300,500),pjte(300,500),pjtm(300,500)
      SAVE  /hijjet1/

      common/hijjet2/nsg,njsg(900),iasg(900,3),k1sg(900,100)
     &         ,k2sg(900,100),pxsg(900,100),pysg(900,100),pzsg(900,100)
     &         ,pesg(900,100),pmsg(900,100)
      SAVE  /hijjet2/
C 
      common/histrng/nfp(300,15),pp(300,15),nft(300,15),pt(300,15)
      SAVE  /histrng/

      common/ranseed/nseed                                     
      SAVE  /ranseed/                                          

      nseed=0                                                 

      do k=1,50
        dndpt(k)=0.
        dndy(k) =0.
      enddo

      frame='CMS'
      
      write(6,*)'===================================================='
      write(6,*)'  Calculation of transverse momentum and rapidity   '
      write(6,*)'       distributions of charged particles in        '
      write(6,*)'hh-, hA- and AA-collisions at fixed impact parameter'
      write(6,*)'    Calculation will be performed in CM System      '
      write(6,*)'===================================================='

      write(6,*)
      write(6,*)'Enter the energy per NN-collision (GeV)'
      read(5,*)  efrm

      write(6,*)   
      write(6,*)'Enter a type of the "projectile" particle'
      write(6,*)
      write(6,*)' P proton,            PBAR anti-proton,'
      write(6,*)' N neutron,           NBAR anti-neutron,'
      write(6,*)' PI+ - positive pion, PI- negative pion,'
      write(6,*)' K+ positive kaon,    K- negative kaon'
      write(6,*)
      write(6,*)' A - nucleus --------------------------'

      read(5,1) proj
1     format(a8)
    
      if(proj.ne.'A') then
        iap=1
        if(proj.eq.'P'   ) izp= 1
        if(proj.eq.'PBAR') izp=-1
        if(proj.eq.'N'   ) izp= 0
        if(proj.eq.'NBAR') izp= 0
        if(proj.eq.'PI+' ) izp= 1
        if(proj.eq.'PI-' ) izp=-1
        if(proj.eq.'K+'  ) izp= 1
        if(proj.eq.'K-'  ) izp=-1
      else
        write(6,*)
        write(6,*)'Enter mass number and charge of the proj. nucleus'
        read(5,*)  iap, izp
      endif

      write(6,*)
      write(6,*)'Enter a type of the "target" particle (same notations)'
      read(5,1) targ

      if(targ.ne.'A') then
        iat=1
        if(targ.eq.'P'   ) izt= 1
        if(targ.eq.'PBAR') izt=-1
        if(targ.eq.'N'   ) izt= 0
        if(targ.eq.'NBAR') izt= 0
        if(targ.eq.'PI+' ) izt= 1
        if(targ.eq.'PI-' ) izt=-1
        if(targ.eq.'K+'  ) izt= 1
        if(targ.eq.'K-'  ) izt=-1
      else
        write(6,*)
        write(6,*)'Enter mass number and charge of the target nucleus'
        read(5,*)  iat, izt
      endif
      
      write(6,*)
      write(6,*)'Enter number of events'

      read(5,*)  n_event

      CALL hijset(efrm,frame,proj,targ,iap,izp,iat,izt)

      write(6,*)' Simulation of interactions with'             
      write(6,*)                                               
      write(6,*)' Proj = ',proj,' and  Targ = ',targ            
      write(6,*)' IAP  =',iap  ,'            IAT  =',iat       
      write(6,*)' IZP  =',izp  ,'            IZT  =',izt       
      write(6,*)                                               
      write(6,*)' Reference frame -   ',frame                  
      write(6,*)' ENERGY            ',efrm,' (GeV)'                     
      write(6,*)' Number of generated events -',n_event       
      write(6,*)                                               

      if(proj.eq.'A'.or.targ.eq.'A') then
        write(6,*)'Enter Min. and Max. values of impact parameter (fm)'
        read(5,*)  bmin, bmax
      else
        bmin=0.0
        bmax=0.0
      endif

      DO 2000 j=1,n_event

        write(6,*)' Event # ',j,' ------------------------------'

        CALL hijing(frame,bmin,bmax)
C
C....calculate rapidity and transverse momentum distributions of
C....produced charged particles:  

        DO 1000 i=1,natt

C....exclude beam nucleons as produced particles:

          if(katt(i,2).EQ.0 .OR. katt(i,2).EQ.10) go to 1000

C....select charged particles only:

          IF(luchge(katt(i,1)) .EQ. 0) go to 1000

          ptr=sqrt(patt(i,1)**2+patt(i,2)**2)
          IF (ptr .GE. 10.0) go to 100

          ipt=1+ptr/0.2
          dndpt(ipt)=dndpt(ipt)+1.0/float(n_event)/0.2/2.0/ptr

100       y=0.5*log((patt(i,4)+patt(i,3))/(patt(i,4)-patt(i,3)))
          IF(abs(y) .GE. 10.0) go to 1000

          iy=1+abs(y)/0.2
          dndy(iy)=dndy(iy)+1.0/float(n_event)/0.2/2.0
1000    CONTINUE
2000  CONTINUE

C....print out the rapidity and transverse momentum distributions:

      do k=1,50
        WRITE(6,2)0.2*(k-1),dndpt(k),dndy(k)
2       format(1x,f5.1,2(2x,e11.4))
      enddo

      stop
      END    

      FUNCTION ran(NSEED)                                      
      ran=rlu(nseed)                                           
      RETURN                                                   
      END