diff --git a/doc/memos/rrf/08011-Fourier-Averaged.pdf b/doc/memos/rrf/08011-Fourier-Averaged.pdf new file mode 100644 index 00000000..c95586d5 Binary files /dev/null and b/doc/memos/rrf/08011-Fourier-Averaged.pdf differ diff --git a/doc/memos/rrf/08011-RRF-Asym-Fourier-Averaged-Ghost.pdf b/doc/memos/rrf/08011-RRF-Asym-Fourier-Averaged-Ghost.pdf new file mode 100644 index 00000000..85876e84 Binary files /dev/null and b/doc/memos/rrf/08011-RRF-Asym-Fourier-Averaged-Ghost.pdf differ diff --git a/doc/memos/rrf/08011-RRF-Histo-Fourier-Averaged-Ghost.pdf b/doc/memos/rrf/08011-RRF-Histo-Fourier-Averaged-Ghost.pdf new file mode 100644 index 00000000..74eb41cc Binary files /dev/null and b/doc/memos/rrf/08011-RRF-Histo-Fourier-Averaged-Ghost.pdf differ diff --git a/doc/memos/rrf/08019-Fourier-Averaged.pdf b/doc/memos/rrf/08019-Fourier-Averaged.pdf new file mode 100644 index 00000000..baee392f Binary files /dev/null and b/doc/memos/rrf/08019-Fourier-Averaged.pdf differ diff --git a/doc/memos/rrf/HAL-9500-t0.pdf b/doc/memos/rrf/HAL-9500-t0.pdf new file mode 100644 index 00000000..90fe700d Binary files /dev/null and b/doc/memos/rrf/HAL-9500-t0.pdf differ diff --git a/doc/memos/rrf/PSI_Logo_narrow_blau.jpg b/doc/memos/rrf/PSI_Logo_narrow_blau.jpg new file mode 100644 index 00000000..aa253979 Binary files /dev/null and b/doc/memos/rrf/PSI_Logo_narrow_blau.jpg differ diff --git a/doc/memos/rrf/PSI_Logo_wide_blau.pdf b/doc/memos/rrf/PSI_Logo_wide_blau.pdf new file mode 100644 index 00000000..13482e93 Binary files /dev/null and b/doc/memos/rrf/PSI_Logo_wide_blau.pdf differ diff --git a/doc/memos/rrf/rrf-notes.pdf b/doc/memos/rrf/rrf-notes.pdf new file mode 100644 index 00000000..89b20118 Binary files /dev/null and b/doc/memos/rrf/rrf-notes.pdf differ diff --git a/doc/memos/rrf/rrf-notes.tex b/doc/memos/rrf/rrf-notes.tex new file mode 100644 index 00000000..74dd9903 --- /dev/null +++ b/doc/memos/rrf/rrf-notes.tex @@ -0,0 +1,265 @@ +\documentclass[twoside]{article} + +\usepackage[english]{babel} +\usepackage{a4} +\usepackage{amssymb} +\usepackage{graphicx} +\usepackage{fancyhdr} +\usepackage{array} +\usepackage{float} +\usepackage{hyperref} +\usepackage{xspace} +\usepackage{rotating} +\usepackage{dcolumn} +\usepackage{url} + +\setlength{\topmargin}{10mm} +\setlength{\topmargin}{-13mm} +% \setlength{\oddsidemargin}{0.5cm} +% \setlength{\evensidemargin}{0cm} +\setlength{\oddsidemargin}{1cm} +\setlength{\evensidemargin}{1cm} +\setlength{\textwidth}{14.5cm} +\setlength{\textheight}{23.8cm} + +\def\mathbi#1{\textbf{\em #1}} + +\pagestyle{fancyplain} +\addtolength{\headwidth}{0.6cm} +\fancyhead{}% +\fancyhead[RE,LO]{\bf RRF Fits}% +\fancyhead[LE,RO]{\thepage} +\cfoot{--- Andreas Suter -- \today ---} +\rfoot{\includegraphics[width=6.4cm]{PSI_Logo_wide_blau.pdf}} + +\DeclareMathAlphabet{\bi}{OML}{cmm}{b}{it} + +\newcommand{\mean}[1]{\langle #1 \rangle} +\newcommand{\lem}{LE-$\mu$SR\xspace} +\newcommand{\musr}{$\mu$SR\xspace} +\newcommand{\etal}{\emph{et al.\xspace}} +\newcommand{\ie}{\emph{i.e.\xspace}} +\newcommand{\eg}{\emph{e.g.\xspace}} + +\newcolumntype{d}[1]{D{.}{.}{#1}} + +\begin{document} + +% Header info -------------------------------------------------- +\thispagestyle{empty} +\noindent +\begin{tabular}{@{\hspace{-0.7cm}}l@{\hspace{6cm}}r} +\noindent\includegraphics[width=3.4cm]{PSI_Logo_narrow_blau.jpg} & + {\Huge\sf Memorandum} +\end{tabular} +% +\vskip 1cm +% +\begin{tabular}{@{\hspace{-0.5cm}}ll@{\hspace{4cm}}ll} +Datum: & \today & & \\[3ex] +Von: & Andreas Suter & An: & \\ +Telefon: & +41\, (0)56\, 310\, 4238 & & \\ +Raum: & WLGA / 119 & cc: & \\ +e-mail: & \verb?andreas.suter@psi.ch? & & \\ +\end{tabular} +% +\vskip 0.3cm +\noindent\hrulefill +\vskip 1cm +% +\section{Rotating Reference Frame Fits} + +High transverse field \musr (HTF-\musr) experiments will typically lead to +rather large data sets since it is necessary to follow the high frequencies +present in the positron decay histograms. Currently the HAL-9500 instrument +at PSI \cite{HAL9500} is operated with $2\time 8$ positron detector, with a +typical number of $\sim 4\times 10^5$ histogram bins. To analyze HTF-\musr +data on rather slugish computer hardware is a challenge. In the last millennium +the people invented the rotating reference frame transformation (RRF) \cite{Riseman90} +to reduce to data sets to be handled. + +Here I will shortly describe the ways how it is implemented in \textsc{Musrfit}, +and why it should be avoided to be used altogether. The starting point of all +is given by the positron decay spectrum which formally takes the form + +\begin{equation}\label{eq:positron_decay_spectrum} + N^{(j)}(t) = N_0^{(j)} \exp(-t / \tau_\mu) \, \left[ 1 + A^{(j)}(t) \right] + N_{\rm bkg}^{(j)}, +\end{equation} + +\noindent where $(j)$ is the index of the positron counter, $N_0$ gives the scale +of recorded positrons, $\tau_\mu$ is the muon lifetime, $A(t)$ the asymmetry, and $N_{\rm bkg}$ +describes the background due to uncorrelated events. + +The idea behind the RRF is twofolded: (i) try to extract $A(t)$, and (ii) shift the high frequency +data set $A(t)$ to lower frequencies such that the number of necessary bins needed can be +reduced (packing / rebinning), and hence the overall number of bins is much smaller. + +As I will try to explain, this is not for free, and there are problems arising from +this kind of data treatment. + +\subsection{Single Histogram RRF Implementation} + +In a first step the asymmetry needs to be determined. This is done the following way: +\begin{enumerate} + \item Determine the background, $N_{\rm bkg}$, at times before $t_0$ ($t_0$ is the time of the muon + implantation). Hopefully the background before and after $t_0$ is equal, which is not always + the case. + \item Multiple the background corrected histogram with $\exp(+t / \tau_\mu)$, this is leading to +\end{enumerate} + +\begin{equation}\label{eq:Mt} + M(t) \equiv \left[ N(t) - N_{\rm bkg} \right]\,\exp(+t / \tau_\mu) = N_0 \left[ 1 + A(t) \right]. +\end{equation} + +\begin{enumerate} + \setcounter{enumi}{2} + \item In order to extract $A(t)$ from $M(t)$, $N_0$ needs to be determined, which is almost the + most tricky part here. The idea is simple: since $A(t)$ is dominated by high frequency signals, + proper averaging over $M(t)$ should allow to determine $N_0$, assuming that $\langle A(t) \rangle = 0$. + Is this assumption always true? \emph{No!} For instance it is \emph{not} true if the averaging is preformed + over incomplete periodes of a single assumed signal. Another case where it will fail is if multiple signals + with too far apart frequencies is present, as \eg in the case of muonium. Said all this, let's come + back and try to determine $N_0$: +\end{enumerate} + +\begin{equation}\label{eq:N0estimate} + N_0 = \sum_{k=0}^{N_{\rm avg}} w_k M(t_k), +\end{equation} + +\noindent where $N_{\rm avg}$ is determined such that $N_{\rm avg} \Delta t \simeq 1 \mu$s ($\Delta t$ being the +time resolution. $1 \mu$s means averaging over many cyles). In order to get a good estimate for $N_{\rm avg}$, $N(t)$ is Fourier +transformed, and from this power spectrum the frequency with the largest amplitude is determined, $\nu_{\rm 0}$. From $\nu_0$, +$\Delta t$, the number of cycles fitting into $1 \mu$s can be determined, and from this and the time resolution $N_{\rm avg}$ +can be calculated. The weight $w_k$ is given by: + +\begin{equation} + w_k = \frac{\left[\Delta M(t_k)\right]^{-2}}{\sum_{j=0}^{N_{\rm avg}} \left[\Delta M(t_j)\right]^{-2}}, +\end{equation} + +\noindent where + +\begin{equation} + \Delta M(t) = \left[ \left(\frac{\partial M}{\partial N} \Delta N\right)^2 + + \left(\frac{\partial M}{\partial N_{\rm bkg}} \Delta N_{\rm bkg}\right)^2 \right]^{1/2} + \simeq \exp(+t/\tau_\mu) \sqrt{N(t)}. +\end{equation} + +\noindent The error estimate on $N_0$ is then + +\begin{equation}\label{eq:N0-rrf} + \Delta N_0 = \sigma_{N_0} = \sqrt{\sum_k w_k^2 \Delta M(t_k)^2}. +\end{equation} + +\noindent Having estimated $N_0$, the asymmetry can be extracted as: + +\begin{equation} + A(t) = M(t) / N_0 - 1. +\end{equation} + +\begin{enumerate} + \setcounter{enumi}{3} + \item Now the actual RRF transformation can take place: $A_{\rm rrf}(t) = 2 \times A(t) \cos(\omega_{\rm rrf} t + \phi_{\rm rrf})$. + The factor of 2 is introduced to conserve the asymmetry amplitude. The idea is the following: Fourier transform theory tells as that +\end{enumerate} + +\begin{equation} + {\cal F}\left\{2 \times A(t) \cos(\omega_{\rm rrf} t + \phi_{\rm rrf}) \right\} = + {\cal F}\left\{A(t)\right\}(\omega-\omega_{\rm rrf}) + {\cal F}\left\{A(t)\right\}(\omega+\omega_{\rm rrf}), +\end{equation} + +\noindent \ie that the Fourier spectrum of $A(t)$ is shifted down and up by $\omega - \omega_{\rm rrf}$ and $\omega + \omega_{\rm rrf}$, respectively. +In order to get rid of the high frequency part (${\cal F}\left\{A(t)\right\}(\omega+\omega_{\rm rrf})$), $A_{\rm rrf}(t)$ will be heavily over-binned, +\ie + +\begin{enumerate} + \setcounter{enumi}{4} + \item Do the rrf packing: $A_{\rm rrf}(t) \rightarrow \langle A_{\rm rrf}(t) \rangle_p$. Packing itself is a filtering of data! Especially this + kind of filter is dispersive \cite{King89}, \ie that it potentially is leading to line shape distortions. For symmetric, rather narrow + lines, this is unlikely to be a problem. However, this might be quite different for complex line shapes as in the case of vortex lattices. +\end{enumerate} + +\noindent The property $\langle A_{\rm rrf}(t) \rangle_p$ is what is fitted. The error on this property is estimated the following way: (i) the +unpacked error of $A(t)$ is: + +\begin{equation} + \Delta A(t) \simeq \frac{\exp(+t/\tau_\mu)}{N_0}\,\left[ N(t) + \left( \frac{N(t)-N_{\rm bkg}}{N_0} \right)^2 \Delta N_0^2 \right]^{1/2}, +\end{equation} + +\noindent and form this the packed $A_{\rm rrf}(t)$ error can be calculated. + +\subsection{Asymmetry RRF Implementation} + +\begin{enumerate} + \item In order to circumvent the difficulties to estimate $N_0$ the asymmetry of the starting positron histograms is formed. For details + see \cite{musrfit_userManual15}. For this, positron detectors geometrically under $180^\circ$ are used. However, due the the spiraling + of the positron in sufficiently high magnetic fields, and the uncertainties of the $t_0$'s, the geometrical phase might \emph{not} + correspond to the positron signal phase! At $B=9$T the uncertainty in $t_0$ by one channel leads to a phase shift of + $\gamma_\mu B \Delta t \cdot (180 / \pi) = 1.7^\circ$. Fig.\ref{fig:hal-9500-t0} shows the $t_0$-region of a typical HAL-9500 spectrum. + It shows that it is very hard to get the absolut value of $t_0$ right. + \item Carry out the RRF transformation $A_{\rm rrf}(t) = 2 \times A(t) \cos(\omega_{\rm rrf} t + \phi_{\rm rrf})$. + \item Do the rrf packing. +\end{enumerate} + +\begin{figure}[h] + \centering + \includegraphics[width=0.5\textwidth]{HAL-9500-t0.pdf} + \caption{The $t_0$ region of a typical HAL-9500 spectrum. The broad black hump with the green line, is the ``prompt'' peak. + It is \emph{not} straight forward how to define $t_0$.}\label{fig:hal-9500-t0} +\end{figure} + + +\section{Discussion} + +Both RRF transformation sketched above have weak points which makes it hard to estimate systematic errors. Both methods will fail at too +low fields of $\lesssim 1$T. The only and single purpose of the RRF transformation is slughish computer power! We developed GPU based fitting +which overcomes \emph{all} this uncontrolled weaknesses and henceforth RRF could be omitted altogether. I still added it for the time being, +since strong GPU/CPU hardware is still a bit costly and therefore not affordable to everyone. + +In order to give a feeling about what might go ``wrong'' with the RRF, I was running a couple of test cases. The chosen asymmetry is + +\begin{equation}\label{eq:asymmetry-simulation} + A^{(j)}(t) = A_0^{(j)} \sum_{k=1}^3 f_k \exp\left[-0.5\cdot (\sigma_k t)^2 \right] \cos(\gamma_\mu B_k t + \phi^{(j)}), +\end{equation} + +\noindent with values found in Tab.\ref{tab:asymmetry-parameters}. For the simulation 4 positron detector signals were generated with +$A_0^{(j)} = \left\{ 0.2554, 0.2574, 0.2576, 0.2566 \right\}$. The further ingredients were: $N_0^{(j)} = \left\{ 27.0, 25.3, 25.6, 26.9 \right\}$, +$N_{\rm bkg}^{(j)} = \left\{ 0.055, 0.060, 0.069, 0.064 \right\}$, and $ \phi^{(j)} = \left\{ 5.0, 95.0, 185.0, 275.0 \right\}$. + +\begin{table}[h] + \centering + \begin{tabular}{c|c|c|c} + $k$ & $f_k$ & $\sigma_k$ & $B_k$ \\ + & & ($1/\mu$s) & (T) \\ \hline\hline + 1 & 0.5 & 7 & $1$ or $9$ \\ + 2 & 0.2 & 0.75 & $1.02$ or $9.02$ \\ + 3 & 0.3 & 0.25 & $1.06$ or $9.06$ + \end{tabular} + \caption{Parameters used in Eq.(\ref{eq:asymmetry-simulation}).}\label{tab:asymmetry-parameters} +\end{table} + +\begin{figure}[h] + \centering + \includegraphics[width=0.45\textwidth]{08011-Fourier-Averaged.pdf} \quad + \includegraphics[width=0.45\textwidth]{08011-RRF-Histo-Fourier-Averaged-Ghost.pdf} \\ + \includegraphics[width=0.45\textwidth]{08011-RRF-Asym-Fourier-Averaged-Ghost.pdf} + \caption{Averaged Fourier power spectra. Top left: from single histogram fit for the 1T data set. + Top right: from single histogram RRF fit. Bottom: from asymmetry RRF fit. Both RRF + sets show ghost lines.}\label{fig:08011-Fourier} +\end{figure} + +Figure.\ref{fig:08011-Fourier} shows the averaged Fourier power spectra for the simulated data sets at 1T. Both RRF transformation are +showing ghost lines, even for optimally chosen RRF rebinning. At higher fields this is less pronounced. The ghost lines have various origins +such as aliasing effects due to the RRF packing not prefectly suppressing the high frequency part of $A_{\rm rrf}(t)$, leakage of the RRF frequency +for not sufficently precise known $N_0$ (see Eq.(\ref{eq:N0-rrf})) for single histogram RRF fits, etc. + +Fits of simulated data as described above (see Eq.(\ref{eq:asymmetry-simulation}), with fields 0.5, 1.0, 3.0, 5.0, 7.0, and 9.0T) show that above +about 1T the model parameters of the RRF fits are acceptable, but the error bars are typically about a factor 3 larger compared to single histogram +fits. The asymmetries of the RRF fits are ``substantially'' too small. The $\chi^2$ values are close to meaningless for the RRF fits, since they +are strongly depending on the RRF packing, time interval chosen, etc. + +To summaries: RRF fits can be used for online analysis if no GPU accelerator is available, but \emph{must not} be used for any final analysis! + +\bibliographystyle{amsplain} +\bibliography{rrf} + +\end{document} \ No newline at end of file diff --git a/doc/memos/rrf/rrf.bib b/doc/memos/rrf/rrf.bib new file mode 100644 index 00000000..6106f488 --- /dev/null +++ b/doc/memos/rrf/rrf.bib @@ -0,0 +1,29 @@ +@Article{ Riseman90, + title = {{The Rotating Reference Frame Transformation in $\mu$SR}}, + author = "T. M. Riseman and J. H. Brewer", + journal = "Hyperfine Interact.", + volume = "65", + year = "1990", + pages = "1107" +} + +@Misc{ HAL9500, + title = {{HAL-9500}}, + author = "", + note = {{\verb!https://www.psi.ch/smus/hal-9500!}}, +} + +@book{ King89, + Address = {New York}, + Author = {R. King and M. Ahmadi and R. Gorgui-Naguib and A. Kwabwe and M. Azimi-Sajadi}, + Edition = {1st}, + Publisher = {Plenum Press}, + Title = {Digital Filtering in one and two Dimensions}, + Year = 1989 +} + +@Misc{ musrfit_userManual15, + title = {{Musrfit User Guide}}, + author = "", + note = {{\verb!https://intranet.psi.ch/MUSR/MusrFit#A_5.3_Asymmetry_Fit!}}, +}