syjkrn.f

Go to the documentation of this file.

CDECK  ID>, SYJKRN.
**************************************************************************
      SUBROUTINE syjkrn(MD,NJETRQ,THRESH,IMODE,
     +  ycuts,iexam,frac,iexam4,npar,idim,ppar,mxjet,
     +  njet,iass,jdim,pjets,y34,y45,ierr)
* Kernel Routine for SY-Jetfindings Scheme or Traditional Jetfindings
* Authour   : Satoru Yamashita
* Creat     : 5-SEP-98
*           :
* Inputs:
*   MD      ; mode (see below)
*   NJETRQ  : fixed number of jets ; =0 variable -- use YCUT etc...
*   THRESH  : threshold (GeV) for the primary jet finding
*   IMODE   : Kernel Jet finding mode
*      1- 6 : YKERN
*       IF(IMODE.EQ.1) THEN
*          CM = 'JADE E0'
*        ELSEIF(IMODE.EQ.2) THEN
*          CM = 'JADE P '
*        ELSEIF(IMODE.EQ.3) THEN
*          CM = 'JADE P0'
*        ELSEIF(IMODE.EQ.4) THEN
*          CM = 'JADE E '
*        ELSEIF(IMODE.EQ.5) THEN
*          CM = 'DURHAM '
*        ELSEIF(IMODE.EQ.6) THEN
*          CM = 'GENEVA '
*        ELSE
*          WRITE(6,281) IMODE
* 281      FORMAT(/,' ### YKERN: IMODE =',I3,' INVALID; SET TO 1 ###')
*          IMODE = 1
*          CM = 'JADE E0'
*        ENDIF
*
*
*
*         7 : cambridge
*      9-11 : LUCLUS
*        12 : cone
*   YCUTS   : (1) ycut or xcut or Cone Radius (2) EPS
*   IEXAM   : Particle re-association method 1-6
*            EXAM(1) = THETA
*            EXAM(2) = PART(4)*PJSTMP(4,J)-PART(1)*PJSTMP(1,J)
*     +                -PART(2)*PJSTMP(2,J)-PART(3)*PJSTMP(3,J)
*            EXAM(3) = DURHAM
*            EXAM(4) = E0JADE
*            EXAM(5) = Geneva Scheme
*            EXAM(6) = EJADE   ! it makes too bad cross talk
*
*   FRAC    : core rate
*   IEXAM4  : Jet merge method 1-8
*            EXAM(1,N) = PL(4,I)*PL(4,J)
*            EXAM(2,N) = PL(4,I)*PL(4,J)-PL(1,I)*PL(1,J)
*     +           -PL(2,I)*PL(2,J)-PL(3,I)*PL(3,J)
*            EXAM(3,N) = DURHAM(I,J)
*            EXAM(4,N) = E0JADE(I,J)
*            EXAM(5,N) = EJADE(I,J)
*            EXAM(6,N) = ANGLE(I,J)
*            EXAM(7,N) = ECOS(I,J)
*            EXAM(8,N) = PL(6,I)*PL(6,J)
*   NPAR    : Number of "particles"
*   IDIM    : first dimension of PPAR array
*   PPAR    : Particle Momenta
*   MXJET   : Maximum Number of Jets (buffer size for PJETS = MXJETS*JDIM)
*   JDIM    : first dimension of PJETS array
*
* Outputs:
*   NJET    : Number of Jets
*   IASS    : Jet assosiation for "particle"s
*   PJETS   : Jet momenta
*   Y34     : YHI (sometimes not filled)
*   Y45     : YLO (sometimes not filled)
*   IERR    : Error flag  0==All O.K.
*     -9999: abnormal error
*      -999: boundary error
*       -99: input mismatch
*      Others : error from various internal calls
* ==========================================================================
*   MD
* There are 11 variations
* MD   Name      : Final NJ    : 1st process  : 2nd process : 3rd process
* 0a. Traditional: variable    :  fixed Ycuts :     -       :   -
* 0b. Traditional: fixed Njets :  fixed Njets :     -       :   -
* 0c. Just Merge : fixed Njets :  fixed Ycuts : do merge    :   -
* 1a. LMODE1     : variable    :  fixed Ycuts : do SYJETF1  :   -
* 1b. LMODE1     : fixed final :  fixed Ycuts : do SYJETF1  : do merge
* 1c. LMODE1     : fixed final :  fixed Njets : do SYJETF1  :   -
* 1d. LMODE1     : fixed final :  fixed Ycuts : do merge    : do SYJETF1
* 2a. LMODE2     : variable    :  fixed Ycuts : do SYJETF2  :   -
* 2b. LMODE2     : fixed final :  fixed Ycuts : do SYJETF2  : do merge
* 2c. LMODE2     : fixed final :  fixed Njets : do SYJETF2  :   -
* 2d. LMODE2     : fixed final :  fixed Ycuts : do merge    : do SYJETF2
* ============================================================================
*  MD & NJETRQ &  IMODE  & YCUTS & IEXAM &  Frac   &  IEXAM4 &     ISYMOD
*  0a &   -    &   1-12  &  > 0  &   -   &   -     &     -   &        1
*  0b &  > 0   &   1-12  &   -   &   -   &   -     &     -   &        2
*  0c &  > 0   &   1-12  &  > 0  &   -   &   -     &    1-8  &        3
*  1a &   -    &   1-12  &  > 0  &  1-6  & 0.0-1.0 &     -   &       11
*  1b &  > 0   &   1-12  &  > 0  &  1-6  & 0.0-1.0 &    1-8  &       12
*  1c &  > 0   &   1-12  &   -   &  1-6  & 0.0-1.0 &     -   &       13
*  1d &  > 0   &   1-12  &  > 0  &  1-6  & 0.0-1.0 &    1-8  &       14
*  2a &   -    &   1-12  &  > 0  &  1-6  & 0.0-1.0 &     -   &       21
*  2b &  > 0   &   1-12  &  > 0  &  1-6  & 0.0-1.0 &    1-8  &       22
*  2c &  > 0   &   1-12  &   -   &  1-6  & 0.0-1.0 &     -   &       23
*  2d &  > 0   &   1-12  &  > 0  &  1-6  & 0.0-1.0 &    1-8  &       24
* ============================================================================
****************************************************************************
      implicit none
      CHARACTER*(*) md
      INTEGER idim,jdim,mxjet
      INTEGER npar
      REAL    ppar(idim,*)
      REAL    pjets(jdim,*)
      INTEGER njetrq,imode,iexam,iexam4,njet,iass(*),ierr
      REAL    thresh,ycuts(*),frac,y34,y45

      INTEGER mxwrk
      parameter(mxwrk = 300)
      INTEGER npar0
      REAL    ppar0(5,mxwrk)
      INTEGER    iuse(mxwrk)
      REAL    ycut,r0,eps0

      INTEGER i,j
      INTEGER ij(mxwrk),ijn(mxwrk),ijj(mxwrk)
      REAL    pj6(6,50),pjn6(6,50),ptmp6(6,50)
      REAL    pord(50),rmax
      INTEGER iiii(50),ii
      REAL    ydhi,ydlo
      INTEGER imode0
      INTEGER isymod

*-----Initialysation for outputs
      ierr = -9999
      njet = 0
      DO i=1,npar
         iass(i) = -1
      ENDDO
      CALL vzero(pjets,mxjet*jdim)
      y34 = 0.0
      y45 = 0.0
      isymod = 0

*-----Check boundary
      IF(md(1:2).EQ.'0a'.OR.md(1:2).EQ.'0A') isymod = 1
      IF(md(1:2).EQ.'0b'.OR.md(1:2).EQ.'0B') isymod = 2
      IF(md(1:2).EQ.'0c'.OR.md(1:2).EQ.'0C') isymod = 3
      IF(md(1:2).EQ.'1a'.OR.md(1:2).EQ.'1A') isymod = 11
      IF(md(1:2).EQ.'1b'.OR.md(1:2).EQ.'1B') isymod = 12
      IF(md(1:2).EQ.'1c'.OR.md(1:2).EQ.'1C') isymod = 13
      IF(md(1:2).EQ.'1d'.OR.md(1:2).EQ.'1D') isymod = 14
      IF(md(1:2).EQ.'2a'.OR.md(1:2).EQ.'2A') isymod = 21
      IF(md(1:2).EQ.'2b'.OR.md(1:2).EQ.'2B') isymod = 22
      IF(md(1:2).EQ.'2c'.OR.md(1:2).EQ.'2C') isymod = 23
      IF(md(1:2).EQ.'2d'.OR.md(1:2).EQ.'2D') isymod = 24
*      write(*,*)'symode selected: ',ISYMOD
*
      IF(isymod.EQ.0) THEN
         WRITE(*,'(A,A)')'unknown MODE detected:'//md(1:2)
*         CALL REPORT('SYJKRN',1,'E')
         ierr  = -999
         RETURN
      ENDIF
      IF(idim.LT.4) THEN
         WRITE(*,'(A,I5)')'too small IDIM',idim
*         CALL REPORT('SYJKRN',2,'E')
         ierr  = -999
         RETURN
      ENDIF
      IF(npar.GT.mxwrk) THEN
         WRITE(*,'(A,I5)')'too many particles',npar
*         CALL REPORT('SYJKRN',3,'E')
         ierr  = -999
         RETURN
      ENDIF
      IF(imode.LE.0.OR.imode.GT.12) THEN
         WRITE(*,'(A,I5)')'unknown Jet find mode',imode
*         CALL REPORT('SYJKRN',4,'E')
         ierr = -999
         RETURN
      ENDIF

*-----check inputs
      ierr = -99
      IF(isymod.EQ.1) THEN
         IF(ycuts(1).LE.0.0) THEN
            WRITE(*,'(A,I6,E12.3)')'ycuts small?',isymod,ycuts(1)
*            CALL REPORT('SYJKRN',5,'E')
            RETURN
         ENDIF
      ELSEIF(isymod.EQ.2) THEN
         IF(njetrq.LE.0) THEN
            WRITE(*,'(A,I6,I7)')'No NjetRQ',isymod,njetrq
*            CALL REPORT('SYJKRN',6,'E')
            RETURN
         ENDIF
      ELSEIF(isymod.EQ.3) THEN
         IF(njetrq.LE.0) THEN
            WRITE(*,'(A,I6,I7)')'No NjetRQ',isymod,njetrq
*            CALL REPORT('SYJKRN',7,'E')
            RETURN
         ENDIF
         IF(ycuts(1).LE.0.0) THEN
            WRITE(*,'(A,I6,E12.3)')'ycuts small?',isymod,ycuts(1)
*            CALL REPORT('SYJKRN',8,'E')
            RETURN
         ENDIF
         IF(iexam4.LE.0.OR.iexam4.GT.8) THEN
            WRITE(*,'(A,I6,I7)')'No merge way',isymod,iexam4
*            CALL REPORT('SYJKRN',9,'E')
            RETURN
         ENDIF
      ELSEIF(isymod.EQ.11.OR.isymod.EQ.21) THEN
         IF(ycuts(1).LE.0.0) THEN
            WRITE(*,'(A,I6,E12.3)')'ycuts small?',isymod,ycuts(1)
*            CALL REPORT('SYJKRN',10,'E')
            RETURN
         ENDIF
         IF(iexam.LE.0.OR.iexam.GT.6) THEN
            WRITE(*,'(A,I6,I7)')'No re-assign way',isymod,iexam
*            CALL REPORT('SYJKRN',11,'E')
            RETURN
         ENDIF
      ELSEIF(isymod.EQ.12.OR.isymod.EQ.22) THEN
         IF(njetrq.LE.0) THEN
            WRITE(*,'(A,I6,I7)')'No NjetRQ',isymod,njetrq
*            CALL REPORT('SYJKRN',12,'E')
            RETURN
         ENDIF
         IF(ycuts(1).LE.0.0) THEN
            WRITE(*,'(A,I6,E12.3)')'ycuts small?',isymod,ycuts(1)
*            CALL REPORT('SYJKRN',13,'E')
            RETURN
         ENDIF
         IF(iexam4.LE.0.OR.iexam4.GT.8) THEN
            WRITE(*,'(A,I6,I7)')'No merge way',isymod,iexam4
*            CALL REPORT('SYJKRN',14,'E')
            RETURN
         ENDIF
         IF(iexam.LE.0.OR.iexam.GT.6) THEN
            WRITE(*,'(A,I6,I7)')'No re-assign way',isymod,iexam
*            CALL REPORT('SYJKRN',15,'E')
            RETURN
         ENDIF
      ELSEIF(isymod.EQ.13.OR.isymod.EQ.23) THEN
         IF(njetrq.LE.0) THEN
            WRITE(*,'(A,I6,I7)')'No NjetRQ',isymod,njetrq
*            CALL REPORT('SYJKRN',16,'E')
            RETURN
         ENDIF
         IF(iexam.LE.0.OR.iexam.GT.6) THEN
            WRITE(*,'(A,I6,I7)')'No re-assign way',isymod,iexam
*            CALL REPORT('SYJKRN',17,'E')
            RETURN
         ENDIF
      ELSEIF(isymod.EQ.14.OR.isymod.EQ.24) THEN
         IF(njetrq.LE.0) THEN
            WRITE(*,'(A,I6,I7)')'No NjetRQ',isymod,njetrq
*            CALL REPORT('SYJKRN',18,'E')
            RETURN
         ENDIF
         IF(ycuts(1).LE.0.0) THEN
            WRITE(*,'(A,I6,E12.3)')'ycuts small?',isymod,ycuts(1)
*            CALL REPORT('SYJKRN',19,'E')
            RETURN
         ENDIF
         IF(iexam4.LE.0.OR.iexam4.GT.8) THEN
            WRITE(*,'(A,I6,I7)')'No merge way',isymod,iexam4
*            CALL REPORT('SYJKRN',20,'E')
            RETURN
         ENDIF
         IF(iexam.LE.0.OR.iexam.GT.6) THEN
            WRITE(*,'(A,I6,I7)')'No re-assign way',isymod,iexam
*            CALL REPORT('SYJKRN',21,'E')
            RETURN
         ENDIF
      ENDIF

*-----Preparation
      npar0 = npar
      DO i=1,npar0
         CALL ucopy(ppar(1,i),ppar0(1,i),min(5,idim))
         IF(idim.LT.5) THEN
            ppar0(5,i)=ppar0(4,i)**2-
     +           ppar0(1,i)**2-ppar0(2,i)**2-ppar0(3,i)**2
            ppar0(5,i)=sqrt(max(0.0,ppar0(5,i)))
         ENDIF
      ENDDO
*      IF(MOD(IMODE0,100).EQ.12) THEN
      IF (imode.EQ.12) THEN    ! bugfix JS
         ycut = 0.0
         r0   = ycuts(1)
         eps0 = ycuts(2)
      ELSE
         ycut = ycuts(1)
         r0   = 0.0
         eps0 = 0.0
      ENDIF

*-----Check threshold
      DO i=1,npar0
         IF(ppar0(4,i).GT.thresh) THEN
            iuse(i) = 1
         ELSE
            iuse(i) = 0
         ENDIF
      ENDDO

*###### S T A R T #######
*-----Traditional Jet findings variable Njet
      IF(isymod.EQ.1) THEN
         imode0 = imode
         CALL syjetf0(imode0,ycut,eps0,r0,0,npar0,ppar0,iuse,
     +        njet,ij,ptmp6,ydhi,ydlo,ierr)
         IF(ierr.NE.0) goto 999
         CALL syjcp6(npar0,ppar0,njet,ij,6,pj6)
*-----Traditional Jet findings fixed Njet
      ELSEIF(isymod.EQ.2) THEN
         imode0 = imode+100
         CALL syjetf0(imode0,ycut,eps0,r0,njetrq,npar0,ppar0,iuse,
     +        njet,ij,ptmp6,y34,y45,ierr)
         IF(ierr.NE.0) goto 999
         CALL syjcp6(npar0,ppar0,njet,ij,6,pj6)
*-----Traditional Jet findings variable Njet and then merge if Njet>Njetrq
      ELSEIF(isymod.EQ.3) THEN
         imode0 = imode
         CALL syjetf0(imode0,ycut,eps0,r0,0,npar0,ppar0,iuse,
     +        njet,ijn,ptmp6,ydhi,ydlo,ierr)
         IF(ierr.NE.0) njet = 0
         IF(njet.LT.njetrq) THEN
            imode0 = imode + 100
            CALL syjetf0(imode0,ycut,eps0,r0,njetrq,npar0,ppar0,iuse,
     +           njet,ijn,ptmp6,y34,y45,ierr)
            IF(ierr.NE.0) goto 999
         ENDIF
         CALL syjcp6(npar0,ppar0,njet,ijn,6,pjn6)
         CALL symrgj(njetrq,iexam4,njet,npar0,ijn,pjn6,ij,pj6,ierr)  ! merge
         IF(ierr.NE.0) goto 999
         njet = njetrq
*-----Re-assignment Jet findings variable Njet ; no merge
      ELSEIF(isymod.EQ.11.OR.isymod.EQ.21) THEN
         imode0 = imode
         CALL syjetf0(imode0,ycut,eps0,r0,0,npar0,ppar0,iuse,
     +        njet,ijn,ptmp6,ydhi,ydlo,ierr)
         IF(ierr.NE.0) goto 999
         IF(isymod.LT.20) THEN
*          3 times iteration
            CALL syjetf1(iexam,njet,npar0,ppar0,ijn,frac,ij ,pj6)
            CALL syjetf1(iexam,njet,npar0,ppar0,ij ,frac,ijn,pj6)
            CALL syjetf1(iexam,njet,npar0,ppar0,ijn,frac,ij ,pj6)
         ELSE  ! 21
*          3 times iteration
            CALL syjetf2(iexam,njet,npar0,ppar0,ijn,frac,ij ,pj6)
            CALL syjetf2(iexam,njet,npar0,ppar0,ij ,frac,ijn,pj6)
            CALL syjetf2(iexam,njet,npar0,ppar0,ijn,frac,ij ,pj6)
         ENDIF
*-----Re-assignment Jet findings variable Njet ;  merge if njet>njetrq
      ELSEIF(isymod.EQ.12.OR.isymod.EQ.22) THEN
         imode0 = imode
         CALL syjetf0(imode0,ycut,eps0,r0,0,npar0,ppar0,iuse,
     +        njet,ijn,ptmp6,ydhi,ydlo,ierr)
         IF(ierr.NE.0) njet = 0
         IF(njet.LT.njetrq) THEN
            imode0 = imode + 100
            CALL syjetf0(imode0,ycut,eps0,r0,njetrq,npar0,ppar0,iuse,
     +           njet,ijn,ptmp6,y34,y45,ierr)
            IF(ierr.NE.0) goto 999
         ENDIF
         IF(isymod.LT.20) THEN
*          3 times iteration
            CALL syjetf1(iexam,njet,npar0,ppar0,ijn,frac,ijj,pjn6)
            CALL syjetf1(iexam,njet,npar0,ppar0,ijj,frac,ijn,pjn6)
            CALL syjetf1(iexam,njet,npar0,ppar0,ijn,frac,ijj,pjn6)
         ELSE  ! 21
*          3 times iteration
            CALL syjetf2(iexam,njet,npar0,ppar0,ijn,frac,ijj,pjn6)
            CALL syjetf2(iexam,njet,npar0,ppar0,ijj,frac,ijn,pjn6)
            CALL syjetf2(iexam,njet,npar0,ppar0,ijn,frac,ijj,pjn6)
         ENDIF
*        Merge
         CALL symrgj(njetrq,iexam4,njet,npar0,ijj,pjn6,ij,pj6,ierr)
         IF(ierr.NE.0) goto 999
         njet = njetrq
*-----Re-assignment Jet findings variable Njet ; no merge
      ELSEIF(isymod.EQ.13.OR.isymod.EQ.23) THEN
         imode0 = imode + 100
         CALL syjetf0(imode0,ycut,eps0,r0,njetrq,npar0,ppar0,iuse,
     +        njet,ijn,ptmp6,y34,y45,ierr)
         IF(ierr.NE.0) goto 999
         IF(isymod.LT.20) THEN
*          3 times iteration
            CALL syjetf1(iexam,njet,npar0,ppar0,ijn,frac,ij ,pj6)
            CALL syjetf1(iexam,njet,npar0,ppar0,ij ,frac,ijn,pj6)
            CALL syjetf1(iexam,njet,npar0,ppar0,ijn,frac,ij ,pj6)
         ELSE  ! 21
*          3 times iteration
            CALL syjetf2(iexam,njet,npar0,ppar0,ijn,frac,ij ,pj6)
            CALL syjetf2(iexam,njet,npar0,ppar0,ij ,frac,ijn,pj6)
            CALL syjetf2(iexam,njet,npar0,ppar0,ijn,frac,ij ,pj6)
         ENDIF
*-----Re-assignment Jet findings variable Njet ;  merge if njet>njetrq
      ELSEIF(isymod.EQ.14.OR.isymod.EQ.24) THEN
         imode0 = imode
         CALL syjetf0(imode0,ycut,eps0,r0,0,npar0,ppar0,iuse,
     +        njet,ijn,ptmp6,ydhi,ydlo,ierr)
         IF(ierr.NE.0) njet = 0
         IF(njet.LT.njetrq) THEN
            imode0 = imode + 100
            CALL syjetf0(imode0,ycut,eps0,r0,njetrq,npar0,ppar0,iuse,
     +           njet,ijn,ptmp6,y34,y45,ierr)
            IF(ierr.NE.0) goto 999
         ENDIF
         CALL syjcp6(npar0,ppar0,njet,ijn,6,pjn6)
         CALL symrgj(njetrq,iexam4,njet,npar0,ijn,pjn6,ijj,pjn6,ierr)
         IF(ierr.NE.0) goto 999
         njet = njetrq
         IF(isymod.LT.20) THEN
*          3 times iteration
            CALL syjetf1(iexam,njet,npar0,ppar0,ijj,frac,ij ,pj6)
            CALL syjetf1(iexam,njet,npar0,ppar0,ij ,frac,ijj,pj6)
            CALL syjetf1(iexam,njet,npar0,ppar0,ijj,frac,ij ,pj6)
         ELSE  ! 21
*          3 times iteration
            CALL syjetf2(iexam,njet,npar0,ppar0,ijj,frac,ij ,pj6)
            CALL syjetf2(iexam,njet,npar0,ppar0,ij ,frac,ijj,pj6)
            CALL syjetf2(iexam,njet,npar0,ppar0,ijj,frac,ij ,pj6)
         ENDIF
      ELSE
         WRITE(*,'(A,I10)')'unknown ISYMOD',isymod
*         CALL REPORT('SYJKRN',22,'E')
         ierr = -9999
         goto 999
      ENDIF
*=========== Energy Ordering ================
      DO i=1,npar0
         ijj(i)=ij(i)
      ENDDO
      DO i=1,njet
         pord(i) = pj6(4,i)
      ENDDO
      DO i=1,njet
         rmax = -1.e+09
         ii=0
         DO j=1,njet
            IF(pord(j).GT.rmax) THEN
               rmax = pord(j)
               ii = j
            ENDIF
         ENDDO
         iiii(i)  = ii
         pord(ii) = -1.e+09
      ENDDO
      DO i=1,npar0
         iass(i) = 0
         DO j=1,njet
            IF(iiii(j).EQ.ijj(i)) THEN
               iass(i) = j
            ENDIF
         ENDDO
      ENDDO
      DO i=1,min(njet,mxjet)
         DO j=1,min(6,jdim)
            pjets(j,i)=pj6(j,iiii(i))
         ENDDO
      ENDDO

      ierr = 0
*---- ALL DONE : GOOD
*      WRITE(*,'(A,I10,I10)')'complete',NJET,ISYMOD
*      CALL REPORT('SYJKRN',23,'I')

      RETURN

 999  CONTINUE
*----- error, so all initialyze again.
      WRITE(*,'(A,I10,I10,I6)')'Error Found',ierr,njet,isymod
*      CALL REPORT('SYJKRN',24,'E')
      njet = 0
      DO i=1,npar
         iass(i)=0
      ENDDO
      CALL vzero(pjets,mxjet*jdim)
      RETURN
      END

*********************************************************************
      SUBROUTINE syjcp6(NPAR,PPAR,NJ,IASS,JDIM,PJ6)
* Just Re-calculate Jet Momenta
*  Author ; Satoru Yamashita
*  Creat  ; 10-May-98
*********************************************************************
      implicit none
      INTEGER npar,iass(*),jdim,nj
      REAL    ppar(5,*),pj6(jdim,*),rr,qq
      INTEGER i
      CALL vzero(pj6,nj*jdim)
      DO i=1,npar
         IF(iass(i).GT.0.AND.iass(i).LE.nj) THEN
            pj6(1,iass(i)) = pj6(1,iass(i)) + ppar(1,i)
            pj6(2,iass(i)) = pj6(2,iass(i)) + ppar(2,i)
            pj6(3,iass(i)) = pj6(3,iass(i)) + ppar(3,i)
            pj6(4,iass(i)) = pj6(4,iass(i)) + ppar(4,i)
         ENDIF
      ENDDO
      IF(jdim.GE.5) THEN
         DO i=1,nj
            rr = pj6(1,i)**2+pj6(2,i)**2+pj6(3,i)**2
            qq = pj6(4,i)**2-rr
            pj6(5,i)=sqrt(max(qq,0.))
            IF(jdim.GE.6) pj6(6,i) = sqrt(rr)
         ENDDO
      ENDIF
      END