Merge branch 'devel-rrf'

This commit is contained in:
suter_a 2016-01-22 15:25:53 +01:00
commit 55f70cd526
70 changed files with 4252 additions and 176 deletions

Binary file not shown.

Binary file not shown.

Binary file not shown.

Binary file not shown.

After

Width:  |  Height:  |  Size: 107 KiB

Binary file not shown.

BIN
doc/memos/rrf/rrf-notes.pdf Normal file

Binary file not shown.

265
doc/memos/rrf/rrf-notes.tex Normal file
View File

@ -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}

29
doc/memos/rrf/rrf.bib Normal file
View File

@ -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!}},
}

View File

@ -15,12 +15,14 @@ h_sources = \
../include/PMusrT0.h \
../include/PPrepFourier.h \
../include/PRunAsymmetry.h \
../include/PRunAsymmetryRRF.h \
../include/PRunBase.h \
../include/PRunDataHandler.h \
../include/PRunListCollection.h \
../include/PRunNonMusr.h \
../include/PRunMuMinus.h \
../include/PRunSingleHisto.h \
../include/PRunSingleHistoRRF.h \
../include/PStartupHandler.h \
../include/PTheory.h
@ -59,12 +61,14 @@ cpp_sources = \
PMusrT0.cpp \
PPrepFourier.cpp \
PRunAsymmetry.cpp \
PRunAsymmetryRRF.cpp \
PRunBase.cpp \
PRunDataHandler.cpp \
PRunListCollection.cpp \
PRunNonMusr.cpp \
PRunMuMinus.cpp \
PRunSingleHisto.cpp \
PRunSingleHistoRRF.cpp \
PStartupHandler.cpp \
PTheory.cpp

View File

@ -8,7 +8,7 @@
***************************************************************************/
/***************************************************************************
* Copyright (C) 2007-2014 by Andreas Suter *
* Copyright (C) 2007-2016 by Andreas Suter *
* andreas.suter@psi.ch *
* *
* This program is free software; you can redistribute it and/or modify *

View File

@ -8,7 +8,7 @@
***************************************************************************/
/***************************************************************************
* Copyright (C) 2007-2014 by Andreas Suter *
* Copyright (C) 2007-2016 by Andreas Suter *
* andreas.suter@psi.ch *
* *
* This program is free software; you can redistribute it and/or modify *
@ -77,12 +77,16 @@ Double_t PFitterFcn::operator()(const std::vector<Double_t>& par) const
if (fUseChi2) { // chi square
value += fRunListCollection->GetSingleHistoChisq(par);
value += fRunListCollection->GetSingleHistoRRFChisq(par);
value += fRunListCollection->GetAsymmetryChisq(par);
value += fRunListCollection->GetAsymmetryRRFChisq(par);
value += fRunListCollection->GetMuMinusChisq(par);
value += fRunListCollection->GetNonMusrChisq(par);
} else { // max likelihood
value += fRunListCollection->GetSingleHistoMaximumLikelihood(par);
value += fRunListCollection->GetSingleHistoRRFMaximumLikelihood(par);
value += fRunListCollection->GetAsymmetryMaximumLikelihood(par);
value += fRunListCollection->GetAsymmetryRRFMaximumLikelihood(par);
value += fRunListCollection->GetMuMinusMaximumLikelihood(par);
value += fRunListCollection->GetNonMusrMaximumLikelihood(par);
}

View File

@ -8,7 +8,7 @@
***************************************************************************/
/***************************************************************************
* Copyright (C) 2007-2015 by Andreas Suter *
* Copyright (C) 2007-2016 by Andreas Suter *
* andreas.suter@psi.ch *
* *
* This program is free software; you can redistribute it and/or modify *

View File

@ -8,7 +8,7 @@
***************************************************************************/
/***************************************************************************
* Copyright (C) 2007-2015 by Andreas Suter *
* Copyright (C) 2007-2016 by Andreas Suter *
* andreas.suter@psi.ch *
* *
* This program is free software; you can redistribute it and/or modify *

View File

@ -8,7 +8,7 @@
***************************************************************************/
/***************************************************************************
* Copyright (C) 2007-2014 by Andreas Suter *
* Copyright (C) 2007-2016 by Andreas Suter *
* andreas.suter@psi.ch *
* *
* This program is free software; you can redistribute it and/or modify *

View File

@ -8,7 +8,7 @@
***************************************************************************/
/***************************************************************************
* Copyright (C) 2007-2014 by Andreas Suter *
* Copyright (C) 2007-2016 by Andreas Suter *
* andreas.suter@psi.ch *
* *
* This program is free software; you can redistribute it and/or modify *

View File

@ -8,7 +8,7 @@
***************************************************************************/
/***************************************************************************
* Copyright (C) 2009-2015 by Bastian M. Wojek / Andreas Suter *
* Copyright (C) 2009-2016 by Bastian M. Wojek / Andreas Suter *
* andreas.suter@psi.ch *
* *
* This program is free software; you can redistribute it and/or modify *

View File

@ -8,7 +8,7 @@
***************************************************************************/
/***************************************************************************
* Copyright (C) 2007-2015 by Andreas Suter *
* Copyright (C) 2007-2016 by Andreas Suter *
* andreas.suter@psi.ch *
* *
* This program is free software; you can redistribute it and/or modify *
@ -311,6 +311,11 @@ Int_t PMsrHandler::ReadMsrFile()
if (!CheckAddRunParameters())
result = PMUSR_MSR_SYNTAX_ERROR;
// check that if RRF settings are present, the RUN block settings do correspond
if (result == PMUSR_SUCCESS)
if (!CheckRRFSettings())
result = PMUSR_MSR_SYNTAX_ERROR;
if (result == PMUSR_SUCCESS) {
CheckMaxLikelihood(); // check if the user wants to use max likelihood with asymmetry/non-muSR fit (which is not implemented)
CheckLegacyLifetimecorrection(); // check if lifetimecorrection is found in RUN blocks, if yes transfer it to PLOT blocks
@ -348,7 +353,7 @@ Int_t PMsrHandler::WriteMsrLogFile(const Bool_t messages)
{
const UInt_t prec = 6; // default output precision for float/doubles
UInt_t neededPrec = 0;
UInt_t neededWidth = 9;
UInt_t neededWidth = 9;
Int_t tag, lineNo = 0, number;
Int_t runNo = -1, addRunNo = 0;
@ -612,9 +617,15 @@ Int_t PMsrHandler::WriteMsrLogFile(const Bool_t messages)
case MSR_FITTYPE_SINGLE_HISTO:
fout << left << "fittype" << MSR_FITTYPE_SINGLE_HISTO << " (single histogram fit)" << endl;
break;
case MSR_FITTYPE_SINGLE_HISTO_RRF:
fout << left << "fittype" << MSR_FITTYPE_SINGLE_HISTO_RRF << " (single histogram RRF fit)" << endl;
break;
case MSR_FITTYPE_ASYM:
fout << left << "fittype" << MSR_FITTYPE_ASYM << " (asymmetry fit)" << endl ;
break;
case MSR_FITTYPE_ASYM_RRF:
fout << left << "fittype" << MSR_FITTYPE_ASYM_RRF << " (asymmetry RRF fit)" << endl ;
break;
case MSR_FITTYPE_MU_MINUS:
fout << left << "fittype" << MSR_FITTYPE_MU_MINUS << " (mu minus fit)" << endl ;
break;
@ -624,6 +635,27 @@ Int_t PMsrHandler::WriteMsrLogFile(const Bool_t messages)
default:
break;
}
} else if (sstr.BeginsWith("rrf_freq", TString::kIgnoreCase) && (fGlobal.GetFitType() == MSR_FITTYPE_SINGLE_HISTO_RRF)) {
fout.width(16);
fout << left << "rrf_freq ";
fout.width(8);
neededPrec = LastSignificant(fGlobal.GetRRFFreq(fGlobal.GetRRFUnit().Data()),10);
fout.precision(neededPrec);
fout << left << fixed << fGlobal.GetRRFFreq(fGlobal.GetRRFUnit().Data());
fout << " " << fGlobal.GetRRFUnit();
fout << endl;
} else if (sstr.BeginsWith("rrf_phase", TString::kIgnoreCase) && (fGlobal.GetFitType() == MSR_FITTYPE_SINGLE_HISTO_RRF)) {
fout.width(16);
fout << "rrf_phase ";
fout.width(8);
fout << left << fGlobal.GetRRFPhase();
fout << endl;
} else if (sstr.BeginsWith("rrf_packing", TString::kIgnoreCase) && (fGlobal.GetFitType() == MSR_FITTYPE_SINGLE_HISTO_RRF)) {
fout.width(16);
fout << "rrf_packing ";
fout.width(8);
fout << left << fGlobal.GetRRFPacking();
fout << endl;
} else if (sstr.BeginsWith("data")) {
fout.width(16);
fout << left << "data";
@ -749,9 +781,15 @@ Int_t PMsrHandler::WriteMsrLogFile(const Bool_t messages)
case MSR_FITTYPE_SINGLE_HISTO:
fout << left << "fittype" << MSR_FITTYPE_SINGLE_HISTO << " (single histogram fit)" << endl;
break;
case MSR_FITTYPE_SINGLE_HISTO_RRF:
fout << left << "fittype" << MSR_FITTYPE_SINGLE_HISTO_RRF << " (single histogram RRF fit)" << endl;
break;
case MSR_FITTYPE_ASYM:
fout << left << "fittype" << MSR_FITTYPE_ASYM << " (asymmetry fit)" << endl ;
break;
case MSR_FITTYPE_ASYM_RRF:
fout << left << "fittype" << MSR_FITTYPE_ASYM_RRF << " (asymmetry RRF fit)" << endl ;
break;
case MSR_FITTYPE_MU_MINUS:
fout << left << "fittype" << MSR_FITTYPE_MU_MINUS << " (mu minus fit)" << endl ;
break;
@ -1107,9 +1145,15 @@ Int_t PMsrHandler::WriteMsrLogFile(const Bool_t messages)
case MSR_PLOT_SINGLE_HISTO:
fout << "PLOT " << fPlots[plotNo].fPlotType << " (single histo plot)" << endl;
break;
case MSR_PLOT_SINGLE_HISTO_RRF:
fout << "PLOT " << fPlots[plotNo].fPlotType << " (single histo RRF plot)" << endl;
break;
case MSR_PLOT_ASYM:
fout << "PLOT " << fPlots[plotNo].fPlotType << " (asymmetry plot)" << endl;
break;
case MSR_PLOT_ASYM_RRF:
fout << "PLOT " << fPlots[plotNo].fPlotType << " (asymmetry RRF plot)" << endl;
break;
case MSR_PLOT_MU_MINUS:
fout << "PLOT " << fPlots[plotNo].fPlotType << " (mu minus plot)" << endl;
break;
@ -1604,9 +1648,15 @@ Int_t PMsrHandler::WriteMsrFile(const Char_t *filename, map<UInt_t, TString> *co
case MSR_FITTYPE_SINGLE_HISTO:
fout << left << "fittype" << MSR_FITTYPE_SINGLE_HISTO << " (single histogram fit)" << endl;
break;
case MSR_FITTYPE_SINGLE_HISTO_RRF:
fout << left << "fittype" << MSR_FITTYPE_SINGLE_HISTO_RRF << " (single histogram RRF fit)" << endl;
break;
case MSR_FITTYPE_ASYM:
fout << left << "fittype" << MSR_FITTYPE_ASYM << " (asymmetry fit)" << endl ;
break;
case MSR_FITTYPE_ASYM_RRF:
fout << left << "fittype" << MSR_FITTYPE_ASYM_RRF << " (asymmetry RRF fit)" << endl ;
break;
case MSR_FITTYPE_MU_MINUS:
fout << left << "fittype" << MSR_FITTYPE_MU_MINUS << " (mu minus fit)" << endl ;
break;
@ -1618,6 +1668,30 @@ Int_t PMsrHandler::WriteMsrFile(const Char_t *filename, map<UInt_t, TString> *co
}
}
// RRF related stuff
if ((fGlobal.GetRRFFreq(fGlobal.GetRRFUnit().Data()) > 0.0) && (fGlobal.GetFitType() == MSR_FITTYPE_SINGLE_HISTO_RRF)) {
fout.width(16);
fout << left << "rrf_freq ";
fout.width(8);
fout << left << fGlobal.GetRRFFreq(fGlobal.GetRRFUnit().Data());
fout << " " << fGlobal.GetRRFUnit();
fout << endl;
}
if ((fGlobal.GetRRFPhase() != 0.0) && (fGlobal.GetFitType() == MSR_FITTYPE_SINGLE_HISTO_RRF)) {
fout.width(16);
fout << "rrf_phase ";
fout.width(8);
fout << left << fGlobal.GetRRFPhase();
fout << endl;
}
if ((fGlobal.GetRRFPacking() != -1) && (fGlobal.GetFitType() == MSR_FITTYPE_SINGLE_HISTO_RRF)) {
fout.width(16);
fout << "rrf_packing ";
fout.width(8);
fout << left << fGlobal.GetRRFPacking();
fout << endl;
}
// data range
if ((fGlobal.GetDataRange(0) != -1) || (fGlobal.GetDataRange(1) != -1) || (fGlobal.GetDataRange(2) != -1) || (fGlobal.GetDataRange(3) != -1)) {
fout.width(16);
@ -1765,9 +1839,15 @@ Int_t PMsrHandler::WriteMsrFile(const Char_t *filename, map<UInt_t, TString> *co
case MSR_FITTYPE_SINGLE_HISTO:
fout << left << "fittype" << MSR_FITTYPE_SINGLE_HISTO << " (single histogram fit)" << endl;
break;
case MSR_FITTYPE_SINGLE_HISTO_RRF:
fout << left << "fittype" << MSR_FITTYPE_SINGLE_HISTO_RRF << " (single histogram RRF fit)" << endl;
break;
case MSR_FITTYPE_ASYM:
fout << left << "fittype" << MSR_FITTYPE_ASYM << " (asymmetry fit)" << endl ;
break;
case MSR_FITTYPE_ASYM_RRF:
fout << left << "fittype" << MSR_FITTYPE_ASYM_RRF << " (asymmetry RRF fit)" << endl ;
break;
case MSR_FITTYPE_MU_MINUS:
fout << left << "fittype" << MSR_FITTYPE_MU_MINUS << " (mu minus fit)" << endl ;
break;
@ -2108,9 +2188,15 @@ Int_t PMsrHandler::WriteMsrFile(const Char_t *filename, map<UInt_t, TString> *co
case MSR_PLOT_SINGLE_HISTO:
fout << "PLOT " << fPlots[i].fPlotType << " (single histo plot)" << endl;
break;
case MSR_PLOT_SINGLE_HISTO_RRF:
fout << "PLOT " << fPlots[i].fPlotType << " (single histo RRF plot)" << endl;
break;
case MSR_PLOT_ASYM:
fout << "PLOT " << fPlots[i].fPlotType << " (asymmetry plot)" << endl;
break;
case MSR_PLOT_ASYM_RRF:
fout << "PLOT " << fPlots[i].fPlotType << " (asymmetry RRF plot)" << endl;
break;
case MSR_PLOT_MU_MINUS:
fout << "PLOT " << fPlots[i].fPlotType << " (mu minus plot)" << endl;
break;
@ -2796,8 +2882,8 @@ Bool_t PMsrHandler::HandleGlobalEntry(PMsrLines &lines)
TString str;
TObjArray *tokens = 0;
TObjString *ostr = 0;
Int_t ival;
Double_t dval;
Int_t ival = 0;
Double_t dval = 0.0;
UInt_t addT0Counter = 0;
// since this routine is called, a GLOBAL block is present
@ -2828,7 +2914,9 @@ Bool_t PMsrHandler::HandleGlobalEntry(PMsrLines &lines)
if (str.IsDigit()) {
Int_t fittype = str.Atoi();
if ((fittype == MSR_FITTYPE_SINGLE_HISTO) ||
(fittype == MSR_FITTYPE_SINGLE_HISTO_RRF) ||
(fittype == MSR_FITTYPE_ASYM) ||
(fittype == MSR_FITTYPE_ASYM_RRF) ||
(fittype == MSR_FITTYPE_MU_MINUS) ||
(fittype == MSR_FITTYPE_NON_MUSR)) {
global.SetFitType(fittype);
@ -2839,6 +2927,55 @@ Bool_t PMsrHandler::HandleGlobalEntry(PMsrLines &lines)
error = true;
}
}
} else if (iter->fLine.BeginsWith("rrf_freq", TString::kIgnoreCase)) {
if (tokens->GetEntries() < 3) {
error = true;
} else {
ostr = dynamic_cast<TObjString*>(tokens->At(1));
str = ostr->GetString();
if (str.IsFloat()) {
dval = str.Atof();
if (dval <= 0.0)
error = true;
}
if (!error) {
ostr = dynamic_cast<TObjString*>(tokens->At(2));
str = ostr->GetString();
global.SetRRFFreq(dval, str.Data());
if (global.GetRRFFreq(str.Data()) == RRF_FREQ_UNDEF)
error = true;
}
}
} else if (iter->fLine.BeginsWith("rrf_packing", TString::kIgnoreCase)) {
if (tokens->GetEntries() < 2) {
error = true;
} else {
ostr = dynamic_cast<TObjString*>(tokens->At(1));
str = ostr->GetString();
if (str.IsDigit()) {
ival = str.Atoi();
if (ival > 0) {
global.SetRRFPacking(ival);
} else {
error = true;
}
} else {
error = true;
}
}
} else if (iter->fLine.BeginsWith("rrf_phase", TString::kIgnoreCase)) {
if (tokens->GetEntries() < 2) {
error = true;
} else {
ostr = dynamic_cast<TObjString*>(tokens->At(1));
str = ostr->GetString();
if (str.IsFloat()) {
dval = str.Atof();
global.SetRRFPhase(dval);
} else {
error = true;
}
}
} else if (iter->fLine.BeginsWith("data", TString::kIgnoreCase)) { // data
if (tokens->GetEntries() < 3) {
error = true;
@ -3113,7 +3250,9 @@ Bool_t PMsrHandler::HandleRunEntry(PMsrLines &lines)
if (str.IsDigit()) {
Int_t fittype = str.Atoi();
if ((fittype == MSR_FITTYPE_SINGLE_HISTO) ||
(fittype == MSR_FITTYPE_SINGLE_HISTO_RRF) ||
(fittype == MSR_FITTYPE_ASYM) ||
(fittype == MSR_FITTYPE_ASYM_RRF) ||
(fittype == MSR_FITTYPE_MU_MINUS) ||
(fittype == MSR_FITTYPE_NON_MUSR)) {
param.SetFitType(fittype);
@ -4027,7 +4166,9 @@ Bool_t PMsrHandler::HandlePlotEntry(PMsrLines &lines)
error = true;
break;
case MSR_PLOT_SINGLE_HISTO: // like: runs 1 5 13
case MSR_PLOT_SINGLE_HISTO_RRF:
case MSR_PLOT_ASYM:
case MSR_PLOT_ASYM_RRF:
case MSR_PLOT_NON_MUSR:
case MSR_PLOT_MU_MINUS:
rl = new PStringNumberList(iter1->fLine.Data());
@ -4439,8 +4580,10 @@ Bool_t PMsrHandler::HandlePlotEntry(PMsrLines &lines)
cerr << endl << ">> [use_fit_ranges [ymin ymax]]";
cerr << endl << ">> [view_packing n]";
cerr << endl;
cerr << endl << ">> where <plot_type> is: 0=single histo asym,";
cerr << endl << ">> where <plot_type> is: 0=single histo,";
cerr << endl << ">> 1=RRF single histo,";
cerr << endl << ">> 2=forward-backward asym,";
cerr << endl << ">> 3=forward-backward RRF asym,";
cerr << endl << ">> 4=mu minus single histo,";
cerr << endl << ">> 8=non muSR.";
cerr << endl << ">> <run_list> is the list of runs, e.g. runs 1 3";
@ -5193,6 +5336,53 @@ Bool_t PMsrHandler::CheckRunBlockIntegrity()
fRuns[i].SetPacking(1);
}
break;
case PRUN_SINGLE_HISTO_RRF:
// check that there is a forward parameter number
if (fRuns[i].GetForwardHistoNo() == -1) {
cerr << endl << ">> PMsrHandler::CheckRunBlockIntegrity(): **ERROR** in RUN block number " << i+1;
cerr << endl << ">> forward parameter number not defined. Necessary for single histogram RRF fits." << endl;
return false;
}
if ((fRuns[i].GetNormParamNo() > static_cast<Int_t>(fParam.size())) && !fFourierOnly) {
// check if forward histogram number is a function
if (fRuns[i].GetNormParamNo() - MSR_PARAM_FUN_OFFSET > static_cast<Int_t>(fParam.size())) {
cerr << endl << ">> PMsrHandler::CheckRunBlockIntegrity(): **ERROR** in RUN block number " << i+1;
cerr << endl << ">> forward histogram number " << fRuns[i].GetNormParamNo() << " is larger than the number of fit parameters (" << fParam.size() << ").";
cerr << endl << ">> Consider to check the manual ;-)" << endl;
return false;
}
}
// check fit range
if (!fRuns[i].IsFitRangeInBin() && !fFourierOnly) { // fit range given as times in usec (RUN block)
if ((fRuns[i].GetFitRange(0) == PMUSR_UNDEFINED) || (fRuns[i].GetFitRange(1) == PMUSR_UNDEFINED)) { // check fit range in RUN block
if (!fGlobal.IsFitRangeInBin()) { // fit range given as times in usec (GLOBAL block)
if ((fGlobal.GetFitRange(0) == PMUSR_UNDEFINED) || (fGlobal.GetFitRange(1) == PMUSR_UNDEFINED)) { // check fit range in GLOBAL block
cerr << endl << "PMsrHandler::CheckRunBlockIntegrity(): **ERROR** in RUN block number " << i+1;
cerr << endl << " Fit range is not defined. Necessary for single histogram fits." << endl;
return false;
}
}
}
}
// check number of T0's provided
if ((fRuns[i].GetT0BinSize() > fRuns[i].GetForwardHistoNoSize()) &&
(fGlobal.GetT0BinSize() > fRuns[i].GetForwardHistoNoSize())) {
cerr << endl << ">> PMsrHandler::CheckRunBlockIntegrity(): **ERROR** in RUN block number " << i+1;
cerr << endl << ">> Found " << fRuns[i].GetT0BinSize() << " T0 entries. Expecting only " << fRuns[i].GetForwardHistoNoSize() << ". Needs to be fixed." << endl;
cerr << endl << ">> In GLOBAL block: " << fGlobal.GetT0BinSize() << " T0 entries. Expecting only " << fRuns[i].GetForwardHistoNoSize() << ". Needs to be fixed." << endl;
return false;
}
// check that RRF frequency is given
if (fGlobal.GetRRFUnitTag() == RRF_UNIT_UNDEF) {
cerr << endl << ">> PMsrHandler::CheckRunBlockIntegrity(): **ERROR** no RRF frequency found in the GLOBAL block." << endl;
return false;
}
// check that RRF packing is given
if (fGlobal.GetRRFPacking() == -1) {
cerr << endl << ">> PMsrHandler::CheckRunBlockIntegrity(): **ERROR** no RRF packing found in the GLOBAL block." << endl;
return false;
}
break;
case PRUN_ASYMMETRY:
// check alpha
if ((fRuns[i].GetAlphaParamNo() == -1) && !fFourierOnly) {
@ -5245,6 +5435,62 @@ Bool_t PMsrHandler::CheckRunBlockIntegrity()
fRuns[i].SetPacking(1);
}
break;
case PRUN_ASYMMETRY_RRF:
// check alpha
if ((fRuns[i].GetAlphaParamNo() == -1) && !fFourierOnly) {
cerr << endl << ">> PMsrHandler::CheckRunBlockIntegrity(): **ERROR** in RUN block number " << i+1;
cerr << endl << ">> alpha parameter number missing which is needed for an asymmetry RRF fit.";
cerr << endl << ">> Consider to check the manual ;-)" << endl;
return false;
}
// check that there is a forward parameter number
if (fRuns[i].GetForwardHistoNo() == -1) {
cerr << endl << ">> PMsrHandler::CheckRunBlockIntegrity(): **ERROR** in RUN block number " << i+1;
cerr << endl << ">> forward histogram number not defined. Necessary for asymmetry RRF fits." << endl;
return false;
}
// check that there is a backward parameter number
if (fRuns[i].GetBackwardHistoNo() == -1) {
cerr << endl << ">> PMsrHandler::CheckRunBlockIntegrity(): **ERROR** in RUN block number " << i+1;
cerr << endl << ">> backward histogram number not defined. Necessary for asymmetry RRF fits." << endl;
return false;
}
// check fit range
if (!fRuns[i].IsFitRangeInBin()) { // fit range given as times in usec
if ((fRuns[i].GetFitRange(0) == PMUSR_UNDEFINED) || (fRuns[i].GetFitRange(1) == PMUSR_UNDEFINED)) {
if ((fGlobal.GetFitRange(0) == PMUSR_UNDEFINED) || (fGlobal.GetFitRange(1) == PMUSR_UNDEFINED)) {
cerr << endl << ">> PMsrHandler::CheckRunBlockIntegrity(): **ERROR** in RUN block number " << i+1;
cerr << endl << ">> Fit range is not defined, also NOT present in the GLOBAL block. Necessary for asymmetry RRF fits." << endl;
return false;
}
}
}
// check number of T0's provided
if ((fRuns[i].GetT0BinSize() > 2*fRuns[i].GetForwardHistoNoSize()) &&
(fGlobal.GetT0BinSize() > 2*fRuns[i].GetForwardHistoNoSize())) {
cerr << endl << ">> PMsrHandler::CheckRunBlockIntegrity(): **ERROR** in RUN block number " << i+1;
cerr << endl << ">> Found " << fRuns[i].GetT0BinSize() << " T0 entries. Expecting only " << 2*fRuns[i].GetForwardHistoNoSize() << " in forward. Needs to be fixed." << endl;
cerr << endl << ">> In GLOBAL block: " << fGlobal.GetT0BinSize() << " T0 entries. Expecting only " << 2*fRuns[i].GetForwardHistoNoSize() << ". Needs to be fixed." << endl;
return false;
}
if ((fRuns[i].GetT0BinSize() > 2*fRuns[i].GetBackwardHistoNoSize()) &&
(fGlobal.GetT0BinSize() > 2*fRuns[i].GetBackwardHistoNoSize())) {
cerr << endl << ">> PMsrHandler::CheckRunBlockIntegrity(): **ERROR** in RUN block number " << i+1;
cerr << endl << ">> Found " << fRuns[i].GetT0BinSize() << " T0 entries. Expecting only " << 2*fRuns[i].GetBackwardHistoNoSize() << " in backward. Needs to be fixed." << endl;
cerr << endl << ">> In GLOBAL block: " << fGlobal.GetT0BinSize() << " T0 entries. Expecting only " << 2*fRuns[i].GetBackwardHistoNoSize() << ". Needs to be fixed." << endl;
return false;
}
// check that RRF frequency is given
if (fGlobal.GetRRFUnitTag() == RRF_UNIT_UNDEF) {
cerr << endl << ">> PMsrHandler::CheckRunBlockIntegrity(): **ERROR** no RRF frequency found in the GLOBAL block." << endl;
return false;
}
// check that RRF packing is given
if (fGlobal.GetRRFPacking() == -1) {
cerr << endl << ">> PMsrHandler::CheckRunBlockIntegrity(): **ERROR** no RRF packing found in the GLOBAL block." << endl;
return false;
}
break;
case PRUN_MU_MINUS:
// needs eventually to be implemented
break;
@ -5631,6 +5877,103 @@ void PMsrHandler::CheckMaxLikelihood()
}
}
//--------------------------------------------------------------------------
// CheckRRFSettings (public)
//--------------------------------------------------------------------------
/**
* <p>Make sure that if RRF settings are found in the GLOBAL section, the fit types
* in the RUN blocks correspond.
*/
Bool_t PMsrHandler::CheckRRFSettings()
{
Bool_t result = true;
Int_t fittype = fGlobal.GetFitType();
// first set of tests: if RRF parameters are set, check if RRF fit is chosen.
if (fGlobal.GetRRFFreq(fGlobal.GetRRFUnit().Data()) != RRF_FREQ_UNDEF) {
if (fittype != -1) { // check if GLOBAL fittype is set
if ((fittype != MSR_FITTYPE_SINGLE_HISTO_RRF) &&
(fittype != MSR_FITTYPE_ASYM_RRF)) {
cerr << endl << ">> PMsrHandler::CheckRRFSettings: **ERROR** found GLOBAL fittype " << fittype << " and";
cerr << endl << ">> RRF settings in the GLOBAL section. This is NOT compatible. Fix it first.";
result = false;
}
} else { // GLOBAL fittype is NOT set
for (UInt_t i=0; i<fRuns.size(); i++) {
fittype = fRuns[i].GetFitType();
if ((fittype != MSR_FITTYPE_SINGLE_HISTO_RRF) &&
(fittype != MSR_FITTYPE_ASYM_RRF)) {
cerr << endl << ">> PMsrHandler::CheckRRFSettings: **ERROR** found RUN with fittype " << fittype << " and";
cerr << endl << ">> RRF settings in the GLOBAL section. This is NOT compatible. Fix it first.";
result = false;
break;
}
}
}
} else {
if (fGlobal.GetRRFPacking() != -1) {
cerr << endl << ">> PMsrHandler::CheckRRFSettings: **WARNING** found in the GLOBAL section rrf_packing, without";
cerr << endl << ">> rrf_freq. Doesn't make any sense. Will drop rrf_packing";
cerr << endl << endl;
fGlobal.SetRRFPacking(-1);
}
if (fGlobal.GetRRFPhase() != 0.0) {
cerr << endl << ">> PMsrHandler::CheckRRFSettings: **WARNING** found in the GLOBAL section rrf_phase, without";
cerr << endl << ">> rrf_freq. Doesn't make any sense. Will drop rrf_phase";
cerr << endl << endl;
fGlobal.SetRRFPhase(0.0);
}
}
// if not a RRF fit, done at this point
if ((fittype != MSR_FITTYPE_SINGLE_HISTO_RRF) &&
(fittype != MSR_FITTYPE_ASYM_RRF)) {
return true;
}
// second set of tests: if RRF fit is chosen, do I find the necessary RRF parameters?
fittype = fGlobal.GetFitType();
if ((fittype == MSR_FITTYPE_SINGLE_HISTO_RRF) ||
(fittype == MSR_FITTYPE_ASYM_RRF)) { // make sure RRF freq and RRF packing are set
if (fGlobal.GetRRFFreq(fGlobal.GetRRFUnit().Data()) == RRF_FREQ_UNDEF) {
cerr << endl << ">> PMsrHandler::CheckRRFSettings: **ERROR** RRF fit chosen, but";
cerr << endl << ">> no RRF frequency found in the GLOBAL section! Fix it.";
return false;
}
if (fGlobal.GetRRFPacking() == -1) {
cerr << endl << ">> PMsrHandler::CheckRRFSettings: **ERROR** RRF fit chosen, but";
cerr << endl << ">> no RRF packing found in the GLOBAL section! Fix it.";
return false;
}
} else { // check single runs for RRF
UInt_t rrfFitCounter = 0;
for (UInt_t i=0; i<fRuns.size(); i++) {
fittype = fRuns[i].GetFitType();
if ((fittype == MSR_FITTYPE_SINGLE_HISTO_RRF) ||
(fittype == MSR_FITTYPE_ASYM_RRF)) { // make sure RRF freq and RRF packing are set
rrfFitCounter++;
}
}
if (rrfFitCounter != fRuns.size()) {
cerr << endl << ">> PMsrHandler::CheckRRFSettings: **ERROR** #Runs (" << fRuns.size() << ") != # RRF fits found (" << rrfFitCounter << ")";
cerr << endl << ">> This is currently not supported.";
return false;
}
if (fGlobal.GetRRFFreq(fGlobal.GetRRFUnit().Data()) == RRF_FREQ_UNDEF) {
cerr << endl << ">> PMsrHandler::CheckRRFSettings: **ERROR** RRF fit chosen, but";
cerr << endl << ">> no RRF frequency found in the GLOBAL section! Fix it.";
return false;
}
if (fGlobal.GetRRFPacking() == -1) {
cerr << endl << ">> PMsrHandler::CheckRRFSettings: **ERROR** RRF fit chosen, but";
cerr << endl << ">> no RRF packing found in the GLOBAL section! Fix it.";
return false;
}
}
return result;
}
//--------------------------------------------------------------------------
// GetGroupingString (public)
//--------------------------------------------------------------------------

View File

@ -8,7 +8,7 @@
***************************************************************************/
/***************************************************************************
* Copyright (C) 2007-2014 by Andreas Suter *
* Copyright (C) 2007-2016 by Andreas Suter *
* andreas.suter@psi.ch *
* *
* This program is free software; you can redistribute it and/or modify *
@ -34,6 +34,8 @@ using namespace std;
#include <boost/algorithm/string.hpp>
using namespace boost;
#include "TMath.h"
#include "PMusr.h"
//+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
@ -708,6 +710,10 @@ void PRawRunData::SetTempError(const UInt_t idx, const Double_t errTemp)
PMsrGlobalBlock::PMsrGlobalBlock()
{
fGlobalPresent = false;
fRRFFreq = RRF_FREQ_UNDEF; // rotating reference frequency in units given by fRRFUnitTag. Only needed for fittype 1
fRRFUnitTag = RRF_UNIT_UNDEF; // RRF unit tag. Default: undefined
fRRFPhase = 0.0;
fRRFPacking = -1; // undefined RRF packing/rebinning
fFitType = -1; // undefined fit type
for (UInt_t i=0; i<4; i++) {
fDataRange[i] = -1; // undefined data bin range
@ -720,6 +726,135 @@ PMsrGlobalBlock::PMsrGlobalBlock()
fPacking = -1; // undefined packing/rebinning
}
//--------------------------------------------------------------------------
// GetRRFFreq (public)
//--------------------------------------------------------------------------
/**
* <p> get RRF frequency value in specific units. If units is unknown, RRF_UNDEF_FREQ will be returned.
*
* \param unit unit string in which the units shall be given
*/
Double_t PMsrGlobalBlock::GetRRFFreq(const char *unit)
{
Double_t freq = 0.0;
// check that the units given make sense
TString unitStr = unit;
Int_t unitTag = RRF_UNIT_UNDEF;
if (!unitStr.CompareTo("MHz", TString::kIgnoreCase))
unitTag = RRF_UNIT_MHz;
else if (!unitStr.CompareTo("Mc", TString::kIgnoreCase))
unitTag = RRF_UNIT_Mcs;
else if (!unitStr.CompareTo("T", TString::kIgnoreCase))
unitTag = RRF_UNIT_T;
else
return RRF_FREQ_UNDEF;
// calc the conversion factor
if (unitTag == fRRFUnitTag)
freq = fRRFFreq;
else if ((unitTag == RRF_UNIT_MHz) && (fRRFUnitTag == RRF_UNIT_Mcs))
freq = fRRFFreq/TMath::TwoPi();
else if ((unitTag == RRF_UNIT_MHz) && (fRRFUnitTag == RRF_UNIT_T))
freq = fRRFFreq*1e4*GAMMA_BAR_MUON; // 1e4 need for T -> G since GAMMA_BAR_MUON is given in MHz/G
else if ((unitTag == RRF_UNIT_Mcs) && (fRRFUnitTag == RRF_UNIT_MHz))
freq = fRRFFreq*TMath::TwoPi();
else if ((unitTag == RRF_UNIT_Mcs) && (fRRFUnitTag == RRF_UNIT_T))
freq = fRRFFreq*1e4*TMath::TwoPi()*GAMMA_BAR_MUON; // 1e4 need for T -> G since GAMMA_BAR_MUON is given in MHz/G
else if ((unitTag == RRF_UNIT_T) && (fRRFUnitTag == RRF_UNIT_MHz))
freq = fRRFFreq/GAMMA_BAR_MUON*1e-4; // 1e-4 need for G -> T since GAMMA_BAR_MUON is given in MHz/G
else if ((unitTag == RRF_UNIT_T) && (fRRFUnitTag == RRF_UNIT_Mcs))
freq = fRRFFreq/(TMath::TwoPi()*GAMMA_BAR_MUON)*1e-4; // 1e-4 need for G -> T since GAMMA_BAR_MUON is given in MHz/G
return freq;
}
//--------------------------------------------------------------------------
// SetRRFFreq (public)
//--------------------------------------------------------------------------
/**
* <p> set RRF frequency value in specific units. If units is unknown, 0.0 will be set.
*
* \param RRF frequency value
* \param unit unit string in which the units shall be given
*/
void PMsrGlobalBlock::SetRRFFreq(Double_t freq, const char *unit)
{
// check that the units given make sense
TString unitStr = unit;
Int_t unitTag = RRF_UNIT_UNDEF;
if (!unitStr.CompareTo("MHz", TString::kIgnoreCase))
unitTag = RRF_UNIT_MHz;
else if (!unitStr.CompareTo("Mc", TString::kIgnoreCase))
unitTag = RRF_UNIT_Mcs;
else if (!unitStr.CompareTo("T", TString::kIgnoreCase))
unitTag = RRF_UNIT_T;
else {
cerr << endl << ">> PMsrGlobalBlock::SetRRFFreq: **ERROR** found undefined RRF unit '" << unit << "'!";
cerr << endl << ">> Will set RRF frequency to 0.0." << endl;
fRRFFreq = 0.0;
fRRFUnitTag = RRF_UNIT_UNDEF;
}
fRRFFreq = freq;
fRRFUnitTag = unitTag;
}
//--------------------------------------------------------------------------
// GetRRFUnit (public)
//--------------------------------------------------------------------------
/**
* <p> returns RRF frequency unit.
*/
TString PMsrGlobalBlock::GetRRFUnit()
{
TString unit;
switch (fRRFUnitTag) {
case RRF_UNIT_UNDEF:
unit = TString("??");
break;
case RRF_UNIT_kHz:
unit = TString("kHz");
break;
case RRF_UNIT_MHz:
unit = TString("MHz");
break;
case RRF_UNIT_Mcs:
unit = TString("Mc");
break;
case RRF_UNIT_G:
unit = TString("G");
break;
case RRF_UNIT_T:
unit = TString("T");
break;
default:
unit = TString("??");
break;
}
return unit;
}
//--------------------------------------------------------------------------
// SetRRFPacking (public)
//--------------------------------------------------------------------------
/**
* <p> set RRF packing.
*
* \param RRF packing
*/
void PMsrGlobalBlock::SetRRFPacking(Int_t pack)
{
if (pack <= 0) {
cerr << endl << "PMsrGlobalBlock::SetRRFPacking: **WARNING** found RRF packing <= 0. Likely doesn't make any sense." << endl;
fRRFPacking = -1; // set to undefined
}
fRRFPacking = pack;
}
//--------------------------------------------------------------------------
// GetDataRange (public)
//--------------------------------------------------------------------------

View File

@ -8,7 +8,7 @@
***************************************************************************/
/***************************************************************************
* Copyright (C) 2007-2014 by Andreas Suter *
* Copyright (C) 2007-2016 by Andreas Suter *
* andreas.suter@psi.ch *
* *
* This program is free software; you can redistribute it and/or modify *
@ -75,7 +75,7 @@ PMusrCanvasPlotRange::PMusrCanvasPlotRange()
void PMusrCanvasPlotRange::SetXRange(Double_t xmin, Double_t xmax)
{
if (xmin > xmax) {
cerr << endl << "PMusrCanvasPlotRange::SetXRange: **WARNING** xmin > xmax, will swap them." << endl;
cerr << endl << ">> PMusrCanvasPlotRange::SetXRange(): **WARNING** xmin > xmax, will swap them." << endl;
fXmin = xmax;
fXmax = xmin;
} else {
@ -97,7 +97,7 @@ void PMusrCanvasPlotRange::SetXRange(Double_t xmin, Double_t xmax)
void PMusrCanvasPlotRange::SetYRange(Double_t ymin, Double_t ymax)
{
if (ymin > ymax) {
cerr << endl << "PMusrCanvasPlotRange::SetYRange: **WARNING** ymin > ymax, will swap them." << endl;
cerr << endl << ">> PMusrCanvasPlotRange::SetYRange(): **WARNING** ymin > ymax, will swap them." << endl;
fYmin = ymax;
fYmax = ymin;
} else {
@ -405,37 +405,76 @@ void PMusrCanvas::SetMsrHandler(PMsrHandler *msrHandler)
}
// check if RRF data are present
if ((fMsrHandler->GetMsrPlotList()->at(0).fRRFPacking > 0) &&
(fMsrHandler->GetMsrPlotList()->at(0).fRRFFreq != 0.0)) {
if (((fMsrHandler->GetMsrPlotList()->at(0).fRRFPacking > 0) &&
(fMsrHandler->GetMsrPlotList()->at(0).fRRFFreq != 0.0)) ||
(fMsrHandler->GetMsrGlobal()->GetRRFPacking() > 0 &&
fMsrHandler->GetMsrGlobal()->GetRRFUnit().CompareTo("??"))) {
fRRFLatexText = new TLatex();
fRRFLatexText->SetNDC(kTRUE);
fRRFLatexText->SetTextFont(62);
fRRFLatexText->SetTextSize(0.03);
fRRFText = new TString("RRF: ");
if (fMsrHandler->GetMsrPlotList()->at(0).fRRFUnit == RRF_UNIT_kHz) {
*fRRFText += TString("#nu_{RRF} = ");
*fRRFText += fMsrHandler->GetMsrPlotList()->at(0).fRRFFreq;
*fRRFText += TString(" (kHz)");
} else if (fMsrHandler->GetMsrPlotList()->at(0).fRRFUnit == RRF_UNIT_MHz) {
*fRRFText += TString("#nu_{RRF} = ");
*fRRFText += fMsrHandler->GetMsrPlotList()->at(0).fRRFFreq;
*fRRFText += TString(" (MHz)");
} else if (fMsrHandler->GetMsrPlotList()->at(0).fRRFUnit == RRF_UNIT_Mcs) {
*fRRFText += TString("#omega_{RRF} = ");
*fRRFText += fMsrHandler->GetMsrPlotList()->at(0).fRRFFreq;
*fRRFText += TString(" (Mc/s)");
} else if (fMsrHandler->GetMsrPlotList()->at(0).fRRFUnit == RRF_UNIT_G) {
*fRRFText += TString("B_{RRF} = ");
*fRRFText += fMsrHandler->GetMsrPlotList()->at(0).fRRFFreq;
*fRRFText += TString(" (G)");
} else if (fMsrHandler->GetMsrPlotList()->at(0).fRRFUnit == RRF_UNIT_T) {
*fRRFText += TString("B_{RRF} = ");
*fRRFText += fMsrHandler->GetMsrPlotList()->at(0).fRRFFreq;
*fRRFText += TString(" (T)");
Int_t rrfUnitTag = -1;
Double_t rrfFreq = 0.0;
if (fMsrHandler->GetMsrPlotList()->at(0).fRRFPacking > 0) { // RRF single histo PLOT
fRRFText = new TString("RRF: ");
rrfUnitTag = fMsrHandler->GetMsrPlotList()->at(0).fRRFUnit;
rrfFreq = fMsrHandler->GetMsrPlotList()->at(0).fRRFFreq;
TString rrfFreqStr("");
rrfFreqStr.Form("%.5g", rrfFreq);
if (rrfUnitTag == RRF_UNIT_kHz) {
*fRRFText += TString("#nu_{RRF} = ");
*fRRFText += rrfFreq;
*fRRFText += TString(" (kHz)");
} else if (rrfUnitTag == RRF_UNIT_MHz) {
*fRRFText += TString("#nu_{RRF} = ");
*fRRFText += rrfFreqStr;
*fRRFText += TString(" (MHz)");
} else if (rrfUnitTag == RRF_UNIT_Mcs) {
*fRRFText += TString("#omega_{RRF} = ");
*fRRFText += rrfFreqStr;
*fRRFText += TString(" (Mc/s)");
} else if (rrfUnitTag == RRF_UNIT_G) {
*fRRFText += TString("B_{RRF} = ");
*fRRFText += rrfFreqStr;
*fRRFText += TString(" (G)");
} else if (rrfUnitTag == RRF_UNIT_T) {
*fRRFText += TString("B_{RRF} = ");
*fRRFText += rrfFreqStr;
*fRRFText += TString(" (T)");
}
*fRRFText += TString(", RRF packing = ");
*fRRFText += fMsrHandler->GetMsrPlotList()->at(0).fRRFPacking;
} else { // RRF single histo FIT
fRRFText = new TString("RRF: ");
rrfUnitTag = fMsrHandler->GetMsrGlobal()->GetRRFUnitTag();
rrfFreq = fMsrHandler->GetMsrGlobal()->GetRRFFreq(fMsrHandler->GetMsrGlobal()->GetRRFUnit().Data());
TString rrfFreqStr("");
rrfFreqStr.Form("%.5g", rrfFreq);
if (rrfUnitTag == RRF_UNIT_kHz) {
*fRRFText += TString("#nu_{RRF} = ");
*fRRFText += rrfFreqStr;
*fRRFText += TString(" (kHz)");
} else if (rrfUnitTag == RRF_UNIT_MHz) {
*fRRFText += TString("#nu_{RRF} = ");
*fRRFText += rrfFreqStr;
*fRRFText += TString(" (MHz)");
} else if (rrfUnitTag == RRF_UNIT_Mcs) {
*fRRFText += TString("#omega_{RRF} = ");
*fRRFText += rrfFreqStr;
*fRRFText += TString(" (Mc/s)");
} else if (rrfUnitTag == RRF_UNIT_G) {
*fRRFText += TString("B_{RRF} = ");
*fRRFText += rrfFreqStr;
*fRRFText += TString(" (G)");
} else if (rrfUnitTag == RRF_UNIT_T) {
*fRRFText += TString("B_{RRF} = ");
*fRRFText += rrfFreqStr;
*fRRFText += TString(" (T)");
}
*fRRFText += TString(", RRF packing = ");
*fRRFText += fMsrHandler->GetMsrGlobal()->GetRRFPacking();
}
*fRRFText += TString(", RRF packing = ");
*fRRFText += fMsrHandler->GetMsrPlotList()->at(0).fRRFPacking;
}
}
@ -609,7 +648,7 @@ void PMusrCanvas::UpdateDataTheoryPad()
// first check that plot number is smaller than the maximal number of runs
if ((Int_t)plotInfo.fRuns[i] > (Int_t)runs.size()) {
fValid = false;
cerr << endl << "PMusrCanvas::UpdateDataTheoryPad: **ERROR** run plot number " << (Int_t)plotInfo.fRuns[i] << " is larger than the number of runs " << runs.size();
cerr << endl << ">> PMusrCanvas::UpdateDataTheoryPad(): **ERROR** run plot number " << (Int_t)plotInfo.fRuns[i] << " is larger than the number of runs " << runs.size();
cerr << endl;
return;
}
@ -620,7 +659,7 @@ void PMusrCanvas::UpdateDataTheoryPad()
}
if (fitType == -1) {
fValid = false;
cerr << endl << "PMusrCanvas::UpdateDataTheoryPad: **ERROR** plottype = " << fPlotType;
cerr << endl << ">> PMusrCanvas::UpdateDataTheoryPad(): **ERROR** plottype = " << fPlotType;
cerr << ", fittype = " << runs[runNo].GetFitType() << "(RUN block)/";
cerr << "fittype = " << globalBlock->GetFitType() << "(GLOBAL block). However, they have to correspond!";
cerr << endl;
@ -643,7 +682,19 @@ void PMusrCanvas::UpdateDataTheoryPad()
if (!data) { // something wrong
fValid = false;
// error message
cerr << endl << "PMusrCanvas::UpdateDataTheoryPad: **ERROR** couldn't obtain run no " << runNo << " for a single histogram plot";
cerr << endl << ">> PMusrCanvas::UpdateDataTheoryPad(): **ERROR** couldn't obtain run no " << runNo << " for a single histogram plot";
cerr << endl;
return;
}
// handle data
HandleDataSet(i, runNo, data);
break;
case MSR_FITTYPE_SINGLE_HISTO_RRF:
data = fRunList->GetSingleHistoRRF(runNo, PRunListCollection::kRunNo);
if (!data) { // something wrong
fValid = false;
// error message
cerr << endl << ">> PMusrCanvas::UpdateDataTheoryPad(): **ERROR** couldn't obtain run no " << runNo << " for a single histogram RRF plot";
cerr << endl;
return;
}
@ -655,7 +706,19 @@ void PMusrCanvas::UpdateDataTheoryPad()
if (!data) { // something wrong
fValid = false;
// error message
cerr << endl << "PMusrCanvas::UpdateDataTheoryPad: **ERROR** couldn't obtain run no " << runNo << " for a asymmetry plot";
cerr << endl << ">> PMusrCanvas::UpdateDataTheoryPad(): **ERROR** couldn't obtain run no " << runNo << " for a asymmetry plot";
cerr << endl;
return;
}
// handle data
HandleDataSet(i, runNo, data);
break;
case MSR_FITTYPE_ASYM_RRF:
data = fRunList->GetAsymmetryRRF(runNo, PRunListCollection::kRunNo);
if (!data) { // something wrong
fValid = false;
// error message
cerr << endl << ">> PMusrCanvas::UpdateDataTheoryPad(): **ERROR** couldn't obtain run no " << runNo << " for a asymmetry RRF plot";
cerr << endl;
return;
}
@ -667,7 +730,7 @@ void PMusrCanvas::UpdateDataTheoryPad()
if (!data) { // something wrong
fValid = false;
// error message
cerr << endl << "PMusrCanvas::UpdateDataTheoryPad: **ERROR** couldn't obtain run no " << runNo << " for a mu minus single histogram plot";
cerr << endl << ">> PMusrCanvas::UpdateDataTheoryPad(): **ERROR** couldn't obtain run no " << runNo << " for a mu minus single histogram plot";
cerr << endl;
return;
}
@ -679,7 +742,7 @@ void PMusrCanvas::UpdateDataTheoryPad()
if (!data) { // something wrong
fValid = false;
// error message
cerr << endl << "PMusrCanvas::UpdateDataTheoryPad: **ERROR** couldn't obtain run no " << runNo << " for a none musr data plot";
cerr << endl << ">> PMusrCanvas::UpdateDataTheoryPad(): **ERROR** couldn't obtain run no " << runNo << " for a none musr data plot";
cerr << endl;
return;
}
@ -699,7 +762,7 @@ void PMusrCanvas::UpdateDataTheoryPad()
default:
fValid = false;
// error message
cerr << endl << "PMusrCanvas::UpdateDataTheoryPad: **ERROR** wrong plottype tag?!";
cerr << endl << ">> PMusrCanvas::UpdateDataTheoryPad(): **ERROR** wrong plottype tag?!";
cerr << endl;
return;
break;
@ -822,13 +885,15 @@ void PMusrCanvas::UpdateInfoPad()
else
tstr = *runs[runNo].GetRunName() + TString(","); // run_name
// histo info (depending on the fittype
if (runs[runNo].GetFitType() == MSR_FITTYPE_SINGLE_HISTO) {
if ((runs[runNo].GetFitType() == MSR_FITTYPE_SINGLE_HISTO) ||
(runs[runNo].GetFitType() == MSR_FITTYPE_SINGLE_HISTO_RRF)) {
tstr += TString("h:");
TString grouping;
fMsrHandler->GetGroupingString(runNo, "forward", grouping);
tstr += grouping;
tstr += TString(",");
} else if (runs[runNo].GetFitType() == MSR_FITTYPE_ASYM) {
} else if ((runs[runNo].GetFitType() == MSR_FITTYPE_ASYM) ||
(runs[runNo].GetFitType() == MSR_FITTYPE_ASYM_RRF)) {
tstr += TString("h:");
TString grouping;
fMsrHandler->GetGroupingString(runNo, "forward", grouping);
@ -1409,7 +1474,7 @@ void PMusrCanvas::SaveGraphicsAndQuit(Char_t *fileName, Char_t *graphicsFormat)
}
if (idx == -1) {
cerr << endl << "PMusrCanvas::SaveGraphicsAndQuit **ERROR**: fileName (" << fileName << ") is invalid." << endl;
cerr << endl << ">> PMusrCanvas::SaveGraphicsAndQuit(): **ERROR** fileName (" << fileName << ") is invalid." << endl;
return;
}
@ -1439,7 +1504,7 @@ void PMusrCanvas::SaveGraphicsAndQuit(Char_t *fileName, Char_t *graphicsFormat)
void PMusrCanvas::ExportData(const Char_t *fileName)
{
if (fileName == 0) { // path file name NOT provided, generate a default path file name
cerr << endl << ">> PMusrCanvas::ExportData **ERROR** NO path file name provided. Will do nothing." << endl;
cerr << endl << ">> PMusrCanvas::ExportData(): **ERROR** NO path file name provided. Will do nothing." << endl;
return;
}
@ -2094,7 +2159,7 @@ void PMusrCanvas::ExportData(const Char_t *fileName)
// open output data-file
fout.open(fileName, iostream::out);
if (!fout.is_open()) {
cerr << endl << ">> PMusrCanvas::ExportData: **ERROR** couldn't open file " << fileName << " for writing." << endl;
cerr << endl << ">> PMusrCanvas::ExportData(): **ERROR** couldn't open file " << fileName << " for writing." << endl;
return;
}
@ -2433,7 +2498,7 @@ void PMusrCanvas::InitMusrCanvas(const Char_t* title, Int_t wtopx, Int_t wtopy,
canvasName += fPlotNumber;
fMainCanvas = new TCanvas(canvasName.Data(), title, wtopx, wtopy, ww, wh);
if (fMainCanvas == 0) {
cerr << endl << "PMusrCanvas::PMusrCanvas: **PANIC ERROR**: Couldn't invoke " << canvasName.Data();
cerr << endl << ">> PMusrCanvas::PMusrCanvas(): **PANIC ERROR** Couldn't invoke " << canvasName.Data();
cerr << endl;
return;
}
@ -2479,7 +2544,7 @@ void PMusrCanvas::InitMusrCanvas(const Char_t* title, Int_t wtopx, Int_t wtopy,
// title pad
fTitlePad = new TPaveText(0.0, YTITLE, 1.0, 1.0, "NDC");
if (fTitlePad == 0) {
cerr << endl << "PMusrCanvas::PMusrCanvas: **PANIC ERROR**: Couldn't invoke fTitlePad";
cerr << endl << ">> PMusrCanvas::PMusrCanvas(): **PANIC ERROR** Couldn't invoke fTitlePad";
cerr << endl;
return;
}
@ -2491,7 +2556,7 @@ void PMusrCanvas::InitMusrCanvas(const Char_t* title, Int_t wtopx, Int_t wtopy,
// data/theory pad
fDataTheoryPad = new TPad("dataTheoryPad", "dataTheoryPad", 0.0, YINFO, XTHEO, YTITLE);
if (fDataTheoryPad == 0) {
cerr << endl << "PMusrCanvas::PMusrCanvas: **PANIC ERROR**: Couldn't invoke fDataTheoryPad";
cerr << endl << ">> PMusrCanvas::PMusrCanvas(): **PANIC ERROR** Couldn't invoke fDataTheoryPad";
cerr << endl;
return;
}
@ -2501,7 +2566,7 @@ void PMusrCanvas::InitMusrCanvas(const Char_t* title, Int_t wtopx, Int_t wtopy,
// parameter pad
fParameterPad = new TPaveText(XTHEO, 0.5, 1.0, YTITLE, "NDC");
if (fParameterPad == 0) {
cerr << endl << "PMusrCanvas::PMusrCanvas: **PANIC ERROR**: Couldn't invoke fParameterPad";
cerr << endl << ">> PMusrCanvas::PMusrCanvas(): **PANIC ERROR** Couldn't invoke fParameterPad";
cerr << endl;
return;
}
@ -2512,7 +2577,7 @@ void PMusrCanvas::InitMusrCanvas(const Char_t* title, Int_t wtopx, Int_t wtopy,
// theory pad
fTheoryPad = new TPaveText(XTHEO, 0.1, 1.0, 0.5, "NDC");
if (fTheoryPad == 0) {
cerr << endl << "PMusrCanvas::PMusrCanvas: **PANIC ERROR**: Couldn't invoke fTheoryPad";
cerr << endl << ">> PMusrCanvas::PMusrCanvas(): **PANIC ERROR** Couldn't invoke fTheoryPad";
cerr << endl;
return;
}
@ -2524,7 +2589,7 @@ void PMusrCanvas::InitMusrCanvas(const Char_t* title, Int_t wtopx, Int_t wtopy,
// info pad
fInfoPad = new TLegend(0.0, 0.0, 1.0, YINFO, "NDC");
if (fInfoPad == 0) {
cerr << endl << "PMusrCanvas::PMusrCanvas: **PANIC ERROR**: Couldn't invoke fInfoPad";
cerr << endl << ">> PMusrCanvas::PMusrCanvas(): **PANIC ERROR** Couldn't invoke fInfoPad";
cerr << endl;
return;
}
@ -2929,9 +2994,9 @@ void PMusrCanvas::HandleDataSet(UInt_t plotNo, UInt_t runNo, PRunData *data)
Double_t dval = (startFitRange - data->GetDataTimeStart())/data->GetDataTimeStep();
if (dval < 0.0) { // make sure that startBin >= 0
startBin = 0;
cerr << endl << "PMusrCanvas::HandleDataSet() **WARNING** found startBin data < 0 for 'use_fit_range', will set it to 0" << endl << endl;
cerr << endl << ">> PMusrCanvas::HandleDataSet(): **WARNING** found startBin data < 0 for 'use_fit_range', will set it to 0" << endl << endl;
} else if (dval >= (Double_t)data->GetValue()->size()) { // make sure that startBin <= length of data vector
cerr << endl << ">> PMusrCanvas::HandleDataSet() **WARNING** found startBin data=" << (UInt_t)dval << " >= data vector size=" << data->GetValue()->size() << " for 'use_fit_range',";
cerr << endl << ">> PMusrCanvas::HandleDataSet(): **WARNING** found startBin data=" << (UInt_t)dval << " >= data vector size=" << data->GetValue()->size() << " for 'use_fit_range',";
cerr << endl << ">> will set it to data vector size" << endl << endl;
startBin = data->GetValue()->size();
} else {
@ -2945,9 +3010,9 @@ void PMusrCanvas::HandleDataSet(UInt_t plotNo, UInt_t runNo, PRunData *data)
dval = (endFitRange - data->GetDataTimeStart())/data->GetDataTimeStep();
if (dval < 0.0) { // make sure that endBin >= 0
endBin = 0;
cerr << endl << "PMusrCanvas::HandleDataSet() **WARNING** found endBin data < 0 for 'use_fit_range', will set it to 0" << endl << endl;
cerr << endl << ">> PMusrCanvas::HandleDataSet(): **WARNING** found endBin data < 0 for 'use_fit_range', will set it to 0" << endl << endl;
} else if (dval >= (Double_t)data->GetValue()->size()) { // make sure that endBin <= length of data vector
cerr << endl << ">> PMusrCanvas::HandleDataSet() **WARNING** found endBin data=" << (UInt_t)dval << " >= data vector size=" << data->GetValue()->size() << " for 'use_fit_range',";
cerr << endl << ">> PMusrCanvas::HandleDataSet(): **WARNING** found endBin data=" << (UInt_t)dval << " >= data vector size=" << data->GetValue()->size() << " for 'use_fit_range',";
cerr << endl << ">> will set it to data vector size" << endl << endl;
endBin = data->GetValue()->size();
} else {
@ -2960,9 +3025,9 @@ void PMusrCanvas::HandleDataSet(UInt_t plotNo, UInt_t runNo, PRunData *data)
Double_t dval = (fMsrHandler->GetMsrPlotList()->at(fPlotNumber).fTmin[runNo] - data->GetDataTimeStart())/data->GetDataTimeStep();
if (dval < 0.0) { // make sure that startBin >= 0
startBin = 0;
cerr << endl << "PMusrCanvas::HandleDataSet() **WARNING** found startBin data < 0 for 'sub_ranges', will set it to 0" << endl << endl;
cerr << endl << ">> PMusrCanvas::HandleDataSet(): **WARNING** found startBin data < 0 for 'sub_ranges', will set it to 0" << endl << endl;
} else if (dval >= (Double_t)data->GetValue()->size()) { // make sure that startBin <= length of data vector
cerr << endl << ">> PMusrCanvas::HandleDataSet() **WARNING** found startBin data=" << (UInt_t)dval << " >= data vector size=" << data->GetValue()->size() << " for 'sub_ranges',";
cerr << endl << ">> PMusrCanvas::HandleDataSet(): **WARNING** found startBin data=" << (UInt_t)dval << " >= data vector size=" << data->GetValue()->size() << " for 'sub_ranges',";
cerr << endl << ">> will set it to data vector size" << endl << endl;
startBin = data->GetValue()->size();
} else {
@ -2972,9 +3037,9 @@ void PMusrCanvas::HandleDataSet(UInt_t plotNo, UInt_t runNo, PRunData *data)
dval = (fMsrHandler->GetMsrPlotList()->at(fPlotNumber).fTmax[runNo] - data->GetDataTimeStart())/data->GetDataTimeStep();
if (dval < 0.0) { // make sure that endBin >= 0
endBin = 0;
cerr << endl << "PMusrCanvas::HandleDataSet() **WARNING** found endBin data < 0 for 'sub_ranges', will set it to 0" << endl << endl;
cerr << endl << ">> PMusrCanvas::HandleDataSet(): **WARNING** found endBin data < 0 for 'sub_ranges', will set it to 0" << endl << endl;
} else if (dval >= (Double_t)data->GetValue()->size()) { // make sure that endtBin <= length of data vector
cerr << endl << ">> PMusrCanvas::HandleDataSet() **WARNING** found endBin data=" << (UInt_t)dval << " >= data vector size=" << data->GetValue()->size() << " for 'sub_ranges',";
cerr << endl << ">> PMusrCanvas::HandleDataSet(): **WARNING** found endBin data=" << (UInt_t)dval << " >= data vector size=" << data->GetValue()->size() << " for 'sub_ranges',";
cerr << endl << ">> will set it to data vector size" << endl << endl;
endBin = data->GetValue()->size();
} else {
@ -3061,9 +3126,9 @@ void PMusrCanvas::HandleDataSet(UInt_t plotNo, UInt_t runNo, PRunData *data)
Double_t dval = (startFitRange - data->GetDataTimeStart())/data->GetTheoryTimeStep();
if (dval < 0.0) { // make sure that startBin >= 0
startBin = 0;
cerr << endl << "PMusrCanvas::HandleDataSet() **WARNING** found startBin theory < 0 for 'use_fit_range', will set it to 0" << endl << endl;
cerr << endl << ">> PMusrCanvas::HandleDataSet(): **WARNING** found startBin theory < 0 for 'use_fit_range', will set it to 0" << endl << endl;
} else if (dval >= (Double_t)data->GetTheory()->size()) { // make sure that startBin <= length of theory vector
cerr << endl << ">> PMusrCanvas::HandleDataSet() **WARNING** found startBin theory=" << (UInt_t)dval << " >= theory vector size=" << data->GetTheory()->size() << " for 'use_fit_range',";
cerr << endl << ">> PMusrCanvas::HandleDataSet(): **WARNING** found startBin theory=" << (UInt_t)dval << " >= theory vector size=" << data->GetTheory()->size() << " for 'use_fit_range',";
cerr << endl << ">> will set it to theory vector size" << endl << endl;
startBin = data->GetTheory()->size();
} else {
@ -3077,9 +3142,9 @@ void PMusrCanvas::HandleDataSet(UInt_t plotNo, UInt_t runNo, PRunData *data)
dval = (endFitRange - data->GetDataTimeStart())/data->GetTheoryTimeStep();
if (dval < 0.0) { // make sure that endBin >= 0
endBin = 0;
cerr << endl << "PMusrCanvas::HandleDataSet() **WARNING** found endBin theory < 0 for 'use_fit_range', will set it to 0" << endl << endl;
cerr << endl << ">> PMusrCanvas::HandleDataSet(): **WARNING** found endBin theory < 0 for 'use_fit_range', will set it to 0" << endl << endl;
} else if (dval >= (Double_t)data->GetTheory()->size()) { // make sure that endBin <= length of theory vector
cerr << endl << ">> PMusrCanvas::HandleDataSet() **WARNING** found endBin theory=" << (UInt_t)dval << " >= theory vector size=" << data->GetTheory()->size() << " for 'use_fit_range',";
cerr << endl << ">> PMusrCanvas::HandleDataSet(): **WARNING** found endBin theory=" << (UInt_t)dval << " >= theory vector size=" << data->GetTheory()->size() << " for 'use_fit_range',";
cerr << endl << ">> will set it to theory vector size" << endl << endl;
endBin = data->GetTheory()->size();
} else {
@ -3095,9 +3160,9 @@ void PMusrCanvas::HandleDataSet(UInt_t plotNo, UInt_t runNo, PRunData *data)
Double_t dval = (fMsrHandler->GetMsrPlotList()->at(fPlotNumber).fTmin[runNo] -data->GetDataTimeStart())/data->GetTheoryTimeStep();
if (dval < 0.0) { // make sure that startBin >= 0
startBin = 0;
cerr << endl << "PMusrCanvas::HandleDataSet() **WARNING** found startBin theory < 0 for 'sub_ranges', will set it to 0" << endl << endl;
cerr << endl << ">> PMusrCanvas::HandleDataSet(): **WARNING** found startBin theory < 0 for 'sub_ranges', will set it to 0" << endl << endl;
} else if (dval >= (Double_t)data->GetTheory()->size()) { // make sure that startBin <= length of theory vector
cerr << endl << ">> PMusrCanvas::HandleDataSet() **WARNING** found startBin theory=" << (UInt_t)dval << " >= theory vector size=" << data->GetTheory()->size() << " for 'sub_ranges',";
cerr << endl << ">> PMusrCanvas::HandleDataSet(): **WARNING** found startBin theory=" << (UInt_t)dval << " >= theory vector size=" << data->GetTheory()->size() << " for 'sub_ranges',";
cerr << endl << ">> will set it to theory vector size" << endl << endl;
startBin = data->GetTheory()->size();
} else {
@ -3107,9 +3172,9 @@ void PMusrCanvas::HandleDataSet(UInt_t plotNo, UInt_t runNo, PRunData *data)
dval = (fMsrHandler->GetMsrPlotList()->at(fPlotNumber).fTmax[runNo] -data->GetDataTimeStart())/data->GetTheoryTimeStep();
if (dval < 0.0) { // make sure that endBin >= 0
endBin = 0;
cerr << endl << "PMusrCanvas::HandleDataSet() **WARNING** found endBin theory < 0 for 'sub_ranges', will set it to 0" << endl << endl;
cerr << endl << ">> PMusrCanvas::HandleDataSet(): **WARNING** found endBin theory < 0 for 'sub_ranges', will set it to 0" << endl << endl;
} else if (dval >= (Double_t)data->GetTheory()->size()) { // make sure that endtBin <= length of theory vector
cerr << endl << ">> PMusrCanvas::HandleDataSet() **WARNING** found endBin theory=" << (UInt_t)dval << " >= theory vector size=" << data->GetTheory()->size() << " for 'sub_ranges',";
cerr << endl << ">> PMusrCanvas::HandleDataSet(): **WARNING** found endBin theory=" << (UInt_t)dval << " >= theory vector size=" << data->GetTheory()->size() << " for 'sub_ranges',";
cerr << endl << ">> will set it to theory vector size" << endl << endl;
endBin = data->GetTheory()->size();
} else {
@ -3418,7 +3483,7 @@ void PMusrCanvas::HandleFourier()
// calculate fourier transform of the data
PFourier fourierData(fData[i].data, fFourier.fUnits, startTime, endTime, fFourier.fDCCorrected, fFourier.fFourierPower);
if (!fourierData.IsValid()) {
cerr << endl << "**SEVERE ERROR** PMusrCanvas::HandleFourier: couldn't invoke PFourier to calculate the Fourier data ..." << endl;
cerr << endl << ">> PMusrCanvas::HandleFourier(): **SEVERE ERROR** couldn't invoke PFourier to calculate the Fourier data ..." << endl;
return;
}
fourierData.Transform(fFourier.fApodization);
@ -3458,7 +3523,7 @@ void PMusrCanvas::HandleFourier()
Int_t powerPad = (Int_t)round(log((endTime-startTime)/fData[i].theory->GetBinWidth(1))/log(2))+3;
PFourier fourierTheory(fData[i].theory, fFourier.fUnits, startTime, endTime, fFourier.fDCCorrected, powerPad);
if (!fourierTheory.IsValid()) {
cerr << endl << "**SEVERE ERROR** PMusrCanvas::HandleFourier: couldn't invoke PFourier to calculate the Fourier theory ..." << endl;
cerr << endl << ">> PMusrCanvas::HandleFourier(): **SEVERE ERROR** couldn't invoke PFourier to calculate the Fourier theory ..." << endl;
return;
}
fourierTheory.Transform(fFourier.fApodization);
@ -3578,7 +3643,7 @@ void PMusrCanvas::HandleDifferenceFourier()
// calculate fourier transform of the data
PFourier fourierData(fData[i].diff, fFourier.fUnits, startTime, endTime, fFourier.fDCCorrected, fFourier.fFourierPower);
if (!fourierData.IsValid()) {
cerr << endl << "**SEVERE ERROR** PMusrCanvas::HandleFourier: couldn't invoke PFourier to calculate the Fourier diff ..." << endl;
cerr << endl << ">> PMusrCanvas::HandleFourier(): **SEVERE ERROR** couldn't invoke PFourier to calculate the Fourier diff ..." << endl;
return;
}
fourierData.Transform(fFourier.fApodization);
@ -4648,14 +4713,14 @@ void PMusrCanvas::PlotData(Bool_t unzoom)
fDataTheoryPad->SetLogy(1);
// set x-axis label
fHistoFrame->GetXaxis()->SetTitle("time (#mus)");
fHistoFrame->GetXaxis()->SetTitle("Time (#mus)");
// set y-axis label
TString yAxisTitle;
PMsrRunList *runList = fMsrHandler->GetMsrRunList();
switch (fPlotType) {
case MSR_PLOT_SINGLE_HISTO:
if (runList->at(0).IsLifetimeCorrected()) { // lifetime correction
yAxisTitle = "asymmetry";
yAxisTitle = "Asymmetry";
} else { // no liftime correction
if (fScaleN0AndBkg)
yAxisTitle = "N(t) per nsec";
@ -4663,8 +4728,12 @@ void PMusrCanvas::PlotData(Bool_t unzoom)
yAxisTitle = "N(t) per bin";
}
break;
case MSR_PLOT_SINGLE_HISTO_RRF:
case MSR_PLOT_ASYM_RRF:
yAxisTitle = "RRF Asymmetry";
break;
case MSR_PLOT_ASYM:
yAxisTitle = "asymmetry";
yAxisTitle = "Asymmetry";
break;
case MSR_PLOT_MU_MINUS:
yAxisTitle = "N(t) per bin";
@ -6106,6 +6175,11 @@ void PMusrCanvas::PlotAverage(Bool_t unzoom)
break;
}
// check if RRF and if yes show a label
if ((fRRFText != 0) && (fRRFLatexText != 0)) {
fRRFLatexText->DrawLatex(0.1, 0.92, fRRFText->Data());
}
fDataTheoryPad->Update();
fMainCanvas->cd();

View File

@ -8,7 +8,7 @@
***************************************************************************/
/***************************************************************************
* Copyright (C) 2007-2014 by Andreas Suter *
* Copyright (C) 2007-2016 by Andreas Suter *
* andreas.suter@psi.ch *
* *
* This program is free software; you can redistribute it and/or modify *

View File

@ -8,7 +8,7 @@
***************************************************************************/
/***************************************************************************
* Copyright (C) 2007-2015 by Andreas Suter *
* Copyright (C) 2007-2016 by Andreas Suter *
* andreas.suter@psi.ch *
* *
* This program is free software; you can redistribute it and/or modify *
@ -181,7 +181,7 @@ void PPrepFourier::DoBkgCorrection()
if ((fBkgRange[0] != -1) && (fBkgRange[1] != -1)) { // background range is given
// make sure that the bkg range is ok
for (UInt_t i=0; i<fRawData.size(); i++) {
if ((fBkgRange[0] >= fRawData[i].rawData.size()) || (fBkgRange[1] >= fRawData[i].rawData.size())) {
if ((fBkgRange[0] >= (Int_t)fRawData[i].rawData.size()) || (fBkgRange[1] >= (Int_t)fRawData[i].rawData.size())) {
cerr << endl << "PPrepFourier::DoBkgCorrection() **ERROR** bkg-range out of data-range!";
return;
}

View File

@ -8,7 +8,7 @@
***************************************************************************/
/***************************************************************************
* Copyright (C) 2007-2014 by Andreas Suter *
* Copyright (C) 2007-2016 by Andreas Suter *
* andreas.suter@psi.ch *
* *
* This program is free software; you can redistribute it and/or modify *
@ -87,7 +87,7 @@ PRunAsymmetry::PRunAsymmetry(PMsrHandler *msrInfo, PRunDataHandler *rawData, UIn
fPacking = fMsrInfo->GetMsrGlobal()->GetPacking();
}
if (fPacking == -1) { // this should NOT happen, somethin is severely wrong
cerr << endl << ">> PRunAsymmetry::PRunAsymmetry: **SEVERE ERROR**: Couldn't find any packing information!";
cerr << endl << ">> PRunAsymmetry::PRunAsymmetry(): **SEVERE ERROR**: Couldn't find any packing information!";
cerr << endl << ">> This is very bad :-(, will quit ...";
cerr << endl;
fValid = false;
@ -342,7 +342,7 @@ void PRunAsymmetry::SetFitRangeBin(const TString fitRange)
Int_t pos = 2*(fRunNo+1)-1;
if (pos + 1 >= tok->GetEntries()) {
cerr << endl << ">> PRunSingleHisto::SetFitRangeBin(): **ERROR** invalid FIT_RANGE command found: '" << fitRange << "'";
cerr << endl << ">> PRunAsymmetry::SetFitRangeBin(): **ERROR** invalid FIT_RANGE command found: '" << fitRange << "'";
cerr << endl << ">> will ignore it. Sorry ..." << endl;
} else {
// handle fgb+n0 entry
@ -370,7 +370,7 @@ void PRunAsymmetry::SetFitRangeBin(const TString fitRange)
fFitEndTime = (fGoodBins[1] - offset - fT0s[0]) * fTimeResolution;
}
} else { // error
cerr << endl << ">> PRunSingleHisto::SetFitRangeBin(): **ERROR** invalid FIT_RANGE command found: '" << fitRange << "'";
cerr << endl << ">> PRunAsymmetry::SetFitRangeBin(): **ERROR** invalid FIT_RANGE command found: '" << fitRange << "'";
cerr << endl << ">> will ignore it. Sorry ..." << endl;
}

File diff suppressed because it is too large Load Diff

View File

@ -8,7 +8,7 @@
***************************************************************************/
/***************************************************************************
* Copyright (C) 2007-2014 by Andreas Suter *
* Copyright (C) 2007-2016 by Andreas Suter *
* andreas.suter@psi.ch *
* *
* This program is free software; you can redistribute it and/or modify *

View File

@ -8,7 +8,7 @@
***************************************************************************/
/***************************************************************************
* Copyright (C) 2007-2014 by Andreas Suter *
* Copyright (C) 2007-2016 by Andreas Suter *
* andreas.suter@psi.ch *
* *
* This program is free software; you can redistribute it and/or modify *

View File

@ -8,7 +8,7 @@
***************************************************************************/
/***************************************************************************
* Copyright (C) 2007-2014 by Andreas Suter *
* Copyright (C) 2007-2016 by Andreas Suter *
* andreas.suter@psi.ch *
* *
* This program is free software; you can redistribute it and/or modify *
@ -58,12 +58,24 @@ PRunListCollection::~PRunListCollection()
}
fRunSingleHistoList.clear();
for (UInt_t i=0; i<fRunSingleHistoRRFList.size(); i++) {
fRunSingleHistoRRFList[i]->CleanUp();
fRunSingleHistoRRFList[i]->~PRunSingleHistoRRF();
}
fRunSingleHistoRRFList.clear();
for (UInt_t i=0; i<fRunAsymmetryList.size(); i++) {
fRunAsymmetryList[i]->CleanUp();
fRunAsymmetryList[i]->~PRunAsymmetry();
}
fRunAsymmetryList.clear();
for (UInt_t i=0; i<fRunAsymmetryRRFList.size(); i++) {
fRunAsymmetryRRFList[i]->CleanUp();
fRunAsymmetryRRFList[i]->~PRunAsymmetryRRF();
}
fRunAsymmetryRRFList.clear();
for (UInt_t i=0; i<fRunMuMinusList.size(); i++) {
fRunMuMinusList[i]->CleanUp();
fRunMuMinusList[i]->~PRunMuMinus();
@ -106,11 +118,21 @@ Bool_t PRunListCollection::Add(Int_t runNo, EPMusrHandleTag tag)
if (!fRunSingleHistoList[fRunSingleHistoList.size()-1]->IsValid())
success = false;
break;
case PRUN_SINGLE_HISTO_RRF:
fRunSingleHistoRRFList.push_back(new PRunSingleHistoRRF(fMsrInfo, fData, runNo, tag));
if (!fRunSingleHistoRRFList[fRunSingleHistoRRFList.size()-1]->IsValid())
success = false;
break;
case PRUN_ASYMMETRY:
fRunAsymmetryList.push_back(new PRunAsymmetry(fMsrInfo, fData, runNo, tag));
if (!fRunAsymmetryList[fRunAsymmetryList.size()-1]->IsValid())
success = false;
break;
case PRUN_ASYMMETRY_RRF:
fRunAsymmetryRRFList.push_back(new PRunAsymmetryRRF(fMsrInfo, fData, runNo, tag));
if (!fRunAsymmetryRRFList[fRunAsymmetryRRFList.size()-1]->IsValid())
success = false;
break;
case PRUN_MU_MINUS:
fRunMuMinusList.push_back(new PRunMuMinus(fMsrInfo, fData, runNo, tag));
if (!fRunMuMinusList[fRunMuMinusList.size()-1]->IsValid())
@ -147,8 +169,12 @@ void PRunListCollection::SetFitRange(const TString fitRange)
{
for (UInt_t i=0; i<fRunSingleHistoList.size(); i++)
fRunSingleHistoList[i]->SetFitRangeBin(fitRange);
for (UInt_t i=0; i<fRunSingleHistoRRFList.size(); i++)
fRunSingleHistoRRFList[i]->SetFitRangeBin(fitRange);
for (UInt_t i=0; i<fRunAsymmetryList.size(); i++)
fRunAsymmetryList[i]->SetFitRangeBin(fitRange);
for (UInt_t i=0; i<fRunAsymmetryRRFList.size(); i++)
fRunAsymmetryRRFList[i]->SetFitRangeBin(fitRange);
for (UInt_t i=0; i<fRunMuMinusList.size(); i++)
fRunMuMinusList[i]->SetFitRangeBin(fitRange);
for (UInt_t i=0; i<fRunNonMusrList.size(); i++)
@ -169,8 +195,12 @@ void PRunListCollection::SetFitRange(const PDoublePairVector fitRange)
{
for (UInt_t i=0; i<fRunSingleHistoList.size(); i++)
fRunSingleHistoList[i]->SetFitRange(fitRange);
for (UInt_t i=0; i<fRunSingleHistoRRFList.size(); i++)
fRunSingleHistoRRFList[i]->SetFitRange(fitRange);
for (UInt_t i=0; i<fRunAsymmetryList.size(); i++)
fRunAsymmetryList[i]->SetFitRange(fitRange);
for (UInt_t i=0; i<fRunAsymmetryRRFList.size(); i++)
fRunAsymmetryRRFList[i]->SetFitRange(fitRange);
for (UInt_t i=0; i<fRunMuMinusList.size(); i++)
fRunMuMinusList[i]->SetFitRange(fitRange);
for (UInt_t i=0; i<fRunNonMusrList.size(); i++)
@ -198,6 +228,27 @@ Double_t PRunListCollection::GetSingleHistoChisq(const std::vector<Double_t>& pa
return chisq;
}
//--------------------------------------------------------------------------
// GetSingleHistoRRFChisq (public)
//--------------------------------------------------------------------------
/**
* <p>Calculates chi-square of <em>all</em> single histogram RRF runs of a msr-file.
*
* <b>return:</b>
* - chi-square of all single histogram RRF runs of the msr-file
*
* \param par fit parameter vector
*/
Double_t PRunListCollection::GetSingleHistoRRFChisq(const std::vector<Double_t>& par) const
{
Double_t chisq = 0.0;
for (UInt_t i=0; i<fRunSingleHistoRRFList.size(); i++)
chisq += fRunSingleHistoRRFList[i]->CalcChiSquare(par);
return chisq;
}
//--------------------------------------------------------------------------
// GetAsymmetryChisq (public)
//--------------------------------------------------------------------------
@ -219,6 +270,27 @@ Double_t PRunListCollection::GetAsymmetryChisq(const std::vector<Double_t>& par)
return chisq;
}
//--------------------------------------------------------------------------
// GetAsymmetryRRFChisq (public)
//--------------------------------------------------------------------------
/**
* <p>Calculates chi-square of <em>all</em> asymmetry RRF runs of a msr-file.
*
* <b>return:</b>
* - chi-square of all asymmetry RRF runs of the msr-file
*
* \param par fit parameter vector
*/
Double_t PRunListCollection::GetAsymmetryRRFChisq(const std::vector<Double_t>& par) const
{
Double_t chisq = 0.0;
for (UInt_t i=0; i<fRunAsymmetryRRFList.size(); i++)
chisq += fRunAsymmetryRRFList[i]->CalcChiSquare(par);
return chisq;
}
//--------------------------------------------------------------------------
// GetMuMinusChisq (public)
//--------------------------------------------------------------------------
@ -299,9 +371,15 @@ Double_t PRunListCollection::GetSingleHistoChisqExpected(const std::vector<Doubl
case PRUN_SINGLE_HISTO:
expectedChisq = fRunSingleHistoList[subIdx]->CalcChiSquareExpected(par);
break;
case PRUN_SINGLE_HISTO_RRF:
expectedChisq = fRunSingleHistoRRFList[subIdx]->CalcChiSquareExpected(par);
break;
case PRUN_ASYMMETRY:
expectedChisq = fRunAsymmetryList[subIdx]->CalcChiSquareExpected(par);
break;
case PRUN_ASYMMETRY_RRF:
expectedChisq = fRunAsymmetryRRFList[subIdx]->CalcChiSquareExpected(par);
break;
case PRUN_MU_MINUS:
expectedChisq = fRunMuMinusList[subIdx]->CalcChiSquareExpected(par);
break;
@ -336,16 +414,17 @@ Double_t PRunListCollection::GetSingleRunChisq(const std::vector<Double_t>& par,
return chisq;
}
Int_t subIdx = 0;
Int_t type = fMsrInfo->GetMsrRunList()->at(idx).GetFitType();
if (type == -1) { // i.e. not forun in the RUN block, try the GLOBAL block
if (type == -1) { // i.e. not found in the RUN block, try the GLOBAL block
type = fMsrInfo->GetMsrGlobal()->GetFitType();
}
// count how many entries of this fit-type are present up to idx
UInt_t subIdx = 0;
for (UInt_t i=0; i<idx; i++) {
if (fMsrInfo->GetMsrRunList()->at(i).GetFitType() == type)
subIdx++;
subIdx = idx;
} else { // found in the RUN block
// count how many entries of this fit-type are present up to idx
for (UInt_t i=0; i<idx; i++) {
if (fMsrInfo->GetMsrRunList()->at(i).GetFitType() == type)
subIdx++;
}
}
// return the chisq of the single run
@ -353,9 +432,15 @@ Double_t PRunListCollection::GetSingleRunChisq(const std::vector<Double_t>& par,
case PRUN_SINGLE_HISTO:
chisq = fRunSingleHistoList[subIdx]->CalcChiSquare(par);
break;
case PRUN_SINGLE_HISTO_RRF:
chisq = fRunSingleHistoRRFList[subIdx]->CalcChiSquare(par);
break;
case PRUN_ASYMMETRY:
chisq = fRunAsymmetryList[subIdx]->CalcChiSquare(par);
break;
case PRUN_ASYMMETRY_RRF:
chisq = fRunAsymmetryRRFList[subIdx]->CalcChiSquare(par);
break;
case PRUN_MU_MINUS:
chisq = fRunMuMinusList[subIdx]->CalcChiSquare(par);
break;
@ -376,7 +461,7 @@ Double_t PRunListCollection::GetSingleRunChisq(const std::vector<Double_t>& par,
* <p>Calculates log max-likelihood of <em>all</em> single histogram runs of a msr-file.
*
* <b>return:</b>
* - chi-square of all single histogram runs of the msr-file
* - log max-likelihood of all single histogram runs of the msr-file
*
* \param par fit parameter vector
*/
@ -390,6 +475,27 @@ Double_t PRunListCollection::GetSingleHistoMaximumLikelihood(const std::vector<D
return mlh;
}
//--------------------------------------------------------------------------
// GetSingleHistoRRFMaximumLikelihood (public)
//--------------------------------------------------------------------------
/**
* <p>Calculates log max-likelihood of <em>all</em> single histogram RRF runs of a msr-file.
*
* <b>return:</b>
* - log max-likelihood of all single histogram runs of the msr-file
*
* \param par fit parameter vector
*/
Double_t PRunListCollection::GetSingleHistoRRFMaximumLikelihood(const std::vector<Double_t>& par) const
{
Double_t mlh = 0.0;
for (UInt_t i=0; i<fRunSingleHistoRRFList.size(); i++)
mlh += fRunSingleHistoRRFList[i]->CalcMaxLikelihood(par);
return mlh;
}
//--------------------------------------------------------------------------
// GetAsymmetryMaximumLikelihood (public)
//--------------------------------------------------------------------------
@ -412,6 +518,28 @@ Double_t PRunListCollection::GetAsymmetryMaximumLikelihood(const std::vector<Dou
return mlh;
}
//--------------------------------------------------------------------------
// GetAsymmetryRRFMaximumLikelihood (public)
//--------------------------------------------------------------------------
/**
* <p> Since it is not clear yet how to handle asymmetry fits with max likelihood
* the chi square will be used!
*
* <b>return:</b>
* - chi-square of all asymmetry RRF runs of the msr-file
*
* \param par fit parameter vector
*/
Double_t PRunListCollection::GetAsymmetryRRFMaximumLikelihood(const std::vector<Double_t>& par) const
{
Double_t mlh = 0.0;
for (UInt_t i=0; i<fRunAsymmetryRRFList.size(); i++)
mlh += fRunAsymmetryRRFList[i]->CalcChiSquare(par);
return mlh;
}
//--------------------------------------------------------------------------
// GetMuMinusMaximumLikelihood (public)
//--------------------------------------------------------------------------
@ -419,7 +547,7 @@ Double_t PRunListCollection::GetAsymmetryMaximumLikelihood(const std::vector<Dou
* <p>Calculates log max-likelihood of <em>all</em> mu minus runs of a msr-file.
*
* <b>return:</b>
* - chi-square of all mu minus runs of the msr-file
* - log max-likelihood of all mu minus runs of the msr-file
*
* \param par fit parameter vector
*/
@ -493,9 +621,15 @@ UInt_t PRunListCollection::GetNoOfBinsFitted(const UInt_t idx) const
case PRUN_SINGLE_HISTO:
result = fRunSingleHistoList[subIdx]->GetNoOfFitBins();
break;
case PRUN_SINGLE_HISTO_RRF:
result = fRunSingleHistoRRFList[subIdx]->GetNoOfFitBins();
break;
case PRUN_ASYMMETRY:
result = fRunAsymmetryList[subIdx]->GetNoOfFitBins();
break;
case PRUN_ASYMMETRY_RRF:
result = fRunAsymmetryRRFList[subIdx]->GetNoOfFitBins();
break;
case PRUN_MU_MINUS:
result = fRunMuMinusList[subIdx]->GetNoOfFitBins();
break;
@ -526,9 +660,15 @@ UInt_t PRunListCollection::GetTotalNoOfBinsFitted() const
for (UInt_t i=0; i<fRunSingleHistoList.size(); i++)
counts += fRunSingleHistoList[i]->GetNoOfFitBins();
for (UInt_t i=0; i<fRunSingleHistoRRFList.size(); i++)
counts += fRunSingleHistoRRFList[i]->GetNoOfFitBins();
for (UInt_t i=0; i<fRunAsymmetryList.size(); i++)
counts += fRunAsymmetryList[i]->GetNoOfFitBins();
for (UInt_t i=0; i<fRunAsymmetryRRFList.size(); i++)
counts += fRunAsymmetryRRFList[i]->GetNoOfFitBins();
for (UInt_t i=0; i<fRunMuMinusList.size(); i++)
counts += fRunMuMinusList[i]->GetNoOfFitBins();
@ -558,7 +698,7 @@ PRunData* PRunListCollection::GetSingleHisto(UInt_t index, EDataSwitch tag)
switch (tag) {
case kIndex:
if ((index < 0) || (index >= fRunSingleHistoList.size())) {
cerr << endl << "PRunListCollection::GetSingleHisto: **ERROR** index = " << index << " out of bounds";
cerr << endl << ">> PRunListCollection::GetSingleHisto(): **ERROR** index = " << index << " out of bounds";
cerr << endl;
return 0;
}
@ -581,6 +721,49 @@ PRunData* PRunListCollection::GetSingleHisto(UInt_t index, EDataSwitch tag)
return data;
}
//--------------------------------------------------------------------------
// GetSingleHistoRRF (public)
//--------------------------------------------------------------------------
/**
* <p>Get a processed single histogram RRF data set.
*
* <b>return:</b>
* - pointer to the run data set (processed data) if data set is found
* - null pointer otherwise
*
* \param index msr-file run index
* \param tag kIndex -> data at index, kRunNo -> data of given run no
*/
PRunData* PRunListCollection::GetSingleHistoRRF(UInt_t index, EDataSwitch tag)
{
PRunData *data = 0;
switch (tag) {
case kIndex:
if ((index < 0) || (index >= fRunSingleHistoRRFList.size())) {
cerr << endl << ">> PRunListCollection::GetSingleHistoRRF(): **ERROR** index = " << index << " out of bounds";
cerr << endl;
return 0;
}
fRunSingleHistoRRFList[index]->CalcTheory();
data = fRunSingleHistoRRFList[index]->GetData();
break;
case kRunNo:
for (UInt_t i=0; i<fRunSingleHistoRRFList.size(); i++) {
if (fRunSingleHistoRRFList[i]->GetRunNo() == index) {
data = fRunSingleHistoRRFList[i]->GetData();
break;
}
}
break;
default: // error
break;
}
return data;
}
//--------------------------------------------------------------------------
// GetAsymmetry (public)
//--------------------------------------------------------------------------
@ -601,7 +784,7 @@ PRunData* PRunListCollection::GetAsymmetry(UInt_t index, EDataSwitch tag)
switch (tag) {
case kIndex: // called from musrfit when dumping the data
if ((index < 0) || (index > fRunAsymmetryList.size())) {
cerr << endl << "PRunListCollection::GetAsymmetry: **ERROR** index = " << index << " out of bounds";
cerr << endl << ">> PRunListCollection::GetAsymmetry(): **ERROR** index = " << index << " out of bounds";
cerr << endl;
return 0;
}
@ -624,6 +807,49 @@ PRunData* PRunListCollection::GetAsymmetry(UInt_t index, EDataSwitch tag)
return data;
}
//--------------------------------------------------------------------------
// GetAsymmetryRRF (public)
//--------------------------------------------------------------------------
/**
* <p>Get a processed asymmetry RRF data set.
*
* <b>return:</b>
* - pointer to the run data set (processed data) if data set is found
* - null pointer otherwise
*
* \param index msr-file run index
* \param tag kIndex -> data at index, kRunNo -> data of given run no
*/
PRunData* PRunListCollection::GetAsymmetryRRF(UInt_t index, EDataSwitch tag)
{
PRunData *data = 0;
switch (tag) {
case kIndex: // called from musrfit when dumping the data
if ((index < 0) || (index > fRunAsymmetryRRFList.size())) {
cerr << endl << ">> PRunListCollection::GetAsymmetryRRF(): **ERROR** index = " << index << " out of bounds";
cerr << endl;
return 0;
}
fRunAsymmetryRRFList[index]->CalcTheory();
data = fRunAsymmetryRRFList[index]->GetData();
break;
case kRunNo: // called from PMusrCanvas
for (UInt_t i=0; i<fRunAsymmetryRRFList.size(); i++) {
if (fRunAsymmetryRRFList[i]->GetRunNo() == index) {
data = fRunAsymmetryRRFList[i]->GetData();
break;
}
}
break;
default: // error
break;
}
return data;
}
//--------------------------------------------------------------------------
// GetMuMinus (public)
//--------------------------------------------------------------------------
@ -644,7 +870,7 @@ PRunData* PRunListCollection::GetMuMinus(UInt_t index, EDataSwitch tag)
switch (tag) {
case kIndex:
if ((index < 0) || (index > fRunMuMinusList.size())) {
cerr << endl << "PRunListCollection::GetMuMinus: **ERROR** index = " << index << " out of bounds";
cerr << endl << ">> PRunListCollection::GetMuMinus(): **ERROR** index = " << index << " out of bounds";
cerr << endl;
return 0;
}
@ -686,7 +912,7 @@ PRunData* PRunListCollection::GetNonMusr(UInt_t index, EDataSwitch tag)
switch (tag) {
case kIndex:
if ((index < 0) || (index > fRunNonMusrList.size())) {
cerr << endl << "PRunListCollection::GetNonMusr: **ERROR** index = " << index << " out of bounds";
cerr << endl << ">> PRunListCollection::GetNonMusr(): **ERROR** index = " << index << " out of bounds";
cerr << endl;
return 0;
}

View File

@ -8,7 +8,7 @@
***************************************************************************/
/***************************************************************************
* Copyright (C) 2007-2014 by Andreas Suter *
* Copyright (C) 2007-2016 by Andreas Suter *
* andreas.suter@psi.ch *
* *
* This program is free software; you can redistribute it and/or modify *

View File

@ -8,7 +8,7 @@
***************************************************************************/
/***************************************************************************
* Copyright (C) 2007-2014 by Andreas Suter *
* Copyright (C) 2007-2016 by Andreas Suter *
* andreas.suter@psi.ch *
* *
* This program is free software; you can redistribute it and/or modify *

View File

@ -8,7 +8,7 @@
***************************************************************************/
/***************************************************************************
* Copyright (C) 2007-2014 by Andreas Suter *
* Copyright (C) 2007-2016 by Andreas Suter *
* andreas.suter@psi.ch *
* *
* This program is free software; you can redistribute it and/or modify *

File diff suppressed because it is too large Load Diff

View File

@ -8,7 +8,7 @@
***************************************************************************/
/***************************************************************************
* Copyright (C) 2007-2014 by Andreas Suter *
* Copyright (C) 2007-2016 by Andreas Suter *
* andreas.suter@psi.ch *
* *
* This program is free software; you can redistribute it and/or modify *

View File

@ -8,7 +8,7 @@
***************************************************************************/
/***************************************************************************
* Copyright (C) 2007-2014 by Andreas Suter *
* Copyright (C) 2007-2016 by Andreas Suter *
* andreas.suter@psi.ch *
* *
* This program is free software; you can redistribute it and/or modify *

View File

@ -8,7 +8,7 @@
***************************************************************************/
/***************************************************************************
* Copyright (C) 2007-2014 by Andreas Suter *
* Copyright (C) 2007-2016 by Andreas Suter *
* andreas.suter@psi.ch *
* *
* This program is free software; you can redistribute it and/or modify *

View File

@ -8,7 +8,7 @@
***************************************************************************/
/***************************************************************************
* Copyright (C) 2007-2014 by Andreas Suter *
* Copyright (C) 2007-2016 by Andreas Suter *
* andreas.suter@psi.ch *
* *
* This program is free software; you can redistribute it and/or modify *

View File

@ -8,7 +8,7 @@
***************************************************************************/
/***************************************************************************
* Copyright (C) 2007-2014 by Andreas Suter *
* Copyright (C) 2007-2016 by Andreas Suter *
* andreas.suter@psi.ch *
* *
* This program is free software; you can redistribute it and/or modify *

View File

@ -8,7 +8,7 @@
***************************************************************************/
/***************************************************************************
* Copyright (C) 2007-2014 by Andreas Suter *
* Copyright (C) 2007-2016 by Andreas Suter *
* andreas.suter@psi.ch *
* *
* This program is free software; you can redistribute it and/or modify *

View File

@ -8,7 +8,7 @@
***************************************************************************/
/***************************************************************************
* Copyright (C) 2007-2015 by Andreas Suter *
* Copyright (C) 2007-2016 by Andreas Suter *
* andreas.suter@psi.ch *
* *
* This program is free software; you can redistribute it and/or modify *

View File

@ -1,5 +1,14 @@
/***************************************************************************
* Copyright (C) 2007-2015 by Andreas Suter *
PFourierCanvas.h
Author: Andreas Suter
e-mail: andreas.suter@psi.ch
***************************************************************************/
/***************************************************************************
* Copyright (C) 2007-2016 by Andreas Suter *
* andreas.suter@psi.ch *
* *
* This program is free software; you can redistribute it and/or modify *

View File

@ -8,7 +8,7 @@
***************************************************************************/
/***************************************************************************
* Copyright (C) 2007-2015 by Andreas Suter *
* Copyright (C) 2007-2016 by Andreas Suter *
* andreas.suter@psi.ch *
* *
* This program is free software; you can redistribute it and/or modify *

View File

@ -8,7 +8,7 @@
***************************************************************************/
/***************************************************************************
* Copyright (C) 2007-2014 by Andreas Suter *
* Copyright (C) 2007-2016 by Andreas Suter *
* andreas.suter@psi.ch *
* *
* This program is free software; you can redistribute it and/or modify *

View File

@ -8,7 +8,7 @@
***************************************************************************/
/***************************************************************************
* Copyright (C) 2007-2014 by Andreas Suter *
* Copyright (C) 2007-2016 by Andreas Suter *
* andreas.suter@psi.ch *
* *
* This program is free software; you can redistribute it and/or modify *

View File

@ -8,7 +8,7 @@
***************************************************************************/
/***************************************************************************
* Copyright (C) 2007-2014 by Andreas Suter *
* Copyright (C) 2007-2016 by Andreas Suter *
* andreas.suter@psi.ch *
* *
* This program is free software; you can redistribute it and/or modify *

View File

@ -8,7 +8,7 @@
***************************************************************************/
/***************************************************************************
* Copyright (C) 2009-2014 by Bastian M. Wojek / Andreas Suter *
* Copyright (C) 2009-2016 by Bastian M. Wojek / Andreas Suter *
* andreas.suter@psi.ch *
* *
* This program is free software; you can redistribute it and/or modify *

View File

@ -8,7 +8,7 @@
***************************************************************************/
/***************************************************************************
* Copyright (C) 2007-2014 by Andreas Suter *
* Copyright (C) 2007-2016 by Andreas Suter *
* andreas.suter@psi.ch *
* *
* This program is free software; you can redistribute it and/or modify *
@ -102,6 +102,7 @@ class PMsrHandler
virtual Bool_t CheckFuncs();
virtual Bool_t CheckHistoGrouping();
virtual Bool_t CheckAddRunParameters();
virtual Bool_t CheckRRFSettings();
virtual void CheckMaxLikelihood();
virtual void GetGroupingString(Int_t runNo, TString detector, TString &groupingStr);

View File

@ -8,7 +8,7 @@
***************************************************************************/
/***************************************************************************
* Copyright (C) 2007-2015 by Andreas Suter *
* Copyright (C) 2007-2016 by Andreas Suter *
* andreas.suter@psi.ch *
* *
* This program is free software; you can redistribute it and/or modify *
@ -52,10 +52,12 @@ typedef struct { char a[7]; } __float128; // needed since cint doesn't know it
#define PMUSR_MSR_FILE_WRITE_ERROR -7
#define PMUSR_DATA_FILE_READ_ERROR -8
#define PRUN_SINGLE_HISTO 0
#define PRUN_ASYMMETRY 2
#define PRUN_MU_MINUS 4
#define PRUN_NON_MUSR 8
#define PRUN_SINGLE_HISTO 0
#define PRUN_SINGLE_HISTO_RRF 1
#define PRUN_ASYMMETRY 2
#define PRUN_ASYMMETRY_RRF 3
#define PRUN_MU_MINUS 4
#define PRUN_NON_MUSR 8
// muon life time in (us), see PRL99, 032001 (2007)
#define PMUON_LIFETIME 2.197019
@ -92,17 +94,21 @@ typedef struct { char a[7]; } __float128; // needed since cint doesn't know it
//-------------------------------------------------------------
// msr fit type tags
#define MSR_FITTYPE_SINGLE_HISTO 0
#define MSR_FITTYPE_ASYM 2
#define MSR_FITTYPE_MU_MINUS 4
#define MSR_FITTYPE_NON_MUSR 8
#define MSR_FITTYPE_SINGLE_HISTO 0
#define MSR_FITTYPE_SINGLE_HISTO_RRF 1
#define MSR_FITTYPE_ASYM 2
#define MSR_FITTYPE_ASYM_RRF 3
#define MSR_FITTYPE_MU_MINUS 4
#define MSR_FITTYPE_NON_MUSR 8
//-------------------------------------------------------------
// msr plot type tags
#define MSR_PLOT_SINGLE_HISTO 0
#define MSR_PLOT_ASYM 2
#define MSR_PLOT_MU_MINUS 4
#define MSR_PLOT_NON_MUSR 8
#define MSR_PLOT_SINGLE_HISTO 0
#define MSR_PLOT_SINGLE_HISTO_RRF 1
#define MSR_PLOT_ASYM 2
#define MSR_PLOT_ASYM_RRF 3
#define MSR_PLOT_MU_MINUS 4
#define MSR_PLOT_NON_MUSR 8
//-------------------------------------------------------------
// map and fun offsets for parameter parsing
@ -133,11 +139,14 @@ typedef struct { char a[7]; } __float128; // needed since cint doesn't know it
//-------------------------------------------------------------
// RRF related tags
#define RRF_UNIT_kHz 0
#define RRF_UNIT_MHz 1
#define RRF_UNIT_Mcs 2
#define RRF_UNIT_G 3
#define RRF_UNIT_T 4
#define RRF_UNIT_UNDEF -1
#define RRF_UNIT_kHz 0
#define RRF_UNIT_MHz 1
#define RRF_UNIT_Mcs 2
#define RRF_UNIT_G 3
#define RRF_UNIT_T 4
#define RRF_FREQ_UNDEF 1.0e10
//-------------------------------------------------------------
/**
@ -542,6 +551,11 @@ class PMsrGlobalBlock {
virtual ~PMsrGlobalBlock() {}
virtual Bool_t IsPresent() { return fGlobalPresent; }
virtual Double_t GetRRFFreq(const char *unit);
virtual TString GetRRFUnit();
virtual Int_t GetRRFUnitTag() { return fRRFUnitTag; }
virtual Double_t GetRRFPhase() { return fRRFPhase; }
virtual Int_t GetRRFPacking() { return fRRFPacking; }
virtual Int_t GetFitType() { return fFitType; }
virtual Int_t GetDataRange(UInt_t idx);
virtual UInt_t GetT0BinSize() { return fT0.size(); }
@ -555,6 +569,9 @@ class PMsrGlobalBlock {
virtual Int_t GetPacking() { return fPacking; }
virtual void SetGlobalPresent(Bool_t bval) { fGlobalPresent = bval; }
virtual void SetRRFFreq(Double_t freq, const char *unit);
virtual void SetRRFPhase(Double_t phase) { fRRFPhase = phase; }
virtual void SetRRFPacking(Int_t pack);
virtual void SetFitType(Int_t ival) { fFitType = ival; }
virtual void SetDataRange(Int_t ival, Int_t idx);
virtual void SetT0Bin(Double_t dval, Int_t idx=-1);
@ -566,10 +583,14 @@ class PMsrGlobalBlock {
private:
Bool_t fGlobalPresent; ///< flag showing if a GLOBAL block is present at all.
Int_t fFitType; ///< fit type: 0=single histo fit, 2=asymmetry fit, 4=mu^- single histo fit, 8=non muSR fit
Int_t fDataRange[4]; ///< data bin range (fit type 0, 2, 4)
PDoubleVector fT0; ///< t0 bins (fit type 0, 2, 4). if fit type 0 -> f0, f1, f2, ...; if fit type 2, 4 -> f0, b0, f1, b1, ...
vector<PDoubleVector> fAddT0; ///< addt0 bins (fit type 0, 2, 4). if fit type 0 -> f0, f1, f2, ...; if fit type 2, 4 -> f0, b0, f1, b1, ...
Double_t fRRFFreq; ///< RRF frequency given in units of (MHz, Mc, T)
Int_t fRRFUnitTag; ///< RRF unit tag
Double_t fRRFPhase; ///< RRF phase in (°)
Int_t fRRFPacking; ///< RRF packing
Int_t fFitType; ///< fit type: 0=single histo fit, 1=single histo RRF fit, 2=asymmetry fit, 4=mu^- single histo fit, 8=non muSR fit
Int_t fDataRange[4]; ///< data bin range (fit type 0, 1, 2, 4)
PDoubleVector fT0; ///< t0 bins (fit type 0, 1, 2, 4). if fit type 0 -> f0, f1, f2, ...; if fit type 2, 4 -> f0, b0, f1, b1, ...
vector<PDoubleVector> fAddT0; ///< addt0 bins (fit type 0, 1, 2, 4). if fit type 0 -> f0, f1, f2, ...; if fit type 2, 4 -> f0, b0, f1, b1, ...
Bool_t fFitRangeInBins; ///< flag telling if fit range is given in time or in bins
Double_t fFitRange[2]; ///< fit range in (us)
Int_t fFitRangeOffset[2]; ///< if fit range is given in bins it can have the form fit fgb+n0 lgb-n1. This variable holds the n0 and n1.

View File

@ -8,7 +8,7 @@
***************************************************************************/
/***************************************************************************
* Copyright (C) 2007-2014 by Andreas Suter *
* Copyright (C) 2007-2016 by Andreas Suter *
* andreas.suter@psi.ch *
* *
* This program is free software; you can redistribute it and/or modify *

View File

@ -8,7 +8,7 @@
***************************************************************************/
/***************************************************************************
* Copyright (C) 2007-2014 by Andreas Suter *
* Copyright (C) 2007-2016 by Andreas Suter *
* andreas.suter@psi.ch *
* *
* This program is free software; you can redistribute it and/or modify *

View File

@ -8,7 +8,7 @@
***************************************************************************/
/***************************************************************************
* Copyright (C) 2007-2014 by Andreas Suter *
* Copyright (C) 2007-2016 by Andreas Suter *
* andreas.suter@psi.ch *
* *
* This program is free software; you can redistribute it and/or modify *

View File

@ -8,7 +8,7 @@
***************************************************************************/
/***************************************************************************
* Copyright (C) 2007-2014 by Andreas Suter *
* Copyright (C) 2007-2016 by Andreas Suter *
* andreas.suter@psi.ch *
* *
* This program is free software; you can redistribute it and/or modify *

View File

@ -8,7 +8,7 @@
***************************************************************************/
/***************************************************************************
* Copyright (C) 2007-2015 by Andreas Suter *
* Copyright (C) 2007-2016 by Andreas Suter *
* andreas.suter@psi.ch *
* *
* This program is free software; you can redistribute it and/or modify *

View File

@ -8,7 +8,7 @@
***************************************************************************/
/***************************************************************************
* Copyright (C) 2007-2014 by Andreas Suter *
* Copyright (C) 2007-2016 by Andreas Suter *
* andreas.suter@psi.ch *
* *
* This program is free software; you can redistribute it and/or modify *

View File

@ -0,0 +1,81 @@
/***************************************************************************
PRunAsymmetryRRF.h
Author: Andreas Suter
e-mail: andreas.suter@psi.ch
***************************************************************************/
/***************************************************************************
* Copyright (C) 2007-2016 by Andreas Suter *
* andreas.suter@psi.ch *
* *
* This program is free software; you can redistribute it and/or modify *
* it under the terms of the GNU General Public License as published by *
* the Free Software Foundation; either version 2 of the License, or *
* (at your option) any later version. *
* *
* This program is distributed in the hope that it will be useful, *
* but WITHOUT ANY WARRANTY; without even the implied warranty of *
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
* GNU General Public License for more details. *
* *
* You should have received a copy of the GNU General Public License *
* along with this program; if not, write to the *
* Free Software Foundation, Inc., *
* 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. *
***************************************************************************/
#ifndef _PRUNASYMMETRYRRF_H_
#define _PRUNASYMMETRYRRF_H_
#include "PRunBase.h"
//---------------------------------------------------------------------------
/**
* <p>Class handling the asymmetry fit.
*/
class PRunAsymmetryRRF : public PRunBase
{
public:
PRunAsymmetryRRF();
PRunAsymmetryRRF(PMsrHandler *msrInfo, PRunDataHandler *rawData, UInt_t runNo, EPMusrHandleTag tag);
virtual ~PRunAsymmetryRRF();
virtual Double_t CalcChiSquare(const std::vector<Double_t>& par);
virtual Double_t CalcChiSquareExpected(const std::vector<Double_t>& par);
virtual Double_t CalcMaxLikelihood(const std::vector<Double_t>& par);
virtual void CalcTheory();
virtual UInt_t GetNoOfFitBins();
virtual void SetFitRangeBin(const TString fitRange);
protected:
virtual void CalcNoOfFitBins();
virtual Bool_t PrepareData();
virtual Bool_t PrepareFitData();
virtual Bool_t PrepareViewData(PRawRunData* runData, UInt_t histoNo[2]);
private:
UInt_t fAlphaBetaTag; ///< \f$ 1 \to \alpha = \beta = 1\f$; \f$ 2 \to \alpha \neq 1, \beta = 1\f$; \f$ 3 \to \alpha = 1, \beta \neq 1\f$; \f$ 4 \to \alpha \neq 1, \beta \neq 1\f$.
UInt_t fNoOfFitBins; ///< number of bins to be be fitted
Int_t fRRFPacking; ///< RRF packing for this particular run. Given in the GLOBAL-block.
PDoubleVector fForward; ///< forward histo data
PDoubleVector fForwardErr; ///< forward histo errors
PDoubleVector fBackward; ///< backward histo data
PDoubleVector fBackwardErr; ///< backward histo errors
Int_t fGoodBins[4]; ///< keep first/last good bins. 0=fgb, 1=lgb (forward); 2=fgb, 3=lgb (backward)
Bool_t SubtractFixBkg();
Bool_t SubtractEstimatedBkg();
virtual Bool_t GetProperT0(PRawRunData* runData, PMsrGlobalBlock *globalBlock, PUIntVector &forwardHisto, PUIntVector &backwardHistoNo);
virtual Bool_t GetProperDataRange(PRawRunData* runData, UInt_t histoNo[2]);
virtual void GetProperFitRange(PMsrGlobalBlock *globalBlock);
};
#endif // _PRUNASYMMETRYRRF_H_

View File

@ -8,7 +8,7 @@
***************************************************************************/
/***************************************************************************
* Copyright (C) 2007-2014 by Andreas Suter *
* Copyright (C) 2007-2016 by Andreas Suter *
* andreas.suter@psi.ch *
* *
* This program is free software; you can redistribute it and/or modify *

View File

@ -8,7 +8,7 @@
***************************************************************************/
/***************************************************************************
* Copyright (C) 2007-2015 by Andreas Suter *
* Copyright (C) 2007-2016 by Andreas Suter *
* andreas.suter@psi.ch *
* *
* This program is free software; you can redistribute it and/or modify *

View File

@ -8,7 +8,7 @@
***************************************************************************/
/***************************************************************************
* Copyright (C) 2007-2014 by Andreas Suter *
* Copyright (C) 2007-2016 by Andreas Suter *
* andreas.suter@psi.ch *
* *
* This program is free software; you can redistribute it and/or modify *
@ -37,7 +37,9 @@ using namespace std;
#include "PMsrHandler.h"
#include "PRunDataHandler.h"
#include "PRunSingleHisto.h"
#include "PRunSingleHistoRRF.h"
#include "PRunAsymmetry.h"
#include "PRunAsymmetryRRF.h"
#include "PRunMuMinus.h"
#include "PRunNonMusr.h"
@ -58,7 +60,9 @@ class PRunListCollection
virtual void SetFitRange(const TString fitRange);
virtual Double_t GetSingleHistoChisq(const std::vector<Double_t>& par) const;
virtual Double_t GetSingleHistoRRFChisq(const std::vector<Double_t>& par) const;
virtual Double_t GetAsymmetryChisq(const std::vector<Double_t>& par) const;
virtual Double_t GetAsymmetryRRFChisq(const std::vector<Double_t>& par) const;
virtual Double_t GetMuMinusChisq(const std::vector<Double_t>& par) const;
virtual Double_t GetNonMusrChisq(const std::vector<Double_t>& par) const;
@ -66,7 +70,9 @@ class PRunListCollection
virtual Double_t GetSingleRunChisq(const std::vector<Double_t>& par, const UInt_t idx) const;
virtual Double_t GetSingleHistoMaximumLikelihood(const std::vector<Double_t>& par) const;
virtual Double_t GetSingleHistoRRFMaximumLikelihood(const std::vector<Double_t>& par) const;
virtual Double_t GetAsymmetryMaximumLikelihood(const std::vector<Double_t>& par) const;
virtual Double_t GetAsymmetryRRFMaximumLikelihood(const std::vector<Double_t>& par) const;
virtual Double_t GetMuMinusMaximumLikelihood(const std::vector<Double_t>& par) const;
virtual Double_t GetNonMusrMaximumLikelihood(const std::vector<Double_t>& par) const;
@ -74,12 +80,16 @@ class PRunListCollection
virtual UInt_t GetTotalNoOfBinsFitted() const;
virtual UInt_t GetNoOfSingleHisto() const { return fRunSingleHistoList.size(); } ///< returns the number of single histogram data sets present in the msr-file
virtual UInt_t GetNoOfSingleHistoRRF() const { return fRunSingleHistoRRFList.size(); } ///< returns the number of single histogram RRF data sets present in the msr-file
virtual UInt_t GetNoOfAsymmetry() const { return fRunAsymmetryList.size(); } ///< returns the number of asymmetry data sets present in the msr-file
virtual UInt_t GetNoOfAsymmetryRRF() const { return fRunAsymmetryRRFList.size(); } ///< returns the number of asymmetry RRF data sets present in the msr-file
virtual UInt_t GetNoOfMuMinus() const { return fRunMuMinusList.size(); } ///< returns the number of mu minus data sets present in the msr-file
virtual UInt_t GetNoOfNonMusr() const { return fRunNonMusrList.size(); } ///< returns the number of non-muSR data sets present in the msr-file
virtual PRunData* GetSingleHisto(UInt_t index, EDataSwitch tag=kIndex);
virtual PRunData* GetSingleHistoRRF(UInt_t index, EDataSwitch tag=kIndex);
virtual PRunData* GetAsymmetry(UInt_t index, EDataSwitch tag=kIndex);
virtual PRunData* GetAsymmetryRRF(UInt_t index, EDataSwitch tag=kIndex);
virtual PRunData* GetMuMinus(UInt_t index, EDataSwitch tag=kIndex);
virtual PRunData* GetNonMusr(UInt_t index, EDataSwitch tag=kIndex);
@ -94,10 +104,12 @@ class PRunListCollection
PMsrHandler *fMsrInfo; ///< pointer to the msr-file handler
PRunDataHandler *fData; ///< pointer to the run-data handler
vector<PRunSingleHisto*> fRunSingleHistoList; ///< stores all processed single histogram data
vector<PRunAsymmetry*> fRunAsymmetryList; ///< stores all processed asymmetry data
vector<PRunMuMinus*> fRunMuMinusList; ///< stores all processed mu-minus data
vector<PRunNonMusr*> fRunNonMusrList; ///< stores all processed non-muSR data
vector<PRunSingleHisto*> fRunSingleHistoList; ///< stores all processed single histogram data
vector<PRunSingleHistoRRF*> fRunSingleHistoRRFList; ///< stores all processed single histogram RRF data
vector<PRunAsymmetry*> fRunAsymmetryList; ///< stores all processed asymmetry data
vector<PRunAsymmetryRRF*> fRunAsymmetryRRFList; ///< stores all processed asymmetry RRF data
vector<PRunMuMinus*> fRunMuMinusList; ///< stores all processed mu-minus data
vector<PRunNonMusr*> fRunNonMusrList; ///< stores all processed non-muSR data
};
#endif // _PRUNLISTCOLLECTION_H_

View File

@ -8,7 +8,7 @@
***************************************************************************/
/***************************************************************************
* Copyright (C) 2007-2014 by Andreas Suter *
* Copyright (C) 2007-2016 by Andreas Suter *
* andreas.suter@psi.ch *
* *
* This program is free software; you can redistribute it and/or modify *

View File

@ -8,7 +8,7 @@
***************************************************************************/
/***************************************************************************
* Copyright (C) 2007-2014 by Andreas Suter *
* Copyright (C) 2007-2016 by Andreas Suter *
* andreas.suter@psi.ch *
* *
* This program is free software; you can redistribute it and/or modify *

View File

@ -8,7 +8,7 @@
***************************************************************************/
/***************************************************************************
* Copyright (C) 2007-2014 by Andreas Suter *
* Copyright (C) 2007-2016 by Andreas Suter *
* andreas.suter@psi.ch *
* *
* This program is free software; you can redistribute it and/or modify *

View File

@ -0,0 +1,83 @@
/***************************************************************************
PRunSingleHistoRRF.h
Author: Andreas Suter
e-mail: andreas.suter@psi.ch
***************************************************************************/
/***************************************************************************
* Copyright (C) 2007-2016 by Andreas Suter *
* andreas.suter@psi.ch *
* *
* This program is free software; you can redistribute it and/or modify *
* it under the terms of the GNU General Public License as published by *
* the Free Software Foundation; either version 2 of the License, or *
* (at your option) any later version. *
* *
* This program is distributed in the hope that it will be useful, *
* but WITHOUT ANY WARRANTY; without even the implied warranty of *
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
* GNU General Public License for more details. *
* *
* You should have received a copy of the GNU General Public License *
* along with this program; if not, write to the *
* Free Software Foundation, Inc., *
* 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. *
***************************************************************************/
#ifndef _PRUNSINGLEHISTORRF_H_
#define _PRUNSINGLEHISTORRF_H_
#include "PRunBase.h"
/**
* <p>Class handling single histogram fit type.
*/
class PRunSingleHistoRRF : public PRunBase
{
public:
PRunSingleHistoRRF();
PRunSingleHistoRRF(PMsrHandler *msrInfo, PRunDataHandler *rawData, UInt_t runNo, EPMusrHandleTag tag);
virtual ~PRunSingleHistoRRF();
virtual Double_t CalcChiSquare(const std::vector<Double_t>& par);
virtual Double_t CalcChiSquareExpected(const std::vector<Double_t>& par);
virtual Double_t CalcMaxLikelihood(const std::vector<Double_t>& par);
virtual void CalcTheory();
virtual UInt_t GetNoOfFitBins();
virtual void SetFitRangeBin(const TString fitRange);
protected:
virtual void CalcNoOfFitBins();
virtual Bool_t PrepareData();
virtual Bool_t PrepareFitData(PRawRunData* runData, const UInt_t histoNo);
virtual Bool_t PrepareViewData(PRawRunData* runData, const UInt_t histoNo);
private:
static const Double_t fN0EstimateEndTime = 1.0; ///< end time in (us) over which N0 is estimated. Should eventually be estimated automatically ...
UInt_t fNoOfFitBins; ///< number of bins to be fitted
Double_t fBackground; ///< needed if background range is given (units: 1/bin)
Int_t fRRFPacking; ///< RRF packing for this particular run. Given in the GLOBAL-block.
Int_t fGoodBins[2]; ///< keep first/last good bins. 0=fgb, 1=lgb
PDoubleVector fForward; ///< forward histo data
PDoubleVector fM; ///< vector holding M(t) = [N(t)-N_bkg] exp(+t/tau). Needed to estimate N0.
PDoubleVector fMerr; ///< vector holding the error of M(t): M_err = exp(+t/tau) sqrt(N(t)).
PDoubleVector fW; ///< vector holding the weight needed to estimate N0, and errN0.
PDoubleVector fAerr; ///< vector holding the errors of estimated A(t)
virtual Bool_t GetProperT0(PRawRunData* runData, PMsrGlobalBlock *globalBlock, PUIntVector &histoNo);
virtual Bool_t GetProperDataRange();
virtual void GetProperFitRange(PMsrGlobalBlock *globalBlock);
virtual Double_t GetMainFrequency(PDoubleVector &data);
virtual Double_t EstimateN0(Double_t &errN0, Double_t freqMax);
virtual Bool_t EstimateBkg(UInt_t histoNo);
};
#endif // _PRUNSINGLEHISTORRF_H_

View File

@ -8,7 +8,7 @@
***************************************************************************/
/***************************************************************************
* Copyright (C) 2007-2014 by Andreas Suter *
* Copyright (C) 2007-2016 by Andreas Suter *
* andreas.suter@psi.ch *
* *
* This program is free software; you can redistribute it and/or modify *

View File

@ -8,7 +8,7 @@
***************************************************************************/
/***************************************************************************
* Copyright (C) 2007-2014 by Andreas Suter *
* Copyright (C) 2007-2016 by Andreas Suter *
* andreas.suter@psi.ch *
* *
* This program is free software; you can redistribute it and/or modify *

View File

@ -8,7 +8,7 @@
***************************************************************************/
/***************************************************************************
* Copyright (C) 2007-2014 by Andreas Suter *
* Copyright (C) 2007-2016 by Andreas Suter *
* andreas.suter@psi.ch *
* *
* This program is free software; you can redistribute it and/or modify *

View File

@ -8,7 +8,7 @@
***************************************************************************/
/***************************************************************************
* Copyright (C) 2007-2014 by Andreas Suter *
* Copyright (C) 2007-2016 by Andreas Suter *
* andreas.suter@psi.ch *
* *
* This program is free software; you can redistribute it and/or modify *

View File

@ -8,7 +8,7 @@
***************************************************************************/
/***************************************************************************
* Copyright (C) 2007-2014 by Andreas Suter *
* Copyright (C) 2007-2016 by Andreas Suter *
* andreas.suter@psi.ch *
* *
* This program is free software; you can redistribute it and/or modify *

View File

@ -8,7 +8,7 @@
***************************************************************************/
/***************************************************************************
* Copyright (C) 2007-2014 by Andreas Suter *
* Copyright (C) 2007-2016 by Andreas Suter *
* andreas.suter@psi.ch *
* *
* This program is free software; you can redistribute it and/or modify *

View File

@ -175,10 +175,10 @@ void musrfit_dump_ascii(char *fileName, PRunListCollection *runList)
// go through the run list, get the data and dump it in a file
int runCounter = 0;
PRunData *data;
// single histos
unsigned int size = runList->GetNoOfSingleHisto();
PRunData *data;
if (size > 0) {
for (unsigned int i=0; i<size; i++) {
data = runList->GetSingleHisto(i);
@ -190,6 +190,19 @@ void musrfit_dump_ascii(char *fileName, PRunListCollection *runList)
}
}
// single histos
size = runList->GetNoOfSingleHistoRRF();
if (size > 0) {
for (unsigned int i=0; i<size; i++) {
data = runList->GetSingleHistoRRF(i);
if (data) {
// dump data
musrfit_write_ascii(fln, data, runCounter);
runCounter++;
}
}
}
// asymmetry
size = runList->GetNoOfAsymmetry();
if (size > 0) {
@ -203,6 +216,19 @@ void musrfit_dump_ascii(char *fileName, PRunListCollection *runList)
}
}
// asymmetry RRF
size = runList->GetNoOfAsymmetryRRF();
if (size > 0) {
for (unsigned int i=0; i<size; i++) {
data = runList->GetAsymmetryRRF(i);
if (data) {
// dump data
musrfit_write_ascii(fln, data, runCounter);
runCounter++;
}
}
}
// muMinus
size = runList->GetNoOfMuMinus();
if (size > 0) {
@ -308,10 +334,10 @@ void musrfit_dump_root(char *fileName, PRunListCollection *runList)
// go through the run list, get the data and dump it in a file
int runCounter = 0;
PRunData *data;
// single histos
unsigned int size = runList->GetNoOfSingleHisto();
PRunData *data;
if (size > 0) {
for (unsigned int i=0; i<size; i++) {
data = runList->GetSingleHisto(i);
@ -323,6 +349,19 @@ void musrfit_dump_root(char *fileName, PRunListCollection *runList)
}
}
// single histo RRF
size = runList->GetNoOfSingleHistoRRF();
if (size > 0) {
for (unsigned int i=0; i<size; i++) {
data = runList->GetSingleHistoRRF(i);
if (data) {
// dump data
musrfit_write_root(f, fln, data, runCounter);
runCounter++;
}
}
}
// asymmetry
size = runList->GetNoOfAsymmetry();
if (size > 0) {
@ -336,6 +375,19 @@ void musrfit_dump_root(char *fileName, PRunListCollection *runList)
}
}
// asymmetry RRF
size = runList->GetNoOfAsymmetryRRF();
if (size > 0) {
for (unsigned int i=0; i<size; i++) {
data = runList->GetAsymmetryRRF(i);
if (data) {
// dump data
musrfit_write_root(f, fln, data, runCounter);
runCounter++;
}
}
}
// muMinus
size = runList->GetNoOfMuMinus();
if (size > 0) {

View File

@ -669,6 +669,7 @@ Int_t main(Int_t argc, Char_t *argv[])
}
switch (fitType) {
case MSR_FITTYPE_SINGLE_HISTO:
case MSR_FITTYPE_SINGLE_HISTO_RRF:
case MSR_FITTYPE_MU_MINUS:
if ((runList->at(i).GetRunNameSize() == 1) && (runList->at(i).GetForwardHistoNoSize() == 1)) { // no addruns / no grouping
// feed necessary data
@ -796,6 +797,7 @@ Int_t main(Int_t argc, Char_t *argv[])
}
break;
case MSR_FITTYPE_ASYM:
case MSR_FITTYPE_ASYM_RRF:
if ((runList->at(i).GetRunNameSize() == 1) && (runList->at(i).GetForwardHistoNoSize() == 1)) { // no addruns / no grouping
// feed necessary data forward
musrT0Data.InitData();

View File

@ -179,14 +179,20 @@ void analyticFakeData(const TString type, UInt_t runNo)
phase.push_back((5.0 + 2.0*rand.Rndm())*TMath::Pi()/180.0 + TMath::TwoPi()/noOfHistos * (Double_t)i);
const Double_t gamma = 0.0000135538817; // gamma/(2pi)
Double_t bb0 = 90000.0; // field in Gauss
Double_t rate0 = 1.0/1000.0; // in 1/ns
Double_t bb0 = 5000.0; // field in Gauss
Double_t rate0 = 7.0/1000.0; // in 1/ns
Double_t frac0 = 0.5;
Double_t bb1 = bb0 + 200.0; // field in Gauss
Double_t rate1 = 0.75/1000.0; // in 1/ns
Double_t frac1 = 0.2;
Double_t bb2 = bb0 + 600.0; // field in Gauss
Double_t rate2 = 0.25/1000.0; // in 1/ns
// fake function parameters header info: only for test purposes
cout << ">> write fake header for TMusrRoot" << endl;
if (type.CompareTo("TLemRunHeader")) {
TDoubleVector dvec;
header->Set("FakeFct/Def", "N0 exp(-t/tau_mu) [1 + A exp(-1/2 (t sigma)^2) cos(gamma_mu B t + phi)] + bkg");
header->Set("FakeFct/Def", "N0 exp(-t/tau_mu) [1 + sum_{k=0}^2 frac_k A_0 exp(-1/2 (t sigma_k)^2) cos(gamma_mu B_k t + phi)] + bkg");
for (UInt_t i=0; i<noOfHistos; i++)
dvec.push_back(n0[i]);
header->Set("FakeFct/N0", dvec);
@ -202,10 +208,20 @@ void analyticFakeData(const TString type, UInt_t runNo)
for (UInt_t i=0; i<noOfHistos; i++)
dvec.push_back(phase[i]);
header->Set("FakeFct/phase", dvec);
prop.Set("B", bb0, "G");
header->Set("FakeFct/B", prop);
prop.Set("lambda", rate0, "1/usec");
header->Set("FakeFct/lambda", prop);
prop.Set("B0", bb0, "G");
header->Set("FakeFct/B0", prop);
prop.Set("rate0", 1.0e3*rate0, "1/usec");
header->Set("FakeFct/rate0", prop);
header->Set("FakeFct/frac0", frac0);
prop.Set("B1", bb1, "G");
header->Set("FakeFct/B1", prop);
prop.Set("rate1", 1.0e3*rate1, "1/usec");
header->Set("FakeFct/rate1", prop);
header->Set("FakeFct/frac1", frac1);
prop.Set("B2", bb2, "G");
header->Set("FakeFct/B2", prop);
prop.Set("rate2", 1.0e3*rate2, "1/usec");
header->Set("FakeFct/rate2", prop);
}
cout << ">> create histo objects" << endl;
@ -228,7 +244,10 @@ void analyticFakeData(const TString type, UInt_t runNo)
histo[i]->SetBinContent(j+1, bkg[i]);
} else {
time = (Double_t)(j-t0[i])*timeResolution;
dval = (Double_t)n0[i]*TMath::Exp(-time/tau)*(1.0+a0[i]*TMath::Exp(-0.5*TMath::Power(time*rate0,2))*TMath::Cos(TMath::TwoPi()*gamma*bb0*time+phase[i]))+(Double_t)bkg[i];
dval = (Double_t)n0[i]*TMath::Exp(-time/tau)*(1.0+
frac0*a0[i]*TMath::Exp(-0.5*TMath::Power(time*rate0,2))*TMath::Cos(TMath::TwoPi()*gamma*bb0*time+phase[i]) +
frac1*a0[i]*TMath::Exp(-0.5*TMath::Power(time*rate1,2))*TMath::Cos(TMath::TwoPi()*gamma*bb1*time+phase[i]) +
(1.0-frac0-frac1)*a0[i]*TMath::Exp(-0.5*TMath::Power(time*rate2,2))*TMath::Cos(TMath::TwoPi()*gamma*bb2*time+phase[i]))+(Double_t)bkg[i];
histo[i]->SetBinContent(j+1, dval);
}
}