*$ CREATE USBSUW.FOR
*COPY USBSUW
*
*=== Usbsuw ===========================================================*
*
      PROGRAM USBSUW

      CHARACTER FILE*232,RUNTIT*80,RUNTIM*32,CHSTAT*10
*
*----------------------------------------------------------------------*
*                                                                      *
*     Copyright (C) 1989-2017      by  Alberto Fasso`, Alfredo Ferrari *
*     All Rights Reserved.                    & Paola Sala             *
*                                                                      *
*                                                                      *
*     USer Binning scoring SUm and Weighting:                          *
*                                                                      *
*     Created on 11 September 1989 by    Alfredo Ferrari & Paola Sala  *
*                                                   Infn - Milan       *
*                                                                      *
*     Last change on  03-Jan-17    by    Alfredo Ferrari               *
*                                                                      *
*----------------------------------------------------------------------*
*
*=== usrbin ===========================================================*
*
*----------------------------------------------------------------------*
*     Module USRBIN:                                                   *
*     A. Ferrari & A. Fasso': user defined binnings                    *
*                                                                      *
*          Last change A. Ferrari 12-feb-1997                          *
*                                                                      *
*                                                                      *
*     Up to MXUSBN user defined binnings are allowed                   *
*            lusbin = logical flag, .true. if at least 1 user defined  *
*                     binning is used                                  *
*            lusevt = logical flag, .true. if at least 1 user defined  *
*                     binning event by event is used                   *
*            lustkb = logical flag, .true. if at least 1 user defined  *
*                     track-length binning is used. Track length bin-  *
*                     ning are recognized as binnings where accurate   *
*                     deposition along the track is requested but      *
*                     for generalized particles other than 208/211     *
*            nusrbn = number of user defined binnings used             *
*            itusbn = type of binning: 0 = cartesian, .ne. 0 = RZ      *
*            idusbn = distribution to be scored: the usual values are  *
*                     allowed.                                         *
*            titusb = binning name                                     *
*            ipusbn = logical unit to print the results on: formatted  *
*                     if > 0, unformatted if < 0                       *
*            kbusbn = initial location in blank common of the consi-   *
*                     dered binning (for index 1, real*8 address for   *
*                     lsngbn = .false., real*4 address for lsngbn      *
*                            = .true.)                                 *
*            nxbin  = number of x (r for RZ) intervals                 *
*            nybin  = number of y (1 for RZ) intervals                 *
*            nzbin  = number of z intervals                            *
*         xlow/high = minimum and maximum x (r   for R-Phi-Z)          *
*         ylow/high = minimum and maximum y (phi for R-Phi-Z)          *
*         zlow/high = minimum and maximum z                            *
*            dxusbn = x (r) bin width                                  *
*            dyusbn = y bin width                                      *
*            dzusbn = z bin width                                      *
*            tcusbn = time cut-off (seconds) for this binning          *
*            bkusbn = 1st Birk's law parameter for this binning        *
*                     (meaningful only for energy scoring)             *
*            b2usbn = 2nd Birk's law parameter for this binning        *
*                     (meaningful only for energy scoring)             *
*            xaxusb = x-axis offset for R-Z binning (not possible for  *
*                     R-Phi-Z), just for back-compatibility, now it    *
*                     should be done through rototraslations           *
*            yaxusb = y-axis offset for R-Z binning (not possible for  *
*                     R-Phi-Z), just for back-compatibility, now it    *
*                     should be done through rototraslations           *
*            levtbn = logical flag for binning to be printed at the    *
*                     end of each event only                           *
*            lntzer = logical flag for printing only non zero cells    *
*            ltrkbn = logical flag for flagging track-length binnings  *
*            lsngbn = logical flag for single precision storing of     *
*                     data                                             *
*            krtnbn = flag of a possible pre-defined rotation (if >0)  *
*                     to be applied to the binnings before scoring     *
*                                                                      *
*----------------------------------------------------------------------*
*
      PARAMETER ( MXUSBN = 400 )
      LOGICAL LUSBIN, LEVTBN, LNTZER, LUSEVT, LUSTKB, LTRKBN, LSNGBN
      CHARACTER*10 TITUSB
      COMMON /UUSRBN/ XLOW  (MXUSBN), XHIGH (MXUSBN), YLOW  (MXUSBN),
     &                YHIGH (MXUSBN), ZLOW  (MXUSBN), ZHIGH (MXUSBN),
     &                DXUSBN(MXUSBN), DYUSBN(MXUSBN), DZUSBN(MXUSBN),
     &                TCUSBN(MXUSBN), BKUSBN(MXUSBN), B2USBN(MXUSBN),
     &                XAXUSB(MXUSBN), YAXUSB(MXUSBN),
     &                NXBIN (MXUSBN), NYBIN (MXUSBN), NZBIN (MXUSBN),
     &                ITUSBN(MXUSBN), IDUSBN(MXUSBN), KBUSBN(MXUSBN),
     &                IPUSBN(MXUSBN), KRTNBN(MXUSBN), LEVTBN(MXUSBN),
     &                LNTZER(MXUSBN), LSNGBN(MXUSBN), LTRKBN(MXUSBN),
     &                NUSRBN, LUSBIN, LUSEVT, LUSTKB
      COMMON /USRCH/  TITUSB(MXUSBN)
      PARAMETER ( MXDUMM = 60000000 )
      PARAMETER ( PIPIPI = 3.141592653589793D+00 )
      DIMENSION GMSTOR (MXDUMM), GESTOR (MXDUMM),
     &          GDSTOR (MXDUMM), GBSTOR (MXDUMM),
     &          JB(MXUSBN)
      DOUBLE PRECISION GDSTOR, GBSTOR, WCTOT , GBHELP
*
      DO 5555 I = 1, MXDUMM
         GDSTOR (I) = 0.D+00
         GBSTOR (I) = 0.D+00
         GMSTOR (I) = 0.E+00
         GESTOR (I) = 0.E+00
 5555 CONTINUE
      NCTOT  = 0
      MCTOT  = 0
      WCTOT  = 0.D+00
      MBATCH = 0
 1111 CONTINUE
      WRITE(*,'('' Type the input file: '',$)')
      READ (*,'(A)')FILE
      LQ     = LNNBLN(FILE)
      IF (LQ .LE. 0) GO TO 3000
      OPEN (UNIT=1,FILE=FILE,STATUS='OLD',FORM='UNFORMATTED',ERR=1111)
      IRECRD = 0
      READ (1,ERR=100,END=1111) 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
      KLAST  = 0
      KMAX   = 0
*  +-------------------------------------------------------------------*
*  |  No statistics inside (or old fashioned file):
      IF ( NBATCH .LE. 1 ) THEN
         DO 700 IB = 1, 10000
            NB = IB
            READ (1,ERR=400,END=1000) MB, TITUSB(NB), ITUSBN(NB),
     &                      IDUSBN(NB),
     &                      XLOW(NB),XHIGH(NB),NXBIN(NB),
     &                      DXUSBN(NB),YLOW(NB),
     &                      YHIGH(NB),NYBIN(NB),DYUSBN(NB),
     &                      ZLOW(NB),ZHIGH(NB),NZBIN(NB),
     &                      DZUSBN(NB),LNTZER(NB),BKUSBN(NB),
     &                      B2USBN(NB),TCUSBN(NB)
            IRECRD = IRECRD + 1
            GO TO 500
 400        CONTINUE
*           BACKSPACE 1
            REWIND 1
            DO IR = 1, IRECRD
               READ (1)
            END DO
            READ (1,ERR=450,END=1000) MB, TITUSB(NB), ITUSBN(NB),
     &                        IDUSBN(NB),
     &                        XLOW(NB),XHIGH(NB),NXBIN(NB),
     &                        DXUSBN(NB),YLOW(NB),
     &                        YHIGH(NB),NYBIN(NB),DYUSBN(NB),
     &                        ZLOW(NB),ZHIGH(NB),NZBIN(NB),
     &                        DZUSBN(NB)
            IRECRD = IRECRD + 1
            LNTZER(NB)=.FALSE.
            TCUSBN(NB)=1.E+38
            BKUSBN(NB)=0.E+00
            B2USBN(NB)=0.E+00
            GO TO 500
 450        CONTINUE
*           BACKSPACE 1
            REWIND 1
            DO IR = 1, IRECRD
               READ (1)
            END DO
            READ (1,ERR=460) CHSTAT,IJX
            IRECRD = IRECRD + 1
            IF ( CHSTAT .EQ. 'STATISTICS' ) GO TO 1000
 460        CONTINUE
            STOP 'UNKNOWN RECORD TYPE'
 500        CONTINUE
            JB(IB)=NB
*           JB(IB)=MB
            K0 = KLAST + 1
            K1 = NXBIN(NB) * NYBIN(NB) * NZBIN(NB) + K0 -1
            KLAST = K1
            KMAX  = MAX ( KMAX, KLAST )
            READ (1) (GMSTOR(J), J = K0, K1)
            IRECRD = IRECRD + 1
            DO J = K0, K1
               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
 700     CONTINUE
 1000    CONTINUE
         IB = IB - 1
*  |
*  +-------------------------------------------------------------------*
*  |  Treat the cases with statistics inside (when NBATCH.GT.1)
      ELSE
         DO IB = 1, 1000
            NB = IB
            READ (1,ERR=1600) MB, TITUSB(NB), ITUSBN(NB),
     &                      IDUSBN(NB),
     &                      XLOW(NB),XHIGH(NB),NXBIN(NB),
     &                      DXUSBN(NB),YLOW(NB),
     &                      YHIGH(NB),NYBIN(NB),DYUSBN(NB),
     &                      ZLOW(NB),ZHIGH(NB),NZBIN(NB),
     &                      DZUSBN(NB),LNTZER(NB),BKUSBN(NB),
     &                      B2USBN(NB),TCUSBN(NB)
            IRECRD = IRECRD + 1
            JB(IB)=NB
*           JB(IB)=MB
            K0 = KLAST + 1
            K1 = NXBIN(NB) * NYBIN(NB) * NZBIN(NB) + K0 -1
            KLAST = K1
            KMAX  = MAX ( KMAX, KLAST )
            READ (1) (GMSTOR(J), J = K0, K1)
            IRECRD = IRECRD + 1
         END DO
1600     CONTINUE
*        BACKSPACE 1
         REWIND 1
         DO IR = 1, IRECRD
            READ (1)
         END DO
         READ (1,ERR=9999) CHSTAT,IJX
         IRECRD = IRECRD + 1
         IF ( CHSTAT .NE. 'STATISTICS' ) GO TO 9999
         IB   = IB - 1
         KLAST= 0
         KMAX = 0
         DO NB = 1, IB
            K0 = KLAST + 1
            K1 = NXBIN(NB) * NYBIN(NB) * NZBIN(NB) + K0 -1
            KLAST = K1
            KMAX  = MAX ( KMAX, KLAST )
            READ (1) (GESTOR(J), J = K0, K1)
            IRECRD = IRECRD + 1
            DO J = K0, K1
               GESTOR (J) = (GESTOR(J)*GMSTOR(J))**2
               GESTOR (J) = GESTOR(J) * DBLE(NBATCH-1) + GMSTOR(J)**2
               GDSTOR (J) = GDSTOR (J) + DBLE (GMSTOR(J))
     &                    * ( DBLE (NCASE) + DBLE (MCASE) * 1.D+09 )
               GBSTOR (J) = GBSTOR (J) + DBLE (GESTOR(J))
     &                    * ( DBLE (NCASE) + DBLE (MCASE) * 1.D+09 )
*  |  |  |  ???? Useless !!!
*              GMSTOR (J) = GMSTOR(J)
*    &                    * ( DBLE (NCASE) + DBLE (MCASE) * 1.D+09 )
*              GESTOR (J) = GESTOR(J)
*    &                    * ( DBLE (NCASE) + DBLE (MCASE) * 1.D+09 )
            END DO
         END DO
      END IF
*  |
*  +-------------------------------------------------------------------*
      CLOSE (UNIT=1)
      GO TO 1111
 3000 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
*  +-------------------------------------------------------------------*
*  |  Real sum of Mbatch> 1 runs:
      IF ( MBATCH - 1 .GT. 0 ) THEN
         DO J = 1, KMAX
            GBHELP    = SQRT (GBSTOR(J))
            GBSTOR(J) = ( ( GBHELP - GDSTOR (J) )
     &                  * ( GBHELP + GDSTOR (J) ) )
     &               / DBLE (MBATCH-1)
            IF ( GBSTOR (J) .LT. 1.D-12 * GDSTOR(J)**2 )
     &           GBSTOR (J) = 0.D+00
            IF ( ABS (GDSTOR (J)) .GT. 0.D+00 )
     &          GBSTOR (J) = SQRT ( GBSTOR (J) ) / ABS (GDSTOR (J))
*    &          GBSTOR (J) = MIN ( SQRT(GBSTOR(J))/GDSTOR(J),0.99D+00 )
         END DO
*  |
*  +-------------------------------------------------------------------*
*  |  Only one run available, no statistics is possible
      ELSE
         DO J = 1, KMAX
            IF ( GDSTOR (J) .GT. 0.D+00 ) THEN
               GBSTOR (J) = 0.99D+00
            END IF
         END DO
      END IF
*  |
*  +-------------------------------------------------------------------*
      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)//'.bnn'
            LQ   = LQ+4
         END IF
      END IF
*  |
*  +-------------------------------------------------------------------*
      OPEN (UNIT=1,FILE=FILE,STATUS='UNKNOWN',FORM='UNFORMATTED')
* End_UNIX_seq
      WRITE (*,*) ' Total Batches ', MBATCH
      WRITE (1) RUNTIT,RUNTIM,SNGL(WCTOT),NCTOT,MCTOT,MBATCH
      KLAST = 0
      DO JJ = 1, IB
         NB = JB (JJ)
         WRITE (1) NB, TITUSB(NB), ITUSBN(NB),
     &             IDUSBN(NB),
     &             XLOW(NB),XHIGH(NB),NXBIN(NB),
     &             DXUSBN(NB),YLOW(NB),
     &             YHIGH(NB),NYBIN(NB),DYUSBN(NB),
     &             ZLOW(NB),ZHIGH(NB),NZBIN(NB),
     &             DZUSBN(NB),LNTZER(NB),BKUSBN(NB),
     &             B2USBN(NB),TCUSBN(NB)  
         K0 = KLAST + 1
         K1 = NXBIN(NB) * NYBIN(NB) * NZBIN(NB) + K0 -1
         KLAST = K1
         WRITE (1) ( SNGL (GDSTOR(J)), J = K0, K1)
      END DO
*  Now write the statistical data, first a record to flag them
      WRITE (1) 'STATISTICS',IB
      KLAST = 0
      DO JJ = 1, IB
         NB = JB(JJ)
         K0 = KLAST + 1
         K1 = NXBIN(NB) * NYBIN(NB) * NZBIN(NB) + K0 -1
         KLAST = K1
         WRITE (1) ( SNGL (GBSTOR(J)), J = K0, K1 )
      END DO
      CLOSE (UNIT=1)
      STOP
 9999 CONTINUE
      STOP 'UNKNOWN RECORD TYPE'
      END

*$ CREATE LNNBLN.FOR
*COPY LNNBLN
*
*=== Lnnbln ===========================================================*
*
      INTEGER FUNCTION LNNBLN (CARD)

*
*----------------------------------------------------------------------*
*                                                                      *
*     Copyright (C) 1992-2005      by    Alfredo Ferrari & Paola Sala  *
*     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       *
*                                                                      *
*     Last change on 13-jun-01     by    Alfredo Ferrari               *
*                                                                      *
*                                                                      *
*----------------------------------------------------------------------*
*
      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

*$ CREATE LLOCCR.FOR
*COPY LLOCCR
*
*=== Lloccr ===========================================================*
*
      INTEGER FUNCTION LLOCCR ( CARD, CWSH )

*
*----------------------------------------------------------------------*
*                                                                      *
*     Copyright (C) 1992-2016      by    Alfredo Ferrari & Paola Sala  *
*     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       *
*                                                                      *
*     Last change on  05-Jun-16    by    Alfredo Ferrari               *
*                                                                      *
*                                                                      *
*----------------------------------------------------------------------*
*
      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

