example_3.f

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

C****************************************************************************
C          Program # 2 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

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

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

      common/himain2/katt(130000,4),patt(130000,4)
      SAVE  /himain2/
C
      dimension gb(101), xb(101), dndp(50)

      common/ranseed/nseed                                     
      SAVE  /ranseed/                                          

      nseed=0                                                 

      do i=1,50
         dndp(i)=0.
      enddo

      frame='CMS'
      
      write(6,*)'===================================================='
      write(6,*)'  Calculation of transverse momentum distribution   '
      write(6,*)'        of charged pions in minimum bias            '
      write(6,*)'            hh-, hA- and AA-collisions              '
      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 per each value of impact',
     ,' parameter (e.a. 10)'

      read(5,*)  n_event

C....initialize HIJING for requested interactions
      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 per B interval -',n_event       
      write(6,*)                                               

C....set BMIN=0 and BMAX=R_A+R_B

      bmin=0.0
      bmax=hipr1(34)+hipr1(35)

C....calculate the Glauber probability and its integrated value:

      dip=(bmax-bmin)/100.0
      gbtot=0.0
      DO 100 i=1,101
         xb(i)=bmin+(i-1)*dip
         ov=profile(xb(i))
         gb(i)=xb(i)*(1.0-exp(-hint1(12)*ov))
         gbtot=gbtot+gb(i)
100   CONTINUE
      write(6,*)'Inelastic X-section (mb) ',gbtot*dip*10.*6.28
  
      pause
C....generating N_EVENT for each of 100 impact parameter intervals:

      nont=0
      gnorm=gbtot
      
      DO 300 ib=1,100
         b1=xb(ib)
         b2=xb(ib+1)

C....normalized Glauber probability:

         w_gb=(gb(ib)+gb(ib+1))/2.0/gbtot

         DO 200 ie=1,n_event
            CALL hijing(frame,b1,b2)

C....count number of events without any interaction
C....and renormalize the total Glauber prabability:

            IF(np+nt .EQ. 0) THEN
              nont=nont+1
              gnorm=gnorm-gb(ib)/float(n_event)
              go to 200
            ENDIF

C....calculate pt distribution of charged pions:

            DO 150 k=1,natt
C....select charged pions only:
               IF(abs(katt(k,1)) .NE. 211) go to 150

C....calculate pt:
               ptr=sqrt(patt(k,1)**2+patt(k,2)**2)

C....calculate pt distribution and weight with normalized
C....Glauber probability to get minimum bias results:

               IF(ptr .GE. 10.0) go to 150
               ipt=1+ptr/0.2
               dndp(ipt)=dndp(ipt)+w_gb/float(n_event)/0.2
150         CONTINUE
200      CONTINUE
300   CONTINUE   

C....renormalize the distribution by the renormalized Glauber
C....probability which excludes the events without any interaction:

      IF(nont.NE.0) THEN
        DO 400 i=1,50
           dndp(i)=dndp(i)*gbtot/gnorm
           if(dndp(i).ne.0.) write(6,350)0.2*(i-1),dndp(i)
350        format(1x,f5.1,2x,e11.4)
400     CONTINUE
      ENDIF

      stop
      END    

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