Fitting

[HFitting]

The fitting routines in HBOOK are based on the Minuit package [bib-MINUIT]. Minuit is conceived as a tool to find the minimum value of a multi-parameter function and analyze the shape of the function around the minimum. The principal application is foreseen for statistical analysis, working on chisquare or log-likelihood functions, to compute the best-fit parameter values and uncertainties, including correlations between the parameters. It is especially suited to handle difficult problems, including those which may require guidance in order to find the correct solution.

2 In the case of chi minimization, the final fitted parameter values

2 correspond to the minimum of the chi function as defined below:

             2    n                                           2
          chi  = i=1 (((C(I)-F(X(I), A ,A ,. .., A ))/(E(I))))
                                      1  2        k

where the following definitions are used:

n
Number of channels (cells) in the 1-dimensional or 2-dimensional histogram ( HFITH, HFITHN) or number of points of the distribution ( HFITV)
C(I)
Contents of channel (cell) I
E(I)
Error on C(I)
X(I)
Coordinate(s) of centre of channel (cell) I or coordinate vector of point I of the distribution
A ..A
Parameters of the parametric function. 1 k
F
Parametric function to be fitted to the data points.

Remarks:

  1. If no errors are specified by the user, they are taken to be the square root of the bin contents. In this case:
    1. empty channels are ignored by the fitting procedure;
    2. If at least one of C(I)<0, then E(I) is set to 1 for all channels (cells).
    3. If errors are correctly defined via HBARX or HPAKE, then all channels will be taken into account.

      The minimization algorithm requires the calculation of derivatives of the function with respect to the parameters and this is normally done numerically by the fitting routine. If the analytical expression of the derivatives is known, the fit can be speeded up by making use of this information (see option D in the control flag of the various fitting routines).

      One and two-dimensional distributions

         +------------------------------------------------------------------+
         |CALL  HFITH (ID,FUN,CHOPT,NP,*PARAM*,STEP,PMIN,PMAX,SIGPAR*,CHI2*) |
         +------------------------------------------------------------------+
                                        

      Action: Fits a given parametric function to the contents of a 1- or 2-dimensional histogram, and optionally superimposes it to the 1-dimensional histogram when editing.

      Input parameters:
      ID
      Histogram identifier.
      FUN
      Parametric function (to be declared EXTERNAL, see [more info]) Can be defined interactively in the interactive version (see PAW manual)
      CHOPT
      Character variable specifying the desired options.
      'B'
      Some or all parameters are bounded. The arrays STEP, PMIN and PMAX must be specified. By default all parameters vary freely.
      'D'
      The user provides the derivatives of the function analytically with the user routine HDERIV. By default derivatives are computed numerically.
      'E'
      Perform a detailed error analysis using the MINUIT routines MINOS
      'F'
      Force storing of the result of the fit bin by bin with the histogram for an any-dimensional histogram.
      'L'
      Use the logaritmic Likelihood fitting method. By default the chisquared method is used.
      'M'
      Invoke interactive MINUIT.
      'N'
      The results of the fit are not stored bin by bin with the histogram. By default the function is calculated at the centre of each bin in the specified range.
      'Q'
      Quiet mode. No printing.
      'R'
      Fit a Restricted area of a 1-D or 2-D histogram. The limits for the fit are given using in the vector IQUEST is in the communication common /QUEST/IQUEST(100), as follows: IFXLOW = IQUEST(11) specifies the lower limit in X, IFXUP = IQUEST(12) specifies the upper limit in X, IFYLOW = IQUEST(13) specifies the lower limit in Y (2-D), IFYUP = IQUEST(14) specifies the upper limit in Y (2-D).
      'U'
      User function value is taken from /HCFITD/FITPAD(24),FITFUN, see section [more info]. All calculations are performed in double precision.
      'V'
      Verbose mode. Results are printed after each iteration. By default only final results are printed.
      'W'
      Set the event weights equal to one. By default weights are taken according to statistical errors.
      NP
      Number of parameters (NP<=25).
      PARAM
      Array of dimension NP with initial values for the parameters.
      STEP
      Array of dimension NP with initial step sizes for the parameters ('B' option only).
      PMIN
      Array of dimension NP with the lower bounds for the parameters ('B' option only).
      PMAX
      Array of dimension NP with the upper bounds for the parameters ('B' option only).
      Output parameters:
      PARAM
      Array of dimension NP with the final fitted values of the parameters.
      SIGPAR
      Array of dimension NP with the standard deviations on the final fitted values of the parameters.
      CHI2
      Chisquared of the fit.

      Fitting one-dimensional histograms with special functions

        +--------------------------------------------------------------------+
        | CALL  HFITHN (ID,CHFUN,CHOPT,NP,*PAR*,STEP,PMIN,PMAX,SIGPAR*,CHI2*) |
        +--------------------------------------------------------------------+
                                        

      Action: Fits the given special function to the contents of a one-dimensional histogram, and optionally superimposes it to the histogram when editing.

      Input parameters:
      ID
      Histogram identifier.
      CHFUN
      Character variable specifying the desired parametric function. Possible keywords are:
      G
      to fit gaussian PAR(1)*exp(-0.5*((x-PAR(2))/PAR(3))**2). The first parameter PAR(1) corresponds to the normalization, the second parameter PAR(2) corresponds to the mean value, while the third parameter PAR(3) corresponds to the half width of the gaussian.
      E
      to fit exponential exp(PAR(1)+PAR(2)*x)
      Pn
      to fit polynomyal PAR(1)+PAR(2)*x+PAR(3)*x**2......+PAR(n+1)*x**n
      Any combination of these keywords with the 2 operators + or * is allowed, e.g. 'p4+g', a combination of a 4th degree polynomial and a gaussian, which needs eight parameters or p2*g+g, a second degree polynomyal and 2 gaussians, needing 9 parameters. The order of the parameters in PAR must correspond to the order of the basic functions. For example, in the first case above, PAR(1:5) apply to the polynomial of degree 4 and PAR(6:8) to the gaussian while in the second case PAR(1:3) apply to the polynomial of degree 2, PAR(4:6) to the first gaussian and PAR(7:9) to the second gaussian. Blanks are not allowed in the expression.
      CHOPT
      'B'
      Some or all parameters are bounded. The arrays STEP, PMIN and PMAX must be specified. By default all parameters vary freely.
      'D'
      The user is assumed to compute derivatives analytically using the routine HDERIV. By default, derivatives are computed numerically.
      'E'
      Perform a detailed error analysis using the MINOS
      'F'
      Force storing of the result of the fit bin by bin with the histogram.
      'L'
      Use the logaritmic Likelihood fitting method. By default the chisquared method is used.
      'M'
      Invoke interactive MINUIT.
      'N'
      The results of the fit are not stored bin by bin with the histogram. By default the function is calculated at the centre of each bin in the specified range.
      'Q'
      Quiet mode. No printing.
      'R'
      Fit a Restricted area of the 1-D histogram. IFTLOW = IQUEST(11) specifies the lower limit of the minimization domain, IFTUP = IQUEST(12) specifies the upper limit of the minimization domain.
      'V'
      Verbose mode. Results are printed after each iteration. By default only final results are printed.
      'W'
      Set event weights to one. By default weights are taken according to statistical errors.
      NP
      Number of parameters.
      PAR
      Array of dimension NP with initial values for the parameters.
      STEP
      Array of dimension NP with initial step sizes for the parameters ('B' option only).
      PMIN
      Array of dimension NP with the lower bounds for the parameters ('B' option only).
      PMAX
      Array of dimension NP with the upper bounds for the parameters ('B' option only).
      Output parameters:
      PAR
      Array of dimension NP with the final fitted values of the parameters.
      SIGPAR
      Array of dimension NP with the standard deviations on the final fitted values of the parameters.
      CHI2
      Chisquared of the fit.

      Fitting one or multi-demensional arrays

      +---------------------------------------------------------------------------------+
      |CALL  HFITV (N,NDIM,NVAR,X,Y,EY,FUN,CHOPT,NP,*PAR*,STEP,PMIN,PMAX,SIGPAR*, CHI2*) |
      +---------------------------------------------------------------------------------+
                                        

      Action: Fits a given parametric function to a number of value pairs with associated errors.

      Input parameters:
      N
      Number of points to be fitted.
      NDIM
      Declared first dimension of array X.
      NVAR
      Dimension of the distribution.
      X
      Array of dimension N containing the X-coordinates of the points to be fitted.
      Y
      Array of dimension N containing the Y-coordinates of the points to be fitted.
      EY
      Array of dimension N containing the errors on the Y-coordinates of the points to be fitted.
      FUN
      Parametric function (to be declared EXTERNAL)
      CHOPT
      Character variable specifying the desired options.
      'B'
      Some or all parameters are bounded. The arrays STEP, PMIN and PMAX must be specified. By default all parameters vary freely.
      'D'
      The user provides the derivatives of the function analytically with the user routine HDERIV. By default derivatives are computed numerically.
      'E'
      Perform a detailed error analysis using the MINUIT routines
      'L'
      Use the logaritnic Likelihood fitting method. By default the chisquared method is used.
      'M'
      Invoke interactive MINUIT.
      'Q'
      Quiet mode. No printing.
      'U'
      User function value is taken from /HCFITD/FITPAD(24),FITFUN, see section [more info]. All calculations are performed in double precision.
      'V'
      Verbose mode. Results are printed after each iteration. By default only final results are printed.
      'W'
      Set the event weights equal to one. By default weights are taken according to statistical errors. MINOS
      NP
      Number of parameters (NP<=25).
      PAR
      Array of dimension NP with initial values for the parameters.
      STEP
      Array of dimension NP with initial step sizes for the parameters ('B' option only).
      PMIN
      Array of dimension NP with the lower bounds for the parameters ('B' option only).
      PMAX
      Array of dimension NP with the upper bounds for the parameters ('B' option only).
      Output parameters:
      PAR
      Array of dimension NP with the final fitted values of the parameters.
      SIGPAR
      Array of dimension NP with the standard deviations on the final fitted values of the parameters.
      CHI2
      Chisquared of the fit.

      Results of the fit

      [sec:resultfit]

      When the fit is complete, the parameters and the errors on the parameters are returned to the calling program. By default (unless option N is specified), the fitted parameters, the errors on these parameters, their names (see below), the chi-squared and the number of degrees of freedom of the fit are stored in a data structure associated to the histogram ID when routines HFITH and HFITHN are invoked.

      1. For HFITH the value of the fitted function in each histogram channel is also stored in the histogram data structure. The parameters are given the default names: P1, P2,...,Pn.
      2. For HFITHN the type of the function being fitted is stored instead of the fitted values for each channel.

        The information stored in the associated data structure is used during the printing/plotting phase. In particular it is written when the histogram is stored in a file with HROUT.

        Naming the parameters of a fit

                             +------------------------------+
                             | CALL  HFINAM (ID,CHPNAM,NPAR) |
                             +------------------------------+
                                          

        Action: Modify the names of the fitted parameters in the data structure associated with the given histogram.

        Input parameter:
        ID
        Histogram identifier.
        CHPNAM
        CHARACTER array specifying the name(s) to be given to the parameter(s) of the fit (truncated to 8 characters).
        NPAR
        Number of parameters specified.

        Retrieving the fitted parameters

            +----------------------------------------------------------------+
            | CALL  HGFIT (ID,NFPAR*,NPFITS*,FITCHI*,FITPAR*,FITSIG*,FITNAM*) |
            +----------------------------------------------------------------+
                                          

        Action: Retrieve the fitted parameters values from the data structure associated with a given histogram.

        Input parameter:
        ID
        Histogram identifier.
        Output parameters:
        NFPAR
        Number of parameters.
        NPFITS
        Number of data points used in the fit. 2
        FITCHI
        chi of the fit.
        FITPAR
        Array of dimension at least NFPAR, containing the value of the fitted parameters.
        FITSIG
        Array of dimension at least NFPAR, containing the standard deviations of values of the fitted parameters.
        FITNAM
        Character array of dimension at least NFPAR, containing the names of the fitted parameters (truncated to 8 characters). See also HFINAM above.

        The user parametric function

        [sec:userparfun]

        When routines FUN, specified as an argument to those routines, is called during the minimization procedure.

        If the U option is specified with these routines, the user function is calculated in double precision (in the case of HFITHN with the predefined functions (G,E,Pn) double precision is always used). In this case you must reference the common block /HCFITD/, which contains the parameter values in double precision, and store the resulting function value in variable FITFUN, as shown in the example below.

                   Example of user function in double precision
                                          

              FUNCTION UFIT(X)
        *    The dimension of DPAR  ||  MUST be 24!
        *                           VV
              DOUBLE PRECISION DPAR(24),FITFUN
              COMMON/HCFITD/DPAR,FITFUN
              FITFUN = DPAR(1)+DPAR(2)*SQRT(X) +DPAR(3)/X
              UFIT=FITFUN
              END
        

        Even is you do not want to use a double precision function value (i.e. you do not specify the U option), you should still compute the fit function as accurately as possible, using double precision variables in strategic places. This function should also be protected against possible arithmetic exception conditions, like under or overflows and negative arguments for square roots or logarithms.

        User specified derivatives

                                  +--------------------+
                                  |CALL  HDERIV (DERIV) |
                                  +--------------------+
                                          

        Action: User provided subroutine, which calculates the derivatives of the parameters of the function being fitted. This routine must be called from the user function HFITH or HFITN.

        Input parameter:
        DERIV
        Array containing the derivatives of the function being fitted.