*
* $Id: pahsmo.F,v 1.1.1.1 1996/03/01 11:38:40 mclareni Exp $
*
* $Log: pahsmo.F,v $
* Revision 1.1.1.1 1996/03/01 11:38:40 mclareni
* Paw
*
*
#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