PartonShowerVeto/src/ktclusdble.f

0001 C-----------------------------------------------------------------------
0002 C-----------------------------------------------------------------------
0003 C-----------------------------------------------------------------------
0004 C     KTCLUS: written by Mike Seymour, July 1992.
0005 C     Last modified November 2000.
0006 C     Please send comments or suggestions to Mike.Seymour@rl.ac.uk
0007 C
0008 C     This is a general-purpose kt clustering package.
0009 C     It can handle ee, ep and pp collisions.
0010 C     It is loosely based on the program of Siggi Bethke.
0011 C
0012 C     The time taken (on a 10MIP machine) is (0.2microsec)*N**3
0013 C     where N is the number of particles.
0014 C     Over 90 percent of this time is used in subroutine KTPMIN, which
0015 C     simply finds the minimum member of a one-dimensional array.
0016 C     It is well worth thinking about optimization: on the SPARCstation
0017 C     a factor of two increase was obtained simply by increasing the
0018 C     optimization level from its default value.
0019 C
0020 C     The approach is to separate the different stages of analysis.
0021 C     KTCLUS does all the clustering and records a merging history.
0022 C     It returns a simple list of the y values at which each merging
0023 C     occured. Then the following routines can be called to give extra
0024 C     information on the most recently analysed event.
0025 C     KTCLUR is identical but includes an R parameter, see below.
0026 C     KTYCUT gives the number of jets at each given YCUT value.
0027 C     KTYSUB gives the number of sub-jets at each given YCUT value.
0028 C     KTBEAM gives same info as KTCLUS but only for merges with the beam
0029 C     KTJOIN gives same info as KTCLUS but for merges of sub-jets.
0030 C     KTRECO reconstructs the jet momenta at a given value of YCUT.
0031 C     It also gives information on which jets at scale YCUT belong to
0032 C     which macro-jets at scale YMAC, for studying sub-jet properties.
0033 C     KTINCL reconstructs the jet momenta according to the inclusive jet
0034 C     definition of Ellis and Soper.
0035 C     KTISUB, KTIJOI and KTIREC are like KTYSUB, KTJOIN and KTRECO,
0036 C     except that they only apply to one inclusive jet at a time,
0037 C     with the pt of that jet automatically used for ECUT.
0038 C     KTWICH gives a list of which particles ended up in which jets.
0039 C     KTWCHS gives the same thing, but only for subjets.
0040 C     Note that the numbering of jets used by these two routines is
0041 C     guaranteed to be the same as that used by KTRECO.
0042 C
0043 C     The collision type and analysis type are indicated by the first
0044 C     argument of KTCLUS. IMODE=<TYPE><ANGLE><MONO><RECOM> where
0045 C     TYPE:  1=>ee, 2=>ep with p in -z direction, 3=>pe, 4=>pp
0046 C     ANGLE: 1=>angular kt def., 2=>DeltaR, 3=>f(DeltaEta,DeltaPhi)
0047 C            where f()=2(cosh(eta)-cos(phi)) is the QCD emission metric
0048 C     MONO:  1=>derive relative pseudoparticle angles from jets
0049 C            2=>monotonic definitions of relative angles
0050 C     RECOM: 1=>E recombination scheme, 2=>pt scheme, 3=>pt**2 scheme
0051 C
0052 C     There are also abbreviated forms for the most common combinations:
0053 C     IMODE=1 => E scheme in e+e-                              (=1111)
0054 C           2 => E scheme in ep                                (=2111)
0055 C           3 => E scheme in pe                                (=3111)
0056 C           4 => E scheme in pp                                (=4111)
0057 C           5 => covariant E scheme in pp                      (=4211)
0058 C           6 => covariant pt-scheme in pp                     (=4212)
0059 C           7 => covariant monotonic pt**2-scheme in pp        (=4223)
0060 C
0061 C     KTRECO no longer needs to reconstruct the momenta according to the
0062 C     same recombination scheme in which they were clustered. Its first
0063 C     argument gives the scheme, taking the same values as RECOM above.
0064 C
0065 C     Note that unlike previous versions, all variables which hold y
0066 C     values have been named in a consistent way:
0067 C     Y()  is the output scale at which jets were merged,
0068 C     YCUT is the input scale at which jets should be counted, and
0069 C          jet-momenta reconstructed etc,
0070 C     YMAC is the input macro-jet scale, used in determining whether
0071 C          or not each jet is a sub-jet.
0072 C     The original scheme defined in our papers is equivalent to always
0073 C     setting YMAC=1.
0074 C     Whenever a YCUT or YMAC variable is used, it is rounded down
0075 C     infinitesimally, so that for example, setting YCUT=Y(2) refers
0076 C     to the scale where the event is 2-jet, even if rounding errors
0077 C     have shifted its value slightly.
0078 C
0079 C     An R parameter can be used in hadron-hadron collisions by
0080 C     calling KTCLUR instead of KTCLUS.  This is as suggested by
0081 C     Ellis and Soper, but implemented slightly differently,
0082 C     as in M.H. Seymour, LU TP 94/2 (submitted to Nucl. Phys. B.).
0083 C     R**2 multiplies the single Kt everywhere it is used.
0084 C     Calling KTCLUR with R=1 is identical to calling KTCLUS.
0085 C     R plays a similar role to the jet radius in a cone-type algorithm,
0086 C     but is scaled up by about 40% (ie R=0.7 in a cone algorithm is
0087 C     similar to this algorithm with R=1).
0088 C     Note that R.EQ.1 must be used for the e+e- and ep versions,
0089 C     and is strongly recommended for the hadron-hadron version.
0090 C     However, R values smaller than 1 have been found to be useful for
0091 C     certain applications, particularly the mass reconstruction of
0092 C     highly-boosted colour-singlets such as high-pt hadronic Ws,
0093 C     as in M.H. Seymour, LU TP 93/8 (to appear in Z. Phys. C.).
0094 C     Situations in which R<1 is useful are likely to also be those in
0095 C     which the inclusive reconstruction method is more useful.
0096 C
0097 C     Also included is a set of routines for doing Lorentz boosts:
0098 C     KTLBST finds the boost matrix to/from the cm frame of a 4-vector
0099 C     KTRROT finds the rotation matrix from one vector to another
0100 C     KTMMUL multiplies together two matrices
0101 C     KTVMUL multiplies a vector by a matrix
0102 C     KTINVT inverts a transformation matrix (nb NOT a general 4 by 4)
0103 C     KTFRAM boosts a list of vectors between two arbitrary frames
0104 C     KTBREI boosts a list of vectors between the lab and Breit frames
0105 C     KTHADR boosts a list of vectors between the lab and hadronic cmf
0106 C       The last two need the momenta in the +z direction of the lepton
0107 C       and hadron beams, and the 4-momentum of the outgoing lepton.
0108 C
0109 C     The main reference is:
0110 C       S. Catani, Yu.L. Dokshitzer, M.H. Seymour and B.R. Webber,
0111 C         Nucl.Phys.B406(1993)187.
0112 C     The ep version was proposed in:
0113 C       S. Catani, Yu.L. Dokshitzer and B.R. Webber,
0114 C         Phys.Lett.285B(1992)291.
0115 C     The inclusive reconstruction method was proposed in:
0116 C       S.D. Ellis and D.E. Soper,
0117 C         Phys.Rev.D48(1993)3160.
0118 C
0119 C-----------------------------------------------------------------------
0120 C-----------------------------------------------------------------------
0121 C-----------------------------------------------------------------------
0122       SUBROUTINE KTCLUS(IMODE,PP,NN,ECUT,Y,*)
0123       IMPLICIT NONE
0124 C---DO CLUSTER ANALYSIS OF PARTICLES IN PP
0125 C
0126 C   IMODE   = INPUT  : DESCRIBED ABOVE
0127 C   PP(I,J) = INPUT  : 4-MOMENTUM OF Jth PARTICLE: I=1,4 => PX,PY,PZ,E
0128 C   NN      = INPUT  : NUMBER OF PARTICLES
0129 C   ECUT    = INPUT  : DENOMINATOR OF KT MEASURE. IF ZERO, ETOT IS USED
0130 C   Y(J)    = OUTPUT : VALUE OF Y FOR WHICH EVENT CHANGES FROM BEING
0131 C                        J JET TO J-1 JET
0132 C   LAST ARGUMENT IS LABEL TO JUMP TO IF FOR ANY REASON THE EVENT
0133 C   COULD NOT BE PROCESSED (MOST LIKELY DUE TO TOO MANY PARTICLES)
0134 C
0135 C   NOTE THAT THE MOMENTA ARE DECLARED DOUBLE PRECISION,
0136 C   AND ALL OTHER FLOATING POINT VARIABLES ARE DECLARED DOUBLE PRECISION
0137 C
0138       INTEGER IMODE,NN
0139       DOUBLE PRECISION PP(4,*)
0140       DOUBLE PRECISION ECUT,Y(*),ONE
0141       ONE=1
0142       CALL KTCLUR(IMODE,PP,NN,ONE,ECUT,Y,*999)
0143       RETURN
0144  999  RETURN 1
0145       END
0146 C-----------------------------------------------------------------------
0147       SUBROUTINE KTCLUR(IMODE,PP,NN,R,ECUT,Y,*)
0148       IMPLICIT NONE
0149 C---DO CLUSTER ANALYSIS OF PARTICLES IN PP
0150 C
0151 C   IMODE   = INPUT  : DESCRIBED ABOVE
0152 C   PP(I,J) = INPUT  : 4-MOMENTUM OF Jth PARTICLE: I=1,4 => PX,PY,PZ,E
0153 C   NN      = INPUT  : NUMBER OF PARTICLES
0154 C   R       = INPUT  : ELLIS AND SOPER'S R PARAMETER, SEE ABOVE.
0155 C   ECUT    = INPUT  : DENOMINATOR OF KT MEASURE. IF ZERO, ETOT IS USED
0156 C   Y(J)    = OUTPUT : VALUE OF Y FOR WHICH EVENT CHANGES FROM BEING
0157 C                        J JET TO J-1 JET
0158 C   LAST ARGUMENT IS LABEL TO JUMP TO IF FOR ANY REASON THE EVENT
0159 C   COULD NOT BE PROCESSED (MOST LIKELY DUE TO TOO MANY PARTICLES)
0160 C
0161 C   NOTE THAT THE MOMENTA ARE DECLARED DOUBLE PRECISION,
0162 C   AND ALL OTHER FLOATING POINT VARIABLES ARE DECLARED DOUBLE PRECISION
0163 C
0164       INTEGER NMAX,IM,IMODE,TYPE,ANGL,MONO,RECO,N,I,J,NN,
0165      &     IMIN,JMIN,KMIN,NUM,HIST,INJET,IABBR,NABBR
0166       PARAMETER (NMAX=512,NABBR=7)
0167       DOUBLE PRECISION PP(4,*)
0168       DOUBLE PRECISION R,ECUT,Y(*),P,KT,ETOT,RSQ,KTP,KTS,KTPAIR,KTSING,
0169      &     KTMIN,ETSQ,KTLAST,KTMAX,KTTMP
0170       LOGICAL FIRST
0171       CHARACTER TITLE(4,4)*10
0172 C---KT RECORDS THE KT**2 OF EACH MERGING.
0173 C---KTLAST RECORDS FOR EACH MERGING, THE HIGHEST ECUT**2 FOR WHICH THE
0174 C   RESULT IS NOT MERGED WITH THE BEAM (COULD BE LARGER THAN THE
0175 C   KT**2 AT WHICH IT WAS MERGED IF THE KT VALUES ARE NOT MONOTONIC).
0176 C   THIS MAY SOUND POINTLESS, BUT ITS USEFUL FOR DETERMINING WHETHER
0177 C   SUB-JETS SURVIVED TO SCALE Y=YMAC OR NOT.
0178 C---HIST RECORDS MERGING HISTORY:
0179 C   N=>DELETED TRACK N, M*NMAX+N=>MERGED TRACKS M AND N (M<N).
0180       COMMON /KTCOMM/ETOT,RSQ,P(9,NMAX),KTP(NMAX,NMAX),KTS(NMAX),
0181      &  KT(NMAX),KTLAST(NMAX),HIST(NMAX),NUM
0182       DIMENSION INJET(NMAX),IABBR(NABBR)
0183       DATA FIRST,TITLE,IABBR/.TRUE.,
0184      &     'e+e-      ','ep        ','pe        ','pp        ',
0185      &     'angle     ','DeltaR    ','f(DeltaR) ','**********',
0186      &     'no        ','yes       ','**********','**********',
0187      &     'E         ','Pt        ','Pt**2     ','**********',
0188      &     1111,2111,3111,4111,4211,4212,4223/
0189 C---CHECK INPUT
0190       IM=IMODE
0191       IF (IM.GE.1.AND.IM.LE.NABBR) IM=IABBR(IM)
0192       TYPE=MOD(IM/1000,10)
0193       ANGL=MOD(IM/100 ,10)
0194       MONO=MOD(IM/10  ,10)
0195       RECO=MOD(IM     ,10)
0196       IF (NN.GT.NMAX.OR.NN.LT.1.OR.(NN.LT.2.AND.TYPE.EQ.1))
0197      &     CALL KTWARN('KTCLUS',100,*999)
0198       IF (TYPE.LT.1.OR.TYPE.GT.4.OR.ANGL.LT.1.OR.ANGL.GT.4.OR.
0199      &    MONO.LT.1.OR.MONO.GT.2.OR.RECO.LT.1.OR.RECO.GT.3)
0200      &     CALL KTWARN('KTCLUS',101,*999)
0201       IF (FIRST) THEN
0202          WRITE (6,'(/,1X,54(''*'')/A)')
0203      &   ' KTCLUS: written by Mike Seymour, July 1992.'
0204          WRITE (6,'(A)')
0205      &   ' Last modified November 2000.'
0206          WRITE (6,'(A)')
0207      &   ' Please send comments or suggestions to Mike.Seymour@rl.ac.uk'
0208          WRITE (6,'(/A,I2,2A)')
0209      &   '       Collision type =',TYPE,' = ',TITLE(TYPE,1)
0210          WRITE (6,'(A,I2,2A)')
0211      &   '     Angular variable =',ANGL,' = ',TITLE(ANGL,2)
0212          WRITE (6,'(A,I2,2A)')
0213      &   ' Monotonic definition =',MONO,' = ',TITLE(MONO,3)
0214          WRITE (6,'(A,I2,2A)')
0215      &   ' Recombination scheme =',RECO,' = ',TITLE(RECO,4)
0216          IF (R.NE.1) THEN
0217          WRITE (6,'(A,F5.2)')
0218      &   '     Radius parameter =',R
0219          IF (TYPE.NE.4) WRITE (6,'(A)')
0220      &   ' R.NE.1 is strongly discouraged for this collision type!'
0221          ENDIF
0222          WRITE (6,'(1X,54(''*'')/)')
0223          FIRST=.FALSE.
0224       ENDIF
0225 C---COPY PP TO P
0226       N=NN
0227       NUM=NN
0228       CALL KTCOPY(PP,N,P,(RECO.NE.1))
0229       ETOT=0
0230       DO 100 I=1,N
0231          ETOT=ETOT+P(4,I)
0232  100  CONTINUE
0233       IF (ETOT.EQ.0) CALL KTWARN('KTCLUS',102,*999)
0234       IF (ECUT.EQ.0) THEN
0235          ETSQ=1/ETOT**2
0236       ELSE
0237          ETSQ=1/ECUT**2
0238       ENDIF
0239       RSQ=R**2
0240 C---CALCULATE ALL PAIR KT's
0241       DO 210 I=1,N-1
0242          DO 200 J=I+1,N
0243             KTP(J,I)=-1
0244             KTP(I,J)=KTPAIR(ANGL,P(1,I),P(1,J),KTP(J,I))
0245  200     CONTINUE
0246  210  CONTINUE
0247 C---CALCULATE ALL SINGLE KT's
0248       DO 230 I=1,N
0249          KTS(I)=KTSING(ANGL,TYPE,P(1,I))
0250  230  CONTINUE
0251       KTMAX=0
0252 C---MAIN LOOP
0253  300  CONTINUE
0254 C---FIND MINIMUM MEMBER OF KTP
0255       CALL KTPMIN(KTP,NMAX,N,IMIN,JMIN)
0256 C---FIND MINIMUM MEMBER OF KTS
0257       CALL KTSMIN(KTS,NMAX,N,KMIN)
0258 C---STORE Y VALUE OF TRANSITION FROM N TO N-1 JETS
0259       KTMIN=KTP(IMIN,JMIN)
0260       KTTMP=RSQ*KTS(KMIN)
0261       IF ((TYPE.GE.2.AND.TYPE.LE.4).AND.
0262      &     (KTTMP.LE.KTMIN.OR.N.EQ.1))
0263      &     KTMIN=KTTMP
0264       KT(N)=KTMIN
0265       Y(N)=KT(N)*ETSQ
0266 C---IF MONO.GT.1, SEQUENCE IS SUPPOSED TO BE MONOTONIC, IF NOT, WARN
0267       IF (KTMIN.LT.KTMAX.AND.MONO.GT.1) CALL KTWARN('KTCLUS',1,*999)
0268       IF (KTMIN.GE.KTMAX) KTMAX=KTMIN
0269 C---IF LOWEST KT IS TO A BEAM, THROW IT AWAY AND MOVE LAST ENTRY UP
0270       IF (KTMIN.EQ.KTTMP) THEN
0271          CALL KTMOVE(P,KTP,KTS,NMAX,N,KMIN,1)
0272 C---UPDATE HISTORY AND CROSS-REFERENCES
0273          HIST(N)=KMIN
0274          INJET(N)=KMIN
0275          DO 400 I=N,NN
0276             IF (INJET(I).EQ.KMIN) THEN
0277                KTLAST(I)=KTMAX
0278                INJET(I)=0
0279             ELSEIF (INJET(I).EQ.N) THEN
0280                INJET(I)=KMIN
0281             ENDIF
0282  400     CONTINUE
0283 C---OTHERWISE MERGE JETS IMIN AND JMIN AND MOVE LAST ENTRY UP
0284       ELSE
0285          CALL KTMERG(P,KTP,KTS,NMAX,IMIN,JMIN,N,TYPE,ANGL,MONO,RECO)
0286          CALL KTMOVE(P,KTP,KTS,NMAX,N,JMIN,1)
0287 C---UPDATE HISTORY AND CROSS-REFERENCES
0288          HIST(N)=IMIN*NMAX+JMIN
0289          INJET(N)=IMIN
0290          DO 600 I=N,NN
0291             IF (INJET(I).EQ.JMIN) THEN
0292                INJET(I)=IMIN
0293             ELSEIF (INJET(I).EQ.N) THEN
0294                INJET(I)=JMIN
0295             ENDIF
0296  600     CONTINUE
0297       ENDIF
0298 C---THATS ALL THERE IS TO IT
0299       N=N-1
0300       IF (N.GT.1 .OR. N.GT.0.AND.(TYPE.GE.2.AND.TYPE.LE.4)) GOTO 300
0301       IF (N.EQ.1) THEN
0302          KT(N)=1D20
0303          Y(N)=KT(N)*ETSQ
0304       ENDIF
0305       RETURN
0306  999  RETURN 1
0307       END
0308 C-----------------------------------------------------------------------
0309       SUBROUTINE KTYCUT(ECUT,NY,YCUT,NJET,*)
0310       IMPLICIT NONE
0311 C---COUNT THE NUMBER OF JETS AT EACH VALUE OF YCUT, FOR EVENT WHICH HAS
0312 C   ALREADY BEEN ANALYSED BY KTCLUS.
0313 C
0314 C   ECUT    = INPUT : DENOMINATOR OF KT MEASURE. IF ZERO, ETOT IS USED
0315 C   NY      = INPUT : NUMBER OF YCUT VALUES
0316 C   YCUT(J) = INPUT : Y VALUES AT WHICH NUMBERS OF JETS ARE COUNTED
0317 C   NJET(J) =OUTPUT : NUMBER OF JETS AT YCUT(J)
0318 C   LAST ARGUMENT IS LABEL TO JUMP TO IF FOR ANY REASON THE EVENT
0319 C   COULD NOT BE PROCESSED
0320 C
0321 C   NOTE THAT ALL FLOATING POINT VARIABLES ARE DECLARED DOUBLE PRECISION
0322 C
0323       INTEGER NY,NJET(NY),NMAX,HIST,I,J,NUM
0324       PARAMETER (NMAX=512)
0325       DOUBLE PRECISION YCUT(NY),ETOT,RSQ,P,KT,KTP,KTS,ETSQ,ECUT,KTLAST,
0326      &     ROUND
0327       PARAMETER (ROUND=0.99999D0)
0328       COMMON /KTCOMM/ETOT,RSQ,P(9,NMAX),KTP(NMAX,NMAX),KTS(NMAX),
0329      &  KT(NMAX),KTLAST(NMAX),HIST(NMAX),NUM
0330       IF (ETOT.EQ.0) CALL KTWARN('KTYCUT',100,*999)
0331       IF (ECUT.EQ.0) THEN
0332          ETSQ=1/ETOT**2
0333       ELSE
0334          ETSQ=1/ECUT**2
0335       ENDIF
0336       DO 100 I=1,NY
0337          NJET(I)=0
0338  100  CONTINUE
0339       DO 210 I=NUM,1,-1
0340          DO 200 J=1,NY
0341             IF (NJET(J).EQ.0.AND.KT(I)*ETSQ.GE.ROUND*YCUT(J)) NJET(J)=I
0342  200     CONTINUE
0343  210  CONTINUE
0344       RETURN
0345  999  RETURN 1
0346       END
0347 C-----------------------------------------------------------------------
0348       SUBROUTINE KTYSUB(ECUT,NY,YCUT,YMAC,NSUB,*)
0349       IMPLICIT NONE
0350 C---COUNT THE NUMBER OF SUB-JETS AT EACH VALUE OF YCUT, FOR EVENT WHICH
0351 C   HAS ALREADY BEEN ANALYSED BY KTCLUS.
0352 C   REMEMBER THAT A SUB-JET IS DEFINED AS A JET AT Y=YCUT WHICH HAS NOT
0353 C   YET BEEN MERGED WITH THE BEAM AT Y=YMAC.
0354 C
0355 C   ECUT    = INPUT : DENOMINATOR OF KT MEASURE. IF ZERO, ETOT IS USED
0356 C   NY      = INPUT : NUMBER OF YCUT VALUES
0357 C   YCUT(J) = INPUT : Y VALUES AT WHICH NUMBERS OF SUB-JETS ARE COUNTED
0358 C   YMAC    = INPUT : Y VALUE USED TO DEFINE MACRO-JETS, TO DETERMINE
0359 C                       WHICH JETS ARE SUB-JETS
0360 C   NSUB(J) =OUTPUT : NUMBER OF SUB-JETS AT YCUT(J)
0361 C   LAST ARGUMENT IS LABEL TO JUMP TO IF FOR ANY REASON THE EVENT
0362 C   COULD NOT BE PROCESSED
0363 C
0364 C   NOTE THAT ALL FLOATING POINT VARIABLES ARE DECLARED DOUBLE PRECISION
0365 C
0366       INTEGER NY,NSUB(NY),NMAX,HIST,I,J,NUM
0367       PARAMETER (NMAX=512)
0368       DOUBLE PRECISION YCUT(NY),YMAC,ETOT,RSQ,P,KT,KTP,KTS,ETSQ,ECUT,
0369      &     KTLAST,ROUND
0370       PARAMETER (ROUND=0.99999D0)
0371       COMMON /KTCOMM/ETOT,RSQ,P(9,NMAX),KTP(NMAX,NMAX),KTS(NMAX),
0372      &  KT(NMAX),KTLAST(NMAX),HIST(NMAX),NUM
0373       IF (ETOT.EQ.0) CALL KTWARN('KTYSUB',100,*999)
0374       IF (ECUT.EQ.0) THEN
0375          ETSQ=1/ETOT**2
0376       ELSE
0377          ETSQ=1/ECUT**2
0378       ENDIF
0379       DO 100 I=1,NY
0380          NSUB(I)=0
0381  100  CONTINUE
0382       DO 210 I=NUM,1,-1
0383          DO 200 J=1,NY
0384             IF (NSUB(J).EQ.0.AND.KT(I)*ETSQ.GE.ROUND*YCUT(J)) NSUB(J)=I
0385             IF (NSUB(J).NE.0.AND.KTLAST(I)*ETSQ.LT.ROUND*YMAC)
0386      &          NSUB(J)=NSUB(J)-1
0387  200     CONTINUE
0388  210  CONTINUE
0389       RETURN
0390  999  RETURN 1
0391       END
0392 C-----------------------------------------------------------------------
0393       SUBROUTINE KTBEAM(ECUT,Y,*)
0394       IMPLICIT NONE
0395 C---GIVE SAME INFORMATION AS LAST CALL TO KTCLUS EXCEPT THAT ONLY
0396 C   TRANSITIONS WHERE A JET WAS MERGED WITH THE BEAM JET ARE RECORDED
0397 C
0398 C   ECUT    = INPUT : DENOMINATOR OF KT MEASURE. IF ZERO, ETOT IS USED
0399 C   Y(J)    =OUTPUT : Y VALUE WHERE Jth HARDEST JET WAS MERGED WITH BEAM
0400 C   LAST ARGUMENT IS LABEL TO JUMP TO IF FOR ANY REASON THE EVENT
0401 C   COULD NOT BE PROCESSED
0402 C
0403 C   NOTE THAT ALL FLOATING POINT VARIABLES ARE DECLARED DOUBLE PRECISION
0404 C
0405       INTEGER NMAX,HIST,NUM,I,J
0406       PARAMETER (NMAX=512)
0407       DOUBLE PRECISION ETOT,RSQ,P,KT,KTP,KTS,ECUT,ETSQ,Y(*),KTLAST
0408       COMMON /KTCOMM/ETOT,RSQ,P(9,NMAX),KTP(NMAX,NMAX),KTS(NMAX),
0409      &  KT(NMAX),KTLAST(NMAX),HIST(NMAX),NUM
0410       IF (ETOT.EQ.0) CALL KTWARN('KTBEAM',100,*999)
0411       IF (ECUT.EQ.0) THEN
0412          ETSQ=1/ETOT**2
0413       ELSE
0414          ETSQ=1/ECUT**2
0415       ENDIF
0416       J=1
0417       DO 100 I=1,NUM
0418          IF (HIST(I).LE.NMAX) THEN
0419             Y(J)=ETSQ*KT(I)
0420             J=J+1
0421          ENDIF
0422  100  CONTINUE
0423       DO 200 I=J,NUM
0424          Y(I)=0
0425  200  CONTINUE
0426       RETURN
0427  999  RETURN 1
0428       END
0429 C-----------------------------------------------------------------------
0430       SUBROUTINE KTJOIN(ECUT,YMAC,Y,*)
0431       IMPLICIT NONE
0432 C---GIVE SAME INFORMATION AS LAST CALL TO KTCLUS EXCEPT THAT ONLY
0433 C   TRANSITIONS WHERE TWO SUB-JETS WERE JOINED ARE RECORDED
0434 C   REMEMBER THAT A SUB-JET IS DEFINED AS A JET AT Y=YCUT WHICH HAS NOT
0435 C   YET BEEN MERGED WITH THE BEAM AT Y=YMAC.
0436 C
0437 C   ECUT    = INPUT : DENOMINATOR OF KT MEASURE. IF ZERO, ETOT IS USED
0438 C   YMAC    = INPUT : VALUE OF Y USED TO DEFINE MACRO-JETS
0439 C   Y(J)    =OUTPUT : Y VALUE WHERE EVENT CHANGED FROM HAVING
0440 C                         N+J SUB-JETS TO HAVING N+J-1, WHERE N IS
0441 C                         THE NUMBER OF MACRO-JETS AT SCALE YMAC
0442 C   LAST ARGUMENT IS LABEL TO JUMP TO IF FOR ANY REASON THE EVENT
0443 C   COULD NOT BE PROCESSED
0444 C
0445 C   NOTE THAT ALL FLOATING POINT VARIABLES ARE DECLARED DOUBLE PRECISION
0446 C
0447       INTEGER NMAX,HIST,NUM,I,J
0448       PARAMETER (NMAX=512)
0449       DOUBLE PRECISION ETOT,RSQ,P,KT,KTP,KTS,ECUT,ETSQ,Y(*),YMAC,KTLAST,
0450      &     ROUND
0451       PARAMETER (ROUND=0.99999D0)
0452       COMMON /KTCOMM/ETOT,RSQ,P(9,NMAX),KTP(NMAX,NMAX),KTS(NMAX),
0453      &  KT(NMAX),KTLAST(NMAX),HIST(NMAX),NUM
0454       IF (ETOT.EQ.0) CALL KTWARN('KTJOIN',100,*999)
0455       IF (ECUT.EQ.0) THEN
0456          ETSQ=1/ETOT**2
0457       ELSE
0458          ETSQ=1/ECUT**2
0459       ENDIF
0460       J=1
0461       DO 100 I=1,NUM
0462          IF (HIST(I).GT.NMAX.AND.ETSQ*KTLAST(I).GE.ROUND*YMAC) THEN
0463             Y(J)=ETSQ*KT(I)
0464             J=J+1
0465          ENDIF
0466  100  CONTINUE
0467       DO 200 I=J,NUM
0468          Y(I)=0
0469  200  CONTINUE
0470       RETURN
0471  999  RETURN 1
0472       END
0473 C-----------------------------------------------------------------------
0474       SUBROUTINE KTRECO(RECO,PP,NN,ECUT,YCUT,YMAC,PJET,JET,NJET,NSUB,*)
0475       IMPLICIT NONE
0476 C---RECONSTRUCT KINEMATICS OF JET SYSTEM, WHICH HAS ALREADY BEEN
0477 C   ANALYSED BY KTCLUS. NOTE THAT NO CONSISTENCY CHECK IS MADE: USER
0478 C   IS TRUSTED TO USE THE SAME PP VALUES AS FOR KTCLUS
0479 C
0480 C   RECO     = INPUT : RECOMBINATION SCHEME (NEED NOT BE SAME AS KTCLUS)
0481 C   PP(I,J)  = INPUT : 4-MOMENTUM OF Jth PARTICLE: I=1,4 => PX,PY,PZ,E
0482 C   NN       = INPUT : NUMBER OF PARTICLES
0483 C   ECUT     = INPUT : DENOMINATOR OF KT MEASURE. IF ZERO, ETOT IS USED
0484 C   YCUT     = INPUT : Y VALUE AT WHICH TO RECONSTRUCT JET MOMENTA
0485 C   YMAC     = INPUT : Y VALUE USED TO DEFINE MACRO-JETS, TO DETERMINE
0486 C                        WHICH JETS ARE SUB-JETS
0487 C   PJET(I,J)=OUTPUT : 4-MOMENTUM OF Jth JET AT SCALE YCUT
0488 C   JET(J)   =OUTPUT : THE MACRO-JET WHICH CONTAINS THE Jth JET,
0489 C                        SET TO ZERO IF JET IS NOT A SUB-JET
0490 C   NJET     =OUTPUT : THE NUMBER OF JETS
0491 C   NSUB     =OUTPUT : THE NUMBER OF SUB-JETS (EQUAL TO THE NUMBER OF
0492 C                        NON-ZERO ENTRIES IN JET())
0493 C   LAST ARGUMENT IS LABEL TO JUMP TO IF FOR ANY REASON THE EVENT
0494 C   COULD NOT BE PROCESSED
0495 C
0496 C   NOTE THAT THE MOMENTA ARE DECLARED DOUBLE PRECISION,
0497 C   AND ALL OTHER FLOATING POINT VARIABLES ARE DECLARED DOUBLE PRECISION
0498 C
0499       INTEGER NMAX,RECO,NUM,N,NN,NJET,NSUB,JET(*),HIST,IMIN,JMIN,I,J
0500       PARAMETER (NMAX=512)
0501       DOUBLE PRECISION PP(4,*),PJET(4,*)
0502       DOUBLE PRECISION ECUT,P,KT,KTP,KTS,ETOT,RSQ,ETSQ,YCUT,YMAC,KTLAST,
0503      &     ROUND
0504       PARAMETER (ROUND=0.99999D0)
0505       COMMON /KTCOMM/ETOT,RSQ,P(9,NMAX),KTP(NMAX,NMAX),KTS(NMAX),
0506      &  KT(NMAX),KTLAST(NMAX),HIST(NMAX),NUM
0507 C---CHECK INPUT
0508       IF (RECO.LT.1.OR.RECO.GT.3) THEN
0509         PRINT *,'RECO=',RECO
0510         CALL KTWARN('KTRECO',100,*999)
0511       ENDIF
0512 C---COPY PP TO P
0513       N=NN
0514       IF (NUM.NE.NN) CALL KTWARN('KTRECO',101,*999)
0515       CALL KTCOPY(PP,N,P,(RECO.NE.1))
0516       IF (ECUT.EQ.0) THEN
0517          ETSQ=1/ETOT**2
0518       ELSE
0519          ETSQ=1/ECUT**2
0520       ENDIF
0521 C---KEEP MERGING UNTIL YCUT
0522  100  IF (ETSQ*KT(N).LT.ROUND*YCUT) THEN
0523          IF (HIST(N).LE.NMAX) THEN
0524             CALL KTMOVE(P,KTP,KTS,NMAX,N,HIST(N),0)
0525          ELSE
0526             IMIN=HIST(N)/NMAX
0527             JMIN=HIST(N)-IMIN*NMAX
0528             CALL KTMERG(P,KTP,KTS,NMAX,IMIN,JMIN,N,0,0,0,RECO)
0529             CALL KTMOVE(P,KTP,KTS,NMAX,N,JMIN,0)
0530          ENDIF
0531          N=N-1
0532          IF (N.GT.0) GOTO 100
0533       ENDIF
0534 C---IF YCUT IS TOO LARGE THERE ARE NO JETS
0535       NJET=N
0536       NSUB=N
0537       IF (N.EQ.0) RETURN
0538 C---SET UP OUTPUT MOMENTA
0539       DO 210 I=1,NJET
0540          IF (RECO.EQ.1) THEN
0541             DO 200 J=1,4
0542                PJET(J,I)=P(J,I)
0543  200        CONTINUE
0544          ELSE
0545             PJET(1,I)=P(6,I)*COS(P(8,I))
0546             PJET(2,I)=P(6,I)*SIN(P(8,I))
0547             PJET(3,I)=P(6,I)*SINH(P(7,I))
0548             PJET(4,I)=P(6,I)*COSH(P(7,I))
0549          ENDIF
0550          JET(I)=I
0551  210  CONTINUE
0552 C---KEEP MERGING UNTIL YMAC TO FIND THE FATE OF EACH JET
0553  300  IF (ETSQ*KT(N).LT.ROUND*YMAC) THEN
0554          IF (HIST(N).LE.NMAX) THEN
0555             IMIN=0
0556             JMIN=HIST(N)
0557             NSUB=NSUB-1
0558          ELSE
0559             IMIN=HIST(N)/NMAX
0560             JMIN=HIST(N)-IMIN*NMAX
0561             IF (ETSQ*KTLAST(N).LT.ROUND*YMAC) NSUB=NSUB-1
0562          ENDIF
0563          DO 310 I=1,NJET
0564             IF (JET(I).EQ.JMIN) JET(I)=IMIN
0565             IF (JET(I).EQ.N) JET(I)=JMIN
0566  310     CONTINUE
0567          N=N-1
0568          IF (N.GT.0) GOTO 300
0569       ENDIF
0570       RETURN
0571  999  RETURN 1
0572       END
0573 C-----------------------------------------------------------------------
0574       SUBROUTINE KTINCL(RECO,PP,NN,PJET,JET,NJET,*)
0575       IMPLICIT NONE
0576 C---RECONSTRUCT KINEMATICS OF JET SYSTEM, WHICH HAS ALREADY BEEN
0577 C   ANALYSED BY KTCLUS ACCORDING TO THE INCLUSIVE JET DEFINITION. NOTE
0578 C   THAT NO CONSISTENCY CHECK IS MADE: USER IS TRUSTED TO USE THE SAME
0579 C   PP VALUES AS FOR KTCLUS
0580 C
0581 C   RECO     = INPUT : RECOMBINATION SCHEME (NEED NOT BE SAME AS KTCLUS)
0582 C   PP(I,J)  = INPUT : 4-MOMENTUM OF Jth PARTICLE: I=1,4 => PX,PY,PZ,E
0583 C   NN       = INPUT : NUMBER OF PARTICLES
0584 C   PJET(I,J)=OUTPUT : 4-MOMENTUM OF Jth JET AT SCALE YCUT
0585 C   JET(J)   =OUTPUT : THE JET WHICH CONTAINS THE Jth PARTICLE
0586 C   NJET     =OUTPUT : THE NUMBER OF JETS
0587 C   LAST ARGUMENT IS LABEL TO JUMP TO IF FOR ANY REASON THE EVENT
0588 C   COULD NOT BE PROCESSED
0589 C
0590 C   NOTE THAT THE MOMENTA ARE DECLARED DOUBLE PRECISION,
0591 C   AND ALL OTHER FLOATING POINT VARIABLES ARE DECLARED DOUBLE PRECISION
0592 C
0593       INTEGER NMAX,RECO,NUM,N,NN,NJET,JET(*),HIST,IMIN,JMIN,I,J
0594       PARAMETER (NMAX=512)
0595       DOUBLE PRECISION PP(4,*),PJET(4,*)
0596       DOUBLE PRECISION P,KT,KTP,KTS,ETOT,RSQ,KTLAST
0597       COMMON /KTCOMM/ETOT,RSQ,P(9,NMAX),KTP(NMAX,NMAX),KTS(NMAX),
0598      &  KT(NMAX),KTLAST(NMAX),HIST(NMAX),NUM
0599 C---CHECK INPUT
0600       IF (RECO.LT.1.OR.RECO.GT.3) CALL KTWARN('KTINCL',100,*999)
0601 C---COPY PP TO P
0602       N=NN
0603       IF (NUM.NE.NN) CALL KTWARN('KTINCL',101,*999)
0604       CALL KTCOPY(PP,N,P,(RECO.NE.1))
0605 C---INITIALLY EVERY PARTICLE IS IN ITS OWN JET
0606       DO 100 I=1,NN
0607          JET(I)=I
0608  100  CONTINUE
0609 C---KEEP MERGING TO THE BITTER END
0610       NJET=0
0611  200  IF (N.GT.0) THEN
0612          IF (HIST(N).LE.NMAX) THEN
0613             IMIN=0
0614             JMIN=HIST(N)
0615             NJET=NJET+1
0616             IF (RECO.EQ.1) THEN
0617                DO 300 J=1,4
0618                   PJET(J,NJET)=P(J,JMIN)
0619  300           CONTINUE
0620             ELSE
0621                PJET(1,NJET)=P(6,JMIN)*COS(P(8,JMIN))
0622                PJET(2,NJET)=P(6,JMIN)*SIN(P(8,JMIN))
0623                PJET(3,NJET)=P(6,JMIN)*SINH(P(7,JMIN))
0624                PJET(4,NJET)=P(6,JMIN)*COSH(P(7,JMIN))
0625             ENDIF
0626             CALL KTMOVE(P,KTP,KTS,NMAX,N,JMIN,0)
0627          ELSE
0628             IMIN=HIST(N)/NMAX
0629             JMIN=HIST(N)-IMIN*NMAX
0630             CALL KTMERG(P,KTP,KTS,NMAX,IMIN,JMIN,N,0,0,0,RECO)
0631             CALL KTMOVE(P,KTP,KTS,NMAX,N,JMIN,0)
0632          ENDIF
0633          DO 400 I=1,NN
0634             IF (JET(I).EQ.JMIN) JET(I)=IMIN
0635             IF (JET(I).EQ.N) JET(I)=JMIN
0636             IF (JET(I).EQ.0) JET(I)=-NJET
0637  400     CONTINUE
0638          N=N-1
0639          GOTO 200
0640       ENDIF
0641 C---FINALLY EVERY PARTICLE MUST BE IN AN INCLUSIVE JET
0642       DO 500 I=1,NN
0643 C---IF THERE ARE ANY UNASSIGNED PARTICLES SOMETHING MUST HAVE GONE WRONG
0644          IF (JET(I).GE.0) CALL KTWARN('KTINCL',102,*999)
0645          JET(I)=-JET(I)
0646  500  CONTINUE
0647       RETURN
0648  999  RETURN 1
0649       END
0650 C-----------------------------------------------------------------------
0651       SUBROUTINE KTISUB(N,NY,YCUT,NSUB,*)
0652       IMPLICIT NONE
0653 C---COUNT THE NUMBER OF SUB-JETS IN THE Nth INCLUSIVE JET OF AN EVENT
0654 C   THAT HAS ALREADY BEEN ANALYSED BY KTCLUS.
0655 C
0656 C   N       = INPUT : WHICH INCLUSIVE JET TO USE
0657 C   NY      = INPUT : NUMBER OF YCUT VALUES
0658 C   YCUT(J) = INPUT : Y VALUES AT WHICH NUMBERS OF SUB-JETS ARE COUNTED
0659 C   NSUB(J) =OUTPUT : NUMBER OF SUB-JETS AT YCUT(J)
0660 C   LAST ARGUMENT IS LABEL TO JUMP TO IF FOR ANY REASON THE EVENT
0661 C   COULD NOT BE PROCESSED
0662 C
0663 C   NOTE THAT ALL FLOATING POINT VARIABLES ARE DECLARED DOUBLE PRECISION
0664 C
0665       INTEGER N,NY,NSUB(NY),NMAX,HIST,I,J,NUM,NM
0666       PARAMETER (NMAX=512)
0667       DOUBLE PRECISION YCUT(NY),ETOT,RSQ,P,KT,KTP,KTS,KTLAST,ROUND,EPS
0668       PARAMETER (ROUND=0.99999D0)
0669       COMMON /KTCOMM/ETOT,RSQ,P(9,NMAX),KTP(NMAX,NMAX),KTS(NMAX),
0670      &  KT(NMAX),KTLAST(NMAX),HIST(NMAX),NUM
0671       DATA EPS/1D-6/
0672       DO 100 I=1,NY
0673          NSUB(I)=0
0674  100  CONTINUE
0675 C---FIND WHICH MERGING CORRESPONDS TO THE NTH INCLUSIVE JET
0676       NM=0
0677       J=0
0678       DO 110 I=NUM,1,-1
0679         IF (HIST(I).LE.NMAX) J=J+1
0680         IF (J.EQ.N) THEN
0681           NM=I
0682           GOTO 120
0683         ENDIF
0684  110  CONTINUE
0685  120  CONTINUE
0686 C---GIVE UP IF THERE ARE LESS THAN N INCLUSIVE JETS
0687       IF (NM.EQ.0) CALL KTWARN('KTISUB',100,*999)
0688       DO 210 I=NUM,1,-1
0689          DO 200 J=1,NY
0690             IF (NSUB(J).EQ.0.AND.RSQ*KT(I).GE.ROUND*YCUT(J)*KT(NM))
0691      &          NSUB(J)=I
0692             IF (NSUB(J).NE.0.AND.ABS(KTLAST(I)-KTLAST(NM)).GT.EPS)
0693      &          NSUB(J)=NSUB(J)-1
0694  200     CONTINUE
0695  210  CONTINUE
0696       RETURN
0697  999  RETURN 1
0698       END
0699 C-----------------------------------------------------------------------
0700       SUBROUTINE KTIJOI(N,Y,*)
0701       IMPLICIT NONE
0702 C---GIVE SAME INFORMATION AS LAST CALL TO KTCLUS EXCEPT THAT ONLY
0703 C   MERGES OF TWO SUB-JETS INSIDE THE Nth INCLUSIVE JET ARE RECORDED
0704 C
0705 C   N       = INPUT : WHICH INCLUSIVE JET TO USE
0706 C   Y(J)    =OUTPUT : Y VALUE WHERE JET CHANGED FROM HAVING
0707 C                         J+1 SUB-JETS TO HAVING J
0708 C   LAST ARGUMENT IS LABEL TO JUMP TO IF FOR ANY REASON THE EVENT
0709 C   COULD NOT BE PROCESSED
0710 C
0711 C   NOTE THAT ALL FLOATING POINT VARIABLES ARE DECLARED DOUBLE PRECISION
0712 C
0713       INTEGER NMAX,HIST,NUM,I,J,N,NM
0714       PARAMETER (NMAX=512)
0715       DOUBLE PRECISION ETOT,RSQ,P,KT,KTP,KTS,Y(*),KTLAST,EPS
0716       COMMON /KTCOMM/ETOT,RSQ,P(9,NMAX),KTP(NMAX,NMAX),KTS(NMAX),
0717      &  KT(NMAX),KTLAST(NMAX),HIST(NMAX),NUM
0718       DATA EPS/1D-6/
0719 C---FIND WHICH MERGING CORRESPONDS TO THE NTH INCLUSIVE JET
0720       NM=0
0721       J=0
0722       DO 100 I=NUM,1,-1
0723         IF (HIST(I).LE.NMAX) J=J+1
0724         IF (J.EQ.N) THEN
0725           NM=I
0726           GOTO 105
0727         ENDIF
0728  100  CONTINUE
0729  105  CONTINUE
0730 C---GIVE UP IF THERE ARE LESS THAN N INCLUSIVE JETS
0731       IF (NM.EQ.0) CALL KTWARN('KTIJOI',100,*999)
0732       J=1
0733       DO 110 I=1,NUM
0734          IF (HIST(I).GT.NMAX.AND.ABS(KTLAST(I)-KTLAST(NM)).LT.EPS) THEN
0735             Y(J)=RSQ*KT(I)/KT(NM)
0736             J=J+1
0737          ENDIF
0738  110  CONTINUE
0739       DO 200 I=J,NUM
0740          Y(I)=0
0741  200  CONTINUE
0742       RETURN
0743  999  RETURN 1
0744       END
0745 C-----------------------------------------------------------------------
0746       SUBROUTINE KTIREC(RECO,PP,NN,N,YCUT,PSUB,NSUB,*)
0747       IMPLICIT NONE
0748 C---RECONSTRUCT KINEMATICS OF SUB-JET SYSTEM IN THE Nth INCLUSIVE JET
0749 C   OF AN EVENT THAT HAS ALREADY BEEN ANALYSED BY KTCLUS
0750 C
0751 C   RECO     = INPUT : RECOMBINATION SCHEME (NEED NOT BE SAME AS KTCLUS)
0752 C   PP(I,J)  = INPUT : 4-MOMENTUM OF Jth PARTICLE: I=1,4 => PX,PY,PZ,E
0753 C   NN       = INPUT : NUMBER OF PARTICLES
0754 C   N        = INPUT : WHICH INCLUSIVE JET TO USE
0755 C   YCUT     = INPUT : Y VALUE AT WHICH TO RECONSTRUCT JET MOMENTA
0756 C   PSUB(I,J)=OUTPUT : 4-MOMENTUM OF Jth SUB-JET AT SCALE YCUT
0757 C   NSUB     =OUTPUT : THE NUMBER OF SUB-JETS
0758 C   LAST ARGUMENT IS LABEL TO JUMP TO IF FOR ANY REASON THE EVENT
0759 C   COULD NOT BE PROCESSED
0760 C
0761 C   NOTE THAT THE MOMENTA ARE DECLARED DOUBLE PRECISION,
0762 C   AND ALL OTHER FLOATING POINT VARIABLES ARE DECLARED DOUBLE PRECISION
0763 C
0764       INTEGER NMAX,RECO,NUM,NN,NJET,NSUB,JET,HIST,I,J,N,NM
0765       PARAMETER (NMAX=512)
0766       DOUBLE PRECISION PP(4,*),PSUB(4,*)
0767       DOUBLE PRECISION ECUT,P,KT,KTP,KTS,ETOT,RSQ,YCUT,YMAC,KTLAST
0768       COMMON /KTCOMM/ETOT,RSQ,P(9,NMAX),KTP(NMAX,NMAX),KTS(NMAX),
0769      &  KT(NMAX),KTLAST(NMAX),HIST(NMAX),NUM
0770       DIMENSION JET(NMAX)
0771 C---FIND WHICH MERGING CORRESPONDS TO THE NTH INCLUSIVE JET
0772       NM=0
0773       J=0
0774       DO 100 I=NUM,1,-1
0775          IF (HIST(I).LE.NMAX) J=J+1
0776          IF (J.EQ.N) THEN
0777             NM=I
0778             GOTO 110
0779          ENDIF
0780  100  CONTINUE
0781  110  CONTINUE
0782 C---GIVE UP IF THERE ARE LESS THAN N INCLUSIVE JETS
0783       IF (NM.EQ.0) CALL KTWARN('KTIREC',102,*999)
0784 C---RECONSTRUCT THE JETS AT THE APPROPRIATE SCALE
0785       ECUT=SQRT(KT(NM)/RSQ)
0786       YMAC=RSQ
0787       CALL KTRECO(RECO,PP,NN,ECUT,YCUT,YMAC,PSUB,JET,NJET,NSUB,*999)
0788 C---GET RID OF THE ONES THAT DO NOT END UP IN THE JET WE WANT
0789       NSUB=0
0790       DO 210 I=1,NJET
0791          IF (JET(I).EQ.HIST(NM)) THEN
0792             NSUB=NSUB+1
0793             DO 200 J=1,4
0794                PSUB(J,NSUB)=PSUB(J,I)
0795  200        CONTINUE
0796          ENDIF
0797  210  CONTINUE
0798       RETURN
0799  999  RETURN 1
0800       END
0801 C-----------------------------------------------------------------------
0802       SUBROUTINE KTWICH(ECUT,YCUT,JET,NJET,*)
0803       IMPLICIT NONE
0804 C---GIVE A LIST OF WHICH JET EACH ORIGINAL PARTICLE ENDED UP IN AT SCALE
0805 C   YCUT, TOGETHER WITH THE NUMBER OF JETS AT THAT SCALE.
0806 C
0807 C   ECUT     = INPUT : DENOMINATOR OF KT MEASURE. IF ZERO, ETOT IS USED
0808 C   YCUT     = INPUT : Y VALUE AT WHICH TO DEFINE JETS
0809 C   JET(J)   =OUTPUT : THE JET WHICH CONTAINS THE Jth PARTICLE,
0810 C                        SET TO ZERO IF IT WAS PUT INTO THE BEAM JETS
0811 C   NJET     =OUTPUT : THE NUMBER OF JETS AT SCALE YCUT (SO JET()
0812 C                        ENTRIES WILL BE IN THE RANGE 0 -> NJET)
0813 C   LAST ARGUMENT IS LABEL TO JUMP TO IF FOR ANY REASON THE EVENT
0814 C   COULD NOT BE PROCESSED
0815 C
0816 C   NOTE THAT ALL FLOATING POINT VARIABLES ARE DECLARED DOUBLE PRECISION
0817 C
0818       INTEGER JET(*),NJET,NTEMP
0819       DOUBLE PRECISION ECUT,YCUT
0820       CALL KTWCHS(ECUT,YCUT,YCUT,JET,NJET,NTEMP,*999)
0821       RETURN
0822  999  RETURN 1
0823       END
0824 C-----------------------------------------------------------------------
0825       SUBROUTINE KTWCHS(ECUT,YCUT,YMAC,JET,NJET,NSUB,*)
0826       IMPLICIT NONE
0827 C---GIVE A LIST OF WHICH SUB-JET EACH ORIGINAL PARTICLE ENDED UP IN AT
0828 C   SCALE YCUT, WITH MACRO-JET SCALE YMAC, TOGETHER WITH THE NUMBER OF
0829 C   JETS AT SCALE YCUT AND THE NUMBER OF THEM WHICH ARE SUB-JETS.
0830 C
0831 C   ECUT     = INPUT : DENOMINATOR OF KT MEASURE. IF ZERO, ETOT IS USED
0832 C   YCUT     = INPUT : Y VALUE AT WHICH TO DEFINE JETS
0833 C   YMAC     = INPUT : Y VALUE AT WHICH TO DEFINE MACRO-JETS
0834 C   JET(J)   =OUTPUT : THE JET WHICH CONTAINS THE Jth PARTICLE,
0835 C                        SET TO ZERO IF IT WAS PUT INTO THE BEAM JETS
0836 C   NJET     =OUTPUT : THE NUMBER OF JETS AT SCALE YCUT (SO JET()
0837 C                        ENTRIES WILL BE IN THE RANGE 0 -> NJET)
0838 C   NSUB     =OUTPUT : THE NUMBER OF SUB-JETS AT SCALE YCUT, WITH
0839 C                        MACRO-JETS DEFINED AT SCALE YMAC (SO ONLY NSUB
0840 C                        OF THE JETS 1 -> NJET WILL APPEAR IN JET())
0841 C   LAST ARGUMENT IS LABEL TO JUMP TO IF FOR ANY REASON THE EVENT
0842 C   COULD NOT BE PROCESSED
0843 C
0844 C   NOTE THAT ALL FLOATING POINT VARIABLES ARE DECLARED DOUBLE PRECISION
0845 C
0846       INTEGER NMAX,JET(*),NJET,NSUB,HIST,NUM,I,J,JSUB
0847       PARAMETER (NMAX=512)
0848       DOUBLE PRECISION P1(4,NMAX),P2(4,NMAX)
0849       DOUBLE PRECISION ECUT,YCUT,YMAC,ZERO,ETOT,RSQ,P,KTP,KTS,KT,KTLAST
0850       COMMON /KTCOMM/ETOT,RSQ,P(9,NMAX),KTP(NMAX,NMAX),KTS(NMAX),
0851      &  KT(NMAX),KTLAST(NMAX),HIST(NMAX),NUM
0852       DIMENSION JSUB(NMAX)
0853 C---THE MOMENTA HAVE TO BEEN GIVEN LEGAL VALUES,
0854 C   EVEN THOUGH THEY WILL NEVER BE USED
0855       DATA ((P1(J,I),I=1,NMAX),J=1,4),ZERO
0856      &  /NMAX*1,NMAX*0,NMAX*0,NMAX*1,0/
0857 C---FIRST GET A LIST OF WHICH PARTICLE IS IN WHICH JET AT YCUT
0858       CALL KTRECO(1,P1,NUM,ECUT,ZERO,YCUT,P2,JET,NJET,NSUB,*999)
0859 C---THEN FIND OUT WHICH JETS ARE SUBJETS
0860       CALL KTRECO(1,P1,NUM,ECUT,YCUT,YMAC,P2,JSUB,NJET,NSUB,*999)
0861 C---AND MODIFY JET() ACCORDINGLY
0862       DO 10 I=1,NUM
0863         IF (JET(I).NE.0) THEN
0864           IF (JSUB(JET(I)).EQ.0) JET(I)=0
0865         ENDIF
0866  10   CONTINUE
0867       RETURN
0868  999  RETURN 1
0869       END
0870 C-----------------------------------------------------------------------
0871       SUBROUTINE KTFRAM(IOPT,CMF,SIGN,Z,XZ,N,P,Q,*)
0872       IMPLICIT NONE
0873 C---BOOST PARTICLES IN P TO/FROM FRAME GIVEN BY CMF, Z, XZ.
0874 C---IN THIS FRAME CMZ IS STATIONARY,
0875 C                   Z IS ALONG THE (SIGN)Z-AXIS (SIGN=+ OR -)
0876 C                  XZ IS IN THE X-Z PLANE (WITH POSITIVE X COMPONENT)
0877 C---IF Z HAS LENGTH ZERO, OR SIGN=0, NO ROTATION IS PERFORMED
0878 C---IF XZ HAS ZERO COMPONENT PERPENDICULAR TO Z IN THAT FRAME,
0879 C   NO AZIMUTHAL ROTATION IS PERFORMED
0880 C
0881 C   IOPT    = INPUT  : 0=TO FRAME, 1=FROM FRAME
0882 C   CMF(I)  = INPUT  : 4-MOMENTUM WHICH IS STATIONARY IN THE FRAME
0883 C   SIGN    = INPUT  : DIRECTION OF Z IN THE FRAME, NOTE THAT
0884 C                        ONLY ITS SIGN IS USED, NOT ITS MAGNITUDE
0885 C   Z(I)    = INPUT  : 4-MOMENTUM WHICH LIES ON THE (SIGN)Z-AXIS
0886 C   XZ(I)   = INPUT  : 4-MOMENTUM WHICH LIES IN THE X-Z PLANE
0887 C   N       = INPUT  : NUMBER OF PARTICLES IN P
0888 C   P(I,J)  = INPUT  : 4-MOMENTUM OF JTH PARTICLE BEFORE
0889 C   Q(I,J)  = OUTPUT : 4-MOMENTUM OF JTH PARTICLE AFTER
0890 C   LAST ARGUMENT IS LABEL TO JUMP TO IF FOR ANY REASON THE EVENT
0891 C   COULD NOT BE PROCESSED
0892 C
0893 C   NOTE THAT ALL MOMENTA ARE DOUBLE PRECISION
0894 C
0895 C   NOTE THAT IT IS SAFE TO CALL WITH P=Q
0896 C
0897       INTEGER IOPT,I,N
0898       DOUBLE PRECISION CMF(4),SIGN,Z(4),XZ(4),P(4,N),Q(4,N),
0899      &  R(4,4),NEW(4),OLD(4)
0900       IF (IOPT.LT.0.OR.IOPT.GT.1) CALL KTWARN('KTFRAM',200,*999)
0901 C---FIND BOOST TO GET THERE FROM LAB
0902       CALL KTUNIT(R)
0903       CALL KTLBST(0,R,CMF,*999)
0904 C---FIND ROTATION TO PUT BOOSTED Z ON THE (SIGN)Z AXIS
0905       IF (SIGN.NE.0) THEN
0906         CALL KTVMUL(R,Z,OLD)
0907         IF (OLD(1).NE.0.OR.OLD(2).NE.0.OR.OLD(3).NE.0) THEN
0908           NEW(1)=0
0909           NEW(2)=0
0910           NEW(3)=SIGN
0911           NEW(4)=ABS(SIGN)
0912           CALL KTRROT(R,OLD,NEW,*999)
0913 C---FIND ROTATION TO PUT BOOSTED AND ROTATED XZ INTO X-Z PLANE
0914           CALL KTVMUL(R,XZ,OLD)
0915           IF (OLD(1).NE.0.OR.OLD(2).NE.0) THEN
0916             NEW(1)=1
0917             NEW(2)=0
0918             NEW(3)=0
0919             NEW(4)=1
0920             OLD(3)=0
0921 C---NOTE THAT A POTENTIALLY AWKWARD SPECIAL CASE IS AVERTED, BECAUSE IF
0922 C   OLD AND NEW ARE EXACTLY BACK-TO-BACK, THE ROTATION AXIS IS UNDEFINED
0923 C   BUT IN THAT CASE KTRROT WILL USE THE Z AXIS, AS REQUIRED
0924             CALL KTRROT(R,OLD,NEW,*999)
0925           ENDIF
0926         ENDIF
0927       ENDIF
0928 C---INVERT THE TRANSFORMATION IF NECESSARY
0929       IF (IOPT.EQ.1) CALL KTINVT(R,R)
0930 C---APPLY THE RESULT TO ALL THE VECTORS
0931       DO 30 I=1,N
0932         CALL KTVMUL(R,P(1,I),Q(1,I))
0933  30   CONTINUE
0934       RETURN
0935  999  RETURN 1
0936       END
0937 C-----------------------------------------------------------------------
0938       SUBROUTINE KTBREI(IOPT,PLEP,PHAD,POUT,N,P,Q,*)
0939       IMPLICIT NONE
0940 C---BOOST PARTICLES IN P TO/FROM BREIT FRAME
0941 C
0942 C   IOPT    = INPUT  : 0/2=TO BREIT FRAME, 1/3=FROM BREIT FRAME
0943 C                      0/1=NO AZIMUTHAL ROTATION AFTERWARDS
0944 C                      2/3=LEPTON PLANE ROTATED INTO THE X-Z PLANE
0945 C   PLEP    = INPUT  : MOMENTUM OF INCOMING LEPTON IN +Z DIRECTION
0946 C   PHAD    = INPUT  : MOMENTUM OF INCOMING HADRON IN +Z DIRECTION
0947 C   POUT(I) = INPUT  : 4-MOMENTUM OF OUTGOING LEPTON
0948 C   N       = INPUT  : NUMBER OF PARTICLES IN P
0949 C   P(I,J)  = INPUT  : 4-MOMENTUM OF JTH PARTICLE BEFORE
0950 C   Q(I,J)  = OUTPUT : 4-MOMENTUM OF JTH PARTICLE AFTER
0951 C   LAST ARGUMENT IS LABEL TO JUMP TO IF FOR ANY REASON THE EVENT
0952 C   COULD NOT BE PROCESSED (MOST LIKELY DUE TO PARTICLES HAVING SMALLER
0953 C   ENERGY THAN MOMENTUM)
0954 C
0955 C   NOTE THAT ALL MOMENTA ARE DOUBLE PRECISION
0956 C
0957 C   NOTE THAT IT IS SAFE TO CALL WITH P=Q
0958 C
0959       INTEGER IOPT,N
0960       DOUBLE PRECISION PLEP,PHAD,POUT(4),P(4,N),Q(4,N),
0961      &  CMF(4),Z(4),XZ(4),DOT,QDQ
0962 C---CHECK INPUT
0963       IF (IOPT.LT.0.OR.IOPT.GT.3) CALL KTWARN('KTBREI',200,*999)
0964 C---FIND 4-MOMENTUM OF BREIT FRAME (TIMES AN ARBITRARY FACTOR)
0965       DOT=ABS(PHAD)*(ABS(PLEP)-POUT(4))-PHAD*(PLEP-POUT(3))
0966       QDQ=(ABS(PLEP)-POUT(4))**2-(PLEP-POUT(3))**2-POUT(2)**2-POUT(1)**2
0967       CMF(1)=DOT*(         -POUT(1))
0968       CMF(2)=DOT*(         -POUT(2))
0969       CMF(3)=DOT*(    PLEP -POUT(3))-QDQ*    PHAD
0970       CMF(4)=DOT*(ABS(PLEP)-POUT(4))-QDQ*ABS(PHAD)
0971 C---FIND ROTATION TO PUT INCOMING HADRON BACK ON Z-AXIS
0972       Z(1)=0
0973       Z(2)=0
0974       Z(3)=PHAD
0975       Z(4)=ABS(PHAD)
0976       XZ(1)=0
0977       XZ(2)=0
0978       XZ(3)=0
0979       XZ(4)=0
0980 C---DO THE BOOST
0981       IF (IOPT.LE.1) THEN
0982         CALL KTFRAM(IOPT,CMF,PHAD,Z,XZ,N,P,Q,*999)
0983       ELSE
0984         CALL KTFRAM(IOPT-2,CMF,PHAD,Z,POUT,N,P,Q,*999)
0985       ENDIF
0986       RETURN
0987  999  RETURN 1
0988       END
0989 C-----------------------------------------------------------------------
0990       SUBROUTINE KTHADR(IOPT,PLEP,PHAD,POUT,N,P,Q,*)
0991       IMPLICIT NONE
0992 C---BOOST PARTICLES IN P TO/FROM HADRONIC CMF
0993 C
0994 C   ARGUMENTS ARE EXACTLY AS FOR KTBREI
0995 C
0996 C   NOTE THAT ALL MOMENTA ARE DOUBLE PRECISION
0997 C
0998 C   NOTE THAT IT IS SAFE TO CALL WITH P=Q
0999 C
1000       INTEGER IOPT,N
1001       DOUBLE PRECISION PLEP,PHAD,POUT(4),P(4,N),Q(4,N),
1002      &  CMF(4),Z(4),XZ(4)
1003 C---CHECK INPUT
1004       IF (IOPT.LT.0.OR.IOPT.GT.3) CALL KTWARN('KTHADR',200,*999)
1005 C---FIND 4-MOMENTUM OF HADRONIC CMF
1006       CMF(1)=         -POUT(1)
1007       CMF(2)=         -POUT(2)
1008       CMF(3)=    PLEP -POUT(3)+    PHAD
1009       CMF(4)=ABS(PLEP)-POUT(4)+ABS(PHAD)
1010 C---FIND ROTATION TO PUT INCOMING HADRON BACK ON Z-AXIS
1011       Z(1)=0
1012       Z(2)=0
1013       Z(3)=PHAD
1014       Z(4)=ABS(PHAD)
1015       XZ(1)=0
1016       XZ(2)=0
1017       XZ(3)=0
1018       XZ(4)=0
1019 C---DO THE BOOST
1020       IF (IOPT.LE.1) THEN
1021         CALL KTFRAM(IOPT,CMF,PHAD,Z,XZ,N,P,Q,*999)
1022       ELSE
1023         CALL KTFRAM(IOPT-2,CMF,PHAD,Z,POUT,N,P,Q,*999)
1024       ENDIF
1025       RETURN
1026  999  RETURN 1
1027       END
1028 C-----------------------------------------------------------------------
1029       FUNCTION KTPAIR(ANGL,P,Q,ANGLE)
1030       IMPLICIT NONE
1031 C---CALCULATE LOCAL KT OF PAIR, USING ANGULAR SCHEME:
1032 C   1=>ANGULAR, 2=>DeltaR, 3=>f(DeltaEta,DeltaPhi)
1033 C   WHERE f(eta,phi)=2(COSH(eta)-COS(phi)) IS THE QCD EMISSION METRIC
1034 C---IF ANGLE<0, IT IS SET TO THE ANGULAR PART OF THE LOCAL KT ON RETURN
1035 C   IF ANGLE>0, IT IS USED INSTEAD OF THE ANGULAR PART OF THE LOCAL KT
1036       INTEGER ANGL
1037       DOUBLE PRECISION P(9),Q(9),KTPAIR,R,KTMDPI,ANGLE,ETA,PHI,ESQ
1038 C---COMPONENTS OF MOMENTA ARE PX,PY,PZ,E,1/P,PT,ETA,PHI,PT**2
1039       R=ANGLE
1040       IF (ANGL.EQ.1) THEN
1041          IF (R.LE.0) R=2*(1-(P(1)*Q(1)+P(2)*Q(2)+P(3)*Q(3))*(P(5)*Q(5)))
1042          ESQ=MIN(P(4),Q(4))**2
1043       ELSEIF (ANGL.EQ.2.OR.ANGL.EQ.3) THEN
1044          IF (R.LE.0) THEN
1045             ETA=P(7)-Q(7)
1046             PHI=KTMDPI(P(8)-Q(8))
1047             IF (ANGL.EQ.2) THEN
1048                R=ETA**2+PHI**2
1049             ELSE
1050                R=2*(COSH(ETA)-COS(PHI))
1051             ENDIF
1052          ENDIF
1053          ESQ=MIN(P(9),Q(9))
1054       ELSEIF (ANGL.EQ.4) THEN
1055         ESQ=(1d0/(P(5)*Q(5))-P(1)*Q(1)-P(2)*Q(2)-
1056      &P(3)*Q(3))*2D0/(P(5)*Q(5))/(0.0001D0+1d0/P(5)+1d0/Q(5))**2
1057         R=1d0
1058       ELSE
1059          ktpair = 0D0
1060          CALL KTWARN('KTPAIR',200,*999)
1061          STOP
1062       ENDIF
1063       KTPAIR=ESQ*R
1064       IF (ANGLE.LT.0) ANGLE=R
1065  999  END
1066 C-----------------------------------------------------------------------
1067       FUNCTION KTSING(ANGL,TYPE,P)
1068       IMPLICIT NONE
1069 C---CALCULATE KT OF PARTICLE, USING ANGULAR SCHEME:
1070 C   1=>ANGULAR, 2=>DeltaR, 3=>f(DeltaEta,DeltaPhi)
1071 C---TYPE=1 FOR E+E-, 2 FOR EP, 3 FOR PE, 4 FOR PP
1072 C   FOR EP, PROTON DIRECTION IS DEFINED AS -Z
1073 C   FOR PE, PROTON DIRECTION IS DEFINED AS +Z
1074       DOUBLE PRECISION KTSING,P(9)
1075       DOUBLE PRECISION COSTH,R,SMALL
1076       INTEGER ANGL,TYPE
1077       DATA SMALL/1D-4/
1078       IF (ANGL.EQ.1.OR.ANGL.EQ.4) THEN
1079          COSTH=P(3)*P(5)
1080          IF (TYPE.EQ.2) THEN
1081             COSTH=-COSTH
1082          ELSEIF (TYPE.EQ.4) THEN
1083             COSTH=ABS(COSTH)
1084          ELSEIF (TYPE.NE.1.AND.TYPE.NE.3) THEN
1085             ktsing = 0D0
1086             CALL KTWARN('KTSING',200,*999)
1087             STOP
1088          ENDIF
1089          R=2*(1-COSTH)
1090 C---IF CLOSE TO BEAM, USE APPROX 2*(1-COS(THETA))=SIN**2(THETA)
1091          IF (R.LT.SMALL) R=(P(1)**2+P(2)**2)*P(5)**2
1092          KTSING=P(4)**2*R
1093       ELSEIF (ANGL.EQ.2.OR.ANGL.EQ.3) THEN
1094          KTSING=P(9)
1095       ELSE
1096          ktsing = 0D0
1097          CALL KTWARN('KTSING',201,*999)
1098          STOP
1099       ENDIF
1100  999  END
1101 C-----------------------------------------------------------------------
1102       SUBROUTINE KTPMIN(A,NMAX,N,IMIN,JMIN)
1103       IMPLICIT NONE
1104 C---FIND THE MINIMUM MEMBER OF A(NMAX,NMAX) WITH IMIN < JMIN <= N
1105       INTEGER NMAX,N,IMIN,JMIN,KMIN,I,J,K
1106 C---REMEMBER THAT A(X+(Y-1)*NMAX)=A(X,Y)
1107 C   THESE LOOPING VARIABLES ARE J=Y-2, I=X+(Y-1)*NMAX
1108       DOUBLE PRECISION A(*),AMIN
1109       K=1+NMAX
1110       KMIN=K
1111       AMIN=A(KMIN)
1112       DO 110 J=0,N-2
1113          DO 100 I=K,K+J
1114             IF (A(I).LT.AMIN) THEN
1115                KMIN=I
1116                AMIN=A(KMIN)
1117             ENDIF
1118  100     CONTINUE
1119          K=K+NMAX
1120  110  CONTINUE
1121       JMIN=KMIN/NMAX+1
1122       IMIN=KMIN-(JMIN-1)*NMAX
1123       END
1124 C-----------------------------------------------------------------------
1125       SUBROUTINE KTSMIN(A,NMAX,N,IMIN)
1126       IMPLICIT NONE
1127 C---FIND THE MINIMUM MEMBER OF A
1128       INTEGER N,NMAX,IMIN,I
1129       DOUBLE PRECISION A(NMAX)
1130       IMIN=1
1131       DO 100 I=1,N
1132          IF (A(I).LT.A(IMIN)) IMIN=I
1133  100  CONTINUE
1134       END
1135 C-----------------------------------------------------------------------
1136       SUBROUTINE KTCOPY(A,N,B,ONSHLL)
1137       IMPLICIT NONE
1138 C---COPY FROM A TO B. 5TH=1/(3-MTM), 6TH=PT, 7TH=ETA, 8TH=PHI, 9TH=PT**2
1139 C   IF ONSHLL IS .TRUE. PARTICLE ENTRIES ARE PUT ON-SHELL BY SETTING E=P
1140       INTEGER I,N
1141       DOUBLE PRECISION A(4,N)
1142       LOGICAL ONSHLL
1143       DOUBLE PRECISION B(9,N),ETAMAX,SINMIN,EPS
1144       DATA ETAMAX,SINMIN,EPS/10,0,1D-6/
1145 C---SINMIN GETS CALCULATED ON FIRST CALL
1146       IF (SINMIN.EQ.0) SINMIN=1/COSH(ETAMAX)
1147       DO 100 I=1,N
1148          B(1,I)=A(1,I)
1149          B(2,I)=A(2,I)
1150          B(3,I)=A(3,I)
1151          B(4,I)=A(4,I)
1152          B(5,I)=SQRT(A(1,I)**2+A(2,I)**2+A(3,I)**2)
1153          IF (ONSHLL) B(4,I)=B(5,I)
1154          IF (B(5,I).EQ.0) B(5,I)=1D-10
1155          B(5,I)=1/B(5,I)
1156          B(9,I)=A(1,I)**2+A(2,I)**2
1157          B(6,I)=SQRT(B(9,I))
1158          B(7,I)=B(6,I)*B(5,I)
1159          IF (B(7,I).GT.SINMIN) THEN
1160             B(7,I)=A(4,I)**2-A(3,I)**2
1161             IF (B(7,I).LE.EPS*B(4,I)**2.OR.ONSHLL) B(7,I)=B(9,I)
1162             B(7,I)=LOG((B(4,I)+ABS(B(3,I)))**2/B(7,I))/2
1163          ELSE
1164             B(7,I)=ETAMAX+2
1165          ENDIF
1166          B(7,I)=SIGN(B(7,I),B(3,I))
1167          IF (A(1,I).EQ.0 .AND. A(2,I).EQ.0) THEN
1168             B(8,I)=0
1169          ELSE
1170             B(8,I)=ATAN2(A(2,I),A(1,I))
1171          ENDIF
1172  100  CONTINUE
1173       END
1174 C-----------------------------------------------------------------------
1175       SUBROUTINE KTMERG(P,KTP,KTS,NMAX,I,J,N,TYPE,ANGL,MONO,RECO)
1176       IMPLICIT NONE
1177 C---MERGE THE Jth PARTICLE IN P INTO THE Ith PARTICLE
1178 C   J IS ASSUMED GREATER THAN I. P CONTAINS N PARTICLES BEFORE MERGING.
1179 C---ALSO RECALCULATING THE CORRESPONDING KTP AND KTS VALUES IF MONO.GT.0
1180 C   FROM THE RECOMBINED ANGULAR MEASURES IF MONO.GT.1
1181 C---NOTE THAT IF MONO.LE.0, TYPE AND ANGL ARE NOT USED
1182       INTEGER ANGL,RECO,TYPE,I,J,K,N,NMAX,MONO
1183       DOUBLE PRECISION P(9,NMAX),KTP(NMAX,NMAX),KTS(NMAX),PT,PTT,
1184      &     KTMDPI,KTUP,PI,PJ,ANG,KTPAIR,KTSING,ETAMAX,EPS
1185       KTUP(I,J)=KTP(MAX(I,J),MIN(I,J))
1186       DATA ETAMAX,EPS/10,1D-6/
1187       IF (J.LE.I) CALL KTWARN('KTMERG',200,*999)
1188 C---COMBINE ANGULAR MEASURES IF NECESSARY
1189       IF (MONO.GT.1) THEN
1190          DO 100 K=1,N
1191             IF (K.NE.I.AND.K.NE.J) THEN
1192                IF (RECO.EQ.1) THEN
1193                   PI=P(4,I)
1194                   PJ=P(4,J)
1195                ELSEIF (RECO.EQ.2) THEN
1196                   PI=P(6,I)
1197                   PJ=P(6,J)
1198                ELSEIF (RECO.EQ.3) THEN
1199                   PI=P(9,I)
1200                   PJ=P(9,J)
1201                ELSE
1202                   CALL KTWARN('KTMERG',201,*999)
1203                   STOP
1204                ENDIF
1205                IF (PI.EQ.0.AND.PJ.EQ.0) THEN
1206                   PI=1
1207                   PJ=1
1208                ENDIF
1209                KTP(MAX(I,K),MIN(I,K))=
1210      &              (PI*KTUP(I,K)+PJ*KTUP(J,K))/(PI+PJ)
1211             ENDIF
1212  100     CONTINUE
1213       ENDIF
1214       IF (RECO.EQ.1) THEN
1215 C---VECTOR ADDITION
1216          P(1,I)=P(1,I)+P(1,J)
1217          P(2,I)=P(2,I)+P(2,J)
1218          P(3,I)=P(3,I)+P(3,J)
1219 c         P(4,I)=P(4,I)+P(4,J) ! JA
1220          P(5,I)=SQRT(P(1,I)**2+P(2,I)**2+P(3,I)**2)
1221          P(4,I)=P(5,I) ! JA (Massless scheme)
1222          IF (P(5,I).EQ.0) THEN
1223             P(5,I)=1
1224          ELSE
1225             P(5,I)=1/P(5,I)
1226          ENDIF
1227       ELSEIF (RECO.EQ.2) THEN
1228 C---PT WEIGHTED ETA-PHI ADDITION
1229          PT=P(6,I)+P(6,J)
1230          IF (PT.EQ.0) THEN
1231             PTT=1
1232          ELSE
1233             PTT=1/PT
1234          ENDIF
1235          P(7,I)=(P(6,I)*P(7,I)+P(6,J)*P(7,J))*PTT
1236          P(8,I)=KTMDPI(P(8,I)+P(6,J)*PTT*KTMDPI(P(8,J)-P(8,I)))
1237          P(6,I)=PT
1238          P(9,I)=PT**2
1239       ELSEIF (RECO.EQ.3) THEN
1240 C---PT**2 WEIGHTED ETA-PHI ADDITION
1241          PT=P(9,I)+P(9,J)
1242          IF (PT.EQ.0) THEN
1243             PTT=1
1244          ELSE
1245             PTT=1/PT
1246          ENDIF
1247          P(7,I)=(P(9,I)*P(7,I)+P(9,J)*P(7,J))*PTT
1248          P(8,I)=KTMDPI(P(8,I)+P(9,J)*PTT*KTMDPI(P(8,J)-P(8,I)))
1249          P(6,I)=P(6,I)+P(6,J)
1250          P(9,I)=P(6,I)**2
1251       ELSE
1252          CALL KTWARN('KTMERG',202,*999)
1253          STOP
1254       ENDIF
1255 C---IF MONO.GT.0 CALCULATE NEW KT MEASURES. IF MONO.GT.1 USE ANGULAR ONES.
1256       IF (MONO.LE.0) RETURN
1257 C---CONVERTING BETWEEN 4-MTM AND PT,ETA,PHI IF NECESSARY
1258       IF (ANGL.NE.1.AND.RECO.EQ.1) THEN
1259          P(9,I)=P(1,I)**2+P(2,I)**2
1260          P(7,I)=P(4,I)**2-P(3,I)**2
1261          IF (P(7,I).LE.EPS*P(4,I)**2) P(7,I)=P(9,I)
1262          IF (P(7,I).GT.0) THEN
1263             P(7,I)=LOG((P(4,I)+ABS(P(3,I)))**2/P(7,I))/2
1264             IF (P(7,I).GT.ETAMAX) P(7,I)=ETAMAX+2
1265          ELSE
1266             P(7,I)=ETAMAX+2
1267          ENDIF
1268          P(7,I)=SIGN(P(7,I),P(3,I))
1269          IF (P(1,I).NE.0.AND.P(2,I).NE.0) THEN
1270             P(8,I)=ATAN2(P(2,I),P(1,I))
1271          ELSE
1272             P(8,I)=0
1273          ENDIF
1274       ELSEIF (ANGL.EQ.1.AND.RECO.NE.1) THEN
1275          P(1,I)=P(6,I)*COS(P(8,I))
1276          P(2,I)=P(6,I)*SIN(P(8,I))
1277          P(3,I)=P(6,I)*SINH(P(7,I))
1278          P(4,I)=P(6,I)*COSH(P(7,I))
1279          IF (P(4,I).NE.0) THEN
1280             P(5,I)=1/P(4,I)
1281          ELSE
1282             P(5,I)=1
1283          ENDIF
1284       ENDIF
1285       ANG=0
1286       DO 200 K=1,N
1287          IF (K.NE.I.AND.K.NE.J) THEN
1288             IF (MONO.GT.1) ANG=KTUP(I,K)
1289             KTP(MIN(I,K),MAX(I,K))=
1290      &           KTPAIR(ANGL,P(1,I),P(1,K),ANG)
1291          ENDIF
1292  200  CONTINUE
1293       KTS(I)=KTSING(ANGL,TYPE,P(1,I))
1294  999  END
1295 C-----------------------------------------------------------------------
1296       SUBROUTINE KTMOVE(P,KTP,KTS,NMAX,N,J,IOPT)
1297       IMPLICIT NONE
1298 C---MOVE THE Nth PARTICLE IN P TO THE Jth POSITION
1299 C---ALSO MOVING KTP AND KTS IF IOPT.GT.0
1300       INTEGER I,J,N,NMAX,IOPT
1301       DOUBLE PRECISION P(9,NMAX),KTP(NMAX,NMAX),KTS(NMAX)
1302       DO 100 I=1,9
1303          P(I,J)=P(I,N)
1304  100  CONTINUE
1305       IF (IOPT.LE.0) RETURN
1306       DO 110 I=1,J-1
1307          KTP(I,J)=KTP(I,N)
1308          KTP(J,I)=KTP(N,I)
1309  110  CONTINUE
1310       DO 120 I=J+1,N-1
1311          KTP(J,I)=KTP(I,N)
1312          KTP(I,J)=KTP(N,I)
1313  120  CONTINUE
1314       KTS(J)=KTS(N)
1315       END
1316 C-----------------------------------------------------------------------
1317       SUBROUTINE KTUNIT(R)
1318       IMPLICIT NONE
1319 C   SET R EQUAL TO THE 4 BY 4 IDENTITY MATRIX
1320       DOUBLE PRECISION R(4,4)
1321       INTEGER I,J
1322       DO 20 I=1,4
1323         DO 10 J=1,4
1324           R(I,J)=0
1325           IF (I.EQ.J) R(I,J)=1
1326  10     CONTINUE
1327  20   CONTINUE
1328       END
1329 C-----------------------------------------------------------------------
1330       SUBROUTINE KTLBST(IOPT,R,A,*)
1331       IMPLICIT NONE
1332 C   PREMULTIPLY R BY THE 4 BY 4 MATRIX TO
1333 C   LORENTZ BOOST TO/FROM THE CM FRAME OF A
1334 C   IOPT=0 => TO
1335 C   IOPT=1 => FROM
1336 C
1337 C   LAST ARGUMENT IS LABEL TO JUMP TO IF A IS NOT TIME-LIKE
1338 C
1339       INTEGER IOPT,I,J
1340       DOUBLE PRECISION R(4,4),A(4),B(4),C(4,4),M
1341       DO 10 I=1,4
1342         B(I)=A(I)
1343  10   CONTINUE
1344       M=B(4)**2-B(1)**2-B(2)**2-B(3)**2
1345       IF (M.LE.0) CALL KTWARN('KTLBST',100,*999)
1346       M=SQRT(M)
1347       B(4)=B(4)+M
1348       M=1/(M*B(4))
1349       IF (IOPT.EQ.0) THEN
1350         B(4)=-B(4)
1351       ELSEIF (IOPT.NE.1) THEN
1352         CALL KTWARN('KTLBST',200,*999)
1353         STOP
1354       ENDIF
1355       DO 30 I=1,4
1356         DO 20 J=1,4
1357           C(I,J)=B(I)*B(J)*M
1358           IF (I.EQ.J) C(I,J)=C(I,J)+1
1359  20     CONTINUE
1360  30   CONTINUE
1361       C(4,4)=C(4,4)-2
1362       CALL KTMMUL(C,R,R)
1363       RETURN
1364  999  RETURN 1
1365       END
1366 C-----------------------------------------------------------------------
1367       SUBROUTINE KTRROT(R,A,B,*)
1368       IMPLICIT NONE
1369 C   PREMULTIPLY R BY THE 4 BY 4 MATRIX TO
1370 C   ROTATE FROM VECTOR A TO VECTOR B BY THE SHORTEST ROUTE
1371 C   IF THEY ARE EXACTLY BACK-TO-BACK, THE ROTATION AXIS IS THE VECTOR
1372 C   WHICH IS PERPENDICULAR TO THEM AND THE X AXIS, UNLESS THEY ARE
1373 C   PERPENDICULAR TO THE Y AXIS, WHEN IT IS THE VECTOR WHICH IS
1374 C   PERPENDICULAR TO THEM AND THE Y AXIS.
1375 C   NOTE THAT THESE CONDITIONS GUARANTEE THAT IF BOTH ARE PERPENDICULAR
1376 C   TO THE Z AXIS, IT WILL BE USED AS THE ROTATION AXIS.
1377 C
1378 C   LAST ARGUMENT IS LABEL TO JUMP TO IF EITHER HAS LENGTH ZERO
1379 C
1380       DOUBLE PRECISION R(4,4),M(4,4),A(4),B(4),C(4),D(4),AL,BL,CL,DL,EPS
1381 C---SQRT(2*EPS) IS THE ANGLE IN RADIANS OF THE SMALLEST ALLOWED ROTATION
1382 C   NOTE THAT IF YOU CONVERT THIS PROGRAM TO SINGLE PRECISION, YOU WILL
1383 C   NEED TO INCREASE EPS TO AROUND 0.5E-4
1384       PARAMETER (EPS=0.5D-6)
1385       AL=A(1)**2+A(2)**2+A(3)**2
1386       BL=B(1)**2+B(2)**2+B(3)**2
1387       IF (AL.LE.0.OR.BL.LE.0) CALL KTWARN('KTRROT',100,*999)
1388       AL=1/SQRT(AL)
1389       BL=1/SQRT(BL)
1390       CL=(A(1)*B(1)+A(2)*B(2)+A(3)*B(3))*AL*BL
1391 C---IF THEY ARE COLLINEAR, DON'T NEED TO DO ANYTHING
1392       IF (CL.GE.1-EPS) THEN
1393         RETURN
1394 C---IF THEY ARE BACK-TO-BACK, USE THE AXIS PERP TO THEM AND X AXIS
1395       ELSEIF (CL.LE.-1+EPS) THEN
1396         IF (ABS(B(2)).GT.EPS) THEN
1397           C(1)= 0
1398           C(2)=-B(3)
1399           C(3)= B(2)
1400 C---UNLESS THEY ARE PERPENDICULAR TO THE Y AXIS,
1401         ELSE
1402           C(1)= B(3)
1403           C(2)= 0
1404           C(3)=-B(1)
1405         ENDIF
1406 C---OTHERWISE FIND ROTATION AXIS
1407       ELSE
1408         C(1)=A(2)*B(3)-A(3)*B(2)
1409         C(2)=A(3)*B(1)-A(1)*B(3)
1410         C(3)=A(1)*B(2)-A(2)*B(1)
1411       ENDIF
1412       CL=C(1)**2+C(2)**2+C(3)**2
1413       IF (CL.LE.0) CALL KTWARN('KTRROT',101,*999)
1414       CL=1/SQRT(CL)
1415 C---FIND ROTATION TO INTERMEDIATE AXES FROM A
1416       D(1)=A(2)*C(3)-A(3)*C(2)
1417       D(2)=A(3)*C(1)-A(1)*C(3)
1418       D(3)=A(1)*C(2)-A(2)*C(1)
1419       DL=AL*CL
1420       M(1,1)=A(1)*AL
1421       M(1,2)=A(2)*AL
1422       M(1,3)=A(3)*AL
1423       M(1,4)=0
1424       M(2,1)=C(1)*CL
1425       M(2,2)=C(2)*CL
1426       M(2,3)=C(3)*CL
1427       M(2,4)=0
1428       M(3,1)=D(1)*DL
1429       M(3,2)=D(2)*DL
1430       M(3,3)=D(3)*DL
1431       M(3,4)=0
1432       M(4,1)=0
1433       M(4,2)=0
1434       M(4,3)=0
1435       M(4,4)=1
1436       CALL KTMMUL(M,R,R)
1437 C---AND ROTATION FROM INTERMEDIATE AXES TO B
1438       D(1)=B(2)*C(3)-B(3)*C(2)
1439       D(2)=B(3)*C(1)-B(1)*C(3)
1440       D(3)=B(1)*C(2)-B(2)*C(1)
1441       DL=BL*CL
1442       M(1,1)=B(1)*BL
1443       M(2,1)=B(2)*BL
1444       M(3,1)=B(3)*BL
1445       M(1,2)=C(1)*CL
1446       M(2,2)=C(2)*CL
1447       M(3,2)=C(3)*CL
1448       M(1,3)=D(1)*DL
1449       M(2,3)=D(2)*DL
1450       M(3,3)=D(3)*DL
1451       CALL KTMMUL(M,R,R)
1452       RETURN
1453  999  RETURN 1
1454       END
1455 C-----------------------------------------------------------------------
1456       SUBROUTINE KTVMUL(M,A,B)
1457       IMPLICIT NONE
1458 C   4 BY 4 MATRIX TIMES 4 VECTOR: B=M*A.
1459 C   ALL ARE DOUBLE PRECISION
1460 C   IT IS SAFE TO CALL WITH B=A
1461 C   FIRST SUBSCRIPT=ROWS, SECOND=COLUMNS
1462       DOUBLE PRECISION M(4,4),A(4),B(4),C(4)
1463       INTEGER I,J
1464       DO 20 I=1,4
1465         C(I)=0
1466         DO 10 J=1,4
1467           C(I)=C(I)+M(I,J)*A(J)
1468  10     CONTINUE
1469  20   CONTINUE
1470       DO 30 I=1,4
1471         B(I)=C(I)
1472  30   CONTINUE
1473       END
1474 C-----------------------------------------------------------------------
1475       SUBROUTINE KTMMUL(A,B,C)
1476       IMPLICIT NONE
1477 C   4 BY 4 MATRIX MULTIPLICATION: C=A*B.
1478 C   ALL ARE DOUBLE PRECISION
1479 C   IT IS SAFE TO CALL WITH C=A OR B.
1480 C   FIRST SUBSCRIPT=ROWS, SECOND=COLUMNS
1481       DOUBLE PRECISION A(4,4),B(4,4),C(4,4),D(4,4)
1482       INTEGER I,J,K
1483       DO 30 I=1,4
1484         DO 20 J=1,4
1485           D(I,J)=0
1486           DO 10 K=1,4
1487             D(I,J)=D(I,J)+A(I,K)*B(K,J)
1488  10       CONTINUE
1489  20     CONTINUE
1490  30   CONTINUE
1491       DO 50 I=1,4
1492         DO 40 J=1,4
1493           C(I,J)=D(I,J)
1494  40     CONTINUE
1495  50   CONTINUE
1496       END
1497 C-----------------------------------------------------------------------
1498       SUBROUTINE KTINVT(A,B)
1499       IMPLICIT NONE
1500 C---INVERT TRANSFORMATION MATRIX A
1501 C
1502 C   A = INPUT  : 4 BY 4 TRANSFORMATION MATRIX
1503 C   B = OUTPUT : INVERTED TRANSFORMATION MATRIX
1504 C
1505 C   IF A IS NOT A TRANSFORMATION MATRIX YOU WILL GET STRANGE RESULTS
1506 C
1507 C   NOTE THAT IT IS SAFE TO CALL WITH A=B
1508 C
1509       DOUBLE PRECISION A(4,4),B(4,4),C(4,4)
1510       INTEGER I,J
1511 C---TRANSPOSE
1512       DO 20 I=1,4
1513         DO 10 J=1,4
1514           C(I,J)=A(J,I)
1515  10     CONTINUE
1516  20   CONTINUE
1517 C---NEGATE ENERGY-MOMENTUM MIXING TERMS
1518       DO 30 I=1,3
1519         C(4,I)=-C(4,I)
1520         C(I,4)=-C(I,4)
1521  30   CONTINUE
1522 C---OUTPUT
1523       DO 50 I=1,4
1524         DO 40 J=1,4
1525           B(I,J)=C(I,J)
1526  40     CONTINUE
1527  50   CONTINUE
1528       END
1529 C-----------------------------------------------------------------------
1530       FUNCTION KTMDPI(PHI)
1531       IMPLICIT NONE
1532 C---RETURNS PHI, MOVED ONTO THE RANGE [-PI,PI)
1533       DOUBLE PRECISION KTMDPI,PHI,PI,TWOPI,THRPI,EPS
1534       PARAMETER (PI=3.14159265358979324D0,TWOPI=6.28318530717958648D0,
1535      &     THRPI=9.42477796076937972D0)
1536       PARAMETER (EPS=1D-15)
1537       KTMDPI=PHI
1538       IF (KTMDPI.LE.PI) THEN
1539         IF (KTMDPI.GT.-PI) THEN
1540           GOTO 100
1541         ELSEIF (KTMDPI.GT.-THRPI) THEN
1542           KTMDPI=KTMDPI+TWOPI
1543         ELSE
1544           KTMDPI=-MOD(PI-KTMDPI,TWOPI)+PI
1545         ENDIF
1546       ELSEIF (KTMDPI.LE.THRPI) THEN
1547         KTMDPI=KTMDPI-TWOPI
1548       ELSE
1549         KTMDPI=MOD(PI+KTMDPI,TWOPI)-PI
1550       ENDIF
1551  100  IF (ABS(KTMDPI).LT.EPS) KTMDPI=0
1552       END
1553 C-----------------------------------------------------------------------
1554       SUBROUTINE KTWARN(SUBRTN,ICODE,*)
1555 C     DEALS WITH ERRORS DURING EXECUTION
1556 C     SUBRTN = NAME OF CALLING SUBROUTINE
1557 C     ICODE  = ERROR CODE:    - 99 PRINT WARNING & CONTINUE
1558 C                          100-199 PRINT WARNING & JUMP
1559 C                          200-    PRINT WARNING & STOP DEAD
1560 C-----------------------------------------------------------------------
1561       INTEGER ICODE
1562       CHARACTER*6 SUBRTN
1563       WRITE (6,10) SUBRTN,ICODE
1564    10 FORMAT(/' KTWARN CALLED FROM SUBPROGRAM ',A6,': CODE =',I4/)
1565       IF (ICODE.LT.100) RETURN
1566       IF (ICODE.LT.200) RETURN 1
1567       STOP
1568       END
1569 C-----------------------------------------------------------------------
1570 C-----------------------------------------------------------------------
1571 C-----------------------------------------------------------------------