hijset.f

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

C
C
C
  SUBROUTINE hijset(EFRM,FRAME,PROJ,TARG,IAP,IZP,IAT,IZT)
  CHARACTER frame*8,proj*8,targ*8,eframe*8                 ! Ilya Seluzhenkov 
  DOUBLE PRECISION  dd1,dd2,dd3,dd4
  common/histrng/nfp(300,15),pp(300,15),nft(300,15),pt(300,15)
  SAVE  /histrng/
  common/hijcrdn/yp(3,300),yt(3,300)
  SAVE  /hijcrdn/
  common/hiparnt/hipr1(100),ihpr2(50),hint1(100),ihnt2(50)
  SAVE  /hiparnt/
        common/hijdat/hidat0(10,10),hidat(10)
  SAVE  /hijdat/
        common/ludat1/mstu(200),paru(200),mstj(200),parj(200)
  SAVE  /ludat1/
  EXTERNAL fnkick,fnkick2,fnstru,fnstrum,fnstrus
  CALL title
  ihnt2(1)=iap
  ihnt2(2)=izp
  ihnt2(3)=iat
  ihnt2(4)=izt
  ihnt2(5)=0
  ihnt2(6)=0
C
  hint1(8)=max(ulmass(2112),ulmass(2212))
  hint1(9)=hint1(8)
C
  IF(proj.NE.'A') THEN
    IF(proj.EQ.'P') THEN
        ihnt2(5)=2212
    ELSE IF(proj.EQ.'PBAR') THEN 
        ihnt2(5)=-2212
    ELSE IF(proj.EQ.'PI+') THEN
        ihnt2(5)=211
    ELSE IF(proj.EQ.'PI-') THEN
        ihnt2(5)=-211
    ELSE IF(proj.EQ.'K+') THEN
        ihnt2(5)=321
    ELSE IF(proj.EQ.'K-') THEN
        ihnt2(5)=-321
    ELSE IF(proj.EQ.'N') THEN
        ihnt2(5)=2112
    ELSE IF(proj.EQ.'NBAR') THEN
        ihnt2(5)=-2112
    ELSE
        WRITE(6,*) proj, 'wrong or unavailable proj name'
        stop
    ENDIF
    hint1(8)=ulmass(ihnt2(5))
  ENDIF
  IF(targ.NE.'A') THEN
    IF(targ.EQ.'P') THEN
        ihnt2(6)=2212
    ELSE IF(targ.EQ.'PBAR') THEN 
        ihnt2(6)=-2212
    ELSE IF(targ.EQ.'PI+') THEN
        ihnt2(6)=211
    ELSE IF(targ.EQ.'PI-') THEN
        ihnt2(6)=-211
    ELSE IF(targ.EQ.'K+') THEN
        ihnt2(6)=321
    ELSE IF(targ.EQ.'K-') THEN
        ihnt2(6)=-321
    ELSE IF(targ.EQ.'N') THEN
        ihnt2(6)=2112
    ELSE IF(targ.EQ.'NBAR') THEN
        ihnt2(6)=-2112
    ELSE
        WRITE(6,*) targ,'wrong or unavailable targ name'
        stop
    ENDIF
    hint1(9)=ulmass(ihnt2(6))
  ENDIF

C...Switch off decay of pi0, K0S, Lambda, Sigma+-, Xi0-, Omega-.

  print *, 'IHPR2(12) = ', ihpr2(12)

  IF(ihpr2(12).GT.0) THEN
     CALL lugive('MDCY(C111,1)=0')
     CALL lugive('MDCY(C310,1)=0')
     CALL lugive('MDCY(C3122,1)=0;MDCY(C-3122,1)=0')
     CALL lugive('MDCY(C3112,1)=0;MDCY(C-3112,1)=0')
     CALL lugive('MDCY(C3212,1)=0;MDCY(C-3212,1)=0')
     CALL lugive('MDCY(C3222,1)=0;MDCY(C-3222,1)=0')
     CALL lugive('MDCY(C3312,1)=0;MDCY(C-3312,1)=0')
     CALL lugive('MDCY(C3322,1)=0;MDCY(C-3322,1)=0')
     CALL lugive('MDCY(C3334,1)=0;MDCY(C-3334,1)=0')

     IF (ihpr2(12).EQ.1) THEN
        CALL lugive('MDCY(C411,1)=0;MDCY(C-411,1)=0')
        CALL lugive('MDCY(C421,1)=0;MDCY(C-421,1)=0')
        CALL lugive('MDCY(C431,1)=0;MDCY(C-431,1)=0')
        CALL lugive('MDCY(C511,1)=0;MDCY(C-511,1)=0')
        CALL lugive('MDCY(C521,1)=0;MDCY(C-521,1)=0')
        CALL lugive('MDCY(C531,1)=0;MDCY(C-531,1)=0')
     ENDIF
  ENDIF

  mstu(12)=0
  mstu(21)=1
  IF(ihpr2(10).EQ.0) THEN
    mstu(22)=0
    mstu(25)=0
    mstu(26)=0
  ENDIF
  mstj(12)=ihpr2(11)
  parj(21)=hipr1(2)
  parj(41)=hipr1(3)
  parj(42)=hipr1(4)
C     ******** set up for jetset
  IF(frame.EQ.'LAB') THEN
     dd1=efrm
     dd2=hint1(8)
     dd3=hint1(9)
     hint1(1)=sqrt(hint1(8)**2+2.0*hint1(9)*efrm+hint1(9)**2)
     dd4=dsqrt(dd1**2-dd2**2)/(dd1+dd3)
     hint1(2)=dd4
     hint1(3)=0.5*dlog((1.d0+dd4)/(1.d0-dd4))
     dd4=dsqrt(dd1**2-dd2**2)/dd1
     hint1(4)=0.5*dlog((1.d0+dd4)/(1.d0-dd4))
     hint1(5)=0.0
     hint1(6)=efrm
     hint1(7)=hint1(9)
  ELSE IF(frame.EQ.'CMS') THEN
     hint1(1)=efrm
     hint1(2)=0.0
     hint1(3)=0.0
     dd1=hint1(1)
     dd2=hint1(8)
     dd3=hint1(9)
     dd4=dsqrt(1.d0-4.d0*dd2**2/dd1**2)
     hint1(4)=0.5*dlog((1.d0+dd4)/(1.d0-dd4))
     dd4=dsqrt(1.d0-4.d0*dd3**2/dd1**2)
     hint1(5)=-0.5*dlog((1.d0+dd4)/(1.d0-dd4))
     hint1(6)=hint1(1)/2.0
     hint1(7)=hint1(1)/2.0
  ENDIF
C   ********define Lorentz transform to lab frame
c
C   ********calculate the cross sections involved with
C     nucleon collisions.
  IF(ihnt2(1).GT.1) THEN
    CALL hijwds(ihnt2(1),1,rmax)
    hipr1(34)=rmax
C     ********set up Wood-Sax distr for proj.
  ENDIF
  IF(ihnt2(3).GT.1) THEN
    CALL hijwds(ihnt2(3),2,rmax)
    hipr1(35)=rmax
C     ********set up Wood-Sax distr for  targ.
  ENDIF
C
C
  i=0
20  i=i+1
  IF(i.EQ.10) go to 30
  IF(hidat0(10,i).LE.hint1(1)) go to 20
30  IF(i.EQ.1) i=2
  DO 40 j=1,9
     hidat(j)=hidat0(j,i-1)+(hidat0(j,i)-hidat0(j,i-1))
     &     *(hint1(1)-hidat0(10,i-1))/(hidat0(10,i)-hidat0(10,i-1))
40  CONTINUE
  hipr1(31)=hidat(5)
  hipr1(30)=2.0*hidat(5)
C
C
  CALL hijcrs
C
  IF(ihpr2(5).NE.0) THEN
    CALL hifun(3,0.0,36.0,fnkick)
C   ********booking for generating pt**2 for pt kick
  ENDIF
  CALL hifun(7,0.0,6.0,fnkick2)
  CALL hifun(4,0.0,1.0,fnstru)
  CALL hifun(5,0.0,1.0,fnstrum)
  CALL hifun(6,0.0,1.0,fnstrus)
C   ********booking for x distribution of valence quarks
  eframe='Ecm'
  IF(frame.EQ.'LAB') eframe='Elab'
  WRITE(6,100) eframe,efrm,proj,ihnt2(1),ihnt2(2),
     &               targ,ihnt2(3),ihnt2(4) 
100 FORMAT(//10x,'****************************************
     &  **********'/
     &  10x,'*',48x,'*'/
     &  10x,'*         HIJING has been initialized at         *'/
     &  10x,'*',13x,a4,'= ',f10.2,' GeV/n',13x,'*'/
     &  10x,'*',48x,'*'/
     &  10x,'*',8x,'for ',
     &  a4,'(',i3,',',i3,')',' + ',a4,'(',i3,',',i3,')',7x,'*'/
     &  10x,'**************************************************')
  RETURN
  END