Reduce amor new files

[jokasimr]

Fixes #38 Fixes #41

~This should not be merged in its current state.~ ~The workflow does not give the same result as the program Jochen has.~ ~Some parts of the notebook should be removed and placed elsewhere in the project, such as the geometry definition.~

Edit:

The workflow still does not match the PSI software results exactly, but at this point it is close enough, and we will fix the rest in a follow-up

[jokasimr]

Here's a description of the basic idea of the procedure. There are two different ways to implement the reference correction, I'm not sure which one is the best to use. Possibly the best is just to use what Jochen is using now, and I think he's using the option that is not implemented in the PR. So this might be something we want to change.

Model the intensity obtained from the sample measurement as $$I^{sample}(z, \lambda)= I^{reference}(z, \lambda) R(Q(z, \lambda))$$ where $I^{reference}$ is the intensity obtained when $R(Q)=1$, $z$ is here not the coordinate in the global coordinate system, but instead the logical $z$-coordinate of the detector.

Two options for how to compute $R(Q)$

1. I think this is what Jochen does. Using $$R(Q(z, \lambda)) = \frac{I^{sample}(z, \lambda)}{I^{reference}(z, \lambda)}$$ then for every $(z, \lambda)$ we have a very noisy measurement of $R(Q(z, \lambda))$. To obtain a good estimate for $R(Q)$ we must combine measurements of $R(Q(z, \lambda))$ in a way that takes the noise into consideration. Directly histogramming $R(Q(z, \lambda))$ in $Q$ gives a bad result because we give as much weight to measurements with lots of noise as we do to measurements with little noise.
   • Possible solution: Scale by the variances of the measurements, this is the maximum likelihood estimate of $R(Q)$ assuming normal distributed noise. Specifically: $$R(Qi) \approx \frac{\int{Q(z, \lambda) \in Qi} \frac{R(Q(z, \lambda))}{Var[R(Q(z, \lambda))]} \ dz \ d\lambda}{\int{Q(z, \lambda) \in Q_i} \frac{1}{Var[R(Q(z, \lambda))]} \ dz \ d\lambda}$$
   • When might the above not work?
     1. If the assumption of normal distributed noise fails, that is, if the number of counts is low.
     2. If the variance estimate is not good.
     3. This is the approach currently used in the workflow implementation. Assuming $R(Q)$ varies slowly compared to the coarseness of the binning: $$I(Qi) = \int{Q(z, \lambda) \in Q_i} I^{sample}(z, \lambda) \ dz \ d\lambda \approx R(Qi) \int{Q(z, \lambda) \in Q_i} I^{reference}(z, \lambda) \ dz \ d\lambda $$ $$R(Qi) \approx \frac{ \int{Q(z, \lambda) \in Qi} I^{sample}(z, \lambda) \ dz \ d\lambda } { \int{Q(z, \lambda) \in Q_i} I^{reference}(z, \lambda) \ dz \ d\lambda } $$

[jokasimr]

Update: Comparison between method 1 and method 2 mentioned above.

The results seem quite similar, so probably okay to leave it as is right now.

Implementation of method 1, for the record

da = pl.compute(FootprintCorrectedData[Sample])\
    .group('z_index').bin(wavelength=pl.compute(WBins))
n = pl.compute(NormalizationFactor)
R = da / sc.values(n)
q = n.coords['Q']

def compute_ioq(R, q, relaxation):
    vR = R.copy(deep=True)
    vR.data = sc.variances(vR.data)

    A = sc.ones(dims=('z_index',), shape=(14*32,)) * (pl.compute(WBins)[1:] - pl.compute(WBins)[:-1])
    # Add small relaxation factor to variances, increases noise in regions with few counts, but makes the inverse of the variances more well behaved
    vR = A / (vR.hist() + sc.scalar(relaxation))
    vR.coords['Q'] = q

    return (R * vR).bins.concat().hist(Q=pl.compute(QBins)) / vR.flatten(to='Q').hist(Q=pl.compute(QBins))

f = sc.plot({
    'method1(0.001)': compute_ioq(R, q, 0.001),
    'method1(0.01)': compute_ioq(R, q, 0.01),
    'method1(0.1)': compute_ioq(R, q, 0.1),
    'method1(0.5)': compute_ioq(R, q, 0.5),
    'method2': pl.compute(NormalizedIofQ).hist(),
}, norm='log')
f

[SimonHeybrock]

The results seem quite similar, so probably okay to leave it as is right now.

There appear to be significant differences for some Q. Is it the scientist's judgement that this is irrelevant?

[jokasimr]

Is it the scientist's judgement that this is irrelevant?

No it is not.

[jokasimr]

Comparison between Jochens workflow and essreflectometry for all data files:

~The shift in $Q$ is likely due to a tof correction that is missing in essreflectometry.~ Probably done, might be another smaller correction left. Other differences might be due to differences in masking. Other differences might be due to bugs.

[jokasimr]

There are some things that I'm aware needs improvement, mainly, docstrings and tests.

[nvaytet]

@SimonHeybrock I think we would like you to check that all your comments have been addressed, but it is not required for you to go through everything again.

[nvaytet]

I added some 'reasonable' x positions for now. We are ignoring x in the coordinate transform, so I think we are good now :+1:

[jokasimr]

I added some 'reasonable' x positions for now. We are ignoring x in the coordinate transform, so I think we are good now 👍

Not sure what you mean by we're ignoring the x coord. It's used in the lambda computation.

Adding some reasonable dummy values is fine for me. I'll add an issue for figuring out the real values.

[nvaytet]

Not sure what you mean by we're ignoring the x coord. It's used in the lambda computation.

Ah it's ignored in the theta computation, but not in the wavelength computation.