PYTHIA/jetset/luthru.F

   1
   2 C*********************************************************************
   3
   4       SUBROUTINE LUTHRU(THR,OBL)
   5
   6 C...Purpose: to perform thrust analysis to give thrust, oblateness
   7 C...and the related event axes.
   8       COMMON/LUJETS/N,K(4000,5),P(4000,5),V(4000,5)
   9       COMMON/LUDAT1/MSTU(200),PARU(200),MSTJ(200),PARJ(200)
  10       COMMON/LUDAT2/KCHG(500,3),PMAS(500,4),PARF(2000),VCKM(4,4)
  11       SAVE /LUJETS/,/LUDAT1/,/LUDAT2/
  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(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(K(I,2)).EQ.0)
  24      &  GOTO 100
  25       ENDIF
  26       IF(N+NP+MSTU(44)+15.GE.MSTU(4)-MSTU(32)-5) THEN
  27         CALL LUERRM(11,'(LUTHRU:) no more memory left in LUJETS')
  28         THR=-2.
  29         OBL=-2.
  30         RETURN
  31       ENDIF
  32       NP=NP+1
  33       K(N+NP,1)=23
  34       P(N+NP,1)=P(I,1)
  35       P(N+NP,2)=P(I,2)
  36       P(N+NP,3)=P(I,3)
  37       P(N+NP,4)=SQRT(P(I,1)**2+P(I,2)**2+P(I,3)**2)
  38       P(N+NP,5)=1.
  39       IF(ABS(PARU(42)-1.).GT.0.001) P(N+NP,5)=P(N+NP,4)**(PARU(42)-1.)
  40       PS=PS+P(N+NP,4)*P(N+NP,5)
  41   100 CONTINUE
  42
  43 C...Very low multiplicities (0 or 1) not considered.
  44       IF(NP.LE.1) THEN
  45         CALL LUERRM(8,'(LUTHRU:) too few particles for analysis')
  46         THR=-1.
  47         OBL=-1.
  48         RETURN
  49       ENDIF
  50
  51 C...Loop over thrust and major. T axis along z direction in latter case.
  52       DO 320 ILD=1,2
  53       IF(ILD.EQ.2) THEN
  54         K(N+NP+1,1)=31
  55         PHI=ULANGL(P(N+NP+1,1),P(N+NP+1,2))
  56         MSTU(33)=1
  57         CALL LUDBRB(N+1,N+NP+1,0.,-PHI,0D0,0D0,0D0)
  58         THE=ULANGL(P(N+NP+1,3),P(N+NP+1,1))
  59         CALL LUDBRB(N+1,N+NP+1,-THE,0.,0D0,0D0,0D0)
  60       ENDIF
  61
  62 C...Find and order particles with highest p (pT for major).
  63       DO 110 ILF=N+NP+4,N+NP+MSTU(44)+4
  64       P(ILF,4)=0.
  65   110 CONTINUE
  66       DO 160 I=N+1,N+NP
  67       IF(ILD.EQ.2) P(I,4)=SQRT(P(I,1)**2+P(I,2)**2)
  68       DO 130 ILF=N+NP+MSTU(44)+3,N+NP+4,-1
  69       IF(P(I,4).LE.P(ILF,4)) GOTO 140
  70       DO 120 J=1,5
  71       P(ILF+1,J)=P(ILF,J)
  72   120 CONTINUE
  73   130 CONTINUE
  74       ILF=N+NP+3
  75   140 DO 150 J=1,5
  76       P(ILF+1,J)=P(I,J)
  77   150 CONTINUE
  78   160 CONTINUE
  79
  80 C...Find and order initial axes with highest thrust (major).
  81       DO 170 ILG=N+NP+MSTU(44)+5,N+NP+MSTU(44)+15
  82       P(ILG,4)=0.
  83   170 CONTINUE
  84       NC=2**(MIN(MSTU(44),NP)-1)
  85       DO 250 ILC=1,NC
  86       DO 180 J=1,3
  87       TDI(J)=0.
  88   180 CONTINUE
  89       DO 200 ILF=1,MIN(MSTU(44),NP)
  90       SGN=P(N+NP+ILF+3,5)
  91       IF(2**ILF*((ILC+2**(ILF-1)-1)/2**ILF).GE.ILC) SGN=-SGN
  92       DO 190 J=1,4-ILD
  93       TDI(J)=TDI(J)+SGN*P(N+NP+ILF+3,J)
  94   190 CONTINUE
  95   200 CONTINUE
  96       TDS=TDI(1)**2+TDI(2)**2+TDI(3)**2
  97       DO 220 ILG=N+NP+MSTU(44)+MIN(ILC,10)+4,N+NP+MSTU(44)+5,-1
  98       IF(TDS.LE.P(ILG,4)) GOTO 230
  99       DO 210 J=1,4
 100       P(ILG+1,J)=P(ILG,J)
 101   210 CONTINUE
 102   220 CONTINUE
 103       ILG=N+NP+MSTU(44)+4
 104   230 DO 240 J=1,3
 105       P(ILG+1,J)=TDI(J)
 106   240 CONTINUE
 107       P(ILG+1,4)=TDS
 108   250 CONTINUE
 109
 110 C...Iterate direction of axis until stable maximum.
 111       P(N+NP+ILD,4)=0.
 112       ILG=0
 113   260 ILG=ILG+1
 114       THP=0.
 115   270 THPS=THP
 116       DO 280 J=1,3
 117       IF(THP.LE.1E-10) TDI(J)=P(N+NP+MSTU(44)+4+ILG,J)
 118       IF(THP.GT.1E-10) TDI(J)=TPR(J)
 119       TPR(J)=0.
 120   280 CONTINUE
 121       DO 300 I=N+1,N+NP
 122       SGN=SIGN(P(I,5),TDI(1)*P(I,1)+TDI(2)*P(I,2)+TDI(3)*P(I,3))
 123       DO 290 J=1,4-ILD
 124       TPR(J)=TPR(J)+SGN*P(I,J)
 125   290 CONTINUE
 126   300 CONTINUE
 127       THP=SQRT(TPR(1)**2+TPR(2)**2+TPR(3)**2)/PS
 128       IF(THP.GE.THPS+PARU(48)) GOTO 270
 129
 130 C...Save good axis. Try new initial axis until a number of tries agree.
 131       IF(THP.LT.P(N+NP+ILD,4)-PARU(48).AND.ILG.LT.MIN(10,NC)) GOTO 260
 132       IF(THP.GT.P(N+NP+ILD,4)+PARU(48)) THEN
 133         IAGR=0
 134         SGN=(-1.)**INT(RLU(0)+0.5)
 135         DO 310 J=1,3
 136         P(N+NP+ILD,J)=SGN*TPR(J)/(PS*THP)
 137   310   CONTINUE
 138         P(N+NP+ILD,4)=THP
 139         P(N+NP+ILD,5)=0.
 140       ENDIF
 141       IAGR=IAGR+1
 142       IF(IAGR.LT.MSTU(45).AND.ILG.LT.MIN(10,NC)) GOTO 260
 143   320 CONTINUE
 144
 145 C...Find minor axis and value by orthogonality.
 146       SGN=(-1.)**INT(RLU(0)+0.5)
 147       P(N+NP+3,1)=-SGN*P(N+NP+2,2)
 148       P(N+NP+3,2)=SGN*P(N+NP+2,1)
 149       P(N+NP+3,3)=0.
 150       THP=0.
 151       DO 330 I=N+1,N+NP
 152       THP=THP+P(I,5)*ABS(P(N+NP+3,1)*P(I,1)+P(N+NP+3,2)*P(I,2))
 153   330 CONTINUE
 154       P(N+NP+3,4)=THP/PS
 155       P(N+NP+3,5)=0.
 156
 157 C...Fill axis information. Rotate back to original coordinate system.
 158       DO 350 ILD=1,3
 159       K(N+ILD,1)=31
 160       K(N+ILD,2)=96
 161       K(N+ILD,3)=ILD
 162       K(N+ILD,4)=0
 163       K(N+ILD,5)=0
 164       DO 340 J=1,5
 165       P(N+ILD,J)=P(N+NP+ILD,J)
 166       V(N+ILD,J)=0.
 167   340 CONTINUE
 168   350 CONTINUE
 169       CALL LUDBRB(N+1,N+3,THE,PHI,0D0,0D0,0D0)
 170
 171 C...Calculate thrust and oblateness. Select storing option.
 172       THR=P(N+1,4)
 173       OBL=P(N+2,4)-P(N+3,4)
 174       MSTU(61)=N+1
 175       MSTU(62)=NP
 176       IF(MSTU(43).LE.1) MSTU(3)=3
 177       IF(MSTU(43).GE.2) N=N+3
 178
 179       RETURN
 180       END