example_4.f

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

C****************************************************************************
C          Program # 3 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/ranseed/nseed                                     
      SAVE  /ranseed/                                          

      CHARACTER *1 key

      nseed=0                                                 

C....switch on triggered jet production:
      ihpr2(3)=1

      frame='CMS'
      
      write(6,*)'===================================================='
      write(6,*)' Simulation of events with triggered hard processes '
      write(6,*)'          in 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
      
      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

      if(proj.eq.'A'.or.targ.eq.'A') then
        write(6,*)
        write(6,*)' Would you like to take jet quenching into account?'
        write(6,*)'                 Y - Yes, N - No?   '

        read(5,2) key
2       format(a1)

        if(key.eq.'Y'.or.key.eq.'y') then
          ihpr2(4)=1 
        else
          ihpr2(4)=0 
        endif
      endif

C....set the pt range of the triggered jets:
      write(6,*)
      write(6,*)' Set the Pt of triggered jets (in GeV/c)'
      read(5,*)   pt_trigger

      write(6,*)
      write(6,*)'Enter number of events'
      read(5,*)  n_events

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',n_events       
      write(6,*)' Triggered jets Pt =>',pt_trigger,' (GeV/c)'
      write(6,*)                                               

      DO 100 i_event=1,n_events

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

        hipr1(10)=-pt_trigger

        CALL hijing(frame,bmin,bmax)

C....print out flavor code of the first jet:
        write(6,*)'Flavor code of the first jet: ',ihnt2(9)

C....and its four momentum:
        write(6,*)'     Px           Py           Pz',
     ,'            E     (GeV/c, GeV)'
        write(6,3)hint1(21),hint1(22),hint1(23),hint1(24)
3       format(4(e11.4,2x))        
        write(6,*)

C....print out flavor code of the second jet:
        write(6,*)'Flavor code of the second jet:',ihnt2(10)

C....and its four momentum:
        write(6,*)'     Px           Py           Pz',
     ,'            E     (GeV/c, GeV)'
        write(6,3)hint1(31),hint1(32),hint1(33),hint1(34)
        write(6,*)

        pause
100   CONTINUE
 
      stop
      END    

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