PHOS/shaker/luthru.f

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