The MINUIT-based mSR fitting programme MINFIT

Ivan Reid, PSI mSR Facility

Introduction

These notes detail the current state of the MINUIT-basedfitting programme MINFIT. Look here to have more information on the available versions.

Functions vs. Signals

In my terminology, I make a distinction between a function and a signal. A function is a complete description of a set of mSR data, i.e.

Y(t) = b + Ne^-t/t ( 1 + f(t) ),

where, of course, N is the normalisation to time t₀, b is the time-independent background, and t is the muon lifetime (2.197037 ms). However, in general, the interesting part of Y(t), i.e. the part f(t), is made up of several components:

f(t) =

å
m

a_m g_m(t),

and it is these components a_m g_m(t) (including the individual amplitudes a_m) which I call signals .

It is apparently the traditionat PSI to define a fitting function (using MuFit) as the sum of signals, whereas at TRIUMF the practice (at least in the ``chemistry'' MINUIT) is to define the full functions. The TRIUMF approach was to have a separate programme for each function, whereas MuFit is a more general-purpose package.

I have combined several aspects of these different approaches into MINFIT; there is only one programme, which fits to a choice of functions. One of these functions is rather special, and allows a choice of several signals to be summed into a tailor-made compound function. In general, I believe that the explicit functionapproach is somewhat more efficient, while the compound-function method does have more flexibility. A quick check showed that an explicit function fit takes about two-thirds of the time taken by the COMPOUNDfunction. Thus, I would recommend that if any particular combination of signals is used on a routine basis, a full function should be implemented for that application.

I have also written MINFIT in such a way that the fit can be applied to several histogramsat once, but now constrain each histogram being fitted to have the same length and packing factor. Whether a fit obtained with differing packing factors has any meaning, I leave open for debate...

To implement the ``functions'' and ``signals'', I have made extensive use of VAX FORTRAN STRUCTUREs , as well as several other VAX FORTRAN features. This has the potentially annoying side-effect of limiting portability, but most features appear to have their counterparts in the new FORTRAN '90 standard so this will cease to be a problem in the future. The structuresare defined in FCN_DEF.INC and SIGNAL.INC, respectively. The relevant portion of FCN_DEF.INC (used for defined functions) is

        structure /MINUIT_FCN/
        character*8     basename        ! base name for .con files etc
        character*8     main_fcn        ! name of the main fcn
        character*8     util_fcn        ! name of utility fcn
        character*64    description     ! description of function
        integer*4       fcn_ptr/0/      ! address, to be filled in
        integer*4       util_ptr/0/     ! ditto
        integer*4       getpar_ptr /0/  ! allow different get_params
        integer*4       base_pars       ! number of base parameters
        integer*4       hist_pars       ! number of additional params
                                        ! for each histogram
        logical         grad_calc       ! .TRUE. if fcn/util calculates
                                        ! partial derivatives
        end structure

The main things of interest to the user are basename , which is used to form default parameter-file names , as well as the interactive prompt, and the base_pars and hist_pars . The base_pars are those parameters which are shared by all histograms, such as field or frequency, relaxation factors, and so on, while the hist_pars are those which vary between histograms, such as normalisation, background, amplitude, phase, and so on. The ordering of the parametersis the base parameters first, followed by the set of histogram parameters for each fitted histogram in turn. ``Initial guess'' valuesof the parameters (including their names) may be read from a configuration (.CON) or parameter (.PAR) file (see later).

The heart of SIGNAL.INC (used to build up the COMPOUND function) is

        PARAMETER (MAXSIGS=8,MAXSIGPARS=8)

        structure /SIGNAL/
        character*32    description     ! short description of function
        integer*4       signal_fcn_ptr/0/ ! address, to be filled in
        integer*2       sig_share /0/, sig_uniq /0/
        integer*2       sig_hist /0/, sig_hist_share /0/
        logical*2       grad_OK         ! .TRUE. if fcn/util calculates
                                        ! partial derivatives
        character*2     psi_code /' '/  ! Old MuFit codes
        character*4     sig_par_name(MAXSIGPARS) /MAXSIGPARS*' '/
                                ! short name for signal parameter
        end structure

Here again, the salient point is the set of parameter countssig_share, sig_uniq, sig_hist, and sig_hist_share. A parameter listis built up (using the short names defined in sig_par_name) for each signal in turn. Firstly come the sig_share parameters, which are global parameters applying to all histograms, but which may be shared between different signals, as e.g., field in a muon-plus-muonium fit. The sig_uniq parameters are, of course, unshareable global parameters such as relaxation factors. The sig_hist parameters are those which vary from histogram to histogram (e.g. amplitudes), while the sig_hist_share parameters are those which vary from histogram to histogram but may be shared between signals (e.g. signal phases in a longitudinal-field experiment, assuming the t₀sare accurately known!). These parameters are grouped in sets of nhist, where nhist is the number of histograms being fitted so that, for example, the amplitudes for all histograms precede the phases for the histograms. The parameter list is either built up from manual input or read from a parameter (.PAR) file. Finally come the set of (normalisation, background) pairs for each histogram being fitted, as these only enter the calculation of c² and gradients after all the component signals have been summed (see FCN.TEX).

Available functions

STANDARD: This is a ``standard'' low-fieldfunction
LORENTZ: A function for a single Lorentzian relaxing muon signal:
GAUSSIAN: A function for a single Gaussian relaxing muon signal:
EXPROOT: A function for a single muon signal with exponential-square-root relaxation:
POWER: A function for a single muon signal with exponential-power relaxation:
ABRAGAM: A function for a single muon signal with Abragam relaxation:
COMPOUND: This is the generic name for functions which are built up from the signals described below:

Available signals

The signals currently available match the above functions, plus several longitudinal/zero-field signals culled from MuFit. I have retained the two-letter codes given to the signals in MuFit, extending the codes as necessary to cover new signals. Note that most oscillating signals occur in two versions, one for muon signals in which the field is given, and the other as an oscillation in which the frequency is explicitly given. This involves just a few more lines of code, plus an additional record in the signal ``database''. The labels below give the internal routine name and the MuFit mnemonic codes for each signal.

stdsig - MU: The low-field muonium-plus-muon signal as in the STANDARD function above:
muonsig - EM: A muon signal with Lorentzian relaxation:
gausssig - GM: A muon signal with Gaussian relaxation:
besselssig - BM: A muon signal with Bessel oscillation:
exprtsig - RM: A muon signal with exponential-square-root relaxation:
powersig - PM: A muon signal with exponential-power relaxation:
abrgmsig - AM: A muon signal with Abragam relaxation:
lorosc - E: An oscillation with Lorentzian relaxation:
gaussosc - G: An oscillation with Gaussian relaxation:
besselosc - B: An damped Bessel oscillation:
exprtosc - R: An oscillation with exponential-square-root relaxation:
powerosc - P: An oscillation with exponential-power relaxation:
abrgmosc - A: An oscillation with Abragam relaxation:
lorentzzer - KL
kubogauss - KG
lorcoszer - E0
gaucoszer - G0

There are a couple of other signals in MuFitthat I have not yet implemented within MINFIT. If anyone wants these (or any others) added to the list, please contact me with a nicely-written(!) formula for the signal.

Defining relations between amplitude parameters

It is sometimes desired to define ``relations'' between the amplitudes of different signals, especially in zero/longitudinal-field experiments. E.g. one might have determined from a transverse-field experiment what the total amplitude is for each telescope, and want to restrict the fit so that the signal amplitudes add up to this figure.

MINFIT contains the facility to define such relations in the COMPOUND function,with the following restrictions:

ONLY amplitudes may be defined as relations (although, curiously, parameters other than amplitudes are allowed to be included in the relation...).
A maximum of eight relations may be defined in any given function.
A maximum of four ``coefficients'' may be included in any relation.

A relation is definedby a string consisting of one to four coefficient/parameter pairs, where each pair has the format

                        <sign>[coefficient][*][<P>n]

where <> indicates a required item and [] an optional item, although at least one of the [coefficient] or [<P>n] must be present to make any sense. I.e. the sign is required (unless it's the first pair in the string, when a `+' may be omitted) and the multiplication sign is always optional. If a coefficient is omitted, it is taken to be 1.0, and if the parameter indication (Pn) is missing, the coefficient is taken to be a constant. Leading and trailing blanks are ignored, except that blanks in the coefficient part may be interpreted as zeros if there is no explicit decimal point.

A few examplesshould serve to illustrate the format:

To indicate that the sum of two signal amplitudes is to be constant, the second can be given as .2-P3, where the first amplitude is parameter 3 in the function.
Likewise -P3-P7-P11+.68, given for a fourth amplitude, indicates that the sum of four amplitudes is to be a constant.
The relation 0.2 P5 indicates that the amplitude is to be one-fifth of parameter 5.

Relations are entered manuallyby including an 'R' (or `r') as the first character in the response for a sig_hist type parameter (currently only amplitudes), followed by the relation. In a .PAR parameter file, they are indicated by an `R' in column 1, and the relation string is read from columns 30 to 80. This is to leave space for the PARFILE command to write the derived parameter name and valuein the ``normal'' format when writing a .PAR file. It also leads to the restriction on the number of coefficients - although the relation is written as compactly as possible, the general format for a coefficient/parameter pair is (SP,1PG10.3,'P',I2.2), which will in fact overflow the available 51 characters if four pairs need to be written out in full (even more so if parameter 100 is included in the relation!). I would anticipate that most relations would be entered by editing existing .PAR files, as it is less ``user-friendly'' to enter them manually.

There are a few things to be aware of when using relations. Firstly, the parameter defined by a relation does not exist, as far as MINUIT is concerned. Therefore it does not appear if you issue commands such as SHOW PARAMETER. To see the numerical value, use the SHOW RELATION command. You will also find the value in a .PAR file written with a PARFILE command, or in a .TXT file produced with the WRITE command. Since it is not a parameter, a relation has no parameter number, so be sure when editing .PAR files that real parameters are numbered consecutively!

Secondly, defining a relation destroys the ability to calculate directly the partial derivatives, so that fits will take somewhat longer because MINUIT has to use differences to get the gradients. This will be more noticeable for fits with many parameters.

Preparations for MINFIT

Firstly, note that there are two versions of MINFIT; the second, called MINFITV , was compiled to support the VAX-9410's Vector Processor, and can give a slight performance increase ( ~ 10%). The vector version will run only on PSICLC, restricting its use to batch runs for most people. The results are indistinguishable between the two versions.

MINFIT contains some rudimentary plotting routines to allow interactive users to inspect the progress of the fit (suggestions for enhancement of this facility are welcomed; see description of the PLOT command below). This uses the graphics package GRAPHX , so this must be IMPORTedbefore MINFIT can be used.

MINFIT logicals

To help locate various files, MINFIT uses VMS logicals which need to be set before it is run. These are

MINFIT_DATA:: MINFIT_DATA is the name of the directory in which the data resides. e.g. to analyse data taken on LTF while working on the cluster one would issue the DCL command
MINFIT_NAME:: MINFIT_NAME contains a pattern to be used by MINFIT to generate data-file names. The pattern is of the form
MINFIT_HELP:: MINFIT_HELP should point to a directory containing the MINFIT help library MINFIT_HELP.HLB. On the cluster, use the MINFIT directory:

MINFIT limitations

Wherever possible, I have used virtual memory allocationto avoid having limitations hard-wired into MINFIT. However, in some instances it has been unavoidable to impose some limits. Current limitations are:

MINFIT uses the standard MINUIT distribution (even though I maintain my own private copy of it). This limits the fitting function to 100 parameters, of which a maximum of 35 may be variable at any given time. MINUIT indicates how many parameters are allowed, and how many may be variable, in the text marked DIMENSIONS when MINFIT is initialised.
The maximum number of histograms is set at the limit imposed by the PSI data formats, i.e. 16.
The maximum number of bins fitted, per histogram, is essentially unlimited, but the PSI data format uses an INTEGER*2 variable to hold the histogram size. I interpret this as unsigned (i.e. any negative value has 65,536 added to it) and, in addition, consider a zero value to be 65,536. size.
I have arbitrarily limited the COMPOUND function to a maximum of eight signals. I can't see that one could have a meaningful fit with this many signals; if so many are necessary, another method of analysis is probably indicated!
Similarly, each signal is limited to a maximum of eight formal parameters (the sum of sig_share, sig_uniq, sig_hist and sig_hist_share). The actual number of parameters per signal may exceed eight, of course, since the histogram parameters are (or may be) duplicated for each histogram fitted.
See the earlier section ``Defining relations between amplitude parameters'' for the restrictions imposed on relations.

Initial guess files for functions: the .CON files

To ease the entering of initial values for the various pre-configured functions, MINFIT will attempt to read these parameters from an initial-guess file basename.CON, where basename is the function name as enumerated above. Should this file not be found in the default directory, the user is prompted for another file name. The extension .CON is the default, and need not be explicitly given.

The .CON file contains parameters in the form

<parno> <parname> <value>,<uncertainty>,<lower bound>,<upper bound>

where <parno> is the parameter number, <parname> is the name of the parameter, and the other elements are self-explanatory. The format under which these files are read (with CARRIAGECONTROL='LIST') is (I3,1X,A10,4F20.0). This means that just the first two elements must be formatted as specified if the numeric elements are separated by commas, which over-ride the F20.0 specifications - the last element should contain a decimal point, however, to ensure the desired interpretation! If bounds are omitted then the parameter is considered unbounded(see the MINUIT documentfor a discussion of the effects of bounding on a parameter). If just one bound is given, the other will be take to be 0.0, possibly accompanied by a MINUIT warning that the bounds were reversed. If the uncertainty (initial guess at the parameter error) is also omitted, or set to zero, the parameter is considered to be constant.

Parameters should be entered in order (i.e. the <parno>s should be in strictly ascending order). All the base parameters must be entered first, followed by as many sets of histogram parameters as are necessary for the number of histograms being fitted. It is preferable to have several sets of histogram parameters in the file, as an error occurs if there are too few parameters; it is not an error to have more parameters than are necessary, as any excess parameters are ignored. This may lead to ``bad'' initial values of, say, phases, if one fit is made to histograms 1 and 2, and then another to histograms 3 and 4, since MINFIT takes the histogram parameters strictly in order without regard to which histogram is actually being fitted. In cases like this, it may be preferable to set up different parameter (.PAR) files for each case - see below.

MINFIT passes the parameter namesfor the base parameters directly through to MINUIT (does anyone have a logical explanation for the A10 format?). For the histogram parameters, each parameter name is stripped of trailing blanks and trimmed back to at most 7 characters, and the string ``_NN'' appended, where NN is the number of the histogram the parameter is associated with, before being passed to MINUIT. Thus if the name in the .CON file was ``Phase'', and histogram 3 was being fitted, MINUIT would receive the parameter name ``Phase_03'', whereas ``Asymmetry'' would be transformed to ``Asymmet_03''.

Example .CON files for all the functions may be found in the MINFIT directory.

Initial guess files for signals

Because the COMPOUND function may be constituted in a multitude of ways, there is no sensible way to make up a .CON file for this function. However, a facility is available to read in the complete set of fitting parameters, together with such things as histogram limits and binning, and so on, for all functions including COMPOUND. A discussion of this is deferred until after manual data entry is described - see the description of .PAR files below.

Running MINFIT

Assuming GRAPHXhas been IMPORTed, and any necessary .CON files are available in the default directory, MINFIT is started by RUNning the appropriate executable file. MINUIT outputs its status on start-up, including details of how many total and variable parameters are allowed, and it's idea of what the machine accuracy is. The first question asked by MINFIT is

Read fit parameters from a file <Y>?

Note that wherever possible, I have indicated default answers to question prompts within angle brackets <...>. Here the default is to use a parameter file. Answering ``N'' will result in prompts for manual data entry.

Manual data entry

Run parameters

The first task of MINFIT is to establish which run is to be fitted:

Enter run number:

to which one naturally replies with the desired number. This is inserted into the format given by MINFIT_NAMEand the file searchedfor in the current directory . If it is not found there, the directory indicated by MINFIT_DATAis then searched. Should this also be unsuccessful, MINFIT terminates with suitable grumblings about the sanityof the operator.

Assuming the file is successfully found, MINFIT then allocates sufficient virtual memoryand reads in the data. The data are maintained in memory to facilitate changes of parameters such as binning range.

Next, the user is prompted for the desired histograms:

Which histogram(s) do you want to fit (^Z quits)?

Here, the user enters the number(s) of the histogram(s) to be fitted. MINFIT can be aborted at this point by typing Ctrl-Z. Note that there is no default on the histogram number(s). MINFIT will re-prompt if an invalid histogram number is given. The histograms which MINFIT thinks you want are then displayed, giving you a chance to change the selection.

For the first selected histogram, the user is now prompted for the binningto be used:

Enter [first bin] [,last bin] [,packing factor]<NN,MMMM,1>:

where NN and MMMM are the first and last good bins as given in the run header. Note that the histogram bin numbers are taken to start from zero, not one! As indicated, the defaults are to take the stated good bin range, with a packing factor of one. The square brackets mean that the item is optional, so the above should be taken to mean

just entering a carriage-return takes the default values.
entering ``,,n'' takes the default bin range and packs by n.
entering ``,m'' takes the bin range from the first good bin to bin m and bins by one.
entering ``l'' takes the bin range from bin l to the last good bin, packing by one.
entering ``l,m,n'' takes the bin range from l to m and packs by n (naturally!).

and so on... The commas are necessary delimiters when taking default values for items followed by non-default values.

For subsequent histograms, only the first bin can be entered, as MINFIT now forces all fitted histograms to have the same length and packing. The default is again the first good bin, and the default end bin and packing are shown for information.

The next prompt is for the t₀ value:

Time-zero is at bin MMM.0: Enter replacement:

where MMM is the time-zero bin given in the run header. You can take this as the default value by pressing <return>, or enter a new value. Conceptually, the zero value is at the mid-point of the bin, e.g. if binning by one with 2.5 ns bins with t₀ at MMM.0, the time for bin MMM+N is 2.5NN ns, or if binning by two, the time for the bin (MMM+N)+(MMM+N+1) is 2.5N+1.25 ns.

Function parameters

If a predefined function is being fitted, MINFIT now reads in the initial parameters from the .CON file, which is transparent to the user unless there is a problem with the file.

Signal parameters

If the COMPOUND function is being used, you must now define the signals to be used and give values for their parameters. This process cycles around a loop for each signal.

Firstly a list of short descriptions of the available signals is presented, along with their MuFitmnemonics. A signal is selected by entering either the signal's listed number or its mnemonic. Just pressing <return> terminates signal selection. It is an error to enter no signal at all and MINFIT will die a terrible deathif you do so (I know, I just tried it...).

For each signal parameter in turn, MINFIT builds up a parameter name from the defined 4-letter sig_par_name by adding ``_SX'' to base parameters and ``_SXHYY'' to histogram parameters, where X is the signal number (in the fit, not the list of signals!) and YY is the histogram number. If the parameter is ``shareable''and a parameter of the same type has already been defined, you will be asked if you want to map it to a previously defined parameter. The default is ``no''. If you answer ``yes'', a list of previously-defined parameters is presented for you to choose from.

If a parameter is not to be shared, MINFIT prompts for value, error, lower, and upper bounds for the parameter. These should be entered as comma-separated values. If no bounds are given, the parameter is taken as unbounded. If no error is given, or if it is given as zero, the parameter will be taken to be fixed.

The parameters are entered in the ordersig_share parameters, followed by sig_uniq parameters, and then the sig_hist and sig_hist_share parameters in multiplicity depending upon the number of histograms.

After all the signals have been defined and their parameters given, the normalisation and background for each histogram are then entered. Here, you need only enter a value, as an uncertainty of 5% is assumed. You can give an explicit uncertainty, however. Note that it is not possible to enter a zero uncertainty to make one of these parameters constant; you will need to use the FIX command if you are perverseenough to want to do this. These parameters are assumed unbounded ; to make them bounded requires the use of a CHANGE or SET LIMITS command.

Automatic data entry - the .PAR file

To avoid all the boring rigmarole described above, it is possible to use a parameter file to enter the run and fit parameters. This file is used if the user answers ``yes'' to the initial MINFIT question. With the above description of manual data entry completed, we can finally proceed to a description of the parameter file, which has a default extension of .PAR.

You will be prompted for the name of the .PAR file, as no default can be established. As with the .CON files, .PAR files are assumed to be in the default directory unless the name explicitly contains a different specification.

The content of the .PAR file closely follows the sequence of manual data entry. Note that if at any time an error occurs while reading the parameter file, MINFIT reverts to manual data entry.. The first line should contain the name of the function to be used (format A8). The second line gives a run number, which is presented as a default run number, allowing the user to choose another run if required.

The next line contains the histogramsto be fitted, in a comma-separated list. At this point, you are asked whether these histograms are OK. If not, input will revert to manual. Next in the file comes a line for each selected histogram in turn, giving first bin, last bin, packing factor, and t₀. As usual, the format is over-specified to allow the use of commas as delimiters. For multiple histograms, the length and packing of the first histogram are used subsequently, to force all histograms to be equal in length.

For all but the COMPOUND function, the function parametersthen follow, in exactly the same format as in the .CON file, and they are interpreted in the same manner. A sample function .PAR file follows:

LORENTZ
   1951
 1, 2, 4,
   93, 4040,  1,   90.0000
   92, 4040,  1,   89.0000
   94, 4040,  1,   91.0000
  1 Field (G)  100.05551    , 2.08382677E-02
  2 MuonRelxn 1.07541315E-05,  1.74050095E-03
  3 Normal_01  1969.7874    ,   1.7836021
  4 Bckgrn_01  20.871749    ,  0.21864081
  5 MuonAsy_01 0.17452410   ,  1.25174937E-03
  6 MuonPha_01 178.21454    ,  0.41584137
  7 Normal_02  2704.7217    ,   2.1024913
  8 Bckgrn_02  29.217813    ,  0.25660021,0,500.0
  9 MuonAsy_02 0.17892908   ,  1.10932507E-03,0,0.35
 10 MuonPha_02 88.514591    ,  0.36190107
 11 Normal_04  2819.3928    ,   2.1489219
 12 Bckgrn_04  35.222845    ,  0.26944879,0,500.0
 13 MuonAsy_04 0.18047479   ,  1.11025710E-03,0,0.35
 14 MuonPha_04 -2.3555467   ,  0.35636249

For the COMPOUNDfunction, the next line gives the number of signalsto be defined. The parametersfor each signal are grouped together in turn, preceded by a line giving the MuFit mnemoniccodefor the particular signal. Parameter orderfollows that for manual input, and again the format is identical to the entries in the .CON files. Parameter numbering is in ascending order, except for those parameters which are mapped to previous parameters- these have the parameter number for the actual parameter, and the parameter name and values, if given, are ignored. Finally, the normalisation and background are given for each histogram in turn. If an error occurs at any time while parameters are being read, the signal definitions are cleared and MINFIT reverts to manual data entry.

This is a sample .PAR file for the COMPOUND function:

COMPOUND
   5762
 1, 2,
  255, 8000,  1,  250.0000
  255, 8000,  1,  250.0000
 4
G 
  1 Freq_S1     35.300000    ,  .10, 30.0,40.
  2 Sigm_S1     14.700000    ,  1.0,  10.0    ,  20.0    
  3 Asym_S1H01 0.20100000    , 2.0E-02, 0.000E+00, 0.50000000    
  4 Asym_S1H02 0.21000000    , 2.0E-02, 0.000E+00, 0.50000000    
  5 Phas_S1H01  50.000000    ,  2.0, 0.000E+00,  100.    
  6 Phas_S1H02 -135.20000    ,  2.0, -200.    , -60.0    
G 
  7 Freq_S2     31.600000    ,  2.0,  20.0    ,  40.0    
  8 Sigm_S2     26.000000    ,  2.0,  15.0    ,  40.0    
  9 Asym_S2H01 5.00000000E-02, 1.0E-02, 0.000E+00, 0.10000000    
 10 Asym_S2H02 6.00000000E-02, 1.0E-02, 0.000E+00, 0.10000000    
  5 Phas_S1H01  50.000000    ,  2., 0.000E+00,  100.    
  6 Phas_S1H02 -135.20000    ,  2., -200.    , -60.0    
G 
 11 Freq_S3     37.700000    , 0.1,  33.0    ,  43.0    
 12 Sigm_S3    0.95000000    , 0.1, 0.000E+00,  2.00    
 13 Asym_S3H01 8.80000000E-03, 1.0E-03, 0.000E+00, 2.000E-02
 14 Asym_S3H02 6.30000000E-03, 1.00, 0.000E+00, 2.000E-02
  5 Phas_S1H01  50.000000    ,  2.0, 0.000E+00,  100.    
  6 Phas_S1H02 -135.20000    ,  2.0, -200.    , -60.0    
E
 15 Freq_S4     322.,.1,315.,330.
 16 Relax_S4    0.,.01,-1,1.
 17 Asym_S4H01  .0233,.0007,0.,.05
 18 Asym_S4H02  .0201,.0007,0.,.05
 19 Phas_S4H01  134.1,2,0,400.
 20 Phas_S4H02  125.9,1,0,400.
 21 Normaln_01  3070.0000    ,  10.0
 22 Bckgrnd_01  14.000000    , 0.700
 23 Normaln_02  2972.0000    ,  10.0
 24 Bckgrnd_02  17.000000    , 0.850

By the way, if these examples look a little bit hacked around, that's because they are. We will see later that you can produce nice regular .PAR files with the PARFILE command,but the flexible input format means that you don't have to stay with something so pretty.

FIXed parameters and .PAR files

To enable parameters which have been FIXed to be saved in the .PAR files and recovered as fixed (as opposed to constant parameters,which you cannot change), the sign of the error value, as written in the file, is used as a flag.Normal parameters have positive errors, fixed parameters have negative errors and constant parameters have an error of 0.0 of course. When a parameter is read from a .PAR file and the error is found to be negative, it is set as a normal parameter (with positive error) and then FIXed.

MINFIT commands

Once everything has been initialised, MINFIT enters its command interpreter, showing as a promptthe name of the selected function. Most of the commands are direct MINUIT commands,so you should read carefully Chapter 4 of the MINUIT documentation.There are several commands which are not passed on (e.g. END, EXIT, RETURN), as they cause MINUIT to execute a FCN call with IFLAG=3, which clears the virtual memory allocated to data storage, causing a crashon the next FCN call. For the same reason, you should avoid executing a CALLFCN 3 command.

Note that the commands MIGRAD,MINIMIZE,MINOSand SEEKmay be entered with a zero first argument (e.g. MIGR ,5,6,7), in which case a default value of 1000 function calls is passed on to MINUIT.

As in the MINUIT documentation, commands are described in the following format:

<command-string> <arg1> [arg2] etc.

<command-string>: is one of the commands listed below,
<argi>: are numerical values of required arguments, if any.
[argi]: are numerical values of optional arguments, if any.

The arguments (if any) are separated from each other and from the command by one or more blanksor a comma. Commands may be given in upper or lower case, and may be abbreviated,usually to three characters. The shortest recognized abbreviations are indicated by the capitalized part of the commands listed below (not in strict alphabetical order!).

Where a parameter number is required, it may be replaced by the parameter name, subject to some restrictions. The name is textually replaced by the parameter number before the command line is parsed, so the name must be exactly given, except for trailing blanks; the comparison is case-insensitive, however. To allow for the textual replacement, the parameter name must be long enough to be overwritten by the number, and it also should not be a superstring of another parameter, nor a substring of a MINFIT command.

?
A question mark (?) is interpreted as the command HELP.

[command] ! [string]
An exclamation mark (!) is used to introduce a comment. Anything on an input line to the right of an exclamation mark (including the exclamation mark itself) is discarded before the input is parsed into a command. This allows command sequences in batch files, for example, to be documented in a minimal fashion. See also the COMMENT command, which inserts comments into the output. Note that an exclamation mark within a COMMENT command is discarded before the COMMENT is echoed.

@[filename]
A command starting with an ``at'' sign (@) designates an indirect command file from which further commands are to be read. The default extension is .DAT(for compatability with TRIUMF MINUIT).

Commands are then taken from the named file until:

a further @-command is encountered, upon which the current file is closed and the newly-named file is opened to read more commands - if the file cannot be opened, input reverts to SYS$INPUT (this allows an indirect command file to be terminated by @GOBBLEDYGOOK or simply @); or,
the end-of-file is encountered, after which control passes back to SYS$INPUT.

Note that although indirect command files can be nested, they are not ``stacked''; i.e. a file is forgotten when a further file is executed and input does not pass back to the first file when the second ends.

Indirect command files can be used to store frequently-used command sequences, and can also be helpful when running MINFIT in batch.

CHAnge [parameter] [value] [error] [lower bound] [upper bound]
The reason for this command is that MINUIT offers no opportunity to change a parameter's value and limits at the same time, so if the new value is outside the existing limits, one must first change the limits, then the parameter, then the limits again!

When entered with no parameters, CHANGE enters an ``interactive'' mode, displaying all the parameters and prompting for a parameter to change. Once a parameter is selected the value, error, and limits can be entered. Any option left out will be left unchanged. Note that to retain the value as unchanged, you must enter ,error..., as entering 0,error... will set the value to zero. Note also that it's not possible to set the error to zero! The interactive mode is exited by entering a null value for the parameter to be changed.

When entered with just a parameter number, CHANGE does not display the parameter list, but goes straight to the prompt for new values.

When entered with values as well as a parameter number, CHANGE will set the parameter and return directly to the command prompt. Note that you must enter a new value in this mode, although error and limits can retain their old values as defaults.

COMment <comment string >
This command writes the entire command line to the output device (terminal or log-file). Its purpose is to allow a minimal documentation of the output (compared to the ! command, which documents the input), so that later scanning of the output can see what has been done. This is mainly useful with batch jobs, using indirect command files and conditional commands (which are not yet implemented...).

e.g. Comment: MIGRAD didn't converge, so I'll relax limits.

FUNction
This allows the user to choose a new function to be fitted to the existing data file. Currently, new parameters must be entered manually.

NEWfile
This allows the user to read in a new data file. There is now a default option to retain the same histograms and binning, and the current fit parameters when fitting the new file. This facilitates analysing a series of runs with similar properties.

REBin
This allows the user to choose new binning for the current histogram(s). The current parameters are retained.

PARfile
This allows the user to write out the current fit settings into a .PAR parameter file, to be used as input into a later fit session. The default file-name is the current function, and the default extension is .PAR.

PLOt [histogram] [mode] [start bin] [end bin] [packing] [binflag]
This allows the user to see the data and the fitted function. The form of the graph is determined by the current PLOTLEVEL (see SET PLOTLEVEL).

There are four plot modes:

Raw data: The data is plotted, along with the current fit.
Asymmetry: The current background and the muon decay are extracted from the data and the current normalisation used to arrive at the asymmetry:

p(t) = (d(t)-b)/(Ne^-t/t)-1.

The current fit is also plotted.

Residuals: The relative difference between the data and the current fit is plotted;

p(t) = ( d(t)-Y(t) ) /Y(t).

Logarithmic: This is a raw data plot with logarithmic count axis.

The start and end bins are relative to the bins being fitted, i.e. the first fitted bin is bin 1. The packing is also the number of fitted bins summed together for display, independently of the binning that these may have from the run data. The binflag controls whether the X-axis is given as time since t0 or as the bin number in the original histogram (this is to allow t0values to be more readily determined). A non-zero value of binflag leads to the bin numbers being shown; zero means to plot as a function of time.

All parameters retain their last value as the default. The initial defaults are the lowest numbered histogram being fitted, mode 1, first and last fitted bin, binning by one, and plotting vs. time. If the specified end bin is too large, it will be truncated to the last fitted bin. Initially, one would most likely want to look at the data and the fit with only a reasonable number of plotted points, with a command such as PLOT ,,,,5.

SET PLOtlevel <level >
This sets the level of plot complexity and determines what is included on the graphs produced by the PLOT command. You can examine the current level with the SHOw PLOtlevel command. The defined levels are:

0 (default) Plots specified/default histogram.

1 Plots all fitted histograms.

2 Plots specified/default histogram plus parameter values.

3 Plots all fitted histograms plus parameter values.

LISt <something >
This is an extension of the MINUIT SHOW command, allowing the data to be written to a file instead of displayed on the screen. The initial use of LIST will prompt for a file name; the default extension is .FIT.

SHOw FUNction
This command shows what function is currently being used to fit data. If it's the COMPOUND function, a list of the signals being used is also given.

SHOw ALL
This command gives a listing of the run parameters (as for SHOW RUN) plus the c2 values (as for SHOW CHISQUARE) as well as a listing of the current values of the fit parameters and any relations which are defined. The SHOW FUNCTION output is also given.

SHOw CHIsquare
This command shows the currentc2 value, as well as the overall reduced c2 value and that for each histogram.

SHOw RELation
This command will write out the name, value, and definition of any parameters which have been defined using relations.

SHOw RUN
This command writes out the current run number, title and subtitles, together with information on the histograms currently being fitted.

SHOw TIMe
This command writes out the elapsed CPU and real time since MINFIT was started. This is mainly useful for performance comparisons, between different computer models, or different fitting functions.

LOOsen [factor]
This command causes MINFIT to expand or ``loosen'' the error estimate on all variable parameters by the given factor (default is 5.0). This may be useful when you have changed a parameter after a fit and want the other parameters to have reasonably large step sizes in the new ``initial guess''. The default factor of 5.0 may be a bit large for the new version of MINUIT, which doesn't seem as prone as older versions to keep reducing the step size to ridiculously small values - experiment with different values to get a feel for how much expanding (or contracting) the errors changes the number of function calls needed to arrive at a converged fit.

If the expanded error would overlap a limit on a bounded parameter, the error is reduced to avoid the overlap. If the expanded error is incompatible with limits (i.e. 2*error > (lim2 - lim1)) the error is left unchanged.

This command comes from the TRIUMF fitting programmes. Donald Arseneaunotes ``LOOSEN is not a magical solution to being stuck away from the true minimum. If MINOS keeps finding new minima, it is more likely that the covariance matrix is bad than the errors are too small. In such a situation, try HESSE.'' One could also perhaps try increasing the strategy level (see SET STRATEGY).

WRIte
This command allows the user to write the current parameters out to a MuFit style .TXT file. It functions almost identically to the corresponding option in MuFit, but I have made a few extensions. MINFIT tries to write the MINOS errors . If they are unavailable, it will use the parabolic error estimate. If that has not been calculated, it uses the parameter error estimate (which should be the same as the parabolic error if any fitting has been done). If the MINOS erros aren't available, the error is written as in MuFit. If they are available, the positive error is written as usual, preceded by a + sign instead of the usual ± sign, and the negative error is written on the next line, preceded by a - sign which is aligned with its partner above. This second line is completely blank except for the negative MINOS errors, and is not included if there are no errors to be given.

To accommodate run numbers greater than 9999, I have also made the run-number in a .TXT file an I5 format instead of I4, reducing the following blanks so that the other fields stay in the same columns.

I have also introduced the suffix _shr, which is used for parameters which are defined as being sharedwith another. _def is now used only for those parameters defined by a relation. As an illustration, here is a sample output from a .TXT file, where the two outputs were written before and after MIGRAD and MINOS commands. I've had to change the ± characters to asterisks because of limitations of LATEX's verbatim environment, and I've also chopped off the histogram identification in column 132.

 Run 1956   PMS8-32    8.000K     100G       -              1add   93   92 4040 4040
    HI2=1.028         B0=20.76*0.22     N0=1940.5*1.8
    HI2=1.015         B0=29.09*0.26     N0=2671.2*2.1
                   GM A1=10.4*0.8       P1=181.4*1.4      F1=103.6*0.6      S1=0.65*0.04
                   GM A1=11.2*0.8       P1=86.9*0.8
                   GM A2=6.58_def       P2=179.9*1.8      F2=88.7*0.5       S2=0.434*0.030
                   GM A2=5.75_def       P2=86.9_shr
 Run 1956   PMS8-32    8.000K     100G       -              1add   93   92 4040 4040
    HI2=1.028         B0=20.76*0.22     N0=1940.5*1.8
    HI2=1.015         B0=29.09*0.26     N0=2671.2*2.1
                   GM A1=10.4+0.8       P1=181.4*1.4      F1=103.6+0.6      S1=0.65+0.04
                             -0.8                                 -0.6             -0.04
                   GM A1=11.2+0.8       P1=86.9*0.8
                             -0.8
                   GM A2=6.58_def       P2=179.9*1.8      F2=88.7+0.5       S2=0.434+0.031
                                                                 -0.5               -0.030
                   GM A2=5.75_def       P2=86.9_shr

QUIT
Leaves MINFIT.

Indirect command files

As noted above, command input can be switched to an indirect command file with the @<filename> command.Input is taken from the file (default extension .DAT) until the file ends, input is switched to another file by a further @<filename> command, or an invalid @-command is read. At end-of-file or upon encountering an invalid @-command, input reverts to SYS$INPUT (keyboard or batch file). Note that the indirect files are not ``stacked'' - at the end of a second indirect file, input reverts directly to SYS$INPUT, not to the calling file.

This facility is intended for frequently-executed command sequences, or to facilitate the control of batch jobs.

Setting parameter limits

You should read section 5.3 of the MINUIT documentation for a discussion on the effects of setting parameter limits. In general, it is unavoidable to set error limits to guide MINUIT in the right direction. This is especially true of some relaxation parameters - it is obvious that the relaxation factor in the EXPROOT function cannot be allowed to become less than zero, or the programme will crash. What is perhaps not so obvious is that it may be necessary to constrain even the lin a Lorentzian signal, since if MINUIT decides to try a value which is too far negative, the exp function will crash with a floating-point overflow.

In general one finds that parameters with imposed limits show a parabolic error estimate which is less than the MINOS error (I don't recall ever seeing it significantly larger).How much less is some function of the extent of the limited region and how symmetrically the limits are placed around the parameter value. In general it is wise to try to obtain the MINOS errorson the parameters of interest, especially for relaxation parameters whose errors are often non-symmetric.

Batch MINFIT (partially obsolete...)

Since MINFIT avoids using the dreaded VAXLIB input routines, it is quite easy to run it in batch.

In general, the easiest way is just to set up a .COM file consisting of all the DCL commands necessary to set up MINFIT (importing GRAPHX, setting logicals, etc.) and start it up, and then just incorporate the MINFIT commands into the command file at the appropriate place (without a leading $-sign, of course). It is also possible to put all the MINFIT commands into a separate file which is defined as input by using the DEFINE/USER SYS$INPUT myfile.name command immediately before running MINFIT.

It may also be advantageous to use indirect command files (see above) to store frequently-used sequences of commands.

An example of a batch command file is the following:

$       import graphx
$       set def disk_244_dat0:[reid.minfit]
$       @txdata
$! txdata.com sets up the appropriate logical definitions
$       r minfit
Y
COMPOUND


SET PRINT 0
MIGRAD
SHOW CHI
MINOS
QUIT
$       exit

An example of the second method is:

$       import graphx
$       set def [reid.minfit]
$       @txdata
$       define/user sys$input r5762.inp
$       r minfit
$       exit

where R5762.INP contains the data

Y
COMPOUND


! That was the prelims; now we get to the MINFIT command interpreter
! Note the use of comments -- the first character must be a !
SET PRINT 0    ! Get more printout
SET STRATEGY 2
! Try to be a bit more rigorous, since it's a difficult fit
MIGRAD         ! Do the fit
SHOW CHI       ! Look at the reduced chi-squares
MINOS          ! Get the full minos errors
QUIT

Note that in both these examples, blank lines are used to accept the default run number and histograms.

For batch jobs, the use of .PARfiles is recommended, and a judicious use of the NEWRUN command should enable a sequence of runs to be analysed sequentially. Currently, however, MINUIT (and hence MINFIT) has no strategy for recovering from errors such as MIGRAD termination without convergence (e.g. due to a parameter coming up against a limit). It is recommended that a suitable set of starting parameters be found interactivelyso that convergence is obtained, and batch jobs can then be used for the more tedious calculations such as obtaining MINOS errors.

It has been noticed, however, that it is often advantageous to SET STRATEGY 2 when fitting, especially for complex fits, or when MINOS errors are required. Setting strategy to 2 seems to require more FCN calls in MIGRAD, but is offset by quicker subsequent convergence in MINOS. It is especially recommended for batch jobs. The following is an illustration of using a batch job with indirect command files to fit a series of runs. This particular file had great problems with the last run in the sequence until the strategy level was increased. First, the command file:

$ import graphx
$ def/nolog minfit_data psm020::dub0:[musr.exp.td_musr.dlog]
$ def/nolog minfit_name deltat_imp_NNNN.dat
$ set def disk_244_src0:[reid.minfit]
$ minfit:==run disk_244_dat0:[reid.minfit]minfit
$ if f$getsyi("NODE_HWTYPE") .eqs. "9XXX" then -
        minfit:==run disk_244_dat0:[reid.minfit]minfitv
$! Note the use of F$GETSYI to determine if we can use the Vector version.
$ minfit
Y
R19xx

Y
set strategy 2  ! Improves behaviour
! (A sum of 2 Gaussians may not be the best fcn for these runs!).
set print 0     ! More printout (e.g. eliminates the need for SHOW MINOS)
@R19xx ! Which see...
new
y
1956
loosen 2.
! Expand errors a bit after the last fit (test the effect of this on your fits)
@R19xx
new
y
1957
loosen 2.
@R19xx
new
y
1958
loosen 2.
@R19xx
new
y
1959
loosen 2.
@R19xx
new
y
1960
loosen 2.
@R19xx
quit ! All done
$ exit

The include file R19XX.DAT contains:

migrad
sho run ! We now get some gratuitous space before indirect commands...
sho chi
sho tim
mino 2000,1,2,3,4   ! MINUIT only MINOSses a maximum of 6 directly-specified
mino 2000,7,8,9,10  ! parameters, so do it in two bits.
sho tim

Fitting TRIUMF data

MINFIT also now understands the TRIUMF data format, so it is no longer necessary to convert to PSI format with TRI2PSI. The recognition is automatic, based on the fact that PSI files have variable-length records, while the TRIUMF files have fixed-length records. It is only necessary to make sure that MINFIT_DATA and MINFIT_NAME are appropriately defined. In fact, it is possible to analyse data directly from TRIUMF (although the initial file read is slightly slow) by defining, e.g.

$ DEFINE MINFIT_DATA ERICH::PRV37:[M15.1991]
$ DEFINE MINFIT_NAME NNNNNN.TRI

Note, however, that pre-1985 data archived on WRM0: is currently (June 1991) stored in IBM format, as variable-length records (although they are all, in fact, 512 bytes).MINFIT currently does not understand this particular file format.

section*Fitting RAL data It is also possible to use MINFIT to work with RAL files (these are mainly the same as PSI files, but with a larger number of histograms). Please ask me or Wolfgang Templ for further details.

Transferring PSI data files

Because PSI data files are variable-length records, they are easily corrupted by transferral via non-VAX means. In particular, ftp transfer will always corrupt a PSI file. The problem:

ftp, coming from the Unix world, knows only two types of files: pure binary, and pure ASCII. In binary files, there are no delimiters. In ASCII files, record delimiters are <LF> characters (a.k.a. newlines).
MS-DOS also knows the two file types, but takes the more usual record delimiter of <CR><LF> for ASCII files. Binary files have no record
VMS has many different file types, and the interpretation can depend on the mode specified on opening the file. As with all files, ultimately the storage is a binary file, allocated in sectors or blocks (typically 512 bytes). The logical transformation into records depends on information within the file header and the file itself. The most usual file structure for ASCII and data files is ``variable'' in which each record is preceded by two bytes giving the length of the rest of the record. This can be further complicated by various utilities or filing services (e.g. RMS) embedding further control information in the data. FORTRAN compilers (both VMS and MS-DOS flavours, at least) are very fond of doing this.

For the particular mode chosen for PSI data-files, FORTRAN does add extra control information. The file is opened specifying that the records should be ``segmented'' with a ``blocksize'' of 2048 bytes. This means that the data for a logical record (i.e. all the data items specified in a single Fortran unformatted write statement) are split into ``segments'' of not more than 2044 (it's less than 2048 because of the extra control information) bytes and each segment is written as a variable length record (i.e. with the 2-byte header). The control information consists of two bytes per record indicating whether this segment is the first, last, middle or only segment of the logical record. This allows the logical record to be reconstructed from the segments.

So in the case of the PSI header block, the 1024-byte header has two bytes added (saying ``this is the only segment in this segmented record'') and is then written to disk as a 1026-byte variable-length record (occupying 1028 bytes in total). A 4K histogram, recorded as one logical record, is split up into eight 2042-byte segments and one 48-byte segment, which are written to disk (with the control information added) as eight 2044-byte variable-length records and one 50-byte variable-length record. Note that the maximum length written per record is only 2046 rather than the 2048 specified in the blocksize parameter. This is probably because, if the file were being written to magtape according to ANSI standards rather than to disk, one segment would become one block on magtape and ANSI block headers consist of a four byte length field rather than the two bytes used in disk files.

There is an extra complication with histograms longer than 4096 bins. In this case, the histogram is split into KDAFHI ``logical records'', each of length LENDAF bins, where LENDAF is 4096 or less. The last ``logical record'' of the histogram is padded with zeros, if necessary, so that all ``logical records'' are the same length (KDAFHI and LENDAF are stored in the file header).

If one ftps a PSI data file in ASCII mode, then the VAX side of the transfer apparently opens the file in such a way that the record-length bytes are interpreted and removed before passing the file contents to the programme. The segment control bytes are however still transferred, (presumably because the FORTRAN-77 I/O system is not being used). Additionally, the DOS side of the transfer assumes that the data are coming from a UNIX system, so every time it encounters an incoming <LF> (0x0A), it adds a <CR> (0x0D) to conform to DOS line delimiters. This obviously screws up the data, to the point that recovery, while possible, is a tedious programming exercise requiring knowledge of the original record lengths (in principle calculable).
If one ftps a PSI data file in binary mode, then the VAX side of the transfer apparently opens the file in such a way that it is interpreted as a pure binary file and all information, including record-length bytes and control bytes, are transferred. No information is added or deleted (so one receives the exact data one sees if one does a binary DUMP of the file under VMS). Thus one needs a filter programme to read in the record length information, then read that many bytes, discard the control bytes and write the rest to an output file, continuing until the end-of-file is achieved.

Such a programme, called PSICVT, has been written in Turbo C, and is available in the usual places.When ftping to the PC, make sure you use BINARY mode! (usage of other transfer methods such as Kermit and Xmodem will have similar results). Then one can use PSICVT to reconstruct a flat binary image of the original data. In general, non-VAX transferral methods are to be avoided. You should restrict yourself to DECNet, magnetic tape, or the jNet SEND/FILE/VMSDUMP command to transfer PSI-format files between VAXen, to ensure that the record structure is maintained. This is not a problem with TRIUMF format files, as they are true binary files and may be transferred at will, so long as a binary switch is activated when available.

File translated from T_EX by T_TH, version 2.67.
On 21 Aug 2000, 11:11.