HIJING/hipyset1_35/lujmas_hijing.F

   1 * $Id$
   2
   3 C*********************************************************************
   4
   5       SUBROUTINE LUJMAS_HIJING(PMH,PML)
   6
   7 C...Purpose: to determine, approximately, the two jet masses that
   8 C...minimize the sum m_H|2 + m_L|2, a la Clavelli and Wyler.
   9 #include "lujets_hijing.inc"
  10 #include "ludat1_hijing.inc"
  11 #include "ludat2_hijing.inc"
  12       DIMENSION SM(3,3),SAX(3),PS(3,5)
  13
  14 C...Reset.
  15       NP=0
  16       DO 110 J1=1,3
  17       DO 100 J2=J1,3
  18   100 SM(J1,J2)=0.
  19       DO 110 J2=1,4
  20   110 PS(J1,J2)=0.
  21       PSS=0.
  22
  23 C...Take copy of particles that are to be considered in mass analysis.
  24       DO 150 I=1,N
  25       IF(K(I,1).LE.0.OR.K(I,1).GT.10) GOTO 150
  26       IF(MSTU(41).GE.2) THEN
  27         KC=LUCOMP_HIJING(K(I,2))
  28         IF(KC.EQ.0.OR.KC.EQ.12.OR.KC.EQ.14.OR.KC.EQ.16.OR.
  29      &  KC.EQ.18) GOTO 150
  30         IF(MSTU(41).GE.3.AND.KCHG(KC,2).EQ.0.AND.LUCHGE_HIJING(K(I,2))
  31      $       .EQ.0)GOTO 150
  32       ENDIF
  33       IF(N+NP+1.GE.MSTU(4)-MSTU(32)-5) THEN
  34          CALL LUERRM_HIJING(11
  35      $        ,'(LUJMAS_HIJING:) no more memory left in LUJETS_HIJING')
  36         PMH=-2.
  37         PML=-2.
  38         RETURN
  39       ENDIF
  40       NP=NP+1
  41       DO 120 J=1,5
  42   120 P(N+NP,J)=P(I,J)
  43       IF(MSTU(42).EQ.0) P(N+NP,5)=0.
  44       IF(MSTU(42).EQ.1.AND.K(I,2).NE.22) P(N+NP,5)=PMAS(101,1)
  45       P(N+NP,4)=SQRT(P(N+NP,5)**2+P(I,1)**2+P(I,2)**2+P(I,3)**2)
  46
  47 C...Fill information in sphericity tensor and total momentum vector.
  48       DO 130 J1=1,3
  49       DO 130 J2=J1,3
  50   130 SM(J1,J2)=SM(J1,J2)+P(I,J1)*P(I,J2)
  51       PSS=PSS+(P(I,1)**2+P(I,2)**2+P(I,3)**2)
  52       DO 140 J=1,4
  53   140 PS(3,J)=PS(3,J)+P(N+NP,J)
  54   150 CONTINUE
  55
  56 C...Very low multiplicities (0 or 1) not considered.
  57       IF(NP.LE.1) THEN
  58          CALL LUERRM_HIJING(8
  59      $        ,'(LUJMAS_HIJING:) too few particles for analysis')
  60         PMH=-1.
  61         PML=-1.
  62         RETURN
  63       ENDIF
  64       PARU(61)=SQRT(MAX(0.,PS(3,4)**2-PS(3,1)**2-PS(3,2)**2-PS(3,3)**2))
  65
  66 C...Find largest eigenvalue to matrix (third degree equation).
  67       DO 160 J1=1,3
  68       DO 160 J2=J1,3
  69   160 SM(J1,J2)=SM(J1,J2)/PSS
  70       SQ=(SM(1,1)*SM(2,2)+SM(1,1)*SM(3,3)+SM(2,2)*SM(3,3)-SM(1,2)**2-
  71      &SM(1,3)**2-SM(2,3)**2)/3.-1./9.
  72       SR=-0.5*(SQ+1./9.+SM(1,1)*SM(2,3)**2+SM(2,2)*SM(1,3)**2+SM(3,3)*
  73      &SM(1,2)**2-SM(1,1)*SM(2,2)*SM(3,3))+SM(1,2)*SM(1,3)*SM(2,3)+1./27.
  74       SP=COS(ACOS(MAX(MIN(SR/SQRT(-SQ**3),1.),-1.))/3.)
  75       SMA=1./3.+SQRT(-SQ)*MAX(2.*SP,SQRT(3.*(1.-SP**2))-SP)
  76
  77 C...Find largest eigenvector by solving equation system.
  78       DO 170 J1=1,3
  79       SM(J1,J1)=SM(J1,J1)-SMA
  80       DO 170 J2=J1+1,3
  81   170 SM(J2,J1)=SM(J1,J2)
  82       SMAX=0.
  83       DO 180 J1=1,3
  84       DO 180 J2=1,3
  85       IF(ABS(SM(J1,J2)).LE.SMAX) GOTO 180
  86       JA=J1
  87       JB=J2
  88       SMAX=ABS(SM(J1,J2))
  89   180 CONTINUE
  90       SMAX=0.
  91       DO 190 J3=JA+1,JA+2
  92       J1=J3-3*((J3-1)/3)
  93       RL=SM(J1,JB)/SM(JA,JB)
  94       DO 190 J2=1,3
  95       SM(J1,J2)=SM(J1,J2)-RL*SM(JA,J2)
  96       IF(ABS(SM(J1,J2)).LE.SMAX) GOTO 190
  97       JC=J1
  98       SMAX=ABS(SM(J1,J2))
  99   190 CONTINUE
 100       JB1=JB+1-3*(JB/3)
 101       JB2=JB+2-3*((JB+1)/3)
 102       SAX(JB1)=-SM(JC,JB2)
 103       SAX(JB2)=SM(JC,JB1)
 104       SAX(JB)=-(SM(JA,JB1)*SAX(JB1)+SM(JA,JB2)*SAX(JB2))/SM(JA,JB)
 105
 106 C...Divide particles into two initial clusters by hemisphere.
 107       DO 200 I=N+1,N+NP
 108       PSAX=P(I,1)*SAX(1)+P(I,2)*SAX(2)+P(I,3)*SAX(3)
 109       IS=1
 110       IF(PSAX.LT.0.) IS=2
 111       K(I,3)=IS
 112       DO 200 J=1,4
 113   200 PS(IS,J)=PS(IS,J)+P(I,J)
 114       PMS=(PS(1,4)**2-PS(1,1)**2-PS(1,2)**2-PS(1,3)**2)+
 115      &(PS(2,4)**2-PS(2,1)**2-PS(2,2)**2-PS(2,3)**2)
 116
 117 C...Reassign one particle at a time; find maximum decrease of m|2 sum.
 118   210 PMD=0.
 119       IM=0
 120       DO 220 J=1,4
 121   220 PS(3,J)=PS(1,J)-PS(2,J)
 122       DO 230 I=N+1,N+NP
 123       PPS=P(I,4)*PS(3,4)-P(I,1)*PS(3,1)-P(I,2)*PS(3,2)-P(I,3)*PS(3,3)
 124       IF(K(I,3).EQ.1) PMDI=2.*(P(I,5)**2-PPS)
 125       IF(K(I,3).EQ.2) PMDI=2.*(P(I,5)**2+PPS)
 126       IF(PMDI.LT.PMD) THEN
 127         PMD=PMDI
 128         IM=I
 129       ENDIF
 130   230 CONTINUE
 131
 132 C...Loop back if significant reduction in sum of m|2.
 133       IF(PMD.LT.-PARU(48)*PMS) THEN
 134         PMS=PMS+PMD
 135         IS=K(IM,3)
 136         DO 240 J=1,4
 137         PS(IS,J)=PS(IS,J)-P(IM,J)
 138   240   PS(3-IS,J)=PS(3-IS,J)+P(IM,J)
 139         K(IM,3)=3-IS
 140         GOTO 210
 141       ENDIF
 142
 143 C...Final masses and output.
 144       MSTU(61)=N+1
 145       MSTU(62)=NP
 146       PS(1,5)=SQRT(MAX(0.,PS(1,4)**2-PS(1,1)**2-PS(1,2)**2-PS(1,3)**2))
 147       PS(2,5)=SQRT(MAX(0.,PS(2,4)**2-PS(2,1)**2-PS(2,2)**2-PS(2,3)**2))
 148       PMH=MAX(PS(1,5),PS(2,5))
 149       PML=MIN(PS(1,5),PS(2,5))
 150
 151       RETURN
 152       END