HIJING/hipyset1_35/luthru_hijing.F

   1 * $Id$
   2
   3 C*********************************************************************
   4
   5       SUBROUTINE LUTHRU_HIJING(THR,OBL)
   6
   7 C...Purpose: to perform thrust analysis to give thrust, oblateness
   8 C...and the related event axes.
   9 #include "lujets_hijing.inc"
  10 #include "ludat1_hijing.inc"
  11 #include "ludat2_hijing.inc"
  12       DIMENSION TDI(3),TPR(3)
  13
  14 C...Take copy of particles that are to be considered in thrust analysis.
  15       NP=0
  16       PS=0.
  17       DO 100 I=1,N
  18       IF(K(I,1).LE.0.OR.K(I,1).GT.10) GOTO 100
  19       IF(MSTU(41).GE.2) THEN
  20         KC=LUCOMP_HIJING(K(I,2))
  21         IF(KC.EQ.0.OR.KC.EQ.12.OR.KC.EQ.14.OR.KC.EQ.16.OR.
  22      &  KC.EQ.18) GOTO 100
  23         IF(MSTU(41).GE.3.AND.KCHG(KC,2).EQ.0.AND.LUCHGE_HIJING(K(I,2))
  24      $       .EQ.0)GOTO 100
  25       ENDIF
  26       IF(N+NP+MSTU(44)+15.GE.MSTU(4)-MSTU(32)-5) THEN
  27          CALL LUERRM_HIJING(11
  28      $        ,'(LUTHRU_HIJING:) no more memory left in LUJETS_HIJING')
  29         THR=-2.
  30         OBL=-2.
  31         RETURN
  32       ENDIF
  33       NP=NP+1
  34       K(N+NP,1)=23
  35       P(N+NP,1)=P(I,1)
  36       P(N+NP,2)=P(I,2)
  37       P(N+NP,3)=P(I,3)
  38       P(N+NP,4)=SQRT(P(I,1)**2+P(I,2)**2+P(I,3)**2)
  39       P(N+NP,5)=1.
  40       IF(ABS(PARU(42)-1.).GT.0.001) P(N+NP,5)=P(N+NP,4)**(PARU(42)-1.)
  41       PS=PS+P(N+NP,4)*P(N+NP,5)
  42   100 CONTINUE
  43
  44 C...Very low multiplicities (0 or 1) not considered.
  45       IF(NP.LE.1) THEN
  46          CALL LUERRM_HIJING(8
  47      $        ,'(LUTHRU_HIJING:) too few particles for analysis')
  48         THR=-1.
  49         OBL=-1.
  50         RETURN
  51       ENDIF
  52
  53 C...Loop over thrust and major. T axis along z direction in latter case.
  54       DO 280 ILD=1,2
  55       IF(ILD.EQ.2) THEN
  56         K(N+NP+1,1)=31
  57         PHI=ULANGL_HIJING(P(N+NP+1,1),P(N+NP+1,2))
  58         CALL LUDBRB_HIJING(N+1,N+NP+1,0.,-PHI,0D0,0D0,0D0)
  59         THE=ULANGL_HIJING(P(N+NP+1,3),P(N+NP+1,1))
  60         CALL LUDBRB_HIJING(N+1,N+NP+1,-THE,0.,0D0,0D0,0D0)
  61       ENDIF
  62
  63 C...Find and order particles with highest p (pT for major).
  64       DO 110 ILF=N+NP+4,N+NP+MSTU(44)+4
  65   110 P(ILF,4)=0.
  66       DO 150 I=N+1,N+NP
  67       IF(ILD.EQ.2) P(I,4)=SQRT(P(I,1)**2+P(I,2)**2)
  68       DO 120 ILF=N+NP+MSTU(44)+3,N+NP+4,-1
  69       IF(P(I,4).LE.P(ILF,4)) GOTO 130
  70       DO 120 J=1,5
  71   120 P(ILF+1,J)=P(ILF,J)
  72       ILF=N+NP+3
  73   130 DO 140 J=1,5
  74   140 P(ILF+1,J)=P(I,J)
  75   150 CONTINUE
  76
  77 C...Find and order initial axes with highest thrust (major).
  78       DO 160 ILG=N+NP+MSTU(44)+5,N+NP+MSTU(44)+15
  79   160 P(ILG,4)=0.
  80       NC=2**(MIN(MSTU(44),NP)-1)
  81       DO 220 ILC=1,NC
  82       DO 170 J=1,3
  83   170 TDI(J)=0.
  84       DO 180 ILF=1,MIN(MSTU(44),NP)
  85       SGN=P(N+NP+ILF+3,5)
  86       IF(2**ILF*((ILC+2**(ILF-1)-1)/2**ILF).GE.ILC) SGN=-SGN
  87       DO 180 J=1,4-ILD
  88   180 TDI(J)=TDI(J)+SGN*P(N+NP+ILF+3,J)
  89       TDS=TDI(1)**2+TDI(2)**2+TDI(3)**2
  90       DO 190 ILG=N+NP+MSTU(44)+MIN(ILC,10)+4,N+NP+MSTU(44)+5,-1
  91       IF(TDS.LE.P(ILG,4)) GOTO 200
  92       DO 190 J=1,4
  93   190 P(ILG+1,J)=P(ILG,J)
  94       ILG=N+NP+MSTU(44)+4
  95   200 DO 210 J=1,3
  96   210 P(ILG+1,J)=TDI(J)
  97       P(ILG+1,4)=TDS
  98   220 CONTINUE
  99
 100 C...Iterate direction of axis until stable maximum.
 101       P(N+NP+ILD,4)=0.
 102       ILG=0
 103   230 ILG=ILG+1
 104       THP=0.
 105   240 THPS=THP
 106       DO 250 J=1,3
 107       IF(THP.LE.1E-10) TDI(J)=P(N+NP+MSTU(44)+4+ILG,J)
 108       IF(THP.GT.1E-10) TDI(J)=TPR(J)
 109   250 TPR(J)=0.
 110       DO 260 I=N+1,N+NP
 111       SGN=SIGN(P(I,5),TDI(1)*P(I,1)+TDI(2)*P(I,2)+TDI(3)*P(I,3))
 112       DO 260 J=1,4-ILD
 113   260 TPR(J)=TPR(J)+SGN*P(I,J)
 114       THP=SQRT(TPR(1)**2+TPR(2)**2+TPR(3)**2)/PS
 115       IF(THP.GE.THPS+PARU(48)) GOTO 240
 116
 117 C...Save good axis. Try new initial axis until a number of tries agree.
 118       IF(THP.LT.P(N+NP+ILD,4)-PARU(48).AND.ILG.LT.MIN(10,NC)) GOTO 230
 119       IF(THP.GT.P(N+NP+ILD,4)+PARU(48)) THEN
 120         IAGR=0
 121         SGN=(-1.)**INT(RLU_HIJING(0)+0.5)
 122         DO 270 J=1,3
 123   270   P(N+NP+ILD,J)=SGN*TPR(J)/(PS*THP)
 124         P(N+NP+ILD,4)=THP
 125         P(N+NP+ILD,5)=0.
 126       ENDIF
 127       IAGR=IAGR+1
 128   280 IF(IAGR.LT.MSTU(45).AND.ILG.LT.MIN(10,NC)) GOTO 230
 129
 130 C...Find minor axis and value by orthogonality.
 131       SGN=(-1.)**INT(RLU_HIJING(0)+0.5)
 132       P(N+NP+3,1)=-SGN*P(N+NP+2,2)
 133       P(N+NP+3,2)=SGN*P(N+NP+2,1)
 134       P(N+NP+3,3)=0.
 135       THP=0.
 136       DO 290 I=N+1,N+NP
 137   290 THP=THP+P(I,5)*ABS(P(N+NP+3,1)*P(I,1)+P(N+NP+3,2)*P(I,2))
 138       P(N+NP+3,4)=THP/PS
 139       P(N+NP+3,5)=0.
 140
 141 C...Fill axis information. Rotate back to original coordinate system.
 142       DO 300 ILD=1,3
 143       K(N+ILD,1)=31
 144       K(N+ILD,2)=96
 145       K(N+ILD,3)=ILD
 146       K(N+ILD,4)=0
 147       K(N+ILD,5)=0
 148       DO 300 J=1,5
 149       P(N+ILD,J)=P(N+NP+ILD,J)
 150   300 V(N+ILD,J)=0.
 151       CALL LUDBRB_HIJING(N+1,N+3,THE,PHI,0D0,0D0,0D0)
 152
 153 C...Select storing option. Calculate thurst and oblateness.
 154       MSTU(61)=N+1
 155       MSTU(62)=NP
 156       IF(MSTU(43).LE.1) MSTU(3)=3
 157       IF(MSTU(43).GE.2) N=N+3
 158       THR=P(N+1,4)
 159       OBL=P(N+2,4)-P(N+3,4)
 160
 161       RETURN
 162       END