]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - HIJING/hijing1_36/hijing.F
Update from Jiri: enables the task to create a very simple offline L1 jet trigger...
[u/mrichter/AliRoot.git] / HIJING / hijing1_36 / hijing.F
index 875a266cd8ce3bade1b606336a19b250d0d26d80..07db9ceb025dd87e8454cdf3b87c414cd853e03f 100644 (file)
@@ -154,10 +154,12 @@ C****************************************************************
        SUBROUTINE HIJING(FRAME,BMIN0,BMAX0)
        CHARACTER FRAME*8
        DIMENSION SCIP(300,300),RNIP(300,300),SJIP(300,300),JTP(3),
-     &                 IPCOL(90000),ITCOL(90000)
+     &                 IPCOL(90000),ITCOL(90000),SCIP2(300,300)
+#define BLANKET_SAVE
 #include "hiparnt.inc"
 C
 #include "hijcrdn.inc"
+#include "hijglbr.inc"
 #include "himain1.inc"
 #include "himain2.inc"
 #include "histrng.inc"
@@ -185,6 +187,15 @@ C
        IF(IHNT2(1).LE.1) GO TO 14
        DO 10 KP=1,IHNT2(1)
 5      R=HIRND(1)
+c
+        if(IHNT2(1).EQ.2) then
+           rnd1=max(RLU_HIJING(NSEED),1.0e-20)
+           rnd2=max(RLU_HIJING(NSEED),1.0e-20)
+           rnd3=max(RLU_HIJING(NSEED),1.0e-20)
+           R=-0.5*(log(rnd1)*4.38/2.0+log(rnd2)*0.85/2.0
+     &          +4.38*0.85*log(rnd3)/(4.38+0.85))
+        endif
+c
        X=RLU_HIJING(0)
        CX=2.0*X-1.0
        SX=SQRT(1.0-CX*CX)
@@ -205,6 +216,13 @@ C                  ********two neighbors cannot be closer than
 C                              HIPR1(29)
 8      CONTINUE
 10     CONTINUE
+c*******************************
+        if(IHNT2(1).EQ.2) then
+           YP(1,2)=-YP(1,1)
+           YP(2,2)=-YP(2,1)
+           YP(3,2)=-YP(3,1)
+        endif
+c********************************
        DO 12 I=1,IHNT2(1)-1
        DO 12 J=I+1,IHNT2(1)
        IF(YP(3,I).GT.YP(3,J)) GO TO 12
@@ -226,6 +244,15 @@ C******************************
        IF(IHNT2(3).LE.1) GO TO 24
        DO 20 KT=1,IHNT2(3)
 15     R=HIRND(2)
+c
+         if(IHNT2(3).EQ.2) then
+            rnd1=max(RLU_HIJING(NSEED),1.0e-20)
+            rnd2=max(RLU_HIJING(NSEED),1.0e-20)
+            rnd3=max(RLU_HIJING(NSEED),1.0e-20)
+            R=-0.5*(log(rnd1)*4.38/2.0+log(rnd2)*0.85/2.0
+     &           +4.38*0.85*log(rnd3)/(4.38+0.85))
+         endif
+c
        X=RLU_HIJING(0)
        CX=2.0*X-1.0
        SX=SQRT(1.0-CX*CX)
@@ -246,6 +273,13 @@ C                  ********two neighbors cannot be closer than
 C                              HIPR1(29)
 18     CONTINUE
 20     CONTINUE
+c**********************************
+         if(IHNT2(3).EQ.2) then
+            YT(1,2)=-YT(1,1)
+            YT(2,2)=-YT(2,1)
+            YT(3,2)=-YT(3,1)
+         endif
+c*********************************
        DO 22 I=1,IHNT2(3)-1
        DO 22 J=I+1,IHNT2(3)
        IF(YT(3,I).LT.YT(3,J)) GO TO 22
@@ -281,8 +315,16 @@ C                  ********Initialize for a new event
        N01=0
        N10=0
        N11=0
+        NELT=0
+        NINT=0
+        NELP=0
+        NINP=0
        NSG=0
        NCOLT=0
+        NPSPECP=0
+        NNSPECP=0
+        NPSPECT=0
+        NNSPECT=0 
 
 C****  BB IS THE ABSOLUTE VALUE OF IMPACT PARAMETER,BB**2 IS 
 C       RANDOMLY GENERATED AND ITS ORIENTATION IS RANDOMLY SET 
@@ -300,6 +342,7 @@ C
           SCIP(JP,JT)=-1.0
           B2=(YP(1,JP)+BBX-YT(1,JT))**2+(YP(2,JP)+BBY-YT(2,JT))**2
           R2=B2*HIPR1(40)/HIPR1(31)/0.1
+          SCIP2(JP,JT)=R2
 C              ********mb=0.1*fm, YP is in fm,HIPR1(31) is in mb
           RRB1=MIN((YP(1,JP)**2+YP(2,JP)**2)
      &          /1.2**2/REAL(IHNT2(1))**0.6666667,1.0)
@@ -334,15 +377,99 @@ C                 ********perform elastic collisions
           RNIP(JP,JT)=RANTOT
           SJIP(JP,JT)=HINT1(18)
           NCOLT=NCOLT+1
+           if (R2.GT.2.D0) THEN
+              write (8,*) R2
+           ENDIF
           IPCOL(NCOLT)=JP
           ITCOL(NCOLT)=JT
 70     CONTINUE
+
+c *** cl glauber ***
+        npart=0
+        xmeana=0D0
+        ymeana=0D0
+        xmeanb=0D0
+        ymeanb=0D0
+        xmeanp=0D0
+        ymeanp=0D0
+        xm2=0D0
+        ym2=0D0
+        xym=0D0
+
+       IF(NCOLT>0) THEN 
+           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(SCIP2(JP,JT).LT.2.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(SCIP2(JP,JT).LT.2.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
+
+           IF (npart.GT.0) THEN
+              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,
+     1             ' ',eccrp,' ',eccpart
+           end if
+        end if
+
 C              ********total number interactions proj and targ has
 C                              suffered
+
        IF(NCOLT.EQ.0) THEN
           NLOP=NLOP+1
            IF(NLOP.LE.20.OR.
      &           (IHNT2(1).EQ.1.AND.IHNT2(3).EQ.1)) GO TO 60
+
+
            RETURN
        ENDIF
 C               ********At large impact parameter, there maybe no
@@ -498,6 +625,43 @@ C          ********conduct soft scattering between JP and JT
 
 200    CONTINUE
 
+c
+c**************************
+c
+       DO 201 JP=1,IHNT2(1)
+c           write(6,*) JP, NFP(JP,3), NFP(JP,4), NFP(JP,5)
+           IF(NFP(JP,5).GT.2) THEN
+              NINP=NINP+1
+           ELSE IF(NFP(JP,5).EQ.2.OR.NFP(JP,5).EQ.1) THEN
+              NELP=NELP+1
+           ENDIF
+
+           IF(NFP(JP,5).LE.2) THEN
+              IF (NFP(JP,3) .EQ. 2212) THEN
+                 NPSPECP = NPSPECP + 1
+              ELSE IF (NFP(JP,3) .EQ. 2112) THEN
+                 NNSPECP = NNSPECP + 1
+              ENDIF
+           ENDIF
+ 201    continue
+       DO 202 JT=1,IHNT2(3)
+           IF(NFT(JT,5).GT.2) THEN
+              NINT=NINT+1
+           ELSE IF(NFT(JT,5).EQ.2.OR.NFT(JT,5).EQ.1) THEN
+              NELT=NELT+1
+           ENDIF
+
+           IF(NFT(JT,5).LE.2) THEN
+              IF (NFT(JT,3) .EQ. 2212) THEN
+                 NPSPECT = NPSPECT + 1
+              ELSE IF (NFT(JT,3) .EQ. 2112) THEN
+                 NNSPECT = NNSPECT + 1
+              ENDIF
+           ENDIF
+ 202    continue
+c     
+c*******************************
+
 C********perform jet quenching for jets with PT>HIPR1(11)**********
 
        IF((IHPR2(8).NE.0.OR.IHPR2(3).NE.0).AND.IHPR2(4).GT.0.AND.
@@ -520,6 +684,7 @@ C********N_STR+1 is the number of strings in fragmentation
 C********the number of strings before a line is stored in K(I,4)
 C********IDSTR is id number of the string system (91,92 or 93)
 C
+c
         IF(IHPR2(20).NE.0) THEN
           DO 360 ISG=1,NSG
                CALL HIJFRG(ISG,3,IERROR)
@@ -539,10 +704,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
@@ -561,11 +726,15 @@ C
                   KATT(NATT,1)=K(I,2)
                   KATT(NATT,2)=20
                   KATT(NATT,4)=K(I,1)
-                  IF(K(I,3).EQ.0 .OR. K(K(I,3),2).EQ.IDSTR) THEN
-                     KATT(NATT,3)=0
-                  ELSE
-                     KATT(NATT,3)=NATT-I+K(I,3)+N_STR-K(K(I,3),4)
-                  ENDIF
+                   IF(K(I,3).EQ.0 ) THEN
+                      KATT(NATT,3)=0
+                   ELSE
+                      IF(K(K(I,3),2).EQ.IDSTR) THEN
+                         KATT(NATT,3)=0
+                      ELSE
+                         KATT(NATT,3)=NATT-I+K(I,3)+N_STR-K(K(I,3),4)
+                      ENDIF
+                   ENDIF
 C       ****** identify the mother particle
                   PATT(NATT,1)=P(I,1)
                   PATT(NATT,2)=P(I,2)
@@ -575,7 +744,7 @@ C       ****** identify the mother particle
                    VATT(NATT,2)=V(I,2)
                    VATT(NATT,3)=V(I,3)
                    VATT(NATT,4)=V(I,4)
-                   
+
                   EATT=EATT+P(I,4)
 360       CONTINUE
 C              ********Fragment the q-qbar jets systems *****
@@ -598,28 +767,24 @@ C                 ********check errors
 C
                N_ST=1
                IDSTR=92
+
+               NFTP=NFP(J_JTP,5)
+               IF(NTP.EQ.2) NFTP=10+NFT(J_JTP,5)
+
                IF(IHPR2(21).EQ.0) THEN
                   CALL LUEDIT_HIJING(2)
-               ELSE
+               ELSE IF (NFTP.EQ. 3 .OR. NFTP .EQ. 13) THEN
 381               N_ST=N_ST+1
-C     PH Problem if N_ST > 9000, some additional protection added
-C     PH 9000 is the first dimension of LUJETS/K(...)
-           IF (N_ST.GT.9000) THEN
-                     call LULIST_HIJING(1)
-                     WRITE(6,*) 'error occured, repeat the event'
-              GOTO 50
-                  ENDIF             
                   IF(K(N_ST,2).LT.91.OR.K(N_ST,2).GT.93) GO TO  381
                   IDSTR=K(N_ST,2)
                   N_ST=N_ST+1
-               ENDIF
+                ENDIF
                IF(FRAME.EQ.'LAB') THEN
                        CALL HIBOOST
                ENDIF
 C              ******** boost back to lab frame(if it was in)
 C
-               NFTP=NFP(J_JTP,5)
-               IF(NTP.EQ.2) NFTP=10+NFT(J_JTP,5)
+
                N_STR=0
                DO 390 I=N_ST,N
                   IF(K(I,2).EQ.IDSTR) THEN
@@ -631,7 +796,7 @@ C
                   KATT(NATT,1)=K(I,2)
                   KATT(NATT,2)=NFTP
                   KATT(NATT,4)=K(I,1)
-                  IF(K(I,3).EQ.0 .OR. K(K(I,3),2).EQ.IDSTR) THEN
+                   IF(K(I,3).EQ.0 .OR. K(K(I,3),2).EQ.IDSTR) THEN
                      KATT(NATT,3)=0
                   ELSE
                      KATT(NATT,3)=NATT-I+K(I,3)+N_STR-K(K(I,3),4)
@@ -641,12 +806,11 @@ C       ****** identify the mother particle
                   PATT(NATT,2)=P(I,2)
                   PATT(NATT,3)=P(I,3)
                   PATT(NATT,4)=P(I,4)
-                   VATT(NATT,1)=V(I,1)
-                   VATT(NATT,2)=V(I,2)
-                   VATT(NATT,3)=V(I,3)
-                   VATT(NATT,4)=V(I,4)
-
                   EATT=EATT+P(I,4)
+                   VATT(NATT,1)=V(I,1)  
+                   VATT(NATT,2)=V(I,2)  
+                   VATT(NATT,3)=V(I,3)  
+                   VATT(NATT,4)=V(I,4)
 390            CONTINUE 
 400       CONTINUE
 C              ********Fragment the q-qq related string systems
@@ -661,11 +825,10 @@ C         ********Fragment the q-qq related string systems
                PATT(NATT,2)=PDR(I,2)
                PATT(NATT,3)=PDR(I,3)
                PATT(NATT,4)=PDR(I,4)
-                VATT(NATT,1)=V(I,1)
-                VATT(NATT,2)=V(I,2)
-                VATT(NATT,3)=V(I,3)
+                VATT(NATT,1)=V(I,1)  
+                VATT(NATT,2)=V(I,2)  
+                VATT(NATT,3)=V(I,3)  
                 VATT(NATT,4)=V(I,4)
-
                EATT=EATT+PDR(I,4)
 450    CONTINUE
 C                      ********store the direct-produced particles
@@ -678,5 +841,6 @@ C
 C              call LULIST_HIJING(1)
                GO TO 50
        ENDIF
+
        RETURN
        END