Offset in Asymmetry

[AnthonyLim23]

Expected behavior

MUSR data should osciallte around a value of zero for the asymmetry.

Actual behavior

There is a non-zero offset in the data

Steps to reproduce the behavior

Load a file (e.g. MUSR00062260.nxs) via the muon analysis. The data will be offset. This is easily visable if you fit a flat backgound and exposc functions to the data. The flat background will have a non-zero value.

Platforms affected

All

[AnthonyLim23]

The offset is numeric only when the flatBackground and the GausOsc functions are used together. If just the flatBackground is fitted to the data then it returns a value of zero (except it is visible that there is an offset). Further investigation has shown that this is due to the error bars being used in the fitting routine. Fitting a flatBackground without error bars improves the answer but the value is still larger than the error of the value. This has been tested with the below code:

#!/usr/bin/env python
#

from math import *
import matplotlib.pylab as plt
import numpy as np
from scipy.optimize import curve_fit

def flat(x,C):
    return C

X,Y,E=np.loadtxt('musr.dat',delimiter=',',unpack=True)
popt,pcov = curve_fit(flat,X,Y)
print "without errors ", popt[0]

popt,pcov = curve_fit(flat,X,Y,sigma=E)
print "with errors    ", popt[0]

plt.plot(X,E)
#plt.plot(X,Y,'k')
#plt.errorbar(X,Y,yerr=E)
plt.show()

[AnthonyLim23]

Fitting a GausOsc + flatBackground during the normalization (use the value for the flatBackground as the normalization constant) does give a value for zero for the intercept, in most cases. However, the method fails for 0 field e.g. 57980 (musr). Using just a flatBackground to determine the normalization constant does remove it if the scenario is the same (e.g. no errorbars are used in both the normalization and fit of the output data). This suggests that the problem is that the method we are using is inappropriate for the data, as the current methods do what they are supposed to.

[AnthonyLim23]

Talked to the scientists and got the following request:

On loading / grouping the data, generate the "normalised counts" by multiplying by exp(t/muon lifetime) and dividing by the number of good frames (perhaps also by the number of detectors in the group). Probably store this in another workspace (called something like MUSR01234__NormCounts ).

Do an initial guess of "N0" by summing the raw counts in the group, over the time range to be used for fitting, and dividing by the integral of the muon lifetime exponential over that range.

To plot the data, generate an Asymmetry workspace with Asymmetry = (NormCounts/N0)-1. Perhaps use a child algorithm for this as it's needed again later.

For fitting: have a composite function called something like "TFAsymmetry" which takes a numeric parameter N0 and also sums the user's choice of fit function(s) included within it, just as the implied CompositeFunction does now. The collection of muon and other functions that might be included within do not need to be changed at all. This function would be automatically entered in the function dialog, not to be deleted, and its parameter N0 initialised to the guess above.

Alternatively TFAsymmetry could be "hidden" from the user with only a line for N0 appearing somewhere. This might make it easier to add/remove it if the user changes the fitted group to one where it is no longer appropriate (see below).

TFAsymmetry calculates N0 * (1 + f0(x,params...) + f1(x,params...) + ...)

The fitting is done to the workspace NormCounts so that workspace acquires spectra "Calc" and "Diff".

Then re-run the asymmetry calculation child algorithm, this time using the N0 resulting from the fit. All three spectra get the same treatment. Plot the Asymmetry workspace (data, calc, diff).

Should N0 be added to the workspace NormCounts as a log value, and updated following the fit, so it need not be carried around elsewhere in the interface and given as an explicit parameter to the child algorithm?

Steve mentions - we need to treat "forward-backward" asymmetry data just as it is now and therefore there may need to be some way for the fitting section to know what sort of group it has been presented with, and how to treat it. It's possible to change the selected workspace (group) once in the Fitting tab, perhaps between one type and the other, and this should not cause problems. Perhaps use the naming scheme (there's "Group" or "Pair" in the name now)?

Another suggestion is to duplicate the Data Analysis tab and have one for single group (TF) data and one for paired F-B asymmetry...

One could alternatively have NormCounts be the actual counts (not lifetime corrected), and the composite function TFAsymmetry and the child algorithm to get asymmetry do that correction too. That might make a subtle difference when fitting the "zero" bins at the end of a low stats spectrum, or make it easier to include a Background term, constant counts per unit time.

[AnthonyLim23]

I have spoken to James Lord and the solution that the scientist would like is below:

The dead-time corrected and grouped raw data is denoted by W, such that W_i is then ith entry in the data and has a corresponding time t_i.

Then the normalized counts, C_{norm} is given by

C_{norm} = W*exp(t/tau)/f, [1]

where tau is the muon lifetime and f is the number of good frames. The first estimate at the normalization constant, N_0^{(0)} is given by

N0^{(0)} = Delta sum{i : a \le t_i \le b} W_i / (\int_a^b exp(-t/tau) dt), [2]

where a is the earliest time in the analysis and b is the latest, Delta is the size of the time step. The bottom of the summation "i : a \le t_i \le b" means "i such that a <=t_i<=b". Equation 2 can be simplified to

N0^{(0)} = Delta sum{i : a \le t_i \le b} W_i /(tau* [exp(-a/tau) - exp(-b/tau) ]). [3]

The first attempt at calculating the asymmetry, A_{sym}^{(0)} is then given by

A{sym}^{(0)} = ( C{norm} / N_0^{(0)} ) -1. [4]

For more accurate asymmetry, A_(sym)^{(1)} the following additional steps are required. The equation

C_{norm} \approx N0^{(1)} * [ 1 + \sum{j=1}^m f_j(t,{\lambda_j})], [5]

where N0^{(1)} is a more accurate normalization constant, M is the total number of additional functions in the fit to C{norm}, f_j is the jth additional function and {\lambda_j} is the set of parameters for f_j. Then the new asymmetry is given by

A{sym}^{(1)} = ( C{norm} / N_0^{(1)} ) -1. [6]