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

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> 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> 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) //-------------------------------------------------------------------------- diff --git a/src/classes/PMusr.cpp b/src/classes/PMusr.cpp index 478db109..a84f4765 100644 --- a/src/classes/PMusr.cpp +++ b/src/classes/PMusr.cpp @@ -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 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) +//-------------------------------------------------------------------------- +/** + *

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) +//-------------------------------------------------------------------------- +/** + *

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) +//-------------------------------------------------------------------------- +/** + *

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) +//-------------------------------------------------------------------------- +/** + *

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) //-------------------------------------------------------------------------- diff --git a/src/classes/PMusrCanvas.cpp b/src/classes/PMusrCanvas.cpp index e34486bf..ec99e93a 100644 --- a/src/classes/PMusrCanvas.cpp +++ b/src/classes/PMusrCanvas.cpp @@ -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(); diff --git a/src/classes/PMusrT0.cpp b/src/classes/PMusrT0.cpp index dfedb61d..c13414fc 100644 --- a/src/classes/PMusrT0.cpp +++ b/src/classes/PMusrT0.cpp @@ -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 * diff --git a/src/classes/PPrepFourier.cpp b/src/classes/PPrepFourier.cpp index 1b542c5b..a28bc29f 100644 --- a/src/classes/PPrepFourier.cpp +++ b/src/classes/PPrepFourier.cpp @@ -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[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; } diff --git a/src/classes/PRunAsymmetry.cpp b/src/classes/PRunAsymmetry.cpp index c04d2a73..a2572b32 100644 --- a/src/classes/PRunAsymmetry.cpp +++ b/src/classes/PRunAsymmetry.cpp @@ -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; } diff --git a/src/classes/PRunAsymmetryRRF.cpp b/src/classes/PRunAsymmetryRRF.cpp new file mode 100644 index 00000000..7e7f21df --- /dev/null +++ b/src/classes/PRunAsymmetryRRF.cpp @@ -0,0 +1,1511 @@ +/*************************************************************************** + + PRunAsymmetryRRF.cpp + + 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. * + ***************************************************************************/ + +#ifdef HAVE_CONFIG_H +#include "config.h" +#endif + +#ifdef HAVE_GOMP +#include +#endif + +#include + +#include +#include +using namespace std; + +#include +#include +#include + +#include "PMusr.h" +#include "PRunAsymmetryRRF.h" + +//-------------------------------------------------------------------------- +// Constructor +//-------------------------------------------------------------------------- +/** + *

Constructor + */ +PRunAsymmetryRRF::PRunAsymmetryRRF() : PRunBase() +{ + fNoOfFitBins = 0; + fRRFPacking = -1; + + // the 2 following variables are need in case fit range is given in bins, and since + // the fit range can be changed in the command block, these variables need to be accessible + fGoodBins[0] = -1; + fGoodBins[1] = -1; +} + +//-------------------------------------------------------------------------- +// Constructor +//-------------------------------------------------------------------------- +/** + *

Constructor + * + * \param msrInfo pointer to the msr-file handler + * \param rawData raw run data + * \param runNo number of the run within the msr-file + * \param tag tag showing what shall be done: kFit == fitting, kView == viewing + */ +PRunAsymmetryRRF::PRunAsymmetryRRF(PMsrHandler *msrInfo, PRunDataHandler *rawData, UInt_t runNo, EPMusrHandleTag tag) : PRunBase(msrInfo, rawData, runNo, tag) +{ + // the 2 following variables are need in case fit range is given in bins, and since + // the fit range can be changed in the command block, these variables need to be accessible + fGoodBins[0] = -1; + fGoodBins[1] = -1; + + fRRFPacking = fMsrInfo->GetMsrGlobal()->GetRRFPacking(); + if (fRRFPacking == -1) { // this should NOT happen, somethin is severely wrong + cerr << endl << ">> PRunAsymmetryRRF::PRunAsymmetryRRF(): **SEVERE ERROR**: Couldn't find any RRF packing information!"; + cerr << endl << ">> This is very bad :-(, will quit ..."; + cerr << endl; + fValid = false; + return; + } + + // check if alpha and/or beta is fixed -------------------- + + PMsrParamList *param = msrInfo->GetMsrParamList(); + + // check if alpha is given + if (fRunInfo->GetAlphaParamNo() == -1) { // no alpha given + cerr << endl << ">> PRunAsymmetryRRF::PRunAsymmetryRRF(): **ERROR** no alpha parameter given! This is needed for an asymmetry fit!"; + cerr << endl; + fValid = false; + return; + } + // check if alpha parameter is within proper bounds + if ((fRunInfo->GetAlphaParamNo() < 0) || (fRunInfo->GetAlphaParamNo() > (Int_t)param->size())) { + cerr << endl << ">> PRunAsymmetryRRF::PRunAsymmetryRRF(): **ERROR** alpha parameter no = " << fRunInfo->GetAlphaParamNo(); + cerr << endl << ">> This is out of bound, since there are only " << param->size() << " parameters."; + cerr << endl; + fValid = false; + return; + } + // check if alpha is fixed + Bool_t alphaFixedToOne = false; + if (((*param)[fRunInfo->GetAlphaParamNo()-1].fStep == 0.0) && + ((*param)[fRunInfo->GetAlphaParamNo()-1].fValue == 1.0)) + alphaFixedToOne = true; + + // check if beta is given + Bool_t betaFixedToOne = false; + if (fRunInfo->GetBetaParamNo() == -1) { // no beta given hence assuming beta == 1 + betaFixedToOne = true; + } else if ((fRunInfo->GetBetaParamNo() < 0) || (fRunInfo->GetBetaParamNo() > (Int_t)param->size())) { // check if beta parameter is within proper bounds + cerr << endl << ">> PRunAsymmetryRRF::PRunAsymmetryRRF(): **ERROR** beta parameter no = " << fRunInfo->GetBetaParamNo(); + cerr << endl << ">> This is out of bound, since there are only " << param->size() << " parameters."; + cerr << endl; + fValid = false; + return; + } else { // check if beta is fixed + if (((*param)[fRunInfo->GetBetaParamNo()-1].fStep == 0.0) && + ((*param)[fRunInfo->GetBetaParamNo()-1].fValue == 1.0)) + betaFixedToOne = true; + } + + // set fAlphaBetaTag + if (alphaFixedToOne && betaFixedToOne) // alpha == 1, beta == 1 + fAlphaBetaTag = 1; + else if (!alphaFixedToOne && betaFixedToOne) // alpha != 1, beta == 1 + fAlphaBetaTag = 2; + else if (alphaFixedToOne && !betaFixedToOne) // alpha == 1, beta != 1 + fAlphaBetaTag = 3; + else + fAlphaBetaTag = 4; + + // calculate fData + if (!PrepareData()) + fValid = false; +} + +//-------------------------------------------------------------------------- +// Destructor +//-------------------------------------------------------------------------- +/** + *

Destructor. + */ +PRunAsymmetryRRF::~PRunAsymmetryRRF() +{ + fForward.clear(); + fForwardErr.clear(); + fBackward.clear(); + fBackwardErr.clear(); +} + +//-------------------------------------------------------------------------- +// CalcChiSquare (public) +//-------------------------------------------------------------------------- +/** + *

Calculate chi-square. + * + * return: + * - chisq value + * + * \param par parameter vector iterated by minuit2 + */ +Double_t PRunAsymmetryRRF::CalcChiSquare(const std::vector& par) +{ + Double_t chisq = 0.0; + Double_t diff = 0.0; + Double_t asymFcnValue = 0.0; + Double_t a, b, f; + + // calculate functions + for (Int_t i=0; iGetNoOfFuncs(); i++) { + fFuncValues[i] = fMsrInfo->EvalFunc(fMsrInfo->GetFuncNo(i), *fRunInfo->GetMap(), par); + } + + // calculate chi square + Double_t time(1.0); + Int_t i, N(static_cast(fData.GetValue()->size())); + + // In order not to have an IF in the next loop, determine the start and end bins for the fit range now + Int_t startTimeBin = static_cast(ceil((fFitStartTime - fData.GetDataTimeStart())/fData.GetDataTimeStep())); + if (startTimeBin < 0) + startTimeBin = 0; + Int_t endTimeBin = static_cast(floor((fFitEndTime - fData.GetDataTimeStart())/fData.GetDataTimeStep())) + 1; + if (endTimeBin > N) + endTimeBin = N; + + // Calculate the theory function once to ensure one function evaluation for the current set of parameters. + // This is needed for the LF and user functions where some non-thread-save calculations only need to be calculated once + // for a given set of parameters---which should be done outside of the parallelized loop. + // For all other functions it means a tiny and acceptable overhead. + asymFcnValue = fTheory->Func(time, par, fFuncValues); + + #ifdef HAVE_GOMP + Int_t chunk = (endTimeBin - startTimeBin)/omp_get_num_procs(); + if (chunk < 10) + chunk = 10; + #pragma omp parallel for default(shared) private(i,time,diff,asymFcnValue,a,b,f) schedule(dynamic,chunk) reduction(+:chisq) + #endif + for (i=startTimeBin; iFunc(time, par, fFuncValues); + break; + case 2: // alpha != 1, beta == 1 + a = par[fRunInfo->GetAlphaParamNo()-1]; + f = fTheory->Func(time, par, fFuncValues); + asymFcnValue = (f*(a+1.0)-(a-1.0))/((a+1.0)-f*(a-1.0)); + break; + case 3: // alpha == 1, beta != 1 + b = par[fRunInfo->GetBetaParamNo()-1]; + f = fTheory->Func(time, par, fFuncValues); + asymFcnValue = f*(b+1.0)/(2.0-f*(b-1.0)); + break; + case 4: // alpha != 1, beta != 1 + a = par[fRunInfo->GetAlphaParamNo()-1]; + b = par[fRunInfo->GetBetaParamNo()-1]; + f = fTheory->Func(time, par, fFuncValues); + asymFcnValue = (f*(a*b+1.0)-(a-1.0))/((a+1.0)-f*(a*b-1.0)); + break; + default: + asymFcnValue = 0.0; + break; + } + diff = fData.GetValue()->at(i) - asymFcnValue; + chisq += diff*diff / (fData.GetError()->at(i)*fData.GetError()->at(i)); + } + + return chisq; +} + +//-------------------------------------------------------------------------- +// CalcChiSquareExpected (public) +//-------------------------------------------------------------------------- +/** + *

Calculate expected chi-square. Currently not implemented since not clear what to be done. + * + * return: + * - chisq value == 0.0 + * + * \param par parameter vector iterated by minuit2 + */ +Double_t PRunAsymmetryRRF::CalcChiSquareExpected(const std::vector& par) +{ + return 0.0; +} + +//-------------------------------------------------------------------------- +// CalcMaxLikelihood (public) +//-------------------------------------------------------------------------- +/** + *

NOT IMPLEMENTED!! + * + * \param par parameter vector iterated by minuit2 + */ +Double_t PRunAsymmetryRRF::CalcMaxLikelihood(const std::vector& par) +{ + cout << endl << "PRunAsymmetryRRF::CalcMaxLikelihood(): not implemented yet ..." << endl; + + return 1.0; +} + +//-------------------------------------------------------------------------- +// GetNoOfFitBins (public) +//-------------------------------------------------------------------------- +/** + *

Calculate the number of fitted bins for the current fit range. + * + * return: number of fitted bins. + */ +UInt_t PRunAsymmetryRRF::GetNoOfFitBins() +{ + CalcNoOfFitBins(); + + return fNoOfFitBins; +} + +//-------------------------------------------------------------------------- +// SetFitRangeBin (public) +//-------------------------------------------------------------------------- +/** + *

Allows to change the fit range on the fly. Used in the COMMAND block. + * The syntax of the string is: FIT_RANGE fgb[+n00] lgb[-n01] [fgb[+n10] lgb[-n11] ... fgb[+nN0] lgb[-nN1]]. + * If only one pair of fgb/lgb is given, it is used for all runs in the RUN block section. + * If multiple fgb/lgb's are given, the number N has to be the number of RUN blocks in + * the msr-file. + * + *

nXY are offsets which can be used to shift, limit the fit range. + * + * \param fitRange string containing the necessary information. + */ +void PRunAsymmetryRRF::SetFitRangeBin(const TString fitRange) +{ + TObjArray *tok = 0; + TObjString *ostr = 0; + TString str; + Ssiz_t idx = -1; + Int_t offset = 0; + + tok = fitRange.Tokenize(" \t"); + + if (tok->GetEntries() == 3) { // structure FIT_RANGE fgb+n0 lgb-n1 + // handle fgb+n0 entry + ostr = (TObjString*) tok->At(1); + str = ostr->GetString(); + // check if there is an offset present + idx = str.First("+"); + if (idx != -1) { // offset present + str.Remove(0, idx+1); + if (str.IsFloat()) // if str is a valid number, convert is to an integer + offset = str.Atoi(); + } + fFitStartTime = (fGoodBins[0] + offset - fT0s[0]) * fTimeResolution; + + // handle lgb-n1 entry + ostr = (TObjString*) tok->At(2); + str = ostr->GetString(); + // check if there is an offset present + idx = str.First("-"); + if (idx != -1) { // offset present + str.Remove(0, idx+1); + if (str.IsFloat()) // if str is a valid number, convert is to an integer + offset = str.Atoi(); + } + fFitEndTime = (fGoodBins[1] - offset - fT0s[0]) * fTimeResolution; + } else if ((tok->GetEntries() > 3) && (tok->GetEntries() % 2 == 1)) { // structure FIT_RANGE fgb[+n00] lgb[-n01] [fgb[+n10] lgb[-n11] ... fgb[+nN0] lgb[-nN1]] + 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 << ">> will ignore it. Sorry ..." << endl; + } else { + // handle fgb+n0 entry + ostr = (TObjString*) tok->At(pos); + str = ostr->GetString(); + // check if there is an offset present + idx = str.First("+"); + if (idx != -1) { // offset present + str.Remove(0, idx+1); + if (str.IsFloat()) // if str is a valid number, convert is to an integer + offset = str.Atoi(); + } + fFitStartTime = (fGoodBins[0] + offset - fT0s[0]) * fTimeResolution; + + // handle lgb-n1 entry + ostr = (TObjString*) tok->At(pos+1); + str = ostr->GetString(); + // check if there is an offset present + idx = str.First("-"); + if (idx != -1) { // offset present + str.Remove(0, idx+1); + if (str.IsFloat()) // if str is a valid number, convert is to an integer + offset = str.Atoi(); + } + fFitEndTime = (fGoodBins[1] - offset - fT0s[0]) * fTimeResolution; + } + } else { // error + cerr << endl << ">> PRunSingleHisto::SetFitRangeBin(): **ERROR** invalid FIT_RANGE command found: '" << fitRange << "'"; + cerr << endl << ">> will ignore it. Sorry ..." << endl; + } + + // clean up + if (tok) { + delete tok; + } +} + +//-------------------------------------------------------------------------- +// CalcNoOfFitBins (protected) +//-------------------------------------------------------------------------- +/** + *

Calculate the number of fitted bins for the current fit range. + */ +void PRunAsymmetryRRF::CalcNoOfFitBins() +{ + // In order not having to loop over all bins and to stay consistent with the chisq method, calculate the start and end bins explicitly + Int_t startTimeBin = static_cast(ceil((fFitStartTime - fData.GetDataTimeStart())/fData.GetDataTimeStep())); + if (startTimeBin < 0) + startTimeBin = 0; + Int_t endTimeBin = static_cast(floor((fFitEndTime - fData.GetDataTimeStart())/fData.GetDataTimeStep())) + 1; + if (endTimeBin > static_cast(fData.GetValue()->size())) + endTimeBin = fData.GetValue()->size(); + + if (endTimeBin > startTimeBin) + fNoOfFitBins = endTimeBin - startTimeBin; + else + fNoOfFitBins = 0; +} + +//-------------------------------------------------------------------------- +// CalcTheory (protected) +//-------------------------------------------------------------------------- +/** + *

Calculate theory for a given set of fit-parameters. + */ +void PRunAsymmetryRRF::CalcTheory() +{ + // feed the parameter vector + std::vector par; + PMsrParamList *paramList = fMsrInfo->GetMsrParamList(); + for (UInt_t i=0; isize(); i++) + par.push_back((*paramList)[i].fValue); + + // calculate functions + for (Int_t i=0; iGetNoOfFuncs(); i++) { + fFuncValues[i] = fMsrInfo->EvalFunc(fMsrInfo->GetFuncNo(i), *fRunInfo->GetMap(), par); + } + + // calculate asymmetry + Double_t asymFcnValue = 0.0; + Double_t a, b, f; + Double_t time; + for (UInt_t i=0; isize(); i++) { + time = fData.GetDataTimeStart() + (Double_t)i*fData.GetDataTimeStep(); + switch (fAlphaBetaTag) { + case 1: // alpha == 1, beta == 1 + asymFcnValue = fTheory->Func(time, par, fFuncValues); + break; + case 2: // alpha != 1, beta == 1 + a = par[fRunInfo->GetAlphaParamNo()-1]; + f = fTheory->Func(time, par, fFuncValues); + asymFcnValue = (f*(a+1.0)-(a-1.0))/((a+1.0)-f*(a-1.0)); + break; + case 3: // alpha == 1, beta != 1 + b = par[fRunInfo->GetBetaParamNo()-1]; + f = fTheory->Func(time, par, fFuncValues); + asymFcnValue = f*(b+1.0)/(2.0-f*(b-1.0)); + break; + case 4: // alpha != 1, beta != 1 + a = par[fRunInfo->GetAlphaParamNo()-1]; + b = par[fRunInfo->GetBetaParamNo()-1]; + f = fTheory->Func(time, par, fFuncValues); + asymFcnValue = (f*(a*b+1.0)-(a-1.0))/((a+1.0)-f*(a*b-1.0)); + break; + default: + asymFcnValue = 0.0; + break; + } + fData.AppendTheoryValue(asymFcnValue); + } + + // clean up + par.clear(); +} + +//-------------------------------------------------------------------------- +// PrepareData (protected) +//-------------------------------------------------------------------------- +/** + *

Prepare data for fitting or viewing. What is already processed at this stage: + * - get all needed forward/backward histograms + * - get time resolution + * - get start/stop fit time + * - get t0's and perform necessary cross checks (e.g. if t0 of msr-file (if present) are consistent with t0 of the data files, etc.) + * - add runs (if addruns are present) + * - group histograms (if grouping is present) + * - subtract background + * + * Error propagation for \f$ A_i = (f_i^{\rm c}-b_i^{\rm c})/(f_i^{\rm c}+b_i^{\rm c})\f$: + * \f[ \Delta A_i = \pm\frac{2}{(f_i^{\rm c}+b_i^{\rm c})^2}\left[ + * (b_i^{\rm c})^2 (\Delta f_i^{\rm c})^2 + + * (\Delta b_i^{\rm c})^2 (f_i^{\rm c})^2\right]^{1/2}\f] + * + * return: + * - true if everthing went smooth + * - false, otherwise. + */ +Bool_t PRunAsymmetryRRF::PrepareData() +{ + // keep the Global block info + PMsrGlobalBlock *globalBlock = fMsrInfo->GetMsrGlobal(); + + // get forward/backward histo from PRunDataHandler object ------------------------ + // get the correct run + PRawRunData *runData = fRawData->GetRunData(*(fRunInfo->GetRunName())); + if (!runData) { // run not found + cerr << endl << ">> PRunAsymmetryRRF::PrepareData(): **ERROR** Couldn't get run " << fRunInfo->GetRunName()->Data() << "!"; + cerr << endl; + return false; + } + + // collect histogram numbers + PUIntVector forwardHistoNo; + PUIntVector backwardHistoNo; + for (UInt_t i=0; iGetForwardHistoNoSize(); i++) { + forwardHistoNo.push_back(fRunInfo->GetForwardHistoNo(i)); + + if (!runData->IsPresent(forwardHistoNo[i])) { + cerr << endl << ">> PRunAsymmetryRRF::PrepareData(): **PANIC ERROR**:"; + cerr << endl << ">> forwardHistoNo found = " << forwardHistoNo[i] << ", which is NOT present in the data file!?!?"; + cerr << endl << ">> Will quit :-("; + cerr << endl; + // clean up + forwardHistoNo.clear(); + backwardHistoNo.clear(); + return false; + } + } + for (UInt_t i=0; iGetBackwardHistoNoSize(); i++) { + backwardHistoNo.push_back(fRunInfo->GetBackwardHistoNo(i)); + + if (!runData->IsPresent(backwardHistoNo[i])) { + cerr << endl << ">> PRunAsymmetryRRF::PrepareData(): **PANIC ERROR**:"; + cerr << endl << ">> backwardHistoNo found = " << backwardHistoNo[i] << ", which is NOT present in the data file!?!?"; + cerr << endl << ">> Will quit :-("; + cerr << endl; + // clean up + forwardHistoNo.clear(); + backwardHistoNo.clear(); + return false; + } + } + if (forwardHistoNo.size() != backwardHistoNo.size()) { + cerr << endl << ">> PRunAsymmetryRRF::PrepareData(): **PANIC ERROR**:"; + cerr << endl << ">> # of forward histograms different from # of backward histograms."; + cerr << endl << ">> Will quit :-("; + cerr << endl; + // clean up + forwardHistoNo.clear(); + backwardHistoNo.clear(); + return false; + } + + // keep the time resolution in (us) + fTimeResolution = runData->GetTimeResolution()/1.0e3; + cout.precision(10); + cout << endl << ">> PRunAsymmetryRRF::PrepareData(): time resolution=" << fixed << runData->GetTimeResolution() << "(ns)" << endl; + + // get all the proper t0's and addt0's for the current RUN block + if (!GetProperT0(runData, globalBlock, forwardHistoNo, backwardHistoNo)) { + return false; + } + + // keep the histo of each group at this point (addruns handled below) + vector forward, backward; + forward.resize(forwardHistoNo.size()); // resize to number of groups + backward.resize(backwardHistoNo.size()); // resize to numer of groups + for (UInt_t i=0; iGetDataBin(forwardHistoNo[i])->size()); + backward[i].resize(runData->GetDataBin(backwardHistoNo[i])->size()); + forward[i] = *runData->GetDataBin(forwardHistoNo[i]); + backward[i] = *runData->GetDataBin(backwardHistoNo[i]); + } + + // check if addrun's are present, and if yes add data + // check if there are runs to be added to the current one + if (fRunInfo->GetRunNameSize() > 1) { // runs to be added present + PRawRunData *addRunData; + for (UInt_t i=1; iGetRunNameSize(); i++) { + // get run to be added to the main one + addRunData = fRawData->GetRunData(*(fRunInfo->GetRunName(i))); + if (addRunData == 0) { // couldn't get run + cerr << endl << ">> PRunAsymmetryRRF::PrepareData(): **ERROR** Couldn't get addrun " << fRunInfo->GetRunName(i)->Data() << "!"; + cerr << endl; + return false; + } + + // add forward run + UInt_t addRunSize; + for (UInt_t k=0; kGetDataBin(forwardHistoNo[k])->size(); + for (UInt_t j=0; jGetDataBin(forwardHistoNo[k])->size(); j++) { // loop over the bin indices + // make sure that the index stays in the proper range + if ((j+(Int_t)fAddT0s[i-1][2*k]-(Int_t)fT0s[2*k] >= 0) && (j+(Int_t)fAddT0s[i-1][2*k]-(Int_t)fT0s[2*k] < addRunSize)) { + forward[k][j] += addRunData->GetDataBin(forwardHistoNo[k])->at(j+(Int_t)fAddT0s[i-1][2*k]-(Int_t)fT0s[2*k]); + } + } + } + + // add backward run + for (UInt_t k=0; kGetDataBin(backwardHistoNo[k])->size(); + for (UInt_t j=0; jGetDataBin(backwardHistoNo[k])->size(); j++) { // loop over the bin indices + // make sure that the index stays in the proper range + if ((j+(Int_t)fAddT0s[i-1][2*k+1]-(Int_t)fT0s[2*k+1] >= 0) && (j+(Int_t)fAddT0s[i-1][2*k+1]-(Int_t)fT0s[2*k+1] < addRunSize)) { + backward[k][j] += addRunData->GetDataBin(backwardHistoNo[k])->at(j+(Int_t)fAddT0s[i-1][2*k+1]-(Int_t)fT0s[2*k+1]); + } + } + } + } + } + + // set forward/backward histo data of the first group + fForward.resize(forward[0].size()); + fBackward.resize(backward[0].size()); + for (UInt_t i=0; iGetDataBin(forwardHistoNo[i])->size(); j++) { // loop over the bin indices + // make sure that the index stays within proper range + if ((j+fT0s[2*i]-fT0s[0] >= 0) && (j+fT0s[2*i]-fT0s[0] < runData->GetDataBin(forwardHistoNo[i])->size())) { + fForward[j] += forward[i][j+(Int_t)fT0s[2*i]-(Int_t)fT0s[0]]; + } + } + } + + // group histograms, add all the remaining backward histograms of the group + for (UInt_t i=1; iGetDataBin(backwardHistoNo[i])->size(); j++) { // loop over the bin indices + // make sure that the index stays within proper range + if ((j+fT0s[2*i+1]-fT0s[1] >= 0) && (j+fT0s[2*i+1]-fT0s[1] < runData->GetDataBin(backwardHistoNo[i])->size())) { + fBackward[j] += backward[i][j+(Int_t)fT0s[2*i+1]-(Int_t)fT0s[1]]; + } + } + } + + // subtract background from histogramms ------------------------------------------ + if (fRunInfo->GetBkgFix(0) == PMUSR_UNDEFINED) { // no fixed background given + if (fRunInfo->GetBkgRange(0) >= 0) { // background range given + if (!SubtractEstimatedBkg()) + return false; + } else { // no background given to do the job, try to estimate it + fRunInfo->SetBkgRange(static_cast(fT0s[0]*0.1), 0); + fRunInfo->SetBkgRange(static_cast(fT0s[0]*0.6), 1); + fRunInfo->SetBkgRange(static_cast(fT0s[1]*0.1), 2); + fRunInfo->SetBkgRange(static_cast(fT0s[1]*0.6), 3); + cerr << endl << ">> PRunAsymmetryRRF::PrepareData(): **WARNING** Neither fix background nor background bins are given!"; + cerr << endl << ">> Will try the following:"; + cerr << endl << ">> forward: bkg start = " << fRunInfo->GetBkgRange(0) << ", bkg end = " << fRunInfo->GetBkgRange(1); + cerr << endl << ">> backward: bkg start = " << fRunInfo->GetBkgRange(2) << ", bkg end = " << fRunInfo->GetBkgRange(3); + cerr << endl << ">> NO WARRANTY THAT THIS MAKES ANY SENSE! Better check ..."; + cerr << endl; + if (!SubtractEstimatedBkg()) + return false; + } + } else { // fixed background given + if (!SubtractFixBkg()) + return false; + } + + UInt_t histoNo[2] = {forwardHistoNo[0], backwardHistoNo[0]}; + + // get the data range (fgb/lgb) for the current RUN block + if (!GetProperDataRange(runData, histoNo)) { + return false; + } + + // get the fit range for the current RUN block + GetProperFitRange(globalBlock); + + // everything looks fine, hence fill data set + Bool_t status; + switch(fHandleTag) { + case kFit: + status = PrepareFitData(); + break; + case kView: + status = PrepareViewData(runData, histoNo); + break; + default: + status = false; + break; + } + + // clean up + forwardHistoNo.clear(); + backwardHistoNo.clear(); + + return status; +} + +//-------------------------------------------------------------------------- +// SubtractFixBkg (private) +//-------------------------------------------------------------------------- +/** + *

Subtracts a fixed background from the raw data. The background is given + * in units of (1/bin); for the Asymmetry representation (1/ns) doesn't make too much sense. + * The error propagation is done the following way: it is assumed that the error of the background + * is Poisson like, i.e. \f$\Delta\mathrm{bkg} = \sqrt{\mathrm{bkg}}\f$. + * + * Error propagation: + * \f[ \Delta f_i^{\rm c} = \pm\left[ (\Delta f_i)^2 + (\Delta \mathrm{bkg})^2 \right]^{1/2} = + * \pm\left[ f_i + \mathrm{bkg} \right]^{1/2}, \f] + * where \f$ f_i^{\rm c} \f$ is the background corrected histogram, \f$ f_i \f$ the raw histogram + * and \f$ \mathrm{bkg} \f$ the fix given background. + * + * return: + * - true + * + */ +Bool_t PRunAsymmetryRRF::SubtractFixBkg() +{ + Double_t dval; + for (UInt_t i=0; iGetBkgFix(0); + + // keep the error, and make sure that the bin is NOT empty + if (fBackward[i] != 0.0) + dval = TMath::Sqrt(fBackward[i]); + else + dval = 1.0; + fBackwardErr.push_back(dval); + fBackward[i] -= fRunInfo->GetBkgFix(1); + } + + return true; +} + +//-------------------------------------------------------------------------- +// SubtractEstimatedBkg (private) +//-------------------------------------------------------------------------- +/** + *

Subtracts the background which is estimated from a given interval (typically before t0). + * + * The background corrected histogramms are: + * \f$ f_i^{\rm c} = f_i - \mathrm{bkg} \f$, where \f$ f_i \f$ is the raw data histogram, + * \f$ \mathrm{bkg} \f$ the background estimate, and \f$ f_i^{\rm c} \f$ background corrected + * histogram. The error on \f$ f_i^{\rm c} \f$ is + * \f[ \Delta f_i^{\rm c} = \pm \sqrt{ (\Delta f_i)^2 + (\Delta \mathrm{bkg})^2 } = + * \pm \sqrt{f_i + (\Delta \mathrm{bkg})^2} \f] + * The background error \f$ \Delta \mathrm{bkg} \f$ is + * \f[ \Delta \mathrm{bkg} = \pm\frac{1}{N}\left[\sum_{i=0}^N (\Delta f_i)^2\right]^{1/2} = + * \pm\frac{1}{N}\left[\sum_{i=0}^N f_i \right]^{1/2},\f] + * where \f$N\f$ is the number of bins over which the background is formed. + * + * return: + * - true + */ +Bool_t PRunAsymmetryRRF::SubtractEstimatedBkg() +{ + Double_t beamPeriod = 0.0; + + // check if data are from PSI, RAL, or TRIUMF + if (fRunInfo->GetInstitute()->Contains("psi")) + beamPeriod = ACCEL_PERIOD_PSI; + else if (fRunInfo->GetInstitute()->Contains("ral")) + beamPeriod = ACCEL_PERIOD_RAL; + else if (fRunInfo->GetInstitute()->Contains("triumf")) + beamPeriod = ACCEL_PERIOD_TRIUMF; + else + beamPeriod = 0.0; + + // check if start and end are in proper order + UInt_t start[2] = {fRunInfo->GetBkgRange(0), fRunInfo->GetBkgRange(2)}; + UInt_t end[2] = {fRunInfo->GetBkgRange(1), fRunInfo->GetBkgRange(3)}; + for (UInt_t i=0; i<2; i++) { + if (end[i] < start[i]) { + cout << endl << "PRunAsymmetryRRF::SubtractEstimatedBkg(): end = " << end[i] << " > start = " << start[i] << "! Will swap them!"; + UInt_t keep = end[i]; + end[i] = start[i]; + start[i] = keep; + } + } + + // calculate proper background range + for (UInt_t i=0; i<2; i++) { + if (beamPeriod != 0.0) { + Double_t timeBkg = (Double_t)(end[i]-start[i])*fTimeResolution; // length of the background intervall in time + UInt_t fullCycles = (UInt_t)(timeBkg/beamPeriod); // how many proton beam cylces can be placed within the proposed background intervall + // correct the end of the background intervall such that the background is as close as possible to a multiple of the proton cylce + end[i] = start[i] + (UInt_t) ((fullCycles*beamPeriod)/fTimeResolution); + cout << "PRunAsymmetryRRF::SubtractEstimatedBkg(): Background " << start[i] << ", " << end[i] << endl; + if (end[i] == start[i]) + end[i] = fRunInfo->GetBkgRange(2*i+1); + } + } + + // check if start is within histogram bounds + if ((start[0] < 0) || (start[0] >= fForward.size()) || + (start[1] < 0) || (start[1] >= fBackward.size())) { + cerr << endl << ">> PRunAsymmetryRRF::SubtractEstimatedBkg(): **ERROR** background bin values out of bound!"; + cerr << endl << ">> histo lengths (f/b) = (" << fForward.size() << "/" << fBackward.size() << ")."; + cerr << endl << ">> background start (f/b) = (" << start[0] << "/" << start[1] << ")."; + return false; + } + + // check if end is within histogram bounds + if ((end[0] < 0) || (end[0] >= fForward.size()) || + (end[1] < 0) || (end[1] >= fBackward.size())) { + cerr << endl << ">> PRunAsymmetryRRF::SubtractEstimatedBkg(): **ERROR** background bin values out of bound!"; + cerr << endl << ">> histo lengths (f/b) = (" << fForward.size() << "/" << fBackward.size() << ")."; + cerr << endl << ">> background end (f/b) = (" << end[0] << "/" << end[1] << ")."; + return false; + } + + // calculate background + Double_t bkg[2] = {0.0, 0.0}; + Double_t errBkg[2] = {0.0, 0.0}; + + // forward + for (UInt_t i=start[0]; i<=end[0]; i++) + bkg[0] += fForward[i]; + errBkg[0] = TMath::Sqrt(bkg[0])/(end[0] - start[0] + 1); + bkg[0] /= static_cast(end[0] - start[0] + 1); + cout << endl << ">> estimated forward histo background: " << bkg[0]; + + // backward + for (UInt_t i=start[1]; i<=end[1]; i++) + bkg[1] += fBackward[i]; + errBkg[1] = TMath::Sqrt(bkg[1])/(end[1] - start[1] + 1); + bkg[1] /= static_cast(end[1] - start[1] + 1); + cout << endl << ">> estimated backward histo background: " << bkg[1] << endl; + + // correct error for forward, backward + Double_t errVal = 0.0; + for (UInt_t i=0; i 0.0) + errVal = TMath::Sqrt(fForward[i]+errBkg[0]*errBkg[0]); + else + errVal = 1.0; + fForwardErr.push_back(errVal); + if (fBackward[i] > 0.0) + errVal = TMath::Sqrt(fBackward[i]+errBkg[1]*errBkg[1]); + else + errVal = 1.0; + fBackwardErr.push_back(errVal); + } + + // subtract background from data + for (UInt_t i=0; iSetBkgEstimated(bkg[0], 0); + fRunInfo->SetBkgEstimated(bkg[1], 1); + + return true; +} + +//-------------------------------------------------------------------------- +// PrepareFitData (protected) +//-------------------------------------------------------------------------- +/** + *

Take the pre-processed data (i.e. grouping and addrun are preformed, background correction already carried out) + * and form the asymmetry for fitting. + */ +Bool_t PRunAsymmetryRRF::PrepareFitData() +{ + // transform raw histo data. At this point, the raw data are already background corrected. + + // 1st: form the asymmetry of the original data + + // forward and backward detectors might have different fgb-t0 offset. Take the maximum of both. + Int_t fgbOffset = fGoodBins[0]-(Int_t)fT0s[0]; + if (fgbOffset < fGoodBins[2]-(Int_t)fT0s[1]) + fgbOffset = fGoodBins[2]-(Int_t)fT0s[1]; + // last good bin (lgb) is the minimum of forward/backward lgb + Int_t lgb_offset = fGoodBins[1]-(Int_t)fT0s[0]+fgbOffset; + if (lgb_offset < fGoodBins[3]-(Int_t)fT0s[1]+fgbOffset) + lgb_offset = fGoodBins[3]-(Int_t)fT0s[1]+fgbOffset; + + Int_t fgb = (Int_t)fT0s[0]+fgbOffset; + Int_t lgb = fgb + lgb_offset; + Int_t dt0 = (Int_t)fT0s[0]-(Int_t)fT0s[1]; + + PDoubleVector asym; + PDoubleVector asymErr; + Double_t asymVal, asymValErr; + Double_t ff, bb, eff, ebb; + for (Int_t i=fgb; i 0.0) + asymValErr = 2.0/pow((ff+bb),2.0)*sqrt(bb*bb*eff*eff+ff*ff*ebb*ebb); + else + asymValErr = 1.0; + asymErr.push_back(asymValErr); + } + + // 2nd: a_rrf = a * 2*cos(w_rrf*t + phi_rrf) + PMsrGlobalBlock *globalBlock = fMsrInfo->GetMsrGlobal(); + Double_t wRRF = globalBlock->GetRRFFreq("Mc"); + Double_t phaseRRF = globalBlock->GetRRFPhase()*TMath::TwoPi()/180.0; + + Double_t startTime = fTimeResolution * (Double_t)fgbOffset; + Double_t time=0.0; + for (UInt_t i=0; iTake the pre-processed data (i.e. grouping and addrun are preformed) and form the asymmetry for view representation. + * Before forming the asymmetry, the following checks will be performed: + * -# check if view packing is whished. + * -# check if data range is given, if not try to estimate one. + * -# check that data range is present, that it makes any sense. + * -# check that 'first good bin'-'t0' is the same for forward and backward histogram. If not adjust it. + * -# pack data (rebin). + * -# if packed forward size != backward size, truncate the longer one such that an asymmetry can be formed. + * -# calculate the asymmetry: \f$ A_i = (\alpha f_i^c-b_i^c)/(\alpha \beta f_i^c+b_i^c) \f$ + * -# calculate the asymmetry errors: \f$ \delta A_i = 2 \sqrt{(b_i^c)^2 (\delta f_i^c)^2 + (\delta b_i^c)^2 (f_i^c)^2}/(f_i^c+b_i^c)^2\f$ + * -# calculate the theory vector. + * + * \param runData raw run data needed to perform some crosschecks + * \param histoNo histogram number (within a run). histoNo[0]: forward histogram number, histNo[1]: backward histogram number + */ +Bool_t PRunAsymmetryRRF::PrepareViewData(PRawRunData* runData, UInt_t histoNo[2]) +{ + // feed the parameter vector + std::vector par; + PMsrParamList *paramList = fMsrInfo->GetMsrParamList(); + for (UInt_t i=0; isize(); i++) + par.push_back((*paramList)[i].fValue); + + // first get start data, end data, and t0 + Int_t start[2] = {fGoodBins[0], fGoodBins[2]}; + Int_t end[2] = {fGoodBins[1], fGoodBins[3]}; + Int_t t0[2] = {(Int_t)fT0s[0], (Int_t)fT0s[1]}; + + // check if the data ranges and t0's between forward/backward are compatible + Int_t fgb[2]; + if (start[0]-t0[0] != start[1]-t0[1]) { // wrong fgb aligning + if (abs(start[0]-t0[0]) > abs(start[1]-t0[1])) { + fgb[0] = start[0]; + fgb[1] = t0[1] + start[0]-t0[0]; + cerr << endl << ">> PRunAsymmetryRRF::PrepareViewData(): **WARNING** needed to shift backward fgb from "; + cerr << start[1] << " to " << fgb[1] << endl; + } else { + fgb[0] = t0[0] + start[1]-t0[1]; + fgb[1] = start[1]; + cerr << endl << ">> PRunAsymmetryRRF::PrepareViewData(): **WARNING** needed to shift forward fgb from "; + cerr << start[0] << " to " << fgb[0] << endl; + } + } else { // fgb aligning is correct + fgb[0] = start[0]; + fgb[1] = start[1]; + } + start[0] = fgb[0]; + start[1] = fgb[1]; + + // make sure that there are equal number of bins in forward and backward + UInt_t noOfBins0 = runData->GetDataBin(histoNo[0])->size()-start[0]; + UInt_t noOfBins1 = runData->GetDataBin(histoNo[1])->size()-start[1]; + if (noOfBins0 > noOfBins1) + noOfBins0 = noOfBins1; + end[0] = start[0] + noOfBins0; + end[1] = start[1] + noOfBins0; + + // check if start, end, and t0 make any sense + for (UInt_t i=0; i<2; i++) { + // 1st check if start is within proper bounds + if ((start[i] < 0) || (start[i] > (Int_t)runData->GetDataBin(histoNo[i])->size())) { + cerr << endl << ">> PRunAsymmetryRRF::PrepareViewData(): **ERROR** start data bin doesn't make any sense!"; + cerr << endl; + return false; + } + // 2nd check if end is within proper bounds + if ((end[i] < 0) || (end[i] > (Int_t)runData->GetDataBin(histoNo[i])->size())) { + cerr << endl << ">> PRunAsymmetryRRF::PrepareViewData(): **ERROR** end data bin doesn't make any sense!"; + cerr << endl; + return false; + } + // 3rd check if t0 is within proper bounds + if ((t0[i] < 0) || (t0[i] > (Int_t)runData->GetDataBin(histoNo[i])->size())) { + cerr << endl << ">> PRunAsymmetryRRF::PrepareViewData(): **ERROR** t0 data bin doesn't make any sense!"; + cerr << endl; + return false; + } + } + + // check if forward and backward histo have the same size, otherwise take the minimum size + UInt_t noOfBins = fForward.size(); + if (noOfBins > fBackward.size()) { + noOfBins = fBackward.size(); + } + + // form asymmetry including error propagation + Double_t asym, error; + Double_t f, b, ef, eb, alpha = 1.0, beta = 1.0; + + // get the proper alpha and beta + switch (fAlphaBetaTag) { + case 1: // alpha == 1, beta == 1 + alpha = 1.0; + beta = 1.0; + break; + case 2: // alpha != 1, beta == 1 + alpha = par[fRunInfo->GetAlphaParamNo()-1]; + beta = 1.0; + break; + case 3: // alpha == 1, beta != 1 + alpha = 1.0; + beta = par[fRunInfo->GetBetaParamNo()-1]; + break; + case 4: // alpha != 1, beta != 1 + alpha = par[fRunInfo->GetAlphaParamNo()-1]; + beta = par[fRunInfo->GetBetaParamNo()-1]; + break; + default: + break; + } + + PDoubleVector asymVec, asymErr; + Int_t dtBin = start[1]-start[0]; + for (Int_t i=start[0]; iGetMsrGlobal(); + Double_t wRRF = globalBlock->GetRRFFreq("Mc"); + Double_t phaseRRF = globalBlock->GetRRFPhase()*TMath::TwoPi()/180.0; + Double_t startTime=fTimeResolution*((Double_t)start[0]-t0[0]); + Double_t time = 0.0; + for (UInt_t i=0; iGetNoOfFuncs(); i++) { + fFuncValues[i] = fMsrInfo->EvalFunc(fMsrInfo->GetFuncNo(i), *fRunInfo->GetMap(), par); + } + + // calculate theory + UInt_t size = runData->GetDataBin(histoNo[0])->size(); + Double_t factor = 1.0; + if (fData.GetValue()->size() * 10 > runData->GetDataBin(histoNo[0])->size()) { + size = fData.GetValue()->size() * 10; + factor = (Double_t)runData->GetDataBin(histoNo[0])->size() / (Double_t)size; + } + fData.SetTheoryTimeStart(fData.GetDataTimeStart()); + fData.SetTheoryTimeStep(fTimeResolution*factor); + for (UInt_t i=0; iFunc(time, par, fFuncValues); + if (fabs(dval) > 10.0) { // dirty hack needs to be fixed!! + dval = 0.0; + } + fData.AppendTheoryValue(dval); + } + + // clean up + par.clear(); + + return true; +} + +//-------------------------------------------------------------------------- +// GetProperT0 (private) +//-------------------------------------------------------------------------- +/** + *

Get the proper t0 for the single histogram run. + * -# the t0 vector size = number of detectors (grouping) for forward. + * -# initialize t0's with -1 + * -# fill t0's from RUN block + * -# if t0's are missing (i.e. t0 == -1), try to fill from the GLOBAL block. + * -# if t0's are missing, try t0's from the data file + * -# if t0's are missing, try to estimate them + * + * \param runData pointer to the current RUN block entry from the msr-file + * \param globalBlock pointer to the GLOBLA block entry from the msr-file + * \param forwardHistoNo histogram number vector of forward; forwardHistoNo = msr-file forward + redGreen_offset - 1 + * \param backwardHistoNo histogram number vector of backwardward; backwardHistoNo = msr-file backward + redGreen_offset - 1 + * + * return: + * - true if everthing went smooth + * - false, otherwise. + */ +Bool_t PRunAsymmetryRRF::GetProperT0(PRawRunData* runData, PMsrGlobalBlock *globalBlock, PUIntVector &forwardHistoNo, PUIntVector &backwardHistoNo) +{ + // feed all T0's + // first init T0's, T0's are stored as (forward T0, backward T0, etc.) + fT0s.clear(); + fT0s.resize(2*forwardHistoNo.size()); + for (UInt_t i=0; iGetT0BinSize(); i++) { + fT0s[i] = fRunInfo->GetT0Bin(i); + } + + // fill in the missing T0's from the GLOBAL block section (if present) + for (UInt_t i=0; iGetT0BinSize(); i++) { + if (fT0s[i] == -1) { // i.e. not given in the RUN block section + fT0s[i] = globalBlock->GetT0Bin(i); + } + } + + // fill in the missing T0's from the data file, if not already present in the msr-file + for (UInt_t i=0; iGetT0Bin(forwardHistoNo[i]) > 0.0) { + fT0s[2*i] = runData->GetT0Bin(forwardHistoNo[i]); + fRunInfo->SetT0Bin(fT0s[2*i], 2*i); + } + } + for (UInt_t i=0; iGetT0Bin(backwardHistoNo[i]) > 0.0) { + fT0s[2*i+1] = runData->GetT0Bin(backwardHistoNo[i]); + fRunInfo->SetT0Bin(fT0s[2*i+1], 2*i+1); + } + } + + // fill in the T0's gaps, i.e. in case the T0's are NEITHER in the msr-file and NOR in the data file + for (UInt_t i=0; iGetT0BinEstimated(forwardHistoNo[i]); + fRunInfo->SetT0Bin(fT0s[2*i], 2*i); + + cerr << endl << ">> PRunAsymmetryRRF::GetProperT0(): **WARRNING** NO t0's found, neither in the run data nor in the msr-file!"; + cerr << endl << ">> run: " << fRunInfo->GetRunName()->Data(); + cerr << endl << ">> will try the estimated one: forward t0 = " << runData->GetT0BinEstimated(forwardHistoNo[i]); + cerr << endl << ">> NO WARRANTY THAT THIS OK!! For instance for LEM this is almost for sure rubbish!"; + cerr << endl; + } + } + for (UInt_t i=0; iGetT0BinEstimated(backwardHistoNo[i]); + fRunInfo->SetT0Bin(fT0s[2*i+1], 2*i+1); + + cerr << endl << ">> PRunAsymmetryRRF::GetProperT0(): **WARRNING** NO t0's found, neither in the run data nor in the msr-file!"; + cerr << endl << ">> run: " << fRunInfo->GetRunName()->Data(); + cerr << endl << ">> will try the estimated one: backward t0 = " << runData->GetT0BinEstimated(backwardHistoNo[i]); + cerr << endl << ">> NO WARRANTY THAT THIS OK!! For instance for LEM this is almost for sure rubbish!"; + cerr << endl; + } + } + + // check if t0 is within proper bounds + for (UInt_t i=0; i (Int_t)runData->GetDataBin(forwardHistoNo[i])->size())) { + cerr << endl << ">> PRunAsymmetryRRF::GetProperT0(): **ERROR** t0 data bin (" << fT0s[2*i] << ") doesn't make any sense!"; + cerr << endl << ">> forwardHistoNo " << forwardHistoNo[i]; + cerr << endl; + return false; + } + if ((fT0s[2*i+1] < 0) || (fT0s[2*i+1] > (Int_t)runData->GetDataBin(backwardHistoNo[i])->size())) { + cerr << endl << ">> PRunAsymmetryRRF::PrepareData(): **ERROR** t0 data bin (" << fT0s[2*i+1] << ") doesn't make any sense!"; + cerr << endl << ">> backwardHistoNo " << backwardHistoNo[i]; + cerr << endl; + return false; + } + } + + // check if addrun's are present, and if yes add the necessary t0's + if (fRunInfo->GetRunNameSize() > 1) { // runs to be added present + PRawRunData *addRunData; + fAddT0s.resize(fRunInfo->GetRunNameSize()-1); // resize to the number of addruns + for (UInt_t i=1; iGetRunNameSize(); i++) { + // get run to be added to the main one + addRunData = fRawData->GetRunData(*(fRunInfo->GetRunName(i))); + if (addRunData == 0) { // couldn't get run + cerr << endl << ">> PRunAsymmetryRRF::GetProperT0(): **ERROR** Couldn't get addrun " << fRunInfo->GetRunName(i)->Data() << "!"; + cerr << endl; + return false; + } + + // feed all T0's + // first init T0's, T0's are stored as (forward T0, backward T0, etc.) + fAddT0s[i-1].clear(); + fAddT0s[i-1].resize(2*forwardHistoNo.size()); + for (UInt_t j=0; jGetAddT0BinSize(i); j++) { + fAddT0s[i-1][j] = fRunInfo->GetAddT0Bin(i, j); + } + + // fill in the T0's from the data file, if not already present in the msr-file + for (UInt_t j=0; jGetT0Bin(forwardHistoNo[j]) > 0.0) { + fAddT0s[i-1][2*j] = addRunData->GetT0Bin(forwardHistoNo[j]); + fRunInfo->SetAddT0Bin(fAddT0s[i-1][2*j], i-1, 2*j); + } + } + for (UInt_t j=0; jGetT0Bin(backwardHistoNo[j]) > 0.0) { + fAddT0s[i-1][2*j+1] = addRunData->GetT0Bin(backwardHistoNo[j]); + fRunInfo->SetAddT0Bin(fAddT0s[i-1][2*j+1], i-1, 2*j+1); + } + } + + // fill in the T0's gaps, i.e. in case the T0's are NOT in the msr-file and NOT in the data file + for (UInt_t j=0; jGetT0BinEstimated(forwardHistoNo[j]); + fRunInfo->SetAddT0Bin(fAddT0s[i-1][2*j], i-1, 2*j); + + cerr << endl << ">> PRunAsymmetryRRF::GetProperT0(): **WARRNING** NO t0's found, neither in the run data nor in the msr-file!"; + cerr << endl << ">> run: " << fRunInfo->GetRunName(i)->Data(); + cerr << endl << ">> will try the estimated one: forward t0 = " << addRunData->GetT0BinEstimated(forwardHistoNo[j]); + cerr << endl << ">> NO WARRANTY THAT THIS OK!! For instance for LEM this is almost for sure rubbish!"; + cerr << endl; + } + } + for (UInt_t j=0; jGetT0BinEstimated(backwardHistoNo[j]); + fRunInfo->SetAddT0Bin(fAddT0s[i-1][2*j+1], i-1, 2*j+1); + + cerr << endl << ">> PRunAsymmetryRRF::GetProperT0(): **WARRNING** NO t0's found, neither in the run data nor in the msr-file!"; + cerr << endl << ">> run: " << fRunInfo->GetRunName(i)->Data(); + cerr << endl << ">> will try the estimated one: backward t0 = " << runData->GetT0BinEstimated(backwardHistoNo[j]); + cerr << endl << ">> NO WARRANTY THAT THIS OK!! For instance for LEM this is almost for sure rubbish!"; + cerr << endl; + } + } + } + } + + return true; +} + +//-------------------------------------------------------------------------- +// GetProperDataRange (private) +//-------------------------------------------------------------------------- +/** + *

Get the proper data range, i.e. first/last good bin (fgb/lgb). + * -# get fgb/lgb from the RUN block + * -# if fgb/lgb still undefined, try to get it from the GLOBAL block + * -# if fgb/lgb still undefined, try to estimate them. + * + * \param runData raw run data needed to perform some crosschecks + * \param histoNo histogram number (within a run). histoNo[0]: forward histogram number, histNo[1]: backward histogram number + * + * return: + * - true if everthing went smooth + * - false, otherwise. + */ +Bool_t PRunAsymmetryRRF::GetProperDataRange(PRawRunData* runData, UInt_t histoNo[2]) +{ + // first get start/end data + Int_t start[2] = {fRunInfo->GetDataRange(0), fRunInfo->GetDataRange(2)}; + Int_t end[2] = {fRunInfo->GetDataRange(1), fRunInfo->GetDataRange(3)}; + // check if data range has been provided in the RUN block. If not, try the GLOBAL block + if (start[0] == -1) { + start[0] = fMsrInfo->GetMsrGlobal()->GetDataRange(0); + } + if (start[1] == -1) { + start[1] = fMsrInfo->GetMsrGlobal()->GetDataRange(2); + } + if (end[0] == -1) { + end[0] = fMsrInfo->GetMsrGlobal()->GetDataRange(1); + } + if (end[1] == -1) { + end[1] = fMsrInfo->GetMsrGlobal()->GetDataRange(3); + } + + Double_t t0[2] = {fT0s[0], fT0s[1]}; + Int_t offset = (Int_t)(10.0e-3/fTimeResolution); // needed in case first good bin is not given, default = 10ns + + // check if data range has been provided, and if not try to estimate them + if (start[0] < 0) { + start[0] = (Int_t)t0[0]+offset; + fRunInfo->SetDataRange(start[0], 0); + cerr << endl << ">> PRunAsymmetryRRF::GetProperDataRange(): **WARNING** data range (forward) was not provided, will try data range start = t0+" << offset << "(=10ns) = " << start[0] << "."; + cerr << endl << ">> NO WARRANTY THAT THIS DOES MAKE ANY SENSE."; + cerr << endl; + } + if (start[1] < 0) { + start[1] = (Int_t)t0[1]+offset; + fRunInfo->SetDataRange(start[1], 2); + cerr << endl << ">> PRunAsymmetryRRF::GetProperDataRange(): **WARNING** data range (backward) was not provided, will try data range start = t0+" << offset << "(=10ns) = " << start[1] << "."; + cerr << endl << ">> NO WARRANTY THAT THIS DOES MAKE ANY SENSE."; + cerr << endl; + } + if (end[0] < 0) { + end[0] = runData->GetDataBin(histoNo[0])->size(); + fRunInfo->SetDataRange(end[0], 1); + cerr << endl << ">> PRunAsymmetryRRF::GetProperDataRange(): **WARNING** data range (forward) was not provided, will try data range end = " << end[0] << "."; + cerr << endl << ">> NO WARRANTY THAT THIS DOES MAKE ANY SENSE."; + cerr << endl; + } + if (end[1] < 0) { + end[1] = runData->GetDataBin(histoNo[1])->size(); + fRunInfo->SetDataRange(end[1], 3); + cerr << endl << ">> PRunAsymmetryRRF::GetProperDataRange(): **WARNING** data range (backward) was not provided, will try data range end = " << end[1] << "."; + cerr << endl << ">> NO WARRANTY THAT THIS DOES MAKE ANY SENSE."; + cerr << endl; + } + + // check if start, end, and t0 make any sense + // 1st check if start and end are in proper order + for (UInt_t i=0; i<2; i++) { + if (end[i] < start[i]) { // need to swap them + Int_t keep = end[i]; + end[i] = start[i]; + start[i] = keep; + } + // 2nd check if start is within proper bounds + if ((start[i] < 0) || (start[i] > (Int_t)runData->GetDataBin(histoNo[i])->size())) { + cerr << endl << ">> PRunAsymmetryRRF::GetProperDataRange(): **ERROR** start data bin doesn't make any sense!"; + cerr << endl; + return false; + } + // 3rd check if end is within proper bounds + if ((end[i] < 0) || (end[i] > (Int_t)runData->GetDataBin(histoNo[i])->size())) { + cerr << endl << ">> PRunAsymmetryRRF::GetProperDataRange(): **ERROR** end data bin doesn't make any sense!"; + cerr << endl; + return false; + } + // 4th check if t0 is within proper bounds + if ((t0[i] < 0) || (t0[i] > (Int_t)runData->GetDataBin(histoNo[i])->size())) { + cerr << endl << ">> PRunAsymmetryRRF::GetProperDataRange(): **ERROR** t0 data bin doesn't make any sense!"; + cerr << endl; + return false; + } + } + + // check that start-t0 is the same for forward as for backward, otherwise take max(start[i]-t0[i]) + if (fabs(static_cast(start[0])-t0[0]) > fabs(static_cast(start[1])-t0[1])){ + start[1] = static_cast(t0[1] + static_cast(start[0]) - t0[0]); + end[1] = static_cast(t0[1] + static_cast(end[0]) - t0[0]); + cerr << endl << ">> PRunAsymmetryRRF::GetProperDataRange **WARNING** needed to shift backward data range."; + cerr << endl << ">> given: " << fRunInfo->GetDataRange(2) << ", " << fRunInfo->GetDataRange(3); + cerr << endl << ">> used : " << start[1] << ", " << end[1]; + cerr << endl; + } + if (fabs(static_cast(start[0])-t0[0]) < fabs(static_cast(start[1])-t0[1])){ + start[0] = static_cast(t0[0] + static_cast(start[1]) - t0[1]); + end[0] = static_cast(t0[0] + static_cast(end[1]) - t0[1]); + cerr << endl << ">> PRunAsymmetryRRF::GetProperDataRange **WARNING** needed to shift forward data range."; + cerr << endl << ">> given: " << fRunInfo->GetDataRange(0) << ", " << fRunInfo->GetDataRange(1); + cerr << endl << ">> used : " << start[0] << ", " << end[0]; + cerr << endl; + } + + // keep good bins for potential latter use + fGoodBins[0] = start[0]; + fGoodBins[1] = end[0]; + fGoodBins[2] = start[1]; + fGoodBins[3] = end[1]; + + return true; +} + +//-------------------------------------------------------------------------- +// GetProperFitRange (private) +//-------------------------------------------------------------------------- +/** + *

Get the proper fit range. There are two possible fit range commands: + * fit given in (usec), or + * fit fgb+offset_0 lgb-offset_1 given in (bins), therefore it works the following way: + * -# get fit range assuming given in time from RUN block + * -# if fit range in RUN block is given in bins, replace start/end + * -# if fit range is NOT given yet, try fit range assuming given in time from GLOBAL block + * -# if fit range in GLOBAL block is given in bins, replace start/end + * -# if still no fit range is given, use fgb/lgb. + * + * \param globalBlock pointer to the GLOBAL block information form the msr-file. + */ +void PRunAsymmetryRRF::GetProperFitRange(PMsrGlobalBlock *globalBlock) +{ + // set fit start/end time; first check RUN Block + fFitStartTime = fRunInfo->GetFitRange(0); + fFitEndTime = fRunInfo->GetFitRange(1); + // if fit range is given in bins (and not time), the fit start/end time can be calculated at this point now + if (fRunInfo->IsFitRangeInBin()) { + fFitStartTime = (fGoodBins[0] + fRunInfo->GetFitRangeOffset(0) - fT0s[0]) * fTimeResolution; // (fgb+n0-t0)*dt + fFitEndTime = (fGoodBins[1] - fRunInfo->GetFitRangeOffset(1) - fT0s[0]) * fTimeResolution; // (lgb-n1-t0)*dt + // write these times back into the data structure. This way it is available when writting the log-file + fRunInfo->SetFitRange(fFitStartTime, 0); + fRunInfo->SetFitRange(fFitEndTime, 1); + } + if (fFitStartTime == PMUSR_UNDEFINED) { // fit start/end NOT found in the RUN block, check GLOBAL block + fFitStartTime = globalBlock->GetFitRange(0); + fFitEndTime = globalBlock->GetFitRange(1); + // if fit range is given in bins (and not time), the fit start/end time can be calculated at this point now + if (globalBlock->IsFitRangeInBin()) { + fFitStartTime = (fGoodBins[0] + globalBlock->GetFitRangeOffset(0) - fT0s[0]) * fTimeResolution; // (fgb+n0-t0)*dt + fFitEndTime = (fGoodBins[1] - globalBlock->GetFitRangeOffset(1) - fT0s[0]) * fTimeResolution; // (lgb-n1-t0)*dt + // write these times back into the data structure. This way it is available when writting the log-file + globalBlock->SetFitRange(fFitStartTime, 0); + globalBlock->SetFitRange(fFitEndTime, 1); + } + } + if ((fFitStartTime == PMUSR_UNDEFINED) || (fFitEndTime == PMUSR_UNDEFINED)) { + fFitStartTime = (fGoodBins[0] - fT0s[0]) * fTimeResolution; // (fgb-t0)*dt + fFitEndTime = (fGoodBins[1] - fT0s[0]) * fTimeResolution; // (lgb-t0)*dt + cerr << ">> PRunSingleHisto::GetProperFitRange(): **WARNING** Couldn't get fit start/end time!" << endl; + cerr << ">> Will set it to fgb/lgb which given in time is: " << fFitStartTime << "..." << fFitEndTime << " (usec)" << endl; + } +} diff --git a/src/classes/PRunBase.cpp b/src/classes/PRunBase.cpp index 7be37fe1..ffa15c86 100644 --- a/src/classes/PRunBase.cpp +++ b/src/classes/PRunBase.cpp @@ -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 * diff --git a/src/classes/PRunDataHandler.cpp b/src/classes/PRunDataHandler.cpp index bebc1cf0..b95c46fb 100644 --- a/src/classes/PRunDataHandler.cpp +++ b/src/classes/PRunDataHandler.cpp @@ -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 * diff --git a/src/classes/PRunListCollection.cpp b/src/classes/PRunListCollection.cpp index 4082e91e..c9c5e56f 100644 --- a/src/classes/PRunListCollection.cpp +++ b/src/classes/PRunListCollection.cpp @@ -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; iCleanUp(); + fRunSingleHistoRRFList[i]->~PRunSingleHistoRRF(); + } + fRunSingleHistoRRFList.clear(); + for (UInt_t i=0; iCleanUp(); fRunAsymmetryList[i]->~PRunAsymmetry(); } fRunAsymmetryList.clear(); + for (UInt_t i=0; iCleanUp(); + fRunAsymmetryRRFList[i]->~PRunAsymmetryRRF(); + } + fRunAsymmetryRRFList.clear(); + for (UInt_t i=0; iCleanUp(); 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; iSetFitRangeBin(fitRange); + for (UInt_t i=0; iSetFitRangeBin(fitRange); for (UInt_t i=0; iSetFitRangeBin(fitRange); + for (UInt_t i=0; iSetFitRangeBin(fitRange); for (UInt_t i=0; iSetFitRangeBin(fitRange); for (UInt_t i=0; iSetFitRange(fitRange); + for (UInt_t i=0; iSetFitRange(fitRange); for (UInt_t i=0; iSetFitRange(fitRange); + for (UInt_t i=0; iSetFitRange(fitRange); for (UInt_t i=0; iSetFitRange(fitRange); for (UInt_t i=0; i& pa return chisq; } +//-------------------------------------------------------------------------- +// GetSingleHistoRRFChisq (public) +//-------------------------------------------------------------------------- +/** + *

Calculates chi-square of all single histogram RRF runs of a msr-file. + * + * return: + * - chi-square of all single histogram RRF runs of the msr-file + * + * \param par fit parameter vector + */ +Double_t PRunListCollection::GetSingleHistoRRFChisq(const std::vector& par) const +{ + Double_t chisq = 0.0; + + for (UInt_t i=0; iCalcChiSquare(par); + + return chisq; +} + //-------------------------------------------------------------------------- // GetAsymmetryChisq (public) //-------------------------------------------------------------------------- @@ -219,6 +270,27 @@ Double_t PRunListCollection::GetAsymmetryChisq(const std::vector& par) return chisq; } +//-------------------------------------------------------------------------- +// GetAsymmetryRRFChisq (public) +//-------------------------------------------------------------------------- +/** + *

Calculates chi-square of all asymmetry RRF runs of a msr-file. + * + * return: + * - chi-square of all asymmetry RRF runs of the msr-file + * + * \param par fit parameter vector + */ +Double_t PRunListCollection::GetAsymmetryRRFChisq(const std::vector& par) const +{ + Double_t chisq = 0.0; + + for (UInt_t i=0; iCalcChiSquare(par); + + return chisq; +} + //-------------------------------------------------------------------------- // GetMuMinusChisq (public) //-------------------------------------------------------------------------- @@ -299,9 +371,15 @@ Double_t PRunListCollection::GetSingleHistoChisqExpected(const std::vectorCalcChiSquareExpected(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& 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; iGetMsrRunList()->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; iGetMsrRunList()->at(i).GetFitType() == type) + subIdx++; + } } // return the chisq of the single run @@ -353,9 +432,15 @@ Double_t PRunListCollection::GetSingleRunChisq(const std::vector& 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& par, *

Calculates log max-likelihood of all single histogram runs of a msr-file. * * return: - * - 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::vectorCalculates log max-likelihood of all single histogram RRF runs of a msr-file. + * + * return: + * - log max-likelihood of all single histogram runs of the msr-file + * + * \param par fit parameter vector + */ +Double_t PRunListCollection::GetSingleHistoRRFMaximumLikelihood(const std::vector& par) const +{ + Double_t mlh = 0.0; + + for (UInt_t i=0; iCalcMaxLikelihood(par); + + return mlh; +} + //-------------------------------------------------------------------------- // GetAsymmetryMaximumLikelihood (public) //-------------------------------------------------------------------------- @@ -412,6 +518,28 @@ Double_t PRunListCollection::GetAsymmetryMaximumLikelihood(const std::vector Since it is not clear yet how to handle asymmetry fits with max likelihood + * the chi square will be used! + * + * return: + * - chi-square of all asymmetry RRF runs of the msr-file + * + * \param par fit parameter vector + */ +Double_t PRunListCollection::GetAsymmetryRRFMaximumLikelihood(const std::vector& par) const +{ + Double_t mlh = 0.0; + + for (UInt_t i=0; iCalcChiSquare(par); + + return mlh; +} + //-------------------------------------------------------------------------- // GetMuMinusMaximumLikelihood (public) //-------------------------------------------------------------------------- @@ -419,7 +547,7 @@ Double_t PRunListCollection::GetAsymmetryMaximumLikelihood(const std::vectorCalculates log max-likelihood of all mu minus runs of a msr-file. * * return: - * - 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; iGetNoOfFitBins(); + for (UInt_t i=0; iGetNoOfFitBins(); + for (UInt_t i=0; iGetNoOfFitBins(); + for (UInt_t i=0; iGetNoOfFitBins(); + for (UInt_t i=0; iGetNoOfFitBins(); @@ -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) +//-------------------------------------------------------------------------- +/** + *

Get a processed single histogram RRF data set. + * + * return: + * - 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; iGetRunNo() == 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) +//-------------------------------------------------------------------------- +/** + *

Get a processed asymmetry RRF data set. + * + * return: + * - 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; iGetRunNo() == 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; } diff --git a/src/classes/PRunMuMinus.cpp b/src/classes/PRunMuMinus.cpp index a792f59b..bdcf7049 100644 --- a/src/classes/PRunMuMinus.cpp +++ b/src/classes/PRunMuMinus.cpp @@ -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 * diff --git a/src/classes/PRunNonMusr.cpp b/src/classes/PRunNonMusr.cpp index e8c5f3a8..2b188141 100644 --- a/src/classes/PRunNonMusr.cpp +++ b/src/classes/PRunNonMusr.cpp @@ -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 * diff --git a/src/classes/PRunSingleHisto.cpp b/src/classes/PRunSingleHisto.cpp index 7d8d14ee..2881a32d 100644 --- a/src/classes/PRunSingleHisto.cpp +++ b/src/classes/PRunSingleHisto.cpp @@ -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 * diff --git a/src/classes/PRunSingleHistoRRF.cpp b/src/classes/PRunSingleHistoRRF.cpp new file mode 100644 index 00000000..55f7ce02 --- /dev/null +++ b/src/classes/PRunSingleHistoRRF.cpp @@ -0,0 +1,1205 @@ +/*************************************************************************** + + PRunSingleHistoRRF.cpp + + 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. * + ***************************************************************************/ + +#ifdef HAVE_CONFIG_H +#include "config.h" +#endif + +#ifdef HAVE_GOMP +#include +#endif + +#include +#include +#include +using namespace std; + +#include +#include +#include +#include + +#include "PMusr.h" +#include "PFourier.h" +#include "PRunSingleHistoRRF.h" + +//-------------------------------------------------------------------------- +// Constructor +//-------------------------------------------------------------------------- +/** + *

Constructor + */ +PRunSingleHistoRRF::PRunSingleHistoRRF() : PRunBase() +{ + fNoOfFitBins = 0; + fBackground = 0; + fRRFPacking = -1; + + // the 2 following variables are need in case fit range is given in bins, and since + // the fit range can be changed in the command block, these variables need to be accessible + fGoodBins[0] = -1; + fGoodBins[1] = -1; +} + +//-------------------------------------------------------------------------- +// Constructor +//-------------------------------------------------------------------------- +/** + *

Constructor + * + * \param msrInfo pointer to the msr-file handler + * \param rawData raw run data + * \param runNo number of the run within the msr-file + * \param tag tag showing what shall be done: kFit == fitting, kView == viewing + */ +PRunSingleHistoRRF::PRunSingleHistoRRF(PMsrHandler *msrInfo, PRunDataHandler *rawData, UInt_t runNo, EPMusrHandleTag tag) : PRunBase(msrInfo, rawData, runNo, tag) +{ + fNoOfFitBins = 0; + + PMsrGlobalBlock *global = msrInfo->GetMsrGlobal(); + + if (!global->IsPresent()) { + cerr << endl << ">> PRunSingleHistoRRF::PRunSingleHistoRRF(): **SEVERE ERROR**: no GLOBAL-block present!"; + cerr << endl << ">> For Single Histo RRF the GLOBAL-block is mandatory! Please fix this first."; + cerr << endl; + fValid = false; + return; + } + + if (!global->GetRRFUnit().CompareTo("??")) { + cerr << endl << ">> PRunSingleHistoRRF::PRunSingleHistoRRF(): **SEVERE ERROR**: no RRF-Frequency found!"; + cerr << endl; + fValid = false; + return; + } + + fRRFPacking = global->GetRRFPacking(); + if (fRRFPacking == -1) { + cerr << endl << ">> PRunSingleHistoRRF::PRunSingleHistoRRF(): **SEVERE ERROR**: no RRF-Packing found!"; + cerr << endl; + fValid = false; + return; + } + + // the 2 following variables are need in case fit range is given in bins, and since + // the fit range can be changed in the command block, these variables need to be accessible + fGoodBins[0] = -1; + fGoodBins[1] = -1; + + if (!PrepareData()) { + cerr << endl << ">> PRunSingleHistoRRF::PRunSingleHistoRRF(): **SEVERE ERROR**: Couldn't prepare data for fitting!"; + cerr << endl << ">> This is very bad :-(, will quit ..."; + cerr << endl; + fValid = false; + } +} + +//-------------------------------------------------------------------------- +// Destructor +//-------------------------------------------------------------------------- +/** + *

Destructor + */ +PRunSingleHistoRRF::~PRunSingleHistoRRF() +{ + fForward.clear(); +} + +//-------------------------------------------------------------------------- +// CalcChiSquare (public) +//-------------------------------------------------------------------------- +/** + *

Calculate chi-square. + * + * return: + * - chisq value + * + * \param par parameter vector iterated by minuit2 + */ +Double_t PRunSingleHistoRRF::CalcChiSquare(const std::vector& par) +{ + Double_t chisq = 0.0; + Double_t diff = 0.0; + + // calculate functions + for (Int_t i=0; iGetNoOfFuncs(); i++) { + Int_t funcNo = fMsrInfo->GetFuncNo(i); + fFuncValues[i] = fMsrInfo->EvalFunc(funcNo, *fRunInfo->GetMap(), par); + } + + // calculate chi square + Double_t time(1.0); + Int_t i, N(static_cast(fData.GetValue()->size())); + + // In order not to have an IF in the next loop, determine the start and end bins for the fit range now + Int_t startTimeBin = static_cast(ceil((fFitStartTime - fData.GetDataTimeStart())/fData.GetDataTimeStep())); + if (startTimeBin < 0) + startTimeBin = 0; + Int_t endTimeBin = static_cast(floor((fFitEndTime - fData.GetDataTimeStart())/fData.GetDataTimeStep())) + 1; + if (endTimeBin > N) + endTimeBin = N; + + // Calculate the theory function once to ensure one function evaluation for the current set of parameters. + // This is needed for the LF and user functions where some non-thread-save calculations only need to be calculated once + // for a given set of parameters---which should be done outside of the parallelized loop. + // For all other functions it means a tiny and acceptable overhead. + time = fTheory->Func(time, par, fFuncValues); + + #ifdef HAVE_GOMP + Int_t chunk = (endTimeBin - startTimeBin)/omp_get_num_procs(); + if (chunk < 10) + chunk = 10; + #pragma omp parallel for default(shared) private(i,time,diff) schedule(dynamic,chunk) reduction(+:chisq) + #endif + for (i=startTimeBin; iat(i) - fTheory->Func(time, par, fFuncValues); + chisq += diff*diff / (fData.GetError()->at(i)*fData.GetError()->at(i)); + } + + return chisq; +} + +//-------------------------------------------------------------------------- +// CalcChiSquareExpected (public) +//-------------------------------------------------------------------------- +/** + *

Calculate expected chi-square. + * + * return: + * - chisq value + * + * \param par parameter vector iterated by minuit2 + */ +Double_t PRunSingleHistoRRF::CalcChiSquareExpected(const std::vector& par) +{ + Double_t chisq = 0.0; + Double_t diff = 0.0; + Double_t theo = 0.0; + + // calculate functions + for (Int_t i=0; iGetNoOfFuncs(); i++) { + Int_t funcNo = fMsrInfo->GetFuncNo(i); + fFuncValues[i] = fMsrInfo->EvalFunc(funcNo, *fRunInfo->GetMap(), par); + } + + // calculate chi square + Double_t time(1.0); + Int_t i, N(static_cast(fData.GetValue()->size())); + + // In order not to have an IF in the next loop, determine the start and end bins for the fit range now + Int_t startTimeBin = static_cast(ceil((fFitStartTime - fData.GetDataTimeStart())/fData.GetDataTimeStep())); + if (startTimeBin < 0) + startTimeBin = 0; + Int_t endTimeBin = static_cast(floor((fFitEndTime - fData.GetDataTimeStart())/fData.GetDataTimeStep())) + 1; + if (endTimeBin > N) + endTimeBin = N; + + // Calculate the theory function once to ensure one function evaluation for the current set of parameters. + // This is needed for the LF and user functions where some non-thread-save calculations only need to be calculated once + // for a given set of parameters---which should be done outside of the parallelized loop. + // For all other functions it means a tiny and acceptable overhead. + time = fTheory->Func(time, par, fFuncValues); + + #ifdef HAVE_GOMP + Int_t chunk = (endTimeBin - startTimeBin)/omp_get_num_procs(); + if (chunk < 10) + chunk = 10; + #pragma omp parallel for default(shared) private(i,time,diff) schedule(dynamic,chunk) reduction(+:chisq) + #endif + for (i=startTimeBin; i < endTimeBin; ++i) { + time = fData.GetDataTimeStart() + (Double_t)i*fData.GetDataTimeStep(); + theo = fTheory->Func(time, par, fFuncValues); + diff = fData.GetValue()->at(i) - theo; + chisq += diff*diff / theo; + } + + return chisq; +} + +//-------------------------------------------------------------------------- +// CalcMaxLikelihood (public) +//-------------------------------------------------------------------------- +/** + *

Calculate log maximum-likelihood. See http://pdg.lbl.gov/index.html + * + * return: + * - log maximum-likelihood value + * + * \param par parameter vector iterated by minuit2 + */ +Double_t PRunSingleHistoRRF::CalcMaxLikelihood(const std::vector& par) +{ + // not yet implemented + + return 0.0; +} + +//-------------------------------------------------------------------------- +// CalcTheory (public) +//-------------------------------------------------------------------------- +/** + *

Calculate theory for a given set of fit-parameters. + */ +void PRunSingleHistoRRF::CalcTheory() +{ + // feed the parameter vector + std::vector par; + PMsrParamList *paramList = fMsrInfo->GetMsrParamList(); + for (UInt_t i=0; isize(); i++) + par.push_back((*paramList)[i].fValue); + + // calculate functions + for (Int_t i=0; iGetNoOfFuncs(); i++) { + fFuncValues[i] = fMsrInfo->EvalFunc(fMsrInfo->GetFuncNo(i), *fRunInfo->GetMap(), par); + } + + // calculate theory + UInt_t size = fData.GetValue()->size(); + Double_t start = fData.GetDataTimeStart(); + Double_t resolution = fData.GetDataTimeStep(); + Double_t time; + for (UInt_t i=0; iFunc(time, par, fFuncValues)); + } + + // clean up + par.clear(); +} + +//-------------------------------------------------------------------------- +// GetNoOfFitBins (public) +//-------------------------------------------------------------------------- +/** + *

Calculate the number of fitted bins for the current fit range. + * + * return: number of fitted bins. + */ +UInt_t PRunSingleHistoRRF::GetNoOfFitBins() +{ + CalcNoOfFitBins(); + + return fNoOfFitBins; +} + +//-------------------------------------------------------------------------- +// SetFitRangeBin (public) +//-------------------------------------------------------------------------- +/** + *

Allows to change the fit range on the fly. Used in the COMMAND block. + * The syntax of the string is: FIT_RANGE fgb[+n00] lgb[-n01] [fgb[+n10] lgb[-n11] ... fgb[+nN0] lgb[-nN1]]. + * If only one pair of fgb/lgb is given, it is used for all runs in the RUN block section. + * If multiple fgb/lgb's are given, the number N has to be the number of RUN blocks in + * the msr-file. + * + *

nXY are offsets which can be used to shift, limit the fit range. + * + * \param fitRange string containing the necessary information. + */ +void PRunSingleHistoRRF::SetFitRangeBin(const TString fitRange) +{ + TObjArray *tok = 0; + TObjString *ostr = 0; + TString str; + Ssiz_t idx = -1; + Int_t offset = 0; + + tok = fitRange.Tokenize(" \t"); + + if (tok->GetEntries() == 3) { // structure FIT_RANGE fgb+n0 lgb-n1 + // handle fgb+n0 entry + ostr = (TObjString*) tok->At(1); + str = ostr->GetString(); + // check if there is an offset present + idx = str.First("+"); + if (idx != -1) { // offset present + str.Remove(0, idx+1); + if (str.IsFloat()) // if str is a valid number, convert is to an integer + offset = str.Atoi(); + } + fFitStartTime = (fGoodBins[0] + offset - fT0s[0]) * fTimeResolution; + + // handle lgb-n1 entry + ostr = (TObjString*) tok->At(2); + str = ostr->GetString(); + // check if there is an offset present + idx = str.First("-"); + if (idx != -1) { // offset present + str.Remove(0, idx+1); + if (str.IsFloat()) // if str is a valid number, convert is to an integer + offset = str.Atoi(); + } + fFitEndTime = (fGoodBins[1] - offset - fT0s[0]) * fTimeResolution; + } else if ((tok->GetEntries() > 3) && (tok->GetEntries() % 2 == 1)) { // structure FIT_RANGE fgb[+n00] lgb[-n01] [fgb[+n10] lgb[-n11] ... fgb[+nN0] lgb[-nN1]] + Int_t pos = 2*(fRunNo+1)-1; + + if (pos + 1 >= tok->GetEntries()) { + cerr << endl << ">> PRunSingleHistoRRF::SetFitRangeBin(): **ERROR** invalid FIT_RANGE command found: '" << fitRange << "'"; + cerr << endl << ">> will ignore it. Sorry ..." << endl; + } else { + // handle fgb+n0 entry + ostr = (TObjString*) tok->At(pos); + str = ostr->GetString(); + // check if there is an offset present + idx = str.First("+"); + if (idx != -1) { // offset present + str.Remove(0, idx+1); + if (str.IsFloat()) // if str is a valid number, convert is to an integer + offset = str.Atoi(); + } + fFitStartTime = (fGoodBins[0] + offset - fT0s[0]) * fTimeResolution; + + // handle lgb-n1 entry + ostr = (TObjString*) tok->At(pos+1); + str = ostr->GetString(); + // check if there is an offset present + idx = str.First("-"); + if (idx != -1) { // offset present + str.Remove(0, idx+1); + if (str.IsFloat()) // if str is a valid number, convert is to an integer + offset = str.Atoi(); + } + fFitEndTime = (fGoodBins[1] - offset - fT0s[0]) * fTimeResolution; + } + } else { // error + cerr << endl << ">> PRunSingleHistoRRF::SetFitRangeBin(): **ERROR** invalid FIT_RANGE command found: '" << fitRange << "'"; + cerr << endl << ">> will ignore it. Sorry ..." << endl; + } + + // clean up + if (tok) { + delete tok; + } +} + +//-------------------------------------------------------------------------- +// CalcNoOfFitBins (protected) +//-------------------------------------------------------------------------- +/** + *

Calculate the number of fitted bins for the current fit range. + */ +void PRunSingleHistoRRF::CalcNoOfFitBins() +{ + // In order not having to loop over all bins and to stay consistent with the chisq method, calculate the start and end bins explicitly + Int_t startTimeBin = static_cast(ceil((fFitStartTime - fData.GetDataTimeStart())/fData.GetDataTimeStep())); + if (startTimeBin < 0) + startTimeBin = 0; + Int_t endTimeBin = static_cast(floor((fFitEndTime - fData.GetDataTimeStart())/fData.GetDataTimeStep())) + 1; + if (endTimeBin > static_cast(fData.GetValue()->size())) + endTimeBin = fData.GetValue()->size(); + + if (endTimeBin > startTimeBin) + fNoOfFitBins = endTimeBin - startTimeBin; + else + fNoOfFitBins = 0; +} + +//-------------------------------------------------------------------------- +// PrepareData (protected) +//-------------------------------------------------------------------------- +/** + *

Prepare data for fitting or viewing. What is already processed at this stage: + * -# get proper raw run data + * -# get all needed forward histograms + * -# get time resolution + * -# get t0's and perform necessary cross checks (e.g. if t0 of msr-file (if present) are consistent with t0 of the data files, etc.) + * -# add runs (if addruns are present) + * -# group histograms (if grouping is present) + * + * return: + * - true if everthing went smooth + * - false, otherwise. + */ +Bool_t PRunSingleHistoRRF::PrepareData() +{ + Bool_t success = true; + + // keep the Global block info + PMsrGlobalBlock *globalBlock = fMsrInfo->GetMsrGlobal(); + + // get the proper run + PRawRunData* runData = fRawData->GetRunData(*fRunInfo->GetRunName()); + if (!runData) { // couldn't get run + cerr << endl << ">> PRunSingleHistoRRF::PrepareData(): **ERROR** Couldn't get run " << fRunInfo->GetRunName()->Data() << "!"; + cerr << endl; + return false; + } + + // collect histogram numbers + PUIntVector histoNo; // histoNo = msr-file forward + redGreen_offset - 1 + for (UInt_t i=0; iGetForwardHistoNoSize(); i++) { + histoNo.push_back(fRunInfo->GetForwardHistoNo(i)); + + if (!runData->IsPresent(histoNo[i])) { + cerr << endl << ">> PRunSingleHistoRRF::PrepareData(): **PANIC ERROR**:"; + cerr << endl << ">> histoNo found = " << histoNo[i] << ", which is NOT present in the data file!?!?"; + cerr << endl << ">> Will quit :-("; + cerr << endl; + histoNo.clear(); + return false; + } + } + + // keep the time resolution in (us) + fTimeResolution = runData->GetTimeResolution()/1.0e3; + cout.precision(10); + cout << endl << ">> PRunSingleHisto::PrepareData(): time resolution=" << fixed << runData->GetTimeResolution() << "(ns)" << endl; + + // get all the proper t0's and addt0's for the current RUN block + if (!GetProperT0(runData, globalBlock, histoNo)) { + return false; + } + + // keep the histo of each group at this point (addruns handled below) + vector forward; + forward.resize(histoNo.size()); // resize to number of groups + for (UInt_t i=0; iGetDataBin(histoNo[i])->size()); + forward[i] = *runData->GetDataBin(histoNo[i]); + } + + // check if there are runs to be added to the current one + if (fRunInfo->GetRunNameSize() > 1) { // runs to be added present + PRawRunData *addRunData; + for (UInt_t i=1; iGetRunNameSize(); i++) { + + // get run to be added to the main one + addRunData = fRawData->GetRunData(*fRunInfo->GetRunName(i)); + if (addRunData == 0) { // couldn't get run + cerr << endl << ">> PRunSingleHistoRRF::PrepareData(): **ERROR** Couldn't get addrun " << fRunInfo->GetRunName(i)->Data() << "!"; + cerr << endl; + return false; + } + + // add forward run + UInt_t addRunSize; + for (UInt_t k=0; kGetDataBin(histoNo[k])->size(); + for (UInt_t j=0; jGetDataBin(histoNo[k])->size(); j++) { // loop over the bin indices + // make sure that the index stays in the proper range + if ((j+(Int_t)fAddT0s[i-1][k]-(Int_t)fT0s[k] >= 0) && (j+(Int_t)fAddT0s[i-1][k]-(Int_t)fT0s[k] < addRunSize)) { + forward[k][j] += addRunData->GetDataBin(histoNo[k])->at(j+(Int_t)fAddT0s[i-1][k]-(Int_t)fT0s[k]); + } + } + } + } + } + + // set forward histo data of the first group + fForward.resize(forward[0].size()); + for (UInt_t i=0; iGetDataBin(histoNo[i])->size(); j++) { // loop over the bin indices + // make sure that the index stays within proper range + if ((j+fT0s[i]-fT0s[0] >= 0) && (j+fT0s[i]-fT0s[0] < runData->GetDataBin(histoNo[i])->size())) { + fForward[j] += forward[i][j+(Int_t)fT0s[i]-(Int_t)fT0s[0]]; + } + } + } + + // get the data range (fgb/lgb) for the current RUN block + if (!GetProperDataRange()) { + return false; + } + + // get the fit range for the current RUN block + GetProperFitRange(globalBlock); + + // do the more fit/view specific stuff + if (fHandleTag == kFit) + success = PrepareFitData(runData, histoNo[0]); + else if (fHandleTag == kView) + success = PrepareViewData(runData, histoNo[0]); + else + success = false; + + // cleanup + histoNo.clear(); + + return success; +} + +//-------------------------------------------------------------------------- +// PrepareFitData (protected) +//-------------------------------------------------------------------------- +/** + *

Take the pre-processed data (i.e. grouping and addrun are preformed) and form the RRF histogram for fitting. + * The following steps are preformed: + * -# get fit start/stop time + * -# check that 'first good data bin', 'last good data bin', and 't0' make any sense + * -# check how the background shall be handled, i.e. fitted, subtracted from background estimate data range, or subtacted from a given fixed background. + * -# estimate N0 + * -# RRF transformation + * -# packing (i.e rebinning) + * + * return: + * - true, if everything went smooth + * - false, otherwise + * + * \param runData raw run data handler + * \param histoNo forward histogram number + */ +Bool_t PRunSingleHistoRRF::PrepareFitData(PRawRunData* runData, const UInt_t histoNo) +{ + // keep the raw data for the RRF asymmetry error estimate for later + PDoubleVector rawNt; + for (UInt_t i=0; i freqMax=" << freqMax << " (MHz)" << endl; + + // "optimal packing" + Double_t optNoPoints = 8; + if (freqMax < 271.0) // < 271 MHz, i.e ~ 2T + optNoPoints = 5; + cout << "info> optimal packing: " << (Int_t)(1.0 / (fTimeResolution*(freqMax - fMsrInfo->GetMsrGlobal()->GetRRFFreq("MHz"))) / optNoPoints); + + // initially fForward is the "raw data set" (i.e. grouped histo and raw runs already added) to be fitted. This means fForward = N(t) at this point. + + // 1) check how the background shall be handled + // subtract background from histogramms ------------------------------------------ + if (fRunInfo->GetBkgFix(0) == PMUSR_UNDEFINED) { // no fixed background given + if (fRunInfo->GetBkgRange(0) >= 0) { + if (!EstimateBkg(histoNo)) + return false; + } else { // no background given to do the job, try estimate + fRunInfo->SetBkgRange(static_cast(fT0s[0]*0.1), 0); + fRunInfo->SetBkgRange(static_cast(fT0s[0]*0.6), 1); + cerr << endl << ">> PRunSingleHistoRRF::PrepareFitData(): **WARNING** Neither fix background nor background bins are given!"; + cerr << endl << ">> Will try the following: bkg start = " << fRunInfo->GetBkgRange(0) << ", bkg end = " << fRunInfo->GetBkgRange(1); + cerr << endl << ">> NO WARRANTY THAT THIS MAKES ANY SENSE! Better check ..."; + cerr << endl; + if (!EstimateBkg(histoNo)) + return false; + } + } else { // fixed background given + for (UInt_t i=0; iGetBkgFix(0); + } + fBackground = fRunInfo->GetBkgFix(0); + } + // here fForward = N(t) - Nbkg + + Int_t t0 = (Int_t)fT0s[0]; + + // 2) N(t) - Nbkg -> exp(+t/tau) [N(t)-Nbkg] + Double_t startTime = fTimeResolution * ((Double_t)fGoodBins[0] - (Double_t)t0); + + Double_t time_tau=0.0; + Double_t exp_t_tau=0.0; + for (Int_t i=fGoodBins[0]; i 0.0) + fW.push_back(1.0/(fMerr[i]*fMerr[i])); + else + fW.push_back(0.0); + } + // now fForward = exp(+t/tau) [N(t)-Nbkg] = M(t) + + // 3) estimate N0 + Double_t errN0 = 0.0; + Double_t n0 = EstimateN0(errN0, freqMax); + + // 4a) A(t) = exp(+t/tau) [N(t)-Nbkg] / N0 - 1.0 + for (Int_t i=fGoodBins[0]; i<=fGoodBins[1]; i++) { + fForward[i] = fForward[i] / n0 - 1.0; + } + + // 4b) error estimate of A(t): errA(t) = exp(+t/tau)/N0 sqrt( N(t) + ([N(t)-N_bkg]/N0)^2 errN0^2 ) + for (Int_t i=fGoodBins[0]; i<=fGoodBins[1]; i++) { + time_tau = (startTime + fTimeResolution * (i - fGoodBins[0])) / PMUON_LIFETIME; + exp_t_tau = exp(time_tau); + fAerr.push_back(exp_t_tau/n0*sqrt(rawNt[i]+pow(((rawNt[i]-fBackground)/n0)*errN0,2.0))); + } + + // 5) rotate A(t): A(t) -> 2* A(t) * cos(wRRF t + phiRRF), the factor 2.0 is needed since the high frequency part is suppressed. + PMsrGlobalBlock *globalBlock = fMsrInfo->GetMsrGlobal(); + Double_t wRRF = globalBlock->GetRRFFreq("Mc"); + Double_t phaseRRF = globalBlock->GetRRFPhase()*TMath::TwoPi()/180.0; + Double_t time = 0.0; + for (Int_t i=fGoodBins[0]; i<=fGoodBins[1]; i++) { + time = startTime + fTimeResolution * ((Double_t)i - (Double_t)fGoodBins[0]); + fForward[i] *= 2.0*cos(wRRF * time + phaseRRF); + } + + // 6) RRF packing + Double_t dval=0.0; + for (Int_t i=fGoodBins[0]; i<=fGoodBins[1]; i++) { + if (fRRFPacking == 1) { + fData.AppendValue(fForward[i]); + } else { // RRF packing > 1 + if (((i-fGoodBins[0]) % fRRFPacking == 0) && (i != fGoodBins[0])) { // fill data + dval /= fRRFPacking; + fData.AppendValue(dval); + // reset dval + dval = 0.0; + } + dval += fForward[i]; + } + } + + // 7) estimate packed RRF errors (see log-book p.204) + // the error estimate of the unpacked RRF asymmetry is: errA_RRF(t) \simeq exp(t/tau)/N0 sqrt( [N(t) + ((N(t)-N_bkg)/N0)^2 errN0^2] ) + dval = 0.0; + // the packed RRF asymmetry error + for (Int_t i=fGoodBins[0]; i<=fGoodBins[1]; i++) { + if (((i-fGoodBins[0]) % fRRFPacking == 0) && (i != fGoodBins[0])) { // fill data + fData.AppendErrorValue(sqrt(2.0*dval)/fRRFPacking); // the factor 2.0 is needed since the high frequency part is suppressed. + dval = 0.0; + } + dval += fAerr[i-fGoodBins[0]]*fAerr[i-fGoodBins[0]]; + } + + // set start time and time step + fData.SetDataTimeStart(fTimeResolution*((Double_t)fGoodBins[0]-(Double_t)t0+(Double_t)(fRRFPacking-1)/2.0)); + fData.SetDataTimeStep(fTimeResolution*fRRFPacking); + + CalcNoOfFitBins(); + + return true; +} + +//-------------------------------------------------------------------------- +// PrepareViewData (protected) +//-------------------------------------------------------------------------- +/** + *

Take the pre-processed data (i.e. grouping and addrun are preformed) and form the histogram for viewing + * with life time correction, i.e. the exponential decay is removed. + *

The following steps are preformed: + * -# check if view packing is whished. + * -# check that 'first good data bin', 'last good data bin', and 't0' makes any sense + * -# transform data sets (see below). + * -# calculate theory + * + * return: + * - true, if everything went smooth + * - false, otherwise + * + * \param runData raw run data handler + * \param histoNo forward histogram number + */ +Bool_t PRunSingleHistoRRF::PrepareViewData(PRawRunData* runData, const UInt_t histoNo) +{ + // -------------- + // prepare data + // -------------- + + // prepare RRF single histo + PrepareFitData(runData, histoNo); + + // check for view packing + Int_t viewPacking = fMsrInfo->GetMsrPlotList()->at(0).fViewPacking; + if (viewPacking > 0) { + if (viewPacking < fRRFPacking) { + cerr << ">> PRunSingleHistoRRF::PrepareViewData(): **WARNING** Found View Packing (" << viewPacking << ") < RRF Packing (" << fRRFPacking << ")."; + cerr << ">> Will ignore View Packing." << endl; + } else { + // STILL MISSING + } + } + + // -------------- + // prepare theory + // -------------- + + // feed the parameter vector + std::vector par; + PMsrParamList *paramList = fMsrInfo->GetMsrParamList(); + for (UInt_t i=0; isize(); i++) + par.push_back((*paramList)[i].fValue); + + // calculate functions + for (Int_t i=0; iGetNoOfFuncs(); i++) { + fFuncValues[i] = fMsrInfo->EvalFunc(fMsrInfo->GetFuncNo(i), *fRunInfo->GetMap(), par); + } + + // check if a finer binning for the theory is needed + UInt_t size = fForward.size(); + Double_t factor = 1.0; + Double_t time = 0.0; + Double_t theoryValue = 0.0; + + if (fData.GetValue()->size() * 10 > fForward.size()) { + size = fData.GetValue()->size() * 10; + factor = (Double_t)fForward.size() / (Double_t)size; + } + fData.SetTheoryTimeStart(fData.GetDataTimeStart()); + fData.SetTheoryTimeStep(fTimeResolution*factor); + + // calculate theory + for (UInt_t i=0; iFunc(time, par, fFuncValues); + if (fabs(theoryValue) > 10.0) { // dirty hack needs to be fixed!! + theoryValue = 0.0; + } + fData.AppendTheoryValue(theoryValue); + } + + return true; +} + +//-------------------------------------------------------------------------- +// GetProperT0 (private) +//-------------------------------------------------------------------------- +/** + *

Get the proper t0 for the single histogram run. + * -# the t0 vector size = number of detectors (grouping) for forward. + * -# initialize t0's with -1 + * -# fill t0's from RUN block + * -# if t0's are missing (i.e. t0 == -1), try to fill from the GLOBAL block. + * -# if t0's are missing, try t0's from the data file + * -# if t0's are missing, try to estimate them + * + * \param runData pointer to the current RUN block entry from the msr-file + * \param globalBlock pointer to the GLOBLA block entry from the msr-file + * \param histoNo histogram number vector of forward; histoNo = msr-file forward + redGreen_offset - 1 + * + * return: + * - true if everthing went smooth + * - false, otherwise. + */ +Bool_t PRunSingleHistoRRF::GetProperT0(PRawRunData* runData, PMsrGlobalBlock *globalBlock, PUIntVector &histoNo) +{ + // feed all T0's + // first init T0's, T0's are stored as (forward T0, backward T0, etc.) + fT0s.clear(); + fT0s.resize(histoNo.size()); + for (UInt_t i=0; iGetT0BinSize(); i++) { + fT0s[i] = fRunInfo->GetT0Bin(i); + } + + // fill in the T0's from the GLOBAL block section (if present) + for (UInt_t i=0; iGetT0BinSize(); i++) { + if (fT0s[i] == -1) { // i.e. not given in the RUN block section + fT0s[i] = globalBlock->GetT0Bin(i); + } + } + + // fill in the T0's from the data file, if not already present in the msr-file + for (UInt_t i=0; iGetT0Bin(histoNo[i]) > 0.0) { + fT0s[i] = runData->GetT0Bin(histoNo[i]); + fRunInfo->SetT0Bin(fT0s[i], i); // keep value for the msr-file + } + } + } + + // fill in the T0's gaps, i.e. in case the T0's are NOT in the msr-file and NOT in the data file + for (UInt_t i=0; iGetT0BinEstimated(histoNo[i]); + fRunInfo->SetT0Bin(fT0s[i], i); // keep value for the msr-file + + cerr << endl << ">> PRunSingleHistoRRF::GetProperT0(): **WARRNING** NO t0's found, neither in the run data nor in the msr-file!"; + cerr << endl << ">> run: " << fRunInfo->GetRunName()->Data(); + cerr << endl << ">> will try the estimated one: forward t0 = " << runData->GetT0BinEstimated(histoNo[i]); + cerr << endl << ">> NO WARRANTY THAT THIS OK!! For instance for LEM this is almost for sure rubbish!"; + cerr << endl; + } + } + + // check if t0 is within proper bounds + for (UInt_t i=0; iGetForwardHistoNoSize(); i++) { + if ((fT0s[i] < 0) || (fT0s[i] > (Int_t)runData->GetDataBin(histoNo[i])->size())) { + cerr << endl << ">> PRunSingleHistoRRF::GetProperT0(): **ERROR** t0 data bin (" << fT0s[i] << ") doesn't make any sense!"; + cerr << endl; + return false; + } + } + + // check if there are runs to be added to the current one. If yes keep the needed t0's + if (fRunInfo->GetRunNameSize() > 1) { // runs to be added present + PRawRunData *addRunData; + fAddT0s.resize(fRunInfo->GetRunNameSize()-1); // resize to the number of addruns + for (UInt_t i=1; iGetRunNameSize(); i++) { + + // get run to be added to the main one + addRunData = fRawData->GetRunData(*fRunInfo->GetRunName(i)); + if (addRunData == 0) { // couldn't get run + cerr << endl << ">> PRunSingleHistoRRF::GetProperT0(): **ERROR** Couldn't get addrun " << fRunInfo->GetRunName(i)->Data() << "!"; + cerr << endl; + return false; + } + + // feed all T0's + // first init T0's, T0's are stored as (forward T0, backward T0, etc.) + fAddT0s[i-1].resize(histoNo.size()); + for (UInt_t j=0; jGetT0BinSize(); j++) { + fAddT0s[i-1][j] = fRunInfo->GetAddT0Bin(i-1,j); // addRunIdx starts at 0 + } + + // fill in the T0's from the data file, if not already present in the msr-file + for (UInt_t j=0; jGetT0Bin(histoNo[j]) > 0.0) { + fAddT0s[i-1][j] = addRunData->GetT0Bin(histoNo[j]); + fRunInfo->SetAddT0Bin(fAddT0s[i-1][j], i-1, j); // keep value for the msr-file + } + } + + // fill in the T0's gaps, i.e. in case the T0's are NOT in the msr-file and NOT in the data file + for (UInt_t j=0; jGetT0BinEstimated(histoNo[j]); + fRunInfo->SetAddT0Bin(fAddT0s[i-1][j], i-1, j); // keep value for the msr-file + + cerr << endl << ">> PRunSingleHistoRRF::GetProperT0(): **WARRNING** NO t0's found, neither in the run data nor in the msr-file!"; + cerr << endl << ">> run: " << fRunInfo->GetRunName(i)->Data(); + cerr << endl << ">> will try the estimated one: forward t0 = " << addRunData->GetT0BinEstimated(histoNo[j]); + cerr << endl << ">> NO WARRANTY THAT THIS OK!! For instance for LEM this is almost for sure rubbish!"; + cerr << endl; + } + } + + // check if t0 is within proper bounds + for (UInt_t j=0; jGetForwardHistoNoSize(); j++) { + if ((fAddT0s[i-1][j] < 0) || (fAddT0s[i-1][j] > (Int_t)addRunData->GetDataBin(histoNo[j])->size())) { + cerr << endl << ">> PRunSingleHistoRRF::GetProperT0(): **ERROR** addt0 data bin (" << fAddT0s[i-1][j] << ") doesn't make any sense!"; + cerr << endl; + return false; + } + } + } + } + + return true; +} + +//-------------------------------------------------------------------------- +// GetProperDataRange (private) +//-------------------------------------------------------------------------- +/** + *

Get the proper data range, i.e. first/last good bin (fgb/lgb). + * -# get fgb/lgb from the RUN block + * -# if fgb/lgb still undefined, try to get it from the GLOBAL block + * -# if fgb/lgb still undefined, try to estimate them. + * + * return: + * - true if everthing went smooth + * - false, otherwise. + */ +Bool_t PRunSingleHistoRRF::GetProperDataRange() +{ + // get start/end data + Int_t start; + Int_t end; + start = fRunInfo->GetDataRange(0); + end = fRunInfo->GetDataRange(1); + + // check if data range has been given in the RUN block, if not try to get it from the GLOBAL block + if (start < 0) { + start = fMsrInfo->GetMsrGlobal()->GetDataRange(0); + } + if (end < 0) { + end = fMsrInfo->GetMsrGlobal()->GetDataRange(1); + } + + // check if data range has been provided, and if not try to estimate them + if (start < 0) { + Int_t offset = (Int_t)(10.0e-3/fTimeResolution); + start = (Int_t)fT0s[0]+offset; + fRunInfo->SetDataRange(start, 0); + cerr << endl << ">> PRunSingleHistoRRF::GetProperDataRange(): **WARNING** data range was not provided, will try data range start = t0+" << offset << "(=10ns) = " << start << "."; + cerr << endl << ">> NO WARRANTY THAT THIS DOES MAKE ANY SENSE."; + cerr << endl; + } + if (end < 0) { + end = fForward.size(); + fRunInfo->SetDataRange(end, 1); + cerr << endl << ">> PRunSingleHistoRRF::GetProperDataRange(): **WARNING** data range was not provided, will try data range end = " << end << "."; + cerr << endl << ">> NO WARRANTY THAT THIS DOES MAKE ANY SENSE."; + cerr << endl; + } + + // check if start and end make any sense + // 1st check if start and end are in proper order + if (end < start) { // need to swap them + Int_t keep = end; + end = start; + start = keep; + } + // 2nd check if start is within proper bounds + if ((start < 0) || (start > (Int_t)fForward.size())) { + cerr << endl << ">> PRunSingleHistoRRF::GetProperDataRange(): **ERROR** start data bin (" << start << ") doesn't make any sense!"; + cerr << endl; + return false; + } + // 3rd check if end is within proper bounds + if ((end < 0) || (end > (Int_t)fForward.size())) { + cerr << endl << ">> PRunSingleHistoRRF::GetProperDataRange(): **ERROR** end data bin (" << end << ") doesn't make any sense!"; + cerr << endl; + return false; + } + + // keep good bins for potential later use + fGoodBins[0] = start; + fGoodBins[1] = end; + + return true; +} + +//-------------------------------------------------------------------------- +// GetProperFitRange (private) +//-------------------------------------------------------------------------- +/** + *

Get the proper fit range. There are two possible fit range commands: + * fit given in (usec), or + * fit fgb+offset_0 lgb-offset_1 given in (bins), therefore it works the following way: + * -# get fit range assuming given in time from RUN block + * -# if fit range in RUN block is given in bins, replace start/end + * -# if fit range is NOT given yet, try fit range assuming given in time from GLOBAL block + * -# if fit range in GLOBAL block is given in bins, replace start/end + * -# if still no fit range is given, use fgb/lgb. + * + * \param globalBlock pointer to the GLOBAL block information form the msr-file. + */ +void PRunSingleHistoRRF::GetProperFitRange(PMsrGlobalBlock *globalBlock) +{ + // set fit start/end time; first check RUN Block + fFitStartTime = fRunInfo->GetFitRange(0); + fFitEndTime = fRunInfo->GetFitRange(1); + // if fit range is given in bins (and not time), the fit start/end time can be calculated at this point now + if (fRunInfo->IsFitRangeInBin()) { + fFitStartTime = (fGoodBins[0] + fRunInfo->GetFitRangeOffset(0) - fT0s[0]) * fTimeResolution; // (fgb+n0-t0)*dt + fFitEndTime = (fGoodBins[1] - fRunInfo->GetFitRangeOffset(1) - fT0s[0]) * fTimeResolution; // (lgb-n1-t0)*dt + // write these times back into the data structure. This way it is available when writting the log-file + fRunInfo->SetFitRange(fFitStartTime, 0); + fRunInfo->SetFitRange(fFitEndTime, 1); + } + if (fFitStartTime == PMUSR_UNDEFINED) { // fit start/end NOT found in the RUN block, check GLOBAL block + fFitStartTime = globalBlock->GetFitRange(0); + fFitEndTime = globalBlock->GetFitRange(1); + // if fit range is given in bins (and not time), the fit start/end time can be calculated at this point now + if (globalBlock->IsFitRangeInBin()) { + fFitStartTime = (fGoodBins[0] + globalBlock->GetFitRangeOffset(0) - fT0s[0]) * fTimeResolution; // (fgb+n0-t0)*dt + fFitEndTime = (fGoodBins[1] - globalBlock->GetFitRangeOffset(1) - fT0s[0]) * fTimeResolution; // (lgb-n1-t0)*dt + // write these times back into the data structure. This way it is available when writting the log-file + globalBlock->SetFitRange(fFitStartTime, 0); + globalBlock->SetFitRange(fFitEndTime, 1); + } + } + if ((fFitStartTime == PMUSR_UNDEFINED) || (fFitEndTime == PMUSR_UNDEFINED)) { + fFitStartTime = (fGoodBins[0] - fT0s[0]) * fTimeResolution; // (fgb-t0)*dt + fFitEndTime = (fGoodBins[1] - fT0s[0]) * fTimeResolution; // (lgb-t0)*dt + cerr << ">> PRunSingleHistoRRF::GetProperFitRange(): **WARNING** Couldn't get fit start/end time!" << endl; + cerr << ">> Will set it to fgb/lgb which given in time is: " << fFitStartTime << "..." << fFitEndTime << " (usec)" << endl; + } +} + +//-------------------------------------------------------------------------- +// GetMainFrequency (private) +//-------------------------------------------------------------------------- +/** + *

gets the maximum frequency of the given data. + * + * \param raw data set. + */ +Double_t PRunSingleHistoRRF::GetMainFrequency(PDoubleVector &data) +{ + Double_t freqMax = 0.0; + + // create histo + Double_t startTime = (fGoodBins[0]-fT0s[0]) * fTimeResolution; + Int_t noOfBins = fGoodBins[1]-fGoodBins[0]+1; + TH1F *histo = new TH1F("data", "data", noOfBins, startTime-fTimeResolution/2.0, startTime+fTimeResolution/2.0+noOfBins*fTimeResolution); + for (Int_t i=fGoodBins[0]; i<=fGoodBins[1]; i++) { + histo->SetBinContent(i-fGoodBins[0]+1, data[i]); + } + + // Fourier transform + PFourier *ft = new PFourier(histo, FOURIER_UNIT_FREQ); + ft->Transform(F_APODIZATION_STRONG); + + // find frequency maximum + TH1F *power = ft->GetPowerFourier(); + Double_t maxFreqVal = 0.0; + for (Int_t i=1; iGetNbinsX(); i++) { + // ignore dc part at 0 frequency + if (iGetNbinsX()-1) { + if (power->GetBinContent(i)>power->GetBinContent(i+1)) + continue; + } + // check for maximum + if (power->GetBinContent(i) > maxFreqVal) { + maxFreqVal = power->GetBinContent(i); + freqMax = power->GetBinCenter(i); + } + } + + // clean up + if (power) + delete power; + if (ft) + delete ft; + if (histo) + delete histo; + + return freqMax; +} + +//-------------------------------------------------------------------------- +// EstimateN0 (private) +//-------------------------------------------------------------------------- +/** + *

Estimate the N0 for the given run. + * + * \param errN0 + */ +Double_t PRunSingleHistoRRF::EstimateN0(Double_t &errN0, Double_t freqMax) +{ + // endBin is estimated such that the number of full cycles (according to the maximum frequency of the data) + // is approximately the time fN0EstimateEndTime. + Int_t endBin = (Int_t)round(fN0EstimateEndTime / fTimeResolution * ceil(freqMax)/freqMax); + + Double_t n0 = 0.0; + Double_t wN = 0.0; + for (Int_t i=0; i PRunSingleHistoRRF::EstimateN0(): N0=" << n0 << "(" << errN0 << ")" << endl; + + return n0; +} + +//-------------------------------------------------------------------------- +// EstimatBkg (private) +//-------------------------------------------------------------------------- +/** + *

Estimate the background for a given interval. + * + * return: + * - true, if everything went smooth + * - false, otherwise + * + * \param histoNo forward histogram number of the run + */ +Bool_t PRunSingleHistoRRF::EstimateBkg(UInt_t histoNo) +{ + Double_t beamPeriod = 0.0; + + // check if data are from PSI, RAL, or TRIUMF + if (fRunInfo->GetInstitute()->Contains("psi")) + beamPeriod = ACCEL_PERIOD_PSI; + else if (fRunInfo->GetInstitute()->Contains("ral")) + beamPeriod = ACCEL_PERIOD_RAL; + else if (fRunInfo->GetInstitute()->Contains("triumf")) + beamPeriod = ACCEL_PERIOD_TRIUMF; + else + beamPeriod = 0.0; + + // check if start and end are in proper order + UInt_t start = fRunInfo->GetBkgRange(0); + UInt_t end = fRunInfo->GetBkgRange(1); + if (end < start) { + cout << endl << "PRunSingleHistoRRF::EstimatBkg(): end = " << end << " > start = " << start << "! Will swap them!"; + UInt_t keep = end; + end = start; + start = keep; + } + + // calculate proper background range + if (beamPeriod != 0.0) { + Double_t timeBkg = (Double_t)(end-start)*fTimeResolution; // length of the background intervall in time + UInt_t fullCycles = (UInt_t)(timeBkg/beamPeriod); // how many proton beam cylces can be placed within the proposed background intervall + // correct the end of the background intervall such that the background is as close as possible to a multiple of the proton cylce + end = start + (UInt_t) ((fullCycles*beamPeriod)/fTimeResolution); + cout << endl << "PRunSingleHistoRRF::EstimatBkg(): Background " << start << ", " << end; + if (end == start) + end = fRunInfo->GetBkgRange(1); + } + + // check if start is within histogram bounds + if ((start < 0) || (start >= fForward.size())) { + cerr << endl << ">> PRunSingleHistoRRF::EstimatBkg(): **ERROR** background bin values out of bound!"; + cerr << endl << ">> histo lengths = " << fForward.size(); + cerr << endl << ">> background start = " << start; + cerr << endl; + return false; + } + + // check if end is within histogram bounds + if ((end < 0) || (end >= fForward.size())) { + cerr << endl << ">> PRunSingleHistoRRF::EstimatBkg(): **ERROR** background bin values out of bound!"; + cerr << endl << ">> histo lengths = " << fForward.size(); + cerr << endl << ">> background end = " << end; + cerr << endl; + return false; + } + + // calculate background + Double_t bkg = 0.0; + + // forward + for (UInt_t i=start; i(end - start + 1); + + fBackground = bkg; // keep background (per bin) + + cout << endl << "info> fBackground=" << fBackground << endl; + + fRunInfo->SetBkgEstimated(fBackground, 0); + + return true; +} diff --git a/src/classes/PStartupHandler.cpp b/src/classes/PStartupHandler.cpp index 3f9b425e..47414ed3 100644 --- a/src/classes/PStartupHandler.cpp +++ b/src/classes/PStartupHandler.cpp @@ -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 * diff --git a/src/classes/PTheory.cpp b/src/classes/PTheory.cpp index 08d7f176..efa1c7d8 100644 --- a/src/classes/PTheory.cpp +++ b/src/classes/PTheory.cpp @@ -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 * diff --git a/src/classes/PUserFcn.cpp b/src/classes/PUserFcn.cpp index d0ed32a2..287d63aa 100644 --- a/src/classes/PUserFcn.cpp +++ b/src/classes/PUserFcn.cpp @@ -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 * diff --git a/src/classes/PUserFcnBase.cpp b/src/classes/PUserFcnBase.cpp index 2257613e..d5ec6525 100644 --- a/src/classes/PUserFcnBase.cpp +++ b/src/classes/PUserFcnBase.cpp @@ -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 * diff --git a/src/include/PFitter.h b/src/include/PFitter.h index 8445468c..c9a09014 100644 --- a/src/include/PFitter.h +++ b/src/include/PFitter.h @@ -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 * diff --git a/src/include/PFitterFcn.h b/src/include/PFitterFcn.h index 1c2d846f..94686e00 100644 --- a/src/include/PFitterFcn.h +++ b/src/include/PFitterFcn.h @@ -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 * diff --git a/src/include/PFourier.h b/src/include/PFourier.h index 3e120b7c..bff94450 100644 --- a/src/include/PFourier.h +++ b/src/include/PFourier.h @@ -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 * diff --git a/src/include/PFourierCanvas.h b/src/include/PFourierCanvas.h index 096097c1..1229b34c 100644 --- a/src/include/PFourierCanvas.h +++ b/src/include/PFourierCanvas.h @@ -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 * diff --git a/src/include/PFourierCanvasLinkDef.h b/src/include/PFourierCanvasLinkDef.h index 8ceb2d47..af10deb8 100644 --- a/src/include/PFourierCanvasLinkDef.h +++ b/src/include/PFourierCanvasLinkDef.h @@ -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 * diff --git a/src/include/PFunction.h b/src/include/PFunction.h index b4b60934..9efd48ff 100644 --- a/src/include/PFunction.h +++ b/src/include/PFunction.h @@ -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 * diff --git a/src/include/PFunctionGrammar.h b/src/include/PFunctionGrammar.h index 10c0c9f0..a03bcc4c 100644 --- a/src/include/PFunctionGrammar.h +++ b/src/include/PFunctionGrammar.h @@ -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 * diff --git a/src/include/PFunctionHandler.h b/src/include/PFunctionHandler.h index cd68390e..d8cc78b8 100644 --- a/src/include/PFunctionHandler.h +++ b/src/include/PFunctionHandler.h @@ -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 * diff --git a/src/include/PMsr2Data.h b/src/include/PMsr2Data.h index 292dc8c9..7974d99a 100644 --- a/src/include/PMsr2Data.h +++ b/src/include/PMsr2Data.h @@ -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 * diff --git a/src/include/PMsrHandler.h b/src/include/PMsrHandler.h index fbc6b98b..1dab817b 100644 --- a/src/include/PMsrHandler.h +++ b/src/include/PMsrHandler.h @@ -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); diff --git a/src/include/PMusr.h b/src/include/PMusr.h index 9074a34b..96df487f 100644 --- a/src/include/PMusr.h +++ b/src/include/PMusr.h @@ -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 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 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. diff --git a/src/include/PMusrCanvas.h b/src/include/PMusrCanvas.h index 19241149..735af4f4 100644 --- a/src/include/PMusrCanvas.h +++ b/src/include/PMusrCanvas.h @@ -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 * diff --git a/src/include/PMusrCanvasLinkDef.h b/src/include/PMusrCanvasLinkDef.h index cf6b7dd3..c9c73b15 100644 --- a/src/include/PMusrCanvasLinkDef.h +++ b/src/include/PMusrCanvasLinkDef.h @@ -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 * diff --git a/src/include/PMusrT0.h b/src/include/PMusrT0.h index 91acf97f..b5a47c2f 100644 --- a/src/include/PMusrT0.h +++ b/src/include/PMusrT0.h @@ -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 * diff --git a/src/include/PMusrT0LinkDef.h b/src/include/PMusrT0LinkDef.h index a14248fd..d7e26cc4 100644 --- a/src/include/PMusrT0LinkDef.h +++ b/src/include/PMusrT0LinkDef.h @@ -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 * diff --git a/src/include/PPrepFourier.h b/src/include/PPrepFourier.h index 12f2b5f4..dd40429f 100644 --- a/src/include/PPrepFourier.h +++ b/src/include/PPrepFourier.h @@ -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 * diff --git a/src/include/PRunAsymmetry.h b/src/include/PRunAsymmetry.h index 6d18983b..1bb04701 100644 --- a/src/include/PRunAsymmetry.h +++ b/src/include/PRunAsymmetry.h @@ -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 * diff --git a/src/include/PRunAsymmetryRRF.h b/src/include/PRunAsymmetryRRF.h new file mode 100644 index 00000000..d0aefd7d --- /dev/null +++ b/src/include/PRunAsymmetryRRF.h @@ -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" + +//--------------------------------------------------------------------------- +/** + *

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& par); + virtual Double_t CalcChiSquareExpected(const std::vector& par); + virtual Double_t CalcMaxLikelihood(const std::vector& 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_ diff --git a/src/include/PRunBase.h b/src/include/PRunBase.h index 3e81b2be..37fad18b 100644 --- a/src/include/PRunBase.h +++ b/src/include/PRunBase.h @@ -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 * diff --git a/src/include/PRunDataHandler.h b/src/include/PRunDataHandler.h index 2618b1ef..f75c1428 100644 --- a/src/include/PRunDataHandler.h +++ b/src/include/PRunDataHandler.h @@ -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 * diff --git a/src/include/PRunListCollection.h b/src/include/PRunListCollection.h index 8a741416..00b1112e 100644 --- a/src/include/PRunListCollection.h +++ b/src/include/PRunListCollection.h @@ -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& par) const; + virtual Double_t GetSingleHistoRRFChisq(const std::vector& par) const; virtual Double_t GetAsymmetryChisq(const std::vector& par) const; + virtual Double_t GetAsymmetryRRFChisq(const std::vector& par) const; virtual Double_t GetMuMinusChisq(const std::vector& par) const; virtual Double_t GetNonMusrChisq(const std::vector& par) const; @@ -66,7 +70,9 @@ class PRunListCollection virtual Double_t GetSingleRunChisq(const std::vector& par, const UInt_t idx) const; virtual Double_t GetSingleHistoMaximumLikelihood(const std::vector& par) const; + virtual Double_t GetSingleHistoRRFMaximumLikelihood(const std::vector& par) const; virtual Double_t GetAsymmetryMaximumLikelihood(const std::vector& par) const; + virtual Double_t GetAsymmetryRRFMaximumLikelihood(const std::vector& par) const; virtual Double_t GetMuMinusMaximumLikelihood(const std::vector& par) const; virtual Double_t GetNonMusrMaximumLikelihood(const std::vector& 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 fRunSingleHistoList; ///< stores all processed single histogram data - vector fRunAsymmetryList; ///< stores all processed asymmetry data - vector fRunMuMinusList; ///< stores all processed mu-minus data - vector fRunNonMusrList; ///< stores all processed non-muSR data + vector fRunSingleHistoList; ///< stores all processed single histogram data + vector fRunSingleHistoRRFList; ///< stores all processed single histogram RRF data + vector fRunAsymmetryList; ///< stores all processed asymmetry data + vector fRunAsymmetryRRFList; ///< stores all processed asymmetry RRF data + vector fRunMuMinusList; ///< stores all processed mu-minus data + vector fRunNonMusrList; ///< stores all processed non-muSR data }; #endif // _PRUNLISTCOLLECTION_H_ diff --git a/src/include/PRunMuMinus.h b/src/include/PRunMuMinus.h index 5c13ce68..6c731dd6 100644 --- a/src/include/PRunMuMinus.h +++ b/src/include/PRunMuMinus.h @@ -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 * diff --git a/src/include/PRunNonMusr.h b/src/include/PRunNonMusr.h index 7b4685e5..3bb7ce97 100644 --- a/src/include/PRunNonMusr.h +++ b/src/include/PRunNonMusr.h @@ -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 * diff --git a/src/include/PRunSingleHisto.h b/src/include/PRunSingleHisto.h index 4e9347f3..3980cff1 100644 --- a/src/include/PRunSingleHisto.h +++ b/src/include/PRunSingleHisto.h @@ -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 * diff --git a/src/include/PRunSingleHistoRRF.h b/src/include/PRunSingleHistoRRF.h new file mode 100644 index 00000000..63b063eb --- /dev/null +++ b/src/include/PRunSingleHistoRRF.h @@ -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" + +/** + *

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& par); + virtual Double_t CalcChiSquareExpected(const std::vector& par); + virtual Double_t CalcMaxLikelihood(const std::vector& 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_ diff --git a/src/include/PStartupHandler.h b/src/include/PStartupHandler.h index f91dc679..9125a67b 100644 --- a/src/include/PStartupHandler.h +++ b/src/include/PStartupHandler.h @@ -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 * diff --git a/src/include/PStartupHandlerLinkDef.h b/src/include/PStartupHandlerLinkDef.h index dd099371..036977e7 100644 --- a/src/include/PStartupHandlerLinkDef.h +++ b/src/include/PStartupHandlerLinkDef.h @@ -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 * diff --git a/src/include/PTheory.h b/src/include/PTheory.h index 7baa44bc..85f79ddf 100644 --- a/src/include/PTheory.h +++ b/src/include/PTheory.h @@ -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 * diff --git a/src/include/PUserFcn.h b/src/include/PUserFcn.h index 0c5c1fa6..dacb4041 100644 --- a/src/include/PUserFcn.h +++ b/src/include/PUserFcn.h @@ -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 * diff --git a/src/include/PUserFcnBase.h b/src/include/PUserFcnBase.h index 558c8c41..a920727c 100644 --- a/src/include/PUserFcnBase.h +++ b/src/include/PUserFcnBase.h @@ -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 * diff --git a/src/include/PUserFcnBaseLinkDef.h b/src/include/PUserFcnBaseLinkDef.h index 5f78c800..843c143a 100644 --- a/src/include/PUserFcnBaseLinkDef.h +++ b/src/include/PUserFcnBaseLinkDef.h @@ -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 * diff --git a/src/musrfit.cpp b/src/musrfit.cpp index 060534f0..f5f5ac44 100644 --- a/src/musrfit.cpp +++ b/src/musrfit.cpp @@ -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; iGetSingleHisto(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; iGetSingleHistoRRF(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; iGetAsymmetryRRF(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; iGetSingleHisto(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; iGetSingleHistoRRF(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; iGetAsymmetryRRF(i); + if (data) { + // dump data + musrfit_write_root(f, fln, data, runCounter); + runCounter++; + } + } + } + // muMinus size = runList->GetNoOfMuMinus(); if (size > 0) { diff --git a/src/musrt0.cpp b/src/musrt0.cpp index d1a118c9..272ebe18 100644 --- a/src/musrt0.cpp +++ b/src/musrt0.cpp @@ -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(); diff --git a/src/tests/analyticFakeData/analyticFakeData.C b/src/tests/analyticFakeData/analyticFakeData.C index 2d80c0ac..65a954da 100644 --- a/src/tests/analyticFakeData/analyticFakeData.C +++ b/src/tests/analyticFakeData/analyticFakeData.C @@ -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; iSet("FakeFct/N0", dvec); @@ -202,10 +208,20 @@ void analyticFakeData(const TString type, UInt_t runNo) for (UInt_t i=0; iSet("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); } }