example_5.f

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

 
C****************************************************************************
C The program for checking conservation laws fulfillment in 
C                   the HIJING model.
C           Written by V.Uzhinsky, CERN, Oct. 2003
C***************************************************************************

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

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

      common/himain2/katt(130000,4),patt(130000,4)
      SAVE  /himain2/

C         ********information of produced particles
C
      CHARACTER *1 key
      CHARACTER *12 fname

      parameter(nwpawc=50000)
      common/pawc/hmemory(nwpawc)

      common/ranseed/nseed                                    
      SAVE  /ranseed/                                         
      nseed=0                                                 

      write(6,*)'====================================================='
      write(6,*)'The program for checking conservation laws fulfillment'
      write(6,*)'              in the HIJING model                    '
      write(6,*)'  Only hadronic (hh) interactions are considered     '
      write(6,*)'====================================================='
      write(6,*)'        You can work in Lab. or CM systems           '
      write(6,*)' Would you like to use CM system? Y - Yes, N - No?   '
    
      read(5,1001) key
1001  format(a1)      
      
      if(key.eq.'Y'.or.key.eq.'y') then
        write(6,*)'CM system is used --------------------------------'
        frame='CMS'
      else
        write(6,*)'LAB system is used -------------------------------'      
        frame='LAB'
      endif
      
      write(6,*)'Enter the corresponding energy (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'      

      read(5,1002) proj
1002  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,1002) 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_events                                             

      write(6,*)'Enter FILENAME for HBOOK output'
      read(5,1003) fname
1003  format(a12)

c------------------------------------------------------------------
      CALL hlimit(nwpawc)

      efrml=efrm*(9./10.)
      efrmh=efrm*(11./10.)     
      CALL hbook1(10,' Energy conservation',100,efrml,efrmh,0.)

      pxl=-5.
      pxh= 5.
      CALL hbook1(20,' Px momentum conservation',100,pxl,pxh,0.)

      pyl=-5.         
      pyh= 5.         
      CALL hbook1(30,' Py momentum conservation',100,pyl,pyh,0.)

      if(frame.eq.'CMS') then
        pzl=-5.
        pzh= 5.
      else
        if(proj.eq.'P'.or.proj.eq.'PBAR'.or.
     &     proj.eq.'N'.or.proj.eq.'NBAR'.or.proj.eq.'A') then
         plab=sqrt(efrm**2-0.88)
        endif

        if(proj.eq.'PI+'.or.proj.eq.'PI-') then
         plab=sqrt(efrm**2-0.0196)
        endif

        if(proj.eq.'K+'.or.proj.eq.'K-') then
         plab=sqrt(efrm-0.25)
        endif

        pzl=plab*(3./4.)
        pzh=plab*(5./4.)
      endif
 
      CALL hbook1(40,' Pz momentum conservation',100,pzl,pzh,0.)

      CALL hbook1(50,' Charge conservation'       ,100,-5.,5.,0.)

      CALL hbook1(60,' Baryon number conservation',100,-5.,5.,0.)

      CALL hbook1(70,' Lepton number conservation',100,-5.,5.,0.)

      CALL hbook1(80,' Azimuthal isotropy'      ,100,-3.2,3.2,0.)
C----------------------------------------------------------------

      CALL hijset(efrm,frame,proj,targ,iap,izp,iat,izt)
C     ********Initialize HIJING

      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_events      
      WRITE(6,*)                                              

        

      bmin=0.0
      bmax=0.0

      DO 2000 i_event=1,n_events
C
        WRITE(6,*)' Event # ',i_event,' ------------------'
C
        CALL hijing(frame,bmin,bmax)
C
        sum_e =0.
        sum_px=0.
        sum_py=0. 
        sum_pz=0.
        sum_q =0.
        sum_b =0.
        sum_l =0. 
      
        DO 3000 i=1,natt
                                                               
          if(katt(i,2).eq. 0) go to 3000 ! reject non-interacting projectile
          if(katt(i,2).eq. 1) go to 3000 ! reject elastic scattering
          if(katt(i,2).eq.10) go to 3000 ! reject non-interacting target
          if(katt(i,2).eq.11) go to 3000 ! reject elastic scattering

          ich=luchge(katt(i,1))/3                              

          amass=ulmass(katt(i,1))                             

          sum_e =sum_e + patt(i,4)
          sum_px=sum_px+ patt(i,1)
          sum_py=sum_py+ patt(i,2)
          sum_pz=sum_pz+ patt(i,3)
          sum_q =sum_q + ich

          if(iabs(katt(i,1)).gt.2000) then
            sum_b=sum_b+iabs(katt(i,1))/katt(i,1)
          endif

          if((iabs(katt(i,1)).gt.10.).and.
     &       (iabs(katt(i,1)).lt.20.))           then
            sum_l=sum_l+iabs(katt(i,1))/katt(i,1)
          endif

          phi=ulangl(patt(i,1),patt(i,2))
          CALL hf1(80,phi,1.)

3000    CONTINUE

        CALL hf1(10,sum_e ,1.) 
        CALL hf1(20,sum_px,1.)
        CALL hf1(30,sum_py,1.)
        CALL hf1(40,sum_pz,1.)
        CALL hf1(50,sum_q ,1.)
        CALL hf1(60,sum_b ,1.)
        CALL hf1(70,sum_l ,1.)
2000  CONTINUE

      CALL hrput(0,fname,'N')
      END

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