Correct npart calculation
authormorsch <morsch@f7af4fe6-9843-0410-8265-dc069ae4e863>
Wed, 12 Sep 2012 08:02:48 +0000 (08:02 +0000)
committermorsch <morsch@f7af4fe6-9843-0410-8265-dc069ae4e863>
Wed, 12 Sep 2012 08:02:48 +0000 (08:02 +0000)
C. Loizides

HIJING/hijing1_36/hijing.F
HIJING/hijing1_36/hijset.F

index 7da5a65..89ad5c9 100644 (file)
@@ -379,6 +379,88 @@ C                  ********perform elastic collisions
           IPCOL(NCOLT)=JP
           ITCOL(NCOLT)=JT
 70     CONTINUE
+
+#if 1
+c *** cl glauber ***
+       IF(NCOLT.EQ.0) THEN 
+           goto 60
+        end if
+
+        npart=0
+        xmeana=0D0
+        ymeana=0D0
+        xmeanb=0D0
+        ymeanb=0D0
+        xmeanp=0D0
+        ymeanp=0D0
+        xm2=0D0
+        ym2=0D0
+        xym=0D0
+
+       DO 1110 JP=1,IHNT2(1)
+           YP(1,JP)=YP(1,JP)+BBX
+           YP(2,JP)=YP(2,JP)+BBY
+           xmeana=xmeana+YP(1,JP)
+           ymeana=ymeana+YP(2,JP)
+           DO 1120 JT=1,IHNT2(3)
+              IF(SCIP(JP,JT).GT.-1.0D0) THEN
+                 npart=npart+1
+                 xmeanp=xmeanp+YP(1,JP)
+                 ymeanp=ymeanp+YP(2,JP)
+                 xm2=xm2+YP(1,JP)*YP(1,JP)
+                 ym2=ym2+YP(2,JP)*YP(2,JP)
+                 xym=xym+YP(1,JP)*YP(2,JP)
+                 goto 1110
+              end if
+ 1120      continue
+ 1110   continue
+
+        DO 1130 JT=1,IHNT2(3)
+           xmeanb=xmeanb+YT(1,JT)
+           ymeanb=ymeanb+YT(2,JT)
+           DO 1140 JP=1,IHNT2(1)
+              IF(SCIP(JP,JT).GT.-1.0D0) THEN
+                 npart=npart+1
+                 xmeanp=xmeanp+YT(1,JT)
+                 ymeanp=ymeanp+YT(2,JT)
+                 xm2=xm2+YT(1,JT)*YT(1,JT)
+                 ym2=ym2+YT(2,JT)*YT(2,JT)
+                 xym=xym+YT(1,JT)*YT(2,JT)
+                 goto 1130
+              end if
+ 1140      continue
+ 1130   continue
+
+        xmeana=xmeana/IHNT2(1)
+        ymeana=ymeana/IHNT2(1)
+        xmeanb=xmeanb/IHNT2(3)
+        ymeanb=ymeanb/IHNT2(3)
+        xmeanp=xmeanp/npart
+        ymeanp=ymeanp/npart
+        xm2=xm2/npart
+        ym2=ym2/npart
+        xym=xym/npart
+
+        sx2=xm2-xmeanp*xmeanp
+        sy2=ym2-ymeanp*ymeanp
+        sxy=xym-xmeanp*ymeanp
+
+        delx=xmeanb-xmeana
+        dely=ymeanb-ymeana
+        dtmp=delx**2+dely**2
+        bbtrue=sqrt(dtmp)
+        dnumt=(sy2-sx2)*(delx**2-dely**2)-4D0*sxy*delx*dely
+        ddent=(sy2+sx2)*bbtrue**2
+        eccrp=dnumt/ddent
+        dtmp=(sy2-sx2)*(sy2-sx2)+4D0*sxy*sxy
+        eccpart=sqrt(dtmp)/(sx2+sy2)
+        eccmc=(sy2-sx2)/(sy2+sx2)
+        write(*,*),'HOUT: ',bb,' ',bbtrue,' ',ncolt,' ',npart,' ',eccrp,
+     1  ' ',eccpart
+        RETURN
+#endif
+
+
 C              ********total number interactions proj and targ has
 C                              suffered
        IF(NCOLT.EQ.0) THEN
@@ -621,10 +703,10 @@ C
                IF(IHPR2(21).EQ.0) THEN
                   CALL LUEDIT_HIJING(2)
                ELSE
-351               N_ST=N_ST+1
+351                N_ST=N_ST+1
                   IF(K(N_ST,2).LT.91.OR.K(N_ST,2).GT.93) GO TO  351
-                  IDSTR=K(N_ST,2)
-                  N_ST=N_ST+1
+                   IDSTR=K(N_ST,2)
+                   N_ST=N_ST+1
                ENDIF
 C
                IF(FRAME.EQ.'LAB') THEN
index 0c2d73a..38c8651 100644 (file)
@@ -215,5 +215,10 @@ C          ********booking for x distribution of valence quarks
      & 10X,'*',8X,'for ',
      & A4,'(',I3,',',I3,')',' + ',A4,'(',I3,',',I3,')',7X,'*'/
      & 10X,'**************************************************')
+
+       write (*,*),'Average jet cross section (mb): " ', HINT1(10)
+       write (*,*),'Jet cross section (mb): " ', HINT1(11)
+       write (*,*),'Inelastic NN cross section (mb): ', HINT1(12)
+        write (*,*),'Total NN cross section (mb): ',HINT1(13) 
        RETURN
        END