*
* $Id$
*
* $Log$
* Revision 1.1  1996/03/01 11:38:39  mclareni
* Initial revision
*
*
#include "paw/pilot.h"
*CMZ :  2.00/01 20/12/92  17.34.35  by  John Allison
*-- Author :    John Allison   08/09/92
      SUBROUTINE PAHSMO (ID, CHOPT, SENSIT, SMOOTH,
     + NPAR, CHI2, NDF, IERSMO)
      CHARACTER CHOPT*(*)
      INTEGER ID, NPAR, NDF, IERSMO
      REAL SENSIT, SMOOTH
* PAW histogram and ntuple smoothing routine.
* Services extended HISTOGRAM/SMOOTH command.
* Includes multiquadric, 353QH and spline smoothing.
*
* Input paramaters:
*   ID     = histogram ID.
*   CHOPT  = option string.  The following charcacters are recognised:
*              0 or 1: overwrite histogram contents.
*                   2: do not overwrite histogram contents.
*                   M: multiquadric smoothing.
*                   N: do not plot result of fit.
*                   Q: 353QH smoothing.
*                   S: spline smoothing.
*                   V: verbose (default for all except 1-D histogram).
*                   F: write Fortran77 function to HQUADF.DAT (multiquadric
*                        only)
*
*                 multiquadric       or           spline
*                 ------------                    ------
*   SENSIT =  sensitivity parameter  or   scale factor for no. of knots
*   SMOOTH =  smoothing parameter    or   additive constant for degree
*
* Output parameters.
*   NPAR   = no. of parameters.
*   CHI2   = chi-squared of fit.
*   NDF    = no. of degrees of freedom of fit.
*   IERSMO = 0 if all's OK.

#include "hbook/hcbook.inc"
#include "hbook/hcbits.inc"
#include "hbook/hcunit.inc"

      CHARACTER*80 CHOPTQ
      LOGICAL VPRINT, FFILE
      LOGICAL HEXIST
      INTEGER ISTORE, IOPT
      INTEGER NKNOTS, IKX, NX, NY
      REAL PROB

* Clear output parameters.
      NPAR   = 0
      CHI2   = 0.
      NDF    = 0
      IERSMO = 0

* Decode option string.
      CHOPTQ = CHOPT
      LCHOPQ = LENOCC (CHOPTQ)
      ISTORE = 2
      IOPT   = 1
      FFILE = .FALSE.
      IF (INDEX (CHOPT, '0') .NE. 0) ISTORE = 0
      IF (INDEX (CHOPT, '1') .NE. 0) ISTORE = 1
      IF (INDEX (CHOPT, '2') .NE. 0) ISTORE = 2
      IF (INDEX (CHOPT, 'M') .NE. 0) IOPT = 1
      IF (INDEX (CHOPT, 'Q') .NE. 0) IOPT = 2
      IF (INDEX (CHOPT, 'S') .NE. 0) IOPT = 3
      IF (INDEX (CHOPT, 'F') .NE. 0) FFILE = .TRUE.
      IF (INDEX (CHOPT, 'V') .NE. 0 .OR. I1 .EQ. 0) THEN
         VPRINT = .TRUE.
         CHOPTQ = CHOPTQ (1: LCHOPQ) // 'V'
      ELSE
         VPRINT = .FALSE.
      END IF

      IF (IOPT .EQ. 1) THEN
* Multiquadric smoothing.
         IF (FFILE) THEN
            CALL PALUNF(60,3,LUNP)
            CALL KUOPEN(LUNP,'HQUADF.DAT','UNKNOWN',ISTAT)
            CALL HSETPR('PLUN',REAL(LUNP))
         END IF
         CALL HQUAD (ID, CHOPTQ, 0, SENSIT, SMOOTH, NPAR, CHI2, NDF,
     +   FMIN, FMAX, IERSMO)
         IF (FFILE) THEN
            CALL HSETPR('PLUN',0.)
            CALL PACLOS(LUNP)
         END IF
* HQUAD's chi-squared is a genuine one, NOT chi-squared per degree of freedom.
         IF (IERSMO .NE. 0) THEN
            WRITE (LOUT, 10000)
            IF (I4 .EQ. 0) WRITE (LOUT, 10100)
         ELSE
            IF (VPRINT) THEN
               WRITE (LOUT, 10200) NPAR, CHI2, NDF
               WRITE (LOUT, 10300) FMIN, FMAX
            END IF
            IF (NDF .GT. 0) THEN
               CPROB = PROB (CHI2, NDF)
               IF (CPROB .LT. 1.E-10) WRITE (LOUT, 10400)
            END IF
         END IF
      END IF
      IF (IOPT .EQ. 2) THEN
* 353QH smoothing.
         IF (I1 .NE. 0) THEN
            CALL HSMOOF (ID, ISTORE, CHI2)
* HSMOOF's chi-squared is a genuine one, NOT chi-squared per degree of freedom,
*   but no. of parameters is a mystery.
            IF (VPRINT) WRITE (LOUT, 10500) CHI2, IQ (LCID + KNCX)
         ELSE
            WRITE (LOUT, 10600)
            IERSMO = 1
            GO TO 10
         END IF
      END IF
      IF (IOPT .EQ. 3) THEN
* Spline smoothing (cubic spline only, no. of knots variable; for more
*   freedom use SPLINE command).
         NKNOTS = 10. * SENSIT
         IKX = SMOOTH + 2
         IF (I1 .NE. 0) THEN
            NX = IQ (LCID + KNCX)
            NPAR = NKNOTS - IKX - 1
            NDF = NX - NPAR
            CALL HSPLI1 (ID, ISTORE, NKNOTS, IKX, CHI2)
* HSPLI1's chi-squared is NOT a genuine one, it is chi-squared per degree of
*   freedom.
            IF (VPRINT) THEN
               WRITE (LOUT, 10700) NKNOTS, IKX, NPAR, CHI2 * NDF, NDF
               WRITE (LOUT, 10800)
            END IF
         ELSE IF (I230 .NE. 0) THEN
            IF (ISTORE .EQ. 0 .OR. ISTORE .EQ. 1) THEN
               NX = IQ (LCID + KNCX)
               NY = IQ (LCID + KNCY)
               NPAR = (NKNOTS - IKX - 1) ** 2
               NDF = NX * NY - NPAR
               CHI2 = 0.
               CALL HSPLI2 (ID, NKNOTS, NKNOTS, IKX, IKX)
               IF (VPRINT) WRITE (LOUT, 10900) NKNOTS, IKX, NPAR
            ELSE
               WRITE (LOUT, 11000)
               IERSMO = 1
               GO TO 10
            ENDIF
            WRITE (LOUT, 10800)
         ELSE
            WRITE (LOUT, 11100)
            IERSMO = 1
            GO TO 10
         ENDIF
      END IF

* Redraw histogram with fitted function.
      IF (HEXIST (ID).AND.INDEX(CHOPT,'N').EQ.0)THEN
         CALL HPLOT (ID, ' ', ' ', 0)
      ENDIF

   10 CONTINUE

10000 FORMAT (1X, 'Unable to do multiquadric smoothing.')
10100 FORMAT (3X, 'Try another method - OPTION = Q or S.')
10200 FORMAT (1X, 'Multiquadric smoothing with', I4,  ' parameters.'/
     +3X,'  Chi-squared', G12.5, ' for', I7, ' degrees of freedom.')
10300 FORMAT (3X, 'Min/max event density:', 2G12.5)
10400 FORMAT (3X, 'Chi-squared probability is very low.'/
     +3X, 'Try a larger sensitivity parameter (e.g., SMOOTH id ! 1.5)'/
     +3X, 'Are the data genuinely random?  Are the bin contents ',
     +'independent?'/
     +3X, 'Are the errors correct?  (This method assumes data are'/
     +3X, 'randomly drawn from their parent probability distribution.)')
10500 FORMAT (1X, '353QH smoothing.' / '  Chi-squared',
     +G12.5, ' for', I7, ' bins.')
10600 FORMAT (1X, '353QH smoothing not available for 2-D or ntuples.')
10700 FORMAT (1X, 'Spline smoothing with', I3, ' knots, degree', I2,
     +' and', I4, ' parameters.'/
     +3X, 'Chi-squared', G12.5, ' for', I7, ' degrees of freedom.')
10800 FORMAT (3X,
     +'(Note: the SPLINE command gives you more flexibility.)')
10900 FORMAT (1X, 'Spline smoothing with', I3, ' knots, degree', I2,
     +' and', I4, ' parameters.')
11000 FORMAT (3X, '2-D spline smoothing routine HSPLI2 always ',
     +'overwrites.'/
     +3X, 'Specify option 1S to confirm.')
11100 FORMAT (1X, 'Spline smoothing not available for ntuples.')
*
      END
