      SUBROUTINE HDPM
 
C-----------------------------------------------------------------------
C  H(ADRONIC) D(UAL) P(ARTON) M(ODEL)
C
C  GENERATOR OF HADRONIC COLLISION INSPIRED BY DUAL PARTON MODEL
C  THIS SUBROUTINE IS CALLED FROM SDPM
C-----------------------------------------------------------------------
 
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
*KEEP,CONST.
      COMMON /CONST/   PI,PI2,OB3,TB3,ENEPER
      DOUBLE PRECISION PI,PI2,OB3,TB3,ENEPER
*KEEP,DPMFLG.
      COMMON /DPMFLG/  NFLAIN,NFLDIF,NFLPI0,NFLCHE,NFLPIF,NFRAGM
      INTEGER          NFLAIN,NFLDIF,NFLPI0,NFLCHE,NFLPIF,NFRAGM
*KEEP,ELADPM.
      COMMON /ELADPM/  ELMEAN,ELMEAA,IELDPM,IELDPA
      DOUBLE PRECISION ELMEAN(37),ELMEAA(37)
      INTEGER          IELDPM(37,13),IELDPA(37,13)
*KEEP,ELASTY.
      COMMON /ELASTY/  ELAST,IELIS,IELHM,IELNU,IELPI
      DOUBLE PRECISION ELAST
      INTEGER          IELIS(20),IELHM(20),IELNU(20),IELPI(20)
*KEEP,INDICE.
      COMMON /INDICE/  NNUCN,NKA0,NHYPN,NETA,NETAS,NPIZER,
     *                 NNC,NKC,NHC,NPC,NCH,NNN,NKN,NHN,NET,NPN
      INTEGER          NNUCN(2:3),NKA0(2:3),NHYPN(2:3),NETA(2:3,1:4),
     *                 NETAS(2:3),NPIZER(2:3),
     *                 NNC,NKC,NHC,NPC,NCH,NNN,NKN,NHN,NET,NPN
*KEEP,INTER.
      COMMON /INTER/   AVCH,AVCH3,DC0,DLOG,DMLOG,ECMDIF,ECMDPM,ELAB,
     *                 FNEUT,FNEUT2,GNU,PLAB,POSC2,POSC3,POSN2,POSN3,
     *                 RC3TO2,S,SEUGF,SEUGP,SLOG,SLOGSQ,SMLOG,
     *                 WIDC2,WIDC3,WIDN2,WIDN3,YCM,YY0,ZN,
     *                 IDIF,ITAR
      DOUBLE PRECISION AVCH,AVCH3,DC0,DLOG,DMLOG,ECMDIF,ECMDPM,ELAB,
     *                 FNEUT,FNEUT2,GNU,PLAB,POSC2,POSC3,POSN2,POSN3,
     *                 RC3TO2,S,SEUGF,SEUGP,SLOG,SLOGSQ,SMLOG,
     *                 WIDC2,WIDC3,WIDN2,WIDN3,YCM,YY0,ZN
      INTEGER          IDIF,ITAR
*KEEP,ISTA.
      COMMON /ISTA/    IFINET,IFINNU,IFINKA,IFINPI,IFINHY
      INTEGER          IFINET,IFINNU,IFINKA,IFINPI,IFINHY
*KEEP,LEPAR.
      COMMON /LEPAR/   LEPAR1,LEPAR2,LASTPI,NRESPC,NRESPN,NCPLUS
      INTEGER          LEPAR1,LEPAR2,LASTPI,NRESPC,NRESPN,NCPLUS
*KEEP,MULT.
      COMMON /MULT/    EKINL,MSMM,MULTMA,MULTOT
      DOUBLE PRECISION EKINL
      INTEGER          MSMM,MULTMA(37,13),MULTOT(37,13)
*KEEP,NEWPAR.
      COMMON /NEWPAR/  EA,PT2,PX,PY,TMAS,YR,ITYP,
     *                 IA1,IA2,IB1,IB2,IC1,IC2,ID1,ID2,IE1,IE2,IF1,IF2,
     *                 IG1,IG2,IH1,IH2,II1,II2,IJ1,NTOT
      DOUBLE PRECISION EA(3000),PT2(3000),PX(3000),PY(3000),TMAS(3000),
     *                 YR(3000)
      INTEGER          ITYP(3000),
     *                 IA1,IA2,IB1,IB2,IC1,IC2,ID1,ID2,IE1,IE2,IF1,IF2,
     *                 IG1,IG2,IH1,IH2,II1,II2,IJ1,NTOT
*KEEP,PAM.
      COMMON /PAM/     PAMA,SIGNUM
      DOUBLE PRECISION PAMA(6000),SIGNUM(6000)
*KEEP,PARPAR.
      COMMON /PARPAR/  CURPAR,SECPAR,PRMPAR,OUTPAR,C,
     *                 E00,E00PN,PTOT0,PTOT0N,THICKH,ITYPE,LEVL
      DOUBLE PRECISION CURPAR(14),SECPAR(14),PRMPAR(14),OUTPAR(14),
     *                 C(50),E00,E00PN,PTOT0,PTOT0N,THICKH
      INTEGER          ITYPE,LEVL
*KEEP,PARPAE.
      DOUBLE PRECISION GAMMA,COSTHE,PHI,H,T,X,Y,CHI,BETA,GCM,ECM
      EQUIVALENCE      (CURPAR(2),GAMMA),  (CURPAR(3),COSTHE),
     *                 (CURPAR(4), PHI ),  (CURPAR(5), H    ),
     *                 (CURPAR(6), T   ),  (CURPAR(7), X    ),
     *                 (CURPAR(8), Y   ),  (CURPAR(9), CHI  ),
     *                 (CURPAR(10),BETA),  (CURPAR(11),GCM  ),
     *                 (CURPAR(12),ECM )
*KEEP,RANDPA.
      COMMON /RANDPA/  FAC,U1,U2,RD,NSEQ,ISEED,KNOR
      DOUBLE PRECISION FAC,U1,U2
      REAL             RD(3000)
      INTEGER          ISEED(103,10),NSEQ
      LOGICAL          KNOR
*KEEP,RATIOS.
      COMMON /RATIOS/  RPI0R,RPIER,RPEKR,RPEKNR,PPICH,PPINCH,PPNKCH,
     *                 ISEL,NEUTOT,NTOTEM
      DOUBLE PRECISION RPI0R,RPIER,RPEKR,RPEKNR,PPICH,PPINCH,PPNKCH
      INTEGER          ISEL,NEUTOT,NTOTEM
*KEEP,RESON.
      COMMON /RESON/   RDRES,RESRAN,IRESPAR
      REAL             RDRES(2),RESRAN(1000)
      INTEGER          IRESPAR
 
*KEEP,REST.
      COMMON /REST/    CONTNE,TAR,LT
      DOUBLE PRECISION CONTNE(3),TAR
      INTEGER          LT
*KEEP,RUNPAR.
      COMMON /RUNPAR/  FIXHEI,THICK0,HILOECM,HILOELB,
     *                 STEPFC,NRRUN,NSHOW,PATAPE,MONIIN,
     *                 MONIOU,MDEBUG,NUCNUC,
     *                 CETAPE,
     *                 SHOWNO,ISHW,NOPART,NRECS,NBLKS,MAXPRT,NDEBDL,
     *                 N1STTR,MDBASE,
     *                 DEBDEL,DEBUG,FDECAY,FEGS,FIRSTI,FIXINC,FIXTAR,
     *                 FIX1I,FMUADD,FNKG,FPRINT,FDBASE
     *                ,GHEISH,GHESIG
      COMMON /RUNPAC/  DSN,HOST,USER
      DOUBLE PRECISION FIXHEI,THICK0,HILOECM,HILOELB
      REAL             STEPFC
      INTEGER          NRRUN,NSHOW,PATAPE,MONIIN,MONIOU,MDEBUG,NUCNUC,
     *                 SHOWNO,ISHW,NOPART,NRECS,NBLKS,MAXPRT,NDEBDL,
     *                 N1STTR,MDBASE
      INTEGER          CETAPE
      CHARACTER*79     DSN
      CHARACTER*20     HOST,USER
 
      LOGICAL          DEBDEL,DEBUG,FDECAY,FEGS,FIRSTI,FIXINC,FIXTAR,
     *                 FIX1I,FMUADD,FNKG,FPRINT,FDBASE
     *                ,GHEISH,GHESIG
*KEND.
 
C-----------------------------------------------------------------------
 
      IF ( DEBUG ) WRITE(MDEBUG,444) (CURPAR(I),I=1,9)
  444   FORMAT(' HDPM  : CURPAR=',1P,9E10.3)
 
C  SET ANTI-LEADER TO PROTON OR NEUTRON; TARGET IS ALWAYS NUCLEON
      CALL RMMAR( RD,1,1 )
      IF ( RD(1) .LE. CONTNE(LT) ) THEN
        ITAR = 13
      ELSE
        ITAR = 14
      ENDIF
 
 
C  CALCULATE LAB AND CM ENERGY
      IF ( ITYPE .NE. 1 ) THEN
        ELAB  = PAMA(ITYPE) * GAMMA
        PLAB  = ELAB * BETA
        S     = PAMA(ITYPE)**2 + PAMA(ITAR)**2 + 2.D0*PAMA(ITAR)*ELAB
      ELSE
C  FOR GAMMA-INDUCED REACTION TAKE PI(0) AS LEADING PARTICLE
        ITYPE = 7
        ELAB  = GAMMA
        PLAB  = ELAB
        S     = PAMA(ITAR)**2 + 2.D0*PAMA(ITAR)*ELAB
      ENDIF
 
      ECMDPM  = SQRT(S)
      IF ( DEBUG ) WRITE(MDEBUG,*) 'HDPM  : ITYPE,ELAB,PLAB,S,ECMDPM=',
     *             ITYPE,SNGL(ELAB),SNGL(PLAB),SNGL(S),SNGL(ECMDPM)
 
C  LN(S), LN(S)**2 AND RAPIDITY OF C. M. SYSTEM IN LAB
      SLOG    = LOG(S)
      SLOGSQ  = SLOG**2
      SMLOG   = LOG( 2.D0 * PAMA(ITAR) * ELAB )
      ELABLG  = LOG(ELAB)
      EPLUSP  = ELAB + PLAB
*     YCM = 0.5D0 * LOG( (ELAB+PAMA(ITAR)+PLAB)/(ELAB+PAMA(ITAR)-PLAB) )
      YCM = 0.5D0 * LOG( (EPLUSP**2     +PAMA(ITAR)*EPLUSP)/
     *                   (PAMA(ITYPE)**2+PAMA(ITAR)*EPLUSP) )
      IF ( DEBUG ) WRITE(MDEBUG,*)'HDPM  : SLOG,SLOGSQ,YCM=',
     *                    SNGL(SLOG),SNGL(SLOGSQ),SNGL(YCM)
 
C-----------------------------------------------------------------------
C  RETURN POINT IF CALCULATION OF PARTICLES GOES WRONG
   1  CONTINUE
 
      IF ( ITYPE .NE. 7 ) THEN
C  CHOOSE NUMBER OF INTERACTIONS IN TARGET
        CALL TARINT
      ELSE
C  FOR GAMMA-INDUCED REACTIONS TAKE ALWAYS ONE COLLISION
        GNU = 1.D0
      ENDIF
 
C-----------------------------------------------------------------------
C  NO DIFFRACTION IF
C  OR       THE NUMBER OF INTERACTIONS IN TARGET IS CHOSEN RANDOMLY
C     AND   MORE THAN ONE INTERACTION TAKES PLACE
C  OR       PRIMARY PARTICLE IS GAMMA (PI0)
C  NOW NFLDIF DECIDES WHETHER DIFFRACTIVE PROCESS POSSIBLE OR NOT
      IF ( ( NFLAIN.EQ.0 .AND. GNU.GT.1.D0 .AND. NFLDIF.EQ.0 )
     *      .OR. ( ITYPE .EQ. 7 ) ) THEN
        IDIF = 0
      ELSE
C  SET DIFFRACTION FLAG IF RANDOM NUMBER < PROBABILITY
        CALL RMMAR( RD,1,1 )
C  IDIF IS 0 : NO DIFFRACTION ; IDIF IS 1 : DIFFRACTION
C  DIFFRACTION RISES WITH ENERGY AND SATURATES AT 10000 GEV
C  ### DAS TUT ES ABER NICHT: ES IST KONSTANT 0.15 (SIEHE DPFUNC) !!!!
        IF ( RD(1) .GT. DPFUNC(ECMDPM) ) THEN
          IDIF = 0
        ELSE
          IDIF = 1
        ENDIF
      ENDIF
 
 
C  PRINTOUT FOR DEBUG
      IF ( DEBUG ) THEN
        WRITE(MDEBUG,*) '    DIFFRACTIVE INTERACTION (0/1) = ',IDIF
      ENDIF
 
C  SET COUNTER FOR REPEAT TO 0
      NREPRD = 0
 
C  GENERATION OF INTERACTION
 1919 CONTINUE
 
C  FLAG TO CHECK NUMBER OF SECONDARIES;
C  IS CHANGED TO 1 IF SECONDARY MULTIPLICITY IS LOW
      ISEL   = 0
C  SET LEADING PARTICLE TO INCOMING PARTICLE AND ANTI-LEADER TO NUCLEON
C  (AS IT COMES FROM TARGET NUCLEUS) BOTH MAY BE CHANGED BY LEPACX
      LEPAR1 = ITYPE
      LEPAR2 = ITAR
 
      IF ( IDIF .EQ. 0 ) THEN
C-----------------------------------------------------------------------
C  NON SINGLE DIFFRACTIVE PROCESS STARTS HERE
 
        CALL NSD
C  CHARGE EXCHANGE ENABLED? EXCHANGE LEADER AND ANTI-LEADER
        LASTPI = 0
        NRESPC = 0
        NRESPN = 0
        NCPLUS = 0
        IF ( NFLCHE .EQ. 0 ) THEN
          CALL LEPACX( ECMDPM,ELABLG,LEPAR1,1 )
          CALL LEPACX( ECMDPM,ELABLG,LEPAR2,2 )
        ENDIF
 1921   CONTINUE
        CALL RNEGBI( NCH,AVCH,ECMDPM )
C  NCH IS # OF ALL CHARGED PARTICLES INCLUDING EXCESS FROM TARGET
        IF ( NCH .LT. 1 ) THEN
          IF ( LEPAR1 .LT. 50  .OR.  LEPAR2 .LT. 50 ) THEN
            NREPRD = NREPRD + 1
            IF ( NREPRD .GT. 10 ) GOTO 1
            GOTO 1921
          ELSE
C  INTERACTION IS ONLY RESONANCE PRODUCTION
            ISEL = 1
          ENDIF
        ENDIF
C  WIDTH PLATEAU FOR CLUSTERS AND FOR CALCULATION OF CENTR.RAP.DENSITY
        DELRAP = 0.6722D0 * (2.95D0 + 0.0302D0 * SLOG)
C  SET RSLOG FOR CALCULATION OF PARTICLE RATIOS
        RSLOG  = SLOG
C  AVERAGE TRANSVERSE MOMENTUM
        CALL AVEPT( ECMDPM,SLOG )
 
      ELSE
C-----------------------------------------------------------------------
C  SINGLE DIFFRACTIVE PROCESS STARTS HERE
 
 1920   CONTINUE
        CALL DIFRAC( NRETDF )
        IF ( NRETDF .EQ. 1 ) GOTO 1
C  CHARGE EXCHANGE ENABLED?  EXCHANGE CHARGE OF DIFFRACTING PARTICLE
        LASTPI = 0
        NRESPC = 0
        NRESPN = 0
        NCPLUS = 0
        IF ( NFLCHE .EQ. 0 ) THEN
          IF ( YY0 .GT. 0.D0 ) THEN
C  PROJECTILE DIFFRACTION
            CALL LEPACX( ECMDIF,DMLOG,LEPAR1,1 )
          ELSE
C  TARGET DIFFRACTION
            CALL LEPACX( ECMDIF,DMLOG,LEPAR2,2 )
          ENDIF
        ENDIF
C  FLUCTUATION OF MULTIPLICITY ACCORDING TO NEG.BIN. DISTRIBUTION
        CALL RNEGBI( NCH,AVCH,ECMDIF )
C  REPEAT CALCULATION AS SOMETHING WENT WRONG
        IF ( NCH .LT. 1 ) THEN
          IF ( (YY0 .GT. 0.D0  .AND.  LEPAR1 .LT. 50)  .OR.
     *         (YY0 .LT. 0.D0  .AND.  LEPAR2 .LT. 50)        ) THEN
            NREPRD = NREPRD + 1
            IF ( NREPRD .GT. 10 ) GOTO 1
            GOTO 1920
          ELSE
C  DIFFRACTIVE INTERACTION IS ONLY RESONANCE PRODUCTION
            ISEL = 1
          ENDIF
        ENDIF
C  SET RSLOG FOR CALCULATION OF PARTICLE RATIOS
        RSLOG = DLOG
C  HERE WE USE ECMDPM, BECAUSE THE MOMENTUM TRANSFER IS DEPENDENT
C  ON THE ENERGY OF THE TOTAL SYSTEM AND NOT ON THE DIFFRACTING MASS
        CALL AVEPT( ECMDPM,SLOG )
 
      ENDIF
 
C-----------------------------------------------------------------------
C  NOW FOR NSD AND DIFFRACTIVE PROCESSES
 
C  IN CASE OF LOW MULTIPLICITY SET FLAG ISEL
      IF ( NCH .LE. 2 ) ISEL=1
C  FNCH IS FLUCTUATING TOT.NUMBER OF CHARGED PARTICLES FOR ALL 3 STRINGS
      FNCH  = DBLE(NCH)
C  RATIO  ALL CHARGED PARTICLES  WITH FLUCTUATION/WITHOUT FLUCTUATION
      XZ    = FNCH / AVCH
C  FNCH3 IS FLUCTUATING NUMBER OF CHARGED PARTICLES FOR 3RD STRING
      FNCH3 = XZ * AVCH3
C  FNCH2 IS FLUCTUATING NUMBER OF CHARGED PARTICLES 1ST AND 2ND STRING
      FNCH2 = FNCH - FNCH3
C  RC3TO2 IS RATIO (CHARGED 3RD STRING)/(CHARGED 1ST AND 2ND STRING)
      IF ( FNCH2 .NE. 0.D0 ) THEN
        RC3TO2 = FNCH3 / FNCH2
      ELSE
        RC3TO2 = 0.D0
      ENDIF
      IF ( DEBUG ) WRITE(MDEBUG,*) '   FNCH,FNCH2,FNCH3,RC3TO2=',
     *          SNGL(FNCH),SNGL(FNCH2),SNGL(FNCH3),SNGL(RC3TO2)
 
C  IS NUMBER OF NEUTRALS FLUCTUATING AS NUMBER OF CHARGED ?
      IF ( NFLPIF .EQ. 0 .OR. IDIF .EQ. 1 .OR. ECMDPM .LT. 60.D0 ) THEN
C  SET NUMBER OF GAMMAS ACCORDING TO NEG. BIN. VARIABLE XZ
C  AS NUMBER OF NEUTRALS FLUCTUATES AS CHARGED.
        SEUGF = SEUGP * XZ
        ZG    = XZ
      ELSE
C  NFLPIF IS 1 MEANS: # OF PI(0) FLUCTUATES AS MEASURED AT COLLIDER
        IF ( ECMDPM .LT. 200.D0 ) THEN
          SEUGF = SEUGP * XZ
*         SEUGF = (0.0786D0*SLOG-0.010D0)*FNCH2 + (0.391D0*SLOG+0.305D0)
        ELSE
C  DETERMINE NEW NUMBER OF GAMMAS WITH FLUCTUATION AROUND SEUGP*XZ
          AGR   = EXP(-XZ)
          DGR   = SEUGP * XZ * (0.9823D0 - 0.3756D0 * AGR)
          SGS   = DGR * (0.14718D0 + 2.53213D0 * AGR)
 723      CONTINUE
          SEUGF = 0.88D0 * RANNOR(DGR,SGS)
          IF ( SEUGF .LT. 1.D0 ) GOTO 723
        ENDIF
C  SET NEGATIVE BINOMIAL VARIABLE ZG FOR GAMMAS
        ZG  = SEUGF / SEUGP
      ENDIF
      SEUGF = MAX( 1.D0, SEUGF )
      IF ( DEBUG ) WRITE(MDEBUG,*) 'HDPM  :XZ,ZG,SEUGF=',
     *      SNGL(XZ),SNGL(ZG),SNGL(SEUGF)
 
C-----------------------------------------------------------------------
C  RATIO ALL-NUCLEON/ALL-CHARGED
C  PARAMETRISATION FROM UA5, NUCL. PHYS. B291 (1987) 445 EQ.(2.4)
      RNUCCH = MAX( 0.D0, -0.008D0 + 0.00865D0 * RSLOG )
C  NUMBER FOR DIRECT NEUTRON/ANTINEUTRON PRODUCTION 1ST AND 2ND STRING
C  MULTIPLY BY 0.5 BECAUSE RATIO RNUCCH GIVES (ALL-NUCL)/(ALL-CHARGED)
C  AND HERE ONLY THE NEUTRON-ANTINEUTRONS ARE COUNTED
      FNUCN  = 0.5D0 * RNUCCH * FNCH2
C  RATIO (ALL CHARGED SIGMAS)/(ALL CHARGED) IS 1/3 OF ALL STRANGE BARYON
C  PARAMETRISATION FORM UA5, NUCL. PHYS. B291 (1987) 445 EQ.(2.5)
      RHYPCH = MAX( 0.D0, (-0.007D0 + 0.0028D0 * RSLOG) * OB3 )
C  NEUTRAL STRANGE BARYONS ARE DOUBLE OF CHARGED STRANGE BARYONS
      FHYPN  = 2.D0 * RHYPCH * FNCH2
C  CORRECT NUMBER OF GAMMAS FROM NEUTRAL HYPERON DECAY S0-->L+GAMMA
      SEUGFC = MAX( 0.D0, SEUGF - 0.5D0 * FHYPN )
C  RATIO CHARGED-KAON/CHARGED PIONS
C  PARAMETRISATION FROM UA5, NUCL. PHYS. B291 (1987) 445 EQ.(2.7)
      RKPI   = MAX (0.D0, 0.024D0 + 0.0062D0 * RSLOG )
C  RKCH IS RATIO (CHARGED-KAON)/(ALL-CHARGED) DERIVED FROM RKPI;
C  THE FACTOR 0.5 IN FRONT OF RNUCCH IS BECAUSE ONLY HALF OF NUCLEONS
C  ARE P/PBAR. THE 1.17 IS AN APROXIMATE CONVERSION FACTOR FROM
C  (ALL-NUCL)/(ALL-CHARGED) TO (ALL-NUCL)/(CHARGED-PI), WHICH IS A BIT
C  ENERGY DEPENDENT (1.14 ...1.21) SEE GEICH-GIMBEL TABLE 7.1
      RKCH   = RKPI / (1.D0 + RKPI + (0.5D0*RNUCCH+RHYPCH) * 1.17D0)
C  K0/K0-BAR  FOR 1ST AND 2ND STRING
C  NEUTRAL KAONS ARE PRODUCED WITH THE SAME RATE AS CHARGED KAONS
      FKA0   = RKCH * FNCH2
C  RATIO ETA/PI(0) IS ASSUMED TO BE INDEPENDENT OF ENERGY = 0.19
C  SEE: ANSORGE ET AL. (UA5-COLLABORATION) Z.PHYS.C43(1989)75
*     RETPI0 = 0.19D0
C  RATIO ETA/PI(0) IS ASSUMED TO BE DEPENDENT ON ENERGY
C  SEE: GEICH-GIMBEL,INT.J.MOD.PHYS.A4(1989)1527 TAB.7.1
C  HECK'S FIT:   RETPI0 IS 0.06 + 0.006*RSLOG + 0.0011*RSLOG**2
      RETPI0 = 0.06D0 + 0.006D0 * RSLOG + 0.0011D0 * RSLOG**2
C  AUXIL1 IS FRACTION OF PI(0)/(PI(0)+ETA)
      AUXIL1 = 1.D0 / (1.D0 + RETPI0)
C  NUMBER OF GAMMAS FROM PI(0) IS 2, FROM ETA IS 3.216 IN AVERAGE;
C  AUXIL2 IS NUMBER OF GAMMA-PRODUCING-PARTICLES: PI(0) AND ETA
      AUXIL2 = SEUGFC / ( AUXIL1 * 2.D0 + (1.D0 - AUXIL1) * 3.216D0 )
      FETA   = (1.D0 - AUXIL1) * AUXIL2
      FPI0   =         AUXIL1  * AUXIL2
C  CORRECT FPI0 BY DECAYS OF STRANGE BARYONS; NEUTRAL: FHYPN*0.357
C  CHARGED: 0.5*FNCH2*RHYPCH*0.5157; IT YIELDS FHYPN*(0.357+0.12893)
      FPI0   = MAX( 0.D0, FPI0 - FHYPN * 0.486D0 )
C  SUMMED NEUTRAL PARTICLES FOR 1ST AND 2ND STRING
      FNEUT2 = FNUCN + FKA0 + FHYPN + FETA + FPI0
C  NEUTRAL PARTICLES FROM 3RD STRING
      FNEUT3 = RC3TO2 * FNEUT2
C  TOTAL NUMBER OF NEUTRALS
      FNEUT  = FNEUT2 + FNEUT3
      NEUTOT = NINT( FNEUT )
C  CALCULATE TOTAL NUMBER OF PARTICLES TO BE CREATED
      NTOTEM = NCH + NEUTOT
      IF ( DEBUG ) WRITE(MDEBUG,*)
     *     '    FNUCN,FKA0,FHYPN,FETA,FPI0,FNEUT2,FNEUT3,NTOTEM=',
     * SNGL(FNUCN),SNGL(FKA0),SNGL(FHYPN),SNGL(FETA),SNGL(FPI0),
     * SNGL(FNEUT2),SNGL(FNEUT3),NTOTEM
C  LIMIT OF SECONDARIES PRODUCED (GIVEN BY SIZE OF ARRAY : 3000)
C  LIMIT IS ARRAY SIZE - SIZE OF LARGEST TARGET NUCLEUS(40)
      IF ( NTOTEM .GE. 2956 ) THEN
        WRITE(MONIOU,*) 'HDPM  : REJECT EVENT WITH ',NTOTEM,
     *                  ' SECONDARIES'
        GOTO 1
      ENDIF
C  SPECIAL TREATMENT IF MULTIPLICITY IS TOO LOW
      IF ( NTOTEM .LE. 3 ) ISEL = 1
 
C  FRACTION OF THE VARIOUS NEUTRAL PARTICLES (NN, K(0), L+S0 AS PAIRS)
C  NORMALIZE WITH THE SUM OF ALL NEUTRAL PARTICLES
      FNORML = 1.D0 / ( 0.5D0 * (FNUCN+FKA0+FHYPN) + FETA + FPI0 )
      RNUCNR = FNUCN * FNORML * 0.5D0
      RKA0R  = FKA0  * FNORML * 0.5D0
      RHYPNR = FHYPN * FNORML * 0.5D0
      RETAR  = FETA  * FNORML
      RPI0R  = FPI0  * FNORML
C  CUMULATED RATIOS (NN, K(0), LAMBDA+SIGMA0 AS PAIRS)
      RPIER  = RPI0R + RETAR
      RPEKR  = RPIER + RKA0R
      RPEKNR = RPEKR + RNUCNR
C THEN THE REMAINDER (1-RPEKNR) MUST BE NEUTRAL HYPERON PAIRS
      IF ( DEBUG ) WRITE(MDEBUG,*)
     *         '    RPI0R,RETAR,RKA0R,RNUCNR,RHYPNR,FNORML=',
     *  SNGL(RPI0R),SNGL(RETAR),SNGL(RKA0R),SNGL(RNUCNR),SNGL(RHYPNR),
     *  SNGL(FNORML)
 
C  PROBABILITY TO PRODUCE CHARGED PIONS IS PROBABILITY NOT TO PRODUCE
C  CHARGED KAONS OR PROTONS OR CHARGED HYPERONS, WHERE PROTON/ANTIPROTON
C  IS HALF OF (ALL-NUCL)/(ALL-CHARGED)
      AUXIL  = RKCH + 0.5D0 * RNUCCH + RHYPCH
      AUXIL3 = 1.D0 - AUXIL
C  RENORMALIZATION AS P/P_BAR, K+-, AND HYPERONS ARE PRODUCED IN PAIRS
C  AUXIL2 IS INVERSE OF NORMALISATION
      AUXIL2 = 1.D0 / (1.D0 - 0.5D0 * AUXIL)
C  CUMULATED PROBABILITIES (PP, K+-, SIGMA+- AS PAIRS)
      PPICH  = AUXIL3 * AUXIL2
      PPINCH =  PPICH + 0.25D0 * RNUCCH * AUXIL2
      PPNKCH =  PPINCH + 0.5D0 * RKCH * AUXIL2
C THEN THE REMAINDER (1-PPNKCH) MUST BE CHARGED HYPERON PAIRS
      IF ( DEBUG ) WRITE(MDEBUG,*) '   PPICH,PPINCH,PPNKCH=',
     *         SNGL(PPICH),SNGL(PPINCH),SNGL(PPNKCH)
 
C  NOW SELECT HOW MANY PARTICLES OF EACH TYPE ARE PRODUCED
      CALL PARNUM( INUMFL )
      IF ( INUMFL .NE. 0 ) GOTO 1919
 
C  DEFINE PARTICLE NUMBERS WHERE SPECIAL RAPIDITY IS CALCULATED
C  FOR PARTICLES FROM TARGET (THIRD STRING)
      PPP = RC3TO2 / (1.D0+RC3TO2)
C  NUMBER OF PARTICLES IN PROTON ANTIPROTON PAIRS FROM TARGET
      ITA = NINT(PPP * 2.D0 * NNC)
C  NUMBER OF PARTICLES IN K+ K- PAIRS FROM TARGET
      ITB = NINT(PPP * 2.D0 * NKC)
C  NUMBER OF PARTICLES IN SIGMA+ SIGMA- PAIRS FROM TARGET
      ITC = NINT(PPP * 2.D0 * NHC)
C  NUMBER OF PI+ PI- FROM TARGET
      ITD = NINT(PPP * NPC )
C  CALCULATE BOUNDARIES
      IA1 = 2
      IA2 = IA1 + ITA
      IB1 = IA1 + 2 * NNC
      IB2 = IB1 + ITB
      IC1 = IB1 + 2 * NKC
      IC2 = IC1 + ITC
      ID1 = IC1 + 2 * NHC
      ID2 = ID1 + ITD
      IE1 = ID1 +     NPC
C  NUMBER OF PARTICLES IN NEUTRON ANTINEUTRON PAIRS FROM TARGET
      IE2 = IE1 + 2 * NNUCN(3)
      IF1 = IE1 + 2 * NNN
C  NUMBER OF PARTICLES IN K0S K0L PAIRS FROM TARGET
      IF2 = IF1 + 2 * NKA0(3)
      IG1 = IF1 + 2 * NKN
C  NUMBER OF PARTICLES IN NEUTRAL HYPERON PAIRS FROM TARGET
      IG2 = IG1 + 2 * NHYPN(3)
      IH1 = IG1 + 2 * NHN
C  NUMBER OF ETA FROM TARGET
      IH2 = IH1 + NETAS(3)
      II1 = IH1 +     NET
C  NUMBER OF PI(0) FROM TARGET
      II2 = II1 + NPIZER(3)
      IJ1 = II1 +     NPN
      IF ( DEBUG ) THEN
        WRITE(MDEBUG,*) '   CHARGED FROM TARGET:',ITA,ITB,ITC,ITD
        WRITE(MDEBUG,*) '   NEUTRAL FROM TARGET:',
     *     2*NNUCN(3),2*NKA0(3),2*NHYPN(3),NETAS(3),NPIZER(3)
        WRITE(MDEBUG,*) '   NTOTEM,IJ1=',NTOTEM,IJ1
      ENDIF
C  REDEFINE TOTAL NUMBER OF SECONDARY PARTICLES : NTOTEM
C  BY CHARGE EXCHANGE AND RESONANCE FORMATION THIS NUMBER MAY BE ALTERED
      NTOTEM = IJ1 - 2
 
C-----------------------------------------------------------------------
C  RATIO OF RAPIDITY DENSITY TO MEAN PSEUDORAPIDITY IN CENTER
C  PARAMETRISATION SEE CAPDEVIELLE, J.PHYS.G:NUCL.PHYS.15(1989)909,EQ.6
      IF ( XZ .LT. 1.5D0 ) THEN
        RDS = (0.24396D0 + 0.70150424D0 * XZ)**2
      ELSE
        RDS = (0.55685D0 + 0.48664753D0 * XZ)**2
      ENDIF
C  CALCULATE NOW: DN/DY AT Y = 0; DC0 IS AVERAGE PSEUDORAPIDITY DENSITY
C  TRAP IS RATIO (RAPID.DENS.)/(PSEUDORAP.DENS.) IN CENTER OF RAPIDITY
      TRAP = 1.25D0
      IF ( IDIF .EQ. 0  .AND.  ECMDPM .GT. 19.4D0 )
     *                 TRAP = MAX( 1.D0, 1.28852D0 - 0.0065D0 * SMLOG )
      DCN2 = DC0 * RDS * TRAP
      IF ( DEBUG ) WRITE(MDEBUG,*) '   RDS,TRAP,DCN2=',
     *                SNGL(RDS),SNGL(TRAP),SNGL(DCN2)
C  AMPLITUDE OF GAUSSIAN 1ST AND 2ND STRING
      ATG2 = FNCH2 / (5.0132566D0 * WIDC2)
C  NEW DEFINITION OF POSITION BASED ON SEMI INCLUSIVE DATA
      SQ2  = 2.D0 * ATG2 / DCN2
C  FINAL POSITION OF GAUSSIAN;  WIDTH WIDC2 IS UNCHANGED
      IF ( SQ2 .GT. 1.D0 )  POSC2 = WIDC2 * SQRT( 2.D0*LOG(SQ2) )
C  DENSITY OF CHARGED IN EXCESS FROM TARGET IN CENTER OF RAPIDITY
      DCN3 = 0.5D0 * (GNU - 1.D0) * DCN2
      IF (DEBUG) WRITE(MDEBUG,*) '   SQ2,POSC2,DCN3=',
     *          SNGL(SQ2),SNGL(POSC2),SNGL(DCN3)
      IF ( DCN3 .GT. 0.D0 ) THEN
C  AMPLITUDE 3RD GAUSSIAN
        ATG3 = FNCH3 / (5.0132566D0 * WIDC3)
C  AMPLITUDE IS DIVIDED BY DENSITY FOR GETTING CENTER OF 3RD GAUSSIAN
        SQ3  = 2.D0 * ATG3 / DCN3
C  CHECK IF ADDITIVE MULTIPLICITY IS TOO LOW
        IF ( SQ3 .GT. 1.D0 )  POSC3 = WIDC3 * SQRT( 2.D0*LOG(SQ3) )
        IF (DEBUG) WRITE(MDEBUG,*)'   SQ3,POSC3=',SNGL(SQ3),SNGL(POSC3)
      ENDIF
 
C  NFLPI0 .EQ. 0 MEANS TREAT PI(0) RAPIDITY ACCORDING TO COLLIDER DATA
      IF ( NFLPI0 .EQ. 0 ) THEN
C  RATIO OF RAPIDITY DENSITY TO MEAN PSEUDORAPIDITY AT CENTER WITH Z<1.5
        IF ( ZG .LT. 1.5D0 ) THEN
          RDG = (0.24396D0 + 0.70150424D0 * ZG)**2
        ELSE
          RDG = (0.55685D0 + 0.48664753D0 * ZG)**2
        ENDIF
C  GAMMAS USE RATIO TRAG TO CALCULATE RATIO OF RAPIDITY TO
C  PSEUDO RAPIDITY DENSITY IN CENTER (TRAG = 1.1 * 0.5 ).
C  FACTOR 0.5 COMES FROM RATIO NEUTRAL/CHARGED, AS WE USE DC0, WHICH
C  IS AVERAGE PSEUDORAPIDITY DENSITY FOR CHARGED PIONS
        TRAG = 0.55D0
        IF ( IDIF .EQ. 0 ) THEN
          IF ( ECMDPM .GT. 19.4D0 )
     *               TRAG = MAX( 0.4D0, 0.6658D0 - 0.01954D0 * SMLOG )
          IF     ( ECMDPM .LE.  50.D0 ) THEN
            DCG = DC0 * RDG * TRAG
          ELSEIF ( ECMDPM .LE. 200.D0 ) THEN
            DCG = DC0 * RDG * TRAG * (1.D0 + 0.18D0 * LOG(ECMDPM/50.D0))
          ELSE
            DCG = DC0 * RDG * TRAG * 1.25D0
          ENDIF
        ELSE
          DCG = DC0 * RDG * TRAG
        ENDIF
C  DEFINE WIDTH OF STRINGS FOR NEUTRAL PIONS AND ETAS
        WIDN2 = WIDC2 * MIN( 1.D0, 1.12275D0 - 0.0208D0 * RSLOG )
C  NEW DEFINITION OF CENTER OF GAUSSIAN BASED ON SEMI INCLUSIVE DATA
C  USING AMPLITUDE OF THE GAUSSIAN FOR NEUTRALS
        AUXIL = 2.D0 / (5.0132566D0 * WIDN2 * DCG)
C  TOTAL MULTIPLICITY USED FOR 1ST AND 2ND STRING OF PI(0) AND ETA
C  IS GIVEN BY THEIR NUMBERS. ANALOGOUS FOR 3RD STRING
        SP2   = DBLE ( NPIZER(2)+NETAS(2)) * AUXIL
C  FINAL CENTER OF GAUSSIANS FOR PI(0) AND ETA (WIDC2 IS UNCHANGED)
        IF ( SP2 .GT. 1.D0 ) THEN
          POSN2 = WIDN2 * SQRT( 2.D0 * LOG(SP2) )
        ELSE
          POSN2 = POSC2
        ENDIF
        WIDN3 = WIDN2
        SP3   = DBLE(NPIZER(3)+NETAS(3)) * AUXIL
        IF ( SP3 .GT. 1.D0 ) THEN
          POSN3 = WIDN3 * SQRT( 2.D0 * LOG(SP3) )
        ELSE
          POSN3 = POSC3
        ENDIF
      ELSE
C  NFLPI0 .EQ. 1 MEANS RAPIDITY OF PI(0) AND ETA SAME AS THAT OF CHARGED
        POSN2 = POSC2
        WIDN2 = WIDC2
        POSN3 = POSC3
        WIDN3 = WIDC3
      ENDIF
      IF ( DEBUG ) WRITE(MDEBUG,*)
     *           '   ZG,RDG,DCG,SP2,SP3,POSN2,POSN3,WIDN2 =',
     *   SNGL(ZG),SNGL(RDG),SNGL(DCG),SNGL(SP2),SNGL(SP3),SNGL(POSN2),
     *   SNGL(POSN3),SNGL(WIDN2)
 
C-----------------------------------------------------------------------
      NREPR1 = 0
C  RETURN POINT. NUMBERS OF PARTICLES REMAIN UNCHANGED FOR NEXT TRY,
C  BUT INDIVIDUAL RAPIDITIES GET NEW VALUES.
C  START FROM BEGINNING IF NO MATCH AFTER 20 TRIES
  30  CONTINUE
      NREPR1 = NREPR1 + 1
      IF ( NREPR1 .GT. 20 ) THEN
        IF ( IDIF .EQ. 1  .AND.  NREPRD .LE. 10 ) GOTO 1919
        GOTO 1
      ENDIF
 
C  FOR TOTAL NUMBER OF PARTICLES ADD 2 FOR LEADER AND ANTILEADER
      NTOT = NTOTEM + 2
 
C  PRODUCTION OF INDIVIDUAL RAPIDITIES FOR ALL SECONDARY PARTICLES
      CALL PARRAP
CC    IF ( DEBUG ) THEN
CC      WRITE (MDEBUG,*) '   RAPIDITIES:'
CC      WRITE (MDEBUG,134) (I,YR(I), I=3,NTOT)
C134    FORMAT('   ',1P, (1X, I4, 5X, G13.6 ))
CC    ENDIF
 
 
C  CALCULATION OF CENTRAL RAPIDITY WITHOUT (ANTI)LEADER
C  MULTIPLICITY IN CENTER OF RAPIDITY DISTRIBUTION
      IZN = 0.D0
      IF ( IDIF .EQ. 0 ) THEN
        DO 111  I = 3,NTOT
          IF ( ABS(YR(I)) .LT. DELRAP ) IZN = IZN + 1
 111    CONTINUE
        IF ( IZN .LT. 1 ) THEN
          IF ( ISEL .EQ. 0 ) GOTO 30
C  IN CASE OF FEW PARTICLES, SET PARTICLE NUMBER IN PLATEAU TO 1
          IZN = 1
        ENDIF
C  CENTRAL RAPIDITY DENSITY FOR CHARGED PARTICLES
        IF ( NTOTEM .GE. 1 ) THEN
          ZNC = MAX( 1.1D0, DBLE(NCH)*IZN/(DBLE(NTOTEM)*2.D0*DELRAP) )
        ELSE
          ZNC = 1.1D0
        ENDIF
      ELSE
C  DIFFRACTION: SHIFT RAPIDITIES + TAKE CENT.RAP.DENS. FROM PARAMETRISAT
        DO 112  I = 3,NTOT
          YR(I) = YR(I) + YY0
 112    CONTINUE
        ZNC = MAX( 1.1D0, DCN2 )
      ENDIF
 
C  ZN ACCOUNTS FOR THE RISE OF PT WITH CENTRAL RAP.DENSITY. THE FORMULA
C  IS A FIT TO UA1 VALUES OF ARNISON ET AL, PHYS.LETT.B118(1982)167
C  REGARD, THAT OUR ZN IS DEFINED DIFFERENT FROM LITERATURE N BY 1
C - - - - - -
C  MODIFICATION AFTER J.N. CAPDEVIELLE, (DEC.96)
      IF ( ECMDPM .LE. 500.D0 ) THEN
        ZN = MAX( 1.00001D0, 3.669D0 / ZNC**0.435D0 + 6.4D0 )
      ELSE
C  TAKE INTO ACCOUNT THE RESULTS OF UA1/MIMI EXPERIMENT
        IF ( ZNC .GE. 3.D0 ) THEN
          ZN = 1.D0 /(ZNC*0.004111D0 + 0.145D0)+3.D0
        ELSE
C  FOR ROCH < 3.00 (MIMI)   (TO BE USED IN  PTRAM)
          ZM = 0.0033D0 * (ZNC-1.56D0)**2 + 0.406D0
          ZN = 2.64D0/ZM + 3.D0
        ENDIF
      ENDIF
C - - - - - -
C  NOW SET PARTICLE TYPE AND TRANSV. MOMENTA FOR NEW PARTICLES IN PPARAM
C  SET ALSO TRANSVERSE MASS FOR ALL PARTICLES (INCL. LEADER+ANTILEADER)
      CALL PPARAM
 
      IF ( IDIF .EQ. 0 ) THEN
C  NOW SET THE RAPIDITY OF THE ANTILEADER ACCORDING TO THE DISTRIBUTION
C  IN THE FEYNMAN X VARIABLE; SET THE RAPIDITY OF LEADER TO CONSUME
C  THE REMAINDER OF ENERGY
        CALL LEDENY( LEDEFL )
        IF ( LEDEFL .NE. 0 ) THEN
          IF ( DEBUG ) WRITE(MDEBUG,*) '   LEDEFL=',LEDEFL
          GOTO 30
        ENDIF
 
C  CALCULATE FOR SINGLE COLLISION SYSTEM C.M. ENERGY + RAPIDITY SHIFT
        IF ( GNU .LE. 1.D0 ) THEN
          JGNU   = 0.D0
          DYGNU  = 0.D0
          ECMJAD = ECMDPM
        ELSE
C  MULTIPLE COLLISION IN TARGET
          JGNU   = NINT(GNU-1.D0)
C  ADD ADDITIONALLY INTERACTING
C  TARGET NUCLEONS TO GET CORRECT JADACH FILTERING
C  CHOSE RANDOMLY WHETHER PROTON OR NEUTRON
          CALL RMMAR( RD,JGNU,1 )
          IPR = 0
          INE = 0
          TARMAS = PAMA(ITYP(2))
          DO 114  I = 1,JGNU
            NTOT       = NTOT + 1
            IF ( RD(I) .LE. .5D0 ) THEN
              ITYP(NTOT) = 13
              INE = INE + 1
            ELSE
              ITYP(NTOT) = 14
              IPR = IPR + 1
            ENDIF
            TMAS(NTOT) = PAMA(ITYP(NTOT))
            TARMAS     = TARMAS + TMAS(NTOT)
            EA(NTOT)   = TMAS(NTOT)
            PX(NTOT)   = 0.D0
            PY(NTOT)   = 0.D0
            PT2(NTOT)  = 0.D0
 114      CONTINUE
 
C  CALCULATE C.M. ENERGY + RAPIDITY SHIFT
*         YCMGNU = 0.5D0 * LOG( (ELAB+TARMAS+PLAB)/(ELAB+TARMAS-PLAB) )
          YCMGNU = 0.5D0 * LOG( (EPLUSP**2     +TARMAS*EPLUSP)/
     *                          (PAMA(ITYPE)**2+TARMAS*EPLUSP) )
          DYGNU  = YCM - YCMGNU
 
C  CALCULATE RAPIDITIES OF ADDITIONALLY INTERACTING TARGET NUCLEONS
C  IN THE CM SYSTEM OF NUCLEON-NUCLEON SYSTEM
          DO 115  I = NTOT-JGNU+1,NTOT
            YR(I) = - YCM
 115      CONTINUE
C  SHIFT RAPIDITIES INTO CM SYSTEM OF GNU+1 MASSES
          DO 113  I = 1,NTOT
            YR(I) = YR(I) + DYGNU
 113      CONTINUE
 
C  CENTER OF MASS ENERGY OF 1 PROJECTILE AND GNU TARGET NUCLEONS TO
C  BE USED IN THE JADACH FILTER.
          ECMJAD = SQRT( PAMA(ITYPE)**2 + TARMAS**2 + 2.D0*TARMAS*ELAB )
 
        ENDIF
 
      ELSE
C  IN CASE OF DIFFRACTION SET THE RAPIDITY OF LEADER AND ANTILEADER
C  IN SUBROUTINE LEADDF
        DYGNU  = 0.D0
        ECMJAD = ECMDPM
        CALL LEADDF( IFLGLD )
        IF ( IFLGLD .NE. 0 ) THEN
          IF ( DEBUG ) WRITE(MDEBUG,*) '   IFLGLD=',IFLGLD
          GOTO 30
        ENDIF
      ENDIF
 
C  CORRECT THE RAPIDITIES TO CONSERVE LONGITUDINAL MOMENTA AND ENERGY
C  USING THE ALGORITHM OF JADACH (SIMPLIFIED VERSION BY R. ATTALLAH)
      CALL JADACH( ECMJAD,JADFLG )
      IF ( JADFLG .NE. 0 ) THEN
        IF ( DEBUG ) WRITE(MDEBUG,*) '   JADFLG=', JADFLG
        IF ( JADFLG .GT. 0 ) GOTO 30
        IF ( JADFLG .LT. 0 ) THEN
          NREPRD = NREPRD + 1
          IF ( NREPRD .GT. 10 ) GOTO 1
          GOTO 1919
        ENDIF
      ENDIF
 
 
C  CALCULATE LAB ENERGIES OF SECONDARY PARTICLES FROM THE RAPIDITIES
C  INCLUDING THE ADDITIONAL TARGET NUCLEONS
      ETOT = 0.D0
      DO 135  I = 1,NTOT
        YR(I) = YR(I) + YCM - DYGNU
        EA(I) = TMAS(I) * COSH( YR(I) )
        ETOT  = ETOT + EA(I)
 135  CONTINUE
 
      IF ( DEBUG ) WRITE(MDEBUG,136)
     *             (I,ITYP(I),PX(I),PY(I),YR(I),EA(I),I=1,NTOT)
 136  FORMAT('  NO ITYP     PX          PY            YR          EA'/
     *      (' ',I4,I3,1X,1P,4G13.6) )
 
C-----------------------------------------------------------------------
C  LOOP OVER ALL SECONDARY PARTICLES AND THE LEADING PARTICLE
C  PUT THEM ON THE STACK
      DO 139  LK = 5,8
        SECPAR(LK) = CURPAR(LK)
 139  CONTINUE
 
C  PROCESS LOOP
      DO 140  J = 1,NTOT
C  REJECTION OF BACKWARD GOING PARTICLES
        IF ( YR(J) .LE. 0.D0 ) THEN
          IF ( DEBUG ) WRITE(MDEBUG,*)'HDPM  : YR REJECT PARTICLE ',J
          GOTO 140
        ENDIF
C  CALCULATE THE PROPERTIES OF ALL SECONDARIES
C  PARTICLE TYPE
        SECPAR(1) = ITYP(J)
C  CALCULATE GAMMA FACTOR
        SECPAR(2) = EA(J) / PAMA(ITYP(J))
C  TOTAL MOMENTUM SQUARED
        PTM = EA(J)**2 - PAMA(ITYP(J))**2
        IF ( PT2(J) .GT. PTM ) THEN
          IF ( DEBUG ) WRITE(MDEBUG,*)'HDPM  : PT REJECT PARTICLE ',J
          GOTO 140
        ENDIF
C  EMISSION ZENITH ANGLE AGAINST TRAJECTORY OF PROJECTILE
        IF ( PTM .EQ. 0.D0 ) THEN
          COSTET = 1.D0
        ELSE
          COSTET = SQRT( 1.D0 - PT2(J) / PTM )
        ENDIF
C  EMISSION AZIMUTH ANGLE
        IF ( PX(J) .NE. 0.D0  .OR.  PY(J) .NE. 0.D0 ) THEN
          PHIJ = ATAN2( PY(J), PX(J) )
        ELSE
          PHIJ = 0.D0
        ENDIF
        CALL ADDANG( COSTHE,PHI, COSTET,PHIJ, SECPAR(3),SECPAR(4) )
        IF ( SECPAR(3) .LT. C(29) ) THEN
C  OMIT UPWARD GOING PARTICLES
          IF (DEBUG) WRITE(MDEBUG,*)'HDPM  : ANGLE REJECT PARTICLE ',J
          GOTO 140
        ENDIF
C  PUT SECONDARY PARTICLES ON STACK, IF NOT GOING UPWARDS
        IF ( J .GT. 2 ) THEN
          CALL TSTACK
        ELSE
C  PUT LEADER OR ANTI-LEADER ON STACK, IF NOT GOING UPWARDS
          IF ( ITYP(J) .GT. 50 ) THEN
C  LEADER OR ANTI LEADER ARE RESONANCES AND DECAY
            IRESPAR = IRESPAR + 1
            IF ( IRESPAR .GE. 1000 ) THEN
              WRITE(MONIOU,*) 'STACK OF RESDEC RANDOM NUMBERS FULL'
              IRESPAR = 999
            ENDIF
            RESRAN(IRESPAR) = RDRES(J)
C  COUNTER FOR ENERGY-MULTIPLICITY MATRIX
            MSMM = MSMM + 1
          ENDIF
          CALL TSTACK
 
C  CALCULATE ELASTICITY FROM ENERGY OF LEADER (REST OF RESONANCE DECAY)
          IF ( J. EQ. 1 ) ELASTI = SECPAR(2)*PAMA(NINT(SECPAR(1)))/ELAB
        ENDIF
C  COUNTERS FOR FIRST INTERACTION
        IF ( FIRSTI ) THEN
          IF     ( SECPAR(1) .EQ.  7.D0  .OR.  SECPAR(1) .EQ.  8.D0
     *       .OR.  SECPAR(1) .EQ.  9.D0                           ) THEN
            IFINPI = IFINPI + 1
          ELSEIF ( SECPAR(1) .EQ. 13.D0  .OR.  SECPAR(1) .EQ. 14.D0
     *       .OR.  SECPAR(1) .EQ. 15.D0 .OR. SECPAR(1) .EQ. 25.D0 ) THEN
            IFINNU = IFINNU + 1
          ELSEIF ( SECPAR(1) .EQ. 10.D0  .OR.  SECPAR(1) .EQ. 11.D0
     *       .OR.  SECPAR(1) .EQ. 12.D0 .OR. SECPAR(1) .EQ. 16.D0 ) THEN
            IFINKA = IFINKA + 1
          ELSEIF ( SECPAR(1) .GE. 71.D0 .AND. SECPAR(1) .LE. 74.D0) THEN
            IFINET = IFINET + 1
          ELSEIF ((SECPAR(1) .GE. 18.D0 .AND. SECPAR(1) .LE. 24.D0)
     *       .OR. (SECPAR(1) .GE. 26.D0 .AND. SECPAR(1) .LE. 32.D0))THEN
            IFINHY = IFINHY + 1
          ENDIF
        ENDIF
 
 140  CONTINUE
 
C  COUNTER FOR ENERGY-MULTIPLICITY MATRIX
      MSMM = MSMM + NTOT - 2
 
C  FILL ELASTICITY IN MATRICES
      MEL = MIN ( 1.D0+10.D0*      MAX( 0.D0, ELASTI ) , 11.D0 )
      MEN = MIN ( 4.D0+ 3.D0*LOG10(MAX( .1D0, EKINL  )), 37.D0 )
      IELDPM(MEN,MEL) = IELDPM(MEN,MEL) + 1
      IELDPA(MEN,MEL) = IELDPA(MEN,MEL) + 1
      IF ( ELASTI .LT. 1.D0 ) THEN
        ELMEAN(MEN) = ELMEAN(MEN) + ELASTI
        ELMEAA(MEN) = ELMEAA(MEN) + ELASTI
      ENDIF
 
      IF ( FIRSTI ) THEN
        ELAST  = ELASTI
        FIRSTI = .FALSE.
      ENDIF
 
      IF ( DEBUG ) WRITE(MDEBUG,*) 'HDPM  : ELAST=',SNGL(ELASTI),
     *                            SNGL(ETOT),SNGL(ELAB)
 
      RETURN
      END
