*
*=== Ustsuw ===========================================================*
*
      PROGRAM USTSUW
*
*----------------------------------------------------------------------*
*                                                                      *
*     Copyright (C) 2003-2019:  CERN & INFN                            *
*     Copyright (C) 2020-    :  CERN                                   *
*     All Rights Reserved.                                             *
*                                                                      *
*     USer Tracklength or collision scoring SUm and Weighting:         *
*                                                                      *
*     Created on   07 July 1990    by    Alfredo Ferrari & Paola Sala  *
*                                                   Infn - Milan       *
*     Modified on  26 June 2021                                        *
*                                                                      *
*----------------------------------------------------------------------*
*
      CHARACTER FILE*232,FILE1*232,RUNTIT*80,RUNTIM*32,CHSTAT*10,
     &          CHCTOT*15
*
*=== usrtrc ==========================================================*
*
*---------------------------------------------------------------------*
*     Module usrtrc:                                                  *
*     A. Ferrari : user defined track-length & collision density      *
*                  estimators                                         *
*                                                                     *
*     Up to MXUSTC user defined track or coll are allowed             *
*            lusrtc = logical flag, .true. if at least 1 user defined *
*                     track or coll is used                           *
*            lustrk = logical flag, .true. if at least 1 user defined *
*                     track is used                                   *
*            luscll = logical flag, .true. if at least 1 user defined *
*                     coll is used                                    *
*            nusrtc = number of user defined track & coll used        *
*            nustrk = number of user defined track used               *
*            nuscll = number of user defined coll used                *
*            iustrk = list of track-length estimators                 *
*            iuscll = list of collision density track estimators      *
*                     2 = linear energy, coll, -1 = logarithmic ener- *
*                         gy, track, -2 = logarithmic energy, coll    *
*            itustc = type of binning: 1 = linear energy, track       *
*                     2 = linear energy, coll, -1 = logarithmic ener- *
*                         gy, track, -2 = logarithmic energy, coll    *
*            idustc = distribution to be scored: the usual values are *
*                     allowed                                         *
*            nrustc = region                                          *
*            vusrtc = volume (cm**3) of the detector                  *
*            llnutc = no low energy neutron scoring if false, yes if  *
*                     true                                            *
*            titutc = track or coll name                              *
*            ipustc = logical unit to print the results on: formatted *
*                     if > 0, unformatted if < 0                      *
*            igmutc = maximum low energy neutron group to be scored   *
*            kbustc = initial location in blank common of the consi-  *
*                     dered track or coll (real*8 address)            *
*                     I T   I S   F O R  A D D R E S S  OF  1 !!!!!   *
*            netcbn = number of energy intervals                      *
*     etclow/etchgh = minimum and maximum energies                    *
*            detcbn = energy bin width                                *
*                                                                     *
*---------------------------------------------------------------------*
*
      PARAMETER ( MXUSTC = 2500 )
      LOGICAL LUSRTC, LUSTRK, LUSCLL, LLNUTC
      CHARACTER*10 TITUTC
      COMMON /USRTCC/ ETCLOW(MXUSTC), ETCHGH(MXUSTC), DETCBN(MXUSTC),
     &                VUSRTC(MXUSTC),
     &                IUSTRK(MXUSTC), IUSCLL(MXUSTC), NETCBN(MXUSTC),
     &                NRUSTC(MXUSTC), ITUSTC(MXUSTC), IDUSTC(MXUSTC),
     &                KBUSTC(MXUSTC), IPUSTC(MXUSTC), IGMUTC(MXUSTC),
     &                LLNUTC(MXUSTC), NUSRTC, NUSTRK, NUSCLL, LUSRTC,
     &                LUSTRK, LUSCLL
      COMMON /USTCH/  TITUTC(MXUSTC)

      COMMON /TCHELP/ ATCLOW(MXUSTC), ATCHGH(MXUSTC), DATCBN(MXUSTC),
     &                NATCBN(MXUSTC)

      PARAMETER ( PIPIPI = 3.141592653589793D+00 )
      PARAMETER ( TWOPIP = 2.D+00 * PIPIPI )

      REAL,DIMENSION(:),ALLOCATABLE :: GMSTOR, EPGMAX, EPGMCK, ACCUMH,
     &                                 ENGMAX
      DOUBLE PRECISION,DIMENSION(:),ALLOCATABLE :: GDSTOR, GBSTOR,
     &                                             ACCUMG, ACUMCU
      DOUBLE PRECISION ACCUMP, COSTHE, AFACT , TOTTOT, TOTERR, DSCORE,
     &                 DSCOR2, DELTAE, DOMEGA, ACCUMT, ACCUME, WCTOT ,
     &                 GBHELP, TOTHLP
      DIMENSION TOTTOT(MXUSTC), TOTERR(MXUSTC), KFUSTC (MXUSTC),
     &          JX (MXUSTC)
      LOGICAL LSTATI
*
*     Note that all addresses are for zero index (KBUSTC,KFUSTC)
*
      DO 4444 I = 1, MXUSTC
         TOTTOT (I) = 0.D+00
         TOTERR (I) = 0.D+00
 4444 CONTINUE
      NCTOT   = 0
      MCTOT   = 0
      WCTOT   = 0.D+00
      MBATCH  = 0
      MAXSIZE = 0
      MNETCBN = 0
      MIGMUTC = 0
      IFILE   = 0
      LSTATI  = .FALSE.
      LQOLD   = 0
 1111 CONTINUE
      WRITE(*,'('' Type the input file: '',$)')
      READ (*,'(A)')FILE
      LQ     = LNNBLN(FILE)
      IF (LQ .LE. 0) GO TO 2000
      OPEN (UNIT=1,FILE=FILE,STATUS='OLD',FORM='UNFORMATTED',ERR=1111)
      IFILE  = IFILE + 1
      IRECRD = 0
      READ (1,ERR=100) RUNTIT,RUNTIM,WEIPRI,NCASE,MCASE,NBATCH
      IRECRD = IRECRD + 1
      WRITE(*,*)RUNTIT
      WRITE(*,*)RUNTIM
      WRITE(*,*)WEIPRI
      WRITE(*,*)NCASE
      WRITE(*,*)MCASE
      WRITE(*,*)NBATCH
      GO TO 200
 100  CONTINUE
      REWIND 1
      IRECRD = 0
      READ (1,ERR=110) RUNTIT,RUNTIM,WEIPRI,NCASE,NBATCH
      IRECRD = IRECRD + 1
      WRITE(*,*)RUNTIT
      WRITE(*,*)RUNTIM
      WRITE(*,*)WEIPRI
      WRITE(*,*)NCASE
      MCASE  = 0
      GO TO 200
 110  CONTINUE
      REWIND 1
      IRECRD = 0
      READ (1,ERR=120) RUNTIT,RUNTIM,WEIPRI,NCASE
      IRECRD = IRECRD + 1
      WRITE(*,*)RUNTIT
      WRITE(*,*)RUNTIM
      WRITE(*,*)WEIPRI
      WRITE(*,*)NCASE
      MCASE  = 0
      NBATCH = 1
      GO TO 200
 120  CONTINUE
      REWIND 1
      IRECRD = 0
      READ (1,ERR=150) RUNTIT,RUNTIM,NCASE
      IRECRD = IRECRD + 1
      WRITE(*,*)RUNTIT
      WRITE(*,*)RUNTIM
      WRITE(*,*)NCASE
      WEIPRI = REAL(NCASE)
      MCASE  = 0
      NBATCH = 1
      GO TO 200
 150  CONTINUE
      REWIND 1
      IRECRD = 0
      WRITE(*,*)' How many particle for this run?'
      READ(*,*)NCASE
      WEIPRI = REAL(NCASE)
      MCASE  = 0
      NBATCH = 1
 200  CONTINUE
      NCTOT  = NCTOT + NCASE
      MCTOT  = MCTOT + MCASE
      IF ( NCTOT .GE. 1000000000 ) THEN
         NCTOT = NCTOT - 1000000000
         MCTOT = MCTOT + 1
      END IF
      WCTOT  = WCTOT  + DBLE (WEIPRI)
      MBATCH = MBATCH + NBATCH

*  +-------------------------------------------------------------------*
*  | First pass to determine array size needed
      IF (MAXSIZE .EQ. 0) THEN
         IRKEEP = IRECRD
         DO IX = 1, MXUSTC
            NTC = IX
            READ (1,ERR=310,END=390)
     &                   MTC, TITUTC(NTC), ITUSTC(NTC), IDUSTC(NTC),
     &                   NRUSTC(NTC), VUSRTC (NTC), LLNUTC (NTC),
     &                   ETCLOW(NTC), ETCHGH(NTC),
     &                   NETCBN(NTC), DETCBN(NTC)
**   &                   ATCLOW(NTC), ATCHGH(NTC),
**   &                   NATCBN(NTC), DATCBN(NTC)
            IRECRD = IRECRD + 1
            IF ( LLNUTC (NTC) ) THEN
               READ (1) IGMUTC(NTC)
               IRECRD = IRECRD + 1
            ELSE
               IGMUTC (NTC) = 0
            END IF
            GO TO 350

 310        CONTINUE
*           BACKSPACE 1
            REWIND 1
            DO IR = 1, IRECRD
               READ (1)
            END DO
            READ (1,ERR=320,END=390)
     &                   MTC, TITUTC(NTC), ITUSTC(NTC), IDUSTC(NTC),
     &                   NRUSTC(NTC), VUSRTC (NTC),
     &                   ETCLOW(NTC), ETCHGH(NTC),
     &                   NETCBN(NTC), DETCBN(NTC)
**   &                   ATCLOW(NTC), ATCHGH(NTC),
**   &                   NATCBN(NTC), DATCBN(NTC)
            IRECRD = IRECRD + 1
            IGMUTC (NTC) = 0
            GO TO 350

 320        CONTINUE
*           BACKSPACE 1
            REWIND 1
            DO IR = 1, IRECRD
               READ (1)
            END DO
            READ (1,ERR=330) CHSTAT,IJX
            IRECRD = IRECRD + 1
            IF ( CHSTAT .EQ. 'STATISTICS') GO TO 390

 330        CONTINUE
            STOP 'UNKNOWN RECORD TYPE'

 350        CONTINUE

            NATCBN(NTC)=1

            READ (1)
            IRECRD = IRECRD + 1

            MAXSIZE = MAXSIZE + ( 2 + NATCBN(NTC) ) *
     &                   ( NETCBN(NTC) + IGMUTC(NTC) )
            MNETCBN = MAX ( MNETCBN, NETCBN(NTC) )
            MIGMUTC = MAX ( MIGMUTC, IGMUTC(NTC) )

         END DO
 390     CONTINUE

         ALLOCATE( GMSTOR(MAXSIZE) )
         ALLOCATE( GDSTOR(MAXSIZE) )
         ALLOCATE( GBSTOR(MAXSIZE) )

         GDSTOR = 0.D+00
         GBSTOR = 0.D+00

         ALLOCATE( ENGMAX(MIGMUTC + 1) )
         ALLOCATE( EPGMAX(MNETCBN + MIGMUTC + 1) )
         ALLOCATE( EPGMCK(MNETCBN + MIGMUTC + 1) )
         ALLOCATE( ACCUMG(MNETCBN + MIGMUTC + 1) )
         ALLOCATE( ACUMCU(MNETCBN + MIGMUTC + 1) )
         ALLOCATE( ACCUMH(MNETCBN + MIGMUTC + 1) )

         REWIND 1
         DO IR = 1, IRKEEP
             READ (1)
         END DO
         IRECRD = IRKEEP

      END IF
*  |
*  +-------------------------------------------------------------------*

      KLAST  = 0
      KMAX   = 0
*  +-------------------------------------------------------------------*
*  |                            Read from binary file the results!!!
      DO 6600 IX = 1, MXUSTC
         NTC = IX
         READ (1,ERR=6400,END=1000)
     &                MTC, TITUTC(NTC), ITUSTC(NTC), IDUSTC(NTC),
     &                NRUSTC(NTC), VUSRTC (NTC), LLNUTC (NTC),
     &                ETCLOW(NTC), ETCHGH(NTC),
     &                NETCBN(NTC), DETCBN(NTC)
**   &                ATCLOW(NTC), ATCHGH(NTC),
**   &                NATCBN(NTC), DATCBN(NTC)
         IRECRD = IRECRD + 1
         IF ( LLNUTC (NTC) ) THEN
            READ (1) IGMUTC(NTC), (ENGMAX (IG), IG=1,IGMUTC(NTC)+1)
            IRECRD = IRECRD + 1
         ELSE
            IGMUTC (NTC) = 0
         END IF
         GO TO 6500
6400     CONTINUE
*        BACKSPACE 1
         REWIND 1
         DO IR = 1, IRECRD
            READ (1)
         END DO
         READ (1,ERR=6450,END=1000)
     &                MTC, TITUTC(NTC), ITUSTC(NTC), IDUSTC(NTC),
     &                NRUSTC(NTC), VUSRTC (NTC),
     &                ETCLOW(NTC), ETCHGH(NTC),
     &                NETCBN(NTC), DETCBN(NTC)
**   &                ATCLOW(NTC), ATCHGH(NTC),
**   &                NATCBN(NTC), DATCBN(NTC)
         IRECRD = IRECRD + 1
         LLNUTC (NTC) = .FALSE.
         IGMUTC (NTC) = 0
         GO TO 6500
6450     CONTINUE
*        BACKSPACE 1
         REWIND 1
         DO IR = 1, IRECRD
            READ (1)
         END DO
         READ (1,ERR=6460,END=1000) CHSTAT,IJX
         IRECRD = IRECRD + 1
         IF ( CHSTAT .EQ. 'STATISTICS') THEN
            LSTATI = .TRUE.
            GO TO 1000
         END IF
6460     CONTINUE
         STOP 'UNKNOWN RECORD TYPE'
6500     CONTINUE
******
         ATCLOW(NTC)=0.D+00
         ATCHGH(NTC)=1.D+00
         DATCBN(NTC)=1.D+00
         NATCBN(NTC)=1
******
         JX (IX) = NTC
***************************************************************
*********** Note KBUSTC is used for 0 address *****************
***************************************************************
         KBUSTC (NTC) = KLAST
         K0 = KLAST + 1
         K1 = NETCBN (NTC) * NATCBN (NTC) + K0 - 1
         K2 = K1 + IGMUTC (NTC) * NATCBN (NTC)
         KLAST = K2
         READ (1) (GMSTOR(J), J = K0, K2)
         IRECRD = IRECRD + 1
         DO J = K0, K2
            GDSTOR (J) = GDSTOR (J) + DBLE (GMSTOR(J))
     &                 * ( DBLE (NCASE) + DBLE (MCASE) * 1.D+09 )
            GBSTOR (J) = GBSTOR (J) + DBLE (GMSTOR(J))**2
     &                 * ( DBLE (NCASE) + DBLE (MCASE) * 1.D+09 )
         END DO
         KMAX = MAX ( KMAX, K2 )
6600  CONTINUE
*  |
*  +-------------------------------------------------------------------*
1000  CONTINUE
      IX = IX-1
      CLOSE (UNIT=1)
*  +-------------------------------------------------------------------*
*  |
      IF ( MBATCH .EQ. 1 .OR. IFILE .EQ. 1 ) THEN
         DO IIX = 1, IX
            NTC = JX (IIX)
            KFUSTC (NTC) = KLAST
            KLAST = KLAST + 2 * ( NETCBN (NTC) + IGMUTC (NTC) )
         END DO
         FILE1 = FILE
         LQOLD = LQ
      END IF
*  |
*  +-------------------------------------------------------------------*
      KLAST = 0
*  +-------------------------------------------------------------------*
*  |
      DO 7000 IIX = 1, IX
         NTC = JX (IIX)
         ACCUMP = 0.D+00
         DO IG = 1, NETCBN (NTC) + IGMUTC (NTC)
            ACCUMG (IG) = 0.D+00
            ACUMCU (IG) = 0.D+00
         END DO
         AFACT  = DBLE (NCASE) + DBLE (MCASE) * 1.D+09
         K0     = KLAST + 1
         K1     = NETCBN (NTC) * NATCBN (NTC) + K0 - 1
         K2     = K1 + IGMUTC (NTC) * NATCBN (NTC)
         KLAST  = K2
         ENERG0 = ETCLOW (NTC)
*  |  +----------------------------------------------------------------*
*  |  |
         IF ( ITUSTC (NTC) .LT. 0 ) THEN
            EFACT  = DETCBN (NTC)
*  |  |
*  |  +----------------------------------------------------------------*
*  |  |
         ELSE
            DELTAE = DBLE (DETCBN (NTC))
         END IF
*  |  |
*  |  +----------------------------------------------------------------*
         ENEOLD = ENERG0
*  |  +----------------------------------------------------------------*
*  |  |
*        IF ( ABS ( ITUSTC (NTC) ) .GE. 2 ) THEN
         IF ( ABS ( ITUSTC (NTC) ) .GE. 3 ) THEN
            OMEGA0 = ATCLOW (NTC)
            GFACT  = DATCBN (NTC)
*  |  |
*  |  +----------------------------------------------------------------*
*  |  |
         ELSE
            OMEGA0 = ATCLOW (NTC)
            DOMEGA = DBLE (DATCBN (NTC))
         END IF
*  |  |
*  |  +----------------------------------------------------------------*
         YGMAX = -1.E+30
         YGMIN = 1.E+30
*  |  +----------------------------------------------------------------*
*  |  |
         DO 6540 IE = 1, NETCBN (NTC) + IGMUTC (NTC)
*  |  |  +-------------------------------------------------------------*
*  |  |  |
            IF ( IE .LE. NETCBN (NTC) ) THEN
*  |  |  |  +----------------------------------------------------------*
*  |  |  |  |
               IF ( ITUSTC (NTC) .LT. 0 ) THEN
                  ENERGE = ENERG0 * EFACT**IE
                  DELTAE = DBLE (ENERGE - ENEOLD)
*  |  |  |  |
*  |  |  |  +----------------------------------------------------------*
*  |  |  |  |
               ELSE
                  ENERGE = ETCLOW (NTC) + DETCBN (NTC) * FLOAT (IE)
               END IF
*  |  |  |  |
*  |  |  |  +----------------------------------------------------------*
*  |  |  |
*  |  |  +-------------------------------------------------------------*
*  |  |  |
            ELSE
               IF ( IE-NETCBN (NTC) .EQ. 1) ENEOLD = ENGMAX (1)
               ENERGE = ENGMAX (IE-NETCBN(NTC)+1)
               DELTAE = DBLE(ENEOLD - ENERGE)
            END IF
*  |  |  |
*  |  |  +-------------------------------------------------------------*
            ENEOLD = ENERGE

            IF (DELTAE.LE. 0.D+00)WRITE(*,*)' **** DELTAE',IE,NTC

*  |  |  +-------------------------------------------------------------*
*  |  |  |
*           IF ( ABS (ITUSTC (NTC)) .GE. 2 ) THEN
            IF ( ABS (ITUSTC (NTC)) .GE. 3 ) THEN
               OMEOLD = 0.0
*  |  |  |
*  |  |  +-------------------------------------------------------------*
*  |  |  |
            ELSE
               OMEOLD = OMEGA0
            END IF
*  |  |  |
*  |  |  +-------------------------------------------------------------*
*  |  |  +-------------------------------------------------------------*
*  |  |  |
            DO 6520 IA = 1, NATCBN (NTC)
*  |  |  |  +----------------------------------------------------------*
*  |  |  |  |
*              IF ( ABS (ITUSTC (NTC)) .GE. 2 ) THEN
               IF ( ABS (ITUSTC (NTC)) .GE. 3 ) THEN
                  OMEGAA = OMEGA0 * GFACT**(IA - 1)
                  DOMEGA = DBLE (OMEGAA - OMEOLD)
*  |  |  |  |
*  |  |  |  +----------------------------------------------------------*
*  |  |  |  |
               ELSE
                  OMEGAA = OMEGA0 + DOMEGA * FLOAT (IA)
               END IF
*  |  |  |  |
*  |  |  |  +----------------------------------------------------------*
               OMEOLD = OMEGAA

               IF (DOMEGA.LE. 0.D+00)WRITE(*,*)' **** DOMEGA',IE,NTC

*  |  |  |  +----------------------------------------------------------*
*  |  |  |  |
               IF ( IE .LE. NETCBN (NTC) ) THEN
                  IEA = IE + NETCBN (NTC) * ( IA - 1 ) + K0 - 1
*  |  |  |  |
*  |  |  |  +----------------------------------------------------------*
*  |  |  |  |
               ELSE
                  IG  = IE - NETCBN (NTC)
                  IEA = IG + IGMUTC (NTC) * ( IA - 1 ) + K1
               END IF
*  |  |  |  |
*  |  |  |  +----------------------------------------------------------*
               DSCORE = DBLE ( GMSTOR (IEA) ) * ( DELTAE * DOMEGA )
               DSCOR2 = DBLE ( GMSTOR (IEA) ) * DOMEGA
               ACCUMP = ACCUMP + DSCORE
               ACUMCU (IE) = ACUMCU (IE) + DSCORE
               ACCUMG (IE) = ACCUMG (IE) + DSCOR2
               YGMAX = MAX ( YGMAX, GMSTOR (IEA) )
               YGMIN = MIN ( YGMIN, GMSTOR (IEA) )
6520        CONTINUE
*  |  |  |
*  |  |  +-------------------------------------------------------------*
6540     CONTINUE
*  |  |
*  |  +----------------------------------------------------------------*
         TOTTOT (NTC) = TOTTOT (NTC) + ACCUMP * AFACT
         TOTERR (NTC) = TOTERR (NTC) + ACCUMP * ACCUMP * AFACT

         ACCUMT=0.D0
         DO IE = 1,NETCBN(NTC)+IGMUTC(NTC)
            ACCUMT = ACCUMT + ACUMCU (IE)
         END DO
         IF ( ABS ( ACCUMT - ACCUMP ) .GT.
     &        1.D-13 * ACCUMP ) WRITE (*,*)' *** ERROR 0', NTC,
     &        ACCUMT, ACCUMP
         ACCUME=0.D0
         DO IE = 1,NETCBN(NTC)
            ACCUME = ACCUME + ACUMCU (IE)
         END DO

*  |  +----------------------------------------------------------------*
*  |  |
         DO IE = NETCBN (NTC) - 1, 1, -1
            ACUMCU (IE) = ACUMCU (IE) + ACUMCU (IE+1)
         END DO
*  |  |
*  |  +----------------------------------------------------------------*
        IF ( NETCBN(NTC) .GT. 0 )
     &      ACUMCU (NETCBN(NTC)+1) = ACUMCU (NETCBN(NTC)+1) + ACUMCU (1)
*  |  +----------------------------------------------------------------*
*  |  |
         DO IE = NETCBN (NTC) + 2, NETCBN (NTC) + IGMUTC (NTC)
            ACUMCU (IE) = ACUMCU (IE) + ACUMCU (IE-1)
         END DO
*  |  |
*  |  +----------------------------------------------------------------*

         IF ( LLNUTC (NTC) .AND.
     &      ABS ( ACUMCU ( NETCBN (NTC) + IGMUTC (NTC) ) - ACCUMP ) .GT.
     &      1.D-13 * ACCUMP ) WRITE (*,*)' *** ERROR', NTC,
     &      ACUMCU ( NETCBN (NTC) + IGMUTC (NTC) ), ACCUMP
*  |  +----------------------------------------------------------------*
*  |  |
         DO IE = 1, NETCBN (NTC) + IGMUTC (NTC)
            KK0 = KFUSTC (NTC) + IE
            GDSTOR (KK0) = GDSTOR (KK0) + ACCUMG (IE) * AFACT
            GBSTOR (KK0) = GBSTOR (KK0) + ACCUMG (IE) * ACCUMG (IE) *
     &                     AFACT
            KK0 = KK0 + NETCBN (NTC) + IGMUTC (NTC)
            GDSTOR (KK0) = GDSTOR (KK0) + ACUMCU (IE) * AFACT
            GBSTOR (KK0) = GBSTOR (KK0) + ACUMCU (IE) * ACUMCU (IE) *
     &                     AFACT
         END DO
*  |  |
*  |  +----------------------------------------------------------------*
         KMAX = MAX ( KMAX, KK0 )
7000  CONTINUE
*  |
*  +-------------------------------------------------------------------*
      GO TO 1111
2000  CONTINUE
      DO J = 1, KMAX
         GDSTOR (J) = GDSTOR (J)
     &              / ( DBLE (NCTOT) + DBLE (MCTOT) * 1.D+09 )
         GBSTOR (J) = GBSTOR (J)
     &              / ( DBLE (NCTOT) + DBLE (MCTOT) * 1.D+09 )
      END DO
      DO IIX = 1, IX
         NTC = JX (IIX)
         TOTTOT (NTC) = TOTTOT (NTC)
     &                / ( DBLE (NCTOT) + DBLE (MCTOT) * 1.D+09 )
         TOTERR (NTC) = TOTERR (NTC)
     &                / ( DBLE (NCTOT) + DBLE (MCTOT) * 1.D+09 )
      END DO
*  +-------------------------------------------------------------------*
*  |  Real sum of Mbatch runs:
      IF ( MBATCH - 1 .GT. 0  .AND. ( .NOT. LSTATI .OR. IFILE .GT. 1 )
     &   ) THEN
         DO J = 1, KMAX
            GBHELP     = SQRT (GBSTOR(J))
            GBSTOR (J) = ( ( GBHELP - GDSTOR (J) )
     &                   * ( GBHELP + GDSTOR (J) ) )
     &                 / DBLE (MBATCH-1)
            IF ( ABS (GBSTOR(J)) .LT. 1.D-12 * GDSTOR(J)**2 )
     &           GBSTOR (J) = 0.D+00
            IF ( GDSTOR (J) .GT. 0.D+00 ) THEN
               GBSTOR (J) = MIN ( SQRT ( GBSTOR (J) ) / GDSTOR(J),
     &                            0.99D+00 )
            END IF
         END DO
         DO IIX = 1, IX
            NTC = JX (IIX)
            TOTHLP = SQRT ( TOTERR (NTC) )
            TOTERR (NTC) = ( ( TOTHLP - TOTTOT (NTC) )
     &                     * ( TOTHLP + TOTTOT (NTC) ) )
     &                   / DBLE (MBATCH-1)
            IF ( ABS ( TOTERR (NTC)) .LT. 1.D-12 * TOTTOT (NTC)**2 )
     &              TOTERR (NTC) = 0.D+00
            IF ( TOTTOT (NTC) .GT. 0.D+00 ) THEN
               TOTERR (NTC) = MIN ( SQRT ( TOTERR(NTC) ) / TOTTOT (NTC),
     &                              0.99D+00 )
            END IF
         END DO
*  |
*  +-------------------------------------------------------------------*
*  |  Only one run, check if it is already a sum with statistics:
      ELSE IF ( LSTATI ) THEN
         OPEN ( UNIT=11, FILE=FILE1, STATUS='OLD', FORM='UNFORMATTED' )
         READ (11)
         DO 7050 IIX = 1, IX
            NTC = JX (IIX)
            READ (11)
            IF ( LLNUTC (NTC) ) READ (11)
            READ (11)
 7050    CONTINUE
         READ (11) CHSTAT,IJX
         IF ( CHSTAT .NE. 'STATISTICS' .OR. IJX .NE. IX ) STOP 'STATI'
         KLAST = 0
         DO 7400 IIX = 1, IX
            NTC = JX (IIX)
            READ (11) TOTA, TOTB
            IF ( ABS (TOTA-TOTTOT(NTC)) .GT. 1.D-05*TOTA )
     &         WRITE (*,*)' *** Error Tota:',TOTA,TOTTOT(NTC), NTC
            TOTERR (NTC) = TOTB
            READ (11) NETCCK,IGMUCK,ETCLCK,(EPGMCK(IG),
     &                IG=1,NETCCK+IGMUCK+1)
            IF ( NETCCK .NE. NETCBN (NTC) )
     &         WRITE (*,*)' *** Error Netcck:',NETCCK,NETCBN(NTC)
            IF ( IGMUCK .NE. IGMUTC (NTC) )
     &         WRITE (*,*)' *** Error Igmuck:',IGMUCK,IGMUTC(NTC)
            IF ( ABS (ETCLCK-ETCLOW(NTC)) .GT. 1.E-05*ETCLCK )
     &         WRITE (*,*)' *** Error Etclck:',ETCLCK,ETCLOW(NTC)
            KK0 = 1
            KK1 = NETCBN (NTC)
            KK2 = NETCBN (NTC) + IGMUTC (NTC)
            READ (11) ( ACCUMH (J), J = KK1, KK0,-1 ),
     &                ( ACCUMH (J), J = KK1+1, KK2 )
            DO 7100 J = KK0, KK1
               IF ( ABS (ACCUMH(J)-GDSTOR(J+KFUSTC(NTC))) .GT.
     &              1.D-05*ACCUMH(J) )
     &         WRITE (*,*)' *** Error Gdstor:',ACCUMH(J),
     &                      GDSTOR(J+KFUSTC(NTC))
 7100       CONTINUE
            READ (11) ( ACCUMH (J), J = KK1, KK0,-1 ),
     &                ( ACCUMH (J), J = KK1+1, KK2 )
            DO 7150 J = KK0, KK2
               GBSTOR (J+KFUSTC(NTC)) = ACCUMH (J)
 7150       CONTINUE
            KFUHLP = KFUSTC (NTC) + NETCBN (NTC) + IGMUTC (NTC)
            READ (11) ( ACCUMH (J), J = KK1, KK0,-1 ),
     &                ( ACCUMH (J), J = KK1+1, KK2 )
            DO 7200 J = KK0, KK2
               IF ( ABS (ACCUMH(J)-GDSTOR(J+KFUHLP)) .GT.
     &              1.D-05*ACCUMH(J) )
     &         WRITE (*,*)' *** Error Gdstor_2:',ACCUMH(J),
     &                      GDSTOR(J+KFUHLP)
 7200       CONTINUE
            READ (11) ( ACCUMH (J), J = KK1, KK0,-1 ),
     &                ( ACCUMH (J), J = KK1+1, KK2 )
            DO 7250 J = KK0, KK2
               GBSTOR (J+KFUHLP) = ACCUMH (J)
 7250       CONTINUE
            K0 = KLAST + 1
            K1 = NETCBN (NTC) * NATCBN (NTC) + K0 - 1
            K2 = K1 + IGMUTC (NTC) * NATCBN (NTC)
            KLAST = K2
            READ (11) ( GMSTOR(J), J = K0, K2 )
            DO 7300 J = K0, K2
               GBSTOR (J) = DBLE (GMSTOR(J))
 7300       CONTINUE
 7400    CONTINUE
         KLAST = 0
         CLOSE ( UNIT=11 )
*  |
*  +-------------------------------------------------------------------*
*  |  Only one plain run, no possibility of making statistics:
      ELSE
         DO J =1, KMAX
            IF (GDSTOR(J) .GT. 0.D+00) THEN
               GBSTOR(J)=0.99D+00
            END IF
         END DO
         DO IIX = 1, IX
            NTC = JX (IIX)
            IF (TOTTOT(NTC) .GT. 0.D+00) THEN
               TOTERR(NTC)=0.99D+00
            END IF
         END DO
      END IF
*  |
*  +-------------------------------------------------------------------*
*  +-------------------------------------------------------------------*
*  |  If 1 run already with statistics, no need to create a new output
      IF ( .NOT. LSTATI .OR. IFILE .GT. 1 ) THEN
         RUNTIM='      ***** Sum file *****'
         WRITE(*,'('' Type the output file name: '',$)')
         READ (*,'(A)')FILE
         LQ = LNNBLN(FILE)
* Start_VAX_seq
*        OPEN (UNIT=1,FILE=FILE,STATUS='NEW',FORM='UNFORMATTED')
* End_VAX_seq
* Start_UNIX_seq
         LQS = LLOCCR(FILE,'/') + 1
*  |  +----------------------------------------------------------------*
*  |  |  Be sure to add the default suffix if no one is specified:
         IF ( LQ .GE. LQS ) THEN
            LQD = INDEX (FILE(LQS:LQ),'.')
            IF ( LQD .EQ. LQ-LQS+1 ) THEN
               LQD = 0
               LQ  = LQ - 1
            END IF
            IF ( LQD .LE. 0 ) THEN
               FILE = FILE(1:LQ)//'.trk'
               LQ   = LQ+4
            END IF
         END IF
*  |  |
*  |  +----------------------------------------------------------------*
         OPEN (UNIT=1,FILE=FILE,STATUS='UNKNOWN',FORM='UNFORMATTED')
* End_UNIX_seq
         WRITE (1) RUNTIT,RUNTIM,SNGL(WCTOT),NCTOT,MCTOT,MBATCH
*  |
*  +-------------------------------------------------------------------*
*  |
      ELSE
         FILE = FILE1
         LQ   = LQOLD
      END IF
*  |
*  +-------------------------------------------------------------------*
      KLAST = 0
*  +-------------------------------------------------------------------*
*  |
      DO 6800 JJ=1,IX
         NTC=JX(JJ)
         IF ( .NOT. LSTATI .OR. IFILE .GT. 1 )
     &   WRITE (1)    NTC, TITUTC(NTC), ITUSTC(NTC), IDUSTC(NTC),
     &                NRUSTC(NTC), VUSRTC (NTC), LLNUTC (NTC),
     &                ETCLOW(NTC), ETCHGH(NTC),
     &                NETCBN(NTC), DETCBN(NTC)
**   &                ATCLOW(NTC), ATCHGH(NTC),
**   &                NATCBN(NTC), DATCBN(NTC)
         IF ( LLNUTC (NTC) .AND. ( .NOT. LSTATI .OR. IFILE .GT. 1 ))THEN
            WRITE (1) IGMUTC(NTC), (ENGMAX (IG), IG=1,IGMUTC(NTC)+1)
         END IF
         K0 = KLAST + 1
         K1 = NETCBN (NTC) * NATCBN (NTC) + K0 - 1
         K2 = K1 + IGMUTC (NTC) * NATCBN (NTC)
         KLAST = K2
         IF ( .NOT. LSTATI .OR. IFILE .GT. 1 )
     &   WRITE (1) (SNGL(GDSTOR(J)), J = K0, K2)
6800  CONTINUE
*  |
*  +-------------------------------------------------------------------*
*  Now write the statistical data, first a record to flag them
      IF ( .NOT. LSTATI .OR. IFILE .GT. 1 ) WRITE (1) 'STATISTICS',IX
*  Then the errors concerning the tracklength/collision data
      KLAST = 0
      INDX = INDLST ( FILE, '.' )
      IF ( INDX .GT. 0 ) THEN
* Start_UNIX_seq
         OPEN (UNIT=3,FILE=FILE(1:INDX-1)//'_sum.lis',STATUS='UNKNOWN')
* End_UNIX_seq
* Start_VAX_seq
*        OPEN (UNIT=3,FILE=FILE(1:INDX-1)//'_sum.lis',STATUS='NEW',
*    &         CARRIAGECONTROL='LIST')
* End_VAX_seq
      ELSE
* Start_UNIX_seq
         OPEN (UNIT=3,FILE=FILE(1:LQ)//'_sum.lis',STATUS='UNKNOWN')
* End_UNIX_seq
* Start_VAX_seq
*        OPEN (UNIT=3,FILE=FILE(1:LQ)//'_sum.lis',STATUS='NEW',
*    &         CARRIAGECONTROL='LIST')
* End_VAX_seq
      END IF
      IF ( INDX .GT. 0 ) THEN
* Start_UNIX_seq
         OPEN (UNIT=4,FILE=FILE(1:INDX-1)//'_tab.lis',STATUS='UNKNOWN')
* End_UNIX_seq
* Start_VAX_seq
*        OPEN (UNIT=4,FILE=FILE(1:INDX-1)//'_tab.lis',STATUS='NEW',
*    &         CARRIAGECONTROL='LIST')
* End_VAX_seq
      ELSE
* Start_UNIX_seq
         OPEN (UNIT=4,FILE=FILE(1:LQ)//'_tab.lis',STATUS='UNKNOWN')
* End_UNIX_seq
* Start_VAX_seq
*        OPEN (UNIT=4,FILE=FILE(1:LQ)//'_tab.lis',STATUS='NEW',
*    &         CARRIAGECONTROL='LIST')
* End_VAX_seq
      END IF
      WRITE (3,*)' **** ',FILE(1:LQ),' ****'
      WRITE (3,*)
      WRITE (3,*)
      WRITE (CHCTOT,'(F15.0)') DBLE (NCTOT) + DBLE (MCTOT) * 1.D+09
      WRITE (3,'(A,2X,A14)')'  Total primaries run:',CHCTOT
      WRITE (3,*)
      WRITE (3,*)'   Total weight of the primaries run:',SNGL(WCTOT)
      WRITE (3,*)
*  +-------------------------------------------------------------------*
*  |
      DO 6850 JJ = 1, IX
         NTC    = JX (JJ)
*        WRITE (2,*)'  N: ',NTC,'  ',TITUTC(NTC)
         IF ( .NOT. LSTATI .OR. IFILE .GT. 1 )
     &   WRITE (1) SNGL(TOTTOT(NTC)),SNGL(TOTERR(NTC))
*        WRITE (2,*)'    Tot. response', SNGL(TOTTOT(NTC)),
*    &              '  +/-',SNGL(TOTERR(NTC))*100.0,' %'
         ENERG0 = ETCLOW (NTC)
*  |  +----------------------------------------------------------------*
*  |  |
         IF ( ITUSTC (NTC) .LT. 0 ) THEN
            EFACT  = DETCBN (NTC)
*  |  |
*  |  +----------------------------------------------------------------*
*  |  |
         ELSE
            DELTAE = DBLE (DETCBN (NTC))
         END IF
*  |  |
*  |  +----------------------------------------------------------------*
         ENEOLD = ENERG0
*  |  +----------------------------------------------------------------*
*  |  |
         DO IE = 1, NETCBN (NTC)
            IG = NETCBN (NTC) - IE + 1
*  |  |  +-------------------------------------------------------------*
*  |  |  |
            IF ( ITUSTC (NTC) .LT. 0 ) THEN
               ENERGE = ENERG0 * EFACT**IE
               DELTAE = DBLE (ENERGE - ENEOLD)
*  |  |  |
*  |  |  +-------------------------------------------------------------*
*  |  |  |
            ELSE
               ENERGE = ETCLOW (NTC) + DETCBN (NTC) * FLOAT (IE)
            END IF
*  |  |  |
*  |  |  +-------------------------------------------------------------*
            EPGMAX (IG) = ENERGE
            ENEOLD = ENERGE
         END DO
*  |  |
*  |  +----------------------------------------------------------------*
*  |  +----------------------------------------------------------------*
*  |  |
         DO IG = NETCBN (NTC) + 1, NETCBN (NTC) + IGMUTC (NTC)
            EPGMAX (IG) = ENGMAX (IG-NETCBN(NTC))
         END DO
*  |  |
*  |  +----------------------------------------------------------------*
*  |  +----------------------------------------------------------------*
*  |  |
         IF (LLNUTC (NTC)) THEN
            EPGMAX (NETCBN(NTC)+IGMUTC(NTC)+1) = ENGMAX (IGMUTC(NTC)+1)
            IF (IDUSTC(NTC) .NE. 8)EPGMAX (NETCBN(NTC)+1) = ETCLOW (NTC)
*  |  |
*  |  +----------------------------------------------------------------*
*  |  |
         ELSE
            EPGMAX (NETCBN(NTC)+1) = ETCLOW (NTC)
         END IF
*  |  |
*  |  +----------------------------------------------------------------*
         IF ( .NOT. LSTATI .OR. IFILE .GT. 1 )
     &   WRITE (1) NETCBN(NTC),IGMUTC(NTC),ETCLOW(NTC),(EPGMAX(IG),
     &             IG=1,NETCBN(NTC)+IGMUTC(NTC)+1)
         KK0 = KFUSTC (NTC) + 1
         KK1 = KFUSTC (NTC) + NETCBN (NTC)
         KK2 = KFUSTC (NTC) + NETCBN (NTC) + IGMUTC (NTC)

*        ACCUMP=0.D+00
*        DO IE=1,NETCBN(NTC)
*           IG = NETCBN (NTC) - IE + 1
*           J=KFUSTC(NTC)+IE
*           ACCUMP=ACCUMP+GDSTOR(J)*(EPGMAX(IG)-EPGMAX(IG+1))
*        END DO
*        WRITE(2,*) ACCUMP
*        DO IE=NETCBN(NTC)+1,NETCBN(NTC)+IGMUTC(NTC)
*           IG = IE
*           J=KFUSTC(NTC)+IE
*           ACCUMP=ACCUMP+GDSTOR(J)*(EPGMAX(IG)-EPGMAX(IG+1))
*        END DO
*        WRITE(2,*) ACCUMP

         WRITE (3,*)
         WRITE (3,*) ' Detector n: ', NTC,'  ',TITUTC (NTC)
         WRITE (3,*) '    (Volume:       ',VUSRTC(NTC),' cmc,'
         WRITE (3,*) '     distr. scored:',IDUSTC(NTC),'   ,'
         WRITE (3,*) '     from reg.',NRUSTC(NTC),' )'
         WRITE (3,*)
         WRITE (3,10000)'    Tot. response (p/cmq/pr)',
     &                   SNGL(TOTTOT(NTC)),'  +/-',
     &                   SNGL(TOTERR(NTC))*100.0,' %'
         WRITE (3,10000)'    ( --> Track l. (cm/pr)    ',
     &                   SNGL(TOTTOT(NTC)*VUSRTC(NTC)),'  +/-',
     &                   SNGL(TOTERR(NTC))*100.0,' % )'
         WRITE (3,*)
         WRITE (3,*)
     &         ' **** Different. Fluxes as a function of energy ****',
     &         ' ****      (integrated over solid angle)        ****'
         WRITE (3,*)
         WRITE (3,*)
         WRITE (3,*) ' Energy boundaries (GeV):'
         WRITE (3,*)
         WRITE (3,10002) ( EPGMAX (J), J=1,NETCBN(NTC) )
         WRITE (3,10003) ' Lowest boundary   (GeV):',
     &                     EPGMAX (NETCBN(NTC)+1)
         WRITE (3,*)
         WRITE (3,*) ' Flux (Part/GeV/cmq/pr):'
         WRITE (3,*)
         WRITE (3,10001) ( SNGL (GDSTOR(J)),' +/- ',
     &                     SNGL (GBSTOR(J)*100.D+00),' % ',
     &                     J = KK1, KK0,-1 )
         IF ( LLNUTC (NTC) ) THEN
            IF (IDUSTC(NTC) .NE. 8)
     &         EPGMAX (NETCBN(NTC)+1) = ENGMAX (1)
            WRITE (3,*)
            WRITE (3,*) ' Energy boundaries (GeV):'
            WRITE (3,*)
            WRITE (3,10002) ( EPGMAX (J), J=NETCBN(NTC)+1,
     &                        NETCBN(NTC)+IGMUTC(NTC) )
            WRITE (3,10003) ' Lowest boundary   (GeV):',
     &                        EPGMAX (NETCBN(NTC)+IGMUTC(NTC)+1)
            WRITE (3,*)
            WRITE (3,*) ' Flux (Part/GeV/cmq/pr):'
            WRITE (3,*)
            WRITE (3,10001) ( SNGL (GDSTOR(J)),' +/- ',
     &                        SNGL (GBSTOR(J)*100.D+00),' % ' ,
     &                        J = KK1+1, KK2 )
            IF (IDUSTC(NTC) .NE. 8)
     &         EPGMAX (NETCBN(NTC)+1) = ETCLOW (NTC)
         END IF
         IF ( .NOT. LSTATI .OR. IFILE .GT. 1 ) THEN
            WRITE (1) ( SNGL (GDSTOR(J)), J = KK1, KK0,-1 ),
     &                ( SNGL (GDSTOR(J)), J = KK1+1, KK2 )
            WRITE (1) ( SNGL (GBSTOR(J)), J = KK1, KK0,-1 ),
     &                ( SNGL (GBSTOR(J)), J = KK1+1, KK2 )
         END IF
*****   here to write on "table" file
         IF (NTC .GT. 1) THEN
*****       write 2 empty lines so gnuplot can use it with 'index nnn'
            WRITE (4,*)
            WRITE (4,*)
         END IF
         WRITE (4,*) '# Detector n: ', NTC,'  ',TITUTC (NTC)
         WRITE (4,*) '# N. of energy intervals' ,NETCBN(NTC)+IGMUTC(NTC)
         IF ( LLNUTC (NTC) ) THEN
            IF (IDUSTC(NTC) .NE. 8)
     &         EPGMAX (NETCBN(NTC)+1) = ENGMAX (1)
            KK = KK2 + 1
            DO 6001 J=NETCBN(NTC)+IGMUTC(NTC), NETCBN(NTC)+1,-1
               KK = KK - 1
               WRITE(4,9501) EPGMAX (J+1), EPGMAX (J),
     &              SNGL (GDSTOR(KK)), SNGL (GBSTOR(KK)*100.D+00)
 6001       CONTINUE
            IF (IDUSTC(NTC) .NE. 8)
     &         EPGMAX (NETCBN(NTC)+1) = ETCLOW (NTC)
         END IF
         KK = KK0 - 1
         DO 6002 J=NETCBN(NTC),1,-1
            KK = KK + 1
            WRITE(4,9501) EPGMAX (J+1), EPGMAX (J),
     &              SNGL (GDSTOR(KK)), SNGL (GBSTOR(KK)*100.D+00)
 6002    CONTINUE
******  end write on tab file
         KK0 = KFUSTC (NTC) + NETCBN (NTC) + IGMUTC (NTC) + 1
         KK1 = KFUSTC (NTC) + 2 * NETCBN (NTC) + IGMUTC (NTC)
         KK2 = KFUSTC (NTC) + 2 * ( NETCBN (NTC) + IGMUTC (NTC) )
         WRITE (3,*)
         WRITE (3,*)
         WRITE (3,*)
     &         ' **** Cumulative Fluxes as a function of energy ****',
     &         ' ****      (integrated over solid angle)        ****'
         WRITE (3,*)
         WRITE (3,*)
         WRITE (3,*) ' Energy boundaries (GeV):'
         WRITE (3,*)
         WRITE (3,*)
         WRITE (3,10002) ( EPGMAX (J), J=1,NETCBN(NTC) )
         WRITE (3,10003) ' Lowest boundary   (GeV):',
     &                     EPGMAX (NETCBN(NTC)+1)
         WRITE (3,*)
         WRITE (3,*) ' Cumul. Flux (Part/cmq/pr):'
         WRITE (3,*)
         WRITE (3,10001) ( SNGL (GDSTOR(J)),' +/- ',
     &                     SNGL (GBSTOR(J)*100.D+00),' % ',
     &                     J = KK1, KK0,-1 )
         IF ( LLNUTC (NTC) ) THEN
            IF (IDUSTC(NTC) .NE. 8)
     &         EPGMAX (NETCBN(NTC)+1) = ENGMAX (1)
            WRITE (3,*)
            WRITE (3,*) ' Energy boundaries (GeV):'
            WRITE (3,*)
            WRITE (3,10002) ( EPGMAX (J), J=NETCBN(NTC)+1,
     &                        NETCBN(NTC)+IGMUTC(NTC) )
            WRITE (3,10003) ' Lowest boundary   (GeV):',
     &                        EPGMAX (NETCBN(NTC) + IGMUTC(NTC) + 1)
            WRITE (3,*)
            WRITE (3,*)
            WRITE (3,*) ' Cumul. Flux (Part/cmq/pr):'
            WRITE (3,*)
            WRITE (3,10001) ( SNGL (GDSTOR(J)),' +/- ',
     &                        SNGL (GBSTOR(J)*100.D+00),' % ' ,
     &                        J = KK1+1, KK2 )
            IF (IDUSTC(NTC) .NE. 8)
     &         EPGMAX (NETCBN(NTC)+1) = ETCLOW (NTC)
         END IF
         IF ( .NOT. LSTATI .OR. IFILE .GT. 1 ) THEN
            WRITE (1) ( SNGL (GDSTOR(J)), J = KK1, KK0,-1 ),
     &                ( SNGL (GDSTOR(J)), J = KK1+1, KK2 )
            WRITE (1) ( SNGL (GBSTOR(J)), J = KK1, KK0,-1 ),
     &                ( SNGL (GBSTOR(J)), J = KK1+1, KK2 )
         END IF
         K0 = KLAST + 1
         K1 = NETCBN (NTC) * NATCBN (NTC) + K0 - 1
         K2 = K1 + IGMUTC (NTC) * NATCBN (NTC)
         KLAST = K2
         IF ( .NOT. LSTATI .OR. IFILE .GT. 1 )
     &   WRITE (1) (SNGL(GBSTOR(J)), J = K0, K2)
6850  CONTINUE
*  |
*  +-------------------------------------------------------------------*
      IF ( .NOT. LSTATI .OR. IFILE .GT. 1 ) CLOSE (UNIT=1)
10000 FORMAT ((1X,A,1X,1P,G14.7,A,1X,G14.7,A))
10001 FORMAT ((2(1X,1P,G14.7,A,1X,G14.7,A)))
10002 FORMAT ((5(1X,1P,G14.7)))
10003 FORMAT ((1X,A,1X,1P,G14.7))
*10004 FORMAT ((1X,A,1X,1P,G14.7,1X,G14.7))
 9501 FORMAT (4(1P,1X,G14.7))

      DEALLOCATE(GMSTOR)
      DEALLOCATE(GDSTOR)
      DEALLOCATE(GBSTOR)
      DEALLOCATE(ENGMAX)
      DEALLOCATE(EPGMAX)
      DEALLOCATE(EPGMCK)
      DEALLOCATE(ACCUMG)
      DEALLOCATE(ACUMCU)
      DEALLOCATE(ACCUMH)

      STOP
      END

*
*=== Lnnbln ===========================================================*
*
      INTEGER FUNCTION LNNBLN (CARD)

*
*----------------------------------------------------------------------*
*                                                                      *
*     Copyright (C) 2003-2019:  CERN & INFN                            *
*     All Rights Reserved.                                             *
*                                                                      *
*     Return the length of character variable CARD up to the last non  *
*     blank character                                                  *
*                                                                      *
*     Created on 21 September 1992 by    Alfredo Ferrari & Paola Sala  *
*                                                   Infn - Milan       *
*                                                                      *
*----------------------------------------------------------------------*
*
      CHARACTER CARD*(*)
*
      LENGTH = LEN (CARD)
      DO 100 LQ = LENGTH, 1, -1
         IF ( CARD (LQ:LQ).NE.' ') GO TO 200
  100 CONTINUE
      LQ     = 0
  200 CONTINUE
      LNNBLN = LQ
      RETURN
*=== End of function Lnnbln ===========================================*
      END

*
*=== Lloccr ===========================================================*
*
      INTEGER FUNCTION LLOCCR ( CARD, CWSH )

*
*----------------------------------------------------------------------*
*                                                                      *
*     Copyright (C) 2003-2019:  CERN & INFN                            *
*     All Rights Reserved.                                             *
*                                                                      *
*     Return the length of character variable CARD up to the last      *
*     (included) occurrence of character CWSH (0 if not found)         *
*                                                                      *
*     Created on 21 September 1992 by    Alfredo Ferrari & Paola Sala  *
*                                                   Infn - Milan       *
*                                                                      *
*----------------------------------------------------------------------*
*
      CHARACTER CARD*(*), CWSH*1
*
      LENGTH = LEN (CARD)
      DO 100 LQ = LENGTH, 1, -1
         IF ( CARD (LQ:LQ) .EQ. CWSH ) GO TO 200
  100 CONTINUE
      LQ     = 0
  200 CONTINUE
      LLOCCR = LQ
      RETURN
*=== End of function Lloccr ===========================================*
      END
