HBOOK routines

The following three routines are intended to help with problems of the above type. Subroutine HMCMLL uses MINUIT to perform the log-likelihood maximisation and return a set of fractions. Functions HMCINI and HMCLNL are provided for those who wish to perform the fit themselves.

N.B. Real parameters and functions are all REAL*8 (Double precision on most machines).

CALL HMCMLL (IDD,IDM,IDW,NSRC,CHOPT,IFIX,FRAC,FLIM,START, STEP,UP,PAR*,DPAR*)

Action: Fits the given Monte Carlo distributions to the data distribution, using a binned maximum likelihood fit which includes the effect of both data and Monte Carlo statistics, and allows weights to be provided for each Monte Carlo distribution. The data, Monte Carlo and weight distributions must be presented in identically binned 1 dimensional histograms. Weight distributions which just change the shape of the Monte Carlo spectra (i.e. not efficiency distributions) must be normalised so that

$\sum$_iw_jia_ji= ∑_ia_ji= N_j.

The best estimate of the fraction of each Monte Carlo distribution present in the data distribution is returned, with an error estimate where required.

Input parameters:

IDD

Data histogram identifier.

IDM

Array of dimension NSRC containing Monte Carlo histogram identifiers.

IDW

Array of dimension NSRC containing weight histogram identifiers. ('W' option only).

NSRC

Number of Monte Carlo sources. Must be greater than 1.

CHOPT

'F'

Fix one or more of the fractions. Default is for all fractions to vary freely.

'L'

Set limits on the fractions as given in FLIM. Default is no limits.

'W'

Use the weight histograms provided. For non existent weight histograms, and if the 'W' option is not requested, a dummy weight histogram in which all entries are 1 is booked.

'S'

Scan the likelihood function with respect to each fit parameter, before and after the fit. If the 'N' option is specified, the function will only be scanned once for each parameter.

'N'

Do not perform the fit.

'P'

Use the parameter start points and initial step sizes provided in START and STEP. If the 'P' option is not specified then the start point for each free parameter is

$\{1\; -\; \Sigma$_fNSRC-N_f}

, where $\Sigma$_f is the sum of the fixed fractions, and $N$_f is the number of fixed fractions; and the initial step size is 0.01.

'E'

Perform a detailed error analysis using the MINUIT routines HESSE and MINOS.

IFIX

Array of dimension NSRC containing '1' if a parameter is to be fixed in the fit, '0' otherwise. ('F' option only).

FRAC

Array of dimension NSRC with the values at which parameters are to be fixed. ('F' option only).

FLIM

Array of dimension (2,NSRC) with the lower, then upper limits on the parameters. ('L' option only.)

START

Array of dimension NSRC with the start values for each parameter. ('P' option only - if 'P' option is chosen then default start values are used if values in START are negative).

STEP

Array of dimension NSRC with initial step values for each parameter. ('P' option only - if 'P' option is chosen then default step values are used if values in STEP are negative).

UP

UP value for the error estimate ('E' option only). Default 0.5 (if user supplies negative or zero value for UP when 'E' option is chosen). See the Minuit manual for definition of UP.

Output parameters:

PAR

Array of dimension NSRC with the final fitted values of the parameters.

DPAR

Array of dimension NSRC with the errors on the final fitted values of the parameters.

CALL HMCINI (IDDATA,IDMC,IDWT,NSRC,CHOPT,IERR)

Action: Initialisation routine for function HMCLNL, needs to be called each time a new set of histograms is introduced (generally once at the beginning of each fit). Performs some error checking and sets up a dummy weight histogram if necessary.

Input parameters:

IDDATA

Data histogram identifier.

IFIXMC

Array of dimension NSRC containing '1' if a

IDMC

Array of dimension NSRC containing Monte Carlo histogram identifiers.

IDWT

Array of dimension NSRC containing weight histogram identifiers. ('W' option only).

NSRC

Number of Monte Carlo sources.

CHOPT

'W'

Use the weight histograms provided. For non existent weight histograms, and if the W option is not requested, a dummy weight histogram in which all entries are 1 is booked.

IERR

Error flag - set to 1 if the parameters sent to HMCINI were not usable (e.g. incompatibility between data and MC histograms, number of MC sources less than 2), 0 otherwise.

HMCLNLVARIAB = HMCLNL(FRAC)

Action: HMCLNL is a double precision function giving the log likelihood (including effect of both data and Monte Carlo statistics) that the data distribution arose from a distribution given by combining the Monte Carlo distributions, weighted by the weights provided, using the fractions given in FRAC. HMCINI must be called before this function may be used.

Input parameters:

FRAC

Array of dimension NSRC containing the fraction of each Monte Carlo distribution you wish to assume is in the data distribution, in order to calculate the log likelihood.

Example of use of HMCMLL

       PROGRAM MCSTAT
 
       PARAMETER (NPAWC=10000)
       COMMON/PAWC/H(NPAWC)
 
       INTEGER NMCSRC
       PARAMETER (NMCSRC=2)
 
* declarations for HMCMLL
       INTEGER IDDATA,IDMC(NMCSRC),IDWT(NMCSRC),IFIXMC(NMCSRC)
       DOUBLE PRECISION FLIM(2,NMCSRC),FSTART(NMCSRC),
     + FSTEP(NMCSRC),FUP,FRAC(NMCSRC),ANS(NMCSRC),
     + DANS(NMCSRC)
 
       DATA IDDATA,IDMC,IDWT/10,1301,1302,1351,1352/
 
       CALL HLIMIT(NPAWC)
 
* Set 'UP' to correspond to 70% confidence level
* for a 2 parameter fit.
       FUP=1.2
 
* Read in your data, Monte Carlo and Weight histograms
       CALL HRGET(10,'HMCHIS.PAW',' ')
       DO 20 JSRC=1,NMCSRC
          CALL HRGET(IDMC(JSRC),'HMCHIS.PAW',' ')
          CALL HRGET(IDWT(JSRC),'HMCHIS.PAW',' ')
20     CONTINUE
 
* perform log likelihood maximisation and error analysis,
* using user weights + setting limits on the fractions.
       WRITE(6,1000)
1000   FORMAT(' ** Fit with error analysis - use user weights
     +and limits',/)
       FLIM(1,1)=0.0
       FLIM(2,1)=1.0D0
       FLIM(1,2)=0.0
       FLIM(2,2)=1.0D0
       CALL HMCMLL(IDDATA,IDMC,IDWT,NMCSRC,'EWL',IFIXMC,
     + FRAC,FLIM,FSTART,FSTEP,FUP,ANS,DANS)
 
       STOP
       END

{\scriptsize
 ** Fit with error analysis - use user weightsand limits
 
  MINUIT RELEASE 93.08  INITIALIZED.   DIMENSIONS 100/ 50  EPSMAC=  0.56E-16
 MCMLL: You have  2 free fractions and  0 fixed
 
 PARAMETER DEFINITIONS:
    NO.   NAME         VALUE      STEP SIZE      LIMITS
     1 'P1        '   0.50000      0.10000E-01    0.00000E+00   1.0000
     2 'P2        '   0.50000      0.10000E-01    0.00000E+00   1.0000
 **********
 **    1 **CALL FCN   1.000
 **********
 **********
 **    2 **SET ERR  0.5000
 **********
 **********
 **    3 **MIGRAD
 **********
 
 FIRST CALL TO USER FUNCTION AT NEW START POINT, WITH IFLAG=4.
 START MIGRAD MINIMIZATION.  STRATEGY 1.  CONVERGENCE WHEN EDM .LT. 0.50E-04
 
 FCN=  -8317.393     FROM MIGRAD    STATUS=INITIATE       8 CALLS       10 TOTAL
                     EDM= unknown      STRATEGY= 1      NO ERROR MATRIX
 
  EXT PARAMETER               CURRENT GUESS       STEP         FIRST
  NO.   NAME        VALUE          ERROR          SIZE      DERIVATIVE
   1      P1       0.50000       0.10000E-01   0.20001E-01   -8.1370
   2      P2       0.50000       0.10000E-01   0.20001E-01    8.7598
                               ERR DEF=  0.50
 
 MIGRAD MINIMIZATION HAS CONVERGED.
 
 MIGRAD WILL VERIFY CONVERGENCE AND ERROR MATRIX.
 COVARIANCE MATRIX CALCULATED SUCCESSFULLY
 
 FCN=  -8319.763     FROM MIGRAD    STATUS=CONVERGED     45 CALLS       47 TOTAL
                     EDM=  0.66E-08    STRATEGY= 1      ERROR MATRIX ACCURATE
 
  EXT PARAMETER                                   STEP         FIRST
  NO.   NAME        VALUE          ERROR          SIZE      DERIVATIVE
   1      P1       0.64665       0.75183E-01   0.25797E-02  -0.77581E-03
   2      P2       0.35335       0.70935E-01   0.24335E-02  -0.10448E-02
                               ERR DEF=  0.50
 
 EXTERNAL ERROR MATRIX.    NDIM=  50    NPAR=  2    ERR DEF=  0.50
  0.570E-02-0.460E-02
 -0.460E-02 0.507E-02
 
 PARAMETER  CORRELATION COEFFICIENTS
       NO.  GLOBAL     1     2
        1  0.85488  1.000-0.855
        2  0.85488 -0.855 1.000
 **********
 **    4 **SET ERR   1.200
 **********
 MCMLL: SET UP VALUE TO  1.20
 MCMLL: FOR  2 FREE PARAMETERS
 **********
 **    5 **MINOS
 **********
 
 MINUIT TASK: *** HBOOK New ll maximisation
 
 FCN=  -8319.763     FROM MINOS     STATUS=SUCCESSFUL    72 CALLS      119 TOTAL
                     EDM=  0.66E-08    STRATEGY= 1      ERROR MATRIX ACCURATE
 
  EXT PARAMETER                  PARABOLIC         MINOS ERRORS
  NO.   NAME        VALUE          ERROR      NEGATIVE      POSITIVE
   1      P1       0.64665       0.11579      -0.11144       0.12598
   2      P2       0.35335       0.10932      -0.11560       0.10891
                               ERR DEF=   1.2
}

Figure: Monte Carlo distributions (left) and data distribution (right).